膜结构初始形态确定的研究
发布时间:2021年12月22日 点击数:2339
膜结构是近几十年来发展起来的一种采用新型材料的空间结构、全新建筑结构形式[1,2].膜结构设计内容和分析方法要比传统结构复杂.这种体系可在荷载作用下产生较大位移, 因而在计算中必须考虑几何非线性.膜结构的计算包括初始形态分析、受力分析和剪裁分析等内容, 它们构成所谓的膜结构全过程分析.膜结构的设计大致可按如下步骤进行:
(1) 总体设计:根据建筑的功能要求等确定膜结构的结构形式, 定出控制点的坐标, 选择膜材.
(2) 找形分析:根据边界条件及初始荷载要求确定相应的初始张力和结构形状.
(3) 荷载分析:验算结构在各种荷载下的变形及应力变化是否满足使用要求.
(4) 裁剪分析:根据膜结构形状, 确定将空间膜曲面用平面膜材表示的裁剪式样.
膜结构作为一种柔性体系, 仅当存在适当预应力时才具有确定的形状, 而且其几何形状是随支承一边界条件和预应力分布状态而变化的;因而膜结构设计的首要内容就是所谓的“形状生成”, 即初始形态分析, 借此来确定支承、边界条件、预应力状态、曲面形状这一综合系统在满足使用要求前提下的优化组合.
本文主要讨论的是张拉膜结构的找形分析.用于膜结构初始形态确定和荷载状态分析的主要方法有:力密度法、动力松弛法和非线性有限元法等.目前百富策略白菜网较广的方法是非线性有限元法, 这是由于非线性有限元法具有较高的计算精度, 并且随着计算机运行速度的大幅度提高, 也在一定程度上提高了该方法的计算速度.由于膜结构在荷载作用下一般处于大变形状态, 所以对于该类柔性结构的有限元计算需考虑结构具有的几何非线性的特点.
1 非线性有限元及计算方法的建立
1.1 位形、应变张量
分别用0xi, txi, t+△txi (i=1, 2, 3) [3] 描述物体内各点在时间0, t, t+△t的位形内的坐标.
左上标表示什么时刻物体的位形.物体位形的变化看作是从0xi, txi的一种数学上的一一对应变换.对于某一固定时刻, 这种变换表示成
txi=txi (0x1, 0x2, 0x3) (i=1, 2, 3) (1)
逆变换0xi=0xi (tx1, tx2, tx3)
(i=, 1, 2, 3)
dtxi=(∂txi∂0xj)d0xj‚d0xi=(∂0xi∂txj)dtxjdtxi=(∂txi∂0xj)d0xj‚d0xi=(∂0xi∂txj)dtxj
(2)
引用符号:t0xi, j=∂txi∂0xj‚0txi‚j=∂0xi∂txj∂txi∂0xj‚t0xi‚j=∂0xi∂txj
左下标表示该量对什么时刻位形的坐标求导数, 右下标“, ”前后的符号表示i量对j量求偏导数的坐标号.因此 (2) 式简写为:
dtxi=t0xi, jd0xj, d0xi=0txi, jdtxj (3)
任意P, Q两点之间在时刻0和时刻t的距离0ds, tds表示:
(0ds) 2=d0xid0xi=0txi, m0txi, ndtxmdtxn
(tds) 2=dtxidtxi=t0xi, mt0xi, nd0xmd0xn
变形表示
(tds) 2- (0ds) 2=t0xi, mt0xi, nd0xmd0xn
-0txi, m0txi, ndtxmdtxn
以变形后txi表示, 由 (3) 式:
(tds) 2- (0ds) 2
=t0xi, mt0xi, n0txi, mdtxi0txn, jdtxj
-0txi, m0txi, ndtxmdtxn
=t0xi, mt0xi, n0txi, mdtxi0txn, jdtxj
-0txi, m0txi, ndtxmdtxn
=t0xi, m0txm, it0xi, n0txn, jdtxidtxj
-0txi, mt0xi, ndtxmdtxn
= (δij-0txk, i0txk, j) dtxidtxj
=2ttttεijdtxidtxj
其中:
ttεij=12(δij−0txk‚i0txk‚j) (4)ttεij=12(δij-t0xk‚it0xk‚j)(4)
ttεij分别称为Almansi应变张量.
由 (3) 式:ttεij=0txk, i0txl, j0tεkl, 用tui, t+Δtui (i=1, 2, 3) 表示各质点在时间t, Δt 的位移, 有
txi=0xi+tui, t+Δtxi=0xi+t+Δtui
或tui=txi-0xi, t+Δtui=t+Δtxi-0xi (5)
求导数两种形式:
ttui, j=ttxi, j-0txi, j=δi, , j-0txi, j
t00tui, j=t00txi, j-0000xi, j=t00txi, j-δi, j
或0txi, j=δi, , j-ttui, j, 0txi, j=δi, j+t00tui, j, 代入 (4) 式
ttεij=12(ttui‚j+ttui‚j−ttuk‚ittuk‚j)ttεij=12(ttui‚j+ttui‚j-ttuk‚ittuk‚j)
1.2 应力
Cauchy应力张量以tτij表示, 代表真实应力. (又称欧拉应力张量) .以dT表示面力, 以dS表示面积.所以
dtTi=tτijtνjtdS
其中tνj是面积微元tdS上法线的方向余弦.
变形前后的面力的关系:
Kirchhoff规定:变形前后的面力与物体内各点的变换:
d0xi=0txi, jdtxj
以相同的规律相联系.即:
d0Ti=0txi, jdtTj
以此定义变形前的应力张量:
d0Ti=t0Si, j0νj0dS=0txi, , jdtTj
t0Sij称为第二类Piola—Kirchhoff应力张量.
左上标t表示应力张量是属于变形后 (时刻t) 位形的, 左下标0表示该量是在变形前位形 (时刻0) 内度量的.
为了得到tτij, t0Si, j的关系, 先确定tνjtdS和0νj0dS关系.由质量守恒定律可以得到[3]
t0xi‚γ0νj0dS=tρ0ρtνγtdS0txi‚γ0νj0dS=tρ0ρtνγtdS
所以:
d0Ti=t0Sji0νj0dS=0txi, αdtTα
=0txi, αtτα βtνβtdS
=0txi‚αtταβ0txj‚β0νj0dS0ρtρt0xi‚αtταβt0xj‚β0νj0dS0ρtρ
=0ρtρ0txi‚α0txj‚βtταβ0νj0dS0ρtρt0xi‚αt0xj‚βtταβ0νj0dS
t0Sji=0ρtρ0txi‚α0txj‚βtταβ0tSji=0ρtρt0xi‚αt0xj‚βtταβ
由 δij= 0txi, αt0xα, j
δij= t0xi, α0txα, j
tτji=tρ0ρt0xi‚αt0xj‚βt0Sβαtτji=tρ0ρ0txi‚α0txj‚β0tSβα
1.3 虚位移原理
由式 (5) 时间t到时间t+△t的位移增量为:
ui=t+△tui-tui
=t+△txi-0xi- (txi-0xi)
=t+△txi-txi
时间t+△t位形内物体的平衡条件及力边界条件相等效的虚位移原理为:
∫t+△tVt+△tτijδt+△teijt+△tdV=t+△tW (6)
t+△tW=∫t+△tSt+△tt+△ttkδukt+△tdS+
∫t+△tVt+△tρt+△tt+△tfkδukt+△tdV
t+△tW是时间t+△t位形的外载荷的虚功, t , f是面载荷, 体载荷.δt+△teij=δ12(t+△tui‚j+t+△tuj‚i)δt+△teij=δ12(t+△tui‚j+t+△tuj‚i)是无穷小应变的变分.以上是参考的时间t+△t位形.
更新拉格朗日格式 (Updated Lagrange Formulation 简称U.L. 格式) 是以时间t的位形作为参考位形.由于求解过程中参考位形是不断改变的, 固称为更新拉格朗日格式.
假定物体的面载荷, 体载荷是保守的, 不依赖于物体位形.则有:
t+△tt+△ttkt+△tdS=t+△tttktdS
即 t+△tt+△ttk=t+△tttk
t+△tρt+△tt+△tfkt+△tdV=t+△tρt+△ttftkdV
即 t+△tt+△tfk=t+△ttfk
t+△t tSij=tρt+△tρt+△tτrs tt+△txi‚r tt+△txj,stt+△tSij=tρt+△tρt+△tτrst+△ttxi‚rt+△ttxj,s
t+△ttεij=t+△terst+△ttxr, it+△ttxs, j
上几式代入 (6) 式转换为更新拉格朗日格式:
∫tVt+△ttSijδt+△ttεjitdV=t+△tW (6b)
引入下列增量分解:
应力增量:t+△ttSij=tτij+tSij
应变增量:t+△ttεij=tεij=teij+tηij
并且teij=12(tui‚j+tuj‚i)‚ tηij=12tuk‚ituk‚j并且teij=12(tui‚j+tuj‚i)‚tηij=12tuk‚ituk‚j
将上几式代入 (6b) 得:
∫tVtSijδtεijtdV+∫tVtτijδtηijtdV
=t+△tW-∫tVtτijδteijtdV (7)
(7) 式给出关于位移增量ui的非线性方程.
物理方程:tSij=Dijkltεkl
Dijkl 是线性本构张量.
应变增量:t+△ttεij=tεij=teij+tηij
代入得式 (7) :
∫tVtDijklteklδteijtdV+ (∫tVtDijklteklδtηij+tηijδtekl+tηijδtηijtdV) +∫tVtτijδtηijtdV
=t+△tW-∫tVtτijδteijtdV
略去扩号中的高阶项:
∫tVtDijklteklδteijtdV+∫tVtτijδtηijtdV
=t+△tW-∫tVtτijδteijtdV (8)
1.4 有限元方程
结点以小写字母k表示, n是单元的结点数.单元坐标和位移用结点值插值表示:
0xi=∑k=1nNk0xkitxi=∑k=1nNktxki0xi=∑k=1nΝk0xiktxi=∑k=1nΝktxik
t+△txi=∑k=1nNkt+△txki (i=1‚2‚3)t+△txi=∑k=1nΝkt+△txik(i=1‚2‚3)
tui=∑k=1nNktukitui=∑k=1nΝktuik
ui=∑k=1nNkuki (i=1‚2‚3) (9)ui=∑k=1nΝkuik(i=1‚2‚3)(9)
txki是结点k在时间t的i方向坐标分量, tuki是结点k在时间t的i方向位移分量, Nk 是结点k的形函数.
如果考虑一个单元, 以 (9) 式求 (8) 式中的位移导数项, 再集成, 导出矩阵方程:
(ttKL+ttKLN) u=t+△tQ-ttF (10)
其中:ttKL=∑e∫tνttBTLtDttBLtdVttΚL=∑e∫tνttBLΤtDttBLtdV
ttKNL=∑e∫tνttBTNLtτttBNLtdVttΚΝL=∑e∫tνttBΝLΤtτttBΝLtdV
ttF=∑e∫tνttBTLtτˆtdVttF=∑e∫tνttBLΤtτ^tdV
t+△tQ是结点载荷向量, ttBL和ttBNL分别是线性应变teij和非线性应变tηij与位移的转换矩阵, tD是材料本构矩阵, tτ和t是Cauchy应力矩阵和向量.∑e∑e单元集成.
采用修正Newton-Raphson迭代[3]式 (10) 为:
(ttKL+ttKNL ) △u (l) =t+△tQ-t+△tt+△tF (l) , (l=0, 1, 2, …) (11)
对于三角学单元可计算得:
ttBL=12A⎡⎣⎢y2−y30x3−x20x3−x2y2−y3000y3−y1 x1−x3 x1−x3y3−y1000y1−y2 x2−x1 x2−x1y1−y2000⎤⎦⎥ttBNL=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢y2−y3x3−x2000000y2−y3x3−x2000000y2−y3x3−x2y3−y1x1−x3000000y3−y1x1−x3000000y3−y1x1−x3y1−y2x2−x1000000y1−y2x2−x1000000y1−y2x2−x1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥ttBL=12A[y2-y300y3-y10y1-y200x3-x20x1-x30x2-x10x3-x2y2-y30x1-x3y3-y10x2-x1y1-y20]ttBΝL=[y2-y300y3-y100y1-y200x3-x200x1-x300x2-x1000y2-y300y3-y100y1-y200x3-x200x1-x300x2-x1000y2-y300y3-y100y1-y200x3-x200x1-x300x2-x1]
tτ=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢tσxtτxy0000tτxytσy000000tσxtτxy0000tτxytσy000000tσxtτxy0000tτxytσy⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥tτ=[tσxtτxy0000tτxytσy000000tσxtτxy0000tτxytσy000000tσxtτxy0000tτxytσy]
t= (tσxtσytτxy)
tD=⎡⎣⎢D11D210D12D22000D33⎤⎦⎥tD=[D11D120D21D22000D33]
D11=Ex/ (1-γxyγyx)
D22=Ey/ (1-γxyγyx)
D12=Eyγyx (1-γxyγyx)
D21=Exγxy (1-γxyγyx)
D22=Gxy
2 算例及结果分析
薄膜结构的设计中, 在找到初始形态之前, 并不能准确确定薄膜结构的初始形状和与之对应的预张力分布状态.也就是说, 这时有两个未知数:一是初始形状, 二是预张力分布状态.这时我们会给定其中的一个而求解另一个, 从而产生两种方法.第一种方法是将初始形状作为已知数, 把初始预张力作为外荷载施加在结构上, 求解达到平衡时的状态.此解法求解比较直接方便, 普通的静力计算程序就可以适用, 并可以得到与设计者给定的曲面形状相近的结果.第二种方法是将初始预张力的分布状态作为已知数, 而把与之对应的平衡的形状作为未知数来求解.此目的是要得到具有给定初始预张力分布的一个平衡的初始形态, 为了保证最终得到的预张力分布即为初始假定的预张力, 我们将舍弃变形协调条件和材料本构关系.也就是设计给定控制点的位置, 以及预张力的分布状态, 然后寻求在此条件下的与之对应的平衡的曲面形状.此时式 (11) 为:
ttKNL△u (l) =-t+△tt+△tF (l) ,
(l=1, 2, 3, …) (12)
算例1:图1 为一四点支承, 四边有索加劲的正方形平面膜, 膜厚为 1mm, 膜的预应力为20N/mm2, Et=2 550N/cm, Gt=800N/cm, γ=0.6 .四边索的预拉力为30 kN.用fortran语言编制相应的计算程序, matlab语言编制图形网格, 通过动态连接, 最后的平衡状态如图2, 图3所示.
