董馨, 劉小民

(西安交通大學(xué)能源與動力工程學(xué)院, 710049, 西安)

水平集方法具有優(yōu)秀的運動界面捕捉能力[1]。任曉輝和封建湖提出了一種不用求解復(fù)雜Hamilton-Jacobi方程的水平集方法,在對Heaviside函數(shù)正則化處理中,基于獲得的形狀導(dǎo)數(shù)和拓?fù)鋵?dǎo)數(shù)的相關(guān)信息加快收斂速度[2];段獻(xiàn)葆和秦新強結(jié)合了傳統(tǒng)的形狀靈敏度分析與改進(jìn)的水平集方法,發(fā)展了一種無需網(wǎng)格重新剖分的水平集方法,減少了運算工作量[3];張彬和劉小民提出了一種可以保持界面位置不動的新型隱式重新初始化的水平集方法,這種新型的水平集方法只需使時間步長滿足原始的CFL條件,保證了界面附近水平集函數(shù)節(jié)點值符號不變,從而實現(xiàn)高效的優(yōu)化[4];Dunning等基于水平集拓?fù)鋬?yōu)化方法發(fā)展了適用于非結(jié)構(gòu)化三維網(wǎng)格的快速行進(jìn)法和迎風(fēng)格式,增加水平集方法的魯棒性和高效性[5];Duan等提出了一種自適應(yīng)網(wǎng)格方法來求解不可壓縮Navier-Stokes的拓?fù)鋬?yōu)化問題,將隱式表示的材料分布信息作為設(shè)計變量,并將計算量集中在界面附近,從而顯著降低了計算代價[6]。Shintaro等提出了一種基于水平集方法的自由灰度拓?fù)鋬?yōu)化方法,在結(jié)構(gòu)邊界處運用可適應(yīng)網(wǎng)格,它通過移動歐拉網(wǎng)格的節(jié)點來保持水平集函數(shù),從而改善網(wǎng)格質(zhì)量[7]。但是他們沒有對水平集函數(shù)演化中的參數(shù)和時間步長參數(shù)進(jìn)行選擇分析和研究,在保證優(yōu)化精度的前提下,合理的參數(shù)可使改進(jìn)的水平集方法實現(xiàn)高效計算,充分發(fā)揮該方法的優(yōu)點。

通過對改進(jìn)的水平集方法[8]中的參數(shù)進(jìn)行選擇和分析,研究優(yōu)化過程中迭代和收斂參數(shù)選擇對優(yōu)化進(jìn)程的影響,使計算能夠快速、穩(wěn)定收斂,達(dá)到最優(yōu),同時保證優(yōu)化結(jié)果的精度。改進(jìn)的水平集方法采用了保持界面位置不動的新型隱式重新初始化水平集方法,結(jié)合無需樣條參數(shù)化網(wǎng)格重新劃分的水平集優(yōu)化方法。

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

Ω代表流體域,它是一個具有連續(xù)邊界Γ=?Ω的開區(qū)域,建立二維不可壓縮Stokes方程如下

式中:u為流體的流動速度;p為流體壓強;f為體積力;ν為動力黏性系數(shù)。ΓS是流體與固體間的交界面;ΓD是除ΓS外的Dirichlet邊界,有Γ=ΓS∪ΓD,u0是邊界ΓD上已知的速度分布。

本文研究的優(yōu)化問題可描述為

(2)

式中:u滿足狀態(tài)式(1),D是優(yōu)化問題的工作區(qū)域。

水平集函數(shù)φ(x,t)滿足

水平集函數(shù)對運動界面進(jìn)行追蹤的基本思想如下:用高一維空間中的水平集函數(shù)去描述低一維空間中的運動邊界,那么水平集函數(shù)的零等值線即為所追蹤的運動邊界。

水平集函數(shù)演化過程中區(qū)域界面始終表達(dá)為

φ(x,t)=0

(4)

對上式兩端同時對t求物質(zhì)導(dǎo)數(shù),可得

(5)

取Vn=x′(t)為界面處的水平集法向速度,水平集演化方程為

(6)

2 數(shù)值方法

2.1 優(yōu)化步驟

基于改進(jìn)水平集方法的流體拓?fù)鋬?yōu)化執(zhí)行過程如下,其中n代表第n次優(yōu)化迭代。

步驟1初始化原始流體區(qū)域,在設(shè)計區(qū)域內(nèi)將水平集函數(shù)初始化為符號距離函數(shù);

步驟2對流體區(qū)域進(jìn)行網(wǎng)格重新劃分;

步驟3求解流動控制方程和優(yōu)化共軛方程;

步驟4求解水平集法向速度Vn;

步驟5將Vn由界面擴(kuò)展到整個設(shè)計區(qū)域;

步驟6執(zhí)行面積約束;

步驟7求解水平集方程,實現(xiàn)水平集函數(shù)演化;

步驟8將水平集函數(shù)重新初始化為符號距離函數(shù)。

2.2 時間步長

在水平集函數(shù)求解過程中,為了防止數(shù)值振蕩,本文方程的離散采用一階Essentially Non-Oscillatory(ENO)格式[9-10]。下面給出上述二維問題的離散形式

式中:x和y分別代表網(wǎng)格的橫向和縱向坐標(biāo)。下標(biāo)i代表x方向第i個網(wǎng)格節(jié)點,下標(biāo)j代表y方向第j個網(wǎng)格節(jié)點。

(8)

+=

-=

其中x+=max(x,0)

式中:S(φ)是符號函數(shù);ε是光滑參數(shù),(εi,j)2采用中心差分進(jìn)行數(shù)值離散,離散形式為

其中時間步長滿足下式

2.3 水平集函數(shù)重新初始化

隱式重新初始化方法是通過迭代求解下式的偏微分方程來實現(xiàn)的[11-12]。

φ(x,0)=φ0(x)

(15)

為了數(shù)值求解的需要,對S(φ0)光滑化

ε的取值為

(17)

3 數(shù)值算例

圖1 拓?fù)鋬?yōu)化的工作區(qū)域與邊界條件

圖2給出了翼型優(yōu)化達(dá)到最優(yōu)時的流線分布情況,從圖2中可以看出本文結(jié)果與文獻(xiàn)[15]結(jié)果吻合良好。

(a)γ=0.8 (b)γ=0.85

(c)γ=0.9 (d)γ=0.95圖2 最優(yōu)形狀及其流線分布

3.1 水平集函數(shù)進(jìn)化次數(shù)的分析

在水平集函數(shù)演化過程中,每次優(yōu)化迭代時進(jìn)行多次水平集函數(shù)進(jìn)化會帶來一定的偏差,因此進(jìn)行水平集函數(shù)重新初始化,以消除水平集函數(shù)在演化過程中偏離符號距離函數(shù)的現(xiàn)象。

引入?yún)?shù)m定義為每次迭代時水平集函數(shù)的進(jìn)化次數(shù)。分別選取m值為10、8、6、4、2和1這6組數(shù)據(jù)進(jìn)行研究,如圖3、圖4所示,分別為流體區(qū)面積比和其波動值隨著迭代步數(shù)n的變化。為了實現(xiàn)在優(yōu)化前期加速計算,在優(yōu)化需要精確結(jié)果的地方提高準(zhǔn)確度,使計算程序能夠達(dá)到自動選擇合理參數(shù)的目的,將水平集函數(shù)進(jìn)化次數(shù)m與迭代步數(shù)n的關(guān)系分別采用分段線性函數(shù)(式(18))、三次函數(shù)(式(19))和冪函數(shù)(式(20))描述如下

(a)m取單一值

(b)m取函數(shù)圖3 不同m取值流體區(qū)面積比隨迭代步數(shù)的變化

圖4 不同m取值時流體區(qū)面積比 的波動值隨迭代步數(shù)的變化

由圖3可看出,流體區(qū)面積由初始面積逐漸上升至約束面積,之后基本保持不變直到迭代收斂。當(dāng)m取值分別為10和8時,迭代收斂后在面積比為0.95附近呈現(xiàn)鋸齒狀的波動,原因是較大的m值導(dǎo)致結(jié)構(gòu)參數(shù)(面積比和目標(biāo)函數(shù))變化較大,面積增加速度太快而導(dǎo)致流場面積大于目標(biāo)面積,較大面積對應(yīng)著較小目標(biāo)函數(shù),約束的施加使得面積反向變化,面積減小,目標(biāo)函數(shù)增大。繼續(xù)演化,面積又增大到大于目標(biāo)面積,如此反復(fù),形成了鋸齒狀演化曲線。參數(shù)m分別取值為6、4、2和1時都可以達(dá)到穩(wěn)定收斂,收斂速度依次下降,收斂時的迭代步數(shù)與計算時間見表1。

表1 m取單一值時的收斂指標(biāo)

圖4中面積比的波動值,即迭代第n步與n-1步的面積比差值的絕對值。m取10和8時的波動值較大,且收斂之后的波動也較大。而其他參數(shù)的波動值都隨著迭代步數(shù)大體呈下降趨勢,尤其在收斂時有近乎豎直的下降梯度,收斂之后波動值接近0并呈現(xiàn)穩(wěn)定趨勢。

目標(biāo)函數(shù)隨迭代步數(shù)的收斂如圖5所示,目標(biāo)函數(shù)值隨著迭代次數(shù)的增加而下降,達(dá)到0.85時幾乎不再隨著迭代次數(shù)發(fā)生變化,此時翼型繞流能量耗散值最小,翼型結(jié)構(gòu)為最優(yōu)拓?fù)浣Y(jié)構(gòu)。當(dāng)m取值為10和8時,收斂到最優(yōu)結(jié)果時產(chǎn)生了波動,收斂情況不理想。線性函數(shù)、三次函數(shù)與冪函數(shù)均在第13步時達(dá)到收斂,其中三次函數(shù)在優(yōu)化過程中收斂速度最快,這是由于它在收斂前期有較大的m,收斂后m逐漸減小為1,使得收斂趨于穩(wěn)定。

(a)m取單一值

(b)m取函數(shù)圖5 不同m取值下目標(biāo)函數(shù)隨迭代步數(shù)的變化

3.2 時間步長的分析

時間步長的選擇對于程序的收斂速度和計算速度起著至關(guān)重要的作用。時間步長過大會導(dǎo)致不滿足courant-Friedrichs-Lewy(CFL)條件,程序計算不穩(wěn)定,最后達(dá)不到收斂指標(biāo);而時間步長過小又會導(dǎo)致計算耗時、收斂速度太慢。因此,選取合適的時間步長尤為重要。

本文計算中時間步長設(shè)置為:Δt=kmin(Δx,Δy),其中k為常數(shù),為滿足CFL條件,k≤0.5。k分別選取0.5、0.2和0.1時進(jìn)行計算研究,流體區(qū)面積比和波動值隨迭代步數(shù)的關(guān)系分別如圖6和圖7所示,迭代收斂步數(shù)和計算時間見表2。為了實現(xiàn)優(yōu)化前期快速達(dá)到最優(yōu)值,優(yōu)化后期穩(wěn)定收斂,使計算程序能夠自動選擇合理參數(shù),將描述時間步長的參數(shù)k與迭代步數(shù)n的關(guān)系分別表達(dá)為線性函數(shù)和指數(shù)函數(shù),公式如下

k=-0.006n+0.55,n≤99

(21)

k=0.833e-0.03n

(22)

由圖6可看出,流體區(qū)面積由初始面積逐漸上升,之后由于加入了面積約束,保持面積比基本不變直到迭代收斂。圖7中的波動值大體呈現(xiàn)下降趨勢,最后接近0穩(wěn)定收斂,隨著k取值的減小,波動更平緩但收斂更慢。

(a)k取單一值

(b)k取函數(shù)圖6 不同k取值流體區(qū)面積比隨迭代步數(shù)的變化

圖7 不同k取值時流體區(qū)面積比的波動值隨迭代步數(shù)的變化

圖8為目標(biāo)函數(shù)隨迭代步數(shù)的變化,目標(biāo)函數(shù)隨著迭代次數(shù)的增加而下降,最后逐漸穩(wěn)定保持目標(biāo)函數(shù)不變,能量耗散數(shù)為0.85,此時可得到最優(yōu)拓?fù)浣Y(jié)構(gòu);線性函數(shù)與指數(shù)函數(shù)分別在第18步與第14步達(dá)到收斂,由于它們隨著迭代步數(shù)變化,因此在計算前期可以達(dá)到快速收斂的目的,在后期防止因為時間步長過大引起的數(shù)值不穩(wěn)定現(xiàn)象。

(a)k取單一值

(b)k取函數(shù)圖8 不同k取值下目標(biāo)函數(shù)值隨迭代步數(shù)的變化

k迭代步數(shù)迭代時間/s0.5185 4150.24314 7100.18227 860

(a)參數(shù)改進(jìn)后的流體區(qū)面積比

(b)參數(shù)改進(jìn)后的目標(biāo)函數(shù)圖9 參數(shù)改進(jìn)后的面積比及目標(biāo)函數(shù)隨迭代步數(shù)的變化

綜合選取兩個參數(shù)對優(yōu)化過程的影響規(guī)律,流體區(qū)面積比及目標(biāo)函數(shù)隨迭代步數(shù)的變化情況如圖9所示。從圖9a中可以看出,流體區(qū)面積隨迭代步數(shù)的增加逐漸上升至面積約束值,達(dá)到穩(wěn)定收斂,收斂速度最慢的為只改進(jìn)k的方案,其次為只改進(jìn)m的方案,最快的為同時改進(jìn)m和k的選取方式。由圖9b可得,目標(biāo)函數(shù)在優(yōu)化前期下降速度最快,接近收斂時目標(biāo)函數(shù)變化逐漸減小,直至穩(wěn)定收斂。同時改進(jìn)m和k的選取方式是在優(yōu)化之初不需要精確的收斂數(shù)值的情況下,增加每一次迭代時水平集函數(shù)的進(jìn)化次數(shù)從而提高計算速度,同時時間步長增加使迭代速度加快,隨著計算進(jìn)行,由于需要獲取精確的收斂數(shù)值,迭代速度減慢。改進(jìn)前與改進(jìn)后的收斂指標(biāo)見表3。參數(shù)選取方式如下:改進(jìn)前m=1,k=0.1;只改進(jìn)m的選取方式為三次函數(shù),k=0.1不變;只改進(jìn)k取指數(shù)函數(shù),m=1;改進(jìn)m選取為三次函數(shù),k選取指數(shù)函數(shù)。經(jīng)過合理選擇m和k,收斂速度提高96%。

表3 參數(shù)改進(jìn)前后的收斂指標(biāo)

4 結(jié) 論

針對低雷諾數(shù)流動條件下二維翼型拓?fù)鋬?yōu)化,引入了水平集函數(shù)進(jìn)化次數(shù)m,研究了不同m時面積比與迭代步數(shù)的變化關(guān)系,并分析了收斂指標(biāo)與波動值;研究了不同時間步長時面積比與迭代步數(shù)的變化關(guān)系,揭示了波動值與迭代步數(shù)的變化關(guān)系,并列出了收斂指標(biāo)。獲得的主要結(jié)論如下。

(1)在相同時間步長下,隨著水平集函數(shù)進(jìn)化次數(shù)m的增大,收斂速度變快,面積比的波動值變大。當(dāng)m>6,具有收斂速度快和面積比跳躍性大的特點,但達(dá)到最優(yōu)面積比后收斂變得不穩(wěn)定,無法得到穩(wěn)定的最優(yōu)形狀?？稍趦?yōu)化初期選用較大m,根據(jù)迭代步數(shù)的增長減小該值。

(2)在相同水平集函數(shù)進(jìn)化次數(shù)下,當(dāng)時間步長增大時,收斂速度和波動值都隨著迭代步數(shù)的增大而增大,達(dá)到最優(yōu)面積比時穩(wěn)定收斂。

(3)針對所需優(yōu)化精確結(jié)果和收斂速度的不同,可在不同面積比下選取適當(dāng)m和k,在最初優(yōu)化時采用較大的水平集進(jìn)化次數(shù)和時間步長,根據(jù)迭代步數(shù)增加逐漸采用較小的水平集進(jìn)化次數(shù)和時間步長,同時實現(xiàn)快速迭代和較容易優(yōu)化收斂的結(jié)果。