刘海峰,李正光

(1. 西安交通大学 数学与统计学院,西安 710049; 2. 吉林大学 数学学院,长春 130012)

结构静力重分析预处理共轭梯度法的一种有效改进

刘海峰1,李正光2

(1. 西安交通大学 数学与统计学院,西安 710049; 2. 吉林大学 数学学院,长春 130012)

考虑预处理共轭梯度法在结构静力重分析中的改进. 利用增量刚度矩阵相对结构修改后刚度矩阵稀疏的特点,修改处理共轭梯度法的实施过程,以达到降低计算量、 节省计算时间的目的. 数值算例验证了该方法的有效性.

预处理共轭梯度法； 静力重分析； 稀疏性； Cholesky分解

在飞行器设计、 桥梁隧道和能源与动力等领域,有许多大型复杂的工程结构需要进行设计. 结构设计通常需要进行多次修改,每次修改都需要对结构进行分析,计算量较大. 为改善这种情况,人们提出了结构静力重分析问题. 结构静力重分析的目的是尽可能多地利用初始分析的相关信息,设计高效算法计算修改后结构的响应,以降低计算量,节省计算时间,加快结构设计进程[1].

目前,重分析方法主要分为两类： 逼近重分析法和精确重分析法. 逼近重分析法又分为迭代逼近法、 组合逼近法、 多点逼近法和单点逼近法[2]. 其中最有效的逼近重分析法是迭代逼近法中的预处理共轭梯度法,该算法用给定的预处理矩阵作用结构修改后的平衡方程,使系数矩阵的条件数大幅度减小,从而达到加快收敛的目的. 此外,预处理共轭梯度法还可以处理结构节点数目发生变化的重分析问题[3]. 文献[4]给出了预处理共轭梯度法在重分析问题中的最新研究进展.

精确重分析法又称为直接重分析法,在不考虑计算机舍入误差的情况下,可计算出结构修改后位移向量的精确解. 这类重分析方法一般都建立在Sherman-Morrison-Woodbury(SMW)低秩校正公式基础上. 文献[5]将直接重分析法应用到裂缝增长计算中,取得了较好的效果.

本文考虑预处理共轭梯度法在结构静力重分析问题中的实施问题,利用增量刚度矩阵比结构修改后的刚度矩阵更稀疏的特点,修改了预处理共轭梯度法的执行过程,从而降低了计算量,节省了计算时间.

1 问题的数学描述

设初始结构自由度数目为m,给定其刚度矩阵K0∈m×m和荷载向量R0∈m,则位移向量x0可使用如下平衡方程：

计算,其中K0是一个对称正定矩阵. 根据初始分析,它的Cholesky分解

假设结构修改后,其刚度矩阵变为

其中ΔK称为增量刚度矩阵. 则相应的平衡方程变为

其中:R∈m表示修改后结构的荷载; 刚度矩阵K∈m×m也是一个对称正定矩阵. 静力重分析的目的是尽量多地利用结构修改前的初始信息计算修改后结构的位移向量x,以降低计算量.

2 静力重分析中的预处理共轭梯度法

在逼近重分析法中,预处理共轭梯度法是最有效的方法. 在结构修改程度较小的情况下,选取初始结构的刚度矩阵做预处理矩阵,使用少数几步迭代即可得到精度较高的近似解.

算法1[6-8]

假设初始猜测x=0

k=0

r=R

δ0=rTr

t=‖R‖2

SolveK0z=rforz

ρk+1=rTz

k=k+1

ifk=1

p=z

else

p=z+βp

end

w=Kp

x=x+αp

r=r-αw

δk=rTr

end

其中ε是事先给定的相对误差精度. 算法1的计算量主要体现在两个步骤： 解方程K0z=r和计算矩阵向量乘积w=Kp. 由于K0的Cholesky分解已知,使用前回代法和后回代法即可对K0z=r进行求解. 下面研究w=Kp的计算问题.

在第k步迭代中,分别记w,p,z,r为wk,pk,zk,rk,则由式(3)在第k步迭代中,有

因此,w=Kp可写成另一种形式：

w=r+ΔKz+βw.

由于在静力重分析问题中,增量刚度矩阵ΔK的非零元素个数通常远少于矩阵K的非零元素个数,所以ΔKz的计算量远小于w=Kp的计算量. 因此可将算法1的执行过程进行修改,得到如下算法.

算法2

假设初始猜测x=0

k=0

r=R

δ0=rTr

t=‖R‖2

SolveK0z=rforz

ρk+1=rTz

k=k+1

ifk=1

p=z

else

p=z+βp

end

w=r+ΔKz+βw

x=x+αp

r=r-αw

δk=rTr

end

算法1称为标准的预处理共轭梯度法,算法2称为改进的预处理共轭梯度法.

3 数值算例

下面用数值算例验证算法2的有效性. 电脑配置为Pentium 4,3.2 GHz CPU,2 Gb RAM. 使用软件为MATLAB 7.2.0.232(R2006a).

考虑如图1所示的桥体结构,桥长660 m,宽24 m,每个桥墩高24 m,6个桥塔均高出桥面36 m. 该桥可离散成一个具有4 420个节点和4 626个单元的有限单元模型. 除18个约束节点外,每个节点有6个自由度,因此整个结构的自由度为26 412. 结构单元有3种类型： 索单元(180个),梁单元(486个),板单元(3 960个). 索和梁单元使用材料的弹性模量均为E1=2×1011Pa,梁单元使用的材料Poisson比为υ1=0.3； 板单元使用材料的弹性模量为E2=3×1010Pa,Poisson比为υ2=0.2. 索单元的横截面积为0.031 416 m2,梁单元的横截面为1 m×1 m,板厚0.2 m. 桥面每个节点受垂直向下的力,大小为P=1×104N. 为加固此桥体结构,将索单元的横截面积增加到0.045 239 m2,梁单元的横截面变为1.2 m×1.2 m. 分别使用算法1和算法2对修改后的结构进行计算,相对误差精度ε=10-6. 由算法1和算法2计算桥结构的计算时间分别为1.095 313 s和0.814 062 5 s. 本文所提出方法的计算时间约为标准预处理共轭梯度法的75%.

图1 桥结构Fig.1 Structure of the bridge

综上所述,本文对静力重分析中预处理共轭梯度法的执行问题进行了研究,利用静力重分析问题中增量刚度矩阵ΔK比结构修改后刚度矩阵K稀疏的特点,将预处理共轭梯度法的执行过程进行修改,节约了计算时间,保留了预处理共轭梯度法的全部优点. 同时,由于本文提出的算法不需要存储K,因此所需的存储空间也小于标准预处理共轭梯度法.

[1]Kassim A M A,Topping B H V. Static Reanalysis: A Review [J]. Journal of Structural Engineering,1987,113(5): 1029-1045.

[2]Kirsch U. Reanalysis of Structures [M]. Dordrecht: Springer,2008: 121-158.

[3]Wu B S,Lim C W,Li Z G. A Finite Element Algorithm for Reanalysis of Structures with Added Degrees of Freedom [J]. Finite Elements in Analysis and Design,2004,40(13/14): 1791-1801.

[4]Voormeeren S,Rixen D. Updating Component Reduction Bases of Static and Vibration Modes Using Preconditioned Iterative Techniques [J]. Computer Methods in Applied Mechanics and Engineering,2013,253(1): 39-59.

[5]Pais M J,Yeralan S N,Davis T A,et al. An Exact Reanalysis Algorithm Using Incremental Cholesky Factorization and Its Application to Crack Growth Modeling [J]. International Journal for Numerical Methods in Engineering,2012,91(12): 1358-1364.

[6]Golub G H,Loan C F,van. Matrix Computations [M]. 3rd ed. Baltimore: The Johns Hopkins University Press,1996: 532-544.

[7]李正光,吴柏生. 关于结构刚度修改的一个新算法 [J]. 吉林大学学报： 理学版,2003,41(1): 36-39. (LI Zhengguang,WU Baisheng. A New Algorithm for Modifications of Structural Stiffness [J]. Journal of Jilin University: Science Edition,2003,41(1): 36-39.)

[8]Kirsch U,Kocvara M,Zowe J. Accurate Reanalysis of Structures by a Preconditioned Conjugate Gradient Method [J]. International Journal for Numerical Methods in Engineering,2002,55(2): 233-251.

(责任编辑： 赵立芹)

EfficientImprovementofPreconditionedConjugateGradientMethodforStructuralStaticReanalysis

LIU Haifeng1,LI Zhengguang2
(1.SchoolofMathematicsandStatistics,Xi’anJiaotongUniversity,Xi’an710049,China;
2.CollegeofMathematics,JilinUniversity,Changchun130012,China)

This paper deals with the implementation of preconditioned conjugate gradient method (PCG) for structural static reanalysis. The implementation of PCG is revised by exploring sparsity of incremental stiffness matrices so that the computational cost and time can be reduced. Numerical example is given to show the effectiveness of the method.

PCG method; static reanalysis; sparsity; Cholesky factorization

2013-12-27.

刘海峰(1982—),男,汉族,博士,讲师,从事计算数学的研究,E-mail： lhflkcc@126.com. 通信作者： 李正光(1975—),男,汉族,博士,教授,从事计算力学的研究,E-mail： lizg@jlu.edu.cn.

国家自然科学基金(批准号： 51005096).

O342

A

1671-5489(2014)05-0971-04