Test case of Stationary Heat Equation in Three Dimensions

1. Introduction

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

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.

Static Heat equation

\[\text{(Heat)} \left\{ \begin{matrix} - 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

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"
                }
            }
        }
    }

7. Weak Formulation

We obtain :

Weak Formulation of Static Heat Equation

\[ \text{(Weak Heat Static)} \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\)

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

\(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

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

10. Numeric Parameters

This section show the parameters used to compute the simulation.

Size of mesh : \(1 \,mm\)
Element type : \(P1\)
Solver : automatic
Number of CPU core : \(8\)

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

11. Result

We obtain :

Temperature \(T\) \((K)\) on geometry

Temperature \(T\) \((K)\) on radial cut

\(T\) \((K)\) on \(Or\) axis

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 and exact solution (\(||T_{num}-T_{exact}||_{L^2(\Omega)}\)) for different mesh sizes.

Table 1. Table of data
hsize (m)	\(\Vert T_{num} - T_{exact} \Vert_{L^2} (K)\)
1e-2	1.658663e-01
3.16e-3	1.9202810377953374e-02
1e-3	1.9148176684869753e-03

Our implementation of static heat equation converges.