伞形膜结构在冰雹荷载作用下的非线性动力响应

张作亮，刘长江，郑周练，刘坚，张梦佳; ZHANG Zuoliang; LIU Changjiang; ZHENG Zhoulian; LIU Jian; ZHANG Mengjia

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

伞形膜结构在冰雹荷载作用下的非线性动力响应 PDF

- ORCID：
张作亮 ^1,2
✉
- ORCID：
刘长江 ¹
✉
- ORCID：
郑周练 ³
✉
- ORCID：
刘坚 ¹
- ORCID：
张梦佳 ¹

1. 广州大学土木工程学院，广州 510006； 2. 国网四川省电力公司南充供电公司，四川南充 637000； 3. 重庆大学土木工程学院，重庆 400045

中图分类号： O327

最近更新：2023-04-17

DOI：10.11835/j.issn.2096-6717.2021.117

摘要

对伞形张拉膜结构在自然灾害等外荷载冲击下的振动反应进行研究，为实际工程膜结构建筑受灾破坏分析提供依据。基于冯·卡门大挠度理论和薄板振动理论，建立伞形张拉膜结构在冰雹荷载作用下的非线性动力响应控制方程。采用贝塞尔函数及欧拉方程求得其振型函数，再运用KBM摄动法和迦辽金法对控制方程求解，得到其动力响应的近似解析解。利用Mathmatics软件对方程进行求解并对求解后的数据进行整理，通过Origin绘图软件画出膜面受迫振动后的时程位移图。采用ANSYS软件建立伞形张拉膜结构的有限元模型；再通过LS-DYNA显示动力学模块对冰雹冲击膜面的动力过程进行模拟，得到膜面动力响应的结果。同理，画出数值模拟结果的位移图，将整理后的理论与数值分析数据结果进行对比分析，结果表明，理论解析和数值分析结果吻合。

关键词

膜结构; 冰雹荷载; 非线性动力响应; 时程位移

膜结构自重轻、刚度小，常应用于大跨度空间结构，对外荷载（例如冰雹、雨、风等荷载）作用较敏感，故在外荷载的作用下破坏事故也时常发生，杨庆山等^[

1]对膜结构的破坏进行了分析。海南三亚市有着“美丽之冠”的地标性建筑膜结构曾于2005年和2013年两次受到严重破坏，破坏形式均为顶部膜材撕裂，造成巨大损失；1974年，在暴风作用下，美国斯波坎博览会馆拉索发生移动，造成膜面大面积松弛，经维护后又受到暴风雨的冲击，致使膜面撕裂；2002年，在经历两次台风袭击后，韩国济州岛体育场膜结构被严重撕裂破坏，造成巨大经济损失。膜结构对外荷载的作用比较敏感，灾害天气中，冰雹是常见的外荷载。Prein等^{[参考文献 2

百度学术}2]分析表明，冰雹直径可达十几厘米甚至更大，对地面农作物、建筑和车辆等造成破坏。冰雹荷载对膜结构产生较大振幅的振动，从而引起其松弛变形，降低失稳临界风速。当再次遭受暴雨、风和雪荷载作用时，膜面极易发生积水、积雪，造成安全隐患，或在风荷载作用下直接发生结构破坏。有必要对膜结构在冰雹动力荷载作用下的响应进行深入研究，揭示其致灾机理和科学规律，以防止或减少工程事故的发生。

膜结构的解析计算理论起源于板壳解析计算理论。Liu等^[

3]利用有限元分析软件模拟出冰雹冲击马鞍型膜面后的动力响应，并分析膜面在外荷载作用下的松弛情况。Sun等^{[参考文献 4

百度学术}4]基于风洞实验，给出了椭圆形支撑膜结构的风压分布特征。Kouchakzadeh等^{[参考文献 5

百度学术}5]对两种连接交叉层叠锥壳的振动响应进行了研究，采用唐纳德和哈密顿原理中的薄壳理论得到其基本控制方程。随后，Bagheri等^{[参考文献 6

百度学术}6]分析了由两个圆锥壳组成的连接壳整体的振动，采用唐纳德理论和哈密顿原理建立了一般运动方程和相关边界条件，然后进行求解。Nejati等^{[参考文献 7

百度学术}7]研究了碳纳米管（CNT）增强旋转功能梯度（FG）截头壳在厚度方向上的静态和自由振动行为。Liu等^{[参考文献 8

百度学术}8]运用冯·卡门大挠度理论和达朗贝尔原理建立了矩形张拉膜的自由振动偏微分方程，基于摄动法和伽辽金算法对受迫振动控制方程进行求解，最后通过数值模拟进行验证和分析。Colliers等^{[参考文献 9

百度学术}9]提出了一种集成压力传感器双曲薄壳风洞模型的原型化方法，并在薄壳壳面上得到了更真实的风压分布，便于计算风压对薄壳结构的动力影响。膜结构在冰雹荷载作用下的研究较少，Anghileri等^{[参考文献 10

百度学术}10]利用LSTC-LS-Dyna建立了3种冰雹数值模型：光滑粒子流体动力学模型、任意拉格朗日欧拉模型和有限元模型。Kamarian等^{[参考文献 11

百度学术}11]对碳纳米管增强复合材料锥壳的自由振动进行了研究分析，在等效FIBer假设的基础上，采用EShelbymori-Tanaka方法对纳米复合圆锥壳的材料性能进行了估算。Li等^{[参考文献 12

百度学术}12]对矩形预应力张拉膜在均匀冲击荷载作用下的动力响应进行了理论和试验研究；在薄板理论的基础上推导了横向自由振动的运动方程，结合伽辽金等方法求解，得到了自由振动的频率、速度、振幅和加速度等运动力学参数的解析解。Zhang等^{[参考文献 13

百度学术}13]对钢/短玻璃纤维增强热塑性复合板进行了低速冲击试验，对冲击响应进行研究和分析，提出了一种基于复合材料动静态试验数据的有限元模型。Zhu等^{[参考文献 14

百度学术}14]提出了一种基于刚性—完全塑性法的简单公式，以检验刚性尖锐物在任意位置对矩形板的动态响应。通过查阅文献发现，膜结构的研究通常是以矩形膜为对象，缺少对伞形膜的研究。矩形膜的研究主要是从理论上分析荷载对膜结构的冲击响应，在膜材设计上运用较多，建立的方程主要基于二维空间，而在实际工程中，膜结构建筑通常带有翘曲，常见为伞形膜结构。考虑工程实际，基于膜建筑的三维空间结构推导出伞形膜结构在冰雹荷载下的动力响应方程，主要用于对工程中翘曲膜面的分析和指导。

利用冯·卡门大挠度理论，结合达朗贝尔原理和薄板振动理论，建立伞形张拉膜结构在冰雹荷载作用下的非线性动力响应控制方程；运用KBM摄动法和迦辽金法求解伞形膜结构动力响应的近似解析解，并进行相应的数值分析，进一步验证近似解析解的合理性。

1 解析理论

1.1 控制方程和边界条件

目前对膜结构的理论分析主要基于二维结构，通过笛卡尔坐标系推导出动力响应方程，对其进一步深化并结合工程实际情况对三维结构即翘曲模型进行理论推导。翘曲膜面控制方程的建立较矩形膜面更复杂，选用极坐标建立控制方程，从而进行翘曲膜面在冲击荷载作用下的动力响应分析和理论推导。

周边固支的伞形薄膜厚为h，高为H，半径为r，膜面预张力分别为N_s和N_θ，如图1所示。

（a）

（b）

图1 伞形薄膜立面图和俯视图

Fig. 1 Elevation and top view of umbrella film

由板壳振动理论可知，伞形薄膜位移边界条件为

\{\begin{matrix} w |_{s = a} = 0 \\ |w| < \infty \end{matrix}

，

\frac{\partial w}{\partial s} = 0

（1）

式中：w为伞形薄膜发生振动时的横向位移，即挠度；s为伞形薄膜斜边上的点到顶点的距离。膜结构的基本方程包含弹性曲面微分方程和平衡微分方程及协调方程。通过弹性力学基本理论得到膜结构的协调方程为

s \frac{d^{2}}{d s^{2}} (s N_{s}) + \frac{d (r N_{θ})}{d r} - λ^{2} N_{s} = - E_{θ} h (\frac{\partial w}{\partial s} + \frac{1}{2} {(\frac{\partial w}{\partial s})}^{2})

，

λ^{2} = \frac{E_{θ}}{E_{s}} = \frac{μ_{θ}}{μ_{s}}

（2）

式中：E_s、E_θ分别为伞形薄膜径向和环向的弹性模量；μ_s、μ_θ为径向和环向的泊松比；N_s和N_θ为径向和环向的应力，伞形薄膜的应力函数设为

\{\begin{matrix} N_{s} = h \frac{1}{s} \frac{\partial ϕ}{s} \\ N_{θ} = h \frac{\partial^{2} ϕ}{\partial s^{2}} \end{matrix}

（3）

将应力函数式（3）代入伞形膜结构的协调方程式（2），得

s \frac{d^{3} ϕ}{d s^{3}} + \frac{d^{2} ϕ}{d s^{2}} - λ^{2} \frac{1}{s} \frac{d ϕ}{d s} = - E_{θ} (\frac{\partial w (s, t)}{\partial s} + \frac{1}{2} {(\frac{\partial w (s, t)}{\partial s})}^{2})

（4）

由于伞形膜结构的受力分布对称，因此，将运动微分方程进行化简，得

\frac{c o t γ}{\partial s} + \frac{\partial^{2} w}{\partial s^{2}} - \frac{1}{s} \frac{\partial w}{\partial s} = ρ \frac{\partial^{2} w}{N \partial t^{2}}

（5）

伞形膜结构的振型函数为W=Rcos （nθ），通过贝塞尔函数求解其参数，得到振型解为

\begin{array}{l} W = (\frac{c x^{2}}{a^{2}} + \frac{d x^{4}}{a^{4}} + \frac{e x^{6}}{a^{6}} + \frac{f x^{8}}{a^{8}} + \frac{g x^{10}}{a^{10}} + \frac{k x^{12}}{a^{12}}) c o s (n θ) \end{array}

（6）

采用欧拉方程对上述高阶变系数微分方程进行降次，转化为齐次线性微分方程，求解得到应力函数为

\{\begin{array}{l} \frac{d ϕ}{d s} = E_{θ} X c o s (n θ) T^{2} + E_{θ} Y c o s (n θ) T \\ \frac{d^{2} ϕ}{d s^{2}} = E_{θ} X' c o s (n θ) T^{2} + E_{θ} Y' c o s (n θ) T \end{array}

（7）

式中：X与Y是关于r的函数。

利用达朗贝尔原理和大挠度理论列出伞形薄膜在阻尼冲击荷载下的运动方程

\begin{array}{l} ρ \frac{\partial^{2} w}{\partial t^{2}} + c \frac{\partial w}{\partial t} + \frac{1}{s} \frac{\partial w}{\partial s} (N_{θ} + N_{0 θ}) - \\ \frac{1}{s} (N_{θ} + N_{0 θ}) - \frac{\partial^{2} w}{\partial s^{2}} (N_{s} + N_{0 s}) = \\ F (t) δ (s - s_{0}) \end{array}

（8）

式中：N_0θ为伞形膜面施加的初始预张力； $δ (s - s_{0})$ 为狄拉克函数；ρ为薄膜的材料密度；c为薄膜的阻尼；F（t ）为冰雹冲击膜面时的冲击力。由薄板振动理论可知挠度w

w = W T (t) c o s (n θ)

（9）

式中：挠度w是关于时间t的函数，将式（3）、式（7）、式（9）带入式（8）中化简，求得

\begin{array}{l} ρ W T ″ + c W T' + \frac{1}{s} W' T (\frac{d^{2} ϕ}{d s^{2}} + N_{0 θ}) - \\ \frac{1}{s} (\frac{d^{2} ϕ}{d s^{2}} + N_{0 θ}) - W ″ Τ (\frac{1}{s} \frac{d ϕ}{d s} + Ν_{0 s}) = \\ F (t) δ (s - s_{0}) \end{array}

(10)

1.2 控制方程求解

由于方程式（10）非常复杂，难以求解其精确解，摄动法是求解非线性方程近似解的有效方法。采用KBM摄动法对控制方程组进行近似求解。首先运用伽辽金法分别对等式两边进行积分，得

\begin{array}{l} \int_{0}^{a} (ρ W^{2} T^{″} + c W T^{'} + \frac{1}{s} W^{'} W T^{3} E_{θ} X^{'} + \\ \frac{1}{s} W^{'} W T^{2} E_{θ} Y^{'} + \frac{1}{s} W^{'} W T N_{0 θ} - \frac{W}{s} E_{θ} X^{'} T^{2} - \\ \frac{W}{s} E_{θ} Y^{'} T + \frac{W}{s} N_{0 θ} - W^{″} W T^{3} \frac{1}{s} E_{θ} X - \\ W^{″} W T^{2} \frac{1}{s} E_{θ} Y - W^{″} W T N_{0 θ}) d s = \\ \int_{0}^{a} F (t) δ (s - s_{0}) d s W (r_{0}) \end{array}

（11）

合并同类项为

\begin{array}{l} \int_{0}^{a} [ρ W^{2} T^{″} + c W T^{'} + \\ (\frac{1}{s} W^{'} W N_{0 θ} - \frac{W}{s} E_{θ} Y^{'} - W^{″} W N_{0 θ}) T - \\ (\frac{W}{s} E_{θ} X^{'} + W^{″} W \frac{1}{s} E_{θ} Y) T^{2} - \\ W^{″} W T^{3} \frac{1}{s} E_{θ} X + \frac{W}{s} N_{0 θ}] d s = \\ \int_{0}^{a} F (t) δ (s - s_{0}) d s W (r_{0}) \end{array}

（12）

将式（3）、式（6）、式（7）代入式（12）计算化简，得

\begin{array}{l} a ρ c o s {(n θ)}^{2} A_{1} T^{″} + a c c o s (n θ) A_{2} T^{'} + \\ (N c o s (n θ) A_{3} - E_{θ} c o s {(n θ)}^{2} A_{4}) T - \end{array}

\begin{array}{l} \frac{1}{a^{2}} (A_{5} E_{θ} c o s {(n θ)}^{2} + q c o s {(n θ)}^{3} E_{θ} A_{6}) T^{2} - \\ \frac{c o s {(n θ)}^{3} E_{θ}}{a^{3}} A_{7} T^{3} = \int_{0}^{a} F (t) δ (s - s_{0}) d s W (r_{0}) \end{array}

（13）

将式（13）化简为二阶变系数微分方程

\begin{array}{l} T^{″} + \frac{c A_{2} T^{'}}{ρ c o s {(n θ)}^{2} A_{1}} + (\frac{N A_{3}}{a ρ c o s (n θ) A_{1}} - \frac{E_{θ} A_{4}}{a ρ A_{1}}) T - \\ (\frac{A_{5} E_{θ} + c o s {(n θ)}^{3} E_{θ} A_{6}}{q a^{3} ρ A_{1}}) T^{2} - \frac{c o s (n θ) E_{θ}}{a^{4} ρ A_{1}} A_{7} T^{3} = \\ \frac{\int_{0}^{a} F (t) δ (s - s_{0}) d s W (r_{0})}{a ρ c o s {(n θ)}^{2} A_{1}} \end{array}

(14)

F（t）为冰雹冲击膜面时的撞击力。

F (t) = - M W (x_{0}, y_{0}) T^{″}

（15）

将式（15）代入式（14）中，化简为齐次方程式

\begin{array}{l} T ″ + \frac{a c A_{2} c o s (n θ)}{A_{1} a ρ c o s {(n θ)}^{2} + M W^{2}} T' + \\ (\frac{N A_{3} c o s (n θ)}{A_{1} a ρ c o s {(n θ)}^{2} + M W^{2}} - \frac{A_{4} E_{θ} c o s {(n θ)}^{2}}{A_{1} a ρ c o s {(n θ)}^{2} + M W^{2}}) \cdot \\ T - (\frac{A_{5} E_{θ} c o s {(n θ)}^{2} + A_{6} E_{θ} c o s {(n θ)}^{3}}{A_{1} a^{3} ρ c o s {(n θ)}^{2} + a^{2} M W^{2}}) T^{3} - \\ \frac{A_{7} E_{θ} c o s {(n θ)}^{3}}{A_{1} a^{4} ρ c o s {(n θ)}^{2} + a^{3} M W^{2}} T^{3} = 0 \end{array}

（16）

式（16）是关于时间T(t)的非线性变系数齐次微分方程，采用KBM摄动法求解该方程的近似解析解。

首先设摄动小参数ε=h/a $≪$ 1，同时令χ=χ（τ）=T（t），则

χ ″ + {ω_{0}}^{2} χ = ε (α_{1} χ^{3} + α_{2} χ^{2} + α_{3} χ')

（17）

其中

\begin{array}{l} {ω_{0}}^{2} = \frac{A_{3} N c o s (n θ)}{A_{1} a ρ c o s^{2} (n θ) + Μ W^{2} s_{0}} - \\ \frac{A_{4} E_{θ} c o s^{2} (n θ)}{A_{1} a ρ c o s^{2} (n θ) + a q M W^{2} s_{0}} \end{array}

α_{1} = \frac{A_{7} E_{θ} c o s^{3} (n θ)}{q h A_{1} a^{3} ρ c o s^{2} (n θ) + q h a^{2} M W_{s_{0}}^{2}}

α_{2} = \frac{A_{5} E_{θ} c o s^{2} (n θ) + Α_{6} E_{θ} q c o s^{3} (n θ)}{q h A_{1} a^{2} ρ c o s^{2} (n θ) + q h a M W_{s_{0}}^{2}}

α_{3} = - \frac{a^{2} c A_{2} c o s (n θ)}{h A_{1} a ρ c o s^{2} (n θ) + h Μ W_{s_{0}}^{2}}

基于KBM摄动法的基本思想，令 $f (χ, χ ″) = α_{1} χ^{3} + α_{2} χ^{2} + α_{3} χ'$ ，并设 $χ = a c o s ψ$ 。式中： $a$ 和 $ψ$ 是关于时间t的函数， $a$ 和 $ψ$ 的值由式（18）、式（19）确定。

\frac{d a}{d t} = - \frac{ε}{ω_{0}} A_{0} (a)

（18）

\frac{d ψ}{d t} = ω_{0} - \frac{ε}{a ω_{0}} C_{0} (a)

（19）

式中：A₀（a）和C₀（a）为

\begin{array}{l} A_{0} (a) = = \frac{1}{2 π} \int_{0}^{2 π} s i n (α_{1} {(a c o s ψ)}^{2} + \\ α_{2} {(α c o s ψ)}^{2} - α_{3} (a w_{0} s i n ψ) d ψ = \\ \frac{1}{2 π} \int_{0}^{2 π} s i n ψ f (a c o s ψ - a ω_{0} s i n ψ) d ψ = - \frac{1}{2} a α_{3} ω_{0} \end{array}

\begin{array}{l} C_{0} (a) = \int_{0}^{2 π} c o s (α_{1} {(a c o s ψ)}^{2} + α_{2} {(α c o s ψ)}^{2} - \\ α_{3} (a w_{0} s i n ψ) d ψ = \frac{1}{2 π} \int_{0}^{2 π} c o s ψ f (a c o s ψ - \\ a ω_{0} s i n ψ) d α = \frac{3}{8} α_{1} a^{3} \end{array}

（20）

\frac{d a}{d t} = \frac{ε α_{3} a}{2}

（21）

\frac{d ψ}{d t} = ω_{0} - \frac{ε}{a ω_{0}} C_{0} (a)

（22）

利用分离变量法求解式（21），得

a = D e^{0.5 α_{2} ε t}

（23）

同理，采用分离变量法求解式（22），得

ψ = (ω_{0} - \frac{3 α_{1} a^{2} ε}{8 ω_{0}}) t + ψ_{0}

(24)

将式（23）和式（24）代入式（17），求得

χ = χ (t) = D e^{0.5 α_{3} ε t} c o s ((ω_{0} - \frac{3 α_{1} ε D^{2} e^{α_{3} ε t}}{8 ω_{0}}) t + ψ_{0})

（25）

式中： $ψ_{0}$ 为初相位；D为振幅，其值均由初始条件确定，在（25）式中，令

ω = ω_{0} - \frac{3 α_{1} ε D^{2} e^{α_{3} ε τ}}{8 ω_{0}}

（26）

式（26）为薄膜有阻尼受迫振动频率，表明在冰雹的冲击作用下，伞形薄膜有阻尼受迫振动的频率与时间t和阻尼都有关。考虑薄膜和冰雹组成的系统，忽略能量耗散，由动量守恒

m v_{0} = m v_{0}' + ρ v_{0}' \int_{0}^{2 π} \int_{0}^{a} W s d θ d s

（27）

式中： $v_{0}'$ 为t=0时刻膜面的初始速度； $v_{0}$ 为冰雹冲击膜面时的初始速度，求得

v_{0}' = \frac{m v_{0}}{m + ρ π a^{2} (\frac{b}{2} + \frac{c}{3} + \frac{d}{4} + \frac{e}{5} + \frac{f}{6} + \frac{g}{7} + \frac{k}{8})}

（28）

假设薄膜受到冰雹的冲击，并且在t=0时刻薄膜与冰雹具有同样的初始速度v₀，因此，设初始条件为

{\frac{d w (s, t)}{d t}}_{t = 0} = W (s_{0}) {\frac{d T (t)}{d t}}_{t = 0} = v_{0}'

（29）

χ' {(t)}_{t = 0} = {\frac{d T (t)}{d t}}_{t = 0} = \frac{v_{0}^{'}}{W s_{0}}

（30）

将式（29）、式（30）代入式(25)，得

0 = D c o s ψ_{0}

（31）

因D不等于零，因此 $c o s ψ_{0} = 0$ ，则， $ψ_{0} = k π /$ 2（k=1,3,5,7…），取 $ψ_{0} = π / 2$ ，代入式(25)并对其求一阶偏导为

\begin{array}{l} \frac{d χ (t)}{d t} = \frac{D}{2} α_{2} ε e^{0.5 α_{2} ε t} c o s ((ω_{0} - \frac{3 α_{1} ε D^{2} e^{α_{2} ε t}}{8 ω_{0}}) t + \frac{π}{2}) - \\ D e^{0.5 α_{2} ε t} s i n ((ω_{0} - \frac{3 α_{1} ε D^{2} e^{α_{2} ε t}}{8 ω_{0}}) t + \frac{π}{2}) \cdot \\ (ω_{0} - \frac{3 α_{1} ε D^{2} e^{α_{2} ε t}}{8 ω_{0}} - \frac{3 α_{1} α_{2} ε^{2} D^{2} e^{α_{2} ε t}}{8 ω_{0}} t) \end{array}

（32）

将式（30）代入式（32）进行求解

D (\frac{3 α_{1} ε D^{2}}{8 ω_{0}} - ω_{0}) = \frac{v_{0}'}{W r_{0}}

（33）

求解式（33）中D的方程式，得到3个根，一个实数根及两个复数根，实数根即为膜面的振幅

D = - {\frac{2 (- \frac{1}{3} (v_{0}' ω_{0} + W s_{0} {ω_{0}}^{2}))}{(W s_{0} ε α_{1})^{\frac{1}{3}}}}^{\frac{1}{3}}

（34）

其中Ws₀为

\begin{array}{l} W s_{0} = \frac{b {s_{0}}^{2}}{a^{2}} + \frac{c {s_{0}}^{4}}{a^{4}} + \frac{d {s_{0}}^{6}}{a^{6}} + \frac{e {s_{0}}^{8}}{a^{8}} + \\ \frac{f {s_{0}}^{10}}{a^{10}} + \frac{g {s_{0}}^{12}}{a^{12}} + \frac{k {s_{0}}^{14}}{a^{14}} \end{array}

（35）

解得薄膜受迫振动的频率为

ω = ω_{0} - \frac{3 α_{1} ε D^{2} e^{α_{3} ε τ}}{8 ω_{0}}

（36）

1.3 算例分析

取伞形膜面环向和径向的弹性模量分别为E_s=1.4×10⁶ kN/m²，E_θ=0.9×10⁶ kN/m²，薄膜密度为ρ=1.7 kg/m²，薄膜厚度为h=1 mm，阻尼系数c=120 N·s/m³，膜面高度H=0.4 m，顶点到周边的距离s=8 m，膜面预张力N=10 kN/m²。设冰雹为理想球体，半径分别为1.7、2.5、3.0、4.5、6.0 cm，并假设冰雹表面圆滑，不考虑空气等因素，冰雹与膜面的接触位置距离顶点2a/3。

根据式(34)、式（25）对冰雹冲击膜面时发生的时程位移进行分析，画出前3阶的时程位移图。

将s和h代入到（3）式中，可以得到应力函数表达式，再将应力函数结果代入协调方程式（4）中，简化得到关于振型函数的方程式。将振型中W=cos（nθ）中n取1，求其一阶振型函数，将质量M、薄膜密度ρ、弹性模量E、膜面预张力N等代入w₀和α₁中，求解振幅D。最后通过Mathmatics软件对式（25）进行计算，并画出时间t=0.1 s内的一阶时程位移曲线图。同理，求出不同粒径冰雹冲击膜面后各振动阶数的时程位移，如图2所示，最后相加求得冰雹冲击膜面的总时程位移变化，如图3所示。

（a） d=1.7 cm

（b） d=2.5 cm

（c） d=3.0 cm

（d） d=4.5 cm

（e） d=6.0 cm

图2 不同粒径冰雹冲击膜面后前3阶的时程位移

Fig. 2 Third-order time-history displacement of hail impinging film with different particle sizes

（a） d=1.7 cm

（b） d=2.5 cm

（c） d=3.0 cm

（d） d=4.5 cm

（e） d=6.0 cm

图3 膜面振动总位移时程曲线

Fig. 3 Time-history curve of total displacement of membrance surface vibration

由冰雹冲击膜面后的时程位移图可知：

1）在冰雹冲击作用下，在与冰雹接触后的第1个周期内，膜面产生的位移最大，随着时间的增加，不断衰减，在t=0.1 s内，振幅急剧下降，位移趋于零，薄膜在平衡位置小幅度来回振动，当t=∞时，位移为零，膜面回到平衡状态。反映出膜结构自身的阻尼特性及能量守恒定律。

2）相同粒径冰雹冲击下，随着振动阶数增加，膜面振动频率也随之增大，位移随之减小，振动周期随时间变化，反映了薄膜结构振动的非线性。

3）随着冰雹粒径的增大，冲击速度越大，薄膜结构的振幅越大。

取直径d为1.7、2.5、3.0、4.5、6.0 cm的冰雹，代入式(32)，绘出撞击点( $s = 2 a / 3$ )的时程总位移曲线图（图3）进行分析。

从图3可以看出，膜面总位移时程呈明显非线性。振幅和周期随时间变化，膜面冲击点在初始振动时位移最大，随时间逐渐变小，当t=∞时，位移为零，膜面回到平衡位置。这是由于薄膜结构自身的阻尼特性，在振动过程中由于能量损耗而产生的变化过程。冰雹的直径大小也直接影响膜面振幅的大小，随着冰雹直径的增大，动能随之增加，导致膜面产生的位移增大，如图4所示。

图4 薄膜冲击点位移随冰雹直径的变化

Fig. 4 Displacement of film impact point varies withhail diameter

由图4可知，薄膜冲击点的位移与冰雹直径的大小成正相关，冰雹直径越大，产生的位移就越大。

2 有限元分析

2.1 材料属性

首先使用ANSYS进行隐式找形分析，膜面单元选择壳单元Shell41，边索采用杆单元Link10，冰雹采用实体单元Solid185，如表1所示。隐式分析中，采用“降温法”对膜材施加温度荷载模拟预张力，然后在LS-DYNA模块中建立冰雹冲击膜面的动力学模型。为了简化计算，定义冰雹为刚体并赋予初速度。根据隐—显式顺序求解要求，将壳单元Shell41转换为显式壳单元Shell163，杆单元Link10转换为显式杆单元LINK167，实体单元Solid185转换为显式实体单元SOLID164。最后输出k文件进行求解。

表1 伞形张拉膜材料属性

Table 1 Properties of umbrella tension film material

单元	h/mm	ρ/(g·m^-3)	E_s/MPa	E_θ/ MPa	泊松比
SHELL41	1	950	1.4×10⁶	0.9×10⁶	0.34
SOLID185		900	53	53	0.30
Link10		7 850	1.5×10⁵	53	0.30

2.2 伞形膜面动力响应分析

按照解析理论算例参数进行伞形膜结构在冰雹冲击下的数值模拟。使用LS-Prepost对求解结果进行后处理，冰雹冲击膜面的动态过程如图5所示。

（a）冲击前应力分布

（b）接触时应力分布

（c）冲击时应力集中

（d）冲向后应力开始扩散

（e）冰雹与膜面开始分离

（f）冰雹与膜面完全分离

图5 膜面冲击点的位移变化过程

Fig. 5 Process of displacement changes of impact point on membrane surface

由图5可知：

1）在受到冲击之前，膜面预张力分布较为均匀，由于边界条件因素，膜面顶部出现应力集中。

2）当冰雹与膜面撞击之后，粒径小的冰雹并不会立刻被膜面弹起，在与膜面发生相对转动之后，沿垂直膜面方向弹起。

3）当冰雹与膜面接触时，接触面的竖向位移最大。随着进一步接触，膜面以产生弹性变形的方式吸收冰雹的动能，膜面的位移首先径向传开，当膜面冲击点位移达到最大时，冰雹的动能则完全转化为膜面的弹性势能。

4）膜面位移最大时，弹性势能最大，其中一部分能量又转化为冰雹的动能，将冰雹以一定的初速度弹射出去。

5）当位移达到最大后，振动能量向膜面纬向发展。

6）振动能量纬向扩散一周后，又回到冲击点位置，影响其振动位移，因此可以得出，能量在伞形膜结构上的传递相互影响，导致膜面动力响应的非线性。

通过后处理软件LS-Prepost提取冲击点位移时程数据，绘出粒径为1.7、2.5、3.0、4.5、6.0 cm的位移时程图，如图6所示。

（a） d=1.7 cm

（b） d=2.5 cm

（c） d=3.0 cm

（d） d=4.5 cm

（e） d=6.0 cm

图6 膜面在不同粒径冰雹冲击下的时程位移

Fig. 6 Time-history displacement of the film surface under hail impacut with different particle size

由图6可知，从直径为1.7 cm的冰雹冲击膜面后的时程位移图可以看出，受到冰雹冲击后，膜面迅速达到最大位移，此时膜面具备的弹性势能最大，膜面冲击点的速度为零，随后位移开始震荡衰减。

与解析理论相似，在达到最大位移后，在阻尼的作用下，膜面冲击点位移逐渐衰减，当时间t足够大时，膜面恢复平衡状态。

膜面冲击点的最大位移随着冰雹直径的增大而增大，在直径6.0 cm时，冰雹位移最大，约为21 mm。因此，冰雹给膜结构带来的影响不可忽略。

不同直径冰雹冲击膜面的最大位移数据对比如图7所示。选取膜面测点D（即解析算例中膜面离顶点距离2a/3的点）的动力响应，对比分析膜面冲击点的最大位移、速度、振动频率的解析解和数值解，如图7所示。

图7 不同粒径冰雹冲击膜面后位移的变化趋势

Fig. 7 Variation trend of displacement after hail impact with different particle size

从图3、图4、图6、图7和表2可以看出，两种结果趋于一致，即随着冰雹直径的增大，位移也变大。在冰雹冲击膜面后速度变为零时，膜面冲击点位移最大，随后开始衰减，当时间足够大时，膜面又恢复到平衡位置。通过理论和数值模拟分析可以看出，数值模拟结果与理论结果的平均误差为8.79%，因此，认为解析理论与数值模拟结果吻合较好，可以应用于工程中。

表2 理论与数值的对比分析

Table 2 Comparison analysis of theoretical and numerical results

冰雹粒径/cm	理论分析位移/mm	数值分析位移/mm	绝对差	相对差/%
1.7	2.28	2.58	0.30	13.15
2.5	4.34	4.13	0.21	4.83
3.0	6.19	6.69	0.50	8.07
4.5	14.03	14.96	0.93	6.62
6.0	19.25	21.42	2.17	11.27

将四面拼接的伞形膜面在冰雹冲击荷载下的试验数据^[

15]与张拉膜的理论及数值模拟结果进行对比分析。试验数据见表3。

表3 不同粒径冰球冲击膜面的试验数据分析

Table 3 Experimental data analysis of ice hockey with different sizes impacting film surface

冰球粒径/cm	位移峰值/mm
1.7	1.99
2.5	3.54
3.0	5.45
4.5	11.03
6.0	16.81

通过数据对比可以看出，相对试验数据，张拉膜得到的位移偏大，这是由于试验中采用的模型为四面拼接而成的膜面，如图7所示，而张拉膜采用的模型表面为一整块且底边为圆形，相较于四面伞形膜面，其受到的张拉力较小，因此，得到的数据相对偏大，但与试验得出的数据变化趋势一致。并且，理论和数值模拟结果与试验结果的相对差较小，因此，可以得出，本文的数值模拟和理论分析结果是合理的。

图8 四面拼接的伞形膜面

Fig. 8 Umbrella shaped membrane surface with four sides splicing

3 结论

采用解析和数值方法研究伞形膜结构在冰雹冲击荷载作用下的动力响应问题，主要结论如下：

1）由于柔性膜材无法承受弯矩，在冰雹冲击作用下，膜面通过产生形变来抵抗冰雹的冲击荷载，导致冲击点位置产生较大的竖向位移。二者的冲击过程是冰雹的动能转化为膜材的弹性势能的过程。

2）在冰雹冲击问题中，通常认为冰雹与接触面的作用时间为0.001 s。而当冰雹冲击柔性膜面时，二者相互作用时间随着冰雹直径的增大而增大。这是因为较大粒径的冰雹具有更大的质量与末速度，在冲击的能量转换过程中，膜面需要产生更大的弹性势能来吸收冰雹的动能。因此，更长时间的能量转换过程对应的接触时间变长。

3）在膜面冲击点达到最大位移后，由于阻尼等影响开始衰减，当时间t足够大时，膜面恢复到平衡位置。

4）通过不同粒径冰雹冲击膜面后的对比可得出，冰雹直径越大，膜面冲击点产生的位移也越大。本次试验中最大粒径冰雹冲击膜面后的位移约为20 mm，当冰雹比较大时，膜面将产生较大位移，因此，冰雹等灾害对于薄膜建筑的影响是不可忽视的。

5）通过理论和数值分析结果的对比可以发现，其位移相对差较小，理论和数值模拟结果在变化规律上趋于一致，通过数值模拟分析冰雹冲击膜面后的振动响应可行。

参考文献

杨庆山, 姜忆南. 张拉索—膜结构分析与设计[M]. 北京: 科学出版社, 2004. [百度学术]

YANG Q S, JIANG Y N. Analysis and design of tensioned cable-membrane structures [M]. Beijing: Science Press, 2004. (in Chinese) [百度学术]

PREIN A F, HOLLAND G J. Global estimates of damaging hail hazard [J]. Weather and Climate Extremes, 2018, 22: 10-23. [百度学术]

LIU C J, YANG S P, LIU J, et al. Dynamic response of saddle membrane structure under hail impact [J]. Engineering Structures, 2020, 214: 110597. [百度学术]

SUN X Y, ARJUN K, WU Y. Investigation on wind tunnel experiment of oval-shaped arch-supported membrane structures [J]. Journal of Wind Engineering and Industrial Aerodynamics, 2020, 206: 104371. [百度学术]

KOUCHAKZADEH M A, SHAKOURI M. Free vibration analysis of joined cross-ply laminated conical shells [J]. International Journal of Mechanical Sciences, 2014, 78: 118-125. [百度学术]

BAGHERI H, KIANI Y, ESLAMI M R. Free vibration of joined conical-conical shells [J]. Thin-Walled Structures, 2017, 120: 446-457. [百度学术]

NEJATI M, ASANJARANI A, DIMITRI R, et al. Static and free vibration analysis of functionally graded conical shells reinforced by carbon nanotubes [J]. International Journal of Mechanical Sciences, 2017, 130: 383-398. [百度学术]

LIU C J, DENG X W, ZHENG Z L. Nonlinear wind-induced aerodynamic stability of orthotropic saddle membrane structures [J]. Journal of Wind Engineering and Industrial Aerodynamics, 2017, 164: 119-127. [百度学术]

COLLIERS J, MOLLAERT M, DEGROOTE J, et al. Prototyping of thin shell wind tunnel models to facilitate experimental wind load analysis on curved canopy structures [J]. Journal of Wind Engineering and Industrial Aerodynamics, 2019, 188: 308-322. [百度学术]

ANGHILERI M, CASTELLETTI L M L, INVERNIZZI F, et al. A survey of numerical models for hail impact analysis using explicit finite element codes [J]. International Journal of Impact Engineering, 2005, 31(8): 929-944. [百度学术]

KAMARIAN S, SALIM M, DIMITRI R, et al. Free vibration analysis of conical shells reinforced with agglomerated Carbon Nanotubes [J]. International Journal of Mechanical Sciences, 2016, 108/109: 157-165. [百度学术]

LI D, ZHENG Z L, TIAN Y, et al. Stochastic nonlinear vibration and reliability of orthotropic membrane structure under impact load [J]. Thin-Walled Structures, 2017, 119: 247-255. [百度学术]

ZHANG Y B, SUN L Y, LI L J, et al. Experimental and numerical investigations on low-velocity impact response of high strength steel/composite hybrid plate [J]. International Journal of Impact Engineering, 2019, 123: 1-13. [百度学术]

ZHU L, SHI S Y, JONES N. Dynamic response of stiffened plates under repeated impacts [J]. International Journal of Impact Engineering, 2018, 117: 113-122. [百度学术]

LIU C J, WANG F, DENG X W, et al. Hailstone-induced dynamic responses of pretensioned umbrella membrane structure [J]. Advances in Structural Engineering, 2021, 24(1): 3-16. [百度学术]