段獻(xiàn)葆， 黨 妍， 秦 玲

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

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

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

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

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

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

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

使得

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

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

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

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

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

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

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

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

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

2.1 靈敏度分析

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

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

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

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

引理1[8]考慮Sα關(guān)于邊界S的一個(gè)小擾動(dòng)

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

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

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

證明 對(duì)于目標(biāo)泛函

有

由引理1，可得

因此，

若?u是光滑的，則有

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

因此，

由泰勒展開，可得

而u｜Γ=0，所以有

設(shè)s表示切線分量，則有

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

因此，本工作使用

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

2.2 變分水平集方法

水平集方法的基本思想是將低維界面表示為一個(gè)函數(shù)在高維空間中的零水平集。對(duì)于開集Ω(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ū)域。

對(duì)方程(15)兩邊關(guān)于時(shí)間進(jìn)行微分，并應(yīng)用鏈?zhǔn)椒▌t，可得到Hamilton-Jacobi 型方程，

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

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

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

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

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

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

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

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

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

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

式中：φ(x)是水平集函數(shù)；h是到邊界距離的支撐（水平集函數(shù)的零水平集）?？梢钥闯觯志W(wǎng)格單元是否標(biāo)記為1 還是0 取決于粗網(wǎng)格單元到界面的距離。

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

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

自適應(yīng)水平集方法的具體步驟如下。

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

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

4 數(shù)值算例

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

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

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

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

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

5 結(jié)束語

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