引言

流固耦合广泛存在于土木工程、机械工程、海洋工程、航空航天工程和生物医学工程等众多工程领域^[1-9]. 以土木工程为例,超高层建筑、高耸结构、空间结构、大跨桥梁及长径线缆等钝体结构在强风作用下风致振动是典型流固耦合现象. 此类结构具有质量轻、柔性大、阻尼小和自振频率低等特点,流固耦合效应突出,风损风毁现象已屡见不鲜. 例如,1940年11月7日美国华盛顿州塔科马窄桥在较低风速(约为设计风速的50%)下发生颤振失稳、进而垮塌. 因此,准确考虑流固耦合作用对保障工程结构安全意义重大.

由于两相介质在交界面处存在耦合作用,一方变化会激发另一方的强烈响应,故解析求解流固耦合异常困难,目前多需依赖数值计算. 流固耦合数值仿真已成为当前工程计算与设计中的重要需求,引起学者和工程师的广泛关注. 一般而言,分区弱耦合算法效率较高,常用于气动弹性问题^[6]. 然而,分区弱耦合算法难以克服强附加质量效应,易导致数值失稳,故多用于高质量比刚体流致振动^[10]. 对于低质量比、大位移或柔性体等复杂情形,可采用分区强耦合算法进行求解.

分区强耦合算法在每时间步内依次迭代计算各单场方程直至收敛,具备良好的物理守恒性和程序模块性^[11],因而备受青睐. 其中,固定点法是一种广泛运用的强耦合方法. Le Tallec和ouro^[1]提出了固定点法的预条件最速下降松弛技术,兼具数值稳定与能量守恒,已用于复杂工业问题. Küttler和Wall^[13]根据界面位移构造Aitken动态松弛因子加速固定点迭代,极大扩展了该技术的应用范围. Bai 等^[14]利用ANSYS-CFX 软件实施固定点迭代计算湍流下三维桥梁断面颤振问题. 由于优化了网格运动,该方案具有工程应用潜力. Sun等^[15]也基于分区强耦合的模块化优势,运用商业软件计算桥梁结构风致振动问题. 笔者及其合作者提出新型混合界面耦合条件,采用固定点法数值模拟各类钝体流致振动问题^[16-18]. 有关流固耦合数值方法的最新研究进展可参考文献[19]. 需指出,笔者提出一种部分强耦合迭代算法,进一步提升计算效率^[17].

本文基于任意拉格朗日--欧拉(arbitrary Lagran- gian-Eulerian, ALE)有限元方法,灵活选择不同物理模型及单场求解手段,提出一种精确高效的分区强耦合算法,数值模拟大位移流固耦合问题,并揭示相关流致振动现象,为流固耦合数值仿真研究提供新思路.

1 流体模型

1.1 控制方程

考虑移动网格情形,不可压缩黏性流体的质量守恒律和动量守恒律可写为ALE描述下的无量纲化Navier-Stokes(NS)方程

$\nabla \cdot {\pb u} = {\bf 0} $(1)

$\dfrac{\partial {\pb u }}{\partial t} + {\pb c } \cdot \nabla {\r {\bf u}} +\nabla \cdot {\pb \siga }^{\r F} - {\pb f }^{\r F} = {\bf 0}$ (2)

式中,$\nabla $为梯度算子,${\pb u}$和$p$为流体速度和压力,对流速度${\pb c}={\pb u } - {\pbw }$,${\pb w}$为网格速度,${\pb \siga }^{\r F}$为流体应力张量,${\pb f}^{\r F}$为流体体力项. 流体本构方程表述为

$ {\pb \siga }^{\r F} = - p{\pb I} + \dfrac{1}{Re}\left( {\nabla {\pb u} + \left( {\nabla {\pb u} } \right)^{\r T} } \right) $ (3)

式中,${\pb I}$为单位矩阵,雷诺数${Re} = \rho^{\r F} UD/\u $,$\rho^{\r F}$为流体密度,$U$和$D$分别为来流特征速度和特征长度,$\u $为流体黏度,上标T表示转置.

初始条件和边界条件指定如下

${\pb u} = {\pb u }^0 \,, \ \ p = p^0$ (4)

${\pb u} = {\pb g }^{\r F} \,,{\pb \siga }^{\r F} \cdot {\pb n} = {\pb h }^{\r F}$ (5)

1.2求解过程

运用半隐式特征线分裂(characteristic-based split, CBS)算法^[0]求解NS方程(1)和(),主要步骤说明如下:

(a) 计算辅助速度场

$\babl {\pb u }^ * - {\pb u }^n = \Delta t\Bigg \{ - {\pb c } \cdot \nabla {\pb u } + \dfrac{1}{Re}\nabla^{\pb u } + \[] \qquad \dfrac{\Delta t}{}\left[ {{\pb c} \cdot \nabla \left( {{\pb c} \cdot \nabla{\pb u}} \right)} \right]\Bigg \}^n \eal\eqno(6)$

(b) 更新压力场

$\nabla ^p^{n + 1} = \dfrac{1}{\Delta t}\nabla \cdot {\pb u}^ * \eqno(7)$

(c) 修正速度场

${\pb u}^{n + 1} - {\pb u}^ * = \Delta t\left[ {\nabla p^{n + 1} -\dfrac{\Delta t}{}\left( {{\pb c } \cdot \nabla ^p} \right)^n} \right] \eqno(8) $

CBS算法的稳定系数为$(\Delta t)^{} /$,与单元尺寸无关,故动网格计算略为省时. 半隐式CBS方法的稳定性较好且对时间步$\Delta t$限制较小. 因CBS算法可对速度和压力采用等阶/低阶插值,选用三节点三角形(T3)单元进行有限元离散.

2 固体模型

2.1 刚体动力学

刚体位移记为${\pb d} = \{ d_{1}, d_{}, \theta \}^{\r T}$,其中$d_{i}$ ($i= 1, $)和 $\theta $ 分别表示平动和转动分量,如图1所示. 假定每个方向自由度相互解耦,则刚体运动方程的无量纲分量形式可写为

$ \ddot {d}_1 + 4\pi f_{{\r R}1} \zeta _1 \dot {d}_1 + \left( {\pi f_{{\r R}1} }\right)^d_1 = \dfrac{C_{\r D} }{_1^ * } (9)$

$ \ddot {d}_ + 4\pi f_{{\r R}} \zeta _ \dot {d}_ + \left( {\pi f_{{\r R}} }\right)^d_ = \dfrac{C_{\r L} }{_^ * } (10)$

$ \ddot {\theta } + 4\pi f_{{\r R}\theta } \zeta _\theta \dot {\theta } + \left( {\pi f_{{\r R}\theta } } \right)^\theta = \dfrac{C_{\r } }{_\theta ^ * } (11)$

式中,换算频率$f_{{\r R}i} = f_{{\r N}i} D / U$,转动换算频率$f_{{\r R}\theta } = f_{{\r N}\theta } D / U$,自然频率$f_{{\r N}i} = \sqrt { k_i / _i } / (\pi) $,转动自然频率$f_{N\theta } =\sqrt {k_\theta / I_\theta } / (\pi )$,阻尼比$\zeta _i = {c_i } / ( \sqrt {_i k_i })$,转动阻尼比$\zeta _\theta = {c_\theta } / (\sqrt {I_\theta k_\theta } ) $,质量比$_i^\ast = {_i } /(\rho ^FD^)$和 $_\theta ^\ast = {I_\theta } /(\rho ^{\r F}D^4)$,$_i $, $c_i $和$k_i$分别为平动方向的质量、阻尼和刚度,$I_\theta $, $c_\theta $和$k_\theta$分别为扭转方向的惯性矩、阻尼和刚度,$C_{\r D} $, $C_{\r L} $和$C_{\r }$为流体阻力系数、升力系数和冲量矩系数.

刚体平面运动还应满足协调条件^[1]. 该条件阐述了边界位移、速度和加速度与重心变量之间的几何关系(参见图1).

2.2. 弹性动力学

几何非线性固体的无量纲运动方程可写为如下形式

$\ddot{\pb d} - \dfrac{1}{^ * }\nabla \cdot {\pb \siga}^{\r S} - {\pb f }^{\r S} = {\bf 0 } \eqno(12)$

式中,位移${\pb d} = \{d_{1}, d_{}\}^{\r T}$,质量比$^{\ast } =\rho^{\r S}/ \rho^{\r F}$,${\pb \siga}^{\r S}$为Cauchy应力,${\pb f}^{\r S}$为结构体力项. St. Venant-Kirchhoff材料本构关系可写为

${\pb S} = {\pb C}:{\pb E}\,,{\pb E} = {\pb F}^{\r T}{\pb F} - {\pb I} \eqno(13)$

式中,${\pb S}$是第二Piola-Kirchhoff应力,${\pb C}$为弹性本构张量,${\pb E}$表示Green-Lagrangian应变,${\pb F}$为变形梯度. 第二Piola-Kirchhoff应力经构形转换可表示为Cauchy应力

${\pb S} = J{\pb F}^{ - 1}{\pb \siga }^{\r S}{\pb F}^{ - {\r T}} \eqno(14)$

其中,$J ={\r det} ({\pb F})$. 初始条件和边界条件如下

${\pb d} = {\pb d }^0\,,\dot{\pb d} = \dot{\pb d}^0\,, {\pb \siga }^{\r S} \cdot {\pb n} = {\pb h}^{\r S} \eqno(15)$

本文基于总体拉格朗日格式运用修正牛顿$\!$-$\!$-$\!$拉普森法^{[ ]}求解非线性方程(1). 另外,对变形梯度的处理详见下文.

2 .3 光滑有限元法

传统有限元法组装的结构刚度矩阵偏刚. 为此,这里采用光滑有限元方法^[3-4]求解式(1)线性化后的增量方程. 光滑处理能有效“软化”刚度矩阵,令有限元解更接近理论解,亦未引入新自由度,其基本原理叙述如下.

根据梯度光滑运算,任意场变量 X 的梯度可近似表示为

$\tilde {\nabla }X\left( {{\pb x}_{\r C} } \right) = \dint_{\varOega _{\r C} } {\nabla X\left( {\pb x}\right)} \Phi \left( {{\pb x} - {\pb x}_{\r C} } \right){\r d}\varOega \eqno(16)$

式中,${\pb x}$代表空间坐标,$\varOega_{\r C}$为光滑域,$\varPhi $为Heaviside型光滑核函数.

由高斯散度定理可知,式(16)可进一步写为

$\babl \tilde {\nabla }X\left( {{\pb x}_{\r C} } \right) = \dint_{\varGaa _{\r C} } {X\left( {\pb x} \right)} {\pb n}\left( {\pb x} \right)\varPhi \left( {{\pb x} - {\pb x}_{\r C} } \right){\r d}\varGaa - \[] \qquad \dint_{\varOega _{\r C} } {X\left( {\pb x} \right)} \nabla \varPhi \left( {{\pb x} - {\pb x}_{\r C} } \right){\r d}\varOega \end{array}\eqno(17) $

其中,$\varGaa_{C}$为$\varOega_{C}$边界,${\pb n}$是$\varGaa _{C}$外法向. 根据非负性和归一性要求,光滑核函数$\varPhi $可定义为

$\varPhi \left( {{\pb x} - {\pb x}_{\r C} } \right) = \left\{ \!\! \begin{array}{lc} {0,} & {\pb x} \notin \varOega _{\r C} \[] 1 / A_{\r C} \,, & {\pb x} \in \varOega _{\r C} \end{array} \right. \eqno(18)$

式中,$A_{\r C}$为光滑域面积.

将式(18)代入式(17)中并注意到常数导数为零,则标量$X$的梯度可写为

$\tilde {\nabla }X\left( {{\pb x}_{\r C} } \right) =\dfrac{1}{A_{\r C} }\dint_{\varGaa _{\r C} } {X\left( {\pb x}\right)} {\pb n}\left( {\pb x} \right){\r d}\varGaa \eqno(19)$

用一点高斯积分即可准确计算上述线积分

$\tilde {\nabla }X\left( {{\pb x}_{\r C} } \right) = \dsu_{i = 1}^{ nsd} \dfrac{1}{A_{{\r C}i} }X\left({{\pb x}_i^{{\r GP}} } \right){\pb n}\left( { {\pb x}_i^{{\r GP}} } \right) l_i \eqno(20)$

其中,${nsd}$为每单元内光滑域个数,${\pb x}^{\r GP}$为高斯点坐标,$l_{i}$为边长.

光滑域构造及其形函数值见后文.

2.4 时间积分

结构运动方程可写为如下一般形式

${\pb }\ddot {\pb d} + {\pb C}\dot {\pb d} + {\pb K \pb d} = {\pb F}_{{\r EX}} \eqno(21)$

式中,${\pb }$, ${\pb C}$和${\pb K}$分别表示结构质量矩阵、阻尼矩阵和刚度矩阵,${\pb F}_{\r EX}$为外力向量. 时间积 分方案选用复合隐式时间积分法^[5-7],其基本思想是：将一时间步分为两子步,分别运用Newark-$\beta $法^[8]和三点向后欧拉差分.

主要原理简述如下.

记$[t, t+\Delta t] =\Delta t[n, n + 1] = \Delta t[n,n +\alpha ] \cup \Delta t[n + \alpha $, $n + 1]$,其中0＜a＜1. 在第一个子步内将结构运动方程于 t^n+α

${\pb }\ddot {\pb d}^{n + \alpha } + {\pb C}\dot {\pb d}^{n + \alpha } +

{\pb K \pb d}^{n + \alpha } = {\pb F}_{{\r EX}}^{n + \alpha } \eqno(22)$

对上式应用Newark-$\beta $法,可得

$\dot{\pb d}^{n + \alpha } = \dot{\pb d}^n + \left( {1 - \gaa } \right) \alpha \Delta t \ddot{\pb d}^{n + \alpha } (23)$

$ {\pb d }^{n + \alpha } = {\pb d }^n + \alpha \Delta t \dot{\pb d}^n + \left({1 / - \beta } \right)\left( {\alpha \Delta t} \right)^\ddot{\pb d}^n + \beta \left( {\alpha \Delta t} \right)^\ddot{\pb d}^{n + \alpha } (24)$

其中, $\beta \geqslant 1/4$, $\gaa \geqslant 1/$. 由式(4)可得加速度

$\ddot{\pb d}^{n + \alpha } = \dfrac{1}{\beta \left( {\alpha \Delta t} \right)^}\left( {{\pb d }^{n - \alpha } - {\pb d }^n} \right) - \dfrac{1}{\beta \alpha \Delta t } \dot{\pb d}^n -\left( {\dfrac{1}{\beta } - 1} \right) \ddot{\pb d}^n (25)$

将式(3) $\si$式(5)代入式(),那么

$\left( {\dfrac{1}{\beta \left( {\alpha \Delta t} \right)^}{\pb } +\dfrac{\gaa }{\beta \alpha \Delta t }{\pb C} +{\pb K }} \right){\pb d }^{n + \alpha } = {\pb F }_{{\r EX}}^{n + \alpha } + {\pb }\left[ {\dfrac{1}{\beta \left( {\alpha \Delta t} \right)^} {\pb d}^n + \dfrac{1}{\beta \alpha \Delta t }\dot{\pb d}^n + \left( {\dfrac{1}{\beta } - 1} \right) \ddot{\pb d}^n}\right] + {\pb C }\left[ {\dfrac{\gaa }{\beta \left( {\alpha \Delta t} \right)^}{\pb d }^n + \left( {\dfrac{\gaa }{\beta } - 1} \right) \dot {\pb d}^n + \left( {\dfrac{\gaa }{\beta } - 1} \right) {\alpha \Delta t} \ddot{\pb d}^n} \right] (26)$

同理,在第二个子步内将结构运动方程于$t^{n + 1}$时刻离散

${\pmb M}\ddot{\pmb d }^{n + 1} + {\pmb C}\dot{\pmb d }^{n + 1} + {\pmb K}{\pmb d }^{n + 1} = {\pmb F }_{{\rm EX}}^{n + 1} \eqno(27)$

因任意变量的时间导数$\dot {Q}^{n + 1}$可根据$Q^n$, $Q^{n + \alpha }$和$Q^{n + 1}$近似获得^[29]

$\ddot {Q}^{n + 1} = c_1 Q^n + c_2 Q^{n + \alpha } + c_3 Q^{n + 1} \eqno(28)$

其中,$c_1 = \dfrac{1 - \alpha }{\alpha \Delta t}$,$c_2 = \dfrac{1}{\left( {1 - \alpha } \right)\alpha \Delta t}$和$c_2 = \dfrac{2 - \alpha }{\left( {1 - \alpha } \right)\Delta t}$. 由式(28)可知

$\dot{\pmb d }^{n + 1} = c_1 {\pmb d }^n + c_2 {\pmb d }^{n + \alpha } + c_3 {\pmb d }^{n + 1} \hfill (29)$

$ \ddot{\pmb d }^{n + 1} = c_1 \dot{\pmb d }^n + c_2 \dot{\pmb d} ^{n + \alpha } + c_3\dot {\pmb d }^{n + 1} \hfill (30)$

将式(29)和式(30)代入式(27)中,可得

$ \babl \left( {c_3 c_3 {\pmb M} + c_3 {\pmb C} + {\pmb K}} \right){\pmb d}^{n + \alpha } = {\pmb F}_{{\rm EX}}^{n + \alpha } -{\pmb M}\left( {c_1 c_3 {\pmb d}^n + c_2 c_3 {\pmb d}^{n + \alpha } + c_1 \dot{\pmb d }^n +c_2 \dot{\pmb d }^{n + \alpha }} \right) + {\pmb C}\left( {c_1 {\pmb d}^n + c_2 {\pmb d}^{n + \alpha }} \right) \end{array}\eqno(31) $

式中, $\beta = 1/4$, $\gamma = 1/2$和 $\alpha = 1/2$. $n + \alpha $ 时刻外力由$n$和$n +1$时刻值线性插值获得; 弹性体的${\pmb K}^{n + \alpha }$由文献[30]所建议方法获得.

表面上,复合隐式时间积分法的计算时间近两倍于Newmark-$\beta $法,但它能适应更大时间步,故总耗时并未显著延长^[25].

一般而言,结构时间步大于流体时间步,上述特性也将更利于流体子循环^[31].

3 网格更新

这里采用简单高效的子块移动技术^[3]结合正交-半扭转弹簧近似法^[33]更新流体动网格,其基本思想是通过粗糙背 景网格(即子块)运动插值获得精细网格的节点坐标,如图所示. 相比单纯使用弹簧近似法,该手段可大幅提高网格更新效率,具体过程描述如下:

(a) 提取粗、细两层网格信息;

(b) 确定每个子块内的细网格节点数;

(c) 计算每个动网格节点的插值系数;

(d) 由结构运动更新界面节点位置;

(e) 如粗网格有内部节点则调用正交-半扭转弹簧近似法,否则跳过此步;

(f) 对子块位置进行插值获得细网格新位置;

(g) 检查两层网格质量.

同时,为满足几何守恒律,将一质量源项^[34]引入压力泊松方程(7)可得

$\nabla ^2p^{n + 1} = \dfrac{1}{\Delta t}\nabla \cdot {\pmb u}^ * +S_{{\rm MST}}^{n + 1} \eqno(32)$

其中

$S_{{\rm MST}}^{n + 1} = \dfrac{1}{2A_{\rm e}^{n + 1} }\left( \left| \begin{array}{cc} {w_1^2 - w_1^1 } & {w_2^2 - w_2^1 } {w_1^3 - w_1^1 } & {w_2^3 - w_2^1 }\end{array} \right| \right)^{n + 1} \eqno(33) $

式中,$A_{\rm e}$为单元面积,网格速度上标表示局部节点编号,下标表示坐标轴方向. 由式(33)易知质量源项在欧拉网格上退化为零. 引入质量源项可避免构造复杂的网格速度格式^[35].

4 分区强耦合算法

4.1 界面条件

在流固交界面$\Sigma $处须满足速度连续性和应力平衡性,即

${\pmb u} = {\pmb d }\,, \ \ {\pmb t }^{\rm F} = {\pmb t}^{\rm S} \eqno(34)$

其中,流体拖拽力和结构拖拽力写为

$ {\pmb t}^{\rm F} = {\pmb \sigma }^{\rm F} \cdot {\pmb n}^{\rm S}$,${\pmb t}^{\rm S} = {\pmb \sigma }^{\rm S} \cdot {\pmb n}^{\rm S} $ (35)

式中,${\pmb n}^{\rm S}$为界面法向向量,其方向由结构域指向流体域. 此外,还需保障界面处的几何连续性

${\pmb x} = {\pmb d} \,, \ \ {\pmb w} = \dot {\pmb d} \eqno(36)$

上述界面条件可通过构建高阶速度增量和高阶应力增量^{[16-18,36-38]}进行修正,在此不再赘述.

4.2 强耦合算法

固定点法配合Aitken's $\Delta ^{2}$松弛^[13]易实施、精度与效率较为理想,已见诸于各类流固耦合问题模拟.

本文采用此手段耦合各场,主要步骤描述如下:

(a) 预测界面位移^[39]

$\hat{\pmb d }_{\Sigma ,k}^{n + 1} = {\pmb d }_\Sigma ^n + \Delta t\left( {1.5 \dot{\pmb d }_\Sigma ^n - 0.5 \dot{\pmb d }_\Sigma ^{n - 1} } \right) \eqno(37) $

(b) 开始迭代并令$k \leftarrow k + 1$;

(c) 更新流体动网格并计算网格速度

${\pmb w}_k^{n + 1} = \dfrac{\hat {\pmb d}_k^{n + 1} - {\pmb d}^n}{\Delta t} \eqno(38) $

(d) 计算质量源项;

(e) 求解流体NS方程并获得流体激力;

(f) 求解结构运动方程;

(g) 计算界面残差,检查是否收敛, 若收敛则进入下一时间步,否则进行第(h)步;

(h) 重新计算Aitken因子^[9],并松弛界面位移

$\hat{\pmb d }_{\Sigma ,k}^{n + 1} = \lambda _k^{n + 1} {\pmb d}_{\Sigma ,k}^{n + 1} + \left( {1 - \lambda _k^{n + 1} } \right) \hat{\pmb d }_{\Sigma ,k - 1}^{n + 1} \eqno(39)$

(i) 返回至第(b)步.

上述分区强耦合算法的流程图如图3所示.

4.3 Aitken加速技术

在每一迭代步中,动态松弛因子按以下递归公式计算^[13]

$\lambda _k^{n + 1} = \left\{ \!\! \begin{array}{ll} {\max \left( {\lambda _{{\rm MAX}} ,\lambda ^n} \right) }\,, & {k = 1} \lambda _{k - 1}^{n + 1} \dfrac{{ \pmb e\pmb r \pmb r }_k^{\rm T} \left( { {\pmb e\pmb r \pmb r }_{k - 1} - { \pmb e\pmb r \pmb r }_k } \right)}{\left| { {\pmb e\pmb r \pmb r }_k - {\pmb e\pmb r \pmb r}_{k - 1} } \right|} \,, & {k \geqslant 2} \end{array} \!\! \right.$ (40)

式中,$k$代表迭代步,${\pmb e\pmb r \pmb r} $是界面位移残差, $\lambda _{\rm MAX} = 0.1$和 $\lambda_{0} =0.5$.

5 数值算例

5.1 桥梁截面颤振

H型桥梁截面在均匀来流下经历竖向平动和转动,如图4所示. 相关参数设置如下^[40]：流体密度$\rho^{\rm F} =1.25$, 流体黏性$\mu = 0.1$,结构质量$m_{2} = 3\,000$,阻尼因子$c_{2} = 100$,弹簧刚度$k_{2} =2\,000$,转动惯量$I_{\theta } = 25\,300$,转动阻尼因子$c_{\theta } = 2\,200$,转动弹簧刚度$k_{\theta } =40\,000$,入口速度$U = 10$,特征长度$D = 12$和雷诺数${Re} = 1\,500$.

在图4中动网格域限定为A2. 网格划分和子块划分如图5所示,其中计算域包含6\,486个T3单元和3\,329个节点. 时间步长取$\Delta t = 0.02$.

验网格无关性,选用1(6 486个T3单元和3 39个节点)和(11 714个T3单元和5 953个节点)进行对比计算,结果列 于表1. 由表1可知,两组结果相差很小,故本文方法不受网格影响. 综合考虑效率与精度,选取1进行后续模拟.

两方向振幅与振动频率列于表. 由表可知,本文结果稍大于文献[40]结果;桥梁截面经历较大 扭转振动且两方向振 动频率相同.图6所示为两方向位移时程. 由该图可知,本文位移时程曲线相当平滑,而文献[41]的位移曲线则略有起伏,说明本文方法表现稳定. 综合表1和图6可知,结构振动频率非常接近其转动自然频率($f_{\rm N\theta } =0.200\,1$),截面扭转强烈而竖向平动微弱. 此时结构运动由扭转振动控 制,出现明显扭转颤振现象.图7进一步给出某典型时刻涡量场.

5.2 节流阀瓣涡激振动

在均匀流动下,槽道内节流阀瓣发生涡激振动,问题定义如图8所示. 系统参数设定为^[42]：流体密度 $\rho^{F} = 1.0$,流体黏性$\mu = 0.001$,结构密度$\rho^{S} = 1000$ (Case A)和$\rho ^{S} = 62.5$ (Case B),杨氏模量$E = 60\,000$,泊松比$\upsilon = 0.45$,入口速度$U = 1$,立杆高度$D = 1$和雷诺数${Re} = 1000$.

动网格域取为非对称域A2,测点置于立杆左上角,见图8. 流体域划分为2\,786个T3单元和1\,523个节点,固体域 则均匀划分为40 × 2个Q4单元. 流场网格和子块划分如图9所示. 由稳定性条件^[23],将一个Q4单元分为4个光滑子域,形函数构造如图10所示. 时间步长取为$\Delta t = 0.005$.

两组不同流体网格验证算法的网格敏感性,即1( 786个T3单元和1 53个节点)和(5 006个T3单元和 644个节点). 两种工况下测点最大水平位移列于表3. 可见,两组网格结果相差较小,再次验证了本文算法不受网格划分影响. 下文计算均基于1网格.

图11给出两种工况下测点水平位移时程曲线,并与文献[42]进行对比.

解释如下：结构密度较大时,结构在黏性流体中的振动主要受惯性影响,流体作用相当于阻尼振子,逐渐消耗结构动能.

Case A中结构振动的衰减特征正说明了这点. 在Case B中,文献[42]的时程曲线仍可维持相当振幅,而本文结构振动则衰减相当迅速. 文献[43]同样报道了这一现象.

系统参数与数值方法的不同导致了计算结果的差异. 为清楚说明涡激振动现象,图12给出Case A中某典型时刻水平速度与压力云图.

6 结论

本文基于ALE有限元公式提出一种分区强耦合方法,数值模拟了大位移流固耦合问题. 采用半隐式CBS算法求解不可压缩NS方程, 可进行速度--压力等阶插值且省时;运用新颖的复合隐式时间积分法对结构方程进行时间积分,能选取较大时间步. 针对几何非线性,选用精确高效的光滑有限元法计算结构增量方程. 动网格更新采用子块移动技术高效结合正交--半扭转弹簧近似法,同时引入质量源项保证流体分步算法满足几何守恒律. 多场耦合经固定点迭代配以Aitken加速技术,简单高效. 本文强耦合算法能灵活选取不同单场技术,兼具程序模块性、精度与效率优点. 从原理看,较高的效率主要源于流场计算、网格更新和耦合迭代等. 数值算例表明,计算结果与已有数据吻合较好,成功揭示了典型流致振动现象.

The authors have declared that no competing interests exist.

---

参考文献

文献选项

• 原文顺序
• 文献年度倒序
• 文中引用次数倒序
• 被引期刊影响因子

[1]	王英学 , 高波 , 任文强 . 高速铁路隧道缓冲结构气动载荷与结构应力特性分析 . 力学学报, 2017 , 49 ( 1 ): 48 - 54 [本文引用: 2] ( Wang Yingxue , Gao Bo , Ren Wenqiang . Aerodynamic load and structure stress analysis on hood of high-speed railway tunnel . Chinese Journal of Theoretical and Applied Mechanics, 2017 , 49 ( 1 ): 48 - 54 (in Chinese)) [本文引用: 2]
[2]	刘晶波 , 宝鑫 , 谭辉 等 . 波动问题中流体介质的动力人工边界 . 力学学报, 2017 , 49 ( 6 ): 1418 - 1427 ( Liu Jingbo , Bao Xin , Tan Hui , et al . Dynamical artificial boundary for fluid medium in wave motion problems . Chinese Journal of Theoretical and Applied Mechanics, 2017 , 49 ( 6 ): 1418 - 1427 (in Chinese))
[3]	杜特专 , 王一伟 , 黄晨光 等 . 航行体水下发射流固耦合效应分析 . 力学学报, 2017 , 49 ( 4 ): 782 - 792 [本文引用: 1] ( Du Tezhuan , Wang Yiwei , Huang Chenguang , et al . Study on coupling effects of underwater launched vehicle . Chinese Journal of Theoretical and Applied Mechanics, 2017 , 49 ( 4 ): 782 - 792 (in Chinese)) [本文引用: 1]
[4]	陈伟民 , 付一钦 , 郭双喜 等 . 海洋柔性结构涡激振动的流固耦合机理和响应 . 力学进展, 2017 , 47 ( 1 ): 25 - 91 ( Chen Weimin , Fu Yiqin , Guo Shuangxi , et al . Fluid-solid coupling and dynamic response of vortex-induced vibration of slender ocean cylinders . Advances in Mechanics, 2017 , 47 ( 1 ): 25 - 91 (in Chinese))
[5]	He T , Zhang K , Wang T . AC-CBS-based partitioned semi-implicit coupling algorithm for fluid-structure interaction using stabilized second-order pressure scheme . Communications in Computational Physics, 2017 , 21 ( 5 ): 1449 - 1474
[6]	周岱 , 何涛 , 涂佳黄 . 流固耦合分析的一种改进CBS 有限元算法 . 力学学报, 2012 , 44 ( 3 ): 494 - 504 [本文引用: 1] ( Zhou Dai , He Tao , Tu Jiahuang . A modified CBS finite element approach for fluid-structure interaction . Chinese Journal of Theoretical and Applied Mechanics, 2012 , 44 ( 3 ): 494 - 504 (in Chinese)) [本文引用: 1]
[7]	陈威霖 , 及春宁 , 徐万海 . 并列双圆柱流致振动的不对称振动和对称性迟滞研究 . 力学学报, 2015 , 47 ( 5 ): 731 - 739 ( Chen Weilin , Ji Chunning , Xu Wanhai . Numerical investigation on the asymmetric vibration and symmetry hysteresis of flow-induced vibration of two side-by-side cylinders . Chinese Journal of Theoretical and Applied Mechanics, 2015 , 47 ( 5 ): 731 - 739 (in Chinese))
[8]	孙旭 , 张家忠 , 黄必武 . 弹性薄膜类流固耦合问题的CBS有限元分析 . 力学学报, 2013 , 45 ( 5 ): 787 - 791 [本文引用: 1] ( Sun Xu , Zhang Jiazhong , Huang Biwu . An application of the cbs scheme in the fluid-membrane interaction . Chinese Journal of Theoretical and Applied Mechanics, 2013 , 45 ( 5 ): 787 - 791 (in Chinese)) [本文引用: 1]
[9]	刘赵淼 , 南斯琦 , 史艺 . 中等严重程度冠状动脉病变模型的血流动力学参数分析 . 力学学报, 2017 , 49 ( 6 ): 1058 - 1064 [本文引用: 1] ( Liu Zhaomiao , Nan Siqi , Shi Yi . Hemodynamic parameters analysis for coronary artery stenosis of intermediate severity model . Chinese Journal of Theoretical and Applied Mechanics, 2017 , 49 ( 6 ): 1058 - 1064 (in Chinese)) [本文引用: 1]
[10]	周岱 , 何涛 , 刘光众 . 基于改进CIBC方法的方柱驰振预测 . 水动力学研究与进展A辑, 2013 , 28 ( 4 ): 391 - 398 [本文引用: 1] ( Zhou Dai , He Tao , Liu Guangzhong . Numerical analysis of free rotation gallop of a rectangular cylinder using improved CIBC method . Chinese Journal of Hydrodynamics, 2013 , 28 ( 4 ): 391 - 398 (in Chinese)) [本文引用: 1]
[11]	He T . Partitioned coupling strategies for fluid-structure interaction with large displacement: Explicit, implicit and semi-implicit schemes . Wind and Structures, 2015 , 20 ( 3 ): 423 - 448 [本文引用: 1]
[12]	Le Tallec P , Mouro J . Fluid structure interaction with large structural displacements . Computer Methods in Applied Mechanics and Engineering, 2001 , 190 ( 24-25 ): 3039 - 3067
[13]	Küttler U , Wall W . Fixed-point fluid-structure interaction solvers with dynamic relaxation . Computational Mechanics, 2008 , 43 ( 1 ): 61 - 72 [本文引用: 3]
[14]	Bai Y , Sun D , Lin J . Three dimensional numerical simulations of long-span bridge aerodynamics, using block-iterative coupling and DES . Computers& Fluids 2010 , 39 ( 9 ): 1549 - 1561 [本文引用: 1]
[15]	Sun D , Owen JS , Wright NG . Application of the k-w turbulence model for a wind-induced vibration study of 2D bluff bodies . Journal of Wind Engineering and Industrial Aerodynamics 2009 , 97 : 77 - 87 [本文引用: 1]
[16]	He T , Zhang K . Combined interface boundary condition method for fluid-structure interaction: Some improvements and extensions . Ocean Engineering, 2015 , 109 : 243 - 255
[17]	He T . A CBS-based partitioned semi-implicit coupling algorithm for fluid-structure interaction using MCIBC method . Computer Methods in Applied Mechanics and Engineering, 2016 , 298 : 252 - 278 [本文引用: 1]
[18]	He T , Zhang K . An overview of the combined interface boundary condition method for fluid-structure interaction . Archives of Computational Methods in Engineering, 2017 , 24 ( 4 ): 891 - 934
[19]	何涛 . 流固耦合数值方法研究概述与浅析 . 振动与冲击, 2018 , 37 ( 4 ): 184 - 190 [本文引用: 1] ( He Tao . Numerical solution techniques for fluid-structure interaction simulations: A brief review and discussion . Journal of Vibration and Shock, 2018 , 37 ( 4 ): 184 - 190 (in Chinese)) [本文引用: 1]
[20]	Zienkiewicz OC , Nithiarasu P , Codina R , et al . The characteristic-based-split procedure: An efficient and accurate algorithm for fluid problems . International Journal for Numerical Methods in Fluids, 1999 , 31 : 359 - 392
[21]	Normura T , Hughes TJR . An arbitrary Lagrangian-Eulerian finite element method for interaction of fluid and a rigid body . Computer Methods in Applied Mechanics and Engineering, 1992 , 95 ( 1 ): 115 - 138
[22]	Bathe KJ , Ramm E , Wilson EL . Finite element formulations for large deformation dynamic analysis . International Journal for Numerical Methods in Engineering, 1975 , 9 ( 2 ): 353 - 386
[23]	Dai KY , Liu GR . Free and forced vibration analysis using the smoothed finite element method (SFEM) . Journal of Sound and Vibration, 2007 , 85 : 803 - 820 [本文引用: 1]
[24]	Hamrani A , Habib SH , Belaidi I . CS-IGA: A new cell-based smoothed isogeometric analysis for 2D computational mechanics problems . Computer Methods in Applied Mechanics and Engineering, 2017 , 315 : 671 - 690
[25]	Bathe KJ , Noh G. Insight into an implicit time integration scheme for structural dynamics. Computers& Structures , 2012 , 98- 99 : 1 - 6 [本文引用: 1]
[26]	He T . On a partitioned strong coupling algorithm for modeling fluid-structure interaction . International Journal of Applied Mechanics, 2015 , 7 : 1550021
[27]	Zhang J . Accuracy of a composite implicit time integration scheme for structural dynamics . International Journal for Numerical Methods in Engineering, 2017 , 109 ( 3 ): 368 - 406
[28]	Newmark NM . A method of computation for structural dynamics . Journal of Engineering Mechanics-ASCE, 1959 , 85 ( 3 ): 67 - 94
[29]	Collatz L . The Numerical Treatment of Differential Equations. New York: Springer-Verlag, 1960 [本文引用: 1]
[30]	Kuhl D , Crisfield MA . Energy-conserving and decaying algorithms in nonlinear structural dynamics . International Journal for Numerical Methods in Engineering, 1999 , 45 : 569 - 599 [本文引用: 1]
[31]	Piperno S , Farhat C , Larrouturou B . Partitioned procedures for the transient solution of coupled aroelastic problems Part I: Model problem, theory and two-dimensional application . Computer Methods in Applied Mechanics and Engineering, 1995 , 124 ( 1 ): 79 - 112 [本文引用: 1]
[32]	Lefrançois E . A simple mesh deformation technique for fluid-structure interaction based on a submesh approach . International Journal for Numerical Methods in Engineering, 2008 , 75 : 1085 - 1101
[33]	Markou GA , Mouroutis ZS , Charmpis DC , et al . The ortho-semi-torsional (OST) spring analogy method for 3D mesh moving boundary problems . Computer Methods in Applied Mechanics and Engineering, 2007 , 196 : 747 - 765 [本文引用: 1]
[34]	Jan YJ , Sheu TWH . Finite element analysis of vortex shedding oscillations from cylinders in the straight channel . Computational Mechanics, 2004 , 33 : 81 - 94 [本文引用: 1]
[35]	Etienne S , Garon A , Pelletier D . Perspective on the geometric conservation law and finite element methods for ALE simulations of incompressible flow . Journal of Computational Physics, 2009 , 228 : 2313 - 2333 [本文引用: 1]
[36]	He T , Zhou D , Bao Y . Combined interface boundary condition method for fluid-rigid body interaction . Computer Methods in Applied Mechanics and Engineering, 2012 ,223- 224 : 81 - 102
[37]	He T , Zhou D , Han Z , et al . Partitioned subiterative coupling schemes for aeroelasticity using combined interface boundary condition method . International Journal of Computational Fluid Dynamics, 2014 , 28 : 272 - 300
[38]	He T . A partitioned implicit coupling strategy for incompressible flow past an oscillating cylinder . International Journal of Computational Methods, 2015 , 12 ( 2 ): 1550012
[39]	Piperno S . Explicit/implicit fluid/structure staggered procedures with a structural predictor and fluid subcycling for 2D inviscid aeroelastic simulations . International Journal for Numerical Methods in Fluids, 1997 , 25 : 1207 - 1226 [本文引用: 1]
[40]	Filippini G , Dalcin L , Nigro N , et al . Fluid-rigid body interaction by PETSc-FEM driven by Python. Mecánica Computacional, 2008 , XXVII : 489 - 504 [本文引用: 2]
[41]	Dettmer W , Perié D . A computational framework for fluid-rigid body interaction: Finite element formulation and applications . Computer Methods in Applied Mechanics and Engineering, 2006 , 195 ( 13-16 ): 1633 - 1666 [本文引用: 1]
[42]	Olivier M , Dumas G , Morissette JF . A fluid-structure interaction solver for nano-air-vehicle flapping wings//Proceedings of the 19th AIAA Computational Fluid Dynamics, San Antonio, USA, 2009 : 1 - 15 [本文引用: 5]
[43]	Baiges J , Codina R . The fixed-mesh ALE approach applied to solid mechanics and fluid-structure interaction problems . International Journal for Numerical Methods in Engineering, 2010 , 81 : 1529 - 1557 [本文引用: 1]