基于CFD的玻璃溫室環(huán)境數(shù)值模擬及優(yōu)化分析

賈鶴鳴，韓駿騁，張森，孫康健，李瑤

東北林業(yè)大學(xué) 機(jī)電工程學(xué)院，黑龍江 哈爾濱 150040

溫室植物的產(chǎn)量取決于溫室內(nèi)環(huán)境的溫度、空氣濕度等小氣候要素的適宜度[1]，尤其是較好的溫度和氣流分布，對(duì)于植物生長(zhǎng)有著極其重要的作用。在所有環(huán)境因子中，適宜的溫度是植物生長(zhǎng)的首要條件，過(guò)低的溫度甚至?xí)斐蓽厥覂?nèi)植物的死亡。因此，溫室通風(fēng)口的設(shè)置和風(fēng)速都會(huì)對(duì)溫室內(nèi)氣流的分布情況有較大的影響。程秀花等[1]提出室內(nèi)氣流的速度和分布情況較易受到植物的影響；汪鄭邦[2]認(rèn)為植物溫室的氣流分布對(duì)作物生長(zhǎng)有較大的影響，不均勻的氣流分布容易在作物周圍形成換流，影響作物的蒸騰和光合作用；劉文合等[3]的研究表明通過(guò)對(duì)溫室進(jìn)行升溫，可以使得溫室內(nèi)氣溫顯著升高，但對(duì)溫室內(nèi)土壤深處的溫度影響不大；周偉等[4]通過(guò)試驗(yàn)證明了通風(fēng)口的設(shè)置和風(fēng)速都會(huì)對(duì)溫室內(nèi)氣流的分布情況有較大的影響，同時(shí)溫室的氣流分布對(duì)植物蒸騰和光合作用有較大的影響。因此，如何合理優(yōu)化溫室結(jié)構(gòu)，得到適宜的風(fēng)速和溫度并使兩者均勻性提高是當(dāng)前急需解決的重要問(wèn)題。

1 材料和方法

1.1 溫室模型

試驗(yàn)在東北林業(yè)大學(xué)進(jìn)行，模擬溫室建立在室內(nèi)，受外界環(huán)境影響因素較小。溫室的長(zhǎng)×寬×高為2 m×0.83 m×1.8 m，由于溫室是密閉的，因此忽略溫室的門，并對(duì)溫室四周及頂部進(jìn)行簡(jiǎn)化，即將溫室的頂部和四周均視為以浮法玻璃為材料?？紤]到植物的冠層阻力，將植物模型簡(jiǎn)化為長(zhǎng)1.50 m、寬0.6 m、高0.5 m的多孔介質(zhì)，置于溫室正中央。暖氣扇規(guī)格為0.4 m×0.4 m，置于溫室西側(cè)，離地面1.3 m高處。采用Gambit進(jìn)行以溫室西北角為零點(diǎn)的溫室模型的建立和網(wǎng)格化，得到的玻璃溫室模型如圖1所示。

圖1 玻璃溫室模型

1.2 監(jiān)測(cè)點(diǎn)布置

在溫室中選取Z=0.415 m處的截面，均勻分布6個(gè)監(jiān)測(cè)點(diǎn)，高度分別為Y=0.6 m和1.2 m，將監(jiān)測(cè)點(diǎn)依次記為P1,P2,···,P6。在平面X=1 m處，共建立4個(gè)監(jiān)測(cè)點(diǎn)，高度分別為0.6 m和1.2 m，記為P7、P8、P9、P10。通過(guò)將在監(jiān)測(cè)點(diǎn)處得到的實(shí)際值與模擬值進(jìn)行對(duì)比驗(yàn)證。采用風(fēng)速儀和溫度計(jì)對(duì)溫室內(nèi)風(fēng)速和溫度進(jìn)行測(cè)量，溫室環(huán)境模擬值由Fluent軟件進(jìn)行模擬，監(jiān)測(cè)點(diǎn)的布置如圖2、3所示。

圖2 Z=0.415 m截面監(jiān)測(cè)點(diǎn)布置（Y=0.6 m和Y=1.2 m）

圖3 Z=0.415 m截面監(jiān)測(cè)點(diǎn)布置（X=1 m）

1.3 數(shù)學(xué)模型

1.3.1 CFD模型

當(dāng)Fluent模擬溫室環(huán)境時(shí)，使用標(biāo)準(zhǔn)k-ε模型[5]將溫室內(nèi)的氣流視為遵循質(zhì)量守恒方程、動(dòng)量守恒方程和能量守恒方程的湍流。

動(dòng)量守恒方程指出，所有的流動(dòng)在實(shí)際情況下都應(yīng)該遵循動(dòng)量守恒定律，定律用公式表示為：

式中的u、v、w分別為在x、y、z方向上溫室內(nèi)氣流流速的分量

動(dòng)量守恒用公式表現(xiàn)為：

式中：μeあ代表有效黏度；β代表空氣熱膨脹系數(shù)；g代 表重力加速度；T代表空氣溫度；Tref代表空氣的參考溫度。

能量守恒用公式表示為

式中：q代 表熱源；代表比熱容； λeあ代表有效導(dǎo)熱系數(shù)[6]。

在處理封閉室內(nèi)因溫差引起的浮升力項(xiàng)時(shí)采用Boussinesq假設(shè)[7]。Boussinesq假設(shè)是對(duì)流體流動(dòng)現(xiàn)象的一種簡(jiǎn)化，假設(shè)認(rèn)為：流體中的黏性耗散忽略不計(jì)；除密度外其他物性為常數(shù)；對(duì)密度僅考慮動(dòng)量方程中與體積力有關(guān)的項(xiàng)，其余各項(xiàng)中的密度亦作為常數(shù)。重力項(xiàng)中的密度如若采用作參考溫度，則在不同溫度T下溫室內(nèi)部的空氣密度ρ可以用數(shù)學(xué)方程表示為：

式中：β是溫室內(nèi)氣流的熱膨脹系數(shù)；T為溫室內(nèi)的溫度；T0為溫室內(nèi)部空氣的參考溫度。

流經(jīng)多孔介質(zhì)的空氣可通過(guò)Darcy-Forchheimer理論描述為

式中：Sφ代表動(dòng)量源項(xiàng)；Kp代表多孔介質(zhì)的滲透性系數(shù)；CF代表非線性動(dòng)量損失系數(shù)；μ代表空氣的動(dòng)力黏度；ρ代表空氣密度；u代表空氣流速[8]。

將植物視為各向同性的多孔介質(zhì)[9]，并建立其與溫室內(nèi)部風(fēng)速的數(shù)學(xué)關(guān)系模型：

式中：ILAV代 表葉面積指數(shù)；CD是植物冠層的阻力系數(shù)。

對(duì)于多孔介質(zhì)動(dòng)量源項(xiàng)，使用冪律模型，自定義經(jīng)驗(yàn)系數(shù)C0為0.395，C1為0.2。植物通過(guò)呼吸作用產(chǎn)生的熱量與溫室中氣流進(jìn)行熱交換，將植物通過(guò)呼吸作用產(chǎn)生的熱量視為常值，在進(jìn)行數(shù)值模擬時(shí)將植物區(qū)域看作“體積熱源”。

1.3.2 計(jì)算域與網(wǎng)格劃分

在使用Gambit劃分網(wǎng)格時(shí)將整個(gè)玻璃溫室作為計(jì)算域，采用非結(jié)構(gòu)化四面體網(wǎng)格(混合TGrid網(wǎng)格類型以及Tet/Hybrid混合單元)將溫室整體劃分為四面體混合網(wǎng)格單元[10]，并局部加密暖風(fēng)扇出風(fēng)口和植物區(qū)域，以適應(yīng)氣流和溫度可能出現(xiàn)的較大變化[11]。

1.3.3 邊界條件

因?yàn)闇厥颐荛]性較好且置于室內(nèi)，不考慮太陽(yáng)輻射和漏風(fēng)的影響，對(duì)模型進(jìn)行簡(jiǎn)化。將暖風(fēng)扇側(cè)的邊界條件設(shè)為速度入口條件，吹出水平方向25 ℃暖風(fēng)，平均速度為1 m/s；覆蓋層全部簡(jiǎn)化為玻璃，設(shè)為壁面條件，其與整個(gè)溫室進(jìn)行能量交換的方式主要為對(duì)流傳熱和傳導(dǎo)[10]，初始溫度設(shè)為20 ℃；為簡(jiǎn)化能量交換過(guò)程，在進(jìn)行模擬時(shí)不考慮其他的傳熱過(guò)程。將擁有大的熱容量的土壤的邊界條件設(shè)為壁面[12]，初始溫度設(shè)為20 ℃。

1.3.4 求解運(yùn)算

在進(jìn)行模擬時(shí)，為修正速度，選用SIMPLEC算法：梯度項(xiàng)選用最小二乘法；壓力項(xiàng)采用標(biāo)準(zhǔn)算法；能量和動(dòng)量采用二階迎風(fēng)；為加快收斂速度，黏性和湍流選用一階迎風(fēng)。

2 FLUENT仿真結(jié)果與驗(yàn)證

2.1 FLUENT仿真結(jié)果

用Fluent進(jìn)行模擬的結(jié)果如圖4、5所示。在吹入熱風(fēng)后，暖風(fēng)口周圍區(qū)域的氣體流速變化明顯，且存在較大梯度。氣流從暖風(fēng)扇水平送入溫室后，略微向下方運(yùn)動(dòng)，在溫室內(nèi)部的東面、暖風(fēng)扇正下方和溫室底部形成渦流，且溫度與風(fēng)速有著直接的關(guān)系。

圖5 Y=0.5 m和Y=0.1 m處溫度云圖

選取溫室內(nèi)植物冠層平面Y=0.5 m和植物內(nèi)部平面Y=0.1 m的速度云圖（圖4）進(jìn)行分析。當(dāng)風(fēng)速為0.30～1.00 m/s時(shí)，植物生長(zhǎng)較好[13]。在使用Fluent進(jìn)行模擬計(jì)算后，在距離地面0.5 m處的平均空氣流速為0.152 m/s，最低速度為0.009 m/s，低于最佳風(fēng)速；而植物冠層部分沒(méi)有處于最適風(fēng)速的區(qū)域出現(xiàn)。在平面Y=0.1 m處氣流平均流速為0.195 m/s，最小流速為0.013 m/s，低于適宜風(fēng)速的最小值0.3 m/s，最大風(fēng)速為0.552 m/s，雖略大于0.3 m/s，但植物內(nèi)部達(dá)到適宜植物生長(zhǎng)的風(fēng)速的部分僅占整體的20%。綜上所述，此時(shí)的氣體流動(dòng)不適于植物的生長(zhǎng)。

選取溫室內(nèi)植物表面Y=0.5 m和植物內(nèi)部Y=0.1 m的溫度云圖（圖5）進(jìn)行分析。在暖風(fēng)口附近溫度最高，溫度與其他區(qū)域相比有很大的梯度。在2個(gè)平面內(nèi)最高溫度可達(dá)24.45 ℃，平均溫度23.17 ℃，溫差的最大值為1.28 ℃，此溫度較適于大部分植物的生長(zhǎng)。

2.2 對(duì)比驗(yàn)證

在試驗(yàn)期間比較10個(gè)監(jiān)測(cè)點(diǎn)的實(shí)際測(cè)量值和使用Fluent模擬得出的計(jì)算值，其結(jié)果如圖6、7所示。

圖6 風(fēng)速在監(jiān)測(cè)點(diǎn)模擬值與實(shí)測(cè)值的比較

圖7 溫度在監(jiān)測(cè)點(diǎn)模擬值與實(shí)測(cè)值比較

根據(jù)圖6，對(duì)比不同監(jiān)測(cè)點(diǎn)的風(fēng)速，所有觀測(cè)點(diǎn)的絕對(duì)誤差的平均值為0.008 m/s，相對(duì)誤差的平均值為5.54%，絕對(duì)誤差和相對(duì)誤差的最大值出現(xiàn)在監(jiān)測(cè)點(diǎn)P6處，分別為0.027 m/s和24.228%。雖然溫室內(nèi)風(fēng)速的模擬值小于測(cè)量值，但大部分情況下，模擬結(jié)果與實(shí)際情況一致。因此可以認(rèn)為模型的速度計(jì)算值可以較好地體現(xiàn)溫室內(nèi)部氣體流速的情況。

根據(jù)圖7，對(duì)比不同監(jiān)測(cè)點(diǎn)的溫度，10個(gè)監(jiān)測(cè)點(diǎn)絕對(duì)誤差均值為0.44 ℃，相對(duì)誤差均值為4%，P10處出現(xiàn)最大絕對(duì)誤差1.66 ℃和最大相對(duì)誤差2.78%，模擬值與從監(jiān)測(cè)點(diǎn)測(cè)出的值差別不大，出現(xiàn)誤差的原因可能是溫室模型密封不夠好而造成漏風(fēng)，導(dǎo)致室溫略低。結(jié)合儀器位置和人員干擾對(duì)實(shí)驗(yàn)結(jié)果造成的影響[14-15]，從整體來(lái)看，溫度的實(shí)際值和計(jì)算值誤差很小，因此可以將溫度的模擬值視作溫室內(nèi)溫度的實(shí)際值。

3 溫室內(nèi)環(huán)境的優(yōu)化設(shè)計(jì)

3.1 優(yōu)化方案設(shè)計(jì)

溫室內(nèi)的速度場(chǎng)、溫度場(chǎng)的均勻性受暖風(fēng)扇的影響較大，因此通過(guò)在溫室中改變進(jìn)風(fēng)口的速度和改變暖風(fēng)口的位置來(lái)設(shè)計(jì)出3種改進(jìn)方案，以滿足植物生長(zhǎng)的需要。方案一，升高暖風(fēng)口速度，即暖風(fēng)口速度從1 m/s升高到2 m/s；方案二，下移暖風(fēng)口的位置，將暖風(fēng)口從1.3 m處下移到1 m處，如圖8所示；方案三，下移暖風(fēng)口的位置至0.7 m處，如圖9所示。將原溫室作為方案零來(lái)進(jìn)行對(duì)比，溫室、植物和暖風(fēng)扇的尺寸視為定量不做改變，且方案一、方案二和方案三的網(wǎng)格劃分與方案零相同，其余邊界條件的設(shè)置也相同。通過(guò)對(duì)不同方案使用Fluent進(jìn)行模擬計(jì)算，比較植物冠層表面和內(nèi)部的風(fēng)速及溫度分布情況，依此來(lái)選擇擁有適宜風(fēng)速溫度的最適宜植物生長(zhǎng)的方案。

圖8 暖風(fēng)扇離地1 m時(shí)的玻璃溫室模型

圖9 暖風(fēng)扇離地0.7 m時(shí)的玻璃溫室模型

3.2 優(yōu)化模擬分析

3.2.1 速度場(chǎng)模擬分析

經(jīng)過(guò)改良后的方案氣流分布模式基本相同：在溫室上部和底部都有渦流的產(chǎn)生，暖風(fēng)扇口附近的氣流流速最大，并隨著流動(dòng)方向減弱，直到在東側(cè)壁面附近開(kāi)始沿溫室頂部返回，溫室西側(cè)的植物表面部分和溫室壁面附近風(fēng)速較小，且會(huì)在植物內(nèi)部的渦流處出現(xiàn)整個(gè)植物區(qū)域中風(fēng)速的最大值。

在方案一中，從圖10速度云圖上看，選取的2個(gè)平面空氣流速平均值為0.695 m/s，均高于方案零的0.152 m/s，且滿足植物生長(zhǎng)的最佳風(fēng)速。風(fēng)速在南北兩側(cè)墻處達(dá)到最大，為1.562 m/s，遠(yuǎn)高于最適宜植物生長(zhǎng)的風(fēng)速1 m/s，這樣的速度容易對(duì)植物造成損傷，此時(shí)植物表面風(fēng)速位于0.3～1.0 m/s的部分約占整個(gè)植物冠層的63%。在植物內(nèi)部Y=0.1 m處最大風(fēng)度為1.581 m/s，平均風(fēng)速0.928 m/s，處在最佳氣流流速范圍內(nèi)的區(qū)域約占27.9%，剩下的大部分區(qū)域風(fēng)速均大于1 m/s，不利于植物生長(zhǎng)，由此可知不能僅增大風(fēng)速來(lái)優(yōu)化溫室內(nèi)部的氣流分布情況。

圖10 方案一速度場(chǎng)模擬

從圖11方案二的云圖可知，y=0.5 m處氣流的流速均值為0.284 m/s，最大為0.632 m/s；y=0.1 m處平均風(fēng)速為0.423 m/s，最小風(fēng)速為0 m/s，最大風(fēng)速為0.754 m/s，均小于方案一。處于適宜風(fēng)速的占25%，雖然相對(duì)于方案零適宜風(fēng)速的占比有所增加，氣流的均勻性與方案一相比較好，但適宜風(fēng)速的占比仍小于方案一。

圖11 方案二速度場(chǎng)模擬

分析圖12方案三的云圖，選取分析的平面平均氣體流速為0.383 m/s，最高流速可達(dá)0.834 m/s，均符合最佳風(fēng)速。在植物部分的最低風(fēng)速為0.07 m/s，大于方案零和方案二，但仍低于0.3 m/s。這時(shí)的植物冠層風(fēng)速達(dá)到適宜風(fēng)速的占70%，高于其他方案，且均勻性最好。

圖12 方案三速度場(chǎng)模擬

3.2.2 溫度場(chǎng)模擬分析

分析方案一、方案二和方案三的溫度情況（圖13～15），并將其與速度云圖進(jìn)行比對(duì)，發(fā)現(xiàn)氣流速度較大的區(qū)域溫度上升較多。

圖13 方案一溫度場(chǎng)模擬

圖14 方案二溫度場(chǎng)模擬

圖15 方案三溫度場(chǎng)模擬

3個(gè)方案的最高溫度分別為24.49、24.45和24.72 ℃，平均溫度為 23.93、23.15 和 23.49 ℃，均高于方案零，且植物區(qū)域的溫度分布均勻性優(yōu)于方案零。就溫度而言，方案一的均值和最大值都為最高，其次為方案三。

通過(guò)模擬分析實(shí)驗(yàn)中設(shè)計(jì)的3種溫室的優(yōu)化方案，可以發(fā)現(xiàn)：1）方案一的氣流均勻性稍差，且有較多的區(qū)域風(fēng)速高于最佳風(fēng)速，不利于植物生長(zhǎng)；2）方案二雖然略優(yōu)于方案零，但對(duì)于氣流場(chǎng)和溫度場(chǎng)并不能提供有效的優(yōu)化；3）方案三的氣流均勻性最好，不存在流速大于1 m/s的區(qū)域，從而不會(huì)對(duì)植物的正常生長(zhǎng)造成威脅，且風(fēng)速大于方案零和方案二?？傮w來(lái)講，方案三的暖風(fēng)扇位置更適宜植物的生長(zhǎng)。

4 結(jié)論

1）溫度數(shù)值模擬結(jié)果與實(shí)驗(yàn)實(shí)測(cè)結(jié)果吻合良好，表明建立的CFD模型有效，模擬玻璃溫室內(nèi)環(huán)境分布是可行的；

2）設(shè)計(jì)的3種結(jié)構(gòu)優(yōu)化方案的風(fēng)速值及均勻性均優(yōu)于方案零。綜合考慮，方案三氣流循環(huán)模式最佳，氣流均勻性最好，適宜植物生長(zhǎng)。

本文僅基于風(fēng)速和暖風(fēng)扇位置進(jìn)行溫室的優(yōu)化，在之后的研究中，可考慮將添加排氣扇和增加暖風(fēng)扇的數(shù)量作為溫室的優(yōu)化方式，考慮植物在陽(yáng)光下的光合作用產(chǎn)生的熱量和外部輻射等會(huì)影響溫室內(nèi)部氣流和溫度的條件，探究溫室內(nèi)部氣流與溫度的分布情況。