CN116432482A - 一种铜管水平连铸温度场在线计算方法 - Google Patents

一种铜管水平连铸温度场在线计算方法 Download PDF

Info

Publication number
CN116432482A
CN116432482A CN202310197215.6A CN202310197215A CN116432482A CN 116432482 A CN116432482 A CN 116432482A CN 202310197215 A CN202310197215 A CN 202310197215A CN 116432482 A CN116432482 A CN 116432482A
Authority
CN
China
Prior art keywords
heat
temperature
node
heat transfer
nodes
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202310197215.6A
Other languages
English (en)
Inventor
刘劲松
李旺
梁海成
崔海涛
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shenyang Ligong University
Original Assignee
Shenyang Ligong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shenyang Ligong University filed Critical Shenyang Ligong University
Priority to CN202310197215.6A priority Critical patent/CN116432482A/zh
Publication of CN116432482A publication Critical patent/CN116432482A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/14—Pipes
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/08—Thermal analysis or thermal optimisation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Investigating Or Analyzing Materials Using Thermal Means (AREA)

Abstract

一种铜管水平连铸温度场在线计算方法，包括：连铸温度场在线计算，几何模型建立及离散化，材料属性设定与实现，边界条件设定与实现，面属性设定与实现，运动传热模型建立，相变潜热模型建立，平衡方程组求解。本发明优点在于：围绕铜管水平连铸稳态传热的有限差分法建模，讨论了几何模型的建立及离散化、材料属性、边界条件、界面属性的设定与实现等问题。为描述连铸过程中相对运动及相变潜热对于传热的影响，通过物质状态变化的微元分析建立了相应数学模型。鉴于铜管水平连铸模型涉及非线性热平衡方程，提出了一种非线方程组的迭代解法，水平连铸稳态温度场在线计算程序，确保在线计算程序的传热计算可靠性。

Description

一种铜管水平连铸温度场在线计算方法
技术领域
本发明涉及铜管水平连铸温度场计算领域，特别涉及一种铜管水平连铸温度场在线计算方法。
背景技术
铜管水平连铸包含流体流动、传热、晶体形核及生长等诸多物理过程，其中传热过程处于核心地位，支配着其它物理过程的进行；此外，在线应用要求计算速度尽能能快，故仅保留传热计算，其它物理过程对传热的影响采用适当的数学近似进行描述。
发明内容
本发明的目的是为了克服上述问题，特提供了铜管水平连铸温度场在线计算方法。
基于有限差分法求解稳态传热问题的基本建模原理；围绕铜管水平连铸稳态传热的有限差分法建模，讨论了几何模型的建立及离散化、材料属性、边界条件、界面属性的设定与实现等问题；为描述连铸过程中相对运动及相变潜热对于传热的影响，通过物质状态变化的微元分析建立了相应数学模型；鉴于铜管连铸模型涉及非线性热平衡方程，介绍了一种非线方程组的迭代解法。
本发明提供了一种铜管水平连铸温度场在线计算方法，具体包括以下内容：
连铸温度场在线计算：
建模原理
在实验研究的基础上，傅立叶于1822年建立了“热的解析理论”，奠定了热传导的理论基础[70]。实验观察发现，单位时间内通过指定截面的导热量，正比于沿该截面法向的温度变化率和截面面积，从而有如下经验公式成立：
Figure BDA0004107600070000011
式中，Qn是指定截面法向上的热流，单位W；K是比例常数，即特定材料的热导率，单位W/m/℃；A是指定截面的面积，单位m2；dT/dn是沿指定截面法向的温度变化率，即温度梯度，单位℃/m。该式即为傅立叶定律的微分表达形式。
以图1(a)所示的平板传热模型为例，设平板为有限大，其厚度均匀、内部连续、材质均一，前侧表面保持恒定温度Tconst(℃)，右侧表面暴露于Tatm(℃)恒温空气中、换热系数h(W/m2/℃)。
图中未标注边界条件的外表面皆保持绝热，在指定边界条件的持续作用下，具有任意初始温度分布的平板将逐渐改变其温度分布，足够长时间后其温度分布不再随时间变化，此时板内通过的热流保持恒定，称这一状态为稳态传热，称此时的温度分布为稳态温度场。
当工艺参数稳定，铜管水平连铸的理想状态就是稳态传热，即连铸结晶器内温度分布不随时间改变。设图1(a)中平板厚度足够小，如此能将其稳态温度场的求解简化为二维问题，简化后模型如图1(b)所示。
将图1(b)所示二维模型划分为横纵密度相等的网格，如图2所示；若能求得全部网格节点处的稳态温度值，也就近似求得了稳态温度场；当网格密度足够大时，即便在局部节点之间进行简单线性插值，也能足够精确地得出二维传热区域内任意位置的稳态温度值。按照节点所处位置在传热区域内部抑或边界、所处边界处的边界条件类型，能将全部节点分为6类；在图2中，按照从左到右、从下到上的顺序对节点进行编号，选择编号为4、22、25、28、43、49的6个节点作为上述6类节点的典型。接下来分别讨论6个典型节点处的热量流动情况，并尝试建立起典型节点与相邻节点之间的关系，进而建立起全部节点之间的关系，最终以已知的边界条件为线索求得各未知节点温度。
图2中，节点4所处位置为传热区域的下边界，该边界已设置恒温边界条件，故有：
T4＝Tconst (2)
式中T4为节点4处的温度。本申请文件中，无特殊说明时数字下标皆表示节点编号。
图2中，节点(22)所处位置为传热区域的左边界，该边界已设置绝热边界条件，故节点(22)只受编号为15、23、29的3个相邻节点影响。为考察节点(22)处热量流动情况，在其周围建立起虚拟的控制体积，控制体积由节点(22)所属的两个网眼的各1/4构成，其边界以虚线表示。根据式(1)，单位时间从节点(22)控制体积的右边界流入控制体积的热量为：
Figure BDA0004107600070000031
式中微分项的下角标e表示右侧，即east；类似的，下文中出现的下角标w、n、s分别表示左侧、上侧、下侧。假设节点间温度线性分布，将式中微分项用差商代替能得：
Figure BDA0004107600070000032
式中K为平板材料的热导率；Δ为网眼边长；A为节点(22)控制体积的右边界面积，A＝Δ·Th，其中，Th为平板厚度。同理能得单位时间从节点(22)控制体积的上边界和下边界流入控制体积的热量分别为：
Figure BDA0004107600070000033
Figure BDA0004107600070000034
由于稳态传热时节点(22)处温度保持不变，由此能知单位时间内节点(22)控制体积的净流入热量为0，故有如下热平衡方程成立：
Q22n+Q22s+Q22e＝0 (7)
整理能得：
T15+2T23+T29-4T22＝0 (8)
图2中，节点(25)所处位置为传热区域的内部，故节点(25)只受编号为18、26、32、24的4个相邻节点影响。根据“单位时间内节点控制体积净流入热量为0”的结论，能得节点(25)与4个相邻节点的温度关系为：
T18+T26+T32+T24-4T25＝0 (9)
图2中，节点(28)所处位置为传热区域的右边界，该边界已设置对流换热边界条件，故节点(28)除了受编号为27、21、35的3个相邻节点影响，还受到对流换热边界条件影响。根据牛顿冷却定律，单位时间从节点(28)控制体积的右边界流入控制体积的热量为：
Q28e＝hA(Tatm-T28) (10)
根据式(1)，有：
Figure BDA0004107600070000041
Figure BDA0004107600070000042
Figure BDA0004107600070000043
根据“单位时间内节点控制体积净流入热量为0”的结论，有：
Q28w+Q28s+Q28e+Q28n＝0 (14)
整理能得：
P27T27+P21T21+P35T35+P28T28+C28＝0 (15)
设材料热导率、对流换热系数为常数，则式中P和C为常数。
图2中，节点(43)所处位置为传热区域的左上角点，其所属的两个边界均已设置绝热边界条件，故节点(43)只受编号为36、44的2个相邻节点影响。根据“单位时间内节点控制体积净流入热量为0”的结论，能得节点(43)与2个相邻节点的温度关系为：
T36+T44-2T43＝0 (16)
图2中，节点(49)所处位置为传热区域的右上角点，其所属的传热区域上边界已设置绝热边界条件，其所属的传热区域右边界已设置对流换热边界条件，基于上文结论不难得出：
Figure BDA0004107600070000044
Figure BDA0004107600070000051
Figure BDA0004107600070000052
Q49w+Q49s+Q49e＝0 (20)
整理能得：
P48T48+P42T42+P49T49+C49＝0 (21)
至此，图2中6个典型节点与各自相邻节点的温度关系已全部找到，由此易知全部节点与各自相邻节点的温度关系。分析6个典型节点与各自相邻节点的温度关系，发现其数学形式都属于线性方程；由于在物理范畴这些关系同时成立，故能在数学范畴将其联立，得到能表示全部节点之间温度关系的数学形式，即线性方程组；该线性方程组具有唯一解，该解即是图1中传热模型的稳态温度场。
几何模型建立及离散化：
铜管水平连铸涉及的关键机械结构如图2.1所示，能见一冷水中发挥主要换热作用的，是与铜套相接触的筒形水冷区域，为提高建模效率能忽略其他部分；石墨模进液口以左的部分浸没在浇铸温度的铜液中，内部温度均匀等于铜液浇铸温度，为提高建模效率能忽略这一部分。综上，图2.1中结构能简化为轴对称结构、从3维降为2维，如图3所示。
对简化后的几何模型进行离散处理，如图4所示，图中点阵即同图2中的网格节点。图3中两个筒形水冷区域厚度方向的温度分布不是主要研究对象，故能将两个筒形水冷区域离散成一维点阵，如图4中“一冷水”、“二冷水”两个点阵区域所示。由图3能见，石墨芯棒带有一定锥度，故图4中“铜液/铸坯”以及“石墨芯棒”两个区域的点阵是非正交排列的。鉴于石墨筒外部的铸坯区域温度梯度较小，为其设置了较小的轴向离散密度。为便于界面换热的建模，令图4中相邻点阵在分界处相互重合。
材料属性设定与实现：
图3所示简化模型包含6个传热区域，其中铜液/铸坯的材料为TP2磷脱氧铜；铜套的材料为某牌号的镍铜，其热物理属性与TP2铜相近，能认为铜套的材料也为TP2铜；此外，石墨筒和石墨芯棒的材料为同一牌号石墨、一冷水和二冷水的材料都能认为是纯水，故当前模型共涉及3种材料，分别是石墨、TP2铜、纯水，其热物理属性的设定参考所述采用商业有限元软件建立的铜管水平连铸流动及传热过程的离线数值模型对应设置即能。
将密度描述为温度的函数，其在热平衡方程中的实现比较简单，将在本章第6节讨论；其他温度相关热物理属性的实现方法与之相似，但在当前模型中热导率的温度相关实现却需要特殊关注。对于图5所示的一维传热模型，将其材料热导率描述为温度的函数，图中w、e分别位于节点1、2及节点(2)、3的中点，w、e之间的材料区域即为节点(2)的控制体积。基于本章第1节所述原理，对于节点(2)有如下热平衡方程成立：
Figure BDA0004107600070000061
式中Kw、Ke分别为w、e处的材料热导率，其中Ke往往被简单取值为K2与K3的算术平均，但这样的简单取值并不准确。
基于热导率在相邻控制体积之间发生阶跃变化的假设，能知图5中节点(2)和e之间的材料热导率为K2，e和节点3之间的材料热导率为K3，从而能将e处热流密度qe表示为：
Figure BDA0004107600070000062
式中最右侧分式的分母部分表示节点(2)、3之间的总热阻。由上式不难得出：
Figure BDA0004107600070000063
显然，从准确表示e处热流密度的角度出发，Ke的取值应该是K2与K3的调和平均，而非算术平均。虽然这一结论是基于热导率阶跃变化假设得出的，然而实践表明，即便在热导率连续变化情况下，调和平均取值方式也要比算术平均更有优势。
边界条件设定与实现：
图3中1处位于石墨筒进液端，相邻表面完全浸没于浇铸温度的铜液；为简化模型，于该处设置恒温边界条件近似取代铜液的加热作用，温度值设定为铜液浇铸温度即能。该边界条件在热平衡方程中的实现参考式(2)。
图3中2处为石墨筒与炉壁的接触面，热量经该接触面由石墨筒流向炉壁；为简化模型，于该处设置对流换热边界条件近似取代炉壁的冷却作用；流体温度参考炉壁温度实测值设定为1050℃，换热系数按第1节参数校核结果设定为3000W/m/℃。该边界条件在热平衡方程中的实现参考式(10)。
图3中3、4处位于一冷区和二冷区之间，表面温度较高，必须考虑辐射换热；该处也有氮气、水汽通过，对流换热同样不容忽视，故能于该处设置“对流换热+辐射换热”混合边界条件。其中对流换热边界条件的实现参考式(10)，辐射换热边界条件依据斯特藩-玻尔兹曼定律(Stefan-Boltzmann law)实现：
Q＝ε1A1σF12(T1 4-T2 4) (25)
式中Q为表面1与表面2之间的辐射传热功率；ε1为表面1的发射率；A1为表面1的面积；ζ为斯特藩-玻尔兹曼常数；F12为表面1到表面2的辐射角系数；T1、T2分别为表面1、表面2的温度。当前边界条件中，对流换热的流体温度按氮气吹入温度和二冷来水温度的算术平均设定，换热系数通过实测温度进行校准；由于换热表面成对平行且相互贴近，辐射换热的角系数能设定为1；式(25)中T1对应图3中3、4处的温度；T2对应着周围零件表面温度，为提高计算效率未将这些零件纳入模型，故按实测温度将T2近似设定为常数。
图3中5、6处位于石墨芯棒末端、暴露于空气之中，然而该处空间被铜管内壁遮挡、几乎完全封闭，故能忽略对流换热而只于该处设置辐射换热边界条件，辐射换热边界条件的实现参考上文。
界面属性设定与实现：
图6所示的一维传热模型包含W、E两个传热区域，W区域的右端点(节点w1处)与E区域的左端点(节点e1处)之间有热导率为K、厚度为δ的界面薄层。当δ→0或K＞＞max(KW,KE)时，界面薄层法向热阻δ/K能忽略不计，
Figure BDA0004107600070000071
而一般情况下，界面热阻显著影响连铸传热过程，故此有必要研究界面热阻在热平衡方程中的实现。
对于图6中传热模型，考察节点w1、e1之间的热流，根据式(1)有：
Figure BDA0004107600070000081
式中Q为节点w1、e1之间的热流，A为界面面积。令h＝K/δ，则式(26)变为：
Figure BDA0004107600070000082
能见式(27)与式(10)形式相同，故在建立节点w1处的热平衡方程时，可将其右侧的模型部分等效为对流换热边界条件，流体温度设定为
Figure BDA0004107600070000083
换热系数设定为K/δ。同理，在建立节点e1处的热平衡方程时，可将其左侧的模型部分等效为对流换热边界条件，流体温度设定为
Figure BDA0004107600070000084
换热系数设定为K/δ。虽然出现在“边界条件”中的
Figure BDA0004107600070000085
未知，但基于模型中其他节点处的已知边界条件，仍可通过联立求解两传热区域共m+n个热平衡方程得到所有未知温度的值。
运动传热模型建立：
铜管水平连铸过程中的冷却水流动、铜液流动、铸坯运动在本质上影响着稳态温度场，一般的稳态传热问题并不涉及这些运动因素，经典传热理论亦无专门介绍，故有必要对其进行研究，并建立起相应的数学模型。
图7所示为铸坯传热区域的一个局部，因铸坯被匀速向右牵出，故在单位时间内总有一定质量的铸坯从左向右通过节点(2)控制体积。稳态传热时，节点(2)控制体积左右边界之间存在一定温差，由此推知铸坯在通过控制体积过程中必然发生了相应程度的热量损失，将这部分热量表示为Qtransport，则：
Qtransport＝ρvACp(T2e-T2w) (28)
式中ρ为铸坯的密度，将密度描述为温度的函数时，能取ρ的值为ρ(T2)；v为铸坯牵引速度；A为节点(2)控制体积左右边界的面积；Cp为铸坯的比热容；T2w、T2e分别为节点(2)控制体积左右边界处的温度。
式(28)中，ρvA即等于单位时间内通过节点(2)控制体积的铸坯总质量，“质量×比热容×温度变化量”即为该质量的铸坯改变一定温度所涉及的热量变化。基于节点间温度线性分布假设，式(28)有如下等价形式：
Figure BDA0004107600070000091
将式(29)加入节点(2)的热平衡方程即在数学关系层面实现了铸坯运动对传热的影响；若将相应的流体流动近似为沿主流方向的刚体平移，也能采用同样的方式描述冷却水流动及铜液流动对传热的影响。
相变潜热模型建立：铜管水平连铸过程中相变潜热的释放同样在本质上影响着稳态温度场。为便于后续讨论，首先阐明固相率的一些相关概念。金属凝固过程中其固液两相之间存在一个分界面(固液界面)，纯金属在温度降至其熔点以下后即完全凝固，故其固液界面比较分明；而一般的合金在温度降至其熔点(液相线)以下后，仍需继续降至更低温度(固相线)才能完全凝固，故其固液两相之间存在一定宽度的过渡区域，该区域内金属介于固液两相之间，故称该区域为糊状区。假设金属在凝固过程中其固相体积分数(即固相率fs，0≤fs≤1)与温度成线性关系，则能按下式计算fs：
Figure BDA0004107600070000092
图8所示为铜液传热区域的一个局部，因铸坯被匀速向右牵出，故在单位时间内总有一定质量的铜液从左向右流经节点(2)控制体积。稳态传热时，节点(2)控制体积左右边界之间存在一定温差，故也存在一定固相率之差，由此推知铜液在流经控制体积过程中必然发生了相应程度的相变潜热释放，将这部分热量表示为Qlatent，则：
Figure BDA0004107600070000093
式中ρ为铸坯的密度；v为铸坯牵引速度；A为节点(2)控制体积左右边界的面积；L为铜液的结晶潜热；fs2w、fs2e分别为节点(2)控制体积左右边界处的固相率。
式(31)中，ρvA即等于单位时间内流经节点(2)控制体积的铜液总质量，“质量×结晶潜热×凝固分数”即为该质量的铜液凝固一定分数所涉及的相变潜热释放。基于节点间温度线性分布假设、固相率与温度成线性关系假设，式(31)有如下等价形式：
Figure BDA0004107600070000101
将式(32)加入节点(2)的热平衡方程即在数学关系层面实现了相变潜热释放对传热的影响。值得注意的，当式(32)中v取值为0时Qlatent等于0，然而这并不意味着实际凝固过程中无相变潜热释放。在铸坯牵引速度为0、边界条件保持不变的情况下：非稳态传热时，固液界面仍未停止向液相中推进，此时相变仍继续发生，从而有相变潜热释放；稳态传热时，固液界面停止向液相中推进，此时能认为无相变发生，从而无相变潜热释放。由此能见，式(32)中v取值为0时Qlatent等于0是有物理意义的，正确反映了特殊情况下的潜热释放情况。
热平衡方程组求解：
对于图1所示的平板传热模型，其所涉及的热平衡方程均为线性方程，故其热平衡方程组的求解采用一般的线性方程组解法即能。然而对于铜管连铸模型，其所涉及的热平衡方程不全是线性的，其中部分方程的系数与未知的节点温度有关，也有部分方程包含非线性的源项；在此情况下，迭代法是热平衡方程组的唯一能行解法。
采用迭代法求解非线性热平衡方程组的过程包括如下几步：开始求解时，对所有节点估计一组T值；将估计的T值代入热平衡方程系数表达式，得到名义上的线性方程组；解名义上的线性方程组，得到新一组T值；以新得到的一组T值作为相较更好的估计值，重复执行2、3、4步，直到这种进一步的重复计算(称之为迭代)不再引起T值显著变化时为止。
这种最终T值不再显著变化的状态称之为迭代的收敛，尽管这个解是由线性方程组求解得到的，但它实际上就是非线性热平衡方程组的正确解。然而也能能经过迭代不会收敛到一个解，T值稳定地飘移或是以一个不断增大的幅度震荡，这一与收敛相对立的过程称之为发散。对于铜管连铸模型，遵循前文所述的建模原理就基本能够保证收敛；为了消除类似于辐射边界条件的非线性元素对迭代收敛的不利影响，能针对这些非线性元素应用一些线性化及松弛的方法。
在所有的迭代法当中，最简单的一种是高斯-赛德尔(Gauss-Seidel)逐点计算法。该法首先要为所有节点估计一组T值，然后开始按一定的顺序逐个更新每个节点的T值。T值的更新是依据式(9)那样的代数关系，即当前节点处的热平衡方程；不同的是，此处方程中所引用的相邻节点的T值必须是其最新值，也即T的初始估计值或最近更新值；重复多轮更新所有节点的T值，即能逐渐获得方程组的收敛解。高斯-赛德尔法的最大缺点是其收敛速度太慢，当模型节点数量较大时尤为明显；由于本方法是以一次迭代前进一个节点间距的速率来传递边界条件信息到模型内部的，故其收敛速度慢是能理解的。
为了克服高斯-赛德尔法收敛速度慢的劣势，将其单次迭代更新一个节点改为更新一行节点。如图9所示，将图中圆点标识处作为当前节点行，以相邻节点行的最新T值作为已知条件，建立当前行各个节点处的热平衡方程并进行联立求解，将所得新一组T值作为当前行的最新T值，如此便完成了一次迭代。显然，上述迭代步中要进行求解的方程组是三对角型方程组，该类型方程组能采用TDMA追赶法[71]进行求解，算法实现比较简单。对于图4所示的传热模型，能采用上述的逐行更新方式在横纵两个方向交替地扫描，如此即能将所有的边界条件信息传递到模型内部，最终得到收敛解。逐行更新的高斯-赛德尔迭代的收敛速度更快，这是因为不管沿行有多少个节点，在行两端的边界条件信息都能立即传递到模型内部。
在线计算程序开发：
基于前文工作开发了铜管水平连铸稳态温度场在线计算程序，程序结构如图10所示。在“系统参数加载”环节：(1)程序访问系统盘根目录下的“root.txt”文件，从中读出用户事先设定的工作目录所在位置。工作目录用于存储从上位机传来的实时工艺参数、程序运行过程中产生的临时文件以及计算结果文件。(2)程序访问工作目录下的“inter.txt”文件，从中读出用户事先设定的工艺参数采集节拍。工艺参数采集节拍包含两个信息：采集次数N和采集间歇Δt。从上位机传来的数据更新频率较快且存在一定扰动，为减小数据扰动，每次计算前采集N组工艺参数，相邻两次采集间隔Δt秒，最后最将N组工艺参数的算术平均作为本次计算的输入，这就是程序在“工艺参数加载”环节执行的操作。
在“模型参数加载”环节，程序访问工作目录下的“modata.txt”文件，从中读出用户事先设定的模型几何尺寸、几何离散参数、材料热物理属性、边界条件及界面换热系数。在“几何离散化”环节，程序根据用户设定的模型几何尺寸和几何离散参数在各传热区域划分网格，并按从1开始依次递增1的正整数序列对网格节点进行编号。为便于后续引用，将节点编号按节点在其所属传热区域内的空间位置存入一个二维数组，数组的行列数等于当前传热区域内网格节点的行列数；为每个传热区域建立一个这样的数组，如此在后续建模过程中需要对某一节点进行操作时，就能方便地通过其所属的传热区域以及其在所属传热区域内的位置引用其编号。
基于前文所述建模原理，能够建立起模型中所有典型节点处的热平衡方程。在“热平衡系数表构建”环节，程序根据典型节点热平衡方程的系数形式计算出所有同类型节点处热平衡方程的系数，并将系数向量按节点在其所属传热区域内的空间位置存入一个二维数组，数组的行列数等于当前传热区域内网格节点的行列数；为每个传热区域建立一个这样的数组，如此在后续的方程组迭代求解过程中能直接引用这些系数而不必重复计算，能极大提高程序的时间效率。
在“定义固相率”、“定义相变潜热”及“定义运动传热”环节，将数学模型定义为相应的函数形式，如此在构建热平衡方程组时即能方便地引用这些函数，而不必为同样的计算过程重复编写代码，从而极大提高程序的空间效率。在“热平衡方程组构建”环节，程序基于之前四个环节的准备工作建立起整个模型的热平衡方程组。
在“方程组解初始化”环节，程序为所有节点估计一组T值；以这组T值为起点，程序在“方程组求解”环节采用迭代法求解已经构建好的热平衡方程组。由图10能见，当前程序在整体上是一个无限循环结构，第一次循环中，程序在“方程组解初始化”环节为同一个传热区域内的节点估计相同的T值，不同传热区域的T值按经验大致设定；后续循环中，程序将上次循环得到的热平衡方程组的收敛解作为“方程组解初始化”环节T的估计值，由于在实际生产中前后两次计算对应的工艺参数一般相差不大，能预计前后两次计算的收敛解也相差不大，故这种“解初始化更新”的技巧能减少后续循环中获得收敛解所需的迭代次数，从而提高程序的时间效率。
在“控制台输出”环节，程序将本次计算开始的系统时间、所采用的工艺参数、计算收敛情况、获得收敛解所经过的迭代次数以及壁钟时间等信息输出到控制台窗口，从而实时地向操作人员显示程序运行状态；此外，相同的信息在“系统日志更新”环节被追加写入工作目录下的“log.txt”文件，以备程序运行意外中断时对故障进行诊断。
如图11、12所示，在“绘制温度场云图”、“绘制固相率云图”环节，程序分别将温度场计算收敛解以及对应的固相率分布绘制为彩色图像，不同颜色对应场量值的不同水平，以便操作人员更直观地查看计算结果。为使温度场计算结果的定量查看更加便捷，在温度场计算结果彩图中添加了带有温度值的等温线。
为验证本方案模型计算结果准确性，将其稳态温度计算结果与采用商业有限元软件建立的铜管水平连铸流动及传热过程的离线数值模型进行对比，温度对比的位置选定在石墨筒厚度中面，对比结果如图13所示。能见当前模型的温度计算结果与采用商业有限元软件建立的铜管水平连铸流动及传热过程的离线数值模型高度吻合，说明当前模型的传热计算是可靠的。
本发明与现有技术相比，其优点在于：
本发明所述的铜管水平连铸温度场在线计算方法，围绕铜管水平连铸稳态传热的有限差分法建模，讨论了几何模型的建立及离散化、材料属性、边界条件、界面属性的设定与实现等问题。为描述连铸过程中相对运动及相变潜热对于传热的影响，通过物质状态变化的微元分析建立了相应数学模型。
鉴于铜管水平连铸模型涉及非线性热平衡方程，提出了一种非线方程组的迭代解法。开发了水平连铸稳态温度场在线计算程序，并将其计算结果与采用商业有限元软件建立的铜管水平连铸流动及传热过程的离线数值模型稳态温度场计算结果进行了对比，验证了本章在线计算程序的传热计算可靠性。
附图说明
下面结合附图及实施方式对本发明作进一步详细的说明：
图1平板传热模型；
图2平板传热二维离散模型；
1-6：边界条件；7：固液界面。A：一冷水，B：铜套，C：石墨筒，
D：二冷水，E：铜液，F：铸坯，G：石墨芯棒。
图3铜管水平连铸轴对称简化模型；
图4铜管水平连铸二维离散模型；
图5材料热导率随温度变化的一维传热模型；
图6含接触换热的一维传热模型；
图7运动传热模型；
图8相变潜热模型；
图9逐行更新的高斯-赛德尔迭代；
图10铜管水平连铸稳态温度场在线计算程序框图；
图11在线程序温度场计算结果；
图12在线程序固相率分布计算结果；
图13本方案模型与采用商业有限元软件建立的铜管水平连铸流动及传热过程的离线数值模型稳态温度计算结果对比。
具体实施方式
下面将结合具体的实施方案对本发明进行进一步的解释，但并不局限本发明，说明书附图所绘示的结构、比例、大小等，均仅用以配合说明书所揭示的内容，以供熟悉此技术的人士了解与阅读，并非用以限定本发明能实施的限定条件，故不具技术上的实质意义，任何结构的修饰、比例关系的改变或大小的调整，在不影响本发明所能产生的功效及所能达成的目的下，均应仍落在本发明所揭示的技术内容所能涵盖的范围内。同时，本说明书中所引用的如“上”、“下”、“前”、“后”、“中间”等用语，亦仅为便于叙述的明了，而非用以限定本发明能实施的范围，其相对关系的改变或调整，在无实质变更技术内容下，当亦视为本发明能实施的范畴。
请参阅图
本发明提供了一种铜管水平连铸温度场在线计算方法，具体包括以下内容：
连铸温度场在线计算：
建模原理
在实验研究的基础上，傅立叶于1822年建立了“热的解析理论”，奠定了热传导的理论基础[70]。实验观察发现，单位时间内通过指定截面的导热量，正比于沿该截面法向的温度变化率和截面面积，从而有如下经验公式成立：
Figure BDA0004107600070000141
式中，Qn是指定截面法向上的热流，单位W；K是比例常数，即特定材料的热导率，单位W/m/℃；A是指定截面的面积，单位m2；dT/dn是沿指定截面法向的温度变化率，即温度梯度，单位℃/m。该式即为傅立叶定律的微分表达形式。
以图1(a)所示的平板传热模型为例，设平板为有限大，其厚度均匀、内部连续、材质均一，前侧表面保持恒定温度Tconst(℃)，右侧表面暴露于Tatm(℃)恒温空气中、换热系数h(W/m2/℃)。
图中未标注边界条件的外表面皆保持绝热，在指定边界条件的持续作用下，具有任意初始温度分布的平板将逐渐改变其温度分布，足够长时间后其温度分布不再随时间变化，此时板内通过的热流保持恒定，称这一状态为稳态传热，称此时的温度分布为稳态温度场。
当工艺参数稳定，铜管水平连铸的理想状态就是稳态传热，即连铸结晶器内温度分布不随时间改变。设图1(a)中平板厚度足够小，如此能将其稳态温度场的求解简化为二维问题，简化后模型如图1(b)所示。
将图1(b)所示二维模型划分为横纵密度相等的网格，如图2所示；若能求得全部网格节点处的稳态温度值，也就近似求得了稳态温度场；当网格密度足够大时，即便在局部节点之间进行简单线性插值，也能足够精确地得出二维传热区域内任意位置的稳态温度值。按照节点所处位置在传热区域内部抑或边界、所处边界处的边界条件类型，能将全部节点分为6类；在图2中，按照从左到右、从下到上的顺序对节点进行编号，选择编号为4、22、25、28、43、49的6个节点作为上述6类节点的典型。接下来分别讨论6个典型节点处的热量流动情况，并尝试建立起典型节点与相邻节点之间的关系，进而建立起全部节点之间的关系，最终以已知的边界条件为线索求得各未知节点温度。
图2中，节点4所处位置为传热区域的下边界，该边界已设置恒温边界条件，故有：
T4＝Tconst (2)
式中T4为节点4处的温度。本申请文件中，无特殊说明时数字下标皆表示节点编号。
图2中，节点(22)所处位置为传热区域的左边界，该边界已设置绝热边界条件，故节点(22)只受编号为15、23、29的3个相邻节点影响。为考察节点(22)处热量流动情况，在其周围建立起虚拟的控制体积，控制体积由节点(22)所属的两个网眼的各1/4构成，其边界以虚线表示。根据式(1)，单位时间从节点(22)控制体积的右边界流入控制体积的热量为：
Figure BDA0004107600070000161
式中微分项的下角标e表示右侧，即east；类似的，下文中出现的下角标w、n、s分别表示左侧、上侧、下侧。假设节点间温度线性分布，将式中微分项用差商代替能得：
Figure BDA0004107600070000162
式中K为平板材料的热导率；Δ为网眼边长；A为节点(22)控制体积的右边界面积，A＝Δ·Th，其中，Th为平板厚度。同理能得单位时间从节点(22)控制体积的上边界和下边界流入控制体积的热量分别为：
Figure BDA0004107600070000163
Figure BDA0004107600070000164
由于稳态传热时节点(22)处温度保持不变，由此能知单位时间内节点(22)控制体积的净流入热量为0，故有如下热平衡方程成立：
Q22n+Q22s+Q22e＝0 (7)
整理能得：
T15+2T23+T29-4T22＝0 (8)
图2中，节点(25)所处位置为传热区域的内部，故节点(25)只受编号为18、26、32、24的4个相邻节点影响。根据“单位时间内节点控制体积净流入热量为0”的结论，能得节点(25)与4个相邻节点的温度关系为：
T18+T26+T32+T24-4T25＝0 (9)
图2中，节点(28)所处位置为传热区域的右边界，该边界已设置对流换热边界条件，故节点(28)除了受编号为27、21、35的3个相邻节点影响，还受到对流换热边界条件影响。根据牛顿冷却定律，单位时间从节点(28)控制体积的右边界流入控制体积的热量为：
Q28e＝hA(Tatm-T28) (10)
根据式(1)，有：
Figure BDA0004107600070000171
Figure BDA0004107600070000172
Figure BDA0004107600070000173
根据“单位时间内节点控制体积净流入热量为0”的结论，有：
Q28w+Q28s+Q28e+Q28n＝0 (14)
整理能得：
P27T27+P21T21+P35T35+P28T28+C28＝0 (15)
设材料热导率、对流换热系数为常数，则式中P和C为常数。
图2中，节点(43)所处位置为传热区域的左上角点，其所属的两个边界均已设置绝热边界条件，故节点(43)只受编号为36、44的2个相邻节点影响。根据“单位时间内节点控制体积净流入热量为0”的结论，能得节点(43)与2个相邻节点的温度关系为：
T36+T44-2T43＝0 (16)
图2中，节点(49)所处位置为传热区域的右上角点，其所属的传热区域上边界已设置绝热边界条件，其所属的传热区域右边界已设置对流换热边界条件，基于上文结论不难得出：
Figure BDA0004107600070000174
Figure BDA0004107600070000181
Figure BDA0004107600070000182
Q49w+Q49s+Q49e＝0 (20)
整理能得：
P48T48+P42T42+P49T49+C49＝0 (21)
至此，图2中6个典型节点与各自相邻节点的温度关系已全部找到，由此易知全部节点与各自相邻节点的温度关系。分析6个典型节点与各自相邻节点的温度关系，发现其数学形式都属于线性方程；由于在物理范畴这些关系同时成立，故能在数学范畴将其联立，得到能表示全部节点之间温度关系的数学形式，即线性方程组；该线性方程组具有唯一解，该解即是图1中传热模型的稳态温度场。
几何模型建立及离散化：
铜管水平连铸涉及的关键机械结构如图2.1所示，能见一冷水中发挥主要换热作用的，是与铜套相接触的筒形水冷区域，为提高建模效率能忽略其他部分；石墨模进液口以左的部分浸没在浇铸温度的铜液中，内部温度均匀等于铜液浇铸温度，为提高建模效率能忽略这一部分。综上，图2.1中结构能简化为轴对称结构、从3维降为2维，如图3所示。
对简化后的几何模型进行离散处理，如图4所示，图中点阵即同图2中的网格节点。图3中两个筒形水冷区域厚度方向的温度分布不是主要研究对象，故能将两个筒形水冷区域离散成一维点阵，如图4中“一冷水”、“二冷水”两个点阵区域所示。由图3能见，石墨芯棒带有一定锥度，故图4中“铜液/铸坯”以及“石墨芯棒”两个区域的点阵是非正交排列的。鉴于石墨筒外部的铸坯区域温度梯度较小，为其设置了较小的轴向离散密度。为便于界面换热的建模，令图4中相邻点阵在分界处相互重合。
材料属性设定与实现：
图3所示简化模型包含6个传热区域，其中铜液/铸坯的材料为TP2磷脱氧铜；铜套的材料为某牌号的镍铜，其热物理属性与TP2铜相近，能认为铜套的材料也为TP2铜；此外，石墨筒和石墨芯棒的材料为同一牌号石墨、一冷水和二冷水的材料都能认为是纯水，故当前模型共涉及3种材料，分别是石墨、TP2铜、纯水，其热物理属性的设定参考所述采用商业有限元软件建立的铜管水平连铸流动及传热过程的离线数值模型对应设置即能。
将密度描述为温度的函数，其在热平衡方程中的实现比较简单，将在本章第6节讨论；其他温度相关热物理属性的实现方法与之相似，但在当前模型中热导率的温度相关实现却需要特殊关注。对于图5所示的一维传热模型，将其材料热导率描述为温度的函数，图中w、e分别位于节点1、2及节点(2)、3的中点，w、e之间的材料区域即为节点(2)的控制体积。基于本章第1节所述原理，对于节点(2)有如下热平衡方程成立：
Figure BDA0004107600070000191
式中Kw、Ke分别为w、e处的材料热导率，其中Ke往往被简单取值为K2与K3的算术平均，但这样的简单取值并不准确。
基于热导率在相邻控制体积之间发生阶跃变化的假设，能知图5中节点(2)和e之间的材料热导率为K2，e和节点3之间的材料热导率为K3，从而能将e处热流密度qe表示为：
Figure BDA0004107600070000192
式中最右侧分式的分母部分表示节点(2)、3之间的总热阻。由上式不难得出：
Figure BDA0004107600070000193
显然，从准确表示e处热流密度的角度出发，Ke的取值应该是K2与K3的调和平均，而非算术平均。虽然这一结论是基于热导率阶跃变化假设得出的，然而实践表明，即便在热导率连续变化情况下，调和平均取值方式也要比算术平均更有优势。
边界条件设定与实现：
图3中1处位于石墨筒进液端，相邻表面完全浸没于浇铸温度的铜液；为简化模型，于该处设置恒温边界条件近似取代铜液的加热作用，温度值设定为铜液浇铸温度即能。该边界条件在热平衡方程中的实现参考式(2)。
图3中2处为石墨筒与炉壁的接触面，热量经该接触面由石墨筒流向炉壁；为简化模型，于该处设置对流换热边界条件近似取代炉壁的冷却作用；流体温度参考炉壁温度实测值设定为1050℃，换热系数按第1节参数校核结果设定为3000W/m/℃。该边界条件在热平衡方程中的实现参考式(10)。
图3中3、4处位于一冷区和二冷区之间，表面温度较高，必须考虑辐射换热；该处也有氮气、水汽通过，对流换热同样不容忽视，故能于该处设置“对流换热+辐射换热”混合边界条件。其中对流换热边界条件的实现参考式(10)，辐射换热边界条件依据斯特藩-玻尔兹曼定律(Stefan-Boltzmann law)实现：
Q＝ε1A1σF12(T1 4-T2 4) (25)
式中Q为表面1与表面2之间的辐射传热功率；ε1为表面1的发射率；A1为表面1的面积；ζ为斯特藩-玻尔兹曼常数；F12为表面1到表面2的辐射角系数；T1、T2分别为表面1、表面2的温度。当前边界条件中，对流换热的流体温度按氮气吹入温度和二冷来水温度的算术平均设定，换热系数通过实测温度进行校准；由于换热表面成对平行且相互贴近，辐射换热的角系数能设定为1；式(25)中T1对应图3中3、4处的温度；T2对应着周围零件表面温度，为提高计算效率未将这些零件纳入模型，故按实测温度将T2近似设定为常数。
图3中5、6处位于石墨芯棒末端、暴露于空气之中，然而该处空间被铜管内壁遮挡、几乎完全封闭，故能忽略对流换热而只于该处设置辐射换热边界条件，辐射换热边界条件的实现参考上文。
界面属性设定与实现：
图6所示的一维传热模型包含W、E两个传热区域，W区域的右端点(节点w1处)与E区域的左端点(节点e1处)之间有热导率为K、厚度为δ的界面薄层。
对于图6中传热模型，考察节点w1、e1之间的热流，根据式(1)有：
Figure BDA0004107600070000211
式中Q为节点w1、e1之间的热流，A为界面面积。令h＝K/δ，则式(26)变为：
Figure BDA0004107600070000212
运动传热模型建立：
铜管水平连铸过程中的冷却水流动、铜液流动、铸坯运动在本质上影响着稳态温度场，一般的稳态传热问题并不涉及这些运动因素，经典传热理论亦无专门介绍，故有必要对其进行研究，并建立起相应的数学模型。
图7所示为铸坯传热区域的一个局部，因铸坯被匀速向右牵出，故在单位时间内总有一定质量的铸坯从左向右通过节点(2)控制体积。稳态传热时，节点(2)控制体积左右边界之间存在一定温差，由此推知铸坯在通过控制体积过程中必然发生了相应程度的热量损失，将这部分热量表示为Qtransport，则：
Qtransport＝ρvACp(T2e-T2w) (28)
式中ρ为铸坯的密度，将密度描述为温度的函数时，能取ρ的值为ρ(T2)；v为铸坯牵引速度；A为节点(2)控制体积左右边界的面积；Cp为铸坯的比热容；T2w、T2e分别为节点(2)控制体积左右边界处的温度。
式(28)中，ρvA即等于单位时间内通过节点(2)控制体积的铸坯总质量，“质量×比热容×温度变化量”即为该质量的铸坯改变一定温度所涉及的热量变化。基于节点间温度线性分布假设，式(28)有如下等价形式：
Figure BDA0004107600070000213
将式(29)加入节点(2)的热平衡方程即在数学关系层面实现了铸坯运动对传热的影响；若将相应的流体流动近似为沿主流方向的刚体平移，也能采用同样的方式描述冷却水流动及铜液流动对传热的影响。
相变潜热模型建立：铜管水平连铸过程中相变潜热的释放同样在本质上影响着稳态温度场。为便于后续讨论，首先阐明固相率的一些相关概念。金属凝固过程中其固液两相之间存在一个分界面(固液界面)，纯金属在温度降至其熔点以下后即完全凝固，故其固液界面比较分明；而一般的合金在温度降至其熔点(液相线)以下后，仍需继续降至更低温度(固相线)才能完全凝固，故其固液两相之间存在一定宽度的过渡区域，该区域内金属介于固液两相之间，故称该区域为糊状区。假设金属在凝固过程中其固相体积分数(即固相率fs，0≤fs≤1)与温度成线性关系，则能按下式计算fs：
Figure BDA0004107600070000221
图8所示为铜液传热区域的一个局部，因铸坯被匀速向右牵出，故在单位时间内总有一定质量的铜液从左向右流经节点(2)控制体积。稳态传热时，节点(2)控制体积左右边界之间存在一定温差，故也存在一定固相率之差，由此推知铜液在流经控制体积过程中必然发生了相应程度的相变潜热释放，将这部分热量表示为Qlatent，则：
Figure BDA0004107600070000222
式中ρ为铸坯的密度；v为铸坯牵引速度；A为节点(2)控制体积左右边界的面积；L为铜液的结晶潜热；
Figure BDA0004107600070000223
分别为节点(2)控制体积左右边界处的固相率。
式(31)中，ρvA即等于单位时间内流经节点(2)控制体积的铜液总质量，“质量×结晶潜热×凝固分数”即为该质量的铜液凝固一定分数所涉及的相变潜热释放。基于节点间温度线性分布假设、固相率与温度成线性关系假设，式(31)有如下等价形式：
Figure BDA0004107600070000224
将式(32)加入节点(2)的热平衡方程即在数学关系层面实现了相变潜热释放对传热的影响。值得注意的，当式(32)中v取值为0时Qlatent等于0，然而这并不意味着实际凝固过程中无相变潜热释放。在铸坯牵引速度为0、边界条件保持不变的情况下：非稳态传热时，固液界面仍未停止向液相中推进，此时相变仍继续发生，从而有相变潜热释放；稳态传热时，固液界面停止向液相中推进，此时能认为无相变发生，从而无相变潜热释放。由此能见，式(32)中v取值为0时Qlatent等于0是有物理意义的，正确反映了特殊情况下的潜热释放情况。
热平衡方程组求解：
对于图1所示的平板传热模型，其所涉及的热平衡方程均为线性方程，故其热平衡方程组的求解采用一般的线性方程组解法即能。然而对于铜管连铸模型，其所涉及的热平衡方程不全是线性的，其中部分方程的系数与未知的节点温度有关，也有部分方程包含非线性的源项；在此情况下，迭代法是热平衡方程组的唯一能行解法。
采用迭代法求解非线性热平衡方程组的过程包括如下几步：开始求解时，对所有节点估计一组T值；将估计的T值代入热平衡方程系数表达式，得到名义上的线性方程组；解名义上的线性方程组，得到新一组T值；以新得到的一组T值作为相较更好的估计值，重复执行2、3、4步，直到这种进一步的重复计算(称之为迭代)不再引起T值显著变化时为止。
这种最终T值不再显著变化的状态称之为迭代的收敛，尽管这个解是由线性方程组求解得到的，但它实际上就是非线性热平衡方程组的正确解。然而也能能经过迭代不会收敛到一个解，T值稳定地飘移或是以一个不断增大的幅度震荡，这一与收敛相对立的过程称之为发散。对于铜管连铸模型，遵循前文所述的建模原理就基本能够保证收敛；为了消除类似于辐射边界条件的非线性元素对迭代收敛的不利影响，能针对这些非线性元素应用一些线性化及松弛的方法。
在所有的迭代法当中，最简单的一种是高斯-赛德尔(Gauss-Seidel)逐点计算法。该法首先要为所有节点估计一组T值，然后开始按一定的顺序逐个更新每个节点的T值。T值的更新是依据式(9)那样的代数关系，即当前节点处的热平衡方程；不同的是，此处方程中所引用的相邻节点的T值必须是其最新值，也即T的初始估计值或最近更新值；重复多轮更新所有节点的T值，即能逐渐获得方程组的收敛解。高斯-赛德尔法的最大缺点是其收敛速度太慢，当模型节点数量较大时尤为明显；由于本方法是以一次迭代前进一个节点间距的速率来传递边界条件信息到模型内部的，故其收敛速度慢是能理解的。
为了克服高斯-赛德尔法收敛速度慢的劣势，将其单次迭代更新一个节点改为更新一行节点。如图9所示，将图中圆点标识处作为当前节点行，以相邻节点行的最新T值作为已知条件，建立当前行各个节点处的热平衡方程并进行联立求解，将所得新一组T值作为当前行的最新T值，如此便完成了一次迭代。显然，上述迭代步中要进行求解的方程组是三对角型方程组，该类型方程组能采用TDMA追赶法[71]进行求解，算法实现比较简单。对于图4所示的传热模型，能采用上述的逐行更新方式在横纵两个方向交替地扫描，如此即能将所有的边界条件信息传递到模型内部，最终得到收敛解。逐行更新的高斯-赛德尔迭代的收敛速度更快，这是因为不管沿行有多少个节点，在行两端的边界条件信息都能立即传递到模型内部。
在线计算程序开发：
基于前文工作开发了铜管水平连铸稳态温度场在线计算程序，程序结构如图10所示。在“系统参数加载”环节：(1)程序访问系统盘根目录下的“root.txt”文件，从中读出用户事先设定的工作目录所在位置。工作目录用于存储从上位机传来的实时工艺参数、程序运行过程中产生的临时文件以及计算结果文件。(2)程序访问工作目录下的“inter.txt”文件，从中读出用户事先设定的工艺参数采集节拍。工艺参数采集节拍包含两个信息：采集次数N和采集间歇Δt。从上位机传来的数据更新频率较快且存在一定扰动，为减小数据扰动，每次计算前采集N组工艺参数，相邻两次采集间隔Δt秒，最后最将N组工艺参数的算术平均作为本次计算的输入，这就是程序在“工艺参数加载”环节执行的操作。
在“模型参数加载”环节，程序访问工作目录下的“modata.txt”文件，从中读出用户事先设定的模型几何尺寸、几何离散参数、材料热物理属性、边界条件及界面换热系数。在“几何离散化”环节，程序根据用户设定的模型几何尺寸和几何离散参数在各传热区域划分网格，并按从1开始依次递增1的正整数序列对网格节点进行编号。为便于后续引用，将节点编号按节点在其所属传热区域内的空间位置存入一个二维数组，数组的行列数等于当前传热区域内网格节点的行列数；为每个传热区域建立一个这样的数组，如此在后续建模过程中需要对某一节点进行操作时，就能方便地通过其所属的传热区域以及其在所属传热区域内的位置引用其编号。
基于前文所述建模原理，能够建立起模型中所有典型节点处的热平衡方程。在“热平衡系数表构建”环节，程序根据典型节点热平衡方程的系数形式计算出所有同类型节点处热平衡方程的系数，并将系数向量按节点在其所属传热区域内的空间位置存入一个二维数组，数组的行列数等于当前传热区域内网格节点的行列数；为每个传热区域建立一个这样的数组，如此在后续的方程组迭代求解过程中能直接引用这些系数而不必重复计算，能极大提高程序的时间效率。
在“定义固相率”、“定义相变潜热”及“定义运动传热”环节，将数学模型定义为相应的函数形式，如此在构建热平衡方程组时即能方便地引用这些函数，而不必为同样的计算过程重复编写代码，从而极大提高程序的空间效率。在“热平衡方程组构建”环节，程序基于之前四个环节的准备工作建立起整个模型的热平衡方程组。
在“方程组解初始化”环节，程序为所有节点估计一组T值；以这组T值为起点，程序在“方程组求解”环节采用迭代法求解已经构建好的热平衡方程组。由图10能见，当前程序在整体上是一个无限循环结构，第一次循环中，程序在“方程组解初始化”环节为同一个传热区域内的节点估计相同的T值，不同传热区域的T值按经验大致设定；后续循环中，程序将上次循环得到的热平衡方程组的收敛解作为“方程组解初始化”环节T的估计值，由于在实际生产中前后两次计算对应的工艺参数一般相差不大，能预计前后两次计算的收敛解也相差不大，故这种“解初始化更新”的技巧能减少后续循环中获得收敛解所需的迭代次数，从而提高程序的时间效率。
在“控制台输出”环节，程序将本次计算开始的系统时间、所采用的工艺参数、计算收敛情况、获得收敛解所经过的迭代次数以及壁钟时间等信息输出到控制台窗口，从而实时地向操作人员显示程序运行状态；此外，相同的信息在“系统日志更新”环节被追加写入工作目录下的“log.txt”文件，以备程序运行意外中断时对故障进行诊断。
如图11、12所示，在“绘制温度场云图”、“绘制固相率云图”环节，程序分别将温度场计算收敛解以及对应的固相率分布绘制为彩色图像，不同颜色对应场量值的不同水平，以便操作人员更直观地查看计算结果。为使温度场计算结果的定量查看更加便捷，在温度场计算结果彩图中添加了带有温度值的等温线。
本发明未尽事宜为公知技术。
尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，能理解在不脱离本发明的原理和精神的情况下能对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。

Claims (10)

1.一种铜管水平连铸温度场在线计算方法，其特征在于：所述的铜管水平连铸温度场在线计算方法包括：连铸温度场在线计算，几何模型建立及离散化，材料属性设定与实现，边界条件设定与实现，面属性设定与实现，运动传热模型建立，相变潜热模型建立，平衡方程组求解。
2.根据权利要求1所述的铜管水平连铸温度场在线计算方法，其特征在于：所述的连铸温度场在线计算：
建模原理：有如下经验公式成立：
Figure FDA0004107600060000011
式中，Qn是指定截面法向上的热流，单位W；K是比例常数，即特定材料的热导率，单位W/m/℃；A是指定截面的面积，单位m2；dT/dn是沿指定截面法向的温度变化率，即温度梯度，单位℃/m；该式即为傅立叶定律的微分表达形式；
平板传热模型，设平板为有限大，其厚度均匀、内部连续、材质均一，前侧表面保持恒定温度Tconst(℃)，右侧表面暴露于Tatm(℃)恒温空气中、换热系数h(W/m2/℃)；
边界条件的外表面皆保持绝热，在指定边界条件的持续作用下，具有任意初始温度分布的平板将逐渐改变其温度分布，足够长时间后其温度分布不再随时间变化，此时板内通过的热流保持恒定，称这一状态为稳态传热，称此时的温度分布为稳态温度场；
当工艺参数稳定，铜管水平连铸的理想状态就是稳态传热，即连铸结晶器内温度分布不随时间改变；
若能求得全部网格节点处的稳态温度值，也就近似求得了稳态温度场；当网格密度足够大时，即便在局部节点之间进行简单线性插值，也能足够精确地得出二维传热区域内任意位置的稳态温度值；按照节点所处位置在传热区域内部抑或边界、所处边界处的边界条件类型，能将全部节点分为6类；
所处位置为传热区域的下边界，该边界已设置恒温边界条件，故有：
T4＝Tconst (2)
式中T4为节点4处的温度；本申请文件中，无特殊说明时数字下标皆表示节点编号；
传热区域的左边界，该边界已设置绝热边界条件，单位时间从控制体积的右边界流入控制体积的热量为：
Figure FDA0004107600060000021
式中微分项的下角标e表示右侧，即east；类似的，下文中出现的下角标w、n、s分别表示左侧、上侧、下侧；假设节点间温度线性分布，将式中微分项用差商代替能得：
Figure FDA0004107600060000022
式中K为平板材料的热导率；Δ为网眼边长；A为节点(22)控制体积的右边界面积，A＝Δ·Th，其中，Th为平板厚度；同理能得单位时间从节点(22)控制体积的上边界和下边界流入控制体积的热量分别为：
Figure FDA0004107600060000023
Figure FDA0004107600060000024
由于稳态传热时节点(22)处温度保持不变，由此能知单位时间内节点(22)控制体积的净流入热量为0，故有如下热平衡方程成立：
Q22n+Q22s+Q22e＝0 (7)
整理能得：
T15+2T23+T29-4T22＝0 (8)
节点(25)所处位置为传热区域的内部，故节点(25)只受编号为18、26、32、24的4个相邻节点影响；根据“单位时间内节点控制体积净流入热量为0”的结论，能得节点(25)与4个相邻节点的温度关系为：
T18+T26+T32+T24-4T25＝0 (9)
节点(28)所处位置为传热区域的右边界，该边界已设置对流换热边界条件，故节点(28)除了受编号为27、21、35的3个相邻节点影响，还受到对流换热边界条件影响；根据牛顿冷却定律，单位时间从节点(28)控制体积的右边界流入控制体积的热量为：
Q28e＝hA(Tatm-T28) (10)
根据式(1)，有：
Figure FDA0004107600060000031
Figure FDA0004107600060000032
Figure FDA0004107600060000033
根据“单位时间内节点控制体积净流入热量为0”的结论，有：
Q28w+Q28s+Q28e+Q28n＝0 (14)
整理能得：
P27T27+P21T21+P35T35+P28T28+C28＝0 (15)
设材料热导率、对流换热系数为常数，则式中P和C为常数；
节点(43)所处位置为传热区域的左上角点，其所属的两个边界均已设置绝热边界条件，故节点(43)只受编号为36、44的2个相邻节点影响；根据“单位时间内节点控制体积净流入热量为0”的结论，能得节点(43)与2个相邻节点的温度关系为：
T36+T44-2T43＝0 (16)
节点(49)所处位置为传热区域的右上角点，其所属的传热区域上边界已设置绝热边界条件，其所属的传热区域右边界已设置对流换热边界条件，基于上文结论不难得出：
Figure FDA0004107600060000034
Figure FDA0004107600060000041
Figure FDA0004107600060000042
Q49w+Q49s+Q49e＝0 (20)
整理能得：
P48T48+P42T42+P49T49+C49＝0 (21)
6个典型节点与各自相邻节点的温度关系已全部找到，由此易知全部节点与各自相邻节点的温度关系；分析6个典型节点与各自相邻节点的温度关系，发现其数学形式都属于线性方程；由于在物理范畴这些关系同时成立，故能在数学范畴将其联立，得到能表示全部节点之间温度关系的数学形式，即线性方程组；该线性方程组具有唯一解，该解即是图1中传热模型的稳态温度场。
3.据权利要求1所述的铜管水平连铸温度场在线计算方法，其特征在于：所述的几何模型建立及离散化：
冷水中发挥主要换热作用的，是与铜套相接触的筒形水冷区域，为提高建模效率能忽略其他部分；石墨模进液口以左的部分浸没在浇铸温度的铜液中，内部温度均匀等于铜液浇铸温度，为提高建模效率能忽略这一部分；
对简化后的几何模型进行离散处理，将两个筒形水冷区域离散成一维点阵，石墨芯棒带有一定锥度，鉴于石墨筒外部的铸坯区域温度梯度较小，为其设置了较小的轴向离散密度；为便于界面换热的建模，令相邻点阵在分界处相互重合。
4.根据权利要求1所述的铜管水平连铸温度场在线计算方法，其特征在于：所述的材料属性设定与实现：
将密度描述为温度的函数，其在热平衡方程中的实现，对于一维传热模型，将其材料热导率描述为温度的函数，图中w、e分别位于节点1、2及节点(2)、3的中点，w、e之间的材料区域即为节点(2)的控制体积；对于节点(2)有如下热平衡方程成立：
Figure FDA0004107600060000043
式中Kw、Ke分别为w、e处的材料热导率，其中Ke往往被简单取值为K2与K3的算术平均，但这样的简单取值并不准确；
基于热导率在相邻控制体积之间发生阶跃变化的假设，能知图5中节点(2)和e之间的材料热导率为K2，e和节点3之间的材料热导率为K3，从而能将e处热流密度qe表示为：
Figure FDA0004107600060000051
式中最右侧分式的分母部分表示节点(2)、3之间的总热阻；由上式不难得出：
Figure FDA0004107600060000052
显然，从准确表示e处热流密度的角度出发，Ke的取值应该是K2与K3的调和平均，而非算术平均。
5.根据权利要求1所述的铜管水平连铸温度场在线计算方法，其特征在于：所述的边界条件设定与实现：
石墨筒进液端，相邻表面完全浸没于浇铸温度的铜液；为简化模型，于该处设置恒温边界条件近似取代铜液的加热作用，温度值设定为铜液浇铸温度即能；该边界条件在热平衡方程中的实现参考式(2)；
石墨筒与炉壁的接触面，热量经该接触面由石墨筒流向炉壁；为简化模型，于该处设置对流换热边界条件近似取代炉壁的冷却作用；流体温度参考炉壁温度实测值设定为1050℃，换热系数设定为3000W/m/℃；该边界条件在热平衡方程中的实现参考式(10)；
一冷区和二冷区之间，表面温度较高，必须考虑辐射换热；该处也有氮气、水汽通过，对流换热同样不容忽视，故能于该处设置“对流换热+辐射换热”混合边界条件；其中对流换热边界条件的实现参考式(10)，辐射换热边界条件依据斯特藩-玻尔兹曼定律(Stefan-Boltzmann law)实现：
Q＝ε1A1σF12(T1 4-T2 4) (25)
式中Q为表面1与表面2之间的辐射传热功率；ε1为表面1的发射率；A1为表面1的面积；ζ为斯特藩-玻尔兹曼常数；F12为表面1到表面2的辐射角系数；T1、T2分别为表面1、表面2的温度。
6.据权利要求1所述的铜管水平连铸温度场在线计算方法，其特征在于：所述的界面属性设定与实现：
一维传热模型包含W、E两个传热区域，W区域的右端点(节点w1处)与E区域的左端点(节点e1处)之间有热导率为K、厚度为δ的界面薄层；
对于传热模型，考察节点w1、e1之间的热流，根据式(1)有：
Figure FDA0004107600060000061
式中Q为节点w1、e1之间的热流，A为界面面积；令h＝K/δ，则式(26)变为：
Figure FDA0004107600060000062
7.据权利要求1所述的铜管水平连铸温度场在线计算方法，其特征在于：所述的运动传热模型建立：铜管水平连铸过程中的冷却水流动、铜液流动、铸坯运动在本质上影响着稳态温度场，并建立起相应的数学模型；稳态传热时，左右边界之间存在一定温差，由此推知铸坯在通过控制体积过程中必然发生了相应程度的热量损失，将这部分热量表示为Qtransport，则：
Qtransport＝ρvACp(T2e-T2w) (28)
式中ρ为铸坯的密度，将密度描述为温度的函数时，能取ρ的值为ρ(T2)；v为铸坯牵引速度；A为节点(2)控制体积左右边界的面积；Cp为铸坯的比热容；T2w、T2e分别为节点(2)控制体积左右边界处的温度；
式(28)中，ρvA即等于单位时间内通过节点(2)控制体积的铸坯总质量，“质量×比热容×温度变化量”即为该质量的铸坯改变一定温度所涉及的热量变化；基于节点间温度线性分布假设，式(28)有如下等价形式：
Figure FDA0004107600060000063
8.根据权利要求1所述的铜管水平连铸温度场在线计算方法，其特征在于：所述的相变潜热模型建立：
铜管水平连铸过程中相变潜热的释放同样在本质上影响着稳态温度场；假设金属在凝固过程中其固相体积分数(即固相率fs，0≤fs≤1)与温度成线性关系，则能按下式计算fs：
Figure FDA0004107600060000071
因铸坯被匀速向右牵出，故在单位时间内总有一定质量的铜液从左向右流经节点(2)控制体积；稳态传热时，节点(2)控制体积左右边界之间存在一定温差，故也存在一定固相率之差，由此推知铜液在流经控制体积过程中必然发生了相应程度的相变潜热释放，将这部分热量表示为Qlatent，则：
Figure FDA0004107600060000072
式中ρ为铸坯的密度；v为铸坯牵引速度；A为节点(2)控制体积左右边界的面积；L为铜液的结晶潜热；
Figure FDA0004107600060000073
分别为节点(2)控制体积左右边界处的固相率；
式(31)中，ρvA即等于单位时间内流经节点(2)控制体积的铜液总质量，“质量×结晶潜热×凝固分数”即为该质量的铜液凝固一定分数所涉及的相变潜热释放；基于节点间温度线性分布假设、固相率与温度成线性关系假设，式(31)有如下等价形式：
Figure FDA0004107600060000074
将式(32)加入节点(2)的热平衡方程即在数学关系层面实现了相变潜热释放对传热的影响。
9.根据权利要求1所述的铜管水平连铸温度场在线计算方法，其特征在于：所述的热平衡方程组求解：
平板传热模型，其所涉及的热平衡方程均为线性方程，故其热平衡方程组的求解采用一般的线性方程组解法即能；然而对于铜管连铸模型，其所涉及的热平衡方程不全是线性的，其中部分方程的系数与未知的节点温度有关，也有部分方程包含非线性的源项；在此情况下，迭代法是热平衡方程组的唯一能行解法；
采用迭代法求解非线性热平衡方程组的过程包括如下几步：开始求解时，对所有节点估计一组T值；将估计的T值代入热平衡方程系数表达式，得到名义上的线性方程组；解名义上的线性方程组，得到新一组T值；以新得到的一组T值作为相较更好的估计值，重复执行2、3、4步，直到这种进一步的重复计算，称之为迭代，不再引起T值显著变化时为止。
10.根据权利要求1所述的铜管水平连铸温度场在线计算方法，其特征在于：所述的在线计算程序开发：
在系统参数加载环节：(1)程序访问系统盘根目录下的“root.txt”文件，从中读出用户事先设定的工作目录所在位置；工作目录用于存储从上位机传来的实时工艺参数、程序运行过程中产生的临时文件以及计算结果文件；(2)程序访问工作目录下的“inter.txt”文件，从中读出用户事先设定的工艺参数采集节拍；工艺参数采集节拍包含两个信息：采集次数N和采集间歇Δt；从上位机传来的数据更新频率较快且存在一定扰动，为减小数据扰动，每次计算前采集N组工艺参数，相邻两次采集间隔Δt秒，最后最将N组工艺参数的算术平均作为本次计算的输入，这就是程序在“工艺参数加载”环节执行的操作；
在“模型参数加载”环节，程序访问工作目录下的“modata.txt”文件，从中读出用户事先设定的模型几何尺寸、几何离散参数、材料热物理属性、边界条件及界面换热系数；在“几何离散化”环节，程序根据用户设定的模型几何尺寸和几何离散参数在各传热区域划分网格，并按从1开始依次递增1的正整数序列对网格节点进行编号；为便于后续引用，将节点编号按节点在其所属传热区域内的空间位置存入一个二维数组，数组的行列数等于当前传热区域内网格节点的行列数；为每个传热区域建立一个这样的数组，如此在后续建模过程中需要对某一节点进行操作时，就能方便地通过其所属的传热区域以及其在所属传热区域内的位置引用其编号；
建立起模型中所有典型节点处的热平衡方程，在热平衡系数表构建环节，程序根据典型节点热平衡方程的系数形式计算出所有同类型节点处热平衡方程的系数，并将系数向量按节点在其所属传热区域内的空间位置存入一个二维数组，数组的行列数等于当前传热区域内网格节点的行列数；为每个传热区域建立一个这样的数组，如此在后续的方程组迭代求解过程中能直接引用这些系数而不必重复计算，能极大提高程序的时间效率；
在“定义固相率”、“定义相变潜热”及“定义运动传热”环节，将数学模型定义为相应的函数形式，如此在构建热平衡方程组时即能方便地引用这些函数，而不必为同样的计算过程重复编写代码，从而极大提高程序的空间效率；在“热平衡方程组构建”环节，程序基于之前四个环节的准备工作建立起整个模型的热平衡方程组；
在“方程组解初始化”环节，程序为所有节点估计一组T值；以这组T值为起点，程序在“方程组求解”环节采用迭代法求解已经构建好的热平衡方程组；当前程序在整体上是一个无限循环结构，第一次循环中，程序在“方程组解初始化”环节为同一个传热区域内的节点估计相同的T值，不同传热区域的T值按经验大致设定；后续循环中，程序将上次循环得到的热平衡方程组的收敛解作为“方程组解初始化”环节T的估计值，由于在实际生产中前后两次计算对应的工艺参数一般相差不大，能预计前后两次计算的收敛解也相差不大，故这种“解初始化更新”的技巧能减少后续循环中获得收敛解所需的迭代次数，从而提高程序的时间效率；
在“控制台输出”环节，程序将本次计算开始的系统时间、所采用的工艺参数、计算收敛情况、获得收敛解所经过的迭代次数以及壁钟时间等信息输出到控制台窗口，从而实时地向操作人员显示程序运行状态；此外，相同的信息在“系统日志更新”环节被追加写入工作目录下的“log.txt”文件，以备程序运行意外中断时对故障进行诊断；
在“绘制温度场云图”、“绘制固相率云图”环节，程序分别将温度场计算收敛解以及对应的固相率分布绘制为彩色图像，不同颜色对应场量值的不同水平，以便操作人员更直观地查看计算结果；为使温度场计算结果的定量查看更加便捷，在温度场计算结果彩图中添加了带有温度值的等温线。
CN202310197215.6A 2023-03-03 2023-03-03 一种铜管水平连铸温度场在线计算方法 Pending CN116432482A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310197215.6A CN116432482A (zh) 2023-03-03 2023-03-03 一种铜管水平连铸温度场在线计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310197215.6A CN116432482A (zh) 2023-03-03 2023-03-03 一种铜管水平连铸温度场在线计算方法

Publications (1)

Publication Number Publication Date
CN116432482A true CN116432482A (zh) 2023-07-14

Family

ID=87089769

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310197215.6A Pending CN116432482A (zh) 2023-03-03 2023-03-03 一种铜管水平连铸温度场在线计算方法

Country Status (1)

Country Link
CN (1) CN116432482A (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116933185A (zh) * 2023-07-26 2023-10-24 常州润来科技有限公司 一种铜管精整复绕的评估方法及系统
CN117236145A (zh) * 2023-11-16 2023-12-15 山东工商学院 一种封闭母线接头温度预测与控制方法
CN118504472A (zh) * 2024-07-18 2024-08-16 中国人民解放军国防科技大学 静态水冷却效能的等效热容预测方法、装置、设备及介质

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116933185A (zh) * 2023-07-26 2023-10-24 常州润来科技有限公司 一种铜管精整复绕的评估方法及系统
CN117236145A (zh) * 2023-11-16 2023-12-15 山东工商学院 一种封闭母线接头温度预测与控制方法
CN118504472A (zh) * 2024-07-18 2024-08-16 中国人民解放军国防科技大学 静态水冷却效能的等效热容预测方法、装置、设备及介质
CN118504472B (zh) * 2024-07-18 2024-09-24 中国人民解放军国防科技大学 静态水冷却效能的等效热容预测方法、装置、设备及介质

Similar Documents

Publication Publication Date Title
CN116432482A (zh) 一种铜管水平连铸温度场在线计算方法
Hu et al. Mathematical modelling of solidification and melting: a review
Dong et al. Determination of interfacial heat-transfer coefficient during investment-casting process of single-crystal blades
Wang et al. A comprehensive numerical model for melting with natural convection
Zhao et al. Front-tracking finite element method for dendritic solidification
Byun et al. Frost modeling under cryogenic conditions
Bangian-Tabrizi et al. An optimization strategy for the inverse solution of a convection heat transfer problem
Bartrons et al. A finite volume method to solve the frost growth using dynamic meshes
CN106874591B (zh) 一种方坯加热过程温度分布的计算方法
Wang et al. Prediction of the main characteristics of the shell and tube bundle latent heat thermal energy storage unit using a shell and single-tube unit
Dai et al. On the optimal heat source location of partially heated energy storage process using the newly developed simplified enthalpy based lattice Boltzmann method
Hadała et al. Experimental identification and a model of a local heat transfer coefficient for water–air assisted spray cooling of vertical low conductivity steel plates from high temperature
Kumar et al. Modeling moving-boundary problems of solidification and melting adopting an arbitrary Lagrangian-Eulerian approach
Fredman et al. Model for temperature profile estimation in the refractory of a metallurgical ladle
Zhan et al. An improved equivalent heat capacity method to simulate and optimize latent thermal energy storage units
Srisuma et al. Analytical solutions for the modeling, optimization, and control of microwave-assisted freeze drying
Hamzaoui et al. A 2D1/2 model for natural convection and solidification in a narrow enclosure
Huang et al. Moisture transport and diffusive instability during bread baking
Xu et al. Inverse method with heat and entropy transport in solidification processing of materials
Hetmaniok et al. Solution of the direct alloy solidification problem including the phenomenon of material shrinkage
Celentano et al. Modeling natural convection with solidification in mould cavities
Virzi Finite element analysis of the thermal history for Czochralski growth of large diameter silicon single crystals
Wu et al. A novel algorithm for solving the classical Stefan problem
Heng et al. Efficient reconstruction of local heat fluxes in pool boiling experiments by goal-oriented adaptive mesh refinement
Berdnikov et al. Numerical simulation of crystal growth processes by means of horizontal unidirectional crystallization from melts with different Prandtl numbers

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination