充气膜结构是一种典型的曲面张力空间结构，广泛应用于大型体育场馆、工作车间、仓库等建筑上，具有极大的发展空间和前景[1]。然而，充气膜结构作为一种柔性结构，对风十分敏感[2]。同时，结构的流固耦合效应和脉动风效应对结构的风压分布和风致响应具有较大影响[3]，这在大多数的实际设计中并未同时考虑，此外，国内对于截椭球形大跨空间结构的抗风设计相关研究也相对较少。

在诸多充气膜结构的风振特性的研究方法中，计算流体动力学技术(CFD)因适应性强、周期短、费用低等优势而得到广泛应用[4]。本文基于CFD技术和计算结构动力学技术(CSD)研究了截椭球形充气膜结构的风荷载响应特性。首先，利用线性滤波器法通过Matlab得到了考虑脉动风的风速时程曲线，并编制了UDF将其导入到Flunet。其次，利用Workbench17.0平台数值模拟并得到了风致响应稳定后不同风向角、高度、内压下刚性结构和柔性结构的风压系数及位移响应分布并对得到的相关结果进行了分析。

一般情况下，工程上大多将常态风假设为不可压缩粘性低速流体，并需要遵循质量守恒定律和能量守恒定律。同时，为了考虑并简化流体的湍流运动，人们给出了几种不同的湍流数值模拟方法。其中，求解时均化Reynolds方程的Reynolds平均法(RANS)计算量较小，因而在工程上得到了广泛的应用[5]。本文采用了RNG k-ε模型，因为RNG k-ε模型在修正后能更好的处理弯曲流线流动[6]。该湍流模型的控制方程如下：

∂ρk∂t+∂ρkui∂xi=∂∂xj(αkμeff∂k∂xj)+Gk+ρε∂ρk∂t+∂ρkui∂xi=∂∂xj(αkμeff∂k∂xj)+Gk+ρε(1)

∂∂tρε+∂∂xi(ρεui)=∂∂xj(αεμeff∂ε∂xj)+C∗1εεkGk−C2ερε2k∂∂tρε+∂∂xi(ρεui)=∂∂xj(αεμeff∂ε∂xj)+C1ε*εkGk−C2ερε2k(2)

其中，ρρ为流体密度，k为湍动能，μiμi为湍动粘度，μeffμeff为湍动能有效粘性系数，εε为湍动能耗散率，GkGk为由平均速度梯度而导致的湍动能产生项，其他系数项一般取：C∗1ε=C1ε−η(1−η/η0)1+βη3C1ε*=C1ε−η(1−η/η0)1+βη3，η=(2Eij⋅Eij)1/2kεη=(2Eij⋅Eij)1/2kε，Eij=12[∂ui∂xj+∂uj∂xi]Eij=12[∂ui∂xj+∂uj∂xi]，C1ε=1.42C1ε=1.42，C2ε=1.68C2ε=1.68，η0=4.377η0=4.377，β=0.012β=0.012，αk=αε=1.39αk=αε=1.39。

对于充气膜结构这类柔性大变形结构，需要考虑结构和流体之间的相互耦合作用。本文采用的是分区弱耦合算法，即在每一时间步内分别独立求解流体和结构的控制方程，然后将流体域中求得的风荷载通过耦合交界面传递给结构域，得到结构的变形，然后再将结构的变形通过耦合交界面传递给流体域，往复进行计算，直至结果趋于稳定。Glück等人[7]较早地采用结构和流体的商业计算软件实现了薄膜结构的流固耦合数值模拟。Nayer等人[8]根据文献[9]气弹模型试验的基础上，基于大涡模拟的分区耦合算法进行了在不同雷诺数下充气截球形膜结构流固耦合效应的数值模拟研究，并根据涡脱落频率和斯特劳哈尔数对涡脱落进行了分类。武岳等人[10]、陈怡然等人[11]和李如地等人[12]基于分区弱耦合算法分别对二维张拉膜结构屋面、山脉形膜结构及典型细长索膜结构进行了流固耦合数值模拟，得到了结构的位移响应时程。

建立2500m×1000m×400m的长方体为流场域，结构纵剖面与风场剖面的面积阻塞率最大值仅为1%，能够满足阻塞率小于3%的要求[13]，并且结构距边界入口5倍长度。根据我国规范，该工程处于B类地貌，基本风速为20m/s，相当于基本风压0.4，梯度风高度350m，标准高度10米，地面粗糙度0.16。流体域入口采用指数率平均风剖面，相关公式如下：

Uˉˉˉ(z)=Uˉˉˉ1(zz1)αUˉ(z)=Uˉ1(zz1)α(3)

k=1.5(uˉ×I)2k=1.5(uˉ×I)2(4)

ω=0.090.75k0.5/Lω=0.090.75k0.5/L(5)

其中z1z1和Uˉˉˉ1Uˉ1分别为参考高度及其基本风速，uˉuˉ为平均风速，ω为耗散率，I为湍流强度，L为湍流积分尺度。

将Matlab中按照达文波风速谱模拟生成的脉动风速时程通过UDF输入到Fluent环境中，参考的线性滤波器法原理详见文献[14]。风场高度10米处的风速时程如图1所示。同时，将该风速谱进行FFT变换，与Davenport进行对比，如图2所示。可以发现，两者功率谱较为吻合，说明了该方法具有可行性。

流体域两侧及顶部采用滑移对称边界条件，底部采用无滑移壁面边界条件，出口采用自由出流边界条件[15]。模型采用四面体网格划分，同时要注意对膜面及附近流场局部进行加密。通过试算和网格无关性检验，确定该处最小网格尺寸0.5m，增长率1.2，网格总数达385万个左右，底部网格划分情况详见图3。

求解器设置方面，本文对流体控制方程采用SIMPLEC算法，选用二阶迎风格式对控制方程进行离散，同时模拟过程中迭代计算的收敛标准为无量纲连续性残差达到 10-5并且结构监测点风压系数趋于稳定，前后迭代误差在5%内。

实际工程中充气膜结构膜材采用1.0mm厚的PVDF涂层膜材和0.8mm厚PVC保温内膜，膜材表面附加有外包PE的Φ10钢索。为了简化模型，我们根据文献[2]给出的方法对模型进行等效处理。本文假定索与膜面无相对滑移，将膜材和索等效为均匀材料，等效后模型在均布荷载作用下所产生的位移相同。最后得到模型的等效密度为2500kg/m3，膜厚1mm，等效弹性模量为1×1091×109Pa，泊松比为0.34。结构底面长轴150m，短轴100m，底边采用固定约束，网格采用shell181划分，结构具体尺寸如图4所示。

将结构与流场的交界面定义为流固耦合交界面，该面作为流体域和结构域数据传递的交互面。

充气膜结构的刚度主要由刚度决定[16]，同时，相关风洞试验[17,18]也证实了结构的高度和风向角也会对结构表面风压分布产生较大影响。因此，本文主要的研究变量为风向角、高度和内压。其中，风向角为变量时，高度取30m，内压取400Pa，风向角分别取0°、15°、30°、45°、60°和90°；高度为变量时，风向角取0°，内压取400Pa，高度分别取26m、30m和34m；内压为变量时，风向角取45°，高度取30m，内压分别取400Pa、500Pa、600Pa和700Pa。

本章以某鞍形张拉膜为例，基于上述方法数值模拟得到了不同风向角下的风压系数分布，并与同济大学TJ2风洞的试验结果[19]进行对比，来说明这一方法的可行性。

该鞍形张拉膜具体参数和模型见表1和图5，流体域尺寸大小为1600m×800m×80m，模型中心距离流场入口处 200m，阻塞率仅为 0.5%。流场域如图6所示。网格选用四面体单元，膜面及附近流场局部需要加密。流场边界条件、湍流模型和求解器设置与前文提到的一致。

本小节给出了B类地貌、15m/s基本风速下，风向角为0°、45°和90°时数值模拟和风洞试验方法下结构的风压系数分布情况详见图7—9。

从图7—9可以看出，在三个不同风向角下数值模拟和风洞试验得到的风压系数分布大致相同，仅在数值上和极少部分区域分布上有差异，这主要可能是风洞试验相似性原理所产生的误差。

因为本文主要关注的是风向角、高度和内压变化对结构表面风压分布的影响，所以给出了结构平均化的风压系数分布图及汇总表。由于篇幅限制，本文仅给出部分结果。

表2给出了不同风向角下充气膜结构柔性模型的风压系数极值汇总表。图10和图11给出了0°刚性模型以及0°、45°、60°和90°柔性模型的风压系数分布图。图12给出了各风向角下柔性模型风压系数沿长轴与短轴的分布情况。

从图10可以发现：相同风向角下刚性和柔性模型的表面风压分布大致相同，仅在数值上有差异，刚性模型中正压和负压系数极值分别为0.853和-0.738，柔性模型中正压和负压系数极值分别为0.784和-0.763；同时结合表2发现结构总体上受负压极值控制，并且考虑流固耦合效应的柔性模型负压系数极值更大，风压分布更均匀[2]。

从图10(b)和图11可以发现，随着风向角的变化，结构风压的大小和分布都发生了较大变化。不同风向角下等值线呈椭圆状分布，风压带大致沿风向垂直方向分布，并且风向角0°和90°时沿结构的中轴线对称分布；大部分区域以负压为主，最大负压出现在结构顶部中心区，此处容易产生强烈的流动分离；最大正压出现在迎风面底部，同时由于流动再附，在背风面也会存在正压区域。随着风向角的增大，风压分布带逐渐变得狭长。风向角90°时结构正压和负压系数极值达到最大，分别为1.35和-2.01，说明该工况下结构最容易因受到较大上拔力而发生破坏。

从图12可以发现，各风向角下风压系数值大致沿短轴与长轴呈对称抛物线形分布。结构顶部中心处负压最大，风压梯度变化最大；风向角0°和45°时结构沿短轴与长轴两端30m和25m范围内存在正压；随着风向角增大，风压系数沿短轴与长轴的变化幅度增大，90°时变化最大。同时，风向角0°和45°时风压系数沿长轴的变化幅度更大，60°和90°时风压系数沿短轴变化幅度更大。

图13给出了26m和34m柔性模型的风压系数分布图，图14给出了26m和34m给出了各风向角下柔性模型风压系数沿短轴与长轴的分布情况。

从图10(b)和图13可以看出，随着结构高度的增加，结构顶部负压系数极值分布带逐渐变得狭长，负压系数极值逐渐增大，26m、30m和34m时分别为-0.691、-0.763和-0.865。从图14可以看出，风压系数值大致沿短轴与长轴呈对称抛物线形分布。其中结构顶部中心处负压最大，风压梯度变化也最大；结构沿长轴两端25m范围内存在正压，沿短轴方向都为负压；随着高度增加，风压系数沿短轴与长轴的变化幅度逐渐增大，34m时变化最大，并且结构沿长轴的风压系数变化幅度更明显。

图15给出了500Pa、600Pa和700Pa柔性模型的风压系数分布图，图16给出了各内压下柔性模型在短轴与长轴的分布情况。

从图11(a)和图15可以看出，随着内压的增加，负压系数极值带逐渐变得狭长，正压区位于结构的迎风面和背风面。而从不同内压下结构正负压系数极值变化趋势发现，随着内压的增加，结构正压和负压系数极值先增加后分别保持在1.41和-1.45左右不变。

从图16可以看出，风压系数值大致沿短轴与长轴呈对称抛物线形分布。其中结构顶部中心负压最大，风压梯度变化也较大；结构沿长轴两端25m范围内存在正压，沿短轴方向都为负压；随着内压增加，内压500Pa—700Pa时风压系数沿短轴与长轴的变化幅度较小，说明内压的增加不一定能有效地减少风压系数极值，该工况下内压宜采用400Pa。

本文研究了考虑流固耦合与脉动风效应下不同风向角和结构参数对截椭球形充气膜结构的风压系数的影响，得到了以下结论。

(1)不同风向角下，结构风压系数分布变化较大，分布带基本都沿垂直风向呈椭圆状分布；结构正压区分布范围较小，主要集中在迎风面和背风面；结构负压区分布较大，最大负压区在结构顶部中心位置，在风向角为90°时有最大负压系数极值；各风向角下风压系数值大致沿短轴与长轴呈对称抛物线形分布，并且随着风向角的增大变化幅度增大，0°和45°时风压系数沿长轴的变化幅度更大，60°和90°时风压系数沿短轴变化幅度更大。

(2)不同高度下，结构表面风压分布大致相同，顶部负压系数极值分布带逐渐变得狭长，并随着高度的增加负压极值增大；随着高度增加，风压系数沿短轴与长轴的变化幅度逐渐增大，并且结构沿长轴的风压系数变化幅度更明显。

(3) 随着内压的增加，结构负压系数极值先增加后保持在-1.45左右不变，说明内压的增加不一定能有效地减少风压系数极值，在其他条件相同的情况下内压宜采用400Pa。