一种适用于四面体网格的P1 Galerkin粘性离散化的无参数弱边界条件处理方法

《Computers & Fluids》：A parameter-free weak boundary condition procedure for P 1 Galerkin viscous discretization for tetrahedral meshes

　　科西莫·塔尔西亚·莫里斯科（Cosimo Tarsia Morisco）与西川弘明（Hiroaki Nishikawa）共同研究，隶属于GammaO团队，该团队属于INRIA-ONERA联合项目——位于法国萨克莱（Saclay）的INRIA中心，巴黎萨克莱大学（Université Paris-Saclay），邮编F-91120。

**摘要**
在本文中，我们提出了一种适用于四面体网格的P1 Galerkin粘性离散化的无参数弱边界条件处理方法。边界顶点的残差通过边界通量进行闭合，随后与内部残差一同求解，以获得所有顶点的数值解。边界通量直接从弱形式的方程中推导出来，并以线性精确的方式离散化，以保持二阶精度。这种方法利用在边界顶点指定的边界条件值，通过对相邻四面体元素进行有限元梯度计算来实现。所提出的方法可以将罚项分解为一致的无参数项，并在边界面上循环高效地实现。此外，该方法在同一框架内也实现了诺伊曼（Neumann）边界条件。实验表明，与强边界条件方法相比，弱边界条件不仅能够产生更高质量的解，还能在基于伴随法的网格调整（mesh adaptation）过程中加快收敛速度。

**引言**
P1 Galerkin粘性离散化方法已被广泛应用于基于有限元方法的实际计算流体动力学（CFD）代码中[1]，以及以顶点为中心的有限体积方法[2]，尤其是在四面体网格上。这种方法能够用紧凑的格式实现二阶精度，并且支持高效的基于边的粘性项实现[3][4]。通常，这种离散化会使用强边界条件，尤其是在固体壁面上施加无滑移边界条件时[2]。例如，边界顶点处的x-动量方程离散化形式为Resjx=uj?0，其中uj表示该点的x方向速度。然而，修改后的残差会导致与内部残差不同的离散方程，这往往需要在算法中采用特殊处理（例如分层聚类[5]），并可能在伴随解中引入不连续行为[6]。尽管强边界条件在许多实际应用中得到使用，但证明其求解器的稳定性较为困难（求解器可能收敛缓慢或根本不收敛[7]）。此外，如果边界层没有得到充分解析，强边界条件也可能导致粘性流动问题的解不准确[8][9]。而且，它不适合在分辨率不足的网格上实现壁函数。另一方面，弱边界条件在分辨率不足的网格上可以无缝工作，并且非常适合用于壁函数的实现，但目前关于弱边界条件的研究非常有限，尤其是对于可压缩纳维-斯托克斯方程（Navier–Stokes equations）的P1 Galerkin粘性离散化。

**文献综述**
关于专门为P1 Galerkin粘性离散化开发的弱边界条件，文献中的相关研究较少。Bazilevs和Hughes[8]提出了一种带有罚项的稳定化有限元公式，用于施加狄利克雷（Dirichlet）边界条件，并将其应用于不可压缩纳维-斯托克斯方程[9]的研究中。他们证明了，在矩形域中使用三维笛卡尔网格时，弱边界条件比强边界条件能够更准确地预测中等至高雷诺数下的不可压缩湍流通道流动。对于非结构化网格，Bazilevs和Hughes的弱边界条件方法已被扩展到可压缩纳维-斯托克斯方程的高阶稳定化有限元公式中[6][10]。如参考文献[6][10]所述，弱边界条件能够在泛函（functional）计算中实现超收敛。NASA的FUN3D代码框架中也开发了一种包含壁函数的弱边界条件方法，该代码使用P1 Galerkin粘性离散化处理非结构化网格[11]。还有其他一些相关工作提出了一些用于弱施加边界条件的罚项[12][13]，但这些方法通常涉及可调参数，这些参数往往依赖于具体问题和网格类型。此外，这些参数还可能与元素片的厚度有关，从而严重影响离散化精度[14]。一种无罚项的非对称Nitsche方法[15]是一个例外，该方法不添加包含参数的罚项，而是通过在变分公式中加入额外项来弱化边界条件的施加。这些额外项不涉及参数，但需要在边界面上积分数值解与边界条件值之间的跳跃。这种方法最初是为不连续Galerkin方法提出的[16][17]。

**本文贡献**
在本文中，我们提出了一种从弱形式方程中推导出弱边界条件的简单方法，该方法与用于无粘性项的有限体积类型离散化方法一致。这是其他研究[6][9][10]中未考虑的独特特点。我们的目标离散化方法是一种混合方法：无粘性项通过定义在顶点周围的对偶体积上的基于边的离散化方法进行离散化，而粘性项则采用P1 Galerkin离散化方法。对于基于边的离散化，可以使用保持精度的边界闭合公式轻松实现弱边界条件（二阶公式见参考文献[18]的附录B，更通用的公式适用于二阶和三阶方法，其高效实现方法见参考文献[19][20]）。这种方法会在边界顶点的残差中加入通过边界对偶面的通量贡献，从而闭合无粘性残差。然而，由于无粘性和粘性项的离散化方法不同，因此无法直接使用相同的边界闭合公式来闭合粘性残差。为了闭合粘性残差，我们需要通过边界对偶面添加粘性通量，并且必须以线性精确的方式来实现，以保持边界顶点的二阶精度。

**主要思想**
我们的目标是从弱形式方程中识别边界积分项，并使用已知在边界顶点的边界条件，通过线性精确算法对这些积分项进行离散化。我们将证明每个边界积分项正好对应于通过边界对偶面的粘性通量。从这个意义上说，该方法可以称为无罚项方法。有趣的是，边界项可以重新表述为通过与边界对偶面平行方向的四面体有限元梯度之和，以及在边界面上的三个顶点处定义的罚项之和。因此，所提出的方法自动引入了罚项，而无需其他基于罚项的方法中常见的参数（如罚项系数或分母中的长度）。该方法与参考文献[15]中的无罚项非对称Nitsche方法相关，尽管该方法也看似无罚项，但实际上包含跳跃项，这一点我们稍后会有讨论。本文介绍了新弱边界条件的推导过程及其在基于伴随法的网格调整中的优势，与强边界条件方法进行了对比。

**节选内容**
**纳维-斯托克斯方程**
稳态可压缩纳维-斯托克斯方程的质量、动量和能量守恒方程表述为：
\[ \div(F_{inv}) - F_{vis} = 0, \]
其中，\(F_{inv}=(\rho u)(\rho u \otimes u)+pI(\rho E + p)u\), \(F_{vis}=0\), \(T_Tu - Q\)。\(F_{inv}\)和\(F_{vis}\)分别代表无粘性和粘性通量，\(I\)是3×3的单位矩阵，\(\rho\)表示密度，\(u\)是速度（笛卡尔坐标下为\(u=(u,v,w)\)，\(p\)是静压，\(E\)是总比能量，定义为流动比内能\(e\)和比动能之和，即\(E=e+\|u\|^2\)。

**离散化**
我们的目标离散化方法是无粘性项采用基于边的离散化，粘性项采用P1离散化。以四面体网格为例，首先描述了带有二阶边界闭合的无粘性离散化方法，之后讨论不带边界闭合的粘性离散化方法。

**弱边界条件**
在本节中，我们提出了一种新的弱边界条件处理方法。首先从弱形式方程中推导出边界积分项，然后描述如何基于线性精确算法离散化这些边界积分项。最后，提供了一组可以在所提出的弱边界条件框架内实现的物理边界条件。

**网格调整**
为了研究所提出的弱边界条件对网格调整的影响，我们考虑了一种基于度量的网格调整方法。在这种方法中，可以使用度量场来规定网格的局部尺寸和各向异性[29]。如果一个四面体的所有边在该度量下长度为单位长度，则称其为单位四面体。类似地，如果一个网格在给定度量下满足单位条件，则该网格由（准）单位元素组成。

**数值结果**
本文研究的流动求解器是Wolf[28]，这是一种适用于非结构化网格的基于顶点的混合有限体积-有限元求解器。在Wolf中，系统通过以下无量纲化群进行无量纲化：
\[ \rho_{\infty}, \|u_{\infty}\|, cv,L_{\ref}, \]
其中，\(\rho_{\infty}\)是自由流密度，\(\|u_{\infty}\|是自由流速度，\(cv\)是理想气体的恒容比热容，\(L_{\ref}\)是用户指定的特征长度尺度。

**结论**
本文提出了一种适用于四面体网格的P1 Galerkin粘性离散化的无参数弱边界条件处理方法。边界残差通过边界通量进行闭合，随后与内部残差一起求解。边界通量直接从弱形式方程中推导出来，并以线性精确的方式离散化。边界闭合过程中使用的有限元梯度是根据边界顶点指定的边界条件值计算的。

**作者贡献声明**
科西莫·塔尔西亚·莫里斯科（Cosimo Tarsia Morisco）负责撰写原始草稿、可视化处理、验证、软件开发、方法论研究、形式分析、数据整理及概念化；西川弘明（Hiroaki Nishikawa）负责撰写原始草稿、验证、方法论研究、形式分析、数据整理及概念化。

**利益冲突声明**
作者声明以下可能被视为潜在利益冲突的财务利益/个人关系：西川弘明指出其研究得到了美国陆军研究办公室（US Army Research Office）的财政支持。如果还有其他作者参与，则他们声明没有已知的可能影响本文工作的财务利益或个人关系。

**致谢**
作者感谢Frederic Alauzet（INRIA, GammaO团队）和Enza Parente（Safran Tech）就本研究提供的有益讨论和建设性意见。Adrien Loseille提供了Feflo.a软件。第二作者感谢美国陆军研究办公室在合同/授权号W911NF-19-1-0429项下的支持。