Multiphysics Learning & Networking

Introduction

The Finite Element Method (FEM) provides a much different way from the FDM and FVM for the discretization of the continuous mechanics description of a physical process. The discretization is made over a domain consisting of subdivisions called elements. These elements are essentially the same as those control volumes in the FVM from a geometric perspective. However, the FEM more emphasizes the grid points as the apexes of the element rather than the centroids, volume, and boundary faces that are emphasized in the FVM. The governing equation is discretized at all the grid points of the element to form a local algebraic equation. Since different elements share grid points, the local algebraic equation will be ensembled into the global governing equation to add up the coefficients from different elements for the same point. Another major difference is that the FEM adopts the weak form of the governing equation instead of the differential form in the FDM or the integral form in the FVM. Usually, either the Galerkin method or variational methods from the calculus of variations are used to derive the weak form. The idea underlying the Galerkin and the variational methods is analogous to the principle of minimum energy in physics. In the FEM, we use polynomials of different orders, which also determine the order of the elements, to approach the true solution. The Galerkin method or variational methods will minimize the error of the approximation by fitting trial functions into the governing equation.

Galerkin Method and Weak Form

Suppose we have a linear differential operator $D$ acting on a function $u$ to produce a function, $p\left(x\right)$, \[D\left(u\left(x\right)\right)=p\left(x\right).\] In the minimal weight residual methods, which the FEM belongs to, we approximate $u$ by a function, $\tilde{u}$, which is a linear combination of basis functions chosen from a linearly independent set: \[u\cong \tilde{u}=\sum _{i=1}^{n}u_{i} \varphi _{i} .\] Substituting the above approximation into the differential operator, $D$, in general, the result of the operation is not $p\left(x\right)$ anymore. Accordingly, an error or residual exists, \[E\left(x\right)=R\left(x\right)=D\left(\tilde{u}\left(x\right)\right)-p\left(x\right)\ne 0.\] The notion in the minimal weight methods is to force the residual to approximate zero in some average sense over the domain, that is, \[\int _{X}R\left(x\right)W_{i} dx =0, i=1,2,...,n\] where the number of weight functions $W_{i} $ is exactly equal to the number of unknown constants $u_{i} $ in $\tilde{u}$. The result is a set of algebraic equations ($n$) for the unknown constants, $u_{i} $. There are five common minimal weight residual methods, which are differentiated by the choices of the $W_{i} $'s. These methods include the collocation method, sub-domain method, least squares method, Galerkin method, and method of moments. In the Galerkin Method, the derivative of the approximating function is used as the weight function: \[W_{i} =\frac{\partial \tilde{u}}{\partial u_{i} } =\varphi _{i} \left(x\right).\] Then we have \[\int _{X}\left[D\left(\tilde{u}\left(x\right)\right)-p\left(x\right)\right]\varphi _{i} \left(x\right)dx =0. \] This is the weak formulation, or weak form, of the original differential equation. $\varphi _{i} \left(x\right)$ is the shape function in the FEM. If the element type is fixed, then all the shape functions for each point will be known prior to the solution process. The outcome of using the weak-form is to convert the governing PDE, which requires the equality of the equation at every infinitesimal point within the domain and on the boundary, to the weak form of the equation, which requires that the above integral-form of a function with the residual of the numerical approximation and the shape functions equals zero. The integral form will bring additional benefits: a relatively lower requirement on the continuity and differentiability of the original differential equation. This weak-form is necessary in the FEM for obtaining the algebraic equation system. However, it is worthwhile to mention that the above way to obtain the weak form can be reinterpreted in a more straightforward way using test functions, which are random functions as follows, \[v\cong \tilde{v}=\sum _{i=1}^{n}v_{i} \varphi _{i} .\] The procedure is to multiply the original PDE with the above formation of the test function and then to integrate their product over the domain. As the test function is a random function and exists in every term of the integral form of the equation, $v_{i} $ can be eliminated and we will obtain the same weak formulation.

Example

In the FEM, the above equation in continuum mechanics will be discretized into an algebraic equation system. Taking the bilinear quadrilateral element (Quadrilateral 4 (Q4)) element in 2D for example, the values of the dependent variable $T$ and test function $v$ can be related to the shape functions $N_{Ti} $ ($N_{T} $) and test functions $v_{i} $ ($v$) of the element in the following way: \[T=\sum _{i=1}^{4}N_{Ti} T_{i} =N_{T} \cdot T,\] \[v=\sum _{i=1}^{4}N_{Ti} v_{i} =N_{T} \cdot v,\] where, \[N_{{\rm T}} =\left[\begin{array}{cccc} {N_{1} } & {N_{2} } & {N_{3} } & {N_{4} } \end{array}\right],\] \[T=\left[\begin{array}{cccc} {T_{1} } & {T_{2} } & {T_{3} } & {T_{4} } \end{array}\right]^{T} .\] The discretized governing equation is formulated as \[\rho C\int _{\Omega }\left(N_{{\rm T}}v\right)^{T}\frac{\partial T}{\partial t} d\Omega=\] \[-\lambda\int _{\Omega }\left(A_{{\rm T}} N_{{\rm T}} v\right)^{T}\cdot\left(A_{{\rm T}} N_{{\rm T}} T\right)d\Omega +\int_{\partial \Omega}\left(N_{{\rm T}}v\right)^{T}\left(-J\cdot n\right)dS+\int_{\Omega }\left(N_{{\rm T}} v\right)^{T}g d\Omega\] where $N$ is the shape function matrix and $A$ is the matrix of the differential operator. The integration will be conducted using the Gauss-Legendre numerical integration. The time differential on the left side of the equation will be handled with an iterative process: \[\left(\frac{\rho C}{\Delta t} \int _{\Omega }N^{T} N d\Omega \right)\left(T^{n+1} -T^{n} \right)=\] \[\left(-\lambda \int _{\Omega }B^{T} B d\Omega \right)\left[\theta T^{n+1} +\left(1-\theta \right)T^{n}\right]+\int_{\partial \Omega }N^{T}\left(-J^{n}\cdot n\right)dS+\int_{\Omega }N^{T} g^{n} d\Omega\] where $T^{n+1} $ and $T^{n} $ are the vectors of unknowns at the next and the current time steps, respectively. The test function is eliminated due to its presence in all terms. Detailed expressions of the matrices in the above equation are as follows, $A_{{\rm T}} =\left[\begin{array}{c} {\frac{\partial }{\partial x} } \\ {\frac{\partial }{\partial y} } \end{array}\right]$(for equation with scalar dependent variables in 2D, independent of element type), \[B_{{\rm T}} =\left[\begin{array}{cccc} {\frac{\partial N_{1} }{\partial x} } & {\frac{\partial N_{2} }{\partial x} } & {\frac{\partial N_{3} }{\partial x} } & {\frac{\partial N_{4} }{\partial x} } \\ {\frac{\partial N_{1} }{\partial y} } & {\frac{\partial N_{2} }{\partial y} } & {\frac{\partial N_{3} }{\partial y} } & {\frac{\partial N_{4} }{\partial y} } \end{array}\right],\] where $N_{i} $' s are the shape functions corresponding to the type of the adopted element. For example, the Q4 element has four functions $N_{1} $, $N_{2} $, $N_{3} $, and $N_{4}$.

There are two different ways to deal with the boundary source term specified by the Neumann boundary condition in the FEM. The first way is to calculate the term with shape functions associated with the nodes on the boundaries when discretizing the weak-form of the governing equation; while the other one is to establish the force (or source) vector matrix separately by multiplying the boundary source with the boundary length (area). For both of them, the boundary force needs to be formulated based on physical mechanisms such as convective heat transfer, which can be predicted with Newton's law of cooling, \[-J\cdot n=h\left(T_{{\rm env}} -T\right), \] where $T_{{\rm env}} $ is the environmental temperature. The second method is employed in this code for the practice problem. When no heat source exists, the discretized governing equation can be expressed as \[K_{{\rm Ta}} \left(T^{n+1} -T^{n} \right)=K_{{\rm Tc}} \left[\theta T^{n+1} +\left(1-\theta \right)T^{n} \right]+F_{{\rm Ts}} ^{n} +F_{{\rm Tb}} ^{n} ,\] where $K_{{\rm Ta}} $ and $ K_{{\rm Tc}} $ are the stiffness matrices relevant to the accumulation (a) and conduction (c), respectively; and $F_{{\rm Ts}} $, $F_{{\rm Tb}} $are the force vectors corresponds to source on boundary surface (s) and body source (b, corresponding to $g$). The term "force" is adopted here to represent the generalized force: source. In the current study, no source (body force) is present, so $F_{{\rm Tb}} $ is excluded. The stiffness matrices are calculated using the following equations: \[K_{{\rm Ta}} =\left(\frac{\rho C}{\Delta t} \int _{\Omega }N_{{\rm T}} ^{T} N_{{\rm T}} d\Omega \right),\] \[K_{{\rm Tc}} =\left(-\lambda \int _{\Omega }B_{{\rm T}} ^{T} B_{{\rm T}} dS\right).\] If the explicit method is adopted, the solution at the current step can be obtained as \[T^{n+1} =C_{{\rm T}} T^{n} +M_{{\rm T}} =\frac{K_{{\rm Ta}} +\left(1-\theta \right)K_{{\rm Tc}} }{K_{{\rm Ta}}-\theta K_{{\rm Tc}}} T^{n} +\frac{F_{{\rm Ts}} }{K_{{\rm Ta}}-\theta K_{{\rm Tc}} }\]