Tutorial 7 – Steady Heat Transfer by Conduction

Posted on July 17, 2017 - Tutorials

Problem Definition

In this tutorial we are going to show how to solve a usual heat conduction problem, which can normally be solved with $\tt{genericHeatSolver2D}$ except for the fact that today we choose our conductive medium to be composed of two distinct materials, such as concrete and styrofoam.

Geometry of the domain is depicted in figure below and has been created in millimeters. Thus, we will need to recalculate it to meters. You will see an additional scaling step in the script.

Both side walls (denoted with Physical Line id = 88) are insulated, bottom wall (id 87) is exposed to ambient air temperature of 253 K and heat transfer coefficient on this boundary equals 24 W/(m^2K) and ambient air temperature around the top wall (id 86) is 293 K with a heat transfer coefficient of just 8 W(m^2K). Thermal conductivity of the styrofoam equals 0.03 W/(mK) and thermal conductivity of the concrete equals 0.35 W/(mK).

Mesh Preparation

Mesh for our problem has been generated in Gmsh. Tutorial 1 introduces into geometry and mesh preparation in Gmsh. The significant difference, when preparing mesh for a simulation of heat transfer in multiple materials is that the solver must be able to determine, which finite elements represent which material. To make this possible, you have to assign in Gmsh two different Physical Surfaces to each of the subdomains. Ids of the physical surfaces created in Gmsh will later be passed to QuickerSim CFD Toolbox thanks to importMeshGmsh function.

Resulting mesh is shown in figure below. Subdomain containing concrete has got physical surface id = 89 and subdomain with styrofoam physical surface id = 90. We will need these values, when specifying distinct thermal conductivities of each material.

Mathematical Model

Mathematical model of a heat conduction problem is Laplace’s (diffusion) equation casted below:

    $$\nabla \cdot (\lambda \nabla T) = 0$$

This partial differential equation must be supplemented with appropriate boundary conditions, which for our problem are as follows:

Insulated boundaries:

    $$\vec n \cdot \nabla T = 0$$

At the top and bottom wall we impose a convection boundary condition, known in heat transfer also as third kind boundary condition or Robin boundary condition. The usual engineering formulation of that condition is as follows:

    $$q_{in} = \alpha(T_{out}-T_{wall})$$

where $q_{in}$ stands for heat flux density entering the solid, $\alpha$ stands for heat transfer coefficient, $T_{out}$ for ambient air or fluid temperature and $T_{wall}$ actual wall temperature.

Contrary to this very popular engineering definition of a heat convection boundary condition, it is recognized as Robin boundary condition described in maths with the following formula:

    $$\vec n \cdot \lambda \nabla T + b_2 \cdot T = b_1 $$

First term of the above sum represents exactly the incoming heat flux density $q_{in}$. We only need to redefine two other terms to have consistent definitions. Please note that few trasformations lead to the observation that choosing:

    $$b_1 = \alpha T_{out}$$

and

    $$b_2 = \alpha$$

yield identical definitions. QuickerSim CFD Toolbox works with the proper definition of Robin boundary condition, so once solving a heat transfer problem, you have to remember to recalculate ambient temperature and heat transfer coefficient. You can also find a corresponding example in help to the $\tt{imposeScalarBoundaryConditions2D}$ function.

Code

In this section we present and comment the code which solves the above problem. First we import the mesh and display it to the screen. You got used to the fact that we work on second order finite elements, if we want to solve a fluid flow problem. That’s not needed in case of heat conduction, so in this example we will be working on linear finite elements. In this case, our original mesh from the mesh file is OK and we do not need to call $\tt{convertMeshToSecondOrder}$ function.

clc;
clear;


% Import mesh, scale from milimeters to meters and display
[p,e,t] = importMeshGmsh('brick.msh');
p = p/1000; % coordinates are stored in p array; Divide by 1000 to scale
displayMesh2D(p,t);

Compare this script with $\tt{genericHeatSolver2D}$, since this one is a simplification of a standard heat conduction solver available in QuickerSim CFD Toolbox. We deal with a linear problem (constant thermal conductivities), so we won’t be iterating. We just once solve system of linear equations, so we omit defining convergence criteria and initialization of the temperature field. However, we need to define appropriate thermal conductivities. Let us check number of elements (elements are stored in t array after importing the mesh – consult it with documentation) and declare in MATLAB a vector of the same size which will store thermal conductivity of each finite element.

nelements = size(t,2);
lambda = 0.35*ones(nelements,1);

By writing the above code we assigned to all elements conductivities of the concrete. We now need to fix our intended mistake and assign proper conductivities to styrofoam (denoted with id = 90). Ids of each element are stored in the last row of the t array. So we assign a lower conductivity for all elements of lambda, for which last row of t has value 90.

lambda(t(end,:)==90) = 0.03;

An optional step is to check, whether we have done it correctly. We would like to display the thermal conductivity field with $\tt{displaySolution2D}$ function. Unfortunately, the latter function requires that a solution (field) value is defined for each node in the mesh. We have defined conductivities for elements. If we want to display conductivity field, we have to approximate it by distributing values from elements to the nodes of these elements. It will be an approximation, since we will obtain averaged values at the concrete/styrofoam interface. Two lines of code below peform the above mentioned visual check. You can also omit this step and you will still have a successful simulation.

lambdaInNodes = moveDataFromElementsToNodes(p,t,lambda);
displaySolution2D(p,t,lambdaInNodes,'Thermal conductivity [W/(mK)]');


% We will now assemble the whole problem matrix and impose convection
% boundary conditions that we have discussed earlier.


% Assemble problem matrix
[D,F] = assembleDiffusionMatrix2D(p,t,lambda);


% Apply boundary conditions
[D,F] = imposeScalarBoundaryCondition2D(p,e,D,F,86,'robin',8*293,8);
[D,F] = imposeScalarBoundaryCondition2D(p,e,D,F,87,'robin',24*253,24);


% Solve equations
T = D\F;


% Display solution
figure(4);
displaySolution2D(p,t,T,'Temperature [deg C]');

Postprocessing

We will only do basic postprocessing of results here by checking the balance of heat flux coming into the domain and coming out. To do that, we first compute the heat flux $\vec{q} = - \lambda \nabla T$. We compute the gradient (the function returns a matrix with two columns – the x-component of heat flux vector is stored in the first column and the y-component in second column). That’s why we have to replicate the values of lambda in each node into two columns (MATLAB will be then able to perform element-by-element multiplication). To copy a 1-column vector $\tt{lambdaInNodes}$ into matrix of two identical vectors, we use the $\tt{repmat}$ function. Eventually, we extract the x-component of $\vec{q}$ and its y-component and pass it to $\tt{boundaryFlux2D}$ function to compute total heat flux through a certain edge. Please note that total heat flux on insulated bounaries (id 88), which should be zero is orders of magnitude lower than that on upper and lower boundary. The latter two are equal to each other within discretization and numerical accuracy and show that the total energy balance is conserved.

q = -repmat(lambdaInNodes,1,2).*solutionGradient2D(p,t,T);
qx = q(:,1);
qy = q(:,2);
Qin = boundaryFlux2D(p,e,[qx;qy],86)
Qout = boundaryFlux2D(p,e,[qx;qy],87)
Qinsulated = boundaryFlux2D(p,e,[qx;qy],88)

Tags