在对流体的数值模拟中,有两种典型的数值方法:拉格朗日方法和欧拉方法.欧拉法能够处理大变形问题,但是由于输运计算会带来误差,同时很难给出精确的材料界面.拉格朗日方法的网格完全跟随流体运动,能够用网格边界清晰地捕捉材料界面;但是当涉及大变形时,网格会产生扭曲,需要结合网格重分和物理量重映来保证模拟的正常实施.

考虑到拉格朗日方法能很好地计算流体界面,在可压缩流体流动的数值模拟中能够很好地处理间断问题,同时非结构三角形网格在处理复杂几何问题时具有明显的优势,本文给出一种使用RKDG方法在非结构网格上求解的拉格朗日格式.本文格式是一种半拉格朗日格式,格式中变量的梯度算子和散度算子采用欧拉坐标表示,时间导数采用拉格朗日坐标表示,从而省去了求解映射雅可比矩阵的步骤,简化了计算过程.在本文格式中,网格节点跟随流体运动,节点速度的求解方法是文献[4]中四边形网格节点速度求解方法向三角形网格的延伸,算法依据文献[5]中的相关内容进行改进,格式选取顶点沿边的熵波速度来代替文献[4]算法中的边法向速度的Roe平均.为了保证模拟的稳定实施,本文利用非结构三角形网格上的HWENO(Hermite weighted essential non-oscillatory)重构来控制间断附近产生的虚假振荡.一些数值算例验证了格式在非结构网格上的二阶精度和鲁棒性.

本文的基本内容如下:第2节给出了欧拉方程组的拉格朗日形式.第3节给出非结构三角形网格上使用RKDG方法求解拉氏(拉格朗日)方程组弱形式的基本过程.第4节给出的算例验证了算法的性能.

在本节,参照文献[6]的做法,我们首先给出了可压缩欧拉方程组的拉格朗日形式,同时给出用于求解的积分弱形式.欧拉框架下的无源可压缩流体力学方程组为

(2.1) $ \begin{equation} \left\{ \begin{array}{ll} { } \frac{\partial (\rho)}{\partial t}+\nabla \cdot (\rho {\bf V} ) = 0, \\ { } \frac{\partial (\rho {\bf V})}{\partial t}+\nabla \cdot (\rho {\bf V} \otimes{\bf V})+\nabla p = 0, \\ { } \frac{\partial (\rho E)}{\partial t}+\nabla \cdot (\rho E{\bf V} )+\nabla \cdot (p {\bf V} ) = 0, \end{array} \right. \end{equation} $

其中,未知量为流体密度$ \rho $,压力$ p $,速度$ {\bf V} = (u, v) $,总能量$ E $,内能为$ {\cal E} = E -\frac{1}{2}\parallel V\parallel^2 $.为保证系统的封闭性,增添气体状态方程: $ p = (\gamma-1)\rho{\cal E}. $

设可压缩流体在初始时刻的区域为$ D $,拉式框架下的坐标为$ (X, Y) $.对于任意时刻,拉式框架下的流体质团对应的欧拉坐标为$ (x, y) $.已知两种坐标的对应关系可由流体质团的轨迹方程表示

$ \frac{{\rm d}x}{{\rm d}t} = u(x, y, t), \ x(X, Y, 0) = X;\ \frac{{\rm d}y}{{\rm d}t} = v(x, y, t), \ y(X, Y, 0) = Y, $

其中$ \frac{\rm d}{{\rm d}t} $表示物质导数, $ (u, v) $为流体质团的速度.理想情况下,如果速度场非常光滑,轨迹方程存在唯一解$ (x(X, Y, t), y(X, Y, t)) $.进而,我们能够定义两种坐标之间存在的映射: $ {\cal A}_t:(X, Y)\longrightarrow(x, y), $其中, $ (x, y) $是轨迹方程的解.当目标质团运动到时刻$ t > 0 $,映射能够将流体区域由拉式坐标映射到欧拉坐标.对于初始时刻拉格朗日区域$ D $,设它在$ t $时刻欧拉框架下对应的区域为: $ {\cal A}_t(D) = {\cal D}. $另外,我们假设:在任意时间$ t > 0 $，映射$ {\cal A}_t $都可逆.则$ {\cal A}_t $的雅各比矩阵可表示为

$ J(X, Y, t) = \left|\begin{array}{cccc} { } \frac{\partial x}{\partial X}{\quad} &{ } \frac{\partial y}{\partial X} \\ { } \frac{\partial x}{\partial Y}{\quad} & { } \frac{\partial y}{\partial Y} \end{array}\right|. $

由于映射可逆,且$ t = 0 $时, $ J = 1; $$ t > 0 $时, $ J > 0 $,易得雅各比矩阵J的时间导数为

(2.2) $ \begin{equation} \frac{{\rm d}J}{{\rm d}t}-J \nabla \cdot {\bf V} = 0, \end{equation} $

其中$ \nabla \cdot {\bf V} = \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y} $表示速度场的散度.

设一个依赖欧拉坐标的流场物理量: $ \Psi = \Psi(x, y, t) $,它的时间导数为

(2.3) $ \begin{equation} \frac{{\rm d}\Psi}{{\rm d}t} = \frac{\partial \Psi}{\partial t}+\frac{\partial \Psi}{\partial x}\cdot\frac{\partial x}{\partial t}+\frac{\partial \Psi}{\partial y}\cdot\frac{\partial y}{\partial t}. \end{equation} $

结合轨迹方程,可得

(2.4) $ \begin{equation} \frac{{\rm d}\Psi}{{\rm d}t} = \frac{\partial \Psi}{\partial t}+\nabla\Psi \cdot {\bf V}, \end{equation} $

其中$ \nabla \Psi = (\frac{\partial \Psi}{\partial x}, \frac{\partial \Psi}{\partial y}) $为物理量梯度.最后,结合(2.2)和(2.4)式得到

(2.5) $ \begin{equation} \frac{{\rm d}(\Psi J)}{{\rm d}t} = J(\frac{\partial \Psi}{\partial t}+\nabla(\Psi {\bf V})). \end{equation} $

将守恒变量$ \Psi = \rho, \rho{\bf V}, \rho E $分别代入(2.5)式,并且结合原欧拉方程组可得方程组对应的拉格朗日形式

(2.6) $ \begin{equation} \left\{ \begin{array}{ll} { } \frac{{\rm d}(\rho J)}{{\rm d}t} = 0, \\ { } \frac{{\rm d}(\rho {\bf V} J)}{{\rm d}t}+J\nabla P = 0, \\ { } \frac{{\rm d}(\rho E J)}{{\rm d}t}+J\nabla\cdot (P{\bf V}) = 0. \end{array} \right. \end{equation} $

对于任意的测试函数$ \Phi(X, Y)\in L_2(D) $,有: $ \frac{{\rm d}(\Phi(X, Y))}{{\rm d}t} = 0. $进而,下列等式显然成立

(2.7) $ \begin{equation} \Phi\frac{{\rm d}(J\Psi)}{{\rm d}t} = \frac{{\rm d}(\Psi J\Phi)}{{\rm d}t}. \end{equation} $

在空间$ L_2(D) $中选取有限元检验函数$ \Phi(X, Y) $,将方程组(2.6)两侧同时乘以$ \Phi(X, Y) $,并在拉氏区域$ D $上积分得

(2.8) $ \begin{equation} \left\{ \begin{array}{ll} { } \frac{\rm d}{{\rm d}t} \int_{D}(\rho J\Phi) {\rm d}(D) = 0, \\ { } \frac{\rm d}{{\rm d}t} \int_{D}(\rho {\bf V} J\Phi) {\rm d}(D)+\int_{D}(J\nabla p \Phi) {\rm d}(D) = 0, \\ { } \frac{\rm d}{{\rm d}t} \int_{D}(\rho E J \Phi) {\rm d}(D)+\int_{D}(J\nabla \cdot(p{\bf V})\Phi) {\rm d}(D) = 0. \end{array} \right. \end{equation} $

由于映射$ {\cal A}_t $在任意时刻都可逆,则: $ \Phi(X, Y) = \Phi({\cal A}^{-1}_t(x, y)) = \Phi\cdot{\cal A}^{-1}_t(x, y), $进而: $ \Phi\cdot{\cal A}^{-1}_t(x, y)\in L_2({\cal D}). $设$ \varphi(x, y) = \Phi\cdot{\cal A}^{-1}_t(x, y) = \Phi(X, Y) $,那么在欧拉框架下的基函数需满足这样的性质

(2.9) $ \begin{equation} \frac{{\rm d}(\varphi(x, y))}{{\rm d}t} = 0. \end{equation} $

已知在某一时刻$ t $目标流体区域欧拉框架下对应的区域为: $ {\cal A}_t(D) = {\cal D} $,且$ J{\rm d}(D) = {\rm d}{\cal D} $,通过代换可得(2.6)式的积分弱形式

(2.10) $ \begin{equation} \left\{ \begin{array}{ll} { } \frac{\rm d}{{\rm d}t} \int_{{\cal D}}(\rho \varphi) {\rm d}{\cal D} = 0, \\ { } \frac{\rm d}{{\rm d}t} \int_{{\cal D}}(\rho {\bf V} \varphi) {\rm d}{\cal D}+\int_{{\cal D}}(\nabla P \varphi) {\rm d}{\cal D} = 0, \\ { } \frac{\rm d}{{\rm d}t} \int_{{\cal D}}(\rho E \varphi) {\rm d}{\cal D}+\int_{{\cal D}}(\nabla \cdot(P{\bf V})\varphi) {\rm d}{\cal D} = 0. \end{array} \right. \end{equation} $

本节将采用间断有限元(DG)方法对弱形式(2.10)进行空间离散化.

首先,对计算空间进行网格划分.初始时刻计算区域为$ D $,将$ D $剖分成包含$ n $个三角网格$ \{D_{i}|i = 1, \cdot\cdot\cdot, n\} $的区域.利用网格速度能够对区域和剖分进行更新,得到新的$ {\cal D} $和剖分$ \{{\cal D}_{i}|i = 1, \cdot\cdot\cdot, n \}$.有限元空间定义为一个分片多项式集合

$ W_h^1 = \{\varphi_h\in L_2({\cal D}):\varphi_h \mid_{{\cal D}_{i}} \in P_1({\cal D}_{i}); 1\leq i\leq n\}, $

其中$ P_1({\cal D}_{i}) $记为$ {\cal D}_{i} $上不超过$ 1 $次的多项式集合.

定义$ \{\lambda_1, \lambda_2, \lambda_3\} (\lambda_1+\lambda_2+\lambda_3 = 1) $作为单元$ {\cal D}_{i} $的面积坐标,并在${\bf \lambda }$$ = (\lambda_1, \lambda_2) $平面内考虑由顶点$ {\bf \lambda }$$ _{\Gamma, 1} = (0, 0) $, $ {\bf \lambda }$$ _{\Gamma, 2} = (1, 0) $, ${\bf \lambda }$$ _{\Gamma, 3} = (0, 1) $确定的标准单元$ \Gamma $.一般单元$ {\cal D}_{i} $和标准单元能够通过如下映射转换: $ ({\cal A}_t)_{i}: (\lambda_1, \lambda_2) \rightarrow (x, y). $在离散的情况下,将每个单元$ {\cal D}_{i} $对应的映射$ ({\cal A}_t)_{i} $组成的集合作为拉格朗日坐标(面积坐标)和欧拉坐标之间的映射$ ({\cal A}_t) $,每个映射$ ({\cal A}_t)_{i} $保持线性和可逆性.我们选取拉格朗日一次节点基函数作为有限元空间的基函数:设$ u_h $为数值解向量$ {\bf U}_h $中任一元素,它能够表达为

$ u_h = \sum \limits_{l = 1}^{3}u_{il}B_{li}(x, y), $

其中$ B_{1i}(x, y), B_{2i}(x, y), B_{3i}(x, y) $是单元$ {\cal D}_{i} $上对应的拉格朗日节点基函数.

使用间断有限元求解(2.10)式的过程如下:寻求一个数值解

$ {\bf U}_h = (\rho_h, (\rho{\bf V})_h, (\rho E)_h)\in W_h^k, $

对所有的检验函数$ \varphi_h\in W_h^k $和单元$ {\cal D}_{i} $满足下列方程组

(3.1) $ \begin{equation} \frac{\rm d}{{\rm d}t} \int_{ {\cal D}_{i} }{\bf U}_h \varphi_{h} {\rm d}{\cal D}_{i}+\int_{{\cal D}_{i}}\nabla \cdot (F({\bf U}_h), G({\bf U}_h)) \varphi_{h} {\rm d}{\cal D}_{i} = 0, \end{equation} $

其中

$ F({\bf U}_h) = \left(\begin{array}{cccc} F_1 \\ F_2 \\ F_3 \\ F_4 \end{array}\right) = \left(\begin{array}{cccc} 0 \\ p_h \\ 0 \\ (pu)_h \end{array}\right), {\quad} G({\bf U}_h) = \left(\begin{array}{cccc} G_1 \\ G_2 \\ G_3 \\ G_4 \end{array}\right) = \left(\begin{array}{cccc} 0 \\ 0 \\ p_h \\ (pv)_h \end{array}\right). $

由格林公式,可得

(3.2) $ \begin{equation} \frac{\rm d}{{\rm d}t} \int_{ {\cal D}_{i} }{\bf U}_h \varphi_{h} {\rm d}{\cal D}_{i} = \int_{{\cal D}_{i}}(F({\bf U}_h), G({\bf U}_h))\cdot\nabla\varphi_{h} {\rm d}{\cal D}_{i} -\int_{\partial{\cal D}_{i}}(F({\bf U}_h), G({\bf U}_h){\bf n})\varphi_{h} {\rm d}\partial{\cal D}_{i}, \end{equation} $

其中$ {\bf n} = (n_1, n_2) $为单元边$ \partial{\cal D}_{i} $的单位外法向量.

由于(3.2)式对所有的检验函数成立,检验函数选取为基函数$ B_{ji}(x, y), j = 1, 2, 3 $可得

(3.3) $ \begin{equation} \frac{\rm d}{{\rm d}t} \sum \limits_{l = 1}^{3}\int_{ {\cal D}_{i} }u_{li}B_{li}B_{ji} {\rm d}{\cal D}_{i} = \int_{{\cal D}_{i}}(F, G)\cdot\nabla B_{ji} {\rm d}{\cal D}_{i} -\int_{\partial{\cal D}_{i}}(\widehat{F({\bf U}_h)}, \widehat{G({\bf U}_h)})\cdot{\bf n}B_{ji}{\rm d}\partial{\cal D}_{i}, \end{equation} $

其中$ (\widehat{F({\bf U}_h)}, \widehat{G({\bf U}_h)})\cdot{\bf n} $为单元边上的数值通量.本文定义一个单元边上的单值函数: $ \widehat{H} = (\widehat{F({\bf U}_h)}, \widehat{G({\bf U}_h)})\cdot{\bf n} $代替数值通量,我们选用文献[2]中的HLLC通量进行数值计算.设

$ (\widehat{m_i})_{lj} = \int_{{\cal D}_{i}} B_{li}(x, y)B_{ji}(x, y){\rm d}{\cal D}_{i}\ (l, j = 1, 2, 3), $

(3.3)式演化为

(3.4) $ \begin{equation} \frac{\rm d}{{\rm d}t} {\bf M_i}\cdot{\bf u}_{hi} = R({\bf U}_{hi}), \end{equation} $

其中$ R({\bf U}_{hi}) $为空间离散算子; $ {\bf u}_{hi} = (u_{hi1}, u_{hi2}, u_{hi3}) $; $ {\bf M_i} = \left((\widehat{m_i})_{lj}\right)_{3\times3} $为质量矩阵.

本节首先给出网格节点移动速度的选取方式.我们对^[4]中的四边形网格速度求解器进行延伸,获得适用于非结构网格的拉格朗日速度.在非结构网格上,假设连接一个顶点$ m $共有$ k_m $条线段$ (k_m = 2, 3, \cdot\cdot\cdot). $以$ k_m = 6 $为例,对这六条线段编号: $ k = 1, \cdot\cdot\cdot, 6. $假设第$ k $条线段的两个端点为$ m $和$ m_k $,规定第$ k $条线段的方向为顶点$ m $指向$ m_k $的方向.我们分别计算左右单元物理量的多项式在$ m $处的值$ (\rho^{k\pm}, (\rho u)^{k\pm}, (\rho v)^{k\pm}, (\rho E)^{k\pm}) $,得到第$ k $条线段左右单元在点$ m $处的速度:

$ (u^{k+}, v^{k+}) = (\frac{(\rho u)^{k+} }{\rho^{k+}}, \frac{(\rho v)^{k+}}{\rho^{k+}}), (u^{k-}, v^{k-}) = ( \frac{(\rho u)^{k-}}{\rho^{k-}}, \frac{(\rho v)^{k-}}{\rho^{k-}}). $

顶点$ m $处的网格速度求解方法如下所示.

步骤1  设$ (n_1^k, n_2^k) $为第$ k $条线段顺时针法向量; $ (t_1^k, t_2^k) = (-n_2^k, n_1^k) $为其切向量;设$ w_t^{k-}, w_t^{k+} $为第$ k $条线段左右两侧的切向速度; $ w_n^{k-}, w_n^{k+} $为其两侧的法向速度.进而,我们计算沿着第$ k $条线段在顶点$ m $处的切向速度

$ w_t^{k} = \frac{1}{2}(w_t^{k-}+w_t^{k+}) , k = 1, 2, \cdots , 6. $

步骤2 在文献[4]中,沿着第$ k $条线段在顶点$ m $处的法向速度$ w_n^{k} $取为线段两侧法向速度的Roe平均.本文改用欧拉框架下HLLC通量(参见文献[5])的熵波速度$ S_M $代替两侧法向速度的Roe平均.相比速度的Roe平均,熵波速度引入的数值粘性更小,更加精确地逼近真实流体在顶点$ m $处的法向速度.我们得到沿着第$ k $条线段在顶点$ m $处的法向速度

$ w_n^{k} = S_M^{k}, k = 1, 2, \cdots , 6, $

其中, $ S_M^{k} $为点$ m $第$ k $条线段上的熵波速度.

步骤3 设$ w_x^{k} $和$ w_y^{k} $为$ m $沿着第$ k $条线段$ x $-$ \!\! $和$ y $-$ \!\! $方向的速度分量.结合第$ k $条线段在$ m $处的法向和切向速度,易得

$ w_x^{k} = (w_n^{k} n_1^k-w_t^{k} n_2^k), w_y^{k} = (w_n^{k} n_2^k+w_t^{k} n_1^k). $

结合所有线段在$ m $的$ x $-$ \!\! $和$ y $-$ \!\! $方向速度分量,可得点$ m $的网格速度

$ u_m^g = \frac{1}{6}(w_x^{1}+w_x^{2}+\cdots +w_x^{6}), v_m^g = \frac{1}{6}(w_y^{1}+w_y^{2}+\cdots +w_y^{6}). $

结合得到的网格移动速度,本文采用二阶Runge-Kutta(RK)方法(参见文献[1])完成(3.3)式在时间方向的离散. DG方法在移动网格上的RK时间离散的具体形式参见文献[8].

对于分片多项式逼近, DG方法将导致在强间断处产生虚假振荡.本文采用一种非结构三角形网格上的HWENO重构(参见文献[3])作为限制器控制虚假振荡.该重构保持了DG方法的紧致性,完美地达到了压制虚假振荡的效果.

我们考虑二维低密度流问题,检测本文格式的精度.初始条件为: $ \rho(x, y, 0) = 1+0.99\cdot \sin(2\pi(x+y)), $$ u(x, y, 0) = 1, v(x, y, 0) = 1, p(x, y, 0) = 1. $本算例的计算区域选取为: $ [0, 1]\times[0, 1]. $算例精确解: $ \rho(x, y, t) = 1+0.99\cdot \sin(2\pi(x+y-2t)), u(x, y, t) = 1, $$ v(x, y, t) = 1, $$ p(x, y, t) = 1. $模拟时间为$ t = 0.1. $表 1给出格式关于密度的误差和精度, $ h({\cal D}(t)) $表示网格中所有单元外接圆的最大直径.从表中可以看出格式能够达到二阶精度.

该激波管问题的初始条件为: $ (\rho_L, u_L, v_L, p_L) = (2.0, 0, 0, 10^9), x \leq 0; (\rho_R, u_R, v_R, p_R) = (0.001, 0, 0, 1), x > 0. $计算区域为$ [-10, 10]\times[0, 0.1] $,比热比系数$ \gamma $为1.4.网格单元边长设为0.005. 图 1给出$ t = 0.0001 $时数值结果,本算例给出$ y = 0.1 $这条线上的数值解.通过比较可以看出本文格式是可行的.

Sedov问题是一个包含由强爆破产生的圆柱对称型强激波传播过程的经典算例,它经常被用于测试拉式格式的鲁棒性.该问题的计算区域设为$ [0, 1.21]\times[0, 1.21] $.初始时刻的网格剖分设为:网格数为3020,节点数为1583,边数为4602.本算例的初始密度为1,初始速度为0,初始时刻原点$ (0, 0) $所在单元比内能为$ 400 $,其余区域为0.边界条件采用固壁条件.该问题在$ t = 1 $时的解析解包含一个以原点为圆心,位于半径1位置的对称圆柱形激波,密度峰值为6. 图 2给出$ t = 1 $时的数值结果.

从图中可以看出格式得到的结果很好地逼近了精确解,同时很好地保持了网格的良好品质和解的对称性.

Salzman问题被广泛应用于拉式算法的健壮性测试中.本算例的初始计算区域为$ \Omega(0) = [0, 1]\times[0, 0.1] $,网格由$ 200\times10 $的三角剖分组成,详情参见文献[7].当$ t = 0 $时盒子内是$ \gamma = \frac{5}{3} $的静止理想气体,气体密度为1,压力为$ 10^{-4} $.问题产生一个向右移动的激波,波后密度为4,水平速度为1.当数值模拟时间达到$ t = 0.6, $激波位置在$ x = 0.8. $图 3给出$ t = 0.6 $时的网格, $ x $-$ \!\! $方向密度与采用文献[4]中的网格节点速度算法得到的结果对比.从图中看出本文结果的激波位置等与解析解吻合,并且本文格式改用熵波速度代替法向速度分量去求解网格速度得到的结果相比文献[4]采用Roe平均得到的结果更加逼近解析解.