段獻(xiàn)葆,李飛飛,秦新強(qiáng)

(西安理工大學(xué)理學(xué)院，西安 710048)

1 引言

大量的工程實(shí)際問題表明，初始設(shè)計(jì)階段選擇一個(gè)好的拓?fù)浠蛐螤羁梢允沟貌牧系扔泻艽蟮墓?jié)省.從Pironneau學(xué)者[1,2]開始，在過去的幾十年中，許多數(shù)學(xué)家和工程師致力于這一領(lǐng)域的研究，并且取得了非常好的成效[3-6].由于流體流動(dòng)及其數(shù)學(xué)描述的復(fù)雜性，直到最近，特別是隨著計(jì)算技術(shù)的發(fā)展，流體中的形狀最優(yōu)控制問題才得到很好的發(fā)展.但是，對(duì)流體力學(xué)中的拓?fù)鋬?yōu)化問題進(jìn)行系統(tǒng)的研究則只有十幾年的時(shí)間[7].

迄今為止，已經(jīng)有很多用來求解形狀最優(yōu)控制問題的算法.其中最出名的，也用得非常成功的算法是由Bendsoe和Kikuchi[8]提出的均質(zhì)化方法(homogenization approach)，該算法把整個(gè)計(jì)算區(qū)域離散成足夠多的微小區(qū)域，然后根據(jù)所求解的物理問題對(duì)這些微小區(qū)域進(jìn)行分析，并按照一定的優(yōu)化算法刪除或增加某些微小區(qū)域，用最終保留下來的微小區(qū)域來描述最優(yōu)拓?fù)浠蛐螤?隨后，Bendsoe[9]對(duì)均質(zhì)化方法加以改進(jìn)，提出了SIMP(solid isotropic microstructure with penalization)方法.SIMP方法屬于變密度法的一種.均質(zhì)化方法是把計(jì)算區(qū)域中的離散點(diǎn)假設(shè)為相對(duì)密度在0到1之間的可變材料，而SIMP方法則是用連續(xù)函數(shù)來表達(dá)區(qū)域離散點(diǎn)處的相對(duì)密度與材料物理屬性之間的對(duì)應(yīng)關(guān)系，通過引入懲罰因子，使中間的密度函數(shù)值趨于0或1，從而取得與均質(zhì)化方法一樣的效果.該方法基于假設(shè)的各向同性材料，以每個(gè)單元的相對(duì)密度作為設(shè)計(jì)變量，不需要引入微結(jié)構(gòu)和均質(zhì)化過程.由于其程序?qū)崿F(xiàn)簡(jiǎn)單，計(jì)算高效穩(wěn)定，并且較好的消除了均質(zhì)化方法的棋盤格現(xiàn)象，SIMP方法在形狀最優(yōu)控制領(lǐng)域取得了非常廣泛的應(yīng)用.

不管是均質(zhì)化方法還是SIMP方法，一般都先要與其他優(yōu)化算法相結(jié)合，如最優(yōu)化準(zhǔn)則方法(OC)[3]，然后進(jìn)行迭代求解.當(dāng)所求解的問題比較簡(jiǎn)單，計(jì)算區(qū)域不太復(fù)雜，方程組規(guī)模較小的情況下，計(jì)算時(shí)間會(huì)比較短，是一個(gè)極為有效的方法.由于流體流動(dòng)的復(fù)雜性，采用這兩種方法時(shí)總計(jì)算量都非常大.其中很大一方面的原因是在實(shí)際優(yōu)化的過程中，一般都是對(duì)整個(gè)計(jì)算區(qū)域進(jìn)行均勻剖分.而在形狀最優(yōu)控制問題中，我們更關(guān)心的是計(jì)算區(qū)域的內(nèi)部邊界處形狀的改變.為此，我們提出了一種基于材料分布的自適應(yīng)網(wǎng)格方法，并將該方法用于求解流體力學(xué)形狀最優(yōu)控制問題.材料分布的信息在使用SIMP方法時(shí)已經(jīng)得到，不需要額外的計(jì)算.所提算法的計(jì)算主要集中在流體的內(nèi)邊界區(qū)域，從而在達(dá)到同樣分辨率的情況下可以有效的減少計(jì)算量.

2 最優(yōu)形狀設(shè)計(jì)問題

在流體力學(xué)最優(yōu)形狀設(shè)計(jì)問題中，一般考慮在控制方程約束下使得流體流動(dòng)的能量耗散達(dá)到最小.設(shè)D為包含所有可能形狀?的工作區(qū)域，是R2中具有Lipschitz連續(xù)邊界的開區(qū)域和p分別表示流體的速度和壓力，f表示體積力.不失一般性，我們假設(shè)流體是在理想介質(zhì)中流動(dòng)，并且符合Darcy定律[7]，即是在點(diǎn)x的局部滲透性常數(shù)的倒數(shù)，可以表示為如下凸函數(shù)

其中是需要優(yōu)化的設(shè)計(jì)變量，也是用來控制介質(zhì)在每個(gè)點(diǎn)處的滲透性.?=0表示固體材料，流體無法通過，而?=1表示流體流動(dòng)的區(qū)域.θ是一個(gè)可以用來控制α(?)的正實(shí)數(shù).在現(xiàn)實(shí)問題中，對(duì)于固體材料比如說墻壁，需要但從數(shù)值計(jì)算的角度只能選擇一個(gè)盡可能大的數(shù)，在本文中最小值為狀態(tài)方程為Stokes問題的最優(yōu)控制問題的一般形式為：求u和?，使得目標(biāo)泛函

在滿足

為求解最優(yōu)控制問題(2)–(8)，需要目標(biāo)泛函(2)關(guān)于區(qū)域的靈敏度分析結(jié)果.這可以用共軛方法得到，即，如果狀態(tài)方程是Stokes問題(3)–(6)，則目標(biāo)泛函(2)關(guān)于?的形狀導(dǎo)數(shù)為[6]

其中ε＞0是一個(gè)小的正實(shí)數(shù)，(P,q)是下面問題的解

在本文中，使用

作為最優(yōu)化準(zhǔn)則方法中的下降方向.

3 最優(yōu)化準(zhǔn)則及自適應(yīng)網(wǎng)格方法

最優(yōu)化準(zhǔn)則方法在上世紀(jì)50年代就已經(jīng)被用于工程結(jié)構(gòu)設(shè)計(jì)，后來數(shù)學(xué)家和工程師把Kuhn-Tucker條件作為最優(yōu)結(jié)構(gòu)滿足的準(zhǔn)則，從而使該方法有了數(shù)學(xué)理論的保障，通用性也得到提高.該方法的最大優(yōu)點(diǎn)是對(duì)設(shè)計(jì)變量修改較大，所以收斂速度快，迭代次數(shù)與結(jié)構(gòu)大小及復(fù)雜程度無關(guān)，并且原理簡(jiǎn)單、直觀、易為工程設(shè)計(jì)人員接受與掌握.最優(yōu)化準(zhǔn)則方法已經(jīng)被用于求解很大一類最優(yōu)形狀設(shè)計(jì)問題.本文所提的算法同樣可以與移動(dòng)漸進(jìn)方法(MMA)[10]或Level set方法[11]等優(yōu)化算法相耦合，用以解決更為復(fù)雜的問題，但為了簡(jiǎn)便且不失一般性，這里采用標(biāo)準(zhǔn)的最優(yōu)化準(zhǔn)則方法[3].優(yōu)化過程中的設(shè)計(jì)變量采用下面的更新準(zhǔn)則

其中為最小相對(duì)密度向量(為了避免奇性，沒有取為0)，ρ＞0是移動(dòng)的范圍，n=1/2是控制靈敏度影響的系數(shù).

本文所提算法的目標(biāo)是在降低計(jì)算量的前提使得最優(yōu)控制問題所得到的形狀有更精確的邊界，所以在網(wǎng)格自適應(yīng)的過程中需要一個(gè)指示函數(shù).這可以從計(jì)算區(qū)域的材料分布信息得到，而這些信息在計(jì)算流體問題的時(shí)候就已經(jīng)求出.我們將用如下函數(shù)作為自適應(yīng)網(wǎng)格的指示函數(shù)

其中在每個(gè)單元上的平均值，h為所有單元的最大邊長(zhǎng).由(15)可以看出，粗網(wǎng)格單元被標(biāo)記為1或0與該單元與邊界的距離有關(guān)，確切的說，只有在邊界附近的單元才被標(biāo)記為1，其他單元被標(biāo)記為0.

4 數(shù)值算例

4.1 算法

在上面理論分析的基礎(chǔ)上，我們可以得到如下算法：

步驟1給定一個(gè)初始的設(shè)計(jì)變量并在整個(gè)計(jì)算區(qū)域上劃分粗網(wǎng)格.令迭代次數(shù)k=1；

步驟2進(jìn)行迭代，直到滿足體積約束條件.第k次迭代由以下幾步構(gòu)成：

步驟2.1所有的網(wǎng)格根據(jù)指示函數(shù)(15)標(biāo)記為1或者0，對(duì)其中標(biāo)記為1的網(wǎng)格進(jìn)行加密；

步驟2.2根據(jù)上一步得到的設(shè)計(jì)變量首先求解Stokes問題(3)–(6)得到速度場(chǎng)；

步驟2.3由(9)求出目標(biāo)泛函的形狀導(dǎo)數(shù)；

步驟2.4由(14)更新設(shè)計(jì)變量得到新的設(shè)計(jì)變量同時(shí)得到流體流動(dòng)的新區(qū)域?k+1.相應(yīng)的得到計(jì)算區(qū)域任一點(diǎn)的局部滲透性常數(shù)α(?(x))；

步驟2.5如果需要的話，輸出計(jì)算結(jié)果.

為了節(jié)省計(jì)算量，如果給出的初始設(shè)計(jì)變量過于簡(jiǎn)單或與可能的最優(yōu)拓?fù)浠蛐螤畈顒e較大，建議先對(duì)優(yōu)化問題求解幾次，然后進(jìn)行網(wǎng)格自適應(yīng)的求解.

4.2 算例

本小節(jié)給出了兩個(gè)經(jīng)典算例用以驗(yàn)證所提算法的可靠性和有效性，這兩個(gè)算例也被很多研究者作為驗(yàn)證算例[6,7,12-16].

在兩個(gè)算例中，固體區(qū)域都是假設(shè)不可滲透的，沒有考慮液體本身的重力，沒有液體流過的邊界滿足無滑移條件.在得到的最優(yōu)化結(jié)果中，黑色區(qū)域表示?=0(固體，液體不可滲透)，白色區(qū)域表示?=1(液體).

算例1在第一個(gè)算例中，我們考慮改變沉浸在液體中物體的形狀，使得流體流動(dòng)的阻力(能量耗散)達(dá)到最小.據(jù)我們所知，這是最經(jīng)典，也是最早被用于流體力學(xué)形狀優(yōu)化的算例[1].同樣的問題曾經(jīng)出現(xiàn)在很多流體力學(xué)形狀或拓?fù)渥顑?yōu)化的文獻(xiàn)中[13-16].

問題的設(shè)計(jì)區(qū)域如圖1所示.考慮Dirchlet問題，區(qū)域外邊界上速度在x方向的分量為1，在y方向的分量為0，沉浸物體的邊界是無滑移條件.經(jīng)過我們的測(cè)試，最終結(jié)果對(duì)設(shè)計(jì)變量的初始值不敏感，所以我們假設(shè)在整個(gè)計(jì)算區(qū)域上的初始值都為0.4.得到最優(yōu)形狀時(shí)流體占整個(gè)計(jì)算區(qū)域的94%，雷諾數(shù)Re=1.這些假設(shè)與其他人工作中的設(shè)定是一致的.

圖1:最小阻力問題的設(shè)計(jì)區(qū)域及邊界條件

最優(yōu)化結(jié)果如圖2所示.兩個(gè)圖形分別為得到最優(yōu)形狀時(shí)對(duì)應(yīng)的自適應(yīng)網(wǎng)格(左)和最優(yōu)形狀(右).從結(jié)果可以看出，最優(yōu)形狀是一個(gè)橄欖球形，與文獻(xiàn)[1]中的分析結(jié)果完全一致.

圖2:算例1的計(jì)算結(jié)果

算例2在第二個(gè)算例中，我們考慮有兩個(gè)入口和兩個(gè)出口的情形，類似的問題同樣被很多學(xué)者研究過[7,12,13,15,16].計(jì)算區(qū)域和邊界條件如圖3所示，有對(duì)稱的兩個(gè)入口和出口.在入口處速度在x方向的分量的大小為拋物狀，出口處的壓力為0，其它邊界上面滿足無滑移條件.區(qū)域的縱橫比為1:1.優(yōu)化的目標(biāo)是使得流體在整個(gè)計(jì)算區(qū)域上面的能量耗散最小，同時(shí)滿足最終流體所占區(qū)域?yàn)檎麄€(gè)計(jì)算區(qū)域的40%.本算例的最終結(jié)果同樣對(duì)設(shè)計(jì)變量的初始值不敏感，所以我們和上一算例一樣假設(shè)在整個(gè)計(jì)算區(qū)域上的初始值都為0.4.

最終的結(jié)果如圖4所示.左邊為最優(yōu)化所得拓?fù)?形狀)相對(duì)應(yīng)的自適應(yīng)網(wǎng)格，右邊最優(yōu)形狀.可以看出本算法所得結(jié)果要優(yōu)于以前文獻(xiàn)中所得到的結(jié)果.同樣，本算例也說明了流體力學(xué)形狀最優(yōu)控制的重要性，因?yàn)槲覀冊(cè)谠O(shè)計(jì)初期是無法確定最終優(yōu)化結(jié)果的拓?fù)浠蛐螤畹?

圖3:算例2的設(shè)計(jì)區(qū)域及邊界條件

圖4:算例2的計(jì)算結(jié)果

5 結(jié)論

本文給出了一種求解流體力學(xué)形狀最優(yōu)控制問題的新算法，該算法非常容易實(shí)施.并且由于運(yùn)算主要集中在流體流動(dòng)的內(nèi)部邊界附近，從而在達(dá)到同樣精度的情況下可以很大的節(jié)省計(jì)算量.盡管最優(yōu)化算法采用的是最優(yōu)性準(zhǔn)則方法，并且只是以Stokes問題為例進(jìn)行討論，本算法完全可以采用其他更好的優(yōu)化算法用以求解更為復(fù)雜的流體力學(xué)最優(yōu)形狀控制問題.提供的數(shù)值算例說明了算法的高效性和穩(wěn)定性.

參考文獻(xiàn)：

[1]Pironneau O.On optimal profiles in Stokes flow[J].Journal of Fluid Mechanics,1973,59(1):117-128

[2]Pironneau O.On optimum design in fluid mechanics[J].Journal of Fluid Mechanics,1974,64(1):97-110

[3]Bends?e M P,Sigmund O.Topology Optimization,Theory,Methods,and Applications[M].New York:Springer-Verlag,2003

[4]Thévenin D,Janiga G.Optimization and Computational Fluid Dynamics[M].Berlin:Springer-Verlag,2008

[5]Plotnikov P,Sokolowski J.Compressible Navier-Stokes Equations:Theory and Shape Optimization[M].Birkhuser:Springer Basel,2012

[6]Mohammadi B,Pironneau O.Applied Shape Optimization for Fluids(2nd Edition)[M].Oxford:Oxford University Press,2010

[7]Borrvall T,Petersson J.Topology optimization of fluid in Stokes flow[J].International Journal for Numerical Methods in Fluid,2003,41(1):77-107

[8]Bends?e M P,Kikuchi N.Generating optimal topologies in structural design using a homogenization method[J].Computer Methods in Applied Mechanics and Engineering,1988,71(2):197-224

[9]Bends?e M P.Optimal shape design as a material distribution problem[J].Structural and Multidisciplinary Optimization,1989,1(4):193-202

[10]Svanberg K.The method of moving asymptotes—a new method for structural optimization[J].International Journal of Numerical Methods in Engineering,1987,24(2):359-373

[11]Osher S,Sethian J A.Level Set Methods and Dynamic Implicit Surface[M].New York:Springer-Verlag,2002

[12]Olesen L H,Okkels F,Bruus H.A high-level programming-language implementation of optimization applied to steady-state Navier-Stokes flow[J].International Journal of Numerical Methods in Engineering,2006,65(7):975-1001

[13]Challis V J,Guest J K.Level set topology optimization of fluids in Stokes flow[J].International Journal of Numerical Methods in Engineering,2009,79(10):1284-1308

[14]Duan X B,Ma Y C,Zhang R.Optimal shape control of fluid flow using variational level set method[J].Physics Letters A,2008,372(9):1374-1379

[15]Guest J K,Prevost J H.Optimization of creeping fluid flows using a Darcy-Stokes finite element[J].International Journal of Numerical Methods in Engineering,2006,66(3):461-484

[16]Pingen G,Evgrafov A,Maute K.Topology optimization of flow domains using the lattice Boltzmann method[J].Structural and Multidisciplinary Optimization,2007,34(6):507-524