这几年薄膜结构在我国的发展势头极为迅猛, 特别是2008年北京奥运会的两个主要建筑物:主体育场“鸟巢”和游泳中心“水立方”, 均采用了薄膜结构, 这更将促进薄膜结构在我国的进一步发展。

由于膜材是只能抗拉不能抗压的柔性材料, 所以膜材的受力状态可以分为纯拉状态、单向褶皱状态 (又称为褶皱状态) 和双向褶皱状态 (又称为松弛状态) 三种。薄膜结构的刚度是由膜材平面内的拉力形成的, 当膜材出现褶皱时, 其受力性能与纯拉状态不同。因此, 为了正确的对薄膜结构进行受力分析, 必须确定膜材的受力状态并采取正确的处理方法。

最早关于褶皱分析的方法是Wagner分析工字型及箱型梁薄腹板所能承担的最大剪力问题时提出的张力场理论 (Tension Field Theory) , Mansfield在其基础上用一个合理的松弛能量密度代替原有的应变能将该理论进行了改进, 后来Pipkin又将其用于各向同性膜材的褶皱分析中[1]。1961年, Stein和Hedgepeth提出了Stein-Hedgepeth理论, 主要用来分析薄膜的局部褶皱, Kang和Adler等在其基础上做了一些改进[2], 使其适用范围更加广泛。另外, Roddeman提出了针对薄膜褶皱的修正变形张量方法, 由于该方法公式复杂, 无法得出薄膜内力和切线刚度矩阵的明确表达式, 后来K.Lu等采用了曲线坐标得到了薄膜内力和切线刚度矩阵的简明表达式[3]。当需要描述褶皱的具体形状, 即褶皱的幅值和波长时, 可以采用采用分叉理论[4]。但是, 以上方法的研究对象一般为卫星天线、雷达、气囊、降落伞等。Fujikake M.等针对张拉薄膜结构采用了修正本构矩阵的方法来研究褶皱问题, 并在ADINA程序中得以实现[5], 但其在列应力与应变关系式时, 将预应力单独写出, 即没有考虑初始预应力所产生的应变, 本文作者认为这与真实情况不符, 不能得到真实的应变, 从而不能正确的判别出膜材的受力状态。

本文作者认为在研究张拉薄膜结构褶皱问题时, 主要关心的是发生褶皱的区域和发生褶皱后结构的受力性能, 而褶皱的幅值和波长并不重要, 并且修正本构矩阵的方法思路清晰, 公式简单, 便于程序实现, 若是还能达到所需的精度应该是一个值得采用的方法。故本文作者在以前工作的基础上, 采用主应力-主应变判别准则和更为合理的修正本构矩阵方法对褶皱问题进行分析, 编写了有限元程序, 并进行了验证和检验。

膜材有纯拉、单向褶皱和双向褶皱三种受力状态, 其判别准则有主应力准则、主应变准则和主应力-主应变准则三种。主应力准则能够正确的判别出膜材是处于纯拉状态还是单向褶皱状态, 但是无法正确的判别出膜材是处于单向褶皱状态还是双向褶皱状态;主应变准则只适用于各向同性膜材, 如ETFE等;主应力-主应变准则能够正确的判别出各种膜材的三种受力状态[4,6]。

设膜材内的主应力为σ1和σ2 (σ2≤σ1) , 主应变为ε1和ε2 (ε2≤ε1) , 可将主应力-主应变准则表述如下:

(1) 当σ2>0时, 膜材处于纯拉状态;

(2) 当σ2≤0且ε1>0时, 膜材处于单向褶皱状态;

(3) 当ε1≤0时, 膜材处于双向褶皱状态。

在有限元分析中, 首先得到的是单元内某一点位移{u v w}, 然后由几何方程[7]:

得到其Green应变e={eX eY eXY}, 进而由物理方程:

即:

得到任一点的Kirchhoff应力{SX SY SXY}。D为偏轴本构矩阵, 它与膜材主轴本构矩阵D的关系为:

其中, ST为主偏轴间的转换矩阵, θ为主轴和偏轴间的夹角 (顺时针为正) , 1E为经向弹性模量, 2E为纬向弹性模量, G12为剪切模量, ν1、ν2分别为经向引起的纬向、纬向引起的经向泊松比。为了得到当前状态下的真实应力, 需要将得到的Kirchhoff应力用下式转换为Cauchy应力σ={σxσyσz}:

式中

由该点的Cauchy应力分量可计算出其主应力:

同样, 应将Green应变用下式转换为Almansi应变ε={εxεyεz}:

并进而计算出其主应变:

这样即可按主应力-主应变准则判别出膜材当前的受力状态。

若判别出某一单元处于单向褶皱或双向褶皱状态则需要对本构矩阵进行修正后再进行后续的计算。下面先给出没有褶皱时膜单元的本构关系, 再分别给出膜材发生单向褶皱和双向褶皱时的修正方法。

由于膜材一般为正交异性材料, 应力应变间存在拉剪耦合作用, 其应力主轴与应变主轴不一定相同[6], 因此主应力与其方向上的应变ε={ε1′ε′2γ1′2}间的关系为:

将其展开后可写为:

由式 (18) 可将γ12表示为:

将其代入式 (17) 得:

由σ2=0可得:

将式 (19) 和式 (21) 代入式 (16) , 可得单向褶皱时主应力σ1与ε1′间的关系:

若令:

则上式可写为:

此时膜单元在第一主应力方向上满足式 (22) , 所以本构关系在其它方向上只产生位移和应变而不产生应力, 为了方便程序编制和满足计算中的收敛要求, 仍将主应力方向上的本构关系写成矩阵形式:

其中η为一任意微小量, 以防计算中出现病态矩阵。

膜单元为双向褶皱时, 应分别令σ2=0和σ1=0, 此时柔性膜材无刚度, 同前, 为了方便程序编制和满足计算中的收敛要求, 主应力方向的本构矩阵为:

在有限元分析中, 需将上面得到的主应力向量和主应力方向上的本构矩阵转换到单元坐标中, 其转换式分别如下[7]:

其中

其中

式中ϕ为主应力方向与x坐标轴间的夹角 (顺时针为正) 。

进行上述修正后即可采用相同的方法继续下一步的荷载分析。

依据前面所说的膜材褶皱判别和处理的方法, 在本研究组开发的索膜结构分析软件CAFA1.0[8,9]的基础上编写了有限元程序, 对各种算例进行了分析。下面首先由经典算例证明程序的正确性, 再给出一些工程算例的计算结果。

图1所示的一矩形膜片, 沿AB、CD两对边作用有均布拉力Ny, 在AC、BD两对边作用有拉力P和弯矩M。当M=0时, 整个膜材处于纯拉状态, 随着M的增大, 膜材底部则会出现沿y方向的单向褶皱, 并且褶皱区域EFDC会随着M的增大而不断变大。下面先给出其理论解, 再与有限元程序的计算结果对比。

拉力P和弯矩M与x方向的应力σx之间的关系为[10]:

膜材内的应力和应变只是y的函数。设膜材为各向同性材料, 其弹性模量为E, 泊松比为v, 厚度为t。在没有褶皱的区域内x方向应力的表达式为σx=Ek (y-b) (式中k为一常数, 由边界条件可确定) , 褶皱区域 (即0<y<b) 内x方向应力的表达式为σx=0。将其代入式 (27) 中可以得到:

由式 (29) 可以得:

将上式代入式 (28) 即可以得出褶皱高度b (0≤b≤h) 与外荷载之间的关系为[11]:

在有限元分析时取E=300MPa, ν=0.3, l=10m, h=3.0m, 厚度为0.001m。为了在有限元分析中引入边界条件, 取其右侧一半进行计算, 根据对称性原则计算简图如图2所示。将O点的x和y方向位移均约束, 其余各点只约束其x方向位移, 沿y向可以自由变形。外荷载Ny、P和M都用节点荷载来等效。保持拉力P=60N不变, 共施加五种弯矩M=0N·m, M=30N·m, M=45N·m, M=60N·m, M=75N·m。计算结果和理论解的对比见图3。

由图3可以看出, 计算值与理论值极为接近, 从而证明了有限元程序的正确性。

平面投影为对角线长6 m的正方形, 四边固定, 最高点A点和C点与最低点B点和D点的高差为1m。有限元模型如图4所示, 经向和纬方向的弹性模量分别为1E=540MPa, 2E=180MPa, 剪切模量为G=20MPa, 两泊松比分别为v12=0.2, v21=0.6, 厚为0.001m, 预应力为σx=σy=0.1k N/m, 找形结果如图5所示。

对其施加均布向下的荷载0.05kN/m2, 经过分析计算可得到如图6所示的褶皱单元分布图, 图中带有+号的单元发生了单向褶皱。

图7和图8为哈尔滨工业大学校园内的一膜结构小品的平面图和立面图。该结构是由六个呈环形分布大小相同的马鞍形单元组成, 形成正六边形的图案。单个单元展开面积为42.02m, 高点距地面9m, 低点距地面5.0m, 外边边长13.6m, 内边边长5.0m。膜材为米乐1000, 取其物理参数为1E=540MPa, 2E=180MPa, 12G=12MPa, 12v=0.2, 21v=0.6, 厚为0.001m。对其一个单元进行分析计算, 预应力为2.0kN/mxyσ=σ=, 内边索预拉力为8.0 k N, 两侧边索预拉力为25.0 k N, 外边索为30.0 k N。找形结果如图9所示。

施加均布向下的荷载1kN/m2, 计算所得的褶皱单元分布如图10所示, 图中标有+号的单元为处于单向褶皱状态的单元。

山西大同的一个膜结构为两个马鞍形相对布置。膜材物理参数为1E=900MPa, 2E=648MPa, G12=31.4MPa, v12=0.018, v21=0.025, 厚为0.001m, 钢索的物理参数为E=2.0×105 MPa, 截面面积A=2.1×10-4m2。左侧膜材单元的预应力为σx=σy=2.0kN/m, 右侧膜材单元的预应力为σx=σy=2.5kN/m, 边索1的预拉力为40kN, 边索2的预拉力为50kN, 边索3的预拉力为75kN。平面尺寸如图11所示, 找形结果图如图12所示, 展开面积为251m2。下面给出其在风荷载作用下的受力分析结果。

大同50年一遇的基本风压为0w=0.55kN/m2, 取其高度变化系数μz=1.0, 风振系数βz=2.0, 采用计算流体动力学方法得到的风荷载体形系数μs如图13所示。由图可以看出, 迎风面的体形系数多为正值, 即受风压力作用, 最大值为0.8;背风面的体形系数多为负值, 即受风吸力的作用, 最大值为-1.8。按此施加风荷载计算结果如图14所示。图中带为+号的为单向褶皱单元, 黑色填充的为双向褶皱单元。

由褶皱图和体形系数图对照可以看出, 褶皱单元的分布与体形系数较大值的分布基本一致, 且由于风吸力比风压力要大, 虽然右侧膜材和钢索单元的预拉力均比左侧单元的大, 但在风压力作用下左侧结构只有单向褶皱的单元, 在风吸力作用下右侧结构还出现了一些双向褶皱的单元。

由上面的算例可以看出, 膜材出现双向褶皱的情况很少, 但是出现单向褶皱的情况还是比较多的, 这进一步说明在薄膜结构分析中考虑膜材褶皱的必要性。

在我国, 有些已建成的膜结构工程出现了撕裂和倒塌现象, 这一方面可能是和风荷载的取值有关, 但对结构正确的分析是保证结构安全的必要条件。膜材褶皱问题是薄膜结构受力分析中无法回避的问题, 要对薄膜结构进行正确的分析计算, 必须要解决这个问题。本文作者在以前工作的基础上, 依据主应力-主应变准则及修正本构矩阵的方法编制了有限元程序, 并通过经典算例证明了程序的正确性。