基于OpenFOAM的豬體空氣阻力系數模擬計算

吳雪飛 李 浩* 施正香

(1.中國農業大學 水利與土木工程學院，北京 100083；2.農業農村部設施農業工程重點實驗室，北京 100083；3.北京市畜禽健康養殖環境工程技術研究中心，北京 100083)

集約化生產中，畜禽動物對舍內氣流的阻擋影響了動物活動區域甚至整個舍內的溫度、濕度和污染物濃度等環境參數[1-5]，同時也阻礙通風系統的調控以及舍內氣流的組織。因此研究畜禽動物氣流阻力變化規律對畜禽舍內環境參數的調控和通風系統的運行有重要意義[6-7]。

氣流流速和流向、豬的體重和姿勢等對豬舍內氣流傳熱的影響已經得到研究[8-11]，但目前關于氣流流速和流向等因素對豬體氣流阻力影響的研究仍然較少。同時由于畜禽動物幾何形狀不規則及動物活動性的影響，使得影響豬群的氣流阻力的因素較多，畜禽動物周圍的氣流變得非常復雜[8]。而單只豬的氣流阻力影響因素相對較少，以單只豬氣流阻力的研究作為豬群阻力研究的基礎。

隨著計算機的快速發展，基于計算流體力學(Computational fluid dynamics，CFD)的數值模擬技術在畜禽舍通風系統中得到越來越多的應用[12]。這些CFD軟件分為商業軟件與開源軟件2種：商業CFD軟件，如ANSYS和Star CCM+，具有便利的圖形化界面，計算精度和穩定性較高，得到廣泛使用，但是商業軟件閉源封裝計算代碼，限制了使用者二次開發，同時高額的證書費用增加了使用成本；開源CFD軟件，如OpenFOAM，因其免費的證書、靈活的程序構架、良好的拓展性等優勢也得到了一定程度的推廣應用[13-14]。近年來，OpenFOAM在農業通風相關的CFD模擬應用中也有涉及[15-16]。然而，OpenFOAM計算穩定性低，軟件操作復雜(直接使用C++語言指令進行操作)，在農業行業的應用依然非常有限[17]， 尤其是針對農業環境領域開發準確性更高的數值模擬算法的潛力沒有得到發揮。因此本研究使用OpenFOAM軟件數值模擬單只豬對氣流的阻擋作用，探討開源軟件OpenFOAM在農業通風研究中的可行性。

綜上，本研究旨在基于開源CFD軟件OpenFOAM數值模擬技術，研究單只豬對氣流的阻擋作用，探究氣流流速和流向對豬體空氣阻力系數的影響規律，以期為豬群阻力研究提供基礎理論依據，為開源CFD軟件在農業建筑通風領域的深化應用提供支撐。

1 研究方法

1.1 幾何模型與計算域

1.1.1單只豬模型和圓柱體的構建

本研究參照體重約為90 kg的育肥豬(圖1(a))進行建模，為了方便建模和減輕網格計算負擔，略去豬的四肢、耳朵、尾巴，得到簡化豬體模型(圖1(b))，該簡化豬體模型尺寸約為：長1.17 m、高0.48 m、寬0.27 m，體表面積1.15 m2，其中簡化豬體模型代替全尺寸豬體模型進行模擬的準確性已經得到證實[10，18-19]。由于活體豬數據采集困難，虛擬風洞中模擬圓柱體作為模型驗證的方法簡單快速[20-21]，因此本研究利用圓柱體阻力系數半經驗公式計算所得的圓柱體阻力系數與模擬所得的阻力系數對比從而驗證模型的準確性。為了確保圓柱體與簡化豬體模型空氣阻力系數的可比性，圓柱體與簡化豬體模型體積相同；同時為保障圓柱體阻力系數半經驗公式的適用性，圓柱的長度與直徑之比為4.0[22]，因此構建直徑為0.37 m，高度1.48 m的圓柱體(圖1(c))。

圖1 實際豬模型，簡化豬體模型和圓柱體關系示意圖Fig.1 Diagram of the relationship between the actual pig model, simplified pig model and cylinder

1.1.2計算域及模型位置確定

本研究構建一個尺寸(長×寬×高)為10.0 m×6.0 m×2.4 m的虛擬風洞，虛擬風洞設置進風口(Inlet)和出風口(Outlet)各1個，壁面4個。豬體中心距離進風口4.0 m，出風口6.0 m，滿足單一物體距離風洞進風口大于物體特征長度的5倍(2.4 m)、距離出風口大于物體特征長度的10倍(4.8 m)的要求，距離壁面約為2.5 m,頂面約2.0 m,豬體的最大風洞阻塞率為1.8%，低于3.0%的設計要求[23]，距離地面0.5 m(豬只站立高度)。

本研究中圓柱體的位置因為試驗目的不同而在虛擬風洞中的位置不同：為了驗證模型的可靠性以及模擬結果的準確性，將圓柱體布置在虛擬風洞高度的中間位置(距離地面1.2 m)，圓柱體與其余壁面的距離與簡化豬體模型相同，從而保證圓柱體空氣阻力系數半經驗公式的測試條件；為檢驗圓柱體作為極度簡化豬體模型的可行性，圓柱體也布置在簡化豬體模型的相同位置進行模擬，以保障其與簡化豬體模型空氣阻力系數的可比性。

1.2 氣流流速與氣流流向

研究不同流速下豬體對氣流的阻擋作用，虛擬風洞進風口的氣流流速v設置6個水平，分別為0.2、0.5、1.0、1.5、2.0和2.5 m/s。氣流流向的改變通過豬體與氣流流向的夾角θ進行表征，θ定義為豬體長軸方向與虛擬風洞進風口氣流流向的夾角。本研究中，風洞氣流流向不發生變化，僅改變豬體的朝向來實現θ的調節，θ設置13個水平，分別為 0°，15°，…，180°，其中θ=0°，90°和180°分別代表豬體頭部垂直迎風，豬體側面與氣流流向垂直和豬體尾部垂直迎風的情況。

1.3 空氣阻力系數定義及驗算

物體受到的氣流阻力由幾何形狀引起的壓差阻力和流體粘性引起的摩擦阻力構成。工程上為了方便計算繞流阻力引入空氣阻力系數(Drag coefficient，Cd)。本研究使用豬體表面平均空氣阻力系數(Cd)量化豬體對氣流阻擋作用的大小，Cd越大，說明豬對氣流的阻擋作用越大?？諝庾枇ο禂涤嬎愎綖椋?/p>

Cd=2F/sv2

(1)

式中：F為物體受到的空氣阻力，N；s為物體迎風面積，m2；v為物體與氣流的相對速度，m/s。

用于模型對比驗證的圓柱體空氣阻力系數半經驗公式由Hongwu Tang等[24]總結已有研究的試驗數據并歸納得到式(2)：

(2)

式中：Re為雷諾數，表征流體流動情況。

使用空氣阻力系數相對誤差(Er)評判圓柱模擬結果的準確性，計算公式為：

(3)

式中：Cd,y為模擬所得的圓柱體空氣阻力系數；Cd,b為圓柱體半經驗公式計算得到的圓柱體空氣阻力系數。

1.4 網格處理

1.4.1網格劃分

OpenFOAM網格的生成主要由BlockMesh和SnappyHexMesh算法聯合處理，主要有4個步驟(圖2)：首先將STL(Stereolithography)格式的三角曲面的目標物生成為OpenFOAM可以識別的曲面網格(圖2(a))；其次，使用BlockMesh算法創建由計算域邊界填充整個區域的初始六面體背景網格(圖2(b))；然后在初始六面體背景網格的基礎上由SnappyHexMesh算法進行網格細化，網格細化可以提高網格精細度，以保證網格能準確計算流場信息，同時，為節約計算資源，網格細化的范圍限制在圍繞豬體展開的2.5 m×2.0 m×2.0 m的長方體空間區域(圖2(c))；最后檢查網格質量，優化網格,得到細化后的網格(圖2(d))。

1.簡化豬體模型；2.六面體背景網格；3.豬模型所在位置；4.網格細化區域；5.細化網格；6.不同級別網格之間的過渡網格；7.背景網格1.Simplified pig model; 2.Hexahedron background meshes; 3.Position of pig model; 4.Mesh refinement region; 5.Refined meshes; 6.Transition meshes between different levels of meshes; 7.Background meshes圖2 OpenFOAM網格生成過程Fig.2 The process of grid generation in OpenFOAM

1.4.2網格數量

對于相同的模型，網格數量主要由網格尺寸決定。為了確定網格尺寸并檢驗網格獨立性，本研究共設置12種網格級別(表1)，這12種網格級別具有相同的網格劃分模式，僅由于背景網格尺寸不同而導致網格細化區域的網格尺寸和總網格數不同。

表1 不同網格級別的網格尺寸及網格總數Table 1 Mesh size of different mesh resolutions and the total number of meshes

1.5 邊界條件

邊界條件直接影響到氣流流動的物理屬性和模擬結果，本研究中虛擬風洞進風口、出風口和4個壁面與豬模型表面均為壁面邊界(wall)，對豬模型近壁面區域和虛擬風洞地面區域使用壁面函數法進行處理(表2)。

表2 不同壁面的邊界條件Table 2 Boundary conditions in different walls

1.6 計算模型與求解

由Menter[25]提出的k-ω SST兩方程湍流模型適用于局部空氣阻力與固體表面之間的氣流模擬，并在相關研究中表現良好[26]，因此本研究使用k-ω SST模型作為模擬的湍流模型。OpenFOAM可以自由選擇求解格式等參數，從而能夠更加具有針對性地求解流場信息，但是容易造成求解的不穩定性，為了檢驗求解格式的穩定性，求解設置如下：數值模擬流體求解使用SIMPLE (Semi-implicit method for pressure-linked equation algorithms)算法，求解器使用OpenFOAM程序中的穩態求解器SimpleFoam，擴散項格式(LaplacianSchemes)設置為Guass linear corrected,表面差分格式(InterpolationSchemes)設置為Linear，梯度格式(GradSchemes)設置為Gauss linear,離散格式(DivSchemes)采用Bounded Guass upwind。

2 結果與討論

2.1 網格劃分對空氣阻力系數的影響

為檢驗網格獨立性，模擬得到v=1.0 m/s的流速下，12種網格處理模式圓柱體空氣阻力系數(表3)。當背景網格尺寸為300～500 mm時，網格數量較少，空氣阻力系數的相對誤差和變化幅度均較大，模擬準確度和穩定性均較差；背景網格尺寸為230～280 mm時，空氣阻力系數的相對誤差和變化幅度均較小，模擬穩定性和準確性均較高，如本研究中第7網格級別的圓柱體模擬所得的空氣阻力系數；背景網格尺寸繼續減小到150～220 mm時，網格數量增加，空氣阻力系數波動較大，說明在較小的網格尺寸下模擬的準確度和穩定性均較差。適當的網格尺寸保證模擬的準確度與穩定性，網格尺寸減小導致網格數量和計算成本的增加但不會提高計算精確度，Hong等[15]使用OpenFAOM模擬農業建筑自然通風時也得到了與本研究類似的結果。

表3 不同網格級別的圓柱體空氣阻力系數模擬結果Table 3 Simulated results of drag coefficient ofcylinder with different mesh levels

2.2 模型驗證與豬體模型對比

網格獨立性檢驗完成后需要驗證模型的準確性，分別在v=0.2、0.5、1.0、1.5、2.0和2.5 m/s流速下計算求解第7網格級別處理模式的圓柱體空氣阻力系數(表4)。圓柱體空氣阻力系數的相對誤差波動較大,原因可能是：由于半經驗公式是基于均勻來流的無限長直立圓柱體進行歸納，流體流動屬于二維，但由于實際圓柱體三維流動的復雜性仍可能導致模擬結果的誤差增加，盡管此次模擬的圓柱體長徑比達到了使用公式所要求的4.0從而減小三維流動的影響，但是不能完全避免三維流動造成的干擾；另外一方面，OpenFOAM作為開源算法，在計算迭代過程中對結果的優化有限，特別是考慮到實際工程中湍流流動復雜，準確的測量結果的獲得難度較大，所以該誤差范圍是合理正常的[27]?；谏鲜隹紤]，使用第7網格級別處理網格既保證網格獨立性，又保證了湍流模型和OpenFOAM計算求解的可靠性，驗證了模擬的準確性。

表4 不同氣流流速下第7網格級別的圓柱體空氣阻力系數模擬結果與相對誤差Table 4 Simulated drag coefficient and the relativeerror of cylinder model with the 7th meshlevel under different air velocities

為檢驗圓柱體作為極度簡化豬體模型的可行性，求解風洞中相同位置處的簡化豬體模型和圓柱體在不同流速下的空氣阻力系數(圖3)。v為0.2～1.5 m/s時，豬的空氣阻力系數Cd,p減小，v為1.5～2.5 m/s時，Cd,p增大，增值較緩(圖3(a))。圓柱的空氣阻力系數Cd,y隨著氣流流速的增大而減小。作為類圓柱體的簡化豬體和圓柱體周圍的氣流流動屬于亞臨界區，空氣阻力系數可以看作定值[24]，所以在v=0.5～2.5 m/s時取平均值得到Cd,p=0.943，Cd,y=1.207，二者大小差別明顯。在變化規律方面，線性擬合圓柱和豬的空氣阻力系數(圖3(b))得到y=0.600 3x+0.223 5,R2=0.603 3，說明圓柱體模型與簡化豬體模型阻力系數相關性不顯著。幾何形狀的差異導致簡化豬體模型(忽略豬體的四肢耳朵等非軀干部分的建模)相較于圓柱體更加能夠表征單體豬體阻力系數變化規律。

2.3 氣流流速和流向對豬體空氣阻力系數的影響

不同流速和豬體與氣流流向的不同夾角下豬體簡化模型空氣阻力系數模擬結果見圖4?？梢姡弘S著豬體迎風面積的增大，空氣阻力系數發生波動，但是存在增大趨勢。θ=45°時，流速較高的情況下Cd,p陡增，這是因為作為類圓柱體的豬體模型在豬體與氣流流向的夾角θ=45°時流場三維效應愈發明顯[28-29]。相同流速下，θ=0°、15°、30°、150°、165°和180°時，豬體氣流阻力較低；θ=45°、60°、120°和135°時，豬體阻力適中；θ=75°、90°和105°時，豬體阻力較大。由于豬舍生產過程中風速以不超過2.0 m/s 為宜，最佳風速在0.2 m/s左右，所以在實際生產中可以根據豬體與氣流流向的夾角范圍劃分阻力區間范圍從而評判豬體對氣流的阻擋作用。

圖4 豬體與氣流流向夾角(θ)和氣流流速(v)變化時豬體空氣阻力系數(Cd,p)的模擬結果Fig.4 Simulated results of drag coefficient of pig (Cd,p) under different airflow angles (θ) and air velocities (v)

當氣流流速v=0.2 m/s時，豬體相同迎風面積，頭部迎風和尾部迎風的情況下Cd,p相差較??；v=1.0或2.5 m/s時，頭部迎風和尾部迎風時Cd,p相差較大。以相同迎風面積，頭部迎風θ=60°和尾部迎風θ=120°為例，v=0.2 m/s時，Cd,p分別為0.87 與0.86，壓力云圖差別較小(圖5(a)、(b))；v=1.0 m/s時，Cd,p分別為0.73和0.85，壓力云圖差別較大(圖5(c)、(d))。這是因為流速較低，表面摩擦阻力是影響阻力的關鍵；流速較高，形狀阻力是影響阻力的關鍵，頭部和尾部迎風，迎風形狀不同則氣流變化明顯，空氣阻力變化較大。

圖5 距離風洞進風口2.4 m處不同豬體與氣流流向夾角(θ)和不同氣流流速(v)下的壓力云圖Fig.5 Pressure contour under different velocities (v) and airflow angles (θ) with 2.4 m from inlet of wind tunnel

3 結 論

本研究基于開源CFD(Computational fluid dynamics)算法OpenFOAM模擬風洞中單只豬對氣流的阻擋作用，比較網格大小對計算結果穩定性和準確性的影響，分析了不同氣流流速和豬體與氣流流向不同的夾角下的豬體表面的空氣阻力系數的變化規律，主要結論如下：

1)開源CFD軟件OpenFOAM能夠適用于畜禽舍空氣阻力的研究，網格處理模式對計算精度和穩定性影響明顯，合適的網格尺寸保證計算結果的穩定性和可靠性，減少網格數，節約計算成本。

2)流速較高時，豬體幾何形狀對氣流流動影響較大，圓柱體作為極度簡化的豬體模型不能夠充分表征豬體的空氣阻力系數的變化規律，二者的相關性較低。

3)氣流流向對豬體空氣阻力系數影響顯著，可以根據豬體與氣流流向的夾角范圍劃分阻力區間范圍從而評判豬體對氣流的阻擋作用的大小。