FESTUNG - 3. 采用 HDG 方法求解对流问题[1]
1. 控制方程
线性对流问题控制方程为
[egin{array}{ll}
partial_t c +
abla cdot f = h, & mathrm{in} ; J imes Omega \
c(x, 0) = c_0(x), & mathrm{on} ; {0} imes Omega \
c(x, t) = c_D, & mathrm{on} ; J imes partial Omega_{in}
end{array}
]
其中 (f = [u^1, u^2]^ op 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}
ight) cdot
abla varphi_{h} mathrm{d} mathbf{x}+int_{partial T_{k}} hat{mathbf{f}}left(lambda_{h}, c_{h}^{-}
ight) cdot oldsymbol{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}^{-} ight)) 表达式为
[hat{mathbf{f}}left(lambda_{h}, c_{h}^{-}
ight) :=mathbf{f}left(lambda_{h}
ight)-alphaleft(lambda_{h}-c_{h}^{-}
ight) oldsymbol{v},
]
(lambda_h) 为计算边界上引入对应的未知解。为了对 (lambda_h) 进行求解，需要引入附加的边界积分方程
[int_{E_{overline{k}} in mathcal{E}_{ ext { int }}} alpha mu_{h}left(2 lambda_{h}-c_{h}^{-}-c_{h}^{+}
ight) mathrm{d} s+int_{E_{overline{k}} in mathcal{E}_{mathrm{bc}}} mu_{h}left(lambda_{h}-c_{partial Omega}left(c_{h}^{-}
ight)
ight) mathrm{d} s=0,
quad
forall mu_{h} in M_{h}.
]
将计算域划分为 (K) 个互不重叠的三角形单元，在单元和边界上分别定义多项式空间
[V_{h} :=left{varphi_{h} : overline{Omega}
ightarrow mathbb{R} : forall T in mathcal{T}_{h},left.varphi_{h}
ight|_{T} in mathbb{P}_{p}(T)
ight}, \
M_{h} :=left{mu_{h} : Gamma
ightarrow mathbb{R} : forall E in mathcal{E}_{h},left.mu_{h}
ight|_{E} in mathbb{P}_{p}(E)
ight}.
]
将 (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}.
]
将上述近似公式代入积分方程中，可得
[egin{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} ackslash 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} ackslash 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} ackslash 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{egin{array}{cc}{sum_{j=1}^{overline{N}} Lambda_{k j} mu_{overline{k} j}} & { ext { on } partial Omega_{ ext { out }}} \ {c_{mathrm{D}}(t, mathbf{x})} & { ext { on } partial Omega_{ ext { in }}}end{array}
ight} }_{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} ackslash partial Omega} mu_{overline{k} j} mu_{overline{k} i} }_{alpha ar{ mathcal{M} }_{mu} } Lambda_{k j}
- underbrace{ alpha sum_{j=1}^{N} int_{partial T_{k} ackslash 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} }_{ ilde{mathcal{M}}_{mu}} Lambda_{k j}
- underbrace{ int_{partial T_{k} cap partial Omega} mu_{overline{k} i}left{egin{array}{cc}{sum_{j=1}^{N} C_{k j} varphi_{k j}} & { ext { on } partial Omega_{ ext { out }}} \ {c_{D}(t, mathbf{x})} & { ext { on } partial Omega_{ ext { in }}}end{array}
ight} }_{- 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}
ight) }_{ar{mathcal{L}}} mathbf{C}+
underbrace{ left(mathcal{S}+mathcal{S}_{mathrm{out}}-alpha mathcal{R}_{mu}
ight) }_{ar{mathcal{M}}} mathbf{Lambda}=
underbrace{ mathcal{H}-mathcal{F}_{varphi, ext { in }} }_{mathcal{B}_{varphi}} \
underbrace{ left(-alpha mathcal{T}-mathcal{K}_{mu, ext { out }}
ight) }_{mathcal{N}} mathbf{C} +
underbrace{ left(alpha ar{mathcal{M}}_{mu}+ ilde{mathcal{M}}_{mu}
ight) }_{mathcal{P}} mathbf{Lambda}=
underbrace{ -mathcal{K}_{mu, ext { in }} }_{mathcal{B}_{mu}}
]
2.2. 时间离散
在获得关于时间的常微分方程，表达式为
[mathcal{M}_{varphi} partial_t mathbf{C} =
mathcal{B}_{varphi} - ar{mathcal{L}} mathbf{C} - ar{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 ar{mathcal{L}}
ight) }_{mathcal{L}} mathbf{C}^{(i)}
+ underbrace{ alpha_{ii} Delta t ar{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)} ight)) 代入，可得 (mathbf{Lambda}^{(i)}) 表达式为
[left( -mathcal{N} mathcal{L}^{-1} mathcal{M} + mathcal{P}
ight) 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. ↩︎