李旭，周洲，薛臣

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

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

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

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

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

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

1 數(shù)值計(jì)算方法

1.1 控制方程

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

(1)

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

(2)

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

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

(3)

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

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

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

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

δ(r)=

(4)

阻力系數(shù)的計(jì)算參考文獻(xiàn)[20]，表達(dá)式為

(5)

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

1.2 反饋力浸入邊界法的改進(jìn)

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

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

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

(6)

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

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

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

(7)

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

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

因此，改進(jìn)的方法將具有更大的適用性。本文改進(jìn)的方法不僅能與含有時(shí)間項(xiàng)的N-S方程進(jìn)行結(jié)合，還能與不含時(shí)間項(xiàng)的N-S方程進(jìn)行迭代求解，本文將對(duì)這一觀點(diǎn)進(jìn)行驗(yàn)證。

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

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

步驟1初始化，讀入物面Lagrange點(diǎn)坐標(biāo)。

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

步驟3待本時(shí)間步迭代收斂，利用最終的力源進(jìn)行本時(shí)間步物體受力的計(jì)算。

步驟4重復(fù)步驟2和步驟3，進(jìn)行下一時(shí)間步的計(jì)算。

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

2 改進(jìn)反饋力浸入邊界法的驗(yàn)證

2.1 靜止圓柱繞流

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

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

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

(8)

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

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

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

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

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

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

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

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

再取圓柱表面的壓力系數(shù)Cp與文獻(xiàn)[29]對(duì)比，如圖5所示?？梢钥闯觯疚母倪M(jìn)方法計(jì)算的壓力系數(shù)與文獻(xiàn)[29]給出的計(jì)算結(jié)果符合很好，表明本文改進(jìn)的反饋力浸入邊界法是有效的。不同于原始方法，本文改進(jìn)的方法能用于不含時(shí)間項(xiàng)N-S方程的迭代求解，且反饋力系數(shù)取值對(duì)收斂結(jié)果沒(méi)有大的影響。由于Simplec計(jì)算求解速度快，本文以后的計(jì)算均采用Simplec算法。

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

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

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

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

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

時(shí)間步長(zhǎng)取為0.002 s，一個(gè)時(shí)間步內(nèi)迭代20次，采用相同的計(jì)算設(shè)置，對(duì)兩種方法的收斂特性進(jìn)行對(duì)比。圖6為阻力系數(shù)的收斂曲線?？梢钥闯觯捎迷挤椒ㄟM(jìn)行隱式計(jì)算時(shí)收斂較慢，特別是在快接近穩(wěn)定值的時(shí)候，最終到500個(gè)時(shí)間步才完全收斂。而本文改進(jìn)的方法則很快收斂，阻力在100個(gè)時(shí)間步后就達(dá)到平穩(wěn)。

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

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

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

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

表2 Re=40時(shí)隱式與顯式結(jié)果對(duì)比

2.2 振蕩圓柱

為驗(yàn)證改進(jìn)方法對(duì)運(yùn)動(dòng)邊界模擬的可行性，選取靜止流體中往復(fù)振動(dòng)的圓柱[29-30]進(jìn)行非定常計(jì)算，采用Simplec算法求解不可壓非定常N-S方程，通過(guò)與文獻(xiàn)對(duì)比來(lái)驗(yàn)證本文方法在運(yùn)動(dòng)物體求解上的有效性。

本算例計(jì)算域的大小取為40D×40D，圓柱將在水平方向進(jìn)行往復(fù)運(yùn)動(dòng)，平衡位置在原點(diǎn)處，運(yùn)動(dòng)方程為

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

(9)

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

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

(10)

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

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

對(duì)于圓柱表面運(yùn)動(dòng)速度的計(jì)算，可直接對(duì)式(9) 求導(dǎo)獲得，然后將該速度代入式(7)中進(jìn)行力源項(xiàng)的計(jì)算。在一個(gè)時(shí)間步內(nèi)，圓柱的位置不會(huì)變化，每迭代一次，流場(chǎng)速度得到更新，反饋力

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

源項(xiàng)就會(huì)進(jìn)行更新，然后再將其添加到動(dòng)量方程中，進(jìn)行下一次迭代，直到該時(shí)間步收斂。

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

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

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

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

圖8 運(yùn)動(dòng)邊界隱式計(jì)算效率對(duì)比Fig.8 Comparison of efficiency for implicit iterations at moving boundary

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

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

2.3 運(yùn)動(dòng)橢圓翼

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

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

(11)

(12)

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

(13)

圖10 橢圓翼運(yùn)動(dòng)路徑[31]Fig.10 Flapping path of ellipse wing[31]

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

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

得到橢圓翼的阻力和升力變化規(guī)律如圖12所示。可以看出，由于橢圓翼流場(chǎng)非定常特性明顯，其氣動(dòng)力變化并不是簡(jiǎn)單的正弦曲線。另外，本文計(jì)算的橢圓翼的升力阻力變化與文獻(xiàn)[14]基本一致，且相比文獻(xiàn)[14]，計(jì)算沒(méi)有明顯的數(shù)值振蕩，表明可以利用本文改進(jìn)的方法進(jìn)行隱式計(jì)算，對(duì)復(fù)雜的非定常運(yùn)動(dòng)進(jìn)行模擬。

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

圖12 撲動(dòng)過(guò)程中的氣動(dòng)力曲線Fig.12 Curves of aerodynamic coefficients during flapping

2.4 靜止圓球

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

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

圓球表面的網(wǎng)格如圖14所示，采用三角形單元來(lái)離散球面，共計(jì)有3 860個(gè)單元。來(lái)流速度分別取15 m/s和25 m/s，對(duì)Re=150和Re=250時(shí)的圓球繞流進(jìn)行計(jì)算，得到流場(chǎng)如圖15所示。可以看出，在Re=150時(shí)，繞圓球的流動(dòng)是軸對(duì)稱的。但在Re=250時(shí)，雖然圓球本身是軸對(duì)稱的，但流動(dòng)的對(duì)稱性減弱，x-y面的流動(dòng)和x-z面的流動(dòng)存在差異，流場(chǎng)只存在面對(duì)稱，即上下對(duì)稱。由于x-z面渦不對(duì)稱，圓球產(chǎn)生了側(cè)力，但流動(dòng)依舊保持穩(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

本文計(jì)算得到圓球的受力與文獻(xiàn)對(duì)比如表4所示?？梢钥闯觯疚挠?jì)算結(jié)果與文獻(xiàn)結(jié)果符合較好，在Re=250時(shí)，圓球產(chǎn)生了側(cè)力(表中用Cc表示)。表明本文改進(jìn)的反饋力浸入邊界方法適合三維物體的計(jì)算。

表4 圓球氣動(dòng)特性對(duì)比

3 結(jié) 論

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

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

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

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