基于惯性理论的膜结构充气过程计算机仿真
发布时间:2021年12月10日 点击数:1997
膜结构是20世纪中期出现的一种新型空间结构形式, 它所能创造的丰富多彩的外部建筑造型是其他建筑结构所不具备的.由于其本身独有的质量轻、强度高、热膨胀系数小、耐太阳辐射等优点, 在世界各地得到广泛使用, 成为一种代表当今建筑技术和材料科学发展水平的新型结构.而充气膜结构作为较早出现的膜结构形式, 是以柔性结构体系来承受屋面风荷载和雪荷载等各种外荷载的作用[1].
对充气膜结构的研究, 主要集中于当膜结构充气完成后, 对其进行静力分析;而在充气膜充气完成前瘪气到充胀的充气过程中, 对膜结构的动态实时分析方面的研究, 则涉及很少.本文在运用杆件惯性理论[2]的基础上对充气造型逐步充气施工过程进行端面分析研究, 并通过程序对充气过程进行了实时模拟.
1 基本模型
考虑到膜结构在完全充气之前, 膜内无应变和应力, 故用铰接链杆机构模型 (图1) 来模拟充气膜结构的充气过程, 此机构模型将充气膜结构截面作基本假定: (1) 充气膜结构截面可划分为有限个二维短杆; (2) 所有二维杆都是刚性的; (3) 杆间连接均为铰接.
对杆单元而言, 不论恒载、活载还是内外气压, 常常是作用在杆单元上的分布荷载、均布荷载或集中荷载[3], 而由于杆件的微小性, 对这些非节点荷载一般可简化为杆件质心处的集中力和弯矩[4,5] (如图2所示) .如果能够实时地找到分布在膜面上的各个结点的受力状态和运动规律, 那么就能实时确定整个膜面的状态.通过铰接链杆机构模型的构造, 我们就可以运用杆件惯性理论来分析各微杆的运动, 从而达到模拟整个充气过程的目的.
2 基于杆件惯性理论的充气膜结构充气模拟
2.1 充气过程膜面运动状态分析
根据基本假定, 整个膜截面可以离散为有限个微小单元杆件, 每一微小单元杆件之间铰接, 可等效为多自由度的机构.
图2表示第i根杆件在某一时刻的运动状态.杆件两端作用力分别为: Nx i1, Ny i1, Nx i2, Ny i2.所有的外力都表示成质心处的等价力: 集中力Fx i, Fy i, 弯矩Mi以及自重mig.由这些力所产生的加速度表示成线加速度ax i, ay i和角加速度εi, 杆件的运动则利用质心处的线速度vx i, vy i和角速度ωi表示.
根据牛顿运动定律, 有式 (1)
miax i=Nx i1+Nx i2+Fx i
miay i=Ny i1+Ny i2+Fy i-mig
Jiεi=[ (Nx i1-Nx i2) sinφi+
(Ny i2-Ny i1) cosφi]li/2+Mi (1)
式中: mi为杆件质量;Ji为转动惯量;li为杆长;φi为杆的倾角.等式右边为结构内力和外力的合力, 写成矩阵形式:
⎡⎣⎢axiayiεi⎤⎦⎥=⎡⎣⎢⎢1/mi0lisinφi2Ji01/mi−licosφi2Ji1/mi0−lisinφi2Ji01/milicosφi2Ji⎤⎦⎥⎥⋅ ⎡⎣⎢⎢⎢⎢Nxi1Nyi1Nxi2Nyi2⎤⎦⎥⎥⎥⎥+⎡⎣⎢Fxi/miFyi/mi−gMi/Ji⎤⎦⎥ (2)[axiayiεi]=[1/mi01/mi001/mi01/milisinφi2Ji-licosφi2Ji-lisinφi2Jilicosφi2Ji]⋅[Νxi1Νyi1Νxi2Νyi2]+[Fxi/miFyi/mi-gΜi/Ji](2)
即为: [Ai]=[Mi]·[Ni]+[Fi] (3)
由于机构中力和运动的关系很复杂, 无法通过积分求出精确解, 因此采用数值解, 即求出杆件在各瞬时的加速度, 由此求出速度、位置, 从而求出杆件运动的轨迹.若要求出杆件的加速度, 必须先求出结构的内力, 为此需利用各杆件的运动协调条件.
2.2 充气过程膜面运动协调条件
膜面各杆件的两端都是相互连接的, 因此同一结点的各杆件端点具有相同的速度和位移, 如图3所示.杆端的位移推导如下: 在经过一个极短的瞬间Δt, 杆件质心处的运动线速度由vx i0, vy i0变为vx i, vy i, 角速度由ωi0变为ωi即:
⎧⎩⎨⎪⎪vxi=vxi0+axiΔtvyi=vyi0+ayiΔtωi=ωi0+εiΔt (4){vxi=vxi0+axiΔtvyi=vyi0+ayiΔtωi=ωi0+εiΔt(4)
而杆端线速度与质心速度的关系如下:
⎡⎣⎢⎢⎢⎢vxi1vyi1vxi2vyi2⎤⎦⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢10100101lisinφi/2−licosφi/2−lisinφi/2licosφi/2⎤⎦⎥⎥⎥⎥⎡⎣⎢vxivyiωi⎤⎦⎥ (5)[vxi1vyi1vxi2vyi2]=[10lisinφi/201-licosφi/210-lisinφi/201licosφi/2][vxivyiωi](5)
因此, 综合式 (2) , (4) , (5) 并考虑到杆端位移与杆端速度的关系 (如Δxi1=vx i1Δt等) 可得到杆端位移与杆端未知力的关系:
⎡⎣⎢⎢⎢⎢Δxi1Δyi1Δxi2Δyi2⎤⎦⎥⎥⎥⎥=(Δt)2⎡⎣⎢⎢⎢⎢10100101lisinφi/2−licosφi/2−lisinφi/2licosφi/2⎤⎦⎥⎥⎥⎥⋅ ⎡⎣⎢⎢1/mi0lisinφi2Ji01/mi−licosφi2Ji1/mi0−lisinφi2Ji01/milicosφi2Ji⎤⎦⎥⎥⋅ ⎡⎣⎢⎢⎢⎢Nxi1Nyi1Nxi2Nyi2⎤⎦⎥⎥⎥⎥+(Δt)2⋅⎡⎣⎢⎢⎢⎢⎢⎢⎢10100101lisinφi2−licosφi2−lisinφi2licosφi2⎤⎦⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢FximiFyimi−gMiJi⎤⎦⎥⎥⎥⎥+ Δt⎡⎣⎢⎢⎢⎢10100101lisinφi/2−licosφi/2−lisinφi/2licosφi/2⎤⎦⎥⎥⎥⎥⎡⎣⎢vxi0vyi0ωi0⎤⎦⎥ (6)[Δxi1Δyi1Δxi2Δyi2]=(Δt)2[10lisinφi/201-licosφi/210-lisinφi/201licosφi/2]⋅[1/mi01/mi001/mi01/milisinφi2Ji-licosφi2Ji-lisinφi2Jilicosφi2Ji]⋅[Νxi1Νyi1Νxi2Νyi2]+(Δt)2⋅[10lisinφi201-licosφi210-lisinφi201licosφi2][FximiFyimi-gΜiJi]+Δt[10lisinφi/201-licosφi/210-lisinφi/201licosφi/2][vxi0vyi0ωi0](6)
可记为: [Δi]=[MiA]·[Ni]+[MiB] (7)
整个膜面机构的运动协调条件如下:
(1) 结点处的杆件的速度、位移相同;如果在第i结点上有j, k两个杆件的1号杆端, 则位移为:
Δxj1=Δxk1, Δyj1=Δyk1 (8)
以此类推, 在铰支座处, 有
Δxi1=0, Δyi1=0 (9)
(2) 在1个结点上, 各杆端结构内力的合力为零, 即
∑Nx=∑Ny=0 (10)
设膜面离散为M个杆件, 这样共有4M个杆端力, 要解出这4M个未知力, 就要有4M个方程, 即要有4M个平衡条件.
另设膜面机构有n个结点, 其中n1个为不动结点.在各个结点上, 杆件的数量不相同, 若用ci表示第i结点处杆件根数, 由于1根杆件有2个端点, 因此, ∑i=1nci=2M∑i=1nci=2Μ
在每个支座结点, 有2个位移为零的条件, 在每个非支座结点处, 有2个合力为零的条件, 合起来为2n个条件.
在非支座结点i上, 有ci根杆件, 相应地有2 (ci-1) 个位移相等的条件, 因此, 整个机构上非支座结点处位移相等的条件总共有
∑i=1n−n12(ci−1)=2∑i=1n−n1ci−2(n−n1)= 4M−2n1−2n+2n1=4M−2n (11)∑i=1n-n12(ci-1)=2∑i=1n-n1ci-2(n-n1)=4Μ-2n1-2n+2n1=4Μ-2n(11)
加上前面支座结点处的2n个条件, 总共正好为4M个条件, 可解.
3 计算机模拟的仿真实现
在利用VC程序实现时, 有两处重点需要注意: 其一是将杆单元的动力矩阵中的元素对应置于机构的总动力矩阵中;其二是多阶稀疏矩阵方程的求解.
在构造总动力矩阵之前, 先将每个杆端进行编号, 对出现几个杆端都相交于同一结点的情况, 可以用两个数组kjd[ ]和arjd[ ]来解决 (其中: kjd[ ]中存放每个结点的杆端的个数, arjd[ ]中存放这些杆端的具体编号) , 具体的程序实现步骤如图4所示.
首先构造初始状态文本, 确定膜面的初始状态.其中包括膜截面线密度, 各个结点的编号、位置坐标, 各根杆件的起始编号, 外部支座约束点的对应结点编号, 膜结构所受的内外气压等荷载值;接下来根据初始状态文本确定各结点和杆件的位置并显示;再生成单元杆的单元动力矩阵, 然后根据机构的运动协调条件将单元动力矩阵中的每一元素对号入座放入总的机构动力矩阵中, 形成总动力矩阵;求解动力矩阵方程, 得到每根杆件的杆端未知力;最后将求得的各杆端未知力代入原来的单元杆的杆端线位移增量与杆端未知力的关系式, 确定每个结点的下一Δt的位置;依次迭代, 从而可以得到整个膜截面的动态运动过程.
4 模拟结果
一圆形充气膜结构, 完全充气时膜截面半径r=0.2 m, 膜弹性模量[6]Ex=Ey=0.8×109 N/m2 (充气过程中可认为Ex=Ey=0) , 泊松比为μ=0.4, 厚度为h=0.5mm, 充气过程加载内压p=0.01MPa, 将最上顶点固定, 对其充气, 迭代步长Δt=0.01s, 得到膜结构充气过程如图5所示.
其中利用该程序实现时, 在迭代次数为0, 550, 1 700, 2 500, 4 000次, 即充气时间分别为5.5, 1.7, 25, 40 s时的膜截面的状态如表1所示.
表1 膜截面随充气时间的运动情况 导出到EXCEL
Table 1Movement of membrane section followingwith inflating time
充气过程 |
(a) | (b) | (c) | (d) | (e) |
|
|||||
迭代步数/次 |
0 | 550 | 1 700 | 2 500 | 4 000 |
充气时间/s |
0 | 5.5 | 1.7 | 25.0 | 40.0 |
膜截面积/m2 |
0.500 0 | 0.605 2 | 0.107 5 | 0.118 8 | 0.122 8 |
|
|||||
|
5 小结
充气膜结构的充气过程是典型的小应变大变形的问题.运用常规的非线性有限元分析理论无法准确分析膜结构在完全充气前的受力状态和膜面形状, 而且操作起来非常困难.采用杆件惯性理论, 可以比较方便地模拟各种形状的充气膜的完全干瘪状态到完全充气状态的施工过程, 避免了一些现有分析软件中对充气前的膜结构无法分析的问题, 为充气膜的施工提供一些参考和帮助.