共轭梯度法在膜结构分析中的应用

高振宇

1 概述

力密度法应用于膜结构分析时建立的方程组如下:

其中，C为结构的拓扑矩阵;Q为力密度矩阵;Px，Py和Pz均为节点荷载，具体推导详见文献[1][2]。随着结构网格划分的增加，总结点数变大，力密度法方程组阶数升高，膜结构分析的最终问题就归结为大型线性方程组的求解上。C矩阵本身是含有大量零元素的稀疏矩阵，并且本文通过推导证明，该方程组具有对称正定的特点，所以，本文提出可利用其方程组的结构特性进行算法优化，以提高计算程序的效率。

一般来说，解线性方程组 Ax=b有两种数值分析方法:1)直接法。但直接法涉及到系数矩阵的分解，当 A为大型稀疏矩阵时，因为所要求的分解因子是稠密的，直接法的效率不高甚至难以实现。2)迭代法。这类方法是生成近似解序列{x(k)}，矩阵 A只在矩阵向量乘法时才用到，而稀疏矩阵的乘法已相当成熟。数值方法的一个基本原则是:求解任一问题都应利用其方程组的结构特性，下面将介绍一种迭代解法——共轭梯度法，其优点正好体现于求解稀疏、对称及正定的方程组的效率上。

2 共轭梯度法在力密度方程组求解中的应用

2.1 力密度方程组的对称正定性

根据式(1)，设任意n为向量p，且力密度矩阵 Q为对角矩阵，其中对角元素 q(i，i)＞0，(i=1，n)，则:

所以方程组系数矩阵(CTQC)为对称正定矩阵。

2.2 共轭梯度法

表1 共轭梯度法与Gauss消去直接法找形结果对比表

表2 共轭梯度法与Gauss消去直接法计算耗时对比表

3 算例分析

根据以上算法，本文编制了相应的程序，并进行了算例分析。

本算例为马鞍形膜结构，采用矩形网格划分，划分间距为0.5 m，膜面预应力为1.0 kN/m，膜网内部力密度值为1.0 kN/m，边索力密度值为14.0 kN/m。结点总数为552，其中4个为固定结点，坐标分别为 99 000 001(12，16，4)，99 000 002(3，10，0)，99 000 003(12，1，4)，99 000 004(20，9，0)。单元总数 1 000，其中“T”单元总数为128。

本文编制的程序计算迭代166次，满足终止准则‖r‖2＜1×e-6。图1，图2分别给出结构网格划分、几何零状态时与找形后的平衡形状的对比，表1，表2分别给出了共轭梯度法和Gauss消去直接法的结点坐标计算结果及其计算耗时的比较。算例分析表明，本文算法精度高、收敛快，是求解力密度方程组的可靠、高效的数值分析方法，为力密度法应用于大型膜结构分析奠定了坚实的基础。

[1] H.J.Schek.The force density methodfor form finding and computation of general network[J].Computer Methods in Applied Mechanics and Engineer，1974(3):115-134.

[2] 王 勇，魏德敏，高振宇.张拉膜结构力密度法混合找形分析[J].计算力学学报，2005(10):629-632.

[3] 郑 超，张建海.预处理共轭梯度法在岩土工程有限元中的应用[J].岩石力学与工程学报，2007(7):2820-2825.

[4] G.H.戈卢布，C.F.范洛恩.矩阵计算[M].袁亚湘，译.北京:科学出版社，2001.

[5] 李庆扬，关 治，白峰杉.数值计算原理[M].第2版.北京:清华大学出版社，2001.