Finite Element Method (FEM)

The implementation follows Lewis et al 2004 and the Galerkin FEM description therein, which is a special case of the weighted residual approach. In the one-dimensional FEM, the temperatures are discretized within elements containing specific node points \(x_i\), at which the temperatures are calculated, resulting in the discretization equation

\[T(x) \approx \overline{T}(x) = \sum_{i=1}^n N_i(x)T_i \quad,\]

where \(n\) is the number of nodes within an element and \(N_i\) are the so-called basis or shape functions, which interpolate the temperature within each element (Lewis et al 2004). MoCSI currently uses linear elements (\(n = 2\)) and thus the shape functions follow a linear interpolation with

\[\begin{split}N_1(x) = \frac{x-x_j}{l_m}\\ N_2(x) = \frac{x_i-x}{l_m} \quad,\end{split}\]

where \(x_i\) is the upper position, \(x_j\) is the lower position and \(l_m = x_j-x_i\) is the length of the \(m\)-th element. Using higher-order polynomials, the discretization allows for higher spatial accuracy, but for the current applications, linear elements have been deemed adequate. While most benefits of the FEM are gained in higher dimensions, the one-dimensional formulation allows for easy incorporation of non-uniform grid spacing and also enables a possible variation of the areas at the nodes within the elements, as shown in the following figure}.

Within the Galerkin method, which belongs to the weighted residual methods and where the weights at the nodes are the corresponding nodal shape functions, the formulation of a solution to a differential equation reads

\[\int_\Omega w_i L\left(\overline{T}\right) \text{d} \Omega = 0 \quad,\]

where \(w_i\) are the weights and \(L\left(\overline{T}\right)\) is the residual when using the approximation to the temperature from Equation 1 (Lewis et al 2004). Thus, for the transient heat conduction equation, we arrive at

\[\int_\Omega N_i(x) \left[ \frac{\partial}{\partial x} \left( k \frac{\partial \overline{T}}{\partial x}\right) + Q - c_p \, \rho \, \frac{\partial \overline{T}}{\partial t} \right] \text{d} \Omega = 0\]

and integrating by parts and rearranging the terms leads to the well known form

\[\boldsymbol{C} \, \frac{\partial \overrightarrow{T}}{\partial t} + \boldsymbol{K} \, \overrightarrow{T} = \overrightarrow{f} \quad,\]

with the capacitance matrix \(\boldsymbol{C}\), the stiffness matrix \(\boldsymbol{K}\), the forcing vector \(\overrightarrow{f}\) and the vector of temperatures at all nodes \(\overrightarrow{T}\) (eq. 6.12-6.16 in Lewis_et_al_2004). For the case of constant physical parameters and linear shape functions, for each element these values equate to

\[\begin{split}\boldsymbol{C_i} = \frac{c_p \, \rho \, A_{m+1/2} \, l_m}{6} \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix} \\ \boldsymbol{K_m} = \frac{k \, A_{m+1/2}}{l_i} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \\ \overrightarrow{f} = \frac{A_{m+1/2} \, Q \, l_m}{2} \begin{bmatrix} 1 \\ 1 \end{bmatrix} \quad, \end{split}\]

with the element length \(l_m\) and the geometric mean area \(A_{m+1/2} = (A_m + A_{m+1}) / 2\). In the case of temperature-dependent parameters, we apply an iterative scheme to solve the equation. For the temporal discretization, MoCSI uses a finite difference fully implicit scheme, so the full discretization reads

\[\boldsymbol{C} \, \overrightarrow{T}^{n+1} + \Delta t \, \boldsymbol{K} \, \overrightarrow{T}^{n+1} = \boldsymbol{C} \, \overrightarrow{T}^{n} + \Delta t \, \overrightarrow{f} \quad .\]