力密度法找形, 是通过给定索网的内力与杆长的比值, 即力密度值, 由静力平衡方程, 得到关于结点坐标的线性方程组, 从而将一个非线性的问题转变成线性问题求解[1,2,3,4]。如果在找形过程中, 考虑支承体系的作用, 结构就可能存在拉索和压杆单元, 这些单元是不能够用力密度来控制内力, 而必须由弹性变形来控制其内力。在用力密度法找形中, 同时考虑弹性的拉索和压杆, 这部分研究报道不多见。本文将采用力密度法理论, 部分单元由力密度控制, 部分单元由弹性控制的找形, 称之为力密度法混合找形。一般情况下找形, 膜网内部单元和边索单元由力密度控制, 起支撑作用的拉索和压杆单元由弹性控制。混合找形的特点是由力密度控制的单元采用线性求解, 由弹性控制的单元采用非线性求解, 两者之间通过连接的结点由位移协调条件和力的平衡条件, 迭代计算逐步求解。

本文论述了张拉膜结构力密度法混合找形的基本理论, 导出了相关的迭代计算公式, 据此编制了相应的计算软件;对工程实例进行了验算, 结果表明, 本文给出的计算结果与德国著名软件easy的计算结果相吻合。

设结构由m个单元n个结点组成, 单元结点坐标差为

uj=xi-xk, vj=yi-yk, wj=zi-zk (1)

其中 i和k分别为单元j的终端结点和始端结点。

单元结点坐标差向量的矩阵形式为

u=Cx, v=Cy, w=Cz (2)

其中 C为单元-结点拓扑关系矩阵, 其形式见文献[3]。

力密度法找形结点平衡方程可表示为

CTQCx=Px, CTQCy=Py, CTQCz=Pz (3)

其中 Px, Py和Pz为结点荷载, Q为力密度矩阵。

其中 Nj和lj分别为第j单元的内力和杆长, qj为第j单元的力密度。

如果单元全部由力密度控制, Q为常数, 其中自由结点数为nf, 固定结点数为ng, 则结点平衡方程为

[CTfCTg]Q[Cf Cg]{xfxg}={0Pxg} (5)[CfΤCgΤ]Q[CfCg]{xfxg}={0Ρxg}(5)

可展开得

CTfQCfxf=-CTfQCgxg (6)

CTgQCfxf+CTgQCgxg=Pxg (7)

其中 xf和xg为自由结点向量和固定结点向量, Cf为单元-自由结点拓扑关系矩阵, Cg为单元-固定结点拓扑关系矩阵。

xg为已知, 求解式 (6) 可得xf;再将xf代入式 (7) 可得到固定结点反力Pxg。

同理可由其他两个方向展开的结点平衡方程。

CTfQCfyf=-CTfQCgyg (8)

CTgQCfyf+CTgQCgyg=Pyg (9)

CTfQCfzf=-CTfQCgzg (10)

CTgQCfzf+CTgQCgzg=Pzg (11)

求解得到另外两个方向的坐标值及固定结点反力。

如果单元为弹性控制, Q不再为常数, 式 (3) 将是非线性方程。第j单元杆长:

lj=(xi−xk)2+(yi−yk)2+(zi−zk)2−−−−−−−−−−−−−−−−−−−−−−−−−−−√=(u2j+v2j+w2j)1/2 (12)lj=(xi-xk)2+(yi-yk)2+(zi-zk)2=(uj2+vj2+wj2)1/2(12)

可以写成向量形式或矩阵形式:

l= (u2+v2+w2) 1/2 (13)

L= (U2+V2+W2) 1/2 (14)

其中 l为单元杆长列向量, L为单元杆长对角矩阵。

力密度矩阵可表示为

Q=H (L-10-L-1) (15)

其中 L-10为单元原长倒数对角矩阵, H为单元弹性刚度对角矩阵, 其对角元素hj为j单元的抗拉刚度。

hj=EjAj (16)

因此, 方程 (3) 可写成下列非线性方程:

⎧⎩⎨⎪⎪⎪⎪Fx=CTH(L−10−L−1)Cx−Px=0Fy=CTH(L−10−L−1)Cy−Py=0Fz=CTH(L−10−L−1)Cz−Pz=0 (17){Fx=CΤΗ(L0-1-L-1)Cx-Ρx=0Fy=CΤΗ(L0-1-L-1)Cy-Ρy=0Fz=CΤΗ(L0-1-L-1)Cz-Ρz=0(17)

采用牛顿迭代法解上述非线性方程组, 对于第i次迭代

⎡⎣⎢ΔxiΔyiΔzi⎤⎦⎥=−K−1i⎡⎣⎢FxiFyiFzi⎤⎦⎥ (18)[ΔxiΔyiΔzi]=-Κi-1[FxiFyiFzi](18)

其中 Ki为第i次迭代的雅可比矩阵,

Ki=⎡⎣⎢CTCTCT⎤⎦⎥[HL−3U2+H(L−10−L−1)HL−3UVHL−3V2+H(L−10−L−1) HL−3UWHL−3VWHL−3W2+H(L−10−L−1)⎤⎦⎥iΚi=[CΤCΤCΤ][ΗL-3U2+Η(L0-1-L-1)ΗL-3UVΗL-3V2+Η(L0-1-L-1)ΗL-3UWΗL-3VWΗL-3W2+Η(L0-1-L-1)]i

⎡⎣⎢CCC⎤⎦⎥ (19)[CCC](19)

则结点坐标的迭代计算公式为

xi+1=xi+Δxi, yi+1=yi+Δyizi+1=zi+Δzi (20)

根据以上建立的计算公式, 本文编制了相应的计算软件, 并对两个工程实例进行了计算和分析。

算例1 马鞍形膜结构

广东中山全球通膜结构小品工程, 膜网格划分间距为0.2 m, 膜面预应力为1 kN/m, 则膜网内部力密度值为1 kN/m, 边索力密度值为14 kN/m。网格划分, 结点总数为846, 单元总数为1426, 其中, 力密度控制单元1414, 弹性控制单元12;迭代6次, 坐标增量误差小于8.72×10-7, 找形结果如图1所示。

表1和表2列出算例1的本文与easy软件找形计算结果的比较。其中表1为固定点反力计算结果;表2为弹性压杆与拉索力的计算结果。

算例2 伞形膜结构

本工程为广州中海名都膜结构小品工程, 膜网格划分径向角度为8°, 环向距离为0.5 m;膜面预应力为1.5 kN/m, 边索力密度值为30 kN/m, 中间边索为70 kN/m。网格划分, 结点总数为676, 单元总数为1198, 其中, 力密度控制单元914, 弹性控制单元284;迭代180次, 坐标增量误差小于6.55×10-6, 找形结果如图2所示。

表3和表4给出算例2的本文与easy软件找形计算结果的比较。其中表3为固定点支反力;表4为弹性压杆与拉索力的值。

由上述计算结果的分析对比可以看出, 本文计算结果与德国著名软件easy计算的结果, 二者误差不到千分之一, 足以证明本文推导的结果正确、可靠;因此, 本文编制的膜结构计算软件完全可以百富策略白菜网于工程设计, 其计算数据格式与easy软件完全兼容。