CN1209704C - 一种铸件充型过程模拟方法 - Google Patents

Abstract

本发明涉及一种铸件充型过程模拟方法。总流程为：首先选择任务类型，如是第一次开始，输入部分文件、控制及物理参数；否则输入上次计算结果，再计算；计算后回主菜单；选前处理任务类型时，判断何种类型几何文件，分别对铸件的构造三维实体几何模型的文件和有限元网格文件读取数据，部分后回主菜单；选后处理任务类型时，分别显示温度场、速度场、压力。它能合理设计充型速度、预测卷气缺陷、可为凝固过程模拟提供准确的初始温度条件。

Description

一种铸件充型过程模拟方法

技术领域

本发明涉及金属材料充型技术，具体地说是一种铸件充型过程模拟方法。

背景技术

充型过程是指液态金属通过浇注系统进入到型腔充满这一短暂的流动过程，它是注塑、铸造等成形过程中的重要环节。如果浇注系统设计不得当或充型速度过大或过小，会引起浇不足、卷气、冷隔等缺陷，严重影响铸件质量。然而，铸件的充型过程是液态金属在不透明砂型中的运动，液态金属流动速度快、温度高，难于直接观察研究，随着计第机及数值技术的发展，可利用计算机对充型过程进行数值模拟，动态再现浇注过程，以指导工艺实践。

发明内容

本发明的目的是提供一种能合理设计充型速度、预测卷气冷隔缺陷、可为凝固过程模拟提供准确的初始温度条件的铸件充型过程模拟方法。

为了实现上述目的，本发明的技术方案是：它是通过建立并求解描述液态金属流动过程的控制方程，得到充型过程中压力场、速度场、温度场以及自由表面的定量变化，从而选择正确的浇注方法、控制成形中缺陷的产生；硬件由计算、前处理、后处理三个任务组成，所述计算任务是完成对铸件充型过程中带有自由表面的Navier-Stokes方程及散热方程的求解，对Navier-Stokes方程的求解借鉴了Los Alamos实验室的报告中二维SOLA-VOF算法，将其推广到了三维，求出流体场之后，再采用帕坦卡教授的求解对流与扩散方程的解法求解散热方程；所述前处理任务是准确地构造铸件的三维几何模型，并进行网格剖分；所述后处理任务的功能是把由数值计算或试验获得的大量数据转换为人的视觉可以感受的图像，把大量的抽象数据有机地组合在一起，展示数据所表现的内容及之间的相互关系，展示场的分布；总流程为：首先选择任务类型，选计算任务类型时，先初始化自主开发的温度场计算类(Ctemperature)，如果是第一次开始，则输入剖分文件、铸造充型过程的控制及物理参数，进行计算后回主菜单；否则输入上次计算结果，再计算，计算后回主菜单；选前处理任务类型时，判断何种类型几何文件，经判断，分别对铸件的构造三维实体几何模型(Constructive Solid Geometry)的文件或有限元网格文件读取数据，进行剖分，保存剖分结果，然后回主菜单；选后处理任务类型时，则分别显示温度场、速度场、压力，回主菜单；

所述计算任务中计算部分是用SOLA求解铸件充型过程压力场和速度场，用VOF确定流动域和自由表面；采用施体-受体(Donor-Acceptor)方法计算铸造充型过程自由表面的扩展；具体流程为：先输入控制、物理性质参数，然后依据D-A算法计算自由表面边界条件，再计算自由表面方向，用t＝t+Δt累加时间值，更新速度场，构造铸件的三维几何模型已充满或者完成，当前计算则结束，否则以上一时刻的速度场、压力场作为估计值，显示计算本时刻的速度值，求出单元的散度D_i，j，k，并判断是否满足连续方程，如果所述判断为肯定结果，则直接求解散热方程，再回到依据Donor-Acceptor(D-A)算法计算自由表面扩展步骤；如果所述判断为否定结果，则修正压力，修正速度后，再重新判断，直到判断结果满足连续方程，至下一步求解散热方程；

所述用SOLA求解铸件充型过程压力场和速度场、用VOF确定流动域和自由表面，其计算步骤为：

(1)由Navier-Stokes方程的离散公式，以初始条件或前一时刻值为基础，计算当前时刻的试算速度；

(2)给每一个单元定义散度D_i，j，k：将所述步骤(1)的试算速度值代入连续性方程离散公式，求出D_i，j，k，若D_i，j，k＝(一般当D_i，j，k小于10^-3，即认为D_i，j，k＝0)，则步骤(1)中所述试算速度值满足连续性方程，即此时的速度场与压力场值既满足动量守恒方程又满足质量守恒方程，当前时间步长计算结束；如整个流场中有任一个单元不能满足连续性方程，即D_i，j，k不等于0，则需要下一步的修正；修正步骤：欲修正速度，须先修正压力；

(3)根据校正压力求出校正后的试算速度；

(4)将校正后的试算速度值代入所述步骤(2)，反复迭代直至所有单元均满足连续性方程；

(5)由体积函数方程确定新的流动域，对表面单元做合理设置；

(6)返回第一步，进入下一时刻计算；

所述计算自由表面扩展设计方法为根据计算单元的上游、下游单元，粗略估计出自由表面的形状，然后，根据不同形状，计算体积函数F在单元的边界面上的通量，从而计算出各个时刻体积函数F值；具体步骤是：

1)对边界上的任意两个单元，区分出施体(Donor)单元和受体(Acceptor)单元，依据所述两个单元边界上的速度的代数符号来确定，Donor单元是上游单元，其体积函数用F_D表示，Acceptor是下游单元，其体积函数用F_A表示；

2)确定出Donor单元和Acceptor单元边界上的体积函数值F_AD；

3)计算边界上的体积流量；

所述自由表面边界条件的设置在SOLA-VOF法中，边界包括自由表面、型壁(即网格边界)、入流边界，结合惯性原理和连续性原理处理自由表面速度，具体是：

1)自由表面分界条件：在所述计算单元内，流体保持原有的方向和流速流动，如果遇到障碍，将出流量均匀地转移到其它可能的出流方向，接着调整出流速度，使新生成的满网格单元能够满足连续性方程；所述对于自由表面的压力条件，满足正应力和切应力为零的条件；

2)型壁边界条件：使用“虚拟单元”的方法，在型壁处，设置虚拟单元，通过设定虚拟单元速度，设定自由滑、无滑等边界条件；

3)入流边界是指充型时浇注速度；

所述前处理任务包括与商品化软件的接口程序及构造实体几何模型造型与剖分两个部分，与商品化软件ProCAST的文件接口程序用于读取有限元网格文件，其具体流程为：由剖分步长计算出剖分步数N_x、N_y、N_z，再计算出单元的中心点坐标P(i，j，k)，判断该点属性，并在输出结果文件中写入对应的属性值；所述前处理任务中的构造实体几何模型造型与剖分是通过基本形体的集合运算来组成复杂形体，常用的集合运算有并、交、差和补运算，在计算机上对铸造充型的形体进行实体造型用数学方法描述出基本形体的点的集合及集合间的运算，据此确定出不同物理属性、不同材质的点的集合，如铸件、砂型和砂芯等；

所述后处理任务主要是以有限差分计算的可视化为处理形式，运用MicrosoftFoundation Class(MFC)提供的图像处理功能，动态显示计算得到的数据，具体表现为：1)具有二维或三维显示功能；2)对图形进行交互缩放的功能；3)动态地显示解随时间的变化情况；4)准确地显示足够多的细节，如温度梯度、等值线等，使人们能直观清晰地定性或定量地获取部分区域或某些点上的某些物理量上的值；5)实现交互实时的随视点位置变更而显示图形；6)交互地修改颜色对照表，并实时改变图形颜色。

本发明具有如下优点：

1.由于现有技术铸件生产中，充型过程是先于凝固过程的，采用本发明，进行模拟计算，能提供对流体场与温度场的耦合模拟，为凝固过程模拟提供了准确的初始温度条件。

2.采用本发明实现的流体场模拟可以为惯性对流和自然对流提供初始速度场，从而能预测卷气缺陷，优化浇注系统设计，以指导工艺实践。

3.本发明采用计算机图形技术把三维数据进行可视化处理，使人们能直接看到场的分布图像，这样，就可以大大加快和加深人们对铸造过程的认识，更为重要的是，它使处理海量数据成为可能，人们可以通过图像来分析数据，提取有用信息。

4.采用本发明得到的模拟结果，能优化浇注工艺，避免了浇不足、卷气等缺陷，优化工艺。

5.本发明采用的前处理任务，能方便准确地构造铸件的三维几何模型和进行网格自动剖分，是开发和推广铸造过程数值模拟软件的基础。

附图说明

图1为本发明程序的总体流程图。

图2为图1中计算部分程序框图。

图3为本发明与商品化软件的接口程序的自动网格剖分程序框图。

图4-a为本发明D-A技术说明中施体(Donor)、受体(Acceptor)单元的定义。

图4-b为本发明D-A技术说明中F_AD取D值的情况。

图4-c为本发明D-A技术说明中F_AD取A值的情况。

图4-d为本发明D-A技术说明中F_AD取A值的另一种情况。

图5为本发明自由表面边界条件说明示意图。

图6-a为本发明托盘件的几何形体正面示意图。

图6-b为本发明托盘件的几何形体反面示意图。

图7为ProCAST产生的托盘件的有限元网格。

图8为本发明剖分结果。

图9为本发明后处理任务输出结果。

图10为Al合金薄壁件的几何三视图。

图11为Al合金薄壁件采用本发明剖分后的结果显示。

图12为本发明在t＝0.763481s所得温度场和流体场结果。

图13为ProCAST在t＝0.89s时所得温度场和流体场结果。

图14为本发明在t＝0.132738s所得温度场和充型流体场结果。

图15为ProCAST在t＝0.1241s时所得温度场和充型流体场结果。

图16为本发明在t＝0.1685s所得温度场和充型流体场结果。

图17为ProCAST在t＝0.1474s时所得温度场和充型流体场结果。

具体实施方式

下面结合附图和实施例对本发明作进一步详细说明。

如图1所示，本软件由计算、前处理、后处理三个任务组成，所述计算任务是完成对铸件充型过程中带有自由表面的Navier-Stokes方程及散热方程的求解，对Navier-Stokes方程的求解借鉴了Los Alamos实验室的报告中二维SOLA-VOF算法，将其推广到了三维，求出流体场之后，再采用S.V.帕坦卡的求解对流与扩散方程的解法求解散热方程；所述前处理任务是准确地构造铸件的三维几何模型，并进行网格剖分；所述后处理任务的功能是把由数值计算或试验获得的大量数据转换为人的视觉可以感受的图像，把大量的抽象数据有机地组合在一起，展示数据所表现的内容及之间的相互关系，展示场的分布；总流程为：首先选择任务类型，选计算任务类型时，先初始化Ctemperature类，如果是第一次开始，则输入剖分文件、铸造充型过程的控制及物理参数，进行计算后回主菜单；否则输入上次计算结果，再计算，计算后回主菜单；选前处理任务类型时，判断何种类型几何文件，经判断，分别对铸件的Constructive SolidGeometry(CSG)造形文件或有限元网格文件读取数据，进行剖分，保存剖分结果，然后回主菜单；选后处理任务类型时，则分别显示温度场、速度场、压力，回主菜单。

如图2所示，所述计算任务中计算部分是用SOLA求解铸件充型过程压力场和速度场，用VOF确定流动域和自由表面；采用Donor-Acceptor(D-A)方法计算铸造充型过程自由表面的扩展；具体流程为：先输入控制、物理性质参数，然后依据D-A算法计算自由表面边界条件，再计算自由表面方向，用t＝t+Δt累加时间值，更新速度场，如果已充满或者完成，当前计算则结束，否则以上一时刻的速度场、压力场作为估计值，显式计算本时刻的速度值，求出单元的散度D_i，j，，k，并判断是否满足连续方程，如果所述判断为肯定结果，则直接求解散热方程，再回到依据D-A算法计算自由表面扩展步骤；如果所述判断为否定结果，则修正压力，修正速度后，再重新判断，直到判断结果满足连续方程，至下一步求解散热方程。

所述用SOLA求解铸件充型过程压力场和速度场、用VOF确定流动域和自由表面，整个计算过程中由速度初值及猜测压力值试算速度场的过程并不参与迭代，因而是一套迭代，其计算步骤为：

(1)由Naver-Stokes方程的离散公式，以初始条件或前一时刻值为基础，计算当前时刻的试算速度；

(2)给每一个单元定义散度D_i，j，k：将所述步骤(1)的试算速度值代入连续性方程离散公式，求出D_i，j，k，所述连续方程离散公式为：

$D_{i,j,k}=\frac{u_{i+\frac{1}{2},j,k}{-u}_{i-\frac{1}{2},j,k}}{\mathrm{\Δ}{x}_i}+\frac{v_{i,j+\frac{1}{2},k}{-v}_{i,j-\frac{1}{2},k}}{\mathrm{\Δ}{y}_j}+\frac{w_{i,j,k+\frac{1}{2}}{-w}_{i,j,k-\frac{1}{2}}}{\mathrm{\Δ}{z}_k};$

若D_i，j，k＝0(一般当D_i，j，k小于10^-3，即认为D_i，j，k＝0)，则步骤(1)中所述试算速度值满足连续性方程，即此时的速度场与压力场值既满足动量守恒方程又满足质量守恒方程，当前时间步长计算结束；如整个流场中有任一个单元不能满足连续性方程，即D_i，j，k不等于0，则需要下一步的修正；修正步骤：欲修正速度，须先修正压力，其修正压力用公式为：

Pⁿ⁺¹＝Pⁿ+δPⁿ：

其中，δPⁿ为压力修正量，其值可以用下式求得：

$\mathrm{\δP}=-D_{i,j,k}/\frac{\∂D_{i,j,k}}{\∂P}$

$\frac{{\∂D}_{i,j,k}}{\∂P}=\frac{\mathrm{\Δt}}{{\mathrm{\Δx}}_i}(\frac{1}{{\mathrm{\ρx}}_{i+\frac{1}{2}}}+\frac{1}{{\mathrm{\ρx}}_{i-\frac{1}{2}}})+\frac{\mathrm{\Δt}}{\mathrm{\Δ}{y}_j}(\frac{1}{\mathrm{\ρ}{y}_{j+\frac{1}{2}}}+\frac{1}{\mathrm{\ρ}{y}_{j-\frac{1}{2}}})+\frac{\mathrm{\Δt}}{{\mathrm{\Δz}}_k}(\frac{1}{{\mathrm{\ρz}}_{k+\frac{1}{2}}}+\frac{1}{\mathrm{\ρ}{z}_{k-\frac{1}{2}}})$

(3)根据校正压力求出校正后的试算速度，利用下式：

$u_{i+\frac{1}{2},j,k}^{n+1}=u_{i+\frac{1}{2},j,k}^n+\frac{\mathrm{\Δt\δ}{P}^n\mathrm{\ω}}{{\mathrm{\ρ\Δx}}_{i+\frac{1}{2}}}$

$u_{i-\frac{1}{2},j,k}^{n+1}=u_{i-\frac{1}{2},j,k}^n-\frac{\mathrm{\Δt\δ}{P}^n\mathrm{\ω}}{{\mathrm{\ρ\Δx}}_{i-\frac{1}{2}}}$

$v_{i,j+\frac{1}{2},k}^{n+1}=v_{i,j+\frac{1}{2},k}^n+\frac{\mathrm{\Δt\δ}{P}^n\mathrm{\ω}}{{\mathrm{\ρ\Δy}}_{j+\frac{1}{2}}}$

$v_{i,j-\frac{1}{2},k}^{n+1}=u_{i,j-\frac{1}{2},k}^n-\frac{\mathrm{\Δt\δ}{P}^n\mathrm{\ω}}{{\mathrm{\ρ\Δy}}_{j-\frac{1}{2}}}$

$w_{i,j,k+\frac{1}{2}}^{n+1}=u_{i,j,k+\frac{1}{2}}^n+\frac{\mathrm{\Δt\δ}{P}^n\mathrm{\ω}}{{\mathrm{\ρ\Δz}}_{k+\frac{1}{2}}}$

$u_{i,j,k-\frac{1}{2}}^{n+1}=u_{i,j,k-\frac{1}{2}}^n+\frac{\mathrm{\Δt\δ}{P}^n\mathrm{\ω}}{{\mathrm{\ρ\Δz}}_{k-\frac{1}{2}}};$

其中，ω为松弛因子(0＜ω＜2)，n与n+1表示校正循环次数；

(4)将校正后的试算速度值代入所述步骤(2)，反复迭代直至所有单元均满足连续性方程；

(5)由体积函数方程确定新的流动域，对表面单元做合理设置；

(6)返回第一步，进入下一时刻计算。

所述计算自由表面扩展设计方法是根据计算单元的上游、下游单元，粗略估计出自由表面的形状，然后，根据不同形状，计算体积函数F在单元的边界面上的通量，从而计算出各个时刻体积函数F值；体积函数控制方程的离散使用了D-A技术，对体积函数F进行求解时，保证了体积函数F的不连续性，且避免在一次计算中，流空单元或流量超过该单元所有的。

具体步骤是：

1)对边界上的任意两个单元，区分出Donor单元和Acceptor单元，依据所述两个单元边界上的速度的代数符号来确定，Donor单元是上游单元，其体积函数用F_D表示，Acceptor是下游单元，其体积函数用F_A表示；

2)确定出Donor单元和Acceptor单元边界上的体积函数值F_AD；如图4-a、图4-b、图4-c、图4-d所示，如果Donor单元是充满流体的单元，则F_AD设置为F_D的值；如果Donor单元是边界单元，则根据Donor单元的流体的流向和边界的取向来定出F_D的值，如果流动方向垂直于界面，则F_AD设置为F_A的值；如果相切，则F_AD设置为F_D的值；如果Acceptor单元是空单元，则F_AD设置为F_A的值；图4-a说明的Donor单元与Acceptor单元的定义，虚线标明了流走的流体体积，图4-b、图4-c、图4-d说明了F_AD的取值，黑色部分标明了实际流走的流体体积；

3)计算边界上的体积流量，计算式如下：

$F_{\mathrm{AD}}{u}_{i+\frac{1}{2},j,k}^{n+1}\mathrm{\Δt}=\mathrm{sgn}\left(u_{i+\frac{1}{2},j,k}^{n+1}\right)\mathrm{Min}\left(F_{\mathrm{AD}}\right|V_x|+\mathrm{CF},F_D\mathrm{\Δ}{x}_D)$

其中： $V_x=u_{i+\frac{1}{2},j,k}^{n+1}\mathrm{\Δt}$

CF＝Max((1.0-F_AD)|V_x|-(1.0-F_D)Δx_D，0.0)

以浇注时的浇注速度为在铸造充型过程中对带有自由表面的Navier-Stokes方程的初始条件，浇注时的液态金属温度为散热方程的初始条件。

所述自由表面边界条件的设置在SOLA-VOF法中，边界包括自由表面、型壁(即网格边界)、入流边界，结合惯性原理和连续性原理处理自由表面速度，具体是：

1)自由表面分界条件：在所述计算单元内，流体保持原有的方向和流速流动，如果遇到障碍，将出流量均匀地转移到其它可能的出流方向，接着调整出流速度，使新生成的满网格单元能够满足连续性方程；所述对于自由表面的压力条件，满足正应力和切应力为零的条件，即

式中，

为自由表面单位法向矢量；考虑三维情况下的自由表面以及压力边界条件，如图5所示，自由表面边界用直平面来近似，由网格中的体积函数F值及自由表面的法向向量可以确定自由表面的位置对于自由表面的压力，通过内部单元与外部压力线性插值的方法获得。即

          p_i，j＝(1-η)·p_N+ηp_s

其中，

$\mathrm{\η}=\frac{d_c}{d}$

为单元中心距离和表面与内部单元距离之比。

2)型壁边界条件：使用“虚拟单元”的方法，在型壁处，设置虚拟单元，通过设定虚拟单元速度，设定自由滑、无滑等边界条件；当液态金属流动时，会在边墙上形成一边界层。在边界层内的流体流动速度远较边界层外低，而此边界层的厚度和流体和流体的性质、流体的速度及边墙性质有关；如果这个边界层的厚度超过计算上所划分的格子大小时，应用无滑移边界；反之，如果这个边界层很薄，那么，采用自由滑移边界；所述无滑移的边界条件，必须符合垂直于边墙的速度为0及平行于边墙的速度为0两个条件；一个自由滑移的边墙边界，符合垂直于边墙的速度必须为0及平行于边墙的剪应力为0两个条件。

3)入流边界是指充型时浇注速度。

所述前处理任务包括与商品化软件的接口程序及CSG造型与剖分两个部分，与商品化软件ProCAST的文件接口程序用于读取由PRECAST创建的有限元网格文件，其文件格式是以数字标题加内容说明来组织格式化文件的，典型的文件格式如prefixd.dat文件，包含对于要模拟的问题的所有信息，包含很多个数据段，其中，每一数据段是以两个标号(LA，LB)开头的，标号表明该段的属性，标号后是该段的具体内容；其中，本发明选用了(1，0)段、(1，1)段、(3，0)段、(4，0)段，详细内容为：(1，0)段，标题段，

包括，LA＝1，

      LB＝0，

      TITLE＝标题信息。

(2，0)段，结点数和单元数段，

包括，LA＝1，

      LB＝1，

      NTNOD＝结点总数，

      NTEL＝单元总数。

(3，0)段，单元信息，

包括，LA＝3，

      LB＝0，

      IEL＝单元数，

      ITYPE＝单元类型，其中，

                 1，8结点长方体单元，

                 2，4结点四面体单元，

                 3，六结点契型单元，

      等等，

N＝材质号，

NODE1＝相关数据，

等等。

(4，0)结点坐标段，

包括，LA＝4，

      LB＝0，

      U1＝坐标单位，

                 1，米，

                 2，厘米，

                 3，毫米，

                 4，英尺，

NODE＝结点号，

X，Y，Z，坐标值；

所述自动剖分是按照数值计算的要求，将计算域划分成微小的六面体单元，在有限差分法中常采用六面体单元，根据铸件的特点，在铸件x、y、z轴方向上分别确定剖分步长Δx_i、Δy_i、Δz_k(i＝1，2，…，N_x；j＝1，2，…，N_y；k＝1，2，…，N_z)，依据几何形体信息，算出名剖分单元的中心点坐标，并以各单元的中心点为代表，确定出单元的物理属性(如铸件、砂型、冷铁)；如图3所示，其具体流程为：由剖分步长计算出剖分步数N_x、N_y、N_z，判断剖分步长N_x、N_y或N_z循环值计算出单元的中心点坐标P(i，j，k)，判断该点属性，然后在输出结果文件中写入对应的属性值，如果剖分步长循环值判断未完，P值未完回剖分步长循环值处继续判断，否则结束。

所述前处理任务中的CSG造型与剖分是通过基本形体的集合运算来组成复杂形体，常用的集合运算有并、交、差和补运算，在计算机上对铸造充型的形体进行实体造型用数学方法描述出基本形体的点的集合及集合间的运算，据此确定出不同物理属性、不同材质的点的集合，如铸件、砂型和砂芯等。

所述后处理任务主要是以有限差分计算的可视化为处理形式，运用MFC提供的图像处理功能，动态显示计算得到的数据，具体表现为：1)具有二维或三维显示功能；2)对图形进行交互缩放的功能；3)动态地显示解随时间的变化情况；4)可视解中应保含浓度过渡的彩色色谱，准确地显示足够多的细节，如温度梯度、等值线等，使人们能直观清晰地定性或定量地获取部分区域或某些点上的某些物理量上的值；5)实现交互实时的随视点位置变更而显示图形，通常是增强解的三维效应的最好办法；6)交互地修改颜色对照表，并实时地改变图形颜色。

实施例1

前处理：采用本发明对铝合金托盘件实现了自动剖分，托盘件的几何形状如图6-a、6-b所示。ProCAST软件进行有限元实体剖分，剖分结果显示如图7所示。使用剖分任务，剖分结果如图8所示。

后处理：利用Visual C++中提供的MFC(Microsoft Foundation Class)库，绘制出如图9所示的本发明后处理任务的一个输出实例。

实施例2

铝合金砂铸工艺模拟：前处理采用CSG法进行几何造型，自动剖分；标准测试的铝合金铸件的几何形状如图10所示；剖分时，采用变步长剖分，最大剖分步长为5mm，最小为1mm，共剖得89870(11×86×95)个单元，剖分时间大约2分钟。铸件剖分结果显示如图11所示；

液态金属从直浇道以0.8m/s的速度浇注，基于MFC的CDC(Class of Device-Context)类，采用本发明后处理显示任务，其本发明计算结果与ProCAST软件的对比见图12、图13；可见，两者的自由表面形态、流动模式十分吻合，证明本发明所采用的数学模型与算法是正确的。对于此铝合金件，伯明翰大学的研究者曾使用X-ray技术，实时观察到不同时刻的液态金属的流动模式，本发明计算的结果与观察到的模式基本符合。

实施例3

IN738合金真空精密铸造模拟：IN738合金是一种高温合金，采用真空精密铸造生产，但由于是薄壁板类件，经常产生浇不足或者疏松等缺陷，多次实验均未能解决。原始工艺一般是浇注温度1400℃，壳温850℃，剖分后共92×87×80＝640320个单元，浇注速度为0.8m/s，IN738各项物性参数如表2；

                                                  表2.IN738合金的物性参数

导热率(w/m/k)	密度(kg/m3)	比热kJ/kg/K	固相线温度(℃)	液相线温度(℃)	动力学粘度系数Pa·s
						20	7800	550	1260	1330	0.01

本发明与ProCAST软件的计算结果对比情况如图14、图15、图16、图17所示，两者吻合的相当好。其中：图14是本发明在t＝0.132s的充型与温度场模拟结果，而图15是ProCAST软件在t＝0.124s时的模拟结果，对比可见，自由表面形态与充满程度基本是一致的。在充型初始的一瞬间，由于激冷作用，铸件表面温度均略有下降；图16、图17是t＝0.168s和t＝0.147s时本发明和ProCAST的模拟结果，在距外表面1/3处，由于液流回转，易出现夹杂等，实验证明也是疏松位置。根据模拟结果，改进浇注方案，将浇注温度提高到1550℃，壳温提高到870℃-900℃，解决了充型问题，避免了浇不足等缺陷。同时由于浇注温度与壳温相差较大，因而在凝固前沿温度梯度大，形成由板件外表面向内顺序凝固方式，有利于补缩，也避免了疏松等缺陷。通过模拟，改进了浇注工艺参数，提高了合格率。

本发明总体框架采用Visual C++提供的文档/视图结构，文档/视图结构实现了数值计算与结果处理的分离；对SOLA-VOF算法用自定义的SOLA-VOF算法类(Csolavof)进行了封装，易于此算法的调试、升级、维护与移植。对散热方程的求解在自定义的Ctemperature类中实现，它继承了Csolavof类的所有成员变量和函数。本发明是标准的Windows 32程序，它的运行是依靠外部发生的事件来驱动的。换句话说，程序等待用户输入，然后做判断并做出适当的处理。这里的输入是由操作系统捕捉到之后，以消息形式(一种数据结构)进入到程序之中。归纳起来，本发明是依靠消息驱动的、使用了面向对象编程技术的、标准的Windows 32位程序。

本发明用面向对象语言Visual C++语言编写，在Win9X、WinNT操作系统上运行。

Claims (7)

1.一种铸件充型过程模拟方法，其特征在于：它是通过建立并求解描述液态金属流动过程的控制方程，得到充型过程中压力场、速度场、温度场以及自由表面的定量变化，从而选择正确的浇注方法、控制成形中缺陷的产生；由计算、前处理、后处理三个任务组成，所述计算任务是完成对铸件充型过程中带有自由表面的Navier-Stokes方程及散热方程的求解，对该方程的求解借鉴了Los Alamos实验室的报告中二维SOLA-VOF算法，将其推广到了三维，求出流体场之后，求解散热方程；所述前处理任务是准确地构造铸件的三维几何模型，并进行网格剖分；所述后处理任务的功能是把由数值计算或试验获得的大量数据转换为人的视觉可以感受的图像，把大量的抽象数据有机地组合在一起，展示数据所表现的内容及之间的相互关系，展示场的分布；总流程为：首先选择任务类型，选计算任务类型时，先初始化自主开发的温度场计算类，如果是第一次开始，则输入剖分文件、铸造充型过程的控制及物理参数，进行计算后回主菜单；否则输入上次计算结果，再计算，计算后回主菜单；选前处理任务类型时，判断何种类型几何文件，经判断，分别对铸件的构造三维实体几何模型的文件或有限元网格文件读取数据，进行剖分，保存剖分结果，然后回主菜单；选后处理任务类型时，则分别显示温度场、速度场、压力，回主菜单。

2.按照权利要求1所述铸件充型过程模拟方法，其特征在于：所述计算任务中计算部分是用SOLA求解铸件充型过程压力场和速度场，用VOF确定流动域和自由表面；采用施体-受体方法计算铸造充型过程自由表面的扩展；具体流程为：先输入控制、物理性质参数，然后依据施体-受体算法计算自由表面边界条件，再计算自由表面方向，用t＝t+Δt累加时间值，更新速度场，构造铸件的三维几何模型已充满或者完成，当前计算则结束，否则以上一时刻的速度场、压力场作为估计值，显示计算本时刻的速度值，求出单元的散度D_i，j，k，并判断是否满足连续方程，如果所述判断为肯定结果，则直接求解散热方程，再回到依据施体-受体算法计算自由表面扩展步骤；如果所述判断为否定结果，则修正压力，修正速度后，再重新判断，直到判断结果满足连续方程，至下一步求解散热方程。

3.按照权利要求2所述铸件充型过程模拟方法，其特征在于：所述用SOLA求解铸件充型过程压力场和速度场、用VOF确定流动域和自由表面，其计算步骤为：

(1)由Navier-Stokes方程的离散公式，以初始条件或前一时刻值为基础，计算当前时刻的试算速度；

(2)给每一个单元定义散度D_i，j，k：将所述步骤(1)的试算速度值代入连续性方程离散公式，求出D_i，j，k，若D_i，j，k＝0，则步骤(1)中所述试算速度满足连续性方程，即此时的速度场与压力场值既满足动量守恒方程又满足质量守恒方程，当前时间步长计算结束；如整个流场中有任一个单元不能满足连续性方程，即D_i，j，k不等于0，则需要下一步的修正；修正步骤：欲修正速度，须先修正压力；

(3)根据校正压力求出校正后的试算速度；

(4)将校正后的试算速度值代入所述步骤(2)，反复迭代直至所有单元均满足连续性方程；

(5)由体积函数方程确定新的流动域，对表面单元做合理设置；

(6)返回第一步，进入下一时刻计算。

4.按照权利要求2所述铸件充型过程模拟方法，其特征在于：所述计算自由表面扩展设计方法为根据计算单元的上游、下游单元，粗略估计出自由表面的形状，然后，根据不同形状，计算体积函数F在单元的边界面上的通量，从而计算出各个时刻体积函数F值；具体步骤是：

1)边界上的任意两个单元，区分出施体单元和受体单元，依据所述两个单元边界上的速度的代数符号来确定，施体单元是上游单元，其体积函数用F_D表示，受体是下游单元，其体积函数用F_A表示；

2)定出施体单元和受体单元边界上的体积函数值F_AD；

3)计算边界上的体积流量。

5.按照权利要求2或4所述铸件充型过程模拟方法，其特征在于：所述自由表面边界条件的设置在SOLA-VOF法中，边界包括自由表面、型壁即网格边界、入流边界，结合惯性源理和连续性原理处理自由表面速度，具体是：

1)自由表面分界条件：在所述计算单元内，流体保持原有的方向和流速流动，如果遇到障碍，将出流量均匀地转移到其它可能的出流方向，接着调整出流速度，使新生成的满网格单元能够满足连续性方程；所述对于自由表面的压力条件，满足正应力和切应力为零的条件；

2)型壁边界条件：使用“虚拟单元”的方法，在型壁处，设置虚拟单元，通过设定虚拟单元速度，设定自由滑、无滑等边界条件；

3)入流边界是指充型时浇注速度。

6.按照权利要求1所述铸件充型过程模拟方法，其特征在于：所述前处理任务包括与商品化软件的接口程序及构造实体几何模型造型与剖分两个部分，与商品化软件ProCAST的文件接口程序用于读取有限元网格文件，其具体流程为：读取所述由剖分步长计算出剖分步数，判断剖分步长N_x、N_y或N_z循环值计算出单元的中心点坐标P(i，j，k)，判断该点属性，然后在输出结果文件中写入对应的属性值，如果剖分步长循环值判断未完，P值未完回剖分步长循环值处继续判断，否则结束；所述前处理任务中的构造实体几何模型造型与剖分是通过基本形体的集合运算来组成复杂形体，常用的集合运算有并、交、差和补运算，在计算机上对铸造充型的形体进行实体造型用数学方法描述出基本形体的点的集合及集合间的运算，据此确定出不同物理属性、不同材质的点的集合，如铸件、砂型和砂芯等。

7.按照权利要求1所述铸件充型过程模拟方法，其特征在于：所述后处理任务主要是以有限差分计算的可视化为处理形式，运用Microsoft Foundation Class提供的图像处理功能，动态显示计算得到的数据，具体表现为：1)具有二维或三维显示功能；2)对图形进行交互缩放的功能；3)动态地显示解随时间的变化情况；4)准确地显示足够多的细节，如温度梯度、等值线，使人们能直观清晰地定性或定量地获取部分区域或某些点上的某些物理量上的值；5)实现交互实时的随视点位置变更而显示图形；6)交互地修改颜色对照表，并实时改变图形颜色。

Images