充气膜结构依靠内部充入的具有一定压力的气体, 使得膜面张紧产生预张力来抵抗外部荷载效应的作用, 常分为气承式充气膜结构 (air supported membrane structure) 和气胀式充气膜结构 (air inflated membrane structure) 。气承式充气膜结构为单层膜结构, 需要将整个建筑结构封闭起来并不断的充入气体来维持内部的压力;1988年建成的东京穹顶是这类结构的典型代表[1]。气胀式充气膜结构大多为双层或多层气囊式结构, 与气承式充气膜结构相比, 这类结构的供压气体仅需充入由上、下层膜面组成的气囊之中, 能够提供较为开阔的空间;典型的建筑如1992年巴塞罗那世博会的德国馆等[2]。近年来, 随着ETFE膜材在建筑中的百富策略白菜网, 气胀式充气膜结构多以双层或多层气枕结构的形式出现;采用拼接的方式将各个气枕固定于建筑上;英国的伊甸园工程, 德国的慕尼黑足球场以及我国的国家游泳中心均采用的是这种新近出现的气枕结构[3]。

充气膜结构的分析设计必须首先进行找形分析。目前常用的数值计算方法是非线性有限单元法[4,5,6], 国内部分研究人员利用有限元软件ANSYS对ETFE气枕结构进行了初始形态分析[7,8]并给出了相应的建议。

为确定充气膜结构找形分析时, 在内压已知情况下膜面张力的大小, 根据泰勒展开公式[9]推导了充气膜结构充气膨胀的体积微分方程, 并由基于U.L.格式的非线性有限单元法[10], 建立了以节点位移和膜面初始张力为未知量的平衡方程, 编制了在给定内部气压及膨胀体积情况下的找形分析程序。文中给出了球面和圆柱面充气膜结构矢高、内压、初始张力相互关系, 以用于检验提出方法的精确性。

对于采用三角形常应变膜单元进行网格划分的充气式膜结构, 在内部气压的作用下膨胀至空间某一平衡位置的体积V可表示如下:

V=∑m(A⋅h)m (1)V=∑m(A⋅h)m(1)

其中, A为第m个单元的平面投影面积, h为第m个单元的高度。

根据海伦公式, 膜单元平面投影面积可由式 (2) ～ (6) 计算得到:

A=p⋅(p−a)⋅(p−b)⋅(p−c)−−−−−−−−−−−−−−−−−−−−−√ (2)A=p⋅(p-a)⋅(p-b)⋅(p-c)(2)

p= (a+b+c) /2 (3)

a=(xj−xi)2+(yj−yi)2−−−−−−−−−−−−−−−−−√ (4)a=(xj-xi)2+(yj-yi)2(4)

b=(xk−xj)2+(yk−yj)2−−−−−−−−−−−−−−−−−−√ (5)b=(xk-xj)2+(yk-yj)2(5)

c=(xi−xk)2+(yi−yk)2−−−−−−−−−−−−−−−−−−√ (6)c=(xi-xk)2+(yi-yk)2(6)

其中, xn, yn (n=i, j, k) 为膜单元三个节点投影面上的坐标。

设zn (n=i, j, k) 为节点高度方向坐标, 三角形膜单元m充气膨胀的高度可按下式计算:

h= (zi+zj+zk) /3 (7)

根据多元函数泰勒展开公式, 充气膜结构膨胀后的体积增量可表示为:

ΔV=∑m(∂V∂xi⋅Δxi+∂V∂yi⋅Δyi+∂V∂zi⋅Δzi+ ∂V∂xj⋅Δxj+∂V∂yj⋅Δyj+∂V∂zj⋅Δzj+ ∂V∂xk⋅Δxk+∂V∂yk⋅Δyk+∂V∂zk⋅Δzk)m (8)ΔV=∑m(∂V∂xi⋅Δxi+∂V∂yi⋅Δyi+∂V∂zi⋅Δzi+∂V∂xj⋅Δxj+∂V∂yj⋅Δyj+∂V∂zj⋅Δzj+∂V∂xk⋅Δxk+∂V∂yk⋅Δyk+∂V∂zk⋅Δzk)m(8)

即

{∂V∂X}T⋅{ΔU}=ΔV (9){∂V∂X}Τ⋅{ΔU}=ΔV(9)

其中, ∂V/∂X为体积微分向量, X为充气膜结构节点坐标向量, ΔU为节点位移向量。

充气膜结构在进行找形分析时仅需考虑几何非线性, 文中采用基于U.L.格式的非线性有限单元法, 由虚功原理得到局部坐标系下膜单元在内部气压作用下的平衡方程为:

[KNL]e·{ΔU}e={P}e-{F}e (10)

式中,

[KNL]e=A·t·[G]T·[M]·[G] (11)

{P}e=∫ANTqdA (12)

{F}e=A·t·[BL]T·{σ} (13)

其中, G为非线性相关矩阵;N为形函数矩阵;q为单位面积气压荷载向量;BL为线性应变位移关系矩阵;t为膜材厚度。

式中, 膜单元应力矩阵M为:

M=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢σx00τxy000σx00τxy000σx00τxyτxy00σy000τxy00σy000τxy00σy⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥ (14)Μ=[σx00τxy000σx00τxy000σx00τxyτxy00σy000τxy00σy000τxy00σy](14)

膜单元应力向量为:

{σ}={σxσyτxy}T (15)

在膜结构的形态分析中, 常采用等张力曲面, 即σx =σy =σ, τxy=0;根据坐标转换矩阵及式 (10) 求得膜单元在整体坐标系下的平衡方程后, 集成得到充气膜结构在整体坐标下的平衡方程为:

[[KNL′]−{P}]⋅{{ΔU} 1σ}=−{F′} (16)[[ΚΝL′]-{Ρ}]⋅{{ΔU}1σ}=-{F′}(16)

式中,

[KNL′]=∑m(A⋅t⋅[G]T⋅[G])m (17)[ΚΝL′]=∑m(A⋅t⋅[G]Τ⋅[G])m(17)

{F′}e=∑m(A⋅t⋅[BL]T⋅⎧⎩⎨⎪⎪110⎫⎭⎬⎪⎪)m (18){F′}e=∑m(A⋅t⋅[BL]Τ⋅{110})m(18)

式中, KNL′为结构总体刚度矩阵;△U为节点位移向量;P为气压作用等效节点荷载向量;F′为等效节点位移相关向量。

式 (16) 为给定找形气压及膨胀体积, 以节点位移和初始张力为待求未知量的平衡方程。由于将初始张力作为未知量, 可用于求解的平衡方程个数较未知量个数相差一个;将式 (9) 引入式 (16) 中以解决该问题, 得到新的平衡方程如下:

[[KNL′]n×n{∂V∂X}Tn−{P}n0]⋅{{ΔU}n 1/σ}={−{F′}n ΔV} (19)[[ΚΝL′]n×n-{Ρ}n{∂V∂X}nΤ0]⋅{{ΔU}n1/σ}={-{F′}nΔV}(19)

其中, 下标n为充气膜结构节点位移个数。

平衡方程式 (19) 的求解, 采用迭代的方法进行找形计算。

在找形分析过程中充气膜结构仅承受内部气压的作用, 且不考虑膜材本构关系的影响。气囊形式或气枕形式的充气膜结构由于上、下层结构对称, 找形分析时可以仅对上层结构进行找形, 即简化为气承式充气膜结构处理。本文的两个算例均采用气承式结构并从平面状态开始找形。

算例1 充气膜结构平面形状为圆形, 其直径为3 m, 周边简支固定, 膜材厚度为0.2×10-3 m, 单元网格划分及节点标号如图1所示。结构内部压力值为200 Pa, 给定充气膨胀后的体积为1.031 m3。找形得到的结果见图2, 膜面矢高为0.2905 m, 膜面的初始张力T=σ×t为0.401×103 N/m。

对于圆形边界, 在内压作用下膜面上张力各处相等的曲面为球面, 其矢高h、内压p及初始张力T有关系如下:

h=(2T/p)−(2T/p)2−r2−−−−−−−−−−√ (20)h=(2Τ/p)-(2Τ/p)2-r2(20)

其中, r为圆形边界的半径。得到矢高之后可由几何关系计算得到曲面坐标的理论值。

由式 (20) 可以得到算例1圆形边界充气膜结构矢高h=0.2911 m, 与数值解非常接近。表1为结点竖向坐标的数值解和理论解, 可以发现数值结果有很好的计算精度。

算例2 平面形状为矩形的充气膜结构长3 m, 宽1.8 m, 沿长度方向简支固定, 沿宽度方向仅可产生竖向位移, 膜材厚度为0.2×10-3 m, 单元网格划分及节点标号如图3所示。结构内部压力值为200 Pa, 给定充气膨胀后的体积为0.604 m3。找形得到的结果见图4, 膜面矢高为0.1680 m, 膜面的初始张力T=σ×t为0.501×103 N/m。

对于算例2的长方形边界, 在内压作用下膜面上张力处处相等的曲面为圆柱面, 其矢高h、内压p及初始张力T有关系如下:

h=(T/p)−(T/p)2−r2−−−−−−−−−√ (21)h=(Τ/p)-(Τ/p)2-r2(21)

其中, 2r为圆柱面平面投影的宽度。得到矢高之后可由几何关系计算得到圆柱面节点坐标的理论值。

由式 (21) 计算得到算例2长方形边界充气膜结构矢高h=0.1673 m, 可以发现求得的数值解与理论解十分接近。表2为节点竖向坐标的数值解和理论解, 通过比较可以看出数值方法计算得到的结果与理论解相差很小。

1) 根据泰勒展开公式推导了充气膜结构充气膨胀的体积微分方程, 基于U.L.格式的非线性有限单元法, 提出了以节点位移和膜面初始张力为未知量的充气膜结构找形分析的体积微分法;

2) 编制了充气膜结构在给定内部气压及膨胀体积情况下的找形分析程序, 给出了球面及圆柱面的矢高、内压、初始张力相互关系以用于校核计算结果;

3) 通过对球面及圆柱面气承式充气膜结构进行找形分析, 将计算结果与理论解进行比较, 验证了本文提出的基于有限元的体积微分法的有效性。