基于Kriging模型的重力壩結(jié)構(gòu)優(yōu)化

王永強(qiáng)，馬孝義，李青峰，柳 燁

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

1 研究背景

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

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

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

2 Kriging代理模型

代理模型作為一種綜合建模技術(shù)，包含試驗(yàn)設(shè)計(jì)和近似方法。用代理模型取代復(fù)雜的數(shù)值模擬分析或昂貴的物理試驗(yàn)可顯著促進(jìn)多學(xué)科、多目標(biāo)設(shè)計(jì)優(yōu)化及概念設(shè)計(jì)等領(lǐng)域的發(fā)展。代理模型方法主要包括多項(xiàng)式響應(yīng)面模型、人工神經(jīng)網(wǎng)絡(luò)、支持向量機(jī)、徑向基函數(shù)模型和Kriging模型等。

2.1 Kriging模型原理

Kriging代理模型是基于統(tǒng)計(jì)理論的插值方法，于1951年由南非地質(zhì)學(xué)家Krige提出，最初用于礦產(chǎn)資源分布預(yù)測(cè)。自20世紀(jì)80年代末，其作為一種新型響應(yīng)面技術(shù)在結(jié)構(gòu)優(yōu)化領(lǐng)域得到廣泛應(yīng)用[12]，它可以對(duì)空間分布的樣本數(shù)據(jù)求出線性無偏內(nèi)插估計(jì)。該模型具有能夠?qū)φ`差做出逐點(diǎn)理論估計(jì)的優(yōu)點(diǎn)，同時(shí)存在相關(guān)函數(shù)需要根據(jù)經(jīng)驗(yàn)確定的缺點(diǎn)。Kriging模型的一般數(shù)學(xué)表達(dá)式為：

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

(1)

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

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為參數(shù)的相關(guān)函數(shù)。常用下式所示的高斯函數(shù)作為相關(guān)函數(shù)：

(5)

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

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

(6)

用各樣本點(diǎn)的響應(yīng)值線性加權(quán)疊加插值計(jì)算待測(cè)點(diǎn)xi的響應(yīng)值：

(7)

(8)

2.2 Kriging模型精度檢驗(yàn)

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

(9)

(10)

3 Kriging代理模型構(gòu)建

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

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

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

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

構(gòu)建壩體安全系數(shù)、橫斷面最大寬度及橫斷面面積代理模型一般需要4個(gè)步驟：①用某種試驗(yàn)設(shè)計(jì)方法產(chǎn)生設(shè)計(jì)變量的樣本點(diǎn)作為輸入數(shù)據(jù)；②用有限元仿真軟件計(jì)算獲得一組輸出數(shù)據(jù)；③用某種擬合方法近似擬合樣本點(diǎn)的輸入輸出關(guān)系，構(gòu)造代理模型；④對(duì)代理模型的擬合精度進(jìn)行檢驗(yàn)和評(píng)價(jià)，并利用其對(duì)新設(shè)計(jì)點(diǎn)的輸出進(jìn)行預(yù)測(cè)。

3.1 樣本點(diǎn)選取

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

3.2 壩體安全系數(shù)、橫斷面最大寬度、橫斷面面積Kriging代理模型構(gòu)建

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

表2 代理模型精度評(píng)估Tab.2 Agent model accuracy assessment

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

4 Kriging模型與遺傳算法結(jié)合的壩體結(jié)構(gòu)優(yōu)化

Kriging模型代替ANSYS有限元仿真分析并結(jié)合遺傳算法對(duì)壩體斷面進(jìn)行優(yōu)化，首先要建立優(yōu)化的數(shù)學(xué)模型。

4.1 設(shè)計(jì)變量選取

將設(shè)計(jì)變量應(yīng)用向量形式表示為Z=[m,n,p,q]，其中：m為上游邊坡系數(shù)，n為下游邊坡系數(shù)，p為上游折點(diǎn)高度與壩高的比值，q為下游折點(diǎn)高度與壩高的比值。

4.2 確定約束函數(shù)

c1=[0,0.60,0.32,0.86]

c2=[0.21,0.80,0.66,0.93]

K≥3；Di≥25 m

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

4.3 建立目標(biāo)函數(shù)

構(gòu)建目標(biāo)函數(shù)時(shí)引入罰函數(shù)[14]，即對(duì)于違反任一約束的情況，根據(jù)約束條件對(duì)目標(biāo)函數(shù)的影響程度，將一個(gè)懲罰項(xiàng)加入其中，引入一個(gè)新目標(biāo)函數(shù)，即：

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

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

4.4 壩體結(jié)構(gòu)優(yōu)化

應(yīng)用MATLAB遺傳算法工具箱編寫基于Kriging模型的優(yōu)化程序，經(jīng)過100次迭代獲得此次優(yōu)化的最佳目標(biāo)值為396.92 m2，此時(shí)各變量值及優(yōu)化初始值見表3。

表3 優(yōu)化前后各變量對(duì)比Tab.3 Comparison of variables before and after optimization

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

5 壩體材料分區(qū)對(duì)壩體應(yīng)力分布的影響

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

5.1 壩體分區(qū)

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

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

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

方案區(qū) 域A1A2A3一123二124

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

5.2 計(jì)算結(jié)果對(duì)比分析

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

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

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

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

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

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

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

表6 壩體分四區(qū)各方案應(yīng)力分析結(jié)果Tab.6 The stress analysis results of the four sections of the dam body

注：括號(hào)內(nèi)正數(shù)表示與壩體均勻未分區(qū)相比應(yīng)力增加量，負(fù)數(shù)代表應(yīng)力減少量。

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

6 結(jié) 語

本文基于Kriging代理模型，結(jié)合遺傳算法和功能梯度材料思想，對(duì)壩體幾何參數(shù)及壩體關(guān)鍵部位材料進(jìn)行優(yōu)化，得到如下結(jié)論。

(1)Kriging模型具有較高精度，基于該模型的優(yōu)化算法克服了優(yōu)化算法與有限元仿真分析直接結(jié)合帶來的計(jì)算量大，耗費(fèi)時(shí)間長(zhǎng)的困難，與后者相比前者大幅度縮短優(yōu)化運(yùn)算時(shí)間。

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

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

□