基于漸進結構優化法的拱壩結構曲線函數優化

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

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

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

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

拱壩優化設計的主要方法是調整初始模型中的各個參數，使拱壩各個截面的受力趨于合理。當前，對拱壩的優化設計多變量、非線性等問題。以往的設計方法多采用準則法和數學規劃法，目前多采用有限單元法，以及由此發展的拓撲優化分析方法。

準則法主要有強度準則、剛度準則、以及能量準則等。具體還可以分為滿應力準則和滿應能準則。

數學規劃法分為線性規劃和非線性規化法兩種。拱壩優化屬于非線性規劃，詳細可分為：列線性規劃法(SLP)；罰函數法；序列二次規劃法(SQP)三種類型。

有限單元法目前已發展成為拱壩體型優化的常用方法。可以綜合考慮地基、地質、地形對拱壩的影響計算壩體的變形和應力，還可以模擬材料間復雜的本構關系，以達到精確近似工程設計的目的。通過市場上大型有限元軟件進行快速求解，并利用軟件的高級功能對其曲線模擬等，既提高了優化效率，又能根據設計要求對拱進行修正。

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

我國在拱壩建設取得顯著進展：拱壩優化不再以經濟指標為單目標進行優化而是考慮應力水平、失效概率、高應力區范圍等安全性指標的多目標優化；拱壩優化已由早期的線性優化已轉變到非線性分析、混合線性優化及拓撲優化等。

結構拓撲優化的主要研究對象是連續體結構。連續體結構的拓撲優化設計出現了不同的求解體系。國內學者主要主要集中在局部應力約束下的強度拓撲優化設計，而國外研究人員多集中研究全局體積約束下的拓撲優化。

1992年，Xie和Steven提出了漸進結構優化算法，主要原理在優化區域以一定的刪除準則逐步刪除低應力單元，使得結構在滿足一定的約束條件下，應力狀態和目標函數達到最優。該法具有刪除單元后無法恢復的缺點；因此，Querin和Steven等對傳統的漸進結構優化算法進行了改進，在刪除單元后應力較高的區域增加單元，使得應力集中現象降低；Kim和Kwak對空間優化問題進行了深入的探討；Huang和Xie提出基于雙方向漸進結構優化的研究方法，對周期結構和具有一種或多種材料的結構的拓撲優化問題進行了拓展。

1 漸進結構優化方法基本思想和優化準則

1.1 漸進結構優化(ESO)方法

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

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

1.2 基于單元生死功能的優化準則

對于靜力設計問題，ESO方法通常采用基于應力準則的優化方法。對于各向同性材料，VonMises應力是目前最為常用的應用準則。

對于平面應力問題VonMises應力σvm定義為：

(1)

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

(2)

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

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

(3)

在改進的刪除率下，又有一些滿足刪除條件的單元，執行新一次單元刪除循環，直到達到沒有單元被刪除，停止迭代達到新的穩態，從而獲得一個最佳結構。理想的最終結構變成一個滿應力的設計即結構的每一點的材料應力達到它的滿強度。

1.3 有限元分析基本流程

(1)選擇準備設計區域，用有限元網格離散全部區域，將準備優化的區域定義不同的單元類型；

(2)對結構施加荷載和邊界約束，選擇靜力分析的方法進行靜力分析；

(4)將以上步驟重復進行，直至(2)式無法滿足且收斂到設計要求為止。此時，對應于RRi的結構已達到穩定狀態。為使優化繼續進行，引進進化率ER，從而下一穩定狀態的刪除率修改為式(3)；

(5)重復步驟2～4步直至結構滿足要求為止。

1.4 拓撲優化軟件實現

本文基于ESO法，利用ANSYS的單元生死對厚拱平面進行拓撲優化設計，采用APDL語言二次開發，編制循環程序設定刪除率，將設計區域內應力相對低的單元或低應變能密度的單元“殺死”。逐步“殺死”低效單元，最終得到滿足應力要求及最有目標的結構模型。

2 厚拱拱壩結構的拓撲優化計算

2.1 拱壩設計資料介紹

本文在對某二維厚拱拱壩進行拓撲優化應用的基礎上，其新意在于二次拱結構形式的分析。基本設計資料為：拱壩形式為單圓心拱壩，如圖1所示，圓心角θ=80°，外徑R=450 m，內徑r=350 m，拱厚T=100 m。厚拱判斷條件：D=0.5(R+r)，T/D=0.25，滿足厚拱(T/D=0.2～0.3)要求。壩體作為均質材料考慮，彈性模量Ec= 30.0 G Pa，泊松比取μc= 0.3，拱壩承受均布水壓力為350 kN/m2。

圖1 厚拱拱壩結構初始尺寸

2.2 拱壩結構拓撲優化計算模型

本例中厚拱結構采用ANSYS單元庫中8節點平面單元plane 82進行平面應力分析，將拱壩結構離散，共劃分10509個節點和5122個單元，其計算模型如圖2所示。

圖2 厚拱拱壩結構有限元計算模型

2.3 拱壩結構拓撲優化結果分析

基于Von Mises屈服準則確定的Von Mises Stress(等效應力)可以快速的確定模型中的最危險區域。

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

圖3 拱壩模型初始化單元等效應力云圖

圖4 拱壩模型初始化單元第一主應力云圖

圖5 拱壩模型初始化節點第一主應力云圖

圖6 迭代40次單元等效應力云圖

圖7 迭代40次單元第一主應力云圖

圖8 迭代40次節點第一主應力云圖

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

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

圖9 拱壩結構體積隨迭代次數變化曲線

圖10 單元等效應力最大值隨迭代次數變化曲線

圖11 節點第一主應最大值隨迭代次數變化曲線

2.4 拱壩結構二次拱擬合分析

為了驗證拱壩結構拓撲優化后，其內部可能存在的二次拱，為此對優化后的結構進行曲線擬合，得出二次拱結構的曲線。現利用Origin對拱壩拓撲優化后結構上的數據進行曲線擬合，其拱壩上下游曲線方程分別為：

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

(4)

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

(5)

上下游拓撲優化曲線與擬合曲線如圖12所示，圖中縱、橫坐標表示距離，單位為m。

圖12 拓撲優化后二次拱曲線

通過擬合的兩條拋物線建立拓撲優化后的二次拱模型，并對其加載及施加位移約束，分析其節點第一主應力及第三主應力，如圖13～14所示。

圖13 擬合結構第一主應力云圖

圖14 擬合結構第三主應力云圖

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

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

上述拓撲優化結果分析得出，單心圓拱壩的二次拱的結構形式呈拋物線，本節在拓撲優化的基礎上，分析厚拱拱壩結構的合理受力形式，采用與單心圓拱壩相同的工程背景，重點分析三心圓形式，對類似拱壩結構的設計具有借鑒和指導作用。

3.1 三心圓拱壩結構受力分析及結構優化

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

圖15 三心圓拱壩結構形式

圖16 三心圓拱應力分析及優化結構的應力

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

4 結 語

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

[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] 周春平，常錦昕. 基于單元生死功能的轉向架構架拓撲優化設計[J]. 計算機仿真, 2010, 27(5):267-271.

[6] 穆春燕. 拓撲優化理論及其在拱壩優化設計中的應用[D]. 南京：河海大學，2006.

[7] 郭中澤. 基于單元特性改變的漸進結構拓撲優化方法及應用研究[D].北京：中國工程物理研究所，2006.

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

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

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

[11] 董桂西,王藏柱.結構優化設計的現狀及展望[J].電力情報, 2000，(1):5-7.