王宁等：基于无网格伽辽金法的连铸坯凝固计算方法 187 KEY WORDS element-free Galerkin method;moving least squares;heat transfer;solidification;continuous casting 结晶器内钢液温度几乎无法直接测量，数值 +Q=pcT 模拟是研究结晶器内铸坯凝固行为的重要途径山 Tlt=0=Tin (1) 随着研究问题的深入，基于网格划分的数值方法 .aT .OT 无法精确重构凝固坯壳形貌和液/固相区，对于裂 -nx+kny=q y 纹和大变形等网格复杂变化的问题甚至无法求 式中，k为导热系数，Wm1.K-l;T和Tn分别为温 解.无网格方法用离散的节点代替网格结构，在自 度和初始浇铸温度，℃；T为温度对于时间的导数； 适应、求解裂纹和大变形等问题体现出显著优势， p为密度，kgm3;c为比热容，Jkg1.K-l;n和n,分 已经应用于流体力学、电磁学、爆炸等-多个领 别为沿x和y方向边界外法线的方向余弦；q为结晶 域，也有少数研究将其应用于连铸过程.Zhang等阿 器/俦坯界面上的热流密度，W,m2,采用Savage 用有限点法和无网格局部彼得罗夫伽辽金方法分 与Pritchard提出的经验公式ll:g=A-BVf,t表示 析了方坯凝固过程中的热弹塑性问题：Vertnik与 时间，S,A和B为与结晶器有关的系数； Sarler忉利用局部径向基函数配点法预测了钢液的 dfs dfs dT 湍流流动和传热行为；Yamasaki等图用基于粒子 O=pL (2) dT dr 的MPS法模拟二冷区复杂的喷洒水流：Vaghefi等 式中，Q为内热源，Jm3s,为固相率，L为钢 用MLPG法分析了坯壳的热弹塑性情况；Hostos 液的凝固潜热，Jkg 等uo用无网格伽辽金方法(element--free Galerkin 将式(2)代入式(1)中得到： method,.EFG)计算圆坯的传热和力学行为.但是 ) a/.ar\ 与有限元、有限差分等模拟方法相比，无网格方法 =pc,其中等效热容c=c-L） dfs 在连铸过程的研究仍然较少，且已有研究均采用 1.2变分处理 均匀布置的节点离散计算区域，暂无针对无网格 取温度变分δT为任意函数，利用分部积分法， 方法节点布置的灵活性及其与传统模拟方法相互 建立能量泛函，得到等效积分弱形式为 验证的研究，而无网格方法最大的优势在于其节 点布置的灵活性，因此探讨多种节点布置方法的 g+tg+actbo-w7dr=0 计算精度是十分必要的 (3) 本文依据EFG法基本原理和铸坯传热凝固控 将求解域Ω离散为有限的N个节点，求解域 制方程，建立了基于EFG法的小方坯凝固过程数 内1时刻任意位置(xy)处的温度T(x,y,)可用邻近 值模型，运用C+语言自行设计和开发了面向对 的n个节点的温度与形函数的乘积加和来表示， 象的铸坯凝固过程无网格数值计算软件.探讨多 得到式(4)： 种节点布置方式下的无网格与有限元方法在计算 精确度、自适应性、网格依赖性等方面的差异，为 Ty0=∑pI, 8T= >Φ8T1(i) 无网格方法求解连铸过程的传热、凝固以及后续 (4) 的应力、多相耦合与界面追踪等问题提供基础 =O(.y.!) 51 .oTi(1) 1铸坯传热/凝固模型的建立 式中，Φ为节点1的形函数，将式(4)代入式(3)，整 依据铸坯的凝固特点，建立模型前作如下假 理得到式(5)： 设：沿拉坯方向传递的热量只占总热量的3%~6%， CNXNTNxI+KNxNTNxI FNxI (5) 忽略铸坯纵向传热，简化为二维非稳态传热过程叫： 式中，Ku=gk.r中，7Φd2,C切=∫QPCet·ΦΦd2, 铸坯顶部和底部绝热；忽略弯月面以上保护渣吸 F1=∫，-q1dr,V为那勃勒算子. 收的热量：忽略结品器的锥度和弧度以及铸坯凝 通过求解式(5)可以得到连铸坯二维横截面 固产生的收缩变形：铸坯热物性参数各向同性 的温度分布，沿浇铸方问，采用以下无条件稳定的 1.1传热控制模型 向后差分法进行离散求解吲，其中△表示时间步长： 连铸坯凝固过程二维非稳态传热方程)： (C++△tK+)T+=△tF+A+C+△T, (6)KEY WORDS    element-free Galerkin method；moving least squares；heat transfer；solidification；continuous casting 结晶器内钢液温度几乎无法直接测量，数值 模拟是研究结晶器内铸坯凝固行为的重要途径[1] . 随着研究问题的深入，基于网格划分的数值方法 无法精确重构凝固坯壳形貌和液/固相区，对于裂 纹和大变形等网格复杂变化的问题甚至无法求 解. 无网格方法用离散的节点代替网格结构，在自 适应、求解裂纹和大变形等问题体现出显著优势， 已经应用于流体力学、电磁学、爆炸等[2−5] 多个领 域，也有少数研究将其应用于连铸过程. Zhang 等[6] 用有限点法和无网格局部彼得罗夫伽辽金方法分 析了方坯凝固过程中的热弹塑性问题；Vertnik 与 Šarler[7] 利用局部径向基函数配点法预测了钢液的 湍流流动和传热行为；Yamasaki 等[8] 用基于粒子 的 MPS 法模拟二冷区复杂的喷洒水流；Vaghefi 等[9] 用 MLPG 法分析了坯壳的热弹塑性情况；Hostos 等 [10] 用无网格伽辽金方 法 (element-free Galerkin method，EFG) 计算圆坯的传热和力学行为. 但是 与有限元、有限差分等模拟方法相比，无网格方法 在连铸过程的研究仍然较少，且已有研究均采用 均匀布置的节点离散计算区域，暂无针对无网格 方法节点布置的灵活性及其与传统模拟方法相互 验证的研究，而无网格方法最大的优势在于其节 点布置的灵活性，因此探讨多种节点布置方法的 计算精度是十分必要的. 本文依据 EFG 法基本原理和铸坯传热凝固控 制方程，建立了基于 EFG 法的小方坯凝固过程数 值模型，运用 C++语言自行设计和开发了面向对 象的铸坯凝固过程无网格数值计算软件. 探讨多 种节点布置方式下的无网格与有限元方法在计算 精确度、自适应性、网格依赖性等方面的差异，为 无网格方法求解连铸过程的传热、凝固以及后续 的应力、多相耦合与界面追踪等问题提供基础. 1    铸坯传热/凝固模型的建立 依据铸坯的凝固特点，建立模型前作如下假 设：沿拉坯方向传递的热量只占总热量的 3% ~ 6%， 忽略铸坯纵向传热，简化为二维非稳态传热过程[11] ； 铸坯顶部和底部绝热；忽略弯月面以上保护渣吸 收的热量；忽略结晶器的锥度和弧度以及铸坯凝 固产生的收缩变形；铸坯热物性参数各向同性. 1.1    传热控制模型 连铸坯凝固过程二维非稳态传热方程[12] ：    ∂ ∂x ( k ∂T ∂x ) + ∂ ∂y ( k ∂T ∂y ) + Q = ρcT˙ T| t=0 = Tin − ( k ∂T ∂x nx +k ∂T ∂y ny ) = q （1） k W·m−1 ·K −1 T Tin T˙ ρ kg ·m−3 J· kg−1 ·K −1 nx ny x y W·m−2 q = A− B √ t t 式中， 为导热系数， ； 和 分别为温 度和初始浇铸温度，℃； 为温度对于时间的导数； 为密度， ；c 为比热容， ； 和 分 别为沿 和 方向边界外法线的方向余弦；q 为结晶 器/铸坯界面上的热流密度， ，采用 Savage 与 Pritchard 提出的经验公式[13] ： , 表示 时间，s，A 和 B 为与结晶器有关的系数； Q = ρL d fs dt = ρL d fs dT dT dt （2） Q J·m−3 ·s −1 fs L J· kg−1 式中， [14] 为内热源， ， 为固相率， 为钢 液的凝固潜热， . ∂ ∂x ( k ∂T ∂x ) + ∂ ∂y ( k ∂T ∂y ) = ρceff ∂T ∂t ceff= ( c− L d fs dT ) 将 式 （ 2） 代 入 式 （ 1） 中 得 到 ： ，其中等效热容 . 1.2    变分处理 取温度变分 δT 为任意函数，利用分部积分法， 建立能量泛函，得到等效积分弱形式为[15] ： w Ω [ k ∂T ∂x ∂δT ∂x +k ∂T ∂y ∂δT ∂y +ρceffT˙δT ] dΩ+ w Γ2 qδTdΓ = 0 （3） Ω T (x, y,t) 将求解域 离散为有限的 N 个节点，求解域 内 t 时刻任意位置 (x,y) 处的温度 可用邻近 的 n 个节点的温度与形函数的乘积加和来表示， 得到式 (4)： T (x, y,t) = ∑n I=1 ΦITI (t), δT = ∑n I=1 ΦIδTI (t) T˙ = ∂T (x, y,t) ∂t = ∑n I=1 ΦI ∂TI (t) ∂t （4） 式中， ΦI 为节点 I 的形函数，将式（4）代入式 (3)，整 理得到式 (5)： CN×NT˙N×1 + KN×NTN×1 = FN×1 （5） KIJ = r Ω k · ∇TΦI ∇ΦJdΩ, CIJ = r Ω ρceff ·ΦT I ΦJdΩ, FI = r Γ2 −q ·ΦIdΓ ∇ 式中， ， 为那勃勒算子. ∆t 通过求解式（5）可以得到连铸坯二维横截面 的温度分布，沿浇铸方向，采用以下无条件稳定的 向后差分法进行离散求解[15] ，其中 表示时间步长： (Ct+∆t + ∆tKt+∆t)Tt+∆t = ∆tFt+∆t +Ct+∆tTt （6） 王    宁等： 基于无网格伽辽金法的连铸坯凝固计算方法 · 187 ·