基于漸進(jìn)結(jié)構(gòu)優(yōu)化法的拱壩結(jié)構(gòu)曲線函數(shù)優(yōu)化

左中杰, 干 泉, 王運興, 翟傳明, 張華棟

( 1.中電投工程研究檢測評定中心， 北京 100840； 2.中機(jī)三勘巖土工程有限公司， 湖北 武漢 430030；3.中建(北京)國際設(shè)計顧問有限公司， 北京 100013)

人類修建拱壩有著悠久的歷史。近年來世界各國著重建設(shè)各種體形的拱壩，然而拱壩建設(shè)是一項大型的水利工程項目，投資巨大，耗時較長，為了減少拱壩建設(shè)費用，縮短建設(shè)周期，開展拱壩的優(yōu)化研究就變得具有現(xiàn)實意義。

盡管拱壩有著安全度高、承載能力強(qiáng)等諸多優(yōu)點，但拱壩的開裂現(xiàn)象幾乎和拱壩自身歷史一樣久遠(yuǎn)。國內(nèi)外由于拱壩開裂而失事或影響運行的事例不勝枚舉。表面性的裂縫不會對大壩的安全性能造成影響，但結(jié)構(gòu)性裂縫，對壩的使用和工程整體性能造成不可估量的破壞。于是對壩體的結(jié)構(gòu)優(yōu)化應(yīng)最大程度地降低裂縫發(fā)生的可能性。

拱壩優(yōu)化設(shè)計的主要方法是調(diào)整初始模型中的各個參數(shù)，使拱壩各個截面的受力趨于合理。當(dāng)前，對拱壩的優(yōu)化設(shè)計多變量、非線性等問題。以往的設(shè)計方法多采用準(zhǔn)則法和數(shù)學(xué)規(guī)劃法，目前多采用有限單元法，以及由此發(fā)展的拓?fù)鋬?yōu)化分析方法。

準(zhǔn)則法主要有強(qiáng)度準(zhǔn)則、剛度準(zhǔn)則、以及能量準(zhǔn)則等。具體還可以分為滿應(yīng)力準(zhǔn)則和滿應(yīng)能準(zhǔn)則。

數(shù)學(xué)規(guī)劃法分為線性規(guī)劃和非線性規(guī)化法兩種。拱壩優(yōu)化屬于非線性規(guī)劃，詳細(xì)可分為：列線性規(guī)劃法(SLP)；罰函數(shù)法；序列二次規(guī)劃法(SQP)三種類型。

有限單元法目前已發(fā)展成為拱壩體型優(yōu)化的常用方法?？梢跃C合考慮地基、地質(zhì)、地形對拱壩的影響計算壩體的變形和應(yīng)力，還可以模擬材料間復(fù)雜的本構(gòu)關(guān)系，以達(dá)到精確近似工程設(shè)計的目的。通過市場上大型有限元軟件進(jìn)行快速求解，并利用軟件的高級功能對其曲線模擬等，既提高了優(yōu)化效率，又能根據(jù)設(shè)計要求對拱進(jìn)行修正。

拓?fù)鋬?yōu)化是最近發(fā)展起來的一種優(yōu)化方法，其中主要有變密度法、均勻法等思想。主要是對拱壩以及基巖進(jìn)行簡化模擬，然后對拱壩給定一個初始體形，設(shè)定要優(yōu)化的變量、目標(biāo)函數(shù)、約束條件，進(jìn)行拓?fù)鋬?yōu)化設(shè)計。目前常用的有限元軟件有NSYS、MSCMARC等，都是比較優(yōu)秀的計算工具。

我國在拱壩建設(shè)取得顯著進(jìn)展：拱壩優(yōu)化不再以經(jīng)濟(jì)指標(biāo)為單目標(biāo)進(jìn)行優(yōu)化而是考慮應(yīng)力水平、失效概率、高應(yīng)力區(qū)范圍等安全性指標(biāo)的多目標(biāo)優(yōu)化；拱壩優(yōu)化已由早期的線性優(yōu)化已轉(zhuǎn)變到非線性分析、混合線性優(yōu)化及拓?fù)鋬?yōu)化等。

結(jié)構(gòu)拓?fù)鋬?yōu)化的主要研究對象是連續(xù)體結(jié)構(gòu)。連續(xù)體結(jié)構(gòu)的拓?fù)鋬?yōu)化設(shè)計出現(xiàn)了不同的求解體系。國內(nèi)學(xué)者主要主要集中在局部應(yīng)力約束下的強(qiáng)度拓?fù)鋬?yōu)化設(shè)計，而國外研究人員多集中研究全局體積約束下的拓?fù)鋬?yōu)化。

1992年，Xie和Steven提出了漸進(jìn)結(jié)構(gòu)優(yōu)化算法，主要原理在優(yōu)化區(qū)域以一定的刪除準(zhǔn)則逐步刪除低應(yīng)力單元，使得結(jié)構(gòu)在滿足一定的約束條件下，應(yīng)力狀態(tài)和目標(biāo)函數(shù)達(dá)到最優(yōu)。該法具有刪除單元后無法恢復(fù)的缺點；因此，Querin和Steven等對傳統(tǒng)的漸進(jìn)結(jié)構(gòu)優(yōu)化算法進(jìn)行了改進(jìn)，在刪除單元后應(yīng)力較高的區(qū)域增加單元，使得應(yīng)力集中現(xiàn)象降低；Kim和Kwak對空間優(yōu)化問題進(jìn)行了深入的探討；Huang和Xie提出基于雙方向漸進(jìn)結(jié)構(gòu)優(yōu)化的研究方法，對周期結(jié)構(gòu)和具有一種或多種材料的結(jié)構(gòu)的拓?fù)鋬?yōu)化問題進(jìn)行了拓展。

1 漸進(jìn)結(jié)構(gòu)優(yōu)化方法基本思想和優(yōu)化準(zhǔn)則

1.1 漸進(jìn)結(jié)構(gòu)優(yōu)化(ESO)方法

1993年，謝憶民和Steven提出了漸進(jìn)結(jié)構(gòu)優(yōu)化方法(Evolutionary Strueture Optimization，簡稱ESO法)。認(rèn)為在設(shè)計區(qū)域內(nèi)，結(jié)構(gòu)上不起作用的材料，即那些低應(yīng)力或低應(yīng)變能密度的材料是低效的，是可以去除的。之后榮見華將這種優(yōu)化方法應(yīng)用于包含應(yīng)力、位移(剛度)、動力學(xué)約束和臨界應(yīng)力領(lǐng)域，取得了很好的成果。

漸進(jìn)結(jié)構(gòu)優(yōu)化方法的基本原理是逐步刪除最不起作用的材料，剩余結(jié)構(gòu)逐步達(dá)到最優(yōu)解。在優(yōu)化迭代中，該方法采用對存在的材料單元，其材料數(shù)編號為非零常數(shù)，而對不存在的材料單元數(shù)編號為零(或使得模型乘以接近為零的系數(shù))。該方法容易利用現(xiàn)有有限元分析軟件，通過多次迭代得以實現(xiàn)其所表現(xiàn)的拓?fù)鋬?yōu)化計算，因此該方法具有較好的通用性。

1.2 基于單元生死功能的優(yōu)化準(zhǔn)則

對于靜力設(shè)計問題，ESO方法通常采用基于應(yīng)力準(zhǔn)則的優(yōu)化方法。對于各向同性材料，VonMises應(yīng)力是目前最為常用的應(yīng)用準(zhǔn)則。

對于平面應(yīng)力問題VonMises應(yīng)力σvm定義為：

(1)

式中σx和σy分別是x和y方向的正應(yīng)力；τxy為剪應(yīng)力

(2)

式中，RRi為當(dāng)前的材料刪除率；一般初始刪除率RR0是人為取定值。在當(dāng)前子迭代中采用相同的RRi值，反復(fù)進(jìn)行子迭代運算，同時刪除滿足條件的單元，直到進(jìn)入應(yīng)力或體積穩(wěn)定狀態(tài)，表明在當(dāng)前迭代步中，已經(jīng)不存在滿足條件需要刪除的單元。在下一次迭代運算中，引進(jìn)一個進(jìn)化率ER，這樣刪除率被改變，即：

RRi+1=RRi+ER，i=0,1,2…

(3)

在改進(jìn)的刪除率下，又有一些滿足刪除條件的單元，執(zhí)行新一次單元刪除循環(huán)，直到達(dá)到?jīng)]有單元被刪除，停止迭代達(dá)到新的穩(wěn)態(tài)，從而獲得一個最佳結(jié)構(gòu)。理想的最終結(jié)構(gòu)變成一個滿應(yīng)力的設(shè)計即結(jié)構(gòu)的每一點的材料應(yīng)力達(dá)到它的滿強(qiáng)度。

1.3 有限元分析基本流程

(1)選擇準(zhǔn)備設(shè)計區(qū)域，用有限元網(wǎng)格離散全部區(qū)域，將準(zhǔn)備優(yōu)化的區(qū)域定義不同的單元類型；

(2)對結(jié)構(gòu)施加荷載和邊界約束，選擇靜力分析的方法進(jìn)行靜力分析；

(4)將以上步驟重復(fù)進(jìn)行，直至(2)式無法滿足且收斂到設(shè)計要求為止。此時，對應(yīng)于RRi的結(jié)構(gòu)已達(dá)到穩(wěn)定狀態(tài)。為使優(yōu)化繼續(xù)進(jìn)行，引進(jìn)進(jìn)化率ER，從而下一穩(wěn)定狀態(tài)的刪除率修改為式(3)；

(5)重復(fù)步驟2～4步直至結(jié)構(gòu)滿足要求為止。

1.4 拓?fù)鋬?yōu)化軟件實現(xiàn)

本文基于ESO法，利用ANSYS的單元生死對厚拱平面進(jìn)行拓?fù)鋬?yōu)化設(shè)計，采用APDL語言二次開發(fā)，編制循環(huán)程序設(shè)定刪除率，將設(shè)計區(qū)域內(nèi)應(yīng)力相對低的單元或低應(yīng)變能密度的單元“殺死”。逐步“殺死”低效單元，最終得到滿足應(yīng)力要求及最有目標(biāo)的結(jié)構(gòu)模型。

2 厚拱拱壩結(jié)構(gòu)的拓?fù)鋬?yōu)化計算

2.1 拱壩設(shè)計資料介紹

本文在對某二維厚拱拱壩進(jìn)行拓?fù)鋬?yōu)化應(yīng)用的基礎(chǔ)上，其新意在于二次拱結(jié)構(gòu)形式的分析?；驹O(shè)計資料為：拱壩形式為單圓心拱壩，如圖1所示，圓心角θ=80°，外徑R=450 m，內(nèi)徑r=350 m，拱厚T=100 m。厚拱判斷條件：D=0.5(R+r)，T/D=0.25，滿足厚拱(T/D=0.2～0.3)要求。壩體作為均質(zhì)材料考慮，彈性模量Ec= 30.0 G Pa，泊松比取μc= 0.3，拱壩承受均布水壓力為350 kN/m2。

圖1 厚拱拱壩結(jié)構(gòu)初始尺寸

2.2 拱壩結(jié)構(gòu)拓?fù)鋬?yōu)化計算模型

本例中厚拱結(jié)構(gòu)采用ANSYS單元庫中8節(jié)點平面單元plane 82進(jìn)行平面應(yīng)力分析，將拱壩結(jié)構(gòu)離散，共劃分10509個節(jié)點和5122個單元，其計算模型如圖2所示。

圖2 厚拱拱壩結(jié)構(gòu)有限元計算模型

2.3 拱壩結(jié)構(gòu)拓?fù)鋬?yōu)化結(jié)果分析

基于Von Mises屈服準(zhǔn)則確定的Von Mises Stress(等效應(yīng)力)可以快速的確定模型中的最危險區(qū)域。

本文漸進(jìn)優(yōu)化方法以單元的等效應(yīng)力為依據(jù)，對滿足刪除條件的單元進(jìn)行選擇、刪除。對拱壩結(jié)構(gòu)初始模型進(jìn)行加載及施加位移約束后，得到相應(yīng)的等效應(yīng)力、主應(yīng)力計算結(jié)果見圖3～5，其中圖3為初始化模型后的單元等效應(yīng)力等值圖，圖4為其結(jié)構(gòu)單元第一主應(yīng)力等值圖，圖5為初始化模型的節(jié)點第一主應(yīng)力。通過對拱壩初始模型設(shè)定刪除率，取RRi=10%，ER=10%，經(jīng)過40次迭代計算后，結(jié)構(gòu)單元刪除為零，得到迭代優(yōu)化40次后的拓?fù)鋬?yōu)化結(jié)果云圖，此時拱壩結(jié)構(gòu)的體積為滿足刪除率的最小體積。其計算結(jié)果如圖6～8所示，文中云圖結(jié)果單位為kPa。

圖3 拱壩模型初始化單元等效應(yīng)力云圖

圖4 拱壩模型初始化單元第一主應(yīng)力云圖

圖5 拱壩模型初始化節(jié)點第一主應(yīng)力云圖

圖6 迭代40次單元等效應(yīng)力云圖

圖7 迭代40次單元第一主應(yīng)力云圖

圖8 迭代40次節(jié)點第一主應(yīng)力云圖

由上述云圖可以得出，拱壩在進(jìn)行迭代的過程中無論是單元等效應(yīng)力，還是單元和節(jié)點的最大主拉應(yīng)力都得到大幅降低，使得拱壩最大拉應(yīng)力處于混凝土的容許拉應(yīng)力狀態(tài)，結(jié)構(gòu)趨于安全。另外還可以得出水平拱圈各區(qū)域主拉應(yīng)力的大致分布情況，以及整體危險點位置。拱圈初始模型中發(fā)生在基巖相接觸的角點處的主拉應(yīng)力最大值，在優(yōu)化過程中混凝土的最大拉應(yīng)力σ1<1.4 MPa，整個水平拱圈較安全，優(yōu)化結(jié)果達(dá)到了預(yù)期的目標(biāo)。

從圖9～11可以計算出，拱壩結(jié)構(gòu)優(yōu)化過程中，結(jié)構(gòu)體積減少33%，單元等效應(yīng)力最大值約減小50%，節(jié)點第一主應(yīng)力最大值減小約為92%。

圖9 拱壩結(jié)構(gòu)體積隨迭代次數(shù)變化曲線

圖10 單元等效應(yīng)力最大值隨迭代次數(shù)變化曲線

圖11 節(jié)點第一主應(yīng)最大值隨迭代次數(shù)變化曲線

2.4 拱壩結(jié)構(gòu)二次拱擬合分析

為了驗證拱壩結(jié)構(gòu)拓?fù)鋬?yōu)化后，其內(nèi)部可能存在的二次拱，為此對優(yōu)化后的結(jié)構(gòu)進(jìn)行曲線擬合，得出二次拱結(jié)構(gòu)的曲線?，F(xiàn)利用Origin對拱壩拓?fù)鋬?yōu)化后結(jié)構(gòu)上的數(shù)據(jù)進(jìn)行曲線擬合，其拱壩上下游曲線方程分別為：

y=4.11×10-16x-2.35×10-3x2+454.68

(4)

y=1.11×1016x-2.48×10-3x2+385.71

(5)

上下游拓?fù)鋬?yōu)化曲線與擬合曲線如圖12所示，圖中縱、橫坐標(biāo)表示距離，單位為m。

圖12 拓?fù)鋬?yōu)化后二次拱曲線

通過擬合的兩條拋物線建立拓?fù)鋬?yōu)化后的二次拱模型，并對其加載及施加位移約束，分析其節(jié)點第一主應(yīng)力及第三主應(yīng)力，如圖13～14所示。

圖13 擬合結(jié)構(gòu)第一主應(yīng)力云圖

圖14 擬合結(jié)構(gòu)第三主應(yīng)力云圖

由云圖中數(shù)據(jù)可知，第一主應(yīng)力的最大值小于壩體材料的抗拉強(qiáng)度(1.4 MPa)，第三主應(yīng)力最大值為3.045 MPa，也小于14 MPa，因此該擬合二次拱壩的主壓應(yīng)力、主拉應(yīng)力滿足材料強(qiáng)度要求。由此可以得出在已經(jīng)開裂的拱壩中存在一定形式的二次拱滿足不使混凝土開裂的要求。

3 其他形式的拱壩合理受力軸線分析

上述拓?fù)鋬?yōu)化結(jié)果分析得出，單心圓拱壩的二次拱的結(jié)構(gòu)形式呈拋物線，本節(jié)在拓?fù)鋬?yōu)化的基礎(chǔ)上，分析厚拱拱壩結(jié)構(gòu)的合理受力形式，采用與單心圓拱壩相同的工程背景，重點分析三心圓形式，對類似拱壩結(jié)構(gòu)的設(shè)計具有借鑒和指導(dǎo)作用。

3.1 三心圓拱壩結(jié)構(gòu)受力分析及結(jié)構(gòu)優(yōu)化

如圖15所示，三心圓的凈寬為單心圓的弦長長度，凈高取250 m。該拱壩承受與單心圓相同荷載，利用Ansys有限元計算，計算分析結(jié)果見圖16。

圖15 三心圓拱壩結(jié)構(gòu)形式

圖16 三心圓拱應(yīng)力分析及優(yōu)化結(jié)構(gòu)的應(yīng)力

云圖中三心圓拱初始模型的節(jié)點第一主應(yīng)力最大值小于1.4 MPa。且經(jīng)結(jié)構(gòu)優(yōu)化后，結(jié)構(gòu)體積減小23.8%，且主拉應(yīng)力較小，不會造成拱壩產(chǎn)生結(jié)構(gòu)性裂縫，滿足結(jié)構(gòu)安全。這說明了三心圓受力軸線也是拱壩工程設(shè)計中需要考慮的結(jié)構(gòu)形式。

4 結(jié) 語

本文主要借助ESO法的基本思想以及對ANSYS單元生死功能的二次開發(fā)，應(yīng)用APDL語言編制循環(huán)程序設(shè)定刪除率，通過對單心圓厚拱拱壩結(jié)構(gòu)的拓?fù)鋬?yōu)化，優(yōu)化結(jié)果滿足安全要求，實現(xiàn)了材料的有效利用。從而表明將拓?fù)鋬?yōu)化方法用于拱壩的體形優(yōu)化設(shè)計中是可行的。同時驗證了拱壩結(jié)構(gòu)材料開裂后二次拱的存在，即在原始模型中存在一個二次拱使得在拱壩高應(yīng)力區(qū)域得到分散，內(nèi)力重分布，且二次拱的結(jié)構(gòu)形式多為拋物線。嘗試分析了三心圓拱壩的受力狀況。

[1] Xie Y M, Steven G P. A simple evolutionary procedure for structural optimization[J]. Computers & Structures， 1993, 49(5): 885-896.

[2] Huang X, Xie Y M. Optimal design of periodic structures using evolutionary topology optimization[J]. Structural & Multidisciplinary Optimization, 2008, 36(6): 597-606.

[3] Huang X, Xie Y M. Bi-directional evolutionary topology optimization of continuum structures with one or multiple materials[J].Computational Mechanics, 2009, 43(3):393-401.

[4] Bendsoe M P. Scale effects in the optimal design of a microstructuredmedium against buckling[J]. Int J Solids and Structures, 1990, (26): 725-741.

[5] 周春平，常錦昕. 基于單元生死功能的轉(zhuǎn)向架構(gòu)架拓?fù)鋬?yōu)化設(shè)計[J]. 計算機(jī)仿真, 2010, 27(5):267-271.

[6] 穆春燕. 拓?fù)鋬?yōu)化理論及其在拱壩優(yōu)化設(shè)計中的應(yīng)用[D]. 南京：河海大學(xué)，2006.

[7] 郭中澤. 基于單元特性改變的漸進(jìn)結(jié)構(gòu)拓?fù)鋬?yōu)化方法及應(yīng)用研究[D].北京：中國工程物理研究所，2006.

[8] 汝乃畢. 拱壩開裂的歷史經(jīng)驗[J].水利學(xué)報，1990, (9)：17-25.

[9] 武 亮,葉文明,何仕華. 基于ANSYS的拱壩等效應(yīng)力分析[J]. 水利水電技術(shù)， 2006, 37(9)：27-29.

[10] 高 健.拱壩開裂約束下體形優(yōu)化設(shè)計研究[J].河海大學(xué)學(xué)報(自然科學(xué)版)，2006, 34(5)：526-529.

[11] 董桂西,王藏柱.結(jié)構(gòu)優(yōu)化設(shè)計的現(xiàn)狀及展望[J].電力情報, 2000，(1):5-7.