FreeFEM Documentation on GitHub

stars - forks

Thermal Conduction

Summary : Here we shall learn how to deal with a time dependent parabolic problem. We shall also show how to treat an axisymmetric problem and show also how to deal with a nonlinear problem

How air cools a plate

We seek the temperature distribution in a plate \((0,Lx)\times(0,Ly)\times(0,Lz)\) of rectangular cross section \(\Omega=(0,6)\times(0,1)\); the plate is surrounded by air at temperature \(u_e\) and initially at temperature \(u=u_0+\frac x L u_1\). In the plane perpendicular to the plate at \(z=Lz/2\), the temperature varies little with the coordinate \(z\); as a first approximation the problem is 2D.

We must solve the temperature equation in \(\Omega\) in a time interval (0,T).

\[\begin{split}\begin{array}{rcl} \partial_t u -\nabla\cdot(\kappa\nabla u) &= 0 & \hbox{ in } \Omega\times(0,T)\\ u(x,y,0) &= u_0+x u_1 &\\ \kappa\frac{\partial u}{\partial n} +\alpha(u-u_e) &= 0 & \hbox{ on } \Gamma\times(0,T) \end{array}\end{split}\]

Here the diffusion \(\kappa\) will take two values, one below the middle horizontal line and ten times less above, so as to simulate a thermostat.

The term \(\alpha(u-u_e)\) accounts for the loss of temperature by convection in air. Mathematically this boundary condition is of Fourier (or Robin, or mixed) type.

The variational formulation is in \(L^2(0,T;H^1(\Omega))\); in loose terms and after applying an implicit Euler finite difference approximation in time; we shall seek \(u^n(x,y)\) satisfying for all \(w\in H^1(\Omega)\):

\[\int_\Omega(\frac{u^n-u^{n-1}}{\delta t} w + \kappa\nabla u^n\nabla w) +\int_\Gamma\alpha(u^n-u_ue)w=0\]
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
// Parameters
func u0 = 10. + 90.*x/6.;
func k = 1.8*(y<0.5) + 0.2;
real ue = 25.;
real alpha=0.25;
real T=5.;
real dt=0.1 ;

// Mesh
mesh Th = square(30, 5, [6.*x,y]);

// Fespace
fespace Vh(Th, P1);
Vh u=u0, v, uold;

// Problem
problem thermic(u, v)
    = int2d(Th)(
          u*v/dt
        + k*(
              dx(u) * dx(v)
            + dy(u) * dy(v)
        )
    )
    + int1d(Th, 1, 3)(
          alpha*u*v
    )
    - int1d(Th, 1, 3)(
          alpha*ue*v
    )
    - int2d(Th)(
          uold*v/dt
    )
    + on(2, 4, u=u0)
    ;

// Time iterations
ofstream ff("thermic.dat");
for(real t = 0; t < T; t += dt){
    uold = u; //equivalent to u^{n-1} = u^n
    thermic; //here the thermic problem is solved
    ff << u(3., 0.5) << endl;
    plot(u);
}

Note

We must separate by hand the bilinear part from the linear one.

Note

The way we store the temperature at point (3, 0.5) for all times in file thermic.dat. Should a one dimensional plot be required, the same procedure can be used. For instance to print \(x\mapsto \frac{\partial u}{\partial y}(x,0.9)\) one would do:

1
2
for(int i = 0; i < 20; i++)
   cout << dy(u)(6.0*i/20.0,0.9) << endl;

Results are shown on Fig. 13 and Fig. 14.

Thermic

Fig. 13 Temperature at \(t=4.9\).

ThermicVST

Fig. 14 Decay of temperature versus time at \(x=3, y=0.5\)

Thermal conduction

Axisymmetry: 3D Rod with circular section

Let us now deal with a cylindrical rod instead of a flat plate. For simplicity we take \(\kappa=1\).

In cylindrical coordinates, the Laplace operator becomes (\(r\) is the distance to the axis, \(z\) is the distance along the axis, \(\theta\) polar angle in a fixed plane perpendicular to the axis):

\[\Delta u = {1\over r}\partial _r(r\partial _r u) + {1\over r^2}\partial ^2_{\theta\theta} u + \partial ^2_{z z}.\]

Symmetry implies that we loose the dependence with respect to \(\theta\); so the domain \(\Omega\) is again a rectangle \(]0,R[\times]0,|[\) . We take the convention of numbering of the edges as in square() (1 for the bottom horizontal …); the problem is now:

\[\begin{split}\begin{array}{rcl} r\partial_t u-\partial _r(r\partial _r u) - \partial _z(r\partial _z u) &= 0 &\hbox{ in } \Omega\\ u(t=0) &= u_0 + \frac z{L_z} (u_1-u)&\\ u|_{\Gamma_4} &= u_0&\\ u|_{\Gamma_2} &= u_1&\\ \alpha(u-u_e) + {\partial u\over \partial n} |_{\Gamma_1\cup\Gamma_3} &= 0& \end{array}\end{split}\]

Note that the PDE has been multiplied by \(r\).

After discretization in time with an implicit scheme, with time steps dt, in the FreeFEM syntax \(r\) becomes \(x\) and \(z\) becomes \(y\) and the problem is:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
problem thermaxi(u, v)
    = int2d(Th)(
          (u*v/dt + dx(u)*dx(v) + dy(u)*dy(v))*x
    )
    + int1d(Th, 3)(
          alpha*x*u*v
    )
    - int1d(Th, 3)(
          alpha*x*ue*v
    )
    - int2d(Th)(
          uold*v*x/dt
    )
    + on(2, 4, u=u0);

Note

The bilinear form degenerates at \(x=0\). Still one can prove existence and uniqueness for \(u\) and because of this degeneracy no boundary conditions need to be imposed on \(\Gamma_1\).

A Nonlinear Problem : Radiation

Heat loss through radiation is a loss proportional to the absolute temperature to the fourth power (Stefan’s Law). This adds to the loss by convection and gives the following boundary condition:

\[\kappa{\partial u\over \partial n} +\alpha(u-u_e) + c[(u + 273)^4 - (u_e+273)^4] = 0\]

The problem is nonlinear, and must be solved iteratively. If \(m\) denotes the iteration index, a semi-linearization of the radiation condition gives

\[{\partial u^{m+1}\over \partial n} + \alpha(u^{m+1}-u_e)+ c(u^{m+1}-u_e) (u^m+u_e +546) ((u^m + 273)^2 + (u_e+273)^2) = 0,\]

because we have the identity \(a^4 - b^4 = (a-b)(a+b)(a^2+b^2)\).

The iterative process will work with \(v=u-u_e\).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
...
// Parameters
real rad=1e-8;
real uek=ue+273;

// Mesh
fespace Vh(Th, P1);
Vh vold, w, v=u0-ue, b;

// Problem
problem thermradia(v, w)
    = int2d(Th)(
          v*w/dt
        + k*(dx(v) * dx(w) + dy(v) * dy(w))
    )
    + int1d(Th, 1, 3)(
          b*v*w
    )
    - int2d(Th)(
          vold*w/dt
    )
    + on(2, 4, v=u0-ue)
    ;

for (real t=0;t<T;t+=dt){
    vold = v;
    for (int m = 0; m < 5; m++){
        b = alpha + rad * (v + 2*uek) * ((v+uek)^2 + uek^2);
        thermradia;
    }
}
vold = v+ue;

// Plot
plot(vold);
Table of content