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{eqnarray} && \p_t u -\n\cdot(\kappa\n u)=0 \hbox{ in } \Omega\times(0,T), \cr&& u(x,y,0)=u_0+x u_1 \cr&& \kappa\frac{\p u}{\p n} +\alpha(u-u_e)=0\hbox{ on } \Gamma\times(0,T). \end{eqnarray}

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\n u^n\n 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{\p u}{\p y}(x,0.9)$ one would do

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

Results are shown on figure 1.

Fig. 1: Temperature at $T=4.9$ Decay of temperature versus time at $x=3, y=0.5$

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}\p _r(r\p _r u) + {1\over r^2}\p ^2_{\theta\theta} u + \p ^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{eqnarray}&& r\p_t u-\p _r(r\p _r u) - \p _z(r\p _z u) = 0 \hbox{ in } \Omega, \cr&& u(t=0) = u_0 + \frac z{L_z} (u_1-u) \cr&& u|_{\Gamma_4} = u_0,\quad u|_{\Gamma_2} = u_1, \quad \alpha(u-u_e) + {\p u\over \p n} |_{\Gamma_1\cup\Gamma_3} = 0. \end{eqnarray}

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$.

The problem is nonlinear, and must be solved iteratively. If $m$ denotes the iteration index, a semi-linearization of the radiation condition gives
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