Galerkin法求解弹性边界条件下圆板的流-固耦合振动特性

陈美霞， 姚仕辉， 谢 坤

(华中科技大学 船舶与海洋工程学院，武汉 430074)

圆板作为常见的结构单元，在工程上有着广泛的应用。在很多情况下，圆板是与流体接触的，如微量泵、蝶阀的圆盘、整体式核反应堆的冷却剂等。流体的存在会显著降低圆板的固有频率，从而改变结构的动力响应特性，因此对圆板和流体的流-固耦合问题进行研究具有重要意义。

Kwak[1]研究了漂浮在无限流域上圆板的振动特性，通过引入无量纲附加质量因子(NAVMI)来评价流体对于圆板的作用，近似求解了圆板的固有频率。后来，Kwak[2]应用Hankel变换计算了嵌在刚性障板中圆板的水动力特性，通过与文献[1]进行比较，发现嵌在刚性障板中的圆板所受到的流体附加质量作用比漂浮在自由液面上的圆板更大。Kwak等[3]求解了漂浮在深度有限的流体上圆板的振动特性，并分析了深度改变对于圆板固有频率的影响。Kwak等通过NAVMI法求解圆板的固有频率时，忽略了流体对圆板振型的影响，因此计算结果有一定误差。Amabili等[4]在考虑流体对结构振型影响的基础上，通过Rayleigh-Ritz法建立控制方程，更精确地求解了结构的固有频率，并与NAVMI法进行对比，结果表明，采用NAVMI法求解低阶频率时误差较小。Bauer[5]研究了弹性薄膜封盖或弹性薄板封盖与容器内流体的耦合振动，结果表明，封盖与流体耦合振动的固有频率会大于没有封盖时流体的晃动频率。Bauer等[6]在Bauer的研究基础上，考虑了流体粘性的影响。Amabili[7]研究了漂浮在容器内流体表面的圆板的振动特性，分析了自由液面表面波对于圆板振动的影响。在文献[7]的基础上，Yousefzadeh等[8]分析了漂浮在自由液面上的功能梯度圆板的振动特性。Cheung等[9]研究了作为开口容器底部结构的弹性圆板的水动力特性，利用叠加原理将自由液面边界条件分离出来，分析了流体自由表面波对圆板振动的影响。Jeong等[10]研究了充满可压缩流体的封闭容器中圆形隔板的振动特性，分析了流体压缩性对于圆板振动特性的影响。Askari等[11]对浸没在开口容器中的圆板的振动特性进行分析，并考虑了液体晃动的影响。Eftekhari[12]引入微分求积法(DQM)求解了圆板的流-固耦合振动问题。Tariverdilo等[13]分别用Fourier-Bessel级数法和变分法分析了圆板的流-固耦合问题，并对比了不同方法的计算结果，但未涉及轴对称模态的计算。Jeong[14]计算了完全相同的两块圆板与流体的耦合振动，分别求解了耦合系统的同相和反相模态。在求解轴对称模态时，Jeong对附加质量矩阵的元素进行了修正，但反相轴对称模态结果的误差仍相对较大。

基于文献[13-14]，本文结合Fourier-Bessel级数展开法和Galerkin法建立控制方程，对与流体接触的边界弹性约束圆板的固有振动特性进行了研究。利用0阶贝塞尔函数的性质，添加附加约束方程，对轴对称模态进行求解。在此基础上分析了弹簧刚度和流体深度对于圆板自由振动的影响。

1 理论分析

1.1 圆板位移的级数表达

如图1所示，圆柱形容器内充满不可压缩的无黏性流体，容器底面和柱面为刚性壁面，顶面覆有弹性圆板。流体密度为ρ0，深度为d；圆板半径为a，板厚为h，圆板边缘设有均匀分布的位移弹簧和转角弹簧，弹簧的分布刚度分别为K和Kφ。

假设圆板满足基尔霍夫薄板理论，则系统的动力平衡方程可表示为

(1)

式中：w，D，ρ分别为板的挠度、弯曲刚度和密度，P为流体压力，▽4为柱坐标系下的双调和算子，其表达式为

(2)

根据模态正交性理论，圆板位移可以展开为真空中振型的级数：

(3)

图1 与流体接触的边界弹性约束圆板Fig.1 Elastically restrained circular plate in contact with fluid

式中：Wnm为圆板在真空中的振型，qnm为广义坐标，n和m分别为节径数和节圆数。Wnm满足式(1)对应的齐次微分方程：

(4)

根据文献[15]可知，圆板在真空中的振型为：

Wnm(r,θ)=[Jn(λnmr)+

CnmIn(λnmr)]cos(nθ)

(5)

式中：Jn和In分别为n阶第一类贝塞尔函数和第一类修正贝塞尔函数。Cnm为与振型有关的系数，λnm为频率系数，满足

(6)

式中：ωnm表示真空中圆板的固有频率。

当圆板材料和尺寸一定时，振型系数Cnm和频率系数λnm取决于弹簧刚度K和Kφ，系数的求解步骤见附录A。由此可知，当圆板边界条件确定时，真空中的振型Wnm为已知函数，只要求解出广义坐标qnm，即可根据式(3)得到圆板在流体中的振型。

1.2 流体速度势函数

假设流体作无旋运动，则可用速度势函数Φ(r,θ,x,t)描述流体运动。结合1.1节中流体不可压缩且无粘性的假设，Φ(r,θ,x,t)满足拉普拉斯方程：

▽2Φ=0

(7)

将速度势函数的时间项和空间项分离，可得[14]

Φ(r,θ,x,t)=iωφ(r,θ,x)eiωt

(8)

式中：φ为空间速度势函数，且同样满足拉普拉斯方程。应用分离变量法可得，对于某一个确定的节径数n，φ的一般解为(参考附录B)：

φ(r,θ,x)=(anx+bn)δn0cosnθ+

Fnscosh(βnsx)]cosnθ

(9)

式中：δn0是Kronecker delta函数，an、bn、Ens、Fns、βns为待求系数，且βns>0。

由于容器柱面和底面为刚性壁面，φ需要满足如下边界条件：

(10)

(11)

将式(9)代入(10)，可以得到关于βns的方程：

(12)

同理，将式(9)代入(11)，可得

a0=0

(13)

Ens=0

(14)

因此，式(9)可以简化为

φ(r,θ,x)=cosnθ[b0δn0+

(15)

1.3 流-固耦合面处的连续性条件

在流-固耦合面，需要满足圆板和流体的速度连续性条件：

(16)

将式(3)、(8)和(15)代入式(16)可得

(17)

其中Rnm(r)=Jn(λnmr)+CnmIn(λnmr)为干圆板振型的径向分量。根据贝塞尔方程本征函数的正交性[16]，在方程(17)左右两边乘以rJn(βnsr)，并在[0,a]积分可得

(18)

当n=0时，常数函数也是贝塞尔方程的本征函数，且对于∀β0s>0，常数函数与J0(β0sr)在[0,a]带权r正交。在方程(17)两边同时乘以r并在[0,a]积分，可得n=0时还需要满足的附加约束方程：

(19)

将式(18)代入(15)，速度势函数φ可以进一步表示为

φ(r,θ,x)=cosnθ[b0δn0+

(20)

其中

(21)

1.4 控制方程求解

假设圆板作小振幅振动，应用线性伯努利方程可得作用于圆板的流体压力：

(22)

此时对于确定的n，平衡方程可以表示为

(23)

应用Galerkin法，在方程(23)两边乘以Wnk(r,θ)并在圆板范围积分，结合真空中振型的正交性[15]可得代数方程：

(24)

分别取m和s的截断数为M和S，可以将式(24)写成矩阵形式。当n>0时，式(24)可表示为

[ρhPn-ω2(ρhZn+ρ0Gn)]qn=0

(25)

其中qn={qn0,qn1,…,qnM}T为系数向量，Pn,Zn,Gn为M+1阶矩阵，各矩阵元素分别为

(26)

(27)

k,m=0,1,2,…,M

(28)

当n=0时，除系数向量q0外b0也是未知量，为求得b0和q0，代数方程组(24)需要和附加约束方程(19)联立求解，表示成矩阵形式分别为

[ρhP0-ω2(ρhZ0+ρ0G0)]q0-ρ0ω2c0b0=0

(29)

l0q0=0

(30)

式中：P0,Z0,G0和q0与n>0的情况有着相同的表达式。c0和l0分别为M+1维的列向量和行向量，其元素分别为

(31)

(32)

此时n=0时的控制方程可以表示为如下矩阵形式：

(33)

式(25)和(33)给出了矩阵形式的控制方程。要使方程存在非平凡解，必须让系数矩阵的行列式为零，由此求解出的ω即为与流体接触的圆板的固有频率；求出对应的非零解即可得广义坐标qnm，代入式(3)可得圆板在流体中的振型。

2 数值计算与结果分析

2.1 方法验证

为验证本文解析方法的正确性，将基于本文解析方法的计算结果与基于有限元法的数值结果进行对比。计算时采用如下参数：圆板半径a=0.2 m，板厚h=0.002 m，密度ρ=2 700 kg/m3，杨氏模量E=69 GPa，泊松比ν=0.3；流体密度ρ0=1 000 kg/m3，深度d=0.02 m；位移弹簧的分布刚度K=103N/m2，转角弹簧的分布刚度Kφ=103N/rad。本文方法采用Matlab实现，截断数M和S均取40。数值结果采用有限元软件Ansys得到，其中圆板采用Shell63单元，流体采用Fluid30单元，网格尺寸取Δ=0.004 m。

对截断数M和S的收敛性分析如表1所示。由表1可知，在保留两位小数时，截断数取40和50的计算结果是一致的，而截断数取30时部分高阶频率尚未收敛，故取截断数M=S=40是合理的。

表1选取不同截断数M和S求得的圆板固有频率

Tab.1Naturalfrequenciesofthecircularplateobtainedbyusingdifferenttruncationnumber

nm固有频率/HzM=S=30M=S=40M=S=5001118.63118.63118.632646.26646.26646.2631 697.141 697.141 697.141011.6111.6011.601305.95305.94305.9421 085.351 085.351 085.3532 409.992 409.992 409.992058.1958.1958.191577.72577.72577.7221 618.651 618.651 618.6533 220.713 220.703 220.7030153.58153.58153.581932.13932.13932.1322 242.582 242.572 242.5734 126.834 126.824 126.82

为验证有限元网格的收敛性，分别对Δ1=0.008 m，Δ2=0.004 m，Δ3=0.002 m三种网格尺寸的有限元模型进行分析，计算结果如表2所示。表中f1、f2、f3为不同模型的固有频率，e1、e2为相对误差，满足e1=(f1-f2)/f2，e2=(f2-f3)/f3。由于有限元网格收敛较慢，在综合考虑计算速度和精度的条件下，本文以相对误差<1%作为有限元网格收敛的评判标准，结合表2的结果，取网格尺寸Δ=Δ2=0.004 m是合理的。

表2采用不同网格尺寸的有限元模型得到的圆板固有频率

Tab.2NaturalfrequenciesofthecircularplateobtainedbyFEMusingdifferentmeshsize

nmΔ1=0.008 mΔ2=0.004 mΔ3=0.002 mf1/Hze1/%f2/Hze2/%f3/Hz01119.230.37118.790.11118.662653.790.88648.100.26646.4131 728.601.441 704.090.421 696.911011.630.1511.610.0411.611308.280.57306.530.17306.0221 101.501.141 089.110.341 085.4532 463.001.722 421.270.512 408.992058.420.2958.250.0958.201583.250.72579.080.22577.8121 644.901.261 624.410.381 618.2833 293.901.803 235.520.553 217.8830154.720.55153.870.15153.631944.291.00934.960.30932.1822 288.701.612 252.420.482 241.7034 244.202.264 150.470.684 122.31注:e1=(f1-f2)/f2,e2=(f2-f3)/f3

分别采用本文方法和有限元法求得的圆板固有频率结果如表3所示。从表中可知，对于n,m=0,1,2,3的模态，本文方法与有限元法计算结果相对误差在1%

表3本文方法与有限元法求得的固有频率结果对比

Tab.3ComparisonofthenaturalfrequenciesobtainedbypresentmethodandFEM

nm有限元法/Hz本文方法/Hz相对误差/%01118.79118.63-0.132648.10646.26-0.2831 704.091 697.14-0.411011.6111.60-0.051306.53305.94-0.1921 089.111 085.35-0.3532 421.272 409.99-0.472058.2558.19-0.111579.08577.72-0.2421 624.411 618.65-0.3533 235.523 220.70-0.4630153.87153.58-0.191934.96932.13-0.3022 252.422 242.57-0.4434 150.474 126.82-0.57

以内，验证了本文方法的正确性。此外，由表3还可以看出随着节径数n和节圆数m的增加，圆板的固有频率变大。需要注意的是，本文假定流体不可压缩，因此流体体积保持不变，故不存在n=m=0的模态。

2.2 圆板与流体的耦合振动特性分析

2.2.1 弹簧刚度对圆板自由振动特性的影响

(a) 自由(K=0, Kφ=0)

(c) 导向(K=0, Kφ=∞)

(d) 固支(K=∞, Kφ=∞)图2 不同边界条件下圆板的无量纲频率Fig.2 Nondimentional frequencies of circular plates with different boundary conditions

流体不仅会降低圆板的固有频率还会改变圆板的振型，图3给出了圆板在真空中和流体中的归一化振型的对比。由图3看出，简支和固支圆板的振型受到流体的影响相对较大，自由边界圆板的振型受到流体影响很小，导向边界圆板在真空中和流体中的振型几乎完全一致。

(a) n=0,m=2自由

(b) n=1,m=2自由

(c) n=0,m=2简支

(d) n=1,m=2简支

(e) n=0,m=2导向

(f) n=1,m=2导向

(g) n=0,m=2固支

(h) n=1,m=2固支图3 圆板在真空中和流体中的归一化振型Fig.3 Normalized mode shapes of the circular plate in vacuo and in fluid

在求解固有频率时，若忽略流体对振型的影响，采用圆板在真空中的振型代替圆板在流体中的振型进行求解，即假定w(r,θ)=Wnm(r,θ)，式(23)可简化为

(34)

当n>0时，式(34)两边乘以Wnm在[0,a]积分可以直接求出固有频率ω，相对于求解矩阵方程(25)的广义特征值问题，将大为减少计算时间。

而对于n=0，式(19)简化为

(35)

此式无法成立，因此在求解轴对称模态(n=0的模态)时，这种简化方法是不适用的。

利用式(34)解出的圆板固有频率如表4所示，表中括号中结果为通过求解式(25)的广义特征值问题得到的固有频率。由表4可知，对于自由和导向边界圆板，忽略流体对振型的影响对圆板固有频率的计算结果影响较小，误差在5%以内；而对于简支和固支圆板，则会产生较大的误差。由图3的结果可知，这是因为自由和导向边界圆板振型受到流体影响更小，采用近似方法计算时假定振型(真空中的振型)与真实振型(流体中的振型)很接近，因此误差较小。从表4还可以看出，对于以上四种边界条件，当周向波数n增加时，近似解法的误差会降低。由此可知，当边界条件为自由、导向或者当所求模态周向波数n较大时，可以忽略流体对振型的影响。

表4 忽略流体对振型影响求得的圆板固有频率Tab.4 Natural frequencies of the circular plate obtained by neglecting the fluid influences on mode shapes

2.2.2 流体深度对圆板自由振动特性的影响

出现上述结果是由于底部刚性壁面会限制流体的垂向运动，且距离刚性底面越近，流体作垂向运动越困难。而圆板受到流体的附加质量取决于流-固耦合层附近流体的运动情况，深度d越小，流-固耦合层附近的流体距离刚性底面越接近，运动越困难，流体随圆板振动时需要圆板提供更大的作用力，反过来说圆板受到流体的作用力也就越大，附加质量也就越大。

(a) n=0

(b) n=1

(c) n=2

(d) n=3图4 流体深度对无量纲频率的影响Fig.4 The influence of fluid depth on nondimentional frequencies

3 总 结

(1) 本文采用Fourier-Bessel级数法和Galerkin法，研究了与流体接触的边界弹性约束圆板的自由振动特性。针对轴对称模态，利用0阶贝塞尔函数的性质，添加附加约束方程进行求解。本文方法与数值法结果吻合良好，验证了本文计算方法的正确性。

(2) 通过改变弹簧刚度研究了几种常见边界条件下圆板的流-固耦合振动特性。结果表明，流体作用会降低圆板的固有频率，且对低阶模态影响更大，此结论对于不同边界条件的圆板均成立。此外，流体还会影响圆板的振型，其中简支和固支边界圆板的振型受流体影响较大，而自由及导向边界圆板的振型受到流体影响较小，针对后两种边界条件，可以用圆板在真空中的振型代替圆板在流体中的振型简化计算，且这种简化具有较高的精度。

(3) 在较低的深度范围内，增大流体深度会显著降低流体对于圆板的附加质量作用；当深度大于1.5倍圆板半径时，流体深度的改变对于圆板自由振动的影响可以忽略。

(36)

(37)

式中：Mr、Nr和Mrθ分别代表单位长度的弯矩、剪力和扭矩，这三种分布内力的表达式为：

(38)

(39)

(40)

式中：ν为圆板的泊松比。将式(38)～(40)和圆板位移表达式(3)代入式(36)和(37)可得：

h1(λnm)+Cnmh2(λnm)=0

(41)

y1(λnm)+Cnmy2(λnm)=0

(42)

式中：h1(λnm),h2(λnm),y1(λnm)和y2(λnm)为关于λnm的函数，其表达式为：

(43)

(44)

(45)

(46)

联立式(41)和(42)并消去Cnm，可得真空中圆板的频率方程：

h1(λnm)y2(λnm)-h2(λnm)y1(λnm)=0

(47)

由式(47)可解得频率系数λnm, 将λnm代入式(41)即可得到振型系数Cnm：

(48)

附录B流场拉普拉斯方程求解

对空间速度势函数φ(r,θ,x)采用分离变量法：

φ(r,θ,x)=R(r)Θ(θ)X(x)

(49)

将式(49)代入拉普拉斯方程，并在方程两边同时除以R(r)×Θ(θ)X(x)，可得

(50)

由于r,θ,x为相互独立的变量，式(50)成立的条件为方程两边等于常数：

(51)

其中μ为常数。可得X(x)满足微分方程：

X″+μX=0

(52)

同时，由式(51)可得

(53)

其中λ为常数。考虑到φ(r,θ,x)在周向以2π为周期，Θ(θ)满足定解条件：

(54)

由方程(54)可求得常数λ和函数Θ(θ)分别为

(55)

在求解自由振动特性时，不妨以振型对称轴作为r轴，则A2=0，可得

Θ(θ)=A1cosnθ

(56)

将式(56)代入式(53)可知R(r)需要满足微分方程：

r2R″+rR′-(μr2+n2)R=0

(57)

针对μ的取值，对以下几种情况进行分类讨论：

1)μ<0，方程(52)的通解为

(58)

此时方程(57)为n阶贝塞尔方程，通解为

(59)

φ(r,θ,x)=Jn(βnsr)cosnθ[Enssinh(βnsx)+

Fnscosh(βnsx)]

(60)

式中：βns>0。

2)μ>0，方程(57)为n阶虚宗量贝塞尔方程，通解为

(61)

由R(0)为有限值知D2=0，同时考虑到柱面刚性壁面条件R′(a)=0得D1=0。因此在μ>0时，φ(r,θ,x)不存在非零解。

3)μ=0，方程(52)的通解为

X(x)=E1x+E2

(62)

此时方程(57)为欧拉型常微分方程，其通解为

(63)

由R(0)为有限值知F2=G2=0，同时考虑到R′(a)=0可得F1=0。因此式μ=0时，速度势函数通解为

φ(r,θ,x)=cosnθ(anx+bn)δn0

(64)

综上可知，对于确定的n，速度势函数的一般解为

φ(r,θ,x)=cosnθ{(anx+bn)δn0+

(65)