Weighted Residual Approach For FEM
In that case, one needs to know the functional for the physical problem. For many engineering problems, one does not know the functional or it is not known as intuitively as in mechanics problems. Instead, the governing equation for the problem would be known. What we want to do now is establish FEM equations based on the governing equation, but without knowing the functional. In this case, it is convenient to use the weighted residual approach to formulate the FEM system equations.
The general form of Eq. (12.1) can be rewritten in the form
where the function f is given as
In general, it is difficult to obtain the exact solution of φ(χ, y) which satisfies Eq. (12.24). Therefore, an approximated solution of φ(χ, y) is sought for, which satisfies Eq. (12.24) in a weighted integral sense, i.e.
where W is the weight function. We hope that the solution of φ(χ, y) that satisfies Eq. (12.26) can be a good approximation of the exact results. This is the basic idea behind the weighted residual approach. This approach is very simple, and can be used in the finite element method to establish the discretized system equations as described below.
To ensure a good approximation, the problem domain is divided into smaller subdomains (elements), as we have done in previous topics, as shown in Figure 12.5. In each element, it is assumed that
where the superscript h indicates that the field variable is approximated, and
in which Ni is the shape function of x and y, and nd is the number of the nodes of the element.
Figure 12.5. Division of problem domain Ω bounded by Γ into elements.
For the triangular element shown in Figure 12.5, nd = 3. In Eq. (12.27), Φ is the field variable at the nodes of the element. There are a number of ways in which the weight function, W, can be chosen when developing the element equations. When the shape functions are used as the weight function, the method is called the Galerkin method, which is one of the most popular methods for developing the FE equation.
Using the Galerkin method, the residual calculated at all the nodes for an element is then evaluated by the equation
Finally, the total residual at each of all the nodes in the problem domain is then assembled and enforced to zero to establish the system equation for the whole system. FEM equations will be developed in the following sections for one- and two-dimensional field problems, using a heat transfer problem as an example.
Id Heat Transfer Problem
One-Dimensional Fin
Consider the one-dimensional fin shown in Figure 12.6. The governing differential equation for a steady-state heat transfer problem for the fin is given by Eq. (12.7). The boundary conditions associated with Eq. (12.7) usually consist of a specified temperature at x = 0
and convection heat loss at the free end, where x = H
where φb is the temperature at the end of the fin and is not known prior to the solution of the problem. Note that the convective heat transfer coefficient in Eq. (12.31) may or may not be the same as that in Eq. (12.7).
Figure 12.6. One-dimensional problem: heat transfer in a thin fin that is divided into n elements.
Using the same finite element technique, the fin is divided into elements as shown in Figure 12.6. In one element, the residual equation, can be obtained using the Galerkin approach as in Eq. (12.29):
where
for heat transfer in thin fins. Note that the minus sign is added to the residual mainly for convenience. Integration by parts is then performed on the first term of the right-hand side of Eq. (12.32), leading to
Using the usual interpolation of the field variable, φΗ, by the shape functions in the 1D case,
and substituting Eq. (12.34) into Eq. (12.33), gives
or
where
is the element matrix of thermal conduction. Matrix
in Eq. (12.36) is defined by
which is the matrix of thermal convection on the circumference of the element. Vector fg is associated with the external heat applied on the element, defined as
Finally, b(e) is defined by
which associates with the gradient of the temperature (or heat flux) at the two ends of the element.
InEq. (12.37),
is the same as the strain matrix in the case of mechanics problems. For linear elements, the shape functions are as follows:
and
Substituting the above equation into Eq. (12.37), the heat conduction matrix is obtained as
In deriving the above equation, D = kA has been used. Compared with Eq. (4.15), it can be seen that kD^ is in fact the analogy of the stiffness matrix of a truss structure. The tensile stiffness coefficient AE/l has been replaced by the heat conductivity coefficient of kA/l.
Similarly, to obtain the convection matrix
corresponding to the heat convection, substitute Eq. (12.42) into Eq. (12.38),
In deriving the above equation, g = hP has been used. Compared with Eq. (4.16), it can be seen that kge) is in fact the analogy of the mass matrix of a truss structure. The total mass of the truss element Apl corresponds to the total heat convection rate of hPl.
The nodal heat vector fQ is obtained when Eq. (12.42) is substituted into Eq. (12.39), giving
In deriving the above equation, Q = q + hPφf (see Eq. (12.8)) has been used. The nodal heat vector consists of the heat supply and the convectional heat input to the fin. The nodal heat vector is the analogy of the nodal force vector for truss elements.
Finally, let us analyse the vector b(e) defined by Eq. (12.40), which is associated with the thermal conditions on the boundaries of the element:
where the subscripts ‘L and ‘R’ stand for the left and right ends of the element. It can be easily proven that at the internal nodes of the fin, b^ and be vanish when the elements are assembled, as illustrated in Figure 12.7. As a result, b needs to be determined only for the nodes on the boundary by using the conditions prescribed. At boundaries where the temperature is prescribed, as shown in Eq. (12.30), b^ or bR are yet to know. In fact, there is no need to know b^ or bR^ in the stage of solving the system equation, as the temperature at the node is already known. The situation is very much the same at the prescribed displacement boundary, where the reaction force is usually unknown at the stage of solving system equations.
However, when there is heat convection at the ends of the fin, as prescribed in Eq. (12.31), bL or bR are obtained using the boundary conditions, since the heat flux there can be calculated.
Figure 12.7. Vector b(e) vanishes at internal points after assembly.
For example, for boundary conditions defined by Eq. (12.31), we have
Since φb is the temperature of the fin at the boundary point, we have φb = φ], which is an unknown. Equation (12.49) can be rewritten as
or
where
and
Note that Eqs. (12.51)-(12.53) are derived assuming that the convective boundary is on node ], which is the right side of the element. If the convective boundary is on the left side of the element, we then have
where
and
Substituting the expressions for b(e) back into Eq. (12.36), we obtain
or in a simplified form of
where
and
Note that and f£e) exist only if the node is on the convective boundary, and they are given by Eqs. (12.52) and (12.53) or Eqs. (12.55) and (12.56), depending on the location of the node. If the boundary is insulated, meaning there is no heat exchange occurring there, both k^ and f(e) vanish, because d^/dx = 0 at such a boundary. After obtaining the element matrices, the residual defined by Eq. (12.58) is assembled, and enforced to equate to zero, which will lead to the following global system equation:
The above equation has the same form as that for a static mechanics problem. The detailed assembly process is described in the examples that follow.
Direct Assembly Procedure
Consider an element equation of residual in the form
or in the expanded form
Consider now two linear one-dimensional elements, as shown in Figure 12.8. We have for element 1,
Figure 12.8. Assembling of two elements demonstrating the direct assembly principle.
and for element 2,
The total residual at a node should be obtained by summarizing all the residuals contributed from all the elements that are connected to the node, and the total residual should vanish to ensure the best satisfaction of the system equations. We then have at node 1,
at node 2,
and at node 3,
Writing Eqs. (12.68), (12.69) and (12.70) in matrix form gives
which is the same as the assembly procedure introduced in Section 3.4.7 and Example 4.2.
