Transferencia del calor

Resumen: En este ejemplo vamos a considerar la ecuación del calor estacionaria escrita en forma de divergencia con un coeficiente de difusión no constante y definido a trozos. Veremos también cómo exportar la triangulación a un archivo externo.

Problema considerado

$\Omega\subset \mathbb{R}^2$ $\Omega = C_0\setminus C_1$ $\partial \Omega = \partial C_0\cup \partial C_1$ $u:\Omega \to \mathbb{R}$ tal que

\begin{matrix} (P) & {\begin{cases} - \nabla \cdot (κ \nabla u) = 0 & en Ω, \\ u = 20 & sobre \partial C_{0}, \\ u = 60 & sobre \partial C_{1}, \end{cases} \end{matrix}

donde

\begin{matrix} (1) & \begin{matrix} κ = {\begin{cases} 1 & en Ω ∖ C_{2}, \\ 3 & en C_{2} . \end{cases} \end{matrix} \end{matrix}

$C_1$ $^\circ$ $C_0$ $^\circ$ $\kappa$ $C_2$ $\Omega$ .

Formulación variacional

Observación: todos los cálculos que se hacen en esta sección son puramente formales, es decir, se supone siempre la regularidad suficiente para que todos ellos estén justificados.

$\eqref{eq.heatex}$ $\widetilde{g}+H_0^1(\Omega)$ $\widetilde{g}\in H^1(\Omega)$ $\widetilde{g}|_{\partial C_0} = 20$ $\widetilde{g}|_{\partial C_1} = 60$ .

$\eqref{eq.heatex}$ $v\in H_0^1(\Omega)$ $\Omega$ $v$ $\partial\Omega$ , se obtiene:

\begin{matrix} \int_{Ω} κ \nabla u \cdot \nabla v d x d y = \int_{Ω} f v d x d y \forall v \in H_{0}^{1} (Ω) . \end{matrix}

$\eqref{eq.heatex}$ $\nabla u \cdot \nabla v$ )

\begin{matrix} (PV) & {\begin{cases} Hallar u \in \tilde{g} + H_{0}^{1} (Ω) tal que \\ \int_{Ω} κ (\partial_{x} u \partial_{x} v + \partial_{y} u \partial_{y} v) d x d y = \int_{Ω} f v d x d y \forall v \in H_{0}^{1} (Ω) . \end{cases} \end{matrix}

Problema aproximado

$P_1$ $\partial\Omega$ $\mathcal{T}_h$ $\overline{\Omega}$ . Consideramos

\begin{matrix} W_{h} = {v_{h} | v_{h} \in C^{0} (\overset{―}{Ω}), v_{h} |_{T} \in P_{1} (T) \forall T \in T_{h}} \subset H^{1} (Ω) \end{matrix}

$\text{dim}\, V_h = M$ $M$ $\mathcal{T}_h$ . Consideramos también

\begin{matrix} V_{h} = {v_{h} | v_{h} \in V_{h}, v_{h} |_{\partial Ω} = 0} \subset H_{0}^{1} (Ω), \end{matrix}

$\mathcal{T}_h$ no $\partial\Omega$ $\eqref{eq.heatechangePV}$ es

\begin{matrix} (PV_{h}) & {\begin{cases} Hallar u_{h} \in {\tilde{g}}_{h} + V_{h} tal que \\ \int_{Ω} κ (\partial_{x} u_{h} \partial_{x} v_{h} + \partial_{y} u_{h} \partial_{y} v_{h}) d x d y = \int_{Ω} f v_{h} d x d y \forall v_{h} \in V_{h}, \end{cases} \end{matrix}

$\widetilde{g}_h\in W_h$ $\widetilde{g}_h(a_i)= \widetilde{g}(a_i)$ $a_i\in \partial\Omega$ $\Omega$ .

Resolución con FreeFEM

$C_0$ $C_1$ $[1,2]\times[-3,3]$ $C_2$ $[-2, -1]\times[-3, 3]$ .

Definición de la frontera del dominio

$C_1$ $C_2$ label = 6 $C_1$ label = 2 $C_2$ label = 90 $C_0$ . Las siguientes instrucciones definen la frontera del dominio.


// Datos
real xa = 1., yb = 3.;
real r = 5.;
int CC0 = 90, CC1 = 6, CC2 = 2;

// Definicion del borde
// Frontera exterior del dominio
border C0(t = 0, 2*pi){x = r*cos(t); y = r*sin(t); label = CC0;};
// Rectangulo C1 
border C11(t=0., 1.){x = xa+t; y = yb; label = CC1;};
border C12(t=0., 1.){x = xa+1; y = yb*(1-2*t); label = CC1;};
border C13(t=0., 1.){x = xa-t+1; y = -yb; label = CC1;};
border C14(t=0., 1.){x = xa; y = yb*(2*t-1); label = CC1;};
// Rectangulo C2 
border C21(t=0., 1.){x = -xa-(1-t); y = yb; label = CC2;};
border C22(t=0., 1.){x = -xa; y = yb*(1-2*t); label = CC2;};
border C23(t=0., 1.){x = -xa-t; y = -yb; label = CC2;};
border C24(t=0., 1.){x = -xa-1; y = yb*(2*t-1); label = CC2;};

Construcción del mallado

mesh $C_2$ $C_2$ $\kappa$ (se verá más adelante).


// Construccion de la triangulacion
mesh Th;
Th = buildmesh( C0(100) + C11(5) + C12(30) + C13(5) + C14(30) + 
                         C21(-5) + C22(-30) + C23(-5) + C24(-30) );
plot( Th , wait = 1);

Definición del espacio de elementos finitos

$\mathcal{T}_h$ $W_h$ $\overline{\Omega}$ $T\in\mathcal{T}_h$ es un polinomio de grado menor o igual que 1.


// Definicion del espacio de EF
fespace Wh(Th, P1);

Definición del problema variacional

$u=20$ $\partial C_0$ $u=60$ $\partial C_1$ bloqueo $W_h$ .


// Declaracion de las variables del problema (elementos de Wh)
Wh u, v;

$\kappa$ $\eqref{kappa}$ $\kappa$ $P_0$ $C_2$ W0 $\kappa$ es un elemento de este espacio y lo representamos gráficamente para comprobar que toma valores constantes.


// Definicion del coeficiente kappa
//  kappa = 1 en Omega \setminus C2
//  kappa = 3 en C2
fespace W0(Th, P0);
W0 kappa = 1 + 2*(x<-1)*(x>-2)*(y<3)*(y>-3);
plot(kappa, wait = 1, fill = true, ps = "kappa.eps");

Definimos de forma habitual el problema variacional usando el tipo problem


// Definicion del problema variacional
problem HeatEx(u, v) =
    int2d(Th)(kappa*(dx(u)*dx(v)+dy(u)*dy(v)))
    + on(CC0, u = 20)
    + on(CC1, u = 60);

Resolución del problema y representación


// Resolucion
HeatEx;

plot(u, C0(100) + C11(5) + C12(30) + C13(5) + C14(30) +
          C21(-5) + C22(-30) + C23(-5) + C24(-30), wait = 1);

// Dibujo de la solucion
plot(u, wait = 1);

Código completo del programa


// Datos
real xa = 1., yb = 3.;
real r = 5.;
int CC0 = 90, CC1 = 6, CC2 = 2;

real[int] colorhsv=[  // modelo de color hsv 
0.6,  0.9,  1.0,
1.0,  0.9,  1.0  ];

// Definicion del borde
// Circunferencia C0
border C0(t = 0, 2*pi){x = r*cos(t); y = r*sin(t); label = CC0;};
// Rectangulo C1 
border C11(t=0., 1.){x = xa+t; y = yb; label = CC1;};
border C12(t=0., 1.){x = xa+1; y = yb*(1-2*t); label = CC1;};
border C13(t=0., 1.){x = xa-t+1; y = -yb; label = CC1;};
border C14(t=0., 1.){x = xa; y = yb*(2*t-1); label = CC1;};
// Rectangulo C2 
border C21(t=0., 1.){x = -xa-(1-t); y = yb; label = CC2;};
border C22(t=0., 1.){x = -xa; y = yb*(1-2*t); label = CC2;};
border C23(t=0., 1.){x = -xa-t; y = -yb; label = CC2;};
border C24(t=0., 1.){x = -xa-1; y = yb*(2*t-1); label = CC2;};

// Triangulacion
mesh Th;
Th = buildmesh(C0(100) + C11(5) + C12(30) + C13(5) + C14(30) + 
                         C21(-5) + C22(-30) + C23(-5) + C24(-30) );
plot(Th, C0(100) + C11(5) + C12(30) + C13(5) + C14(30) +
          C21(-5) + C22(-30) + C23(-5) + C24(-30), wait = 1);

// Espacio de elementos finitos: P1
fespace Wh(Th, P1);
Wh u, v;

// Definicion del coeficiente kappa
//  kappa = 1 en Omega \setminus C2
//  kappa = 3 en C2
fespace W0(Th, P0);
W0 kappa = 1 + 2*(x<-1)*(x>-2)*(y<3)*(y>-3);
plot(kappa, wait = 1, fill = true, ps = "kappa.eps");

// Definicion del problema
problem HeatEx(u, v) =
    int2d(Th)(kappa*(dx(u)*dx(v)+dy(u)*dy(v)))
    + on(CC0, u = 20)
    + on(CC1, u = 60);

// Resolucion
HeatEx;

// Representacion
plot(u, wait = 1, fill = true, hsv = colorhsv);

Complementos

Se puede salvar la triangulación para usarla en otro programa usando el comando savemesh. Se guardan los números de las etiquetas de las fronteras, pero no se guardan los nombres de las variables que hemos usado para definirlas (en este ejemplo CC0 , CC1 y CC2). Para usarlas posteriormente, hay que utilizar expresamente sus valores.


savemesh(Th, "triangulacion.msh");

Más concretamente, para resolver el problema de este ejemplo, una vez guardada la triangulación en un fichero, hay que leerla usando el comando readmesh y la resolución del problema se haría de la forma siguiente:


// Lectura de la triangulacion
mesh Th = readmesh("triangulacion.msh");

// Espacio de elementos finitos: P1
fespace Vh(Th, P1);
Vh u, v;

// Definicion del coeficiente kappa
//  kappa = 1 en Omega \setminus C2
//  kappa = 3 en C2
fespace W0(Th, P0);
W0 kappa = 1 + 2*(x<-1)*(x>-2)*(y<3)*(y>-3);

// Definicion del problema
problem HeatEx(u, v) =
    int2d(Th)(kappa*(dx(u)*dx(v)+dy(u)*dy(v)))
    + on(90, u = 20)
    + on(6,  u = 60);
    
// Resolucion y representacion grafica
HeatEx;
plot(u, wait = 1);

Ejercicio

$\eqref{eq.heatex}$ realizando la resolución manual con varf.

Anna Doubova - Rosa Echevarría - Dpto. EDAN - Universidad de Sevilla