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

二维铸造充型过程数值模拟的特征分数步长法∗鲁统超1，葛亮21山东大学数学与系统科学学院, （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 lf V f V tρρρ−∂0+⋅+⋅=∂r r "" (2.4)其中sf ,lf 为固体和液体对应的分数，且ss l l ff ρρρ+=,1s l f f += 这里将液态金属看作不可压缩流体，并假设铸造充型过程进行很快，液态金属来不及凝固，从而方程简化为u vV 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 Nt ∆=节点,其中01k N =,,⋅⋅⋅⋅⋅,。

采用交错网格来离散整个计算域，压力布置在单元的中心，速度变量布置在单元界面上，如（图1 ）所示- 2 -（图1） 111122221;2;3;4;5i j i j i j i j i j u u v v p ,+,−,,−,+首先令，()(nnu x y u x y t ,=,,))n()(np x y p x y t ,=,,，利用分数步长法处理动量方程：12111112n n n nn n n n x uu u u p u v g u t x y x µn ερρ+,−∂∂∂++=−++∆+∂∂∂ (3.1)121121112n n n u up txερ++,′−∂=−+∆∂ (3.2)12211112n n n nn n n n y vv v v p u v g v t x y y µn ερρ+,−∂∂∂++=−++∆+∆∂∂∂ (3.3)121221112n n nvv p tyερ++,′−∂=−+∆∂ (3.4) 式中1n n pp p +′=+，为截断误差。

其中(3.1),(3.3)为第一步，由时刻的已知状态量，求中间过渡量()n i j O t ε,=∆n t 12n u +,12n v+,第二步由这一过渡量求 1n t+时刻的值，同时计入压力项的变化对方程的影响。

先处理下半个时间步长式(3.2)对x 求偏微分，(3.4)对求偏微分，并将两式叠加得：y 1122112222111()()()n n n n u v u vp p O t t x y t x y x yρ++++′′∂∂∂∂∂∂+−+=−++∆∆∂∂∆∂∂∂∂() 为了保持不可压流特征，对时刻强迫其满足不可压条件，即1n t t +=110n n u v x y++∂∂+=∂∂ 从而可得到方程：- 3 -11()p DU O t tρ′−=−+∆∆ (3.5) 其中1122n n u v DU x y++∂∂=+∂∂ 下面我们给出方程(3.5)在边界上的处理将方程(3.2)，(3.4)投影到边界Γ的外法线而得到：11()()n p n V V n t+∗ΓΓΓ′∂=−−∂∆r r ⋅ 其中V 为 ∗Γr 12n V +r 在边界上的值，并且由[ 6]可知数值解与V ∗Γr 的取值无关，故取从而得到： 1n V V ∗+Γ=r r()p nΓ0′∂=∂下面我们分析下半个时间步的计算与收敛性 由方程(3.2),(3.4),(3.5)121 111()2n n u up O t t x ρ++′−∂=−+∆∆∂121 111()2n n v vp O t tyρ++′−∂=−+∆∂∆1122111()n n u v ()p DU O t t t x yρ++∂∂′−∆=−=−++∆∆∆∂∂ 写为向量形式为： 1211()n n V V p O t tρ++−′=−+∆∆r r "1211()n p div O t V tρ+′−∆=−+∆∆r利用迭代法求解方程，令12n w V+=r ,0q p ′=,t τ=∆下面我们用[ 9 ]中的方法得到下面方程：111()m mm m q q B q divw ττ+++−0−∆+= (3.6)- 4 -11()m mm w w q βττ++−0+=" (3.7)其中1β=，B 为一个预处理算子，我们取其为22(1)B ττ=−−∆+∆ 从而当(01)τ∈,充分小时，B 为正定的算子，从而有下式1212k k k B k ∃,,−∆≤≤−∆把B 代入(3.6)有11()2m mm q q q divw ττ+m +−−∆=∆− (3.8)对于上式，由[ 8 ]当 时，(3.8)的解收敛到(3.6)的解。

m →∞下面我们来证明迭代算法(3.6),(3.7)的收敛性 令表示迭代误差，其中q,w 表示精确解，令mmm y q q r w w =−,=−m11ˆˆˆˆm m m m t t yy rr y y yy y r w rw r ττ++−−=,=,=,=,=,=ˆˆ0t By y divr τ−+= (3.9)ˆ0t r yβτ+=" (3.10)考虑(3.9)的弱形式112210ˆˆ()()()0t()B y B v y v r v v H D τ,+,−,=,∀∈""" (3.11) 取ˆ2v yτ=代入(3.11), 并对(3.10)两端用ˆ2r τ作内积，取两式之和整理得：22222222321ˆˆˆ20B B t B t y y y y r r r ττβτβτβτ||||−||||+||||+||||+||||−||||+||||=22()()B L y By y Ω||||=,另外从(3.10)中得到11ˆˆrr ββ=−"y只要为梯度函数，即有为梯度函数。

r ˆr从而由弱形式(3.11)有11221ˆˆ()()()ˆt t B B y B v r v yv y yv v v ττ|,||,||,|≤+≤||+||||||||||||""""""|||| 即1222212ˆˆ)t B t B r y y k y y τ||||≤||||+||||≤||||+||||1ˆτ (3.12)- 5 -将(3.12)两端平方并乘以2βτλ, 其中0λ>为待定系数222222222221ˆˆˆ(12)2(1)B t B B y k y y r r y τβλττβλβτλβτβτ||||+−||||+−||||+||||+||||≤||||+||||2ˆr 选λ充分小， 使2212010k βλτβλ−>,−>可导出下式22222ˆˆ(1)(1)B B Cy C r y r k τβττβτ′+||||++||||≤||||+||||2其中22(1)0C C λτβλ′=,=−>从而取2211max(1)(1)1k r C C ττ−−′=+,+<得到：2222ˆˆ()B B y r r y r βτβτ||||+||||≤||||+||||得到迭代法收敛。

不同的离散方式只要使B 满足上面条件，就能保证收敛性。

现给出下半步长的全离散格式： 由(3.6)令12w w w ⎛⎞⎜⎟⎜⎟⎜⎟⎜⎟⎝⎠=11221111m mi j i jm mi j i j q q t wwx ρ+,,+,+,,+,−∆=−∆ (3.13)11221122m mi j i jm mi j i j q q t w w y ρ,+,+,,+,,+−∆=−∆ (3.14)111122222111111122(4)2()4m m m m m m mi ji ji j i j i j i j i j m m m mi j i j i j i j t q q q q q q q w w w w t x yρ+,,,+,−,,+,−,+,,−,,,+,,−∆=−−−−−−−∆−+∆∆ (3.15) 下一步离散(3.1),(3.3)： 令()x y φ,==11{}()2V x y t φτ−∂∂=+⋅,∂∂r "从而有- 6 -1()n nu p x y u x µφτρρ∂∂,=−+∂∂1()n nv p x y v y µφτρρ∂∂,=−+∂∂ 用τ-特征方向的向后差商逼近前半个步长：1122()()n n u u x t x x y ττ++∂∂,=,=,∂∂r r1122()()n n v v x t x x y ττ++∂∂,=,=,∂∂r r11ˆ()()()()(n n nux ux u x x y x y tφφτ++∂−,≈,=∂∆r r r1122ˆ()()()()n n nvx vx v x x y tφτ++∂−,=∂∆r r r其中ˆxx V t =−∆r r r 对方程中其他项进行离散有111()n ni j i jn i j p p p x xρρ+,,,−∂−≈−∂∆31111122222112222()(nn n n n n i j i ji ji j i j i j ni j u u u u u u u xyµµρρ+,−,+,+,++,+,−,−+−+∆≈+∆∆2)从而有1211223111112222221211122ˆ122()n n n n i ji j i j i j n nnnnni j i j i j i j i j i j n i ju u p p tx uu uuuue x y ρµρ++,+,+,,+,−,+,+,++,+,−+,−−=−∆∆−+−+++∆∆+ (3.16)1211223111112222221211122ˆ122()n n n ni j i j i j i j n n n nn ni j i j i j i j i j i j n i j v v p p ty vvvvvve xyρµρ+,+,+,+,,+,−+,+,+−,+,+,+−−=−∆∆−+−+++∆∆+ (3.17)- 7 -式中1n i j u +,(1n i j v +,)为()在u v 1n t+时刻(i j )x y ,点的值，ˆn i j u ,(ˆn i j v ,)为()在时刻u v nt ˆi j x ,r的值,ni j p ,为p 在时刻(nt )i j x y ,点的值,12ni j e +,与12ni j e ,+为截断误差。