薄膜结构的膜材厚度一般是微米级或者毫米级别的, 压缩和弯曲刚度很小, 当受到压应力作用时, 薄膜会产生褶皱, 单向或者双向褶皱的扩展都会使局部稳定性丧失, 因此, 研究膜褶皱的力学特性, 对设计膜结构具有重要的意义.笔者编制了判断褶皱产生的有限元程序, 并且给出算例, 验证了程序的可靠性.

目前, 研究褶皱主要有2种理论分析方法, 即张力场理论和屈曲分析理论.张力场理论从膜中存在的应力场分析面内褶皱的区域、方向, 但不能得到具体褶皱面外变形参数;屈曲分析是考虑褶皱与薄膜的弯曲刚度, 结合褶皱变形与屈曲相似, 认为褶皱是膜的局部失稳, 可以预测分析褶皱的具体特征.

Stein and Hedgepeth[1]1961年给出褶皱膜应力场表达式, 给出3种典型模型算例的褶皱区域, Reissner, Mansfield, Epstein Forcinito等基于张力场理论做了不少研究工作, 以张力场为基础的研究方法主要有3种, 可变泊松比方法, 修正的变形梯度方法, 修正的本构矩阵方法.

随着计算机的发展, 研究人员开始结合分叉屈曲理论, 采用有限元数值模拟方法研究褶皱现象, Wong和Pellegrino[2,3,4,5]从2000年开始, 利用ABAQUS, 采用屈曲有限元方法, 考虑膜的微小弯曲刚度, 用薄壳shell单元, 结合后屈曲理论, 分别对矩形膜的剪切模型和正方形膜四角受拉褶皱形状进行分析, 得到褶皱变形的具体褶皱参数 (褶皱的波长和幅值) 的具体表达式.

Tomoshi Miyamura[6]在2000年运用分叉屈曲理论结合有限元, 研究了圆形张拉膜在面内扭矩作用下产生褶皱的模拟.考虑几何非线性, 采用四结点等参元, 分析在一定网格密度下的有限数量的褶皱.

本文采用张力场理论, 编制程序, 对膜的褶皱区域进行判断.

前面提到的Stein和Hedgepeth 理论, 又称褶皱理论, 描述了部分褶皱的膜单元的受力状况, 膜主要有纯拉、单向褶皱和双向褶皱3种受力状态, 因此膜出现褶皱时, 有3种褶皱判断准则, 即主应力准则、主应变准则和主应力-主应变准则.

1) 主应力准则σ2>0, 张紧状态;σ2≤0, σ1>0, 褶皱状态;σ2≤0, σ1≤0, 松弛状态;

2) 主应变准则ε2>0, 张紧状态;ε2≤0, ε1>0, 褶皱状态;ε2≤0, ε1≤0, 松弛状态;

3) 主应力-主应变准则σ2>0, 张紧状态;σ2≤0, ε1>0褶皱状态;σ1≤0, ε1≤0, 松弛状态.

结合膜几何非线性有限元基本理论[7], 在局部坐标系中按平面应力问题分析.利用张力场理论方法, 对无褶皱的膜结构程序[8], 利用主应力-主应变准则, 针对出现单向褶皱状态的膜单元, 修改本构矩阵, 编制程序[9], 判断褶皱单元的出现及褶皱区域.

其中, 在主应力和主应变方向上, 本构关系DˉˉˉDˉ为膜单元在应力主轴坐标系内的本构矩阵.

⎧⎩⎨⎪⎪σ1σ20⎫⎭⎬⎪⎪=⎡⎣⎢Dˉˉˉ11Dˉˉˉ21Dˉˉˉ31Dˉˉˉ12Dˉˉˉ22Dˉˉˉ32Dˉˉˉ13Dˉˉˉ23Dˉˉˉ33⎤⎦⎥⎧⎩⎨⎪⎪ε1ε2γ12⎫⎭⎬⎪⎪{σ1σ20}=[Dˉ11Dˉ12Dˉ13Dˉ21Dˉ22Dˉ23Dˉ31Dˉ32Dˉ33]{ε1ε2γ12}

, (1)

展开后写成 (2) ～ (4) 式

σ1=Dˉˉˉ11ε1+Dˉˉˉ12ε2+Dˉˉˉ13γ12σ1=Dˉ11ε1+Dˉ12ε2+Dˉ13γ12, (2)

σ2=Dˉˉˉ21ε1+Dˉˉˉ22ε2+Dˉˉˉ23γ12σ2=Dˉ21ε1+Dˉ22ε2+Dˉ23γ12, (3)

0=Dˉˉˉ31ε1+Dˉˉˉ32ε2+Dˉˉˉ33γ120=Dˉ31ε1+Dˉ32ε2+Dˉ33γ12, (4)

整理得σ2=[Dˉˉˉ21−Dˉˉˉ23Dˉˉˉ31Dˉˉˉ33]ε1+[Dˉˉˉ22−Dˉˉˉ23Dˉˉˉ32Dˉˉˉ33]ε2σ2=[Dˉ21-Dˉ23Dˉ31Dˉ33]ε1+[Dˉ22-Dˉ23Dˉ32Dˉ33]ε2. (5)

再结合σ2=0得主应变和主应力之间关系式

ε2=Dˉˉˉ21Dˉˉˉ33−Dˉˉˉ23Dˉˉˉ31Dˉˉˉ23Dˉˉˉ32−Dˉˉˉ22Dˉˉˉ33ε1ε2=Dˉ21Dˉ33-Dˉ23Dˉ31Dˉ23Dˉ32-Dˉ22Dˉ33ε1, (6)

σ1=1Dˉˉˉ23Dˉˉˉ32−Dˉˉˉ22Dˉˉˉ33[Dˉˉˉ11(Dˉˉˉ23Dˉˉˉ32−Dˉˉˉ22Dˉˉˉ33)+Dˉˉˉ12(Dˉˉˉ21Dˉˉˉ33−Dˉˉˉ23Dˉˉˉ31)+Dˉˉˉ13(Dˉˉˉ31Dˉˉˉ22−Dˉˉˉ21Dˉˉˉ32)]ε1σ1=1Dˉ23Dˉ32-Dˉ22Dˉ33[Dˉ11(Dˉ23Dˉ32-Dˉ22Dˉ33)+Dˉ12(Dˉ21Dˉ33-Dˉ23Dˉ31)+Dˉ13(Dˉ31Dˉ22-Dˉ21Dˉ32)]ε1. (7)

令a=1Dˉˉˉ23Dˉˉˉ32−Dˉˉˉ22Dˉˉˉ33[Dˉˉˉ11(Dˉˉˉ23Dˉˉˉ32−Dˉˉˉ22Dˉˉˉ33)+Dˉˉˉ12(Dˉˉˉ21Dˉˉˉ33−Dˉˉˉ23Dˉˉˉ31)+Dˉˉˉ13(Dˉˉˉ31Dˉˉˉ22−Dˉˉˉ21Dˉˉˉ32)]a=1Dˉ23Dˉ32-Dˉ22Dˉ33[Dˉ11(Dˉ23Dˉ32-Dˉ22Dˉ33)+Dˉ12(Dˉ21Dˉ33-Dˉ23Dˉ31)+Dˉ13(Dˉ31Dˉ22-Dˉ21Dˉ32)],

则写作σ1=aε1, 因此褶皱区域得本构矩阵可以表示成 (8)

Dˉˉˉσ=⎡⎣⎢a000ηa000ηa⎤⎦⎥Dˉσ=[a000ηa000ηa]

.其中, η为避免奇异矩阵产生的系数.

给出判断褶皱单元的流程图, 如图1所示.

建立矩形膜有限元模型, 由于膜是对称的, 有限元分析时取右半侧, 按照线弹性各向同性材料分析.参数分别取为弹性模量E=900 MPa, 泊松比v=0.3, 长度l=5.0 m, 高度h=3.0 m, 厚度0.1×10-2m.在ANSYS中建立模型, 采用shell63单元, 如图2所示, 矩形左边中点的x, y方向加约束, 其余边只约束x方向.外荷载包括拉力、弯矩都用等效节点荷载表示.保持预拉力不变, 即在整个膜片上初拉应力0.02 MPa不变, 再加载.

增大弯矩M的值, 发现负应力区域逐渐增大, 褶皱区域随着弯矩增加的变化趋势, 图3-5所示为几个有代表性的弯矩作用下的变形图及主应力分布图.由于shell63单元具有一定的抗压刚度, 板的中性面不会上移, 即负应力区域不会大于0.5倍的高度, 即负应力总是在中性层以下的区域, 如图3-5所示为几个特殊弯矩作用下x方向的变形应力图.

针对图6的理论模型[1], 拉力P和弯矩M以及均匀拉力Ny, 产生均布应力σ.文献针对具有褶皱膜, 给出褶皱区域的解析解, 其中边界线y=b代表褶皱区域, 拉力P和弯矩M与x方向的正应力σx关系为

P=t∫h0σxdyM=t∫h0σx(y−h2)dy⎫⎭⎬⎪⎪⎪⎪⎪⎪⎪⎪Ρ=t∫0hσxdyΜ=t∫0hσx(y-h2)dy}

. (9)

并且褶皱区域存在下列关系式:

bh=1−2PkEth2−−−−√2Mh=1−232PkEth2−−−−√⎫⎭⎬⎪⎪bh=1-2ΡkEth22Μh=1-232ΡkEth2}

. (10)

计算结果对比如图7, 本文程序计算结果与理论解比较接近, 说明程序计算膜的褶皱区域变化趋势的正确性, 验证了理论解[1] , 从而验证程序在判断褶皱区域的可靠性.

由于ANSYS软件计算褶皱单元, 是按照薄壳单元计算, 薄壳单元具有一定的刚度, 褶皱区域不会大于板的一半高度, 具有一定的局限性.由于膜相对于薄壳来说, 几乎没有抗弯刚度, 在受压力时中性面发生了上移, 褶皱区域增大到中性层上部.

本文中采用的修改刚度矩阵方法, 即结合主应力-主应变准则, 对出现褶皱的程序单元修改本构关系矩阵, 通过理论和程序分析, 对于研究膜的褶皱区域具有重要的理论指导意义.

虽然本程序在一定程度上能够实现膜的褶皱区域的判断, 但是仍旧有一些数值解收敛问题, 比如η值的选取问题等.

综合目前的数值模拟研究现状, 主要是针对某些具体的特殊形状的平面膜结构进行分析, 没有给出一般荷载和边界条件下的褶皱特征.褶皱对整个膜结构应力分布的影响, 目前也没有具体的计算和分析, 因此, 有必要研究局部褶皱对整体膜结构性能影响.用商业软件分析膜结构, 计算耗时, 出现褶皱时按照无压应力单元处理, 有时不能得到收敛解.因此, 能否建立一个褶皱宏单元模型, 即把褶皱区域看作一个或几个宏单元, 使得计算程序收敛较快, 而且适合不同边界、荷载情况下的褶皱情况, 来研究局部失稳对整体结构的影响.