李旭,周洲,薛臣
西北工業(yè)大學(xué) 航空學(xué)院,西安 710072
工程實際中存在大量的流固耦合問題,物體存在復(fù)雜的運動或變形,這給數(shù)值求解帶來了困難。傳統(tǒng)的貼體網(wǎng)格方法已經(jīng)實現(xiàn)了對此類問題的計算,如動網(wǎng)格方法、嵌套網(wǎng)格方法。但是這些方法都不可避免地需要網(wǎng)格的運動,這就帶來了新的問題。當(dā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ì)量,還降低了復(fù)雜幾何的網(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)成功應(yīng)用到了柔性變形[10-12]、撲翼運動[13-15]、葉柵顫振[16-19]、多相流計算[19]等問題上。
Goldstein的反饋力浸入邊界法含有對時間的積分項,這也是目前國內(nèi)外學(xué)者將該方法與非定常N-S方程求解結(jié)合的原因之一。對控制方程的求解采用的方法包括譜方法[9]、投影法[15,20]、格子玻爾茲曼方法[21]等,多為顯式時間推進。Goldstein的反饋力方法在應(yīng)用時有兩個不便:第一,積分和比例環(huán)節(jié)的反饋系數(shù)需要人為給定;第二,計算存在穩(wěn)定性問題。原始的反饋力浸入邊界法在進行顯式計算時,出于穩(wěn)定性的要求,求解方法有非常嚴(yán)格的時間步長限制[20,22]。顯式計算本來就有穩(wěn)定性要求,采用反饋力浸入邊界法后,受反饋力系數(shù)取值的影響,計算的時間步將必須取得很小,不然會引起發(fā)散,這降低了求解的效率。文獻[21]指出Lee研究了在采用反饋力浸入邊界法時不同時間推進方法的穩(wěn)定性范圍,對比了Adams-Bashforth格式與Runge-Kutta格式的差異,結(jié)果表明Runge-Kutta類方法穩(wěn)定性較強,但CFL(Courant-Friedrichs-Lewy)數(shù)仍不能超過1。Shoele和Zhu[12]將加權(quán)隱式Crank-Nicolson格式與反饋力浸入邊界法結(jié)合,但反饋力源項仍然按照時間步進行更新,為保證格式的有界性,時間步長仍要受到限制。
由此可知,顯式求解制約了反饋力浸入邊界法發(fā)揮作用。與此相對,隱式方法通常是無條件穩(wěn)定的,計算時可采用更大的時間步,更加適合實際的工程計算。因此,將隱式求解方法與反饋力浸入邊界法結(jié)合是一種不錯的選擇。另外,對于一些不可壓的定常流動,往往可直接利用松弛迭代方法對不含時間項的N-S方程進行求解,如果能將反饋力浸入邊界法與這類定常求解方法結(jié)合,將進一步拓展該方法的使用范圍,但目前還未見到這方面研究公開發(fā)表。
基于以上思考,本文對Goldstein的反饋力浸入邊界法進行了改進,將原始方法中的時間積分項改為對迭代次數(shù)的求和。與原始方法相比,改進的方法不再含有與時間相關(guān)的參數(shù)。本文對一系列的定常與非定常算例進行了計算,結(jié)果表明改進后的方法可同時用于定常與非定常N-S方程的求解,這為以后反饋力浸入邊界法的應(yīng)用提供了新的選擇。
非定常不可壓的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)格單元信息交流的關(guān)鍵,其形式并不唯一,本文選用的Delta函數(shù)為
δ(r)=
(4)
阻力系數(shù)的計算參考文獻[20],表達式為
(5)
式中:右邊第1項表示物面上力源的積分,S指物體表面,fx為反饋力法計算物面上一點x方向的動量源項;第2項反映運動加速度的效應(yīng),V表示物體的體積;uc,x表示物體質(zhì)心平動速度在x方向的分量。同理可得物體其他方向的受力。本文所有算例的密度均取為ρ=1.0 kg/m3。需說明的是,本文在采用式(2)計算時,反饋力參數(shù)均取正值。
Goldstein的反饋力浸入邊界法(后文統(tǒng)稱原始方法)與顯式求解方法結(jié)合將嚴(yán)重限制計算的效率,難以在工程應(yīng)用中推廣。而隱式格式通常沒有穩(wěn)定性的限制,將反饋力浸入邊界法與隱式格式結(jié)合,是一種有效的解決辦法,也能更好地發(fā)揮反饋力浸入邊界法的優(yōu)點,這正是本文改進方法的出發(fā)點。
在進行數(shù)值計算時,Goldstein反饋力的計算公式可寫為
β′(u(x,t)-u(xs,t))
(6)
式中:i表示時間步數(shù),t=N′·Δt,N′表示當(dāng)前的時間步數(shù)。即式(3)中對時間的積分項化為了對時間步的求和。通過對不同時間步的速度誤差進行反饋來最終滿足物面邊界條件,這是目前國內(nèi)外學(xué)者計算時采用的通用方式。
從純反饋的原理出發(fā),通過對式(6)進行分析,本文認(rèn)為既然α′和β′需要人為給定,可以對其重新進行定義,令α=α′·Δt,β=β′,則可得
β(u(x,N)-u(xs,N))
(7)
式中:N表示當(dāng)前計算的迭代次數(shù)。式(7)不含時間相關(guān)的參數(shù),積分環(huán)節(jié)轉(zhuǎn)化為不同次迭代速度誤差的求和,而非式(6)中不同時間步速度誤差的求和,這是與原始方法最大的不同之處。相同的是,仍然是通過對流場速度誤差的反饋來計算力源項。相比式(6),改進后的方法有兩方面的好處:第一,由于不含時間相關(guān)參數(shù),可直接與一些定常N-S方程迭代求解方法結(jié)合;第二,在非定常計算時,通過迭代次數(shù)來計算力源項更加適合隱式求解,因為隱式推進一個時間步內(nèi)往往有多次迭代的要求。
進一步來看,方程求解的迭代過程就是流域內(nèi)的網(wǎng)格單元速度不斷調(diào)整,達到收斂,從而與給定邊界條件相適應(yīng)的過程,這正好與控制系統(tǒng)中的反饋作用相類似。當(dāng)實際的輸出量與給定值存在偏差時,系統(tǒng)就不斷通過誤差反饋來減小偏差,從而使得輸出量等于給定值。也就是說,直接利用迭代過程中的速度誤差來計算力源項,并將其反饋給流體以期實現(xiàn)無滑移邊界條件是可行的,反饋原理的本質(zhì)不應(yīng)與N-S方程的具體求解方法有關(guān)。
因此,改進的方法將具有更大的適用性。本文改進的方法不僅能與含有時間項的N-S方程進行結(jié)合,還能與不含時間項的N-S方程進行迭代求解,本文將對這一觀點進行驗證。
Fluent是目前一款應(yīng)用較廣泛的CFD商業(yè)軟件,其基于有限體積法,求解器分為兩大類:密度基和壓力基。密度基求解器定常與非定常計算均按照時間推進求解,壓力基求解器則不同。對于定常計算,壓力基求解器求解的是不含時間項的N-S方程,非定常計算才含有時間項,且為隱式時間推進[24]。本文選用壓力基求解器,基于軟件中的用戶自定義函數(shù)(User-Defined Functions, UDFs),通過一系列定常與非定常算例的計算,驗證本文改進反饋力浸入邊界法的有效性。
各個算例空間離散均采用二階精度。對于非定常計算,時間推進則采用一階隱式。定常計算可直接根據(jù)迭代進行力源的更新。對于非定常隱式計算,改進方法的實施步驟如下:
步驟1初始化,讀入物面Lagrange點坐標(biāo)。
步驟2每個時間步開始,如果物面位置發(fā)生變化,則更新Lagrange點,并在本時間步的迭代過程中保持其坐標(biāo)不變。根據(jù)迭代次數(shù),利用改進方法進行力源項的更新,并將其添加到動量方程中,參與控制方程的求解。
步驟3待本時間步迭代收斂,利用最終的力源進行本時間步物體受力的計算。
步驟4重復(fù)步驟2和步驟3,進行下一時間步的計算。
反饋力源的計算在DEFINE_ADJUST宏中進行,并通過DEFINE_SOURCE宏將其加給動量方程。關(guān)于UDF更多的使用細節(jié),可參考Fluent自帶的UDF Manual進行了解。
對靜止圓柱繞流的計算是驗證數(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。入口采用速度入口,出口采用壓力出口,上下邊界設(shè)置為對稱邊界條件,空間離散采用二階精度。將分別基于Simplec和Couple算法對靜止圓柱繞流進行計算。其中Simplec為分離式求解(經(jīng)典的壓力修正算法),Couple則為耦合式求解,以此來檢驗改進方法對不同求解格式的可靠性。
圖2 計算網(wǎng)格分布Fig.2 Grid distribution in simulation
確定好網(wǎng)格和邊界條件后,就需要對α和β進行取值。雖然反饋原理的本質(zhì)不變,但α和β的取值會影響到反饋的效果。α和β選取過大,會引起計算振蕩,甚至導(dǎo)致顯式計算發(fā)散;選取過小,會使得收斂變慢[9]。正如前文所說,目前并沒有理論的方法確定α和β的最優(yōu)值,這與計算狀態(tài)和求解器本身相關(guān),本文選取的值也是經(jīng)過嘗試確定的。
本文研究了α和β的取值對計算結(jié)果的影響。由于物體靜止,式(5)中的第2項為0?;赟implec算法,對比不同反饋力參數(shù)對計算收斂的影響,結(jié)果如圖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)格上,進行下一次迭代,如此反復(fù)進行,直到計算收斂。
取Simplec計算收斂后流場的流線如圖4所示。由于雷諾數(shù)較小,可以看出圓柱后方形成了兩個對稱、穩(wěn)定的分離渦。表1是計算得到的圓柱繞流特征參數(shù)(阻力系數(shù)CD和無量綱回流區(qū)長度L/D),分別是兩種求解方法與改進的反饋力浸入邊界法結(jié)合(表中用Simplec和Couple表示),并與文獻[26-28]的結(jié)果對比。通過圓柱繞流特征參數(shù)的對比可知,本文計算得到的圓柱阻力系數(shù)和回流區(qū)的長度與文獻符合較好。
再取圓柱表面的壓力系數(shù)Cp與文獻[29]對比,如圖5所示。可以看出,本文改進方法計算的壓力系數(shù)與文獻[29]給出的計算結(jié)果符合很好,表明本文改進的反饋力浸入邊界法是有效的。不同于原始方法,本文改進的方法能用于不含時間項N-S方程的迭代求解,且反饋力系數(shù)取值對收斂結(jié)果沒有大的影響。由于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次,采用相同的計算設(shè)置,對兩種方法的收斂特性進行對比。圖6為阻力系數(shù)的收斂曲線。可以看出,采用原始方法進行隱式計算時收斂較慢,特別是在快接近穩(wěn)定值的時候,最終到500個時間步才完全收斂。而本文改進的方法則很快收斂,阻力在100個時間步后就達到平穩(wěn)。
參考文獻[14],定義CFL數(shù)為
式中:Δx=2.96×10-6m,為圓柱周圍Euler網(wǎng)格尺寸。本文改進方法隱式計算與文獻顯式結(jié)果對比如表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時隱式與顯式結(jié)果對比
為驗證改進方法對運動邊界模擬的可行性,選取靜止流體中往復(fù)振動的圓柱[29-30]進行非定常計算,采用Simplec算法求解不可壓非定常N-S方程,通過與文獻對比來驗證本文方法在運動物體求解上的有效性。
本算例計算域的大小取為40D×40D,圓柱將在水平方向進行往復(fù)運動,平衡位置在原點處,運動方程為
x(t)=-Asin(2πft)
(9)
式中:x為圓柱水平方向位移;A為振動幅值;f為振動頻率。
在本算例中,影響計算結(jié)果的關(guān)鍵參數(shù)為雷諾數(shù)Re和KC(Keulegan-Carpenter)數(shù),分別定義為
(10)
式中:umax為圓柱振動的最大速度,umax=2πfA。本文選取Re=100,KC=5。
參考文獻[8],四周遠場邊界設(shè)置為無剪切的物面條件。由于采用隱式計算,一個周期只取360個時間步,每個時間步內(nèi)迭代次數(shù)為20次。反饋力參數(shù)取值見表3。
對于圓柱表面運動速度的計算,可直接對式(9) 求導(dǎo)獲得,然后將該速度代入式(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中渦量的正負分別用虛線和實線表示,當(dāng)圓柱向遠離平衡位置方向運動時,其后會形成兩個強度相同、方向相反的對稱渦;而當(dāng)圓柱向平衡位置返回時,其將分開之前形成的這對渦,并改變渦的旋轉(zhuǎn)方向,這種現(xiàn)象與文獻[30]一致。
圖7 不同相位下的渦量/壓力等值線Fig.7 Vorticity/pressure contours at different phase angles
取改進方法計算的圓柱一個周期的阻力曲線與文獻進行對比,如圖9所示。可以看出,阻力的最大、最小值并不是在特殊相位(t=0.25nT,n為正整數(shù),T為振動周期)時取得,這與圓柱運動過程中流體的慣性有關(guān),表現(xiàn)出了遲滯特性,計算得到的振動圓柱阻力的變化與文獻[30]符合較好。
圖8 運動邊界隱式計算效率對比Fig.8 Comparison of efficiency for implicit iterations at moving boundary
圖9 阻力系數(shù)變化曲線Fig.9 Time history of drag coefficient
以上結(jié)果表明本文改進的反饋力浸入邊界法可與隱式求解方法結(jié)合,對運動物體進行迭代計算,且收斂效果優(yōu)于原始方法。
振動圓柱只有一個方向的運動,橢圓翼的運動要更加復(fù)雜,其不僅有平動,還有繞自身質(zhì)心的轉(zhuǎn)動,可用其來研究昆蟲翅膀運動時非定常流動機理[31]。本節(jié)將選取靜止流體中的橢圓翼來驗證本文改進方法對物體復(fù)雜運動的模擬能力。
橢圓翼運動如圖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。由于運動比較復(fù)雜,物面取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所示??梢钥闯?,由于橢圓翼流場非定常特性明顯,其氣動力變化并不是簡單的正弦曲線。另外,本文計算的橢圓翼的升力阻力變化與文獻[14]基本一致,且相比文獻[14],計算沒有明顯的數(shù)值振蕩,表明可以利用本文改進的方法進行隱式計算,對復(fù)雜的非定常運動進行模擬。
圖11 一個周期內(nèi)不同時刻的渦量圖Fig.11 Vorticity fields during one period
圖12 撲動過程中的氣動力曲線Fig.12 Curves of aerodynamic coefficients during flapping
為進一步驗證本文方法對三維物體模擬的可靠性,選取靜止圓球進行定常計算,求解不含時間項的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所示。可以看出,在Re=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所示??梢钥闯觯疚挠嬎憬Y(jié)果與文獻結(jié)果符合較好,在Re=250時,圓球產(chǎn)生了側(cè)力(表中用Cc表示)。表明本文改進的反饋力浸入邊界方法適合三維物體的計算。
表4 圓球氣動特性對比
本文對Goldstein的反饋力浸入邊界法進行了研究和改進,將原始方法中含有的對速度誤差的時間積分轉(zhuǎn)化為迭代過程中速度誤差的求和,利用多個定常與非定常算例驗證了改進方法的有效性。得到的主要結(jié)論有:
1) 可直接基于迭代次數(shù)對反饋力源項進行計算。與原始方法相比,由于不含時間相關(guān)參數(shù),改進的反饋力浸入邊界法可與求解定常N-S方程的松弛迭代方法結(jié)合,用于定常流動計算。
2) 隱式時間推進可以避免反饋力浸入邊界法顯式求解時嚴(yán)格的時間步長限制。本文改進的反饋力浸入邊界法適合于隱式時間推進,可實現(xiàn)對非定常N-S方程的求解,且收斂特性要優(yōu)于原始方法。
3) 改進的方法與原始方法均是基于速度誤差反饋的思想求解力源項,但本文改進的反饋力浸入邊界法有更廣泛的適用性。