CN114417675A - 一种异型坯连铸凝固传热的有限元计算方法 - Google Patents

一种异型坯连铸凝固传热的有限元计算方法 Download PDF

Info

Publication number
CN114417675A
CN114417675A CN202210092730.3A CN202210092730A CN114417675A CN 114417675 A CN114417675 A CN 114417675A CN 202210092730 A CN202210092730 A CN 202210092730A CN 114417675 A CN114417675 A CN 114417675A
Authority
CN
China
Prior art keywords
heat transfer
unit
continuous casting
node
finite element
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
CN202210092730.3A
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.)
Northeastern University China
Original Assignee
Northeastern University China
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 Northeastern University China filed Critical Northeastern University China
Priority to CN202210092730.3A priority Critical patent/CN114417675A/zh
Publication of CN114417675A publication Critical patent/CN114417675A/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]
    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B22—CASTING; POWDER METALLURGY
    • B22D—CASTING OF METALS; CASTING OF OTHER SUBSTANCES BY THE SAME PROCESSES OR DEVICES
    • B22D11/00—Continuous casting of metals, i.e. casting in indefinite lengths
    • B22D11/16—Controlling or regulating processes or operations
    • B22D11/22—Controlling or regulating processes or operations for cooling cast stock or mould
    • B22D11/225—Controlling or regulating processes or operations for cooling cast stock or mould for secondary cooling
    • 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
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/20—Drawing from basic elements, e.g. lines or circles
    • G06T11/203—Drawing of straight lines or curves
    • 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)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mechanical Engineering (AREA)
  • Continuous Casting (AREA)

Abstract

本发明公开了一种异型坯连铸凝固传热的有限元计算方法，涉及钢连铸技术领域。该方法包括根据异型坯断面尺寸建立1/2异型坯二维几何模型；对1/2异型坯二维几何模型进行网格划分，建立1/2异型坯二维有限元模型；针对异型坯连铸凝固传热问题，建立1/2异型坯二维有限元凝固传热模型；利用1/2异型坯二维有限元凝固传热模型计算连铸过程中的温度场变化。本发明充分考虑了异型坯复杂的断面形状，采用三角形单元和一般四边形单元，可以更好地模拟实际异型坯横截面；根据冷区不同选择合适的边界条件，更准确地模拟异型坯连铸生产全流程凝固传热，描述了异型坯连铸过程的温度场变化。

Description

一种异型坯连铸凝固传热的有限元计算方法
技术领域
本发明涉及钢连铸技术领域，尤其涉及一种异型坯连铸凝固传热的有限元计算方法。
背景技术
所谓异型坯是指除了方坯、板坯、圆坯、矩形坯以外具有复杂断面的连铸坯，其中应用最为广泛的就是H型钢。H型钢是一种经济断面型材，具有壁薄、翼缘内外侧平行及腿端为直角等特点，力学性能和使用性能十分优越，广泛应用于国民经济建设的各个领域之中，如桥梁、大型船舶、抗震设施建筑以及高层建筑等方面均大量采用H型钢。目前我国H型钢发展的空间很大，尤其大中规格有较大缺口，产品质量也与国外先进水平有一定的差距。异型坯质量缺陷产生原因有很多。相较于横截面相对简单的方坯、板坯等常规铸坯，横截面更为复杂的异型坯在生产过程中还会产生一些特殊质量缺陷。异型坯质量缺陷可以归纳为表面缺陷和内部缺陷两类。表面缺陷主要包括表面纵裂、角部纵裂、表面横裂、星状裂纹、表面夹渣、气孔和凹坑。内部缺陷主要包括中心裂纹、矫直裂纹、中心偏析、中心疏松和三角区裂纹。
表面纵裂表现为在铸坯长度方向上呈纵向分布的裂纹，裂纹深度小于2mm，腹板和翼缘上都有可能存在，以腹板为多，是异型坯的主要质量缺陷之一。表面纵裂产生的主要原因是冷却水配比不当等。异型坯断面上各点散热条件差别较大，冷却不均匀，表面不同位置温度差较大，导致异型坯的表面裂纹发生率常常高于常规铸坯。表面横裂表现为在铸坯长度上呈横向分布的裂纹，裂纹深度小于2mm。高二冷比水量的情况下，铸坯进入矫直区温度低，异型坯表面横裂发生率较高。矫直裂纹表现为断面翼梢处沿腹板厚度方向呈纵向分布的裂纹，裂纹深度一般小于2mm。当热铸坯在矫直时，如果正处于低温脆性区，易产生矫直裂纹。中心偏析表现为断面腹板中心处沿厚度方向呈横向分布的长条形状，长度一般在3～4mm；中心疏松表现为断面上腹板中心线处分布有细微的孔隙，中心偏析和中心疏松往往相伴存在，恶化钢的力学性能，降低了钢的韧性和耐蚀性，严重影响产品质量。当二冷区采用强冷制度时，易造成中心偏析和中心疏松。角区裂纹表现为断面翼缘处沿腹板厚度方向呈纵向分布的裂纹，裂纹深度一般大于2mm。当翼缘板鼓肚、翼缘板冷却速度太快、冷却不对称时，可引起异形坯三角区裂纹。
因此铸坯在结晶器和二冷区的均匀冷却对于减少铸坯表面和内部缺陷起着至关重要的作用。想要实现铸坯的均匀冷却，就要先建立异型坯在线调控系统，快速精准地预测出每个时刻铸坯的温度变化情况，并实时调节冷却水量。而异型坯在线调控系统的关键就在于凝固传热核心算法的确定，即异型坯凝固传热模型的建立。
连铸过程的边界条件比较复杂，目前的各种商业软件无法在短时间内对其凝固传热过程进行计算，而现有的个人开发的在线凝固传热模型主要应用于方坯、板坯等常规铸坯，针对异型坯的并不多。例如文献号为CN 103192048 B的中国专利《一种基于精准热物性参数的连铸坯凝固冷却过程模拟方法》需要依托连铸冷却过程模拟与优化软件，且目前仅针对中厚板坯，无法拓展到异型坯；文献号为CN 109446748 A的中国专利申请《一种模拟连铸圆坯凝固过程的方法》通过建立圆坯凝固传热数学模型，模拟连铸圆坯凝固过程，仅适用于连铸圆坯，采用的数值模拟方法为有限差分法，对不规则区域处理繁琐且对区域的连续性等要求较严；文献号为CN 111199119 A的中国专利申请《连铸异形坯坯头温度模拟方法》主要针对异型坯坯头，未考虑连铸坯整体温度变化情况，采用的数值模拟方法为有限体积法且划分正交网格单元，对异型坯复杂断面模拟较差，计算精度难以保证；文献号为CN110941889 A的中国专利申请《连铸异型坯微观和宏观裂纹萌发和扩展的研究方法》通过ANSYS模拟出异型坯在结晶器中的温度场和应力场从而找出铸坯裂纹最易萌生与扩展的位置，但温度场的模拟是借助ANSYS有限元软件，并未探究其原理，无法应用于异型坯在线调控系统。
综上所述，建立异型坯凝固传热模型，计算连铸过程各时刻温度场对于实时调控冷却水，实现连铸坯均匀冷却，提高铸坯质量有着重要意义，但现有模型还未有针对异型坯连铸全流程的凝固传热计算。
发明内容
针对上述现有技术存在的不足，本发明提供一种异型坯连铸凝固传热的有限元计算方法，旨在预测异型坯凝固过程各时刻温度变化情况，推进异型坯模拟和传热分析并为建立异型坯在线凝固传热模型，优化二冷区配水，减少异型坯表面缺陷，提高异型坯质量提供理论基础。
本发明的技术方案为：
一种异型坯连铸凝固传热的有限元计算方法，该方法包括如下步骤：
步骤1：根据异型坯断面尺寸，利用有限元分析软件建立1/2异型坯二维几何模型；
步骤2：对1/2异型坯二维几何模型进行网格划分，建立1/2异型坯二维有限元模型并保存网格信息；
步骤3：针对异型坯连铸凝固传热问题，建立1/2异型坯二维有限元凝固传热模型；
步骤4：利用1/2异型坯二维有限元凝固传热模型计算连铸过程中的温度场变化。
进一步地，根据所述的异型坯连铸凝固传热的有限元计算方法，所述网格信息包括单元信息、节点信息、边界信息；所述单元信息包括单元编号以及单元中包含的节点；所述节点信息包括节点编号以及节点坐标；所述边界包括上边界即内弧侧，下边界即外弧侧，左边界或者右边界即侧弧侧；所述边界信息包括边界上具有的单元信息和节点信息。
进一步地，根据所述的异型坯连铸凝固传热的有限元计算方法，所述步骤3包括如下步骤：
步骤3.1：构建有限元传热基本方程；
步骤3.2：构建瞬态传热问题的有限元计算模型；
步骤3.3：构建平面三节点三角形单元，并基于该单元针对不同的传热边界条件进行传热分析；
步骤3.4：构建平面四节点四边形单元，并基于该单元针对不同的传热边界条件进行传热分析。
进一步地，根据所述的异型坯连铸凝固传热的有限元计算方法，所述平面平面三节点三角形单元的构建方法为：以逆时针方向将3个节点编码为i、j、m构建平面三节点三角形传热单元；所述平面四节点四边形单元的构建方法为：选择四个顶点作为插值节点，并以逆时针方向将节点编码为i，j，k和l构建典型的平面四节点四边形单元。
进一步地，根据所述的异型坯连铸凝固传热的有限元计算方法，所述基于平面四节点四边形单元针对不同的传热边界条件进行传热分析的内容为：在参数空间中选择正方形作为参考正方形单元；将所述正方形向物理空间中四边形的几何映射函数写成节点形函数的组合形式
Figure BDA0003489694840000031
其中正方形节点形函数为
Figure BDA0003489694840000032
利用有限元法将物理空间任意四边形单元温度刚度矩阵元素对全局坐标函数的计算方式转化为在参数空间中正方形单元对局部坐标的计算方式：
Figure BDA0003489694840000041
物理空间中任意四边形单元热容矩阵元素也在参数空间的正方形单元中进行计算：
Figure BDA0003489694840000042
物理空间中任意四边形单元内热源产生的节点载荷也在参数空间的正方形单元中进行计算：
Figure BDA0003489694840000043
上式中，x，y为笛卡尔坐标系；ξ，η为参数空间的坐标系；xi、xj、xk、xl分别为节点i、j、k、l的横坐标；yi、yj、yk、yl分别为节点i、j、k、l的纵坐标；kmn为单元温度刚度矩阵元素；k为热传导系数；Nm和Nn均为形状函数；a、b为参数空间的积分上限、下限，分别为所述正方形的坐标变化范围；J为雅可比矩阵；|J|为雅可比矩阵行列式；J*为雅可比矩阵的伴随矩阵；cmn代表单元热容矩阵元素；ρ为材料密度；cT为材料比热；pm代表节点载荷；qv为单位体积内材料内部单位体积中产生的热量，称为内热源密度；l为该单元的传热边界长度。
进一步地，根据所述的异型坯连铸凝固传热的有限元计算方法，所述步骤4包括如下步骤：
步骤4.1：读取异型坯钢种信息、连铸机结构参数和连铸工艺参数，计算固相线温度与液相线温度；
步骤4.2：导入1/2异型坯二维有限元模型的网格信息；
步骤4.3：确定时间步长；所述时间步长为每次计算切片移动的时间；所述切片代表的是所述1/2异型坯二维有限元模型；
步骤4.4：判断单元相区，计算单元的物性参数；
步骤4.5：计算单元传热矩阵、节点载荷、单元热容矩阵并集成总传热矩阵、总温度载荷、总热容矩阵；所述集成总传热矩阵就是将各个单元传热矩阵分解为子矩阵，然后将各个节点序列所决定的分块子矩阵按照单元节点编号“对号填充”进总传热矩阵中；所述集成总温度载荷、总热容矩阵均与所述集成总传热矩阵同理；
步骤4.6：求解温度场，也就是求解有限元二维非稳态热传导控制方程
Figure BDA0003489694840000051
同时判断切片位置是否超过空冷区，如果是则计算结束，如果否则切片继续移动，转到步骤4.4。
进一步地，根据所述的异型坯连铸凝固传热的有限元计算方法，所述钢种信息包括钢牌号、化学成分；所述连铸机结构参数包括结晶器高度、二冷区结构参数；所述二冷区结构参数包括二冷区数量、长度、入口位置和出口位置；所述连铸工艺参数包括浇铸温度、拉速、结晶器热流密度、二冷各区水流量及水温、环境温度和异型坯断面尺寸；所述异型坯断面尺寸包括翼梢厚度、宽面宽度、窄面宽度、腹板长度、腹板厚度以及內缘弧半径、翼缘半径、圆角半径。
进一步地，根据所述的异型坯连铸凝固传热的有限元计算方法，所述相区是指液相区、固相区和固液两相区；所述物性参数包括单元密度、单元比热、单元固相率和单元热传导系数。
进一步地，根据所述的异型坯连铸凝固传热的有限元计算方法，所述判断单元相区的方法为：首先通过反距离加权插值法计算所述平面四节点四边形单元的中心温度，以及通过二维线性插值法计算三角形单元的中心温度；最后根据单元中心温度判断单元所处相区。
进一步地，根据所述的异型坯连铸凝固传热的有限元计算方法，采用高斯-赛德尔迭代法求解矩阵方程
Figure BDA0003489694840000052
计算各节点温度。
总体而言，通过本发明所构思的以上技术方案较现有技术具有以下有益效果：本发明提供的一种异型坯连铸凝固传热的有限元计算方法，充分考虑了异型坯复杂的断面形状，采用三角形单元和一般四边形单元，可以更好地模拟实际异型坯横截面；根据冷区不同选择合适的边界条件，更准确地模拟异型坯连铸生产全流程凝固传热。总之，充分考虑了异型坯连铸的实际情况，描述了异型坯连铸过程的温度场变化，从而为异型坯均匀冷却及提高异型坯质量提供支撑和指导。
附图说明
图1为本发明实施例提供的异型坯连铸凝固传热有限元计算方法的流程示意图；
图2为本发明实施例提供的异型坯断面尺寸示意图；
图3为本发明实施例提供的1/2异型坯二维有限元模型示意图；
图4为本发明实施例提供的异型坯在连铸机不同位置的温度分布云图；其中图(a)为结晶器出口处铸坯的温度分布云图；图(b)为二冷1区出口处铸坯的温度分布云图；图(c)为二冷2区出口处铸坯的温度分布云图；图(d)为二冷3区出口处铸坯的温度分布云图；图(e)为二冷4区出口处铸坯的温度分布云图；图(f)为二冷5区出口处铸坯的温度分布云图；
图5为本发明实施例提供的异型坯各特征点温度变化曲线图；
图6为本发明实施例提供的异型坯翼缘表面中心计算温度与实测温度对比图。
具体实施方式
为了便于理解本申请，下面将参照相关附图对本申请进行更全面的描述。附图中给出了本申请的较佳实施方式。但是，本申请可以以许多不同的形式来实现，并不限于本文所描述的实施方式。相反地，提供这些实施方式的目的是使对本申请的公开内容理解的更加透彻全面。
本实施例以国内某钢厂450mm×350mm×90mmQ235钢异型坯连铸过程为例，采用本发明的异型坯连铸凝固传热的有限元计算方法，计算该Q235钢异型坯连铸过程中的各时刻温度变化情况。
本实施例的异型坯连铸凝固传热有限元计算方法，如图1所示，具体包括以下步骤：
步骤1：根据异型坯断面尺寸，利用有限元分析软件建立1/2异型坯二维几何模型；
根据异型坯断面尺寸，确定关键点坐标并导入ANSYS有限元分析软件中，然后在ANSYS中由关键点连线并生成几何面，建立1/2异型坯二维几何模型。所述异型坯断面尺寸包括翼梢厚度、宽面宽度、窄面宽度、腹板长度、腹板厚度以及內缘弧半径、翼缘半径、圆角半径。关键点包括腹板端点、翼梢端点、內缘弧端点、翼缘端点、圆角端点。
本实施例中，异型坯宽面宽度为0.435m，窄面宽度为0.320m，翼梢厚度为0.045m，腹板长度为0.229m，腹板厚度为0.090m，翼缘半径为0.020m，內缘弧半径为0.100m，圆角半径为0.010m。本实施例的异型坯断面尺寸示意图如图2所示。
步骤2：对1/2异型坯二维几何模型进行网格划分，建立1/2异型坯二维有限元模型并保存网格信息；
在ANSYS中确定单元类型和单元形状，并确定单元尺寸，然后对1/2异型坯二维几何模型进行网格划分，建立1/2异型坯二维有限元模型。将网格信息输出为DAT格式文件并保存。但由于异型坯断面形状较为复杂，ANSYS进行区域剖分时无法完全使用四边形单元进行网格划分，存在部分三角形单元。
所述网格信息包括单元信息、节点信息、边界信息。所述单元信息包括单元编号以及单元中包含哪些节点；所述节点信息包括节点编号以及节点坐标；所述边界包括上边界即内弧侧，下边界即外弧侧，左边界或者右边界即侧弧侧，由于为1/2异型坯，无右边界或者无左边界；所述边界信息包括边界上具有的单元信息和节点信息。
本实施例中，单元类型选择PLANE55，单元形状选择四边形单元，单元尺寸选择2mm，由ANSYS有限元软件自动进行区域单元剖分。本实施例的1/2异型坯二维有限元模型如图3所示。
步骤3：针对异型坯凝固传热问题，建立1/2异型坯二维有限元凝固传热模型；
步骤3.1：推导有限元传热基本方程；
传热过程的基本变量就是温度，它是物体中的几何位置以及时间的函数。
根据Fourier传热定律和能量守恒定律，可以建立热传导问题的控制方程，即物体的瞬态温度场T(x,y,z,t)应满足如下方程(1).
Figure BDA0003489694840000071
其中ρ为材料密度，kg/m3；cT为材料比热，J/(kg·K)；κx、κy、κz分别沿x,y,z方向的热传导系数，W/(m·K)；Q(x,y,z,t)为物体内部的热源强度，W/kg。
传热边界条件有三类，即
第一类边界条件S1为Dirichlet条件，在边界上给定温度值
Figure BDA0003489694840000072
第二类边界条件S2为给定热流密度的Neumann条件
Figure BDA0003489694840000073
第三类边界条件S3为给定对流换热的Neumann条件
Figure BDA0003489694840000074
其中nx,ny,nz为边界外法线的方向余弦；
Figure BDA0003489694840000075
为在符合第一类边界条件S1的传热边界上给定的温度，℃；
Figure BDA0003489694840000076
为在符合第二类边界条件S2的传热边界上的给定凝固热流密度，W/m2；
Figure BDA0003489694840000077
为物体与周围介质的对流换热系数，W/(m2·K)；T∞为环境温度；t为时间，s；并且物体Ω的边界为
Figure BDA0003489694840000078
若该问题的初始条件IC为
Figure BDA0003489694840000079
其中，
Figure BDA00034896948400000710
为初始温度，℃。
相应的变分提法为，在满足边界条件S1、S2、S3及初始条件IC的许可温度场中，真实的温度场使以下泛函I取极小值，即
Figure BDA0003489694840000081
在实际问题的处理过程中，边界条件S2和S3事先较难满足，因此，可以将这两个条件耦合进泛函中，即
Figure BDA0003489694840000082
步骤3.2：建立瞬态传热问题的有限元计算模型；
在瞬态传热问题中，单元的温度场将随时间变化，即
Figure BDA0003489694840000083
其中N(x,y,z)为形状函数，这里的节点温度
Figure BDA0003489694840000084
是随时间变化的，即
Figure BDA0003489694840000085
将(8)式代入(7)式中，并对
Figure BDA0003489694840000086
求变分极值，可得到
Figure BDA0003489694840000087
其中
Figure BDA0003489694840000088
Figure BDA0003489694840000089
Figure BDA00034896948400000810
Figure BDA00034896948400000811
其中方程(10)叫做单元传热方程；
Figure BDA00034896948400000812
称为单元传热矩阵；
Figure BDA00034896948400000813
为单元节点温度列阵；
Figure BDA00034896948400000814
为单元热容矩阵；
Figure BDA00034896948400000815
为单元节点等效温度载荷列阵。
步骤3.3：构建平面三节点三角形单元，并基于该单元针对不同的传热边界条件进行传热分析；
以逆时针方向将3个节点编码为i、j、m构建平面三节点三角形传热单元，由3节点i,j,m组成的平面三节点三角形传热单元推导以下三种情形的单元矩阵：
(1)无传热边界，即完全为内部单元；
(2)如果该单元的某条边例如jm边为符合第二类传热边界条件S2的传热边界时：由
Figure BDA0003489694840000091
常数来描述；
(3)如果该单元的某条边例如jm边为符合第三类传热边界条件S3的传热边界时：由
Figure BDA0003489694840000092
常数来描述。
该单元的节点温度列阵为
Figure BDA0003489694840000093
取单元温度场的插值关系为
Figure BDA0003489694840000094
其中形状函数矩阵N为
N＝[Ni Nj Nm] (17)
其中
Figure BDA0003489694840000095
Figure BDA0003489694840000096
Figure BDA0003489694840000097
Figure BDA0003489694840000098
(i,j,m循环)
Figure BDA0003489694840000099
其中A为单元面积，m2。
下面分别就几种不同的传热边界，分别计算相应的传热矩阵
Figure BDA00034896948400000910
和节点等效温度载荷列阵
Figure BDA00034896948400000911
(1)完全为内部单元(无传热边界)
将形状函数表达式(17)代入
Figure BDA00034896948400000912
和
Figure BDA00034896948400000913
的计算公式(11)及(12)中(注意这里仅考虑二维问题)，有
Figure BDA0003489694840000101
Figure BDA0003489694840000102
(2)对于该单元的某条边例如jm边为符合第二类传热边界条件S2的传热边界时
将形状函数表达式(17)和BC(S2)表达式(3)代入
Figure BDA0003489694840000103
和
Figure BDA0003489694840000104
的计算公式(11)及(12)中(注意这里仅考虑二维问题)，有
Figure BDA0003489694840000105
Figure BDA0003489694840000106
其中，L为jm边的长度，m。
(2)对于该单元的某条边例如jm边为符合第三类传热边界条件S3的传热边界时
同样，将形状函数表达式(17)和BC(S2)表达式(4)代入
Figure BDA0003489694840000107
和
Figure BDA0003489694840000108
的计算公式(11)及(12)中，有
Figure BDA0003489694840000109
Figure BDA0003489694840000111
Figure BDA0003489694840000112
步骤3.4：构建平面四节点四边形单元，并基于该单元针对不同的传热边界条件进行传热分析。
选择四个顶点作为插值节点，并以逆时针方向将节点编码为i，j，k和l构建典型的平面四节点四边形单元，单元内部温度场采用双线性插值函数进行计算：
T(x,y)＝a1+a2x+a3y+a4xy (25)
式中，a1、a2、a3、a4为待定系数。
任意平面四节点四边形单元中双线性插值函数的系数计算十分困难，而且传热矩阵、热容矩阵及载荷向量表达式中有些积分项无法得到解析的原函数。所以，有限元法不直接计算任意形状的四边形单元模型，而是采用等参映射方法间接计算四边形单元模型。
本发明在参数空间中选择正方形[-1,1]×[-1,1]作为参考正方形单元。即参数空间中正方形顶点的坐标分别为i(-1,-1)、j(1,-1)、k(1,1)和l(-1,1)。将参数空间中[-1,1]×[-1,1]正方形向物理空间中四边形的几何映射函数写成节点形函数的组合形式，即
Figure BDA0003489694840000113
其中正方形节点形函数为
Figure BDA0003489694840000114
利用有限元法将物理空间任意四边形单元温度刚度矩阵元素对全局坐标函数的计算方式转化为在参数空间中正方形单元对局部坐标的计算方式，即
Figure BDA0003489694840000121
参数空间的积分上下限分别为正方形的坐标变化范围，对于[-1,1]×[-1,1]正方形，则a＝-1，b＝1。
物理空间中任意四边形单元热容矩阵元素也可以在参数空间的正方形单元中进行计算，即
Figure BDA0003489694840000122
物理空间中任意四边形单元内热源产生的节点载荷也在参数空间的正方形单元中进行计算，即
Figure BDA0003489694840000123
其中，qv为单位体积内材料内部单位体积中产生的热量，称为内热源密度，W/m2；J称为雅可比矩阵
Figure BDA0003489694840000124
|J|称为雅可比矩阵行列式
Figure BDA0003489694840000125
J*称为雅可比矩阵的伴随矩阵
Figure BDA0003489694840000126
关于传热边界的处理，由4节点i,j,k,l组成的平面四节点四边形单元，针对不同的传热边界条件进行传热分析，推导以下三种情形的单元矩阵：
(1)无传热边界，即完全为内部单元；
(2)如果该单元的某条边例如ij边为符合第二类传热边界条件S2的传热边界时：由
Figure BDA0003489694840000131
常数来描述；
(3)如果该单元的某条边例如ij边为符合第三类传热边界条件S3的传热边界时：由
Figure BDA0003489694840000132
常数来描述。
(1)无传热边界，即完全为内部单元；
无传热边界，完全为内部单元时，单元温度刚度矩阵元素的计算公式为式(28)，单元热容矩阵元素的计算公式为式(29)，节点载荷的计算公式为式(30)。
(2)如果该单元的某条边例如ij边为符合第二类传热边界条件S2的边界时：由常数来描述；
单元温度刚度矩阵元素的计算公式为式(28)，单元热容矩阵元素的计算公式为式(29)，节点载荷的计算公式如下：
Figure BDA0003489694840000133
(3)如果该单元的某条边例如ij边为符合第三类传热边界条件S3的边界时：由
Figure BDA0003489694840000134
常数来描述。
单元热容矩阵元素的计算公式为式(19)，单元温度刚度矩阵元素和节点载荷的计算公式如下：
Figure BDA0003489694840000135
Figure BDA0003489694840000141
Figure BDA0003489694840000142
上式中，kmn和kmm均为单元温度刚度矩阵元素；k为热传导系数；Nm和Nn均为形状函数；a为参数空间的积分上限；b为参数空间的积分下限；J为雅可比矩阵；|J|为雅可比矩阵行列式；J*为雅可比矩阵的伴随矩阵；x，y为笛卡尔坐标系；ξ，η为参数空间的坐标系；cmn代表单元热容矩阵元素；ρ为材料密度；cT为材料比热；pm代表节点载荷；qv为单位体积内材料内部单位体积中产生的热量，称为内热源密度；L为该单元的传热边界长度；
Figure BDA0003489694840000143
为在符合边界条件S2的边界上的给定凝固热流密度；
Figure BDA0003489694840000144
为物体与周围介质的对流换热系数；T∞为环境温度；t为时间。
步骤4：利用1/2异型坯二维有限元凝固传热模型计算连铸过程中的温度场变化。
步骤4.1：读取异型坯钢种信息、连铸机结构参数和连铸工艺参数，计算固相线温度与液相线温度；
所述钢种信息包括钢牌号、化学成分；所述连铸机结构参数包括结晶器高度、二冷区结构参数；所述二冷区结构参数包括二冷区数量、长度、入口位置和出口位置；所述连铸工艺参数包括浇铸温度、拉速、结晶器热流密度、二冷各区水流量及水温、环境温度和异型坯断面尺寸。所述异型坯断面尺寸包括翼梢厚度、宽面宽度、窄面宽度、腹板长度、腹板厚度以及內缘弧半径、翼缘半径、圆角半径。
本实施例中，钢牌号为Q235，该钢种化学成分如表1所示。
表1国内某钢厂450mm×350mm×90mm Q235异型坯化学成分
化学成分 C Si Mn P S
质量分数(％) 0.19 0.21 0.47 0.025 0.027
本实施列中，异型坯连铸机结晶器高度为0.7m；共有5个二冷区，其中二冷1区长度为0.647m，入口和出口分别距弯月面0.7m和1.347m；2区长度为1.486m，入口和出口分别距弯月面1.347m和2.833m；3区长度为2.337m，入口和出口分别距弯月面2.833m和5.170m；4区长度为2.378m，入口和出口分别距弯月面5.170m和7.548m；5区长度为2.224m，入口和出口分别距弯月面7.548m和9.772m。Q235异型坯的浇铸温度为1545℃，拉速为0.9m/min；结晶器宽面热流密度为247w/m2，窄面热流密度为243.677w/m2。二冷区冷却条件如表2所示。二冷水温为25℃，空冷区环境温度为90℃。
表2国内某钢厂450mm×350mm×90mm Q235异型坯二冷冷却条件
二冷区 1 2 3 4 5
内弧水量(L/min) 36 33 21 13 13
外弧水量(L/min) 36 33 21 13 13
侧弧水量(L/min) 40 37 19 8.4 8.4
总水量(L/min) 112 103 60 34..4 34.4
本实施例中，根据Q235钢的C含量选择相应固液相线温度的计算公式
Figure BDA0003489694840000151
Figure BDA0003489694840000152
其中Tl为液相线温度；Ts为固相线温度；[C]、[Si]、[Mn]、[P]、[S]、[Cr]、[Al]、[Ni]、[V]、[Mo]分别为对应元素的质量分数。
步骤4.2：导入1/2异型坯二维有限元模型的网格信息；
将步骤2导出的DAT格式文件导入计算程序中，并储存起来。所述网格信息包括单元信息、节点信息、边界信息。所述单元信息包括单元编号以及单元中包含哪些节点；所述节点信息包括节点编号以及节点坐标；所述边界信息包括边界上有哪些单元和节点。
步骤4.3：确定时间步长；
所述时间步长即每次计算切片(因二维模型无厚度，称为切片)移动的时间。时间步长乘拉速即为切片每次移动的距离，移动距离累计即为切片所处的位置。
本实施例中，时间步长为0.1s，拉速为0.9m/min。
步骤4.4：判断单元相区，计算单元的物性参数；
四边形单元通过反距离加权插值法计算单元中心温度，三角形单元通过二维线性插值计算单元中心温度，根据单元中心温度判断单元所处相区并计算单元物性参数。
所述相区是指液相区、固相区和固液两相区；
所述物性参数包括单元密度、单元比热、单元固相率和单元热传导系数。
本实施例中，由步骤4.2可知固相线温度为1485℃，液相线温度为1515℃。根据单元温度判断单元相区并采用合适的公式计算单元的物性参数。
步骤4.5：计算单元传热矩阵、节点载荷、单元热容矩阵并集成总传热矩阵、总温度载荷、总热容矩阵；
根据步骤3.2和3.3的公式分别计算各个单元的单元传热矩阵、节点载荷、单元热容矩阵，并对边界上的单元添加边界条件。内弧侧、外弧侧以及侧弧侧根据切片所处连铸机位置的不同根据步骤4.1中读取的结晶器和二冷区冷却情况进行计算。需要注意的是实际生产中喷淋水并非均匀分布，因此计算边界条件时需要考虑边界单元在边界不同位置处的水量。
所述集成总传热矩阵就是将各个单元矩阵分解为子矩阵，然后将各个节点序列所决定的分块子矩阵按照单元节点编号“对号填充”进总传热矩阵中。所述集成总温度载荷、总热容矩阵同理。
步骤4.6：求解温度场，即求解有限元二维非稳态热传导控制方程
Figure BDA0003489694840000161
同时判断切片位置是否超过空冷区，如果是则计算结束，如果否则切片继续移动，转到步骤4.4；
采用高斯-赛德尔迭代法求解矩阵方程
Figure BDA0003489694840000162
计算各节点温度，求解完成后判断切片位置是否超过空冷区，如果已经超过空冷区，计算结束，异型坯连铸凝固传热计算完成，如果没有超过空冷区，则切片移动到下一个位置，转到步骤4.4。
本实施例采用高斯-赛德尔迭代法求解矩阵方程，误差为所有节点两次迭代之差的平方和再开方，误差容限为0.001，为了提高计算效率，采用OpenCV中的并行计算函数parallel_for进行CPU并行计算。
步骤5：对异型坯连铸凝固传热计算结果进行可视化及结果后处理；
异型坯连铸凝固传热计算完成后提取连铸机不同位置处的计算结果，导入Tecplot软件中进行可视化处理，形成温度云图；分别提取腹板中心、內缘中心、窄面中心以及圆角中心的温度变化情况，导入Origin软件中进行可视化处理，形成温度变化曲线。
本实施例提供了国内某钢厂450mm×350mm×90mm Q235异型坯在连铸机不同位置的温度分布云图，如图4所示，其中图(a)为结晶器出口处铸坯的温度分布云图；图(b)为二冷1区出口处铸坯的温度分布云图；图(c)为二冷2区出口处铸坯的温度分布云图；图(d)为二冷3区出口处铸坯的温度分布云图；图(e)为二冷4区出口处铸坯的温度分布云图；图(f)为二冷5区出口处铸坯的温度分布云图。相对应的各特征点温度变化曲线，如图5所示。翼缘表面中心计算温度与实测温度对比，如图6所示。从图6可以看出模拟计算结果与实际温度基本符合。由于异形坯形状特点和冷却特点，随着冷却进行，异形坯中心高温区域逐渐减小，且逐渐向內缘侧移动。由于內缘呈凹状，该处传热效果较差，內缘表面温度高于表面其他位置处温度。异形坯角部位置以二维传热方式，热量散失较快，因此角部温度最低。
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。

Claims (10)

1.一种异型坯连铸凝固传热的有限元计算方法，其特征在于，该方法包括如下步骤：
步骤1：根据异型坯断面尺寸，利用有限元分析软件建立1/2异型坯二维几何模型；
步骤2：对1/2异型坯二维几何模型进行网格划分，建立1/2异型坯二维有限元模型并保存网格信息；
步骤3：针对异型坯连铸凝固传热问题，建立1/2异型坯二维有限元凝固传热模型；
步骤4：利用1/2异型坯二维有限元凝固传热模型计算连铸过程中的温度场变化。
2.根据权利要求1所述的异型坯连铸凝固传热的有限元计算方法，其特征在于，所述网格信息包括单元信息、节点信息、边界信息；所述单元信息包括单元编号以及单元中包含的节点；所述节点信息包括节点编号以及节点坐标；所述边界包括上边界即内弧侧，下边界即外弧侧，左边界或者右边界即侧弧侧；所述边界信息包括边界上具有的单元信息和节点信息。
3.根据权利要求1所述的异型坯连铸凝固传热的有限元计算方法，其特征在于，所述步骤3包括如下步骤：
步骤3.1：构建有限元传热基本方程；
步骤3.2：构建瞬态传热问题的有限元计算模型；
步骤3.3：构建平面三节点三角形单元，并基于该单元针对不同的传热边界条件进行传热分析；
步骤3.4：构建平面四节点四边形单元，并基于该单元针对不同的传热边界条件进行传热分析。
4.根据权利要求3所述的异型坯连铸凝固传热的有限元计算方法，其特征在于，所述平面平面三节点三角形单元的构建方法为：以逆时针方向将3个节点编码为i、j、m构建平面三节点三角形传热单元；所述平面四节点四边形单元的构建方法为：选择四个顶点作为插值节点，并以逆时针方向将节点编码为i，j，k和l构建典型的平面四节点四边形单元。
5.根据权利要求4所述的异型坯连铸凝固传热的有限元计算方法，其特征在于，所述基于平面四节点四边形单元针对不同的传热边界条件进行传热分析的内容为：在参数空间中选择正方形作为参考正方形单元；将所述正方形向物理空间中四边形的几何映射函数写成节点形函数的组合形式
Figure FDA0003489694830000011
其中正方形节点形函数为
Figure FDA0003489694830000021
利用有限元法将物理空间任意四边形单元温度刚度矩阵元素对全局坐标函数的计算方式转化为在参数空间中正方形单元对局部坐标的计算方式：
Figure FDA0003489694830000022
物理空间中任意四边形单元热容矩阵元素也在参数空间的正方形单元中进行计算：
Figure FDA0003489694830000023
物理空间中任意四边形单元内热源产生的节点载荷也在参数空间的正方形单元中进行计算：
Figure FDA0003489694830000024
上式中，x，y为笛卡尔坐标系；ξ，η为参数空间的坐标系；xi、xj、xk、xl分别为节点i、j、k、l的横坐标；yi、yj、yk、yl分别为节点i、j、k、l的纵坐标；kmn为单元温度刚度矩阵元素；k为热传导系数；Nm和Nn均为形状函数；a、b为参数空间的积分上限、下限，分别为所述正方形的坐标变化范围；J为雅可比矩阵；|J|为雅可比矩阵行列式；J*为雅可比矩阵的伴随矩阵；cmn代表单元热容矩阵元素；ρ为材料密度；cT为材料比热；pm代表节点载荷；qv为单位体积内材料内部单位体积中产生的热量，称为内热源密度；l为该单元的传热边界长度。
6.根据权利要求5所述的异型坯连铸凝固传热的有限元计算方法，其特征在于，所述步骤4包括如下步骤：
步骤4.1：读取异型坯钢种信息、连铸机结构参数和连铸工艺参数，计算固相线温度与液相线温度；
步骤4.2：导入1/2异型坯二维有限元模型的网格信息；
步骤4.3：确定时间步长；所述时间步长为每次计算切片移动的时间；所述切片代表的是所述1/2异型坯二维有限元模型；
步骤4.4：判断单元相区，计算单元的物性参数；
步骤4.5：计算单元传热矩阵、节点载荷、单元热容矩阵并集成总传热矩阵、总温度载荷、总热容矩阵；所述集成总传热矩阵就是将各个单元传热矩阵分解为子矩阵，然后将各个节点序列所决定的分块子矩阵按照单元节点编号“对号填充”进总传热矩阵中；所述集成总温度载荷、总热容矩阵均与所述集成总传热矩阵同理；
步骤4.6：求解温度场，也就是求解有限元二维非稳态热传导控制方程
Figure FDA0003489694830000031
同时判断切片位置是否超过空冷区，如果是则计算结束，如果否则切片继续移动，转到步骤4.4。
7.根据权利要求6所述的异型坯连铸凝固传热的有限元计算方法，其特征在于，所述钢种信息包括钢牌号、化学成分；所述连铸机结构参数包括结晶器高度、二冷区结构参数；所述二冷区结构参数包括二冷区数量、长度、入口位置和出口位置；所述连铸工艺参数包括浇铸温度、拉速、结晶器热流密度、二冷各区水流量及水温、环境温度和异型坯断面尺寸；所述异型坯断面尺寸包括翼梢厚度、宽面宽度、窄面宽度、腹板长度、腹板厚度以及內缘弧半径、翼缘半径、圆角半径。
8.根据权利要求6所述的异型坯连铸凝固传热的有限元计算方法，其特征在于，所述相区是指液相区、固相区和固液两相区；所述物性参数包括单元密度、单元比热、单元固相率和单元热传导系数。
9.根据权利要求8所述的异型坯连铸凝固传热的有限元计算方法，其特征在于，所述判断单元相区的方法为：首先通过反距离加权插值法计算所述平面四节点四边形单元的中心温度，以及通过二维线性插值法计算三角形单元的中心温度；最后根据单元中心温度判断单元所处相区。
10.根据权利要求6所述的异型坯连铸凝固传热的有限元计算方法，其特征在于，采用高斯-赛德尔迭代法求解矩阵方程
Figure FDA0003489694830000032
计算各节点温度。
CN202210092730.3A 2022-01-26 2022-01-26 一种异型坯连铸凝固传热的有限元计算方法 Pending CN114417675A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210092730.3A CN114417675A (zh) 2022-01-26 2022-01-26 一种异型坯连铸凝固传热的有限元计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210092730.3A CN114417675A (zh) 2022-01-26 2022-01-26 一种异型坯连铸凝固传热的有限元计算方法

Publications (1)

Publication Number Publication Date
CN114417675A true CN114417675A (zh) 2022-04-29

Family

ID=81276782

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210092730.3A Pending CN114417675A (zh) 2022-01-26 2022-01-26 一种异型坯连铸凝固传热的有限元计算方法

Country Status (1)

Country Link
CN (1) CN114417675A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110941889A (zh) * 2019-08-29 2020-03-31 华北理工大学 连铸异型坯微观和宏观裂纹萌发和扩展的研究方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110941889A (zh) * 2019-08-29 2020-03-31 华北理工大学 连铸异型坯微观和宏观裂纹萌发和扩展的研究方法
CN110941889B (zh) * 2019-08-29 2024-03-22 华北理工大学 连铸异型坯微观和宏观裂纹萌发和扩展的研究方法

Similar Documents

Publication Publication Date Title
Guo et al. Numerical simulation for the tip leakage vortex cavitation
Hefny et al. CFD analysis of pollutant dispersion around buildings: Effect of cell geometry
Davidson Effect of inclined vortex generators on heat transfer enhancement in a three-dimensional channel
Yanke et al. Simulation of slag-skin formation in electroslag remelting using a volume-of-fluid method
Ikeda et al. Single-phase CFD applicability for estimating fluid hot-spot locations in a 5× 5 fuel rod bundle
Chennu et al. Development of heat transfer coefficient and friction factor correlations for offset fins using CFD
CN114417675A (zh) 一种异型坯连铸凝固传热的有限元计算方法
Piper et al. A CFD study of the thermo-hydraulic characteristics of pillow-plate heat exchangers
Eldeeb et al. Pillow plate heat exchanger weld shape optimization using approximation and parallel parameterized CFD and non-uniform rational B-splines
Chatterjee et al. Unsteady mixed convection heat transfer from tandem square cylinders in cross flow at low Reynolds numbers
Zhang et al. Numerical analysis of thermal-hydraulic characteristics on serrated fins with different attack angles and wavelength to fin length ratio
Petrovic et al. Numerical and experimental performance investigation of a heat exchanger designed using topologically optimized fins
Liu et al. Lift-drag characteristics of S-shaped hydrofoil under different cloud cavitation conditions
CN116052806A (zh) 一种不规则断面连铸坯凝固传热的有限体积计算方法
CN111931392B (zh) 等离子加热中间包底吹氩参数优化方法和装置
CN114491855A (zh) 一种圆坯连铸凝固传热的有限元计算方法
CN112983560A (zh) 一种涡轮气冷叶片温度场快速评估方法
CN104573206B (zh) 一种基于有限元热力耦合的型钢断面热形状尺寸设计方法
CN115455760A (zh) 铸坯在线热状态跟踪方法
CN109522677A (zh) 一种用于热轧带钢温度控制的带钢横断面分层计算的方法
CN114491692A (zh) 变截面连铸坯的温度分布及固相率分布计算方法及系统
Li et al. Influences of casting speed and sen depth on fluid flow in the funnel type mold of a thin slab caster
CN111128316A (zh) 一种具有直裂纹或异质拼接材料的热性能分析方法
CN114491926B (zh) 一种粗糙换热面局部钙类污垢沉积分析方法
Xie et al. Numerical analysis of coupled fluid flow, heat transfer and solidification in ultra-thick slab continuous casting mold

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