FESTUNG 模型介绍 - 2. 对流问题隐式求解
1. 控制方程
对流问题的控制方程为
[partial_t C + partial_x u^1 C + partial_y u^2 C = 0, \
egin{array}{cl}
C = C_D & mathrm{on} ; partial Omega_D, \
-
abla C cdot mathbf{v} = g_N & mathrm{on} ; partial Omega_N.
end{array}
]
在边界处包含 Dirichlet 和 Neumann 两种边界条件。
2. 数值离散
2.1. 空间离散
采用间断有限元方法对控制方程进行离散，包括采用一次分部积分并利用 Green-Gauss 公式将方程转化为
[underbrace{ int_{mathcal{T}_k} varphi_{ki} sum_{j=1}^N partial_t C_{kj} varphi_{kj} }_{I}
- underbrace{ int_{mathcal{T}_k} partial_x varphi_{ki} sum_{j=1}^N C_{kj} varphi_{kj} sum_{j=1}^N u^1_{kl} varphi_{kl}
- int_{mathcal{T}_k} partial_y varphi_{ki} sum_{j=1}^N C_{kj} varphi_{kj} sum_{j=1}^N u^1_{kl} varphi_{kl} }_{II} \
+ underbrace{ left {
egin{array}{c}
sum_{n=1}^3 int_{partial mathcal{T}_k} v_{kn}^1 u^1_{kn} varphi_{ki} sum_{j=1}^N C_{kj}^* varphi_{kj}
+ sum_{n=1}^3 int_{partial mathcal{T}_k} v_{kn}^2 u^2_{kn} varphi_{ki} sum_{j=1}^N C_{kj}^* varphi_{kj} quad mathrm{on} ; mathcal{E}_{Omega} \
sum_{n=1}^3 int_{partial mathcal{T}_k} v_{kn}^1 u^1_{kn} varphi_{ki} sum_{j=1}^N C_D^* varphi_{kj}
+ sum_{n=1}^3 int_{partial mathcal{T}_k} v_{kn}^2 u^2_{kn} varphi_{ki} sum_{j=1}^N C_D^* varphi_{kj} quad mathrm{on} ; mathcal{E}_{D}
end{array}
ight }
}_{III} = 0.
]
代入质量矩阵 (mathcal{M}) 等，可将半离散方程化为矩阵形式
[mathcal{M} partial_t C + left( mathcal{G}^1 + mathcal{G}^2 + mathcal{R}
ight) C = mathcal{K}_D + mathcal{K}_N
]
其中 (mathcal{K}_D) 和 (mathcal{K}_N) 为边界条件产生的常数向量。最终得到关于时间的常微分方程组
[partial_t C = mathcal{S}( C, t ), quad mathcal{S}( C, t ) = mathcal{M}^{-1} left[ mathcal{K}_D + mathcal{K}_N - left( mathcal{G}^1 + mathcal{G}^2 + mathcal{R}
ight) C
ight]
]
2.2. 时间离散
在获得关于时间的半离散方程后，可采用 s 步对角隐式 Runge-Kutta 方法(DIRK)计算，其计算流程为
[egin{array}{ll}
mathbf{C}^{(i)} :=mathbf{C}^{n}+Delta t^{n} sum_{j=1}^{i} a_{i j} mathcal{S}left(mathbf{C}^{(j)}, t^{(j)}
ight), & ext { for } quad i=1, ldots, s \
mathbf{C}^{n+1} :=mathbf{C}^{n}+Delta t^{n} sum_{i=1}^{s} b_{i} mathcal{S}left(mathbf{C}^{(i)}, t^{(i)}
ight)
end{array}
]
在 DIRK 方法中，(b_j = a_{sj})，因此计算时不必进行最后一步 (mathbf{C}^{n+1}) 的计算，而只需把 (mathbf{C}^{(s)}) 值赋给 (mathbf{C}^{n+1}) 即可。
将空间算子表达式代入，可得每一步 (mathbf{C}^{(i)}) 的具体计算公式为
[mathbf{C}^{(i)} =mathbf{C}^{n}+Delta t^{n} sum_{j=1}^{i} a_{i j} mathcal{M}^{-1} left[ mathcal{K}_D + mathcal{K}_N - left( mathcal{G}^1 + mathcal{G}^2 + mathcal{R}
ight) C^{(j)}
ight]
]
将与 (mathbf{C}^{(i)}) 相关的项移到等号左侧可得
[left( 1 + Delta t^n a_{ii} mathcal{M}^{-1} mathcal{A}
ight) mathbf{C}^{(i)} = mathbf{C}^{n} + Delta t^n a_{ii} mathcal{M}^{-1} mathcal{V} + sum_{j=1}^{i-1} a_{i j} mathcal{M}^{-1}( mathcal{V} - mathcal{A} mathbf{C}^{(j)} )
]
其中 (mathcal{A} = mathcal{G}^1 + mathcal{G}^2 + mathcal{R})， (mathcal{V} = mathcal{K}_D + mathcal{K}_N)。将 (mathcal{M}) 同时乘以等号两侧，可得
[left( mathcal{M} + Delta t^n a_{ii} mathcal{A}
ight) mathbf{C}^{(i)} =
mathcal{M} mathbf{C}^{n} + Delta t^n a_{ii} mathcal{V} + sum_{j=1}^{i-1} a_{i j} ( mathcal{V} - mathcal{A} mathbf{C}^{(j)} )
]
两侧同乘以 (left( mathcal{M} + Delta t^n a_{ii} mathcal{A}
ight)^{-1})，即可获得 (mathbf{C}^{(i)}) 的数值解。
3. 模型实现
在 FESTUNG 中，solveSubStep.m 实现了 DIRK 方法中每步 (mathbf{C}^{(i)}) 的计算过程。
首先，组装系数矩阵 sysA 与 sysV，分别代表 (mathcal{A}) 与 (mathcal{V})。当求解恒定问题时，由于此时半离散方程不包含时间变化项，因此
[C = mathcal{A}^{-1} mathcal{V}
]
当求解问题为非恒定时，首先计算方程组的右端项 sysR = (mathcal{M} mathbf{C}^{n} + Delta t^n a_{ii} mathcal{V} + sum_{j=1}^{i-1} a_{i j} ( mathcal{V} - mathcal{A} mathbf{C}^{(j)} )) ，其代码为
% Computing the rhs
sysR = (problemData.tau * problemData.A(nSubStep, nSubStep)) * sysV + problemData.globM * problemData.cDiscRK{1};
for j = 1 : nSubStep - 1
sysR = sysR + (problemData.tau * problemData.A(nSubStep, j)) * problemData.rhsRK{j};
end % for
随后，将系数矩阵的逆左乘至右端项 sysR 上，即可获得 (mathbf{C}^{(i)}) 的数值解
% Compute the next step
problemData.cDiscRK{nSubStep + 1} = (problemData.globM + (problemData.tau * problemData.A(nSubStep, nSubStep)) * sysA) sysR;