劉 洋，孫旭曙，黃葉寧

（三峽大學(xué) 水利與環(huán)境學(xué)院，湖北 宜昌 443002）

對大體積混凝土和凍土等溫度敏感材料研究，要求解溫度分布問題。而瞬態(tài)溫度問題求解方法[1-3]，一般用有限差分法推求某個(gè)時(shí)間時(shí)溫度場方程。但有限差分法進(jìn)行時(shí)間域離散時(shí)，數(shù)值解會(huì)在迭代過程出現(xiàn)震蕩情況，導(dǎo)致數(shù)值解精度下降。對有限差分法解的震蕩問題，林金木[4]應(yīng)用障礙迭代法在實(shí)數(shù)子空間求解瞬態(tài)溫度場，提高了解的精度，但較難找出帶約束的定解條件；李元清[5]針對柴油機(jī)瞬態(tài)溫度場問題，提出一種改進(jìn)的Crank-Nicholson法，但是這種方法在實(shí)際運(yùn)用中局限性大。

本文用一種新的方法[6]（Runge-Kutta法），并借助MATLAB分別編寫三階Runge-Kutta法和兩點(diǎn)循環(huán)的Crank-Nicholson法的有限元程序。且對單元熱容矩陣進(jìn)行處理，用集中質(zhì)量矩陣代替單元熱容矩陣。通過兩種方法說明，Runge-Kutta法相對于Crank-Nicholson法，在求解瞬態(tài)溫度場方程的優(yōu)勢，為Runge-Kutta法處理瞬態(tài)溫度場問題的理論提供依據(jù)。

1 數(shù)學(xué)模型及有限元求解方程

1.1 數(shù)學(xué)模型

考慮多孔介質(zhì)的各向同性，忽略水汽遷移、熱量對流的情況下，飽和土壤的二維溫度場方程為：

邊界條件（Dirichlet條件）：

式中T為溫度；T0為初始邊界溫度；C為體積熱容量，取1.68×104J/（m3·n）；λ 為導(dǎo)熱系數(shù)，取50W/（m·n）。

1.2 瞬態(tài)溫度場數(shù)值方法

瞬態(tài)溫度場數(shù)值方法，在空間域上用伽遼金法，在時(shí)間域上用Runge-Kutta法。首先，用伽遼金法對溫度場方程（1）進(jìn)行空間域離散可得：

由式（3）可得：

直接求大型矩陣M-1很困難，這里用到集中質(zhì)量矩陣[2]。假定單元的質(zhì)量集中在結(jié)點(diǎn)上，那么單元熱容矩陣Me轉(zhuǎn)換的單元集中質(zhì)量矩陣是對角線矩陣，且該質(zhì)量矩陣?yán)碚撋鲜钦ǖ摹?/p>

將單元熱容矩陣Me轉(zhuǎn)換為單元集中質(zhì)量矩陣：

接下來，用三階Runge-Kutta[11]法對式（4）進(jìn)行時(shí)間域離散。設(shè)t=n時(shí)，溫度為Tn；t=n+1時(shí)，溫度為Tn+1；其中Tn，1和Tn，2可由插值公式（6）和（7）表示：

式中Δt為時(shí)間步長，等于Tn+1與Tn差值；

則t=n+1，溫度Tn+1求解方程為：

2 二維穩(wěn)定溫度場解析解

矩形區(qū)域內(nèi)，兩邊（x=0，x=a）始終保持0度，另外兩邊（y=0，y=b）的溫度分別為f（x）和g（x）。

應(yīng)用分離變量法[7]，設(shè)T（x，y）=X（x）Y（y）：

其中λ 為常數(shù)，因此可得兩個(gè)常微分方程：

由齊次邊界條件（11）知：T（0，y）=0，T（a，y）=0，即X（0）=X（a）=0。首先，求解常微分方程邊值問題：X″（x）+λX（x）=0，X（0）=X（a）=0的非零解。

（1）當(dāng)λ≤0時(shí)，沒有非平凡解。

（2）當(dāng)λ＞0時(shí)，有非平凡解。

將λ 代入方程（14）可得：

其通解為：

可得出方程（9）滿足齊次邊界條件（11）的一系列特解：

由于方程（9）和邊界條件（11）是齊次的，因此：

仍滿足方程（13）和齊次邊界條件（14）。再應(yīng)用非齊次邊界條件：

則有關(guān)系式：

利用傅里葉系數(shù)公式得：

聯(lián)立式（20）和（21）得：

滿足方程（9）（10）（11）的通解為：

取a=1，b=2，f（x）=-6，g（x）=0，帶入式（25），此時(shí)穩(wěn)定溫度場解析解為：

3 結(jié)果分析

瞬態(tài)溫度場達(dá)到穩(wěn)定 （300次迭代）時(shí)，Crank-Nicholson法和Runge-Kutta法兩種數(shù)值方法得出模擬值與解析值進(jìn)行對比分析。

溫度場分布結(jié)果如圖1，圖2和圖3。

圖1 溫度等值線圖（解析解）

圖2 溫度等值線圖（C-N法）

圖3 溫度等值線圖（R-K法）

由圖1、圖2和圖3可知，數(shù)值解溫度分布情況和解析解溫度分布相同的。說明Crank-Nicholson法有限程序和Runge-Kutta法有限程序是正確的，且和解析值完全吻合。

取若干節(jié)點(diǎn)進(jìn)一步進(jìn)行溫度對比分析，如表1。最后對兩種數(shù)值方法進(jìn)行運(yùn)算效率對比，如表2。

表1 300次迭代兩種數(shù)值方法模擬結(jié)果

表2 300次迭代兩種數(shù)值運(yùn)算時(shí)間對比

從表1可知，瞬態(tài)溫度場運(yùn)行到穩(wěn)定時(shí)，數(shù)值解與解析解的誤差都在5%以內(nèi)，但相對于Crank-Nicholson法，Runge-Kutta法得的節(jié)點(diǎn)溫度精度更高，且提高2%左右。

從表2 中可知，瞬態(tài)溫度場運(yùn)算到穩(wěn)定時(shí)，Crank-Nicholson 法編寫MATLAB 程序運(yùn)算時(shí)長為1.531s，而Runge-Kutta法編寫MATLAB程序運(yùn)算時(shí)長0.991s，時(shí)長縮短0.54s，計(jì)算效率提高35%。顯然，相比Crank-Nicholson法，Runge-Kutta法大幅度提高運(yùn)算效率，節(jié)省了數(shù)值模擬時(shí)間。

4 結(jié)語

對瞬態(tài)溫度場進(jìn)行數(shù)值模擬，用Crank-Nicholson法和Runge-kutta法進(jìn)行時(shí)間域離散，對兩者運(yùn)行結(jié)果和運(yùn)算時(shí)間對比分析發(fā)現(xiàn)：

（1）Runge-Kutta法用于處理瞬態(tài)溫度場分布問題，結(jié)果非常滿意。

（2）不管是精度，還是運(yùn)算效率，Runge-kutta法都有明顯提高。

（3）用集中質(zhì)量矩陣取代單元熱容矩陣，結(jié)果較好。