Stokes 問題形狀優(yōu)化中自適應水平集方法的應用

段獻葆， 黨 妍， 秦 玲

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

由于工業(yè)領域的迫切需要, 自Pironneau[1]開創(chuàng)性的工作以來，流體中的形狀優(yōu)化問題就吸引了很多工程師和數(shù)學家的關注。大量的實例表明，在工業(yè)設計的早期階段，選擇合適的形狀或拓撲結構可以大大地節(jié)省成本。在過去的幾十年里，流體動力學中的形狀優(yōu)化方法已取得了很大進展，許多方法已經(jīng)被引入并應用到工業(yè)設計問題中[2-7]。流體形狀優(yōu)化的細節(jié)和更多內(nèi)容參見文獻[8]。

Osher 等[9]提出的水平集方法(level set method)最初用來數(shù)值跟蹤或捕捉流體動力學中的動態(tài)界面和自由邊界。水平集方法不直接跟蹤界面，可以很容易地處理演化過程中的拓撲變化。水平集方法目前已經(jīng)被廣泛地應用到與形狀優(yōu)化相關的領域，如材料科學、數(shù)字圖像處理、幾何光學和生物學等[10-15]。在經(jīng)典的水平集方法中，水平集函數(shù)在演化的過程中需要保持為距離函數(shù)，但這一性質是很難保持的。通常的處理方法是在演化的過程中對水平集函數(shù)進行重新初始化，而這一過程往往是非常復雜費時的，且對邊界有很強的依賴性。為了解決這個問題，Chan 等[16]提出了變分水平集方法，該方法在水平集函數(shù)演化的過程中不再需要重新初始化，從而可以進一步提高計算效率。

在經(jīng)典的水平集方法中，水平集函數(shù)是定義在整個求解區(qū)域、采用均勻網(wǎng)格的有限差分格式實現(xiàn)，這使得求得較高精度數(shù)值解的時候會增加計算量。而很多時候人們關心的只是水平集函數(shù)的零水平集，為了克服這一局限性，提出了窄帶法[17]和局部水平集等[18]方法。此外，人們也提出了自適應方法來提高求解效率和精度[19]。

本工作的主要目的是研究流體流動引起的形狀優(yōu)化問題。為了提高優(yōu)化過程的計算效率，提出了一種幾何自適應網(wǎng)格方法。該方法具有靈活性、通用性等優(yōu)點，易于與現(xiàn)有的有限元方法軟件集成，非常容易實施。

1 控制方程和形狀優(yōu)化問題

設計算區(qū)域Ω是R2中一個具有Lipschitz 連續(xù)邊界Γ:=?Ω的有界連通區(qū)域。以偏微分方程為狀態(tài)方程的形狀優(yōu)化問題的一般形式如下：

使得

式中：γ(x)(0 ≤γ(x)≤1)是優(yōu)化變量，用來控制當?shù)亟橘|的滲透性；VD是設計領域的體積；β是流體可能占據(jù)的最大區(qū)域。

式中：φ(x)是水平集函數(shù)；γ=0表示固體材料，γ=1表示液體材料。方程(1)中目標泛函J是關于狀態(tài)變量u和設計變量γ的函數(shù)，對于阻力最小化問題，具有如下形式。

式中：ν=1/Re(Re是雷諾數(shù))是運動黏性系數(shù)。

式(2)～(3)構成了形狀優(yōu)化問題中的Dirichlet 型狀態(tài)約束。若狀態(tài)方程為Stokes 方程，則對應的Dirich-let 邊值問題為

式中：函數(shù)u=(u1,u2)；g、f和標量函數(shù)p分別表示流體的速度、邊界上規(guī)定的流速、外力和壓力。假設在理想多孔介質中流動的流體滿足Darcy 定律，即外力f和流體速度成正比。因此，f=P(x)u，P(x)是滲透介質在位置x的滲透率。基于不可壓縮性條件，利用散度定理容易看出函數(shù)g必須滿足相容性條件，

式中：n為單位外法向量。

與Borrvall 等[20]的研究相同，本工作選擇一個凸插值函數(shù)表示流體的滲透性參數(shù)P，

式中：0 ≤γ≤1,q是一個實的正參數(shù)，用于調(diào)整P(γ)的形狀。對于實際問題來說，不能滲透的固壁需要Pmax=∞, 但在實際計算中只能選擇一個有限值。本工作取Pmax=104，Pmin=0。

2 靈敏度分析和變分水平集方法

2.1 靈敏度分析

為了對水平集函數(shù)進行演化，需要得到優(yōu)化問題中目標函數(shù)(式(1))的靈敏度分析信息。設Ωα ?R2是區(qū)域Ω的一個擾動，Ωα=(Id+α)(Ω):={x+α(x),其中α是正的小參數(shù)，使得x ∈Ω}，δΩ=ΩαΩ，邊界為Γα={x+α(x)n(x),?x ∈Γ:=?Ω}。

設w=u(Ωα)-u(Ω)≡uα-u,是u在區(qū)域Ωα的零延拓。對于向量α, 開集Ω(α)是區(qū)域Ω小的擾動，(Id+α)是R2中一個微分同胚。

目標泛函J(Ω)關于區(qū)域Ω的形狀導數(shù)為W1,∞(R2,R2)中的Fr′echet 導數(shù)，

式中：DJ(Ω)是一個連續(xù)線性函數(shù)。

引理1[8]考慮Sα關于邊界S的一個小擾動

式中：α是關于x(s)∈S的函數(shù)(s曲線橫坐標)；λ是一個趨向于0 的正數(shù)。設Ωα=Ω(Sα)，則對于S附近的任意連續(xù)函數(shù)f有

定理1如果約束是Stokes 問題(式(6)～(8))，則目標函數(shù)J(u(γ),γ)關于區(qū)域Ω的形狀導數(shù)為

式中：n為局部單位外法向量。

證明 對于目標泛函

有

由引理1，可得

因此，

若?u是光滑的，則有

將式(6)乘以試驗函數(shù)w，并利用Green 公式，可得

因此，

由泰勒展開，可得

而u｜Γ=0，所以有

設s表示切線分量，則有

把式(11)和(12)代入式(10)，可得

因此，本工作使用

作為水平集函數(shù)的演化速度。

2.2 變分水平集方法

水平集方法的基本思想是將低維界面表示為一個函數(shù)在高維空間中的零水平集。對于開集Ω(x,t)，x ∈R2，t ∈R+，可以通過

定義R2×R+上的函數(shù)φ(x,t)來確定區(qū)域Ω(x,t)。邊界Γ(x,t) :=?Ω(x,t)為零水平集，

在區(qū)域Ω(x,t)外部有φ(x,t)＞0，即水平集函數(shù)φ(x,t)是滿足

的函數(shù)D為字體區(qū)域。

對方程(15)兩邊關于時間進行微分，并應用鏈式法則，可得到Hamilton-Jacobi 型方程，

式中：V= dφ(x,t)/dt是移動邊界的法向速度。式(16)即為水平集函數(shù)的運動方程。在水平集方法中，單位外法向量n和表面曲率κ分別為

在經(jīng)典的水平集方法中，水平集函數(shù)在演化過程中需要保持為符號距離函數(shù)的形式，而重新初始化的過程相當復雜。水平集函數(shù)φ(x,t)是否是一個距離函數(shù)當且僅當它滿足｜?φ(x,t)｜= 1。所以在文獻[16]中積分

被用來描述水平集函數(shù)φ(x,t)在Ω ∈R2中與距離函數(shù)的密切程度。假設水平集函數(shù)φ(x,t)的法向導數(shù)在邊界上為0，

則φ(x,t)在ψ方向的Gateaux 導數(shù)為

在本算法中，水平集函數(shù)將由基于形狀靈敏度的Hamilton-Jacobi 方程求解，

式中：λ是平衡(Δφ(x,t)-κ)影響的一個正參數(shù)。方程(18)的第二項能使水平集函數(shù)的零水平集趨向區(qū)域的邊界。與經(jīng)典的水平集方法相比，使用方程(18)求解水平集函數(shù)時，重新初始化的過程在優(yōu)化過程中可以忽略。

3 自適應網(wǎng)格方法

本工作通過一個數(shù)值阻尼項將參考域中的無滑移條件正則化，該數(shù)值阻尼項對于固體或流體速度較慢的材料來說很大的，而對于流體材料是接近于0。目標函數(shù)被重新表述為流體通過多孔介質時的功率耗散。解用0-1 拓撲表示，不包含中間情形，即為0 的區(qū)域是流體，為1的區(qū)域是固體，通過固體材料的流動是由Darcy 定律控制的。

基于上述分析，下面給出自適應網(wǎng)格法的主要步驟。

（1）對整個計算區(qū)域進行均勻剖分，在該網(wǎng)格上求解水平集函數(shù)。粗網(wǎng)格水平集函數(shù)的取值是對粗網(wǎng)格進行細化的關鍵。使用以下水平集函數(shù)作為細化指示子，

式中：φ(x)是水平集函數(shù)；h是到邊界距離的支撐（水平集函數(shù)的零水平集）。可以看出，粗網(wǎng)格單元是否標記為1 還是0 取決于粗網(wǎng)格單元到界面的距離。

（2）對細化指標標記為1 的粗網(wǎng)格進一步劃分為均勻的細網(wǎng)格。優(yōu)化過程中，粗網(wǎng)格內(nèi)均勻細網(wǎng)格的生成是自動完成的。

（3）對優(yōu)化問題進行求解。

自適應水平集方法的具體步驟如下。

（1）給出初始形狀Ω0和初始化水平集函數(shù)φ0(x)構建統(tǒng)一的粗網(wǎng)格的基礎上最初的形狀。設迭代次數(shù)k=1。

（2）進行迭代，直到達到體積約束。第k次迭代包括以下步驟：①所有粗網(wǎng)格按細化準則(式(19))進行分類，標記為1 的網(wǎng)格進一步劃分為均勻細網(wǎng)格，細化后重新排列單元和節(jié)點編號；②由所得數(shù)據(jù)可以得到γk，也即關于材料的分布，從而可以得到計算區(qū)域，利用細化網(wǎng)格求解Stokes 問題(式(6)～(8))，得到速度場；③由式(9)求出靈敏度分析結果；④通過求解粗網(wǎng)格上的水平集方程(18)更新水平集函數(shù)，由式(5)可以得到新的區(qū)域Ωk+1；⑤若需要，則進行可視化處理。

4 數(shù)值算例

本工作采用一個經(jīng)典的二維形狀優(yōu)化問題驗證所提算法的可行性和有效性，該算例可以在很多流體形狀優(yōu)化的文獻中找到[1-2,15,21-22]。假設固體區(qū)域是不能滲透的并且忽略外力項，固壁邊界采用無滑移邊界條件。在所有數(shù)據(jù)中，灰度圖片顯示最優(yōu)結果，黑色對應的設計值γ= 0(固體)，白色對應γ= 1(流體)。考慮的是浸入外部Stokes 流的物體形狀的優(yōu)化問題，目標是使物體引起的阻力最小。最小阻力問題的設計域如圖1 所示，其中狀態(tài)約束是Stokes方程，上、下及入流邊界上速度在x方向上是1，在y方向上是0。

圖1 最小阻力問題的設計域[1]Fig.1 Design domain and boundary conditions for the minimum drag example[1]

在本算例中，Stokes 問題(式(6)～(8))的求解采用有限元方法，在得到靈敏度分析結果后，通過插值映射到求解Hamilton-Jacobi 方程(式(18))的有限差分網(wǎng)格上。而在求解方程(18)的過程中，需要滿足CFL(Courant-Friedrichs-Lewy) 條件。本工作所采用的時間步長Δt需要滿足Δt≤Δx/｜V ｜max，其中｜V ｜max是由式(13)所得到結果在整個求解區(qū)域內(nèi)的最大值。在求解過程中，對于水平集函數(shù)的演化沒有重新初始化的過程，這也相應地減少了很多求解過程中的迭代，進一步提升了計算效率。本工作使用多種不同的初始形狀和體積約束進行優(yōu)化，發(fā)現(xiàn)結果對初始形狀和體積約束不敏感。本算例采用的初始形狀是一個簡單的圓形，流體體積為20%。最終的液體體積限制在80%和94%，結果在Re= 1 時得到。所有參數(shù)都與已有文獻保持一致。優(yōu)化結果如圖2 所示，其中最優(yōu)形狀對應的自適應網(wǎng)格見圖(a)和(b)，相應的水平集函數(shù)等值線圖見(c)和(d)。水平集函數(shù)的零水平集對應內(nèi)部實心物體的邊界，最優(yōu)形狀如圖2(e)和(f)所示。可以看出，該優(yōu)化過程得到了類似橄欖球的形狀，這符合文獻[1]的理論分析得到的形狀優(yōu)化結果，與文獻[15, 21-22]所得到的數(shù)值模擬結果也完全吻合。

圖2 最小阻力問題的最優(yōu)結果Fig.2 Optimal results of the minimum drag example

本算例的初始網(wǎng)格所用三角元1 332 個，節(jié)點數(shù)是714，求解水平集方程用的差分網(wǎng)格是20×20. 求解時間是516.24 s。為獲得相同的分辨率，本工作也通過均勻的細化網(wǎng)格（所用三角元13 126 個，節(jié)點6 108 個，差分網(wǎng)格是75×75）求解了該問題，最終結果一致，但運行時間1 652.96 s，約為自適應網(wǎng)格方法的3.2 倍，表明本算法的計算效率得到非常大的提升。

5 結束語

本工作給出了一種求解流體力學中形狀優(yōu)化問題的自適應水平集方法。該方法計算量小，數(shù)值實現(xiàn)簡單。使用這種自適應方法，計算量主要集中在界面(邊界)附近，因此，與全區(qū)域均勻網(wǎng)格相比，達到了相同的分辨率時的計算成本大大降低。雖然本工作的討論是以二維Stokes 問題為例，但是該方法可以推廣到求解復雜流體流動的形狀或拓撲優(yōu)化問題。