基于力源擬化邊界思想的工程結構繞流模擬

基于力源擬化邊界思想的工程結構繞流模擬

0流-固耦合工程抗風空氣繞流設計驗證是計算流量值的重要目標，也是其應用價值的體現(xiàn)。但在實際工程中,構件通常具有任意幾何外形并伴有流-固耦合現(xiàn)象。要得到更為精確的CFD數值模擬結果,精細地構造結構模型以及合理地解決動網格求解是必不可少的前提條件。目前數值模擬方法多采用貼體網格構造物體邊界,但受制于計算效率等多方面因素,并不能很好地滿足對結構細部構件的模擬,且對動態(tài)問題求解困難。這也是CFD發(fā)展至今在工程應用方面所遇到的主要瓶頸。因此,發(fā)展出能夠涵蓋任意幾何邊界并適用于流-固耦合求解的CFD方法已成為工程結構抗風研究領域亟待解決的問題。在此過程中,高效性和實用性是2個重要的評判標準。相比于其他CFD方法,浸入邊界法模型概念簡單,最大優(yōu)勢在于能夠很好地適用于流-固耦合問題的求解,通過構造隱含邊界避開動網格所帶來的計算效率的降低。雖然其尚未拓展至結構工程抗風領域,但其已經廣泛應用于醫(yī)學、仿生學等領域中流-固耦合現(xiàn)象的研究,已形成一個較為成熟的體系［１－１１］。鑒于此,將浸入邊界法應用于工程抗風領域具有很強的可行性。然而目前開展此類方法研究的大部分文獻多集中于單一邊界:圓柱,進行繞流模擬分析,還未全面地展開對其他幾何邊界物體繞流模擬研究,且缺少對方法的改進以提高實用性［１２－１５］。因此,全面系統(tǒng)地開展任意邊界物體繞流的綜合研究,并進一步研究分析算法本身作用機理才是實現(xiàn)其工程應用最根本的出發(fā)點。為此,本文選定適用于計算剛性邊界物體繞流的浸入邊界算法(VirtualBoundaryMethod),結合高效求解方法———分步法(FractionalStepMeth-od)及有效的邊界處理方式———交錯網格,重點介紹插值方法的運用,通過邊界點與周圍網格點的數據傳遞,構造物體邊界并分析其形成機理;系統(tǒng)的展開對靜止方柱、圓柱、有風攻角的靜止方柱繞流及二維橋梁斷面的繞流數值模擬,展示本文對任意邊界物體繞流模擬的全面性和可行性。1流體速度界面控制方程早期浸入邊界法模型研究對象均為彈性邊界,應用于剛性邊界計算時,其邊界實現(xiàn)需設定較大的彈簧常數值,這種情況導致流動方程剛性增大,計算效率降低。針對于此,Goldstein將浸入邊界法同反饋循環(huán)控制力源結合,發(fā)展出適用于剛性邊界的VirtualBoundary方法［１４，１６－１７］。這種方法主要特點即在于反饋,對邊界上流體施加反饋力源以使流體速度不斷趨近于物體邊界實際速度,以構造邊界。VirtualBoundary方法基本控制方程組可表達為式中:U為速度;t為時間;ρ為流體密度;P為壓力,v為流體運動粘度;F為力源項,施加在物體邊界點上。通過回饋循環(huán)控制物體剛性邊界的形成為式中:F(XＳ,t)為物體邊界點的坐標;αｆ、βｆ為反饋常數項,其取值直接關系到邊界擬化效率;V(XＳ,t)為物體邊界速度;u(XＳ,t)為邊界點流體速度;u(XＳ,t)-V(XＳ,t)體現(xiàn)邊界處流體速度與邊界速度的差異。表達式正是基于此差異,通過反饋機制回饋出力源使兩者速度接近,以滿足無滑移邊界條件。采用龍哥庫塔或亞當斯-巴什福斯公式差分方程組時,為滿足計算穩(wěn)定性,時間步長限制公式如式(4)所示為式中:k為差分方法中一階常數項?？梢钥闯?過大的反饋系數值αｆ、βｆ,會因計算穩(wěn)定性要求導致時間步長縮小,計算效率降低。此類現(xiàn)象普遍存在于VirtualBoundary方法中。2基于雙線性插值的拉格朗日誤差獲取如式(1)~式(3)所述,在VirtualBoundary算法中邊界的構造依賴于反饋函數,而形成有效反饋函數的關鍵前提是獲得實時的邊界點流體速度,其一般通過邊界識別和數值傳遞2個過程來獲取。模擬首先在計算區(qū)域中辨別出用以代表物體邊界的邊界點,稱為拉格朗日點,其空間位置取決于物體邊界的運動狀態(tài),對于靜止物體,拉格朗日點空間位置不隨時間而變化。其他計算網格點稱之為歐拉點,歐拉點空間位置不隨時間變化。其次實際情況中,具有復雜幾何外形的物體邊界任意穿過計算格子以及流-固耦合產生的邊界運動,導致拉格朗日點不能與歐拉點重合,拉格朗日點流動參數值的獲取需要通過其周圍歐拉點上的數值插值———即構造傳遞函數來完成,以進一步形成作用于邊界點上的反饋力源。浸入邊界法中的傳遞函數最初以Delta函數為原型［１８］。但傳遞函數是浸入邊界法中作為邊界點同周圍網格點進行數據傳遞的橋梁,必須保證便于參與方程離散求解和繞流模擬計算的精度,其不能僅以Delta函數最初的簡單形式應用,必須得到更細化的構造。針對于此,Peskin首先定義浸入邊界法中傳遞函數的廣義表達式為［１８］式中:rｘ=(xｋ-xｉ)/h;rｙ=(yk-yｉ)/h;rｚ=(zｋ-zｉ)/h;h為網格邊距大小;xｋ、yｋ、zｋ為歐拉點即網格點坐標;xｉ、yｉ、zｉ為邊界點即拉格朗日點坐標。構造式(5)中插值函數f(r)的表達式為此種插值函數雖然能夠保證很好的插值精度,但形式較為復雜,編程實現(xiàn)較為繁瑣,不利于擴展至工程應用研究。實際上,插值函數f(r)只須滿足插值函數準則式(8)~式(10),也可以轉化為其他更為簡潔的構造形式。Goldstein等曾選擇譜函數作為其實現(xiàn)數據間傳遞的途徑,但結果表明此方法易造成邊界數值震蕩［１４］。本文借助于Saiki等提出的精度較好且更易于實際應用的雙線性插值方法(BilinearInterpola-tionMethod)構造插值函數,實現(xiàn)本次模擬［１５］。圖1中P為邊界點,(i,j),(i+1,j),(i+1,j+1),(i,j+1)為其周圍4個網格點。邊界點P上流動參數值U(xｓ)通過雙線性插值由周圍網格點上參數值Uｉ，ｊ插值得出式中:xｉ、yｊ為歐拉點即網格點坐標;xｓ、yｓ為邊界點即拉格朗日點坐標。函數d(x)的表達式為3邊界流體速度的彈簧體系對VirtualBoundary算法模型方程式(1)簡化得出式中:q=u-v對應式(4)所表述的速度差異;方程右邊第一項包含反饋系數αｆ的反饋積分,可以形象地看成連接微觀流體單元同物體表面的彈簧力,以使得該流體單元始終依附在擬化邊界上;第二項βｆq等效于周圍流體單元對彈簧運動產生的阻尼,控制響應的振蕩頻率。式(13)進一步揭示出在VirtualBoundary算法模型中,反饋函數實質上構成邊界流體速度的彈簧體系,以滿足無滑移壁面條件。選定方柱在VirtualBoundary算法模型下,展開邊界擬化分析,觀察不同反饋系數αｆ、βｆ組合下邊界形成效率,以邊界點速度累積誤差L作為無滑移邊界實現(xiàn)評判標準,進而選擇合適的反饋系數展開后續(xù)數值模擬為式中:nｂ為邊界點個數;uｉ為邊界點流體速度。下頁圖2為不同常系數組合下,邊界點處速度累積誤差L的分析結果。圖2中前2組常系數組合圖2(a),圖2(b)邊界點速度累積誤差收斂為0大約需9~27s,滿足壁面條件所需時間過長,影響程序模擬效率,不利于對動態(tài)流場的追蹤;增大常系數取值,如圖2(c)所示,誤差收斂所需時間銳減至0.7s;繼續(xù)增大常系數取值圖2(d),誤差收斂所需時間約0.6s。比較圖2(c)與圖2(d),雖然繼續(xù)增大反饋系數取值(如圖2(d)所示)仍能提高邊界實現(xiàn)效率,但增幅較小,而且較大的反饋系數組合導致方程剛性過大,計算效率降低,因此圖2(c)所選取的反饋系數組合更適合展開數值模擬。進一步分析采用圖2(c)所示反饋系數組合所形成的方柱內部流體,如圖3所示。圖3中內部流場速度平均值在整個計算時域內其數值波動處于很低的幅度,且其數值大小恒定在10－３左右,符合物體邊界擬化形成后,內部流體(In-ternalFlow)與外部流場完全隔絕,相互獨立,內部流體處于衰減狀態(tài)的算法理論。4出口邊界條件計算網格采用等距笛卡爾網格,反饋系數組合αf=-160000,βf=-600,流動從左至右,入口邊界條件為:u=1,v=0。計算域上下邊界采用對稱邊界:,v=0,出口邊界條件設置,物體中心位置距離上下邊界H/2,物體距離出口邊界約為15D,以減少邊界效應對計算流域內流體流動的影響。模擬簡單鈍體時,如:方柱,圓柱等,網格間距dX=dY=0.025,Δt=0.00015s,網格總數為800×320,阻塞率B/H設定為0.125。橋梁斷面模擬時,網格間距dX=dY=0.05,Δt=0.00015s,網格總數2500×500,阻塞率B/H設定為0.04。雷諾數Re定義為Re=Du/v,其中D為物體特征長度,u為入流速度,H見圖4。計算區(qū)域和邊界條件如圖4所示。4.1不同雷諾數下阻力系數變化規(guī)律模擬方柱繞流,不同雷諾數方柱繞流形態(tài)如圖5所示,圖5通過速度流線反映不同雷諾數方柱的繞流形態(tài)。圖5中雷諾數Re為10,30時,流動對稱且方柱尾流區(qū)出現(xiàn)流動分離現(xiàn)象,產生上下對稱的穩(wěn)定回流漩渦。隨著雷諾數增大,回流渦長度變大。雷諾數Re增至60時,方柱回流旋渦對稱性消失,方柱后緣流體形成不穩(wěn)定的剪切層,進而形成向下游流動的旋渦。圖6,圖7為不同雷諾數下阻力系數變化規(guī)律及回流區(qū)長度。圖6顯示雷諾數在10~50范圍內,方柱后部回流區(qū)長度呈現(xiàn)線性變化,對應圖5(a),圖5(b)。圖7中在雷諾數較低時,阻力系數Cｄ呈現(xiàn)較高的值,隨著雷諾數升高,阻力系數逐漸降低,在雷諾數為100時降到最低,之后阻力系數Cｄ趨于平穩(wěn)。本文數值模擬中回流渦長度Lｗ及阻力系數Cｄ隨雷諾數變化曲線比較同類文獻[19,20],基本符合。4.2不同雷諾數下圓柱阻力系數變化規(guī)律下頁圖8為圓柱在雷諾數分別為40,50的速度流線圖。如圖8所示,雷諾數為40時,圓柱繞流仍為對稱流動,尾流區(qū)形成對稱分離渦。雷諾數升至50時,圓柱尾流對稱渦破壞,此雷諾數為圓柱繞流的關鍵雷諾數,同文獻結論相同。下頁圖9對應不同雷諾數圓柱阻力系數變化規(guī)律,反映出圓柱阻力系數隨雷諾數增高逐漸降低,Re>100后,阻力系數仍處于減小趨勢。圖10為不同雷諾數下圓柱尾流區(qū)回流區(qū)長度,Re<50時,對稱尾流下的回流區(qū)長度隨雷諾數線性增加;Re>50后,尾流區(qū)對稱性消失,回流區(qū)長度呈現(xiàn)非線性銳減趨勢。圖11為不同雷諾數有風攻角的方柱繞流流線圖,雷諾為40時,有風攻角的方柱尾流對稱渦即消失,應為其迎風面尖銳邊界點所致。在雷諾數為50時,流動在方柱后緣形成較大的速度回旋區(qū)。Re繼續(xù)增加至100,流動在方柱下端出現(xiàn)較小的速度回旋區(qū)。同其他物體繞流相比,有風攻角的方柱在較低雷諾數時尾流對稱性就消失。4.3種攻角下梁完成上述任意邊界鈍體繞流模擬后,繼續(xù)采用本文浸入邊界算法程序,展開Re為2×10３下二維橋梁斷面繞流的初步模擬。下頁圖12為0°攻角下橋梁斷面壓力分布圖及渦度顯示圖,來流在橋梁前緣處產生流動分離,橋梁上表面產生速度回旋區(qū),對應圖12(b),0°攻角下橋梁尾流區(qū)呈現(xiàn)規(guī)律性旋渦脫落。圖13為不同攻角下(θ=0°,3°)流場速度顯示圖。θ=3°時,橋梁斷面尾部速度流場出現(xiàn)擾動,相對于0°攻角,可以觀察到尾部流場出現(xiàn)明顯的流場非對稱性。表1為2種攻角下所得出的氣動力參數,由于橋梁斷面接近于流線體,其阻力系數和升力系數均小于相同雷諾數時方柱所對應數值。結果顯示:斷面攻角增大后,阻力系數Cｄ小幅度增加,升力系數CＬ出現(xiàn)顯著提升,對比2種攻角下的斯托羅哈數Sｔ,可以看出攻角的改變并不影響斷面尾部的旋渦脫落頻率。5數值模擬結果(1)基于VirtualBoundary方法結合較為成熟的雙線性插值編寫出流場模擬程序,分析VirtualBoundary算法模型下物體邊界擬化機理,討論不同反饋系數組合對邊界擬化的效率,以邊界點累積誤差為標準,選擇出合適的反饋系數進行數值模擬并顯示出邊界形成后物體內部流體衰減狀態(tài)。(2)對方柱、圓柱、有風攻角的方柱進行數值模擬,得出不同邊界鈍體尾流對稱分離渦消失所對應的雷諾數,并將模擬所得其他各氣動參數同其他文獻驗證,吻合度均較好,驗證了本算法程序在傳統(tǒng)CFD領域中的可行性和普適性。(3)拓展至