仿真報(bào)告修改后

1、基于matlab對(duì)管式反應(yīng)器內(nèi)傳熱過(guò)程的研究1.熱傳導(dǎo)的基本原理 熱傳導(dǎo)是由物質(zhì)內(nèi)部分子、原子和自由電子等微觀粒子的熱運(yùn)動(dòng)而產(chǎn)生的熱量傳遞現(xiàn)象。熱傳導(dǎo)的機(jī)理非常復(fù)雜，簡(jiǎn)而言之，非金屬固體內(nèi)部的熱傳導(dǎo)是通過(guò)相鄰分子在碰撞時(shí)傳遞振動(dòng)能實(shí)現(xiàn)的；金屬固體的導(dǎo)熱主要通過(guò)自由電子的遷移傳遞熱量；在流體特別是氣體中，熱傳導(dǎo)則是由于分子不規(guī)則的熱運(yùn)動(dòng)引起的。1.1溫度場(chǎng)和等溫面任一瞬間物體或系統(tǒng)內(nèi)各點(diǎn)溫度分布的空間，稱(chēng)為溫度場(chǎng)。在同一瞬間，具有相同溫度的各點(diǎn)組成的面稱(chēng)為等溫面。因?yàn)榭臻g內(nèi)任一點(diǎn)不可能同時(shí)具有一個(gè)以上的不同溫度，所以溫度不同的等溫面不能相交。1.2溫度梯度 從任一點(diǎn)開(kāi)始，沿等溫面移動(dòng)，如圖1所示

2、，因?yàn)樵诘葴孛嫔蠠o(wú)溫度變化，所以無(wú)熱量傳遞；而沿和等溫面相交的任何方向移動(dòng)，都有溫度變化，在與等溫面垂直的方向上溫度變化率最大。將相鄰兩等溫面之間的溫度差t與兩等溫面之間的垂直距離n之比的極限稱(chēng)為溫度梯度，其數(shù)學(xué)定義式為：1.3傅里葉定律 導(dǎo)熱的機(jī)理相當(dāng)復(fù)雜，但其宏觀規(guī)律可以用傅里葉定律來(lái)描述，其數(shù)學(xué)表達(dá)式為： 圖1 溫度梯度與傅里葉定律 傅里葉定律表明：在熱傳導(dǎo)時(shí)，其傳熱速率與溫度梯度及傳熱面積成正比。2.模型描述 如圖所示，一高溫管式反應(yīng)器由內(nèi)徑 5cm，壁厚 2.5cm 的陶瓷材料制成，為實(shí)現(xiàn)高溫條件下的反應(yīng)，需要對(duì)反應(yīng)器進(jìn)行保溫，保溫層厚度為 5cm，陶瓷和保溫層的物性為： 陶瓷：1=

3、2600 kg/m3 c1=1150 J/(kg·K ) k1=3.0 W/(m·K) 保溫層：2=600 kg/m3 c2=200 J/(kg·K ) k2=0.2 W/(m·K) 反應(yīng)器初始溫度是 300 K，內(nèi)表面突然被暴露在 1500K 的反應(yīng)物流中，管壁和反應(yīng)物流間的對(duì)流換熱系數(shù)為 500 W/(m2·K)。請(qǐng)給出保溫層內(nèi)表面溫度隨時(shí)間的變化規(guī)律，以及保溫層內(nèi)溫度分布，討論保溫層厚度的影響。（保溫層外表面視為絕熱） 圖2 管式反應(yīng)器的結(jié)構(gòu)3.構(gòu)建數(shù)學(xué)模型 該問(wèn)題是柱坐標(biāo)系下的一維非穩(wěn)態(tài)導(dǎo)熱問(wèn)題， 比較特殊之處在于該體系分陶瓷層、 保溫

4、層兩部分，共有內(nèi)、中、外三個(gè)邊界條件。 建立柱坐標(biāo)系并傳熱微分方程進(jìn)行化簡(jiǎn)，不難得到如下關(guān)系式： 其中，無(wú)量綱溫度（=1500K,=300K） 根據(jù)題目給出的參數(shù)可校驗(yàn)發(fā)現(xiàn)該圓柱是厚壁物體， 則必須解以上偏微分方程， 考慮用差分方程的方法給出數(shù)值解。 為保證結(jié)果的穩(wěn)定性不受步長(zhǎng)影響，選擇用時(shí)間向后、空間一階向后、空間二階中心差分方程。 以上上標(biāo)n表示時(shí)間,下標(biāo)i表示半徑。t,r分別為時(shí)間、半徑的步長(zhǎng)。 由此對(duì)于給定時(shí)間點(diǎn)， 可以得到一個(gè)大規(guī)模方程組， 自變量為各個(gè) r 位置的無(wú)量綱溫度。 以上即可解決三個(gè)非邊界點(diǎn)以外的所有差分點(diǎn)的列式，下面再給出邊界點(diǎn)的列式。 (1) r=0.0025m 處，

5、陶瓷材料與熱流界面處在該位置可通過(guò)熱平衡法得到其表達(dá)式: 進(jìn)而做差分近似處理得 對(duì)于無(wú)內(nèi)熱源的邊界微元層，我們可以認(rèn)為在足夠小的微元層內(nèi)無(wú)熱量積累。則有： 進(jìn)而做差分近似處理得 這樣，方程求解會(huì)變得簡(jiǎn)單一些。 (2) r=0.0050m 處，陶瓷材料與保溫層界面處，同理的可以由熱量衡算給出其微分表達(dá)式。進(jìn)而做差分近似處理得 此處左邊表達(dá)式用向前差分，右邊為向后差分。 (3) r=0.010m 處，保溫層外界面處，壁面處絕熱，故近似認(rèn)為，由此經(jīng)整理可得到完整的差分方程組。 不妨選取r=0.0005m，t=1s,并用x代替，則有： 方程(1)： 方程(2)方程(50): 方程(51): 方程(52

6、)方程(150): 方程(151): 1,2分別為陶瓷層，保溫層的導(dǎo)熱系數(shù)。 借助 Matlab 對(duì)以上大矩陣做高斯-賽達(dá)爾迭代進(jìn)行求解。4.結(jié)果分析4.1 整體圖像 取時(shí)間 t=1:2000s ，空間半徑 r=0.0025:0.0005:0.0100，將無(wú)量綱溫度轉(zhuǎn)換為對(duì)溫度 T。先從整體上觀察溫度場(chǎng)隨時(shí)間、空間的變化規(guī)律。 圖3 整體溫度場(chǎng)變化圖像從圖上可得到信息： 1）隨著時(shí)間增加，同一空間點(diǎn)（同一薄圓柱殼層）的溫度逐漸增加， 最初為 300K，在 t=2000s時(shí)各空間點(diǎn)均已經(jīng)上升至 1200K 以上。 2）同一時(shí)間，隨著半徑的增大（遠(yuǎn)離熱反應(yīng)物流） ，溫度逐漸降低。以 t=2000s

7、 為例，陶瓷層內(nèi)壁 1476.6K，保溫層外壁 1205.1K。 3）還可看出在兩層導(dǎo)熱介質(zhì)交界處存在溫度梯度突變；余位置溫度梯度連續(xù)，曲線、曲面光滑。 4.2 保溫層內(nèi)表面溫度 Tc 分析 先做出 Tc-t 圖像，再截取部分 Tc-t 觀察其變化規(guī)律。 圖4 Tc-t關(guān)系圖像 表1 Tc-t 變化數(shù)據(jù) t/sTc/Kt/sTc/K1300.000051308.91472300.005052309.57833300.009753310.26824300.014354310.98415300.019055311.72596300.023656312.49337300.028357313.2861

8、8300.033058314.10399300.037859314.946410300.042760315.813311300.047961316.704312300.053562317.618913300.059763318.556814300.066964319.517615300.075465320.500916300.085866321.506317300.098567321.506318300.114468322.533419300.134269324.651020300.158970325.740621300.189371326.850322300.226672327.979423

9、300.271873329.127824300.326174330.294825300.390875331.480126300.466976332.683227300.555877333.903828300.658678335.141429300.776579336.395530300.910780337.665831301.062381338.951932301.232482340.253433301.421983341.569834301.631884342.900835301.863085344.245936302.116386345.604837302.392487346.977138

10、302.692188348.362539303.015989349.760540303.364490351.170841303.738091352.593042304.137192354.026943304.562193355.472044305.013394356.928045305.490895358.394546305.994996359.871447306.525597361.358248307.082998362.854649307.666999364.360350308.2775100365.8750 截取了前100s的數(shù)據(jù)進(jìn)行分析，發(fā)現(xiàn)保溫層內(nèi)壁處溫度變化還是比較緩慢的。經(jīng)過(guò)10

11、0s的時(shí)間溫度僅從300K升至365K左右。 當(dāng)t=2000s時(shí) 該位置的溫度上升至1423.9K,已經(jīng)比較接近熱流溫度。 之所以說(shuō)其溫度變化緩慢，是因?yàn)樘沾蓪觾?nèi)溫度梯度很陡，以t=100s為例給出溫度場(chǎng)圖.看以下圖5可發(fā)現(xiàn),保溫層內(nèi)溫度變化相比陶瓷層十分緩慢，明顯這是因?yàn)楹笳哂懈〉膶?dǎo)熱系數(shù)k。 圖5 t=100s時(shí)整體溫度場(chǎng)4.3 保溫層內(nèi)溫度場(chǎng)分析 圖6 保溫層內(nèi)溫度場(chǎng) 保溫層內(nèi)溫度隨位置、時(shí)間的變化規(guī)律與之前整體圖像規(guī)律一致。但是與前面對(duì)比易知，在保溫層內(nèi)溫度無(wú)論是隨時(shí)間還是隨空間的變化都比陶瓷層內(nèi)平緩。4.4 保溫層厚度的影響 定性的分析我們很容易知道： 1.隨保溫層厚度增加，保溫層

12、外表面溫度會(huì)降低，單位表面積散熱量會(huì)減小 2.但是保溫層變厚后，保溫層外表面積增大，故總的散熱量是否減小則不一定 3.但基于該問(wèn)題， 題目認(rèn)為保溫層是絕熱的， 那么保溫層的厚度如何改變都不會(huì)影響散熱。 下面試取邊界層厚度 d=0.025 m, d=0.060 m（原先 d=0.050 m） 觀察整體溫度場(chǎng)圖像 以及保溫層內(nèi)表面的溫度變化曲線、 t=100s 時(shí)整體溫度場(chǎng)。 將結(jié)果與之前比較。 1）d=0.025 m < d0 圖7 d=0.025m 整體溫度場(chǎng) 顯然可以看出此時(shí)保溫層外壁的溫度變高了，達(dá)到 1400K 以上。圖8 d=0.025m 保溫層內(nèi)表面變溫曲線 與原圖像對(duì)比，總體

13、變化趨勢(shì)不變，但溫度上升變快了，當(dāng)t=2000s時(shí)保溫層內(nèi)表面溫度達(dá)1455.9K,比d=0.0050m 時(shí)高了200K左右。 圖9 d=0.025m時(shí)t=100s 整體圖像 保溫層內(nèi)側(cè)溫度提高，其余基本一致。 2）d=0.060 m > d0圖10 d=0.060m 整體溫度場(chǎng)圖像 對(duì)比之前兩張圖，當(dāng)保溫層變厚，可以看出保溫效果變好，例如 t=2000s 時(shí)保溫層外壁溫度低至1100K以下。 圖11 d=0.060m t=100s 整體圖像 與前兩個(gè)d的情形相比，變化趨勢(shì)不變；與d=0.05m相比，陶瓷層內(nèi)溫度降低不明顯，例如保溫層內(nèi)壁溫度此時(shí)為365.875K，與d=0.05m的數(shù)據(jù)

14、持平。但是因?yàn)槠浔貙幼兒?，其外壁溫度更低，在t=2000s時(shí)外壁僅1035K。 圖12 d=0.060m 保溫層內(nèi)壁溫度變化 對(duì)比發(fā)現(xiàn) 這條曲線與d=0.05m時(shí)的曲線基本一致，t=2000s時(shí)保溫層內(nèi)壁溫度也與之基本持平。 可見(jiàn)此時(shí)增加保溫層厚度對(duì)陶瓷層的溫度場(chǎng)影響不大，但是使保溫層的溫度場(chǎng)“延長(zhǎng)”了，以致外壁溫度更低。故可得出結(jié)論： 1. 原題情形下，t=2000s 時(shí) 保溫層內(nèi)壁 1476.6K，外壁 1205.1K 2. 隨著保溫層厚度增加，保溫層外壁溫度降低，但保溫層內(nèi)壁溫度未必受明顯影響 5.matlab程序dr=0.0005; dt=1; M=2600*1150*( dr)2/

15、 dt; N=600*200*( dr)2/ dt; for i=1:151 r(i)=0.025+(i-1)* dr; end T=ones(151,2500); x=zeros(151,150,2500); A1=eye(151,151); A2=zeros(151,151); alpha1=3/(2600*1150); alpha2=0.2/(600*200); for j=2:50 A1(j,j)=(dr2)/(alpha1* dt)- dr/r(j)+2; end for j=1:150 A2(j,j+1)=-1; end A1(1,1)=500+3/dr; A2(1,2)=-3/d

16、r; A3=zeros(151,151); for j=1:150 A3(j+1,j)= dr /r(j+1)-1; end for j=52:150 A1(j,j)=( dr 2)/(alpha2* dt)- dr /r(j)+2; endA1(51,51)=-(3+0.2); A2(51,52)=0.2; A3(51,50)=3; A=A1+A2+A3;A(151,151)=1; A(151,150)=-1;b=zeros(151,1,2500);L=-1*tril(A,-1); U=-1*triu(A,1); D=diag(diag(A); j=1; x(:,:,1)=1; for n=

17、2:2000 b(1,1,n)=0; %由前一時(shí)間循環(huán)結(jié)果給下一時(shí)間的方程組b賦值 for q=2:50 b(q,1,n)= dr 2/(alpha1* dt) *x(q,j,n-1); end b(51,1,n)=0; for q=52:150 b(q,1,n)= dr2/(alpha2* dt) *x(q,j,n-1); end b(151,1,n)=0;for j=1:149 %設(shè)計(jì)循環(huán)，高斯-賽德?tīng)柕蠼?符合要求時(shí)跳出循環(huán),得到迭代次數(shù) j if j=1 x(:,j+1,n)=inv(D-L)*U*x(:,j,n)+inv(D-L)*b(:,:,n); elseif norm(x(:,j,n)-x(:,j-1,n),inf) >10-5 x(:,j+1,n)=inv(D-L)*