Runge-Kutta法在瞬态温度场的应用 刘洋1¹，孙旭曙2¹，黄叶宁3¹

（1.三峡大学水利与环境学院 湖北 宜昌 443002）

摘要：对瞬态温度场求解常用差分法，但差分法会随着迭代过程而出现震荡，精度下降，效率不高，而Runge-Kutta法是一种特殊的单步法，精度高，效率高，广泛应用于求解常微分问题。本文用MATLAB对瞬态温度场分布问题求解的Crank-Nicholson法和Runge-Kutta法进行编程及实现，通过结果对比发现，Runge-Kutta法的计算精度不仅比Crank-Nicholson法高，且模拟效率也较后者的有显著提高，其中模拟效率提高幅度最大，达到35%。

关键词：瞬态温度场；Runge-Kutta法；Crank-Nicholson法.

0引 言

对大体积混凝土和冻土等对温度敏感材料研究，要求解温度分布问题。而瞬态温度问题求解方法^[1-3]，一般用有限差分法推求某个时间时温度场方程。但有限差分法进行时间域离散时，数值解会在迭代过程出现震荡情况，导致数值解精度下降。对有限差分法解的震荡问题，林金木^[4]应用障碍迭代法在实数子空间求解瞬态温度场，提高了解的精度，但较难找出带约束的定解条件；而本文用另外一种新方法（Runge-Kutta法），李元清^[5]针对柴油机瞬态温度场问题，提出一种改进的Crank-Nicholson法，但是这种方法在实际运用中局限性大。

而本文用一种新的方法^[6]（Runge-Kutta法），并借助MATLAB分别编写三阶Runge-Kutta法和两点循环的Crank-Nicholson法的有限元程序。且对单元热容矩阵进行处理，用集中质量矩阵代替单元热容矩阵。通过两种方法说明，Runge-Kutta法相对于Crank-Nicholson 法，在求解瞬态温度场方程的优势，为Runge-Kutta法处理瞬态温度场问题的理论提供依据。

1 数学模型及有限元求解方程

1.1 数学模型

考虑多孔介质的各向同性，忽略水汽迁移、热量对流的情况下，，饱和土壤的二维温度场方程为：

                                (1)

边界条件（Dirichlet条件）：                                  (2)

式中：T为温度；为初始边界温度；C为体积热容量，取； 为导热系数，取。

1.2瞬态温度场数值方法

瞬态温度场数值方法，在空间域上用伽辽金法，在时间域上用Runge-Kutta法。下面介绍这两种方法。

首先，用伽辽金法对温度场方程（1）进行空间域离散可得：

                       (3)

其中：单元热容矩阵

单元热传导矩阵

单元温度荷载列阵

将方程(3)的左项移右边为：                 （4）                    

直接求大型矩阵很困难，这里用到集中质量矩阵^[2]。假定单元的质量集中在结点上，那么单元热容矩阵转换的单元集中质量矩阵是对角线矩阵，且该质量矩阵理论上是正定的。

将单元热容矩阵转换为单元集中质量矩阵方法如下：

   (5)

因此，单元热容矩阵的逆为。

接下来，用三阶Runge-Kutta^[11] 法对式(4)进行时间域离散。 

设时，温度为，时，温度为，其中和可由插值公式(6)和（7）表示

   (6)

(7)

式中：时间步长为，它等于与差值;；。

那么，温度求解方程为                        (8)

2 二维稳定温度场解析解

矩形区域内，两边（，）始终保持0度，另外两边（，）的温度分别为和。

     （，）      (9)

，                     (10)

 ,                             (11)

应用分离变量法^[7]，设

                                        (12)

其中为常数，因此可得两个常微分方程

                                                (13)

由齐次边界条件(11)知

，，即

首先，求解常微分方程边值问题

                 ，的非零解       (14)

（1）当时，问题(14)没有非平凡解。

（2）当时，问题（14）有非平凡解。此时，，对应的。

接下来，考虑方程                                (15)

将代入方程（15）可得：                      (16)

其通解为： （n=1,2,3,…）                         (17)   

这样就得出方程(9)满足齐次边界条件(11)的一系列特解

     （n=1,2,3,…）            (18)

由于方程(9)和边界条件(11)是齐次的，因此

                               (19)

仍满足方程(13)和齐次边界条件(15)。

再应用非齐次边界条件

                       ，

则有关系式                   

利用傅里叶系数公式得：

                                            (20)

                                      (21）

联立式(20)和(21)得：

      (22)

                       (23）

因此，满足方程(9)(10)(11)的通解为

            (24)

取a=1,b=2,, 带入式(24)，此时稳定温度场解析解为：

      (25)

3 结果分析

瞬态温度场达到稳定（300次迭代）时，Crank-Nicholson法和Runge-Kutta法两种数值方法得出模拟值与解析值进行对比分析。

温度场分布结果如图1，图2和3所示。

图1 温度等值线图（解析解）     图2 温度等值线图（C-N法）        图3 温度等值线图（R-K法）

图1、2和3可知：数值解温度分布情况和解析解温度分布相同的。说明Crank-Nicholson法有限程序和Runge-Kutta法有限程序是正确的，且和解析值是完全吻合。

接下来，取若干节点进行温度对比分析，具体见表1

表1  300次迭代两种数值方法模拟结果

从表1可知，瞬态温度场运行到稳定时，数值解与解析解的误差都在5%以内，但相对于Crank-Nicholson法，Runge-Kutta法得的节点温度精度更高，且提高2%左右。

最后，对两种数值方法进行运算效率对比，见表2

表2  300次迭代两种数值运算时间对比

从表2中可知：瞬态温度场运算到稳定时，Crank-Nicholson法编写MATLAB程序运算时长为1.531s,而Runge-Kutta法编写MATLAB程序运算时长为0.991s,时长缩短0.54s，计算效率提高35%。显然，相比Crank-Nicholson法，Runge-Kutta法大幅度提高运算效率，节省了数值模拟时间。

4 结  论

对瞬态温度场进行数值模拟，用Crank-Nicholson法和Runge-kutta法进行时间域离散，对两者运行结果和运算时间对比分析发现：(1)Runge-Kutta法用于处理瞬态温度场分布问题，结果非常满意。(2)不管是精度，还是运算效率，Runge-kutta法都有明显提高。(3)用集中质量矩阵取代单元热容矩阵，结果较好。

参考文献：

[1]孔祥谦.有限元法在传热学中的应用.北京：科学出版社,1986.

[2]王勖成编著.有限单元法[M].北京：清华大学出版社,2003.

[3]R.Glowinski. Lectures on Numerical Methods for Non-linear Variational Problems. New York: Springer-Verlag, 1980.

[4]林金木.瞬态温度场的新解法[D].湖南：湖南大学学报,1996.

[5]李元清，高世义.有限元求解瞬态温度场的新方法.北京：内燃机学报,1990.

[6] 秦军.Runge-Kutta法在求解微分方程模型中的应用[D].合肥：安徽大学,2010

[7]吴崇试.数学物理方法[M].北京：北京大学出版社,2003,12.