基于狀態(tài)空間法的大型氣候環(huán)境實驗室熱負荷仿真模型

馬建軍，劉海燕，成竹

（中國飛機強度研究所，西安 710065）

基于狀態(tài)空間法的大型氣候環(huán)境實驗室熱負荷仿真模型

馬建軍，劉海燕，成竹

（中國飛機強度研究所，西安 710065）

目的 建立可用于大型氣候環(huán)境實驗室初期設計熱負荷計算的方法和模型，對實驗室的熱負荷進行計算，為制冷系統(tǒng)選型提供依據(jù)。方法 通過分析實驗室的組成結構和熱負荷來源，基于狀態(tài)空間法建立各部分的熱負荷計算方法，在Simulink中搭建以制冷量為輸入?yún)?shù)的實驗室熱負荷計算仿真模型，對實驗室空載降溫過程進行仿真。結果 該仿真模型可以快速對實驗室熱負荷進行計算，并且可以方便地增減熱負荷模塊，實驗室在降溫過程中熱負荷達3500 kW以上，地板結構的熱負荷占到了總熱負荷50%以上。結論 熱負荷計算結果可應用于空調系統(tǒng)和制冷系統(tǒng)的設計和選型，以及控制策略的優(yōu)化。

大型氣候環(huán)境實驗室；熱負荷；狀態(tài)空間法；Simulink

大型氣候環(huán)境實驗室是以飛機和其他大型裝備為試驗對象，在實驗室內模擬飛機或其他裝備在地面遭受的諸如高溫、低溫、濕熱、日照、降雨、降霧、凍雨、大風、降雪等自然氣候環(huán)境[1]，考核飛機及地面保障設備對這些氣候環(huán)境的適應性，為飛機等大型裝備的環(huán)境適應性設計提供參考。大型氣候實驗室熱負荷計算是制冷/加熱系統(tǒng)、空氣處理系統(tǒng)設計的出發(fā)點，決定了實驗室溫度指標能否實現(xiàn)、系統(tǒng)設計、投資和運行成本，因此準確地計算熱負荷非常重要[2—3]。大型氣候環(huán)境實驗室空間巨大、結構復雜、熱負荷源項多，若以穩(wěn)態(tài)或假設空氣溫度線性變化作為邊界條件進行熱負荷計算，將帶來很大誤差，因此需進行合理的瞬態(tài)熱負荷計算。現(xiàn)有的計算方法如有限差分法[4—5]，編程復雜、適應性較差；CFD流固耦合計算[6—8]則耗時過長。狀態(tài)空間法[9]簡化了許多不必要考慮的環(huán)節(jié)，計算簡便，可在 Simulink中快速搭建模型，得到較準確的總的制冷量需求，并可驗證控制策略等。對于大型氣候環(huán)境實驗室，可建立基于狀態(tài)空間法，并聯(lián)合CFD仿真的熱負荷計算方法。利用CFD仿真計算穩(wěn)態(tài)時，不同氣流組織下圍護結構及試件表面的對流換熱系數(shù)[10]，將對流換熱系數(shù)作為狀態(tài)空間的輸入?yún)?shù)進行動態(tài)仿真，確定不同工況下的總熱負荷，以此設計制冷系統(tǒng)。

1　大型氣候環(huán)境實驗室組成結構

大型氣候實驗室為“房中房”結構，外層為鋼結構建筑，內層為環(huán)境室，如圖1所示。其有效容積超過105m3。環(huán)境室側墻及天花板均由200 mm厚的聚胺脂保溫板拼接而成，地面結構由上到下分為三層：300 mm加強混凝土、300 mm泡沫玻璃、300 mm加強混凝土。

圖1 大型氣候實驗室Fig.1 Large climatic environmental test laboratory

實驗室采用的是間接制冷的方式[11]，制冷系統(tǒng)產生的冷量通過載冷劑系統(tǒng)輸送到空氣處理系統(tǒng)對空氣進行冷卻?？諝馓幚硐到y(tǒng)包括新風系統(tǒng)和循環(huán)風系統(tǒng)[12]，如圖2所示。其中新風系統(tǒng)將室外新風最低處理至-25 ℃，然后送入循環(huán)風系統(tǒng)與室內回風混合并進一步降溫，實驗室內所有的熱負荷將通過循環(huán)風系統(tǒng)的空氣換熱器除去。因此，實驗室空氣處理系統(tǒng)、載冷劑系統(tǒng)、制冷系統(tǒng)應將該熱負荷作為設計輸入?yún)?shù)。

圖2 空氣處理系統(tǒng)Fig.2 Air handling system

2　建立熱負荷仿真模型

大型氣候環(huán)境實驗室的熱負荷主要來源為室內空氣、保溫圍護結構、混凝土地面結構、新風、內部鋼結構、照明燈、循環(huán)風機、試驗件等。由于試驗件多種多樣，難以建立通用的計算模型，暫不考慮，下面建立各個部分的熱負荷計算模型。

2.1 室內空氣降溫放熱

室內空氣溫度的變化是引起保溫圍護結構、地面結構等溫度變化的直接原因，同時結構溫度變化產生熱負荷反過來影響室內空氣溫度的變化。同時新風系統(tǒng)持續(xù)不斷地補充新風以保持實驗室壓力，在每個時刻室內的空氣質量均不同，空氣的溫度變化按式（1）計算。

式中：Ta為空氣溫度，K；Q為制冷量，W；q為圍護結構、地面結構等的熱負荷，W；Va為實驗室內部凈體積，m3。

2.2 保溫圍護結構對流傳熱

保溫圍護結構主要為聚胺脂夾心板，將保溫板分成n層，如圖3所示。

圖3 保溫板分層Fig.3 Layering of insulation panel

取n+1個點溫度作為變量，建立n+1維狀態(tài)變量，空氣溫度當作輸入變量，根據(jù)固體傳熱方程建立常微分方程組[13]。

式中：ci為單位比熱容，ci=δiρicpi；δi為節(jié)點i與節(jié)點i+1之間的厚度，m；ρi為密度，kg/m3；cpi為定壓比熱，J/(kg·K)；Ri為熱阻，Ri=λiδi；λi為導熱系數(shù)，W/(m·K)；Ta，ha分別為室內溫度及表面對流換熱系數(shù)，K，W/(m2·K)；Tb，hb分別為室外溫度及表面對流換熱系數(shù)，K，W/(m2·K)。

將式（2）改寫成：

其中A1是(n+1)×(n+1)維的矩陣。

B1是(n+1)×2的矩陣：

u1是2×1維的輸入列向量：

保溫板向空氣傳遞的熱負荷為：

式中：a1為保溫板的面積。

將式（7）寫成如下形式：

式中：Y=[q1]，C1=[a1ha0…0]，D1= [-a1ha0]。

式（3）和式（8）共同組成了保溫板的狀態(tài)空間[14]。表面對流換熱系數(shù)可通過查相關設計手冊、實測或CFD分析獲得。

2.3 地面結構對流傳熱

同聚胺脂保溫板一樣，將地面結構分層，如圖4所示。上層混凝土分成w層，中層泡沫玻璃分成k層，底層混凝土分成n層。?。╳+k+n+1)個點溫度作為狀態(tài)變量，根據(jù)固體導熱微分方程建立地面結構的狀態(tài)方程。

圖4 地面結構分層Fig.4 Layering of floor structure

與保溫板不同的是，底層混凝土下面是土壤結構，內部有通風管，可以認為其溫度是不變的，節(jié)點（w+k+n+1)溫度微分方程為：

地面向空氣傳遞的熱負荷為：

式中：a2為地面結構的表面積。

構建地面結構的傳熱狀態(tài)空間為：

其中T為(w+k+n+1)維列向量：

A2為(w+k+n+1)×(w+k+n+1)維矩陣，B2為(w+k+n+1)×1維矩陣，C2為(w+k+n+1)×1維矩陣，D2為(w+k+n+1)×1維矩陣，Y=q3。

2.4 新風降溫放熱

降溫過程中，空氣密度變小將導致室內形成負壓，造成室外氣體滲入，損壞保溫結構。因此需持續(xù)不斷地為實驗室補充新風，維持室內保持正壓。這是利用一套新風系統(tǒng)實現(xiàn)的，新風系統(tǒng)以質量流量fm˙補充溫度最低至-25 ℃的新風進入循環(huán)風風道，然后在循環(huán)風處理機組內進一步降溫。新風的熱負荷由式（13）計算。

式中：ta為實驗室內空氣溫度，℃；tf為新風溫度，℃。

2.5 循環(huán)風機做功產熱

循環(huán)風機對空氣的做功為：

式中：P為風機全壓，Pa；V˙為空氣流量，m3/s。

在不同的溫度下，空氣的密度不同，風機的全壓按式（15）修正。

式中：P0為標準工況時風機全壓，Pa；B為當?shù)卮髿鈮毫?，Pa；t0為設計工作溫度，℃；t為工作空氣溫度，℃。

2.6 照明燈產熱

實驗室內的照明燈功率恒為q5。

2.7 內部鋼結構降溫放熱

內部鋼結構為實驗室頂部的承力鋼架，其完全處理空氣對流當中，相對表面積較大，因此可假設其溫度變化率與空氣一致，鋼結構的質量為msteel，其向室內空氣傳熱量為：

2.8 仿真模型搭建

在MATLAB中編寫并計算保溫板和地面結構狀態(tài)空間系數(shù)矩陣，在 Simulink[15]中，依據(jù)圖 5所示的系統(tǒng)原理，搭建各個部件的仿真模塊，最終完成輸入?yún)?shù)為制冷量的實驗室熱負荷計算仿真模型，此處記為模型一。

圖5 模型一 輸入?yún)?shù)為制冷量的熱負荷仿真模型Fig.5 Model 1: Thermal load calculation model with cooling capacity as input

3　應用及結果分析

實驗室要求24 h內從常溫降至-55 ℃，降溫速度不大于 3 ℃/min。初步確定新風流量為 24.6 kg/s，循環(huán)空氣流量為300 m3/s，風機標準環(huán)境下總壓為1800 Pa，經CFD氣流組織仿真得到地面對流換熱系數(shù)為 15 W/(m2·K)，保溫結構表面為 10 W/(m2·K)。照明功率為40 kW，內部鋼結構質量為140 t。

假設制冷量隨室內空氣溫度變化呈線性變化，初始溫度+35 ℃時制冷量為4000 kW，-55 ℃時制冷量減少至3000 kW。降溫過程中各個部分的熱負荷見圖6，地面結構和空氣溫度見圖7，保溫板結構溫度見圖8。

從圖 6—8中可以看出，24 h空氣溫度降至-55 ℃?？諝饨禍厮俾氏瓤旌舐?，4 h左右后溫度降至-25 ℃，對應的熱負荷迅速降低。地面結構熱負荷先升后降，最高達2200 kW，降溫結束時仍達1500 kW，占總熱負荷的接近 1/2。地面結構的上層混凝土被凍透，中層泡沫玻璃保溫性能良好，其底面溫度幾乎未發(fā)生改變。保溫板表層是0.8 mm的不銹鋼，傳熱性能良好，因此其表層溫度與空氣溫度幾乎保持一致，且保溫板保溫性能良好，最外層溫度僅有略微下降。

圖6 模型一 各模塊熱負荷曲線Fig.6 Model 1: Thermal load of different module

圖7　模型一 地面結構不同深度溫度曲線Fig.7 Model 1: Temperature of floor structure in different depth

圖8 模型一 保溫板不同深度溫度曲線Fig.8 Model 1: Temperature of insulation panel in different depth

將圖5所示的原理將模型稍做調整，就可變?yōu)槌Ｓ玫囊允覂瓤諝鉁囟惹€為輸入?yún)?shù)的熱負荷計算模型，記為模型二。

圖9 模型二 輸入?yún)?shù)為室內空氣溫度曲線的熱負荷仿真模型Fig.9 Model 2: Thermal load calculation model with air temperature in door as input

通常假設降溫過程中空氣溫度線性變化：空氣初始溫度為+35 ℃，經24 h線性降至-55 ℃。降溫過程中各模塊的熱負荷見圖 10，地面結構溫度見圖11，保溫板溫度曲線見圖12。

圖10 模型二 各模塊熱負荷曲線Fig.10 Model 2: Thermal load of different module

圖11 模型二 地面結構不同深度溫度曲線Fig.11 Model 2: Temperature of floor structure in different depth

圖12 模型二 保溫板不同深度溫度曲線Fig.12 Model 2: Temperature of insulation panel in different depth

從圖10—12中可以看出，總熱負荷和各分項熱負荷均是逐漸遞增的?？偀嶝摵蓮?50 kW逐漸遞增至3800 kW，地板結構熱負荷從初始的0 kW增至2000 kW，超過總熱負荷1/2。兩種計算模型的結果對比見表1。

表1 兩種模型計算結果對比Table 1 Results of the two thermal load calculation models

制冷系統(tǒng)溫度越低，其COP值越低，若以模型二的計算結果作為制冷系統(tǒng)的設計輸入?yún)?shù)，將造成機組選型偏大，增大初投資和運行成本。同時降溫結束后，模型一得到的地表溫度為-33 ℃，比模型二的-24 ℃低，這將更有利于試驗區(qū)的空氣溫度均勻性。

4　結語

文中采用狀態(tài)空間法建立的，以制冷量為輸入?yún)?shù)的大型氣候環(huán)境實驗室熱負荷仿真模型可以快速、準確地計算瞬態(tài)熱負荷，模塊化的設計使模型容易搭建和修改，且可以很容易調整為以空氣溫度曲線為輸入?yún)?shù)的熱負荷仿真模型。降溫過程若假設空氣溫度線性變化，計算得到的熱負荷在低溫段偏大，以此為設計參數(shù)將導致制冷系統(tǒng)選型偏大，同時地表溫度降低緩慢，不利于試驗區(qū)的空氣溫度均勻性。

實驗室建造中，將在地板結構、保溫板結構等中預埋溫度傳感器，根據(jù)實測結果對仿真模型進行修正，使仿真結果更精確?？筛鶕?jù)具體試驗工況編寫相應的熱負荷模型，如降雪、吹風、太陽輻照等，實時預測熱負荷，為制冷系統(tǒng)運行和控制系統(tǒng)控制策略優(yōu)化提供參數(shù)。

[1] 唐虎, 李喜明. 飛機氣候試驗[J]. 裝備環(huán)境工程, 2012, 9(1): 60—65. TANG Hu, LI Xi-ming. Climatic Test of Aircraft[J]. Equipment Environment Engineering, 2012, 9(1): 60—65.

[2] 王浚, 黃本誠, 萬才大, 等. 環(huán)境模擬技術[M]. 北京:國防工業(yè)出版社, 1996. WANG Jun, HUANG Ben-cheng, WAND Cai-da, et al. Environment Simulation Technology[M]. Beijing: National Defense Industry Press, 1996.

[3] 李兆堅. 大型綜合氣候環(huán)境室空調制冷設計[J]. 暖通空調, 1999, 29(5): 57—59. LI Zhao-jian. Air Conditioning and Refrigeration Design of a Large Multi-purpose Environmental Test Chamber[J]. HV&AC, 1999, 29(5): 57—59.

[4] 李兆堅. 環(huán)境實驗室圍護結構不穩(wěn)定傳熱特性分析[J].環(huán)境技術,1992(2): 22—44. LI Zhao-jian. Transit Heat Transfer Analysis of Climatic Test Chamber Insulation Layer[J]. Environmental Techonlogy, 1992(2): 22—24.

[5] 吳志勇. 環(huán)境試驗箱圍護結構傳熱研究[J]. 裝備環(huán)境工程, 2007, 4(5): 51—56. WU Zhi-yong. Study on Heat Transfer of Exterior Pretected Construction of Environmental Test Chamber[J]. Equipment Environmental Engineering, 2007, 4(5): 51—56.

[6] 王怡, 劉國華, 陳靖. 建筑室內環(huán)境數(shù)值模擬中的耦合方法[J]. 建筑科學與工程學報, 2006, 23(4): 79—83. WANG Yi, LIU Guo-hua, CHEN Jing. Coupling Method by Numerical Simulation of Building Indoor Environment [J]. Journal of Architecture and Civil Engineering, 2006, 23(4): 79—83.

[7] ZHAI Zhi-qiang, CHEN Qing-yan, et al. On Approaches to Couple Energy Simulation and Computational Fluid Dynamics Programs[J]. Building and Environment, 2002 (37): 857—864.

[8] 劉海燕, 馬建軍. 環(huán)境實驗室氣流組織與圍護結構耦合傳熱研究[J]. 裝備環(huán)境工程, 2016, 13(2): 44—51. LIU Han-yan, MA Jian-jun. Coupled Heat Transfer Analysis of Airflow and Enclosure Structures of Environment Test Chamber[J]. Equipment Environmental Engineering, 2016, 13(2): 44—51.

[9] 翟喚春, 尚娜, 王浚. 環(huán)境試驗設備的狀態(tài)空間法仿真設計[J]. 北京航空航天大學學報, 2004, 30(2): 109—113.ZHAI Huan-chun, SHANG Na, WANG Jun. Simulation Desing Based on State-space for Environment Test Equipment[J]. Journal of Beijing University of Aeronautics and Astronautics, 2004, 30(2): 109—113.

[10] ZHAI Zhi-qiang, CHEN Qing-yan. Numerical Determination and Treatment of Convective Heat Transfer Coefficient in the Coupled Building Energy and CFD Simulation[J]. Building and Environment, 2004 (39): 1001—1009.

[11] 張惠, 劉海燕, 李喜明, 等. 大型氣候環(huán)境實驗室冷媒的分析選擇[J]. 裝備環(huán)境工程, 2015, 12(2): 104—109. ZHANG Hui, LIU Hai-yan, LI Xi-ming, et al. Analysis and Selection of Refrigerants Used In the Large Climatic Environment Test Facility[J]. Equipment Environmental Engineering, 2015, 12(2): 104—109.

[12] 劉海燕, 馬建軍, 張惠. 大型氣候環(huán)境實驗室空氣處理系統(tǒng)方案討論[J]. 裝備環(huán)境工程, 2014, 11(5): 107—113. LIU Hai-yan, MA jian-jun, ZHANG Hui. Discussionon Design of the Air Handing System in Large Climatic Environ-mentalTestLaboratory[J]. Equipment Environment Engineering, 2014, 11(5): 107—113.

[13] 陶文銓. 數(shù)值傳熱學[M]. 第二版. 西安: 西安交通大學出版社, 2001. TAO Wen-quan. Numerical Heat Transfer[M]. 2nd Edition. Xi’an: Xi’an Jiaotong University Press, 2001.

[14] 胡壽松. 自動控制原理[M]. 北京: 科學出版社, 2013: 419—511. HU Shou-song. Principle of automatic control[M]. Beijing: Science Press, 2013: 419—511.

[15] 王正林, 郭陽寬. MATLAB/SIMULINK與過程控制系統(tǒng)仿真[M]. 北京: 電子工業(yè)出版社, 2012. WANG Zheng-lin, GUO Yang-kuan. MATLAB/ SIMULINK and Simulation of Process Control System [M]. Beijing: Electronic Industry Press, 2012.

A Climatic Test Laboratory Cooling Capacity Calculation Model Based on Space State Method

MA Jian-jun, LIU Hai-yan, CHENG Zhu
(Aircraft Strength Research Institute of China, Xi’an 710065, China)

Objective To establish a thermal-load calculation model applied in large climatic test laboratory, the calculation result will help the design of cooling system. Methods First, every potential thermal source was analyzed. Then thermal calculation method and simulation module for thermal source was founded in Simulink. At last a thermal-load simulation models of cooling capacity dominant was set up and applied in the thermal calculation of chamber cooling process. Results The thermal load calculation models is time efficiency for thermal calculation of large climatic test laboratory and convenient in modules update. Thermal load of large climatic test laboratory during cooling down is larger than 3500 kW, thermal load of the floor is more than 50% of the total thermal load. Conclusion The calculation result will help the design of air-condition system and cooling system, and the optimization of control strategy.

large climatic test laboratory; thermal load; state-space; Simulink

2016-09-28；Revised：2016-09-29

10.7643/ issn.1672-9242.2016.05.019

TJ01

A

1672-9242(2016)05-0115-07

2016-09-28；

2016-09-29

中國航空工業(yè)集團公司技術創(chuàng)新基金（2013F62302）

Fund：Technology Innovation Fund of Aviation Industry Corporation of China（2013F62302）

馬建軍（1989—），男，安徽阜陽人，碩士，工程師，主要研究方向為飛機氣候環(huán)境試驗技術。

Biography：MA Jian-jun（1989—）， Male， from Fuyang， Anhui， Master， Engineer， Research focus: aircraft climatic environmental test technology.