二维铸造充型过程数值模拟的特征分数步长法

二维铸造充型过程数值模拟的特征分数步长法∗

鲁统超1，葛亮2

1山东大学数学与系统科学学院, （250100） 2

山东大学数学与系统科学学院, （250100）

E-mail ：lutc@

摘 要：铸造充型过程的数学模型是包括连续性方程和动量方程的偏微分方程组。本文利用分数步长法将动量方程分裂成两部分，对第一个方程采用特征差分法进行处理，对第二个方程结合连续性方程进行处理后得到压力的 泊松方程，用迭代法进行求解，给出了收敛性分析和稳定性条件。

关键词：分数步长；特征差分；收敛性；迭代。

1. 引 言

铸造生产的实质就是直接将液态金属浇入铸型并在铸型中凝固和冷却，进而得到铸件。液态金属的充型过程是铸件形成的第一个阶段。许多铸造缺陷(如卷气、夹渣、浇不足、冷隔及砂眼等)都是在充型不利的情况下产生的。因此，了解并控制充型过程是获得优质铸件的重要条件。但是，由于充型过程非常复杂，长期以来人们对充型过程的把握和控制主要是建立在大量实验基础上的经验准则。随着计算机的发展，铸件充型过程数值模拟才得到广泛应用。

充型过程流场数值模拟的主控方程均为非线性方程。其计算使用有限差分或有限元等数值方法求解质量守恒方程(连续性方程)和动量守恒 方程即Navier-Stokes 方程，以得出流体运动规律。在以前的研究中，Chorin(1968)和Temam(1969)分别独立的提出投影法。1972年由Minnesota 大学的Patankar 与Spalding 提出了simple 算法，这是一个压力修正算法，在以后的研究中又有simplec 方法,Raithby 提出的simplex 方法, Sheng 等提出的simplet 算法。

本文中利用分数步长法的思想将动量方程分裂成两部分，对第一个方程采用特征差分法求解，对第二个方程结合连续性方程进行处理后得到压力的 泊松方程，我们用迭代法进行求解，给出了收敛性分析和稳定性条件。

2. 问题的数学模型

铸造充型过程的模型主要由连续性方程和动量方程组成。 (a) 流体的动量方程

1x u p V u g u t x µρρ

∂∂=−⋅−++∆∂∂r " (2.1)

∗

本课题得到教育部高等学校博士点基金资助，编号：20030422049

- 1 -

1y v p V v g v t y µρρ

∂∂=−⋅−++∆∂∂r " (2.2)

其中, ，V u ()x y ,∈Ω(0]t J T ∈=,i vj =+r r r

是速度向量，p 为压力，, 为重力

分量，x g y g ρ为平均密度，µ为黏度。

(b) 流体区域的连续性方程

()0V t

ρρ∂+⋅=∂r " (2.3)

[()]()l s l l

f V f V t

ρρρ−∂0

+⋅+⋅=∂r r "" (2.4)

其中s

f ,l

f 为固体和液体对应的分数，且s

s l l f

f ρρρ+=,1s l f f += 这里将液态金属

看作不可压缩流体，并假设铸造充型过程进行很快，液态金属来不及凝固，从而方程简化为

u v

V x y

∂∂⋅=+=∂∂r " (2.5)

(c)速度边界条件

速度边界条件包括型壁边界和流体的自由表面边界。对于型壁边界,与型壁垂直方向应

满足速度为0。对于自由表面边界要求其表面速度满足连续性条件，即

0u v V x y

∂∂⋅=+=∂∂r "

3. 方程求解

对进行剖分,设x,y 方向的步长分别为Ωx ∆,y ∆, 节点i x i x =∆,。时间步长

j y j =∆y T N

t ∆=

节点,其中01k N =,,⋅⋅⋅⋅⋅,。采用交错网格来离散整个计算域，压力布置在单元的中

心，速度变量布置在单元界面上，如（图1 ）所示

- 2 -

（图1） 11112

2

2

2

1;2;3;4;5i j i j i j i j i j u u v v p ,+,−,,−,+

首先令，()(n

n

u x y u x y t ,=,,))n

()(n

p x y p x y t ,=,,，利用分数步长法处理动量方程：

1

2

111112

n n n n

n n n n x u

u u u p u v g u t x y x µn ερρ+,−∂∂∂++=−++∆+∂∂∂ (3.1)

1

2

1121

112n n n u u

p t

x

ερ++,′

−∂=−

+∆∂ (3.2)

1

2

211112

n n n n

n n n n y v

v v v p u v g v t x y y µn ερρ+,−∂∂∂++=−++∆+∆∂∂∂ (3.3)

1

2

1

22

1112

n n n

v

v p t

y

ερ++,′−∂=−

+∆∂ (3.4) 式中1

n n p

p p +′=+，为截断误差。其中(3.1),(3.3)为第一步，由时刻的已知

状态量，求中间过渡量()n i j O t ε,=∆n t 1

2

n u +,1

2

n v

+,第二步由这一过渡量求 1

n t

+时刻的值，同时计入压力项的

变化对方程的影响。

先处理下半个时间步长

式(3.2)对x 求偏微分，(3.4)对求偏微分，并将两式叠加得：

y 1

1

22

11

2222111()()()n n n n u v u v

p p O t t x y t x y x y

ρ++++′′∂∂∂∂∂∂+−+

=−++∆∆∂∂∆∂∂∂∂() 为了保持不可压流特征，对时刻强迫其满足不可压条件，即

1

n t t +=

11

0n n u v x y

++∂∂+=∂∂ 从而可得到方程：

- 3 -