基于格子Boltzmann 的混合過程建模與流致振動特性

尹亞星，王彤，王承彥，張彥康，徐仕承，譚大鵬

(浙江工業大學 機械工程學院，浙江 杭州 310014)

工程流體靜態混合是高壓、高速勢場驅動下的流體介質流經固定機械構件而產生相互交融的過程.靜態混合過程伴隨強剪切、回流反沖、壁面沖擊等物理現象,由流固耦合產生的流致振動較顯著.因此,研究靜態混合過程所涉及的流固耦合振動特性對于航空航天、綠色石化、核電、新能源等工業領域具有重要意義[1-3].在化工強化過程中需要多種高效混合設備,靜態混合器作為其中一種,其共性特征如下：流體介質為工程多相流,固體邊界為薄壁圓柱殼結構；應用過程中流體經過內部混合葉片或孔板后不斷分流、拉伸、碰撞,具有高度非線性動力學特征；流體沖擊壁面,迫使管殼產生振動及噪聲.在絕大多數情況下,這種振動和噪聲是有害的[4].因此,開展流體激勵下靜態混合器流固耦合振動響應特性,探究其流場演化及振動傳播機理,總結其振動響應規律,降低靜態混合過程壁面振動和保證混合安全,在理論和實踐上有著重大研究意義.

在內部流體沖擊作用下靜態混合器中會形成復雜的流致振動,其特征參數與靜態混合管道幾何尺寸、內置元件、流場分布等因素有關[5-6].目前尚未有成熟的理論模型對靜態混合過程流致振動規律進行定量分析,只能在某種假設下或結合實驗研究,得到有限的流場特性[7-9].Haddadi 等[10]采用雷諾平均方程(k-ω)湍流模型對靜態混合器中流場進行研究,最終得出的結果與實驗數據較一致.Murasiewicz 等[11]采用大渦模擬方法得到Kenics混合器3 個分量(軸向、切向和徑向)上速度波動的時間演變,得到靜態混合器內部的穩定速度.孟輝波等[12]通過實驗研究靜態混合器內瞬態壁壓波動特性,發現管壁各點均存在一定程度的壓力波動.Nikolic 等[13]對管道的動力特性和穩定性問題進行研究,得出當內部流體為恒定流時,在較高流速下管道不僅會發生彎曲而且還會發生顫振現象的結論.Leclaire 等[14]采用格子Boltzmann 方法（lattice Boltzmann method，LBM）模擬層流和過渡流中液-液系統的壓降和流速之間的關系,最終仿真結果與實驗結果較吻合.此外權登輝等[15-16]對Kensic 靜態混合器和LSM 靜態混合器進行內部流場的研究,得到混合流場的傳熱和混合性能演變規律.

從上述文獻可以推斷,當前對靜態混合的研究主要集中在兩相動態建模、穩定流場、壁面壓力波動等方面.針對流固耦合界面求解、流致振動特性問題的研究較少.靜態混合內部可以認為是柔性體多個速度分量與剛體產生強烈流固耦合的過程,導致流致振動具有高度非線性特征,這無疑增加了流固耦合求解精度及穩定性.因此需要靜態混合過程流固耦合建模與求解方法,以開展靜態混合器內部流固耦合振動現象的研究.

針對以上問題,本研究提出基于LBM 的多速度分量流固耦合模型,結合大渦模擬（large eddy simulation,LES）湍流模型和流固弱耦合求解方法.基于上述模型,建立靜態混合過程流固耦合動力學模型,并討論內流沖擊與壁面振型的關聯性,揭示靜態混合過程中流體與壁面的相互作用機理.

1 靜態混合數學模型

為了描述高度非線性特征下的靜態混合剛柔接觸過程，須詳細考慮各種速度分量.因此,對于LBM 的流固耦合模型,采用多速度分量的格子模型,并通過完全離散的玻爾茲曼輸運方程進行求解.對于流固耦合振動過程的求解，隨著速度分量的增加,會面臨求解精度及穩定性的挑戰.針對此問題,給出流場-結構場弱耦合求解策略.

1.1 基于LBM-LES 的靜態混合數學模型

靜態混合過程具有復雜邊界條件,在工作過程中具有強剪切、回流反沖以及強烈的壁面碰撞效應,這種現象可以視為多個柔性體與剛體碰撞進而產生強湍流促進混合.

基于介觀的格子-玻爾茲曼方法(LBM)是把時間和空間完全離散,然后將流體離散為若干粒子,將流場劃分為格子,讓這些粒子的分布函數沿網格線運動,并在網格點上根據一定的規則相互碰撞.分布函數的演化在宏觀上反映了流體的運動規律,流場的密度、速度、壓力等力學參數都可以由分布函數的計算得到[17].本研究提出的基于LBM 的靜態混合流固耦合動力學模型,采用的格子模型為D3Q27 速度模型(D 表示維數,Q 表示粒子運動方向的總數),如圖1 所示.該格子模型每個離散節點比其他格子模型具有更高的自由度,更加適用于流體與壁面碰撞過程,能獲得更加精確的壁面振動特性.

圖1 D3Q27 速度模型Fig.1 Velocity model of D3Q27

速度離散后的玻爾茲曼輸運方程如下：

式中：fα=fα(r,t) 表示t時刻r點處α 方向上的粒子分布函數,eα為相應的離散速度,Ωα為碰撞算子,Fα為外力F在α 方向的投影.

為了求解式(1),須對其時間及空間離散,得到完全離散化的玻爾茲曼方程,并對其濾波得到如下公式：

式中：G(r,r′) 為濾波函數.

采用多松弛時間算法[18],則碰撞因子 Ωα可以寫為

相關研究表明,LES 湍流模型能比其他模型更好地描述靜態混合器內部流場[7].根據LES 方法,引入附加黏性系數,即亞格子尺度渦黏系數 υT.為了模擬亞格子湍流,選用wall-adapting localeddy（WALE）渦黏模型,可以更好地反映近壁面區域渦黏系數變化規律,并適合復雜湍流模擬,求解形式如下：

由于靜態混合器的復雜邊界條件,所采用的多松弛時間計算方法的散射算子在中心動量空間實現并獨立松弛,然后執行逆變換以恢復概率分布函數.通過局部宏觀速度離散粒子速度,提高給定速度集的伽利略不變性和數值穩定性.WALE渦黏模型提供了在葉片邊界處局部渦流黏度,可提高近壁面碰撞效應從而促進流固耦合求解精度.

LBM 方法中空間是格子化的,時間也是離散化的,因此,為了模擬真實的物理過程,須考慮量綱的變化問題[20].假設格子上長度的量綱為L,時間量綱為T,質量量綱為G,加上“′”表示對應格子參數,不加則表示對應的宏觀參數.以粒子直徑、流體的運動學黏滯系數以及流體的密度作為基準量.設粒子直徑為D,流體的運動學黏滯系數為ν,流體密度為 ρ,則則宏觀參數到格子參數的變換關系如下：時間t=t′T,長度l=l′L,速度v=壓強P=

1.2 結構域數學模型

采用粒子法與網格法相結合的流固耦合方法,耦合策略采用分區耦合即流體域與結構域采用不同的控制方程.控制結構單元運動變形的動力學方程如下：

式中：M、C、K分別為結構的質量矩陣、瑞利阻尼矩陣和結構剛度矩陣；F(t)為施加在結構上的時域變化外力；y為結構單元節點的位移矢量；α1和α2為系數,與結構的固有頻率及阻尼比相關.

為了得到結構單元節點的位移信息,Newmark-β方法[21]被廣泛應用于式(10)的求解.在t=t+Δt時刻,由于式(10)同時包含未知矢量尚須補充2 組方程才能使式(10)封閉.對結構單元節點的速度和位移進行泰勒展開,以實現對式(10)的求解：

式中：β 和γ 為權重參數,與計算靜態混合器模型穩定性相關,本研究分別取β＝0.25,γ＝0.5.

1.3 LBM 流固耦合方法

基于上述數學模型,提出流固耦合弱耦合求解方法,如圖2 所示.該方法計算流程主要包括初始化邊界條件、LBM 控制方程和LES 方程求解、結構場控制方程求解、流固耦合面分布函數重構等.在分區求解過程中選擇大小合適的計算域,劃分流場網格,進行流固耦合計算初始化,對流場和結構賦初值；求解LBM 流場控制方程和LES 湍流模型,得到流體的壓力和速度.采集流固耦合面上流體作用力,將流體作用力施加在結構上,求解式(10)并提取結構運動的位移、速度數據；根據結構位移更新結構在流場中的位置,根據Yu 等[22]提出的插值反彈格式方法,通過鄰近點的流體粒子分布函數進行插值計算,最終保證流固耦合面上速度、壓力、變形等物理量的守恒.上述求解過程不斷重復,最終輸出結構變形結果和流場計算結果.

圖2 LBM 流固耦合計算流程圖Fig.2 Flowchart of fluid-structure interaction calculation for LBM

2 靜態混合動力學模型

2.1 靜態混合器幾何模型

以Ross LPD 型號靜態混合器為對象實例進行研究.為了簡化建模過程及降低對計算機硬件的要求,將幾何模型進行合理的簡化,簡化后的幾何模型如圖3 所示.該幾何模型只保留了靜態混合器自身主要構件,包含2 個入口和1 個出口與大氣相通,以及心軸與混合葉片.以入口處為原點,建立柱坐標系,分別用x軸表示軸向,r軸表示徑向,θ 軸表示周向描述整個模型.根據相同比例對數值計算進行建模,靜態混合器和流體的幾何參數如表1 所示.表中，L0為長度,R為內徑,h為壁厚,d為心軸直徑,ρ 為靜態混合器密度,E為彈性模量,μ為泊松比,ρw為水的密度,μw為黏度.內部一系列交錯的半橢圓板周期交叉焊接在圓柱心軸上[23],3 種螺旋葉片夾角 αb分別取80°、90°、100°,3 組葉片左旋右旋交錯排布.

表1 靜態混合器和流體物理參數Tab.1 Physics parameters of static mixer and fluid

圖3 靜態混合空間結構示意圖Fig.3 Diagram of static mixed space structure

2.2 邊界條件和初始條件

在流固耦合求解中,流體通道的主入口（入口1）直徑為70 mm,側入口（入口2）直徑為10 mm.流體入口設置為速度入口,出口均為自由出口,入口2 速度恒定為4 m/s,由于主次入口直徑比為7∶1,相同入口速度情況下的體積流量比為49∶1,因此主要考慮入口1 的速度變化對靜態混合器壁面的振動影響.側面采用零摩擦的自由滑移墻面邊界,靜態混合器的邊界條件設置為兩端固支,流固耦合面選取為靜態混合器內壁面以及葉片表面.為了觀察靜態混合器變形隨時間變化的過程,采用非定常求解方法進行計算.計算時間步長為0.00 002 s,時間步為10 000 步,總時間為0.2 s.由于總仿真時間較長,在0.2 s 后,靜態混合器內部流體已達到穩定狀態.

2.3 靜態混合流固耦合模型驗證

基于上述幾何尺寸模型,建立靜態混合過程流固耦合動力學模型,如圖4 所示.格子玻爾茲曼模型使用的格子劃分策略采用八叉樹晶格結構劃分流場格子,然后生成粒子,固體域采用非結構化網格,在流固耦合面處進行節點與格子匹配,保證固體域網格和流體域最大格子尺寸相同,目的是使流體域粒子各速度分量與固體域網格節點之間準確地相互傳遞數據,保證數值模擬的準確性.

圖4 流固耦合動力學模型Fig.4 Hydrodynamic model of fluid structure coupling

考慮到格子大小和網格尺寸會直接影響結果的精度,使用不同的格子和網格密度進行無關性驗證.模擬pout=0 Pa,入口速度vin=4 m/s 的工況下靜態混合過程流固耦合,使用5 種網格尺寸(編號1～5),考察格子尺寸對數值結果的影響.對0.2 s 時x=0.275 m,θ=0°處外壁面振幅以及總仿真時間進行對比,綜合考慮流體域粒子數量和固體域網格數量與振幅和仿真時間的變化趨勢,如表2 所示.表中,M為流體粒子數,N為固體域網格數,Ar為混合器壁面徑向振幅.隨著粒子數量與網格密度增加,幅值逐漸相差較小.編號3、4、5 這3 種網格尺寸,振幅呈現基本一致,表示已經滿足了網格無關性要求,保證了數值計算的準確性和可重復性.本研究最終采用方案4 進行流體域與結構域離散化處理.

表2 網格無關性驗證Tab.2 Grid independence verification

為了驗證該方法的正確性,利用Ansys-Fluent軟件計算相同物理參數下的算例,并將軸向位移時間曲線與LBM 方法計算得到的結果對比.在Fluent 中流體同樣采用LES 模型,在流固耦合的每一步內,先求解流體控制方程,得到速度場、壓力場以及作用在壁面的作用力,然后求解固體域方程,將靜態混合器流固耦合面變形通過Fluent動網格宏進行處理,實現動邊界對流場的反饋作用,從而進行雙向流固耦合.

LBM 方法和Fluent 流固耦合分析后的變形對比，如圖5 所示.圖中,L為軸向距離,dr為混合器壁面徑向位移.通過2 種模擬方法得到0.2 s 時刻在θ=0°沿x軸方向上的位移響應曲線,2 種計算得到的結果大致吻合.由于2 種方法是截然不同的理論公式,在網格劃分過程中，傳統方法是流體網格與固體網格間相互傳遞數據,LBM 方法采用的則是粒子與固體網格之間相互傳遞數據,本研究通過LBM 方法采用27 個速度分量在流固耦合面處傳遞更多方向上的速度、壓力數據.因此采用LBM 流固耦合建模方法,可以更加真實精確地模擬靜態混合器內部的流體沖擊.

圖5 靜態混合器壁面位移曲線對比Fig.5 Comparison of wall displacement curve of static mixer

3 數值計算結果與討論

靜態混合器物理空間內部葉片的存在增加了湍流場的無序性和非線性,為了獲得其中流固耦合特性,須研究靜態混合器整體外壁面位移響應和不同流速、不同葉片夾角條件下的靜態混合器局部外壁面位移響應和幅頻響應.

3.1 靜態混合器流致振動響應分析

在靜態混合過程中,由于內部流體瞬態沖擊,流致振動規律具有高度非線性特征.為了更好地觀察管道內流體對管道振動的影響,在靜態混合器模型壁面上添加監測點,沿管體x軸方向,平均每隔0.025 m 設置1 個監測點,共設置23 個監測點.由于模型具有對稱性,在半周向進行監測點的選取,周向選取相位角分別為θ=0°,30°,60°,···,180°的a～g共7 個監測點.如圖6 所示為靜態混合器在內部流體激勵作用下的形變云圖.可以看出，中間部位變形最大.

圖6 周向監測點示意圖及內部流體激勵作用下靜態混合器形變云圖Fig.6 Schematic diagram of circumferential monitoring points and deformation of static mixer under internal fluid excitation

如圖7 所示為靜態混合器沿x軸方向的振幅均方根.圖中,RA為振幅均方根.可以看出，靜態混合器呈現五階振型模態,中部振幅均方根最大.由于靜態混合器和管道類似屬于細長結構,會在重力作用下發生較為明顯的形變.靜態混合器內介質流經葉片的管道時,由于葉片的阻礙,流體速度發生改變，軸向速度轉變為徑向速度和切向速度,這些加速的力反過來作用在管道上,引起管道的附加振動.

圖7 靜態混合器沿x 軸方向振幅均方根Fig.7 RMS amplitude of static mixer along x axial direction

在數值仿真結果的處理中,為了便于系統比較,選取混合器振幅最大位置(即x=0.275 m)處的振動情況進行研究.對圓周方向上的7 個監測點提取徑向、周向和軸向位移響應,如圖8 所示.圖中，dθ、dx分別為周向、軸向位移.可以看出,在初始階段,位移響應趨勢不太穩定,表明流體在剛沖入時對壁面產生壓力,壁面也會施加給流體一個回彈力,短時間內這種反復的耦合響應導致初始位移失穩,在一定時間后達到穩定的力場,進而形成一個穩定的振動響應.

圖8 圓周方向各監測點的位移響應Fig.8 Displacement response of each monitoring point in circumferential direction

由圖8(b)、(c)可知,當振動響應穩定后,周向和軸向位移響應曲線幾乎完全重合.由圖8(a)的徑向位移響應對比可以看出,a比g的變形大,b比f的變形大,c比e的變形大,表示對稱的兩側所受的力大小不等,主要原因是流體受到重力的影響,混合器壁面上側的流體沖擊力與流體重力相互抵消,混合器壁面下側流體沖擊力與流體重力疊加.同樣，g、a變形相差值大于b、f和c、e兩組.另外,對比3 個方向上的振動位移響應曲線可知,靜態混合器的徑向位移響應最大,其次是周向,軸向的位移響應最小,而且徑向位移的峰值更為密集.這表示靜態混合器內部流體沖擊導致的附加振動對徑向位移響應的影響最大,同時驗證了靜態混合器內部流體速度轉化會導致附加振動且使得靜態混合器壁面振動現象更加復雜.

3.2 不同入口速度對靜態混合器振動特性的影響

在靜態混合器實際應用過程中,常采用不同的入口速度提高混合器的混合效率,而不同流速變化將產生不同的振動響應.為了研究不同流速對靜態混合器壁面沖擊的影響，選取t=0.2 s 的內部流場速度分布以及x=0.275 m,θ=0°處靜態混合器外壁面監測點的壁面位移響應變化情況.如圖9(a)～(c)所示為3 種vin下流場速度分布云圖.圖中,v為流場速度.可以看出,流體與第1 組葉片接觸,流動方向突然改變,形成高速旋轉的流體；之后繼續與第2 組和第3 組快速接觸,必然引起局部湍流渦團的耗散,增加湍流狀態的非線性,這也是引起靜態混合器壁面附加振動的主要原因.同時,對比3 種工況下云圖.可以發現，速度越大,在葉片表面附近速度突變程度越大,且主要發生在第1 組葉片處,表明第1 組葉片比其他2 組更容易受磨損,損壞概率更高.

圖9 不同入口速度下的流場速度分布云圖Fig.9 Flow field velocity distribution nephogram of different inlet velocity

如圖10 所示為不同流速流體激勵下靜態混合器壁面振動響應曲線.圖中,Ar為徑向振幅,f為頻率.由圖10(a)可以看出，入口流速對混合器的位移響應和振動頻率有較大影響,三者振型都不太規則,平均振幅差別顯著,隨著流體速度增大,流體對靜態混合器沖擊振動的影響愈發顯著,流體速度越大,平均振幅越大.

圖10 不同流速流體激勵下靜態混合器壁面振動響應曲線Fig.10 Vibration response curves of static mixer wall under fluid excitation at different flow rates

由圖10(b)可以看出流體入口速度為3、5、7 m/s的工況,三者振動情況較為相似,靜態混合器的振動普遍存在多模態共存的現象,振動情況較分散,這是因為靜態混合器內部不規則漩渦會誘發隨機振動.三者工況的主振動頻率分別為114.4、19.9、14.9 Hz,可見隨著內部流體速度的變化,流速較大的工況使振動更加趨于低頻化,促使靜態混合器在低頻段更容易發生共振現象.對比不同速度下的振幅可知,隨著速度增加,最大振幅差距顯著,7 m/s 時最大振幅為3 m/s 時最大振幅的4.86 倍,5 m/s 時最大振幅為3 m/s 時最大振幅的3.38 倍,表明速度越大,最大振幅越大,隨著速度增大,最大振幅增加的比值有減小的趨勢.靜態混合器的振動過程通常會由多個頻率共同影響,但其振動主要受某一階模態控制,其他模態對靜態混合器的振動影響較小.

3.3 不同混合原件對靜態混合器的振動特性的影響

在靜態混合器設計過程中,混合原件扮演著重要的角色,不同的混合原件對流場的剪切作用不同從而產生不同的振動誘導因素,造成不同程度的疲勞損傷.在實際情況中,主入口速度越大,壓力降越大,能量損失越多.考慮到靜態混合器實際情況,采用偏小的主入口速度,主要研究葉片對振動特性的影響.選取t=0.2 s 時刻入口1 速度4 m/s 內流場速度分布以及x=0.275 m,θ=0°處靜態混合器外壁面監測點的壁面位移響應的變化情況.不同葉片夾角下流場速度分布如圖11 所示.可以發現,隨著葉片夾角的減小,第1 組葉片表面處流速突變越來越明顯,而第3 組葉片處流速突變減弱.這是因為第1 組葉片夾角最先改變流場的徑向剪切流動能量,誘導軸向速度轉化為徑向速度,而軸向速度的減小勢必導致第3 組葉片處軸向速度向徑向速度的轉化減少.管壁附近速度遠遠大于中心部位;葉片夾角越小,混合前、后葉片流速相差越大，葉片夾角較小時混合后的速度波動比混合前起伏更大且穩定時流速相對較高,表明葉片夾角越小,對流體的引流剪切作用越強,可能導致混合器的穩定性下降.

圖11 不同夾角混合原件的流場速度分布云圖Fig.11 Velocity distribution diagram of mixed elements with different angles

如圖12 所示為不同葉片夾角對混合器壁面位移響應的影響.可以看出，靜態混合器中的位移響應有延遲,初始階段周向位移反應較迅速,徑向和軸向則需要一定時間達到穩態;在達到穩態時,角度越小,徑向位移響應越大；對于周向位移,90°時周向位移最小,100°周向位移最大；軸向位移隨著角度變化，夾角越小，位移響應越大.對比3 個方向的位移響應,徑向位移響應最大,其次是周向,軸向位移響應最小.

圖12 不同葉片對混合器壁面位移響應的影響Fig.12 Effect of different blades on wall displacement response of mixer

如圖13 所示為不同葉片對混合器壁面幅頻響應的影響.圖中，Aθ、Ax分別為周向、軸向振幅.可以看出，當改變葉片的夾角時,幅頻響應的低頻段出現幅值較大的峰值,周向方向葉片夾角越小峰值越大,而徑向方向和軸向方向,在葉片夾角為90°時峰值偏小,增大或減小夾角低頻段峰值都會增大,這種現象是由于結構改變造成剛度變化,造成的影響主要表現在幅頻曲線的低頻段.不同葉片夾角其主振動頻率未發生改變,徑向、周向、軸向主振動頻率分別為4.97、4.97、44.77 Hz.原因是徑向主頻率與軸向主頻率頻率受葉片的影響,流體速度經葉片引導,軸向速度轉變為徑向速度和周向速度,從而引發周向和徑向附加振動,振動頻率偏小,整體未出現低頻化現象,表明改變葉片夾角不會促使靜態混合器在低頻段發生共振現象.

圖13 不同葉片對混合器壁面幅頻響應影響Fig.13 Effect of different blades on wall amplitude-frequency response of mixer

4 結論

(1) 對壁面的位移響應分析表明,在混合過程中靜態混合器會在流體重力的作用下呈現五階振型,內部葉片的存在將流場軸向速度轉化為切向速度引發附加振動,使得徑向位移響應最大,而且徑向位移峰值更為密集.

(2) 在葉片相同的工況,增加入口流體速度,葉片附近流場速度突變會越明顯,且主要發生在第1 組葉片處；當入口流速較大時,靜態混合器振動頻率趨向于低頻化,振幅會增大,即加劇了振動響應.

(3) 在流速相同的工況下,葉片夾角越小,對流體的引流剪切作用越強,管壁處速度顯著增大,振動幅值增大.周向位移在葉片夾角為90°時位移響應最小.不管混合葉片怎么改變,頻率都不發生改變,葉片夾角對振動頻率無影響.

基于LBM 流固耦合模型的靜態混合系統具有高度的復雜性,本研究在此方面進行了有益的嘗試.在理論方面可以為基于LBM 的流體-結構耦合求解方法和多速度分量描述系統的建模提供參考；在技術方面,可以為工業過程監測提供更多的技術解決方案,具有較好的工程應用前景.

本研究僅通過數值計算模擬單向流條件下靜態混合器壁面振動響應特點,今后工作重點將放在LBM 多相流流固耦合動力學模型建立上,以實現多種流體的混合計算.除此之外,本研究偏向于理論計算,后續將建立靜態混合實驗平臺對數值計算模型進行驗證.