基于卡方檢驗的計算流體動力學網格無關性分析

謝 宇， 任金波， 黃煌輝， 張 翔

(福建農林大學機電工程學院,福州 350002)

計算流體動力學(computational fluid dynamics, CFD)[1],其技術取得了快速的發展[2-3],數值模擬方法的應用日趨廣泛[4-7],所面對的問題越來越復雜,但數值模擬在不同工程環境下的表現值不一樣[8],因此數值模擬的可靠性驗證是判斷數值模擬可信度的重要方面,計算結果可靠性分析方法主要有基準解比較法、實驗比較法和網格收斂性分析法三種[9-13],網格收斂性分析法中最流行的是Roache[14]提出通過計算網格收斂指數(grid convergence index,GCI) 來評估數值誤差的方法,現階段GCI方法已廣泛應用于數值計算誤差,Craig等[15]利用GCI對厭氧消化器攪拌模型進行誤差分析,Karimia等[16]分析了GCI 中的關鍵參數對旋液分離器數值誤差的影響,鄭秋亞等[17]針對M6機翼的繞流問題,就網格密度對 5 套不同網格的CFD模型進行了GCI分析, Longest等[18]使用GCI分析不同網格樣式對于生物流動系統的流速場等因素的影響,劉厚林等[19]借助GCI對比分析了不同網格類型的離心泵內部流場情況,GCI方法在一定程度上提高模擬的準確性并提供理論支持,但仍存在以下不足:① GCI方法在網格以整數倍形式加密時計算結果可靠,但是在實際操作中(尤其是多維模型時)往往很難做到,同時若網格加密倍數過小,解的變化不明顯；②在多維模型中,不同方向網格加密倍數可能不一樣,須在每個方向先分別計算相應的 GCI,然后再迭加求其和,而在實際操作中很難做到；③截差階數p較難確定；④漸近范圍的判斷沒有絕對值,網格加密到一定程度時計算結果會進入到漸近收斂范圍,再對網格進行加密沒有太大的意義[20-21],因此,確定合適的漸進收斂范圍很重要。

基于此，對刀片切泥制漿數值模擬進行CFD網格無關性的分析,探討在混合網格下網格無關性的網格數量,在允許誤差范圍內借助卡方檢驗驗證網格無關性,再與GCI方法進行對比,最終和試驗結果來判斷是否合理,探究適用于刀片切泥制漿模擬的網格無關性所對應的最優網格數,以分析在非整數加密以及加密方式不一致情況下的網格收斂分析準確性,以期為今后相關研究提供借鑒和啟發。

1 模型與方法

1.1 物理模型

以刀片旋轉切割泥漿來培漿的過程為研究對象,設計泥漿槽長寬高為2 000 mm×2 000 mm×1 600 mm,一側有一臺階,刀片轉速為40 rad/s,泥漿被旋轉的刀片切割飛濺至一側臺階上,來模擬現實中水稻田里泥漿被刀片旋轉打到育秧盤的工作環境,其物理模型如圖1所示。

圖1 刀片切泥制漿模型簡圖Fig.1 Model diagram of blade cutting mud

1.2 建模及網格劃分方法

用Pro/Engineer軟件對模型進行建模,使用CFD前處理軟件(the integrated computer engineering and manufacturing code for computational fluid dynamics, ICEMCFD)盡可能對結構模型采用結構化網格劃分,網絡模型的劃分方法和網格質量對數值模擬結果具有極大影響,網格quality是判定網格質量的重要方面,使大其數值接近1,采用混合網格,泥漿槽部分采用結構網格劃分,刀片部分采用非結構網格劃分,刀片的邊緣和刀面細化處理,如圖2所示。

圖2 刀片區域網格Fig.2 Domain of the blade grid

1.3 計算模型

采用商用軟件Fluent16.1,泥漿設置為賓漢姆模型,近壁面處理采用標準壁面函數；采用混合模型(mixture model)對氣液相互作用進行了表征；邊界均設置為Wall；采用滑移網格(sliding mesh)描述刀片旋轉；使用三維雙精度求解器；壓力-速度耦合方法采用SIMPLE算法及一階迎風格式；湍流模型采用RNGk-ε模型；時間步長取為0.000 1 s,計算時間為1 s,在每個時間步長內設置最代次數為20次,收斂殘差設為0.001。

圖3 泥漿體積分布Fig.3 Domain of the blade grid

2 數值模擬的網格無關性分析

現主要探討模擬計算結果中的泥漿飛濺量、泥漿飛濺速度的網格無關性檢驗。設計9套網格數量對重要部分(刀片區域)和非重要部分(泥漿槽流體區域)進行不同程度的加密處理來進行網格無關性驗證,模擬1 s后,觀察其模擬結果,圖3為刀軸橫截面處泥漿體積分布圖,圖4所示為刀軸橫截面處泥漿飛濺速度圖。 數值模擬數據如表1所示,泥漿飛濺體積量是飛過距離刀軸一側0.2 m橫截面的泥漿體積量,如圖5所示；平均側面飛濺速度是距離刀軸一側0.2 m處泥漿的平均飛濺速度,如圖6所示。由圖5、圖6可知，網格數量對仿真模擬結果具有較大影響,當網格數量接近300×104時,模擬結果變化差異不大。

圖4 泥漿飛濺速度Fig.4 Mud splash velocity

圖5 泥漿飛濺體積量Fig.5 Mud splash volume

圖6 泥漿平均測面飛濺速度Fig.6 Mud average side surface splash velocity

3 卡方檢驗

(1)

(2)

表1 泥漿飛濺收斂項卡方檢驗分析Table 1 Chi-square test analysis of mud splash convergence

4 網格收斂性指數分析

下面引入網格收斂指數(grid convergence index,GCI)對網格獨立性作進一步的討論分析。

網格收斂誤差ε為

(3)

式(3)中,f1、f2分別為細網格收斂解與粗網格收斂解。網格加密比定義為

(4)

式(4)中,hk為每個網格的平均間距,由式(4)計算得到:

(5)

式(5)中，ΔVi為每個網格單元的體積；Nk為每套網格的節點數。網格加密比還可以簡化為

(6)

網格收斂指數GCI定義為

(7)

式(7)中，Fs為安全因子,當使用3套或3套以上網格來估算GCI時,Fs=1.25,P為收斂精度,取P=1.97。

GCI的計算結果如表2、表3所示。從A6開始,泥漿平均測面飛濺速度計算的GCI分別為2.60%、2.68%和1.57%,GCI均小于3%,泥漿飛濺量計算的GCI分別為2.77%、2.19%和1.88%,GCI均小于3%, 滿足其網格收斂準則。

綜上所述,卡方檢驗判定網格獨立性與GCI判定網格收斂性的結果是一致的,此方案在網格接近300×104時數值模擬的計算值網格收斂。

表2 泥漿飛濺量GCI收斂分析Table 2 GCI convergence of mud splash

注：hi為每個網格的平均間距,r為網格加密比，f(A)為泥漿飛濺量，ε為網格收斂誤差，GCI為網格收斂指數。

表3 泥漿平均測面飛濺速度GCI收斂分析Table 3 GCI convergence of mud average side surface splashing velocity

注：f(V)為混漿飛濺速度。

5 試驗驗證

對卡方檢驗的準確性進行進一步的輔助驗證,試驗采用長方體泥漿槽,試驗所用泥漿為福建省福州市閩侯縣南通鎮水稻田泥漿,泥漿高度為刀軸中心,時間為10月,溫度為32 ℃,微耕機功率為4 kW,作業轉速為40 rad/s,進行3次試驗,1 s內泥漿槽泥漿下降的高度乘以槽的長和寬即為泥漿飛濺量,試驗過程如圖7所示,試驗結果表明:1 s內泥漿飛濺量為16.97 L,與A9誤差為8.72%。槽箱泥漿深度不夠、機器無法立即達到所需轉速等原因會引起泥漿飛濺量的減少,所以試驗驗證的數值和卡方檢驗的結果是趨近的。

圖7 泥漿槽泥漿飛濺試驗Fig.7 Experiment of mud splash in mud tank

6 結論

(1)基于較高質量的網格模型和合理的數值模擬方法,建立九種網格方案,對結果影響較大區域進行一定比例加密,泥漿飛濺平均速度和飛濺體積整體上均呈現出顯著的上升趨勢,且上升幅度越來越小,最后都在接近300×104網格左右能得到其網格無關性的結果。

(2)取標準差為一組中網格最少的參數的5%誤差來進行卡方檢驗,泥漿飛濺量和平均側面飛濺速度均從A6開始網格收斂。從A6開始的連續網格,網格收斂指數(GCI)均小于3%,當網格數目接近300萬時,數值模擬的計算值與網格數目無關。實驗驗證1 s泥漿飛濺量與數值模擬誤差為8.72%。

(3)通過網格收斂指數GCI來檢驗卡方檢驗用于網格無關性驗證的方法,得出的結果是一致的,同時和實驗驗證的結果是趨近的,可得在一定的誤差范圍內借助卡方檢驗的方法驗證網無關性是合理的。

(4)GCI方法更注重前后兩項的聯系,而卡方檢驗注重于某一范圍內數據前后的變化和聯系,能夠在不同加密倍數情況下來判斷網格收斂情況,本方法更加適用于找出趨于穩定后的最佳網格數量。相對而言,卡方檢驗也更加復雜細致。對于采用的均值未知的卡方檢驗方法用于網格無關性檢驗的關鍵是標準差的取值,基于不同的模型,標準差取值大小需要進一步研究。