张拉膜结构作为一种新兴结构, 近20年以来在国内外广泛流行, 倍受推崇。尤其在大跨度的体育展览场馆, 更是展现了其不可替代的优越性。同时在国内城市化进程加速的过程中, 做为建筑美和力学美完美统一的膜结构小品也倍受青睐。

张拉膜结构的实现和传统结构有许多不同的地方:首先, 从建筑方案到结构体系的确立是方案设计和结构计算多次反复的过程;其次, 膜结构的找形和荷载分析是整个过程的关键, 找形和荷载分析的速度和准确性决定了方案的进展和成功与否。

传统的膜结构找形和荷载分析的方法有力密度法、动力松弛法和非线性有限元法等。力密度法和动力松弛法是简化的方法, 它们在大规模提高计算速度的同时在一定程度上牺牲了计算的精度。而非线性有限元法因为考虑了几何非线性, 在提高计算精度方面有明显的优势, 但是以牺牲计算速度为代价。考虑如何弃弊从利, 国内外研究者做了大量的工作, 比如论证了三种方法在理论方面是有限元理论的简化或者延伸[1];把力密度法的简化形式结合到非线性有限元法进行找形[2]。本文通过对动力松弛法和非线性有限元法的研究和工程百富策略白菜网, 提出综合这两种方法的混合法, 实现对膜结构的设计和分析。在大大缩短分析周期的同时, 保证了计算的精度。同时提出了保证混合法可行性的检验方法。

膜结构动力松弛法的核心内容是三角形膜单元和杆单元的内力等效方法和动力松弛法的概念[3]。

当用等应力来找形时, 如图1, 膜单元和杆单元之间有如下的内力等效关系:

TmLm=tσ2tanαm,m=i,j,kΤmLm=tσ2tanαm,m=i,j,k (1)

其中:Tm是杆的内力;Lm、αm分别是杆的长度和对应角的大小;t、σ分别是膜单元的厚度和应力。

设在运动的t时刻, 结构的i节点在坐标方向上有残余力:

Rtij=Mi⋅V?tijRijt=Μi⋅V?ijt (2)

其中:Rtijijt为t时刻作用在节点i上j方向的残余力;Mij是i节点的虚拟质量;V?V?tijijt是t时刻i节点在j方向上的加速度。

用中心差分形式表示, 则有:

V?tij=(Vt+Δt/2ij−Vt−Δt/2ij)/ΔtV?ijt=(Vijt+Δt/2-Vijt-Δt/2)/Δt (3)

则t+Δt/2时刻的速度为:

Vt+Δt/2ijijt+Δt/2=Vt−Δt/2ijijt-Δt/2+Δt/Mi·Rtij (4)

可以确定t+Δt/2时刻的i节点坐标是:

xt+Δtijijt+Δt=xtijijt+Δt·Vt+Δt/2ijijt+Δt/2 (5)

利用 (4) 式和 (5) 式可计算结构各节点在新时刻t+Δt的坐标值。而t+Δt时刻节点的残余力是:

Rt+Δt/2ijijt+Δt/2=Pt+Δtijijt+Δt+∑∀n∀n(Xkj−XijL)t+Δtn⋅Tt+Δtn (6)∀n∀n(Xkj-XijL)nt+Δt⋅Τnt+Δt(6)

其中:Pt+Δtijijt+Δt是t+Δt时刻节点i在j方向的外荷载分量;n是与节点i相连的杆单元数;Xt+Δtkjkjt+Δt是与节点i相邻的另一节点k在t+Δt时刻的坐标;Lt+Δtnnt+Δt是t+Δt时刻相对应杆的长度;Tt+Δtnnt+Δt是t+Δt时刻相对应杆的内力。

结构体系每一时刻对应的总动能有:

Ek=12∑Mi∑j=13V2ijEk=12∑Μi∑j=13Vij2 (7)

为保证计算收敛, 节点的虚拟质量按以下公式计算:

Mi=Δt22⋅SimaxΜi=Δt22⋅Simax (8)

其中Simax为节点i的最大可能刚度, 且有

Simax=∑m(TmLm)t+Δt+∑n(TnLn)t+ΔtSimax=∑m(ΤmLm)t+Δt+∑n(ΤnLn)t+Δt (9)

其中:Tm、Lm分别是交于i节点的由膜单元等代成的杆单元的内力和长度;Tn、Ln分别是交于节点的索单元的内力和长度。

找形的过程中如果体系的能量收敛到小于一个要求的小量, 就找到了符合要求的形面, 否者就要重新选择初始参数再次计算。

非线性有限元法是考虑膜结构的大位移特性, 并使用一定的非线性代数方程组的数值解法来对膜结构进行找形和荷载分析[4,5]。

对膜结构进行有限元分析时, 采用三节点三角形单元对膜面进行网格划分, 并且使用如图2的局部坐标系, 位移函数和节点坐标有如下关系:

{u}={1 x y}⋅[Φ]−1⋅{dx}{v}={1 x y}⋅[Φ]−1⋅{dy}{w}={1 x y}⋅[Φ]−1⋅{dz}{u}={1xy}⋅[Φ]-1⋅{dx}{v}={1xy}⋅[Φ]-1⋅{dy}{w}={1xy}⋅[Φ]-1⋅{dz}

(10)

由于膜单元呈现几何非线性, 膜单元的几何关系中应变需要考虑位移的二次项, 有

εx=∂u∂x+12{(∂u∂x)2+(∂v∂x)2+(∂w∂x)2}εy=∂v∂y+12{(∂u∂y)2+(∂v∂y)2+(∂w∂y)2}γxy=∂v∂x+∂u∂y+{∂u∂x∂u∂y+∂v∂x∂v∂y+∂w∂x∂w∂y} (11)εx=∂u∂x+12{(∂u∂x)2+(∂v∂x)2+(∂w∂x)2}εy=∂v∂y+12{(∂u∂y)2+(∂v∂y)2+(∂w∂y)2}γxy=∂v∂x+∂u∂y+{∂u∂x∂u∂y+∂v∂x∂v∂y+∂w∂x∂w∂y}(11)

将式 (10) 代入 (11) , 可整理得到应变的位移一次项和和二次项的表达式:

{ε}=[A]{d}+[12dTBTBd 12dTCTCd dTBTCd]T{ε}=[A]{d}+[12dΤBΤBd12dΤCΤCddΤBΤCd]Τ (12)

膜材由基布和涂层组成。涂层主要起保护和自洁的作用, 基布由纤维织成, 是主要承力物质。基布的力学性能是各向异性的。但由于织物结构有经纬向, 因此在计算时近似假设膜材是正交各向异性的, 应力应变关系为

{σ}=[D]·{ε} (13)

其中, [D]为弹性矩阵, 当局部坐标x, y与经纬向重合时:

[D]=⎡⎣⎢d11d210d12d22000d33⎤⎦⎥[D]=[d11d120d21d22000d33]

根据虚功原理, 得到虚功方程为:

∫∫∫V[ (σ(0)xx(0)+σx) δεx+ (σ(0)yy(0)+σy) δεy

+ (τ(0)xyxy(0)+τxy) δγxy]dV= (f (0) +f) Tδd (14)

其中, σ(0)xx(0)、σ(0)yy(0)、τ(0)xyxy(0)为膜单元的初始应力, f (0) 为节点上的初始力。

利用变分原理, 由式 (12) 可得

δεx=[A1]δd+[dTBTB]δd

δεy=[A2]δd+[dTBTB]δd+[dTCTC]δd

δγxy=[A3]δd+dT[CTB+BTC]δd (15)

将式 (15) 和式 (13) 代入方程 (14) , 整理, 用矩阵表示为:

[kE+kG]d=f-r (16)

其中, Am和hm为膜单元的三角形面积和厚度。

r=AmhmATσ (0) -f (0)

kE=Amhm[ATDA]

kG=Amhm[σ(0)xx(0)BTB+σ(0)yy(0)CTC+τ(0)xyxy(0) (BTC+CTB) ]

同理按照上面的思路, 同样可以得到局部坐标下两节点索单元的几何非线性模型。把膜单元和索单元的局部坐标系下的刚度矩阵转换到整体坐标系下并组合总的刚度矩阵得到以下平衡方程:

K (d) ·d=f (17)

然后可以用修正的Newton-Raphson法解这个非线性方程组获得相关的解。

当用非线性有限元法找形时, 膜和索取为无限弹性, 也就是膜单元和索单元的杨氏模量为0, 而初始的预应力在计算过程中保持不变, 从而求得和预应力相对应的结构形面;在荷载分析过程中, 膜和索的实际材料参数参与工作, 并且因为膜和索是不能受压力的材料, 当计算中主应力出现负值时, 需要相应的对应力和杨氏模量进行处理, 以防止刚度矩阵的奇异。

根据上面的理论阐述, 编制了相应的动力松弛法和非线性有限元的程序来对索膜结构进行分析。

用动力松弛法 (DRM) 找形时, 对计算结果进行三个参数的控制, 首先是能量的控制, 结构的总动能小于一个极小值;其次是节点残余力的控制以满足节点平衡;第三是节点变位的控制, 以保证节点的结果坐标达到要求的精度。用非线性有限元法 (FNEM) 计算时, 只对节点残余力进行控制以保证节点平衡。

算例1:用悬链面对两种方法的结果进行对比。悬链面的解析解为z=−a{ln(x2+y2−−−−−−√+x2+y2−a2−−−−−−−−−−√)−lna}+hz=-a{ln(x2+y2+x2+y2-a2)-lna}+h.图3为悬链面初始网格。其中a=10m, b=50m, h=22.9243m, 悬链面的上下端为固定约束。初始面取为如图3b的圆台面, 环向划分成18等分, 径向划分成6等分。使用FERRARI1202T2的膜材, 厚度为0.8mm, 经纬向预张力都为2.5kN/m。对两种方法都采用最大节点残余力Δfmax≤0.001N的控制。在配制是CPU为PIV1.7GHz, 内存为DDR256M的PC上进行运算。收敛后, 动力松弛法迭代496步, 资源占用CPU99%, 内存810k, 用时2秒。收敛后总动能Ek≤0.2×10-5J, 节点坐标精度达到0.4×10-5m。非线性有限元法迭代5150步, 资源占用CPU99%, 内存7480k, 用时27分30秒。收敛后的结果形面如图4。表1给出用两种方法计算后, 如图3a中对应节点的结果坐标。

表中1、7是约束点。由表1可知对同一个悬链面, 采用相同的初始面、初始参数和残余力精度控制, 找形完成后, 两种方法的对应节点的坐标在毫米的精度上是一样的。两种结果和解析解相比, 能控制在1%的误差内。

结合其他工程算例的比较发现:对同一个索膜结构, 如果初始网格和初始参数一样, 两种方法又采用相同的节点残余力的收敛标准, 找形完后节点坐标在要求的精度范围内能够达到相同的结果。但是动力松弛法的收敛步数只有有限元法的10%左右, 收敛时间更是降低到0.1%。既然如此, 我们可以把动力松弛法的找形结果直接用非线性有限元进行荷载分析, 从而在大大提高设计速度的同时, 确保了荷载分析的精确性, 这就是混合法的概念。

这是杭州某体育训练中心田径场看台的悬臂挑篷膜结构。采用张拉膜与空间钢管结构体系, 建筑平面图如图5。内圈、外圈为同心圆弧, 外圈半径223m, 弧长138.24m, 内圈半径207.5m, 弧长128.63m, 悬挑跨度为17.488m。挑篷共有18榀桁架, 每榀之间的角度不全相同。挑篷端部标高为16.121m, 悬挑端标高为25.096m。挑篷下端部通过通长的联系梁相连。

索膜结构部分使用FERRARI1202T2的膜材, 密度为1.05kg/m2, 经纬向预张力都为2.5kN/m。其中边索和脊索的直径Ø=0.014m, 单位长度索重1.08kg/m。找形时边索索力30kN, 脊索为造型索, 索力为2kN。考虑到每榀膜结构大小是不同的, 对整个膜结构进行初始网格划分, 然后用动力松弛法进行找形, 节点残余力控制Δfmax≤0.001N。初始网格划分时控制三角形网格的边长在1m左右。图6是其中一片膜结构用动力松弛法找形完的结果。找形完后整个结构共有2742个节点, 5168个膜单元, 629个索单元。

在荷载分析之前, 我们可以用称为零荷载检验的方法对动力松弛法的找形结果进行检查。也就是对结果形面在赋予实际的材料常数后, 用非线性有限元法进行只在给定预应力作用下的荷载分析, 观察预应力从分布和位移的情况。对找形结果用这种方法检查, 迭代3步后收敛, 有0.1mm的位移, 膜应力在2.499～2.501kN/m之间, 边索索力都是30.000kN, 脊索索力在1.998～2.001kN之间。所以结构基本上没有位移和应力重分布, 动力松弛法的找形结果可以直接用于非线性有限元的荷载分析。

荷载分析主要计算四种受力:自重、雪荷载和如图7所示沿径向的X向和-X向风荷载。环向因为整个膜结构起伏很小, 不计算这个方向的来风荷载。参考建筑结构荷载规范GB50009-2001中单跨单坡屋面和双跨双坡屋面的积雪分布系数的要求, 和单坡顶盖和前后纵墙半开屋盖的风荷载体型系数的要求, 并且参考日本膜构造建筑物构造设计手册中的风洞试验资料, 确定如图7所示三个不同的区间给定积雪分布系数μr和风荷载体型系数μs如表3。

风荷载的计算依据公式:

wk=βz×μs×μz×w0 (18)

其中:风振系数βz取1.8;风压高度变化系数μz取1.25;基本风压w0取0.4kN/m2。

设计雪载的计算依据公式:

sk=μr×s0 (19)

其中, 基本雪压s0取0.45kN/m2。

参考图6a, 表4给出了在各种荷载作用下的分析结果。

膜材的抗拉强度是σm=112kN/m, 膜的允许应力是σ≤σm8=14kN/mσ≤σm8=14kΝ/m, 上述各种荷载作用下的膜内应力都在许用应力的范围内。直径Ø14索的最小破断力104kN, 而荷载作用下最大索力是54.0kN, 满足要求。同时膜结构只在X向风荷载作用有很小比例面积的褶皱。膜结构的位移也很小。因此通过荷载分析, 膜结构的设计是满足要求的。

取如图5所示的边榀K9和中榀K21给出索膜部分在各种荷载下对钢结构的作用力, 见图8、图9。

由上面的反力图可以看出, 边榀曲梁设计主要考虑如何平衡Y向的侧向力, 曲梁接地支座的设计也是考虑如何平衡这个力所引起的弯矩, 尤其是连接边索处的端部力对弯矩起控制作用, Z向反力在边榀钢结构设计中基本可以忽略。雪荷载和-X向风荷载使这个Y向侧向力重分布, 端部的力变小;X向风荷载使结构承受风吸作用, 端部的Y向侧向力比自重情况下增加了5倍左右, 是边榀的最不利工况。中榀曲梁设计主要考虑如何平衡X向力, 这个力使曲梁受压, 曲梁的接地支座设计则考虑如何平衡这个力所引起的弯矩, 尤其是连接边索处的端部力对弯矩起控制作用, 而钢结构自重产生的弯矩能在一定程度上平衡这个弯矩, Y向反力在中榀钢结构设计中基本可以忽略。雪荷载和-X向风荷载基本上不引起X向反力的变化, 只是Z向反力有一定程度的变化, 而且在接近接地支座处, 对抗弯的设计影响不大。X向风荷载使结构承受风吸作用, 端部的X向反力比自重情况下增加了6倍左右, 是中榀的最不利工况。

以上采用先对膜结构进行荷载分析, 然后计算出反力给钢结构的分开设计的方法, 能够使钢结构的受力性状明确, 容易确定钢结构受力的重点部位。

膜结构找形分析时, 对相同的初始形面, 当动力松弛法和非线性有限元法采用相同的残余力控制标准时, 能够收敛获得符合工程分析要求的相同的最终形面。

索膜结构分析时, 用动力松弛法进行找形, 然后直接用非线性有限元法进行荷载分析的混合法, 能够弃弊从利, 在大大缩短设计周期的同时, 又保证分析的精度, 从而实现了工程的快速选型和下部结构的准确分析。

考虑到实际工程的复杂性, 采用零荷载检验的方法可以对动力松弛法的找形结果用有限元法进行荷载分析可靠性的有效检查。