向量式有限元膜单元及其在膜结构褶皱分析中的应用-项目案例-污水池加盖-反吊膜|膜加盖-除臭加盖-膜结构公司-上海华喜膜结构工程有限公司
网站首页 解决方案 项目案例 新闻动态 膜材介绍 关于华喜 联系方式 EN
首页 > 新闻动态 > 行业动态

向量式有限元膜单元及其在膜结构褶皱分析中的应用

发布时间:2021年9月26日 点击数:1371

0引言

膜结构是一种被广泛使用的柔性结构,只具有很小的抗压或抗弯刚度,面外荷载作用下容易产生压缩应力从而发生褶皱。褶皱不仅会影响结构整体的美观,更会破坏结构的受力性能,甚至导致膜材的撕裂。

膜结构褶皱研究的主要方法之一是Wagner在分析工字形及箱形梁薄腹板最大剪力时提出的张力场理论[1],Miller等[2]将可变泊松比用应变关系代替,得到了褶皱状态时的膜材本构方程。Adler[3]将此褶皱模型嵌入有限元软件ABAQUS,发现薄膜中的压缩应力会使薄膜单元产生负几何刚度项,数值分析时出现刚度矩阵奇异问题。Roddeman等[4]通过修改应变张量使薄膜处于单轴拉伸状态,忽略薄膜的弯曲刚度,考虑褶皱对结构应力状态及力传递过程的影响,该方法在采用Newton-Raphson方法求解联立方程时在褶皱区会发生奇异,迭代很难收敛。Lu等[5]在Roddeman模型中使用曲线坐标得到了褶皱扩展方向。Nakashino等[6]也在此基础上提出了通过修改弹性张量实现对褶皱薄膜的分析。Fujikake等[7]针对张拉膜结构提出修正本构矩阵的方法,并在ADINA程序中实现; 谭锋[8]对该方法进行了改进,考虑初始预应力产生的应变,对膜结构褶皱进行了分析。

膜结构褶皱的另一种研究方法是基于薄壳理论的非线性屈曲分析,其基本思想是将膜结构的褶皱现象看成膜面的局部失稳,可有效模拟褶皱的面外形变,获得褶皱幅度、波长等具体信息,该方法主要应用于空间飞行器、卫星天线等高精度空间结构[9],国内也有学者将其应用到张拉膜材和气枕式充气膜结构的褶皱分析中[10]。由于膜结构在面外荷载作用下产生的大变形及褶皱效应,上述两种传统有限元方法在分析过程中均存在单元刚度矩阵奇异、迭代不易收敛等问题。

向量式有限元法[11,12,13]是美国普渡大学丁承先教授提出的一种基于点值描述和向量力学理论的分析方法。它以质点的运动来描述体系变形,不因其是结构或机构而改变; 其计算流程是逐点逐步循环,计算过程中无单元刚度矩阵,不存在刚度矩阵奇异的问题,无需求解复杂的非线性联立方程组,即也不存在迭代不收敛问题; 且通过引入逆向运动,可消除刚体位移所带来的数值误差。因此非常适合于分析大变位结构和机构问题。已有学者在向量式有限元膜单元开发方面做了一些工作,文献[14-15]分别推导了三节点三角形膜单元和四节点四边形膜单元的基本方程,并进行了算例分析。

作者首先推导三角形膜单元的向量式有限元基本公式,描述运动解析的原理及变形坐标系下单元节点内力的求解方法,编制计算程序并进行算例验证。然后将向量式有限元引入膜结构的褶皱分析中,采用主应力-主应变判别准则和修改本构矩阵的方法来处理膜结构的褶皱问题。最后对一个典型膜结构在面外荷载作用下的力学行为及褶皱发展过程进行跟踪分析。

1 向量式有限元膜单元的基本理论

向量式有限元基本单元由结构离散成的有质量的质点和质点间无质量的单元( 本文为常应变三角形单元) 组成。质点a的运动满足平动微分方程:

 


式中,ma为质点a的质量,¨xa为质点的加速度向量,α为阻尼参数,Fa为质点受到的合力向量,且a=fext_a+ fint_a,其中fext_a为质点的外力向量,fint_a为膜单元传递给质点的内力向量。

每个时间步内,通过中央差分公式数值求解质点运动方程式( 1) 获得时间步末时刻的质点新位置,进而得到该步的质点全位移。将式( 1) 转化为中央差分形式,包括无初始条件( 连续) 和有初始条件( 不连续) 时两种情况,如式( 2) 所示( 省略下标a) 。

 


式中,Δt为时间步长,

1. 1 单元节点位移

向量式有限元法采用逆向运动的方法,从单元节点全位移中扣除刚体位移( 包括刚体平移和刚体转动) 以得到纯变形量。图1所示为三角形膜单元abc在t0- t时段的空间运动图,在t0时刻单元处于abc位置,并于t时刻运动到a'b'c'位置。

图 1 膜单元的空间运动 Fig. 1 Space movement of membrane element

图 1 膜单元的空间运动 Fig. 1 Space movement of membrane element   下载原图


在运动过程中,单元经历的刚体位移包含刚体平移和刚体转动。选择质点a为参考点,则a点的位移向量ua为单元的刚体平移,让单元按向量( - ua)作逆向平移运动,得到各质点扣除刚体平移后的相对位移为:

 


扣除刚体平移后,计算单元的刚体转动,总的刚体转动向量θ可分为两个部分,即

 


式中,θ1表示单元平面外转动向量,θ2表示单元平面内转动向量。

如图2所示,单元扣除刚体平移后,处于a″b″c″位置,此时与t0时刻的单元abc并不处于同一平面内,平面外转动角度θ1( 平面外转轴单位向量nop,0从n转动到n″的角度) 为

 


平面外转轴单位向量

 


其中lop、mop和nop分别为沿x、y和z向的单位方向余弦; 从而有单元平面外转动向量θ1= θ1nop,0

图 2 膜单元的平面外转动 Fig. 2 Out-of-plane rotation of membrane element

图 2 膜单元的平面外转动 Fig. 2 Out-of-plane rotation of membrane element   下载原图


质点i( i = a,b,c) 的平面外逆向刚体转动位移向量可由下式计算得到:

 


式中,平面外逆向转动矩阵R*op( - θ1) =[1 - cos( - θ1) ]×A2op+ sin( - θ1) Aop,单元的边向量r″ai= x″i- x″a= x'i- x'a,

 


在经历平面外的逆向转动之后,单元回到t0时刻所处的平面上,得到各节点坐标为:

 


平面内的逆向转动位移示意见图3。

图 3 膜单元的平面内逆向转动位移示意 Fig. 3 In-plane rotation of membrane element

图 3 膜单元的平面内逆向转动位移示意 Fig. 3 In-plane rotation of membrane element   下载原图


平面内刚体转动角度θ2

 


式中,Δφi( i = a,b,c) 为质点面内转动的角度,根据质点与形心连线向量的转动夹角计算得到: |Δφi| =cos- 1( ei·ei) ,( i = a,b,c) ,|Δφi|∈[0,π]。其中, ,xi和xC分别为单元节点i和形心C的坐标向量。

Δφi的正负可根据下式判断:

 


平面内转轴单位向量

 


即为膜单元在t0时刻的单位法向向量,其中lip,mip,nip分别为沿x,y,z方向的单位方向余弦。

从而有质点i ( i = a,b,c) 的面内转动向量为Δφi= Δφin0,Δφi∈[- π,π],单元平面内转动向量

参考式( 7) 可以得到单元质点在平面内的逆向转动位移为:

 


式中,平面内逆向转动矩阵 ,单元的边向量 ,其中

扣除单元的刚体平移、面外刚体转动和面内刚体转动后,得到单元质点的纯位移向量为:

 


单元节点内力

单元变形满足的虚功方程为

 


其中fi为整体坐标系下膜单元节点i的内力。

求解膜单元节点内力时,向量式有限元采用定义变形坐标系的方法将空间膜元问题转化到平面上,如图4所示定义,变形坐标系的原点在参考节点a,这样节点a的位移为零,即可消除两个未知变量。定义坐标轴 ^x的方向为节点b的变形位移方向,则可以消除节点b的竖向位移分量。

图 4 单元变形坐标系 Fig. 4 Deformation coordinate system of element

图 4 单元变形坐标系 Fig. 4 Deformation coordinate system of element   下载原图


变形坐标系与整体坐标系之间的转换关系为

 


式中 是坐标系间转换矩阵,

变形坐标系下必有 ^xz= 0,因而式( 15) 可简写成

 


在变形坐标系中单元节点的位移向量为

 


由于变形坐标系下 简写成

 


引入传统有限元的形函数,得到单元上任意节点的位移

 


式中, 是单元形函数,形式同传统有限元。

单元节点位移的组合向量可记为

 


由于变形坐标系下必有 ,因而式( 20) 可简写成

 


由单元形函数描述的整个膜单元变形分布,可获得单元的应变分布为

 


式中, ,B的形式同传统有限元,B*为去除B中对应 的行和列后获得的简化形式。

假设材料为线弹性,在得到单元的应变矩阵后,可得到单元的应力矩阵

 


则单元的变形虚功为

 


式中, 为单元的初始应力。

因而有节点内力向量

 


其中, 另外三个分量 可由平面膜单元的静力平衡条件得到。

式( 25) 求得的是单元在变形坐标系下虚拟状态的节点内力分量,需将其转化为t时刻的内力。首先将变形坐标系下的单元节点内力向量扩展为三维向量形式 再通过转换矩阵得到整体坐标系下的内力,然后经过正向运动转换( 包括刚体转动和刚体平移,后者不影响转换) 得到t时刻的内力为

以上所求得的fi是单元节点i由于单元纯变形产生的内力,反向作用于对应单元节点i的质点a上即得到三角形膜单元传给质点a的力fint_a

2 若干问题处理

2. 1 临界时间步长和临界阻尼系数的估值

向量式有限元以质点为研究对象,质点与相邻质点间的相互作用是由多个应力和应变量来描述的,可等效简化成质点为m,弹簧刚度为k的单自由度体系,其自振圆频率为 ,临界阻尼为Ccr=2mω。采用显式中央差分方法求解质点运动方程式时,临界时间步长 临界阻尼系数 选取的实际时间步长需要满足h < hc。

固体内部点与点之间的相互作用宏观上可分为轴向拉伸、剪切和弯曲作用,一般轴向拉伸刚度最大,弯曲刚度最小,对于薄膜单元仅存在轴向拉伸刚度。因而采用轴向拉伸刚度来近似估算临界时间步长hc的下限值。轴力拉伸组件的kt= EA / l,则有

对于临界阻尼系数α 本文假定存在薄膜弯曲刚度D < E/200( 此时已可视为忽略弯曲作用) 来近似估算其下限值。选取的实际阻尼系数需要满足α < αc,且应避免α过大( 过阻尼) 而无振荡效应, α过小而振荡效应显著,因此导致趋于平衡状态非常缓慢。对于薄膜单元,有 对于静力问题,通过引入阻尼系数α来实现最终趋于静力平衡状态,实际选取阻尼系数α值时可取小于αc的接近值; 对于动力问题,无阻尼时直接取α = 0即可。

2. 2 节点应力的插值计算

通过本文上述方法可获得三角形CST薄膜单元的应力( 常应力单元) ,实际需要的通常是单元节点处的应力,因而需要对计算得到的应力进行插值处理。本文中采用较为简单的应力处理方法,即取围绕节点各单元形心处von Mises应力的平均值作为该节点的von Mises应力。这样得到的节点von Mises应力值是围绕该节点的有限区域内的von Mises应力平均值。

3 膜材褶皱的处理

3. 1 褶皱判别准则

目前对膜材褶皱的判别方式有三种[8]: 主应力准则、主应变准则和主应力-主应变准则。主应力判断准则在得到主应力之前无法判断褶皱是否发生,会出现判断不准确的情况; 主应变准则对于存在拉剪耦合作用的正交异性膜材无法使用; 主应力-主应变准则可以准确得到膜材的受力情况,通过主应力判断是否受拉,主应变判断是否松弛,对各向异性材料也适用。

采用主应力-主应变准则来判别以下3种受力状态。设膜材主应力σ1≥σ2,主应变ε1≥ε2,则:

1) 当σ2> 0时,膜材处于纯拉状态,无褶皱产生;

2) 当ε1≤0时,膜材处于松弛状态,此时发生双向褶皱;

3) 除上述两条以外的情况,膜材发生单向褶皱。

3. 2 褶皱处理方法

采用张力场理论中修正本构矩阵的方法[8]处理膜的褶皱问题。对于存在拉剪耦合作用的膜材,其主应力方向上的应力-应变关系为

 


式中,T是主应力方向本构矩阵与弹性主轴方向本构矩阵的转换关系矩阵,可表达为

 


其中β是主应力方向与弹性主轴方向夹角。

单向褶皱时,根据σ2= 0,可以在式 ( 28 ) 中消去ε'1、ε'2,得到

 


 


双向褶皱时,膜材失去刚度,计算中只需将膜材的刚度矩阵置零即可,即令式( 30) 中a0= 0。

3. 3 应力转换计算

向量式有限元中的纯变形是在虚拟位置时的变形坐标系下得到的,而在初始时刻t0的初始应力σ0是整体坐标系下的应力分量,因此在计算过程中需要转换到虚拟状态的变形坐标系分量才能与得到的应力增量进行处理,即

 


经式( 31) 处理后的应力可以与计算得到的应力增量在虚拟状态相加减,得到的总应力再转换回到整体坐标系,并经过正向运动转换( 包括刚体转动和刚体平移,后者不影响转换) 回到t时刻有:

 


 


式中,Rp= RipRop

根据得到的t时刻单元的应力分量,可由材料力学方法计算平面单元的主应力,同理可以求出单元主应变,根据前文述及的主应力-主应变准则可以判断单元是否发生褶皱或松弛,然后由判别得到的情况,根据式( 32) 、( 33) 对单元的弹性矩阵进行修正,并采用修正后的弹性矩阵计算单元内力,得到单元内力之后采用向量式有限元方法计算单元作用在质点上的内力,代入质点运动控制方程即可进行膜结构在荷载作用下的褶皱分析。

4 算例分析

基于上文所推导的三角形膜单元基本公式,首先采用Matlab编制了膜单元的向量式有限元分析程序。下面通过算例验证所推导公式和编制程序的正确性。

4. 1 算例 1

为了说明编制的膜单元程序在大角度刚体转动情况下的计算稳定性,以图5所示方形膜片为例进行分析。膜材质量密度ρ = 1. 0 kg/m3,弹性模量E =104Pa,泊松比υ = 0,膜厚度为0. 01 m,几何尺寸为1m×1 m。采用5个质点4个三角形膜单元来模拟其面内及面外刚体转动,分别通过给5个质点施加一定初速度( 图5) 来激励转动。

图6a、6b分别为面内旋转和面外旋转的图形变化,图6c、6d则为对应的质点1在旋转过程中坐标随时间的变化情况。计算结果符合膜片的实际运动情况,且与文献[14-15]结果一致,说明本文编写的计算程序是正确的,可用于处理膜单元的刚体位移问题。

图 5 平面膜片刚体转动算例 Fig. 5 Numerical example of rigid-body rotation for planar patch

图 5 平面膜片刚体转动算例 Fig. 5 Numerical example of rigid-body rotation for planar patch   下载原图


图 6 平面膜片刚体转动计算结果 Fig. 6 Results of rigid-body rotation for planar patch

图 6 平面膜片刚体转动计算结果 Fig. 6 Results of rigid-body rotation for planar patch   下载原图


4. 2 算例 2

为测试膜单元的收敛性能,对图7所示受内压的薄膜结构进行分析,其初始形状由以下双抛物线曲面函数所确定:

 


取矩形膜片边长m = n = 1 m,矢高f = 0. 1 m,四边固定,膜材质量密度ρ = 106kg / m3,弹性模量为E = 6. 0 GPa,泊松比υ = 0. 267,膜厚度1. 0 mm,薄膜结构划分为200个三角形单元,在薄膜内部施加1. 5 MPa内压。在内压作用下薄膜结构将向上膨胀,为尽快得到收敛解,采取缓慢加载和施加阻尼的方式减除动力震荡效应; 根据本文2. 1节所推导的公式,分别取时间步长h =4. 0×10- 4,阻尼系数α =100。

图8a所示为内压作用下膜的变形图,图8b为内压为1. 0 MPa时膜顶点位移计算的收敛过程,图8c为本文计算结果与文献[16-17]结果的比较。可以看出,本算例有很好的收敛结果,且与文献[16-17]结果一致。

图 7 双抛物线膜结构 Fig. 7 Double-parabolic-surface membrane structure

图 7 双抛物线膜结构 Fig. 7 Double-parabolic-surface membrane structure   下载原图


图 8 内压作用下双抛物线膜结构中心点位移计算结果 Fig. 8 Results of center point displacement of doubleparabolic-surface membrane structure under internal pressure

图 8 内压作用下双抛物线膜结构中心点位移计算结果 Fig. 8 Results of center point displacement of doubleparabolic-surface membrane structure under internal pressure   下载原图


5 膜结构褶皱分析程序编制和算例

在上节膜结构基本程序的基础上,利用本文第3节中给出的方法,即采用主应力-主应变准则和修改本构矩阵方法,进一步编制了可考虑褶皱效应的膜结构分析程序,其流程如图9所示。

为验证程序的可靠性,对图10所示一典型的菱形膜结构算例进行分析,相关参数取值参考文献[8]。菱形膜结构平面投影对角线长6 m,四边固定,高点A、B与低点C、D之间的高差为1 m。AB、CD方向分别为膜材的经向和纬向,经向和纬向的弹性模量分别为E1= 900 MPa,E2= 600 MPa,剪切模量G = 20 MPa,泊松比分别为υ12= 0. 6,υ21= 0. 4,厚度为1 mm,预应力σx= σy= 2 k N / m,作用均布向上荷载为qz= 2. 5 k N / m2。结构划分为200个三角形膜单元,如图10所示。

图 9 流程图 Fig. 9 Flow chart

图 9 流程图 Fig. 9 Flow chart   下载原图


图 10 菱形膜结构 Fig. 10 Diamond membrane structure

图 10 菱形膜结构 Fig. 10 Diamond membrane structure   下载原图


采取逐级加载的方式,跟踪膜结构在外力作用下的变形全过程,并捕捉发生褶皱的区域。图11给出了若干典型荷载步下的膜结构变形及褶皱单元分布情况。当均布荷载很小时,结构变形很小,不会出现褶皱; 当荷载为0. 3 k N/m2时,结构开始出现单向褶皱,此时变形仍较小,发生褶皱的单元也较少; 当荷载达到0. 4 k N/m2时,结构的变形增大,出现的褶皱单元也变多,主要沿高点A、B之间分布; 当荷载达到0. 5 k N/m2时,结构中出现的单向褶皱单元达到最多,数量约占总单元数的46% ,此时结构变形较大,膜面出现不光滑; 当荷载达到1. 0 k N/m2时,结构中单向褶皱的单元数目有所减少,但是出现了若干双向褶皱的单元,表明部分膜面出现松弛,此时膜面变得更不光滑; 当荷载达到1. 5 k N/m2时,结构单向褶皱数目继续减少,双向褶皱的单元消失,结构变形继续变大; 当荷载达到2. 5 k N/m2时,结构中褶皱的单元基本消失,结构在向上荷载作用下,继续发生变形,最终可以达到正高斯曲面的平衡状态。该算例表明,利用本文方法及程序,可以成功捕捉膜结构从褶皱出现、褶皱大范围产生直至褶皱最终消失的全过程,克服了一般有限元方法中当膜面由于褶皱发生局部畸变而形成机构时的计算不收敛问题。

图 11 褶皱单元分布与结构变形 Fig. 11 Wrinkled element distribution and structural deformation

图 11 褶皱单元分布与结构变形 Fig. 11 Wrinkled element distribution and structural deformation   下载原图


在本算例对称膜结构的褶皱分析全过程中,当荷载很小( 图11a) 或荷载很大( 图11f) 时,均仅在边界存在小范围的单向褶皱,此时褶皱位置满足严格的对称性; 当荷载较小( 图11b) 或荷载较大( 图11e)时,单向褶皱范围有所扩大,由于褶皱问题本身的复杂性和计算误差的存在,除个别位置外大部分褶皱位置仍是满足对称的; 当对应褶皱范围很大的荷载( 图11c) 或甚至出现松弛的荷载( 图11d) 时,在基本满足褶皱对称分布的情况下,褶皱和松弛引起的褶皱不对称分布的影响也趋于增大。

上述膜结构褶皱从出现、扩展到最终消失的过程和褶皱位置的( 非) 对称性分布情况与文献[8]的计算结果基本一致,验证了本文方法在膜结构褶皱分析问题中的有效性和准确性。图12给出了荷载-中心节点z向位移dz曲线,并与文献[8]的结果进行了比较,可见,两者的计算结果符合较好。

图 12 荷载-中心节点位移曲线 Fig. 12 Load-displacement curve of central node

图 12 荷载-中心节点位移曲线 Fig. 12 Load-displacement curve of central node   下载原图


图13所示为不同荷载条件下中心节点O的竖向位移收敛情况。图13a为荷载0. 5 k N/m2条件下节点的位移-时间曲线,可以看出,由于褶皱的出现,计算过程中出现在收敛值附近范围内波动而不收敛的现象; 图13b为荷载2. 5 k N/m2条件下节点的位移-时间曲线,由于最终结构重新进入双向张拉状态,最终计算达到平稳的收敛解。

图 13 不同荷载条件下的中心节点位移-时间曲线 Fig. 13 Center nodal displacement-time curves under different loading conditions

图 13 不同荷载条件下的中心节点位移-时间曲线 Fig. 13 Center nodal displacement-time curves under different loading conditions  下载原图


6 结论

1) 文中推导了向量式有限元三角形膜单元的基本方程,详细阐述了通过逆向运动处理膜单元的平面内、外刚体位移从而获得单元节点纯变形的过程,以及进一步通过纯变形得到单元节点内力的求解方法。在此基础上编制了膜单元的计算程序,算例分析表明,所编制的程序可以很好地完成膜结构的刚体运动及弹性变形计算,验证了推导结果的正确性。

2) 进一步将向量式有限元引入膜结构的褶皱分析,采用主应力-主应变准则和修正本构矩阵方法处理膜结构的褶皱问题,编制程序对典型膜结构在荷载作用下的褶皱发展情况进行跟踪分析,成功捕捉了从褶皱出现、褶皱大范围扩展直至褶皱最终消失的全过程,可有效克服传统有限元方法存在的单元刚度矩阵奇异、迭代不易收敛等问题。

专题报道             more...
  • 轨道交通中膜结构的应
    ...

    查看更多

  • 膜结构建筑保温内衬技
    刚查县为青海省海北藏族自治州辖县,青海省措温波高原海滨藏城演艺中心,作为刚查县的标志性建筑,演艺中心为直径50米的圆形建...

    查看更多

  • 膜结构幕墙的应用
    膜结构幕墙是膜结构在建筑外围护结构的应用,具有膜结构的共同特性和优点:膜结构是一种非传统的全新结构方式。...

    查看更多

  • 膜结构屋面的应用
    屋盖是房屋最上部的围护结构,应满足相应的使用功能的要求,为建筑提供适宜的内部空间环境。屋盖也是房屋顶部的承重结构,受到材...

    查看更多

  • 膜结构应用于环保工程
    随着我国国民经济飞速发展和市政基础设施建设全面展开,特别是污水处理厂等环保项目日益增多,其中有相当数量的污水处理厂的厌氧...

    查看更多

  • 膜结构在污水处理厂中
    相当数量的污水处理厂的厌氧池、污泥浓缩池、生物絮凝池等建于居民区、厂区的周边,污水池的环境、风貌及污水臭味等直接影响人们...

    查看更多

关于华喜

硬件实力 质量控制 发展历程 公司简介

软件实力 经营理念  解决方案 联系方式

中国华喜建筑网站

+021-59198545 400-176-6885 dshx@hxmjg99.com www.hxmjg.com 沪ICP备08009856号 使用条款