图1 变截面柔性矩形薄板模型
Dynamic characteristic analysis of flexible rectangular thin plate with variable cross section in airflow
柔性体-流体耦合运动是指柔性体的抗弯刚度一般很小，流体不仅会对柔性体产生作用力影响其运动，使之产生大位移和变形，而且柔性体的变形反过来又会影响流场结构，二者相互作用、相互影响.因此，柔性体与运动流体的耦合运动不仅与来流的速度和方向有关，还受到柔性体的密度和抗弯刚度等材料参数的影响.浸没在稳定均匀流中的弹性柔性薄板的稳定性因其在许多工程实际应用的相关性而受到广泛关注：研究自然界中鱼在水中的游动［1］，植物在风中重构和伏倒［2-3］，鸟类昆虫的飞行和小型潜水器的优化［4-5］，防止薄膜产品在加工时因气流发生颤振因而损坏［6］，治疗呼吸障碍以及分析大动脉中的动脉瘤和人工心脏瓣膜［7-8］；在实际工程中，远距离特高压输电电线和海底运气运油等管道，均要考虑柔性体与流体的耦合作用，防止发生颤振对工程产生危害［9］.
实验方面，Watanabe等［10］在风洞实验中得到了刚度、质量比和张力对柔性体颤振速度的影响；Zhang等［11］在流动肥皂膜中用2条柔性细丝为模型，研究了细丝之间的流体动力学耦合作用，显示柔性体不仅可以在流体中振颤，也会有绷直静止状态；Jia等［12］在实验中用不同柔性体材料，研究了柔性体在风中的运动，证实风速对颤振频率和波数的影响，并用小扰动分析理论对其推导，实验数据和数值模拟对理论结果得到了很好的验证；Kim等［13］在高雷诺数开环风洞中探索了前端自由和后端夹紧柔性体的自激动力学，证实无量纲弯曲刚度、流速和长宽比对柔性体自激振荡的影响；Liao等［14］给出不同边界条件下矩形板与流场耦合的压力场和速度，说明流场对流-固耦合系统动力学行为的影响；Liu和Lu［15］在实验中用不同厚度的板测试了变截面柔性矩形薄板在空气中的自然振动特性，并利用柱的附加质量理论分析了变截面柔性矩形薄板在静态水中的自然振动特性；Shoele 和 Mittal［16］利用非线性特征值分析来研究通道之间柔性板的稳定性，计算了不同约束条件和放置方式下板的临界速度随板的质量比的变化关系；Wang等［17］和王文江［18］研究了静止流体中不同形状和不均匀弦向刚度分布的柔性板的颤振运动，显示柔性体的周围流场环境、自身构型及材料力学性质等因素对其动力推进特性的影响.
在自然界和工程领域中，为了提高材料的利用率，优化结构的重量，更符合自然界模型等需要，因此，在很多情况下要根据外部荷载和边界条件让板的厚度作各种相应变化.本文基于理想流体的Euler方程和 Euler-Bernoulli梁理论［12］，建立了在稳定气流中变截面柔性矩形薄板的振动模型.考虑变截面柔性板的厚度只沿垂直于几何中面方向呈线性分布，用线性小扰动分析法研究了气流中变截面柔性矩形薄板的动力特性.
在Euler-Bernoulli梁模型的基础上，建立两端固定变截面柔性矩形薄板的几何模型和稳定气流中变截面柔性矩形薄板颤振运动模型如图1和2所示.变截面柔性矩形薄板的两端固定在O和A两点（假设OA之间距离足够长，可以减少颤振时两端固定对于变截面柔性矩形薄板起始端和末尾端的影响），且以O点为原点，OA方向为x轴建立坐标系.h0和λ分别为该变截面柔性矩形薄板的初始厚度和截面变化系数，令λ=tanφ，则变截面柔性矩形薄板某一横向位置处的截面厚度h(x)=h0+ λx，稳定气流在变截面柔性矩形薄板上下表面均匀流过，其初始流动速度为U.
图1 变截面柔性矩形薄板模型
图2 变截面柔性矩形薄板在气流中的运动模型
两端固定的变截面柔性矩形薄板服从线性化的Euler-Bernoulli梁理论，轴向气流假设为非粘性不可压缩的理想流体，则由动力学方程可得变截面柔性矩形薄板在稳定气流中颤振的平衡微分方程为
式中：E为变截面柔性矩形薄板的弹性模量，I(x)为横截面惯性矩，ρ1为变截面柔性矩形薄板的线密度，B(x)为横截面面积，w(x ， t)为变截面柔性矩形薄板在稳定气流中颤振的挠度，Δp为变截面柔性矩形薄板上下表面压强的差值（上下表面压强不一致，因而造成变截面柔性矩形薄板在气流中颤振）.
稳定匀速气流的流速（U）为非常低的马赫数.因此，忽略外层的粘性应力并假设不可压缩性，并且流场的质量和动量守恒，由欧拉方程表示为：
式中：u为变截面柔性矩形薄板上颤振时的气流速度，∇为哈密顿算子，ρ2为气流的密度，p为颤振时气流的压强场.
气流在变截面柔性矩形薄板上具有不可滑移和不可渗透性，满足流-固耦合的速度相容条件，则变截面柔性矩形薄板上下表面的颤振速度和气流运动的速度一致，即
式中n是变截面柔性矩形薄板在任意一点垂直于板面的法向量.
假定两端之间距离足够长，则颤振时可以减少流-固耦合运动时固定端对板起始和末尾端的影响.变截面柔性矩形薄板在稳定气流时的法向模态为正弦模态，其复频率为Ω.因此，为了寻找线性化形式问题的解，令Ω=ω+iσ，其中实部（ω）为变截面柔性矩形薄板运动时的频率，虚部（σ）为波扰动时间增长率.当变截面柔性矩形薄板因气流发生颤振时，用线性小扰动分析理论给出各变量的值［19］：
式中：p为变截面柔性矩形薄板颤振运动时流场内的动压强，p0为变截面柔性矩形薄板在未发生颤振时流场内的静压强，U为流场中气流与板接触时的流速，k为变截面柔性矩形薄板颤振时的波数，D为变截面柔性矩形薄板颤振时的振幅，κ为颤振时沿着z方向的空间衰减率，ux和uz分别为变截面柔性矩形薄板颤振运动时气流在x和z方向上的速度.其中P′、U1和W是利用小扰动分析理论加入的小变化量，ρ2U2为压强项中的量纲参考值，W和U1为速度项中的量纲参考值.4个参变量w、ux、uz和p的小扰动线性化都是将表达式中的值分为2种不同状态下的和，即静止状态对应的零阶项与颤振运动时一阶项的和.
将式（5）～（8）代入到式（2）～（4）中，求解未知量得：
在式（11）中，因变截面柔性矩形薄板的颤振运动对气流压强的影响随着z坐标值的增大而不断地削弱.因此，本文设定：当z≥ 0时，k= κ；当z＜ 0时，k=- κ.式（10）和（12）中的正负号表示所在板面的z轴上下区域.
对于变截面柔性矩形薄板在气流中颤振主要因素——Δp的值，把式（10）代入式（5）得到
将式（8）和（13）代入变截面柔性矩形薄板的平衡微分方程（1）中，得到线性小扰动分析理论变换后的控制方程为
将变截面柔性矩形薄板的截面厚度h(x)=h0+λx代入方程（14）中，并将其无量纲化，引入无量纲变量
式中：
为无量纲密度比，
为无量纲复频率
为无量纲颤振波数，
为无量纲横向位置，
为无量纲弹性模量.将无量纲量代入方程（14），得到变截面柔性矩形波动板的无量纲平衡方程
在变截面柔性矩形薄板与流体耦合运动的研究中，文中设
=x/h0，柔性结构的厚度（h0）相对于长度是一个量级比较小的值，当
=1 000时，λ=0.001、0.003和0.005，柔性体两端厚度相比较已经有明显的变化，对结构的整体刚度已有不小的影响；当λ取较大值时，结构整体刚度值会比较大（柔性体的刚度值是个小值），且本文讨论的是柔性薄板，故取 3个典型值 .在方程（16）第 1项中的
，当 λ=0.001、0.003 和 0.005，这项的值相比较
是一个特别小的值，并且带有虚数i，表明此项只表示变化趋势.因此，将此项忽略，则式（16）可写为
式（17）得到一个关于
的二次方程式，求解得
令式（18）根号下的值为0，得到一个关于
的表达式，把此时
的值称为截至波数
，则有
当变截面柔性矩形薄板颤振的
时，由表达式（18）可知，此时表达式中
的值恒为正，表示变截面柔性矩形薄板颤振是不增长也不衰减的稳定的流-固耦合颤振运动；当变截面柔性矩形薄板颤振的
时，由
和式（18）可以得到无量纲频率（
）和无量纲扰动时间增长率（
）的表达式分别为：
为了研究不同
和k对变截面柔性矩形薄板颤振
的影响，本文通过算例展开讨论.取壁板的材料特性参数 
=1000，得到
和
对变截面柔性矩形薄板颤振
的影响如图3所示.当λ=0时，变截面柔性矩形薄板退化为等截面柔性矩形薄板，且与Jia等［12］的无量纲频率一致；当变截面柔性矩形薄板的
时，变截面柔性矩形薄板的
随
的增大而增大，并且在增大到一定的值后会慢慢趋于一个定值，表明变截面柔性矩形薄板在气流中当
到达一定值后，流-固耦合运动会趋于稳定；当
确定时，变截面柔性矩形薄板的
随着λ的增大而减小，这是由于变截面系数增大也使柔性体整体的等效刚度和等效密度增大，从而使
变小；当λ一定时，变截面柔性矩形薄板的
随着
的增大而减小，这也是由于柔性体整体的等效刚度和等效密度增大，从而使
变小 .
图3 不同截面变化系数（λ）下无量纲波数（
）和无量纲横向位置（
）对无量纲频率（
）的影响（a）λ=0；（b）λ=0.001；（c）λ=0.003；（d）λ=0.005
在
=10 和 1000 时，
和 λ对变截面柔性矩形薄板
的影响曲线如图4所示.等截面柔性矩形薄板的
与
无关，只 随着
的 增 长 而 增 大 .变 截 面柔性矩形薄板的
随着
的变大，其影响趋势也在增加，即
愈大，
衰减的愈大 .当变截面柔性矩形薄板颤振的
时，即不稳定范围内
与
无关.
图4 不同无量纲横向位置（
）下无量纲波数（
）和截面变化系数（λ）对无量纲频率（
）的影响（a）
=10；（b）
=1000
、
和 λ对
的影响如图 5所示 .由式（21）可知，当 =0时，变截面柔性矩形薄板的
与 λ无关，所以在图5（b）中，当 λ=0、0.001、0.003和0.005时，4条
与
的关系曲线趋势都是一致的 .由图 5（a）可知，当=0 时，变截面柔性矩形薄板的
随着
的增大而增大，且增大到极大值后会立刻减小为0；λ一定时，
随着
的增大而减小，则可以得到变截面柔性矩形薄板在流-固耦合的颤振运动中抗弯刚度越大，越容易趋于稳定.
图5 无量纲波数（
）、无量纲弹性模量（
）和截面变化系数（λ）对无量纲扰动时间增长率（
）的影响（a）λ=0.001；（b）
=1000
在 
=1000 时，λ、
和
对变截面柔性矩形薄板的影响如图6所示.当λ=0时，变截面柔性矩形薄板退化为等截面柔性矩形薄板，与Jia等［12］的
一致；当变截面柔性矩形薄板颤振的
，λ一定时，随着
的增长，
和
在减小；
一定时，随着
的增长，
在不断的增加，到达极大值之后会立刻衰减为0，即此时变截面柔性矩形薄板颤振趋于稳定.
图6 截面变化系数（λ）、无量纲波数（
）和无量纲弹性模量（
）对无量纲扰动时间增长率（
）的影响（a）λ=0；（b）λ=0.001；（c）λ=0.003；（d）λ=0.005
、
和 λ对
的影响曲线如图 7 所示 .当变截面柔性矩形薄板颤振的
，同一变截面柔性矩形薄板的
在不同
的增长率是不同的，当
=10-4时
=0的
比
=1000 的值大，但随着
的增加，
=1000下
的增加趋势明显比
=0 要大，因而在
=8×10-4附近有一个交点，即2个不同
有一个相同的扰动，但不一定同时发生；当
和
一定时，变截面柔性板
的极大值随着λ的增大而减小，即λ越大（或变截面柔性板的尾部），变截面柔性板越先趋于稳定.
图7 无量纲波数（
）、无量纲横向位置（
）和截面变化系数（λ）对无量纲扰动时间增长率（
）的影响曲线（a）λ=0.003；（b）
=1000，
=1000
变截面柔性矩形薄板颤振的
时，
随着
的增大得到一个极大值，设
对应的
为最不稳定波数，即
时的波数记为
，由式（21）可得
由于研究的材料为柔性体薄板，其抗弯刚度是一个量级特别小的值.因此，
也是一个小值，令
→ 0时，根据式（22）可以求得
的近似值为
将式（23）分别代入式（20）和（21）中，得到
所对应的
和
分别为：
计算得到 
所对应的
近似表达式 .现在将
和
引入物理量纲有：
由式（26）可知，变截面柔性矩形薄板在稳定恒速气流中发生颤振后，最大频率等于动压强项ρ2U与惯性项ρ1(h0+λx)的比值，最大频率与气流流速成正比，与柔性体密度和柔性体厚度均成反比.变截面柔性矩形薄板在稳定恒速气流中发生颤振后，
等于动压强项 ρ2U与刚度项E(h 0+ λx)3的比值 .
与U成正比，与柔性体弹性模量和柔性体厚度均成反比 .由式（26）和（27）可以得到 ωm/km的值与气流流速
也是正比.
在本文的分析中，由于研究的是弹性弯曲，柔性体的重力显然是次主导因素.但由于重力会使变截面柔性矩形薄板发生微小的形变从而打破流-固耦合颤振时的对称平衡状态.设变量速度（Ug）为重力影响的柔性体运动速度比例，由动压强
与单位表面的薄板质量
，得到重力变量速度（Ug）为
由式（28）可知，Ug是一个由动压强和单位表面的薄板质量决定的值.把 Ug代入式（26）和（27）中，可以消去由重力对最不稳定波数和频率的影响，最终得到精确的最不稳定波数和最大频率值.
本文建立了在稳定气流中变截面柔性矩形薄板的颤振模型，考虑变截面柔性板的厚度只沿垂直于几何中面方向呈线性分布，用线性小扰动分析法研究了气流中变截面柔性矩形薄板的动力特性.研究了无量纲波数、薄板的横向位置和不同变截面系数对于变截面柔性矩形薄板的无量纲角频率、无量纲扰动时间增长率的影响，并分析了气流速度对变截面柔性板波数和频率的影响，得到了最不稳定颤振时的波数和频率值，得到如下结论：
（1）变截面柔性矩形薄板的随的增大而增大，并且在增大到一定的值后会慢慢趋于一个定值；变截面柔性矩形薄板的随着λ和的增大而减小.
（2）变截面柔性矩形薄板的随着的增大而增大，且增大到极大值后会立刻减小为0；随着-E的增大而减小，则变截面柔性矩形薄板在流-固耦合的颤振运动中抗弯刚度越大越容易趋于稳定，即在受到扰动后刚度越大的变截面柔性板（或变截面柔性板的尾部）越容易趋于稳定.
（3）变截面柔性矩形薄板在稳定恒速气流中发生颤振后，最大频率等于动压强与惯性项的比值.最大频率与气流流速成正比，与柔性体密度和柔性体厚度均成反比；最不稳定波数等于动压强项与刚度项的比值，最不稳定波数与气流流速成正比，与柔性体弹性模量和柔性体厚度均成反比.
［1］LIAO J C，BEAL D N，LAUDER G V，et al.Fish exploiting vortices decrease muscle activity［J］.Science，2003，302（5650）：1566-1569.
［2］GOSSELIN F，LANGRE E D，MACHADO-ALMEIDA B A.Drag reduction of flexible plates by reconfiguration［J］.Journal of Fluid Mechanics，2010，650：319-341.
［3］ CISONNI J，LUCEY A D，ELLIOTT N S J.Flutter of structurally inhomogeneous cantilevers in laminar channel flow［J］.Journal of Fluids Structures，2019，90：177-189.
［4］王思莹.柔性体与流体耦合运动的数值模拟和实验研究［D］.合肥：中国科学技术大学，2010.
［5］彭泽瑞.仿生自主推进柔性板集群运动的流固耦合数值研究［D］.合肥：中国科学技术大学，2018.
［6］魏恒.均匀来流中柔性体流固耦合数值模拟［D］.合肥：中国科学技术大学，2018.
［7］HUANG L X.Flutter of cantilevered plates in axial flow［J］.Journal of Fluids and Structures，1995，9（2）：127-147.
［8］ PIRNAR J，DOLENC-GROŠELJ L，FAJDIGA I，et al.Computational fluid-structure interaction simulation of airflow in the human upper airway［J］.Journal of Biomechanics，2015，48（13）：3685-3691.
［9］张成瑶.单个或多个柔性体自主推进的流固耦合数值研究［D］.合肥：中国科学技术大学，2020.
［10］WATANABE Y，SUZUKI S，SUGIHARA M，et al.An experimental study of paper flutter［J］.Journal of Fluids and Structures，2002，16（4）：529-542.
［11］ ZHANG J，CHILDRESS S，LIBCHABER A，et al.Flexible filaments in a flowing soap film as a model for one-dimensional flags in a two-dimensional wind［J］.Nature，2000，408（6814）：835-839.
［12］JIA P，ANDREOTTI B，CLAUDIN P.Paper waves in the wind［J］.Physics of Fluids，2015，27（10）：27-63.
［13］KIM D，COSSÉ J，CERDEIRA C H，et al.Flapping dynamics of an inverted flag［J］.Journal of Fluid Mechanics，2013，736：R1-R12.
［14］LIAO C Y，CHEN G W，HSU H W，et al.Theoretical analysis of vibration characteristics of rectangular thin plate fully immersed in fluid with finite dimension［J］.International Journal of Mechanical Sciences，2021，189：105979.
［15］ LIU L F，LU D G.Research on mechanism of the largeamplitude and narrow-band vibration of a flexible flat plate in the rectangular channel［J］.Nuclear Engineering and Design，2012，250：560-566.
［16］SHOELE K，MITTAL R.Flutter instability of a thin flexible plate in a channel［J］.Journalof Fluid Mechanics，2016，786：29-46.
［17］WANG W J，HUANG H B，LU X Y.Interplay of chordwise stiffness and shape on performance of selfpropelled flexible flapping plate［J］.Physics of Fluids，2021，33（9）：1-15.
［18］王文江.变刚度柔性体自主推进特性的流固耦合数值研究［D］.合肥：中国科学技术大学，2021.
［19］高丽敏，李晓军，谢建，等.对转压气机最先失速级的小扰动理论分析［J］.航空学报，2013，34（3）：533-540.