FreeFEM Sample Project (a-2-6-1)

Heat diffusion in an Igloo

by Lukas Waldner

Introduction

Aim of this project is to obtain a solution for a transient heat transport problem using the open source software FreeFEM++, which uses a finite element method for spatial discretisation. The discretisation in time should be realised using a Crank-Nicolson scheme. Both, van-Neumann and Dirichlet boundary conditions have been used for the last case. The FreeFEM++ code should first be tested on a simple geometry which provides the ability to compare the numerical to an analytic solution.

In the next step the tested code was used to calculate a solution for a more complicated case. The requested boundary conditions were applied and computed on this geometry.

Method

The heat diffusion problem in an “Igloo” geometry is solved using the Crank-Nicolson-discretization method in time linked to the spatial discretization with finite elements. For the boundary conditions a combination of Neumann- and Dirichlet-boundary conditions is implemented.

For each case the FreeFEM++ code is structured in four parts:

      - Definition of the mesh containing the boundaries of the problem
      - Definition of the problem parameters
      - Definition of the governing PDE
      - The numeric solver (Crank-Nicolson in time)

To compare the numerical and the analytical solution, the $ \text{RMS} $/$ L2 $-Error is calculated for each case which is defined as:
\begin{equation} \text{RMS}(t) = \frac{1}{A} \sqrt{\iint_{\Omega}\left[T(x,t) - T_{\text{exact}}(x,t)\right]^2 dA} \end{equation} $ T_{\text{exact}} $ is only suitable for simple test-cases. There has been no appropriate analytical solution for the “Igloo”-geometry linked to the chosen boundary conditions. The analytical solution for the test cases is correct for Dirichlet-boundary conditions:
\begin{equation} T_{\text{exact}}(x,t) = \frac{1}{4 \pi \sigma t} \exp\left(\frac{x^2 + y^2}{-4\sigma t}\right) \end{equation} with diffusivity $ \sigma $ and vanishing heat source ($ F(x,t)=0 $). The governing PDE for the heat transfer problems is:
\begin{equation} \frac{\partial T}{\partial t} = \sigma \frac{\partial^2 T}{\partial x^2} + F(x,t). \end{equation} The discretisation of the problem considering FEM and a Crank-Nicolson method in time is defined as:
\begin{equation} \!\begin{aligned} & \iint_{\Omega}\left(T_{n+1} - T_n\right)dA = \\ & \iint_{\Omega} \frac{\sigma \Delta t}{2}\left[\left(\frac{\partial^2 T}{\partial x^2}\right)_n + \left(\frac{\partial^2 T}{\partial y^2}\right)_n + \left(\frac{\partial^2 T}{\partial x^2}\right)_{n+1} + \left(\frac{\partial^2 T}{\partial y^2}\right)_{n+1}\right]dA \\ & + \iint_{\Omega} \frac{\Delta t}{2}\left(F_n + F_{n+1}\right)dA \end{aligned} \end{equation} To analyse the convergence behaviour of the FreeFEM++ code, three test cases are considered:

      (a) Test case Square
      (b) Test case Circle
      (c) Test case Ellipse
      (d) Project case "Igloo"

They are compared to the mentioned analytical solution using the $ \text{RMS} $/$ L2 $-Error.

Results and Discussion

(a) Test case Square

Fig. 1. Square geometry

Fig. 2. Square Mesh

Fig. 3. Solution at $ t = 2 $

Fig. 4. Solution at $ t = 3 $

Convergence study for $ t = 3 $

To study the behaviour of the $ \text{RMS} $/$ L2 $-Error a refining process was implemented in the Free-FEM++ code. The following diagrams show the decreasing of the error due to the finer grid and the decreasing time step ($ dt $). One can also observe a saturation behaviour linked to higher time steps. The first diagram contains a refining parameter (grid parameter $ n $) which was decreased from $ 10 $ to $ 50 $. The amount of triangle elements was used to analyse the convergence behaviour is shown in the second diagram.

Fig. 5. Convergence study

(b) Test case Circle

Fig. 6. Circle geometry

Fig. 7. Circle mesh

Fig. 8. Solution at $ t = 2 $

Fig. 9. Solution at $ t = 3 $

(c) Test case Ellipse

Fig. 10. Ellipse geometry

Fig. 11. Ellipse mesh

Fig. 12. Solution at $ t = 2 $

Fig. 13. Solution at $ t = 3 $

(d) Project case "Igloo"

The geometry is shown in Fig. 14. The following boundary conditions are applied:

      - Red boundary: fixed temperature
      - Yellow boundary: insulation
      - Green boundary: constant heat flux

Initial state is defined by constant temperature in the whole space. In the next step, the temperature field is calculated.

Fig. 14. Igloo geometry

Fig. 15. Igloo mesh

Fig. 16. Solution at $ t = 2 $

Fig. 17. Solution at $ t = 3 $

Remark

An analytical solution consistent with the problem has to be chosen. In the present work, due to the pure diffusive behaviour of the heat equation, the diffusion of a point-wise spot was used as benchmark.

Note: The relevant source codes will be added here later.


Back to Main Menu: Numerical Methods in Fluid Dynamics