随着计算机技术的不断发展和先进数值算法的提出，数值模拟已经成为求解复杂流固耦合（fluid-structure interactions，FSI）问题的有效工具. Peskin^[1]基于非贴体网格提出浸没边界法（immersed-boundary method，IBM）用于模拟血液流经心脏瓣膜问题，核心思想是在Navior-Stokes方程中引入力源项来强化无滑移速度边界条件的施加，通过Delta函数将流固耦合边界力分配到周围流体节点. IBM避免了贴体网格算法中的网格更新操作，为有效处理复杂流固耦合问题指引新的方向，并促进浸没类方法的发展，如改进的浸没边界法^[2-3]、虚拟区域法^[4-7]、浸没有限元法^[8-9]、浸没边界-格子玻尔兹曼法^[10-11]和笛卡尔/浸没边界混合算法^[12-13]等.

Hou等^[14]根据流固耦合力作用范围将浸没类方法分为2类，浸没边界法和浸没域法. 浸没边界法将流固耦合力分配到边界周围流体节点，即在固体域内的部分流体节点也受到耦合力的作用，被称为虚拟流体. 浸没域法在浸没边界法的基础上，将虚拟流体扩展到整个固体域内，以模拟占有实际体积的浸没固体，如虚拟区域法和浸没有限元法. 虚拟区域法是Glowinski等^[4-7]提出的，将包含固体的随时间变化的流体转化为不随时间变化的拓展流体域，应用分布式拉格朗日乘子法在虚拟区域上施加刚体约束，最终得到刚体粒子在流场中的运动过程. 浸没有限元法（immersed finite element method，IFEM）最早由Zhang等^[8-9]提出，将流固耦合力施加到整个虚拟流体域内，分布的体力使虚拟流体保持与固体运动一致的约束条件. 在IFEM基础上，通过改进固体求解器和流固耦合力计算方法，Zhang等^[15-16]提出浸没光滑有限元法，Wang等^[17-19]进一步提出浸没光滑点插值法，能够有效模拟大变形流固耦合问题.

有限元法（finite element method，FEM）常用于求解固体运动或变形，采用低阶单元如简单的三角形或四面体单元会使模型刚度过硬、精度较低，但这类单元在离散复杂几何域时可以自动剖分. Liu等^[20]建立双重弱形式理论，基于梯度光滑技术软化模型刚度，提出光滑点插值法，提高三角形或四面体单元的计算性能，在固体力学分析中得到广泛应用. 光滑有限元方法在所构建的光滑域内采用线性插值方法，而光滑点插值方法弱化场函数的连续性，可以采用高阶插值格式. 边基光滑点插值法（edge-based smoothed point interpolation method，ES-PIM）和点基光滑点插值法（node-based smoothed point interpolation method，NS-PIM）是光滑点插值方法中较为重要的2种方法. ES-PIM能够保证空间和时间上的稳定性，在相同的网格下，相比于有限元法和点基光滑点插值方法，能够产生更加准确的结果，并且具备超收敛和超精度特性^[21]. 点基光滑点插值方法软化模型程度较高，抵抗网格畸变能力较强，能够提供弹性问题能量解的上限^[22-23]. 鲁欢等^[24-25]将模型刚度偏软的点基光滑点插值法和偏硬的有限元方法相结合形成点基局部光滑点插值法^[13]（node-based partly smoothed point interpolation method，NPS-PIM），通过构造适宜的固体模型刚度提高固体求解的准确性. Zhang等^[19]将NPS-PIM与浸没方法相结合求解大变形流固耦合问题.

本研究在浸没方法框架下，分别以FEM、ES-PIM和NPS-PIM作为固体求解器，比较它们在求解流固耦合问题时的计算性能. 不可压缩黏性流体采用半隐式特征线分裂（characteristic-based split，CBS）算法求解，CBS算法基于有限元离散Navior-Stokes方程，能够避免对流项及压力项的数值振荡问题，广泛用于各种复杂流动问题中.

在浸没边界法中，在Navior-Stokes方程中引入体力项以考虑浸没边界对流场的影响：

(1) $\begin{array}{l} {{\rm{\rho }}^{\rm{f}}}({{\partial {{{v}}^{\rm{f}}}} / {\partial {\rm{t}}}}) + {{\rm{\rho }}^{\rm{f}}}({{{v}}^{\rm{f}}} \cdot \nabla ) \cdot {{{v}}^{\rm{f}}} = - \nabla {{{p}}^{\rm{f}}} + {\rm{\mu }}{\nabla ^2}{{{v}}^{\rm{f}}} + {{\rm{\rho }}^{\rm{f}}}{{g}} + {{{f}}^{{\rm{f}},{\rm{FSI}}}}. \end{array} $

式中： ${{\rm{\rho }}^{\rm{f}}}$、 ${{{v}}^{\rm{f}}}$、 ${\rm{\mu}} $、 ${{{p}}^{\rm{f}}}$、 ${{{f}}^{{\rm{f}},{\rm{FSI}}}}$分别为流体的密度、速度、动力黏性系数、压力和流固耦合力， ${{g}}$为重力加速度.

Zhang等^[8]在提出的IFEM中将虚拟流体引入整个固体域，并且假设虚拟流体与真实流体具有相同的物理特性，同时与固体的运动特性保持一致. 将流固耦合力作用于整个虚拟域内，将流固耦合力定义为

(2) $\begin{array}{l} {{{f}}^{{\rm{f}},{\rm{FSI}}}}{\rm{ = }}{{\rm{\rho }}^{\rm{f}}}({{\partial {{{v}}^{\rm{f}}}} / {\partial {\rm{t}}}}) + {{\rm{\rho }}^{\rm{f}}}({{{v}}^{\rm{f}}} \cdot \nabla ) \cdot {{{v}}^{\rm{f}}} +\nabla {{{p}}^{\rm{f}}} - {\rm{\mu }}{\nabla ^2}{{{v}}^{\rm{f}}} - {{\rm{\rho }}^{\rm{f}}}{{g}}. \end{array}$

式（2）适用范围为被固体覆盖的虚拟流体域 ${{{\varOmega }}^{{\rm{fc}}}}$，并且满足条件： ${{{\varOmega }}^{{\rm{s}}}}={{{\varOmega }}^{{\rm{fc}}}}$、 ${{{\varGamma }}^{{\rm{fc}}}} = {{{\varGamma }}^{\rm{s}}}$.其中， ${{{\varOmega }}^{{\rm{s}}}}$为固体域， ${{{\varGamma }}^{{\rm{fc}}}}$、 ${{{\varGamma }}^{{\rm{s}}}}$分别为虚拟流体边界和固体边界.

因此，真实流体和虚拟流体分布于整个计算域，可以基于固定欧拉网格描述整个流体域，采用拉格朗日网格描述固体. 整个流固耦合系统控制方程可以分解为3个部分：

1）不可压缩黏性流体的控制方程：

(3) $\left. \begin{array}{l} {{\rm{\rho}} ^{\rm{f}}}({{\partial {{{v}}^{\rm{f}}}} / {\partial {\rm{t}}}}) + {{\rm{\rho}} ^{\rm{f}}}({{{v}}^{\rm{f}}}\cdot \nabla ) \cdot {{{v}}^{\rm{f}}}{\rm{ = }} - \nabla {{{p}}^{\rm{f}}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{\mu}} {\nabla ^2}{{{v}}^{\rm{f}}} + {{\rm{\rho}} ^{\rm{f}}}{{g}} ,\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\nabla \cdot {{{v}}^{\rm{f}}} = 0 .\\ \end{array} \right\}$

2）固体运动方程：

(4) ${{\rm{\rho}} ^{\rm{s}}}{\ddot {{u}}^{\rm{s}}} = {{\partial {{{P}}^{\rm{s}}}} / {\partial {{{X}}^{\rm{s}}}}} + {{\rm{\rho}} ^{\rm{s}}}{{g}}.$

式中， ${{\rm{\rho}} ^{\rm{s}}}$、 ${{{u}}^{\rm{s}}}$、 ${{{P}}^{\rm{s}}}$分别为固体的密度、位移、第一类皮奥拉-基尔霍夫应力.

3）流固耦合条件：

流固耦合系统的速度条件和力条件分别为

(5) ${{{v}}^{\rm{f}}} = {{{v}}^{\rm{s}}},$

(6) ${{f}}_{}^{{\rm{s}},{\rm{FSI}}} = - {{f}}_{}^{{\rm{f}},{\rm{FSI}}}.$

式中： ${{v}}^{\rm{s}} $、 ${{f}}_{}^{{\rm{s}},{\rm{FSI}}} $分别为固体的速度、流固耦合力. 式（5）、（6）适用于虚拟流体域.

本研究采用CBS算法，通过有限元方法对流体域离散得到求解不可压缩黏性流体控制方程：

(7) $ \begin{split} \!\!M_{IJ}^{\rm{f}}{(^*}v_{Ji}^{\rm{f}}{ - ^n}v_{Ji}^{\rm{f}})/\Delta t = & -{ ^n} C{{_{IJ}^{\rm{f}}}^n}v_{Ji}^{\rm{f}}{ - ^n}F_{Ii}^{\rm{f}} - (\Delta t/2){^n}K{{_{IJ}^{\rm{f}}}^n}v_{Ji}^{\rm{f}} + \\ &{{\;^n}F_{Ii}^{{\rm{f}},t}{{\rm{ + }}^n}F_{Ii}^{{\rm{f}},{{g}}}}{{{\rm{ = }}^*}{\rm{RHS}}_{Ii}^{\rm{f}}{ + ^n}F_{Ii}^{{\rm{f}},{{g}}}} \end{split} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!,$

(8) $^{n}H_{IJ}^{\rm{f}}{}^{n + 1}p_J^{\rm{f}}{\rm{ = }} - {{{\rm{(1}}}/ {\Delta t}})Q_{IJi}^{\rm{f}}{}^*v_{Ji}^{\rm{f}}\;,$

(9) $ \begin{split} M_{IJ}^{\rm{f}}{{({}^{n + 1}v_{Ji}^{\rm{f}}{ - ^n}v_{Ji}^{\rm{f}})} / {\Delta t}}& = M_{IJ}^{\rm{f}}{{({}^*v_{Ji}^{\rm{f}}{ - ^n}v_{Ji}^{\rm{f}})} / {\Delta t}} - G_{IJi}^{\rm{f}}{}^np_J^{\rm{f}}= \\ & {}^{n + 1}{\rm{RHS}}_{Ii}^{\rm{f}} \;. \end{split} $

式中： $\Delta t$为时间步长， $^nv_{Ji}^{\rm{f}}$、 ${}^np_J^{\rm{f}}$分别为单元第J个流体节点在第n步时的速度和压力， ${}^*v_{Ji}^{\rm{f}}$为单元第J个流体节点中间的修正速度， ${}^{n + 1}v_{Ji}^{\rm{f}}$、 ${}^{n + 1}p_J^{\rm{f}}$分别为单元第J个流体节点在第n+1步时的速度和压力， ${}^*{\rm{RHS}}_{Ii}^{\rm{f}}$、 ${}^{n + 1}{\rm{RHS}}_{Ii}^{\rm{f}}$为简化后的替代等式左侧的变量，其他变量的具体计算形式如下：

$ \begin{split} &M_{IJ}^{\rm{f}} = \int\limits_{{\varOmega ^{\rm{f}}}} {{{\rm{\rho}} ^{\rm{f}}}\varPhi _I^{\rm{f}}\varPhi _J^{\rm{f}}{\rm{d}}\varOmega } ,\;{}^nF_{Ii}^{\rm{f}}{\rm{ = }}\int\limits_{{\varOmega ^{\rm{f}}}} {{{(\partial \varPhi _I^{\rm{f}}} / {\partial x_j^{\rm{f}}}})\tau _{ij}^{\rm{f}}{\rm{d}}\varOmega }, \\ &{}^nC_{IJ}^{\rm{f}}{\rm{ = }}\int\limits_{{\varOmega ^{\rm{f}}}} {{{\rm{\rho}} ^{\rm{f}}}{{\varPhi}} _I^{\rm{f}}({{\partial ({}^nv_j^{\rm{f}}{{\varPhi}} _J^{\rm{f}})} / {\partial x_j^{\rm{f}}}}){\rm{d}}\varOmega }, \\ &{}^nH_{IJ}^{\rm{f}}{\rm{ = }}\int\limits_{{\varOmega ^{\rm{f}}}} {{{(\partial {{\varPhi}} _I^{\rm{f}}} / {\partial x_j^{\rm{f}}}})({{\partial {{\varPhi}} _J^{\rm{f}}} / {\partial x_i^{\rm{f}}}}){\rm{d}}\varOmega }, \\ &Q_{IJ}^{\rm{f}}{\rm{ = }}\int\limits_{{\varOmega ^{\rm{f}}}} {{{\rm{\rho}} ^{\rm{f}}}{{\varPhi}} _I^{\rm{f}}{{(\partial {{\varPhi}} _J^{\rm{f}}} / {\partial x_i^{\rm{f}})}}{\rm{d}}\varOmega } ,G_{IJi}^{\rm{f}}{\rm{ = }}\int\limits_{{\varOmega ^{\rm{f}}}} {{{\varPhi}} _I^{\rm{f}}{{(\partial {{\varPhi}} _J^{\rm{f}}} / {\partial x_i^{\rm{f}}}}){\rm{d}}\varOmega }, \\ &{}^nK_{IJ}^{\rm{f}}{\rm{ = }}\int\limits_{{\varOmega ^{\rm{f}}}} {\frac{{\partial ({}^nv_k^{\rm{f}}{{\varPhi}} _I^{\rm{f}})}}{{\partial x_k^{\rm{f}}}}{{\rm{\rho}} ^{\rm{f}}}\frac{{\partial ({}^nv_j^{\rm{f}}{{\varPhi}} _J^{\rm{f}})}}{{\partial x_j^{\rm{f}}}}{\rm{d}}\varOmega },\; {}^nF_{Ii}^{{\rm{f}},t}{\rm{ = }}\int\limits_{{\varGamma ^{\rm{f}}}} {{{\varPhi}} _I^{\rm{f}}{}^n\tau _{ij}^{\rm{f}}n_j^{\rm{f}}{\rm{d}}\varGamma } ,\\ &{}^nF_{Ii}^{{\rm{f}},{{g}}}{\rm{ = }}\int\limits_{{\varOmega ^{\rm{f}}}} {{{\varPhi}} _I^{\rm{f}}{{\rm{\rho}} ^{\rm{f}}}{g_i}{\rm{d}}\varOmega } . \end{split} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!$

式中： ${{\varPhi }}_I^{\rm{f}}$、 ${{\varPhi }}_J^{\rm{f}}$分别为单元第I、J个流体节点的形函数， ${}^n\tau _{ij}^{\rm{f}}$为流体单元偏斜应力， $\varOmega^{\rm f} $为流体域.

CBS算法首先通过式（7）采用典型的伽辽金弱形式离散Navior-Stokes方程，求得中间时刻的速度；速度和压力采用相同的插值形函数 ${{\varPhi}} _I^{\rm{f}}$，通过式（8）得到压力；再通过式（9）得到修正后的速度值，其详细推导过程可参考文献[26]、[27].

固体域某一点的位移 ${{u}}_{}^{\rm{s}}$、速度 ${{v}}_{}^{\rm{s}}$和加速度 ${{a}}_{}^{\rm{s}}$可以通过插值得到：

(10) ${{u}}_{}^{\rm{s}} = \sum\limits_I {{{\varPhi}} _I^{\rm{s}}{{u}}_I^{\rm{s}}} {\rm{ }},\;{\rm{ }}{{v}}_{}^{\rm{s}} = \sum\limits_I {{{\varPhi}} _I^{\rm{s}}{{v}}_I^{\rm{s}}} {\rm{ , }}\;{{a}}_{}^{\rm{s}} = \sum\limits_I {{{\varPhi}} _I^{\rm{s}}{{a}}_I^{\rm{s}}}. $

式中： ${{\varPhi}} _I^{\rm{s}}$为固体插值形函数，在有限元法中采用线性插值，在S-PIM中采用点插值形函数.

采用标准的伽辽金弱形式离散固体运动方程，固体的运动离散方程为

(11) ${{M}}_{IJ}^{\rm{s}}{{a}}_J^{\rm{s}} = {{f}}_I^{{\rm{s,ext}}} - {{f}}_I^{{\rm{s,int}}}.$

式中： ${{M}}_{IJ}^{\rm{s}}$为集中质量阵； ${{f}}_I^{{\rm{s,ext}}}$、 ${{f}}_I^{{\rm{s,int}}}$分别为外力和内力. 其中内力可以表示为

(12) ${{f}}_I^{{\rm{s,int}}}{\rm{ = }}\int\limits_{{}^0\varOmega _i^{{\rm{sd}}}} {{{B}}_I^{}\left\{ {{S}} \right\}} {\rm{d}}\varOmega .$

式中： ${{S}}$为第2类皮奥拉-基尔霍夫应力，其具体求解可以参考文献[28]； $^0{\varOmega}^{\rm{sd}}_i$为光滑固体域； ${{{B}}_I}$为应变-位移矩阵，可以表示为

(13) ${{{B}}_I} = {{{B}}_{I0}}{\rm{ + }}{{{B}}_{I1}}.$

式中： ${{{B}}_{I0}}$为初始状态应变， ${{{B}}_{I1}}$为变形增量.

(14) $\left. {\begin{split} & {{{{B}}_{I0}}{\rm{ = }}\left[ {{{B}}_{I0}^{{J_1}}, \cdots ,{{B}}_{I0}^{{J_{{{{N}}^{\rm{s}}}}}}} \right]} ,\\ & {{{{B}}_{I1}} = \left[ {{{B}}_{I1}^{{J_1}}, \cdots ,{{B}}_{I1}^{J_{{{{N}}^{\rm{s}}}}^{}}} \right]} . \end{split}} \right\} $

式中： ${J_i}$为该域内的第 $i$个支点， ${{{N}}^{\rm{s}}}$为支点总数.

光滑点插值方法（smoothed point interpolation method，S-PIM）与有限元法的主要区别在于光滑域的构造以及采用梯度光滑技术进行应变场的光滑操作. 根据无重叠无间隙的原则构造不同形式的光滑域，形成不同类型且具备不同特性的光滑点插值方法. 如图1所示为边基光滑域的构造格式，由三角形单元边的节点和其相邻单元的中心点连接而成^[21]. 如图2所示为点基光滑域的构造格式，由节点周围单元的中心点及边的中点连接而成^[19]. 在点基光滑域基础上，可以进一步构造点基局部光滑域. 如图3所示为点基局部光滑域构造格式，将点基光滑域分割成2个部分，基于点基局部光滑技术分别采用点基光滑点插值法和有限元法进行计算^[19].

假设三角形单元边长为 $2b$，NPS-PIM的光滑子域长为 $(1 - \alpha )b$，则有限元法的边长为 $2\alpha b$， $\alpha $为边长比. 计算域面积为

(15) $\tilde A_{\rm{L}}^{\text{ε}} = {\alpha ^2}A_{\rm{L}}^{\rm{s}},\;\tilde A_{\rm{L}}^{\rm{s}} = (1 - {\alpha ^2})A_{\rm{L}}^{\rm{s}};\;\alpha \in [0,1.0].$

式中， $A_{\rm{L}}^{{{\rm{s}}}}$、 $\tilde A_{\rm{L}}^\text{ε}$、 $\tilde A_{\rm{L}}^{{{\rm{s}}}}$分别为整个光滑域面积、FEM计算域面积和NS-PIM计算域面积. ${\alpha ^2} = 0$表示整个计算域都采用NS-PIM进行插值计算； ${\alpha ^2} = 1.0$表示整个计算域都采用有限元法进行计算. Zhang等^[19]研究该参数对计算结果的影响，并将 ${\alpha ^2} = 0.4$作为推荐值.

在光滑域内可以采用梯度光滑技术软化模型刚度，光滑后的位移梯度可以表示为

(16) $ \begin{split} \bar u_{i,j}^{\rm{s}}({x_{\rm{L}}}) =& \displaystyle\sum\limits_I {\left[{{(1} / {A_i^{{\rm{sd}}}}})\int_{\varGamma _i^{{\rm{sd}}}} {{{\varPhi}} _I^{\rm{s}}({x^{\rm{s}}})} {{n}}_{}^{{\rm{sd}}}{\rm{d}}\varGamma \right]} {{u}}_I^{\rm{s}} = \\ & \displaystyle\sum\limits_I {{{(1} / {A_i^{{\rm{sd}}}}})} \sum\limits_{{n}} {{W_n}\left( {{{\varPhi}} _I^{\rm{s}}({x^{\rm{s}}}){{n}}_{}^{{\rm{sd}}}} \right)} {{u}}_I^{\rm{s}}. \end{split} $

式中： $A_i^{{\rm{sd}}}$为第i个光滑域的面积， $\varGamma _i^{{\rm{sd}}}$为光滑域边界， $x^{\rm{s}}$为光滑域内节点， $I$为光滑域内节点个数， $n$为光滑域内高斯点个数， ${W_n}$为高斯点的权重系数， ${{n}}_{}^{{\rm{sd}}}$为光滑域边界的外法线矢量， $\bar u_{i,j}^{\rm{s}}({x_{\rm{L}}})$为光滑域内光滑位移梯度.

在采用非贴体网格进行流固耦合信息传递时，须进行数值插值，可以通过搜索算法建立点与单元的对应关系. 在施加流固耦合速度条件时，虚拟流体节点的速度可以从固体域插值得到：

(17) $ \begin{split} {}^{n + 1}{{\hat v}}_{Ii}^{{\rm{fc}}} = \sum\limits_{L = I,J,K} {{}^{n + 1}{{\varPhi}} _L^{\rm{s}}({}^{n + 1}x_I^{\rm{f}})} {}^{n + 1}{{v}}_{Li}^{\rm{s}}{\rm{ }};\; ^{n+1}x_I^{\rm{f}} \in {}^{n + 1}{{{\varOmega}} ^{{\rm{fc}}}}. \end{split} $

式中： ${}^{n + 1}x_I^{\rm{f}}$为固体单元所覆盖的流体节点，I、J、K为固体单元节点， ${}^{n + 1}{{v}}_{Li}^{\rm{s}}$为固体节点速度， ${}^{n + 1}{{\varPhi}} _L^{\rm{s}}({}^{n + 1}x_I^{\rm{f}})$为该固体单元在点 ${}^{n + 1}x_I^{\rm{f}}$处的插值形函数； ${}^{n + 1}\hat {{v}}_{Ii}^{{\rm{fc}}}$为虚拟流体域流体节点速度.

对于作用于流体域的流固耦合力，根据式（2）、（9）可以得到流固耦合力的充分离散形式：

(18) ${}^{n + 1}F_i^{{\rm{f,FSI}}} = M_{IJ}^{\rm{f}}\frac{{{}^{n + 1}v_{Ji}^{\rm{f}}{ - ^n}v_{Ji}^{\rm{f}}}}{{\Delta t}} - {}^{n + 1}{\rm{RHS}}_{Ii}^{\rm{f}}.$

${}^{n + 1}F_i^{{\rm{f,FSI}}}$为作用于流体域的流固耦合力，基于欧拉网格计算得到，而作用于固体上的流固耦合力基于虚拟流体区域内拉格朗日网格计算得到：

(19) $ \begin{split} {}^{n + 1}F_i^{{\rm{s,FSI}}} = & - {}^{n + 1}F_i^{{\rm{fc,FSI}}} \\ {\rm{= }}& - M_{IJ}^{{\rm{fc}}}\frac{{{}^{n + 1}v_{Ji}^{{\rm{fc}}}{ - ^n}v_{Ji}^{{\rm{fc}}}}}{{\Delta t}} - {}^nC_{IJ}^{{\rm{fc}}}{}^nv_{Ji}^{{\rm{fc}}} - \\ &{}^nF_{Ii}^{{\rm{fc}}} - {{(\Delta t} / 2}){}^nK_{IJ}^{{\rm{fc}}}{}^nv_{Ji}^{{\rm{fc}}} + G_{IJi}^{{\rm{fc}}}{}^np_J^{{\rm{fc}}}. \end{split} $

式中：fc表示虚拟流体； ${}^{n + 1}v_{Ii}^{{\rm{fc}}} = {}^{n + 1}v_{Ii}^{\rm{s}},{}^nv_{Ii}^{{\rm{fc}}} = {}^nv_{Ii}^{\rm{s}}$，拉格朗日节点上的压力和剪切力可以通过插值得到.

浸没光滑点插值法求解流固耦合问题的具体流程如下：1）离散流体和固体域，并进行变量初始化；2）施加流固耦合力条件计算固体运动方程；3）根据单元与节点对应关系通过插值施加流固耦合速度条件，并求解流体运动方程；4）根据流体压力和剪切力计算流固耦合力；5）更新变量，并进行下一时间步长的计算.

如图4所示，本算例考虑流域内受到重力作用沉降的圆盘，以验证流固耦合算法的可靠性. 圆盘的直径2R=0.25 cm，流体域的宽度W=2.0 cm，高度H=5.0 cm. 在初始时刻，圆盘中心距离流体域的顶端和左侧距离L=1.0 cm. 重力加速度g=980 cm/s²，流域四周固壁为无滑移边界，流体速度满足 $v_1^{\rm{f}} = v_2^{\rm{f}} = 0$，并且顶端压力边界条件为 ${{{p}}^{\rm{f}}} = 0$. 在离散时流体和固体域均采用三角形单元，流体域划分为16281个节点和32000个单元；固体域划分为101个节点和168个单元. 固体材料参数设置如下：密度 ${{\rm{{\rm{\rho}} }}^{\rm{s}}} = 2$ g/cm³，弹性模量 ${{{E}}^{\rm{s}}} = 1\;000$ g/（cm·s²），泊松比 ${{\rm{\nu }}^{\rm{s}}} = 0.3$. 流体参数如下：密度 ${{\rm{{\rm{\rho}} }}^{\rm{f}}} = 1$g/cm³，动力黏性系数 ${{\rm{{\rm{\mu}} }}^{\rm{f}}} = 1.5$ g/（cm·s）. 本算例采用有限元法求解固体，因为圆盘在运动过程中变形较小，不同固体求解器对结果影响较小；流体和固体均采取相同的时间步长 ${\rm{d}}t = {10^{ - 5}}$ s.

将圆盘下落过程中竖直y方向的速度与经验公式进行数值算法验证对比，经验公式具体形式^[29]为

(20) $\begin{split}v_2^* = &\frac{{\left( {{{\rm{{\rm{\rho}} }}^{\rm{s}}}{\rm{ - }}{{\rm{{\rm{\rho}} }}^{\rm{f}}}} \right)g{R^2}}}{{4{{\rm{{\rm{\mu}} }}^{\rm{f}}}}} [ \ln\; \left( {{W / {2R}}} \right) - 0.915\;7 + \\ & 1.724\;4{\left( {{{2R} / W}} \right)^2} - 1.730\;2{\left( {{{2R} / W}} \right)^4} ] . \end{split}$

如图5（a）所示为该问题的计算简化模型，如图5（b）所示为该模型的有限元网格划分示意图. 在图5中，圆盘最终稳定的速度略大于参考解的速度. 这是由于参考解并未考虑圆盘弹性的作用，在本研究中，圆盘的弹性降低了周围流体对圆盘的作用力，从而增加了固体最终稳定的速度. 故圆盘下落过程中竖直方向的稳定速度略大于经验公式（式（20））的计算结果.

圆盘在下落过程的不同时刻处，所在不可压缩黏性流体中的x方向的速度云图如图6所示. y方向的速度云图如图7所示. 流场的压力云图如图8所示. 可以看出，只有在圆盘周围的流体的速度较大，且随着下落时间的增加，周围流场的速度随圆盘的速度而变化. 图8中的压强场有静压和动压，静压不随着时间变化，动压随着速度变化. 圆盘在下落过程中会带动周围流体运动，从而会导致其周围流场产生动压，故整个流体域在圆盘周围的压力场是变化的，而流体压力场其他部分变化不大. 说明圆盘下落的情况符合实际情况，证明用浸没式光滑点插值方法处理二维不可压缩黏性流体的流固耦合问题的可靠程度较高.

如图9所示，一底端固定的变形梁置于流域内. 流域尺寸为L=4 cm， H=1 cm；梁尺寸为a=0.04 cm，长度b=0.80 cm. 在本算例中忽略重力，流域底端满足速度无滑移条件，顶端满足法向速度 $v_2^{\rm{f}} = 0$；流域左端入口速度满足：

(21) $\left. \begin{array}{c} v_1^{\rm{f}}\left( t \right) = 1.5\left( { - {y^2} + 2y} \right) ,\\ v_2^{\rm{f}}\left( t \right) = 0 . \\ \end{array} \right\}$

流域右端压力边界条件为 ${{{p}}^{\rm{f}}} = 0$. 固体材料属性如下：密度 ${{\rm{\rho}} ^{\rm{s}}} = 7.8\;{\rm{ g/c}}{{\rm{m}}^3}$，弹性模量 ${E^{\rm{s}}} = {10^5}\;{\rm{ g/(cm}} \cdot {{\rm{s}}^2})$，泊松比 ${\nu ^{\rm{s}}} \!\!=\!\! 0.3$. 流体材料属性如下：密度 ${{\rm{{\rm{\rho}} }}^{\rm{f}}}\!\! =\!\! 1.0\;{\rm{ g/c}}{{\rm{m}}^3}$，动力黏性系数 ${{\rm{{\rm{\mu}} }}^{\rm{f}}} = 0.1\;{\rm{ g/(cm}} \cdot {\rm{s}})$.

对于该流固耦合问题，流体域和固体域在离散时均采用三角形单元，固体域设置5组网格，具体的网格信息如表1所示；流体域设置2组网格，具体网格信息如表2所示. 在计算时流体和固体均采取相同的时间步长，前4套固体网格下时间步长 ${\rm{d}}t = {10^{ - 4}}$，第5套固体精细网格下时间步长 ${\rm{d}}t = 5 \times {10^{ - 5}}$.

在计算时分别采用FEM、ES-PIM、NPS-PIM作为固体求解器. 在网格设置为MF1、MS1时，计算不同方法下梁端点的位移时间历程曲线. 由于文献[15]也使用浸没方法，网格划分采用三角形单元，使用41438个节点划分流体域，1457个节点划分固体域，较精密，并且最终得到收敛的结果，本研究将其作为参考解进行比较，结果如图10所示. 图中，u^s为梁端点的位移. 可以看出，模型刚度软化的光滑点插值法结果好于刚度偏硬的FEM的结果；基于点基局部光滑技术的NPS-PIM构造更为适宜的模型刚度，结果比ES-PIM更接近参考解. 如图11所示为流固耦合系统稳定时的固体位移云图. 图中，u_x为固定位移. 可以看出，NPS-PIM的顶端位移结果偏大. 如图12所示为3种方法计算的流场水平速度云图和流线图. 可以看出，在梁上方和后方分别形成高速区域和低速区域；3种方法下的流场云图并无太大区别，梁周围流场的过渡较光滑，整个流域的流场速度云图和流线图较光顺.

为了进一步比较3种方法的计算效率，定义二范数误差：

(22) $e_u^{\rm{s}} = {\left\| {u_i^{{\rm{s,num}}} - u_i^{{\rm{s,ref}}}} \right\|_{{L_2}}}\bigg/{\left\| {u_i^{{\rm{s,ref}}}} \right\|_{{L_2}}}.$

式中： $u_i^{{\rm{s,num}}}$、 $u_i^{{\rm{s,ref}}}$分别为数值解和参考解.

数值解为前4套固体域划分网格尺寸的计算结果；取固体域加密网格情况下，即当固体网格单元尺寸为1/300时，且固体求解器为有限元方法时的计算结果作为参考解.

综合考虑计算误差和CPU计算时间，定义计算效率^[18]为

(23) $\eta = {1 / {(t e_u^{\rm{s}})}}.$

式中：3种方法的计算时间t使用的是各自的实际CPU计算时间；计算误差为式（22）得出的计算误差.

通过归一化处理，令有限元方法的计算效率为1.000，其他2种方法的值取对于有限元方法的相对值. 第1种流体网格对应的不同固体求解器下的计算时间、误差以及计算效率的相对值分别如表3、4、5所示；第2种流体网格下的计算效率如表6所示.

在本算例中，由表3~6中的数据分析可知：ES-PIM和NPS-PIM的计算误差明显小于FEM，且NPS-PIM的计算误差最小；ES-PIM和NPS-PIM的计算时间大于FEM. 选取不同的流体网格，流固耦合模型采用NPS-PIM的计算效率均高于FEM；当流体域网格节点数为10 251时，ES-PIM的计算效率高于FEM；当流体域网格节点数为40 501 时，ES-PIM的计算效率随固体域逐渐加密而低于FEM. NPS-PIM在3种固体求解器中的计算效率最高，并且这种趋势随着固体网格或流体网格的加密而更加明显.

顶腔驱动流体运动并作用于超弹性墙，是经过广泛验证的经典流固耦合算例，许多学者都对此进行了验证^{[15, 18, 30]}. 如图13所示，方腔尺寸如下：L₁=2 m，位于方腔底端的弹性墙长度为L₂=2 m，高度H=0.5 m. 超弹性墙采用Neo-Hookean材料模型，密度 ${{\rm{{\rm{\rho}} }}^{\rm{s}}} = 1.0\;{\rm{ kg/}}{{\rm{m}}^3}$，材料常数 ${C_{10}} = 0.1\;{\rm{ kg/(m}} \cdot {{\rm{s}}^2})$， ${C_{01}} = 0$，体积模量 $K = 0$. 流体材料属性如下：密度 ${{\rm{{\rm{\rho}} }}^{\rm{f}}} = 1.0\;{\rm{ kg/}}{{\rm{m}}^3}$，流体的动力黏性系数 ${{\rm{{\rm{\mu}} }}^{\rm{f}}} = 0.2\;{\rm{ kg/(m}} \cdot {\rm{s}})$. 方腔顶盖的驱动速度满足：

(24) $v_x^{\rm{f}} = \left\{ \begin{array}{*{20}{l}} {\rm{0}}{\rm{.5}} {\sin ^2}\;({\text{π}} {x_1}/0.6),&{x^{\rm{f}}} \in {\rm{[0}}{\rm{,0}}{\rm{.3]}} ;\\ 0.5,&{x^{\rm{f}}} \in {\rm{[0}}{\rm{.3,1}}{\rm{.7]}}; \\ {\rm{0}}{\rm{.5}} {\rm{si}}{{\rm{n}}^{\rm{2}}}\;{\rm{(}}{\text{π}} ({x_1} - 2)/0.6{\rm{)}},&{x^{\rm{f}}} \in {\rm{[1}}{\rm{.7,2}}{\rm{.0]}} . \end{array} \right.$

弹性墙顶部自由运动，其他边界均为固壁，流域其他边界满足不可滑移速度边界条件，底部中点压力设置为参考值0. 对于该流固耦合问题，流体域和固体域在离散时均采用规则三角形单元. 在计算时流体和固体均采取相同的时间步长，前2套固体网格下时间步长 ${\rm{d}}t = {10^{ - 4}}$，第3套精细固体网格下时间步长 ${\rm{d}}t = 5 \times {10^{ - 5}}$.

在计算时分别采用FEM、ES-PIM和NPS-PIM作为固体求解器. 流体域网格采用16641个节点，固体域网格划分信息如表7所示.

数值解为前2套固体域划分网格尺寸的计算结果；参考解取固体域加密网格情况下，即当固体网格单元尺寸为1/100 cm时，固体求解器为有限元方法时的计算结果. 计算误差采用式（22），计算时间为具体的CPU计算耗时，处理结果如表8、9所示；计算效率采用式（23），归一化结果如表10所示.

如图14所示为NPS-PIM计算求得的流线分布和水平速度云图，将弹性墙顶部变形与圆圈表示的参考解^[30]对比，结果较吻合，流域底部弹性墙变形也较光顺. 由表8、9中数据可知，ES-PIM和NPS-PIM的计算误差明显小于FEM，NPS-PIM的计算误差最小；ES-PIM和NPS-PIM的计算时间长于FEM，ES-PIM的计算时间最长. 采用ES-PIM和NPS-PIM的计算效率均高于FEM，NPS-PIM在3种固体求解器中的计算效率最高.

（1）选用适当的固体求解器可以显著提高流固耦合问题的求解精度和计算效率.

（2）相比于有限元方法作为固体求解器时产生过硬的刚度模型，光滑点插值方法通过梯度光滑技术可以软化模型刚度，可以更准确地模拟流固耦合问题固体结构的模型刚度，从而有效提高求解精度，同时提高求解效率.

（3）点基局部光滑点插值法因构造了更为适宜的模型刚度，在流固耦合模型中作为固体求解器的求解精度和计算效率均得到显著提高，在本研究所采用的3种固体求解器中表现最好.