基于Kriging模型的重力壩結構優化

王永強，馬孝義，李青峰，柳 燁

(西北農林科技大學水利與建筑工程學院，陜西 楊凌 712100)

1 研究背景

重力壩是由混凝土或漿砌石修筑的大體積擋水建筑物，其基本剖面是直角三角形，整體由若干壩段組成，主要依靠自重維持穩定。在壩體結構優化設計中常將高精度有限元分析與優化算法相結合進行計算[1,2]，這些工作往往計算量大、效率低，因此使用高效的優化設計方法對壩體結構進行優化很有必要。

為解決計算精度與效率之間的矛盾，有關學者提出了基于試驗設計和近似理論的代理模型方法。代理模型是一種利用已知點的響應信息來預測未知點響應值的近似擬合模型，計算量小且計算結果與高精度仿真模型計算結果相近，同時能過濾掉數值仿真中由于網格劃分、迭代收斂準則等因素造成的數值計算噪聲。楊麗[3]等基于Kriging模型和遺傳算法的齒輪修形減振優化驗證了Kriging模型預測嚙合剛度函數的精確性；鄭少平[4]等將Kriging代理模型技術應用于船舶板架強度和穩定性計算中指出，相比較徑向基函數和多項式響應面，Kriging代理模型對設計過程中起決定作用的控制應力的擬合效果較好，滿足工程精度要求；崔寶珍[5]等基于多項式響應面模型的立柱結構優化指出，代理模型比有限元仿真可顯著減少優化時間；王夢寒等[6]基于Kriging模型與遺傳算法結合的RHCM成型工藝參數優化驗證了該優化策略的可行性和合理性。上述研究表明Krig-ing模型法已成為結構分析和設計領域的研究熱點，該方法已在航天、船舶、汽車、橋梁等領域中被廣泛使用，在重力壩優化設計中的應用則罕見報道。在重力壩的優化設計中，壩踵應力控制是重難點[7]，薛小軍[8]指出，對最合理壩體體型進行材料分區可達到改善壩體應力分布的效果，但未對最合理體型壩體上、下游邊坡系數進一步優化；崔盼[9]指出，壩踵處使用彈性模量較小的混凝土，壩趾處使用彈性模量相對較大的混凝土，形成彈性模量梯度，可大幅度減小壩踵處最大拉應力，但未進一步指明該梯度的值。

本文借鑒前人研究成果，在蔣繪靜[1]的基礎上，構建壩體安全系數、橫斷面最大寬度和橫斷面面積的Kriging代理模型，以壩體上、下游邊坡系數及上、下游折點高度與壩高的比值為變量，以安全系數及橫斷面最大寬度為約束條件，以橫斷面面積最小為目標函數，利用遺傳算法[10]選擇出最優斷面；并在此基礎上，應用功能梯度材料思想[11]對壩體材料進行分區，使壩體應力分布更加合理。該研究可為相關單位進行混凝土重力壩結構優化設計提供參考。

2 Kriging代理模型

代理模型作為一種綜合建模技術，包含試驗設計和近似方法。用代理模型取代復雜的數值模擬分析或昂貴的物理試驗可顯著促進多學科、多目標設計優化及概念設計等領域的發展。代理模型方法主要包括多項式響應面模型、人工神經網絡、支持向量機、徑向基函數模型和Kriging模型等。

2.1 Kriging模型原理

Kriging代理模型是基于統計理論的插值方法，于1951年由南非地質學家Krige提出，最初用于礦產資源分布預測。自20世紀80年代末，其作為一種新型響應面技術在結構優化領域得到廣泛應用[12]，它可以對空間分布的樣本數據求出線性無偏內插估計。該模型具有能夠對誤差做出逐點理論估計的優點，同時存在相關函數需要根據經驗確定的缺點。Kriging模型的一般數學表達式為：

f(x)=g(x)+z(x)

(1)

式中：g(x)為回歸部分，一般為自變量x的多項式，可以是零階、一階或二階的形式；z(x)為隨機過程，其統計特性為：

E[z(x)]=0

(2)

Var[z(x)]=σ2

(3)

Cov[z(xi),z(x)]=σ2R(c,x,xi)

(4)

式中：R(c,x,xi)是以c為參數的相關函數。常用下式所示的高斯函數作為相關函數：

(5)

式中：di表征待測點xi與樣本點之間的距離，即di=|xj-xi|,j=1,…,n；ci表示相關函數在樣本點第i個方向的常數參量。由z(x)統計特性可得：

E[f(x)]=g(x)

(6)

用各樣本點的響應值線性加權疊加插值計算待測點xi的響應值：

(7)

(8)

2.2 Kriging模型精度檢驗

精度檢驗是判斷代理模型好壞的關鍵，若精度符合要求，則所建代理模型可代替高精度有限元仿真分析，否則需從試驗設計、設計變量等方面做出調整，進而重新構建代理模型。常用最大誤差[式(9)]和均方根誤差[式(10)]作為誤差檢驗標準：

(9)

(10)

3 Kriging代理模型構建

基于蔣繪靜[1]的研究成果，以某一典型混凝土重力壩為分析對象，選取其非溢流壩斷面建立二維模型，壩高29.65 m，壩頂寬5 m,如圖1所示，材料參數見表1。壩體承受的荷載有壩體自重、上下游靜水壓力及壩底揚壓力。

圖1 壩體橫斷面示意圖(單位：m)Fig.1 Schematic diagram of the cross section of the dam

材料彈性模量E/GPa泊松比容重/(kN·m-3)抗壓強度/MPa抗拉強度/MPa125.50.172013.401.542200.2525--3414[9]38[9]0.170.172424 ——

注：材料1、3、4分別用于壩體、材料2用于壩基。

構建壩體安全系數、橫斷面最大寬度及橫斷面面積代理模型一般需要4個步驟：①用某種試驗設計方法產生設計變量的樣本點作為輸入數據；②用有限元仿真軟件計算獲得一組輸出數據；③用某種擬合方法近似擬合樣本點的輸入輸出關系，構造代理模型；④對代理模型的擬合精度進行檢驗和評價，并利用其對新設計點的輸出進行預測。

3.1 樣本點選取

建立代理模型的前提是選取一系列能夠表示向量空間任何一部分的樣本點，要求取樣的數目可以是任意的且變量的維數對取樣方法影響較小，從而能夠較準確地反映有限元仿真模型的輸入、輸出關系。拉丁超立方取樣法作為一種特殊的多維分層抽樣方法能夠滿足這些要求。鄧乾旺[13]等人通過和隨機抽樣法對比，指出拉丁超立方所選取的樣本點具有分布均勻且很快能達到收斂的特點。通常情況下，樣本點數目不應少于(k+1)·(k+2)/2(k是變量個數)。本文選取m，n，p，q(分別為上、下游邊坡系數，上、下游折點高度與壩高的比值)作為試驗設計的設計變量，采用拉丁超立方取樣法選取200個樣本點用于模型構建，隨機產生獨立于樣本點的40個檢測點用于模型精度評估。

3.2 壩體安全系數、橫斷面最大寬度、橫斷面面積Kriging代理模型構建

利用有限元仿真軟件ANSYS APDL語言編寫參數化命令流，在MATLAB平臺上編制程序調用命令流文件，獲得樣本點的安全系數、橫斷面最大寬度、橫斷面面積響應值，分別將二階多項式、高斯函數作為回歸多項式、相關函數構建Kriging模型，并用均方根誤差和最大誤差對構建的代理模型精度進行評估(見表2)。

表2 代理模型精度評估Tab.2 Agent model accuracy assessment

如表2所示，安全系數代理模型的均方根誤差為0.24%，用均方根誤差除以仿真平均數可得安全系數誤差為0.04%，工程要求誤差比例不超過20%，最大誤差為4.31%，該誤差滿足精度要求[5]。在ANSYS中，橫斷面最大寬度與斷面面積不受網格數量及劃分方式的影響，且與設計變量有明確的數量關系，故二者代理模型的模擬值與真值完全相同。因此，文中構建的安全系數、橫斷面最大寬度、橫斷面面積代理模型可以近似替代ANSYS仿真模型對壩體結構進行優化設計。

4 Kriging模型與遺傳算法結合的壩體結構優化

Kriging模型代替ANSYS有限元仿真分析并結合遺傳算法對壩體斷面進行優化，首先要建立優化的數學模型。

4.1 設計變量選取

將設計變量應用向量形式表示為Z=[m,n,p,q]，其中：m為上游邊坡系數，n為下游邊坡系數，p為上游折點高度與壩高的比值，q為下游折點高度與壩高的比值。

4.2 確定約束函數

c1=[0,0.60,0.32,0.86]

c2=[0.21,0.80,0.66,0.93]

K≥3；Di≥25 m

式中：c1為設計變量上界；c2為設計變量下界；K為安全系數代理模型；Di為橫斷面最大寬度代理模型。

4.3 建立目標函數

構建目標函數時引入罰函數[14]，即對于違反任一約束的情況，根據約束條件對目標函數的影響程度，將一個懲罰項加入其中，引入一個新目標函數，即：

F=f(A)·e[f(a)+f(b)]

式中：f(A)為橫斷面面積代理模型；f(a)、f(b)為罰函數；F為新目標函數。

4.4 壩體結構優化

應用MATLAB遺傳算法工具箱編寫基于Kriging模型的優化程序，經過100次迭代獲得此次優化的最佳目標值為396.92 m2，此時各變量值及優化初始值見表3。

表3 優化前后各變量對比Tab.3 Comparison of variables before and after optimization

如表3所示，優化后4個設計變量分別減小5.00%、13.75%、15.87%、3.37%，斷面面積減小14.02%。此外，基于Kriging模型的優化算法對壩體橫斷面優化所用時間為78 s，若MATLAB直接調用ANSYS軟件進行優化所用時間為3 657 s，前者顯著縮短運算時間。對優化后的值利用ANSYS進行有限元仿真驗證，可得壩體最大拉應力為0.237 MPa，最大壓應力值為1.511 MPa，其值都小于許用值，抗滑穩定約束由抗剪斷公式驗算，得安全系數為k=6.53>3，滿足抗滑穩定要求[10]，由代理模型計算所得安全系數為k=6.55，滿足精度要求。

5 壩體材料分區對壩體應力分布的影響

基于上述所得驗證模型，根據《混凝土重力壩設計規范》(SL319-2005)[15]，應用功能梯度材料思想在壩體不同區域應用不同強度的混凝土材料，探究合理的分區及混凝土布局方式。

5.1 壩體分區

(1)分區方式1：壩體分兩區，如圖2(a)所示。各區材料布局見表4。

(2)分區方式2：壩體分四區，如圖2(b)所示。各區材料布局見表5，計算分析壩體應力。

圖2 壩體分區示意圖(單位：m)Fig.2 Schematic diagram of dam body division

方案區 域A1A2A3一123二124

注：數字1、2、3、4代表不同混凝土材料，其彈性模量分別為25.5、20、14、38 GPa。

5.2 計算結果對比分析

在基本荷載作用下，首先對壩體均勻未分區(壩體使用材料1，壩基使用材料2)和壩體分兩區進行有限元計算，壩踵-壩趾應力分布如圖3、圖4所示。

如圖3、圖4所示，增大壩體下部混凝土彈性模量，壩踵處主拉應力及壩趾處主壓應力減小，但y向應力增大，不利于壩體穩定；減小壩體下部混凝土彈性模量，壩趾處y向應力減小，但壩踵處主拉應力及壩趾處主壓應力增大。因此，需要探索一種細致的分區方式，使壩體應力分布更加合理。將壩體分四區各方案應力分析結果如表6所示。

表5 壩體分四區各區材料布局Tab.5 The dam body is divided into four districts and the material layout of each district

注：數字1、2、3、4代表不同混凝土材料，其彈性模量分別為25.5、20、14、38 GPa。

圖3 壩踵-壩趾y方向應力分布Fig.3 Stress distribution in the y direction from dam heel to toe

圖4 壩踵-壩趾第一主應力分布Fig.4 Distribution of the first principal stress from dam heel to toe

如表6所示，相比較均勻未分區，A1、A3區均采用材料3時，壩趾第一主應力小幅減小，壩踵及壩趾y方向應力分別減小40.57%和25.33%，但壩踵第一主應力增加12.83%，不利于壩踵處應力分布；A1、A3區均采用材料4時，壩踵及壩趾第一主應力分別減少10.70%和44.00%，但壩趾y向應力增加幅度較大，不利于壩體下游穩定，故方案三和方案四不合理。A1區采用材料3，A3區采用材料4時，壩踵及壩趾第一主應力均減小，但壩趾處y向應力增大；A1區采用材料4，A3區采用材料3時，雖然壩趾處第一主應力及y向應力、壩踵y向應力不同程度減小，但壩踵處第一主應力增加約27.27%，不利于壩體穩定，故方案五、方案六不可行。

表6 壩體分四區各方案應力分析結果Tab.6 The stress analysis results of the four sections of the dam body

注：括號內正數表示與壩體均勻未分區相比應力增加量，負數代表應力減少量。

方案八與方案六情況基本相同。A1區采用材料3，A3區采用材料1時，壩踵、壩趾第一主應力、壩踵、壩趾y向應力分別減小3.74%、40.00%、19.81%、8.33%；故方案七可行。A1區采用材料4，A3區采用材料1時，壩踵第一主應力及y向應力增大，壩趾及第一主應力及y向應力減?。籄1區采用材料1，A3采用材料4時，壩踵及壩趾第一主應力減小，但壩踵及壩趾y向應力增大，故方案九、方案十不可行。壩體分區高度未到達上游折點處，故上游折點處應力在各方案中應力變化不大，該處值在應力分析中不出現負值即可。方案五、方案十與方案七梯度方向相同，但二者梯度值分別為24、12.5 GPa，均大于方案七的11.5 GPa，導致壩體應力分布不合理。

6 結 語

本文基于Kriging代理模型，結合遺傳算法和功能梯度材料思想，對壩體幾何參數及壩體關鍵部位材料進行優化，得到如下結論。

(1)Kriging模型具有較高精度，基于該模型的優化算法克服了優化算法與有限元仿真分析直接結合帶來的計算量大，耗費時間長的困難，與后者相比前者大幅度縮短優化運算時間。

(2)上、下游折點高度與壩高的比值、上、下游邊坡系數及壩體橫斷面面積初始值分別為0.63、0.89、0.20、0.80、461.68 m2，優化后各值分別減小15.87%、3.37%、5.00%、13.75%，斷面面積減小14.02%。

(3)壩踵處采用彈性模量為14 GPa混凝土，壩趾處采用彈性模量為25.5 GPa混凝土且彈性模量梯度不大于11.5 GPa，有利于壩體應力分布。

□