劉德釗 辛宜聰 榮 莉 王朝元 李保明

(1.浙江大學(xué) 生物系統(tǒng)工程與食品科學(xué)學(xué)院，杭州 310058；2.奧胡斯大學(xué) 土木與建筑學(xué)院，丹麥 奧胡斯 8000；3.中國農(nóng)業(yè)大學(xué) 水利與土木工程學(xué)院，北京 100083;4.農(nóng)業(yè)農(nóng)村部設(shè)施農(nóng)業(yè)裝備與信息化重點實驗室，杭州 310058)

良好的豬舍內(nèi)環(huán)境對豬健康生長和生產(chǎn)力的提高至關(guān)重要。計算流體力學(xué)(Computational fluid dynamics, CFD)近年來被廣泛用于畜禽舍內(nèi)氣流、溫濕度和氣體濃度的模擬研究，已成為研究舍內(nèi)環(huán)境的重要手段之一[1-3]。在以往對豬舍內(nèi)環(huán)境的CFD模擬中，為簡化計算豬體模型通常被省略[4-5]。然而，在對豬舍進行模擬時，將豬體存在因素考慮在內(nèi)對準(zhǔn)確模擬舍內(nèi)空氣流動非常重要。因此有些研究嘗試在舍內(nèi)環(huán)境模擬中建立實際尺寸的豬體模型[6-7]，如Seo等[8]對全尺寸商業(yè)豬場中648頭豬進行詳細(xì)建模，采用CFD模擬探究了舍內(nèi)溫度分布和氣流模式。然而，在有大量動物的情況下采用這種方式進行建模，普通臺式機以及專業(yè)的Workstation計算機根本無法滿足如此巨大的計算量要求，同時計算時間會成倍增加。

采用多孔介質(zhì)模型簡化動物區(qū)域可以有效解決這一難題[9-10]。應(yīng)用多孔介質(zhì)模型，需要先獲得模型設(shè)定所需的阻力系數(shù)，即多孔介質(zhì)區(qū)域的入口速度與進出口壓降之間的關(guān)系。Du等[11]求取養(yǎng)雞場的籠養(yǎng)雞區(qū)域的阻力系數(shù)，將其簡化為多孔介質(zhì)模型，探究養(yǎng)雞場內(nèi)氣流、溫度和相對濕度的分布，結(jié)果顯示模擬與實測結(jié)果吻合良好。Cheng等[12]探究蛋雞的幾何形狀、分布以及體重對阻力系數(shù)的影響，為將蛋雞區(qū)域簡化為多孔介質(zhì)模型提供了參考。在養(yǎng)豬場中，直接測量養(yǎng)豬區(qū)域的風(fēng)速以及壓降存在較大困難，因此通常采用CFD模擬獲得不同風(fēng)速下養(yǎng)豬區(qū)域的壓降進而計算多孔介質(zhì)模型所需的阻力系數(shù)[13]。然而，豬體的不同幾何模型假設(shè)對氣流產(chǎn)生的阻力存在差異。在此前的研究中，學(xué)者們通常會將豬體模型進行簡化，省略細(xì)小的部位，例如尾巴、耳朵、四肢，在保證模擬準(zhǔn)確度的情況下節(jié)約計算成本[14]。Li等[15]將省略了細(xì)小部位的豬體模型用于對流換熱損失的模擬研究，可是在實際畜禽場的模擬中，由于動物的實際尺寸和不規(guī)則的形狀，省略細(xì)小部位的動物簡化模型仍然會產(chǎn)生較大的計算量[16]。因此，學(xué)者們開始使用簡單的幾何模型代替動物模型，如圓柱體、橢球體和球體等[17-18]。Mondaca與Choi[16]發(fā)現(xiàn)球體、圓柱體和6個圓柱的裝配體均能較好的模擬奶牛平均熱通量。Li等[19]將實際尺寸的豬模型與圓柱體模型進行對比，結(jié)果表明2種模型的對流換熱系數(shù)具有較強的相關(guān)性，可以采用圓柱模型作為豬的簡化模型。但是該研究僅考慮了豬在站立時的情況，而在實際養(yǎng)殖過程中育肥豬大約80%的時間都處于躺臥狀態(tài)[20]，而且目前仍沒有研究表明何種簡化幾何體可代表豬的躺臥模型以及不同幾何體簡化模型之間阻力系數(shù)有何差異性。因此，研究豬在躺臥姿勢下對氣流的阻礙作用是非常有必要的。

本研究以躺臥的育肥豬作為研究對象，將簡化后的豬體躺臥模型作為參考模型，并選擇半橢球、橢球、半圓柱3種簡單幾何體分別代替簡化豬模型，計算得到的阻力系數(shù)可用于之后的豬舍全尺寸模擬，探究4種模型的阻力系數(shù)之間的差異以及對動物區(qū)域氣流以及壓力分布的影響，為多孔介質(zhì)模型在動物區(qū)域的使用和動物模型的簡化與選擇提供參考。

1 材料與方法

1.1 CFD模擬方法驗證

為確保在對不同躺臥豬模型的探究過程中所采用的模擬方法的有效性，參考Cheng等[12]已發(fā)表文獻中的風(fēng)洞實測值對本研究的CFD模擬進行驗證。Cheng等采用直徑為150 mm的球體代替1.5 kg的雞模型，在長為6 910 mm、寬和高均為500 mm的風(fēng)洞內(nèi)進行試驗，主要測量了雞區(qū)域內(nèi)5條垂直測量線(L2～L6)的風(fēng)速值，再對此風(fēng)洞試驗進行CFD模擬。這與本研究中豬區(qū)域以及計算域的氣流特性相近，且同樣分為X、Y和Z3 個方向分別求取阻力系數(shù)。CFD模擬方法的驗證采用探究不同躺臥豬模型阻力系數(shù)時所用的模擬方法對文獻中的風(fēng)洞試驗進行CFD模擬，參考Cheng等的幾何模型及邊界設(shè)置，幾何模型如圖1所示。采用非結(jié)構(gòu)化四面體網(wǎng)格對計算域進行網(wǎng)格劃分，并選擇3種數(shù)量的網(wǎng)格對驗證模擬進行網(wǎng)格獨立性檢驗，高密度網(wǎng)格數(shù)量為3 574 326，中等密度網(wǎng)格數(shù)量為2 708 984，低密度網(wǎng)格數(shù)量為1 769 593，模擬所得到的雞區(qū)域進出口的壓降依次分別是0.381、0.383 和0.398 Pa，由于高密度網(wǎng)格與中等密度網(wǎng)格設(shè)定下壓降的誤差僅為-0.52%，因此選擇中等密度的網(wǎng)格進行驗證模擬。采用標(biāo)準(zhǔn)k-ε湍流模型和scalable壁面函數(shù)，壓力-速度耦合選用SIMPLE算法，動量、湍動能和湍流耗散率選用二階迎風(fēng)離散格式。將本研究CFD模擬得到的風(fēng)速廓線與文獻中的風(fēng)洞實測值進行對比，計算得到平均相對誤差，用于驗證本研究模擬方法的有效性，平均相對誤差計算公式如下式：

圖1 驗證幾何模型Fig.1 Geometric model of validation

(1)

式中：MRE為平均相對誤差；n為測點數(shù)；xi為第i個點的速度模擬值，m/s；yi為第i個點的速度測量值，m/s。

1.2 豬躺臥模型的阻力系數(shù)

1.2.1豬的躺臥模型

本研究選擇55 kg育肥豬作為研究對象，研究表明在豬的所有躺臥姿勢中，側(cè)躺且四肢伸展姿勢占比高于60%[21]，因此本研究中豬的躺臥模型均基于此種姿勢。豬重量與體尺之間的冪函數(shù)關(guān)系如下[22]：

L=(26.8±0.7)M(0.33±0.01)

(2)

H=(8.7±0.2)M(0.30±0.01)

(3)

W=(7.7±0.2)M(0.37±0.01)

(4)

式中：L為躺臥豬體總長度，cm；H為躺臥豬體高度，cm；W為躺臥豬體寬度，cm；M為豬的體重，kg。

由上述公式計算得到豬主要身體部位的尺寸，結(jié)合實際中豬的形狀進行簡化并建立了躺臥時的簡化豬模型。以簡化豬模型作為參考模型，將簡化豬模型繼續(xù)簡化為3種簡單幾何模型，分別為半橢球模型、橢球模型和半圓柱模型，根據(jù)簡化豬模型的體尺建立半橢球模型、橢球模型和半圓柱模型，再對這3種簡單幾何模型的體尺進行微調(diào)使得這3種模型的體積都與簡化豬模型的體積相同，如圖2所示。考慮到圓柱體放置在地板上時與地板間的空隙較大，不符合豬在躺臥時的實際情況，因此沒有選擇圓柱體模型。豬的尾巴、耳朵和四肢等部位在劃分網(wǎng)格時會顯著增加網(wǎng)格數(shù)量，而這些部位在豬躺臥的情況下對氣流的影響較小，因此在建模時省略。

圖2 豬的4種躺臥模型Fig.2 Four geometric models of lying pigs

1.2.2豬的分布

豬的數(shù)量和在豬欄內(nèi)的分布對氣流有非常大的影響。模擬所選擇的豬場位于浙江省義烏市(29°21′N，119°94′E)，根據(jù)此豬場的實際養(yǎng)殖規(guī)模以及內(nèi)部結(jié)構(gòu)，其中單個豬欄長L=5 680 mm、寬W=4 850 mm，豬欄高度根據(jù)豬站立時的高度H設(shè)置為700 mm[9]。豬欄中實心地板區(qū)域和漏縫地板區(qū)域分別占豬欄總面積的30.99%和69.01%，共有27頭豬。豬的躺臥方向參考豬舍內(nèi)不同時段拍攝的實際照片得到。根據(jù)Bjerg等[9]提出的豬欄內(nèi)豬的分布規(guī)律，實心地板區(qū)域應(yīng)有13頭豬躺臥、2頭豬站立，漏縫地板區(qū)域應(yīng)有9只豬躺臥、3只豬站立。結(jié)合Li等[15]在CFD模擬中采用的豬舍內(nèi)分布和此研究中實際豬和豬欄的大小，本次模擬的豬的分布以半橢球模型分布為例，如圖3所示。本研究的研究對象僅為豬的躺臥模型，因此4種不同躺臥模型下豬的分布完全相同，僅躺臥模型不同，豬欄內(nèi)共有5頭豬站立，均采用簡化豬模型，離地277.9 mm。

圖3 半橢球躺臥豬模型分布Fig.3 Distribution of half ellipsoid lying pig model

1.2.3計算域與邊界條件

計算豬體所在區(qū)域X、Y和Z3個方向上的阻力系數(shù)的計算域如圖4所示。進風(fēng)口設(shè)置為速度入口邊界，風(fēng)速取0.05、0.10、0.15、0.20和0.25 m/s。出風(fēng)口設(shè)置為壓力出口邊界，豬的表面設(shè)置為壁面邊界，其余面均設(shè)置為對稱邊界。

圖4 不同方向的計算域Fig.4 Computational domains of different directions

1.2.4網(wǎng)格劃分與數(shù)值求解

基于CFD驗證模擬中網(wǎng)格劃分經(jīng)驗對計算域進行網(wǎng)格劃分，同樣采用非結(jié)構(gòu)化四面體網(wǎng)格?？刂品匠滩捎没谟邢摅w積的離散方法，即采用雷諾平均N-S方程(RANS)。同時，本模擬研究采用標(biāo)準(zhǔn)k-ε湍流模型和scalable壁面函數(shù)。標(biāo)準(zhǔn)k-ε湍流模型自從被提出之后就變成工程流場計算中主要的湍流模型，scalable壁面函數(shù)能避免計算結(jié)果惡化，對于任意細(xì)化的網(wǎng)格能給出一致的解。此外，求解壓力耦合方程的半隱方法(SIMPLE)用于壓力-速度場修正，而動量、湍動能和湍流耗散率選用二階迎風(fēng)離散格式。

1.2.5阻力系數(shù)的計算

CFD模擬通常將動物所在區(qū)域簡化為多孔介質(zhì)模型進行求解，能夠降低建模以及網(wǎng)格計算的難度。準(zhǔn)確的阻力系數(shù)對多孔介質(zhì)模型的設(shè)置至關(guān)重要，包括黏性阻力系數(shù)與慣性阻力系數(shù)。有限厚度的多孔介質(zhì)模型的阻力系數(shù)如下式Darcy定律所示[13]：

(5)

式中：Δp為壓力降，Pa；1/α為黏性阻力系數(shù)，1/m2；C2為慣性阻力系數(shù)，1/m；μ為空氣動力粘度，Pa/s；ν為風(fēng)速，m/s；ρ為空氣密度，kg/m3；Δn為多孔介質(zhì)區(qū)域厚度，m。

將風(fēng)速與模擬所得的進出口的壓降擬合成常數(shù)項為0的一元二次方程，其二次項與一次項前的系數(shù)通過式(5)可解得慣性阻力系數(shù)與黏性阻力系數(shù)。

1.3 多孔介質(zhì)模型驗證

為了能夠?qū)⑻膳P豬模型簡化為多孔介質(zhì)模型，在求解得到阻力系數(shù)之后還需要采用多孔介質(zhì)模型驗證其準(zhǔn)確性。將養(yǎng)豬區(qū)域設(shè)置為多孔介質(zhì)區(qū)域并設(shè)置阻力系數(shù)，其余設(shè)置均與躺臥豬模型的阻力系數(shù)求解相同，在不同風(fēng)速下進行CFD模擬，將所得驗證壓降值與不同躺臥豬模型的模擬壓降值進行對比，驗證阻力系數(shù)的準(zhǔn)確性。

2 結(jié)果與分析

2.1 CFD模擬方法驗證

圖5為CFD模擬得到的5條不同測量線(L2～L6)的風(fēng)速廓線與文獻中風(fēng)洞實測值的對比。從圖中可以看出，在大多數(shù)情況下模擬值與實測值吻合較好，在L5實測值與模擬值之間的差異性較大，這是因為標(biāo)準(zhǔn)k-ε模型本身不能很好模擬回流區(qū)所造成的。本研究中模擬值與實測值的平均相對誤差為8.41%，Cheng等采用相同的湍流模型所得誤差為8.10%。因此，本研究中采用的模擬方法能夠合理地模擬動物區(qū)域的氣流運動。

圖5 不同測量線的模擬風(fēng)速廓線與實驗測量值Fig.5 Simulated velocities profiles and experimental measurements along different measuring lines

2.2 不同幾何模型對單位長度壓降的影響

圖6為簡化豬模型、半橢球模型、橢球模型和半圓柱模型在X、Y和Z方向上速度與單位長度壓降之間的關(guān)系以及擬合方程。由圖6可以看出，單位長度的壓降隨速度的增加而增大，并且在同一方向上，速度越大，4種模型的單位長度壓降值相差越大。以簡化豬模型作為參考模型，簡化豬模型與半橢球模型之間的相對誤差在X方向上為-21.6%～-21.2%，在Z方向上為-24.5%～-23.5%，在Y方向上為6.9%～7.0%。簡化豬模型與橢球模型之間的相對誤差在X方向上為-2.9%～-1.1%，在Z方向上為-11.3%～-9.0%，在Y方向上為14.7%～13.6%。簡化豬模型與半圓柱模型之間的相對誤差在X方向上為-6.9%～-6.6%，在Z方向上為-11.6%～-10.3%，在Y方向上為30.6%～31.6%。由此可見，在X和Z方向上，簡化豬模型與橢球模型的單位長度壓降值的相對誤差最小，與半橢球模型之間的相對誤差最大，這可能是由于半橢球模型底面與地面完全重合，氣流沿著模型表面流線型邊緣線流動，減少了對空氣的阻力作用，導(dǎo)致單位長度的壓降值最小。在圖6(a)和6(b)中，單位長度的壓降值從大到小依次為簡化豬模型、橢球模型、半圓柱模型，但是在圖6(c)中，半圓柱模型單位長度壓降值最大，這可能是因為半圓柱模型在Y方向上的投影面積最大，而投影面積與阻力成正相關(guān)，投影面積越大，對氣流阻力越大，單位長度壓降越大。

圖6 4種躺臥豬模型在X、Y和Z方向上速度與單位長度壓降之間的關(guān)系Fig.6 Relationship between velocity and pressure drop per unit length of four lying pig models in X, Y and Z directions

2.3 不同幾何模型對阻力系數(shù)的影響

表1為4種模型在不同方向上的慣性阻力系數(shù)與黏性阻力系數(shù)，由表1可知，在X、Y和Z3種方向上，4種模型均在Y方向上的慣性阻力系數(shù)最大，X方向上的慣性阻力系數(shù)最小，這可能是因為4種模型在Y方向上的投影均大于X和Z方向，這與Cheng等[12]研究簡化雞模型阻力系數(shù)的結(jié)果一致。橢球模型與簡化豬模型的慣性阻力系數(shù)在X和Z2個方向上最為接近，其在X、Z和Y方向上的慣性阻力系數(shù)與簡化豬模型之間的相對誤差分別為-4.0%、-12.2%和14.7%。黏性阻力系數(shù)值在表1中無明顯規(guī)律，且將風(fēng)速值代入公式(5)中經(jīng)過計算可知一次項對單位長度壓降的貢獻較小，因此沒有對黏性阻力系數(shù)進行詳細(xì)分析。根據(jù)阻力系數(shù)可知，橢球模型與簡化豬模型最相似，因此在今后對多孔介質(zhì)模型的應(yīng)用研究中，為了降低建模的難度，可以使用橢球模型代替簡化豬模型對阻力系數(shù)進行求解，并且為今后利用多孔介質(zhì)模型模擬全尺寸豬舍的CFD研究提供了可靠的數(shù)據(jù)。本研究中所得阻力系數(shù)的結(jié)果是在假設(shè)豬的體重為55 kg的情況下計算得到的，并不能代表豬在其他體重下的阻力系數(shù)，在后續(xù)的研究中將會對不同體重下豬體模型的阻力系數(shù)進行系統(tǒng)性的研究。

2.4 不同幾何模型對氣流分布的影響

圖7為不同幾何模型在X方向上，風(fēng)速為0.25 m/s 時豬欄內(nèi)y=160 mm平面上的風(fēng)速分布，由圖中可以看出，不同豬體模型會影響風(fēng)速的分布。簡化豬模型、半橢球模型、橢球模型和半圓柱模型在此平面的最大風(fēng)速分別為0.42、0.46、0.44 和0.45 m/s，結(jié)合圖中可知，半橢球與半圓柱對風(fēng)速分布影響較為明顯，模型背風(fēng)面低風(fēng)速區(qū)域與高風(fēng)速區(qū)域明顯多于其他2種模型，這可能是因為半圓柱與半橢球模型底面與豬欄地板重合，對氣流完全阻擋，造成背風(fēng)面處風(fēng)速更低，其在X方向上慣性阻力系數(shù)較低，因此高風(fēng)速區(qū)域也相對較多。橢球模型與簡化豬模型對氣流影響最相似，與表1中阻力系數(shù)結(jié)果一致。然而在Cheng等[12]的研究中，通過分析雞模型與橢球模型的風(fēng)速分布圖，最大風(fēng)速分別為5.78 和4.30 m/s，得出橢球模型不能代替雞模型的結(jié)論，這可能是因為雞模型相對于橢球模型來說更為復(fù)雜，且橢球模型將雞的脖子以上部分全部省略，導(dǎo)致2種模型相差較大，并且雞在雞籠中的分布更為緊密，孔隙率低，因此模型的不同會導(dǎo)致在風(fēng)速分布上的顯著差異。而對于本研究，躺臥時豬的體態(tài)相比雞更為圓潤，豬在躺臥時的高度僅約為站立時的一半，因此上方空間更大，所以橢球模型與豬模型對氣流分布的影響并無顯著差異。此外，Li等[23]采用CFD探究雞的對流換熱系數(shù)，結(jié)果表明，球體模型與雞模型的對流換熱系數(shù)擬合較好，在后續(xù)研究中可以采用球體模型代替雞模型。因此，在簡化模型時需要根據(jù)不同動物的特征以及在舍內(nèi)的分布情況進行簡化，同時模型的簡化程度也受到研究目的的影響。

表1 4種躺臥豬模型在X、Y和Z方向上的阻力系數(shù)Table 1 Resistant coefficients of four lying pig models in X, Y and Z directions

圖7 4種躺臥豬模型在y=160 mm平面上的風(fēng)速分布Fig.7 Velocity distribution of four lying pig models in the plane of y=160 mm

2.5 不同幾何模型對壓力分布的影響

圖8為不同幾何模型在X方向上，風(fēng)速為 0.25 m/s 時豬欄內(nèi)y=160 mm平面上的壓力分布。因為4種豬模型的分布完全相同，所以壓力分布的變化趨勢相同，但是局部壓力值存在差異。簡化豬模型、半橢球模型、橢球模型和半圓柱模型在此平面上的最大壓力分別為0.086、0.068、0.097 和0.081 Pa。從圖中可以看出，豬模型的存在使得在豬的身后產(chǎn)生了氣流分離與低壓區(qū)，較大的壓力值均出現(xiàn)在靠近進風(fēng)口的豬模型的迎風(fēng)面，而在遠(yuǎn)離進風(fēng)口的位置處則出現(xiàn)負(fù)壓。在豬分布密集區(qū)域，4種模型的壓力分布差異不明顯，而在豬模型分布稀疏的區(qū)域，模型迎風(fēng)面產(chǎn)生的較大壓力值差異顯著。這與Bustos-Vanegas等[24]采用多孔介質(zhì)模型模擬奶牛區(qū)域的空氣流動的研究中得到的壓力分布規(guī)律一致，并且提出雖然動物模型與多孔介質(zhì)模型得到的局部壓力相差較大，但是總壓降是相同的，說明多孔介質(zhì)模型仍是適用的。

圖8 4種躺臥豬模型在y=160 mm平面上的壓力分布Fig.8 Pressure distribution of four lying pig models in the plane of y=160 mm

2.6 多孔介質(zhì)模型驗證

圖9為簡化豬模型在X、Y和Z方向上的模擬壓降值與多孔介質(zhì)模型壓降值的對比。由圖可知，簡化豬模型的單位長度壓降值與采用簡化豬模型阻力系數(shù)的多孔介質(zhì)模型的單位長度壓降值吻合度高。以簡化豬模型為例，可以說明采用CFD模擬的方法求解阻力系數(shù)是有效的，同時驗證了阻力系數(shù)的正確性。在之后的研究中可以直接采用多孔介質(zhì)模型代替復(fù)雜的動物模型。

圖9 X、Y和Z方向上簡化豬模型的多孔介質(zhì)模型驗證Fig.9 Validation of porous media model for simplified pig model in X, Y and Z directions

3 結(jié) 論

本研究對簡化豬模型、半橢球模型、橢球模型和半圓柱模型4種躺臥豬模型進行了CFD模擬，探究了不同躺臥豬模型對阻力系數(shù)的影響。得出以下結(jié)論：

1)采用本研究的CFD模擬方法得到的風(fēng)速值與文獻中的風(fēng)速實測值吻合較好，平均相對誤差為8.41%。

2)簡化豬模型與橢球模型單位長度壓降之間的相對誤差在X方向上為-2.2%～-1.1%，在Z方向上為-11.3%～-9.0%，在Y方向上為14.7%～13.6%。在X、Z和Y方向上的慣性阻力系數(shù)相對誤差分別為-4.0%、-12.2%和14.7%，在今后的研究中可以使用橢球模型代替簡化豬模型。

3)不同躺臥豬模型會影響豬欄內(nèi)風(fēng)速和壓力的分布。橢球模型與簡化豬模型欄內(nèi)的氣流分布基本一致，半橢球與半圓柱模型則對氣流分布產(chǎn)生了顯著影響。壓力分布的差異主要體現(xiàn)在躺臥豬模型分布稀疏的區(qū)域。

4)將通過躺臥豬模型求解得到的阻力系數(shù)用于多孔介質(zhì)模型，能夠得到與模擬單位長度壓降值吻合度較高的驗證值，說明多孔介質(zhì)模型對于簡化動物區(qū)域是適用的。