FEM for Heat Transfer Problems (Finite Element Method) Part 3

Worked Example

Example 12.1: Heat transfer along ID fin of rectangular cross-section

The temperature distribution in the fin, as shown in Figure 12.9, is to be calculated using the finite element method. The fin is rectangular in shape, 8 cm long, 0.4 cm wide, and 1 cm thick. Assume that convection heat loss occurs from the right end of the fin.

Analysis of the problem. The fin is divided uniformly into four elements with a total of five nodes. Each element is with a length of I = 2 cm. The system equation should be 5 x 5. At node 1 the temperature is specified, therefore there is no need to calculate k^ and f£e) for element 1 as only the temperature is requested. As nodes 2, 3 and 4 are internal, there is also no need to calculate k^ and f(e). Since heat convection is occurring on node 5, we have to calculate k^ and f(e) using Eqs. (12.52) and (12.53).

Figure 12.9. A one-dimensional fin of a rectangular cross-section.

Data preparation

Solution. The general forms of the element matrices for elements 1, 2 and 3 are

Substituting the data into the above two equations, we have

The general forms of the element matrices for element 4 are

Substituting the data into the above equations, we have

The next step is to assemble the element matrices together to form the global system equations. The assembly is carried out using the direct assembly procedure described in the previous section. We obtain

Note that at node 1 the temperature is to be fixed at a particular value. This requires a heat exchanger (or heat source) there for it to happen. The heat source is unknown at this stage, and we do not need to know it actually in computing the temperature distribution. However, it is required in balancing the equation. We therefore simply note it as Q*. Since the temperature at node 1 is given as 80°C, we actually have four unknown temperatures in a 5 x 5 matrix equation. This implies that we can actually eliminate Q* in the first system equation as indicated below, and still solve the four unknowns with the remaining 4 x 4 matrix equation. Note that the term in the second row, first column (circled term) must also be accounted for in the four remaining system equations:

Solving the 4 x 4 equation gives

Remarks

Formulations and examples given in this section have clearly shown that the heat transfer problem is very similar to the mechanics problem in terms of FEM treatment. The displacement and force correspond to the temperature and heat flux, respectively. We also showed the analogies of some of the element matrices and vectors between heat transfer problems and mechanics problems. However, we did not mention the mechanics counterparts of matrix and vector that come from the heat convection boundary. Since the heat flux on the boundary depends upon the unknown field variable (temperature), the resulting term (k^) needs to be combined with the heat convection matrix. In fact, in mechanics problems, we can also have a similar situation when the structure is supported by elastic supports such as springs. The reaction force on the boundary depends upon the unknown field variable (displacement) at the boundary. For such a mechanics system, we will have an additional matrix corresponding to and a vector corresponding to f^. However, in mechanics problems, it is often more convenient to formulate such an elastic support using an additional element. The stiffness of the support is then assembled to the global stiffness in the same systematic direct assembly procedure. The reaction force of the elastic support is found after the global equations system is solved for the displacements.

We conclude that the Galerkin residual formulation gives exactly the same FE equations as those we obtained using the energy principles discussed in previous topics.

Composite Wall

Consider heat transfer through a composite wall, as shown in Figure 12.10. The governing equation is given by Eq. (12.9). If there is no heat source or sink in a layer (q = 0 within the layer), one linear element is enough for modelling the entire layer (why?)1. For the case shown in Figure 12.10, a total of three elements, one for each layer, should be used if there is no heat source within these layers. This argument holds even if there is no heat source or sink in the interface between the layers.

As for the boundary conditions, at any one or both of the outer surfaces, the temperature or heat convection or heat insulation could be specified. The convective boundary conditions are given as

Figure 12.10. Heat transfer through a composite wall of three layers. One element is sufficient to model a layer if there is no heat supply/sink in the layers.

and

Therefore, all the elements developed in the section for the 1D fin are also valid for the case of the composite wall, except that and fQe) do not exist because the g and Q vanish in the case of a composite wall (see Eq. (12.10)).

The general form of the element stiffness matrix is given by

where k^ is obtained either by Eqs. (12.52) or (12.55), if there is heat convection occurring on the surface, or else k^ = 0 if the surface is insulated. As for a surface with a specified temperature, there is no need to calculate k^. For the force vector, f^, it is obtained either by Eqs. (12.53) or (12.56) if there is heat convection occurring on the surface, or zero if the surface is insulated. In the case where the surface has a specific temperature, it is not necessary to have for calculation of the temperature. It can always be calculated after the temperature field has been found.

Worked Example

Example 12.2: Heat transfer through a composite wall of two layers

Figure 12.11 shows a composite wall consisting of two layers of different materials. Various heat transfer parameters are shown in the figure. The temperature distribution through the composite wall is to be calculated by the finite element method.

Figure 12.11. Heat transfer through a composite wall of two layers.

Analysis of problem. The wall is divided into two elements with a total of three nodes. Hence, the system equation should be 3 x 3. At node 3, the temperature is specified, therefore there is no need to calculate and for element 2. Since the heat convection is occurring on node 1, k^ and f(e have to be calculated using Eqs. (12.55) and (12.56), respectively.

Data preparation. For element (1),

For element (2),

Solution. The element matrices for element (1) are

The element matrix for element (2) is

The next step is to assemble the element matrices of the two elements together to form the global system equations, which leads to

Note again that Q* is as yet unknown, but is required to balance the equation, as the temperature at node 3 is fixed. Having only two unknowns of temperature, the first two equations in the above give

The above matrix equation (actually consisting of two simultaneous equations) is solved, and the solution is given as

Example 12.3: Heat transfer through layers of thin films

Figure 12.12 shows the process of producing thin film layers of different materials on a substrate using physical deposition techniques. A heat supply is provided on the upper surface of the glass substrate. At the stage shown in Figure 12.12, a layer of iron and a layer of platinum have been formed. The thickness of these layers and the thermal conductivities for these materials of these layers are also shown in the figure. Convection heat loss occurs on the lower platinum surface, and the ambient temperature is 150°C. A heat supply is provided to maintain the temperature on the upper surface of the glass substrate at 300°C. The temperature distribution through the thickness of the layers of the thin film system is to be calculated by the finite element method.

Analysis of problem. This problem is actually similar to the previous example on the composite wall. Hence, 1D elements can be used for this purpose. The layers are divided into three elements with a total of four nodes, as shown in Figure 12.12. Since the temperature is specified at node 1, there is no need to calculate and f^ for element 1. k^ and f^ for element 2 is zero, since there is no heat convection occurring at either nodes 2 or 3. Since there is heat convection occurring at node 4, k^ and f(e) have to be calculated using Eqs. (12.52) and (12.53), respectively.

Data preparation. For element (1),

For element (2),

Figure 12.12. Heat transfer during a thin film deposition process.

For element (3),

Solution. The element matrix for element (1) is

The element matrix for element (2) is

And finally, the element matrices for element (3) are

The next step is to assemble the element matrices of the two elements together to form the global system equations, which leads to

Since φ1 is given to be 300°C, we can therefore reduce the above equation to a 3 x 3 matrix to solve for the remaining three unknown temperatures:

The above matrix equation is solved, and the solution is given as

To calculate the heat flux, Q* on the top of the glass substrate, we can now substitute the temperatures into the first equation of the matrix equation in Eq. (12.106) to obtain

2D Heat Transfer Problem

Element Equations

This section deals with heat transfer problems in two-dimensions that is governed by Eq. (12.1). The procedure for obtaining the FEM equations for 2D heat transfer problems is the same as that for the 1D problems described in the previous sections, except that the mathematical manipulation is more tedious due to the additional dimension.

Let us assume that the problem domain is divided into elements, as shown in Figure 12.5. For one element in general, the residual can be obtained by the Galerkin method as

Note that the minus sign is added to the residual mainly for convenience. The integration in Eq. (12.110) for the residual must be evaluated so as to obtain the element matrices, but in this case, the integration is much more complex than the 1D case because the integration is performed over the area of the element. Recall that in the 1D case, the integral is evaluated by parts, but in this 2D case we need to use Gauss’s divergence theorem instead.

Using the product rule for differentiation first, the following expression can be obtained:

The first integral in Eq. (12.110) can then be obtained by

where Ae is the area of the element. Gauss’s divergence theorem can be stated mathematically for this case as

where θ is the angle of the outwards normal on the boundary Te of the element with respect to the x-axis. Equation (12.113) is thus substituted into Eq. (12.112) to obtain

In a similar way, the second integral in Eq. (12.110) can be evaluated to obtain

The two integrals in Eqs. (12.114) and (12.115) are substituted back into the residual in Eq. (12.110) to give

The field variable φ is now interpolated from the nodal variables by shape functions as in Eq. (12.34), which is then substituted into Eq. (12.116) to give

or in matrix form,

where

As in the 1D case, the vector b(e) is related to the derivatives of temperature (heat flux) on the boundaries of the element. It will be evaluated in the next section in detail. For now, Eqs. (12.120)-(12.122) will be evaluated and analysed.

The integral in Eq. (12.120) can be rewritten in the matrix form by defining

and the gradient vector as

where B is the strain matrix given by

Note that to obtain the above equation, the usual shape function given by Eq. (12.28) is utilized. Using Eqs. (12.123)-(12.125), it can be easily verified that

Therefore, the general element ‘stiffness’ matrix for 2D elements given by Eq. (12.120) becomes

which is exactly the same as Eq. (3.71) that is obtained using Hamilton’s principle, except that the matrix of material elasticity is replaced by the matrix of heat conductivity. Note also that Eq. (12.121) kge) is the same as the matrix given by Eq. (3.75) for mechanics problems. We observe again that the Galerkin weighted residual formulation produces the same set of FE equations as those produced by the energy principle.

Triangular Elements

Using shape functions of a triangular element, the field function of temperature, φ, can be interpolated as follows:

where Ni (i = 1, 2, 3) are the three shape functions given by Eq. (7.22), and φi (i = 1, 2, 3) are the nodal values of temperature at the three-nodes of the triangular element shown in Figure 12.13.

Figure 12.13. Linear triangular element.

Note that the strain matrix B is constant for triangular elements, and is given by Eq. (7.38). Using Eq. (12.127), can be evaluated easily as the integrand is a constant matrix if the material constants Dx and Dy do not change within the element:

Expanding the matrix product yields

It is noted that the stiffness matrix is symmetrical.

The matrix, kg defined by Eq. (12.121) can be evaluated as

The above integral is carried out using the following factorial formula Eq. (7.43), and the fact that the area coordinates are the same as the shape functions, just as we have done for the mass matrix.For example,

Using the area coordinates and the factorial formula in Eq. (7.43), the matrix kg is found as

The element force vector fQ defined in Eq. (12.122) also involves the integration of shape functions, and can also be obtained using the factorial formula in Eq. (7.43):

It is assumed that Q is a constant within the element. The value of QA is equally shared by the three nodes of the rectangular element.

Rectangular Elements

Consider now a four-node, rectangular element as shown in Figure 7.8. The field function φ is interpolated over the element as follows:

Note that for rectangular elements, the natural coordinate system is again adopted, as in mechanics problems (Figure 7.8). The shape functions are given by Eq. (7.51), and are known as bilinear shape functions. The strain matrix for the rectangular element is given by Eq. (7.55). Note that for bilinear elements, the strain matrix is no longer constant. Using Eqs. (12.127), (7.51) and (7.55), can be evaluated as

The matrix kife defined by Eq. (12.121) can be evaluated as

which results in

The element force vector, f^, defined in Eq. (12.122) also involves the integration of the shape functions, and with substitution of the shape functions, Eq. (7.51), it can be obtained as

Note that in the above, Q is assumed to be constant within the element. The value of QA is equally shared by the four nodes of the rectangular element.

In this subsection, matrices k^, kg, fQ have been evaluated exactly and explicitly for rectangular elements. In engineering practice, however, it is very rare to use rectangular elements unless the geometry of the problem domain is also a rectangular one. Very often, quadrilateral elements with four nodes and four non-parallel sides are used for the complex geometry of the problem domain. Formulating FEM equations for quadrilateral elements has been detailed in Section 7.4. Note that with quadrilateral elements, it is difficult to work out the exact explicit form of the element matrices.