錢向東吳有奇林荔珊

(1.河海大學(xué)力學(xué)與材料學(xué)院,江蘇南京 210098;2.上海市水務(wù)規(guī)劃設(shè)計(jì)研究院,中國上海 200232)

下限極限分析的子迭代路徑跟蹤內(nèi)點(diǎn)算法

錢向東1,吳有奇1,林荔珊2

(1.河海大學(xué)力學(xué)與材料學(xué)院,江蘇南京 210098;2.上海市水務(wù)規(guī)劃設(shè)計(jì)研究院,中國上海 200232)

采用路徑跟蹤內(nèi)點(diǎn)法求解有限元下限極限分析所對應(yīng)的非線性規(guī)劃問題。在非線性方程組的Newton算法中引入子迭代過程,能夠直接采用位移型有限元的數(shù)據(jù)存儲格式和求解工具,并且大量計(jì)算可以在單元一級完成。改進(jìn)的算法可直接利用現(xiàn)有的位移型有限元程序,實(shí)現(xiàn)過程簡單。算例表明,該算法的效率和精度均可以得到保證。

下限極限分析法;有限元;非線性規(guī)劃;路徑跟蹤內(nèi)點(diǎn)法

下限極限分析法的目標(biāo)是尋求滿足靜力容許條件的最大外荷載,即極限荷載的最大下限[1]。當(dāng)采用有限單元法進(jìn)行求解時,下限極限分析問題最終轉(zhuǎn)化為以荷載因子為目標(biāo)函數(shù)、應(yīng)力向量為變量,滿足平衡條件和屈服條件的約束最優(yōu)化問題。其數(shù)學(xué)表達(dá)式為

式中:β——可變荷載的比例因子;E——系統(tǒng)的平衡矩陣;σ——系統(tǒng)的應(yīng)力向量;Rv——系統(tǒng)的可變荷載向量;Rc——系統(tǒng)的不變荷載向量;fi(σ)——應(yīng)力屈服函數(shù);NF——系統(tǒng)的應(yīng)力表征點(diǎn)數(shù)。

當(dāng)屈服函數(shù)是應(yīng)力的非線性函數(shù)時,最優(yōu)化問題(1)為非線性約束最優(yōu)化問題,早期的解法都是將約束條件線性化[2-3],繼而轉(zhuǎn)化為線性最優(yōu)化問題。對于大型有限元問題,這種線性化可能導(dǎo)致約束數(shù)量的急劇增加。近年來,有關(guān)學(xué)者為了避免因線性化而導(dǎo)致約束數(shù)量急劇增加的問題,都在致力于尋求直接求解非線性最優(yōu)化問題。Herskovits等[4-5]、Zouain等[6]提出一種適用于上、下限分析的二步可行方向迭代算法,由滿足全部等式約束條件的初始點(diǎn)出發(fā),采用Newton迭代和方向偏轉(zhuǎn)技術(shù),搜索可行方向確定滿足不等式約束的前進(jìn)方向和步長。Lyamin等[7]進(jìn)一步發(fā)展此方法,以可行點(diǎn)作為初始點(diǎn),提出了新的方向偏轉(zhuǎn)方法和分步Newton算法,并應(yīng)用于處理二、三維土力學(xué)極限分析問題。Krabbenhoft等[8]采用松弛變量和罰函數(shù)技術(shù)處理不等式約束,提出了一種適合于所有單元模型和屈服條件的路徑跟蹤內(nèi)點(diǎn)法[8-12],該方法不需要調(diào)整迭代方向,通過調(diào)整懲罰因子保證迭代點(diǎn)接近約束邊界而不超越邊界。這些方法最終都?xì)w結(jié)為非線性方程組的分步Newton算法,每步迭代中必須求解一個代數(shù)方程組,其系數(shù)矩陣的階數(shù)和稀疏性均不同于剛度矩陣,解法的效率有待進(jìn)一步研究,且很難利用現(xiàn)有的位移型有限元軟件。

本文將Krabbenhoft等[8]的算法加以改造,在分步Newton算法中引入子迭代過程,使其能夠直接采用位移型有限元的數(shù)據(jù)存儲和求解格式。

1 路徑跟蹤內(nèi)點(diǎn)法

引入松弛變量s≥0和對數(shù)懲罰函數(shù)后,下限極限分析的約束最優(yōu)化問題可以等價地表示為[8]

式中:μ——懲罰因子,采用迭代法計(jì)算時,其數(shù)值在迭代過程中逐漸減少。

由Kuhn-Tucker條件可導(dǎo)出問題(2)的最優(yōu)化條件,相應(yīng)的Lagrange函數(shù)為

式中:λ、v——Lagrange乘子向量和對偶變量。

一階Kuhn-Tucker極值條件要求Lagrange函數(shù)的梯度向量等于0,即

將式(3)代入式(4)后得

式中:e——單位向量。

采用Newton法求解非線性方程組(5),迭代格式為

式中:I——單位矩陣;H——Hessian矩陣,其表達(dá)式為

對于某些屈服準(zhǔn)則,可能會出現(xiàn)H非正定的情況,此時可引入一微小的正的攝動參數(shù)ε,把H改寫為

一般取ε=10-5maxλj,其中δij=1(i=j),δij=0(i≠j)。

以上就是下限極限問題(1)的路徑跟蹤內(nèi)點(diǎn)法。

2 分步Newton算法

考慮到矩陣T的高度稀疏性,為了減少方程組的規(guī)模,以達(dá)到提高計(jì)算效率的目的,展開式(6)中第一式,略去迭代次數(shù)標(biāo)記k,可以得到以下5個方程組:

由式(12)和式(11)有

將式(14)、式(15)代入式(9)有

將式(16)代入式(13)有

將式(17)代入式(10)可得求解Δβ的方程:

從而可以分步解出

分析比較表明,這種分步Newton算法一般比用式(6)直接求解約快10倍[7-8]。

3 子迭代算法

在分步Newton算法中,每一迭代步都必須對矩陣W、EW-1ET求逆。W是塊對角矩陣,求逆過程計(jì)算量不大;但EW-1ET的稀疏性、有效的存儲格式和求逆算法目前還沒有專門的研究。本文提出一種避免求逆矩陣(EW-1ET)-1且能夠直接利用有限元剛度矩陣的子迭代算法。

若取K0為通常意義下的彈性剛度矩陣,則可以構(gòu)造一種子迭代格式計(jì)算式(19)和(20)。繼續(xù)令g=EW-1q+r5,則式(19)為

再令δ=K-1Rv、φ=K-1g,則有Kδ=Rv、Kφ=g,將式(24)代入后有:

采用簡單迭代法可以求解方程(27)和(28),取δ0=0,φ0=0,有

解出δ和φ后,再代入式(26)有

把式(31)代入式(20)后,可以迭代求解Δv,取Δv0=0:

W-1是一個塊對角矩陣,存儲量很少,而K0y和ETy的計(jì)算均可在單元一級完成,因此,不必存儲整體的E矩陣。

4 算 例

如圖1所示,一個正方形板,邊長為2L,上下兩邊均有一長為L/2的切口,左右兩邊受均布荷載σ的作用。該模型材料為各向同性的均質(zhì)材料,服從Von Mises屈服條件,平面應(yīng)變情況下為

取σ0=3 MPa、L=1 m時的極限荷載為σc≈1.13156 MPa[8]。

考慮問題的對稱性,取1/4建立有限元計(jì)算模型,采用4結(jié)點(diǎn)四邊形等參單元模擬。4等分(N=4)的網(wǎng)格如圖2所示。采用2等分(N=2)至38等分(N=38)的10種網(wǎng)格進(jìn)行計(jì)算,結(jié)果見表1。

圖1 切口方板Fig.1 Slotted square block

圖2 N=4時的有限元法網(wǎng)格Fig.2 Element discretization forN=4

Newton算法的控制條件為滿足如下關(guān)系式時迭代結(jié)束:

5 結(jié) 語

極限分析對應(yīng)的非線性規(guī)劃問題是一類典型的約束最優(yōu)化問題,Krabbenhoft等提出的路徑跟蹤內(nèi)點(diǎn)法適用于各種有限元離散和屈服條件,無需對非線性約束條件進(jìn)行線性化。本文對該算法加以改造,在Newton算法中引入子迭代過程,使其能夠直接采用位移型有限元的數(shù)據(jù)存儲格式和求解工具,大量計(jì)算可以在單元一級完成。子迭代算法可直接利用現(xiàn)有的位移型有限元程序,易于實(shí)現(xiàn)。算例表明,子迭代算法的效率和精度均可以得到保證。

表1 極限荷載計(jì)算結(jié)果Table1 Calculated results of limit load σcfor slotted square block

[1]徐秉業(yè),劉信聲.結(jié)構(gòu)塑性極限分析[M].北京:中國建筑工業(yè)出版社,1985:1-11.

[2]BOTTERO A,NEGRE R,PASTOR J,et al.Finite element method and limit analysis theory for soil mechanics problems[J].Computational Methods in Applied Mechanics and Engineering,1980,22:131-149.

[3]SLOAN S W.Lower bound limit analysis using finite elements and linear programming[J].International Journal for Numerical and Analytical Methods in Geomechanics,1988,12:61-77.

[4]HERSKOVITS J.A two-stage feasible directions algorithm for nonlinearly constrained optimization[J].Mathematical Programming,1986,36:19-38.

[5]HERSKOVITS J.Feasible direction interior-point technique for nonlinear optimization[J].Journal of Optimization Theory and Applications,1998,99(1):121-146.

[6]ZOUAIN N,HERSKOVITS J,BORGES L A,et al.An iterative algorithm for limit analysis with nonlinear yield functions[J].International Journal of Solids and Structures,1993,30(10):1397-1417.

[7]LYAMIN A V,SLOAN S W.Lower bound limit analysis using non-linear programming[J].International Journal for Numerical Methods in Engineering,2002,55:573-611.

[8]KRABBENHOFT K,DAMKIDLE L.A general non-linear optimization algorithm for lower bound limit analysis[J].International Journal for Numerical Methods in Engineering,2003,56:165-184.

[9]PASTOR F,LOUTE E.Solving limit analysis problems:an interior-point method[J].Communications in Numerical Methods in Engineering,2005,11:631-642.

[10]SIMON J W,WEICHERT D.Numerical lower bound shakedown analysis of engineering structures[J].Computer Methods in Applied Mechanics and Engineering,2011,41-44:2828-2839.

[11]LI Huaxiang.A nonlinear programming approach to limit analysis of non-associated plastic flow materials[J].Mathematics and Mechanics of Solids,2013,18:524-542.

[12]SIMON J W,H?WERB D,WEICHERTB D.A starting-point strategy for interior-point algorithms for shakedown analysis of engineering structures[J].Engineering Optimization,2014,5:648-668.

Subiterative interior point path-following algorithm for lower bound limit analysis

QIAN Xiangdong1,WU Youqi1,LIN Lishan2
(1.College of Mechanics and Materials,Hohai University,Nanjing 210098,China; 2.Shanghai Water Planning and Design Research Institute,Shanghai 200232,China)

The interior point path-following method was applied to lower bound finite element limit analysis of the nonlinear programming problem.By applying the subiterative process in the Newton algorithm to solving nonlinear equations,the data storage format and solver of the displacement-based finite element method can be used directly, and large amounts of calculation can be performed element by element.The improved algorithm can be simply implemented with the existing displacement-based finite element program.An example shows that the subiteration algorithm can obtain high efficiency and precision.

algorithm for lower bound limit analysis;finite element;nonlinear programming;interior point pathfollowing algorithm

TU311

:A

:1000-1980(2015)03-0244-05

10.3876/j.issn.1000-1980.2015.03.009

2015-01 24

國家自然科學(xué)基金重點(diǎn)項(xiàng)目(11132003)

錢向東(1963—),男,江蘇吳江人,教授,博士,主要從事工程力學(xué)、水工結(jié)構(gòu)工程研究。E-mail:xdqian@hhu.edu.cn