铸件凝固模拟过程中的几种时间步长计算方法探讨_龚文邦

(13)
从 x 、y 、z 方向流入 单元(i , j , k)的 热量总 和为
6
∑qi , 根据热量守恒原理 , 有 :
i =1
∑ ρC p Δx 0 Δy 0 Δz 0
T′0 - T0 Δτ
6
= qi i =1
式中 T 0 ———图 1 单元在当前时刻的温度 ;
(14)
T′0 ———图 1 单元在 t +Δτ时刻的温度 ;
件的时间步长 。
2. 3. 2 边界单元时间步长的计算
当网格单元处于边界时 , 将(12)式代入(14)式 :
∑ ρCp Δx 0 Δy 0 Δz0
T′0 - T0 Δτ
l
= (αc +αr)si(Ti i =1
- T0)
∑6
+ si(T i - T 0 )/ i =1
Δx 0 2λ0
+
Δx i 2λi
+
γ0 , i
(19)
将(19)式整理 , 可得 :
T′0
=
Δτ A
ΔAτ-
B
T 0 +Cຫໍສະໝຸດ Baidu
(20)
式中 B 、C 为正 数 , 根 据计算单 元所处的 位置不
同 , B 、C 有不同的表达式及不同的计算值 。 从(20)式
可知 , 要使 T′0 数值解稳定 , 必须满足 :
(4)
式中 T L ———液相线温度 ; T S ———固相线温度 。
将式(4)代入式(3), 整理得 :
ρcp
+ TL
L -
TS
ΔT
=
λ6(TΔix
-T Δx
0)Δτ
(5)
令 ce =cp +L /(T L - TS ), 为等效比热容(无相变时 , ce =c), 则有 :
ΔT
=ρ6cλe
Ti - T 0 Δx Δx
由式(17)可知 , 要使 T′0 有稳定的解 , 则 T0 的系数必
须 >0 , 即
ΔAτ-
6λ0 Si Δx 0 +λ0 γ0 , i
≥0
解不等式 , 有 :
Δτ≤ρcp Δx 0 Δy 0 Δz 0 (Δx 0 +λ0 γ0 , i)/6λ0 Si (18)
由上式计算出的 Δτ即为内部网格点满足收敛约束条
(15)
式中下标 0 表示图中当前单元(i , j , k), Δy 0 Δz0 =
S1 为传热方向上的单元面积 。
q1 为一垂直于 x 轴的面上传入单元(i , j , k)的热 量 , 同理 , 从垂直于 x 轴的另一面及垂直于 y 轴 、z 轴
的面传入的热量可用通式表示为 :
qi =si(T i - T0 )/
Δτ
(6)
为了使数值解稳定 , 必须满足 :
ΔT ≤|T i - T 0 |
(7)
由此可得 :
Δτ≤ р6cλe Δx Δx
(8)
(8)式即可计算出 F ourier 方程满足数值解稳定
的时间步长 。
2 考虑对流 、辐射换热时单元网格时间步长的计算
铸件 /铸型系统的传热通常包含高温金属及铸型
表面的辐射传热 、液态金属内部及其与铸型间的对流
基金资助项目(2002 年) 作者简介 :龚文邦(1965- ), 湖南岳 阳人 , 研究 员 , 在读博 士 , 主要从 事
材料成形过程控制与工艺研究. Email :g on gwenb @pub lic. w h. hb . cn
ρcp
T t
=
x
λT x
+y
λT y
+z
λT z
+Q′ (1)
式中 ρ———密度 ;
图 1 单元(i , j , k)与相邻 单元的热交换
Fig . 1 Heat ex change between units(i, j , k)and adjoining units
在 Δτ时间内 , 单元(i , j , k)的总热量变化 Q 为 :
Q
=ρC p Δx 0 Δy 0 Δz 0
T′0 - T0 Δτ
2. 3 时间步长的计算 2. 3. 1 内部单元时间步长的计算
当网格单元处于内部时 , 将(16)式代入(14)式 :
ρC p Δx 0 Δy 0 Δz 0
T′0 - T 0 Δτ
=
∑6
S i(T i - T0 )/ i =1
Δx 0 2 λ0
+
Δx i 2 λi
+γ0 ,
i
令 A =ρCp Δx 0 Δy 0 Δz 0 , 假定以均匀网格剖分 , 且
qi ———图 1 单元从第 i 个相邻单元传入的热量
(i =1 ～ 6);
Δx0 、Δy 0 、Δz0 ———图 1 单元的 网格尺寸 , 为均
匀网 格 时 , Δx0(i) = Δy0(i) =
Δz 0(i) 。
qi
=
Δx 0 2 λ0
Δy 0 Δz 0
+2Δλxi-
i1,
1 j,
k
+γ[ 0] [ i- 1 , j , k]
(11)
αr
=
εσ(T
4 i
-
T 40) /S i(Ti
-
T 0) 为辐射换热系数 。
在边界的单元 , 一般同时存在对流和辐射 , 由此传入单
元的热量为 :
qi =qic +qir =(αc +αr )Si(T i - T 0 ) (12)
2. 2 能量守恒形式下的差分方程
网格单元(i , j , k)与相邻单元间的换热见图 1 。
换热 、金属向铸型的热传导 3 种形式 。
2. 1 三种换热模型
对于铸件 内部网格 单元 , 传热模 型依据 F ourier
方程 , 如式(1)。 对于铸件的边界单元 , 单元某个面与空气或铸型
接触 , 则该单元面与空气或铸型之间产生辐射与对流 。
对流换热可用 New t on 冷却定律描述 :
为了在满足稳定求解的前提下 , 尽量加快计算速 度 , 本文根据不同情况经理论推导 , 得到了几种时间步 长的计算方法 。 1 满足 Fo urier 方程收敛条件的时间步长计算
金属液充型后 , 金属与铸型内部传热主要以不稳 定导热方式进行 。其控制方程为 :
收稿日期 :2005-05-28 ; 修订日期 :2006-01-17 基金项目 :国家自然科学基金资助项目(10176009);国家航天支撑技 术
cp ———比热容 , J / kg K ; T ———温度 , K ; t ———时间 , s ; λ———导热系数 , W /m K ;
Q′———内热源 , Q′=рL f s/ T ; L ———熔化潜热 , J / kg ; f s ———固相率 。
按差分格式将方 程(1)离散 , 并 假定 Δx =Δy =
q
=
εσ(T
4 f
-
T
4 w
)
对铸件边界单元有
qir =εσ(T4i - T4o )
(10)
式中 qir ———辐射换热传入的热量 ;
ε———辐射系数 ;
σ———S tefen-Bolt zman 常数 。
为了便于计算 , 将式(10)转换为(9)式的形式 :
qir = αr Si(Ti - T 0)
Δτ+ρQc′p Δτ
(3)
式中 Σ0 、T′0 ———分别为计算单元在当前时刻和下一
《铸造技术》04 /2006
龚文邦等 :铸件凝固模拟过程中的几种时间步长 计算方法探讨
39 9
时刻的温度 。
在凝固区间内 , 假定潜热在整个结晶温度区间均匀释
放 , 则有 :
Q′Δτ=-
ρL
ΔT TL -
Ts
摘要 :为 了尽可能缩短铸件凝固模拟的计算时间 , 从不 同的角 度对凝 固模拟 求解中 的时间步 长的计 算进行 了分析 ;并根 据能 量守恒原理 , 在考虑热传导 、对流 、辐射情况下进行了时间 步长的 推算 , 得出 了差分 网格单元 在不同 条件 、不同 位置的实 用时 间步长计算公式 。 通过实例验证 , 采用准确的时间步长 , 既能 避免计算发散又能提高计算效率 。 关键词 :凝固模拟 ;时间步长 ;计算方法 中图分类号 :TG 244 文献标识码 :A 文章编号:1000-8365(2006)04-0398-03
Δx 0 2 λ0
+
Δx i 2 λi
+γ0 ,
i
(16)
式中 Si ———qi 传入方向上的单元面积 ; T i ———图 1 中网格单元的第 i 个相邻单元当前 时刻的温度 ;
V ol. 27 N o. 4
40 0
F O UN D RY T ECH NO LO G Y
A pr. 2006
γ0 , i ———图 1 单元与第 i 个相邻单元间的热阻 。
q =α(T f - Tw )
对铸件边界单元网格 , 有 :
qic = αc Si(T i - T o)
(9)
式中 qic ———对流传入单元的热量 ; αc ———对流换热系数 ;
Si ———换热面积 ;
T i ———环境热力学温度 ;
T o ———计算单元在当前时刻的热力学温度 。 而辐射换热遵循 S tefen-Bolt zm an 定律 :
Discussion on Algorithm for Time Step in Casting Solidification Simulation
GONG Wen-bang, CHEN Li-l iang, LIU Rui-xiang , ZHOU Jian-xin , LIAO Dun-ming (School of Materials Science and Engineering, Huazhong University of Science and Technology ,Wuhan 430074 , China)
铸件凝固过程是一个非常复杂的物理化学变化过 程 , 由热量传输 、动量传输 、质量传输及相变等一系列 过程耦合而成 , 精确描述凝固过程十分困难 。 因此 , 在 凝固模拟中通常基于“瞬时充型 、初温均布” 的假设而 只求解温度场 。 温度场数值求解中广泛采用有限差分 显式算法 。此方法求解计算过程简便 , 但其稳定性和 收敛性强烈地受时间步长和空间步长的影响 。在适当 的条件下 , 从稳定性可以推出收敛性 , 即稳定是收敛的 充分条件 。因此 , 求解温度场方程之前确定适当的时 间步长 Δτ是稳定 、快速完成求解过程的前提 。一般可 通过试算或根据经验选定一个较小的 Δτ值 。
Δz , 有下式 :
6
∑ T′0 - T 0
Δτ
=
λ ρcp
Ti
i =1
Δx Δx
+ρQc′p
(2)
6
∑ Ti =Tx+Δx +Tx-Δx +Ty+Δy +Ty-Δy +Tz+Δz +Tz-Δz
i=1
a =λ/ρcp 可将(2)式转换为 :
6
∑ T′0 = T 0
+a
i
=1
Ti Δx
- 6T Δx
0
各网格材质特性一致时 , 上述方程可简化为 :
A(T′0 - T 0 ) Δτ
=
6λ0 Si (T i - T0 Δx 0 +λ0 γ0 , i
)
,
经整理 , 可得 :
T′0
=
Δτ A
ΔAτ-
6λ0 Si Δx 0 +λ0 γ0 , i
T0 +ΔAτ
6λ0 Si T i Δx 0 +λγ0 , i
(17)
Abstract :In order to cut calculation time of solidification simulation to as short as possible , the calculation of time step in sol idification simulation is analyzed from different aspect. Time step is calculated by considering the heat conduction, convection and radiation based on the conservation of energy. Practical calculation formula of time step for different grid unit under different condition and place is obtained. The experiments show that :They are onot onty avoid operating radiation but improve fate by correct time step. Key words :Sol idification simulation ;Time step ;Algorithm
铸造技术
39 8
F O UN D RY T ECH NO LO G Y
计算机应用 Computer Applicat ion
V ol. 27 N o. 4 A pr. 2006
铸件凝固模拟过程中的几种时间步长计算方法探讨
龚文邦 , 陈立亮 , 刘瑞祥 , 周建新 , 廖敦明
(华中科技大学材料科学与工程学院 , 湖北 武汉 430074)