李旭，周洲，薛臣

西北工業(yè)大學 航空學院，西安 710072

工程實際中存在大量的流固耦合問題，物體存在復雜的運動或變形，這給數(shù)值求解帶來了困難。傳統(tǒng)的貼體網(wǎng)格方法已經(jīng)實現(xiàn)了對此類問題的計算，如動網(wǎng)格方法、嵌套網(wǎng)格方法。但是這些方法都不可避免地需要網(wǎng)格的運動，這就帶來了新的問題。當物體出現(xiàn)大變形或大位移時，動網(wǎng)格方法計算可能會出現(xiàn)負網(wǎng)格，而嵌套網(wǎng)格方法則容易出現(xiàn)孤兒網(wǎng)格[1]，這都會影響計算的穩(wěn)定性。

浸入邊界法[2-3]采用笛卡爾網(wǎng)格進行計算，物面用一系列離散的Lagrange點表示，通過邊界對周圍流體網(wǎng)格施加力或直接施加速度邊界條件來等效物體對流體的作用。在對運動邊界進行計算時，只需要改變Lagrange點的作用區(qū)域，并不需要流體網(wǎng)格的運動，避免了傳統(tǒng)貼體網(wǎng)格方法網(wǎng)格的變形或運動，不僅保證了網(wǎng)格質(zhì)量，還降低了復雜幾何的網(wǎng)格生成難度。經(jīng)過多年的發(fā)展，目前已經(jīng)出現(xiàn)了許多種不同形式[4-5]的浸入邊界法。Goldstein的反饋力浸入邊界法(又稱虛擬邊界法，Virtual Boundary Method)正是這些方法中的一種[6-7]。該方法基于非定常Navier-Stokes(N-S)方程提出，其吸取了控制理論中反饋的思想，分別采用積分和比例環(huán)節(jié)，利用物面的速度誤差計算力源項，并利用力源不斷對流場進行反饋，最終實現(xiàn)無滑移邊界條件。該方法形式簡單，不用對物體內(nèi)外的流體網(wǎng)格單元進行區(qū)分，也不用對已有求解器進行大的改動，因此方便使用[8-9]，已經(jīng)成功應用到了柔性變形[10-12]、撲翼運動[13-15]、葉柵顫振[16-19]、多相流計算[19]等問題上。

Goldstein的反饋力浸入邊界法含有對時間的積分項，這也是目前國內(nèi)外學者將該方法與非定常N-S方程求解結合的原因之一。對控制方程的求解采用的方法包括譜方法[9]、投影法[15,20]、格子玻爾茲曼方法[21]等，多為顯式時間推進。Goldstein的反饋力方法在應用時有兩個不便：第一，積分和比例環(huán)節(jié)的反饋系數(shù)需要人為給定；第二，計算存在穩(wěn)定性問題。原始的反饋力浸入邊界法在進行顯式計算時，出于穩(wěn)定性的要求，求解方法有非常嚴格的時間步長限制[20,22]。顯式計算本來就有穩(wěn)定性要求，采用反饋力浸入邊界法后，受反饋力系數(shù)取值的影響，計算的時間步將必須取得很小，不然會引起發(fā)散，這降低了求解的效率。文獻[21]指出Lee研究了在采用反饋力浸入邊界法時不同時間推進方法的穩(wěn)定性范圍，對比了Adams-Bashforth格式與Runge-Kutta格式的差異，結果表明Runge-Kutta類方法穩(wěn)定性較強，但CFL(Courant-Friedrichs-Lewy)數(shù)仍不能超過1。Shoele和Zhu[12]將加權隱式Crank-Nicolson格式與反饋力浸入邊界法結合，但反饋力源項仍然按照時間步進行更新，為保證格式的有界性，時間步長仍要受到限制。

由此可知，顯式求解制約了反饋力浸入邊界法發(fā)揮作用。與此相對，隱式方法通常是無條件穩(wěn)定的，計算時可采用更大的時間步，更加適合實際的工程計算。因此，將隱式求解方法與反饋力浸入邊界法結合是一種不錯的選擇。另外，對于一些不可壓的定常流動，往往可直接利用松弛迭代方法對不含時間項的N-S方程進行求解，如果能將反饋力浸入邊界法與這類定常求解方法結合，將進一步拓展該方法的使用范圍，但目前還未見到這方面研究公開發(fā)表。

基于以上思考，本文對Goldstein的反饋力浸入邊界法進行了改進，將原始方法中的時間積分項改為對迭代次數(shù)的求和。與原始方法相比，改進的方法不再含有與時間相關的參數(shù)。本文對一系列的定常與非定常算例進行了計算，結果表明改進后的方法可同時用于定常與非定常N-S方程的求解，這為以后反饋力浸入邊界法的應用提供了新的選擇。

1 數(shù)值計算方法

1.1 控制方程

非定常不可壓的N-S方程組為

(1)

式中：ρ為密度;p為壓力；μ為層流黏性系數(shù)；F(x,t)為流體的力源項，其計算公式為

(2)

其中：f(xs,t)為物體的力源項；δ(x-xs)為Delta函數(shù)[23]，積分前的負號表示物體的力源與流體的力源方向相反。Goldestein的反饋力浸入邊界法計算物體的力源項的表達式為

β′(u(x,t)-u(xs,t))

(3)

式中：u(x,t)為流體計算得到的物面速度，可通過Delta函數(shù)計算得到；u(xs,t)為物面真正的速度。式(3)含有兩個待定參數(shù)：α′和β′，需要人為給定，控制原理如圖1所示。

圖1 反饋力浸入邊界法原理Fig.1 Principle of feedback forcing immersed boundary method

圖1中的比例積分環(huán)節(jié)代表式(3)，通過每個時間步的速度誤差反饋，使|u(x,t)-u(xs,t)|趨于0，從而實現(xiàn)無滑移邊界條件。

Delta函數(shù)是物面節(jié)點與流場網(wǎng)格單元信息交流的關鍵，其形式并不唯一，本文選用的Delta函數(shù)為

δ(r)=

(4)

阻力系數(shù)的計算參考文獻[20]，表達式為

(5)

式中：右邊第1項表示物面上力源的積分，S指物體表面，fx為反饋力法計算物面上一點x方向的動量源項；第2項反映運動加速度的效應，V表示物體的體積；uc,x表示物體質(zhì)心平動速度在x方向的分量。同理可得物體其他方向的受力。本文所有算例的密度均取為ρ=1.0 kg/m3。需說明的是，本文在采用式(2)計算時，反饋力參數(shù)均取正值。

1.2 反饋力浸入邊界法的改進

Goldstein的反饋力浸入邊界法(后文統(tǒng)稱原始方法)與顯式求解方法結合將嚴重限制計算的效率，難以在工程應用中推廣。而隱式格式通常沒有穩(wěn)定性的限制，將反饋力浸入邊界法與隱式格式結合，是一種有效的解決辦法，也能更好地發(fā)揮反饋力浸入邊界法的優(yōu)點，這正是本文改進方法的出發(fā)點。

在進行數(shù)值計算時，Goldstein反饋力的計算公式可寫為

β′(u(x,t)-u(xs,t))

(6)

式中：i表示時間步數(shù)，t=N′·Δt，N′表示當前的時間步數(shù)。即式(3)中對時間的積分項化為了對時間步的求和。通過對不同時間步的速度誤差進行反饋來最終滿足物面邊界條件，這是目前國內(nèi)外學者計算時采用的通用方式。

從純反饋的原理出發(fā)，通過對式(6)進行分析，本文認為既然α′和β′需要人為給定，可以對其重新進行定義，令α=α′·Δt，β=β′，則可得

β(u(x,N)-u(xs,N))

(7)

式中：N表示當前計算的迭代次數(shù)。式(7)不含時間相關的參數(shù)，積分環(huán)節(jié)轉(zhuǎn)化為不同次迭代速度誤差的求和，而非式(6)中不同時間步速度誤差的求和，這是與原始方法最大的不同之處。相同的是，仍然是通過對流場速度誤差的反饋來計算力源項。相比式(6)，改進后的方法有兩方面的好處：第一，由于不含時間相關參數(shù)，可直接與一些定常N-S方程迭代求解方法結合；第二，在非定常計算時，通過迭代次數(shù)來計算力源項更加適合隱式求解，因為隱式推進一個時間步內(nèi)往往有多次迭代的要求。

進一步來看，方程求解的迭代過程就是流域內(nèi)的網(wǎng)格單元速度不斷調(diào)整，達到收斂，從而與給定邊界條件相適應的過程，這正好與控制系統(tǒng)中的反饋作用相類似。當實際的輸出量與給定值存在偏差時，系統(tǒng)就不斷通過誤差反饋來減小偏差，從而使得輸出量等于給定值。也就是說，直接利用迭代過程中的速度誤差來計算力源項，并將其反饋給流體以期實現(xiàn)無滑移邊界條件是可行的，反饋原理的本質(zhì)不應與N-S方程的具體求解方法有關。

因此，改進的方法將具有更大的適用性。本文改進的方法不僅能與含有時間項的N-S方程進行結合，還能與不含時間項的N-S方程進行迭代求解，本文將對這一觀點進行驗證。

Fluent是目前一款應用較廣泛的CFD商業(yè)軟件，其基于有限體積法，求解器分為兩大類：密度基和壓力基。密度基求解器定常與非定常計算均按照時間推進求解，壓力基求解器則不同。對于定常計算，壓力基求解器求解的是不含時間項的N-S方程，非定常計算才含有時間項，且為隱式時間推進[24]。本文選用壓力基求解器，基于軟件中的用戶自定義函數(shù)(User-Defined Functions， UDFs)，通過一系列定常與非定常算例的計算，驗證本文改進反饋力浸入邊界法的有效性。

各個算例空間離散均采用二階精度。對于非定常計算，時間推進則采用一階隱式。定常計算可直接根據(jù)迭代進行力源的更新。對于非定常隱式計算，改進方法的實施步驟如下：

步驟1初始化，讀入物面Lagrange點坐標。

步驟2每個時間步開始，如果物面位置發(fā)生變化，則更新Lagrange點，并在本時間步的迭代過程中保持其坐標不變。根據(jù)迭代次數(shù)，利用改進方法進行力源項的更新，并將其添加到動量方程中，參與控制方程的求解。

步驟3待本時間步迭代收斂，利用最終的力源進行本時間步物體受力的計算。

步驟4重復步驟2和步驟3，進行下一時間步的計算。

反饋力源的計算在DEFINE_ADJUST宏中進行，并通過DEFINE_SOURCE宏將其加給動量方程。關于UDF更多的使用細節(jié)，可參考Fluent自帶的UDF Manual進行了解。

2 改進反饋力浸入邊界法的驗證

2.1 靜止圓柱繞流

對靜止圓柱繞流的計算是驗證數(shù)值方法的經(jīng)典算例。首先，本文對雷諾數(shù)Re=40時靜止圓柱的流動進行了計算，以此來驗證本文改進的反饋力方法在定常N-S方程求解時的有效性。采用Fluent壓力基求解器中兩種不同的定常計算方法，控制方程求解均不含時間項，進行松弛迭代求解。此時由于計算不再是時間推進求解，Goldstein的原始公式不再適用。

圓柱表面均勻分布75個Lagrange點，所在流域采用均勻網(wǎng)格，網(wǎng)格間距Δx的選取采用文獻[25]的方法,計算網(wǎng)格量為29 380。

如圖2所示，計算源項時需要對網(wǎng)格單元進行搜索，為減少搜索的范圍，流體網(wǎng)格被分成了內(nèi)外兩個域。圓柱阻力系數(shù)的定義為

(8)

式中：U∞為自由來流速度，U∞=5.84 m/s；D為圓柱直徑，D=1×10-4m。采用矩形計算域，60D×40D。入口采用速度入口，出口采用壓力出口，上下邊界設置為對稱邊界條件，空間離散采用二階精度。將分別基于Simplec和Couple算法對靜止圓柱繞流進行計算。其中Simplec為分離式求解(經(jīng)典的壓力修正算法)，Couple則為耦合式求解，以此來檢驗改進方法對不同求解格式的可靠性。

圖2 計算網(wǎng)格分布Fig.2 Grid distribution in simulation

確定好網(wǎng)格和邊界條件后，就需要對α和β進行取值。雖然反饋原理的本質(zhì)不變，但α和β的取值會影響到反饋的效果。α和β選取過大，會引起計算振蕩，甚至導致顯式計算發(fā)散；選取過小，會使得收斂變慢[9]。正如前文所說，目前并沒有理論的方法確定α和β的最優(yōu)值，這與計算狀態(tài)和求解器本身相關，本文選取的值也是經(jīng)過嘗試確定的。

本文研究了α和β的取值對計算結果的影響。由于物體靜止，式(5)中的第2項為0?；赟implec算法，對比不同反饋力參數(shù)對計算收斂的影響，結果如圖3所示。

圖3 反饋力參數(shù)對收斂的影響Fig.3 Influence of feedback forcing parameters for iteration

在3種不同取值下，阻力系數(shù)均收斂到了相同值。取α=1.0，β=0.1時，阻力系數(shù)的振蕩比較明顯，即超調(diào)量大；取α=0.01，β=0.001時，阻力系數(shù)波動幅值較小，但達到平穩(wěn)速度較慢。因此，確定本算例反饋力參數(shù)α=0.1，β=0.01。

在計算過程中，每進行一次迭代，反饋力的源項都會得到更新，再將式(7)的力源分布到周圍的流體網(wǎng)格上，進行下一次迭代，如此反復進行，直到計算收斂。

取Simplec計算收斂后流場的流線如圖4所示。由于雷諾數(shù)較小，可以看出圓柱后方形成了兩個對稱、穩(wěn)定的分離渦。表1是計算得到的圓柱繞流特征參數(shù)(阻力系數(shù)CD和無量綱回流區(qū)長度L/D)，分別是兩種求解方法與改進的反饋力浸入邊界法結合(表中用Simplec和Couple表示)，并與文獻[26-28]的結果對比。通過圓柱繞流特征參數(shù)的對比可知，本文計算得到的圓柱阻力系數(shù)和回流區(qū)的長度與文獻符合較好。

再取圓柱表面的壓力系數(shù)Cp與文獻[29]對比，如圖5所示。可以看出，本文改進方法計算的壓力系數(shù)與文獻[29]給出的計算結果符合很好，表明本文改進的反饋力浸入邊界法是有效的。不同于原始方法，本文改進的方法能用于不含時間項N-S方程的迭代求解，且反饋力系數(shù)取值對收斂結果沒有大的影響。由于Simplec計算求解速度快，本文以后的計算均采用Simplec算法。

圖4 圓柱流線(Re=40)Fig.4 Streamlines on circular cylinder at Re=40

表1 Re=40阻力系數(shù)和回流區(qū)長度對比

圖5 靜止圓柱壓力系數(shù)Fig.5 Pressure coefficient for stationary circular cylinder

以上計算求解的是不含時間項的N-S方程，為驗證改進方法在非定常隱式求解上的有效性，本文對Re=40時的圓柱繞流進行了非定常模擬。同樣基于Simplec算法，對含有時間項的N-S方程進行求解，對原始反饋力浸入邊界法與改進方法的收斂效果進行對比。

在進行隱式計算時，原始方法的源項是按照時間步進行計算的，則每個時間步只在最開始進行更新，一個時間步內(nèi)其值保持不變。本文改進方法的源項基于迭代次數(shù)計算，每一次迭代后源項均會得到修正，一個時間步內(nèi)其值一般是變化的。

時間步長取為0.002 s，一個時間步內(nèi)迭代20次，采用相同的計算設置，對兩種方法的收斂特性進行對比。圖6為阻力系數(shù)的收斂曲線。可以看出，采用原始方法進行隱式計算時收斂較慢，特別是在快接近穩(wěn)定值的時候，最終到500個時間步才完全收斂。而本文改進的方法則很快收斂，阻力在100個時間步后就達到平穩(wěn)。

參考文獻[14]，定義CFL數(shù)為

式中：Δx=2.96×10-6m，為圓柱周圍Euler網(wǎng)格尺寸。本文改進方法隱式計算與文獻顯式結果對比如表2所示。由于本文是隱式求解，與顯式計算[8, 14]相比，CFL數(shù)可以取的很大，這也體現(xiàn)了隱式求解的好處。

圖6 靜止邊界隱式計算收斂效率對比Fig.6 Comparison of efficiency for implicit iterations at stationary boundary

從數(shù)值求解上來看，在代數(shù)方程求解的過程中，源項影響的是方程組的常數(shù)項。隱式求解每一次迭代都修正源項，使得源項與速度場的同步性得到改善，能起到加速收斂的作用。也就是說，改進的方法更加適合非定常隱式求解。

表2 Re=40時隱式與顯式結果對比

2.2 振蕩圓柱

為驗證改進方法對運動邊界模擬的可行性，選取靜止流體中往復振動的圓柱[29-30]進行非定常計算，采用Simplec算法求解不可壓非定常N-S方程，通過與文獻對比來驗證本文方法在運動物體求解上的有效性。

本算例計算域的大小取為40D×40D，圓柱將在水平方向進行往復運動，平衡位置在原點處，運動方程為

x(t)=-Asin(2πft)

(9)

式中：x為圓柱水平方向位移；A為振動幅值；f為振動頻率。

在本算例中，影響計算結果的關鍵參數(shù)為雷諾數(shù)Re和KC(Keulegan-Carpenter)數(shù)，分別定義為

(10)

式中：umax為圓柱振動的最大速度，umax=2πfA。本文選取Re=100，KC=5。

參考文獻[8]，四周遠場邊界設置為無剪切的物面條件。由于采用隱式計算，一個周期只取360個時間步，每個時間步內(nèi)迭代次數(shù)為20次。反饋力參數(shù)取值見表3。

對于圓柱表面運動速度的計算，可直接對式(9) 求導獲得，然后將該速度代入式(7)中進行力源項的計算。在一個時間步內(nèi)，圓柱的位置不會變化，每迭代一次，流場速度得到更新，反饋力

表3 不同算例反饋力參數(shù)Table 3 Feedback parameters for different cases

源項就會進行更新，然后再將其添加到動量方程中，進行下一次迭代，直到該時間步收斂。

與2.1節(jié)相同，流場網(wǎng)格分為了兩個域，內(nèi)區(qū)域網(wǎng)格量為40 000，圓柱運動過程都在其中，外區(qū)域網(wǎng)格量為31 824，總計網(wǎng)格為71 824，物面有75個Lagrange點。在計算過程中，只對圓柱所在區(qū)域的網(wǎng)格進行搜索判斷以及反饋力源項的賦值。運動穩(wěn)定后，得到流場不同時刻的壓力及渦量如圖7所示。

圖7中渦量的正負分別用虛線和實線表示，當圓柱向遠離平衡位置方向運動時，其后會形成兩個強度相同、方向相反的對稱渦；而當圓柱向平衡位置返回時，其將分開之前形成的這對渦，并改變渦的旋轉(zhuǎn)方向，這種現(xiàn)象與文獻[30]一致。

圖7 不同相位下的渦量/壓力等值線Fig.7 Vorticity/pressure contours at different phase angles

取改進方法計算的圓柱一個周期的阻力曲線與文獻進行對比，如圖9所示?？梢钥闯?，阻力的最大、最小值并不是在特殊相位(t=0.25nT，n為正整數(shù)，T為振動周期)時取得，這與圓柱運動過程中流體的慣性有關，表現(xiàn)出了遲滯特性，計算得到的振動圓柱阻力的變化與文獻[30]符合較好。

圖8 運動邊界隱式計算效率對比Fig.8 Comparison of efficiency for implicit iterations at moving boundary

圖9 阻力系數(shù)變化曲線Fig.9 Time history of drag coefficient

以上結果表明本文改進的反饋力浸入邊界法可與隱式求解方法結合，對運動物體進行迭代計算，且收斂效果優(yōu)于原始方法。

2.3 運動橢圓翼

振動圓柱只有一個方向的運動，橢圓翼的運動要更加復雜，其不僅有平動，還有繞自身質(zhì)心的轉(zhuǎn)動，可用其來研究昆蟲翅膀運動時非定常流動機理[31]。本節(jié)將選取靜止流體中的橢圓翼來驗證本文改進方法對物體復雜運動的模擬能力。

橢圓翼運動如圖10[31]所示，質(zhì)心的平動和繞質(zhì)心的轉(zhuǎn)動規(guī)律為

(11)

(12)

阻力系數(shù)和升力系數(shù)的定義分別為

(13)

圖10 橢圓翼運動路徑[31]Fig.10 Flapping path of ellipse wing[31]

式中：c為橢圓翼的長軸長，取c=1 m。橢圓翼的短軸長b=0.25 m。一個周期T分為360個時間步，內(nèi)迭代次數(shù)為20。選取計算域為80c×80c，四周遠場邊界采用無剪切壁面，總網(wǎng)格量為164 160。由于運動比較復雜，物面取150個Lagrange點。

計算穩(wěn)定后，得到橢圓翼運動周期內(nèi)不同時刻的渦量圖如圖11所示。在向下運動過程中，橢圓翼前后緣會產(chǎn)生一對反向旋轉(zhuǎn)的渦(t=0.25T)，由于橢圓翼自身的旋轉(zhuǎn)作用，這對渦將逐漸靠近(t=0.44T)。在向上撲動過程中，之前形成的這對渦會從橢圓翼表面脫落(t=0.74T)，形成一對相互作用的偶極子，并向下運動(t=0.99T)。

得到橢圓翼的阻力和升力變化規(guī)律如圖12所示?？梢钥闯觯捎跈E圓翼流場非定常特性明顯，其氣動力變化并不是簡單的正弦曲線。另外，本文計算的橢圓翼的升力阻力變化與文獻[14]基本一致，且相比文獻[14]，計算沒有明顯的數(shù)值振蕩，表明可以利用本文改進的方法進行隱式計算，對復雜的非定常運動進行模擬。

圖11 一個周期內(nèi)不同時刻的渦量圖Fig.11 Vorticity fields during one period

圖12 撲動過程中的氣動力曲線Fig.12 Curves of aerodynamic coefficients during flapping

2.4 靜止圓球

為進一步驗證本文方法對三維物體模擬的可靠性，選取靜止圓球進行定常計算，求解不含時間項的N-S方程。采用Simplec算法，計算域大小為80D×60D×60D，D=0.001 m。

為減少網(wǎng)格量，本文在圓球所處的流域進行了局部網(wǎng)格加密，如圖13所示，加密后總的網(wǎng)格量為282萬。遠場邊界分別采用速度入口和壓力出口。

圓球表面的網(wǎng)格如圖14所示，采用三角形單元來離散球面，共計有3 860個單元。來流速度分別取15 m/s和25 m/s，對Re=150和Re=250時的圓球繞流進行計算，得到流場如圖15所示?？梢钥闯觯赗e=150時，繞圓球的流動是軸對稱的。但在Re=250時，雖然圓球本身是軸對稱的，但流動的對稱性減弱，x-y面的流動和x-z面的流動存在差異，流場只存在面對稱，即上下對稱。由于x-z面渦不對稱，圓球產(chǎn)生了側(cè)力，但流動依舊保持穩(wěn)態(tài)。

圖13 局部網(wǎng)格加密Fig.13 Local refining for computation mesh

圖14 圓球表面的Lagrange網(wǎng)格單元Fig.14 Lagrange mesh of sphere surface

圖15 不同雷諾數(shù)圓球繞流流線Fig.15 Streamlines for flow past sphere at different Re

本文計算得到圓球的受力與文獻對比如表4所示。可以看出，本文計算結果與文獻結果符合較好，在Re=250時，圓球產(chǎn)生了側(cè)力(表中用Cc表示)。表明本文改進的反饋力浸入邊界方法適合三維物體的計算。

表4 圓球氣動特性對比

3 結 論

本文對Goldstein的反饋力浸入邊界法進行了研究和改進，將原始方法中含有的對速度誤差的時間積分轉(zhuǎn)化為迭代過程中速度誤差的求和，利用多個定常與非定常算例驗證了改進方法的有效性。得到的主要結論有：

1) 可直接基于迭代次數(shù)對反饋力源項進行計算。與原始方法相比，由于不含時間相關參數(shù)，改進的反饋力浸入邊界法可與求解定常N-S方程的松弛迭代方法結合，用于定常流動計算。

2) 隱式時間推進可以避免反饋力浸入邊界法顯式求解時嚴格的時間步長限制。本文改進的反饋力浸入邊界法適合于隱式時間推進，可實現(xiàn)對非定常N-S方程的求解，且收斂特性要優(yōu)于原始方法。

3) 改進的方法與原始方法均是基于速度誤差反饋的思想求解力源項，但本文改進的反饋力浸入邊界法有更廣泛的適用性。