基于有限差分法的瞬态温度场计算误差研究_张涛

以双曲线型绝热温升 θ (τ ) = θ 0
和向后差分中点龄期方法求解精度较高，当绝热温升速率值较小，在 0.5 天步长下，基本能够保证工程 应用精度。但当绝热温升速率值较大，且步长进一步加大时，误差也是不容忽视的。如当 L =1，时间 步长为 1 天时，向后差分中点龄期方法计算最终绝热温升将产生 0.07 θ 0 的误差；中点差分方法计算最 终绝热温升将产生 0.14 θ 0 的误差。 事实上，以差分方法计算水化热增量，必然会造成水化热温升计算的误差。其大小与绝热温升变 收稿日期:2003-03-07 录用日期:2003-05-07 作者简介:张涛(1972-),女,工程师,结构工程博士 研究方向:结构有限元分析及设计 1
张涛 1 张宇鑫 2 （1 上海市政工程设计研究院轨道所，上海, 200092；2 上海师范大学建筑工程学院，上海, 201418）
摘要:本文研究了瞬态温度场有限差分计算中时间步长对计算精度的影响，提出了有效的修正方法，并对修正龄期 进行了分析，为数值模拟提供了可靠的理论依据。研究结果表明，文中所提出的方法对减少由于时间步长增大而造成的 计算误差是简明而有效的。 关键词：温度场 时间步长 水化热 计算误差 中图分类号：TU375 文献标识码：A
f (τ n ) + f (τ n − ∆τ n ) ∆τ n ， 2
f (τ n ) ∆τ n 。由于向前差分方法的差分步长容易受计算稳定性限制，现已很少采用。中点差分和向后差分
方法均是隐式解法，计算都是无条件稳定的。目前采用向后差分法较多，中点差分法精度相对较高， 但有时求解会产生振荡。
τ 为例， 后差分法和中点差分法对绝热温升计算结果的精度影响如 L +τ 图 1 所示：直接应用后点龄期计算的向后差分方法计算误差是非常严重的， L =1 时，以 0.5 天为步长， 最终绝热温升计算结果将产生 0.13 θ 0 的误差；以 1 天为步长，将产生 0.27 θ 0 的误差。中点差分方法
１．２ 计算绝热温升（θ０＝１） １．１ １ ０．９ ０．８ ０．７ ０．６ ０．５ ０．４ 向后差分后点龄期计算，步长0.5天 向后差分后点龄期计算，步长1 天 时间（ｄ） 中点差分，步长1天 中点差分,步长0.5天 精确值 向后差分中点龄期计算，步长0.5 向后差分中点龄期计算，步长1天 计算绝热温升（θ０＝１） １．１ １ ０．９ ０．８ ０．７ ０．６ ０．５ ０．４ 精确值 向后差分中点龄期计算，步长 0.5 天 向后差分中点龄期计算，步长 1 天 向后差分后点龄期计算，步长 0.5 天 向后差分后点龄期计算，步长 1 时间（ｄ） 天 中点差分，步长1 天 中点差分，步长0.5 天
4 结论
通过分析，可得出以下主要结论：(1)在用有限差分法求解混凝土瞬态温度场时，求解精度与差分 方法、差分步长及温度变化速率都有很大关系。直接采用向后差分中点龄期计算方法将引起结果的较 大误差。(2)对于浇筑初期的混凝土温度场，因温度变化较快，故由差分步长引起的计算误差较大，初 期修正系数较大，随温度变化逐渐平缓，修正系数最终趋于 1。(3)本文提出的修正方法，能够有效减 小由于增大时间步长给计算带来的误差影响，同时保证绝热温升总量不变。
在瞬态温度场有限差分计算中，时间步长的选定直接影响着计算规模的大小。选定大的计算步长 可以有效减少计算量，但同时会引起较大的计算误差。国内外文献中[1-4]，对温度场计算中的网格数量 和计算次数如何减少进行了一些研究。其中一些算法和理论虽然在不同程度上提高了运算速度，但也 引起了一些负效应，有些模型涉及到过渡单元的使用和反复调整，使并层在程序处理上遇到许多麻烦， 数据交换的自动化程序设计有较大难度，反而大大增加了前处理的工作量，在实际应用上受到了限制。 因而加大时间步长仍是一种提高温度场运算速度的最直接方法。本文将针对时间步长过大造成的误差 进行研究，以寻求简单实用的修正方法，从而有效提高运算速度和精度。
第4期 2003 年 12 月
广 东 交 通 职 业 技 术 学 院 学 报
Journal Of Guang Dong Communication Polytechnic
NO.4 December 2003
文章编号：1671-8496(2003)04-0001-03
基于有限差分法的瞬态温度场计算误差研究
β (T − Ta ) λ
其中 α 为导温系数,θ为绝热温升，T 为温度；Γ1 , Γ2 , Γ3 为一、二、三类边界； T0 为边界给定的温度；q 为 给定的边界热流量； β 为放热系数； λ 为热流系数。 方程(1)的求解通常采用有限元-差分方法，即在空间域用有限元方法进行离散，在时间域用差分法 进行离散。求解总体方程为[6]： 1− s 1 1− s [H ] + {T } + {T } + [R ] [H ] − 1 [R ] {F } + {F } = 0 (2)
Study on Computing Error of Transient Temperature Field Based on Finite Difference Method
ZHANG Tao，2ZHANG Yu-xin ( Shanghai Municipal Engineering Design Institute,Shanghai 200092, China; 2School of Architecture Engineering of Shanghai Normal University, Shanghai 201418, China)
∆τ 2 1 3 1 )∆τ2 =[ f (((m−1)k + )∆τ1) + f (((m−1)k + )∆τ1) +...+ f ((mk− )∆τ1)]∆τ1 (6) 2 2 2 2 为以 ∆τ 2 为时间步长第 m 时间段向后差分计算绝热温升增量修正系数，得到: ∆θm = λBmf (m∆τ 2 −
2 1 1 = ( L + ( m - ) ∆τ 2 ) 2 ∑ 2 1 k ( L + i∆ τ 1 ) 2 i = ( m−1) k + 2 mk 1
(8)
2.2 修正龄期的确定 由于混凝土浇筑初期温度变化较大，时间步长对计算结果影响较大，因而初期的修正尤为重要。 随混凝土温度变化逐渐趋于平缓，时间步长对计算结果的影响逐渐减小，当满足一定的精度要求时，
２５ ２３ ２１ １９ １７ １５ １３ １１ ０ ２０ ４０ ６０ ８０ 时间（ｄ） １００ 修正后 修正前
∆t t
０
t i−1
t i− 0. 5∆t
ti
图 2 误差形成示意图
图 3 水化热温升速率修正计算结果
2 温度场计算误差修正
为了减少计算量，同时保证求解精度，需要对由于时间步长过大造成的误差进行修正。在混凝土 施工期，影响混凝土温度变化的最主要因素是混凝土自身的水化热温升。因此，时间步长主要由水化 热温升速率来确定。 dθ 2.1 水化热升温速率 修正系数 λ m 的确定 dτ 设能满足精度要求的较小步长为 ∆τ 1 ，拟采用的较大步长为 ∆τ 2 ，在第 i 时间段 ∆τ i = τi - τ i -1 。当采用 向后差分中点龄期计算方法时， 由差分法特点可知， 在第 m 时间段有 ∆τ 2 = k∆τ1 (3) τ m = m∆τ 2 = mk∆τ 1 (4)在较小步长 ∆τ 1 下，可认为： ∆θ m = f (τ m − ∆τ1 / 2) ∆τ1 (5)随步长加大，(5)式左右两端差加大，等式条 件无法严格满足。在较大步长 ∆τ 2 下，如果存在系数 λBm 满足：
1 2
3 算例分析
混 凝 土 导 温 系 数 a = 0.01m 2 / d , 导 热 系 数 λ = 226.1kJ / m ⋅ d⋅o C ， 绝 热 温 升 公 式 采 用 双 曲 线 形 式 θ(t ) = 27.45t (1.1597 + t ) 。采用后点差分中点龄期法，计算以 1 天为步长，以 0.1 天步长对计算过程进行 补偿。 温升速率修正后的计算结果如图 3 所示， 其中修正前后对计算结果精度影响的最大值达到 13％。
第4期 2003 年 12 月
广 东 交 通 职 业 技 术 学 院 学 报
Journal Of Guang Dong Communication Polytechnic
NO.4 December 2003
化速率及时间步长有关，变化速率越快，时间步长越大，计算结果误差越大。当时间步长足够小时， 误差将趋于 0，时间差分步长越大，温度变化速率“以直代曲”的累计误差就越大。如图 2 所示。
参考文献: [1] Seatta A,Scotta R,Vitaliani R.Stress analysis of concrete structures subjected to variable thermal loads[J]. Journal of structures engineering,ASCE,1995,(2) [2] S.B.Tatro,E.K.Schrader.Thermal consideration for roller compacted concrete[J].ACI journal,1985，March. [3] Verback G J,Foster C W.Long-time study of cement performance in concrete[J]. Proc.ASTM,50,1950 [4] 王铁梦.大体积混凝土的瞬态温度场和温度收缩应力的计算机仿真[J].工业建筑，1990(1) [5] 钱壬章．俞昌铭.林文贵.传热分析与计算[M].高等教育出版社,1987 [6] J.A. 亚当斯,D.F. 罗杰斯.传热学计算机分析[M]. 科学出版社,1980
2
第4期 2003 年 12 月
广 东 交 通 职 业 技 术 学 院 学 报
Journal Of Guang Dong Communication Polytechnic
NO.4 December 2003
则无需再进行补偿。 设 ε 为允许误差， f (τ ) 为绝热温升变化速率。当Biblioteka Baidu用后点差分中点龄期方法计算，若在第 m 时间 段满足： f (τ m − ∆τ ) < ε (9) 即无需对水化热温升速率进行修正。
s ∆ ôn
n
s
s ∆ ôn

n -1
s
n -1
n
式中：{Tn}和{Tn-1}分别为时间τn 和τn-1 时结点温度向量； ∆ô n = τ n − τ n −1 是时间步长； H 、 R 为与热 学参数及单元形状有关的矩阵； F 是与水化热速率、降温速率及边界条件有关的项。 1.2 差分方法与求解精度 方程（2）中，s 为 0、1/2 和 1 时，分别对应了三种不同的差分方法：向前差分法、中点差分法和 向后差分法。三种情况下， ∆ô n 时间段绝热温升增量分别为 ∆θ = f (τ n -1 ) ∆τ n ，
mk -
则称 λBm
∑ f (i ∆τ 1 )
λBm =
i = ( m− 1) k + 1 2
1 2
1 kf (( m - )∆τ 2 ) 2
(7) 即以步长 ∆τ 2 求解时，只需将第 m 时间步对应的绝热温升上升速率
dθ dτ
用 λBm 进行修正，即能保证以 ∆τ 1 ， ∆τ 2 为步长求解绝热温升增量相同。 入式（7） ，得： λBm
1 温度场理论模型
1.1 三维瞬态温度场的有限差分模型 三维瞬态温度场的控制方程及边界条件为[5]：
TΓ1 = T0 ( x, y, z ), TaΓ =−
2
∂T ∂2T ∂ 2T ∂2T ∂θ = α ∂X2 + ∂Y2 + ∂Z2 + ∂τ ∂τ
=−
Γ3
(1)
1 ∂T q ( x, y, z ), λ ∂n
０
１０ ２０ ３０ ４０ ５０ ６０ ７０ ８０ ９０ １００
０
１０ ２０ ３０ ４０ ５０ ６０ ７０ ８０ ９０ １００
a
双曲线形式绝热温升计算 n=1
b 双曲线形式绝热温升计算 n=2
图 1 两种差分方法计算精度对比
２７ θ 绝热温升（℃） erf ∆θ ≈ ∂θ ∂t ∆t
t− 0. 5∆t