在应用研究方面：Borazjani等[2-3]模拟了血管动脉流动和人工心脏瓣膜流动；Zhu和Peskin[4]模拟了流动肥皂泡薄膜中飘动的柔性细丝，并与实验结果[5]进行了对比；Mittal等[6]将尖锐接触曲线(sharp interface curvib approach)和浸入边界法结合，成功模拟了雷诺数Re=6 300的三维鱼胸鳍绕流，预测了绕流的主要性质。
在算法改进方面：Goldstein等[7]在1993年结合谱方法提出了虚拟边界法(Ghost Boundary Method，GBM)，通过引入反馈力来考虑固体边界对流场的影响，从而使得物面无滑移边界条件得到满足；Mittal等提出了尖锐接触面(sharp interface)方法，并先后应用于复杂构型的不可压缩流[8]、可压缩流[9]的模拟中；基于有限元方法[10-11]的浸入边界法几乎同时被Boffi和Gastaldi[12]、Wang和Liu[13]、Zhang等[14]提出，其中，Boffi和Gastaldi[12]提出的方法用变分公式代替了传统的Dirac-delta函数，Heltai等[15-16]通过流固耦合标准算例验证了该方法的计算精度。
浸入类方法多采用笛卡儿网格，导致网格尺度足够小才能在边界上达到要求的计算精度，但全局均匀加密网格会增加计算量，降低计算效率。因此，国内外研究者多在浸入类方法中引入网格自适应技术解决这一问题。网格自适应技术可分为网格节点移动自适应(r-refinement)、网格局部升阶自适应(p-refinement)、网格局部加密自适应(h-refinement)。考虑到网格节点移动自适应不便处理几何大变形，网格局部升阶自适应无法提高边界几何描述精度，本文采用网格局部加密自适应技术，即通过在局部区域加入新的网格单元，细化局部网格，提高网格的分辨率以更好地捕捉流场特征和固体边界几何特征。如吴泽艳等[17]将网格局部加密自适应应用于简单WENO限制器(Weighted Essentially Non-Oscillatory limiter)-间断Galerkin方法，并实现了对欧拉方程的求解。基于树状结构的分层笛卡儿网格自适应技术最先由Berger和Oliger[18]提出并应用于求解一维和二维的双曲型偏微分方程。Abgrall等[19]改进了一种带自适应结构网格的Level Sets浸入边界法，并利用各向异性的网格自适应弥补了边界模拟的精度。Peng和Zhou[20]基于当地DFD(Domain-Free Discretization)的浸入边界法研究了一种基于流场特征的网格自适应技术，以涡量作为自适应指示因子，但未考虑几何层面的自适应。Xu等[21]基于浸入边界-格子玻尔兹曼方法(Immersed Boundary-Lattice Boltzmann Method, IB-LBM)提出了新的几何自适应方法。该方法针对初始网格为多层不同尺度的笛卡儿网格，以边界位置作为自适应指示因子，在边界运动过程中通过加入/删除新的网格点完成自适应过程，但未考虑流场特征自适应。
对于固体浸没于流体并自由漂动的问题，流固边界不易形成强剪切层，单纯从流场特征或几何特征出发难以保证对此类问题的模拟精度。为了改善这一不足，本文在Boffi和Gastaldi[12]提出的浸入有限元方法基础之上，提出了一种耦合流场特征和几何特征的网格自适应策略，并应用于方腔顶盖驱动流中黏弹性圆盘运动问题。
1 基本算法 浸入有限元方法采用有限元方法求解弱形式控制方程，其积分运算无需直接求解传统的浸入边界法中的Dirac-delta函数，避免了因Dirac-delta函数的分布影响界面插值的精度的问题。计算域中流体区域采用非贴体笛卡儿网格离散，固体区域使用贴体网格提供边界几何信息。
动边界流固耦合问题的控制方程由质量方程、动量方程以及固体的运动方程组成。如图 1所示，整个物理区域为固定于空间的区域Ω，其中固体区域为Bt，流体域为Ω\Bt。固体的运动可由映射ζ:B→Bt, x=ζ(s, t)描述，其中B为初始时刻的固体区域，s∈B, x∈Ω, 时间t∈[0, T]，T为最终时间。
 |
| 图 1 计算域的定义 Fig. 1 Definition of calculation domain |
| 图选项 |
针对不包含边界的区域控制方程可写为
 | (1) |
 | (2) |
 | (3) |
式(1)为质量方程，式(2)为动量方程, 式(3)为固体边界上满足的连续性条件。ρ为质量密度，u为速度矢量，σ为柯西应力张量，b为单位质量的外力密度，n为固体边界的单位法向矢量。对于外边界?Ω，需满足特定的Dirichlet边界条件和Neumann边界条件。
根据网格定义域不同，密度ρ可分为流体域和固体域两部分：
 | (4) |
式中：ρf为流体的密度；ρs为固体的密度。对于不可压缩牛顿流体，ρf为常数。
流体的柯西应力可表达为
 | (5) |
式中：pf为压力；I为单位张量；μf为流体的动力黏性系数；L=gradu。
对于不可压缩超弹性结构体，柯西应力可表示为
 | (6) |
式中：ps为满足不可压缩约定的拉格朗日乘子；σsv和σse分别为柯西应力的黏性部分和弹性部分。
各部分的柯西应力表达式为
 | (7) |
式中：μs为固体的动力黏性系数；F为变形梯度张量，可由当前时刻物质点的当前位置ζ(s, t)定义，即F=?ζ(s, t)/?s；J=det F；Pse为第一Piola-Kirchhoff应力张量。
本文采用INH (Incompressible Neo Hookean)模型作为弹性模型。INH模型的第一Piola-Kirchhoff应力张量可定义为
 | (8) |
式中：Ge为弹性剪切模量。
固体的运动通过位移场描述，位移场定义为
 | (9) |
对式(9)求全导数，可得
 | (10) |
式(1)、式(2)和式(10)组成了本文的控制方程。
因此，该方法中需求解的未知数为u(x, t), p(x, t), w(s, t)(x∈Ω, s∈B, t∈[0, T])。
令V0、Q和Z分别为速度、压力以及位移对应的测试函数空间，并取任意的v∈V0、q∈Q以及z∈Z分别对控制方程式(1)、式(2)和式(10)应用变分原理, 并代入柯西应力的表达式可得其变分方程[15]为
 | (11) |
 | (12) |
 | (13) |
式中：dV和dV′分别为区域Ω和B上的微元；τg为分部积分后的边界项；da为边界?Ω上的微元；J(s, t)为衡量几何体积变化的雅可比矩阵；ρs0的定义为
 | (14) |
在本文的求解器中，对于有限元的选取，离散空间的速度uh和位移wh选取二次连续Lagrange有限元，对于压力ph选取一次间断Legendre有限元。
2 网格自适应技术 网格自适应可以分为3个步骤：①标记网格的加密和粗化；②执行网格自适应，更新网格；③将流场解映射到自适应之后的新网格。文本仅研究步骤①的算法的策略，不讨论步骤②、③的细节。
2.1 基于几何特征的自适应 为了更好地模拟固体区域，需根据固体位置进行几何自适应。本文的几何自适应是对每一时间步固体所在的笛卡儿网格区域进行全部或部分加密。因此，构造合适的几何加密区域是几何自适应的核心。
由第1节介绍的控制方程(式(1)、式(2)和式(10))可知，位移场可在求解完控制方程后得到，因此，针对任一时刻t的固体区域Bt上网格点的坐标ζ(s, t)计算式为
 | (15) |
如图 2所示，Bt为t时刻的固体区域且有Bt的质心为Gt，取gr为几何加密的比例系数且满足0≤gr≤1，将Bt区域以Gt为中心缩放为原来的(1-gr)倍得到Ω2区域，则加密区域为Bt\Ω2区域(图中的Ω1区域)所在的笛卡儿网格区域。当gr=1时，固体所对应的笛卡儿网格区域全部被加密。
 |
| 图 2 几何自适应示意图 Fig. 2 Schematic of adaptive refinement based on geometry feature |
| 图选项 |
2.2 基于流场特征的自适应 基于流场特征的自适应技术通常是在流场状态量变化剧烈的区域对网格局部加密，从而减小流场解的误差。该方法的核心是构造合适的自适应指示因子。对于不可压缩流动来说，压强和速度是一组可以完备描述流场状态的量。本文选取由速度梯度导出的涡量作为网格自适应的指示因子。设任一时间步j的离散区域为Ωh，对于Ωh内的任一网格单元K有自适应指示因子wj(K)，即
 | (16) |
在确定了自适应指示因子之后，还需确定自适应门槛值以保证加密或粗化后网格量可以稳定在预期值。本文根据需加密/粗化区域指示因子的和占指示因子总和的比例确定自适应的门槛值，具体采用的方法如下。
记E为wj(K)的和，即
 | (17) |
则网格自适应范围由式(18)决定：
 | (18a) |
 | (18b) |
式中：ar为加密比例常数，ac为粗化比例常数，两者均介于0和1之间且满足ar+ac < 1；Mr和Mc分别表示将要被加密和粗化的单元的集合，其生成算法步骤如下：
1) 计算每个网格单元K上的指示因子并求和。
2) 对Ωh中所有的自适应指示因子排序。
3) 由空集开始按照指示因子从大到小的顺序将对应的网格单元加入加密集合，直到式(18a)不满足时对应的整个加密集合即为Mr。
4) 由空集开始按照指示因子从小到大的顺序将对应的网格单元加入粗化集合，直到式(18b)不满足时对应的整个粗化集合即为Mc。
相比于常规的自适应门槛值的设置(当网格单元指示因子大于某一特定值时加密，小于某一特定值时粗化)，本文在合适的加密/粗化比例常数下，能保证自适应过程中网格量的波动较小。
2.3 网格悬点的处理 网格局部加密会导致悬点问题，如图 3所示，C2单元自适应加密一次后，其他单元不受影响，但V3、V4、V5和V6成为了悬点。悬点的处理方式与边界条件类似，同样可分为强形式和弱形式，二者的区别在于悬点所在界面(如图 3中的f1界面)的连续性约束强弱不同。相比于弱形式，强形式的约束可在总体矩阵中删除悬点自由度所在的行和列，减小问题的求解规模，提高计算效率。因此，本文采用强形式的约束处理悬点。
 |
| 图 3 网格悬点形成 Fig. 3 Hanging point formation of mesh |
| 图选项 |
当网格自适应出现一条边上多个悬点时，会降低网格的插值精度。为避免多悬点问题的出现，本文的求解器对网格层级较小一侧的网格进行加密。
2.4 两种自适应指标的综合 当几何自适应和流场自适应的标识符(加密标识符、粗化标识符和无标识符)同时存在于同一单元时，需制定合理的标识符设置规则确定该单元最终的标识符。本文采用的标识符设置规则为加密标识符的优先级最高，粗化标识符次之，无标识符的优先级最低。根据这一规则，可得如图 4所示的自适应流程。
 |
| 图 4 网格自适应流程图 Fig. 4 Flowchart of mesh adaptation |
| 图选项 |
3 算例验证 本文以方腔顶盖驱动圆盘运动问题为例，对比研究了均匀加密、流体或几何单指标自适应以及本文耦合自适应策略的计算精度。计算问题的定义参考文献[16]，如图 5所示，方腔为边长为L的正方形，其中的流体和固体均为不可压缩介质且密度相同，故计算中忽略体积力。流体为黏性体，固体为黏弹性体。本算例中设置μf=μs=0.01 Pa·s，ρf=ρs=1.0 kg/m3，圆盘的弹性剪切模量Ge=1.0 Pa。初始时刻半径为R=0.2m，圆心为C的圆盘在方腔顶盖速度为U=1.0 m/s的作用下运动，其中，C的坐标为(0.6, 0.5) m。圆盘与流体的密度相同，因此在运动过程中两者之间无强剪切作用。
 |
| 图 5 方腔顶盖驱动圆盘问题说明 Fig. 5 Illustration of disk entrained in a lid-driven cavity flow |
| 图选项 |
根据方腔区域网格的均匀加密次数以及设置自适应指标的不同共划分为6种网格形式计算。表 1中给出了不同网格形式的对比。计算基准网格是由单位正方形均匀二分5次得到的32×32均匀笛卡儿网格，记为Case1。在基准网格的基础上分别均匀加密1次和2次得到的网格记为Case2和Case3。自适应算例由基准网格开始计算，在各自自适应策略的驱动下对网格进行加密：Case4为流场单指标自适应，Case5为几何单指标自适应，Case6为设置耦合自适应的情况。不同初始加密层级对应的有限元自由度对比(设置自适应的情况为初始时刻的自由度)如表 2所示，固体单元始终保持相同的自由度，而控制体区域的网格随着加密层级的不同而不同。固体区域的网格示意图如图 6所示。
表 1 6种不同形式的网格设置 Table 1 Mesh setting for six cases
| Case | 初始加密层级 | 几何自适应 | 流场自适应 | 流场加密系数 | 流场粗化系数 | 几何加密系数 |
| Case1 | 5 | 否 | 否 | — | — | — |
| Case2 | 6 | 否 | 否 | — | — | — |
| Case3 | 7 | 否 | 否 | — | — | — |
| Case4 | 5 | 否 | 是 | 0.85 | 0.1 | — |
| Case5 | 5 | 是 | 否 | — | — | 1.0 |
| Case6 | 5 | 是 | 是 | 0.6 | 0.35 | 0.4 |
表选项
表 2 初始时刻网格自由度 Table 2 Degrees of freedom of mesh at initial time
| 初始加密层级 | 控制体 | 固体 | 总自由度 | |||
| 网格量 | 自由度 | 网格量 | 自由度 | |||
| 5(Case1, Case4~Case6) | 1 024 | 11 522 | 320 | 2 626 | 14 148 | |
| 6(Case2) | 4 096 | 45 570 | 320 | 2 626 | 48 196 | |
| 7(Case3) | 16 384 | 181 250 | 320 | 2 626 | 183 876 | |
表选项
 |
| 图 6 固体区域网格 Fig. 6 Mesh of solid area |
| 图选项 |
为便于对比各策略的计算精度，通过调整自适应参数使各自适应算例的自由度规模保持一致。如表 1所示，Case4中流场自适应的加密比例系数ar=0.85，粗化比例系数ac=0.1；Case5中对圆盘所在的笛卡儿网格区域全部加密，即设置gr=1.0；Case6中流场自适应的加密比例系数ar=0.6，粗化比例常数ac=0.35，且设置几何部分加密，加密比例系数gr=0.4。图 7展示了不同情况下的计算自由度DoF随时间的变化。由图中可以看出，通过设置不同的自适应比例系数，自适应过程中Case4、Case5和Case6的计算自由度相当。
 |
| 图 7 自由度随时间变化 Fig. 7 Variation of degrees of freedom with time |
| 图选项 |
计算精度通过圆盘区域速度散度2-范数以及圆盘的运动轨迹2个指标衡量。
由于本文求解的为不可压缩固体，因此，理想情况下圆盘区域每个离散单元的体积与总体积均保持不变。因此，本文选取圆盘区域速度散度的2-范数作为误差判断的标准之一。
除了考虑圆盘运动过程中的体积守恒，圆盘的运动轨迹也是人们关心的一个数值模拟结果。由于整个圆盘的运动轨迹难以量化，因此，本文对圆盘上4个点(初始位置为图 6中的P1、P2、P3、P4)的运动轨迹进行分析。另外，由于本算例中圆盘的运动轨迹无理论解，此处以最密网格Case3的运动轨迹作为基准，对比其他情况下的运动轨迹与Case3的误差。
图 8展示了6种情况下圆盘区域速度散度的2-范数div v_norm随时间的变化。表 3展示了6种不同情况下速度散度2-范数的时间平均值。结合图 8和表 3可知，当不设置自适应时，随着网格加密层级的增加，圆盘区域速度散度的2-范数逐渐减小；当设置自适应时，在3种不同自适应的计算自由度相当的情况下，Case5和Case6的圆盘区域速度散度2-范数与Case3相当，然而Case4与Case1相比，仅在6 s附近自适应对结果有所改善，此时，Case4的平均值在10-2量级，而Case5和Case6均在10-3量级，与Case3的平均值差值仅为10-4量级。因此，几何单指标自适应和耦合自适应的计算误差相当，而流场单指标自适应的误差远大于其他2种自适应的情况。这是由于以圆盘区域速度散度2-范数作为误差评判指标时，仅考虑了圆盘各个离散单元在流场中的相对位置，因此基于几何特征的自适应占主导。
 |
| 图 8 圆盘区域速度散度2-范数随时间的变化 Fig. 8 Variation of 2-norm of divergence of velocity in disk area with time |
| 图选项 |
表 3 速度散度2-范数平均值 Table 3 Average value of 2-norm of divergence of velocity
| Case | 速度散度2-范数平均值 |
| Case1 | 0.022 733 |
| Case2 | 0.009 509 |
| Case3 | 0.002 129 |
| Case4 | 0.017 686 |
| Case5 | 0.002 195 |
| Case6 | 0.002 242 |
表选项
图 9~图 12分别为4个点轨迹x坐标和y坐标的误差(xerror为x坐标误差，yerror为y坐标误差)。4个点的轨迹误差均说明3种自适应情况在计算自由度相当的情况下，Case6的误差最小，Case4的误差其次，Case5的误差最大。为了更直观地反映出不同自适应情况下的误差大小，表 4展示了4个点的x坐标和y坐标轨迹误差2-范数的大小(xnorm和ynorm分别为x坐标、y坐标的轨迹误差2-范数)。从表中可以看出，Case4、Case5、Case6的轨迹误差2-范数分别为10-2、10-1、10-3量级。因此，耦合自适应的轨迹计算误差最小，流场单指标自适应的次之，几何单指标自适应的误差大于其他两种自适应的情况。造成几何单指标自适应方法误差大的原因是以圆盘上某点轨迹作为评判指标时考虑了圆盘上的点在流场中的绝对位置，这与流场本身的模拟精度有关，此时基于流场特征的自适应占主导。
 |
| 图 9 P1点轨迹x坐标和y坐标误差随时间的变化 Fig. 9 Variation of error of x and y coordination for point P1 trajectory with time |
| 图选项 |
 |
| 图 10 P2点轨迹x坐标和y坐标误差随时间的变化 Fig. 10 Variation of error of x and y coordination for point P2 trajectory with time |
| 图选项 |
 |
| 图 11 P3点轨迹x坐标和y坐标误差随时间的变化 Fig. 11 Variation of error of x and y coordination for point P3 trajectory with time |
| 图选项 |
 |
| 图 12 P4点轨迹x坐标和y坐标误差随时间的变化 Fig. 12 Variation of error of x and y coordination for point P4 trajectory with time |
| 图选项 |
表 4 四个点的轨迹误差2-范数 Table 4 2-norm of trajectory error for four points
| Case | P1 | P2 | P3 | P4 | |||||||
| xnorm | ynorm | xnorm | ynorm | xnorm | ynorm | xnorm | ynorm | ||||
| Case1 | 0.381 017 | 0.309 065 | 0.175 893 | 0.169 636 | 0.293 602 | 0.247 673 | 0.130 823 | 0.079 848 | |||
| Case2 | 0.135 22 | 0.104 793 | 0.067 734 | 0.067 368 | 0.084 943 | 0.068 042 | 0.042 763 | 0.019 129 | |||
| Case4 | 0.044 911 | 0.033 318 | 0.022 543 | 0.016 775 | 0.041 846 | 0.024 332 | 0.029 854 | 0.025 508 | |||
| Case5 | 0.370 969 | 0.301 401 | 0.197 824 | 0.171 009 | 0.308 501 | 0.266 88 | 0.126 822 | 0.081 096 | |||
| Case6 | 0.004 382 | 0.001 385 | 0.001 874 | 0.002 688 | 0.001 667 | 0.001 309 | 0.002 096 | 0.000 906 | |||
表选项
图 13~图 15分别给出了3种自适应策略下t=3，4，5，6 s这4个时刻的速度云图(图中v为速度矢量的模)以及网格自适应状态示意图。其中，图 13为仅设置流场自适应的情况；图 14为仅设置几何自适应的情况；图 15为设置耦合自适应的情况。从图 13~图 15中可以看到，弹性体在流体驱动下运动以及在流体剪切驱动下变形的现象。顶盖附近是流体的强剪切区，弹性体形变剧烈。从图 13中可以看出，当仅设置流场自适应时，圆盘仅在流体速度场作用下运动，圆盘与周围流体不存在明显的剪切作用，因此，自适应方法无法感知弹性体的存在，加密区域仅为靠近方腔上壁面涡量较大的区域。当仅设置几何自适应时，如图 14所示，自适应加密区域仅为圆盘周围区域，但是流动的模拟精度难以保证，依赖于流场的弹性体运动的模拟精度也就会受到不利影响。当设置流场和几何耦合自适应时，如图 15所示，自适应加密的区域同时涉及了靠近方腔顶盖附近涡量较大的区域和圆盘周围的区域。
 |
| 图 13 仅设置流场自适应情况下3~6 s的速度云图 Fig. 13 Velocity contour at t=(3-6) s when setting mesh adaptation based flow features only |
| 图选项 |
 |
| 图 14 仅设置几何自适应情况下3~6 s的速度云图 Fig. 14 Velocity contour at t=(3-6) s when setting mesh adaptation based geometry features only |
| 图选项 |
 |
| 图 15 设置耦合自适应情况下3~6 s的速度云图 Fig. 15 Velocity contour at t=(3-6) s when setting mesh adaptation based on flow and geometry features |
| 图选项 |
综上所述，当设置耦合自适应后，可同时加密涡量较大的区域和圆盘周围区域，因此相比于流场或几何单指标自适应，耦合自适应对圆盘区域速度散度2-范数和圆盘的运动轨迹误差的模拟精度更高。
4 结论 网格自适应的效果取决于自适应指标和策略的设置。针对浸入液体中黏弹性运动体的模拟，本文基于浸入有限元方法开展了相应的自适应方法研究，具体结论如下：
1) 构造了一种耦合流场涡量和运动固体几何位置的自适应策略，并在流场自适应中采用指示因子总和的比例作为门槛值，在几何自适应中根据网格点到质心的距离场全部/部分加密。
2) 通过方腔驱动圆盘运动算例验证了耦合自适应方法的优势。耦合自适应中，流场自适应保证圆盘运动轨迹的计算精度，几何自适应保证圆盘体积的计算精度。因此，当分别以圆盘区域速度散度2-范数和圆盘上特征点的运动轨迹误差作为评判标准，在流场自适应、几何自适应和耦合自适应3种情况的计算自由度相当时，耦合自适应的结果优于其他2种自适应的情况。
目前，基于浸入有限元的自适应方法研究还在进一步发展之中，尤其针对复杂几何外形的流固耦合问题。因此，本文下一步工作将围绕复杂几何外形的流固耦合问题开展，研究高效、鲁棒的适用于复杂几何外形的网格自适应方法。
参考文献
| [1] | PESKIN C. Flow patterns around heart valves:A numerical method[J]. Journal of Computational Physics, 1972, 10(2): 252-271. |
| [2] | BORAZJANI I, GE L, SOTIROPOULOS F. High-resolution fluid-structure interaction simulations of flow through a bi-leaflet mechanical heart valve in an anatomic aorta[J]. Annual Review of Biomedical Engineering, 2010, 38: 326-344. DOI:10.1007/s10439-009-9807-x |
| [3] | SOTIROPOULOS F, BORAZJANI I. A review of state-of-the-art numerical methods for simulating flow through mechanical heart valves[J]. Medical & Biological Engineering & Computing, 2009, 47(3): 245-256. |
| [4] | ZHU L, PESKIN C S. Simulation of a flapping flexible filament in a flowing soap film by the immersed boundary method[J]. Journal of Computational Physics, 2002, 179(2): 452-468. |
| [5] | ZHANG J, CHILDRESS S, LIBCHABER A, et al. Flexible filaments in a flowing soap film as a model for one-dimensional flags in a two-dimensional wind[J]. Nature, 2000, 408(6814): 835-839. DOI:10.1038/35048530 |
| [6] | MITTAL R, SESHADRI V, UDAYKUMAR H S. Flutter, tumble and vortex induced autorotation[J]. Theoretical and Computational Fluid Dynamics, 2004, 17(3): 165-170. |
| [7] | GOLDSTEIN D, HANDLER R, SIROVICH L. Modeling a no-slip flow boundary with an external force field[J]. Journal of Computational Physics, 1993, 105(2): 354-366. |
| [8] | GHIAS R, MITTAL R, DONG H. A sharp interface immersed boundary method for compressible viscous flows[J]. Journal of Computational Physics, 2007, 225(1): 528-533. |
| [9] | MITTAL R, DONG H. A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[J]. Journal of Computational Physics, 2008, 227(10): 4825-4852. DOI:10.1016/j.jcp.2008.01.028 |
| [10] | 秦望龙, 吕宏强, 伍贻兆. 基于混合网格的高阶间断有限元黏流数值解法[J]. 力学学报, 2013, 45(6): 987-991. QIN W L, LV H Q, WU Y Z. High-order discontinuous Galerkin solution of N-S equations on hybrid mesh[J]. Chinese Journal of Theoretical and Applied Mechanics, 2013, 45(6): 987-991. (in Chinese) |
| [11] | 于剑, 阎超. Navier-Stokes方程间断Galerkin有限元方法研究[J]. 力学学报, 2010, 42(5): 962-970. YU J, YAN C. Study on discontinuous Galerkin method for Navier-Stokes equations[J]. Chinese Journal of Theoretical and Applied Mechanics, 2010, 42(5): 962-970. (in Chinese) |
| [12] | BOFFI D, GASTALDI L. A finite element approach for the immersed boundary method[J]. Computers and Structures, 2003, 81(8-11): 491-501. DOI:10.1016/S0045-7949(02)00404-2 |
| [13] | WANG X, LIU W K. Extended immersed boundary method using FEM and RKPM[J]. Computer Methods in Applied Mechanics and Engineering, 2004, 193(12-14): 1305-1321. DOI:10.1016/j.cma.2003.12.024 |
| [14] | ZHANG L, GERSTENBERGER A, WANG X, et al. Immersed finite element method[J]. Computer Methods in Applied Mechanics and Engineering, 2004, 193(21-22): 2051-2067. DOI:10.1016/j.cma.2003.12.044 |
| [15] | HELTAI L, FRANCESCO C. Variational implementation of immersed finite element methods[J]. Computer Methods in Applied Mechanics and Engineering, 2012, 229-232: 110-127. DOI:10.1016/j.cma.2012.04.001 |
| [16] | ROY S, HELTAI L, COSTANZO F. Benchmarking the immersed finite element method for fluid-structure interaction problems[J]. Computers and Mathematics with Applications, 2015, 69(10): 1167-1188. DOI:10.1016/j.camwa.2015.03.012 |
| [17] | 吴泽艳, 王立峰, 武哲. 基于简单WENO-间断Galerkin的Euler方程自适应计算[J]. 北京航空航天大学学报, 2016, 42(4): 806-814. WU Z Y, WANG L F, WU Z. Adaptive simple WENO limiter-discontinuous Galerkin method for Euler equations[J]. Journal of Beijing University of Aeronautics and Astronautics, 2016, 42(4): 806-814. DOI:10.13700/j.bh.1001-5965.2015.0237 (in Chinese) |
| [18] | BERGER M J, OLIGER J. Adaptive mesh refinement for hyperbolic partial differential equations[J]. Journal of Computational Physics, 1984, 53(3): 484-512. |
| [19] | ABGRALL R, BEAUGENDRE H, DOBRZYNSKI C. An immersed boundary method using unstructured anisotropic mesh adaptation combined with level-sets and penalization techniques[J]. Journal of Computational Physics, 2014, 257: 83-101. DOI:10.1016/j.jcp.2013.08.052 |
| [20] | PENG B, ZHOU C H. An approach of dynamic mesh adaptation for simulating 3-dimensional unsteady moving-immersed-boundary flows[J]. International Journal for Numerical Methods in Fluids, 2018, 87(4): 180-201. DOI:10.1002/fld.4486 |
| [21] | XU L C, TIAN F B, YOUNG J, et al. A novel geometry-adaptive Cartesian grid based immersed boundary-lattice Boltzmann method for fluid-structure interactions at moderate and high Reynolds numbers[J]. Journal of Computational Physics, 2018, 375: 22-56. DOI:10.1016/j.jcp.2018.08.024 |
