引 言

流固耦合现象广泛存在于自然界及各类工程问题中, 如昆虫与鸟类的飞行^[1]、血流和血管壁的相互作用^[2]、高层建筑、高耸结构及桥梁结构的风致振动^[3-4]、飞行器的颤振与运动稳定性^[5]以及医疗器械的设计与优化^[6-7]等．在工程问题中, 由于流场与弹性结构相互作用会诱发许多不良现象 (如静稳定性发散、颤振、极限环振荡及涡激振动等), 这些现象可能引起结构破坏或疲劳损伤, 因此在工程设计与计算中考虑流固耦合作用对于结构安全具有重要意义．

流固耦合涉及到不同物理场的多尺度(时间和空间)非线性耦合, 其快速准确求解一直以来都是力学与工程领域的一个难题^[8]．目前工程上流固耦合问题的求解策略主要分为两类: 分区算法和整体算法^[9]．在分区算法中, 流体和固体区域分别独立地进行求解并通过交界面交换位移和作用力等数据, 这类方法具备良好的程序模块性, 可以直接利用现有的计算流体力学 (CFD) 和计算固体力学 (CSD) 应用软件及程序, 且对计算机内存要求较低, 是目前流固耦合问题的主流求解方法^[10]．

分区算法又分为分区显式(弱耦合)算法^[11-12]和分区隐式(强耦合)算法^[13-14], 其中分区显式算法在每个时间步对流体和固体方程只分别求解一次, 由于存在时间滞后效应, 流固界面上的守恒律无法满足, 因此时间步长必须取得非常小, 并且数值模拟结果有可能收敛到错误的解上. 分区显式算法主要适用于气动弹性等较大质量比的情形^[15]．分区隐式算法则在每时间步内交替迭代求解各单场方程直至收敛以满足界面守恒, 常用的求解方法主要有固定点法^[16]、Newton类方法^[17]和优化方法^[18]．分区隐式算法的这种交替迭代求解方式本质上是一种非线性Gauss-Seidel型算法, 它的收敛速度较慢, 并且容易出现不收敛的情况, 因此通常需要采用一些加速收敛技术如Aiken松弛算法^[19]．此外, 当流体和固体的密度非常接近时, 迭代耦合方法会出现由附加质量效应导致的数值不稳定性^[20]．

为了避免分区算法的这些问题, 近些年来整体算法逐渐受到越来越多的关注^[21-25]．整体算法又称全耦合算法, 其对流固控制方程以及界面守恒条件同时进行求解, 可较好地适应强附加质量效应, 以及不会恶化整体系统和各子系统的求解精度．此外, 由于比分区算法少了一层迭代, 收敛速度相对较快, 且更适合并行计算, 但其缺点是方程具有极强的非线性, 且已失去了子问题控制方程系统矩阵的特点和可利用之处, 不能直接利用现有的流体、固体算法和软件, 需要重新研究求解算法．文献[26-27]对分区算法和整体算法的计算性能进行了详细比较．

流固耦合问题整体求解的计算规模通常非常巨大, 且离散矩阵具有病态性, 因而对方程迭代求解技术要求很高, 目前已经发展了一些有效的方法(主要是Newton类方法)用于其离散非线性方程组及相应线性系统的求解^[28-29], 不过由于三维问题的计算量巨大, 并且现有的算法大多缺乏较好的并行可扩展性, 难以处理大规模问题的数值模拟, 因此目前关于三维问题的数值结果尚不多．

本文将针对三维非定常流固耦合问题, 构造一种适合大规模并行计算的整体求解并行算法．对于流固运动界面, 我们采用任意拉格朗日−欧拉 (arbitrary Lagrangian-Eulerian, ALE) 方法进行描述^[30]. 对于流体、固体及动网格方程, 采用非结构有限元进行统一离散, 然后构造一种结合Newton类非线性求解器, Krylov子空间算法类线性求解器以及基于区域分解的Schwarz预条件子^[31]的Newton-Krylov-Schwarz (NKS) 算法进行高效求解．NKS算法在求解大规模线性及非线性方程组时非常有效, 已经成功地应用于各类大规模科学计算问题, 我们在之前的工作中利用这类方法来构造求解流体动力学方程的并行算法^[32-33], 在数千核的计算规模中取得了良好的收敛性和并行可扩展性, 本文将之拓展到流固耦合问题, 通过标准测试算例验证其计算效果, 并对算法的收敛性、稳定性、鲁棒性和并行可扩展性等数值性能进行测试分析．

1.   控制方程

本文考虑不可压缩流体与弹性体的相互作用问题．下面我们基于ALE方法给出要求解的流固耦合控制方程及相应的边界条件, 对于流体、固体和网格等不同场的变量, 我们分别用上标$ f,\;s,\;m $来表示区分．

1.1   ALE描述及动网格方程

如图1所示, $\varOmega _t^f$和$\varOmega _t^s$分别表示$t$时刻的流体区域和固体区域, 流体和固体区域的交界面为$\varGamma _t^i = \partial \varOmega _t^f \cap \partial \varOmega _t^s$, 除去交界面, 将流体与固体区域各自的边界分别记作$\varGamma _t^f = \partial \varOmega _t^f\backslash \varGamma _t^i$, $\varGamma _t^s = \partial {{\varOmega}} _t^s\backslash \varGamma _t^i$．将初始时刻的区域形状作为参考构形并记作${\overset{\frown} \varOmega^f } = \varOmega _0^f$及${\overset{\frown} \varOmega^s } = \varOmega _0^s$.

图  1  ALE框架

Figure  1.  ALE configuration

下载: 全尺寸图片 幻灯片

对于固体变形, 我们采用Lagrange描述法．对于流体运动, 采用ALE描述方法．定义$\overset{\lower0.5 em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\varOmega } _{}^f$到$\varOmega _t^f$之间的ALE映射

${A_t}:\;\;\overset{\lower0.5 em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\varOmega } _{}^f \to \varOmega _t^f,\;\;{\boldsymbol{x}}({\boldsymbol{X}},t) = {\boldsymbol{X}} + {{\boldsymbol{d}}^m}({\boldsymbol{X}},t),\;\;{\boldsymbol{X}} \in \overset{\lower0.5 em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\varOmega } _{}^f$

其中${{\boldsymbol{d}}^m}({\boldsymbol{X}},t)$表示网格位移, 通过求解给定的网格运动方程得到．为了方便, 这里我们采用调和方程来描述网格位移, 即有

边界条件为

其中${{\boldsymbol{d}}^s}$表示交界面处的固体位移．除了调和方程, 也可以采用其他方程来描述网格位移, 例如弹性方程, 它对大变形的情况具有更好的适应性, 不过与此同时也引入了一些附加的待定参数．我们的测试结果显示, 使用调和方程对于一般的变形已经可以得到良好的计算结果．

1.2   流体方程

ALE坐标下, 不可压缩Navier-Stokes方程可写为

其中$ {\rho ^f} $和${{\boldsymbol{u}}^f}$分别表示流体的密度和运动速度, ${{\boldsymbol{u}}^m} = \dfrac{{\partial {{\boldsymbol{d}}^m}}}{{\partial t}}$表示网格运动速度,${{\boldsymbol{\sigma}} ^f} = - {p^f}{\boldsymbol{I}} +$ ${\mu ^f}\left[ {\nabla {{\boldsymbol{u}}^f} + {{(\nabla {{\boldsymbol{u}}^f})}^{\rm{T}}}} \right]$为流场的Cauchy应力张量, 这里$ {p^f} $为流体压强, $ {\mu ^f} $为流体的动力黏性系数．边界条件给定为

其中$\varGamma _{t,D}^f和\varGamma _{t,N}^f$分别为流体边界$\varGamma _t^f$中的Dirichlet边界部分和Neumann边界部分, ${{\boldsymbol{n}}^f}$为边界面的外法向单位向量, g^f和h^f为给定函数, ${{\boldsymbol{u}}^s} = \dfrac{{\partial {{\boldsymbol{d}}^s}}}{{\partial t}}$为固体在交界面处的速度．

1.3   固体方程

对于固体变形, 我们采用线弹性模型, 在Lagrange描述下的弹性动力学方程为

其中$ {\rho ^s} $和${{\boldsymbol{d}}^s}$分别表示固体的密度和位移, ${{\boldsymbol{\sigma}} ^s} = {\lambda ^s}(\nabla \cdot {{\boldsymbol{d}}^s}){\boldsymbol{I}} + {\mu ^s}\left[ {\nabla {{\boldsymbol{d}}^s} + {{(\nabla {{\boldsymbol{d}}^s})}^{\rm{T}}}} \right]$为固体的应力张量, 这里$ {\lambda ^s} $和$ {\mu ^s} $为弹性体的Lamé系数, 它们与杨氏模量$ {E^s} $和Poisson比$ {\nu ^s} $满足如下关系

为了方便求解, 将式(9)化为如下一阶导数形式

边界条件给定为

式(1), 式(4), 式(5), 式(10)统一构成流固耦合控制方程组, 式(2)~式(3), 式(6)~式(8), 式(11)~式(12)组成边界条件.

2.   数值方法

2.1   有限元离散

针对流固耦合方程, 空间方向我们采用非结构有限元进行离散, 为此我们首先给出控制方程的Galerkin变分形式．

对于动网格方程, 分别定义如下试探函数空间和测试函数空间

则方程(1)的变分形式为: 寻找${{\boldsymbol{d}}^m} \in V_{}^m$, 使得对任意${{\boldsymbol{\phi}} ^m} \in V_0^m$, 都有

注意流固界面处的位移连续性条件在$V_{}^m$中作了限制.

对于固体方程, 类似地定义相应函数空间

则方程(10)的变分形式为: 寻找${{\boldsymbol{d}}^s} \in V_0^s,\;{{\boldsymbol{u}}^s} \in V_0^s$, 使得对任意${{\boldsymbol{\phi}} ^s} \in V_0^s$与${{\boldsymbol{\psi}} ^s} \in V_0^s$, 都有

注意上式左边第三项利用到了流固界面处的应力连续性条件(12)．

最后, 对于流体方程, 定义函数空间

则方程(4)和(5)的变分形式为: 寻找${{\boldsymbol{u}}^f} \in V_{}^f$, ${p^f} \in P$, 使得对任意${{\boldsymbol{\phi}} ^f} \in V_0^f$与${q^f} \in P$, 都有

其中${\boldsymbol{\varepsilon}} ({{\boldsymbol{u}}^f}) = \dfrac{1}{2}\left[ {\nabla {{\boldsymbol{u}}^f} + {{(\nabla {{\boldsymbol{u}}^f})}^{\rm{T}}}} \right]$为应变率张量. 注意流固界面处的速度连续性条件在$V_{}^f$中作了限制．

对于变分式(13)~式(15), 采用非结构有限元进行离散, 其中网格方程和固体方程采用标准的线性有限元, 流体方程采用P1-P1稳定化有限元^[34]. 通过定义$V_{}^m,\;V_0^m,\;V_0^s,\;V_{}^f,\;V_0^f,\;P$等函数空间相应的有限元空间, 然后将式(13)~式(15)中的测试函数取为线性有限元的基函数, 即可得到有限元离散方程. 需要注意的是对于不可压缩流动问题, 当速度和压力采用等阶的P1-P1有限单元时, 其离散方程不满足Ladyzhenskaya-Babuška-Breezi (LBB) inf- sup 条件, 需要添加稳定化项．本文采用流线迎风Petrov-Galerkin (SUPG)方法^[35], 设流体区域的网格为$\varOmega _t^f = \left\{ K \right\}$, 则流体方程的有限元变分形式改为

其中$ B_h^f $为原变分式(15)的左端项, $ {\left( {\; \cdot \;,\;\; \cdot \;} \right)_K} $表示在单元K上的积分, ${\boldsymbol{u}}_h^f,\;p_h^f,\;{\boldsymbol{d}}_h^m$分别为有限元网格上的待解速度场、压力场和网格位移函数, ${\boldsymbol{\phi}} _h^f,\;q_h^f$为相应的有限元基函数, $L_h^f = {\rho ^f}\dfrac{{\partial {\boldsymbol{u}}_h^f}}{{\partial t}} + {\rho ^f}({\boldsymbol{u}}_h^f -$${\boldsymbol{u}}_h^m) \cdot \nabla {\boldsymbol{u}}_h^f + \nabla p_h^f,$${\boldsymbol{U}}_h^f = - {\tau _m}L_h^f$, $ {\tau _c},\;{\tau _m},\;{\tau _b} $为稳定化参数, 分别定义如下

式中${G_{ij}} = \displaystyle\sum\limits_{k = 1}^3 {\dfrac{{\partial {\xi _k}}}{{\partial {x_i}}}\dfrac{{\partial {\xi _k}}}{{\partial {x_j}}}}$, $\dfrac{{\partial {\xi _k}}}{{\partial {x_i}}}$为表征标准单元(参考单元)与物理单元之间坐标转换的Jacobian 矩阵; $\Delta t$为时间离散的步长; 常数${c_1}$一般取为时间方向离散格式的阶数．本文拟采用二阶全隐BDF格式, 故${c_1} = 2$; 常数${c_2} = 36$．

通过有限元离散后, 我们得到如下半离散形式的非线性常微分方程组

其中${\boldsymbol{y}}(t)$为整个计算区域的所有待求未知量, 即由网格位移${\boldsymbol{d}}_h^m$, 固体位移${\boldsymbol{d}}_h^s$及速度${\boldsymbol{u}}_h^s$, 流体速度${\boldsymbol{u}}_h^f$及压力${{p}}_h^f$一起组成的向量, M为质量矩阵, ${\boldsymbol{L}}\left( {{\boldsymbol{y}}(t)} \right)$为除时间导数项外的非线性函数项．注意由于流体区域的网格是随时间变化的, 因此这里质量矩阵也为${\boldsymbol{y}}(t)$的函数．

在时间方向, 采用二阶BDF全隐格式进行离散, 得

因此在每一时间步, 需要求解如下大规模非线性方程组

在启动步中, 先采用半步长的隐式Euler法计算出BDF2格式所需要的初值．

2.2   Newton-Krylov-Schwarz并行算法

在每一时间步, 需要求解大规模非线性方程组(18)．为了高效并行地进行求解, 本文构造一种基于区域分解法的Newton-Krylov-Schwarz (NKS) 并行求解算法, 主要包括非线性方程求解器、线性求解器和预条件子三部分．为了方便推导, 略去上标, 将方程组写为一般形式

采用带线搜索的非精确Newton法来求解非线性方程组(19), 假设第k迭代步的值为$ {{\boldsymbol{y}}^{(k)}} $, 则第k + 1迭代步的值更新为

其中步长$ {\alpha ^{(k)}} $采用线搜索方法^[36]来计算, Newton修正量$\delta {{\boldsymbol{y}}^{(k)}}$通过求解如下线性方程组得到

这里${{\boldsymbol{J}}_k} = \dfrac{{\partial {\boldsymbol{F}}({{\boldsymbol{y}}^{(k)}})}}{{\partial {\boldsymbol{y}}}}$为第k步的Jacobian矩阵, ${\boldsymbol{F}}({{\boldsymbol{y}}^{(k)}})$为第k步的残差．采用广义最小残差法 (generalized minimal residual method, GMRES) 求解大规模线性方程组(21)．GMRES属于一种Krylov子空间迭代法, 它具有占用内存小、收敛速度快的优点, 而且比经典迭代法更适合大规模并行计算．此外由于Newton修正量$\delta {{\boldsymbol{y}}^{(k)}}$实际上无需精确计算, 为了减小计算量, 我们只作近似求解, 迭代收敛条件设定为

其中$ {\eta _k} $为相对收敛误差, 当$ {\eta _k} \to 0 $时, 则上述算法退化为精确Newton法．

对于Jacobian矩阵${{\boldsymbol{J}}_k}$, 一般可以通过无矩阵化方法^[37]或者有限差分方法进行近似求解, 也可以通过解析方法精确计算．在我们之前关于流体问题的测试中, 采用精确的Jacobian矩阵可以使算法具有更好的收敛性和鲁棒性^[32], 因此, 本文对Jacobian矩阵也采用解析表达式进行计算．Jacobian矩阵的主要计算量为流体方程的非线性对流项及稳定化项, 我们注意到在这些项中对网格位移及速度的导数项的量级较小, 因此进行了相应的简化处理, 即忽略对网格位移及速度的导数, 这可以有效降低算法的计算量以及编程的工作量, 并且数值测试表明这种处理几乎不影响算法的收敛性．

最后, 为了加快收敛速度, 对线性方程组(21)进行预条件处理

其中${\boldsymbol{M}}_k^{ - 1}$为右预条件子, 即${{\boldsymbol{J}}_k}$的近似逆矩阵．本文采用基于区域分解方法的限制型加性Schwarz预条件子 (restricted additive Schwarz, RAS)^[31] 进行预处理, 其构造方式如下．首先, 将计算区域$\varOmega$划分成${n_p}$个非重叠的子区域${\varOmega _l}\left( {l = 1,\;2, \cdots ,{n_p}} \right)$, 这里${n_p}$等于并行计算中的MPI进程数; 然后, 通过在各个子区域外围增加$\delta $层相邻区域的网格单元, 将各子区域延拓成重叠型子区域$\varOmega _l^\delta$, 这里$\delta $为整数, 表示子区域间的重叠层数．

对于各个重叠型子区域$\varOmega _l^\delta$, 定义限制算子(矩阵)${\boldsymbol{R}}_l^\delta$, 表示整个计算区域$\varOmega$的待求解未知量到子区域$\varOmega _l^\delta$相应未知量之间的映射, 即有

其中I为与${\boldsymbol{y}}_l^\delta$同阶的单位矩阵．限制算子的转置表示相应的延拓算子, 于是对应各个重叠子区域$\varOmega _l^\delta$的局部Jacobian 矩阵可表示为

记${(\tilde {\boldsymbol{J}}_l^\delta )^{ - 1}}$为${\boldsymbol{J}}_l^\delta$的近似或精确逆矩阵, 即局部预条件子, 那么加性Schwarz 预条件子定义为

而限制型加性Schwarz 预条件子则定义为

其中$ {\boldsymbol{R}}_l^0 $表示整个计算区域$\varOmega$的待求解未知量到非重叠子区域${\varOmega _l}$相应未知量之间的限制算子．相比标准的加性Schwarz 预条件子, 限制型加性Schwarz 预条件子在限制或者延拓中不传递重叠区域对应的分量, 因此节省了通信时间, 并且收敛速度更快．因此本文采用限制型加性Schwarz 预条件子．

在预条件矩阵${\boldsymbol{M}}_k^{ - 1}$的构造过程中, 主要的计算量是求解子区域的近似逆矩阵${(\tilde {\boldsymbol{J}}_l^\delta )^{ - 1}}$．对于大规模并行计算, 预条件子的构建在整套算法的总计算时间中通常占有很大比例, 但是其不需要精确逼近被预处理的问题, 因此这里采用按节点分块的不完全LU分解算法 (point-block ILU) 来求解相关的子区域问题．ILU算法包括不同的填充层级．当填充层级较大时, ILU对Jacobian矩阵的逆有更好的逼近, 因此达到收敛时需要的线性迭代次数会较少, 但是与此同时, 较大的填充层级意味着花费在子区域预条件子构造的时间较多, 因此总的计算时间不一定就少. 为了达到较好的性能, 有个折衷的选择, 需要根据实际情况进行选取．

整个NKS算法的计算流程如下:

(1) 采用上一个时间步的解作为初始值

(2) 对$k = 0,\;1,\;2, \cdots $, 执行以下步骤, 直至收敛:

① 构造${\boldsymbol{F}}({\boldsymbol{y}}_{}^{(k)})$的Jacobian矩阵${\boldsymbol{J}}_k^{}$;

② 采用带有预处理的Krylov 子空间迭代法GMRES求解如下线性方程组

③ 通过线搜索获得步长$ {\alpha _k} $;

④ 更新迭代解

3.   数值算例与算法性能测试

本节首先采用标准测试算例对本文算法进行验证, 然后对算法的收敛性、鲁棒性及并行可扩展性等数值性能进行测试研究. 计算平台为国家超级计算广州中心的天河二号集群, 其计算节点的硬件配置为双路12核Intel(R) Xeon(R) E5 CPU@ 2.4 GHz, 64 G内存. 算法程序基于PETSc软件^[38]编写, 计算网格采用非结构四面体网格, 利用CUBIT软件^[39]生成, 网格并行分区采用ParMETIS软件^[40]实现. 计算过程中, Newton 法的收敛误差限rtol设定为1.0 × 10⁻⁶, 最大迭代次数设定为30; 线性求解器采用重启型GMRES, 相对误差限设为1.0 × 10⁻⁴, 重启次数设为400, 最大迭代次数设为2000; 默认情况下, 区域分解的重叠层数设为2, 子区域问题采用ILU(2)近似求解．

3.1   算法验证

采用文献[22,28]的弹性障碍物绕流标准测试计算模型, 如图2所示, 计算区域为$(0,\;L) \times (0,\;H) \times $ $\left( - \dfrac{W}{2},\;\dfrac{W}{2}\right)$, 其中各几何参数分别为: 槽道的长$L = 1.5\;{\text{m}}$、宽$W = 0.8\;{\text{m}}$、高$H = 0.4\;{\text{m}}$, 内部弹性阻碍物的长$l = 0.2\;{\text{m}}$、宽$w = 0.4\;{\text{m}}$、高$h = $$0.2\;{\text{m}}$, 弹性体与入口的距离为$e = 0.4\;{\text{m}}$; 流体的物理参数为${\rho ^f} = 1.0 \;\times {10^3}\;{\text{kg/}}{{\text{m}}^{\text{3}}}$, 动力黏性系数为${\mu ^f} = $ $1.0\;{\text{Pa}} \cdot {\text{s}}$; 弹性体的物理参数为${\rho ^s} = 1.0$ ${\text{kg/}}{{\text{m}}^{\text{3}}}$, 杨氏模量$ {E^s} = $$1.4 \;{\text{MPa}}$, Poisson比$ {\nu ^s} = 0.4 $, 相应的Lamé系数为${\lambda ^s} = 2 \;{\text{MPa}}$, ${\mu ^s} = 0.5\;{\text{MPa}}$．

图  2  计算区域

Figure  2.  Computational domain

下载: 全尺寸图片 幻灯片

入口处给定入流速度

其中$s(t)$是一个时间方向的光滑因子, 定义为

$\bar v$为平均流速, 这里取$\bar v = 1$．壁面采用无滑移边界条件, 出口处采用无应力条件${\boldsymbol{n}} \cdot {{\boldsymbol{\sigma}} ^f} = 0$; 弹性体的底座处给定Dirichlet边界条件${{\boldsymbol{d}}^s} = {{\boldsymbol{u}}^s} = 0$．

采用四个不同规模的网格进行计算, 网格单元数与自由度总数如表1所示．对于这四个不同网格, 分别采用24, 48, 96, 192个处理器核进行计算, 程序的运行内存, Newton法的平均迭代次数, 每Newton步的平均GMRES迭代次数, 以及平均每时间步的计算时间也同时在表1给出．

表  1  计算网格

Table  1.  Computational meshes

Mesh	1	2	3	4
cells	1.89 × 104	1.30 × 105	1.08 × 106	6.03 × 106
DoF	2.90 × 104	1.78 × 105	1.47 × 106	8.11 × 106
np	24	48	96	192
memory/MB	51.86	141.5	453.8	1140
Newton	2.0	2.0	2.0	2.0
GMRES	11.6	17.0	28.2	40.4
time step/s	1.18	4.15	19.88	77.72

下载: 导出CSV 

| 显示表格

为了验证数值结果的准确性, 我们计算位于弹性体顶点处的两个观测点${{P}}1\;(0.4,\;0.2,\; - 0.2)$及${{P2}}\;(0.5,\; 0.2,\; - 0.2)$的位移, 并与文献[22]的计算结果进行比较, 如表2, 3所示．可以看到本文计算结果与文献结果基本一致, 表明本文算法的求解结果是有效和可靠的．

表  2  观测点${{P}}1\;(0.4,\;0.2,\; - 0.2)$处的位移

Table  2.  Displacements at point ${P}1\;(0.4,\;0.2,\; - 0.2)$

	$d_x^{}$/m	$d_y^{}$/m	$d_z^{}$/m
mesh 1	1.984 × 10−3	5.152 × 10−4	−9.431 × 10−5
mesh 2	1.729 × 10−3	5.052 × 10−4	−6.124 × 10−5
mesh 3	1.602 × 10−3	4.832 × 10−4	−5.775 × 10−5
mesh 4	1.598 × 10−3	4.847 × 10−4	−6.785 × 10−5
Ref. [22]	1.63 × 10−3	5.05 × 10−4	−1.22 × 10−4

下载: 导出CSV 

| 显示表格

表  3  观测点${{P2}}\;(0.5,\;0.2,\; - 0.2)$处的位移

Table  3.  Displacements at point ${{P2}}\;(0.5,\;0.2,\; - 0.2)$

	$d_x^{}$/m	$d_y^{}$/m	$d_z^{}$/m
mesh 1	1.964 × 10−3	−5.377 × 10−4	−4.855 × 10−5
mesh 2	1.700 × 10−3	−4.683 × 10−4	−7.854 × 10−5
mesh 3	1.563 × 10−3	−4.251 × 10−4	−7.185 × 10−5
mesh 4	1.550 × 10−3	−4.240 × 10−4	−7.216 × 10−5
Ref. [22]	1.54 × 10−3	−4.65 × 10−4	−3.13 × 10−5

下载: 导出CSV 

| 显示表格

图3给出了$t = 3\;{\rm{s}}$时的弹性变形及流场速度分布图, 为了验证较大变形的情况, 我们同时也给出了弹性模量${E^s} =350$ kPa及${E^s} = 87.5$ kPa的结果．而当进一步降低弹性体的弹性模量时, 由于弹性体变形过大, 网格出现奇性, 从而导致计算发散. 注意这种发散不是数值格式不稳定性带来, 不能通过减小时间步长来消除．

图  3  不同杨氏模量下的弹性变形及速度分布图

Figure  3.  Deformation and velocity distribution with different Young's modulus

下载: 全尺寸图片 幻灯片

3.2   算法性能测试

首先对算法的收敛性能进行测试．影响算法收敛性能的参数主要有两个: 一是区域分解的重叠层数$\delta $; 二是构造子区域局部预条件子的线性求解器. 分别取$\delta = $1, 2, 然后采用不同填充级的子区域求解器ILU(k)进行计算, 测试结果如表4所示, 其中Newton表示平均每个时间步的非线性迭代次数, GMRES表示平均每个Newton步的线性迭代次数, Time表示平均每个时间步的计算时间(单位为秒)．可以看到, Newton迭代次数基本不受影响, 表明不同填充级的ILU线性求解器均可作为子区域问题的求解器．另外, 由于填充级增加时, ILU对Jacobian矩阵的逆有更好的逼近, 因此可以看到线性迭代次数显著随之减少．不过每时间步的计算时间并没有相应减少, 原因是当填充级增加时, 虽然线性迭代次数显著减少了, 但是花在子区域预条件子构造的时间相应增加了, 因此有个折衷的选择, 在本算例中, 填充级取为$k \leqslant 2$均可, 当$k = 3$时计算时间则显著增加．另外, 由于本算例中Newton法的收敛速度很快, 因此子区域的重叠层数对计算时间影响不大．

表  4  子区域的重叠层数及子问题求解器对算法性能的影响

Table  4.  Impact of the overlapping size and the subsolver on the NKS algorithm

$\delta $	Subsolver	Newton	GMRES	Time/s
1	ILU(0)	2.0	62.3	18.07
1	ILU(1)	2.0	40.0	18.15
1	ILU(2)	2.0	31.2	19.34
1	ILU(3)	2.0	29.2	34.98
2	ILU(0)	2.0	54.3	17.40
2	ILU(1)	2.0	34.4	18.23
2	ILU(2)	2.0	28.2	19.88
2	ILU(3)	2.0	24.5	37.50

下载: 导出CSV 

| 显示表格

接下来我们考察NKS算法的稳定性. 表5给出取不同时间步长时算法的收敛性能．由于时间方向采用了全隐离散格式, 因此算法具有很好的稳定性, 当时间步长$\Delta t$从0.001增大至0.128时, 虽然Newton迭代步数、线性迭代步数和每步的计算时间都略有增加, 但总的来说算法仍能快速收敛, 表明算法可以在精度允许的范围内自由地选择时间步长进行数值模拟．

表  5  不同时间步长对NKS算法收敛性能的影响

Table  5.  Performance of NKS with respect to $\Delta t$

$\Delta t$	Newton	GMRES	Time/s
0.001	2.0	19.5	21.34
0.002	2.0	22.6	21.54
0.004	2.0	26.2	21.97
0.008	2.0	39.0	22.96
0.016	2.0	47.2	23.99
0.032	2.0	55.9	24.25
0.064	2.6	71.8	32.98
0.128	3.0	70.4	38.17

下载: 导出CSV 

| 显示表格

随后对NKS算法的鲁棒性进行测试, 即考察取不同的物理参数对算法收敛性能的影响．表6给出了仅改变流体黏性系数, 而其他参数保持不变时算法收敛性能的测试结果, 表7给出仅改变固体密度时的测试结果, 表8则给出仅改变固体的杨氏模量和Poisson比的测试结果．可以看到在不同的物理参数下, 算法的收敛性能基本保持不变, 表明了算法具有良好的鲁棒性．

表  6  关于流体黏性系数的算法鲁棒性测试

Table  6.  Robustness of the algorithm with respect to the fluid viscosity

${\mu ^f}$	Newton	GMRES	Time/s
1.0	2.0	26.2	21.97
0.1	2.0	26.6	22.04
0.01	2.0	26.8	21.98

下载: 导出CSV 

| 显示表格

表  7  关于固体密度的算法鲁棒性测试

Table  7.  Robustness of the algorithm with respect to the solid density

${\rho ^s}$	Newton	GMRES	Time/s
1	2.0	60.8	25.29
10	2.0	53.2	24.41
100	2.0	45.6	23.75
1000	2.0	26.2	21.97
2000	2.0	22.4	21.13
4000	2.0	14.0	20.85

下载: 导出CSV 

| 显示表格

表  8  关于杨氏模量和Poisson比的算法鲁棒性测试

Table  8.  Robustness of the algorithm with respect to the Young’s modulus and Poisson’s ratio

${E^s}$	${\nu ^s}$	Newton	GMRES	Time/s
1.4 × 106	0.1	2.0	25.0	21.62
1.4 × 106	0.2	2.0	25.7	21.87
1.4 × 106	0.3	2.0	25.8	21.96
1.4 × 106	0.4	2.0	26.2	21.97
1.4 × 106	0.48	2.0	41.4	23.39
1.4 × 106	0.4	2.0	26.2	21.97
7.0 × 105	0.4	2.0	22.9	21.67
3.5 × 105	0.4	2.0	20.2	21.44
1.75 × 105	0.4	2.0	18.0	21.21
8.75 × 104	0.4	2.0	16.8	21.05

下载: 导出CSV 

| 显示表格

最后, 对算法的并行性能进行测试, 如表9所示, 其中n_p表示计算的进程数(等于处理器的核数), speedup为加速比, ideal表示理想加速比, efficiency表示并行效率．可以看到, 随着进程数的增加, 非线性迭代次数保持不变, 线性迭代次数则略有增加, 每迭代步的计算时间则显著减少．当进程数从192 增加到1536时, 算法取得了超线性的加速比, 当进程数增加到3072时, 并行效率为91%．算法的加速比曲线如图4所示, 测试结果表明了本算法具有良好的并行效率和可扩展性．

表  9  算法并行可扩展性测试

Table  9.  Parallel performance and scalability of the algorithm

${n_p}$	192	384	768	1536	3072
Newton	2.0	2.0	2.0	2.0	2.0
GMRES	45.7	47.0	49.3	52.4	53.5
time	83.41	35.07	17.20	9.65	5.73
speedup	1	2.38	4.85	8.64	14.56
ideal	1	2	4	8	16
efficiency	100%	119%	121%	108%	91%

下载: 导出CSV 

| 显示表格

图  4  算法加速比

Figure  4.  Parallel speedup of the algorithm

下载: 全尺寸图片 幻灯片

3.3   与其他算法性能的比较

本节将本文方法与文献中已有的一些算法进行比较. 对不同算法的性能直接进行对比是一项复杂的工作, 一方面由于不同文献通常采用了不同的测试算例, 另一方面算例的计算规模、算法的参数设置以及进行测试的机器也不相同. 本文选取平均每时间步的非线性迭代次数作为算法收敛性能的考量, 另外, 基于不同文献测试算例的问题规模不同, 参照文献[29]的做法, 我们定义一个时间步内“每CPU核在每秒钟内所求解的自由度数”这一指标来比较不同算法的计算效率. 最后, 为了考察算法处理大规模问题的能力, 对算法的并行性能也进行对比. 关于全耦合算法, 目前文献中的测试算例大多都是关于二维问题的, 三维算例较少, 在我们的检索范围内, 找到了文献[22, 28-29]这几种不同的数值方法, 这几种方法均采用Newton法作为非线性求解器, Krylov子空间迭代法作为线性求解器, 它们的主要区别在于预条件算子及其求解算法不同, 表10给出了本文方法与它们在收敛性能与并行性能方面的比较结果. 可以看到, 在不考虑算法收敛参数设定以及测试机器的情况下, 本文方法每时间步的Newton迭代次数最少, 每CPU核在每秒钟内所求解的自由度数仅次于文献[28]的方法, 在并行可扩展性与并行效率上, 本文方法远远高于其他方法, 表明在处理大规模问题上具有较好的潜能.

表  10  本文算法与其他文献算法的性能比较

Table  10.  Comparison of different numerical approaches found in literature

	Solver	Dofs	Newton steps	Ratio*	CPU cores	Parallel efficiency	Source
Ref. [22]	Newton-Block LDU	14 000 000	6 ~ 9	109.9	16 ~ 256	41%	Fig.14, Fig.15
Ref. [28]	Newton-Krylov- Multigrid-Richardson	120 902	4	1139	1	−	Tab. 1, Tab. 12
Ref. [29]	Newton-Multigrid	3 531 304	5.15	450.3	1 ~ 32	31%	Tab. 6, Tab. 9, Fig. 7
this work	Newton-Krylov- Schwarz	8 110 000	2	543.5	192 ~ 3072	91%	Tab. 1, Tab. 9
*Note: Here the “Ratio” indicates the degree of freedom persecond computed with one CPU core for each time step

下载: 导出CSV 

| 显示表格

4.   结 论

本文针对三维流固耦合问题, 提出了一种有限元数值求解的全隐全耦合可扩展并行算法. 为方便处理复杂的计算区域, 采用了三维非结构网格进行空间离散, 为使算法具有更好稳定性和鲁棒性, 采用二阶BDF全隐格式进行时间离散. 对于时空离散后得到的大规模非线性代数系统, 构造了基于区域分解的Newton-Krylov-Schwarz并行算法对流体、固体与动网格方程进行一次性整体求解. 采用弹性障碍物绕流问题的标准测试算例对算法的准确性进行了验证．算法性能测试结果表明本文提出的算法在数千处理器核条件下具有良好的并行可扩展性能. 通过测试算法在不同计算时间步长、不同物理参数下的收敛性能, 验证了算法具有良好的稳定性和鲁棒性. 另外, 本文对动网格方程采用了较为简单的调和方程模型, 数值结果表明其在大变形的情况下计算会发散, 为了避免这一问题, 可以考虑更复杂的弹性方程模型, 这些问题有待进一步深入研究.