劉海峰,李正光

(1. 西安交通大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,西安 710049; 2. 吉林大學(xué) 數(shù)學(xué)學(xué)院,長春 130012)

結(jié)構(gòu)靜力重分析預(yù)處理共軛梯度法的一種有效改進(jìn)

劉海峰1,李正光2

(1. 西安交通大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,西安 710049; 2. 吉林大學(xué) 數(shù)學(xué)學(xué)院,長春 130012)

考慮預(yù)處理共軛梯度法在結(jié)構(gòu)靜力重分析中的改進(jìn). 利用增量剛度矩陣相對結(jié)構(gòu)修改后剛度矩陣稀疏的特點(diǎn),修改處理共軛梯度法的實(shí)施過程,以達(dá)到降低計(jì)算量、 節(jié)省計(jì)算時(shí)間的目的. 數(shù)值算例驗(yàn)證了該方法的有效性.

預(yù)處理共軛梯度法； 靜力重分析； 稀疏性； Cholesky分解

在飛行器設(shè)計(jì)、 橋梁隧道和能源與動(dòng)力等領(lǐng)域,有許多大型復(fù)雜的工程結(jié)構(gòu)需要進(jìn)行設(shè)計(jì). 結(jié)構(gòu)設(shè)計(jì)通常需要進(jìn)行多次修改,每次修改都需要對結(jié)構(gòu)進(jìn)行分析,計(jì)算量較大. 為改善這種情況,人們提出了結(jié)構(gòu)靜力重分析問題. 結(jié)構(gòu)靜力重分析的目的是盡可能多地利用初始分析的相關(guān)信息,設(shè)計(jì)高效算法計(jì)算修改后結(jié)構(gòu)的響應(yīng),以降低計(jì)算量,節(jié)省計(jì)算時(shí)間,加快結(jié)構(gòu)設(shè)計(jì)進(jìn)程[1].

目前,重分析方法主要分為兩類： 逼近重分析法和精確重分析法. 逼近重分析法又分為迭代逼近法、 組合逼近法、 多點(diǎn)逼近法和單點(diǎn)逼近法[2]. 其中最有效的逼近重分析法是迭代逼近法中的預(yù)處理共軛梯度法,該算法用給定的預(yù)處理矩陣作用結(jié)構(gòu)修改后的平衡方程,使系數(shù)矩陣的條件數(shù)大幅度減小,從而達(dá)到加快收斂的目的. 此外,預(yù)處理共軛梯度法還可以處理結(jié)構(gòu)節(jié)點(diǎn)數(shù)目發(fā)生變化的重分析問題[3]. 文獻(xiàn)[4]給出了預(yù)處理共軛梯度法在重分析問題中的最新研究進(jìn)展.

精確重分析法又稱為直接重分析法,在不考慮計(jì)算機(jī)舍入誤差的情況下,可計(jì)算出結(jié)構(gòu)修改后位移向量的精確解. 這類重分析方法一般都建立在Sherman-Morrison-Woodbury(SMW)低秩校正公式基礎(chǔ)上. 文獻(xiàn)[5]將直接重分析法應(yīng)用到裂縫增長計(jì)算中,取得了較好的效果.

本文考慮預(yù)處理共軛梯度法在結(jié)構(gòu)靜力重分析問題中的實(shí)施問題,利用增量剛度矩陣比結(jié)構(gòu)修改后的剛度矩陣更稀疏的特點(diǎn),修改了預(yù)處理共軛梯度法的執(zhí)行過程,從而降低了計(jì)算量,節(jié)省了計(jì)算時(shí)間.

1 問題的數(shù)學(xué)描述

設(shè)初始結(jié)構(gòu)自由度數(shù)目為m,給定其剛度矩陣K0∈m×m和荷載向量R0∈m,則位移向量x0可使用如下平衡方程：

計(jì)算,其中K0是一個(gè)對稱正定矩陣. 根據(jù)初始分析,它的Cholesky分解

假設(shè)結(jié)構(gòu)修改后,其剛度矩陣變?yōu)?/p>

其中ΔK稱為增量剛度矩陣. 則相應(yīng)的平衡方程變?yōu)?/p>

其中:R∈m表示修改后結(jié)構(gòu)的荷載; 剛度矩陣K∈m×m也是一個(gè)對稱正定矩陣. 靜力重分析的目的是盡量多地利用結(jié)構(gòu)修改前的初始信息計(jì)算修改后結(jié)構(gòu)的位移向量x,以降低計(jì)算量.

2 靜力重分析中的預(yù)處理共軛梯度法

在逼近重分析法中,預(yù)處理共軛梯度法是最有效的方法. 在結(jié)構(gòu)修改程度較小的情況下,選取初始結(jié)構(gòu)的剛度矩陣做預(yù)處理矩陣,使用少數(shù)幾步迭代即可得到精度較高的近似解.

算法1[6-8]

假設(shè)初始猜測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的計(jì)算量主要體現(xiàn)在兩個(gè)步驟： 解方程K0z=r和計(jì)算矩陣向量乘積w=Kp. 由于K0的Cholesky分解已知,使用前回代法和后回代法即可對K0z=r進(jìn)行求解. 下面研究w=Kp的計(jì)算問題.

在第k步迭代中,分別記w,p,z,r為wk,pk,zk,rk,則由式(3)在第k步迭代中,有

因此,w=Kp可寫成另一種形式：

w=r+ΔKz+βw.

由于在靜力重分析問題中,增量剛度矩陣ΔK的非零元素個(gè)數(shù)通常遠(yuǎn)少于矩陣K的非零元素個(gè)數(shù),所以ΔKz的計(jì)算量遠(yuǎn)小于w=Kp的計(jì)算量. 因此可將算法1的執(zhí)行過程進(jìn)行修改,得到如下算法.

算法2

假設(shè)初始猜測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稱為標(biāo)準(zhǔn)的預(yù)處理共軛梯度法,算法2稱為改進(jìn)的預(yù)處理共軛梯度法.

3 數(shù)值算例

下面用數(shù)值算例驗(yàn)證算法2的有效性. 電腦配置為Pentium 4,3.2 GHz CPU,2 Gb RAM. 使用軟件為MATLAB 7.2.0.232(R2006a).

考慮如圖1所示的橋體結(jié)構(gòu),橋長660 m,寬24 m,每個(gè)橋墩高24 m,6個(gè)橋塔均高出橋面36 m. 該橋可離散成一個(gè)具有4 420個(gè)節(jié)點(diǎn)和4 626個(gè)單元的有限單元模型. 除18個(gè)約束節(jié)點(diǎn)外,每個(gè)節(jié)點(diǎn)有6個(gè)自由度,因此整個(gè)結(jié)構(gòu)的自由度為26 412. 結(jié)構(gòu)單元有3種類型： 索單元(180個(gè)),梁單元(486個(gè)),板單元(3 960個(gè)). 索和梁單元使用材料的彈性模量均為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. 橋面每個(gè)節(jié)點(diǎn)受垂直向下的力,大小為P=1×104N. 為加固此橋體結(jié)構(gòu),將索單元的橫截面積增加到0.045 239 m2,梁單元的橫截面變?yōu)?.2 m×1.2 m. 分別使用算法1和算法2對修改后的結(jié)構(gòu)進(jìn)行計(jì)算,相對誤差精度ε=10-6. 由算法1和算法2計(jì)算橋結(jié)構(gòu)的計(jì)算時(shí)間分別為1.095 313 s和0.814 062 5 s. 本文所提出方法的計(jì)算時(shí)間約為標(biāo)準(zhǔn)預(yù)處理共軛梯度法的75%.

圖1 橋結(jié)構(gòu)Fig.1 Structure of the bridge

綜上所述,本文對靜力重分析中預(yù)處理共軛梯度法的執(zhí)行問題進(jìn)行了研究,利用靜力重分析問題中增量剛度矩陣ΔK比結(jié)構(gòu)修改后剛度矩陣K稀疏的特點(diǎn),將預(yù)處理共軛梯度法的執(zhí)行過程進(jìn)行修改,節(jié)約了計(jì)算時(shí)間,保留了預(yù)處理共軛梯度法的全部優(yōu)點(diǎn). 同時(shí),由于本文提出的算法不需要存儲(chǔ)K,因此所需的存儲(chǔ)空間也小于標(biāo)準(zhǔn)預(yù)處理共軛梯度法.

[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]李正光,吳柏生. 關(guān)于結(jié)構(gòu)剛度修改的一個(gè)新算法 [J]. 吉林大學(xué)學(xué)報(bào)： 理學(xué)版,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.

(責(zé)任編輯： 趙立芹)

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—),男,漢族,博士,講師,從事計(jì)算數(shù)學(xué)的研究,E-mail： lhflkcc@126.com. 通信作者： 李正光(1975—),男,漢族,博士,教授,從事計(jì)算力學(xué)的研究,E-mail： lizg@jlu.edu.cn.

國家自然科學(xué)基金(批準(zhǔn)號(hào)： 51005096).

O342

A

1671-5489(2014)05-0971-04