一阶中子输运方程稳定有限元方法

发表于 | 分类于

标准的Galerkin有限元方法只适用于椭圆方程，直接应用于双曲方程会产生伪震荡，所以需要一些技巧来消除这些震荡，使解变得稳定。这些技巧主要有两个方面：第一是将一阶的双曲方程变为二阶的椭圆方程，比如最小二乘方法（GLS)、迎风流线扩散方法(SUPG)等；另一种是将离散的格式写为迎风格式，比如间断有限元(DG)等，以下对这些稳定算法进行简要的介绍。

Galerkin

离散过程

输运方程
$$Ω⋅∇ϕ+σϕ = s$$

弱形式
$$(φ,Ω⋅∇ϕ)+(φ,σϕ) = (φ,s)$$

分部积分
$$<φ,Ω⋅nϕ> - (Ω⋅∇φ,ϕ)+(φ,σϕ) = (φ,s)$$

边界条件
$$<φ,Ω⋅nϕ>^{-} = G$$

带边界条件的方程
$$<φ,Ω⋅nϕ>^{+} - (Ω⋅∇φ,ϕ)+(φ,σϕ) = (φ,s) - G$$

等价方程1
$$ -<φ,Ω⋅nϕ>^{-} + (φ,Ω⋅∇ϕ)+(φ,σϕ) = (φ,s)-G$$

等价方程2
$$ [<φ,Ω⋅nϕ>^{+}-<φ,Ω⋅nϕ>^{-}] /2 + [(φ,Ω⋅∇ϕ)- (Ω⋅∇φ,ϕ)]/2+(φ,σϕ) = (φ,s)-G$$

矩阵形式

定义
$$L = (φ,Ω⋅∇φ)$$

$$R = (φ,σφ)$$

$$O = <φ,Ω⋅nφ>$$

$$O^{+} = <φ,Ω⋅nφ>^{+}$$

$$O^{-} = <φ,Ω⋅nφ>^{-}$$

$$b = (φ,s)-G$$

所以Galerkin离散后的方程为
$$[L - O^{-}+R]ϕ = b$$

$$[ -L^{T}+ O^{+}+R]ϕ = b$$

$$[ (L-L^{T})/2+(O^{+}-O^{-})/2+R]ϕ = b$$

SUPG (Streamline Upwinding Petrov Galerkin)

离散过程

流线扩散方法属于Petrov Galerkin 方法，该方法的检验函数与试探函数所在的空间不一样，检验函数所在的空间为 $φ+τΩ⋅∇φ$，所以输运方程的弱形式为
$$(φ,Ω⋅∇ϕ)+(φ,σϕ) + (τΩ⋅∇φ,Ω⋅∇ϕ)+(τΩ⋅∇φ,σϕ) = (φ,s)+ (τΩ⋅∇φ,s)$$

矩阵形式

定义
$$D=(τΩ⋅∇φ,Ω⋅∇ϕ)$$

$$C_{r}^{T} = (τΩ⋅∇φ,σϕ)$$

$$b = (φ,s)+ (τΩ⋅∇φ,s)-G$$

所以SUPG离散之后的方程为
$$[-L^{T}+ O^{+}+R + D + L_{r}^{T}]ϕ = b$$

对于非真空区域，取$τ = \frac{1}{σ}$，则$L_{r}^{T} = L^{T}$，所以此时
$$[ O^{+}+R + D ]ϕ = b$$

GLS (Galerkin Least Squares)

最小二乘公式的推导方式有很多种，可以从变分原理出发，也可以从加权残值法出发。以下从最小二乘泛函出发推导。令 $T = Ω⋅∇+σ$, 定义输运方程的全局残差的二范数

$$R = \int_v \int_{Ω} (T\phi - s)^2dΩdV$$

对输运方程求解就是要使残差最小，即该泛函的一阶变分为零，所以

$$δR = 2 \int_v \int_{Ω} Tδϕ(T\phi - s)dΩdV = 0$$

所以，最小二乘方法的弱形式为

$$(Tδϕ,Tϕ) = (Tδϕ,s)$$

令 $φ = δϕ$，所以

$$(Ω⋅∇φ,Ω⋅∇ϕ) + (σφ,Ω⋅∇ϕ) + (Ω⋅∇φ,σϕ) + (σφ,σϕ)= (Ω⋅∇φ,s)+ (σφ,s)$$

边界条件

最小二乘方法的边界条件不好处理，对于$S_N$方法可以全部使用强制边界条件；而对于$P_N$方法，需要使用弱形式的边界条件，即在原来变分形式中减去边界上的变分的 $c$ 倍，即

$$(Ω⋅∇φ,Ω⋅∇ϕ) + (σφ,Ω⋅∇ϕ) + (Ω⋅∇φ,σϕ) + (σφ,σϕ) - ^{-} = (Ω⋅∇φ,s)+ (σφ,s)$$

其中c是一个常数，c的选取是一个很重要的问题。

SGS (Sub-Grid Scale)

离散过程

该方法将变量分为连续部分与间断部分之和，其中$ ϕ_c$为连续部分，$ϕ_d$为间断部分，最终会形成两个方程
$$ϕ = ϕ_c+ϕ_d$$

其弱形式为
$$(φ_c,Ω⋅∇ϕ)+(φ_c,σϕ) = (φ_c,s) \tag{$1$}$$

$$(φ_d,Ω⋅∇ϕ)+(φ_d,σϕ) = (φ_d,s) \tag{$2$}$$

对方程1的$ϕ$分部积分
$$<φ_c,Ω⋅nϕ> - (Ω⋅∇φ_c,ϕ)+(φ_c,σϕ) = (φ_c,s)$$

对方程2的$ϕ_d$分部积分
$$<φ_d,Ω⋅nϕ_d> - (Ω⋅∇φ_d,ϕ_d)+(φ_d,Ω⋅∇ϕ_c)+(φ_d,σϕ) = (φ_d,s)$$

对于区域边界，有如下边界条件
$$<φ,Ω⋅nϕ_c>^{-} = G$$

$$<φ,Ω⋅nϕ_d> = 0$$

对于内部单元边界
$$<φ_d,Ω⋅nϕ_d>^{-} = 0$$

带入边界条件得
$$<φ_c,Ω⋅nϕ_c>^{+} - (Ω⋅∇φ_c,ϕ)+(φ_c,σϕ) = (φ_c,s) - G$$

$$<φ_d,Ω⋅nϕ_d>^{+} - (Ω⋅∇φ_d,ϕ_d)+(φ_d,Ω⋅∇ϕ_c)+(φ_d,σϕ) = (φ_d,s)$$

矩阵形式

$$[ O^{+}-L^{T}+R]ϕ_c+[ -L^{T}+R]ϕ_d = s_c \tag{$3$}$$

$$[ L+R]ϕ_c+[ O^{+}-L^{T}+R]ϕ_d = s_d \tag{$4$}$$

方程(3)描述整个求解区域，方程(4)描述单个单元。所以对于每个内部单元

$$K_{11} = -L^{T}+R$$

$$K_{12} = -L^{T}+R$$

$$K_{21} = L+R$$

$$K_{22} = O^{+}-L^{T}+R$$

对于边界
$$K_{11} = O^{+}$$

由方程(4)求解$\phi_d$
$$\phi_d=K_{22}^{-1}(s_d - K_{21}\phi_c)$$

然后带入方程(3)中，得
$$[K_{11}-K_{12}K_{22}^{-1}K_{21}]\phi_c = s_c - K_{12}K_{22}^{-1}s_d$$

迭代求解

采用迭代方法求解

  1. 令$\phi_d = 0$，求解$K_{11}\phi_c = s_c$，得到$\phi_c$
  2. 计算每个单元上的误差源 $s_d = s_d - K_{21}\phi_c$
  3. 在每个单元上求解 $K_{22}\phi_d = s_d$，得到 $\phi_d$
  4. 更新源项 $s_c =s_c - K_{12}\phi_d$
  5. 求解方程 $K_{11}\phi_c = s_c$，得到新的 $\phi_c$

感觉迭代 2 次应该就差不多了。

结果

迭代求解失败，方程不满足局部守恒关系。迭代会震荡得更厉害。

SUPG （LS） + SGS

想采用LS方法和SGS方法结合起来求解。因为最小二乘方法中为处理边界条件，公式中出现 $\frac{1}{σ}$，所以不适合于真空边界条件求解。所以单元矩阵在分真空边界上采用最小二乘方法，在真空区域采用SGS方法。

以下方法是我自己提出来的

FSGS (快速SGS)

最终要求解的代数方程组为
$$K\phi = F$$
其中 $F = Es$，所以
$$K = A-BD^{-1}C = -L^{T}+R + O^{+} -(-L^{T}+R)D^{-1}(L+R)$$

$$K = -L^{T}+R+ O^{+} + L^{T}D^{-1}L + L^{T}D^{-1}R -RD^{-1}L - RD^{-1}R$$

$$E = I-(-L^{T}+R)D^{-1} = I+L^{T}D^{-1} - RD^{-1}$$

非真空区域

令 $D=R$，此时
$$K = -L^{T}+R+ O^{+} + L^{T}R^{-1}L + L^{T} -L - R = O^{+} + L^{T}R^{-1}L -L$$

$$E = L^{T}R^{-1}$$

真空区域

令 $D = I$，此时

$$K = O^{+} + L^{T}L -L $$

$$E = L^{T}$$

OWLS (最优加权最小二乘)

不带边界条件连续有限元离散之后的输运方程为
$$(L+R)\phi = s$$

在方程两端同时乘以 $(L+R)^{T}D$，$D$是一个权重矩阵，所以
$$(L+R)^{T}D(L+R)\phi = (L+R)^{T}Ds$$

所以
$$K =L^{T}DL + L^{T}DR+R^{T}DL + R^{T}DR$$

$$E =L^{T}D+ R^{T}D$$

非真空区域

取$D = R^{-1}$，由于矩阵 $R$ 是对称的，所以
$$K =L^{T}R^{-1}L + L^{T}+L + R$$

$$E =L^{T}R^{-1}+ I$$

该方程减去一阶方程$(L+R)\phi = s$，得
$$K =L^{T}R^{-1}L + L^{T}$$

$$E =L^{T}R^{-1}$$

由于 $L^{T} = O-L$，所以
$$K =L^{T}R^{-1}L + O-L$$

$$E =L^{T}R^{-1}$$

带入边界条件后，该式与FSGS是完全一样的

真空区域

取 $D = I$，则
$$K = O^{+} + L^{T}L -L $$

$$E = L^{T}$$

所以，最小二乘方法与FSGS是完全等价的。

为什么要减去一阶方程

不带边界条件的最小二乘方法的弱形式为

$$(Tφ,Tϕ) = (Tφ,s)$$

利用分部积分将检验函数上的导数都转移到求解量上，得
$$<φ,n⋅Ω(Ω⋅∇ϕ+σϕ-s)>-(φ,Ω⋅∇(Ω⋅∇ϕ)-Ω⋅∇s-σ^{2}ϕ+σs) = 0 $$

所以最小二乘泛函的Euler-Lagrange方程为
$$\left\{\begin{matrix}
Ω⋅∇(Ω⋅∇ϕ)-Ω⋅∇s-σ^{2}ϕ+σs=0,\phi \in V\\
Ω⋅∇ϕ+σϕ-s,\phi \in Γ
\end{matrix}\right.$$

所以在一个非真空材料区域内满足SAAF方程，在材料边界满足一阶中子输运方程。如果直接对原来的最小二乘变分求解表示一阶方程在材料边界上严格满足，不存在离散误差，而这是不可能的，在材料交界面上一阶方程不能严格满足。所以原来的最小二乘变分必须减去一阶方程。这样做有很多好处，包括实现了最小二乘方程的全局守恒性，可以直接处理自然边界条件，可以计算内真空问题，解决了传统最小二乘方法的所有问题。