CN108388727A - 一种适用于凝固过程中辐射换热的计算方法 - Google Patents

Abstract

本发明属于铸造中的凝固过程数值模拟热领域，并公开了一种适用于凝固过程中辐射换热的计算方法。该方法包括下列步骤：(a)将辐射换热的物体离散为多个辐射单元，构建以辐射单元为中心的半球作为辐射换热模型，从球心发出的多条辐射射线是将所述半球分割为多个子空间，每条辐射射线代表其相应的子空间的辐射换热；(b)每个辐射单元上中的辐射射线的方向向量；(c)求解每个子空间的辐射权重因子；(d)计算每个子空间的辐射换热和辐射单元的总辐射换热，构建辐射单元辐射换热的边界条件，采用有限元求解的方法求解边界条件。通过本发明可以准确地实时处理铸造凝固过程中的辐射换热，提高凝固过程数值模拟的准确性。

Description

一种适用于凝固过程中辐射换热的计算方法

技术领域

本发明属于铸造中的凝固过程数值模拟热领域，更具体地，涉及一种 适用于凝固过程中辐射换热的计算方法。

背景技术

我国是全球最大的铸件生产国家，铸造对国民经济的发展起着重要的 作用，是基础工业的一项主要支柱，近年来，随着计算机模拟技术的日益 成熟，数值模拟手段在优化铸造宏观过程方面的作用越来越重要。

铸造凝固过程数值模拟往往涉及到辐射换热的计算，而辐射换热边界 相比与其他换热边界(定温、对流、给定热流量)，其复杂之处在于辐射 现象是三维空间分布，使得辐射换热本身就需要进行离散化处理，而且由 于辐射换热不需要空间介质，使得计算辐射换热前必须进行辐射对应关系 的查找即射线追踪。在定向凝固数值模拟中，由于铸件与型壳会进行抽拉 运动以及定向凝固的复杂工艺(辐射挡板、水冷结晶器、液态金属液等)， 使得辐射换热边界的处理更为困难。目前已有辐射换热表征方法往往不进 行射线追踪操作，而简单地将辐射换热等效为对流换热，降低了凝固过程 数值模拟的准确性。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种适用于凝固 过程中辐射换热的计算方法，采用构建以辐射单元为中心的半球作为辐射 换热模型，每个子空间的辐射换热由对应的射线代替，通过追踪各个射线 的投射位置以及投射材质，从而获得射线投射的终点表面单元的温度，从 而准确地计算辐射换热对凝固过程的影响，由此解决凝固过程中的辐射换 热准确性低的技术问题。

为实现上述目的，按照本发明，提供了一种适用于凝固过程中辐射换 热的计算方法，其特征在于，该方法包括下列步骤：

(a)将辐射换热的物体离散为多个辐射单元，对于每个辐射单元，构 建以辐射单元为中心的半球作为辐射换热模型，从球心发出的多条辐射射 线将所述半球分割为多个子空间，每条辐射射线代表其相应的子空间的辐 射换热；

(b)在每个辐射单元上选取该辐射单元的几何中心作为坐标原点，该 辐射单元的法向量方向作为Z轴正方向构建三维空间坐标系，根据三维空 间方向向量的旋转求解该辐射单元上的辐射射线的方向向量；

(c)根据步骤(b)获取的辐射射线的方向向量，获得该辐射射线所 代表的子空间在所述三维空间坐标系中分别与x轴的夹角和z轴的夹角θ， 利用该和θ分别求解每个子空间在其相应的辐射单元中的辐射权重因子；

(d)根据每个子空间的辐射权重因子计算每个子空间的辐射换热和相 应的辐射单元的总辐射换热，根据该总辐射换热和辐射权重因子构建该辐 射单元辐射换热的边界条件，采用有限元求解的方法求解边界条件，由此 获得此获得由辐射引起的热量交换。

优选地，在步骤(b)中，所述辐射射线的方向向量优选采用下列表达 式，

其中，

R(t)是辐射射线的方向向量，P是辐射单元的几何中心，t是辐射射线的 位置自变量，是辐射单元的单位法向向量，是在辐射单元中构建的空间 三维直角坐标系沿x方向的单位向量，是是在辐射单元中构建的空间三维 直角坐标系沿y方向的单位向量。

优选地，在步骤(c)中，所述辐射权重因子优选采用下列表达式，

其中，为射线在x-y平面上的投影与x轴的夹角，是的间隔，θ是 辐射射线与z轴的夹角，Δθ是θ的间隔。

优选地，在步骤(c)中，所述相应的辐射单元的总辐射换热Q_rad优选 采用下列表达式，

其中，Q_j是第j个子空间的辐射换热，n是子空间额总数量，α_j是辐射 权重因子，ε和ε_j分别是辐射射线起点和终点所在表面单元的黑度，S和S_j分 别是辐射射线起点和终点所在表面单元的面积，u和u_j分别是辐射射线起点 和终点所在表面单元的温度，σ为Stefan-Boltzman常数。

优选地，在步骤(d)中，所述边界条件优选采用下列表达式，

其中，

k是导热系数，n_x是辐射射线起点所在表面单元外法向向量x向分量，n_y是辐射射线起点所在表面单元外法向向量方向向量y向分量，n_z是辐射射 线起点所在表面单元外法向向量方向向量z向分量，u_e是辐射射线起点所 在表面单元的平均温度。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够 取得下列有益效果：

1、本发明首先通过对辐射换热物体进行物理上的离散处理获得多个辐 射单元，然后通过显式计算各个射线的等效辐射换热，准确地表征整体辐 射换热效果；

2、本发明通过采用辐射换热边界的有限元计算模型，该边界条件简单 直观，贴切建立的辐射换热模型，使得辐射换热模型数值化，直观地进行 数值求解；

3、本发明通过借助计算图形学知识和几何关系，获得辐射换热处理所 需的三维空间内辐射射线方向向量，可以简便地进行辐射换热所需的物理 离散。

附图说明

图1是按照本发明的优选实施例所构建的辐射换热的计算方法的流程 图；

图2是按照本发明的优选实施例所构建的辐射换热的计算方法的流程 图在三维直角坐标下的辐射换热模型射线的定义；

图3是按照本发明的优选实施例所构建的点和方向向量相对于任意轴 旋转的情形示意图；

图4是按照本发明的优选实施例所构建的显示在平面上的一般旋转 示意图；

图5是按照本发明的优选实施例所构建的三维空间的辐射线方向向量 示意图；

图6是按照本发明的优选实施例所构建的LMC工艺整体几何模型和网 格模型；

图7是按照本发明的优选实施例所构建的采用所提出的辐射换热表征 方法计算的铸件温度场分布(单位：℃)。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图 及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体 实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的 本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可 以相互组合。

图1是按照本发明的优选实施例所构建的辐射换热的计算方法的流程 图，如图1所示，一种适用于凝固过程计算的辐射换热表征方法，该方法 包括下列步骤：

(a)对辐射换热进行物理上的离散处理，即将上半球空间分割为多个 子空间，每个子空间的辐射换热由对应的射线代替，通过追踪各个射线的 投射位置以及投射材质来处理辐射换热。

(b)通过有限元理论，发现辐射换热边界条件与对流换热边界条件有 着相同的形式，获得辐射换热边界的计算模型。

(c)在三维空间内方向向量的旋转基础上，获得辐射换热处理所需的 三维空间内辐射射线方向向量。

该方法的具体步骤如下：

(1)辐射换热物理离散

由于辐射换热是三维空间内密布的，因此需要先对辐射换热进行物理上 的离散处理。如图2所示，对于单个表面辐射单元i而言，认为其与外界的 辐射换热通过单元_i的几何中心向上半球发射射线进行简化计算。因此将上 半球空间分割为多个子空间，每个子空间的辐射换热由对应的射线代替， 通过追踪各个射线的投射位置以及投射材质(型壳、炉体、辐射挡板、水 冷结晶器等)，从而获得射线投射的终点表面单元j及其温度u_j，认为此射线所代表的子空间的投射单元的温度均为u_j，且射线所代表的子空间的能 量全部投射到终点表面单元上，即对应的角系数

(2)辐射射线方向向量的求解

(2.1)三维空间内方向向量的旋转

射线追踪过程中需要确定三维空间内不同射线的方向向量，该方向向 量可以通过相关的旋转操作获得。图3表示了点和方向向量相对一个任意 轴旋转的一般情形。

下面是其中用到的符号的说明：

Q,定义旋转轴的点和单位向量

θ旋转角度

P被旋转的点

T(P)旋转得到的点

被旋转的方向向量

旋转得到的方向向量

垂直于的平面

在上的投影

考虑一个相对单位向量(它是经过Q的一条旋转轴)逆时针方向旋转 并且旋转角度为θ(右手法则)的旋转。向量可以相对于另一个向量分解为 平行和垂直分量，我们在此将投影到得到和注意，由于方向向 量与位置无关，我们可将画成从旋转轴上的Q开始，以使图形更易于理解。

为了有助于理解以下的内容，参见图4，其中显示了垂直于且包含 的平面

据此可以做如下的推论：

由于并且T是线性变换，因此可得：

平行和垂直向量分量的定义如下：

分别将式(29)和式(30)代入式(26)和式(27)可得：

将式(31)和式(32)代入式(28)可得：

(2.2)三维空间内辐射射线方向向量的求解

在三维空间内方向向量的旋转基础上，可以获得辐射换热处理所需的 三维空间内辐射射线方向向量。图5为三维空间的辐射线方向向量示意图。

对于单个表面单元，表面单元的3个有序节点为{V₀,V₁,V₂}，P为表面单 元的几何中心，为表面单元的单位法向向量，则有：

单位向量与平行，单位向量与构成局部直角坐标系，则有：

为表面单元所发射的某一条射线，与表面单元(即)夹角为θ，在 x-y平面的投影与x轴的夹角根据三维空间中方向向量的旋转理论，可 以认为由旋转得到：以为旋转轴旋转得到方向向量 以为旋转轴旋转得到射线方向向量 为在x-y平面的投影。则有：

在此已经求出了对应的射线起点P以及射线方向向量则有：

针对表面单元上半球所辐射的射线，可知实际上 需要将θ和离散化取值，考虑到(见辐射换热处理方式)。将Δθ和 取值为于是θ取3个离散值，取12个离散值，即：

那么易知对于每个表面单元而言，其一共有36条辐射射线，每对取 值对应一条辐射射线的射线方程。

(3)根据求得的辐射射线的方向向量获得每条射线的方向，即在三维 坐标中与x轴夹角与z轴夹角θ；

每条射线所代表的子空间可由x轴夹角水平面间隔z轴夹角θ和 垂直面间隔Δθ决定，结合表面单元i的温度u，可以得到子空间内表面单元 i与外界的辐射换热Q_j，将表面单元i所发出的每条射线所代表的子空间的辐 射换热进行叠加，即可得到表面单元i与外界的总辐射换热Q_rad，设所分割的 子空间个数为n，则有：

其中，ε和ε_j分别为射线起点及射线终点所在表面单元的黑度；S和S_j分 别为射线起点及射线终点所在表面单元的面积，m²；σ为Stefan-Boltzman 常数，J/(s·m²·K⁴)；α_j为辐射权重因子，用来表征射线所携带的辐射占表 面单元i向外总辐射的比值。

对于α_j的求解过程如下：对于定位在局部坐标系x-y平面的表面单元 i，根据Lambert定律，其在空间上的定向辐射力为：

E_θ＝E_n·cosθ (3)

其中，E_θ为与表面单元夹角为θ的射线所具有的定向辐射力， W/(m²·sr)；E_n为法线方向射线的辐射力，W/(m²·sr)。

在球坐标系内，子空间上一个微元的面积为：

其中，r为微元与表面单元i的距离，m；θ为射线与z轴的夹角；为 射线在x-y平面上的投影与x轴的夹角。

每个子空间的总辐射力E_j由每条射线通过能量积分得到：

表面单元i在所有子空间上的总辐射力为：

辐射权重因子α_j表示E_j占表面单元i总辐射力E的比例，可得：

(4)辐射换热有限元处理

辐射换热只是针对表面单元而言，因此换热的处理实际上是边界条件 的处理，相应的辐射换热边界条件为：

将式(1)、式(2)以及式(7)代入式(8)，可得：

令：

于是，式(9)可改写为：

对式(12)进行线性化，将u在u_e附近展开，并取一次项：

其中，u_e为表面单元i的平均温度。于是式(12)可改写为：

令：

于是，式(14)可改写为：

(在边界Γ_r上) (17)

其中：

对式(18)和式(19)仔细分析，针对特定表面单元、特定时刻，两 式中的各个物理量均为已知的，即针对特定时刻的特定表面而言，D₀与D₁为 已知量。需要说明的是：上述公式计算中温度采用绝对温度，所以D₀和D₁的 计算中温度为绝对温度，所以式(17)改写为：

采用有限元的方法，计算上述表达式(20),由此完成辐射单元辐射换 热的计算。

不难发现辐射换热边界条件与对流换热边界条件有着相同的形式，易 知考虑辐射换热的非稳态温度场的有限元求解方程为：

各矩阵的单元集成形式为：

其中，是单元辐射热交换对热传导矩阵的修正；是单元的辐射换 热边界的温度载荷，T是计算域的温度场矩阵，矩阵中的每个元素代表划分 网格后每个节点的温度。其中矩阵的元素由下式给出：

应用实例

为了说明本发明所提出的辐射换热表征方法的实用性，计算一定向凝 固工艺的凝固过程。图6为该定向凝固工艺的整体几何模型和网格模型。 其中，部件包括：铸件、型壳、辐射挡板以及炉体(加热区和冷却区)，其 中型壳网格是在铸件网格的基础上随型生长得到，型壳厚度为5mm，炉体划 分面网格即可。

图7为采用所提出的辐射换热表征方法计算的铸件温度场分布。从模 拟结果可以看出，所提出的辐射换热表征方法可以准确地实时处理铸造凝 固过程中的辐射换热，提高凝固过程数值模拟的准确性。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已， 并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等 同替换和改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种适用于凝固过程中辐射换热的计算方法，其特征在于，该方法包括下列步骤：

(a)将辐射换热的物体离散为多个辐射单元，对于每个辐射单元，构建以辐射单元为中心的半球作为辐射换热模型，从球心发出的多条辐射射线将所述半球分割为多个子空间，每条辐射射线代表其相应的子空间的辐射换热；

(b)在每个辐射单元上选取该辐射单元的几何中心作为坐标原点，该辐射单元的法向量方向作为Z轴正方向构建三维空间坐标系，根据三维空间方向向量的旋转求解该辐射单元上的辐射射线的方向向量；

(c)根据步骤(b)获取的辐射射线的方向向量，获得该辐射射线所代表的子空间在所述三维空间坐标系中分别与x轴的夹角和z轴的夹角θ，利用该和θ分别求解每个子空间在其相应的辐射单元中的辐射权重因子；

(d)根据每个子空间的辐射权重因子计算每个子空间的辐射换热和相应的辐射单元的总辐射换热，根据该总辐射换热和辐射权重因子构建该辐射单元辐射换热的边界条件，采用有限元求解的方法求解该边界条件，由此获得由辐射引起的热量交换。

2.如权利要求1所述的一种适用于凝固过程中辐射换热的计算方法，其特征在于，在步骤(b)中，所述辐射射线的方向向量优选采用下列表达式，

其中，

R(t)是辐射射线的方向向量，P是辐射单元的几何中心，t是辐射射线的位置自变量，是辐射单元的单位法向向量，是在辐射单元中构建的空间三维直角坐标系沿x方向的单位向量，是是在辐射单元中构建的空间三维直角坐标系沿y方向的单位向量。

3.如权利要求1或2所述的一种适用于凝固过程中辐射换热的计算方法，其特征在于，在步骤(c)中，所述辐射权重因子优选采用下列表达式，

其中，为射线在x-y平面上的投影与x轴的夹角，是的间隔，θ是辐射射线与z轴的夹角，Δθ是θ的间隔。

4.如权利要求1-3任一项所述的一种适用于凝固过程中辐射换热的计算方法，其特征在于，在步骤(c)中，所述相应的辐射单元的总辐射换热Q_rad优选采用下列表达式，

其中，Q_j是第j个子空间的辐射换热，n是子空间额总数量，α_j是辐射权重因子，ε和ε_j分别是辐射射线起点和终点所在表面单元的黑度，S和S_j分别是辐射射线起点和终点所在表面单元的面积，u和u_j分别是辐射射线起点和终点所在表面单元的温度，σ为Stefan-Boltzman常数。

5.如权利要求1-4任一项所述的一种适用于凝固过程中辐射换热的计算方法，其特征在于，在步骤(d)中，所述边界条件优选采用下列表达式，

其中，

k是导热系数，n_x是辐射射线起点所在表面单元外法向向量x向分量，n_y是辐射射线起点所在表面单元外法向向量方向向量y向分量，n_z是辐射射线起点所在表面单元外法向向量方向向量z向分量，u_e是辐射射线起点所在表面单元的平均温度。