建筑膜结构中最常用的织物膜材是一种在纤维基体上涂覆树脂涂层而制成的柔性复合材料,其抗压和抗弯能力很弱,主要依靠膜面曲率变化和薄膜力来抵抗外力,荷载作用下膜面构形容易发生较大改变,因此相应的体系平衡条件必须建立在变形后的构形上,并计入大位移引起的二阶或更高阶效应,按照几何非线性分析方法进行求解[1]。同时,织物膜材本身还具有正交异性和非线性拉伸特性,这使得对于膜材本构关系的表征和变形行为的描述也较一般结构更为复杂。很多膜结构分析设计限定于小应变、大位移的情况,将膜材假定为线弹性体,忽略材料非线性而只考虑几何非线性,显然这种简化计算结果较实际情况会有偏差,同时也限制了膜材强度的利用率。
早期针对薄膜大变形的分析主要是基于理论方法,也获得了一些具有简单边界和规则几何形状的薄膜大挠度问题的解析解[2],但一般都要求施加的荷载不能过大而且要符合对称性条件,这显然无法满足实际膜结构分析的需要。因此,目前的分析主要采用数值方法,比如有限单元法、动力松弛法、应力密度法、质点-弹簧模拟法、边界元法、无网格法等,其中以非线性有限元法发展最为成熟。比如,近年来由Onãte等[3]提出的有限单元片法较常应变(C0)和线性应变(C1)单元的精度和效率都有所改善,已在建筑与航天膜结构的分析中获得了应用;Bletizinger等[4]对该方法做了进一步完善,通过引入位形迭代算法使其在找形、荷载分析过程中能够考虑裁剪效应和优化的影响。然而,以连续介质力学理论为基础的上述方法用于处理薄膜的大变形问题仍较为繁琐,通常需借助复杂的非线性迭代算法,并对收敛性和稳定性做特别考虑,特别是需要同时涉及几何、材料等强非线性效应时,求解难度将大大增加。
为了研究机织建筑膜材的非线性拉伸行为,有很多学者根据纱线和涂层材料的几何形态特征和变形机制,从微观角度提出了一些由梁、杆、板等元素构成的膜材内部组织结构作为其力学分析模型[5],较成熟的有Kato栅格模型[6]和Pargana菱形截面模型[7]等,这种方法虽然可以较为精确地描述膜材的拉伸性能,但模拟过程复杂,需要借助相应的力学理论建立对应的函数模型才能用于结构整体分析。Reese等[8]为分析充气膜结构的材料大变形问题,基于由微观力学模型计算得到的大量“数值试验”数据,构建了可表征各向异性膜材非线性拉伸特性的广义应变能函数;Goncalve和Campello[9]曾提出一种比较实用的非线性正交异性弹性材料模型,其中将总应变能分成各向同性和正交异性两个部分,分别对应于膜材的弯曲变形和薄膜变形。不过,文献[8―9]均没有考虑到膜材形面变化较大时经纬线正交异性发生改变的情况。
有限质点法[10]是基于向量式力学的基本理论提出的面向工程结构复杂行为的数值分析方法。与传统分析方法相比,它不再以函数连续体来描述分析对象,而是将求解域离散为相互联系的质点系统,着眼于该系统内各个质点所受的力,结构静、动力问题均通过质点运动方程的时间积分求解,从而获得到结构整体的运动变形状态,其优势表现在分析结构强非线性行为时无需迭代求解和特殊修正。目前该方法已经应用于机构运动、接触碰撞、结构屈曲、倒塌破坏、固体大变形、膜结构找形、结构多尺度精细化等结构复杂行为的分析[11―17],均取得了较好的效果。
本文通过理论推导和数值分析,将有限质点法推广到薄膜结构的大变形行为分析中。文中首先给出薄膜的有限质点法基本理论和控制方程,阐述有限质点法处理几何、材料非线性因素的基本思路;然后根据有限质点的计算特点建立正交异性织物膜材的线性和非线性拉伸本构模型,推导本构矩阵的计算方法,并给出一种实用的材料主轴方向确定方式;最后通过算例的模拟和分析,检验该方法的正确性和有效性。
根据膜材的变形特点,有限质点法将薄膜结构离散为一组具有三维平动自由度的质点集合,而质点间由一系列仅有面内变形的膜元相连。结构上所有的质量、内力和外力作用均由质点承担,通过有限个质点在微小时间段上的运动轨迹来描述整个结构在空间和时间上的连续变形和运动状态,即称之为“点值描述”。
如图1所示,假设一块膜面在外力作用作用下经历了一个连续变形和运动过程,在初始时刻t0=0、中间时刻ta=t-τ以及当前时刻t时的几何构形分别为V0、Va和V,该时间历程中(t0-t)膜面上所有质点运动的时间轨迹被分割成一系列独立的时间片段(t0, t1, t2, …, t),每一个时间段的运动轨迹即是一个途径单元,问题转化为求解膜面质点集合在各个途径单元内的位置、速度等物理量的变化情况。
图1 薄膜结构的“点值描述”模型
Fig.1 Point description model of a membrane structure
膜面质点在每个途径单元内的运动都遵循牛顿第二定律,且运动轨迹独立。对于分布于空间膜面上的任意一个质点α,运动变量可分解为沿空间坐标轴方向的三个平动自由度分量,对应坐标轴方向上的三个质点力。以为参考,点α的全位移向量为
途径单元ta≤t≤tb内质点运动方程式可表达为:
式中:Mα为质点α的质量,包括该点的集中质量及相连薄膜元所提供的等效质量;为直接作用在质点α上的集中外力向量;
分别为与质点α相连接的薄膜元i提供的等效外力向量和内力向量;
为作用在质点α上的阻尼力[15];nc代表与质点相连的薄膜元个数。
有限质点法采用显式时间积分法求解运动方程,避免复杂迭代和收敛问题。将质点的速度和加速度写成位移和时间的中央差分表达式,代入式(1)后即可得到途径单元内各时间步的质点位移、速度等物理量。需要指出,由于显式直接积分算法的收敛性为有条件稳定,故计算时间步长需满足每一步差分运算的误差控制条件。根据稳定性准则[18],薄膜的计算临界步长应满足:
式中:ce为薄膜元内当前波速E、ν和ρ分别为膜材弹性模量、泊松比和密度
为薄膜元的最小特征长度,
其中β对三角形取0,四边形取1,Ae为薄膜元面积,Li (i=1,2,3,4)为边长。
由于有限质点法中控制方程是一种动力形式的平衡方程,要获得所求问题的静态解可通过引入虚拟阻尼耗能和斜坡-平台式缓慢逐级加载两种方式将各质点振幅控制在一个有限区间[10],并持续削弱其振荡效应直至稳定在允许残差范围内。为进一步减少平衡位置附近振荡,加快静力解收敛速度,本文计算中对于与系统动能耗散效率密切相关的阻尼参数采用虚设值,将参数μ取为由系统最小固有频率计算得到的临界值[19],即:
式中:Ccr为临界阻尼系数;ωmin为系统自由振动基频,可由Rayleigh法[20]计算得到。
无论是虚拟阻尼耗能还是缓慢加载方式,实际上都属于对静力平衡态的逼近算法而不可能获得绝对的静止状态,因此本文在静力问题求解中综合选用下列收敛准则作为计算结束的判断条件:1) 残余力收敛准则,即位移收敛准则,即
εf和εd分别为设定的残余力和位移收敛控制精度,取值一般为10-3~10-5。
膜面质点内力仅与相连膜元的纯变形相关,有限质点法通过虚拟逆向运动扣除刚体平移和刚体转动,计算膜元的纯变形和内力。以图2所示膜元为例,在ta和t时刻的位置分别为Va(1a-2a-3a)和V(1-2-3),以Va为参考构形,令膜元做逆向虚拟平移(-Δu1)和转动(-θ)至虚拟位置Vr(1r-2r-3r),求得膜元的纯变形再引入变形坐标
得到膜元面内的非零独立节点变形位移
求出膜元的纯变形后,可借助修正的内插函数由虚功方程求出虚拟位置上的节点等效内力向量再通过正向运动与坐标变换得到对应当前t时刻膜元构形的真实内力,将内力反向集成到相连质点上,最终得到全局坐标下的质点合内力。
图2 膜元纯变形的求解过程
Fig.2 Solution process of the membrane pure deformation
膜面质点外力与内力的计算类似,考虑到分布荷载与节点力之间的等效关系不会因变形和刚体位移而发生改变,可先将分布荷载q(x)=经逆向运动后转移到虚拟位置Vr来求等效关系,然后随膜元经正向运动恢复到t时刻的真实位置V,即得到当前构形下与质点相连膜元上分布荷载的等效外力
薄膜质点内力和外力具体推导过程参见文献[16]。
对于膜面变形中的几何非线性因素,关键问题是如何正确处理内力计算中不同构形上的刚体位移和变形。如1.3节所述,有限质点法采用“虚拟运动”的方式将转动做近似扣除,降低了刚体运动对内力求解的影响,将大位移大转动问题转化为小位移小转动问题,并实现参考构形的正确转变。其中,逆向运动后残余转动为小量,可以直接取Va(1a-2a-3a)为参考构形用微应变和工程应力求解变形和内力。节点力是一组平衡力,经过正向运动后平衡力的值不会改变,仅在方向上有差异,满足了内力的时间状态与空间状态相互对应的要求。
需要指出,有限质点法的内力计算不对结构的几何线性和非线性运动变形加以区别,都采取统一的模式进行处理,这使得几何非线性求解变成一种自然的过程。而采用显式时间积分法求解离散质点运动轨迹的二阶微分方程式(1),也不存在非线性迭代的概念。
有限质点法中通过途径单元的合理配置将内力求解简化为小变形的情况,材料因素在小变形情况下只会影响质点内力的计算,不修改计算的基本框架。程序中每一种材料模型都可以作为单独的计算模块,在需要的时候直接调用即可。而且,由于每个途径单元基础架构均进行更新,因此应变和应力都可以用全位移来计算,避免了使用复杂的变率量以及因无法采用合适的工程试验数据而造成的材料本构描述上的困难。针对膜材的具体做法为:在每个时间步根据单元的应力状态由膜材相应状态下的本构关系计算出应力增量Δσ,与途径单元初始时刻的应力σa叠加后代入虚功方程计算即得到该时间步的单元内力。因此采用有限质点法分析正交异性和非线弹性膜材的大变形问题不会在本质上增加求解难度,关键是本构模型的建立和材料参数的确定,本文3.1节和3.2节将做详细阐述。
一般当膜面仅发生大位移和小应变时,采用近似线性化应力-应变关系引起的误差很小,此时可假定膜材为线弹性体,在变形后仍保持正交异性。但当膜面几何发生较大改变后(一般认为应变超过2%为大变形),膜材的正交异性被破坏,这时必须遵照复合材料非线性本构理论建立膜材大变形模型,以尽可能准确描述膜材的应力-应变关系。
正交异性膜材在平面内有一对互相垂直的对称基准轴(即弹性主轴),其方向与经纬向纱线的排列一致。假定膜材为线弹性体,则根据广义虎克定律可建立织物膜材在径、纬两个主轴方向上的增量本构关系如下[1]:
式中:为膜材在弹性主轴方向上的应力-应变关系矩阵;Ew、Δσw、Δεw分别为经向的弹性模量、应力增量和应变增量;Ef、Δσf、Δεf分别为纬向的弹性模量、应力增量和应变增量;ΔσM和ΔεM分别为材料主轴方向上的应力增量和应变增量;Gwf、Δτwf、Δγwf分别为剪切模量、剪应力增量和剪应变增量;νwf、νfw分别为经向对纬向的泊松比和纬向对经向泊松比。
有限质点法中膜元的纯变形和内力计算都是在变形坐标系下进行的,即应变增量和应力增量
均按变形坐标分量描述,但变形坐标轴方向一般与弹性主轴方向并不一致,这就需要借助主、偏轴方向之间的应力、应变转换关系推导出变形坐标系下的本构关系。
对于正交异性膜材,有限质点法的内力计算中涉及到四种坐标系(如图3所示),分别为全域坐标系o-xyz、膜元局部坐标系o'-x'y'z'、膜元变形坐标系和材料弹性主轴坐标系
。其中,z'轴、
轴及
轴重合,且都沿着膜元的法线方向n或
故膜元内定义的三个坐标系始终位于同一平面内。
图3 膜元线弹性应力-应变计算相关的若干坐标系
Fig.3 The coordinates accounting for the stress-strain relationship defined in a membrane element
通过计算轴之间的方向余弦,可得到变形坐标向量
与弹性主轴坐标向量
之间的变换关系如下:
式中:为变形坐标与弹性主轴坐标的转换矩阵;
别代表
轴与弹性主轴坐标系各坐标轴的方向余弦。考虑到
所以有
和
由可得到变形坐标系
面内的应力增量
和应变增量
与材料弹性主轴坐标系
面内的应力增量
和应变增量
之间分别有转换关系如下:
式中:Tσ为应力变换矩阵;分别变形坐标系和弹性主轴坐标系下的膜元面内应力增量;Tε为应变变换矩阵,且有
分别变形坐标系和弹性主轴坐标系下的膜元面内应变增量。
假设变形坐标系下的膜材本构矩阵为,途径单元初始ta时刻膜面内应力为
,利用弹性主轴方向上的本构关系式
及式(6a)和式(6b),当前时刻t(ta≤t≤tb)以变形坐标分量描述的应力-应变关系可表示为:
需要指出,上述推导膜材在变形坐标下本构计算公式的过程中,确定膜元弹性主轴坐标与局部坐标之间的相对位置关系是一个非常关键的问题,同时也是膜结构大变形分析计算中的一个难点。以往许多文献回避了此问题或者没能提出合理的解决方案,本文这里给出一种实用计算方法。
图4中假设所示膜曲面上的虚线分别指向膜材的经、纬线方向,问题的关键是求出这两簇曲线方向向量与各个膜元的局部坐标向量(ex'、ey'、ez')之间的夹角φ。一般情况下,当膜结构的裁剪样式确定后,在膜片的边缘或角点位置上根据与边线的几何关系就能够确定相应的经、纬线主轴方向,因而至少有一个膜元内的弹性主轴方向是已知的。这里假设膜元1内弹性主轴方向已知,其经、纬线方向单位向量分别记作具体步骤如下:1) 找到膜元1的某一相邻膜元2,所处平面分别记作
和
(单位法线向量分别为N1和N2,交线向量记作i1=-i2如图5所示);2) 在
平面内按右手准则建立一对直角坐标系I1-J1-N1和I2-J2-N2,坐标轴方向分别沿着向量i1、j1、n1及i2、j2、n2,其中j1=n1×i1、j2=n2×i2;3)在由I1-J1所构成的坐标平面内沿
方向作一条直线,与i1、j1方向轴分别交于a点和b点,在i2轴和j2轴坐标方向上分别找出与a点和b点坐标值相同的c点和d点。这样,c-d连线方向就代表所要求的膜元2中的一个弹性主轴
从而实现材料主轴向相邻膜元的传递。
图4 膜材弹性主轴与局部坐标轴的关系
Fig.4 The relationship between the principal axes of elasticity and local axes defined in the membrane
图5 弹性主轴向量在膜元之间的传递方式
Fig.5 The transmission of the principal axes of elasticity between the membrane elements
织物膜材经历大变形后,经、纬向纱线之间的夹角将发生改变,线弹性条件下的正交异性本构模型不再适用。本文借助分段线性化拟合函数和应变能密度函数来近似描述膜材经纬向纱线的应力-应变关系;同时,根据膜材以纱线基层受力为主的特点,假设纱线处于单轴拉伸状态,分别按照膜材经纬向不同的拉伸性能建立等效单轴受拉本构模型,各项参数由膜材单、双轴试验数据确定。
假设当前时刻t薄膜元i内的经线和纬线构成一个平面主轴坐标系轴与局部坐标系x'轴之间的夹角分别为φ1和φ2,可确定
轴的方向向量
不一定正交),而参照式(7)可知
轴与变形坐标系的
轴之间的方向余弦分别为
图6 大变形后膜元内的经纬向坐标轴与应力、应变定义
Fig.6 The warp and fill coordinates with the stress and strain therein after large deformation in a membrane element
膜面内经纬向纱线可认为处于单向受拉状态,则变形坐标系平面内的应变增量
、应力增量
与经纬向应变增量
应力增量
之间有如下变换关系:
式中:为应变转换矩阵;
为应力转换矩阵;且有
适用于
坐标轴正交或非正交的情况[6]。
当前时刻t(ta≤t≤tb)膜材主轴方向上的应力-应变关系式可表示为将式(8)和式(9)代入后可得:
式中,分别为ta时刻变形坐标系下和沿经纬主轴方向的本构矩阵。
按照基本假设,经纬向纱线各自均处于单轴拉伸状态,故本构矩阵这样主轴方向上的应力-应变关系在形式上与单轴状态是一致的,因而可以通过定义两个主轴方向上的等效单轴应力-应变关系来建立膜材的非线性本构方向,即模型。若已知膜材在经纬线方向上的拉伸刚度函数分别为
则当前时刻t(ta≤t≤tb)主轴方向上的应力可写成如下矩阵形式:
需要注意的是,在计算本构矩阵时,其中的拉伸刚度函数
反映的是当前途径单元初始时刻(t=ta)沿膜材经纬向的单轴拉伸力学性质,其大小取决于ta时刻两个材料主轴方向上累积的线应变
但由于膜材经纬主轴方向及其夹角受剪切变形的影响会发生改变,因此当计算进入到下一个途径单元时(t=tb),应对主轴方向及其应变进行更新。关于主轴方向,令tb时刻
轴和
轴与局部坐标系x'、y'轴之间的方向余弦分别为
则由局部坐标系x'o'y'平面内的应变增量
与经纬主轴
方向上线应变增量
之间的换算关系可知夹角
分别满足以下关系,即:
式中,分别为当前途径单元初始时刻(t=ta)经、纬向材料主轴与局部坐标系的x'轴之间的夹角。
和
方向上的累积应变
更新过程如下:1) 根据主轴
与变形坐标轴
之间的方向余弦
确定转换矩阵
将ta时刻主轴方向上线应变
转至变形坐标下,加上ta-tb时段内的增量应变
2) 利用材料主轴向量
由式(8)将
改成以tb时刻的材料主轴坐标分量表示,得到更新后的累积线应变
,即
这里假设膜材在轴方向上都处于单向受力状态,因此
可直接按膜材沿经、纬线方向的应力-应变试验曲线来定义(包括线弹性、非线弹性以及非弹性)。织物膜材的试验曲线可以按照斜率大小进行分段,同一段内曲线上各点斜率基本保持在一个较小的变化范围内。基于此特点,本文首先考虑采用一种简单的多折线分段函数模型(piece-wise model,简称P-W模型)来定义
如图7所示,将膜材经、纬向应力-应变曲线以斜率转折点(1w/1f、2w/2f、3w/3f...、nw/nf)为界划分为n段(具体分段数需结合曲线斜率变化情况与计算精度要求来确定),转折点对应的应变值记作
在每一曲线段内分别进行线性拟合,利用最小二乘法求出拟合直线的斜率,记作
于是,经、纬向分段定义的拉伸刚度函数可表示为:
图7 经纬向应力-应变关系分段线性拟合
Fig.7 The piece-wise linear fitting for the stress-strain relationships in warp and fill directions
另外,在单调加载条件下,也可借助初始各向同性超弹性材料的本构模型来定义(如Murnaghan材料模型、Neo-Hooken模型、Mooney-Rivlin模型、Ogden模型等[18])。这类材料模型的基本特征是本构关系通过应变能密度函数W来表示,而W一般可表示为应变张量不变量I1、I2、I3的函数,基本公式为
本节通过几个简单膜面的几何、材料大变形算例,将该方法的计算结果与相关文献及试验结果进行对比,验证上述理论及所编程序的正确性。
下面以承受竖向荷载作用的马鞍面膜结构为例,考察本文方法求解各向同性和正交异性膜面大变形行为的能力。如图8所示,该马鞍形膜面投影形状为矩形,初始形状按以下函数确定,即:
图8 马鞍形膜结构的几何、边界与荷载条件
Fig.8 Geometry, boundary and loading conditions of the saddle-shaped membrane structure
本例中膜面投影边长取为a=b=1.0 m,1、3为低点,2、4为高点,h1=h3=0,h2=h4=0.3 m,四边简支固定。膜面预张力通过施加边界强迫位移(ux=uy=0.01 m)来实现,在竖向均布荷载(Pz=0 MPa~4.5 MPa)作用下产生向上变形,忽略重力影响。膜材厚度ha=0.003 m,密度ρ=1.05×103 kg/m3,材料本构考虑各向同性和正交异性两种情况:各向同性时弹性模量为E=6.0×103 MPa,泊松比ν=0.3;正交异性时经纬向弹性模量分别为Ew=6.0×103 MPa和Ef=2.5×103 MPa,泊松比分别为νwf=0.72和νfw=0.30,剪切模量Gwf=1.50×102 MPa。当膜材为正交异性时,膜面经纬向布置考虑两种;一种是与边界平行布置(经向沿x方向),另一种是沿对角线布置(经向沿y=x方向)。为检验算法的收敛性,计算时使用两种疏密程度的质点布置方式(沿边线方向质点分布分别为7×7和21×21),都分别用三节点(T3)和四节点(Q4)膜元进行连接,相应的计算模型记作T7×7、T21×21及Q7×7、Q21×21。
图9(a)和图9(b)分别为采用各向同性和正交异性模型(经纬向与边界平行)计算得到的结构中心点M的竖向位移随荷载变化情况。考虑膜材为各向同性情况,本文计算得到的M点荷载-位移曲线的变化趋势与Tsiatas等[21]采用比拟方程法(AEM)的分析结果基本一致;当膜材正交异性时,本文方法结果同样也与ANSYS的计算结果较为接近,两种膜面布置条件下的平均误差仅分别为7.16%、3.39%(T3)和4.97%、2.85%(Q4),由两种膜元都能获得较高的求解精度。
图9 马鞍形膜面中心点M的位移-荷载曲线
Fig.9 Load vs.displacement curve of particle M in the center
进一步考察膜材正交异性对结构整体的影响,图10给出了两种不同布置情况下由T21×21模型计算得到的Pz=4.5 MPa时膜面位移结果,并与各向同性膜材的位移分布情况进行比较。不难看出,考虑正交异性时膜面位移大小及分布与各向同性时相比都发生明显改变,总体上来看刚度较大的膜材经向截面(A-A截面与C-C截面)上的位移要小于刚度较小的膜材纬向截面(B-B截面与D-D截面),特别是当膜材经纬向与边界方向平行布置时,位移不再呈对称分布,膜面发生整体扭剪变形。可见,膜材的正交异性对其受力变形情况有较大影响,分析中应予以充分考虑。
图10 不同材料本构下的马鞍形膜面位移计算结果比较
Fig.10 The displacement comparison of the saddle-shaped membrane considering different constitutive laws
如图11所示,膜片计算模型的几何尺寸和技术参数参照文献[22―23]提供的试验资料确定,试验中采用PVC膜材。其中,单轴拉伸测试采用长条形试件,按照与膜材经向偏离θ角度剪取一组试样(θ=0°~90°,间隔为15°);双轴拉伸测试采用平面十字形切缝试件,分别按照θ=0°和45°进行取样,经纬向应力加载比例为1∶1,加载速率20 (N/m)/s。表1和表2分别列出了按经纬向应力-应变试验数据拟合得到的多折线模型(P-W模型)与Murnaghan模型(M模型)的非线性本构材料参数。计算中对单轴和双轴拉伸试件分别设置964个和1648个离散质点,并分别通过力控制和位移控制方式进行加载,以更接近试验中的加载情况。
图11 单轴与双轴拉伸试件模型
Fig.11 The axial and biaxial tensile test models
表1 P-W模型材料参数
Table 1 The material parameters of the P-W model
材料主轴Tξε˜ξ/MPa、T()()ηε˜η/MPa ε˜ξ、ε˜η范围单轴 双轴 单轴 双轴经向k1108 1932 0~0.0167 0~0.0104 1w k192 276 0.0167~0.1327 0.0104~0.1131 2w k485 1385 ≥0.1327 ≥0.1131 3w k113 192 0~0.0473 0~0.0342 1f纬向k276 462 0.0473~0.0761 0.0342~0.0644 2f k174 167 0.0761~0.2034 0.0644~0.1685 3f k291 875 ≥0.2034 ≥0.1685 4f
表2 Murnaghan模型材料参数
Table 2 The material parameters of the Murnaghan model
材料主轴C1/MPaC2/MPaC3/MPa C4/MPaC5/MPa经向单轴-62.3 4578.4 312.6-13084.2324.9双轴-108.9 6495.8 264.3-27987.4742.6纬向单轴12.8 1374.5 59.8-4873.1-95.2双轴-38.6 2916.9 31.5-3975.8-37.9
图12为本文计算得出的受拉试件应力-应变曲线与试验结果的比较。从图12可以看出,无论是沿膜材主轴方向还是偏轴方向进载,本文方法都能获得与试验曲线趋势一致的结果。其中,采用简单的P-W线性化模型的计算结果与试验值有较高的拟合优度(决定系数R2≥0.9)。图13还给出了十字形膜片在拉伸应力(σ =36 Ν/m)作用下的等效应变分布情况,可以看到中心正方形边缘区域由于所连接的各切条横向应变得到释放(各切条处于单向拉伸状态),导致该处等效应变与内部处于双向受力状态的核心区域相比应变值减小,而在悬臂切缝与核心区域的交界位置及角点附近由于应力集中效应膜材发生局部较大变形,这也与试验中试件破坏位置吻合。算例表明,本文方法通过引入各向异性非线性膜材本构模型可以正确计算典型薄膜试件的单、双向拉伸大变形行为。
图12 单轴和双轴受拉膜片应力-应变关系曲线
Fig.12 The stress-strain curves of the membrane in axial and biaxial tension
图13 双轴受拉膜片中等效应变分布(σ =40 Ν/m)
Fig.13 The Equivalent strain in the biaxial tensile membrane
本算例考察一个充气薄膜圆管在内压和集中力共同作用下的大变形行为(计入膜材的非线性拉伸特性)。薄膜圆管的几何尺寸与荷载条件如图14所示。圆管直径R=21 cm,膜材厚度h=1 cm,质量密度ρ=1.75 g/cm3,经纬向本构关系都按照Neo-Hooken超弹性非线性本构模型定义,其中Lamé常数μ均取为1.0 kg/cm2。薄膜管初始处于零应力状态,充气内压p从0开始缓慢递增至0.045 kg/cm2,然后在圆管外表面沿纵向施加集中线荷载F1,按照维持管内气体总量恒定(即pV=const)和压力恒定两种情况分别模拟结构响应,忽略重力影响。利用对称性条件,计算时可以只取圆管环向1/4部分。计算模型沿环向和纵向采用11×11、21×21、31×31三种质点布置形式对膜面进行离散,对应模型编号①~③。
图14 薄膜圆管的几何尺寸与荷载条件
Fig.14 Geometry and loading conditions of the membrane tube
对于采用Neo-Hooken材料的薄膜圆管,文献[24]给出了充气压力与环向变形关系的理论解:pR0/μh=1-γ-4,其中p为充气内压,γ=R/R0为环向拉伸变形量(R0和R分别为圆管充气前后的半径)。图15将本文方法计算得到的充气压力与环向变形之间的非线性关系与文献结果进行了比较,最大误差为9.4%,特别是采用细化模型③时可以获得与理论解非常接近的结果,显示出本文方法的精确性及良好的收敛性。图16给出了采用模型③计算得到的p=0.045 kg/cm2时圆管截面形状,同时也比较了集中线荷载作用下的膜面变形情况,与文献结果对比显示两者变化趋势保持一致。通过本算例的分析,验证了有限质点法模拟膜结构材料非线性变形行为的有效性。
图15 充气内压与环向变形量关系比较
Fig.15 Inflation vs.pressure of a tube-comparison
图16 集中线荷载作用下的充气圆管截面变形
Fig.16 Deformation of the fully inflated membrane tube subjected to a concentrated linear load
本文将限质点法应用于正交异性膜材大变形问题的研究。分别阐述了有限质点法处理膜材几何非线性和材料非线性因素的基本思路和计算方法,建立了建筑织物膜材的正交异性线弹性模型和经、纬主轴等效单轴受拉非线弹性模型,推导了相应的以变形坐标分量描述的本构矩阵,并提出了一种实用的材料弹性主轴确定方法,物理概念清晰,易于编程实现。通过一系列数值算例对本文方法进行了测试,并将计算结果和相关文献及试验结果进行了对比,表明该方法从理论到程序均能有效求解薄膜的几何和材料大变形问题,且具有较高的精度。
有限质点法突破了传统分析方法在结构非线性问题中遇到的困难,计算中独立求解各质点位移,不需要组集单元的刚度矩阵,也不需要迭代求解控制方程式,步骤清晰,通用性较好,在非线性问题分析中具有优势和灵活性。特别是对于线性与非线性问题不区别对待,都采用统一的计算模式处理,每种材料本构都可以作为单独模块在程序中引入,考虑膜材的正交异性或非线性均不会在本质上增加问题求解的难度。需要指出的是,本文关注的重点是膜结构的整体大变形分析方法,为了确保论述主题明确、避免行文冗长,计算中暂未涉及褶皱因素的影响,此部分工作将在后续论文中介绍。本文工作也为进一步研究薄膜结构的张拉成形、褶皱、破坏等更复杂的非线性行为打下了良好基础。
[1]杨庆山, 姜忆南.张拉索—膜结构分析与设计[M].北京: 科学出版社, 2004: 45―70.Yang Qingshan, Jiang Yinan.Analysis and design of tensioned cable-membrane structure [M].Beijing:Science Press, 2004: 45―70.(in Chinese)
[2]Kao R, Perrone N.Large deflections of axisymmetric circular membranes [J].International Journal of Solids and Structures, 1971, 7(12): 1601―1612.
[3]Onãte E, Flores F G, Marcipar J.Membrane structures formed by low pressure inflatable tubes: New analysis method and recent construction [C]// Onãte E, Kropelin B.Textile Composites and Inflatable Structures II.Netherlands: Springer, 2008: 163―194.
[4]Bletzinger K U, Linhard J, Wuchner R.Advanced numerical methods for the form finding and patterning of membrane structures [C]// Pimenta P D M, Wriggers P.New trends in thin structures: Formulaiton, Optimization and Coupled Problems.Vienna: Springer, 2010: 133―154.
[5]Zhang Y, Xu J, Zhang Q.Advances in mechanical properties of coated fabrics in civil engineering [J].Journal of Industrial Textiles, 2018, 48(1): 255―271.
[6]Kato S, Yoshino T, Minami H.Formulation of constitutive equations for fabric membranes based on the concept of fabric lattice model [J].Engineering Structures, 1999, 21(8): 691―708.
[7]Pargana J B, Lloyd-Smith D, Izzuddin B A.Advanced material model for coated fabrics used in tensioned fabric structures [J].Engineering Structures, 2007, 29(7):1323―1336.
[8]Reese S, Raible T, Wriggers P.Finite element modelling of orthotropic material behaviour in pneumatic membranes [J].International Journal of Solids and Structures, 2001, 38(52): 9525―9544.
[9]Gonçalve F R, Campello E M B.Orthotropic material models for the nonlinear analysis of structural membranes [J].Journal of the Brazilian Society of Mechanical Sciences and Engineering, 2014, 36(4):887―899.
[10]喻莹.基于有限质点法的空间钢结构连续倒塌破坏研究[D].杭州: 浙江大学, 2010.Yu Ying.Progressive collapse of space steel structures based on the finite particle method [D].Hangzhou:Zhejiang University, 2010.(in Chinese)
[11]Yu Y, Luo Y Z.Motion analysis of deployable structures based on the rod hinge element by the finite particle method [J].Proceedings of The Institution of Mechanical Engineers Part G-Journal of Aerospace Engineering,2009, 223(7): 156―171.
[12]俞锋, 尹雄, 罗尧治, 等.考虑接触点摩擦的索滑移行为分析[J].工程力学, 2017, 34(8): 42―50.Yu Feng, Yin Xiong, Luo Yaozhi, et al.Cable sliding analysis considering frictional effect [J].Engineering Mechanics, 2017, 34(8): 42―50.(in Chinese)
[13]张鹏飞, 罗尧治, 杨超.薄壳屈曲问题的有限质点法求解[J].工程力学, 2017, 34(2): 12―20.Zhang Pengfei, Luo Yaozhi, Yang Chao.Buckling analysis of thin shell using the finite particle method [J].Engineering Mechanics, 2017, 34(2): 12―20.(in Chinese)
[14]喻莹, 谭长波, 金林, 等.基于有限质点法的单层球面网壳强震作用下连续倒塌破坏研究[J].工程力学,2016, 33(5): 134―141.Yu Ying, Tan Changbo, Jin Lin, et al.Research on seismic progressive collapse of single-layer reticulated dome using the finite particle method [J].Engineering Mechanics, 2016, 33(5): 134―141.(in Chinese)
[15]罗尧治, 杨超.求解平面固体几何大变形问题的有限质点法[J].工程力学, 2013, 30(4): 260―268.Luo Yaozhi, Yang Chao.The finite particle method for solving geometric large deformation of planar solids [J].Engineering Mechanics, 2013, 30(4): 260―268.(in Chinese)
[16]Yang C, Shen Y, Luo Y.An efficient numerical shape analysis for light weight membrane structures [J].Journal of Zhejiang University Science A, 2014, 15(4): 255―271.
[17]郑延丰, 罗尧治.基于有限质点法的多尺度精细化分析[J].工程力学, 2016, 33(9): 21―29.Zheng Yanfeng, Luo Yaozhi.Multi-scale fine analysis based on the finite particle method [J].Engineering Mechanics, 2016, 33(9): 21―29.(in Chinese)
[18]Belytschko T, Liu W K, Moran B.Nonlinear Finite Elements for Continua and Structures [M].New York:John Wiley & Son, 2000: 544―549.
[19]Rezaiee-Pajand M, Kadkhodayan M, Alamatian J.A new method of fictitious viscous damping determination for the dynamic relaxation method [J].Computers &Structures, 2011, 89(9/10): 783―794.
[20]Oakley D R, Knight Jr N F.Adaptive dynamic relaxation algorithm for non-linear hyperelastic structures Part I.Formulation [J].Computer Method in Applied Mechanics and Engineering, 1995, 126(1/2): 67―89.
[21]Tsiatas G, Katsikadelis J.Large deflection analysis of elastic space membranes [J].International Journal For Numerical Methods in Engineering, 2006, 65(2): 264―294.
[22]易洪雷, 丁辛, 陈守辉.PES/PVC 膜材料拉伸性能的各向异性及破坏准则[J].复合材料学报, 2005, 22(6):98―102.Yi Honglei, Ding Xin, Chen shouhui.Orthotropic behavior and strength criterion of PES/PVC membrane materials under tensile loading [J].Acta Materiae Compositae Sinica, 2005, 22(6): 98―102.(in Chinese)
[23]Ambroziak A, Klosowski P.Mechanical testing of technical woven fabrics [J].Journal of Reinforced Plastics and Composites, 2013, 32(10): 726―739.
[24]Bonet J, Wood R D, Mahaney J.Finite element analysis of air supported membrane structures [J].Computer Methods in Applied Mechanics and Engineering, 2000,190(10): 579―595.
LARGE DEFORMATION ANALYSIS OF ORTHOTROPIC MEMBRANES USING THE FINITE PARTICLE METHOD
杨 超(1986―),男,浙江人,博士后,主要从事大跨度空间结构研究(E-mail: 04tmgcyc@zju.edu.cn);
郑延丰(1987―),男,福建人,博士后,主要从事大跨度空间结构研究(E-mail: yanfeng39@zju.edu.cn).