Ecuación del calor con difusión variable

Resumen: En este ejemplo vamos a considerar un problema relativo a la ecuación del calor en el que el coeficiente de conductividad térmica es una matriz no constante. En la discretización en tiempo se consideran los métodos implícitos de Euler y de Gear.

Problema considerado

Sean Ω=[−Lx,Lx]×[−Ly,Ly]$\Omega = [-L_x,L_x]\times [-L_y,L_y]$ y T>0$T>0$. Consideramos el problema de hallar u:Ω×(0,T)→R$u:\Omega\times(0,T) \to \mathbb{R}$ tal que

donde f:Ω→R$f: \Omega \to \mathbb{R}$, K:Ω→R2×2$K:\Omega \to \mathbb{R}^{2\times2}$ y u0:Ω→R$u_0: \Omega\to \mathbb{R}$ son funciones dadas.

Para calcular una solución aproximada, procedemos como sigue:

1. Discretizamos en tiempo tratando la variable t$t$ con el mismo enfoque que para un problema de Cauchy para una EDO: el operador diferencial con las derivadas respecto de t$t$ se aproxima por un esquema en diferencias. Se obtiene así una sucesión de problemas estacionarios en Ω$\Omega$.

2. Cada problema estacionario se discretiza mediante el método de los elementos finitos sobre una triangulación dada, como en los ejemplos de problemas elípticos considerados con anterioridad.

    

Discretización en tiempo Euler implícito

Consideramos una partición uniforme del intervalo [0,T]$[0,T]$ tomando M+1$M+1$ puntos equidistantes 0=t0<⋯<tM=T$0=t_0 <\dots < t_M=T$, es decir dividimos el intervalo [0,T]$[0,T]$ en M$M$ segmentos de igual longitud Δt$\Delta t$:

Para cada n=1,…,M$n = 1,\dots, M$, particularizamos la ecuación diferencial de (1)$\eqref{eq.heat2D.1}$ y las condiciones de contorno para t=tn$t = t_n$:

Denotamos

Para n=0$n = 0$, se tiene

Para n≥1$n\geq 1$, sustituyendo en (2)$\eqref{eq.heat2D.2}$ la derivada de u$u$ respecto de t$t$ por la aproximación en diferencias

se tiene la siguiente ecuación (aproximada)

Para cada n=1,2,…,n$n= 1,2,\dots,n$, tenemos una EDP estacionaria en la incógnita es un$u^n$ ( ya que un−1$u^{n-1}$ ha sido calculada en la etapa anterior).

Por tanto, el esquema discreto en tiempo es como sigue: dado u0(x,y)=u0(x,y)$u^0(x,y) = u_0(x,y)$, para cada n=1,…,M$n=1,\dots, M$, calcular un:Ω→R$u^n:\Omega \to \mathbb{R}$ tal que

Obsérvese que, para cada n$n$, (3)$\eqref{eq.heat2D.3}$ es un problema que no depende del tiempo, que podemos aproximar usando el método de los elementos finitos razonando como en los ejemplos anteriores.

Formulación débil o variacional

El espacio natural para buscar la solución del problema (3)$\eqref{eq.heat2D.3}$ es H01(Ω)$H^1_0(\Omega)$. Para obtener su formulación variacional, multiplicamos en la ecuación ambos miembros por una función v∈H01(Ω)$v\in H_0^1(\Omega)$ e integramos en Ω$\Omega$. Integrando por partes y teniendo en cuenta que v$v$ vale cero sobre ∂Ω$\partial\Omega$, obtenemos:

Obsérvese que todos los problemas (PV)$\eqref{eq.heat2D.PV}$, para los distintos valores de n=1,2,…,M$n=1,2,\dots,M$, comparten la misma forma bilineal

y que solo difieren en la forma lineal. En la resolución numérica, esto se traduce en el hecho de que todos los sistemas lineales a resolver tienen la misma matriz, que, por lo tanto, solo habrá que construir y factorizar la primera vez.

Código

Resolvemos el problema para: Lx=Ly=10$L_x = L_y=10$, T=5$T = 5$, u0=μ0e−(x2+y2)/2σ0$u_0 = \mu_0 \, e^{-(x^2+y^2)/2\sigma_0}$,

x
//
//  Problem data
//
real L = 10., T = 3.;
real f0 = 1., x1 = -4., x2 = -3., y1 = 2.5, y2 = 3.5;
real mu0 = 1., sigma0 = 0.5;
real sigma = -0.5/sigma0;

func K11 = 1.+x*x;
func K22 = 1.+y*y;
func K12 = 0.1*x*y;
func f = f0*(x1 < x)*(x < x2)*(y1 < y)*(y < y2);
func u0 = mu0 * exp(sigma*(x*x+y*y));

//
//  Discretization data
//
int  nt = 50, nx = 50, ny = 50;
real dt = T/nt, undt = 1./dt;

//
//  Mesh
//
mesh Th = square(nx, ny, [L*(-1+2*x), L*(-1+2*y)]);
plot(Th, wait=true);

//
//  FE space:
//
fespace Vh(Th,P1);
Vh u, v;
Vh fh = f, uold;
plot(fh, fill = true, value = true, 
         cmm = "PDE right hand side", wait = true);

//
//  Problem
//  Option "solver" allow to choose the method to use
problem Heat2DEuler(u , v, solver = Cholesky, init = 1) =
        int2d(Th)(undt * u*v + K11*dx(u)*dx(v) + K22*dy(u)*dy(v)
                             + K12*(dx(u)*dy(v) + dy(u)*dx(v)))
      - int2d(Th)(fh*v) - int2d(Th)(undt*uold*v)
      + on(1,2,3,4, u = 0.);

//
//  Computations and plots
//
u = u0;
plot(u, nbiso = 40, dim = 3, fill=true, 
        cmm = "Initial u, t = 0", wait = true);
real t;
for (int k = 1; k < nt+1; k ++){
    uold = u;
    t = k*dt;
    Heat2DEuler;
    plot(u, nbiso = 40, dim = 3, fill=true, value = true,
          cmm = "Solution al t = " + t, wait=true);
}
//  Use the "prev=true" option in plots with t > 0 
//  to preserve the aspect ratio of previous plot

Distretización en tiempo Gear implícito

Consideramos ahora una discretización en tiempo usando el esquema implícito de Gear de dos pasos:

Sustituyendo en la ecuación diferencial de (2)$\eqref{eq.heat2D.2}$, obtenemos, para n≥2$n\geq2$

Para cada n≥2$n\geq 2$, tenemos una EDP estacionaria en la incógnita es un$u^n$ (un−1$u^{n-1}$ y un−2$u^{n-2}$ han debido ser calculadas en las etapas anteriores). Obviamente, esto no sirve para calcular u1$u^1$, que deberá ser calculada a partir de u0$u^0$ utilizando un esquema de un solo paso, por ejemplo el de Euler.

Por tanto, el algoritmo completo es como sigue:

1. Tomar u0(x,y)=u(x,y,0)$u^0(x,y) = u(x,y,0)$,
2. Calcular u1(x,y)≈u(x,y,t1)$u^1(x,y) \approx u(x,y,t_1)$ solución de

3. Para cada n=2,…,M$n=2,\dots, M$, calcular un≈u(x,y,tn)$u^n \approx u(x,y,t_n)$ solución de

Código

xxxxxxxxxx
//
//  Data
//
real L = 10., T = 3.;
real f0 = 1., x1 = -4., x2 = -3., y1 = 2.5, y2 = 3.5;
real mu0 = 1., sigma0 = 0.5;
real sigma = -0.5/sigma0;

func K11 = 1.+x*x;
func K22 = 1.+y*y;
func K12 = 0.1*x*y;
func f = f0*(x1 < x)*(x < x2)*(y1 < y)*(y < y2);
func u0 = mu0 * exp(sigma*(x*x+y*y));
//  Discretization data
int  nt = 50, nx = 50, ny = 50;
real dt = T/nt;
real dt1 = 1./dt, dt2 = 0.5*dt1, dt3 = 3.*dt2, dt4 = 2*dt1;

//
//  Mesh
//
mesh Th = square(nx, ny, [L*(-1+2*x), L*(-1+2*y)]);
plot(Th, wait=true);

//
//  FE space:
//
fespace Vh(Th,P1);
Vh u, v;
Vh fh = f, uold1, uold2;
// uold1  -->  u_{i-1},  uold2 -->  u_{i-2}
plot(fh, fill = true, value = true, cmm = "PDE right hand side", wait = true);

//
//  Problem for the first step: Euler squeme
//
uold1 = u0;
plot(uold1, nbiso = 40, dim = 3, fill=true, cmm = "Initial u, t = 0", wait = true);
solve Euler(u , v, solver = Cholesky) =
        int2d(Th)(dt1 * u*v + K11*dx(u)*dx(v) + K22*dy(u)*dy(v)
                             + K12*(dx(u)*dy(v) + dy(u)*dx(v)))
      - int2d(Th)(fh*v) - int2d(Th)(dt1*uold1*v)
      + on(1,2,3,4, u = 0.);
plot(u, nbiso = 40, dim = 3, fill=true, cmm = "Solution u, t = "+ dt, wait = true);

//
//  Problem for Gear steps
problem Heat2DGear(u , v, solver = Cholesky, init = 1) =
      int2d(Th)(dt3 * u*v + K11*dx(u)*dx(v) + K22*dy(u)*dy(v)
                           + K12*(dx(u)*dy(v) + dy(u)*dx(v)))
    - int2d(Th)(fh*v) - int2d(Th)(dt4*uold1*v) + int2d(Th)(dt2*uold2*v)
    + on(1,2,3,4, u = 0.);

//
//  Computations of Gear steps
//
real t;
for (int k = 2; k < nt+1; k ++){
    uold2 = uold1; 
    uold1 = u;
    t = k*dt;
    Heat2DGear;
    plot(u, nbiso=40, dim = 3, fill=true, value = true, prev = true,
            cmm = "Solution al t = " + t, wait=true);
}

Ejercicios

Resolver el problema (1)$\eqref{eq.heat2D.1}$, mediante los dos esquemas descritos antes, cambiando las condiciones de contorno por las siguientes:

1. Condición de Neumann homogénea sobre toda la frontera del dominio: (K∇u)⋅n=0$(K\nabla u)\cdot n = 0$ sobre ∂Ω×(0,T)$\partial\Omega\times (0,T)$.

2. Condición de Fourier sobre toda la frontera: (K∇u)⋅n=pu$(K\nabla u)\cdot n = p u$ sobre ∂Ω×(0,T)$\partial\Omega\times (0,T)$, con p=2$p = 2$.

3. Condiciones mixtas Dirichlet-Neumann, tomando ∂Ω=Γ1∪Γ2∪Γ3∪Γ4$\partial\Omega = \Gamma_1\cup\Gamma_2\cup\Gamma_3\cup\Gamma_4$:

   {u=0sobre Γ2×(0,T),(K∇u)⋅n=0 sobre Γ1∪Γ3∪Γ4×(0,T).

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