B Appendix B - The Finite Element Method

The goal of the finite element method [51] is to determine an approximate solution T ^e that satisfies the differential equation and its boundary conditions through the means of a variational or weighted-integral method. The method is illustrated here for a linear set of basis functions.

To begin, the radial energy equation (7) is written in weak Galerkin form:

(99)

The approximate solution T ^e that is sought is one for which a weight function h(r) can be found that satisfies the equation over the entire domain in a weighted-integral sense. An integration by parts on the term that contains the highest order derivative results in the equation

(100)

Note that this integration by parts has reduced (weakened) the polynomial order required of the solution T ^e by shifting some of the differentiation to the weight function. The temperature solution is now required to be only linearly continuous throughout the domain; that is, the first derivative of temperature may now be discontinuous.

(The integration described above is actually performed over the entire spherical domain. Since neither the weight function h(r) nor the quantity in curly brackets depends on angle θ or φ, integration has been performed over both angles, and the resulting factor of 4π has been divided out.)

The domain of the problem is now subdivided into subdomains called finite elements. If a solution can be found that solves the above weighted-integral equation in each individual element without sacrificing the necessary continuity between elements, then the desired approximate solution has been found.

From equation (100), the equation that must hold in each element i bounded by nodes at radii r_i and r_i+1 is

(101)

Note that the boundary condition is applied only for the outermost element, for which i = N − 1, where N is the number of nodes.

A local coordinate system is now defined,

(102)

with L(t) = (r_i+1(t)− r_i(t)) the radial thickness of the element. The average radius of element i is defined as r_av = (r_i + r_i+1)/2. Since within the element neither the weighting function h nor the temperature solution T requires more than a first order polynomial, these functions can be represented as a sum of nodal values multiplied by a set of linear basis functions,

(103)

and

(104)

where h₁ = (1 − z)/2 and h₂ = (1 + z)/2. The heat sink term can be described similarly as

(105)

where _j = ρ_pφ _pB exp(−E/RT_j) is the value of the Arrhenius expression at the local node j. This change of coordinates relates the global numbering of nodes i and i + 1 bounding element i to local node numbers j = 1 and 2 respectively.

Substituting the local coordinate system into equation (101) and applying the outer boundary condition from equation (65) results in the following matrix equation for element i:

(106)

where

(107)

(108)

(109)

(110)

(111)

(112)

After integration, the matrices become

(113)

(114)

(115)

(116)

The next step in setting up the finite element problem is assembling the matrices from individual elements to form a single linear algebra equation to be solved. The two requirements for connecting the entire domain are continuity of temperature and balance of heat flux from one element to the next. The first requirement is satisfied by using the global numbering of nodes to sum matrices of individual elements in the correct positions. The final matrices [B], [C], [A], and [G] are N x N tridiagonal matrices. The second requirement is satisfied automatically by the finite element method as a "natural" boundary condition. To see this, consider again the radial energy equation in weak Galerkin form (99), but this time consider only the portion of the integration carried out over two adjacent elements a and b, with thermal conductivities k_a and k_b. These elements are defined by nodes a1 and a2 and nodes b1 and b2 respectively, where nodes a2 and b1 are identical and shared. Performing an integration by parts on the highest order derivative term in the weak Galerkin equation integrated from node a1 to node b2 results in an integrated part plus boundary terms:

(117)

The integrated parts have been set to 0 through the formulation of the linear algebra equation as already described. Since the boundary terms must also sum to zero, it must be true that at shared nodes the terms must be equal. At node a2/b1 therefore,

(118)

This is precisely the required condition that heat flux is matched from one element to the next.

(119)

where the matrices have been assembled from those for individual elements. This equation is treated numerically by using the standard tools for partial differential equations to set up the matrix equation

(120)

to determine the nodal temperatures at the next timestep t_n+1 = t_n + Δt from the temperatures at the current timestep t_n. For the Crank-Nicolson scheme used in the code for this model, the matrices in this equation are

(121)

(122)

(123)

The superscripts n and n + 1 for vectors {T}, {}, and {K} refer to the application of values of temperature at the previous timestep and those to be determined at this timestep, respectively. Note that {ⁿ⁺¹} and {K ⁿ+1} appear in the RHS of this equation The values in these vectors are determined iteratively at each timestep until the solution converges. The superscripts n and n + 1 for the matrices [A], [B], [C], and [G] indicate that element sizes and material properties from previous and current timesteps respectively are applied.

Since all matrices in this equation are tridiagonal, the solution for each timestep is obtained quickly and easily. The finite element code written for this model was verified by comparison with an analytical solution for heat transport in a sphere [54].

For quadratic basis functions, h₁ and h₂ are replaced by the set

(124)

and the solution is obtained using a sparse matrix solver.