建筑膜结构(下文简称膜结构)作为一种新型结构，其具有质量轻、刚度小、跨度大等优点，现已被广泛百富策略白菜网于大跨度空间结构中[1]。然而其缺点也十分明显，对外界荷载作用十分敏感，易产生振动、大变形以及撕裂破坏等，导致膜结构无法继续适用甚至造成工程事故[2]。目前，针对膜结构在动力荷载作用下的相关研究多集中于强风、暴雨等荷载作用，并且已经取得了诸多突破。1986年Kimoto等[3]对比分析了单向悬挂膜结构在不同条件下的气动稳定性，得出不同风速下具有不同的结构振幅峰值。同年，Uematsu等[4]对鞍形张拉膜结构进行了风洞模试验研究，分析出膜结构振动的变化情况。向阳等[5,6]和武岳等[7]设计制作了鞍形和伞形张拉膜结构的模型并分别进行了风洞试验研究，提出膜结构的振动为随机脉动风场激励下的强迫振动，不同的相关条件均会对结构响应产生一定的影响。武岳等[8,9,10]进行了伞形膜结构风荷载作用下的试验研究，得出膜面位移风振系数随着膜面倾角不同而不同，并给定了相关参数建议。杨庆山等[11,12,13]研究了薄膜结构的风致动力响应、气弹动力效应和失稳临界风速，并编制程序对膜结构非线性风致动力响应进行分析，指出膜材的褶皱对膜结构的失稳有一定影响。Guo等[14]对矩形张拉膜结构在冲击荷载作用下的动力响应进行了试验研究，得到了位移-时间曲线、位移功率谱密度曲线和振动频率。

综上，膜结构在外界荷载作用下的研究已由试验上升到计算机模拟再到理论，逐步走向成熟，同时也反映了膜结构是一种对外荷载十分敏感的柔性结构，易产生较大的变形。冰雹作为常见的自然灾害之一，其对人类活动的影响较大，并且较强的冰雹冲击荷载对柔性建筑物(如膜结构等)可能带来更为严重的破坏。所以，对膜结构在冰雹冲击作用下的动力响应研究显得非常重要。一方面，冰雹冲击作用往往会使膜结构产生较大的振动，加速应力松弛，从而影响膜结构的适用性与安全性；另一方面，目前关于膜结构的设计规范规定在膜结构设计时需考虑风、雪荷载效应，而未将冰雹等竖向冲击荷载纳入荷载效应组合[15]。因此，本文以典型的三维空间膜结构—伞形膜结构作为研究对象，通过不同粒径的冰雹进行冲击试验，得到伞形膜结构的动力响应和张拉松弛率等数据。然后进行数值模拟，并与试验结果对比分析。结果表明，冰雹冲击会造成伞形膜结构产生大变形，并加速应力松弛，甚至导致膜面撕裂破坏。该结果为更精确地研究群冰雹冲击伞形膜结构的动力响应打下基础，也为张拉膜结构的设计、施工和维护提出工程参考依据，以防止膜结构在冰雹冲击作用后，发生褶皱、大变形，进而防止局部撕裂破坏、整体风致失稳等工程事故。

根据CECS 158—2015《膜结构技术规程》,本次试验的膜材采用P类膜材，基本试验参数如表1所示。

试验用的试件模型采用矢跨比h/L为1/3,跨度L为1.5 m, 高度h为0.5 m的空间曲面模型，即伞形模型。试件是由初始平面形态的膜材通过专业膜结构找形分析软件3D3S进行找形分析，并确定最终的伞形膜结构形态特征。在找形分析时也一并做了裁剪分析，裁剪分析是为了将平面膜材搭接成为空间曲面模型而又不产生褶皱，裁剪分析时得到各膜片尺寸，并确定膜片搭接方法。在裁剪分析得到膜面尺寸后，经由专业膜结构公司对其加工，采用热合的方式进行连接，然后将确定的边索穿入膜套中，并制备好角点连接件，如图1所示。

在现代气象学上，冰雹专指直径在0.5 cm以上的固体降水[16]。根据冰雹自然形成过程，其形状多以球形和椭球形出现[17],故本次试验采用硅胶模具制作冰球来模拟冰雹，如图2所示。

本次采用的试验装置分三大类，分别为试验张拉装置、试验加载装置以及数据采集装置。其中试验张拉装置主要采用改进的“十字形螺旋杆式膜材张拉装置”进行膜材张拉，其初始设计平面图、试验装置如图3(a)所示；试验加载装置主要有预应力加载装置、冰雹释放装置组成。预应力加载装置是将膜面四个角点与张拉装置连接，膜面预应力通过边索以及角点张拉共同作用施加。边索穿过膜面预留的模套，通过两头螺栓转动施加边索拉力，测量边索张拉长度从而确定边索拉力。而伞形膜的形态是由膜面预应力以及边索张拉确定，故试件四边边索采用直径为16 mm的6′19IWRC型钢丝索，同时膜面中部固定有直径为12 cm的钢柱支承以提供试验所需伞形膜结构矢跨比，其示意图如图3(b)和图3(c)所示。冰雹释放装置由框架和铁纱网组成，铁纱网平面尺寸为1 800 mm×1 000 mm, 铁纱网方格大小为2.5 mm×2.5 mm。释放冰球时，将铁纱网定位于伞形膜结构上方一定高处，并将冰球中预制的导线穿过铁砂方格网，通过铁砂方格网定位冰球位置，定位完成后将导线束成一端并固定，待冰球静置后解除约束释放导线，从而使冰球做自由落体运动并冲击伞形膜结构，如图3(d)所示；数据采集装置主要有动力响应测量装置ZLDS100、边索拉力测量装置HP-10K、膜面预张力测量装置NK-200等组成，如图3(e)所示。

伞形膜结构为对称结构，试验时取1/4膜面进行试验，测点分别位于膜面跨度1/12,1/8,1/4及1/2处，共计6个测点，测点具体布置如图4所示。

气象资料显示，观测到的冰雹中有80%～90%的冰雹直径小于3.0 cm, 然而也有大直径冰雹出现。故本次试验采用不同粒径的冰球来模拟冰雹进行冲击试验，其粒径分别为1.7 cm, 2.5 cm, 3.0 cm, 4.5 cm, 6.0 cm。

假设冰球冲击膜面之前做自由落体运动，势能转化动能，试验冰球冲击膜面时的冲击速度等于理论计算的末速度，而且通过计算只要高度达到24 m以上，可以认为冰雹末速度稳定在理论计算的末速度，因此根据冰雹末速度公式，冰雹在空气中下落时遇到的阻力为F即[18]

F=CdSρ2V2F=CdSρ2V2 (1)

式中：Cd为阻力系数；r为空气密度，kg/m3;V为冰雹相对于空气的下落速度，m/s; S为冰雹最大截面积，cm2。冰雹达到末速度时，阻力与重力平衡，即

mg=CdSρ2V2mg=CdSρ2V2 (2)

式中：m为冰雹质量，g; g为重力加速度，m/s2;而对于直径d为的圆球形冰雹，cm, 这里的m=π6d3ρim=π6d3ρi,其中ρi为冰雹密度；冰雹最大截面积S=π4d2S=π4d2。将其代入表达式(2)中，整理可得圆形冰球末速度的普遍公式

VS=(43⋅ρiρ⋅gCdd)12VS=(43⋅ρiρ⋅gCdd)12 (3)

为更精确地得到圆形冰雹末速度VS(m/s),从而引入雷诺数Re=ρVdμRe=ρVdμ,其中，μ为空气动力黏性系数，得出

VS=43d2⋅gρiμ(CdRe)−1VS=43d2⋅gρiμ(CdRe)-1 (4)

而直径不同的冰雹在下落时其阻力不尽相同，故将阻力系数Cd分段而论：

当冰雹直径为0.5 cm≤d≤5.0 cm时，阻力系数取值为Cd=0.41d0.185;

当冰雹直径为5.0 cm≤d≤7.5 cm时，阻力系数取值为Cd=0.575-0.020 8(d-4)2;

当冰雹直径为6.4 cm≤d≤16.0 cm时，阻力系数取值为Cd=0.109d0.253。

因此，由式(4)可计算得到不同球形冰雹的末速度，也就得到冰球冲击膜面时的冲击速度，则计算结果如表2所示

试验采用5种不同粒径的冰球，并按不同粒径进行编号，其试验工况如表3所示。

试验中膜面预应力为3 MPa, 边索拉力为7.5 kN,冰球置于伞形膜结构张拉装置正上方25.2 m处的冰雹释放装置中，通过1.7 cm, 2.5 cm, 3.0 cm, 4.5 cm, 6.0 cm 5种不同粒径冰球释放后做自由落体运动来冲击伞形膜面点D,利用数据采集装置测得点D的振动位移时程和速度时程，如图5和图6所示。

分析图5可知，冰球冲击膜面时，膜面冲击点振幅明显增加，且不同粒径冰球冲击的最大振幅也不相同：粒径分别为1.7 cm, 2.5 cm, 3.0 cm, 4.5 cm, 6.0 cm的冰球，其冲击的最大振幅为1.99 mm, 3.54 mm, 7.06 mm, 11.03 mm, 16.81 mm。由冰雹末速度公式计算可知冰球直径越大则速度越大，因此冰球冲击膜面时膜面抵抗冰雹冲击所产生的变形也越大。

由于伞形膜结构具有阻尼作用，故冰球冲击膜面后会有振动衰减现象，但冰球的直径越大衰减也就越缓慢，这是因为直径越大的冰球冲击膜面时其动能较大，膜面振动也越剧烈，消散动能的时间也就越长。

由图6可知，不同粒径的冰球冲击膜面的速度响应也不相同，冰球直径越大，膜面速度响应也越大；试验中直径为6 cm的冰球冲击产生的速度响应最大，数值大小为9.36 cm。

为更直观的表示膜面冲击点D的位移时程、速度时程变化趋势，如图7所示。

由图7所示，膜面测点D受冰球冲击的振幅随着冰球粒径的增大而增大，冰球粒径为1.7～2.5 cm冲击膜面的振幅增长平缓，冰球粒径为2.5～6.0 cm冲击膜面的振幅呈近似线性增长；膜面测点D受冰球冲击的峰值速度随着冰球粒径的增大而增大，冰球粒径为1.7～4.5 cm冲击膜面的峰值速度增长较快，但在冰球粒径4.5～6.0 cm冲击膜面的速度增长放缓。

这里将粒径为3.0 cm冰球置于伞形膜结构张拉装置正上方25.2 m处的冰雹释放 装置中，释放后做自由落体运动冲击伞形膜面点D,通过数据采集装置得到测点A～F共计6个测点时程图，不同测点的位移、速度时程，如图8和图9所示。

由图8可知，冰球冲击点为测点D,则测点D的振幅最大，最大振幅为7.06 mm。在膜面经线方向，C,B,A三个测点随着离点D距离依次增大，在点D受到冰球冲击后，膜面振动传递以及耗能消散有滞后的现象。C、B、A三个测点的最大振幅也呈递减的趋势，最小振幅出现在测点A,振幅大小为3.35 mm, 分析其原因主要是测点A距离冲击点较远，且与钢柱支承较近，冰球冲击所带来的能量主要由钢柱耗散。在膜面纬线方向，E,F两个测点呈对称分布，所以两测点的最大振幅十分接近，又由于两测点距边索较近，受边索阻曳作用明显，其最大振幅小于点D。

由图9可知，冰球冲击点D后，各点振动波在传递时具有滞后现象， F,E,C,B,A测试的经过0.54 s, 0.39 s, 0.46 s, 0.43 s, 0.51 s后达到峰值，这符合振动波传递特点。膜面离点D越近的点C,E,F的最大振幅相对较小，而A,B两点最大振幅相对较大。这是由于C,E,F三点均处于膜面能量耗散较小区域，变形较大，其势能也较大，而动能较小，因此速度时程振幅较小。而A,B两点距离钢柱支承较近，由于钢柱的能量耗散，故势能较小，而动能较大，进而速度时程振幅也较大。

伞形膜结构是由四边边索以及角点张拉共同提供膜面预应力来保持其一定的刚度，在荷载的作用时候，结构才会有足够的承载力以及形态上的完整。由于时变特性，索膜材料会产生应力松弛，而为研究荷载冲击对伞形膜结构的松弛产生怎样的影响，本文主要针对边索的拉力变化进行了实验分析。为检测试验中冰球冲击前后对边索拉力变化，试验中采用数显式推拉力计，在冲击前后对四个边索监测并读出其拉力大小N1,N2。这里定义边索松弛率为(N1−N2)N1×100%(Ν1-Ν2)Ν1×100%,则试验结果如表4和图10所示。

由以上可知，伞形膜结构受到冰球冲击的瞬时，其四边边索拉力均有损失，从而导致伞形膜面预应力降低，其原因是因为冰球冲击荷载较大，使索和膜产生了不可恢复的塑性变形。从而说明冰雹的冲击会加速伞形膜结构的松弛，再遇风雨等其他荷载时会加大膜面的振动变形，导致伞形膜结构撕裂破坏。

本文使用了通用的显式动力学有限元分析软件ANSYS/LS-DYNA,基于隐式-显式顺序解法，模拟了冰雹对膜的冲击过程。

首先，在隐式计算中，膜采用了Shell41单元，索选择LINK10单元，冰雹采用了Solid185单元。而在显式动态分析时，相应地转换为SHELL163单元、LINK167单元、SOLID164单元。网格划分中，采用了映射三角形单元来生成膜表面以及索的网格，映射六面体实体单元来生成冰雹的网格。索网格划分大小为0.5 cm, 冰雹冲击点膜面网格大小为0.25 cm, 冰雹网格大小为1 cm。划分网格完成后，整体模型单元个数为8 544个，膜单元为7 559个，而索单元为121个，冰雹单元为864个。在显式动力学分析中，将表2所示的冲击速度百富策略白菜网于冰雹，并对SHELL和SOLID单元设置了面面自动接触，接触面为冰雹表面、目标面为膜面，并为完全弹性碰撞，恢复系数为1。设置了相应的接触面参数来描述冰雹大变形接触和动态冲击时膜间的复杂相互作用。

模拟时的参数设置选用与现场试验相同的材料属性参数，如表5所示。考虑到实验中膜材搭接处由于搭接造成的膜面厚度增大，而具有脊索效应。所以数值模拟中为了考虑这一因素而在搭接处加入了脊索。同时在进行冰雹冲击伞形膜结构数值模拟时进行了如下假设：

(1) 数值模拟中索、膜均只承受拉力，不承受压力，这与试验伞形膜结构相一致；

(2) 数值模拟中索为理想柔性，索的材料符合胡克定律；

(3) 数值模拟中索、膜在工作范围内均处于弹性阶段；

(4) 数值模拟中膜与四个边索设置为铰连接；

(5) 数值模拟中膜与中部的钢柱、膜角连接点以及脊索设置为刚性连接。

数值模拟中分别对5种不同粒径的冰雹：1.7 cm, 2.5 cm, 3.0 cm, 4.5 cm和6.0 cm冲击伞形膜结构膜面测点D进行模拟分析，通过后处理程序LS-PREPOST读出位移时程。以冰雹粒径6.0 cm为例，数值模拟冲击过程如图11所示。

如图11所示，数值模拟中冰雹冲击过程总结如下：

(1)冰雹冲击膜面测点D之前，膜面受力较为均匀，膜面预应力水平为3 MPa, 边索拉力为7.5 kN,最大应力出现在膜面脊索和膜面中心钢柱处；

(2)冰雹与膜面接触时，膜面通过变形来抵抗冰雹的冲击力，膜面产生由冰雹冲击带来的应力；

(3)当冰雹将所有势能转化到膜面上时，此时冰雹冲击对膜面产生的位移达到最大；

(4)膜面上的动能以波的形式向四周传递，由于伞形膜结构的特点，波首先沿纬向方向传播，再沿经向方向传播，因测点D距离边索较近，故波最快传递到边索上，并通过边索进行耗能，未耗散完的能量又以波的形式传递回点D;

(5)波通过更长时间会传递到钢柱、脊索处，能量也会在此处进行耗散，并将未耗散完的能量以波的形式传回点D,由此会产生干涉现象；

(6)膜面具有预应力，从而在冰雹冲击达到最大位移后，会给予冰雹一定初速度将其“弹”出膜面。

在LS-PREPOST后处理中，读取数值模拟中冰雹粒径为1.7 cm, 2.5 cm, 3.0 cm, 4.5 cm和6.0 cm共计5种粒径冲击伞形膜面测点D的位移时程，以及粒径3 cm冰雹冲击伞形膜面点D时A～F的位移时程，时程曲线如图12和图13所示。

如图12和图13所示，粒径越大冰雹冲击产生的振幅也越大，距离冲击点越近振幅也越大。而距离冲击测点由近及远不同测点在达到最大振幅的时间有所滞后，测点A～F达到最大振幅的时间分别为8.2 ms, 7.9 ms, 7.1 ms, 6.4 ms, 9.8 ms, 9.7 ms, 这种滞后现象也符合波在膜面传递的规律。而波在传递的时候有干涉现象发生，故上图位移时程图在衰减时偶有振幅增大现象出现。最大振幅、速度具体数值如表6和表7所示。

为更方便观测不同粒径冰雹冲击膜面数值模拟动力响应变化规律，如图14所示。

由粒径为1.7 cm冰雹冲击膜面点D位移时程可知，从冰雹接触膜面时刻为9 ms, 冰雹冲击达到最大位移时刻为11.5 ms, 由此可知冰雹冲击膜面达到最大位移作用时间为2.5 ms。根据动量定理可得到冰雹冲击伞形膜结构上的集中冲击荷载公式F=F(τ)/A,这样，计算得出不同粒径冰雹冲击伞形膜面作用力和冲击应力，具体计算结果如表8所示。

变化趋势如图15所示。

数值模拟中主要进行了对不同粒径的冰雹冲击伞形膜结构的模拟，并得出了不同测点的动力响应。为了检验试验的准确性，对比分析中主要分析冰雹冲击膜面测点的最大位移、速度。为方便观测对比分析，如图16所示。

由图16可知，在冰球冲击荷载作用下，试验、数值模拟的位移、速度变化趋势都呈增长趋势，即随着粒径的增大而冲击产生的位移、速度也随着增大。在粒径为1.7 cm, 2.5 cm, 3.0 cm冰雹冲击中，试验与数值模数据十分接近，从而验证了冰雹冲击伞形膜结构数值模拟的准确性，反映出了动力响应规律。然而粒径为4.5 cm, 6.0 cm冰雹冲击产生的位移以及速度出现了较大的误差，而且粒径越大，误差越明显。

试验结果与数值模拟间会相应的出现误差，比如上文所说的粒径4.5 cm, 6.0 cm冰雹冲击试验误差。但是这种误差是由多方面引起，现归纳如下几点：

(1)试验中使用的膜材为聚酯纤维织物为基材表面涂覆聚合物涂层的织物膜材，基材织物由平织方法织成，涂层覆盖于基材表面，而且膜材自身为正交异性材料，在经向和纬向上弹性模量有细微差别；在数值模拟中，ANSYS与LS-DYNA都将膜材视为弹性材料，因此ANSYS中设置的是各向同性材料。因此在膜材材料设置上的不同会使得试验和数值模拟产生误差。

(2)试验是在自然环境下进行，会受到外界各种干扰，特别是冰雹做自由落体运动时受到风、热、空气阻力、空气阻曳力等影响，冲击膜面时的末速度与理论末速度以及冲击点都有微小偏差，而数值模拟中所有工况冲击试验都是在理想条件下进行的。

(3)试验膜材由两片膜材进行热合而成，在热合处材料特性有一定的变化，而数值模拟中膜面相应位置加入了脊索以应对试验膜材热合连接处的脊索效应，因此材料上的不同会使得两者产生误差。

(4)试验中，为方便边索拉力控制和监测，在膜角点和膜张拉点之间设计了拉力传感器和角点连接件，在进行冲击试验时，振动波会沿着膜面将能量传递到膜角点并在膜角点耗散；在数值模拟中将膜面四个角点设为刚性连接，未考虑角点连接件对膜面振动波耗散的影响，而波在传递到角点时会形成反射波，因此数值模拟中较试验中干涉现象会更为严重。

(5)膜材在裁剪、安装过程中不可避免地有误差产生，而且在试验过程中操作误差也会引起偏差。

现场试验中主要采用不同粒径的冰雹冲击伞形膜面不同测点，采集冲击膜面后动力响应数据，之后进行数值模拟分析以得到相应的动力响应数据并进行对比分析，通过分析膜面测点的相关数据可以得出以下主要结论：

(1)冰雹冲击膜面测点时，膜面会产生较大位移，且位移会随着冰雹粒径的增大而增大，会随着预应力的降低而增大。

(2)冰雹冲击荷载会加速伞形膜结构的松弛，实际工程中松弛后，再遇风雨等其他荷载时会加速膜面振动，使膜面产生更大的变形，甚至导致事故发生。

(3)不同粒径冰雹冲击伞形膜面的冲击力不同，粒径为1.7 cm冰雹冲击膜面作用力为11.54 N,而粒径为6.0 cm冰雹冲击膜面作用力高达669.08 N。

(4)通过试验结果与数值模拟结果对比分析，两者误差较小，因此利用有限元分析软件ANSYS与LS-DYNA能够很好的展现冰雹冲击伞形膜结构全过程。

因此，在膜结构设计中，建议将冰雹作为偶然荷载考虑在荷载不利组合中。同时适当增大膜面预应力可提供结构抵抗变形的能力。在膜结构的日常使用阶段，结构受到冰雹等冲击荷载作用后，需及时对膜结构进行松弛检查和二次张拉维护。