謝 宇, 任金波, 黃煌輝, 張 翔
(福建農(nóng)林大學機電工程學院,福州 350002)
計算流體動力學(computational fluid dynamics, CFD)[1],其技術取得了快速的發(fā)展[2-3],數(shù)值模擬方法的應用日趨廣泛[4-7],所面對的問題越來越復雜,但數(shù)值模擬在不同工程環(huán)境下的表現(xiàn)值不一樣[8],因此數(shù)值模擬的可靠性驗證是判斷數(shù)值模擬可信度的重要方面,計算結(jié)果可靠性分析方法主要有基準解比較法、實驗比較法和網(wǎng)格收斂性分析法三種[9-13],網(wǎng)格收斂性分析法中最流行的是Roache[14]提出通過計算網(wǎng)格收斂指數(shù)(grid convergence index,GCI) 來評估數(shù)值誤差的方法,現(xiàn)階段GCI方法已廣泛應用于數(shù)值計算誤差,Craig等[15]利用GCI對厭氧消化器攪拌模型進行誤差分析,Karimia等[16]分析了GCI 中的關鍵參數(shù)對旋液分離器數(shù)值誤差的影響,鄭秋亞等[17]針對M6機翼的繞流問題,就網(wǎng)格密度對 5 套不同網(wǎng)格的CFD模型進行了GCI分析, Longest等[18]使用GCI分析不同網(wǎng)格樣式對于生物流動系統(tǒng)的流速場等因素的影響,劉厚林等[19]借助GCI對比分析了不同網(wǎng)格類型的離心泵內(nèi)部流場情況,GCI方法在一定程度上提高模擬的準確性并提供理論支持,但仍存在以下不足:① GCI方法在網(wǎng)格以整數(shù)倍形式加密時計算結(jié)果可靠,但是在實際操作中(尤其是多維模型時)往往很難做到,同時若網(wǎng)格加密倍數(shù)過小,解的變化不明顯;②在多維模型中,不同方向網(wǎng)格加密倍數(shù)可能不一樣,須在每個方向先分別計算相應的 GCI,然后再迭加求其和,而在實際操作中很難做到;③截差階數(shù)p較難確定;④漸近范圍的判斷沒有絕對值,網(wǎng)格加密到一定程度時計算結(jié)果會進入到漸近收斂范圍,再對網(wǎng)格進行加密沒有太大的意義[20-21],因此,確定合適的漸進收斂范圍很重要。
基于此,對刀片切泥制漿數(shù)值模擬進行CFD網(wǎng)格無關性的分析,探討在混合網(wǎng)格下網(wǎng)格無關性的網(wǎng)格數(shù)量,在允許誤差范圍內(nèi)借助卡方檢驗驗證網(wǎng)格無關性,再與GCI方法進行對比,最終和試驗結(jié)果來判斷是否合理,探究適用于刀片切泥制漿模擬的網(wǎng)格無關性所對應的最優(yōu)網(wǎng)格數(shù),以分析在非整數(shù)加密以及加密方式不一致情況下的網(wǎng)格收斂分析準確性,以期為今后相關研究提供借鑒和啟發(fā)。
以刀片旋轉(zhuǎn)切割泥漿來培漿的過程為研究對象,設計泥漿槽長寬高為2 000 mm×2 000 mm×1 600 mm,一側(cè)有一臺階,刀片轉(zhuǎn)速為40 rad/s,泥漿被旋轉(zhuǎn)的刀片切割飛濺至一側(cè)臺階上,來模擬現(xiàn)實中水稻田里泥漿被刀片旋轉(zhuǎn)打到育秧盤的工作環(huán)境,其物理模型如圖1所示。
圖1 刀片切泥制漿模型簡圖Fig.1 Model diagram of blade cutting mud
用Pro/Engineer軟件對模型進行建模,使用CFD前處理軟件(the integrated computer engineering and manufacturing code for computational fluid dynamics, ICEMCFD)盡可能對結(jié)構模型采用結(jié)構化網(wǎng)格劃分,網(wǎng)絡模型的劃分方法和網(wǎng)格質(zhì)量對數(shù)值模擬結(jié)果具有極大影響,網(wǎng)格quality是判定網(wǎng)格質(zhì)量的重要方面,使大其數(shù)值接近1,采用混合網(wǎng)格,泥漿槽部分采用結(jié)構網(wǎng)格劃分,刀片部分采用非結(jié)構網(wǎng)格劃分,刀片的邊緣和刀面細化處理,如圖2所示。
圖2 刀片區(qū)域網(wǎng)格Fig.2 Domain of the blade grid
采用商用軟件Fluent16.1,泥漿設置為賓漢姆模型,近壁面處理采用標準壁面函數(shù);采用混合模型(mixture model)對氣液相互作用進行了表征;邊界均設置為Wall;采用滑移網(wǎng)格(sliding mesh)描述刀片旋轉(zhuǎn);使用三維雙精度求解器;壓力-速度耦合方法采用SIMPLE算法及一階迎風格式;湍流模型采用RNGk-ε模型;時間步長取為0.000 1 s,計算時間為1 s,在每個時間步長內(nèi)設置最代次數(shù)為20次,收斂殘差設為0.001。
圖3 泥漿體積分布Fig.3 Domain of the blade grid
現(xiàn)主要探討模擬計算結(jié)果中的泥漿飛濺量、泥漿飛濺速度的網(wǎng)格無關性檢驗。設計9套網(wǎng)格數(shù)量對重要部分(刀片區(qū)域)和非重要部分(泥漿槽流體區(qū)域)進行不同程度的加密處理來進行網(wǎng)格無關性驗證,模擬1 s后,觀察其模擬結(jié)果,圖3為刀軸橫截面處泥漿體積分布圖,圖4所示為刀軸橫截面處泥漿飛濺速度圖。 數(shù)值模擬數(shù)據(jù)如表1所示,泥漿飛濺體積量是飛過距離刀軸一側(cè)0.2 m橫截面的泥漿體積量,如圖5所示;平均側(cè)面飛濺速度是距離刀軸一側(cè)0.2 m處泥漿的平均飛濺速度,如圖6所示。由圖5、圖6可知,網(wǎng)格數(shù)量對仿真模擬結(jié)果具有較大影響,當網(wǎng)格數(shù)量接近300×104時,模擬結(jié)果變化差異不大。
圖4 泥漿飛濺速度Fig.4 Mud splash velocity
圖5 泥漿飛濺體積量Fig.5 Mud splash volume
圖6 泥漿平均測面飛濺速度Fig.6 Mud average side surface splash velocity
(1)
(2)
表1 泥漿飛濺收斂項卡方檢驗分析Table 1 Chi-square test analysis of mud splash convergence
下面引入網(wǎng)格收斂指數(shù)(grid convergence index,GCI)對網(wǎng)格獨立性作進一步的討論分析。
網(wǎng)格收斂誤差ε為
(3)
式(3)中,f1、f2分別為細網(wǎng)格收斂解與粗網(wǎng)格收斂解。網(wǎng)格加密比定義為
(4)
式(4)中,hk為每個網(wǎng)格的平均間距,由式(4)計算得到:
(5)
式(5)中,ΔVi為每個網(wǎng)格單元的體積;Nk為每套網(wǎng)格的節(jié)點數(shù)。網(wǎng)格加密比還可以簡化為
(6)
網(wǎng)格收斂指數(shù)GCI定義為
(7)
式(7)中,F(xiàn)s為安全因子,當使用3套或3套以上網(wǎng)格來估算GCI時,Fs=1.25,P為收斂精度,取P=1.97。
GCI的計算結(jié)果如表2、表3所示。從A6開始,泥漿平均測面飛濺速度計算的GCI分別為2.60%、2.68%和1.57%,GCI均小于3%,泥漿飛濺量計算的GCI分別為2.77%、2.19%和1.88%,GCI均小于3%, 滿足其網(wǎng)格收斂準則。
綜上所述,卡方檢驗判定網(wǎng)格獨立性與GCI判定網(wǎng)格收斂性的結(jié)果是一致的,此方案在網(wǎng)格接近300×104時數(shù)值模擬的計算值網(wǎng)格收斂。
表2 泥漿飛濺量GCI收斂分析Table 2 GCI convergence of mud splash
注:hi為每個網(wǎng)格的平均間距,r為網(wǎng)格加密比,f(A)為泥漿飛濺量,ε為網(wǎng)格收斂誤差,GCI為網(wǎng)格收斂指數(shù)。
表3 泥漿平均測面飛濺速度GCI收斂分析Table 3 GCI convergence of mud average side surface splashing velocity
注:f(V)為混漿飛濺速度。
對卡方檢驗的準確性進行進一步的輔助驗證,試驗采用長方體泥漿槽,試驗所用泥漿為福建省福州市閩侯縣南通鎮(zhèn)水稻田泥漿,泥漿高度為刀軸中心,時間為10月,溫度為32 ℃,微耕機功率為4 kW,作業(yè)轉(zhuǎn)速為40 rad/s,進行3次試驗,1 s內(nèi)泥漿槽泥漿下降的高度乘以槽的長和寬即為泥漿飛濺量,試驗過程如圖7所示,試驗結(jié)果表明:1 s內(nèi)泥漿飛濺量為16.97 L,與A9誤差為8.72%。槽箱泥漿深度不夠、機器無法立即達到所需轉(zhuǎn)速等原因會引起泥漿飛濺量的減少,所以試驗驗證的數(shù)值和卡方檢驗的結(jié)果是趨近的。
圖7 泥漿槽泥漿飛濺試驗Fig.7 Experiment of mud splash in mud tank
(1)基于較高質(zhì)量的網(wǎng)格模型和合理的數(shù)值模擬方法,建立九種網(wǎng)格方案,對結(jié)果影響較大區(qū)域進行一定比例加密,泥漿飛濺平均速度和飛濺體積整體上均呈現(xiàn)出顯著的上升趨勢,且上升幅度越來越小,最后都在接近300×104網(wǎng)格左右能得到其網(wǎng)格無關性的結(jié)果。
(2)取標準差為一組中網(wǎng)格最少的參數(shù)的5%誤差來進行卡方檢驗,泥漿飛濺量和平均側(cè)面飛濺速度均從A6開始網(wǎng)格收斂。從A6開始的連續(xù)網(wǎng)格,網(wǎng)格收斂指數(shù)(GCI)均小于3%,當網(wǎng)格數(shù)目接近300萬時,數(shù)值模擬的計算值與網(wǎng)格數(shù)目無關。實驗驗證1 s泥漿飛濺量與數(shù)值模擬誤差為8.72%。
(3)通過網(wǎng)格收斂指數(shù)GCI來檢驗卡方檢驗用于網(wǎng)格無關性驗證的方法,得出的結(jié)果是一致的,同時和實驗驗證的結(jié)果是趨近的,可得在一定的誤差范圍內(nèi)借助卡方檢驗的方法驗證網(wǎng)無關性是合理的。
(4)GCI方法更注重前后兩項的聯(lián)系,而卡方檢驗注重于某一范圍內(nèi)數(shù)據(jù)前后的變化和聯(lián)系,能夠在不同加密倍數(shù)情況下來判斷網(wǎng)格收斂情況,本方法更加適用于找出趨于穩(wěn)定后的最佳網(wǎng)格數(shù)量。相對而言,卡方檢驗也更加復雜細致。對于采用的均值未知的卡方檢驗方法用于網(wǎng)格無關性檢驗的關鍵是標準差的取值,基于不同的模型,標準差取值大小需要進一步研究。