基于OpenFOAM的豬體空氣阻力系數(shù)模擬計(jì)算

吳雪飛 李 浩* 施正香

(1.中國農(nóng)業(yè)大學(xué) 水利與土木工程學(xué)院，北京 100083；2.農(nóng)業(yè)農(nóng)村部設(shè)施農(nóng)業(yè)工程重點(diǎn)實(shí)驗(yàn)室，北京 100083；3.北京市畜禽健康養(yǎng)殖環(huán)境工程技術(shù)研究中心，北京 100083)

集約化生產(chǎn)中，畜禽動(dòng)物對(duì)舍內(nèi)氣流的阻擋影響了動(dòng)物活動(dòng)區(qū)域甚至整個(gè)舍內(nèi)的溫度、濕度和污染物濃度等環(huán)境參數(shù)[1-5]，同時(shí)也阻礙通風(fēng)系統(tǒng)的調(diào)控以及舍內(nèi)氣流的組織。因此研究畜禽動(dòng)物氣流阻力變化規(guī)律對(duì)畜禽舍內(nèi)環(huán)境參數(shù)的調(diào)控和通風(fēng)系統(tǒng)的運(yùn)行有重要意義[6-7]。

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

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

綜上，本研究旨在基于開源CFD軟件OpenFOAM數(shù)值模擬技術(shù)，研究單只豬對(duì)氣流的阻擋作用，探究氣流流速和流向?qū)ωi體空氣阻力系數(shù)的影響規(guī)律，以期為豬群阻力研究提供基礎(chǔ)理論依據(jù)，為開源CFD軟件在農(nóng)業(yè)建筑通風(fēng)領(lǐng)域的深化應(yīng)用提供支撐。

1 研究方法

1.1 幾何模型與計(jì)算域

1.1.1單只豬模型和圓柱體的構(gòu)建

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

圖1 實(shí)際豬模型，簡(jiǎn)化豬體模型和圓柱體關(guān)系示意圖Fig.1 Diagram of the relationship between the actual pig model, simplified pig model and cylinder

1.1.2計(jì)算域及模型位置確定

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

本研究中圓柱體的位置因?yàn)樵囼?yàn)?zāi)康牟煌谔摂M風(fēng)洞中的位置不同：為了驗(yàn)證模型的可靠性以及模擬結(jié)果的準(zhǔn)確性，將圓柱體布置在虛擬風(fēng)洞高度的中間位置(距離地面1.2 m)，圓柱體與其余壁面的距離與簡(jiǎn)化豬體模型相同，從而保證圓柱體空氣阻力系數(shù)半經(jīng)驗(yàn)公式的測(cè)試條件；為檢驗(yàn)圓柱體作為極度簡(jiǎn)化豬體模型的可行性，圓柱體也布置在簡(jiǎn)化豬體模型的相同位置進(jìn)行模擬，以保障其與簡(jiǎn)化豬體模型空氣阻力系數(shù)的可比性。

1.2 氣流流速與氣流流向

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

1.3 空氣阻力系數(shù)定義及驗(yàn)算

物體受到的氣流阻力由幾何形狀引起的壓差阻力和流體粘性引起的摩擦阻力構(gòu)成。工程上為了方便計(jì)算繞流阻力引入空氣阻力系數(shù)(Drag coefficient，Cd)。本研究使用豬體表面平均空氣阻力系數(shù)(Cd)量化豬體對(duì)氣流阻擋作用的大小，Cd越大，說明豬對(duì)氣流的阻擋作用越大。空氣阻力系數(shù)計(jì)算公式為：

Cd=2F/sv2

(1)

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

用于模型對(duì)比驗(yàn)證的圓柱體空氣阻力系數(shù)半經(jīng)驗(yàn)公式由Hongwu Tang等[24]總結(jié)已有研究的試驗(yàn)數(shù)據(jù)并歸納得到式(2)：

(2)

式中：Re為雷諾數(shù)，表征流體流動(dòng)情況。

使用空氣阻力系數(shù)相對(duì)誤差(Er)評(píng)判圓柱模擬結(jié)果的準(zhǔn)確性，計(jì)算公式為：

(3)

式中：Cd,y為模擬所得的圓柱體空氣阻力系數(shù)；Cd,b為圓柱體半經(jīng)驗(yàn)公式計(jì)算得到的圓柱體空氣阻力系數(shù)。

1.4 網(wǎng)格處理

1.4.1網(wǎng)格劃分

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

1.簡(jiǎn)化豬體模型；2.六面體背景網(wǎng)格；3.豬模型所在位置；4.網(wǎng)格細(xì)化區(qū)域；5.細(xì)化網(wǎng)格；6.不同級(jí)別網(wǎng)格之間的過渡網(wǎng)格；7.背景網(wǎng)格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網(wǎng)格生成過程Fig.2 The process of grid generation in OpenFOAM

1.4.2網(wǎng)格數(shù)量

對(duì)于相同的模型，網(wǎng)格數(shù)量主要由網(wǎng)格尺寸決定。為了確定網(wǎng)格尺寸并檢驗(yàn)網(wǎng)格獨(dú)立性，本研究共設(shè)置12種網(wǎng)格級(jí)別(表1)，這12種網(wǎng)格級(jí)別具有相同的網(wǎng)格劃分模式，僅由于背景網(wǎng)格尺寸不同而導(dǎo)致網(wǎng)格細(xì)化區(qū)域的網(wǎng)格尺寸和總網(wǎng)格數(shù)不同。

表1 不同網(wǎng)格級(jí)別的網(wǎng)格尺寸及網(wǎng)格總數(shù)Table 1 Mesh size of different mesh resolutions and the total number of meshes

1.5 邊界條件

邊界條件直接影響到氣流流動(dòng)的物理屬性和模擬結(jié)果，本研究中虛擬風(fēng)洞進(jìn)風(fēng)口、出風(fēng)口和4個(gè)壁面與豬模型表面均為壁面邊界(wall)，對(duì)豬模型近壁面區(qū)域和虛擬風(fēng)洞地面區(qū)域使用壁面函數(shù)法進(jìn)行處理(表2)。

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

1.6 計(jì)算模型與求解

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

2 結(jié)果與討論

2.1 網(wǎng)格劃分對(duì)空氣阻力系數(shù)的影響

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

表3 不同網(wǎng)格級(jí)別的圓柱體空氣阻力系數(shù)模擬結(jié)果Table 3 Simulated results of drag coefficient ofcylinder with different mesh levels

2.2 模型驗(yàn)證與豬體模型對(duì)比

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

表4 不同氣流流速下第7網(wǎng)格級(jí)別的圓柱體空氣阻力系數(shù)模擬結(jié)果與相對(duì)誤差Table 4 Simulated drag coefficient and the relativeerror of cylinder model with the 7th meshlevel under different air velocities

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

2.3 氣流流速和流向?qū)ωi體空氣阻力系數(shù)的影響

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

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

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

圖5 距離風(fēng)洞進(jìn)風(fēng)口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 結(jié) 論

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

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

2)流速較高時(shí)，豬體幾何形狀對(duì)氣流流動(dòng)影響較大，圓柱體作為極度簡(jiǎn)化的豬體模型不能夠充分表征豬體的空氣阻力系數(shù)的變化規(guī)律，二者的相關(guān)性較低。

3)氣流流向?qū)ωi體空氣阻力系數(shù)影響顯著，可以根據(jù)豬體與氣流流向的夾角范圍劃分阻力區(qū)間范圍從而評(píng)判豬體對(duì)氣流的阻擋作用的大小。