FESTUNG - 3. 采用 HDG 方法求解对流问题1
1. 控制方程
线性对流问题控制方程为
\[
\begin{array}{ll}
\partial_t c + \nabla \cdot f = h, & \mathrm{in} \; J \times \Omega \\
c(x, 0) = c_0(x), & \mathrm{on} \; {0} \times \Omega \\
c(x, t) = c_D, & \mathrm{on} \; J \times \partial \Omega_{in}
\end{array}
\]
其中 \(f = [u^1, u^2]^\top c\) 为通量项。
2. 数值离散
2.1. 空间离散
采用间断有限元方法对控制方程进行离散，包括进行一次分部积分并采用 Green-Gauss 公式导出面积分项，可得
\[
\int_{T_{k}} \partial_{t} c_{h} \varphi_{h} \mathrm{d} \mathbf{x}-\int_{T_{k}} \mathbf{f}\left(c_{h}\right) \cdot \nabla \varphi_{h} \mathrm{d} \mathbf{x}+\int_{\partial T_{k}} \hat{\mathbf{f}}\left(\lambda_{h}, c_{h}^{-}\right) \cdot \boldsymbol{v} \varphi_{h}^{-} \mathrm{d} s=\int_{T_{k}} h \varphi_{h} \mathrm{d} \mathbf{x},
\quad \forall \varphi_{h} \in V_{h}
\]
其中面积分中引入近似通量项 \(\hat{\mathbf{f}}\left(\lambda_{h}, c_{h}^{-}\right)\) 表达式为
\[
\hat{\mathbf{f}}\left(\lambda_{h}, c_{h}^{-}\right) :=\mathbf{f}\left(\lambda_{h}\right)-\alpha\left(\lambda_{h}-c_{h}^{-}\right) \boldsymbol{v},
\]
\(\lambda_h\) 为计算边界上引入对应的未知解。为了对 \(\lambda_h\) 进行求解，需要引入附加的边界积分方程
\[
\int_{E_{\overline{k}} \in \mathcal{E}_{\text { int }}} \alpha \mu_{h}\left(2 \lambda_{h}-c_{h}^{-}-c_{h}^{+}\right) \mathrm{d} s+\int_{E_{\overline{k}} \in \mathcal{E}_{\mathrm{bc}}} \mu_{h}\left(\lambda_{h}-c_{\partial \Omega}\left(c_{h}^{-}\right)\right) \mathrm{d} s=0,
\quad
\forall \mu_{h} \in M_{h}.
\]
将计算域划分为 \(K\) 个互不重叠的三角形单元，在单元和边界上分别定义多项式空间
\[
V_{h} :=\left\{\varphi_{h} : \overline{\Omega} \rightarrow \mathbb{R} : \forall T \in \mathcal{T}_{h},\left.\varphi_{h}\right|_{T} \in \mathbb{P}_{p}(T)\right\}, \\
M_{h} :=\left\{\mu_{h} : \Gamma \rightarrow \mathbb{R} : \forall E \in \mathcal{E}_{h},\left.\mu_{h}\right|_{E} \in \mathbb{P}_{p}(E)\right\}.
\]
将 \(c_h\) 和 \(\lambda_h\) 分别用单元和边界上多项式进行近似
\[
c_h = \sum_{j=1}^N c_{kj} \varphi_{kj}, \quad \lambda_h = \sum_{j=1}^{\hat{N}} \Lambda_{\hat{k} j} \mu_{\hat{k} j}.
\]
将上述近似公式代入积分方程中，可得
\[
\begin{array}{r}
\underbrace{ \sum_{j=1}^{N} \int_{T_{k}} \varphi_{k i} \varphi_{k j} }_{\mathcal{M}_{\varphi}} \partial_{t} C_{k j} -
\underbrace{ \sum_{j=1}^{N} \sum_{m=1}^{2} \int_{T_{k}} u^{m}(t, \mathbf{x}) \partial_{x^{m}} \varphi_{k i} \varphi_{k j} }_{- \mathcal{G}^2 - \mathcal{G}^2} C_{k j} \\
+ \underbrace{ \sum_{j=1}^{\overline{N}} \int_{\partial T_{k} \backslash \partial \Omega} \varphi_{k i} (\mathbf{u}(t, \mathbf{x}) \cdot v) \mu_{\overline{k} j} }_{\mathcal{S}} \Lambda_{k j}
- \underbrace{ \alpha \sum_{j=1}^{\overline{N}} \int_{\partial T_{k} \backslash \partial \Omega} \varphi_{k i} \mu_{\overline{k} j} }_{\alpha \mathcal{R}_{\mu}} \Lambda_{k j} \\
+ \alpha \underbrace{ \sum_{j=1}^{N} \int_{\partial T_{k} \backslash \partial \Omega} \varphi_{k i} \varphi_{k j} }_{\alpha \mathcal{R}_{\varphi}} C_{k j} \\
+ \underbrace{ \int_{\partial T_{k} \wedge \partial \Omega} \varphi_{k i}(\mathbf{u}(t, \mathbf{x}) \cdot v)\left\{\begin{array}{cc}{\sum_{j=1}^{\overline{N}} \Lambda_{k j} \mu_{\overline{k} j}} & {\text { on } \partial \Omega_{\text { out }}} \\ {c_{\mathrm{D}}(t, \mathbf{x})} & {\text { on } \partial \Omega_{\text { in }}}\end{array}\right\} }_{\mathcal{S}_{out}, \mathcal{F}_{\varphi, in}}
= \underbrace{ \int_{T_{k}} h \varphi_{k i} \mathrm{dx} }_{\mathcal{H}}
\end{array}
\]
关于 \(\lambda_h\) 的控制方程为
\[
\underbrace{ \alpha \sum_{j=1}^{\overline{N}} \int_{\partial T_{k} \backslash \partial \Omega} \mu_{\overline{k} j} \mu_{\overline{k} i} }_{\alpha \bar{ \mathcal{M} }_{\mu} } \Lambda_{k j}
- \underbrace{ \alpha \sum_{j=1}^{N} \int_{\partial T_{k} \backslash \partial \Omega} \mu_{\overline{k} i} \varphi_{k j} C_{k j} }_{- \alpha \mathcal{T} } \\
+ \underbrace{ \sum_{j=1}^{\overline{N}} \int_{\partial T_{k} \cap \partial \Omega} \mu_{\overline{k} j} \mu_{k i} }_{\tilde{\mathcal{M}}_{\mu}} \Lambda_{k j}
- \underbrace{ \int_{\partial T_{k} \cap \partial \Omega} \mu_{\overline{k} i}\left\{\begin{array}{cc}{\sum_{j=1}^{N} C_{k j} \varphi_{k j}} & {\text { on } \partial \Omega_{\text { out }}} \\ {c_{D}(t, \mathbf{x})} & {\text { on } \partial \Omega_{\text { in }}}\end{array}\right\} }_{- \mathcal{K}_{\mu, out}, \mathcal{K}_{\mu, in}} =0
\]
将离散方程简写为矩阵形式
\[
\mathcal{M}_{\varphi} \partial_{t} \mathbf{C}+
\underbrace{ \left(-\mathcal{G}^{1}-\mathcal{G}^{2}+\alpha \mathcal{R}_{\varphi}\right) }_{\bar{\mathcal{L}}} \mathbf{C}+
\underbrace{ \left(\mathcal{S}+\mathcal{S}_{\mathrm{out}}-\alpha \mathcal{R}_{\mu}\right) }_{\bar{\mathcal{M}}} \mathbf{\Lambda}=
\underbrace{ \mathcal{H}-\mathcal{F}_{\varphi, \text { in }} }_{\mathcal{B}_{\varphi}} \\
\underbrace{ \left(-\alpha \mathcal{T}-\mathcal{K}_{\mu, \text { out }}\right) }_{\mathcal{N}} \mathbf{C} +
\underbrace{ \left(\alpha \bar{\mathcal{M}}_{\mu}+ \tilde{\mathcal{M}}_{\mu}\right) }_{\mathcal{P}} \mathbf{\Lambda}=
\underbrace{ -\mathcal{K}_{\mu, \text { in }} }_{\mathcal{B}_{\mu}}
\]
2.2. 时间离散
在获得关于时间的常微分方程，表达式为
\[
\mathcal{M}_{\varphi} \partial_t \mathbf{C} =
\mathcal{B}_{\varphi} - \bar{\mathcal{L}} \mathbf{C} - \bar{\mathcal{M}} \mathbf{\Lambda} =
\mathcal{R}_{RK}(t, \mathbf{C}, \mathbf{\Lambda}) \\
0 = \mathcal{B}_{\mu} - \mathcal{N} \mathbf{C} - \mathcal{P} \mathbf{\Lambda}
\]
采用 DIRK 格式对其进行求解，其每分部表达式为
\[
\mathcal{M}_{\varphi} \mathbf{C}^{(i)} = \mathcal{M}_{\varphi} \mathbf{C}^n
+ \Delta t \sum_{j=1}^i \alpha_{ij} \mathcal{R}_{RK} (t^{(j)}, \mathbf{C}^{(j)}, \mathbf{\Lambda}^{(j)}) \\
0 = \mathcal{B}_{\mu}^{(i)} - \mathcal{N} \mathbf{C}^{(i)} - \mathcal{P} \mathbf{\Lambda}^{(i)}
\]
将 \(\mathcal{R}_{RK}\) 表达式代入，并将包含 \(\mathbf{C}^{(i)}\) 和 \(\mathbf{\Lambda}^{(i)}\) 的项移到等号左侧，可得
\[
\underbrace{ \left( \mathcal{M}_{\varphi} + \alpha_{ii} \Delta t \bar{\mathcal{L}} \right) }_{\mathcal{L}} \mathbf{C}^{(i)}
+ \underbrace{ \alpha_{ii} \Delta t \bar{\mathcal{M}}^{(i)} }_{\mathcal{M}} \mathbf{\Lambda}^{(i)}
= \underbrace{ \mathcal{M}_{\varphi} \mathbf{C}^n
+ \alpha_{ii} \Delta t \mathcal{B}_{\varphi}^{(i)} +
\Delta t \sum_{j=1}^{i-1} \alpha_{ij} \mathcal{R}_{RK} (t^{(j)}, \mathbf{C}^{(j)}, \mathbf{\Lambda}^{(j)}) }_{\mathcal{Q}} \\
\mathcal{N} \mathbf{C}^{(i)} + \mathcal{P} \mathbf{\Lambda}^{(i)} = \mathcal{B}_{\mu}^{(i)}
\]
即
\[
\mathcal{L} \mathbf{C}^{(i)} + \mathcal{M} \mathbf{\Lambda}^{(i)} = \mathcal{Q} \\
\mathcal{N} \mathbf{C}^{(i)} + \mathcal{P} \mathbf{\Lambda}^{(i)} = \mathcal{B}_{\mu}^{(i)}
\]
为了求解此方程组，将第一个方程中 \(\mathbf{C}^{(i)} = \mathcal{L}^{-1} \left( \mathcal{Q} - \mathcal{M} \mathbf{\Lambda}^{(i)} \right)\) 代入，可得 \(\mathbf{\Lambda}^{(i)}\) 表达式为
\[
\left( -\mathcal{N} \mathcal{L}^{-1} \mathcal{M} + \mathcal{P} \right) \mathbf{\Lambda}^{(i)} = \mathcal{B}_{\mu}^{(i)} - \mathcal{N} \mathcal{L}^{-1} \mathcal{Q}
\]
此时，即可获得每步的 \(\mathbf{\Lambda}^{(i)}\) 与 \(\mathbf{C}^{(i)}\) 数值解。
Jaust A, Reuter B, Aizinger V, Schütz J, Knabner P. FESTUNG: A MATLAB / GNU Octave toolbox for the discontinuous Galerkin method. Part III: Hybridized discontinuous Galerkin (HDG) formulation. ArXiv:170904390 [Math] 2017.↩