自適應網格方法在Stokes問題形狀最優控制中的應用?

段獻葆,李飛飛,秦新強

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

1 引言

大量的工程實際問題表明，初始設計階段選擇一個好的拓撲或形狀可以使得材料等有很大的節省.從Pironneau學者[1,2]開始，在過去的幾十年中，許多數學家和工程師致力于這一領域的研究，并且取得了非常好的成效[3-6].由于流體流動及其數學描述的復雜性，直到最近，特別是隨著計算技術的發展，流體中的形狀最優控制問題才得到很好的發展.但是，對流體力學中的拓撲優化問題進行系統的研究則只有十幾年的時間[7].

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

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

2 最優形狀設計問題

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

其中是需要優化的設計變量，也是用來控制介質在每個點處的滲透性.?=0表示固體材料，流體無法通過，而?=1表示流體流動的區域.θ是一個可以用來控制α(?)的正實數.在現實問題中，對于固體材料比如說墻壁，需要但從數值計算的角度只能選擇一個盡可能大的數，在本文中最小值為狀態方程為Stokes問題的最優控制問題的一般形式為：求u和?，使得目標泛函

在滿足

為求解最優控制問題(2)–(8)，需要目標泛函(2)關于區域的靈敏度分析結果.這可以用共軛方法得到，即，如果狀態方程是Stokes問題(3)–(6)，則目標泛函(2)關于?的形狀導數為[6]

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

在本文中，使用

作為最優化準則方法中的下降方向.

3 最優化準則及自適應網格方法

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

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

本文所提算法的目標是在降低計算量的前提使得最優控制問題所得到的形狀有更精確的邊界，所以在網格自適應的過程中需要一個指示函數.這可以從計算區域的材料分布信息得到，而這些信息在計算流體問題的時候就已經求出.我們將用如下函數作為自適應網格的指示函數

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

4 數值算例

4.1 算法

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

步驟1給定一個初始的設計變量并在整個計算區域上劃分粗網格.令迭代次數k=1；

步驟2進行迭代，直到滿足體積約束條件.第k次迭代由以下幾步構成：

步驟2.1所有的網格根據指示函數(15)標記為1或者0，對其中標記為1的網格進行加密；

步驟2.2根據上一步得到的設計變量首先求解Stokes問題(3)–(6)得到速度場；

步驟2.3由(9)求出目標泛函的形狀導數；

步驟2.4由(14)更新設計變量得到新的設計變量同時得到流體流動的新區域?k+1.相應的得到計算區域任一點的局部滲透性常數α(?(x))；

步驟2.5如果需要的話，輸出計算結果.

為了節省計算量，如果給出的初始設計變量過于簡單或與可能的最優拓撲或形狀差別較大，建議先對優化問題求解幾次，然后進行網格自適應的求解.

4.2 算例

本小節給出了兩個經典算例用以驗證所提算法的可靠性和有效性，這兩個算例也被很多研究者作為驗證算例[6,7,12-16].

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

算例1在第一個算例中，我們考慮改變沉浸在液體中物體的形狀，使得流體流動的阻力(能量耗散)達到最小.據我們所知，這是最經典，也是最早被用于流體力學形狀優化的算例[1].同樣的問題曾經出現在很多流體力學形狀或拓撲最優化的文獻中[13-16].

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

圖1:最小阻力問題的設計區域及邊界條件

最優化結果如圖2所示.兩個圖形分別為得到最優形狀時對應的自適應網格(左)和最優形狀(右).從結果可以看出，最優形狀是一個橄欖球形，與文獻[1]中的分析結果完全一致.

圖2:算例1的計算結果

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

最終的結果如圖4所示.左邊為最優化所得拓撲(形狀)相對應的自適應網格，右邊最優形狀.可以看出本算法所得結果要優于以前文獻中所得到的結果.同樣，本算例也說明了流體力學形狀最優控制的重要性，因為我們在設計初期是無法確定最終優化結果的拓撲或形狀的.

圖3:算例2的設計區域及邊界條件

圖4:算例2的計算結果

5 結論

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

參考文獻：

[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