相较于传统的金属制螺旋桨，柔性材料螺旋桨因其优越的抗腐蚀、高强度、水弹性变形自适应等物理性能成为当今船舶推进研究的一大热点^[1]。迄今为止，对柔性材料螺旋桨的研究主要集中在使用流固耦合方法对其稳态性能进行预报及优化等方面^[2-6]。Young等^[7]提出了基于螺旋桨面元法与有限元的流固耦合方法，使用该方法研究了一系列柔性复合材料螺旋桨的稳态水动力性能，并针对稳态性能使用遗传算法进行了桨叶结构的优化^[8]。文献[9-11]使用计算流体力学(computational fluid dynamics, CFD)与有限元分析(finite element method, FEM)耦合的方法，针对非均匀来流中的柔性复合材料螺旋桨的水弹性性能进行了优化计算。相对于BEM/FEM(boundary element method, BEM)耦合的方法，CFD/FEM耦合方法计算速度较慢，但对于流场细节的预报更精确。由于实验技术的限制，目前对于柔性螺旋桨的实验仍处于宏观水动力性能测试的方面^[12-13]，对其结构变形及振动的实验测量等问题还有待解决。对于柔性材料螺旋桨动态性能的研究，许多文献都有涉及^[14-15]，但由于实验欠缺等原因仍需进一步研究。

本文着重于采用数值方法研究柔性螺旋桨的非稳态特性。以各向同性材料螺旋桨为例，从时域和频域分析了流场和结构的变化规律和频谱分布，揭示了其内在联系和耦合机理，为柔性材料螺旋桨的设计研究提供了参考。

1 物理模型及网格

本文采用的计算模型是一型柔性材料3D打印加工成型的侧斜桨p1374，由MARINTEK设计、加工和试验，其几何参数与柔性材料主要力学参数如文献[16]。相同几何的对比金属模型为铝制刚性桨，两者几何完全相同，考虑到叶片受力大小，可认为铝制桨为刚性体，弹性变形可忽略。两者实桨模型均为4叶右旋，直径D为0.25 m，如图 1所示。

由于在流固耦合计算过程中，几何模型每次迭代都会变形，使用网格变形算法对网格进行更新，为适应叶片的复杂几何外形并保证网格质量，包裹螺旋桨的旋转域以及旋转域附近的部分静止计算域选用多面体网格，远场静止域部分采用拉伸网格。所有固壁面采用统一尺度的边界层网格，边界层设置10层网格，第1层网格厚度根据工作条件设置小于0.01 mm，以保证y+ < 5。叶片梢部、导边等流动变化复杂的区域也均做了加密处理。为尽可能消除远场边界对计算结果的影响，参考Morgut给出的设置^[17]，静止域速度入口到桨盘面距离取16D，静压出口距桨盘32D，静止域自由滑移圆周边界取24D。计算域设置及网格见图 2。

2 数值方法 2.1 流场控制方程

流场计算控制方程为雷诺平均Navier-Stokes(reynolds average navier-stokes, RANS)方程，考虑流动介质为水，其不可压形式为：

$ \frac{\partial \bar{u}}{\partial {{x}_{i}}}=0 $

(1)

$\frac{\partial {{{\bar{u}}}_{i}}}{\partial t}+\frac{\partial {{{\bar{u}}}_{i}}\partial {{{\bar{u}}}_{j}}}{\partial {{x}_{j}}}=-\frac{1}{\rho }\frac{\partial \bar{p}}{\partial {{x}_{i}}}+\text{ }v\frac{{{\partial }^{2}}{{{\bar{u}}}_{i}}}{\partial {{x}_{j}}\partial {{x}_{j}}}+\frac{\partial {{{\bar{\tau }}}_{ij}}}{\partial {{x}_{j}}} $

(2)

式中ν和ρ分别为流体运动粘力系数和密度。应用多参考系模型(multiple reference frame, MRF)来求解旋转计算域稳态流场，在两域交界面做通量传递。非定常流场求解模型采用滑移网格方法，相比于多参考系模型其旋转域网格以ω的角速度在每一步都进行更新，其求解结果更精确可靠。湍流模型选用剪力输运k-ω模型(shear stress transport, SST k-ω)，其固壁面边界层处理方式为增强壁面模型，相较其余湍流模型具有较高的精度，在工程应用中比较广泛^[18]。

2.2 结构计算数值模型

有限元计算模块控制方程为：

$ {\mathit{\boldsymbol{M}}_{\rm{s}}}\left\{ {\ddot u\left( t \right)} \right\} + {\rm{ }}{\mathit{\boldsymbol{C}}_{\rm{s}}}\left\{ {\dot u\left( t \right)} \right\} + {\mathit{\boldsymbol{K}}_s}\left\{ {u\left( t \right)} \right\} = {\rm{ }}{\mathit{\boldsymbol{F}}_{{\rm{ext}}}}\left( t \right) $

(3)

式中：M_s、C_s和K_s分别为质量矩阵、阻尼矩阵和刚度矩阵；F_ext为结构负载。在螺旋桨流固耦合计算中，结构负载主要包括：离心力、科氏力和所受的流体作用力：

$ {\mathit{\boldsymbol{F}}_{{\rm{ext}}}}\left( t \right) = {\mathit{\boldsymbol{F}}_{{\rm{cori}}}} + {\rm{ }}{\mathit{\boldsymbol{F}}_{{\rm{cent}}}} + {\mathit{\boldsymbol{F}}_{\rm{f}}}\left( t \right) $

(4)

式中：流体作用力F_f(t)可以写做水动力压力在流固耦合表面法向上的积分

$ {\mathit{\boldsymbol{F}}_{\rm{f}}} = \int\limits_s {} {\mathit{\boldsymbol{N}}^{\rm{T}}}\left\{ P \right\}{\rm{d}}S $

(5)

式中：结构表面的F_f为流体作用力；N为流固耦合表面的法向量；P为流体作用在物体表面的压力；S为流固耦合表面积。

2.3 流固耦合设置

在双向流固耦合计算中，可以将螺旋桨叶片看做大变形柔性体，流体计算和结构计算使用隐式迭代，每次内迭代过程中交换数据，并采用局部网格重构法生成新的流体网格。考虑到网格变形和滑移运动，采用任意拉格朗日欧拉方法(arbitrary Lagrangian-Eulerian, ALE)来解决计算过程中流固耦合边界上节点不对应的问题。

3 计算结果及分析 3.1 非定常流场结果及分析

本文应用滑移网格模型求解URANS(unsteady reynolds average navier-stokes, URANS)控制方程，从而得到螺旋桨瞬时流场压力分布，采用隐式方法耦合，在每次内迭代完成时，使用ALE方法将其传递给固体结构做有限元变形计算，得到新的模型进而更新流体区域，然后进行下一次内迭代，重复进行直至该时间步收敛。

螺旋桨转速固定为7 rad/s，耦合时间步长取5×10^-4 s，内迭代最大步数为10，为保证计算结果最终收敛仿真总时间设置为5 s，并且使用稳态计算结果对流场进行初始化。图 3给出了6个不同时刻的r/R=0.9处桨剖面压力分布展开图，来流V=0.86 m/s方向为-X。可以看出，在不同时刻下，叶片转动到不同位置，其桨叶剖面展开图中压力分布基本保持一致。最大和最小压力分别出现在压力面和吸力面的导边处，不同时刻压力极值有所变化，但改变幅度不明显。随边处低压区域是由于螺旋桨梢涡和随边泻出涡共同影响，在梢涡及随边泻出涡生成、发展脱落直到耗散的过程中，随边尾流区域压力随之改变，从而产生了一个作用于叶片结构的非定常负载。

图 4给出了桨叶r/R=0.9截面处，螺旋桨旋转一周内的压力极值的波动曲线，其中压力最大值出现在压力面导边处，最小值在导边吸力面处。由于流动的不均匀性，导边附着涡在压力面和吸力面的生成及耗散的频率不同，从而造成了压力面和吸力面的压力波动曲线形式不同。另外可以看出，同一周期内，吸力面压力最小值波动曲线近似呈正弦变化，而压力面压力最大值与桨叶位置关系不规律。其中一个原因是叶片压力面导边处的附着涡结构与吸力面不同；另一个原因是由于桨叶之间的互相干扰，从图 3可以看出，前一个叶片的梢涡形成的诱导涡会对压力面的压力分布产生影响^[19]，从而非定常的涡耗散会与原本涡结构的波动叠加，因此产生了图 4中不同频率压力波动叠加的结果。

桨叶的非定常压力分布可以看做，定常压力分布与非定常波动量的叠加，由于流动速度以及涡结构的不同，叠加后叶片不同区域的压力波动频率也有所差异。而螺旋桨叶片总推力为叶片压力负载在叶片表面上的积分，因此螺旋桨叶片非定常总推力是叶片所有区域的不同频率压力波动的叠加。压力分布在每个时刻都传递给叶片结构，相当于在结构上作用了一个非定常变化的压力分布负载，而当结构改变的时候，流场的非定常波动也随之改变，包括流动速度和涡结构，互相作用直至非定常流动和结构振动共同收敛。

3.2 结构变形及应力分析

结构边界条件为桨叶叶根具有周向旋转自由度无变形，叶片表面为流固耦合面。图 5给出了t=5 s时刻桨叶x和y方向变形分布，可以看出，最大弯曲变形出现在负载最大、厚度最小的叶梢处，同一半径处的变形基本呈均匀分布，考虑到叶片结构采用的材料为各向同性柔性材料，这说明叶片变形主要为向+x的弯曲变形，其扭曲变形不大。而yz方向变形极值出现在叶梢部和导边处，这是因为叶片迎流方向存在流体阻力，从而形成叶片周向负载，因此叶梢部会向后变形，从而拉扯导边，因此yz方向极值在叶片梢部和导边处。

图 6给出了叶片等效应力分布云图，可以看到，由于叶根约束的存在，最大应力出现在叶片根部中间部分。由于叶片流动阻力的存在，其结构负载存在周向分量，因此叶片导边和随边应力较大，而叶梢自由无约束因此最小应力出现在叶梢部分。这说明叶片轴向推力主要影响着叶片根部应力，而导边部分和随边部分主要受周向流体阻力的影响。因此，在进行柔性桨叶片强度校核的时候，应该分2种情况分别进行考虑，从而确定其强度要求及破坏形式。

3.3 流场脉动及结构模态分析

图 7(a)给出了0.5 s内螺旋桨推力波动曲线，曲线采样频率为2 000 Hz。可以看出，螺旋桨推力波动范围在68.0~72.0 N，平均推力为70.5 N。包络曲线波动稳定表明非定常计算收敛，低频包络频率为螺旋桨旋转周期，波峰波谷区域存在由非定常压力波动和结构振动造成的高频波动。为了更进一步分析其动态波动特性，使用快速傅里叶转换(fast Fourier transform, FFT)的方法将时域响应转换到频域，其0波幅为70.27 N，图 7(b)给出其5~1 kHz波谱。

在柔性螺旋桨非定常推力的频域波谱中，0波幅可以看做是受静力负载变形后的刚性螺旋桨在定常流稳态条件下的推力响应，而高阶频率分量则是由流场非定常波动引起的。可以看出，非定常波动的频率分量主要处于40~500 Hz，其中156.25 Hz是低频主要分量，波幅最高为1.389 N。对比柔性螺旋桨的非定常推力与变形后刚性螺旋桨定常计算结果，相同工作条件下，定常计算推力为70.7 N，与柔性桨的直流分量相差很小，而非定常推力高频分量的贡献最大值在1.389 N与变形后刚性螺旋桨推力相比波动幅度不大。因此，在均匀来流下，各向同性柔性材料螺旋桨性能的计算，使用定常方法替代非定常方法得到的推力和扭矩结果是可以接受的。

为进一步研究柔性桨叶片在非定常流场中的振动特性，使用有限元与边界元声学耦合的方法计算了螺旋桨桨叶的干湿模态，其5阶内干湿模态频率见表 1。可以看出，螺旋桨叶片干湿模态频率具有很大的变化，这是因为在流体中，由于流体阻尼的存在从而降低了叶片的自然频率。因此，在设计过程中考虑共振等因素时，需要考虑螺旋桨的工作环境导致的模态频率的改变。由于模态频率的改变，叶片在受到流体非定常负载时会产生振动变形，图 8(a)给出了叶片r/R=0.9处螺距角在1 s内的变化曲线，并使用FFT方法将其转换成频域响应，曲线如图 8(b)所示。

其中叶片在0.9R处的螺距角频域直流分量(0)为22.58°，而初始几何该截面的螺距角为20.15°，这说明柔性螺旋桨叶片在受到静力负载时0.9R处的螺距角变形为两者之差2.43°。而从频谱响应曲线可以看出，贡献最大的2个频率分别为70.1 Hz和140.6 Hz，2者分别与螺旋桨叶片的前两阶湿模态频率接近，这说明在流体非定常负载的作用下，其中处于这2者频率下的流体负载分量与螺旋桨叶片产生了共振，从而导致2者频率下的变形较明显。与之有明显对比的是贡献相对较少的156.3 Hz频率，这也是图 7中所显示的流体非定常作用力。这说明，在受到流体非定常负载的作用下，影响结构变形的主要因素为结构的前2阶湿模态，其次为流体非定常负载的频率分量。

4 结论

1) 影响桨叶表面压力波动的主要因素在于涡结构的生成与耗散，涡结构受螺旋桨叶片几何和工作条件的影响并与结构变形及振动互相干扰。

2) 各向同性柔性材料桨叶的变形主要体现在推力方向的弯曲变形。

3) 桨叶应力分布与桨叶所受轴向推力和周向的流体阻力有关，周向流体阻力决定了桨叶导边和随边的应力大小。

4) 非定常流场的负载波动具有多个频段信息，而桨叶振动频率主要受前2阶湿模态频率的影响，非定常流场负载频率的影响次于模态频率影响。

综上所述，柔性桨流场的非定常特性与结构振动密切相关，而对与各向同性材料的螺旋桨在只考虑推进性能的计算时，准静态的流固耦合方法是合适的。在关心螺旋桨振动特性和非定常流场以及声场的情况下，非定常流场计算是很必要的。由于涡结构是影响柔性桨的振动特性的一大因素，因此在后续工作中可以考虑大涡模拟流固耦合的方法开展更进一步研究。