0.   引　言

置于流体中的结构做变速运动时，会产生附加质量与流体阻尼，进而改变结构的固有频率和阻尼比。研究人员对结构附加质量与流体阻尼的理论分析和试验研究进行了大量的工作^[1-5]。理论分析主要基于势流理论及小振幅假设，其中，势流理论忽略了粘性的影响，因此无法给出流体阻尼，而小振幅假设并不适用于所有情况，特别是对于流致振动这类结构产生大幅振动的情况。随着计算机硬件与数值方法的迅速发展，人们开始采用数值模拟方法研究附加质量与流体阻尼，Yang等^[6]对N-S方程进行线性化，并采用有限元方法求解同心及偏心圆的附加质量与流体阻尼，该方法与基于小振幅假设的理论解相吻合。李刚等^[7-8]采用计算流体动力学（CFD）求解器得到结构在强迫振动下的流体力，将其分解为加速度与速度分量，由此获取附加质量与流体阻尼，并对比有/无粘条件对附加质量的影响。但这些研究并未考虑对于结构振动幅值的影响。

本文基于CFD通用求解器Fluent对同心圆柱结构进行数值模拟，通过用户自定义函数（UDF）设定内圆运动方程，运用嵌套网格技术实现网格运动，采用最小二乘法拟合CFD计算所得的流体力曲线，进而得到附加质量和流体阻尼，最终分析结构振幅及流体粘性对附加质量与流体阻尼的影响。

1.   数学模型

1.1   控制方程

守恒形式的N-S方程^[9]如下：

式中，$ \boldsymbol{\tau} $为粘性应力张量；$ \rho $为密度；$ {\boldsymbol{u}} $为速度矢量；$ t $为时间；$ P $为压力；g为重力加速度。

对式（1）、式（2）进行雷诺时均，其中连续性方程形式不变，动量方程产生额外的雷诺应力$ -\rho \overline{u^{\prime} v^{\prime}} $（u′、$ v^{\prime}$为不同方向上的速度脉动量）。为求解雷诺应力，常用的方法是采用涡粘假设（Boussinesq假设），即将雷诺应力表示为如下形式：

式中，μ_t为湍流粘度，该值由扩散强度所控制；$ U、V $为不同方向上的速度平均量。

本研究使用k-ω 剪切应力输运（SST）湍流模型计算湍流粘度，该模型解决了标准k-ε模型在近壁区域模拟精度不足以及标准$ k {\text{-}}\omega $模型在自由剪切流中性能欠佳的问题。圆柱运动近壁区存在较强剪切层，选用该模型能够较好地捕捉流场中的剪切层结构。

1.2   可压液体

固定密度液体会使壁面运动造成的压力波速度为无穷大，这有可能导致计算发散，因此常采用液体状态方程来描述液体处于平衡状态时压力、体积模量及密度的关系，其中Tait方程就是最常用的半经验液体状态方程^[10]。在Fluent中，Tait方程可简化为如下形式：

式中，$ p_{0} $为参考压力（绝对值）；$ \rho_{0} $为参考密度（$ p_{0} $下的值）；$ K_{0} $为参考体积模量（$ p_{0} $下的值）；$ {n} $为密度指数；$ p $为压力（绝对值）。

1.3   附加质量与流体阻尼计算方法

假定物体在静止的流场中做速度$ V=A \omega \cos (\omega t) $的振荡运动。根据Bearman假设^[1]，物体所受流体力可以表示为$ F \sin (\omega t+\varepsilon) $。将该式展开得到：

式中，$ \dot{x} $为结构速度；$ \ddot{x} $为结构加速度；$ A $为结构振幅；$ \varepsilon $为相位差；$ \omega $为运动角频率；F为流体力振幅。

式（5）右侧第一项与加速度同相，该项为附加质量项，右侧第二项与内圆柱运动反相，为流体阻尼项。定义附加质量m_a为：

流体阻尼c_a为：

利用CFD计算得到流体力曲线，采用最小二乘法拟合以获取流体力振幅和相位差，将其代入式（6）、式（7）即可得到结构的附加质量与流体阻尼。

2.   物理模型

同心圆柱结构可以简化为如图1所示的二维模型，其中内圆直径$ d $为0.0454 m，流体介质为常温水，运动角频率设置为1 rad/s，由此得到运动雷诺数Re，见式（8）。通过改变流体粘性、外圆直径$ D $和结构振幅的取值来获取不同工况下结构的附加质量与流体阻尼，进而开展对比分析。

图  1  同心圆柱几何模型示意图

Figure  1.  Schematic Diagram of Concentric Cylinder Geometry

下载: 全尺寸图片 幻灯片

式中，$ \nu$为运动粘度。

3.   数值计算

3.1   网格划分及无关性验证

本文采用通用的CFD前处理软件ICEM CFD进行网格划分，并采用嵌套网格技术处理网格运动。嵌套网格技术可以在运动过程中维持网格质量，避免传统方法易出现网格质量变差甚至出现负体积的情况，适用于刚体大幅运动的情况。内圆壁面的第一层网格高度取为0.2 mm（此高度保证所有工况下第一层网格的无量纲壁面距离$ y^{+}$小于1），增长率为1.1，并保证插值区域处的网格尺寸相当。图2给出了直径比D/d=2时同心圆柱结构的网格示意图。

图  2  D/d=2时同心圆柱结构网格示意图

Figure  2.  Schematic Diagram of Concentric Cylinder Mesh (D/d=2)

下载: 全尺寸图片 幻灯片

以D/d=2、无量纲振幅A/d=0.1为例开展网格无关性验证，在保持网格总体布置一致的情况下，通过修改前景网格的周向尺寸和背景网格总体尺寸获取不同网格总量下的流体力振幅与相位差，如图3所示。可以看出，当网格总量加密到$ 4.02 \times 10^{4} $后，曲线趋于平缓，继续加密网格对结果影响较小，综合计算资源与精度，选定总量为$ 4.02 \times 10^{4} $的网格开展后续分析。

图  3  流体力振幅与相位差随网格总量变化曲线

Figure  3.  Variation Curves of Fluid Force Amplitude and Phase Difference with the Total Mesh Quantity

下载: 全尺寸图片 幻灯片

3.2   求解器设置

边界条件设置为：内外圆壁面均为无滑移壁面，流体介质为可压液体。

为保证嵌套网格的精度，开启双精度模式，并选择具有算子分裂的压力隐式（PISO）算法求解不可压N-S方程。经过多次尝试，最终选择的时间步长为0.015 s，模拟总时间为20 s。

4.   计算结果及后处理分析

4.1   流场分布

图4给出了D/d=2、A/d=0.1工况下，速度函数在0°、90°、180°和270°相位角时的流场速度云图。根据速度表达式可知，平衡位置处的内圆速度最大，此时其带动着环腔内的流体一起运动，导致流场整体速度最大。

图  4  不同相位角下的速度云图

Figure  4.  Contour of Velocity at Different Phase Angles

下载: 全尺寸图片 幻灯片

4.2   曲线拟合

不同直径比在A/d=0.1下的流体力时程曲线如图5所示，随着直径比的减小，数值计算得到的流体力振荡得越强烈，这是由封闭的壁面边界条件所致。为剔除振荡分量，获取合适的流体力曲线，采用1.3节的理论分析假定拟合的函数为$ F \sin (\omega t+\varepsilon) $，其中流体力振幅和相位差均为未知量，并选取5~20 s的数据作为样本，通过最小二乘法得到拟合的流体力时程曲线。可以看出，拟合曲线与数值计算所得曲线吻合得非常好，拟合优度均大于0.98。

图  5  不同直径比在A/d=0.1下流体力时程曲线

Figure  5.  Fluid Force Curves with Time at A/d=0.1

下载: 全尺寸图片 幻灯片

4.3   无量纲振幅与粘性的影响

图6给出了不同直径比和无量纲振幅条件下的附加质量系数[$ 4m_{\mathrm{a}}/\left(\rho\pi d^2\right) $]和流体阻尼系数[$ 4c\mathrm{_a}/\left(\rho\pi d^2\omega\right) $]曲线，其中，解析解由Chen等^[1]根据线性化的不可压N-S方程[见式（9）、式（10）]得到，该方程基于小振幅假设，忽略了对流项的影响。可以看出，A/d=0.01时的结果与解析解、文献[6]所得结果均相近，说明了数值计算及曲线拟合的准确性。

图  6  附加质量系数与流体阻尼系数随直径比变化图

Figure  6.  Variation of Added Mass Coefficient and Fluid Damping Coefficient with Diameter Ratio

下载: 全尺寸图片 幻灯片

从图6a可以发现，不同无量纲振幅条件下的有/无粘流体的附加质量系数非常接近，即无量纲振幅的大小并不影响附加质量系数；与有粘流体相比，无粘流体的附加质量系数偏低，这是由于粘性的存在会使结构运动“拖拽”更多流体一起做加速运动。图6b给出了有粘条件下的流体阻尼系数（无粘条件下流体阻尼系数接近于0，故未在图中给出），可以发现，有粘条件下的流体阻尼系数随无量纲振幅的增加而增加。

根据文献[11]可知，流体力可以分为压差力和粘性力，流体阻尼也可分为压差阻尼和粘性阻尼。利用UDF分别输出流体力中的压差力和粘性力，并通过曲线拟合最终得到不同无量纲振幅和直径比条件下的压差阻尼系数和粘性阻尼系数，如图7所示。可以看出，相比于压差阻尼系数，粘性阻尼系数在不同无量纲振幅和直径比下的变化较小；同一直径比条件下的总阻尼值随无量纲振幅增加而增加，但是压差阻尼系数占比几乎不变，即两种阻尼等比增减；同一无量纲振幅条件下的压差阻尼系数和粘性阻尼系数均随直径比的减小而增大，此时压差阻尼系数占比也随之增大。

图  7  压差/粘性阻尼系数随直径比和无量纲振幅变化的曲面图

图中标记的数值为压差阻尼系数占比

Figure  7.  Surface Diagram of Viscosity/Pressure Damping Coefficient with Diameter Ratio and Dimensionless Amplitude

下载: 全尺寸图片 幻灯片

4.4   公式修正

Chen^[1]通过线性化的N-S方程得到同心圆柱体运动的流体力的数学表达式，其流体阻尼的解析解为：

式中，M_d为排开介质质量；$\mathrm{Im}(h) $为求解得到的中间变量h的虚部，详细描述见文献[1]。

由4.3节可知，无粘条件下的流体阻尼接近于0，有粘条件下流体阻尼随无量纲振幅变化而变化，即流体阻尼既受粘性影响也受无量纲振幅的影响，然而式（11）并没有考虑无量纲振幅的影响，因此需要对该式进行修正。图8给出了D/d=2、不同无量纲振幅下的流体阻尼系数。

图  8  D/d=2时流体阻尼系数随无量纲振幅变化的曲线

Figure  8.  Variation Curve of Fluid Damping Coefficient with Dimensionless Amplitude (D/d=2)

下载: 全尺寸图片 幻灯片

参考边界层相似理论所得的表达形式^[4]，选取A/d=0.1、0.2和D/d=2的CFD计算值作为样本数据进行拟合，得到流体阻尼修正公式：

式中，$ H(x) $为单位阶跃函数。图9给出了修正公式所得结果与CFD计算结果的对比，两者最大误差不超过5%，说明了修正公式的准确性。需要注意的是，图9并未给出D/d=1.6、A/d=0.3下的流体阻尼系数，这是由于当内圆运动到靠近外圆壁面时，壁面效应影响明显，流体力函数不再近似为一条正弦曲线，即此时的Bearman假设不再成立，附加质量与流体阻尼在单个周期内不再恒定。

图  9  不同无量纲振幅和直径比条件下的修正值与CFD计算值

Figure  9.  Modified Values and CFD Calculated Values under Different Dimensionless Amplitude and Diameter Ratios

下载: 全尺寸图片 幻灯片

5.   结　论

本文基于CFD通用求解器Fluent，通过嵌套网格技术和UDF实现了流场的数值模拟，计算分析了同心圆柱结构的附加质量和流体阻尼，得到如下结论：

（1）基于Bearman假设确定函数形状，并采用最小二乘法对流体力时程曲线进行拟合，得到的结果与解析解非常相近，说明了该方法的准确性。

（2）无粘条件下的流体阻尼系数接近于0，附加质量系数也由于粘性“拖拽”作用的消失而减小。

（3）无量纲振幅对附加质量影响较小，但对流体阻尼影响明显。通过流体力分解发现，压差阻尼和粘性阻尼随无量纲振幅的增减而等比例增减，而直径比的变化会改变两种阻尼的占比。

（4）通过曲线拟合得到了考虑无量纲振幅影响的流体阻尼修正公式，通过对比发现，该修正公式解与CFD计算结果具有较好的一致性。