永凍土層熱傳導(dǎo)問題數(shù)學(xué)建模

1、摘要本文針對永凍土層上關(guān)于路基熱傳導(dǎo)的問題，通過對不同材料層的密度、比熱容、傳熱系數(shù)進(jìn)行研究，建立微分方程模型，利用Matlab與Lingo軟件進(jìn)行求解。問題一，考慮在分析各材料層進(jìn)行后，給出空氣溫度傳入路基規(guī)律，以及各材料層的溫度分布。同時(shí)，已知外界的溫度是關(guān)于時(shí)間的函數(shù)，凍土層的溫度是不變的零下溫度。首先，我們通過中國選礦技術(shù)網(wǎng)以及中國天氣網(wǎng)分別獲得各材料層的密度、比熱容、傳熱系數(shù)等數(shù)據(jù)，和拉薩最近24小時(shí)的溫度數(shù)據(jù)。通過擬合得到溫度與時(shí)間的關(guān)系函數(shù)，建立一維熱傳導(dǎo)方程的微分方程模型。隨后利用向前差分的方法求出方程的近似數(shù)值解，因?yàn)榻缑嫣幍臒醾鲗?dǎo)率處于平衡，且溫度相等，則可以一層一層向下計(jì)

2、算得出各材料層的溫度分布規(guī)律。問題二，考慮在一些設(shè)備的支架不能固定在解凍土層上，必須固定在永凍土層中的情況下，地下土層的解凍位置，并給出解凍砂土與凍結(jié)砂土的分界線。由問題一的求解可以計(jì)算出L4的值，在已知上界x值與下界的0c溫度值后，同問題一的求解方法可以給出解凍砂土與凍結(jié)砂土的分界線。問題三，考慮結(jié)合溫度分布、成本及耐用性，給出各層材料的最佳厚度。結(jié)合鐵路建設(shè)施工保障，我們將耐用性作為出發(fā)點(diǎn)，分別從壓強(qiáng)及壓實(shí)度考慮耐用性的約束條件。由于壓實(shí)度與含水量存在聯(lián)系，而含水量與溫度存在關(guān)系，故建立起壓實(shí)度與溫度的相關(guān)關(guān)系。將成本作為目標(biāo)函數(shù)，壓實(shí)度與壓強(qiáng)的限制作為約束條件，建立線性規(guī)劃模型，由此解出

3、最少成本為46334.72元。問題四，考慮在以上問題的基礎(chǔ)上，結(jié)合我國青藏鐵路永凍土層地基進(jìn)行仿真，并為施工單位提出合理建議。因?yàn)楸绢}的前三問即是在查閱青藏鐵路路基修建相關(guān)數(shù)據(jù)的基礎(chǔ)上進(jìn)行的，故問題四的仿真即已經(jīng)得到相應(yīng)的解決。通過對以上問題的求解進(jìn)行合理性分析，即可對施工單位給出合理建議。為了簡化計(jì)算量，提高求解速度，本題中的微分方程模型使用向前差分的方法求出近似數(shù)值解，而且對模型的可行性及有效性進(jìn)行了一定的分析，所得結(jié)果十分合理。本文的優(yōu)點(diǎn)在于利用差分方法求解一維熱傳導(dǎo)微分方程模型的近似數(shù)值解，使得材料內(nèi)界面的條件處理得較為容易。同時(shí)，在前三個(gè)問題的求解中，將問題背景設(shè)定在青藏鐵路的修建中

4、，在一定程度上對問題四的求解提供了較大的幫助。關(guān)鍵詞：熱傳導(dǎo)問題拋物型方程數(shù)值模擬1問題重述在永凍土地上鋪設(shè)道路、飛機(jī)跑道和某些結(jié)構(gòu)的地基。分析這類地基的結(jié)構(gòu),有瀝青層、混凝土層、干砂層、石子層、絕緣材料層，再下面是濕的沙土層和凍土層(參見圖1)。已知，外界的溫度是關(guān)于時(shí)間的函數(shù)，凍土層的溫度是不變的零下溫度。請回答以下問題：圖1:永凍土地地基結(jié)構(gòu)圖(1)分析空氣溫度傳入路基規(guī)律，各材料層的溫度分布；(2)由于一些設(shè)備的支架不能固定在解凍土層上，必須固定在永凍土層中,因此確定地下土層的解凍位置非常重要，請給出解凍砂土與凍結(jié)砂土的分界線；(3)請給出各層材料的最佳厚度(結(jié)合溫度分布、成本和耐用性

5、)；(4)請結(jié)合我國青藏鐵路永凍土層地基進(jìn)行仿真，并給施工單位給出合理建議。2問題分析本問題要求建立模型和設(shè)計(jì)算法，在永凍土地上鋪設(shè)道路、飛機(jī)跑道和某些結(jié)構(gòu)的地基的問題背景中，結(jié)合溫度分布、成本和耐用性，給出各層材料的最佳厚度，并對我國青藏鐵路永凍土層地基進(jìn)行仿真。由于問題的局限性，我們首先應(yīng)該找到合理的數(shù)據(jù)來對相關(guān)背景進(jìn)行一定的了解，定位模型的求解方法，循序漸進(jìn)地對此問題進(jìn)行合理的求解。問題一中，要求在題中所給的情況下，分析空氣溫度傳入路基規(guī)律，以及各材料層的溫度分布。經(jīng)過分析，我們發(fā)現(xiàn)建立模型所需的數(shù)據(jù)不足，故應(yīng)首先搜集題目中各材料熱傳導(dǎo)相關(guān)的密度、導(dǎo)熱系數(shù)、比熱容以及環(huán)境溫度等數(shù)據(jù)。問題

6、二中，要求確定地下土層的解凍位置，給出解凍砂土與凍結(jié)砂土的分界線，以應(yīng)對由于一些設(shè)備的支架不能固定在解凍土層上，必須固定在永凍土層中的情況。由問題一建立的模型，本問題即是求出溫度為0七時(shí)所對應(yīng)的距離。問題三中，要求結(jié)合溫度分布、成本和耐用性，給出各層材料的最佳厚度。在分析問題一、二后，則是對各目標(biāo)進(jìn)行權(quán)重約束，以求得最佳的厚度。問題四中，要求結(jié)合我國青藏鐵路永凍土層地基進(jìn)行仿真，并給施工單位給出合理建議。3模型假設(shè)1、假設(shè)環(huán)境溫度不會(huì)發(fā)生突變，沒有極端天氣出現(xiàn)；2、假設(shè)路基材料分布均勻；3、假設(shè)考慮壓實(shí)度時(shí)，可以將半凍砂土層以上的四層材料合并為一層;4、假設(shè)對于各層只考慮重力所產(chǎn)生的壓強(qiáng)。4符

7、號說明a(t)外界的溫度關(guān)于時(shí)間的函數(shù)To凍土層的溫度K(i=1,2,5)每個(gè)材料層的熱傳導(dǎo)率q每個(gè)材料層的比熱R每個(gè)材料層的密度u要求解的第i個(gè)材料層的溫度分布x地下某一點(diǎn)距地面的跑離t時(shí)刻值li每層材料的厚度其他局部符號在引用時(shí)給出了具體說明。5模型的建立與求解5.1 模型一的建立與求解問題一中要求分析空氣溫度傳入路基規(guī)律，以及各材料層的溫度分布。我們知道，一切穩(wěn)定的數(shù)學(xué)物理問題都可以用橢圓型微分方程模型來描述，若在穩(wěn)態(tài)之前有一個(gè)過程，則討論這種漸進(jìn)穩(wěn)定的數(shù)學(xué)物理過程可以用拋物型為方程模型來描述1。5.1.1 模型一的建立分析題目可知，地基的結(jié)構(gòu)有瀝青層，混凝土層，干砂石層，防水隔溫層，半

8、凍砂土層，外界的溫度是關(guān)于時(shí)間的函數(shù)a(t),凍土層的溫度是不變的零下溫度To,每個(gè)材料層都有熱傳導(dǎo)率K(i=1,2,，5),問題一則要分析空氣溫度傳入路基的規(guī)律，以及各材料層的溫度分布2。首先，由于路基各材料層是均勻的，所以要分析的熱傳導(dǎo)方程可歸為一維的熱傳導(dǎo)方程研究。路基各材料層則化為一維熱傳導(dǎo)方程為2生(1).:tjG",；x其中，Ui是要求解的第i個(gè)材料層的溫度分布；ki,q,R為此材料層的傳熱系數(shù)，比熱和密度，它們都是已知的常數(shù)，且工0。CiPi由模型假設(shè)可知，層與層的界面處傳導(dǎo)率處于平衡，且溫度相等，即ki=ki+Ui1,Ui=Ui書，x=L,(i=1,2,3,4)(2)

9、x二x基于以上分析，在半凍砂土層假設(shè)為濕砂層(4是其傳熱系數(shù))的情況下，所建立的路基熱傳導(dǎo)微分方程模型為，1-2cuikiuui笈CiPiex2旦=ki書UUi=Ui書，x=L,(i=1,2,3,4)(3)dx改Uo(O,t)=a(t)iU5(L5,t)=To5.1.2 模型一的求解由以上所建立的拋物型微分方程模型進(jìn)行分析，對問題一的求解我們給出了如下的具體步驟：(1)依據(jù)題意，我們在中國選礦技術(shù)網(wǎng)中找到瀝青層，混凝土層，干砂石層，防水隔溫層和半凍砂土層的密度，比熱容和熱導(dǎo)系數(shù)等值，部分難以確定的材料選取了最接近的材料進(jìn)行替代，所找到的數(shù)據(jù)材料如下表。表1各層材料相關(guān)參數(shù)表各材料層干密度耳(K

10、g/m3)熱導(dǎo)系數(shù)九2(w/m*)比熱容C(KJ/(Kg朱)價(jià)格P(元/t)瀝青層(石油瀝青)14000.271.682800混凝土層(c30混凝土)23001.510.92770.5干砂石16000.581.0170防水隔溫層（乳化瀝青膨脹珍珠巖）4000.121.551000砂土層（濕砂土）23501.421.07500以上則為各材料層的密度、熱導(dǎo)系數(shù)、比熱容和價(jià)格（價(jià)格因素在問題三中需要考慮，故在這里一起給出）的數(shù)據(jù)。止匕外，由于溫度的傳導(dǎo)涉及到凍土地帶的環(huán)境溫度，故我們結(jié)合青藏鐵路的修建，選取拉薩最近24小時(shí)內(nèi)的溫度3考慮到中國天氣網(wǎng)只可以給出最近24小時(shí)的較精確溫度，故我們問題的季節(jié)

11、背景選擇在夏秋之交4。表2拉薩24小時(shí)各時(shí)刻溫度表時(shí)亥IJ000102030405060708091011溫度121212131313131313121313時(shí)亥IJ121314151617181920212223溫度151718192122222120161312以上則為拉薩最近24小時(shí)各時(shí)刻的溫度表。（2）通過查閱鐵路建設(shè)國家規(guī)范，我們將瀝青層的厚度設(shè)定為4cm,混凝土層設(shè)定為10cm,干砂石層設(shè)定為7cm,防水隔溫層設(shè)定為10cm5。（3）通過擬合數(shù)據(jù)得到外界溫度關(guān)于時(shí)間的關(guān)系式，根據(jù)問題一中建立的第一種邊界的一維拋物型方程，考慮到解析解不易求得，故在這里我們利用差分的方法求出其近似的數(shù)

12、值解6。差分求解拋物型偏微分方程的過程則是首先對xt平面進(jìn)行網(wǎng)格剖分，分別取h,七為x方向和t方向的步長，用兩族平行直線x=Xk=kh（k=±1,±2：）,節(jié)點(diǎn)為（Xk,tj）,為簡單起見,t=tj=jT（j=0,1,2,），將Xt平面剖分成矩形網(wǎng)格,記（k,j）=（Xk，tj）,u（k,j）=u（Xk,tj）,則在網(wǎng)格內(nèi)點(diǎn)（k,j）處，對一采取向前差商t.2u(k,j1)-u(k,j)公式，對T采用二階中心差商公式，可得到一維熱傳導(dǎo)方程的差分近似：=0(4)kiu(k1)-2u(k)u(k-1)h（4）求出第一層底部的溫度后，令k=1,依據(jù)公式（2）即可求得下一層底部的溫

13、度，以下每一層材料溫度的求解依次循環(huán)以上步驟。5.1.3 模型一的結(jié)果及分析針又t5.1.2中的求解步驟，運(yùn)用Matlab編程求解。針對附錄一（1）中的程序文件，運(yùn)行可得到如下結(jié)果。外界環(huán)境溫度隨時(shí)間的變化函數(shù)為：a(t)=15.11+2.216cos(1.757t)+3.687sin(1.757t)-1.895cos(3.514t)+1.728sin(3.514t)外界溫度隨時(shí)間變化的函數(shù)如下圖:圖1外界溫度隨時(shí)間變化的函數(shù)第一層求解所得的溫度隨距離或時(shí)間變化的圖像分別為:圖2第一層溫度隨距離分布圖時(shí)間變量，圖3第一層溫度隨時(shí)間分布圖其他層溫度隨時(shí)間，隨距離的變化圖像見附錄一（2）0由上可知

14、，當(dāng)固定另一因素時(shí)，溫度隨著與路基表面距離的增加越來越小，其隨時(shí)間的變化則呈現(xiàn)先增加后減少的趨勢。5.2 模型二的建立與求解根據(jù)題意，問題二中要求確定地下土層的解凍位置，給出解凍砂土與凍結(jié)砂土的分界線，以應(yīng)對由于一些設(shè)備的支架不能固定在解凍土層上，必須固定在永凍土層中的情況。由問題一建立的模型，本問題即是求出09時(shí)所對應(yīng)的距離。5.2.1 問題二模型的建立根據(jù)問題一中所建立的模型，即可得出第五層溫度隨著距離變化的圖像關(guān)系。由圖分析，即可給出0P時(shí)所對應(yīng)的距離7。5.2.2 問題二模型的求解問題二中建立的模型用來確定第五層的解凍位置，即在已知上界的x值與下界溫度后，給出解凍砂土與凍結(jié)砂土的分界線

15、。由函數(shù)圖象分析，即可給出0P時(shí)所對應(yīng)的距離。5.2.3 問題二模型的結(jié)果分析針又t5.2.2中的求解步驟，依據(jù)問題一中所得的第四層溫度隨距離的圖像分析進(jìn)行求解。所得圖像如下圖：盤慶分存國嚏E寸間變量t空間變發(fā)H圖4第五層溫度隨距離變化圖像解得解凍砂土與凍結(jié)砂土之間的分界線為路面下120cm-140cm。結(jié)合國家公路施工技術(shù)規(guī)范，分析數(shù)據(jù)可知所得結(jié)果比較合理，故可作為問題而的求解答案。5.2.4 模型二的誤差分析因?yàn)椴罘址ㄊ菍?shù)值解得近似，所以對最后結(jié)果的誤差分析十分重要。依據(jù)附錄二中的程序文件，得出如下的誤差分析圖像：副片藜馬拈混缸七廠括:,然0肥C若J.u025彳n?015D.1DOS0J

16、110.20.3-0.40.5D.&07D.BD.91x圖5數(shù)值解與精確解比較圖像由上可知，數(shù)值解與精確解的最大相對誤差非常小，為0.42877%,且數(shù)值解的求解結(jié)果是穩(wěn)定收斂的，故該結(jié)果可以作為本問題的求解。5.3 模型三的建立與求解由題目可知，問題三需要在問題一、二的求解啟發(fā)中，結(jié)合溫度分布、成本和耐用性，給出各層材料的最佳厚度。在分析問題一、二后，則是對各目標(biāo)進(jìn)行權(quán)重約束，以求得最佳的厚度。5.3.1 問題三模型的建立由于題目中所給相關(guān)數(shù)據(jù)較少，由此我們思考對最佳厚度設(shè)置約束條件。對于耐用性，我們從壓實(shí)度和壓強(qiáng)兩方面進(jìn)行考慮。首先將上面四層材料合為一層以簡化計(jì)算，查閱文獻(xiàn)可知：壓

17、實(shí)度=|乂100%，試驗(yàn)干密度、試驗(yàn)干密度=實(shí)際密度x（1+含水率）而實(shí)際干密度則與溫度、水分等條件相關(guān)，則在此時(shí)將溫度分布考慮進(jìn)模型的建立。由此，需要擬合溫度與水分的關(guān)系8。而對于壓強(qiáng)，我們假設(shè)這里路基的承重來源只有上部材料的重力，則可得出P=Fjgg=VggSSS"g(6)如下推導(dǎo)公式：對上述公式中出現(xiàn)的局部符號說明如下:P表示某層材料所承受的壓強(qiáng)；F表示某層材料所承受的壓力；S表示承壓接觸面的表面積；V表示表面積S對應(yīng)的材料的體積；m表示體積V對應(yīng)的材料的質(zhì)量；:表示材料的密度；l表示材料層的厚度。將目標(biāo)函數(shù)設(shè)定為壓實(shí)度大于95%且能承受最大壓強(qiáng)為2760MPa時(shí)求得的最小成本

18、。用Q表示成本,minQ=1400x2.8li+2300x16712+1600，0.07la+400xII4(7)同時(shí)添加對各層壓強(qiáng)的限制，使各層所承受壓強(qiáng)大于路基建設(shè)所應(yīng)承載的最小壓強(qiáng)值。同時(shí)，分析數(shù)據(jù)知，壓實(shí)度=11(含水率)100%95%(8)基于以上分析，用r表示含水率，以(7)為目標(biāo)函數(shù)，(8)為約束條件建立優(yōu)化模型如下：minQ=14002.8l1+2300167l2+16000.07l3+4001l41400M10l1之700(9)2300M10l2之2760s.t1600m1013至960400M10l4>4001>0.95J+r5.3.2 問題三模型的求解由于問題

19、三是在結(jié)合溫度分布、成本、耐用性多種因素的基礎(chǔ)上對厚度問題進(jìn)行的研究，所以需要考慮各個(gè)因素的相關(guān)數(shù)據(jù)，以及各個(gè)因素之間的相互關(guān)系，具體的求解步驟如下：(1)查閱資料對溫度與水分含量的關(guān)系給出合理的擬合關(guān)系，并得出具體函數(shù)公式；(2)對壓強(qiáng)和壓實(shí)度的影響因素進(jìn)行化簡，得到相應(yīng)的約束關(guān)系；(3)列出壓實(shí)度大于95%時(shí)的目標(biāo)函數(shù)，在Lingo中求解線性規(guī)劃問題即可。5.3.3 問題三模型的結(jié)果分析針又t5.3.1中的模型，按照5.3.2中的求解步驟，首先利用Matlab求解得出溫度隨水分變化的函數(shù)關(guān)系式為：y=-4.139+33.95cos(0.01754x)+3.358sin(0.01754x)(

20、10)所畫出的函數(shù)圖像為：圖6溫度與水分的擬合圖像擬合度數(shù)據(jù)如下：gof=sse:0.6035rsquare:0.9991dfe:5adjrsquare:0.9985rmse:0.3474有數(shù)據(jù)可以看出，溫度與水分的函數(shù)關(guān)系擬合得較為合理。對所列出的目標(biāo)函數(shù)與約束條件，運(yùn)用Lingo編程求解，解得壓實(shí)度大于95%時(shí)的目標(biāo)函數(shù)所取得最小值為46334.72元。此時(shí)，第一、二、三、四層的厚度分別為5cm,12cm,6cm,10cm。同上分析，可以看出此數(shù)據(jù)與實(shí)際施工情況中的數(shù)據(jù)相差較小，故可作為本題的求解。5.4 模型四的建立與求解由題目可知，問題四要求結(jié)合我國青藏鐵路永凍土層地基進(jìn)行仿真，并給施

21、工單位給出合理建議。在對之前問題求解的基礎(chǔ)上，得出合理的分析與建議。5.4.1 問題四模型的建立問題一、二中選取的數(shù)據(jù)即是青藏鐵路路基建設(shè)的相關(guān)數(shù)據(jù)，同時(shí)將外界環(huán)境溫度設(shè)定為拉薩市最近24小時(shí)的溫度。故由此所得的厚度范圍即為對青藏鐵路永凍土層地基建設(shè)的仿真9。5.4.2 問題四模型的求解由上所述原理對問題四進(jìn)行求解，具體求解步驟如下：(1)由以上三個(gè)問題的求解，給出各材料層最佳厚度取值表；(2)通過對所得結(jié)果進(jìn)行分析，給出施工工人的合理化建議。5.4.3 問題四模型的結(jié)果分析由以上問題的求解，可以得出如下結(jié)果：表3各材料層最佳厚度取值表材料層厚度瀝青層(石油瀝青)5cm混凝土層(C30混凝土)

22、12cm干砂石6cm防水隔溫層(乳化瀝青膨脹珍珠巖)10cm通過對數(shù)據(jù)與國家青藏鐵路路基建設(shè)進(jìn)行比對，以上所得求解數(shù)據(jù)基本符合鐵路建設(shè)規(guī)范，故可以作為本問題的仿真結(jié)果。5.4.3問題四模型求解對施工工人的建議(1)由以上求解可以給出各層材料大致的厚度取值，可將此作為一個(gè)基本的施工參照給施工工人一定的理論指導(dǎo)；(2)由于模型建設(shè)中提出了相應(yīng)的假設(shè)，故在實(shí)際問題的施工時(shí)，相關(guān)施工單位可結(jié)合自身施工經(jīng)驗(yàn)對厚度取值進(jìn)行一定的調(diào)整，同時(shí)也應(yīng)進(jìn)行多次項(xiàng)目預(yù)試驗(yàn)，以充分考慮鐵路承重、修建成本以及后發(fā)事故。6模型的評價(jià)與改進(jìn)方向6.1 模型的評價(jià)由于本問題中，模型建立的背景為青藏鐵路的修建，同時(shí)對于各材料的相

23、關(guān)熱學(xué)參數(shù)進(jìn)行了唯一取值，對于此情況下的鐵路路基熱傳導(dǎo)問題可以合理有效地解決。如果修建環(huán)境發(fā)生改變便，則各材料層成分以及環(huán)境溫度均會(huì)發(fā)生改變，故應(yīng)對該模型進(jìn)行改進(jìn)繼而求得新的求解情況。6.1.1 模型的優(yōu)點(diǎn)(1)本模型中，所選取的數(shù)據(jù)均為與各材料層成分接近程度最大的材料的熱學(xué)參數(shù)，建立的線性模型函數(shù)設(shè)置合理，對于題中所給問題可以得到滿意的求解方案。(2)算法設(shè)計(jì)較為合適，計(jì)算量較小，能夠高效運(yùn)行出結(jié)果。(3)利用差分代替微分的思維求解一維熱傳導(dǎo)微分方程模型，可以使材料內(nèi)界面的條件處理得較為容易。(4)本模型在前三個(gè)問題的求解中，將問題背景設(shè)定在青藏鐵路的修建中，在一定程度上對問題四的求解提供了

24、較大的幫助。6.1.2 模型的缺點(diǎn)(1)問題一中，我們所選取的材料數(shù)據(jù)并不能充分考慮各層材料成分，溫度數(shù)據(jù)也僅為一個(gè)地點(diǎn)24小時(shí)內(nèi)的數(shù)據(jù)；(2)問題三中，在考慮壓實(shí)度時(shí)，將半凍砂土層以上材料壓縮為一種材料，在一定程度上簡化了問題模型。6.2 模型的改進(jìn)方向(1)對于所選取的溫度及材料熱學(xué)數(shù)據(jù)，可以選取多組求平均值，以使求解結(jié)果更具代表性；(2)求解一維熱傳導(dǎo)微分方程時(shí)，對于初始條件的選取稍有欠缺，可以結(jié)合實(shí)際情況選擇更好的初值條件。(3)對于差分求解一維熱傳導(dǎo)微分方程，雖在一定程度上減少了求解困難，但其步長選取不太合理，可以進(jìn)一步研究問題模型，以選取更合理地步長，以使求解結(jié)果與真實(shí)值更加接近。

25、參考文獻(xiàn)1張辰熙，于明.季節(jié)凍土層下人工凍土墻溫度場實(shí)驗(yàn)研究J.低溫建筑術(shù).2010(6):144-146.2文斌，吳青柏，蔣觀利，張鵬.模擬退火優(yōu)化算法的凍土熱傳導(dǎo)參數(shù)反分析J.巖土力學(xué).2013(8):15-16.3王建偉，李書民，李振卿.凍土層等效導(dǎo)熱系數(shù)的探討J.防滲技術(shù).2001(3):36-38.4李述訓(xùn)，吳通化.凍土溫度狀況研究方法和應(yīng)用分析J.冰川凍土.2004(4):142-148.5候彬彬，盧德唐，陳永輝，孔祥言.非達(dá)西滲流模型模擬青藏鐵路拋石路基凍土層溫度J.水工滲流研究與應(yīng)用進(jìn)展第五屆全國水利工程滲透學(xué)術(shù)研討會(huì)論文集.2006.6史冊.熱傳導(dǎo)方程有限差分法的MATLAB

26、實(shí)現(xiàn)J.咸陽師范學(xué)院學(xué)報(bào).2009(4):24-26.7張世民.青藏鐵路多年凍土路基熱-力穩(wěn)定性數(shù)值仿真分析J.土木建筑與環(huán)境工程.2012(34):3-7.8穆彥虎，馬巍,牛富俊，李國玉，王大雁，劉永智.青藏鐵路多年凍土區(qū)普通路基熱狀況監(jiān)測分析J.冰川凍土.2014(2):10-13.9田立慧，凌賢長，王力娜,張鋒.青藏鐵路高溫多年凍土區(qū)列車行駛路基長期永久變形數(shù)值模擬研究J.地震工程學(xué)報(bào).2014(4):21-25.附錄問題一外界環(huán)境溫度隨時(shí)間變化的函數(shù)關(guān)系式擬合程序creatFit.mfunctionfitresult,gof=createFit(x,y)%CREATEFIT(X,Y)%

27、Createafit.%Datafor'temperature'fit:%XInput:x%YOutput:y%Output:%fitresult:afitobjectrepresentingthefit.%gof:structurewithgoodness-offitinfo.%另請參閱FIT,CFIT,SFIT.%由MATLAB于17-Aug-201519:27:00自動(dòng)生成%Fit:'temperature'.xData,yData=prepareCurveData(x,y);%Setupfittypeandoptions.ft=fittype('

28、fourier2');opts=fitoptions('Method','NonlinearLeastSquares');opts.Display='Off;opts.Normalize='on'opts.StartPoint=000001.9316882339819;%Fitmodeltodata.fitresult,gof=fit(xData,yData,ft,opts);%Plotfitwithdata.figure('Name','temperature');h=plot(fitresult

29、,xData,yData);legend(h,'yvs.x','temperature','Location','NorthEast');%Labelaxesxlabel('x');ylabel('y');gridonnihe.mx=1:24;y=121213131313131312131315171819212222212016131212;fitresult,gof=createFit(x,y)偏微分方程求解程序functionUxt=PDEParabolicClassicalExplicit(

30、uX,uT,phi,psi1,psi2,M,N,C)%古典顯式格式求解拋物型偏微分方程%Uxt=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C)%方程：u_t=C*u_xx0<=x<=uX,0<=t<=uT%初值條件：u(x,0)=phi(x)%邊值條件：u(0,t)=psi1(t),u(uX,t)=psi2(t)%輸出參數(shù)：U-解矩陣，第一行表示初值，第一列和最后一列表示邊值，第二行表示第2層%x-空間變量%t-時(shí)間變量%輸入?yún)?shù)：uX-空間變量x的取值上限%uT-時(shí)間變量t的取值上限%phi-初值條件，

31、定義為內(nèi)聯(lián)函數(shù)%psi1-邊值條件，定義為內(nèi)聯(lián)函數(shù)%psi2-邊值條件，定義為內(nèi)聯(lián)函數(shù)%M-沿x軸的等分區(qū)間數(shù)%N-沿t軸的等分區(qū)間數(shù)%C-系數(shù)，默認(rèn)情況下C=1%應(yīng)用舉例：%uX=1;uT=0.2;%M=15;N=100;C=1;%phi=inline('sin(pi*x),);psi1=inline(,0');psi2=inline(,0,);%Uxt=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C);%inline('263+15.11+2.216*cos(1.757*t)+3.687*sin(1.7

32、57*t)-1.895*cos(2*1.757*t)+1.728*sin(2*1.757*t)')%設(shè)置參數(shù)C的默認(rèn)值ifnargin=7C=1;end%計(jì)算步長dx=uX/M;%x的步長dt=uT/N;%t的步長x=(0:M)*dx;t=(0:N)*dt;r=C*dt/dx/dx;%步長比r1=1-2*r;ifr>0.5disp('r>0.5,不穩(wěn)定')end%計(jì)算初值和邊值U=zeros(M+1,N+1);fori=1:M+1U(i,1)=phi(x(i);endforj=1:N+1U(1,j)=psi1(t(j);U(M+1,j)=psi2(t(j);

33、end%逐層求解forj=1:Nfori=2:MU(i,j+1)=r*U(i-1,j)+r1*U(i,j)+r*U(i+1,j);endendU=U'%作出圖形mesh(x,t,U);title('溫度分布圖像')ylabel('空間變量x')xlabel('時(shí)間變量t')zlabel('溫度U')return;溫度隨時(shí)間，隨距離變化的圖像對問變量1圖7第二層溫度隨時(shí)間的變化圖像3出溫度紗布剪套圖8第二層溫度隨距離的變化圖像10用度分布圖情111171RdII-I1iin09溫度分布圖溫on?口出oce009r012o.u

34、口徑o1002時(shí)間變量t圖9第二層溫度隨時(shí)間的變化圖像堂間變量圖10第三層溫度隨距離的變化圖像瀛宸分布園藝圖11第四層溫度隨時(shí)間的變化圖像空間變重K圖12第四層溫度隨距離的變化圖像問題二模型二誤差分析程序：%求土壤溫度分布clear;clc;formatshortea=input('請輸入系數(shù)a的值：')；l=input('請輸入長度l的值：)；M=input('請輸入將區(qū)間0,l等分的個(gè)數(shù)M:');ot=input('請輸入時(shí)間增量ot的值：)；n=input('請輸入運(yùn)行次數(shù)n的值：)；ox=l/M;x0=zeros(M+1,1);f

35、orii=1:Mx0(ii+1)=ii*ox;endu=sin(pi*x0/l);%t=0時(shí)u(x,t)的值r=aA2*ot/(ox)A2;forii=1:n%數(shù)據(jù)的輸入B=zeros(M-1,1);%存放系數(shù)矩陣主對角線元素A=zeros(M-2,1);%存放系數(shù)矩陣主對角線元素下方次對角線的元素C=zeros(M-2,1);%存放系數(shù)矩陣主對角線元素上方次對角線的元素S=zeros(M-1,1);%存放右端的常數(shù)項(xiàng)forii=1:M-2B(ii)=1+2*r;A(ii)=-r;C(ii)=-r;S(ii)=u(ii+1,1);endB(M-1)=1+2*r;S(M-1)=u(M,1);u(

36、1,2)=0;u(M+1,2)=0;S(1,1)=S(1,1)+r*u(1,2);S(M-1,1)=S(M-1,1)+r*u(M+1,2);%追趕法S(1)=S(1)/B(1);T=B(1);k=2;whilek=MB(k-1)=C(k-1)/T;T=B(k)-A(k-1)*B(k-1);S(k)=(S(k)-A(k-1)*S(k-1)/T;k=k+1;endk=1;whilek=M-1S(M-1-k)=S(M-1-k)-B(M-1-k)*S(M-k);k=k+1;endu(2:M,2)=S;%把結(jié)果放入矩陣u中u(:,1)=u(:,2);%過河拆橋end%計(jì)算精確值，存放在u的第二列forx

37、=0:Mu(x+1,2)=exp(-(pi*a/l)A2*n*ot)*sin(pi*x*ox/l);end%計(jì)算最大相對誤差ez=zeros(M-1,1);forii=2:Mez(ii-1)=abs(u(ii,1)-u(ii,2)/u(ii,2);endE=max(ez);fprintf('最后時(shí)刻數(shù)值解與精確解分別為：n');disp(u);fprintf(,差分法得到的結(jié)果與正確結(jié)果的最大相對誤差為：,);disp(num2str(E*100)'%');%畫二維圖比較plot(x0,u(:,1),'g-',x0,u(:,2),'m*');legend('數(shù)值解',精確解)xlabel('x'),ylabel('u(x,t)')title('數(shù)值解與精確解比較')問題三溫度含水量擬合程序creatFit.mfunctionfitresult,gof=createFit(x,y)%CREATEFIT(X,Y)%Createafit.%Datafor'u