基于高阶间断伽辽金数值方法的浸没边界法.pdf

文章编号：0258-1825(2023)09-0096-11基于高阶间断伽辽金数值方法的浸没边界法王婷婷，吕宏强*，安慰(南京航空航天大学航空学院，南京210016)摘要：为了降低复杂边界条件下贴体网格生成难度，近年来基于笛卡尔网格的浸没边界法逐渐成为研究此类问题的主要数值方法之一，然而计算精度和计算效率仍然是此类方法目前面临的挑战。与传统的空间 2 阶精度的有限体积格式相比，空间精度为 3 阶或 3 阶以上的高精度方法具有空间精度高、数值分辨率高、数值耗散小的优点，而作为高精度数值方法之一的间断伽辽金方法在浸没边界方面的应用仍较少。本文将高阶间断伽辽金方法的高精度优势与浸没边界法无需贴体网格的优势结合起来，提出了适用于可压缩流动的高精度浸没边界法。其中边界条件采用体积惩罚方式实现，同时采用牛顿法迭代以及 MPI 并行提高计算效率，物面的数据重建采用插值点处的逆距离权重（IDW-IP）方法代替高阶格式下多项式插值方法。本文基于笛卡尔网格测试验证了二维定常和非定常情况下的数值模拟效果，并与传统贴体网格的计算结果进行了对比。关键词：浸没边界法；间断伽辽金；体积惩罚；界面函数；数据重建中图分类号：O357.1；V211.3文献标识码：Adoi：10.7638/kqdlxxb-2022.0138Immersed boundary method based on high-orderdiscontinuous Galerkin methodWANGTingting，LYUHongqiang*，ANWei(College of Aerospace Engineering,Nanjing University of Aeronautics and Astronautics,Nanjing210016,China)Abstract:To decrease the difficulty of generating body-fitted meshes for complex geometries,theimmersedboundary(IB)methodbasedonCartesianmeshhasgraduallybecomeapopularnumericalmethodforstudying such problems in recent years.However,accuracy and efficiency of this method still need to beimproved.Comparedwiththetraditionalfinitevolumeschemewithsecond-orderspatialaccuracy,third-orderorhigher-ordernumericalmethodshavetheadvantagesofhighspatialaccuracy,highnumericalresolutionandlownumericaldissipation.However,thediscontinuousGalerkin(DG)method,asoneofthehigh-precisionnumericalmethods,isstillrarelyappliedwiththeimmersedboundary.Thispaperproposesahigh-orderdiscontinuousGalerkinmethodforcompressibleflowsbycombiningtheadvantageofthehigh-orderdiscontinuousGalerkinmethodandtheimmersedboundarymethod.Theboundaryconditionsinthispaperarerealizedbythevolumepenalizationmethod.NewtonsmethoditerationandMPIareusedtoimprovecomputationalefficiency.Theinversedistanceweightatinterpolationpoint(IDW-IP)insteadofthehigh-orderpolynomialdefinedineachelementtoperformhigh-orderdatareconstructionattheobjectsurface.BasedontheCartesiangrids,numericaltestsareperformedfortwo-dimensionalsteadyandunsteadyflows,andcomparisonswiththoseontraditionalbody-fittedgridsaregiven.Keywords:immersedboundarymethod；discontinuousGalerkin；volumepenalization；maskfunction；datareconstruction收稿日期：2022-09-07；修订日期：2022-10-09；录用日期：2022-10-24；网络出版时间：2022-11-22基金项目：气动噪声控制重点实验室基金（ANCL20190103）作者简介：王婷婷（1997），女，浙江温州人，硕士研究生，研究方向：计算流体力学.E-mail：通信作者：吕宏强*，教授，研究方向：计算流体力学.E-mail：引用格式:王婷婷,吕宏强,安慰.基于高阶间断伽辽金数值方法的浸没边界法J.空气动力学学报,2023,41(9):96106.WANGTT,LYUHQ,ANW.Immersedboundarymethodbasedonhigh-orderdiscontinuousGalerkinmethodJ.ActaAerodynamicaSinica,2023,41(9):96106(inChinese).doi:10.7638/kqdlxxb-2022.0138第41卷第9期空气动力学学报Vol.41，No.92023年9月ACTA AERODYNAMICA SINICASept.，2023 0 引言浸没边界法（immersedboundarymethod,IBM）最初由 Peskin 于 1972 年1提出，用于心脏瓣膜周围流动的模拟。在该方法中1-2，流体在笛卡尔网格上求解，几何体被一组拉格朗日点所取代，通过在控制方程中添加源项实现边界对流体的作用，如 Lee3和Nicolaou4所述，这类 IBM 称为连续强迫 IBM。自从IBM 提出以来，该方法得到了巨大的发展并广泛应用。Mohd-Yusof5和 Fadlun 等6提 出 了 直 接 强 迫IBM，与连续强迫 IBM 不同的是，不使用拉格朗日点并且强迫项在笛卡尔网格上直接进行计算，控制方程仅仅在流体域的计算网格上求解。何跃龙等7引入Choi 等8提出的指数插值法对目标点速度进行插值，建立了适用于可压缩流动的浸没边界法。Yang等9采用大涡模拟与浸没边界相结合的方法，对后向台阶流动进行了数值模拟。Giannenas 和 Laizet10提出了一种基于三次样条重建的浸没边界方法，用于基于笛卡尔网格的湍流中浸没物体的高保真模拟。Kou等11提出一种基于体积惩罚的结合高阶通量重建方法的浸没边界法。尽管 IBM 已得到广泛应用，但采用笛卡尔网格结合高阶方法的探索仍较少。近年来，间断伽辽金（discontinuousGalerkin,DG）方法已成为 CFD 中最有前途的高阶离散技术之一12-13。1997 年，Bassi 和 Rebay 等14-15采用 DG 方法求解了欧拉方程和 N-S 方程，随后 Bassi 等16又将 DG 方法拓展到可压缩流和不可压缩流的 DNS（direct numerical simulation）和 ILES（implicit largeeddysimulation）。Landmann17在其博士论文中发展了并行高阶 DG 方法并对二维 N-S 方程及 RANS（Reynolds-averagedNavier-Stokes）方程进行了求解。DG 方法在国内起步较晚，蔚喜军和周铁18在 2005 年采用龙格-库塔间断伽辽金（Runge-KuttadiscontinuousGalerkin,RKDG）方法求解了二维可压缩欧拉方程。于剑和阎超19通过引入全局提升算子和局部提升算子，发展了求解 N-S 方程的 DG 方法的一般框架。张来平等20发展了 DG/FV 混合方法，计算效率相较于同阶精度的 DG 方法更高。党亚斌等21发展了一种基于解析精确 Jacobian 矩阵的 GMRES隐式方法，用于求解可压缩层流及湍流问题，提高了隐式高精度DG 方法的计算稳定性和计算效率。国外已有研究者将 DG 方法与浸没边界法相结合，但大多采用切割网格的方法。针对基于笛卡尔网格的浸没式几何体，Qin 等22提出了一种求解欧拉方程的 DG 方法。Mller等23提出了一种高阶离散格式，用于求解具有浸没边界的可压缩欧拉和 N-S 方程，在一个由水平集函数隐式定义的域中使用 DG 离散化，采用了一种非侵入式单元聚集技术。虽然切割网格法对于静态边界能够较好的实现，然而切割网格法仅限于低精度，而且小切割单元的处理以及动态边界的处理比较复杂。本文引入 Kou 等11提出的基于体积惩罚的浸没边界实现思路，提出基于惩罚项模式的高阶 DG 浸没边界法，采用典型算例对所提出方法进行测试验证。1 控制方程可压缩黏性流体的控制方程17写成：Ut+Fi(U)=Fv(U,U)(1)UFiFv式中：为守恒变量；为无黏通量；为黏性通量。2 间断伽辽金空间离散因为式（1）中存在 2 阶偏导数项，在求解前必须对方程进行降阶处理。采用 Landmann17使用的混合方法，在计算中把 1 阶偏导数项当作额外变量，引入辅助变量。式（1）可以改写成如下形式：U=0(2)Ut+Fi(U)Fv(U,)=0(3)W采用间断伽辽金有限元法对式（2）和式（3）进行空间离散。在式（2）和式（3）的两边同乘一个任意的测试函数，并将测试函数的高阶表达式代入，在计算域单元边界上的积分中引入逆风格式，得到如下方程组：wejde=wejUh(Uh,U+h)ndwejUhde(4)twejUhde+wejFi(Uh,U+h)ndwejFi(Uh)de=wejFv(Uh,U+h,+)ndwejFv(Uh,)de(5)Uh(Uh,U+h)nUhU+hUh(Uh,U+h)n式（4）中为辅助数值通量，其中和分别表示交界面左单元和右单元的变量值。本文式（4）中辅助数值通量和式（5）中的黏第9期王婷婷等：基于高阶间断伽辽金数值方法的浸没边界法97Fv(Uh,U+h,+)n性数值通量均采用 BR2 格式24。3 隐式时间离散对式（4）和式（5）进行数值积分之后，得到的方程组如下所示：Mdudt=R(u)(6)u=u1,u2,uk,uNeleTNeleukR(u)=R1,R2,Rk,RNeleTRkM式中：为全局自由度矢量，为 全 局 单 元 的 总 数，为 单 元 k 的 自 由 度；为全局残值矢量，为单元 k 的残值矢量；为全局块质量矩阵。tn+1采用隐式时间离散方法，式（6）右端残值项取新的时间步，计算得到：Mdudt=R(un+1)(7)采用向后差分方法处理之后，式（7）可以写成：Mun+1unt=R(un+1)(8)由于求解的是非线性的控制方程，所以对式（8）的右端项进行线性化：R(uk+1)=R(uk)+(Ru)k(uk+1uk)(9)(R/u)ktn+1式中，为需要求解的大型线性方程组的雅克比矩阵。在每个时间步内采用牛顿迭代法进行求解，对于定常算例，只需进行一个牛顿步就可以达到需要的精度：u(0)=u(n)(10)Mt(Ru)k(uk+1uk)=R(uk)Mt(uku0)uk+1=uk+c(uk+1uk)(11)un+1=u(kmax)(12)k=(0,1,kmaxc 0,1tnkmaxtn+1其中：)；为阻尼因子。在一个时间步内即一个牛顿迭代循环内，以第 步的数值解作为初始值进行牛顿迭代，直到式（11）中的残值下降到设定阈值，将第个牛顿步的结果作为新的时间步的数值解并进入下一个时间步进行计算。4 浸没边界法IBM 的基本思想是通过适当的数值处理向非贴体网格施加边界条件。体积惩罚是一种典型的方法，它通过惩罚固体中积分点的速度来施加边界条件。4.1 体积惩罚法体积惩罚法通过在控制方程中引入惩罚源项来施加边界条件。在这种方法中，首先定义一个用于区分流体区域和固体区域的界面函数：(x,t)=1,if x S0,otherwise(13)(x,t)该界面函数用于确定是否对当前积分点施加 IBM力。对于移动的边界，是随时间而变化的。对于高阶笛卡尔网格的体积惩罚方法（图 1），所有被固体区域覆盖的积分点都需要进行惩罚以施加边界条件。固体流体流体点=0固体点=1图 1 高阶方法的体积惩罚示意图Fig.1 Volume penalization for high-order method用 IBM 处理的 N-S 方程可以写成：Ut=R(U)+S(U)(14)S式中，“”为 IBM 源项。us=(us,vs)T考虑固体内部速度狄利克雷边界条件，源项为：S(U)=0usuvsv2(u2s+v2s)2(u2+v2)(15)us=(0,0)T式中，为体积惩罚参数。考虑无滑移边界条件时，设条件。本文惩罚项遵循文献 25 中的公式，以施加绝热壁的边界条件。上述方程也可用于运动物体，其中固体速度根据运动方程进行更新。通过结合算子分裂和隐式时间离散方法，控制方程等式（14）在时间上进行。使用 2 阶 Strang 分裂26方法添加源项，正如 Piquet 等27所讨论的，通过 Strang分裂可以精确计算动量和能量方程的惩罚项。在时间步骤 n，执行以下操作：1）U1Unt1=S(U1),t1=t/2,U1,0=Un(16)2）U2U1t2=R(U2),t2=t,U2,0=U1(17)98空气动力学学报第41卷3）Un+1U2t3=S(Un+1),t3=t/2,Un+10=U2(18)Mt(Ru)kt在数值计算中，将控制方程的变量进行无量纲化处理17。下文中如无特别说明，惩罚参数均为=1103。在时间步骤 n，执行第 2 步，即对式（17）进行牛顿法迭代时，式（11）中的无量纲时间步长=0.1。4.2 边界表示和界面函数固体边界的表示是 IBM 方法的关键之一，其主要用来进行界面函数的定义和气动系数的计算。=0 浸没边界的离散化应足够灵活，以处理复杂的几何形状。在本文中，选择使用一组拉格朗日标记点来表示实体边界，定义为浸没边界点。标记点由线性元素（即二维线段）连接。使用该表示法可以有效地计算几何量，包括曲面法线、用于数据重建的插值模板、曲面距离28-29。浸没边界由离散的点定义，点在网格中的位置用水平集方法可描述为满足，是到浸没边界的有符号距离30。对于圆柱体、球体或机翼等简单几何体，存在解析形状函数，可用于定义界面函数。例如，对于直径为 1的圆，有：(x,t)=1,if x2+y2 0.520,otherwise(19)4.3 曲面数据重建为了从模拟结果中获得相关物理量（如压力、摩擦力）的分布，固体表面物理量的重建是一个重要工作。在高阶框架中，可以直接考虑使用每个单元中定义的高阶多项式来执行数据重建。然而，这种方法有一些潜在的缺点，可能导致结果不准确。这主要是因为在大多数情况下，此类插值方法将涉及浸入固体中的积分点，而这些积分点通常具有非物理值。除多项式插值外，另一种方法是从几个最近的积分点插值数据。与基于高阶多项式的插值相比，其具有更大的灵活性来选择插值点（interpolationpoint,IP），而不涉及具有非物理值的固体点。这里，采用文献 31中描述的适用于高阶框架的插值方法。通常，文献 31中的插值方法基于 IB 点与插值点之间的逆距离。首先从 IB 点周围最近的流体点中选择候选插值点，然后插值保守变量的值和梯度。本文采用插值点处的逆距离权重（inversedistanceweightatinterpolationpoint,IDW-IP）方法进行插值，如图 2 所示。A-模板点B-带点C-IB点IP-插值点n-IB点法向场单元带单元内部单元nABCrIPd1,id2,iV图 2 IDW-IP 方法的插值模板示意图Fig.2 Interpolation stencil for IDW-IP method插值点在这种情况下被定义为一个点，位于IB 点的法线上，在这个点上，属性从周围的单元内插，然后转移到 IB 点。关于插值点，首先要确定它的位置，其次是它的属性。为了确定插值点的位置，使用如下关系式：dIP=ni=1d2,i/d1,ini=11/d1,i(20)d1d2式中：为 A 点到 IB 点法线的垂直距离；为 A 点到IB 点的法线的距离投影；i 为模板中场单元和带单元的固体外积分点。为了确定插值点的压力值，使用以下关系式：PIP=ni=1Pi/d1,ini=11/d1,i(21)在这种方法中，在插值点确定的压力被复制到 IB 点。5 算例测试在本节中，应用不同测试算例验证所提出的基于高阶间断伽辽金数值方法的浸没边界法。5.1 定常圆柱绕流us=(0,0)Tx 5D,10D y 5D,5Dx D,D y D,Dh=0.03D第一个测试案例为二维定常圆柱绕流，雷诺数和马赫数分别设置为 40 和 0.2。壁面采用无滑移绝热壁面边界条件，因此在固体中施加。矩形计算域的大小为、，其中D 是 圆 柱 的 直 径，设 置 为 1。在 正 方 形 区 域、，使用均匀网格。均匀网格尺寸 h 分别选择为 0.03D（网格 1）、0.05D（网格 2）和0.1D（网格 3）。网格 1 如图 3 所示，局部均匀网格大小。在网格 1 情况下，由 0 阶至 3 阶 DG 计算得到的流线如图 4、物面附近马赫数云图如图 5所示。从图第9期王婷婷等：基于高阶间断伽辽金数值方法的浸没边界法99Cp中可以看出，随着阶数的提高，同一网格下计算结果的精度显著提高，尤其是物面附近的数值解更加光滑。在网格 1 情况下，不同阶数计算得到的压力系数与 文 献 11 中 贴 体 网 格 模 拟 结 果 的 对 比 如图 6 所示，图中实线表示文献 11 贴体网格的计算结果，符号表示本文方法数值模拟得到的结果（不同符号代表不同的阶数）。结果表明，压力系数分布随着计算阶数的提高与贴体网格结果具有良好的一致性。CpCp不同网格尺度下压力系数的对比如图 7 所示，图中分别给出了 0 阶至 3 阶 DG、物面附近网格尺度为 0.03D、0.05D、0.1D 时的结果对比。由图 7（a-d）所示，在同 1 阶数下，随着网格尺度的减小曲线更光滑，且总体曲线更接近文献 11 贴体网格的测试结果；无论网格 2 还是网格 3，再次验证了同一网格下，计算阶数越高时精度越高，模拟计算得到的结果与文献结果吻合越好。tCp设置无量纲时间步长=1103。在网格 1 情况下，1 阶 DG 至 3 阶 DG 计算得到的压力系数与文献 11 中贴体网格模拟结果的对比见图 8。结果表明，本文方法计算所得压力系数结果与文献结果具有良好的一致性，且阶数越高与文献结果的吻合度越高。由图 6 和图 8 可以看出，同为网格 1 且惩罚参数42y/Dy/D00 x/Dx/D51024332211002211335(a)全局网格(b)局部网格图 3 计算网格（网格 1）Fig.3 Computational grid(grid 1)42024y/Dx/D246Ma020.220.200.180.160.140.120.100.080.060.040.0200.220.200.180.160.140.120.100.080.060.040.0200.220.200.180.160.140.120.100.080.060.040.020Ma0.220.200.180.160.140.120.100.080.060.040.02042024y/Dx/D2460242024y/Dx/D246Ma0242024y/Dx/D246Ma02(a)0阶DG(b)1阶DG(c)2阶DG(d)3阶DG图 4 网格 1 下不同阶数流线图和马赫数云图Fig.4 Streamlines and Mach number contours calculated by DG with different orders(grid 1)100空气动力学学报第41卷相同时，时间步长的取值会对计算结果产生影响。与无量纲时间步长为 0.1 时相比，当无量纲时间步长设置为与惩罚参数相同时（即 1103），1 阶 DG 至 3 阶DG 得到的计算结果均更光滑且与文献结果吻合更好。ttt=t本文方法惩罚项的添加基于 Brinkman 惩罚法25，该方法中将固体视为渗透率极低的多孔介质。Engels等32研究发现，当使用显式时间离散时，惩罚参数的选择受 CFL 的限制，显式时间步长必须小于惩罚参数；当采用隐式时间离散时，允许惩罚参数 小于时间步长。Labert33对 N-S 方程进行了研究，揭示了惩罚参数与时间步长之间的密切耦合误差在惩罚参数 约等于时间步长时达到饱和。因此在体积惩罚法中通常建议，一种解释为：惩罚项作为速度上的强阻尼，其引入了一个特征时间尺度，该特征时间尺度阶数为，必须用时间推进来解决。1.01.00.50.5000.50.51.01.0y/Dx/D1.01.00.50.5000.50.51.01.0y/Dx/DMa0.220.200.180.160.140.120.100.080.060.040.0200.220.200.180.160.140.120.100.080.060.040.0200.220.200.180.160.140.120.100.080.060.040.0200.220.200.180.160.140.120.100.080.060.040.020Ma1.01.00.50.5000.50.51.01.0y/Dx/DMa1.01.00.50.5000.50.51.01.0y/Dx/DMa(a)0阶DG(b)1阶DG(c)2阶DG(d)3阶DG图 5 网格 1 情况下不同阶数物面附近马赫云图Fig.5 Mach number contours calculated by DG with different orders(grid 1)2.01.51.00.50Cp0.51.01.52.0050100/()1500.03D p30.03D p20.03D pl0.03D p0Jiaqing KOU et al11图 6 无量纲时间步长 0.1 时网格 1 情况下不同阶数压力系数t=0.1Fig.6 Pressure distribution calculated by DGwith different orders(grid 1,)第9期王婷婷等：基于高阶间断伽辽金数值方法的浸没边界法101 5.2 非定常圆柱绕流Ma=0.2 Re=200 x 5D,10D y 5D,5Dx D,D y D,D为了验证对非定常解的有效性，测试二维圆柱周围的低雷诺数涡脱落。当、时，计算域大小为、，其中 D 是圆柱的直径，设置为1。在正方形区域、，使用均匀网格，均匀网格尺寸 h 选为 0.05D（图 9）。StSt=0.193St=0.20St=0.194St=0.194对不同阶数下的圆柱绕流进行数值模拟，1 阶DG 和 3 阶 DG 的马赫数云图分别如图 10 和图 11 所示。不论 1 阶 DG 还是 3 阶 DG 结果，从全流场的图上均可以看到圆柱尾迹中有周期性的旋涡脱落。从物面附近流场图可以看出，1 阶 DG 的物面附近存在阶梯状，3 阶 DG 的物面附近较 1 阶 DG 的更为光滑。计算描述涡脱落的特征参数，并与文献 34-35对比。1 阶 DG 的，3 阶 DG 的，与文献 34 中以及文献 35 中较为吻2.01.51.00.50Cp0.51.01.52.0050100/()1500.03D p00.05D p00.1D p02.01.51.00.50Cp0.51.01.52.0050100/()1500.03D p10.05D p10.1D p12.01.51.00.50Cp0.51.01.52.0050100/()1500.03D p20.05D p20.1D p22.01.51.00.50Cp0.51.01.52.0050100/()1500.03D p30.05D p30.1D p3(a)0阶DG(b)1阶DG(c)2阶DG(d)3阶DGJiaqing KOU et al11Jiaqing KOU et al11Jiaqing KOU et al11Jiaqing KOU et al11t=0.1图 7 不同网格尺度下压力系数对比图（）t=0.1Fig.7 Pressure distributions calculated by different grid resolutions()2.01.51.00.50Cp0.51.01.52.0050100/()1500.03D p30.03D p20.03D plJiaqing KOU et al11图 8 无量纲时间步长 0.001 时网格 1 下不同阶数压力系数t=0.001Fig.8 Pressure distributions calculated by DGwith different orders(grid 1,)102空气动力学学报第41卷合。本文所采用的体积惩罚法对边界条件的实现是间接的，采用的是一种弱约束方式，该方法对边界条件的实现精度与贴体网格相比较低，边界附近的计算精度依赖于网格密度、数值方法的精度以及惩罚参数的设置32-33。y=0图 12 为中轴线上时均流向速度分布。由图可看出，1 阶 DG 和 3 阶 DG 的模拟结果中，圆柱内部的速度为 0，同时沿来流方向在圆柱后方近壁面流速为负值，存在先减后增的趋势，故圆柱后方存在回流区，在远处流速趋于稳定。整个回流区及近尾区的流42y/D00 x/D510245图 9 0.05D非定常计算网格Fig.9 Computational grid for unsteady flow simulation(h=0.05D)42005x/Dx/D10245y/DMa0.280.260.240.220.200.180.160.140.120.100.080.060.040.020Ma0.280.260.240.220.200.180.160.140.120.100.080.060.040.0201.51.52.02.53.01.01.00.50.5000.50.51.51.01.0y/D(a)全流场(b)物面附近流场图 10 Re=200 时圆柱马赫数云图（p1）Fig.10 Mach number contours at Re=200(p1)42005x/Dx/D10245y/DMa0.280.260.240.220.200.180.160.140.120.100.080.060.040.020Ma0.280.260.240.220.200.180.160.140.120.100.080.060.040.0201.51.52.02.53.01.01.00.50.5000.50.51.51.01.0y/D(a)全流场(b)物面附近流场图 11 Re=200 时圆柱马赫数云图（p3）Fig.11 Mach number contours at Re=200(p3)1.00.80.60.4u/u00.200123x/D4560.2p1p3R Pasquetti et al34y=0图 12 沿中轴线的时均流向速度分布y=0Fig.12 Mean streamwise velocity in the direction第9期王婷婷等：基于高阶间断伽辽金数值方法的浸没边界法103速分布与文献 34 结果吻合良好，但在恢复远区，可能由于远处网格尺度较大导致本文计算的结果偏小。总体上，3 阶 DG 的结果相较于 1 阶 DG 更接近文献 34 的结果。x/D=1.0y/D=0图 13 至图 15 分别给出了圆柱后的三个不同剖面、2.0、4.5 上的时均流向速度分布。图中实线为文献 34 计算结果，实心三角和空心方形分别表示 1 阶 DG 和 3 阶 DG 的计算结果。可以看出，1 阶DG 和 3 阶 DG 的计算结果与文献 34 结果在不同剖面处流向速度随坐标的变化规律一致，均关于圆柱中心线（）对称，且在两侧趋于稳定数值。通过三个剖面处的结果对比，可以看出，同一网格下 3 阶DG 的计算结果相较于 1 阶 DG 的计算结果与文献更为吻合。总体来说，本文方法成功模拟了低雷诺数下非定常问题，且阶数越高模拟的结果精度越高。1.51.00.50u/u00213123y/Dx/D=4.5p1p3R Pasquetti et al34x/D=4.5图 15 截面的时均流速分布x/D=4.5Fig.15 Mean streamwise velocity at 6 结论本文提出了适用于可压缩流动数值模拟的结合高阶 DG 的浸没边界法，用体积惩罚方法来实现物面边界条件，避免了切割网格。相较于显式迭代，对惩罚项采用隐式时间离散允许较小的惩罚参数。采用IDW-IP 方法对物面附近的数据进行重建，以获得物面表面压力。研究得出以下结论：1）通过不同的流动模拟，包括静态圆柱定常和非定常的流动，验证了本文方法计算可压缩流动的有效性和准确性。2）高阶 DG 的高精度优势在本文的浸没边界方法中得到了显著体现。随着计算阶数的提高，边界条件的实现精度以及整体的数值精度都得到了显著提高。3）与贴体网格方法相比，IBM 最大的优势在于网格生成十分简单，缺点在于边界条件的实现是采用弱约束方式，因此边界处理的精度不如贴体网格方法，计算结果依赖于网格密度、高阶 DG 的阶数以及惩罚参数 的设置。4）本文方法为未来复杂外形及移动边界等要求高质量贴体网格的计算提供了更为简便的思路。但对移动边界的处理还需进一步研究验证。参 考 文 献：PESKINCS.Flowpatternsaroundheartvalves:anumericalmethodJ.JournalofComputationalPhysics,1972,10(2):252271.doi:10.1016/0021-9991(72)90065-41PESKINCS.NumericalanalysisofbloodflowintheheartJ.JournalofComputationalPhysics,1977,25(3):220252.doi:10.1016/0021-9991(77)90100-02LEEI,CHOIH.Adiscrete-forcingimmersedboundarymethodforthe31.51.00.500.5u/u00213123y/Dx/D=1.0p1p3R Pasquetti et al34x/D=1.0图 13 截面的时均流速分布x/D=1.0Fig.13 Mean streamwise velocity at 1.51.00.50u/u00213123y/Dx/D=2.0p1p3R Pasquetti et al34x/D=2.0图 14 截面的时均流速分布x/D=2.0Fig.14 Mean streamwise velocity at 104空气动力学学报第41卷fluid-structure interaction of an elastic slender bodyJ.Journal ofComputationalPhysics,2015,280:529546.doi:10.1016/j.jcp.2014.09.028NICOLAOU L,JUNG S Y,ZAKI T A.A robust direct-forcingimmersed boundary method with enhanced stability for moving bodyproblemsincurvilinearcoordinatesJ.Computers&Fluids,2015,119:101114.doi:10.1016/pfluid.2015.06.0304MOHD-YUSOFJ.Combinedimmersed-boundary/B-splinemethodsforsimulations of flow in complex geometriesR.NASA ARS/StanfordUniversityCTR,Stanford.5FADLUNEA,VERZICCOR,ORLANDIP,etal.Combinedimmersed-boundaryfinite-differencemethodsforthree-dimensionalcomplexflowsimulationsJ.JournalofComputationalPhysics,2000,161(1):3560.doi:10.1006/jcph.2000.64846何跃龙,李盾,刘帅,等.基于指数插值的浸没边界法在可压缩流模拟中的应用研究J.空气动力学学报,2016,34(4):426432.HEYL,LID,LIUS,etal.Theapplicationofpower-lawinterpolationbasedimmersedboundarymethodincompressibleflowsimulationJ.ActaAerodynamicaSinica,2016,34(4):426432(inChinese).7CHOIJ,OBEROIRC,EDWARDSJR,etal.AnimmersedboundarymethodforcomplexincompressibleflowsJ.JournalofComputationalPhysics,2007,224(2):757784.doi:10.1016/j.jcp.2006.10.0328YANGDD,HESD,SHENL,etal.LargeeddysimulationcoupledwithimmersedboundarymethodforturbulentflowsoverabackwardfacingstepJ.ProceedingsoftheInstitutionofMechanicalEngineers,Part C:Journal of Mechanical Engineering Science,2021,235(15):095440622095489.doi:10.1177/09544062209548929GIANNENAS A E,LAIZET S.A simple and scalable immersedboundary method for high-fidelity simulations of fixed and movingobjectsonaCartesianmeshJ.AppliedMathematicalModelling,2021,99(2):606627.doi:10.1016/j.apm.2021.06.02610KOU JIAQING,JOSHI S,HURTADO-DE-MENDOZA A,et al.Immersedboundarymethodforhigh-orderfluxreconstructionbasedonvolumepenalizationJ.JournalofComputationalPhysics,2022,448:110721.doi:10.1016/j.jcp.