水平集與靈敏度分析相結合的流體形狀識別與優化方法

劉小民，張彬，張卓颯

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

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

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

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

1 數學模型

1．1 問題描述

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

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

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

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

1．2 靈敏度分析

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

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

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

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

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

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

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

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

2 改進的水平集方法

2．1 水平集函數的演化

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

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

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

2．2 法向速度擴展

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

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

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

將式（12）、式（13）結合起來，先在流體域直接求出邊界上較準確的水平集法向速度，然后運用快速行進法將法向速度由流體域向未定義的固體域擴展，最后再運用求解偏微分方程的方法向整個水平集函數求解區域進行擴展。可見，本文方法在不做其他處理的情況下，直接、較準確地求出邊界上的法向速度，且將界面附近的法向速度十分準確地沿界面的法線方向擴展到整個水平集函數求解區域。

3 數值方法

3．1 無需樣條參數化的重新網格劃分方法

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

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

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

式中：ξ、η滿足

這樣，式（14）變為

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

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

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

（2）給出初始區域Ω0，將Ω0初始化為符號距離函數；

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

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

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

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

（7）在結構化網格上運用差分方法求解水平集方程，并進行水平集函數的演化；

（8）3次優化迭代后進行1次水平集函數重新初始化，然后回到步驟（3）進一步循環，直到收斂為止。

3．3 有效性驗證

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

根據雷諾準則，任意2個符合雷諾準則的流動，只要滿足雷諾數相同條件，其流動情況一致。本文研究的流動問題均滿足雷諾準則，所以算例中只給出了流動雷諾數。考慮到對稱性，特征尺度L為1／2進口寬度，特征速度U 為進口處平均速度。將研究問題的工作區域D描述為四邊形區域，流體域用Ω表示。在區域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。優化問題的目標函數為式（3），狀態約束為式（2），流場求解區域為串列雙圓繞流區域，如圖1所示，由定理2給出目標函數的靈敏度分析結果。

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

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

圖2 重新網格劃分方法的網格劃分和速度矢量分布

圖3 浸沒邊界方法的網格劃分和速度矢量分布

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

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

4 數值算例

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可以看出：當迭代20步時，初始形狀發生了拓撲變化，表明本文方法可以有效處理帶有拓撲變化的優化問題，Stokes流動形狀識別過程中得到的最終形狀與目標形狀很接近（見圖6）。目標函數在形狀識別過程中的變化如圖8所示。由圖8可以看出，初始迭代時目標函數變化較大，之后變化減緩，直至保持不變為止，由此獲得最優結果。

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

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

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

圖7 Stokes流動形狀識別過程

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

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

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

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

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

5 結 論

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

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

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

［1］ 夏天翔，姚衛星．連續體結構拓撲優化方法評述 ［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，FEDKIW 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］張彬．水平集方法的改進及其在流體形狀最優控制中的應用 ［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．

［本刊相關文獻鏈接］

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

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

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

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

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

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