劉小民，張彬，張卓颯

（西安交通大學(xué)能源與動(dòng)力工程學(xué)院，710049，西安）

流體形狀識(shí)別是反問題中的一個(gè)重要的研究方向，在生物醫(yī)學(xué)領(lǐng)域的心臟結(jié)構(gòu)反演、血管形狀設(shè)計(jì)及病變檢測(cè)等問題中廣泛應(yīng)用。帶有拓?fù)渥兓男螤钭R(shí)別問題是一類具有特殊目標(biāo)函數(shù)的拓?fù)鋬?yōu)化問題。拓?fù)鋬?yōu)化不依賴于給定結(jié)構(gòu)形狀或拓?fù)涞膬?yōu)化設(shè)計(jì)方法，最早用于解決基于固體力學(xué)原理的結(jié)構(gòu)優(yōu)化設(shè)計(jì)問題，主要包括均勻化法、變密度法、變厚度法、獨(dú)立連續(xù)映射法、漸進(jìn)結(jié)構(gòu)優(yōu)化法和水平集方法等［1］。水平集方法于1988年提出，Osher等在提出水平集方法的同時(shí)給出了求解水平集方程高精度的穩(wěn)定數(shù)值解法［2］。在基于固體力學(xué)的拓?fù)鋬?yōu)化研究領(lǐng)域，最早把水平集方法引入其中的是Sethian和Wiegmann，他們根據(jù)計(jì)算獲得的von Mises應(yīng)力進(jìn)行了水平集方程的求解［3］。將水平集方法與靈敏度分析相結(jié)合，可使水平集方法在處理優(yōu)化問題時(shí)更加簡便、準(zhǔn)確［4－6］。

隨著結(jié)構(gòu)拓?fù)鋬?yōu)化在固體力學(xué)領(lǐng)域的發(fā)展，流體拓?fù)鋬?yōu)化由Borrvall等應(yīng)用變密度方法得以實(shí)現(xiàn)［7］。但是，變密度方法會(huì)產(chǎn)生中間密度材料，這使得基于變密度方法的流體拓?fù)鋬?yōu)化很難得到較為精確的結(jié)果。最近，Duan等運(yùn)用一種無需重新初始化的變分水平集方法分別針對(duì)二維Stokes流動(dòng)及Navier－Stokes流動(dòng)進(jìn)行了拓?fù)鋬?yōu)化，實(shí)現(xiàn)了基于變分水平集方法的流體拓?fù)鋬?yōu)化［8］。Zhou等應(yīng)用傳統(tǒng)水平集方法對(duì)穩(wěn)態(tài)二維及三維Navier－Stokes流動(dòng)實(shí)施了拓?fù)鋬?yōu)化，并主張每一次優(yōu)化都要重新進(jìn)行網(wǎng)格劃分和法向速度擴(kuò)展［9］。Challis等針對(duì)二維和三維Stokes流動(dòng)進(jìn)行了拓?fù)鋬?yōu)化，但由優(yōu)化方法得到的邊界并不光滑［10］。Duan等應(yīng)用變分水平集方法對(duì)Stokes流動(dòng)形狀識(shí)別問題進(jìn)行了研究［11］。Chantalat等在Cartesian網(wǎng)格上運(yùn)用一種罰方法與水平集方法相結(jié)合的方法進(jìn)行了流體形狀識(shí)別和形狀優(yōu)化［12］，但該方法仍然存在優(yōu)化所得邊界不光滑以及不易解決復(fù)雜流動(dòng)等缺點(diǎn)，對(duì)含有梯度項(xiàng)的法向速度的直接求解不夠準(zhǔn)確，也無法直接加入無滑移邊界條件。

本文研究的目的是，通過提出一種新型的法向速度擴(kuò)展方法和無需樣條參數(shù)化的重新網(wǎng)格劃分方法，來改進(jìn)水平集優(yōu)化方法，通過求解Navier－Stokes和Stokes流動(dòng)形狀識(shí)別問題，來驗(yàn)證改進(jìn)水平集優(yōu)化方法的有效性。

1 數(shù)學(xué)模型

1．1 問題描述

設(shè)Ω是一個(gè)具有Lipschitz連續(xù)邊界Γ＝?Ω的有界開區(qū)域，代表流體域，則二維不可壓縮的Navier－Stokes方程如下

式中：Γ＝ΓN∪ΓD∪ΓS，并且ΓN、ΓD和ΓS兩兩之間無重合部分；ν是動(dòng)力黏性系數(shù)；f是體積力；u0是邊界ΓD上已知的速度分布；σ＝ν▽u－pI。如果式（1）中的對(duì)流項(xiàng)（u·▽）u被略去，則二維不可壓縮Stokes方程如下

由此，形狀識(shí)別問題可描述為

式中：u是滿足狀態(tài)約束式（1）或式（2）的速度分布；ud是已知的速度分布；D是優(yōu)化問題的工作區(qū)域，對(duì)于所有可能的Ω均有Ω∈D。

1．2 靈敏度分析

定理1 設(shè)Ω是一個(gè)具有光滑邊界Γ＝?Ω的有界開區(qū)域，h∈W1，∞（R2，R2），其中R是實(shí)數(shù)集，Navier－Stokes流動(dòng)形狀識(shí)別問題的目標(biāo)函數(shù)的Eulerian導(dǎo)數(shù)為［9］

式中：u是式（1）的解；w是如下共軛方程的解。該共軛方程如下

定理2 設(shè)Ω是一個(gè)具有光滑邊界Γ＝?Ω的有界開區(qū)域，h∈W1，∞（R2，R2），其中R是實(shí)數(shù)集，Stokes流動(dòng)形狀識(shí)別問題中目標(biāo)函數(shù)的Eulerian導(dǎo)數(shù)為［13］

式中：u是式（2）的解；w是如下共軛方程的解。該共軛方程如下

綜上，水平集演化方程為下面的Hamilton－Jacobi方程，即

式（8）中的水平集法向速度Vn應(yīng)該由以上的靈敏度分析結(jié)果確定。由定理1可以提取Navier－Stokes流動(dòng)形狀識(shí)別問題的水平集法向速度，即

由定理2可以確定Stokes流動(dòng)形狀識(shí)別問題的水平集函數(shù)演化的法向速度，即

式中：“∶”表示雙點(diǎn)乘運(yùn)算。

2 改進(jìn)的水平集方法

2．1 水平集函數(shù)的演化

水平集方法有許多優(yōu)點(diǎn)［14］：容易描述邊界法向及曲線等幾何特征，可以很好地控制演化過程中的拓?fù)渥兓?，能方便地?yīng)用偏微分方程描述，容易進(jìn)行數(shù)值求解。

水平集函數(shù)φ（x，t）在演化過程中滿足的方程如下

在流動(dòng)形狀識(shí)別問題中只有界面處的Vn有物理意義，為了在整個(gè)工作區(qū)域上求解式（11），應(yīng)進(jìn)行法向速度擴(kuò)展。

2．2 法向速度擴(kuò)展

文獻(xiàn)［15］利用快速行進(jìn)法進(jìn)行了法向速度擴(kuò)展，此方法在近界面處可以得到較準(zhǔn)確的擴(kuò)展結(jié)果，所求解的偏微分方程如下

雖然由文獻(xiàn)［15］可以得到靠近界面處節(jié)點(diǎn)上的法向速度，但是實(shí)施起來比較困難，特別是當(dāng)由本層節(jié)點(diǎn)尋找下一層節(jié)點(diǎn)位置時(shí)，實(shí)施過程很繁瑣。本文將此方法用于流體第一層節(jié)點(diǎn)到固體第一層節(jié)點(diǎn)上的法向速度的擴(kuò)展，對(duì)于整個(gè)水平集函數(shù)求解區(qū)域的擴(kuò)展，可以通過以下方程來實(shí)現(xiàn)［16］

式（13）為全場(chǎng)擴(kuò)展方法，基本思想是：將界面處的q值沿界面法向向兩側(cè)傳播。由式（13）的特點(diǎn)可知，遠(yuǎn)離界面節(jié)點(diǎn)上的q都是以靠近界面處節(jié)點(diǎn)上的q為參考的，q＝Vn即可表示水平集法向速度的擴(kuò)展。

將式（12）、式（13）結(jié)合起來，先在流體域直接求出邊界上較準(zhǔn)確的水平集法向速度，然后運(yùn)用快速行進(jìn)法將法向速度由流體域向未定義的固體域擴(kuò)展，最后再運(yùn)用求解偏微分方程的方法向整個(gè)水平集函數(shù)求解區(qū)域進(jìn)行擴(kuò)展?？梢?，本文方法在不做其他處理的情況下，直接、較準(zhǔn)確地求出邊界上的法向速度，且將界面附近的法向速度十分準(zhǔn)確地沿界面的法線方向擴(kuò)展到整個(gè)水平集函數(shù)求解區(qū)域。

3 數(shù)值方法

3．1 無需樣條參數(shù)化的重新網(wǎng)格劃分方法

本文采用的數(shù)值方法是無需樣條參數(shù)化的重新網(wǎng)格劃分方法，關(guān)鍵是在零水平集邊界上直接進(jìn)行三角網(wǎng)格劃分。不失一般性且便于表述，以二維問題為例描述如下。

當(dāng)已知網(wǎng)格節(jié)點(diǎn)rk－1的坐標(biāo)（xrk－1，yrk－1），求下一個(gè)三角網(wǎng)格節(jié)點(diǎn)rk的坐標(biāo)（xrk，yrk）時(shí)，先確定rk所在的四邊形單元。如果rk－1與rk的距離為Δdk－1，則rk滿足

式中：φrk為rk處的水平集函數(shù)值，rk所在單元的4個(gè)頂點(diǎn)a、d、f、e處的水平集函數(shù)值分別為φa、φd、φf、φe。設(shè)四邊形單元內(nèi)水平集函數(shù)值呈線性分布，則φrk值可由φa、φd、φf、φe進(jìn)行雙線性插值獲得，即

式中：ξ、η滿足

這樣，式（14）變?yōu)?/p>

式（18）可以通過Newton迭代法等進(jìn)行求解。在零水平集邊界上劃分出三角網(wǎng)格的邊界節(jié)點(diǎn)后，就可以運(yùn)用Delaunay等三角形化方法在流體域生成三角網(wǎng)格。

3．2 水平集方法的實(shí)施步驟

（1）給出初始的速度分布ud；

（2）給出初始區(qū)域Ω0，將Ω0初始化為符號(hào)距離函數(shù)；

（3）運(yùn)用無樣條參數(shù)化的重新網(wǎng)格劃分方法在流體域劃分出三角網(wǎng)格；

（4）利用非結(jié)構(gòu)化網(wǎng)格上的有限體積方法求解狀態(tài)方程（對(duì)于Navier－Stokes流動(dòng)形狀識(shí)別問題需要求解式（1），對(duì)于Stokes流動(dòng)形狀識(shí)別問題需要求解式（2））的速度u與壓力p，進(jìn)一步求解共軛方程（對(duì)于Navier－Stokes流動(dòng)形狀識(shí)別問題需要求解式（5），對(duì)于Stokes流動(dòng)形狀識(shí)別問題需要求解式（7））得到w和q；

（5）求取目標(biāo)函數(shù)J（u，Ω），判斷目標(biāo)函數(shù)值是否收斂，如果收斂結(jié)束程序，否則運(yùn)行步驟（6）；

（6）由式（9）（或式（10））求取水平集法向速度，并按照式（13）進(jìn)行法向速度擴(kuò)展；

（7）在結(jié)構(gòu)化網(wǎng)格上運(yùn)用差分方法求解水平集方程，并進(jìn)行水平集函數(shù)的演化；

（8）3次優(yōu)化迭代后進(jìn)行1次水平集函數(shù)重新初始化，然后回到步驟（3）進(jìn)一步循環(huán)，直到收斂為止。

3．3 有效性驗(yàn)證

在水平集拓?fù)鋬?yōu)化中，物理場(chǎng)的求解主要有重新網(wǎng)格劃分和浸沒邊界方法，本文采用重新網(wǎng)格劃分方法。為了驗(yàn)證重新網(wǎng)格劃分方法的有效性，分別采用重新網(wǎng)格劃分和浸沒邊界這2種方法對(duì)Stokes方程進(jìn)行求解。

根據(jù)雷諾準(zhǔn)則，任意2個(gè)符合雷諾準(zhǔn)則的流動(dòng)，只要滿足雷諾數(shù)相同條件，其流動(dòng)情況一致。本文研究的流動(dòng)問題均滿足雷諾準(zhǔn)則，所以算例中只給出了流動(dòng)雷諾數(shù)?？紤]到對(duì)稱性，特征尺度L為1／2進(jìn)口寬度，特征速度U 為進(jìn)口處平均速度。將研究問題的工作區(qū)域D描述為四邊形區(qū)域，流體域用Ω表示。在區(qū)域D的4條邊界ΓD上均給定Dirichlet邊界條件，速度為常量。如果以u(píng)＝（u1，u2）表示速度，那么u＊＝u／U＝（u＊1，u＊2）＝（u1／U，u2／U）在4條邊界上均有u＊1＝1、u＊0＝0，Ω與固體域的交界面則看作具有無滑移邊界條件的邊界ΓS。優(yōu)化問題的目標(biāo)函數(shù)為式（3），狀態(tài)約束為式（2），流場(chǎng)求解區(qū)域?yàn)榇须p圓繞流區(qū)域，如圖1所示，由定理2給出目標(biāo)函數(shù)的靈敏度分析結(jié)果。

圖1 Navier－Stokes流動(dòng)形狀識(shí)別問題的目標(biāo)流場(chǎng)

運(yùn)用重新網(wǎng)格劃分和浸沒邊界方法求解串列雙圓繞流問題時(shí)采用的網(wǎng)格分別如圖2a和圖3a所示。由圖2a可以看出，采用重新網(wǎng)格劃分方法時(shí)需先將雙圓形狀提取出來，然后在流體域進(jìn)行適體網(wǎng)格劃分。由圖3a可以看出，浸沒邊界方法是在整個(gè)設(shè)計(jì)區(qū)域上進(jìn)行網(wǎng)格劃分且不需要提取邊界。上述2種方法的網(wǎng)格劃分都是采用Delaunay方法生成的，在區(qū)域外邊界上的網(wǎng)格具有相同的邊界節(jié)點(diǎn)數(shù)。重新網(wǎng)格劃分和浸沒邊界方法求解雙圓繞流問題時(shí)獲得的速度矢量分布分別如圖2b和圖3b所示。由圖3b可以看出，雙圓附近出現(xiàn)了大范圍的低速區(qū)，這與串列雙圓實(shí)際流動(dòng)狀態(tài)不符。重新網(wǎng)格劃分和浸沒邊界方法求解串列雙圓繞流時(shí)獲得的流線分布分別如圖4和圖5所示。由圖4可以看出，重新網(wǎng)格劃分方法求解雙圓繞流時(shí)在界面附近獲得了合理的流動(dòng)分布，表明重新網(wǎng)格劃分方法能有效提高物理場(chǎng)控制方程的求解精度。由圖5可以看出，圓形界面處的流線沒有沿圓的邊界繞流，而是進(jìn)入了固體域，這與Stokes流動(dòng)的真實(shí)情況不符。

圖2 重新網(wǎng)格劃分方法的網(wǎng)格劃分和速度矢量分布

圖3 浸沒邊界方法的網(wǎng)格劃分和速度矢量分布

圖4 重新網(wǎng)格劃分方法獲得的串列圓形流線分布

圖5 浸沒邊界方法獲得的串列圓形流線分布

4 數(shù)值算例

4．1 Stokes流動(dòng)形狀識(shí)別問題

Stokes流動(dòng)是一種典型流動(dòng)，當(dāng)Re較小時(shí)，Navier－Stokes方程中對(duì)流項(xiàng)（u·▽）u的作用很小，甚至可以忽略不計(jì)，此時(shí)Navier－Stokes流動(dòng)可以當(dāng)作Stokes流動(dòng)進(jìn)行處理。對(duì)于Stokes流動(dòng)形狀識(shí)別問題，目標(biāo)形狀仍然取圖1所示的雙圓形狀，Re＝40。目標(biāo)流場(chǎng)的速度矢量分布ud＝（ud1，ud2），如圖6所示。

圖6 Stokes流動(dòng)形狀識(shí)別問題的目標(biāo)流場(chǎng)

一個(gè)四邊形作為串列雙圓形狀識(shí)別問題的初始形狀如圖7a所示，迭代10、20和30步時(shí)獲得的形狀及速度分布如圖7b～圖7d所示，Stokes流動(dòng)形狀識(shí)別過程得到的最終形狀如圖7e所示。由圖7可以看出：當(dāng)?shù)?0步時(shí)，初始形狀發(fā)生了拓?fù)渥兓?，表明本文方法可以有效處理帶有拓?fù)渥兓膬?yōu)化問題，Stokes流動(dòng)形狀識(shí)別過程中得到的最終形狀與目標(biāo)形狀很接近（見圖6）。目標(biāo)函數(shù)在形狀識(shí)別過程中的變化如圖8所示。由圖8可以看出，初始迭代時(shí)目標(biāo)函數(shù)變化較大，之后變化減緩，直至保持不變?yōu)橹?，由此獲得最優(yōu)結(jié)果。

4．2 Navier－Stokes流動(dòng)形狀識(shí)別問題

對(duì)于一個(gè)目標(biāo)為雙圓Navier－Stokes流動(dòng)形狀識(shí)別問題，其工作域、流體域、邊界條件等與Stokes流動(dòng)形狀識(shí)別問題相同，優(yōu)化問題的目標(biāo)函數(shù)為式（3），狀態(tài)約束為式（1），Re＝40。定理1給出了Navier－Stokes流動(dòng)形狀識(shí)別問題的靈敏度分析結(jié)果。一般情況下，目標(biāo)流場(chǎng)的速度ud＝（ud1，ud2）是通過實(shí)驗(yàn)測(cè)得的，本文主要是驗(yàn)證方法的有效性，所以目標(biāo)流場(chǎng)并未采用實(shí)驗(yàn)測(cè)量結(jié)果，而是通過數(shù)值計(jì)算結(jié)果獲得的目標(biāo)流場(chǎng)的速度分布。

串列雙圓Navier－Stokes流動(dòng)形狀識(shí)別問題的初始形狀如圖9a所示，迭代10、20、30和40步時(shí)獲得的形狀及速度分布如圖9b～圖9e所示，最終獲得的目標(biāo)形狀如圖9f示。由圖9可以看出：迭代40步時(shí)，初始形狀發(fā)生了拓?fù)渥兓砻鞅疚姆椒梢蕴幚韼в型負(fù)渥兓腘avier－Stokes流動(dòng)形狀識(shí)別問題；隨著迭代的繼續(xù)進(jìn)行，最終獲得的形狀與目標(biāo)形狀比較接近。目標(biāo)函數(shù)在形狀識(shí)別過程中的變化如圖10所示。由圖10可以看出：目標(biāo)函數(shù)的變化一直呈下降趨勢(shì)，最終達(dá)到最小值且保持不變。

圖7 Stokes流動(dòng)形狀識(shí)別過程

圖8 Stokes流動(dòng)形狀識(shí)別問題的目標(biāo)函數(shù)變化過程

圖9 Navier－Stokes流動(dòng)形狀識(shí)別過程

圖10 目標(biāo)函數(shù)在形狀識(shí)別過程中的變化

比較圖6和圖9可以看出，當(dāng)Re＝40時(shí)，Stokes和Navier－Stokes流動(dòng)形狀識(shí)別過程存在一定的差異。對(duì)于Navier－Stokes流動(dòng)形狀識(shí)別問題，迭代20步時(shí)尚未出現(xiàn)拓?fù)渥兓ㄒ妶D6c）；對(duì)于Stokes流動(dòng)形狀識(shí)別問題，迭代20步時(shí)初始形狀出現(xiàn)拓?fù)渥兓?，由區(qū)域中間一個(gè)孔洞變成了兩個(gè)孔洞（見圖9c）。由圖6e和圖9f可以看出，相比Navier－Stokes流動(dòng)形狀識(shí)別，Stokes識(shí)別的最終結(jié)果更加接近目標(biāo)形狀，即使在相同Re下，Stokes識(shí)別與Navier－Stokes仍有區(qū)別，這是Navier－Stokes方程中存在對(duì)流項(xiàng)（u·▽）u引起的。相對(duì)于Stokes流動(dòng)，Navier－Stokes流動(dòng)中對(duì)流項(xiàng)對(duì)流動(dòng)狀態(tài)的影響較大，Navier－Stokes流動(dòng)形狀識(shí)別問題中對(duì)應(yīng)的水平集法向速度的變化相對(duì)較大。當(dāng)時(shí)間步長一定時(shí)，為了滿足CFL條件，水平集方程求解需要對(duì)最大法向速度加以限制。最大法向速度相同時(shí)，Navier－Stokes流動(dòng)形狀識(shí)別問題的整體收斂速度較慢，因此在Navier－Stokes問題中，迭代40步后初始形狀才會(huì)發(fā)生拓?fù)渥兓赟tokes問題中，迭代20步后初始形狀即可達(dá)到拓?fù)渥兓?/p>

從算例中還可以看出，本文Stokes和Navier－Stokes流動(dòng)形狀識(shí)別問題最后所獲得的最優(yōu)形狀與實(shí)際形狀存在一些差異。這種差異可能由以下原因引起的：①對(duì)于復(fù)雜的流動(dòng)，梯度類優(yōu)化方法的局部收斂性使得優(yōu)化問題容易陷入局部極值，從以上算例的分析中可以看出，Stokes優(yōu)化問題相比Navier－Stokes優(yōu)化問題的收斂性更差；②本文目標(biāo)流場(chǎng)的形狀識(shí)別難度較高，在靠近流固界面處流動(dòng)速度接近0，而在圓柱內(nèi)部設(shè)定為0，這種速度分布無疑會(huì)給形狀的準(zhǔn)確識(shí)別帶來一定的難度。針對(duì)Navier－Stokes和Stokes流動(dòng)形狀識(shí)別，盡管在收斂速度及與目標(biāo)形狀的吻合程度方面本文還存在一定的偏差，但最終都能獲得與目標(biāo)流場(chǎng)和目標(biāo)形狀吻合較好的結(jié)果，這也表明本文發(fā)展的方法能夠有效處理流體形狀識(shí)別問題。

5 結(jié) 論

（1）本文提出了水平集法向速度擴(kuò)展方法，并將其與無需樣條參數(shù)化的重新網(wǎng)格劃分相結(jié)合，從而得到了一種較為精確的流體拓?fù)鋬?yōu)化方法。運(yùn)用此方法處理流體拓?fù)鋬?yōu)化問題時(shí)，既可直接施加于流體邊界，又可降低樣條參數(shù)化等繁瑣過程。

（2）對(duì)于 Navier－Stokes和Stokes流動(dòng)形狀識(shí)別問題，本文采用改進(jìn)的水平集方法可以有效處理優(yōu)化過程中的拓?fù)渥兓?，可以直接?yōu)化得到光滑的邊界，無需進(jìn)行邊界光滑化處理。

（3）對(duì)于 Navier－Stokes和Stokes流動(dòng)形狀識(shí)別問題，最終收斂的形狀與目標(biāo)形狀有一定的偏差，造成這種偏差的主要原因是梯度類優(yōu)化方法的全局性較差，但可以通過選擇全局性較好的方法來解決。在梯度類優(yōu)化領(lǐng)域，拓?fù)鋵?dǎo)數(shù)方法的全局性較好，用其進(jìn)行形狀識(shí)別可有效減小識(shí)別偏差，或者在本文優(yōu)化結(jié)果的基礎(chǔ)上，運(yùn)用非梯度類的全局優(yōu)化算法進(jìn)一步進(jìn)行優(yōu)化，可以減小最優(yōu)形狀與目標(biāo)形狀的偏離。

［1］ 夏天翔，姚衛(wèi)星．連續(xù)體結(jié)構(gòu)拓?fù)鋬?yōu)化方法評(píng)述 ［J］．航空工程進(jìn)展，2011，2（1）：1－11．XIA Tianxiang，YAO Weixing．A survey of topology optimization of continuum structure［J］．Advances in Aeronautical Science and Engineering，2011，2（1）：1－11．

［2］ OSHER S，F(xiàn)EDKIW R．Level set methods and dynamic implicit surfaces［M］．Berlin，Germany：Springer，2003：1－258．

［3］ SETHIAN J，WIEGMANN A．Structural boundary design via level set and immersed interface methods［J］．Journal of Computational Physics，2000，163（2）：489－528．

［4］ WANG M Y，WANG X，GUO D．A level set method for structural topology optimization ［J］．Computer Methods in Applied Mechanics Engineering，2003，192（1／2）：227－246．

［5］ ALLAIRE G，JOUVE J，TOADER A M．A level set method for shape optimization ［J］．Comptes Rendus de l’Académie des Sciences：Series I，2002，334（12）：1125－1130．

［6］ ALLAIRE G，JOUVE J，TOADER A M．Structural optimization using sensitivity analysis and a level set method［J］．Journal of Computational Physics，2004，194（1）：363－393．

［7］ BORRVALL T，PETERSSON J．Topology optimization of fluids in Stokes flow ［J］．International Journal for Numerical Methods in Fluids，2003，41（1）：77－107．

［8］ DUAN X，MA Y，ZHANG R．Optimal shape control of fluid flow using variational level set method ［J］．Physics Letters：A，2008，372（9）：1374－379．

［9］ ZHOU S，LI Q．A variational level set method for the topology optimization of steady－state Navier－Stokes flow ［J］．Journal of Computational Physics，2008，227（24）：10178－10195．

［10］CHALLIS V J，GUEST J K．Level set topology optimization of fluids in Stokes flow ［J］．International Journal for Numerical Methods in Engineering，2009，79（10）：1284－1308．

［11］DUAN X B，MA Y C，ZHANG R．Shape－topology optimization of Stokes flow via variational level set method［J］．Applied Mathematics and Computation，2008，202（1）：200－209．

［12］CHANTALAT F，BRUNEAU C H，GALUSINSKI C，et al．Level－set，penalization and Cartesian meshes：aparadigm for inverse problems and optimal design［J］．Journal of Computational Physics，2009，228（17）：6291－6315．

［13］張彬．水平集方法的改進(jìn)及其在流體形狀最優(yōu)控制中的應(yīng)用 ［D］．西安：西安交通大學(xué)，2011．

［14］ADALSTEINSSON D，SETHIAN J A．The fast construction of extension velocities in level set methods

［J］．Journal of Computational Physics，1999，148（1）：2－22．

［15］PENG D，MERRIMAN B，OSHER S，et al．A PDE－based fast local level set method［J］．Journal of Computational Physics，1999，155（2）：410－438．

［16］陶文銓．計(jì)算傳熱學(xué)的近代進(jìn)展 ［M］．北京：科學(xué)出版社，2000：1－427．

［本刊相關(guān)文獻(xiàn)鏈接］

王芳梅，范虹．利用改進(jìn)CV模型連續(xù)水平集算法的核磁共振乳腺圖像分割．2014，48（2）：38－43．［doi：10．7652／xjtuxb 201402007］

高燕華，劉玉歡，喻罡．多尺度非參數(shù)化水平集的超聲心動(dòng)圖分割．2013，47（2）：53－57．［doi：10．7652／xjtuxb201302009］

田輝，王琳琳，葉陽輝，等．應(yīng)用改進(jìn)粒子水平集法對(duì)液滴碰撞的數(shù)值研究．2012，46（1）：125－130．［doi：10．7652／xjtuxb 201201023］

田輝，王琳琳，葉陽輝，等．一種改進(jìn)的高效粒子水平集法及其應(yīng) 用．2011，45（11）：34－38．［doi：10．7652／xjtuxb201111 007］

李彥鵬，王煥然．低沖擊能量液滴與球面碰撞沉積特性的數(shù)值研究．2009，43（7）：21－24．［doi：10．7652／xjtuxb200907005］

段現(xiàn)報(bào)，馬逸塵，韓西安．變分水平集方法在Stokes問題形狀識(shí) 別 中 的 應(yīng) 用．2008，42（10）：1313－13．［doi：10．7652／xjtuxb200810025］