滲流問(wèn)題的拓?fù)鋬?yōu)化

周向陽(yáng)，范良志

（武漢紡織大學(xué) 機(jī)械工程與自動(dòng)化學(xué)院，湖北 武漢 430073）

滲流問(wèn)題的拓?fù)鋬?yōu)化

周向陽(yáng)，范良志

（武漢紡織大學(xué) 機(jī)械工程與自動(dòng)化學(xué)院，湖北 武漢 430073）

分析拓?fù)鋬?yōu)化中的密度懲罰函數(shù)插值法SIMP（Solid isotropic material with penalization model），將結(jié)構(gòu)力學(xué)中的拓?fù)鋬?yōu)化方法應(yīng)用到滲流問(wèn)題的拓?fù)鋬?yōu)化設(shè)計(jì)中，用有限元的方法建立了簡(jiǎn)單的理想狀態(tài)下的滲流問(wèn)題的拓?fù)鋬?yōu)化數(shù)學(xué)模型，并采用基于剃度法的數(shù)值解法－優(yōu)化準(zhǔn)則法（OC系列算法），以設(shè)計(jì)具有最小能量損耗情況下的流體最佳流動(dòng)路徑的拓?fù)浞植肌Ｒ砸唤M二維滲流問(wèn)題為例，說(shuō)明了該模型的有效性，為滲流問(wèn)題的優(yōu)化設(shè)計(jì)提供了一種有效的新思路和方法。

拓?fù)鋬?yōu)化；滲流；最小能量損耗；最佳流動(dòng)路徑

結(jié)構(gòu)的拓?fù)鋬?yōu)化設(shè)計(jì)是結(jié)構(gòu)的尺寸優(yōu)化設(shè)計(jì)和形狀優(yōu)化設(shè)計(jì)以后，在結(jié)構(gòu)優(yōu)化領(lǐng)域出現(xiàn)的一種新型的富有挑戰(zhàn)性的研究方向[1]。從有關(guān)文獻(xiàn)看，連續(xù)體結(jié)構(gòu)的拓?fù)鋬?yōu)化設(shè)計(jì)技術(shù)，特別是SIMP方法[1,2]，已經(jīng)成功的用于機(jī)械設(shè)計(jì)[3,4]、MEMS系統(tǒng)[5,6]、材料設(shè)計(jì)[7]等方面。在航空航天、汽車制造等固體結(jié)構(gòu)領(lǐng)域，拓?fù)鋬?yōu)化技術(shù)不僅僅是學(xué)術(shù)上的應(yīng)用，而且成為一種實(shí)用的設(shè)計(jì)工具被廣泛采用。把拓?fù)鋬?yōu)化理論用于流體場(chǎng)流體流動(dòng)設(shè)計(jì)是拓?fù)鋬?yōu)化技術(shù)應(yīng)用的新的熱點(diǎn)研究方向之一。其方法是：在給定邊界條件的設(shè)計(jì)域Ω中，確定哪部分是流體哪部分是非流體，使得滿足規(guī)定流體部分體積比的某個(gè)目標(biāo)最小化。

由于滲流是流體中比較簡(jiǎn)單的情況，在實(shí)際工程中也比較常見(jiàn)，比如閘壩的滲流，農(nóng)田地下排水系統(tǒng)，油氣地下滲流，裂隙巖體滲流以及裂隙排水等，作者首先嘗試將用于固體結(jié)構(gòu)的剛度拓?fù)鋬?yōu)化技術(shù)用于滲流問(wèn)題的優(yōu)化設(shè)計(jì)。本文主要考慮流體在多孔介質(zhì)中的流動(dòng)，例如水在地壩中的流動(dòng)，在管道中或者圍繞固體的流動(dòng)，并且只考慮理想流體在穩(wěn)太、無(wú)旋（即流體質(zhì)點(diǎn)只是平動(dòng)）、不可壓縮（等質(zhì)量密度）和無(wú)粘（沒(méi)有粘性）狀態(tài)下的情況，流體與表面之間的摩擦也被忽略，流體也并不滲透到周圍物體或并不與物體的表面分開(kāi)。取介質(zhì)相對(duì)密度為設(shè)計(jì)變量，流體在整體設(shè)計(jì)域的體積比為約束，以系統(tǒng)的最小能量損耗為目標(biāo)函數(shù)。要說(shuō)明的是：目的在于探討拓?fù)鋬?yōu)化技術(shù)在流體方面的應(yīng)用的可行性和方法，由于本文選取理想滲流，所以要使模型符合實(shí)際的應(yīng)用還有待進(jìn)一步的研究。

1 SIMP方法回顧

SIMP模型主要通過(guò)引入懲罰因子，在材料的彈性模量和單元相對(duì)密度之間建立起一種顯示的非線性對(duì)應(yīng)關(guān)系。它的作用是當(dāng)設(shè)計(jì)變量的值在(0,1)之間時(shí)，對(duì)中間密度值進(jìn)行懲罰，使中間密度值逐漸向0-1兩端聚集，這樣可以使連續(xù)變量的拓?fù)鋬?yōu)化模型能較好地逼近原來(lái)0-1離散變量的優(yōu)化模型。

SIMP材料模型的數(shù)學(xué)表達(dá)形式：

E表示插值以后的彈性模量，0E為初始彈性模量表示單元j的設(shè)計(jì)變量即j單元的相對(duì)密度，為了避免剛度矩陣奇異，給xj一個(gè)大于0的下限值表示j單元初始剛度矩陣，表示第j單元優(yōu)化后的剛度矩陣。 為兩數(shù)學(xué)模型中對(duì)中間密度材料的懲罰因子。為有效壓縮中間密度材料，

以結(jié)構(gòu)的最小柔度設(shè)計(jì)問(wèn)題為例，其拓?fù)鋬?yōu)化模型可表示為：

X表示設(shè)計(jì)變量， C表示結(jié)構(gòu)的柔度，V表示優(yōu)化后的有限元單元體積列向量， V*表示優(yōu)化的目標(biāo)體積即體積約束，F(xiàn)表示力矢量，K表示結(jié)構(gòu)的剛度矩陣，U表示位移矢量。

這里只有一個(gè)約束條件，一般情況下采用優(yōu)化準(zhǔn)則法。它是由目標(biāo)函數(shù)和約束條件構(gòu)成的拉格朗日函數(shù)，在滿足Kuhn-Tucker條件下推導(dǎo)出相應(yīng)的迭代求解公式。在綜合考慮設(shè)計(jì)變量上下限的情況下，可得問(wèn)題（3）式的優(yōu)化準(zhǔn)則法求解公式如下：

式中 為阻尼系數(shù)，引入 的目的是為了確保數(shù)值計(jì)算的穩(wěn)定性和收斂性，有關(guān) 取值范圍的討論詳見(jiàn)文獻(xiàn)[8]。

2 滲流問(wèn)題拓?fù)鋬?yōu)化模型

類似地，可以建立滲流問(wèn)題的拓?fù)鋬?yōu)化模型。考慮在多孔介質(zhì)中的二維流體流動(dòng)，根據(jù)質(zhì)量守恒和達(dá)西定律，對(duì)于常數(shù)滲透系數(shù)流體的流動(dòng)微分方程為：

式中，Kxx和Kyy分別為多孔介質(zhì)在x和y方向的滲透系數(shù)，Φ是流體的水頭或者速度勢(shì)函數(shù)，是單位體積的體積流動(dòng)率。

式中Cx和Cy為表面S2的單位法線矢量的方向余弦，與圖1中所示的相同。式（8）表明，在邊界表面S1上已知邊界流體水頭或者速度勢(shì)ΦB，而式（9）則表明，在垂直于表面S2上的勢(shì)的剃度或者速度為已知。在不滲透的邊界上速度或剃度等于0。

一般來(lái)說(shuō)，基本單元可以有內(nèi)部源或匯，例如來(lái)自泵的源或匯，或來(lái)自河或江的表面邊界流動(dòng)率。為了包括這些影響，單元外載荷包括作用在整個(gè)單元上的均勻內(nèi)源Q和作用在表面上的均勻表面流動(dòng)率源*q，則力矩陣的項(xiàng)分別是：

由控制方程、邊界條件和初始條件，構(gòu)造滲流場(chǎng)的范函表達(dá)式，然后變分求解。在求解區(qū)域上離散單元，選取插值函數(shù)，經(jīng)過(guò)單元分析和組裝，即可得到滲流結(jié)構(gòu)的有限元表達(dá)式

類似的可以建立滲流問(wèn)題的拓?fù)鋬?yōu)化數(shù)學(xué)模型為：

圖1 流體流動(dòng)的邊界條件

式中[B]為聯(lián)系水力剃度與節(jié)點(diǎn)勢(shì)的矩陣，與固體結(jié)構(gòu)有限元中的單元應(yīng)變矩陣相似，其求法也與單元應(yīng)變矩陣一樣。[D]為材料性質(zhì)矩陣，對(duì)于各向同性滲透流體中有則

3 算法和二維算例

目前適用于拓?fù)鋬?yōu)化的優(yōu)化算法主要包括兩類：優(yōu)化準(zhǔn)則法（又稱OC系列算法）[9]、序列線性規(guī)劃或序列二次規(guī)劃算法，后者較為典型的是移動(dòng)漸進(jìn)算法（又稱MMA算法，the method of moving asymptotes）[10]。

優(yōu)化準(zhǔn)則法是一種間接的優(yōu)化方法，它不直接優(yōu)化目標(biāo)函數(shù)，而是基于Kuhn—Tucker條件，通過(guò)構(gòu)造Lagrange函數(shù)來(lái)形成設(shè)計(jì)變量的更新方案，一般適用于大量設(shè)計(jì)變量、單目標(biāo)、單約束條件下問(wèn)題的優(yōu)化。對(duì)于OC系列方法的詳細(xì)介紹見(jiàn)文獻(xiàn)[9]。

MMA方法是一種更高級(jí)的數(shù)學(xué)規(guī)劃算法，能夠適用于單約束情況和多目標(biāo)、多約束情況下問(wèn)題的求解，但其計(jì)算收斂性不夠理想。在實(shí)際工程應(yīng)用中，應(yīng)根據(jù)不同情況選用不同的優(yōu)化求解算法。關(guān)于MMA系列方法的詳細(xì)介紹參考Bruyneel等人的文獻(xiàn)[10]。

這里的問(wèn)題屬于單目標(biāo)單約束的優(yōu)化問(wèn)題，所以選擇優(yōu)化準(zhǔn)則法（OC算法），并在matlab上實(shí)現(xiàn)。

以二維平面問(wèn)題為例，設(shè)計(jì)出在不同的初始邊界條件時(shí)，使系統(tǒng)具有最小能量損耗的最佳流體流動(dòng)路徑。

為了消除拓?fù)鋬?yōu)化數(shù)值計(jì)算中的棋盤格和網(wǎng)格依賴現(xiàn)象，在算法中采用了懲罰項(xiàng)，取懲罰因子p=3.0，加入了阻力項(xiàng)，取阻力系數(shù)η= 0.5，采用了局部敏度過(guò)濾算法[11]。模型離散為四節(jié)點(diǎn)四邊形單元。滲透系數(shù)均為1。拓?fù)浞植紙D中，不可滲透材料部分用黑色表示，流體流動(dòng)部分用白色表示。

算例一：區(qū)域?yàn)?0×30單元，左右邊界分別有兩點(diǎn)水力，速度為1，水流進(jìn)為正，流出為負(fù)。結(jié)構(gòu)尺寸如圖2（a），體積約束比為0.3，優(yōu)化結(jié)果如圖2（b），迭代過(guò)程如圖2（c）。

算例二：區(qū)域?yàn)?0×30單元，左上角下右角兩點(diǎn)水力，速度為1。結(jié)構(gòu)尺寸如圖3（a），體積約束比為0.25，優(yōu)化結(jié)果如圖3（b），迭代過(guò)程如圖3（c）。

算例三：區(qū)域?yàn)?0×30單元，左邊界均勻分布水力，速度為1，右邊界中間節(jié)點(diǎn)勢(shì)為0。結(jié)構(gòu)尺寸如圖4（a）。體積約束比 f為0.5，優(yōu)化結(jié)果如圖4（b），迭代過(guò)程如圖4（c）。體積約束比 f為0.4，優(yōu)化結(jié)果如圖4（d），迭代過(guò)程如圖4（e）。

算例四：區(qū)域?yàn)?0×50單元，左右邊界分別有一點(diǎn)水力，速度為1。結(jié)構(gòu)尺寸如圖5（a），在中間偏下位置有一個(gè)直徑為15的圓形空洞。體積約束比 f為0.55，優(yōu)化結(jié)果如圖5（b），迭代過(guò)程如圖5（c）。體積約束比 f為0.3，優(yōu)化結(jié)果如圖5（d），迭代過(guò)程如圖5（e）。

圖2（a）算例一設(shè)計(jì)域

圖2（b）算例一優(yōu)化結(jié)果

圖2（c）算例一迭代過(guò)程

圖3（a）算例二的設(shè)計(jì)域

圖3（b）算例二優(yōu)化結(jié)果

圖3（c）算例二迭代過(guò)程

圖4（a）算例三的設(shè)計(jì)域

圖4（b）體積比0.5的優(yōu)化結(jié)果

圖4（c）體積比0.5的迭代過(guò)程

圖4（d）體積比0.4的優(yōu)化結(jié)果

圖4（e）體積比0.4的迭代過(guò)程

圖5（a）算例四的設(shè)計(jì)域

圖5（b）體積比0.55的優(yōu)化結(jié)果

圖5（c）體積比0.55的迭代過(guò)程

圖5（d）體積比0.3的優(yōu)化結(jié)果

圖5（e）體積比0.3的迭代過(guò)程

4 結(jié)論

從以上算例我們看出，拓?fù)鋬?yōu)化方法在滲流問(wèn)題的拓?fù)鋬?yōu)化設(shè)計(jì)中的應(yīng)用，沒(méi)有出現(xiàn)象在固體結(jié)構(gòu)設(shè)計(jì)中出現(xiàn)的數(shù)值不穩(wěn)定現(xiàn)象，比如多孔材料、棋盤格、網(wǎng)格依賴性、局部極值等。在算例三和四中出現(xiàn)了一些小島問(wèn)題只與約束體積比有關(guān)，當(dāng)系統(tǒng)采用的約束體積比太小時(shí)才會(huì)出現(xiàn)。算例表明了拓?fù)鋬?yōu)化方法在解滲流問(wèn)題時(shí)的可行性和有效性。

將拓?fù)鋬?yōu)化方法應(yīng)用于滲流問(wèn)題的拓?fù)鋬?yōu)化設(shè)計(jì)是一種新的思想和方法。該方法能夠設(shè)計(jì)出具有最小總體能量損耗情況下的流體最佳流動(dòng)路徑，為滲流結(jié)構(gòu)設(shè)計(jì)提供了一種新的方法和思路。

[1] Bends?e M, Sigmund O. Topology optimization: Theory, Methods, and Applications [M]. New York: Springer, 2003.

[2] Bendsoe M, Sigmund O. Material interpolations in topology optimization [J]. Archive of Applied Mechanics, 1999,69: 635～654.

[3] 石琴，姚成，馬恒永.集裝箱半掛車車架結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)[J].農(nóng)業(yè)機(jī)械學(xué)報(bào),2005,36(1):10～12，17.

[4] 范文杰，徐進(jìn)永，張子達(dá).基于雙向漸進(jìn)結(jié)構(gòu)優(yōu)化法的裝載機(jī)動(dòng)臂拓?fù)鋬?yōu)化[J].農(nóng)業(yè)機(jī)械學(xué)報(bào),2006,37(11):24～27.

[5] 左孔天, 陳立平, 王書亭, 等. 用拓?fù)鋬?yōu)化方法進(jìn)行微型柔性機(jī)構(gòu)的設(shè)計(jì)研究[J]. 中國(guó)機(jī)械工程, 2004, 15(21): 1886-1890.

[6] 羅震, 蒙永立, 郭文德, 等．分布式柔性機(jī)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì)的理論和算法[J]. 機(jī)械工程學(xué)報(bào)，2006，42(10): 27～36

[7] Sigmund O. A new class of extremal composites [J]. J Mech Phys Solid, 2000, 48(2)：397-428.

[8] 左孔天, 陳立平, 王書亭，等. 基于人工材料密度的新型拓?fù)鋬?yōu)化理論和算法研究[J]. 機(jī)械工程學(xué)報(bào), 2004, 40(12):23～29.

[9] Hassani B, Hinton E. Homogenization and Structural Topology Optimization Theory, Practice and Software [M]. London: Springer,1999.

[10] Bruyneel M, Duysinx P, Fleury C. A family of MMA approximations for structural optimization [J]. Structural and Multidiscipline Optimization, 2002, 24: 263～276.

[11] Sigmund O, Petersson J. Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards,mesh-dependencies and local minima [J]. Structural Optimization, 1998,16: 68～75.

Topology Optimization of Seepage Flow Problems

ZHOU Xiang-yang, FAN Liang-zhi
(College of Mechanical Engineering and Automation, Wuhan Textile University, Wuhan Hubei 430073, China)

The SIMP method used in topology optimization is analyzed. The topology optimization method is applied into the topology optimization design of seepage flow. The finite element method is used to model the simple flow problems under the ideal conditions. And the optimization problem is solved with a gradient-based math-programming algorithm (OC) that is driven by analytical sensitivities. The topology optimization model for seepage flow is established, which is used to minimize the energy dissipation in the system and find the optimal flow route. Several two-dimensional examples were tested. The results show the model is effective. A new approach for the design of seepage structure is presented.

Topology Optimization; Seepage Flow; Energy Dissipation Minimization; Optimal Flow Route

TH11

A

1009－5160(2011)06－0079－06

周向陽(yáng)（1977-），女，講師，博士，研究方向：數(shù)字化設(shè)計(jì)與制造.

國(guó)家自然科學(xué)基金（50805109）.