DOI:10.19476/j.ysxb.1004.0609.2000.02.018
铸件三维充型过程耦合数值模拟
周彼德 薛祥 糜忠兰 孙小波 张春晖 马建
哈尔滨工业大学材料科学与工程学院!哈尔滨150001
哈尔滨工业大学材料科学与?
摘 要:
采用改进后的SOLA VOF方法处理自由表面 , 有效地克服了传统的SOLA VOF法中的“施体 受体”法在处理三维自由表面时精度差和时间步长小的缺点。在自由表面处理上采用了新的压力插值函数 , 使三维自由表面的处理变得简单明了。用有限差分法建立了铸件三维充型过程流动与传热的耦合计算数值模型 , 并对大型薄壁铝合金铸件进行了数值计算。后处理中采用插值技术 , 直观地显示了充填过程中自由表面的位置和形状
关键词:
大型薄壁件 ;SOLA-VOF ;流体流动 ;传热 ;
中图分类号: TG248
收稿日期: 1998-12-07
3D coupling numerical simulation of mold filling
Abstract:
The numerical model for 3D fluid flow coupled with heat transfer was established with FDM and the free surfaces were handled by an improved SOLA VOF technique. A pressure interpolation function was adopted as the free surface boundary condition. The example of large thin wall Al alloy casting was calculated with this model. In the post processing, the location and shape of free surfaces during mold filling were displayed visually with a linear interpolation method.
Keyword:
large thin-wall castings; SOLA-VOF; fluid flow; heat transfer;
Received: 1998-12-07
随着计算机及数值计算技术的完善, 利用计算机对充型过程进行数值模拟已成为可能。 充填过程数值模拟即通过建立并求解描述这一过程的微分方程, 得到充型过程中压力场、 速度场、 温度场以及自由表面的定量变化, 从而为选择正确的浇注方法、 控制成形中缺陷的产生奠定科学基础。
虽然充型过程数值模拟日益引起人们的重视, 但由于涉及的控制方程多、 计算量大且迭代结果容易发散, 特别是对自由表面问题的处理难度大, 模拟结果难以进行实验验证, 使充型过程的数值模拟工作受到阻碍。 此外, 由于充型过程伴随着传热现象, 使液体温度降低, 尤其是自由表面前沿, 当液体温度低于液相线温度时, 固相开始析出, 粘度急剧增加, 流动阻力不断增大。 当流股头部固相量达某一临界值时, 该部分停止流动, 从而形成冷隔或浇不足, 导致铸件不能完全充填
[1 ]
。 采用传统的实验方法预测这些缺陷的形成既不科学也不经济, 而用数值计算方法模拟充填过程, 不但可以预测上述缺陷的位置, 为优化铸造工艺奠定基础, 还可为后续温度场的计算提供初始温度分布
[2 ,3 ,4 ]
。 为此, 本论文中采用了改进后的SOLA-VOF法处理自由表面, 建立了流场与温度场耦合的数值模型, 数值计算采用有限差分法对模型进行离散化。
1 铸件三维充型过程的数学模型
液态金属的充填过程可认为是不可压缩的牛顿型流体的非定常流动过程, 流体流动遵循质量守恒、 动量守恒和能量守恒。 据此, 本文模拟采用的数学模型如下:
1) 连续性方程
[6 ]
?
u
?
x
+
?
v
?
y
+
?
w
?
z
=
0
?
?
?
(
1
)
2) 动量守恒方程
[6 ]
?
u
?
t
+
u
?
u
?
x
+
v
?
u
?
y
+
w
?
u
?
z
=
?
?
?
?
γ
(
?
2
u
?
x
2
+
?
2
u
?
y
2
+
?
2
u
?
z
2
)
-
?
p
ρ
?
x
+
g
x
?
?
?
(
2
)
?
v
?
t
+
u
?
v
?
x
+
v
?
v
?
y
+
w
?
v
?
z
=
?
?
γ
(
?
2
v
?
x
2
+
?
2
v
?
y
2
+
?
2
v
?
z
2
)
-
?
p
ρ
?
y
+
g
y
?
?
?
(
3
)
?
w
?
t
+
u
?
w
?
x
+
v
?
w
?
y
+
w
?
w
?
z
=
?
?
γ
(
?
2
w
?
x
2
+
?
2
w
?
y
2
+
?
2
w
?
z
2
)
-
?
p
ρ
?
z
+
g
z
?
?
?
(
4
)
式中 u , v , w —x , y , z 方向的速度分量, m/s; p —压力, Pa; t —时间, s; g x , g y , g z —x , y , z 方向的重力加速度, m/s2 ; ρ —流体密度, kg/m3 ; γ —流体运动粘度, m2 /s。
3) 能量方程
[6 ]
?
Τ
?
t
+
u
?
Τ
?
x
+
v
?
Τ
?
y
+
w
?
Τ
?
z
=
?
?
λ
ρ
C
p
(
?
2
Τ
?
x
2
+
?
2
Τ
?
y
2
+
?
2
Τ
?
z
2
)
?
?
?
(
5
)
式中 T —温度, K; λ —流体导热系数, W/ (m·K) ; C p —流体比热, J/ (g·K) 。
4) 体积函数方程
[7 ]
?
F
?
t
+
u
?
F
?
x
+
v
?
F
?
y
+
w
?
F
?
z
=
0
?
?
?
(
6
)
式中 F —流体体积分数。
为减小几何误差, 提高计算精度, 研究中采用非均匀网格剖分技术, 将研究对象划分为六面体单元。 求解动量守恒方程时采用交错网格技术, 以避免出现棋盘形压力场导致的不正确解。 用VOF解体积函数方程时采用了改进后的算法, 对方程 (6) 进行离散处理
[5 ]
:
F
Ρ
*
-
F
Ρ
n
Δ
t
+
F
e
*
u
e
-
F
w
*
u
w
Δ
x
+
?
?
F
n
*
v
n
-
F
s
*
v
s
Δ
y
+
F
t
*
w
t
-
F
b
*
w
b
Δ
z
=
0
?
?
?
(
7
)
式中 下标P 表示变量F 的位置在所计算单元的中心; 下标e , w , n , s , t , b 分别表示单元的六个面; 上标n 表示时间水平, F * 表示对n +1时刻的F 值的猜测值。 为了避免出现负的F 值或F 值比1大时造成的体积误差, 必须对自由表面单元六个面上的F * 值进行修正, 以便获得F n +1P 值。 这需要一个迭代过程。
在边界条件处理时, 由于三维情况下自由表面往往相互交叉或相互重叠, 导致同一自由表面单元可能存在几个插值单元, 为此, 研究中采用的自由表面边界条件函数为
[5 ,8 ]
:
p
Ρ
=
(
F
Ρ
-
0
.
5
)
1
Ν
∑
n
=
1
Ν
p
n
F
n
+
(
1
.
5
-
F
Ρ
)
p
s
?
?
?
(
8
)
式中 F P —表面单元P 的流体体积分数; F n —相邻满单元的流体体积分数; N —相邻满单元的个数; p n —相邻满单元的压力, Pa; p s —环境压力, Pa; p P —表面单元压力, Pa。
计算过程中, 为避免发散, 时间步长必须满足以下三个不等式:
Δ
t
1
<
min
{
Δ
x
min
|
u
e
|
?
Δ
y
min
|
v
n
|
,
Δ
z
min
|
w
t
|
}
?
?
?
(
9
)
a e +a w +a n +a s +a t +a b <1 (10)
Δt =min (Δt 1 , Δt 2 ) (11)
式 (9) 中的min{…}表示取最小值, Δx min 表示空间步长的最小值,
|
u
e
|
表示x方向的速度分量的绝对值, 其余类似; 式 (10) 中的ae 表示动量离散方程的系数, 其余类似; 式 (11) 中的Δ t2 是由式 (10) 计算出的时间步长。
在以上理论基础上, 开发的铸件三维充型过程耦合数值模拟程序的计算步骤如下:
1) 初始化速度场、 压力场和温度场;
2) 利用初始条件和边界条件或前一时刻的计算结果计算最小时间步长;
3) 求解自由表面的体积函数方程;
4) 用上一时刻的压力场及初始条件、 边界条件解显式动量方程, 得新一时刻的试算速度场;
5) 通过迭代计算修正压力场和相应的试算速度场, 直到满足连续性方程为止;
6) 解隐式能量方程求得温度场;
7) 把修正后的速度场作为当前速度场, 修正后的压力场作为上一时刻的压力场, 回到第1) 步, 直到充填时间达到预定时间。
2 模拟程序的验证
为了验证所开发的模拟程序的可靠性, 本文针对大型薄壁铝合金铸件 (1 000 mm ×480 mm ×7 mm ) 进行了计算。 该铸件如图1所示。 铸型向下倾斜15°浇注。 计算时所采用的物性参数如表1所示, 运动学粘度为0.013 cm 2 /s 。
图1 铝合金平板铸件的FDM网格模型
Fig.1 FDM grid model of aluminum alloy plate casting
充填过程中流体形态变化如图2所示。 所有图中色温对照表中的温标显示的不是温度, 而是将流体体积分数扩大1 000倍后的值, 表示单元的充满程度。 1 000所对应的白色表示满单元, 而0则表示空单元。 由图2可以看出, 由于模拟计算时假设浇注系统与型腔具有一定的斜度, 这样, 当液流到达型腔末端时, 液流失去速度头开始向上翻转形成回流, 最后两股液流在离型腔末端不远处汇合。 由此可见, 在这种充填条件下, 容易发生缺陷的部位在离型腔末端不远处。 另外, 由于倾斜浇注, 沿铸件断面及同一平面上具有较大的温度梯度。
为了进一步检验计算程序, 研究了不同形状的铸件的充填情况, 分别进行了顶注和底注式充填模拟计算, 以便观察计算结果的对称性和自由表面的模糊程度。 图中的色标同样是表示单元内的流体体积分数。
图3 (a) , (b) 分别为顶注式和底注式模拟验证结果。 顶注式的铸件大小为150 mm×150 mm×30 mm, 网格尺寸为10 mm×10 mm×10 mm, 入口在铸件的上方中心, 速度为1 m/s。 图3 (a) 中只显示了最底层的充填情况, 由图可见, 模拟结果的对称性很好。 底注式的大小为30 mm×150 mm×150
表1 模拟计算使用的热物性参数
Table 1 Thermal physical parameters used in simulation
Material
Conductivity / (W·m-1 ·K-1 )
Specific heat capacity / (J·g-1 ·K-1 )
Density / (g·cm-3 )
Latent heat / (J·g-1 )
Liquidus point /℃
Solidus point /℃
ZL101
114
0.98
2.28
554.94
618
571.3
Green sand
1.28
2.303
1.8
-
-
-
图2 薄壁平板件充填流场模拟结果 (t=1.09 s)
Fig.2 Simulated result of mold filling for thin-wall plate (t =1.09 s)
图3 顶注式充填 (a) 和底注式充填 (b) 流场模拟结果
Fig.3 Simulated results of mold filling for top-pouring (a) and bottom-pouring (b) (a) —t =0.38 s; (b) —t =0.18 s
mm, 网格尺寸为10 mm×10 mm×10 mm, 入口速度为4 m/s。 从图3 (b) 可以看出, 由于入口速度太大, 液流到达横浇道末端后回流, 沿内浇道向上直冲到型腔顶部, 然后形成回流。 可见自由表面清晰, 并未产生界面模糊现象。 由此可以看出模拟计算是成功的。
另外, 为了检验所开发的三维充型过程耦合模拟计算程序的准确性, 研究中对铝合金薄壁平板铸件进行了实际浇注, 并用热电偶测量铸件8个点的温度—时间曲线, 然后与模拟计算的结果进行了比较。 二者基本一致。
3 结论
1) 通过对不同形状铸件的模拟计算, 验证了改进后的计算方法具有很好的对称性和清晰的自由表面。 说明了模拟程序具有良好的模拟自由表面的能力和可靠性。
2) 由模拟结果可以看出, 用所开发的三维耦合模拟计算程序可以很好地模拟充型过程, 从而为预测缺陷的位置, 优化铸造工艺设计奠定基础。
参考文献
[1 ] CHANGQing ming (常庆明 ) .大型薄壁铸件三维充型及传热的数值模拟 [D] .Harbin :HarbinInstituteofTechnology, 1 994:1 2 3~ 1 2 4 .
[] [2 ] LIChen xi, (李晨曦 ) .铸造充型过程流场计算机模拟技术述评 [J] .Foundry (铸造 ) , 1 997 (2 ) :47~ 50 .
[3] LeeWB , LuHYandLuiYB .JournalofMaterialsPro cessingTechnology .1 995 , 52 :2 4 8~ 2 69.
[4] UsmaniAS , CrossJTandLewisRW .InternationalJournalforNumericalMethodsinEngineering , 1 992 , 32 :787~ 80 6 .
[5] XUEX .ModelingandSimulationofFluidFlowandHeatTransferinMoldFilling [D] .Denmark :TechnicalUni versityofDenmark , 1 992 :1 4 2~ 1 4 6 .
[6] PatankarSV .NumericalHeatTransferandFluidFlow[M] .McGraw Hill, USA :HemispherePublishingCorp .1 980 .
[7] HirtCWandNicholsBD .Volumeoffluid (VOF) methodforthedynamicsoffreeboundaries [J] .JournalofComputationalPhysics, 1 981 , 39.
[8] XUEXandThorpeWB .Predictionofair entrapmentduringfillingathin sectionindirectsqueezecasting [J] .AFSTransactions, 1 995 , 1 0 3 :743