1、精選優(yōu)質文檔-傾情為你奉上全國第X屆研究生數(shù)學建模競賽題 目 恒溫箱溫度變化的數(shù)學建模摘 要恒溫箱是航空、汽車、家電、科研等領域必備的測試設備，用于測試和確定電工、電子及其他產品及材料進行高溫試驗的溫度環(huán)境變化后的參數(shù)及性能。因此,對溫度進行測量和控制也是科學實驗和工業(yè)生產中經常需要解決的重要問題。隨著計算機技術和控制理論的發(fā)展,以及對產品質量要求的提高,人們對高精度的溫度測量和控制的要求也越來越高。電加熱設備溫度特性復雜,其溫度的測量和控制亦顯得尤為重要和復雜。因此，本問題具有很強的研究背景和實際應用價值。 為描述保溫箱溫度的變化規(guī)律，文章首先對環(huán)境溫度進行三次樣條插值，將環(huán)境溫度的采樣間隔

2、變成和恒溫箱一樣，也就是一分鐘。然后根據(jù)能量平衡原理和熱力學知識列出恒溫箱和隔熱層從室溫開始加熱這一階段的熱量平衡方程組，方程組一共有6個未知量，取78個分鐘的數(shù)據(jù)帶入方程組中得到156個方程式，用matlab求解這個超定方程組，求得6個未知量的值，也就是恒溫箱和隔熱層的各項參數(shù)值。隨后將剩余兩個分鐘的數(shù)據(jù)帶入求得的方程組中進行檢驗，誤差在接受范圍內，說明建模合理。接下來，將熱量平衡代數(shù)方程組變?yōu)楣β势胶馕⒎址匠探M，在考慮恒溫箱控制條件的情況下用龍格庫塔法進行求解，得到接下來恒溫箱和隔熱層的溫度變換規(guī)律，畫出溫度變換曲線。最后對恒溫箱的性能進行了評價。關鍵詞：恒溫箱 隔熱層 熱學平衡方程式 超

3、定方程組 龍格庫塔法參賽密碼 （由組委會填寫）參賽隊號 隊員姓名 陳文哲 管俊 孫鵬偉 XX大學承辦 恒溫箱溫度變化的數(shù)學建模一、問題重述恒溫箱是航空、汽車、家電、科研等領域必備的測試設備，用于測試和確定電工、電子及其他產品及材料進行高溫試驗的溫度環(huán)境變化后的參數(shù)及性能。因此,對溫度進行測量和控制也是科學實驗和工業(yè)生產中經常需要解決的重要問題。隨著計算機技術和控制理論的發(fā)展,以及對產品質量要求的提高,人們對高精度的溫度測量和控制的要求也越來越高。電加熱設備溫度特性復雜,其溫度的測量和控制亦顯得尤為重要和復雜。因此，本問題具有很強的研究背景和實際應用價值。 模型要求熱源即加熱器以恒定的速率加熱恒

4、溫箱。恒溫箱體內裝有一個自動的溫度控制器，在溫度高于68.2華氏度會自動關閉加熱系統(tǒng), 在溫度低于67.8華氏度時會自動打開加熱系統(tǒng)，但溫度控制器讀取溫度有時間間隔，設定溫度控制器讀取溫度有時間間隔為每1分鐘1次。 根據(jù)上述分析，本文擬解決如下問題：（1）、根據(jù)已有的資料數(shù)據(jù)，利用所學知識及相關參考資料，通過研究溫度傳播規(guī)律,確定加熱方案,并建立恒溫箱系統(tǒng)的熱數(shù)學模型。；（2）、建立恒溫箱體的溫度變化模型。并根據(jù)小型試驗數(shù)據(jù)及相關數(shù)學、物理方法，驗證模型準確度和適用性；（3）、對恒溫箱產品的質量優(yōu)劣做出評價；（4）、對建立的模型進行誤差分析，通過模型參數(shù)調節(jié)優(yōu)化模型效果，并根據(jù)修改模型預測最有

5、效數(shù)學模型。根據(jù)誤差原因分析等對未來需要進行的試驗和研究工作提出了一些建議；（5）、根據(jù)相似準則設計相關實驗，對已建立模型進行推廣和擴充。以達到理論-實驗-工程三者的互相驗證。 二、基本假設（1）、環(huán)境溫度不受熱源、恒溫箱和隔熱層的影響；（2）、忽略空氣與恒溫箱及隔熱層之間的對流換熱、熱輻射；（3）、所有的導熱過程僅考慮穩(wěn)態(tài)導熱；（4）、所有物體的物性不隨溫度變化；（5）、環(huán)境溫度始終低于67.8華氏度；（6）、恒溫箱體與箱頂隔熱層溫度始終不低于外界環(huán)境溫度；三、符號約定 恒溫箱比熱容 恒溫箱質量 時間內恒溫箱溫度增量 恒溫箱箱體熱導率 恒溫箱某一時刻溫度 環(huán)境溫度 恒溫箱與外界環(huán)境的接觸表面

6、積 恒溫箱體的厚度 隔熱層某時刻溫度 恒溫箱體和隔熱層的接觸表面積 熱源時間內所釋放的熱量 隔熱層的比熱容 隔熱層的質量 時間隔熱層溫度的增量 隔熱層的熱導率 隔熱層的厚度專心-專注-專業(yè)四、模型的建立與求解基本理論（1）熱傳導只要有溫差,熱量就會自發(fā)地從高溫物體傳到低溫物體,所以傳熱是一種最常見的物理過程。實際上,熱是一種擴散效應,物體將能量通過熱的方式釋放出來,因此熱其實是能量的傳遞。熱的傳遞方式主要包括以下三種:熱傳導、熱對流和熱輻射。在對恒溫箱的數(shù)學建模中，我們忽略熱對流和熱輻射，僅考慮熱傳導的主要作用。熱傳導是指當物體內具有溫差或者不同溫度的物體相互接觸時,在各部分之間沒有發(fā)生相對宏

7、觀位移情況下,通過物質微觀粒子(分子、原子或自由電子)的熱運動而進行的熱量傳遞。導熱是固體內唯一的熱量傳遞方式;在氣體和液體中,盡管也同樣存在著導熱現(xiàn)象。圖2 恒溫箱壁面導熱示意圖如圖所示恒溫箱壁，一塊厚度為，橫截面積為的平板，兩側表面分別維持著均勻溫度和。實驗表明：單位時間內從表面1 傳遞到表面2的熱量與壁面間的溫差以及導熱面積成正比，而與平板的厚度成反比，即方程中的比例系數(shù)反映了材料導熱能力的大小，叫做導熱系數(shù)或者熱導率，單位是。導熱率是一個物性參數(shù)，不同材料的導熱率不同，即使是同一種材料，在不同溫度下也具有不同的熱導率。為了確保恒溫箱內的光譜儀能夠安全可靠的工作,我們需要將恒溫箱內多余的

8、熱量排放到外面，從而保證恒溫箱內的光譜儀能夠獲得適宜的工作溫度，這是系統(tǒng)級熱設計。總的來說,恒溫箱的熱設計可以分為3個主要方面。第一：組件級熱設計，即對元器件級的熱設計。第二：封裝級設計，即對于電子控制模塊、散熱設備、PCB板的熱設計。第三：系統(tǒng)級設計，即對于恒溫箱內膽及外殼的熱設計。系統(tǒng)級的熱設計主要研究電子設備所處環(huán)境的溫度對其產生的影響，環(huán)境溫度是系統(tǒng)級熱分析的重要邊界條件，其熱設計是采取措施控制環(huán)境溫度，使電子設備在適宜的溫度環(huán)境下進行工作。隨著電子技術的發(fā)展，尤其是微電子技術的發(fā)展，電子元器件和設備的尺寸正迅速縮小，而功率卻一直在增大，使得單位體積容納的熱量越來越多，特別是在惡劣環(huán)境

9、下電子產品以及國防電子儀器要求必須達到良好的電磁兼容性能及三防性能等。（2）三次樣條差值將環(huán)境溫度每隔15min采集到的數(shù)據(jù)利用三次樣條差值成每隔1min的數(shù)據(jù)。算法原理：設在區(qū)間上給定個節(jié)點，在節(jié)點處的函數(shù)值為。若函數(shù)滿足如下三條：（1） 在每個子區(qū)間上，是三次多項式；（2） ；（3） 在區(qū)間上，的二階導數(shù)連續(xù)。則稱為函數(shù)在區(qū)間上的三次樣條插值函數(shù)。子區(qū)間上的的表達式為：關于參數(shù)的方程組（三彎矩方程組）：牛頓插值多項式：構造牛頓插值多項式首先列出差商表，進而由差商表寫出牛頓插值多項式。階差商為：零階差商和一階差商：（3）龍格庫塔法算法原理：考慮用函數(shù)f(x,y)在若干點上的函數(shù)值的線性組合來

10、構造近似公式，構造時要求近似公式在(xi,yi)處的Taylor展開式與解y(x)在xi處的Taylor展開式的前面幾項重合，從而使近似公式達到所需要的階數(shù)。既避免求高階導數(shù),又提高了計算方法精度的階數(shù)。或者說,在xi,xi+1這一步內多計算幾個點的斜率值，然后將其進行加權平均作為平均斜率，則可構造出更高精度的計算格式，這就是龍格庫塔（Runge-Kutta）法的基本思想。用類似上述的處理方法，只需在區(qū)間xi,xi+1上用四個點處的斜率加權平均作為平均斜率K*的近似值，構成一系列四階龍格庫塔公式。具有四階精度，即局部截斷誤差是O(h5)。對于高階常微分方程的初值問題：將其轉化為一階微分方程組求

11、解：標準的四級四階龍格庫塔法的向量形式：其分量形式：（4）超定方程的最小二乘解小二乘法廣泛地應用于工程計算中，用最小二乘法消除（平滑）誤差，用最小二乘法從有噪聲的數(shù)據(jù)中提取信號，從海量數(shù)據(jù)中找出數(shù)據(jù)變化的趨勢，。甚至利用簡單函數(shù)計算復雜函數(shù)的近似值，我們并不期望它的近似值多么精確（事實上很多時候也不用很精確），盡管如此還是希望計算出的近似數(shù)據(jù)與原始數(shù)據(jù)之間有相似之處。如果從線性代數(shù)角度來理解最小二乘法，實際上是將一個高維空間的向量投影到低維子空間所涉及的工作。超定方程組的最小二乘解當方程組GX=b的方程數(shù)多于未知數(shù)個數(shù)時，對應的系數(shù)矩陣G的行數(shù)大于列數(shù)，此時方程組被稱為是超定方程組。設G=(g

12、iu)m×n，當m>n時即所謂的高矩陣，絕大多數(shù)情況下，超定方程組沒有古典意義下的解。超定方程組的最小二乘解是一種廣義解，是指使殘差r = b GX 的2-范數(shù)達取極小值的解，即超定方程組的最小二乘解是指正規(guī)方程組GTGX=GTb的解。如果系數(shù)矩陣(GTG)可逆，則正規(guī)方程組有唯一解。此時，最小二乘解可以形式地寫為如下形式X=(GTG)-1GTb兩種常用的方法如下1用對稱矩陣的三角分解法解正規(guī)方程組GTGX = GTb；記A=GTG，則A是對稱矩陣，由三角分解A = L D LT，其中L是下三角矩陣，D是對角矩陣。將這一算法寫過三個過程：解下三角方程組：LY1 = GTb；解對

13、角方程組：DY2 = Y1 ；解上三角方程組：LTY3 = Y22用矩陣的QR分解直接求解超定方程組由QR分解（正交三角分解）G=QR，其中Q是正交矩陣，R是上三角矩陣。將QR分解代入最小二乘解表達式中，得X=( R T QTQR)-1(QR)Tb = R-1QTb由此可知，用分解方法求超定方程組的最小二乘解只需求解上三角方程組RX=QTb模型1 熱源給恒溫箱體加熱，恒溫箱體給箱頂隔熱層加熱；恒溫箱及箱頂隔熱層向外界空氣散熱。此時，恒溫箱未達到68.2華氏度，恒溫箱體溫度高于箱頂隔熱層，如圖圖 3 模型1恒溫箱熱量傳播方向圖模型2恒溫箱達到68.2華氏度后，箱頂隔熱層未達到67.8華氏度，熱源

14、停止加熱；恒溫箱體給箱頂隔熱層加熱，恒溫箱及箱頂隔熱層向外界空氣散熱，如圖圖 4 模型2恒溫箱熱量傳播方向圖模型3恒溫箱及箱頂隔熱層溫度均高于67.8華氏度，熱源停止加熱；恒溫箱溫度低于箱頂隔熱層，箱頂隔熱層給恒溫箱加熱；恒溫箱及箱頂隔熱層向外界空氣散熱，如圖圖 5 模型3恒溫箱熱量傳播方向圖模型4恒溫箱體溫度低于67.8華氏度后，此時箱頂隔熱層溫度高于67.8華氏度，熱源加熱，熱源與箱頂隔熱層給恒溫箱體加熱，直到箱頂隔熱層與恒溫箱體溫度相等。圖 6 模型4恒溫箱熱量傳播方向圖模型建模根據(jù)已測的環(huán)境溫度、恒溫箱體和箱頂隔熱層溫度數(shù)據(jù)，在對每隔15min測量一次的環(huán)境溫度數(shù)據(jù)進行三次樣條差值后得

15、到每隔1min一個數(shù)據(jù)，如圖7所示。初始時刻，恒溫箱體與隔熱層溫度等于環(huán)境溫度。此時熱源給恒溫箱體加熱，恒溫箱體給箱頂隔熱層加熱；恒溫箱及箱頂隔熱層向外界空氣散熱。這段過程，恒溫箱未達到68.2華氏度，恒溫箱體溫度高于箱頂隔熱層。根據(jù)能量守恒，恒溫箱體增加的熱量與恒溫箱對環(huán)境空氣和箱頂隔熱層的散熱之和等于熱源提供的熱量。 （1）恒溫箱體對箱頂隔熱層的散熱等于箱頂隔熱層增加的熱量與箱頂隔熱層對環(huán)境空氣的散熱之和。 （2）式中各變量含義分別為是恒溫箱比熱容，是恒溫箱質量，是時間內恒溫箱溫度增量，是時間內恒溫箱溫度增量，是恒溫箱箱體熱導率，是恒溫箱某一時刻溫度，是環(huán)境溫度，是恒溫箱與外界環(huán)境的接觸表

16、面積，是恒溫箱體的厚度，是隔熱層某時刻溫度，是恒溫箱體和隔熱層的接觸表面積，是熱源時間內所釋放的熱量，是隔熱層的比熱容，是隔熱層的質量，是時間隔熱層溫度的增量，是隔熱層的熱導率，是隔熱層的厚度。將上式進一步簡化為 （3） （4）式中，；定為1min取差值后的80組數(shù)據(jù)中的前78組數(shù)據(jù)，求解超定方程中的A、B、C、D、E五個參數(shù)。A=22.4495B=63.9245C=0.0366D=122.5048E=0.4669將能量守恒方程變換為功率守恒方程： （5） （6）代入參數(shù)值可得： （7） （8）將第79和80組數(shù)據(jù)帶入上面的方程式，絕對誤差分別為0.5624，0.3112和0.602，0.33

17、1，在接受范圍內，驗證了方程式的正確性。圖 7 環(huán)境溫度實際測量與三次采樣根據(jù)式（7）和（8）反復迭代計算恒溫箱與箱頂隔熱層的溫度，計算前80min的數(shù)據(jù)，可以得到恒溫箱體與箱頂隔熱層溫度變化曲線，如圖8和圖9所示。圖 8 前80min恒溫箱溫度擬合曲線圖 9 前80min箱頂隔熱層溫度擬合曲線由圖可知，由建模得到的恒溫箱體與箱頂隔熱層溫度變化，與溫度實際測量曲線擬合得很好。再根據(jù)模型2、3、4討論的情況采用比較恒溫箱、箱頂隔熱層溫度與臨界溫度的高低，再選擇模型，具體程序見附錄。由此可以得到恒溫箱體和箱頂隔熱層溫度的變化曲線。仿真結果見圖10和圖11。圖 10 恒溫箱與箱頂隔熱層溫度變化曲線圖

18、 11 恒溫箱溫度變化曲線由建模結果可知，箱頂隔熱層始終低于68華氏度，箱頂隔熱層與環(huán)境溫度變化趨勢相近。恒溫箱體在剛開始持續(xù)受熱，約在260分鐘加熱到68華氏度，然后在67.8華氏度與68.2華氏度附近來回震蕩。對比圖10和圖11可知，隨著環(huán)境溫度的升高，恒溫箱體溫度震蕩頻率降低。環(huán)境溫度決定了恒溫箱溫度的震蕩頻率，即熱源開關的頻率。恒溫箱箱體溫度趨近于穩(wěn)定后（300分鐘1440分鐘），加熱時間Theat約為317分鐘，不加熱時間約為814分鐘,加熱時間與不加熱時間比值為317/814=0.4017。恒溫箱箱體溫度趨近于穩(wěn)定后（300分鐘1440分鐘），最高溫度68.4457華氏度，最低溫度

19、67.6212華氏度。模型的改進 恒溫箱溫度控制算法設計建模前提是熱源提供的速率一致，可以采用PID控制調節(jié)熱源的加熱速率。使恒溫箱加熱頻率降低，響應速率變快。PID 控制算法簡單、魯棒性能好、可靠性高，被廣泛應用于工業(yè)過程并取得了良好的控制效果。下圖是PID 控制的原理框圖。將給定值與實際輸出值相減，得到控制偏差：將偏差經過比例、積分、微分作用后線性相加，得到PID控制器輸出的控制量，輸入到被控對象進行控制： 其傳遞函數(shù)形式為：式中，是比例系數(shù)，是積分時間常數(shù)，是微分時間常數(shù)。比例環(huán)節(jié)：即時成比例地反應控制系統(tǒng)的偏差信號，偏差一旦產生，控制器立即動作以減小誤差。比例作用大，可以加快調節(jié)，減少

20、誤差。（比例也不能過大，過大的比例會使系統(tǒng)的穩(wěn)定性下降，甚至造成系統(tǒng)的不穩(wěn)定。）積分環(huán)節(jié)：主要作用是消除靜差，提高系統(tǒng)的控制精度。積分作用的強弱取決于積分時間常數(shù)，其值越大，積分作用越弱；反之越強，引入積分調節(jié)可使系統(tǒng)穩(wěn)定性下降，動態(tài)響應變慢。微分環(huán)節(jié)：能反應偏差信號的變化趨勢（變化速率），并能在偏差信號變得太大之前引入一個修正信號，從而減小系統(tǒng)的超調，加快系統(tǒng)的過渡過程，減小調節(jié)時間。但微分作用對噪聲干擾有放大作用。 參考文獻1 李嘉林. 基于TEC制冷加熱管制熱的恒溫箱設計D.復旦大學,2012. 2 陳國初,張琳,郝寧眉等. 恒溫恒液位系統(tǒng)的動態(tài)機理建模與仿真J.青島大學學報(工程技術版

21、)，2003,18(2):46-50.3 葉慶銀,劉斌,臧潤清等.小型恒溫箱恒溫性能的理論及實驗研究J.制冷與空調,2007,7(4):48-50,81.DOI:10.3969/j.issn.1009-8402.2007.04.012.4 傅秦生. 熱工基礎與應用M.北京：機械工業(yè)出版社，2007.5 李乃成，梅立泉. 數(shù)值分析M.北京：科學出版社，2011.附錄三次樣條插值代碼：hwx,grc=textread('hwxwd.txt','%f%f','headerlines',1)wmwd=textread('wmwd.txt'

22、,'%f')figureplot(t1,wmwd)T_interp=0:1:1440;wmwd1=interp1(t1,wmwd,T_interp,'linear');wmwd2=interp1(t1,wmwd,T_interp,'spline');plot(t1(1:7),wmwd(1:7),'o',T_interp(1:90),wmwd1(1:90),'r+',T_interp(1:90),wmwd2(1:90),'g*')溫度曲線擬合代碼：調用函數(shù)：Tair=wmwd2;wdqx(0,1000

23、,34.6256,34.6256,Tair,1,1000)龍格庫塔法數(shù)值積分函數(shù)：function hwxwd,grcwd = wdqx( x0,xn,Tbox0,Tgere0,Tair,h,n )%定義初值%h:步長；n:節(jié)點個數(shù)x=x0:h:xn;y(1,1)=Tbox0;%賦恒溫箱、隔熱層溫度初值y(2,1)=Tgere0;i=1;while i<=nwhile (i<=n)&&(y(1,i)<68.2) K1(:,1)=h*func1(y(1,i),y(2,i),Tair(i),h*func3(y(1,i),y(2,i),Tair(i)' K2

24、(:,1)=h*func1(y(1,i)+K1(1,1)/2,y(2,i)+K1(2,1)/2,Tair(i),h*func3(y(1,i)+K1(1,1)/2,y(2,i)+K1(2,1)/2,Tair(i)' K3(:,1)=h*func1(y(1,i)+K2(1,1)/2,y(2,i)+K2(2,1)/2,Tair(i),h*func3(y(1,i)+K2(1,1)/2,y(2,i)+K2(2,1)/2,Tair(i)' K4(:,1)=h*func1(y(1,i)+K3(1,1),y(2,i)+K3(2,1),Tair(i),h*func3(y(1,i)+K3(1,1),y(2,i)+K3(2,1),Tair(i)' y(:,i+1)=y(:,i)+1/6*(K1+2*K2+2*K3+K4);%加熱循環(huán) i=i+1;end while (i<=n)&&(y(1,i)>67.8) K1(:,1)=h*func2(y(1,i),y(2,i),Tair(i),h*fun