Test case of Transient Heat Equation in Three Dimensions

1. Introduction

This page presents the simulation of temperature on geometry of torus in transient case with electric source term computed with CFPDEs application.

The electric source term isn’t the real physics source term because it varies over time.

2. Run the calculation

The command line to run this case is :

    mpirun -np 32 feelpp_toolbox_coefficientformpdes --config-file=heat.cfg --cfpdes.gmsh.hsize=1e-3

3. Data Files

The case data files are available in Github here :

CFG file - Edit the file
JSON file - Edit the file
GEO file - Edit the file

4. Equation

We solve the heat equation (Static Heat) in the stationary mode.

Heat equation

\[\text{(Heat)} \left\{ \begin{matrix} \rho C_p \, \frac{\partial T}{\partial t} - k \Delta T = Q \text{ on } \Omega_c \\ \frac{\partial T}{\partial \mathbf{n}} = 0 \text{ on } \Gamma_{Nc} \\ k \, \frac{\partial T}{\partial \mathbf{n}} = h \, \left( T - T_c \right) \text{ on } \Gamma_{Rc} \\ \end{matrix} \right.\]

5. Geometry

The geometry is a torus in cartesian coordinates \((x,y,z)\).

Geometry in three dimensions

The geometrical domains are :

Conductor : the torus, composed by conductive material
Interior and Exterior : interior and exterior of ring (or right and left of rectangle in axisymmetric coordinates \((r,z)\)) correspond to \(\Gamma_R\)
Bottom and Upper : up and bottom of the ring (or up and bottom of the rectangle in axisymmetric coordinates \((r,z)\)) correspond to \(\Gamma_N\)

Symbol

Description

value

unit

\(r_{int}\)

interior radius of torus

\(75e-3\)

\(r_{ext}\)

exterior radius of torus

\(100.2e-3\)

\(z_1\)

half-height of torus

\(25e-3\)

6. Boundary Conditions

We impose the boundary conditions :

Neumann : \(\frac{\partial T}{\partial \mathbf{n}} = 0\) on Interior and Exterior
Robin : \(-k \, \frac{\partial T}{\partial \mathbf{n}} = h \, \left( T - T_c \right)\) on Upper and Bottom

We initialize on Conductor : \(T(t=0,r,z)=T_i\).

On JSON file, the boundary conditions are writed :

Boundary conditions on JSON file

    "BoundaryConditions":
    {
        "heat":
        {
            "Robin":
            {
                "Interior":
                {
                    "expr1":"h:h",
                    "expr2":"h*T_c:h:T_c"
                },
                "Exterior":
                {
                    "expr1":"h:h",
                    "expr2":"h*T_c:h:T_c"
                }
            }
        }
    }

On JSON file, the intial conditions are writed : .Initial conditions on JSON file

    "InitialConditions":
    {
        "temperature":
        {
            "Expression":
            {
                "myic":
                {
                    "markers":"Conductor",
                    "expr":"T_i:T_i"
                }
            }
        }
    }

7. Weak Formulation

We obtain :

Weak Formulation of Heat Equation

\[ \text{(Weak Heat)} \int_{\Omega_c}{ \rho C_p \, \frac{\partial T}{\partial t} \, \phi } + \int_{\Omega_c}{ k \, \nabla T \cdot \nabla \phi } = \int_{\Omega_c}{ Q \, \phi } - \int_{\Gamma_{Rc}}{ h \, (T-T_c) \, \phi }\]

8. Parameters

The parameters of problem are :

On Conductor :

Symbol

Description

Value

Unit

\(Q\)

source term

\(\sigma \, \left( \frac{U}{2\pi \, r} \right)\)

\(U\)

electrical potential

\(1\)

\(Volt\)

\(\sigma\)

electrical conductivity

\(58e6\)

\(S/m\)

\(k\)

thermal conductivity

\(380\)

\(W/m/K\)

\(h\)

convective coefficient

\(8e4\)

\(W/m^2/K\)

\(T_c\)

cooling temperature

\(293\)

\(K\)

\(T_i\)

initial temperature

\(293\)

\(K\)

\(\rho\)

density

\(10000\)

\(kg/m^3\)

\(C_p\)

thermal capacity

\(380\)

\(J/K/kg\)

On JSON file, the parmeters are writed :

Parameters on JSON file

    "Parameters":
    {
	    "h":80000,  // W/m2/K
        "T_c":293,  // K
        "T_i":293,  // K
        "U":1,      // V

        // Constants of analytical solve
        "a":1933.10,    // K
        "b":0.40041,    // K
        "rmax":0.0861910719118454,  // m
        "Tmax":364.446  // K
    }

9. Coefficient Form PDEs

We use the application Coefficient Form PDEs. The coefficient associate to Weak Formulation are :

On Conductor :

Coefficient

Description

Expression

\(d\)

damping or mass coefficient

\(\rho \, C_p\)

\(c\)

diffusion coefficient

\(k\)

\(f\)

source term

\(\sigma \left( \frac{U}{2\pi \, r} \right)^2 = \sigma \left( \frac{U}{2\pi \, \sqrt{x^2+y^2}} \right)^2\)

On JSON file, the coefficients are writed :

CFPDEs coefficients on JSON file

    "Materials":
    {
        "Conductor":
        {
            "k":380,        // W/m/K
            "sigma":58e+6,  // S.m-1
            "rho":10000,    // kg/m3
            "Cp":380,       // J/K/kg

            "heat_c":"k:k",
            "heat_f":"sigma*(U/2/pi)*(U/2/pi)/(x*x+y*y):sigma:U:x:y",
            "heat_d":"rho*Cp:rho:Cp"
        }
    }

10. Numeric Parameters

This section show the parameters used to compute the simulation.

Size of mesh : \(1 \, mm\)
Time Parameters :
- Time step : \(0.1 \, s\)
- Initial Time : \(0 \, s\)
- Final Time : \(15 \, s\)
Element type : \(P1\)
Solver : automatic
Number of CPU core : \(32\)

Mesh of Geometry (size of mesh \(10 \, mm\))

11. Results

We can vizualize the result on geometry :

Temperature \(T(K)\)

I plot the maximal temperature \(T_{max}\) across time :

\(T_{max}\) \((K)\) in time

The maximum temperature increase in time and stabilized around \(T_{max}=364 \, K\).

I plot the temperature at \(t=15s\), at this time, we consider the system stationary.

\(T\) \((K)\) on \(Or\) axis at \(t=15s\)

The temperature is constant by \(z\). We can see with the cut, the temperature have a profil of bell. The temperature is less at the cooled interior and exterior than the middle of torus.

12. Validation

12.1. Exact solution

In this section, we compute the exact solution of problem.

To simplify the problem, we pass on cylindric coordinates.

The problem is independant of \(\theta\) and \(z\), thus, the equation (Heat Axis) becomes :

\[ - \frac{1}{r} \frac{\partial \left( r \frac{\partial T}{\partial r} \right)}{\partial r} = \sigma \left( \frac{u}{2\pi \, r} \right)\]

Thus :

\[ T = a \, ln(r) - \frac{\sigma}{2k} \left( \frac{U}{2\pi} \right)^2 \, ln^2(r) + b\]

With \(a\) and \(b\) integrate constants.

With the Robin boundary conditions :

\[ -k \frac{\partial T}{\partial \mathbf{n}} = h(r) \left( T(r) - T_c(r) \right)\]

So :

\[\left\{ \begin{eqnarray*} \left( k+h_{ext} r_{ext} ln(r_{ext}) \right) \, a + h_{ext} r_{ext} b = \frac{\sigma}{2k} \left( \frac{U}{2\pi} \right)^2 ln(r_{ext}) \left( h_{ext} r_{ext} ln(r_{ext})+2k \right) + h_{ext} r_{ext} T_{c \, ext} \\ \left( -k+h_{int} r_{int} ln(r_{int}) \right) \, a + h_{int} r_{int} b = \frac{\sigma}{2k} \left( \frac{U}{2\pi} \right)^2 ln(r_{int}) \left( h_{int} r_{int} ln(r_{int})-2k \right) + h_{int} r_{int} T_{c \, int} \end{eqnarray*} \right.\]

With \(T_{c \, int}=T(r_{int})\), \(T_{c \, ext}=T(r_{ext})\), \(h_{int}=h(r_{int})\) and \(h_{ext}=h(r_{ext})\).

And, with our conditions :

\[\left\{ \begin{eqnarray*} a &=& \frac{\sigma}{2k} \left( \frac{U}{2 \pi} \right)^2 \\ b &=& k \left( \frac{1}{h_{int} r_{int}} + \frac{1}{h_{ext} r_{ext}} \right) + ln(\frac{r_{ext}}{r_{int}}) \end{eqnarray*} \right.\]

The temperature can be right like :

\[ T(r) = - a ln^2(\frac{r}{r_{max}}) + T_{max}\]

With \(T_{max} = T(r_{max})\) the maximal temperature and \(r_{max}\) its input.

And :

\[ r_{max} = e^{\frac{T_{c \, int}-T_{c \, ext} + ac}{2 \, a \, b}}\]

\[ T_{max} = \frac{2 \, a \, k}{h_{int} r_{int} + h_{ext} r_{ext}} \, ln(h_{int} r_{int} + h_{ext} r_{ext}) + \frac{h_{int} r_{int} T_{c \, int} + h_{ext} r_{ext} T_{c \, ext}}{h_{int} r_{int} + h_{ext} r_{ext}} + a \frac{ h_{int} r_{int} ln^2(\frac{r_{int}}{r_{max}}) + h_{ext} r_{ext} ln^2(\frac{r_{ext}}{r_{max}}) }{h_{int} r_{int} + h_{ext} r_{ext}}\]

12.2. Convergence Test

I compute the L2-error between the numeric result at final time \(t=t_f\) and exact solution of static problem (\(||T_{num}(t=t_f)-T_{exact}(t=t_f)||_{L^2(\Omega)}\)) for different mesh sizes.

hsize (m)	\(\Vert T_{num}(t=t_f) - T_{exact}(t=t_f) \Vert_{L^2(\Omega)} (K)\)
0.1	7.1966488756673541e-01
0.0316	6.8845957529152346e-01
0.01	1.6586735241573908e-01
0.00316	1.9206462876028615e-02
0.001	1.9189853233305876e-03

Our implementation of transient heat equation converges.