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

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

2、算得出各材料層的溫度分布規(guī)律。問(wèn)題二，考慮在一些設(shè)備的支架不能固定在解凍土層上，必須固定在永凍土層中的情況下，地下土層的解凍位置，并給出解凍砂土與凍結(jié)砂土的分界線。由問(wèn)題一的求解可以計(jì)算出的值，在已知上界值與下界的溫度值后，同問(wèn)題一的求解方法可以給出解凍砂土與凍結(jié)砂土的分界線。問(wèn)題三，考慮結(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元。問(wèn)題四，考慮在以上問(wèn)題的基礎(chǔ)上，結(jié)合我國(guó)青藏鐵路永凍土層地基進(jìn)行仿真，并為施工單位提出合理建議。因?yàn)楸绢}的前三問(wèn)即是在查閱青藏鐵路路基修建相關(guān)數(shù)據(jù)的基礎(chǔ)上進(jìn)行的，故問(wèn)題四的仿真即已經(jīng)得到相應(yīng)的解決。通過(guò)對(duì)以上問(wèn)題的求解進(jìn)行合理性分析，即可對(duì)施工單位給出合理建議。為了簡(jiǎn)化計(jì)算量，提高求解速度，本題中的微分方程模型使用向前差分的方法求出近似數(shù)值解，而且對(duì)模型的可行性及有效性進(jìn)行了一定的分析，所得結(jié)果十分合理。本文的優(yōu)點(diǎn)在于利用差分方法求解一維熱傳導(dǎo)微分方程模型的近似數(shù)值解，使得材料內(nèi)界面的條件處理得較為容易。同時(shí)，在前三個(gè)問(wèn)題的求解中，將問(wèn)題背景設(shè)定在青藏鐵路的修建中，在一定程

4、度上對(duì)問(wèn)題四的求解提供了較大的幫助。關(guān)鍵詞：熱傳導(dǎo)問(wèn)題 拋物型方程 數(shù)值模擬1 問(wèn)題重述在永凍土地上鋪設(shè)道路、飛機(jī)跑道和某些結(jié)構(gòu)的地基。分析這類(lèi)地基的結(jié)構(gòu)，有瀝青層、混凝土層、干砂層、石子層、絕緣材料層，再下面是濕的沙土層和凍土層（參見(jiàn)圖1）。已知，外界的溫度是關(guān)于時(shí)間的函數(shù)，凍土層的溫度是不變的零下溫度。請(qǐng)回答以下問(wèn)題：L0L1L2L3L4L5瀝青層混凝土層干砂石層防水隔溫層半凍砂土層外界溫度（t）K1K2K3K4K5永凍土層溫度T0圖1：永凍土地地基結(jié)構(gòu)圖（1）分析空氣溫度傳入路基規(guī)律，各材料層的溫度分布；（2）由于一些設(shè)備的支架不能固定在解凍土層上，必須固定在永凍土層中，因此確定地下土層

5、的解凍位置非常重要，請(qǐng)給出解凍砂土與凍結(jié)砂土的分界線；（3）請(qǐng)給出各層材料的最佳厚度（結(jié)合溫度分布、成本和耐用性 ）；（4）請(qǐng)結(jié)合我國(guó)青藏鐵路永凍土層地基進(jìn)行仿真，并給施工單位給出合理建議。2 問(wèn)題分析本問(wèn)題要求建立模型和設(shè)計(jì)算法，在永凍土地上鋪設(shè)道路、飛機(jī)跑道和某些結(jié)構(gòu)的地基的問(wèn)題背景中，結(jié)合溫度分布、成本和耐用性，給出各層材料的最佳厚度，并對(duì)我國(guó)青藏鐵路永凍土層地基進(jìn)行仿真。由于問(wèn)題的局限性，我們首先應(yīng)該找到合理的數(shù)據(jù)來(lái)對(duì)相關(guān)背景進(jìn)行一定的了解，定位模型的求解方法，循序漸進(jìn)地對(duì)此問(wèn)題進(jìn)行合理的求解。問(wèn)題一中，要求在題中所給的情況下，分析空氣溫度傳入路基規(guī)律，以及各材料層的溫度分布。經(jīng)過(guò)分析

6、，我們發(fā)現(xiàn)建立模型所需的數(shù)據(jù)不足，故應(yīng)首先搜集題目中各材料熱傳導(dǎo)相關(guān)的密度、導(dǎo)熱系數(shù)、比熱容以及環(huán)境溫度等數(shù)據(jù)。問(wèn)題二中，要求確定地下土層的解凍位置，給出解凍砂土與凍結(jié)砂土的分界線，以應(yīng)對(duì)由于一些設(shè)備的支架不能固定在解凍土層上，必須固定在永凍土層中的情況。由問(wèn)題一建立的模型，本問(wèn)題即是求出溫度為0時(shí)所對(duì)應(yīng)的距離。問(wèn)題三中，要求結(jié)合溫度分布、成本和耐用性，給出各層材料的最佳厚度。在分析問(wèn)題一、二后，則是對(duì)各目標(biāo)進(jìn)行權(quán)重約束，以求得最佳的厚度。問(wèn)題四中，要求結(jié)合我國(guó)青藏鐵路永凍土層地基進(jìn)行仿真，并給施工單位給出合理建議。3 模型假設(shè)1、假設(shè)環(huán)境溫度不會(huì)發(fā)生突變，沒(méi)有極端天氣出現(xiàn)； 2、假設(shè)路基材料

7、分布均勻；3、假設(shè)考慮壓實(shí)度時(shí)，可以將半凍砂土層以上的四層材料合并為一層；4、假設(shè)對(duì)于各層只考慮重力所產(chǎn)生的壓強(qiáng)。4 符號(hào)說(shuō)明外界的溫度關(guān)于時(shí)間的函數(shù)凍土層的溫度每個(gè)材料層的熱傳導(dǎo)率每個(gè)材料層的比熱每個(gè)材料層的密度要求解的第個(gè)材料層的溫度分布地下某一點(diǎn)距地面的距離時(shí)刻值每層材料的厚度其他局部符號(hào)在引用時(shí)給出了具體說(shuō)明。5 模型的建立與求解5.1模型一的建立與求解問(wèn)題一中要求分析空氣溫度傳入路基規(guī)律，以及各材料層的溫度分布。我們知道，一切穩(wěn)定的數(shù)學(xué)物理問(wèn)題都可以用橢圓型微分方程模型來(lái)描述，若在穩(wěn)態(tài)之前有一個(gè)過(guò)程，則討論這種漸進(jìn)穩(wěn)定的數(shù)學(xué)物理過(guò)程可以用拋物型為方程模型來(lái)描述1。5.1.1模型一的建

8、立分析題目可知，地基的結(jié)構(gòu)有瀝青層，混凝土層，干砂石層，防水隔溫層，半凍砂土層，外界的溫度是關(guān)于時(shí)間的函數(shù)，凍土層的溫度是不變的零下溫度，每個(gè)材料層都有熱傳導(dǎo)率，問(wèn)題一則要分析空氣溫度傳入路基的規(guī)律，以及各材料層的溫度分布2。首先，由于路基各材料層是均勻的，所以要分析的熱傳導(dǎo)方程可歸為一維的熱傳導(dǎo)方程研究。路基各材料層則化為一維熱傳導(dǎo)方程為 （1）其中，是要求解的第個(gè)材料層的溫度分布；，為此材料層的傳熱系數(shù)，比熱和密度，它們都是已知的常數(shù)，且。由模型假設(shè)可知，層與層的界面處傳導(dǎo)率處于平衡，且溫度相等，即 （2）基于以上分析，在半凍砂土層假設(shè)為濕砂層（是其傳熱系數(shù)）的情況下，所建立的路基熱傳導(dǎo)微

9、分方程模型為 （3） 5.1.2模型一的求解由以上所建立的拋物型微分方程模型進(jìn)行分析，對(duì)問(wèn)題一的求解我們給出了如下的具體步驟：（1）依據(jù)題意，我們?cè)谥袊?guó)選礦技術(shù)網(wǎng)中找到瀝青層，混凝土層，干砂石層，防水隔溫層和半凍砂土層的密度，比熱容和熱導(dǎo)系數(shù)等值，部分難以確定的材料選取了最接近的材料進(jìn)行替代，所找到的數(shù)據(jù)材料如下表。表1 各層材料相關(guān)參數(shù)表各材料層干密度熱導(dǎo)系數(shù)比熱容 價(jià)格瀝青層（石油瀝青）14000.271.682800混凝土層（c30混凝土）23001.510.92770.5干砂石16000.581.0170防水隔溫層（乳化瀝青膨脹珍珠巖）4000.121.551000砂土層（濕砂土）23

10、501.421.07500以上則為各材料層的密度、熱導(dǎo)系數(shù)、比熱容和價(jià)格（價(jià)格因素在問(wèn)題三中需要考慮，故在這里一起給出）的數(shù)據(jù)。此外，由于溫度的傳導(dǎo)涉及到凍土地帶的環(huán)境溫度，故我們結(jié)合青藏鐵路的修建，選取拉薩最近24小時(shí)內(nèi)的溫度3?？紤]到中國(guó)天氣網(wǎng)只可以給出最近24小時(shí)的較精確溫度，故我們問(wèn)題的季節(jié)背景選擇在夏秋之交4。表2 拉薩24小時(shí)各時(shí)刻溫度表時(shí)刻000102030405060708091011溫度121212131313131313121313時(shí)刻121314151617181920212223溫度151718192122222120161312以上則為拉薩最近24小時(shí)各時(shí)刻的溫度表。

11、（2）通過(guò)查閱鐵路建設(shè)國(guó)家規(guī)范，我們將瀝青層的厚度設(shè)定為4cm，混凝土層設(shè)定為10cm，干砂石層設(shè)定為7cm，防水隔溫層設(shè)定為10cm5。（3）通過(guò)擬合數(shù)據(jù)得到外界溫度關(guān)于時(shí)間的關(guān)系式，根據(jù)問(wèn)題一中建立的第一種邊界的一維拋物型方程，考慮到解析解不易求得，故在這里我們利用差分的方法求出其近似的數(shù)值解6。差分求解拋物型偏微分方程的過(guò)程則是首先對(duì)平面進(jìn)行網(wǎng)格剖分，分別取，為方向和方向的步長(zhǎng)，用兩族平行直線，將平面剖分成矩形網(wǎng)格，節(jié)點(diǎn)為，為簡(jiǎn)單起見(jiàn)，記，則在網(wǎng)格內(nèi)點(diǎn)處，對(duì)采取向前差商公式，對(duì)采用二階中心差商公式，可得到一維熱傳導(dǎo)方程的差分近似： （4）（4）求出第一層底部的溫度后，令，依據(jù)公式（2）即

12、可求得下一層底部的溫度，以下每一層材料溫度的求解依次循環(huán)以上步驟。5.1.3模型一的結(jié)果及分析針對(duì)5.1.2中的求解步驟，運(yùn)用Matlab編程求解。針對(duì)附錄一（1）中的程序文件，運(yùn)行可得到如下結(jié)果。外界環(huán)境溫度隨時(shí)間的變化函數(shù)為：（5）外界溫度隨時(shí)間變化的函數(shù)如下圖： 圖1 外界溫度隨時(shí)間變化的函數(shù)第一層求解所得的溫度隨距離或時(shí)間變化的圖像分別為：圖2 第一層溫度隨距離分布圖圖3 第一層溫度隨時(shí)間分布圖其他層溫度隨時(shí)間，隨距離的變化圖像見(jiàn)附錄一（2）。由上可知，當(dāng)固定另一因素時(shí)，溫度隨著與路基表面距離的增加越來(lái)越小，其隨時(shí)間的變化則呈現(xiàn)先增加后減少的趨勢(shì)。5.2模型二的建立與求解根據(jù)題意，問(wèn)題

13、二中要求確定地下土層的解凍位置，給出解凍砂土與凍結(jié)砂土的分界線，以應(yīng)對(duì)由于一些設(shè)備的支架不能固定在解凍土層上，必須固定在永凍土層中的情況。由問(wèn)題一建立的模型，本問(wèn)題即是求出0時(shí)所對(duì)應(yīng)的距離。5.2.1問(wèn)題二模型的建立根據(jù)問(wèn)題一中所建立的模型，即可得出第五層溫度隨著距離變化的圖像關(guān)系。由圖分析，即可給出0時(shí)所對(duì)應(yīng)的距離7。5.2.2問(wèn)題二模型的求解問(wèn)題二中建立的模型用來(lái)確定第五層的解凍位置，即在已知上界的值與下界溫度后，給出解凍砂土與凍結(jié)砂土的分界線。由函數(shù)圖象分析，即可給出0時(shí)所對(duì)應(yīng)的距離。5.2.3問(wèn)題二模型的結(jié)果分析針對(duì)5.2.2中的求解步驟，依據(jù)問(wèn)題一中所得的第四層溫度隨距離的圖像分析進(jìn)

14、行求解。所得圖像如下圖：圖4 第五層溫度隨距離變化圖像解得解凍砂土與凍結(jié)砂土之間的分界線為路面下120cm-140cm。結(jié)合國(guó)家公路施工技術(shù)規(guī)范，分析數(shù)據(jù)可知所得結(jié)果比較合理，故可作為問(wèn)題而的求解答案。5.2.4模型二的誤差分析因?yàn)椴罘址ㄊ菍?duì)數(shù)值解得近似，所以對(duì)最后結(jié)果的誤差分析十分重要。依據(jù)附錄二中的程序文件，得出如下的誤差分析圖像：圖5 數(shù)值解與精確解比較圖像由上可知，數(shù)值解與精確解的最大相對(duì)誤差非常小，為0.42877% ，且數(shù)值解的求解結(jié)果是穩(wěn)定收斂的，故該結(jié)果可以作為本問(wèn)題的求解。5.3模型三的建立與求解由題目可知，問(wèn)題三需要在問(wèn)題一、二的求解啟發(fā)中，結(jié)合溫度分布、成本和耐用性，給出

15、各層材料的最佳厚度。在分析問(wèn)題一、二后，則是對(duì)各目標(biāo)進(jìn)行權(quán)重約束，以求得最佳的厚度。5.3.1問(wèn)題三模型的建立由于題目中所給相關(guān)數(shù)據(jù)較少，由此我們思考對(duì)最佳厚度設(shè)置約束條件。對(duì)于耐用性，我們從壓實(shí)度和壓強(qiáng)兩方面進(jìn)行考慮。首先將上面四層材料合為一層以簡(jiǎn)化計(jì)算，查閱文獻(xiàn)可知： 而實(shí)際干密度則與溫度、水分等條件相關(guān)，則在此時(shí)將溫度分布考慮進(jìn)模型的建立。由此，需要擬合溫度與水分的關(guān)系8。而對(duì)于壓強(qiáng)，我們假設(shè)這里路基的承重來(lái)源只有上部材料的重力，則可得出如下推導(dǎo)公式： （6）對(duì)上述公式中出現(xiàn)的局部符號(hào)說(shuō)明如下：表示某層材料所承受的壓強(qiáng)；表示某層材料所承受的壓力；表示承壓接觸面的表面積； 表示表面積對(duì)應(yīng)的

16、材料的體積； 表示體積對(duì)應(yīng)的材料的質(zhì)量；表示材料的密度；表示材料層的厚度。將目標(biāo)函數(shù)設(shè)定為壓實(shí)度大于95%且能承受最大壓強(qiáng)為2760MPa時(shí)求得的最小成本。用表示成本， （7）同時(shí)添加對(duì)各層壓強(qiáng)的限制，使各層所承受壓強(qiáng)大于路基建設(shè)所應(yīng)承載的最小壓強(qiáng)值。同時(shí)，分析數(shù)據(jù)知， （8）基于以上分析，用表示含水率，以（7）為目標(biāo)函數(shù)，（8）為約束條件建立優(yōu)化模型如下： （9）5.3.2問(wèn)題三模型的求解由于問(wèn)題三是在結(jié)合溫度分布、成本、耐用性多種因素的基礎(chǔ)上對(duì)厚度問(wèn)題進(jìn)行的研究，所以需要考慮各個(gè)因素的相關(guān)數(shù)據(jù)，以及各個(gè)因素之間的相互關(guān)系，具體的求解步驟如下：（1）查閱資料對(duì)溫度與水分含量的關(guān)系給出合理的擬

17、合關(guān)系，并得出具體函數(shù)公式；（2）對(duì)壓強(qiáng)和壓實(shí)度的影響因素進(jìn)行化簡(jiǎn)，得到相應(yīng)的約束關(guān)系；（3）列出壓實(shí)度大于95%時(shí)的目標(biāo)函數(shù)，在Lingo中求解線性規(guī)劃問(wèn)題即可。5.3.3問(wèn)題三模型的結(jié)果分析針對(duì)5.3.1中的模型，按照5.3.2中的求解步驟，首先利用Matlab求解得出溫度隨水分變化的函數(shù)關(guān)系式為： （10）所畫(huà)出的函數(shù)圖像為：圖6 溫度與水分的擬合圖像擬合度數(shù)據(jù)如下：gof = sse: 0.6035 rsquare: 0.9991 dfe: 5 adjrsquare: 0.9985 rmse: 0.3474有數(shù)據(jù)可以看出，溫度與水分的函數(shù)關(guān)系擬合得較為合理。對(duì)所列出的目標(biāo)函數(shù)與約束條件

18、，運(yùn)用Lingo編程求解，解得壓實(shí)度大于95%時(shí)的目標(biāo)函數(shù)所取得最小值為46334.72元。此時(shí)，第一、二、三、四層的厚度分別為5cm，12cm，6cm，10cm。同上分析，可以看出此數(shù)據(jù)與實(shí)際施工情況中的數(shù)據(jù)相差較小，故可作為本題的求解。5.4模型四的建立與求解由題目可知，問(wèn)題四要求結(jié)合我國(guó)青藏鐵路永凍土層地基進(jìn)行仿真，并給施工單位給出合理建議。在對(duì)之前問(wèn)題求解的基礎(chǔ)上，得出合理的分析與建議。5.4.1問(wèn)題四模型的建立問(wèn)題一、二中選取的數(shù)據(jù)即是青藏鐵路路基建設(shè)的相關(guān)數(shù)據(jù)，同時(shí)將外界環(huán)境溫度設(shè)定為拉薩市最近24小時(shí)的溫度。故由此所得的厚度范圍即為對(duì)青藏鐵路永凍土層地基建設(shè)的仿真9。5.4.2問(wèn)

19、題四模型的求解由上所述原理對(duì)問(wèn)題四進(jìn)行求解，具體求解步驟如下：（1）由以上三個(gè)問(wèn)題的求解，給出各材料層最佳厚度取值表；（2）通過(guò)對(duì)所得結(jié)果進(jìn)行分析，給出施工工人的合理化建議。5.4.3問(wèn)題四模型的結(jié)果分析由以上問(wèn)題的求解，可以得出如下結(jié)果：表3 各材料層最佳厚度取值表材料層厚度瀝青層（石油瀝青）5cm混凝土層（c30混凝土）12cm干砂石6cm防水隔溫層（乳化瀝青膨脹珍珠巖）10cm通過(guò)對(duì)數(shù)據(jù)與國(guó)家青藏鐵路路基建設(shè)進(jìn)行比對(duì)，以上所得求解數(shù)據(jù)基本符合鐵路建設(shè)規(guī)范，故可以作為本問(wèn)題的仿真結(jié)果。5.4.3問(wèn)題四模型求解對(duì)施工工人的建議（1）由以上求解可以給出各層材料大致的厚度取值，可將此作為一個(gè)基本

20、的施工參照給施工工人一定的理論指導(dǎo)；（2）由于模型建設(shè)中提出了相應(yīng)的假設(shè)，故在實(shí)際問(wèn)題的施工時(shí)，相關(guān)施工單位可結(jié)合自身施工經(jīng)驗(yàn)對(duì)厚度取值進(jìn)行一定的調(diào)整，同時(shí)也應(yīng)進(jìn)行多次項(xiàng)目預(yù)試驗(yàn)，以充分考慮鐵路承重、修建成本以及后發(fā)事故。6 模型的評(píng)價(jià)與改進(jìn)方向6.1 模型的評(píng)價(jià)由于本問(wèn)題中，模型建立的背景為青藏鐵路的修建，同時(shí)對(duì)于各材料的相關(guān)熱學(xué)參數(shù)進(jìn)行了唯一取值，對(duì)于此情況下的鐵路路基熱傳導(dǎo)問(wèn)題可以合理有效地解決。如果修建環(huán)境發(fā)生改變便，則各材料層成分以及環(huán)境溫度均會(huì)發(fā)生改變，故應(yīng)對(duì)該模型進(jìn)行改進(jìn)繼而求得新的求解情況。6.1.1 模型的優(yōu)點(diǎn)（1）本模型中，所選取的數(shù)據(jù)均為與各材料層成分接近程度最大的材料的

21、熱學(xué)參數(shù)，建立的線性模型函數(shù)設(shè)置合理，對(duì)于題中所給問(wèn)題可以得到滿意的求解方案。（2）算法設(shè)計(jì)較為合適，計(jì)算量較小，能夠高效運(yùn)行出結(jié)果。（3）利用差分代替微分的思維求解一維熱傳導(dǎo)微分方程模型，可以使材料內(nèi)界面的條件處理得較為容易。（4） 本模型在前三個(gè)問(wèn)題的求解中，將問(wèn)題背景設(shè)定在青藏鐵路的修建中，在一定程度上對(duì)問(wèn)題四的求解提供了較大的幫助。6.1.2 模型的缺點(diǎn)（1）問(wèn)題一中，我們所選取的材料數(shù)據(jù)并不能充分考慮各層材料成分，溫度數(shù)據(jù)也僅為一個(gè)地點(diǎn)24小時(shí)內(nèi)的數(shù)據(jù)；（2） 問(wèn)題三中，在考慮壓實(shí)度時(shí)，將半凍砂土層以上材料壓縮為一種材料，在一定程度上簡(jiǎn)化了問(wèn)題模型。6.2 模型的改進(jìn)方向（1）對(duì)于所

22、選取的溫度及材料熱學(xué)數(shù)據(jù)，可以選取多組求平均值，以使求解結(jié)果更具代表性；（2）求解一維熱傳導(dǎo)微分方程時(shí)，對(duì)于初始條件的選取稍有欠缺，可以結(jié)合實(shí)際情況選擇更好的初值條件。（3）對(duì)于差分求解一維熱傳導(dǎo)微分方程，雖在一定程度上減少了求解困難，但其步長(zhǎng)選取不太合理，可以進(jìn)一步研究問(wèn)題模型，以選取更合理地步長(zhǎng)，以使求解結(jié)果與真實(shí)值更加接近。參考文獻(xiàn)1 張辰熙,于明.季節(jié)凍土層下人工凍土墻溫度場(chǎng)實(shí)驗(yàn)研究J.低溫建筑 術(shù).2010(6):144-146.2 文斌,吳青柏,蔣觀利,張鵬.模擬退火優(yōu)化算法的凍土熱傳導(dǎo)參數(shù)反分析J.巖 土力學(xué).2013(8):15-16.3 王建偉,李書(shū)民，李振卿.凍土層等效導(dǎo)熱

23、系數(shù)的探討J.防滲技術(shù).2001(3):36-38.4 李述訓(xùn),吳通化.凍土溫度狀況研究方法和應(yīng)用分析J.冰川凍土.2004(4)： 142-148.5 候彬彬,盧德唐,陳永輝,孔祥言.非達(dá)西滲流模型模擬青藏鐵路拋石路基凍土層溫度J.水工滲流研究與應(yīng)用進(jìn)展第五屆全國(guó)水利工程滲透學(xué)術(shù)研討會(huì)論文集.2006.6 史冊(cè).熱傳導(dǎo)方程有限差分法的MATLAB實(shí)現(xiàn)J.咸陽(yáng)師范學(xué)院學(xué)報(bào).2009(4):24-26.7 張世民.青藏鐵路多年凍土路基熱-力穩(wěn)定性數(shù)值仿真分析J.土木建筑與環(huán)境工程.2012(34):3-7.8 穆彥虎,馬巍,牛富俊,李國(guó)玉,王大雁,劉永智.青藏鐵路多年凍土區(qū)普通路基熱狀況監(jiān)測(cè)分析

24、J.冰川凍土.2014(2):10-13.9 田立慧，凌賢長(zhǎng),王力娜,張鋒.青藏鐵路高溫多年凍土區(qū)列車(chē)行駛路基長(zhǎng)期永久變形數(shù)值模擬研究J.地震工程學(xué)報(bào).2014(4):21-25.附錄問(wèn)題一外界環(huán)境溫度隨時(shí)間變化的函數(shù)關(guān)系式擬合程序creatFit.mfunction fitresult, gof = createFit(x, y)%CREATEFIT(X,Y)% Create a fit.% Data for 'temperature' fit:% X Input : x% Y Output: y% Output:% fitresult : a fit object repr

25、esenting the fit.% gof : structure with goodness-of fit info.% 另請(qǐng)參閱 FIT, CFIT, SFIT.% 由 MATLAB 于 17-Aug-2015 19:27:00 自動(dòng)生成% Fit: 'temperature'.xData, yData = prepareCurveData( x, y );% Set up fittype and options.ft = fittype( 'fourier2' );opts = fitoptions( 'Method', 'Non

26、linearLeastSquares' );opts.Display = 'Off'opts.Normalize = 'on'opts.StartPoint = 0 0 0 0 0 1.9316882339819;% Fit model to data.fitresult, gof = fit( xData, yData, ft, opts );% Plot fit with data.figure( 'Name', 'temperature' );h = plot( fitresult, xData, yData );l

27、egend( h, 'y vs. x', 'temperature', 'Location', 'NorthEast' );% Label axesxlabel( 'x' );ylabel( 'y' );grid onnihe.mx=1:24;y=12 12 13 13 13 13 13 13 12 13 13 15 17 18 19 21 22 22 21 20 16 13 12 12;fitresult, gof = createFit(x, y)偏微分方程求解程序function U x t=

28、PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C)%古典顯式格式求解拋物型偏微分方程%U x t=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C)%方程：u_t=C*u_xx 0 <= 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 -

29、時(shí)間變量%輸入?yún)?shù)：uX -空間變量x的取值上限% uT -時(shí)間變量t的取值上限% phi -初值條件，定義為內(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');%U x t=PDEParabolicClassicalExplicit

30、(uX,uT,phi,psi1,psi2,M,N,C);% inline('263+15.11+2.216*cos(1.757*t)+3.687*sin(1.757*t)-1.895*cos(2*1.757*t)+1.728*sin(2*1.757*t)')%設(shè)置參數(shù)C的默認(rèn)值if nargin=7 C=1;end%計(jì)算步長(zhǎng)dx=uX/M;%x的步長(zhǎng)dt=uT/N;%t的步長(zhǎng)x=(0:M)*dx;t=(0:N)*dt;r=C*dt/dx/dx;%步長(zhǎng)比r1=1-2*r;if r > 0.5 disp('r > 0.5,不穩(wěn)定')end%計(jì)算初值和邊值

31、U=zeros(M+1,N+1);for i=1:M+1 U(i,1)=phi(x(i);endfor j=1:N+1 U(1,j)=psi1(t(j); U(M+1,j)=psi2(t(j);end%逐層求解for j=1:N for i=2:M U(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;

32、溫度隨時(shí)間，隨距離變化的圖像圖7 第二層溫度隨時(shí)間的變化圖像圖8 第二層溫度隨距離的變化圖像圖9 第二層溫度隨時(shí)間的變化圖像圖10 第三層溫度隨距離的變化圖像圖11 第四層溫度隨時(shí)間的變化圖像圖12 第四層溫度隨距離的變化圖像問(wèn)題二模型二誤差分析程序：% 求土壤溫度分布clear;clc;format short ea=input(' 請(qǐng)輸入系數(shù)a 的值：');l=input(' 請(qǐng)輸入長(zhǎng)度l 的值：');M=input(' 請(qǐng)輸入將區(qū)間0,l等分的個(gè)數(shù)M：');ot=input(' 請(qǐng)輸入時(shí)間增量ot 的值：');n=input

33、(' 請(qǐng)輸入運(yùn)行次數(shù)n 的值：');ox=l/M;x0=zeros(M+1,1);for ii=1:Mx0(ii+1)=ii*ox;endu=sin(pi*x0/l);%t=0 時(shí)u(x,t)的值r=a2*ot/(ox)2;for ii=1:n%數(shù)據(jù)的輸入B=zeros(M-1,1);%存放系數(shù)矩陣主對(duì)角線元素A=zeros (M-2,1);%存放系數(shù)矩陣主對(duì)角線元素下方次對(duì)角線的元素C=zeros (M-2,1);%存放系數(shù)矩陣主對(duì)角線元素上方次對(duì)角線的元素S=zeros(M-1,1);%存放右端的常數(shù)項(xiàng)for ii=1:M-2B(ii)=1+2*r;A(ii)=-r;C(i

34、i)=-r;S(ii)=u(ii+1,1);endB(M-1)=1+2*r;S(M-1)=u(M,1);u(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;while k=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;while k=M-1S(M-1-k)=S(M-1-k)-B(M-1-k)*S(M-k);k=k+1;endu(2:M,

35、2)=S; %把結(jié)果放入矩陣u 中u(:,1)=u(:,2);% 過(guò)河拆橋end%計(jì)算精確值，存放在u 的第二列for x=0:Mu(x+1,2)=exp(-(pi*a/l)2*n*ot)*sin(pi*x*ox/l);end%計(jì)算最大相對(duì)誤差ez=zeros(M-1,1);for ii=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é)果的最大相對(duì)誤差為：');disp(num2str

36、(E*100) '%');%畫(huà)二維圖比較plot(x0,u(:,1),'g-',x0,u(:,2),'m*');legend(' 數(shù)值解',' 精確解')xlabel('x'),ylabel('u(x,t)')title(' 數(shù)值解與精確解比較')問(wèn)題三溫度含水量擬合程序creatFit.mfunction fitresult, gof = createFit(x, y)%CREATEFIT(X,Y)% Create a fit.% Data for 'untitled fit 1' fit:% X Input : x% Y Output: y% Output:% fitresult : a fit object representing the fit.% gof : structure with goodness-of fit info.% 另請(qǐng)參閱 FIT, CFIT, SFIT.% 由 MATLAB 于 19-Aug-2015 09:26:16 自動(dòng)生成% Fit: