目前，膜结构风振响应机理仍不明确，风致事故仍时有发生。事实上，膜结构抗风分析的难点在于结构与风场之间的流固耦合效应，这种相互作用会加剧流场的分离和旋涡脱落，以致结构吸收的能量大于振动所消耗的能量，从而出现“气弹失稳”现象[1]。而对膜结构流固耦合效应的研究一直是结构风工程的热点和难点问题之一，这其中包括了附加质量和气承刚度对薄膜结构振动频率的影响等，通过试验研究认为流固耦合作用对结构振动频率有较大的影响[2]。因此，对膜结构进行风振响应分析时，应考虑流固耦合效应对膜结构振动频率的影响。

在膜结构振动频率的试验研究方面，Irwin等[3]基于蒙特利尔体育馆气弹模型的风洞试验，给出了一个较为简便的计算膜结构附加质量的估算公式。Sewall等[4]对三角形薄膜进行了试验研究，并提出了附加质量分布模型，通过该模型得到的三角形薄膜，其一阶自振频率与试验结果比较接近。王磊等[5]通过简化模型对圆膜在空气中的自振频率进行了研究，得到的计算频率与试验频率接近，并采用该模型对Sewall的三角形薄膜[4]进行验证，所得到的自振频率比Sewall提出的方法更为精确。孙晓颖等[6]和陈昭庆等[7]分别对开敞式和封闭式单向张拉膜结构进行试验研究，发现当表面流场中的旋涡主频与结构基频接近时，会出现涡激共振，从而导致气弹失稳，但是并没有对膜结构的振动频率进行深入研究。

在理论研究方面，目前的研究主要集中于附加质量对膜结构振动频率的影响。Minami[8]采用薄型机翼理论，推导了矩形平面薄膜在空气中的附加质量。杨庆山[9]和沈世钊等[10]较早地对附加质量进行了研究，认为静止流体中振动物体的附加空气质量与结构表面形状、结构振动方向以及雷诺数有关。毛国栋等[11]根据势流理论，推导了形状较为简单的薄膜在空气中振动的附加质量计算公式。孙晓颖[12]采用涡格法求解了流体在敞开式薄膜表面产生的涡强密度，并根据能量守恒定律得到了敞开式薄膜的附加质量。Zhou等[13]采用边界元方法对静风下薄膜结构进行分析，所得到的自振频率与已有的试验数据较为吻合。Chen等[14]推导了单向张拉薄膜在均匀流场中附加质量的解析表达式，并通过风洞试验加以验证。

上述研究主要针对附加质量对膜结构的振动影响，但是实际环境中膜结构与空气动力相耦合，有诸多气动参数会对结构的振动频率产生影响。针对以上不足，文中采用简化气弹模型对考虑流固耦合效应的薄膜结构振动频率进行研究。采用升力面理论对敞开式薄膜结构气动力进行研究，通过涡格法近似求解开敞式薄膜的振动频率。基于势流理论，建立空气与结构相互耦合的动力平衡方程。通过边界元与有限元方法求解封闭式薄膜结构振动频率。通过对开敞式薄膜和封闭式薄膜在来流风下的振动频率特性进行研究，分析两种膜结构在不同条件下振动频率的变化规律和影响因素。

由于流固耦合的复杂性，文献[15]中提出了一种简化气弹模型，即膜结构在风荷载作用下的动力方程可表示为

式中:Ms、Cs、Ks分别为结构的质量、阻尼系数和刚度矩阵;F为气动力，是时间t、结构位移x (t) 、速度x· (t) 、加速度x·· (t) 的函数。

由于式 (1) 等号两边都含有结构运动项，需要通过简化方法，对气动力方程进行解耦。因此，引入以下两个基本假定:

1) 拟定常假定。该假定认为结构运动只对此时的表面风压产生影响，不会影响流体本身的脉动特性。因此，气动力方程可近似表示为

式中:p (t) 为流体自身脉动产生的气动力;f[x (t) , ]为结构运动产生的附加气动力。

2) 强相关假定。附加气动力在频域上接近结构主振动频率的分量，这是由于受结构运动影响较大，而远离结构主振动频率的部分则可以忽略不计。为此，对式 (2) 进一步解耦，可表示为

式中: 分别为气动力与结构加速度、速度、位移的相关系数。

将式 (3) 代入式 (1) 中，通过整理可得

式中:Ma、Ca、Ka分别为附加质量、气动阻尼和气承刚度。

式 (4) 为简化气弹模型，不考虑具体的流体环境，将流体与结构视为整体，通过数学模型将二者联系在一起。

由结构动力学可知，质量和刚度对结构的振动频率影响较大，而阻尼的影响可以忽略不计。对于开敞式薄膜，可不考虑气承刚度[10]，所以开敞式薄膜的振动频率只需考虑附加质量的影响;而对于封闭式薄膜，由于所采用的分析方法考虑了附加质量、气动阻尼、气承刚度的影响，不再进行近似处理。

如图1所示，对于开敞式薄膜结构，由升力面理论可知，来流在薄膜表面产生的旋涡可以采用附着涡和自由涡表征。附着涡和自由涡沿x向和y向均有变化，沿y向相邻两剖面间拖出的自由涡环量即为这两个剖面上附着涡的环量差[16]。假设附着涡在dxdy面元上的附着涡密度为γm (x, y, t) ，则环量dΓ (x, y, t) =γm (x, y, t) dx，自由涡密度 。由此建立y向自由涡线密度与x向附着涡之间的关系。

定义p1、p2为薄膜结构下表面和上表面的扰动气压，假设空气不可压缩、无黏性和无旋转，速度为V，则根据非定常伯努利方程可得

式 (5) 中的速度分量Vx1=V+u1, Vy1=v1, Vz1=w1, Vx2=V+u2, Vy2=v2, Vz2=w2，其中u1、v1、w1为薄膜下表面扰动速度，u2、v2、w2为薄膜上表面扰动速度，气流的速度势函数 可采用扰动速度势函数 和来流产生的速度势Vx的叠加，则有? 

由环量的定义有如下关系式:

根据绕流物体的边界条件，z向空气扰动速度w为[17]

其中:Zs为薄膜结构的运动曲面函数Zs (x, y, t) =Z0 (x, y) +h (x, y, t) ，其中Z0 (x, y) 为薄膜结构初始应力下的曲面函数，h (x, y, t) 为薄膜的位移。

将式 (6) ～ (9) 代入式 (10) ，可得:

由速度势理论可知，扰动速度势为 ，则有

将式 (6) ～ (9) 、 (11) ～ (13) 代入式 (5) 中，可以得到薄膜表面总压强差，即

由式 (14) 可知，若求解得到附着涡密度γm (x, y, t) ，即可以得到薄膜表面总压强差p。为此，采用涡格法来近似求解γm (x, y, t) 。

将图2所示的整个投影平面可以看成由有限个格子组成，每个格子上均布置有马蹄涡，且其马蹄涡的涡强是常数，但是不同马蹄涡的涡强各不相同，取每个格子3/4弦线处的中点作为控制点，在控制点上求其数值解。

根据毕奥-萨伐尔定理，得到每个马蹄涡对控制点的诱导速度的表达式。由于扰动速度u、v远小于风速V，则根据式 (10) 可知，每个控制点的诱导速度可以由结构位移来表示。由此可以得到结构位移与马蹄涡无量纲涡强γj的代数方程。

为求解γj，假定翼展长度为l，定义无量纲涡强γj=Гj/Vl，则格子翼的平均附着涡密度为

由式 (15) 计算得到的是每个格子翼的平均附着涡密度，将其代入式 (14) ，可以求解得到薄膜结构表面的总压强差。根据能量守恒定律，膜周围的空气在dt时间内的动能负增量等于作用在膜结构上的总气动压力p在dt时间内做的功d W[18]，即:

其中， 为膜结构沿法线方向的速度，Ma为薄膜结构单位面积上的附加质量。

联立式 (16) 、 (17) 可得

采用文献[14]的单向张拉膜结构气弹模型风洞试验，模型尺寸为1.2 m×0.6 m，将膜的长边固定，短边无约束，底部支撑为开敞式。膜材料为乳胶薄膜，面密度ρs=0.413 3 kg/m2，厚度0.4 mm，来流为均匀流场，风向垂直于长边方向。试验风速和预张力见表1。

采用涡格法对上述试验模型进行计算，其中第n阶振型下的结构位移方程为

其中，A为结构振幅，l为薄膜跨度，取0.6 m，第n阶振动频率为 

根据式 (19) 可以求得无量纲涡强γj。根据式 (14) 、 (15) 、 (18) 、 (19) ，可以得到结构的附加空气质量。考虑附加质量的结构振动中，其一阶振动频率与试验的一阶频率见表2。

由表2可知，采用涡格法对上述算例进行计算，结果比较理想，可见，采用涡格法对敞开式薄膜进行附加质量计算具有可靠性。

对于封闭式薄膜，假设其支撑在无限大的刚性平面上，同时假设空气是不可压缩的理想势流，具有无黏性、无旋转的特点。薄膜上表面空气流速为V，下表面流速为零 (图3) 。空气的扰动速度势满足拉普拉斯方程，即。通过边界元直接法，得到薄膜上P点的速度势函数? (P) 在其区域内与边界值? (Q) 的关系式[19,20]，即

其中:S为薄膜表面积;r (P, Q) 表示点P到Q点的距离。

在平面积S上的边界条件为Neumann类型，是薄膜和空气相互耦合的条件，可表示为

其中，w为薄膜的竖向位移，w=w (x, y, t) ;? 为薄膜位移对来流方向的偏导数。

在结构上表面产生的扰动压强[20]为

式中，ρ为空气密度， 为速度势函数?对来流方向的偏导数。

当速度V=0时，通过式 (20) ～ (22) 可以得到膜结构下表面扰动压强p2，薄膜最终所受的总气动压强p=p1 (P, t) +p2 (P, t) 。假定μ代表薄膜的单位面积质量，T1和T2分别代表薄膜在x和y向上的单位长度张力，则膜结构的线性动力方程为

将速度势、气动力和位移的时间变量与空间变量进行分离，有  (其中λ=γ±iω，ω为振动圆频率) ，则式 (20) ～ (23) 可以表示为:

对式 (25) ～ (27) 采用Galerkin法[21]建立有限元方程组。将薄膜表面离散成若干三节点三角形单元，得到以下代数方程:

其中: 分别为单元节点上相关的物理量组成的向量，矩阵B1、B2、B3、B4是由有限元单元矩阵B1e、B2e、B3e、B4e组成，相应表达式分别为:

式中，N=[N1, N2, N3]是单元形函数组成的矩阵，采用高斯积分进行计算。

边界元划分与有限元一致，采用同样的三节点三角形单元对薄膜结构进行离散，采用边界元法[22]对边界积分方程 (24) 建立矩阵方程:

式 (34) 中矩阵A的元素Amn表达式为

其中:Nn (j) 为第n个节点在单元j上的形函数，e为含有节点n的有限元单元的数量，r (P, Q) 为点P与点Q之间的距离，Sj为单元j的面积。

采用高斯积分计算方程式 (36) ，当m=n时，其为1/r类型的奇异积分，则需采用特殊的数值方法计算[23]。

根据式 (28) 、 (29) 、 (35) 可以得到薄膜上表面压强为

其中B5=B1-1B2，当V=0时，可以得到薄膜下表面压强，即有

则薄膜结构总压强为

将式 (39) 代入结构离散动力方程式 (30) ，可得

式中:K为刚度矩阵;C为阻尼矩阵;M为质量矩阵。刚度矩阵K=Ks+Ka，其中薄膜刚度矩阵为Ks=-T1B3-T2B4，空气刚度矩阵为Ka=VρB1B5AB5;质量矩阵M=Ms+Ma，薄膜质量矩阵为Ms=μB1，空气质量矩阵为Ma=2ρB1A;阻尼矩阵C=VρB1 (AB5+B5A) B5AB5。

采用文献[5]的圆形薄膜振动试验进行验证，模型为直径0.3 m的圆形薄膜，放置于真空箱中，试验时利用真空泵改变真空箱内的气压，用以测试不同气压下膜单元的自振特性。由于试验工况较多，选取了环境激励下乳胶薄膜的试验工况进行计算。乳胶薄膜单元四周固定，通过顶升边界的方法均匀施加张力，包括三种不同应力水平，其中σ1=0.297 MPa、σ2=0.178 MPa、σ3=0.092 MPa。具体的试验工况见表3。

采用有限元与边界元相结合的方法对该试验各工况进行计算，图4为乳胶薄膜在环境激励下的各工况下计算频率与试验频率，表4为计算结果与试验结果的相对误差。

由于试验中没有量测得到第三阶振动频率，所以试验频率和误差的第三阶缺省。从表4中可以看出，计算结果与试验大部分结果较为接近，而C4工况的计算误差较大，这是由于试验测得的数据误差较大的缘故，可见，本文所采用的方法对封闭式薄膜的计算结果较为可靠。

采用上述方法对长方形膜和圆形膜两种典型形状薄膜的振动频率进行参数研究。对于长方形薄膜结构，风速分别为0、3、6、9、12、15 m/s，计算工况见表5。

图5为开敞式和封闭式膜结构的振动频率随风速的变化曲线，由图可知，振动频率随风速增大而减小。开敞式和封闭式薄膜的振动频率分别用f1和f2表示，从图5a和图5c可以发现，当张力较大，f1相对较小;而张力T=20 N/m时，f1、f2值接近;当风速V=15 m/s时，f2<f1。令沿来流方向与垂直来流方向的结构尺寸的比值为λ，由图5b和图5d可知，当λ较小时，f2>f1。当λ较大时，长边固定的长方形薄膜的f2相对要小。

对于薄膜结构，影响振动频率的主要因素是附加质量Ma，图6为上述各工况下Ma的变化曲线，由图可知，Ma随风速增大而增大，这是造成膜结构振动频率随风速增大而减小的主要原因。令开敞式薄膜和封闭式薄膜的附加质量分别为Ma1和Ma2，从图6a和图6c可以看出，长方形膜四边固定时，在张力T=100 N/m时，Ma1>Ma2，从而使f1<f2;而在其他情况下，Ma1<Ma2，而在大部分情况下f2>f1，这是由于封闭式薄膜存在气承刚度，使得刚度增加，从而导致f2>f1。在张力T=20 N/m时，因为Ma2>Ma1，使得f2<f1。由图6b和图6d可知，当λ较小时，Ma2相对于Ma1要小，使得f2>f1。当λ较大时，Ma2相对于Ma1较大，而对四边固定的长方形膜，f1和f2的数值接近，这同样是由于封闭式薄膜存在气承刚度的原因;对长边固定的长方形膜，由于Ma2相对于Ma1要大的很多，f2<f1。

对圆形膜结构进行分析，其半径为0.15 m，周边固定。风速V分别为0、3、6、9、12、15 m/s，张力T分别为20、50、100 N/m。对于开敞式薄膜结构，工况X5;对于封闭式薄膜结构，工况Y5。

图7为开敞式和封闭式圆形膜在不同张力下振动频率随风速的变化曲线。由图7可知，f1在各张力下稍大于f2，这是由于圆形膜的λ较大，Ma2相对于Ma1较大 (图8) 。

1) 张力、结构尺寸、风速等因素对结构频率均有影响。风速增大，振动频率降低;当张力较大时，开敞式薄膜的振动频率相对于封闭式薄膜的振动频率要小;而张力较小，两种薄膜结构的振动频率比较接近。当沿来流方向与垂直来流方向的结构尺寸的比值较小时，封闭式薄膜的振动频率大于开敞式薄膜;而比值较大时，长边固定的长方形薄膜会出现封闭式薄膜的振动频率小于开敞式薄膜的情况。

2) 附加质量是影响薄膜结构振动频率的主要因素，附加质量随风速增大而增大，从而导致振动频率下降。封闭式薄膜的附加质量大于开敞式薄膜的附加质量，而开敞式薄膜的振动频率却没有比封闭式的大，这说明气承刚度对封闭式薄膜振动频率的影响不可忽视。