1. Introduction

Thermal transfer is a common and important process in physics and other researchers. Analytical solution can be obtained under certain boundary conditions, but under most realistic conditions, only numerical solutions are accessible. People have introduced several approximate methods of solution to heat conduction problems based on variational principles. With the development in scientific calculation, the finite element method is now commonly applied to field problems.

In the 21st century, the power consumption of building has become a topic in both architecture and physics. It is considered more environmentally friendly and economic to choose more comfortable and power efficient designs. By applying the finite element method to heat transfer calculation of houses, the temperature distribution and power efficiency of different heating systems can be predicted which may help engineers in heat system design.

Precise simulation often involves resolving Navier-Stokes equation, which is skipped in this article by considering a simplified model with given air flow speed. Such simplification is employed as an approach to solve heat convection problem in general heating system. Explicitness and efficiency in such compute method is showcased, without unnecessary omission of certain physical property.

This paper will be divided into four parts. Firstly, general FEM methodology involved in this paper is introduced. The second part of this paper will be the construction of heating system model and corresponding FEM mesh. Then, the specific simplification on the circumfluence in heat convection calculation is illustrated. Finally, there will be comparison on numerical result and discussion on the limitation of such analysis.

2. Literature Review

In order to solve real problems defined on domains with complex boundary condition, the finite element method (FEM) was introduced in the research of heat conduction problem. The application of FEM on field problem including heat conduction started in 1965 [1]. In consequence, time independent finite elements with Gurtin’s variational principle in the study of transient heat conduction problem is formulated. The analysis considered complex solids of arbitrary shape with temperature and heat flux boundary conditions and developed in detail for two-dimensional bodies which are idealized by systems of triangular elements [2,3]. The validation on one-dimensional example of constant heat flux applied to a semi-infinite solid succeeded [4]. Not too long after that, transient linear and non-linear two-dimensional heat conduction problems were solved by the finite element method. Using rectangular prisms in a space-time domain as the finite elements, the finite element method has been found to be stable, convergent to the exact solution [5]. Then, for 2D and 3D transient heat conduction problem, isoparametric finite elements are formulated and are used in a time-stepping solution using the Galerkin approach [6]. More recently, a more generalized finite-element formulation for the space-time domain for heat conduction in solid media is presented [7].

In engineering field, finite element method on heat transfer problem considering heat convection effect was widely applied. Convective heat transfer from internal room surfaces has major effect on the thermal comfort, air movement and heating and cooling loads for the room. Natural convection heat transfer coefficients of a heated wall, a heated floor and a heated ceiling which can be calculated using computational fluid dynamics (CFD) codes [8]. Convection heat transfer and flow calculations for electronic modelling analysis can be done on commercial software Motor-CAD [9]. The finite element method was also capable of solving solidification problem. Different time-stepping schemes were compared for calculations performance with the consistent heat-capacity matrix and with a lumped heat-capacity matrix [10]. Combining proper orthogonal decomposition (POD) with FEM, POD-FEM is a method of reducing the number of degrees of freedom and the overall computing time. Good accuracy of the POD-FEM solver has been observed on solving linear and nonlinear transient heat conduction problems [11]. An approach using a Cartesian reference system to treat the anisotropic thermal conduction of problems for which the solid medium is characterized by a set of tensors of thermal conductivity of different natures (Cartesian and/or cylindrical and/or spherical), with or without phase change has been made [12]. A more specific study can be carried out in terms of specific boundary condition, such as circular enclosure and a sinusoidal cylinder [13], a heated microcantilever and a surrounding air environment [14], etc.

3. Methodology & Result

In different heating systems, heat pipe temperature, heat transfer rate, position, size etc., can be represented by different parameters in the heat transfer equation and its boundary condition. Numerically solve the heat transfer equation through finite element method then it is possible to analyze and compare the heating systems condition. In this section, only a brief introduction of the applied finite element method is included, for more detail, one may infer to [15].

In the finite element method, a computational domain is given, and corresponding mesh is constructed. By the construction of a mesh of an object, its numerical domain is subdivided into a collection of finite elements with variable shapes and sizes, interconnecting with each other in a discrete number of nodes.

Heat transfer system is one of steady state problems involving partial differential equation. The solution of the partial differential equation is approximated in each element by a low-order polynomial. In the form of matrix, global solution is represented in vector \( u \) made up with the solution value in all elements. In term of the heat transfer system, \( u \) is the temperature in the considered domain. Considering only isotropic condition, the primary heat field equation of a heat transfer system (with invariant character of gas and friction heat ignored),

\( \frac{λ}{ρ{C_{p}}}{∇^{2}}u-v\cdot ∇u=\frac{∂u}{∂t}=f in Ω. \) (1)

Where air viscosity is \( λ \) , heat capacity at constant pressure is \( {C_{p}} \) . In steady state, \( \frac{∂u}{∂t}=f=0 \) .

Using Galerkin’s method, the polynomial equation can be represented by matrix equation as:

\( Au=f. \) (2)

The element conduction matrix \( \) is often referred to as the stiffness matrix. It is straightforward to verify \( A \) is a symmetric and positive definite matrix and thus the solution \( u \) exists and unique.

After iterations, multigrid solver gives out the solution of \( Au=f. \) For more details please refer to [16].

The program is constructed on MATLAB based on iFEM*. According to the geometry of the considered domain, cuboid and polygon prism solid tetrahedron meshes are generated. Partial differential equation is discretized by assembling the stiffness matrix according to corresponding low-order polynomial.

*L. Chen. iFEM: an integrated finite element method package in MATLAB. Technical Report, University of California at Irvine, 2009.

In different heating systems, boundary condition varies in value and form. Dirichlet boundary condition is applied to surface where steady temperature is supplied by thermostatic radiator, while on heat dissipating surface like walls, windows and roofs or heat inflowing surface like heated floors Neumann boundary condition is applied.

By subdividing each tetrahedron in a triangulation into eight sub-tetrahedra of equal volume and choosing the ordering of sub-tetrahedron so that recursive application to any initial tetrahedron yields elements of at most three congruence classes, the mesh is refined once in every iteration. Expectation of error in solution drops as iterations time increases. Gradually increase iteration time to improve precision and fineness of the global solution, till the simulation of system temperature is steady and provides enough detail for further evaluations.

In one tetrahedra element, the number 1 element has four nodes 1,2,3,4. Given the temprature in all nodes \( {u_{i}}(i=1,2,3,4) \) , then the temprature at any point \( X \) in the tetrahedron is determined approximately.

\( u(X)=\sum _{i}^{4}{u_{i}}{ϕ_{i}}(X) \) (3)

Where, \( {ϕ_{i}}(X) \) is what so called shape functioin.

\( X=( \begin{array}{c} x \\ y \\ z \end{array} )=(\begin{matrix}{x_{1}} \\ {y_{1}} \\ {z_{1}} \\ \end{matrix} \begin{matrix}{x_{2}} \\ {y_{2}} \\ {z_{2}} \\ \end{matrix} \begin{matrix}{x_{3}} \\ {y_{3}} \\ {z_{3}} \\ \end{matrix} \begin{matrix}{x_{4}} \\ {y_{4}} \\ {z_{4}} \\ \end{matrix})( \begin{array}{c} {ϕ_{1}}(X) \\ {ϕ_{2}}(X) \\ {ϕ_{3}}(X) \\ {ϕ_{4}}(X) \end{array} ) \) (4)

\( {ϕ_{1}}(X)+{ϕ_{2}}(X)+{ϕ_{3}}(X)+{ϕ_{4}}(X)=1 in this tetrahedra element \) (5)

\( {ϕ_{i}}(X)=0 outside this tetrahedra element \) (6)

In a mesh, there are numerous elements, and corresponding sets of shape function to each element. Given (6), in a mesh with N nodes,

\( u(X)=\sum _{i}^{N}{u_{i}}{ϕ_{i}}(X) \) (7)

which is an unique representation of \( u(X) \) .

Generally the simulation of an heat convection involved system requires to solve Navier-Stokes equation, which is a complex and hardware demanding task.

\( \frac{∂v}{∂t}+∇\cdot (vv)=-∇\frac{p}{ρ}+λ{∇^{2}}v+\frac{f}{ρ} \) (8)

In order to get \( v \) without solving Navier-Stokes equation, a cursory but straight-forward simplification of circumfluence is needed. The rotational symmetry of cylinder case gives an explicite example of circumfluence simplification. In case of a cylinder room wraped with heat emitting walls and roof, but the floor is heated evenly. It would be warmer at the center of the floor, and cooler near the walls, thus a simple rotational symmetric coil shaped circumfluence will form.

To keep the condition numerical simple, a rectangular loop is considered as an replacement for the real complex circumfluence.

\( v=\begin{cases} \begin{array}{c} {v_{z}}={v_{{z_{0}}}}(1-2r) ,(z-2r)(z+2r-2) \lt 0 \\ {v_{r}}={v_{{r_{0}}}}(1-z),(z-2r)(z+2r-2) \gt 0 \end{array} \end{cases} \) (9)

Bring the air velocity given by the rectangular loop to equation (1), and apply the boundary condition.

\( \frac{λ}{ρ{C_{p}}}{∇^{2}}u-v\cdot ∇u=\frac{∂u}{∂t}=f in Ω, \)

\( u={g_{D}} in {Γ_{D}}, \frac{∂u}{∂n}={g_{N}} in {Γ_{N}} \) (10)

Multiply a test function \( w \) , integrate over Ω, and use integration by parts to obtain the

corresponding variational formulation:

\( a(u,w)=(f,w) \) (11)

Similar to (7), the test function

\( w=\sum _{i}^{N}{w_{i}}{ϕ_{i}}(X) \) (12)

We then define as isomorphium function space on nodes \( {W_{h}}≅ {R^{N}} \) :

\( w=\sum _{i}^{N}{w_{i}}{ϕ_{i}}(X)⟷w={({w_{1}},…,{w_{N}})^{T}} \) (13)

and call \( w \) the coordinate vector of v relative to the basis \( \lbrace {ϕ_{i}}\rbrace _{i=1}^{N} \) . Following the terminology in elasticity, we reintroduce the stiffness matrix mentioned in(2).

\( A={\lbrace {a_{ij}}\rbrace _{N×N}}, {where a_{ij}}=a({ϕ_{j}},{ϕ_{i}}) \) (14)

and also the load vector \( f \) :

\( f=\lbrace 〈f,{ϕ_{k}}〉 \rbrace _{k=1}^{N}∈{R^{N}} \) (15)

Thus, the matrix product is tantamount to the inner product \( a( \cdot , \cdot ) \) in the variant equation.\

\( a(u,w)=a(\sum _{i}^{N}{u_{i}}{ϕ_{i}}(X),\sum _{i}^{N}{w_{i}}{ϕ_{i}}(X))=\sum _{i}^{N}a({ϕ_{i}},{ϕ_{j}}){u_{i}}{w_{j}}={w^{T}}Au \) (16)

Therefore for any vector \( u∈{R^{N}} \) , \( {u^{T}}Au=a(u,u)≥0 \) , and equals 0 if and only if \( u \) is zero. The solution of \( Au=f \) exists and unique. \( u={A^{-1}}f \) is the numerical solution of the original equation(10).

The height and diameter is 2, heat influx intensity through the floor is 5, heat emmision intensity through the wall and roof is 1. Rectangular circumfluence velocity \( {v_{{z_{0}}}}={v_{{r_{0}}}} \) =2.

When circumfluence is removed, the equation reduce to poisson form. In such case, the simulated is equivalent solid air.

In the cylinder case, to build the initial mesh, is equally divided into sixteen triangular prism segments, then each is divided into four tentrahydron. Under such condition, the heat is sperated more evenly than the case without convection. Thus temprature fluctuation in the whole system is reduced spatially, after reaching the steady state.

Compare the case with and without the air circumluence, the temperature range droped.

4. Discussion

Convection is usually the dominant form of heat transfer in liquids and gases. In the simulation of heat convection system, heated and less-dense fluid move upwards, cool and dense fluid move downwards. Fluid motion is influenced by gravity, pressure and friction. Heat convection calculation is quite complicated and will be even worse with complex boundary condition. Aiming to simplify the problem, a cylinder domain is chosen. Thanks to the axisymmetry, making reasonable simplified assumptions about the air circumfluence is not with much difficulty.

In order to solve the heat convection problem precisely, asset of equations including continuity equation, gas state equation, heat transfer equation and Navia-Stokes equation need to be solved. It is nearly impossible, to solve such a system described by different kinds of parameters: temprature, density, velocity. But by making bold assumption on circumfluence velocity, the unknown parameter to be solved is reduced to only temprature.

Qualitative result can be raised from this research, but no quantitative result is obtained. It needs precise physical constants and material properties to apply this calculation for engineering use.

5. Conclusion

Using the finite element method to solve heat transfer problem in solid medium is common, but the same thing in fluid medium is really complicated. In order to solve heat convection problem and simulate a specific heat system efficiently, Navier-Stokes equation is skiped and air flow velocity of certain circumfluence corresponding to the boundary condition is crassly assumpted.

Based on the finite element method of field problem, a simplified model of heat convection of a room filled with air is solved. The heat transfer in such model is quite different when convection is considered, but we should note that heat convection plays the major role in air heat transfer. To keep the calculation explicit and simple, cylinder axisymmeric domain is chosen, and the heat source is just heated floor with constant heat transfer rate. In the simplified method, air circumfluence is described as a rectangular loop with constant speed and intuitive distribution. The domain is diveded into tetradydron mesh, on which the finite element calculation is based. The result shows that, the circumfluence significantly mitigates the temprature distribution fluctuation in the room, giving a clue that heat convection is important in fluid heat transfer problems. The assumption on circumfluence reduced calculation complexity greatly without affecting qualitive result.

However, some limitations should be noted. First, the idealized circumfluence may not follow physical principles precisely and makes the calculation result diverge from reality. Second, this specific assumption is only reasonable in the corresponding simplified domain and condition. One may found his calculation go wrong, if applied the assumption to a different condition.

After all, this calculation emphasizes explicitness and simplicity, trying to avoid complex calculation needed for heat convection problem. The method might be useful for verifiction of heat convection effect and crude estimate on engineering design.