# 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).

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):

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):

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:

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:

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<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); |