Stabilization for the Finite Element Method

Both the neutron transport equation and the mass transport equations are hyperbolic and are numerically unstable when solved using the finite element method. There are several approaches to mitigate the oscillations that occur in the neutron transport equation and mass transport problems dominated by streaming and advection (respectively). The following sections discuss the approaches taken by Gnat for both sets of governing equations, starting with the neutron transport equation and ending with the advection-diffusion equations.

Neutron Transport

Gnat provides two stabilization schemes for the neutron transport equation, both being described below.

1. The Self-Adjoint Angular Flux Method

The first approached implemented is the self-adjoint angular flux (SAAF) method derived by (TODO: REFERENCE HERE). We can write the transport equation where the collision term is on the left and all remaining terms are on the right:

var element = document.getElementById("moose-equation-d9cb2fcd-f896-4238-afce-591dfcdaf874");katex.render("\\Sigma_{r,\\,g}\\Psi_{g} = \\sum_{g' = 1}^{G}\\int_{4\\pi}\\Sigma_{s,\\, g'\\rightarrow g}f_{g'\\rightarrow g}\\Psi_{g'}\\, d\\hat{\\Omega}' + S_{g} - \\frac{\\partial}{\\partial t}\\frac{\\Psi_{g}}{v_{g}} - \\hat{\\Omega}\\cdot\\vec{\\nabla}\\Psi_{g}", element, {displayMode:true,throwOnError:false});

Inverting the collision operator yields the following equation, which is known as the angular flux equation:

(1)var element = document.getElementById("moose-equation-4544425d-02d3-406b-8bdc-4ec9ab4215fb");katex.render("\\Psi_{g} = \\frac{1}{\\Sigma_{r,\\,g}} \\Big(\\sum_{g' = 1}^{G}\\int_{4\\pi}\\Sigma_{s,\\, g'\\rightarrow g}f_{g'\\rightarrow g}\\Psi_{g'}\\, d\\hat{\\Omega}' + S_{g} - \\frac{\\partial}{\\partial t}\\frac{\\Psi_{g}}{v_{g}} - \\hat{\\Omega}\\cdot\\vec{\\nabla}\\Psi_{g}\\Big)", element, {displayMode:true,throwOnError:false});

Eq. (1) can be substituted into Governing Equations, Eq. (3) and algebraically simplified to yield the SAAF weak form of the neutron transport equation:

(2)var element = document.getElementById("moose-equation-dab1de19-d04e-4e82-b6a3-0db749f6d69d");katex.render("\\Big( \\psi_{j} + \\frac{1}{\\Sigma_{r,\\,g}}\\vec{\\nabla}\\psi_{j}\\cdot\\hat{\\Omega},\\, \\frac{\\partial}{\\partial t}\\frac{\\Psi^{h}_{g}}{v_{g}} \\Big)_{V} + \\Big\\langle \\psi_{j},\\, \\hat{n}\\cdot\\hat{\\Omega}\\Psi^{h}_{g} \\Big\\rangle_{\\Gamma} + \\Big( \\vec{\\nabla}\\psi_{j}\\cdot\\hat{\\Omega},\\, \\frac{1}{\\Sigma_{r,\\,g}}\\hat{\\Omega}\\cdot\\vec{\\nabla}\\Psi_{g}^{h}\\Big)_{V} + \\Big( \\psi_{j},\\, \\Sigma_{r,\\,g}\\Psi^{h}_{g} \\Big)_{V}\\\\ - \\Big( \\psi_{j} + \\frac{1}{\\Sigma_{r,\\,g}}\\vec{\\nabla}\\psi_{j}\\cdot\\hat{\\Omega},\\, \\sum_{g' 1}^{G}\\int_{4\\pi}\\Sigma_{s,\\, g'\\rightarrow g}f_{g'\\rightarrow g}\\Psi^{h}_{g'}\\, d\\Omega' \\Big)_{V} - \\Big( \\psi_{j} + \\frac{1}{\\Sigma_{r,\\,g}}\\vec{\\nabla}\\psi_{j}\\cdot\\hat{\\Omega},\\, S_{g} \\Big)_{V} = 0", element, {displayMode:true,throwOnError:false});

Eq. (2) is second order and therefore stable for continuous finite elements. However, the division by the removal cross-section makes Eq. (2) ill-suited for regions where the removal cross-section is void or close to being a void. To remedy this an alternative form of the angular flux equation is proposed by (TODO: REFERENCE HERE):

(3)var element = document.getElementById("moose-equation-740f8be6-cb00-4e15-9fcc-1641b4a6ebae");katex.render("\\Psi_{g} = \\Big(1 - \\tau_{g}\\Sigma_{r,\\, g}\\Big)\\Psi_{g} + \\tau_{g}\\Big(\\sum_{g' = 1}^{G}\\int_{4\\pi}\\Sigma_{s,\\, g'\\rightarrow g}f_{g'\\rightarrow g}\\Psi_{g'}\\, d\\Omega' + S_{g} - \\frac{\\partial}{\\partial t}\\frac{\\Psi_{g}}{v_{g}} - \\hat{\\Omega}\\cdot\\vec{\\nabla}\\Psi_{g}\\Big)", element, {displayMode:true,throwOnError:false});

where $var element = document.getElementById("moose-equation-fd8b71da-86bc-46d7-a512-a1b6479bd5f8");katex.render("\\tau_{g}", element, {displayMode:false,throwOnError:false});$ is defined as:

(4)var element = document.getElementById("moose-equation-ebfc34e2-c6a6-4990-a375-079d92f836f2");katex.render("\\tau_{g} = \\begin{cases} \\frac{1}{c\\Sigma_{r,\\,g}}, & ch\\Sigma_{r,\\,g} \\geq \\zeta\\\\ \\frac{h}{\\zeta}, & ch\\Sigma_{r,\\,g} < \\zeta \\end{cases}", element, {displayMode:true,throwOnError:false});

h is the maximum vertex separation within a single finite element, c and $var element = document.getElementById("moose-equation-c1fb1381-562e-4672-8cfa-1ce73d179c0a");katex.render("\\zeta", element, {displayMode:false,throwOnError:false});$ are constants to determine the maximum stabilization in a void region. Selecting $var element = document.getElementById("moose-equation-6e41711f-cb51-4d41-b025-6bee0e0e0eda");katex.render("\\zeta = 0", element, {displayMode:false,throwOnError:false});$ collapses Eq. (3) into Eq. (1). Substituting Eq. (3) into Governing Equations, Eq. (3) yields the void-treated SAAF weak form which is implemented in Gnat as the SAAF scheme:

(5)var element = document.getElementById("moose-equation-11d1eec2-4ea5-4cac-974b-2ab897c514de");katex.render("\\Big( \\psi_{j} + \\tau_{g}\\vec{\\nabla}\\psi_{j}\\cdot\\hat{\\Omega},\\, \\frac{\\partial}{\\partial t}\\frac{\\Psi^{h}_{g}}{v_{g}} \\Big)_{V} + \\Big\\langle \\psi_{j},\\, \\hat{n}\\cdot\\hat{\\Omega}\\Psi^{h}_{g} \\Big\\rangle_{\\Gamma} + \\Big( \\vec{\\nabla}\\psi_{j}\\cdot\\hat{\\Omega},\\, \\tau_{g}\\hat{\\Omega}\\cdot\\vec{\\nabla}\\Psi_{g}^{h} + (\\tau_{g}\\Sigma_{r,\\, g} - 1)\\Psi_{g}^{h} \\Big)_{V}\\\\ + \\Big( \\psi_{j},\\, \\Sigma_{r,\\,g}\\Psi^{h}_{g} \\Big)_{V} - \\Big( \\psi_{j} + \\tau_{g}\\vec{\\nabla}\\psi_{j}\\cdot\\hat{\\Omega},\\, \\sum_{g' = 1}^{G}\\int_{4\\pi}\\Sigma_{s,\\, g'\\rightarrow g}f_{g'\\rightarrow g}\\Psi^{h}_{g'}\\, d\\Omega' \\Big)_{V} - \\Big( \\psi_{j} + \\tau_{g}\\vec{\\nabla}\\psi_{j}\\cdot\\hat{\\Omega},\\, S_{g} \\Big)_{V} = 0", element, {displayMode:true,throwOnError:false});

2. The Upwinding Method

The upwinding approach for the neutron transport equation follows from the discontinuous weak form of the neutron transport equation:

(6)var element = document.getElementById("moose-equation-c59466e7-085a-45d4-999c-f9bee6461bbe");katex.render("\\Big( \\psi_{j},\\, \\frac{\\partial}{\\partial t}\\frac{\\Psi^{k}_{g}}{v_{g}} \\Big)_{D^{k}} + \\Big\\langle \\psi_{j},\\, \\hat{n}\\cdot\\hat{\\Omega}\\mathcal{F}^{*}_{g} \\Big\\rangle_{F^{k}} - \\Big( \\vec{\\nabla}\\psi_{j}\\cdot\\hat{\\Omega},\\, \\Psi^{k}_{g} \\Big)_{D^{k}}\\\\ + \\Big( \\psi_{j},\\, \\Sigma_{r,\\,g}\\Psi^{k}_{g} \\Big)_{D^{k}} - \\Big( \\psi_{j},\\, \\sum_{g' = 1}^{G}\\int_{4\\pi}\\Sigma_{s,\\, g'\\rightarrow g}f_{g'\\rightarrow g}\\Psi^{k}_{g'}\\,d\\Omega' \\Big)_{D^{k}} - \\Big( \\psi_{j},\\, S_{g} \\Big)_{D^{k}} = 0", element, {displayMode:true,throwOnError:false});

where:

$var element = document.getElementById("moose-equation-9758a5aa-2584-473b-aaa7-ca6206d4d926");katex.render("\\Psi^{k}_{g}", element, {displayMode:false,throwOnError:false});$ are the trial functions defined on the single finite element $var element = document.getElementById("moose-equation-ee755e7e-371a-4550-a95e-c573e8e0a27f");katex.render("k", element, {displayMode:false,throwOnError:false});$ . It is important to note that the trial functions are discontinuous over the full domain.
$var element = document.getElementById("moose-equation-710cc924-bd7e-49e2-bf56-6a987df3e9c9");katex.render("D^{k}", element, {displayMode:false,throwOnError:false});$ is the volume of the finite element $var element = document.getElementById("moose-equation-72c6bb02-1734-4338-aa00-66050034b210");katex.render("k", element, {displayMode:false,throwOnError:false});$ such that $var element = document.getElementById("moose-equation-e7457ec3-4ee8-42bb-b6ef-352a66513e42");katex.render("D^{k}\\in V", element, {displayMode:false,throwOnError:false});$ .
$var element = document.getElementById("moose-equation-85e0c6e5-f50e-4c35-ba57-dcb942623811");katex.render("F^{k}", element, {displayMode:false,throwOnError:false});$ is the bounding surface of the finite element $var element = document.getElementById("moose-equation-12e4718e-d212-44ce-a5a9-899cb9d9e34f");katex.render("k", element, {displayMode:false,throwOnError:false});$ .
$var element = document.getElementById("moose-equation-6aa047b7-f735-4d46-886d-052072cce3d2");katex.render("\\mathcal{F}^{*}_{g}(\\vec{r}) = \\Psi^{k}_{g}(\\vec{r}),\\, \\vec{r}\\in F^{k}", element, {displayMode:false,throwOnError:false});$ is the interface numerical angular flux.

The upwinding scheme is comparatively simple to the SAAF scheme. For each finite element $var element = document.getElementById("moose-equation-3b5d59dd-5cae-4926-9c60-94148aa67082");katex.render("k", element, {displayMode:false,throwOnError:false});$ the faces which compose the bounding surface $var element = document.getElementById("moose-equation-047bdc57-8ec6-4d64-be8c-088d42279204");katex.render("F^{k}", element, {displayMode:false,throwOnError:false});$ are traversed. Each face is classified as either upwind or downwind according to the following scheme:

Upwind faces are faces where $var element = document.getElementById("moose-equation-3d75c74f-8e98-45ec-a2d1-6553958f5537");katex.render("\\hat{n}\\cdot\\hat{\\Omega} \\geq 0", element, {displayMode:false,throwOnError:false});$ .
Downwpwind faces are faces where $var element = document.getElementById("moose-equation-27d29b34-97c6-4950-ae30-cc62434b2a89");katex.render("\\hat{n}\\cdot\\hat{\\Omega} \\lt 0", element, {displayMode:false,throwOnError:false});$ .

The numerical flux is then equal to the value of the neighbour's trial function at the interface for downwind faces, and the value of the current element's trial function at the interface.

Mass Transport

The stabilization approach taken for the isotope mass transport equation is the SUPG method (TODO: REFERENCE HERE). The SUPG method substitutes the set of continuous test functions for a discontinuous set of test functions which add artificial diffusion in the direction of advection. The derivation of the stabilized SUPG form begins with the discontinuous test functions:

(7)var element = document.getElementById("moose-equation-3996eb72-6ff5-4b4d-8d5c-a0bcb76fb7fd");katex.render("\\psi_{j} + \\tau\\vec{v}\\cdot\\vec{\\nabla}\\psi_{j},\\, \\tau = \\frac{h}{2||v||}", element, {displayMode:true,throwOnError:false});

where $var element = document.getElementById("moose-equation-cec3913d-6e0c-45f3-9133-c3f03176b567");katex.render("\\vec{v}", element, {displayMode:false,throwOnError:false});$ is the velocity and $var element = document.getElementById("moose-equation-9cc43067-6907-4e73-a98a-421bab4cf71a");katex.render("h", element, {displayMode:false,throwOnError:false});$ is the minimum separation between vertices in a single finite element. Rearranging Governing Equations, Eq. (6) such that it equals zero, multiplying by the modified test functions Eq. (7) and integrating over the domain $var element = document.getElementById("moose-equation-479b3c80-67cd-452d-be6d-ab5f3c8638fe");katex.render("V", element, {displayMode:false,throwOnError:false});$ yields the following:

(8)var element = document.getElementById("moose-equation-394da40b-3626-49d1-9c06-f3ea455ac1e8");katex.render("\\Big(\\psi_{j} + \\tau\\vec{v}\\cdot\\vec{\\nabla}\\psi_{j},\\,\\frac{\\partial}{\\partial t}N_{i}\\Big)_{V} + \\Big(\\psi_{j} + \\tau\\vec{v}\\cdot\\vec{\\nabla}\\psi_{j},\\,\\nabla\\cdot\\Big(\\vec{v}N_{i}\\Big)\\Big)_{V} + \\Big(\\psi_{j} + \\tau\\vec{v}\\cdot\\vec{\\nabla}\\psi_{j},\\,\\lambda_{i}N_{i}\\Big)_{V} + \\Big(\\psi_{j} + \\tau\\vec{v}\\cdot\\vec{\\nabla}\\psi_{j},\\,\\sum_{g = 1}^{G}\\sigma_{a,g,i}N_{i}\\Phi_{g}\\Big)_{V}\\\\ - \\Big(\\psi_{j} + \\tau\\vec{v}\\cdot\\vec{\\nabla}\\psi_{j},\\,\\nabla\\cdot\\Big(D_{i}\\vec{\\nabla}N_{i}\\Big)\\Big)_{V} - \\Big(\\psi_{j} + \\tau\\vec{v}\\cdot\\vec{\\nabla}\\psi_{j},\\,\\sum_{i' = 1}^{I}\\sum_{g = 1}^{G}\\sigma_{a,g,i'\\rightarrow i}N_{i'}\\Phi_{g}\\Big)_{V} - \\Big(\\psi_{j} + \\tau\\vec{v}\\cdot\\vec{\\nabla}\\psi_{j},\\,\\sum_{i' = 1}^{I}\\lambda_{i'}f_{i'\\rightarrow i}N_{i'}\\Big)_{V} = 0", element, {displayMode:true,throwOnError:false});

Applying integration by parts and the divergence theorem yields the final weak form for the isotope mass transport equation:

Neutron Transport
Mass Transport