水平集與靈敏度分析相結(jié)合的流體形狀識別與優(yōu)化方法

劉小民，張彬，張卓颯

（西安交通大學能源與動力工程學院，710049，西安）

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

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

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

1 數(shù)學模型

1．1 問題描述

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

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

由此，形狀識別問題可描述為

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

1．2 靈敏度分析

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

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

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

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

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

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

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

式中：“∶”表示雙點乘運算。

2 改進的水平集方法

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

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

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

在流動形狀識別問題中只有界面處的Vn有物理意義，為了在整個工作區(qū)域上求解式（11），應進行法向速度擴展。

2．2 法向速度擴展

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

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

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

將式（12）、式（13）結(jié)合起來，先在流體域直接求出邊界上較準確的水平集法向速度，然后運用快速行進法將法向速度由流體域向未定義的固體域擴展，最后再運用求解偏微分方程的方法向整個水平集函數(shù)求解區(qū)域進行擴展?？梢姡疚姆椒ㄔ诓蛔銎渌幚淼那闆r下，直接、較準確地求出邊界上的法向速度，且將界面附近的法向速度十分準確地沿界面的法線方向擴展到整個水平集函數(shù)求解區(qū)域。

3 數(shù)值方法

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

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

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

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

式中：ξ、η滿足

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

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

3．2 水平集方法的實施步驟

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

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

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

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

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

（6）由式（9）（或式（10））求取水平集法向速度，并按照式（13）進行法向速度擴展；

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

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

3．3 有效性驗證

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

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

圖1 Navier－Stokes流動形狀識別問題的目標流場

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

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

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

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

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

4 數(shù)值算例

4．1 Stokes流動形狀識別問題

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

圖6 Stokes流動形狀識別問題的目標流場

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

4．2 Navier－Stokes流動形狀識別問題

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

串列雙圓Navier－Stokes流動形狀識別問題的初始形狀如圖9a所示，迭代10、20、30和40步時獲得的形狀及速度分布如圖9b～圖9e所示，最終獲得的目標形狀如圖9f示。由圖9可以看出：迭代40步時，初始形狀發(fā)生了拓撲變化，表明本文方法可以處理帶有拓撲變化的Navier－Stokes流動形狀識別問題；隨著迭代的繼續(xù)進行，最終獲得的形狀與目標形狀比較接近。目標函數(shù)在形狀識別過程中的變化如圖10所示。由圖10可以看出：目標函數(shù)的變化一直呈下降趨勢，最終達到最小值且保持不變。

圖7 Stokes流動形狀識別過程

圖8 Stokes流動形狀識別問題的目標函數(shù)變化過程

圖9 Navier－Stokes流動形狀識別過程

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

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

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

5 結(jié) 論

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

（2）對于 Navier－Stokes和Stokes流動形狀識別問題，本文采用改進的水平集方法可以有效處理優(yōu)化過程中的拓撲變化，可以直接優(yōu)化得到光滑的邊界，無需進行邊界光滑化處理。

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

［1］ 夏天翔，姚衛(wèi)星．連續(xù)體結(jié)構(gòu)拓撲優(yōu)化方法評述 ［J］．航空工程進展，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］張彬．水平集方法的改進及其在流體形狀最優(yōu)控制中的應用 ［D］．西安：西安交通大學，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］陶文銓．計算傳熱學的近代進展 ［M］．北京：科學出版社，2000：1－427．

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

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

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

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

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

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

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