高温模拟下冰块的热传递和融化过程

一、科学问题提出  

材料失效，这种现象在物理世界中普遍存在，力学范涛中主要代表物质失去承载能力的现象。进一步解释材料失效，是指材料或结构在外部作用下（如应力、温度、化学腐蚀等）无法继续履行其设计功能的现象。  

失效包括但不限于强度失效、疲劳失效、蠕变失效、腐蚀失效、磨损失效和应力腐蚀失效等。失效与破坏的区别在于失效是功能上的损失，而破坏是物理上的断裂。失效的形式可以是弹性变形失效、塑性变形失效、韧性断裂失效、脆性断裂失效、疲劳断裂失效和腐蚀失效。  

材料失效的原因可能包括设计不当、材料缺陷、铸造缺陷、焊接缺陷和热处理缺陷等。失效机理通常分为断裂、变形、磨损和腐蚀四种类型。影响材料失效的因素包括微观结构和宏观形态的变化，以及材料的力学性能、化学性能和环境因素。  

ABAQUS提供了很多材料模型进行失效模拟；还提供了Umeshmotion子程序和ALE技术，为磨损等材料侵蚀的仿真创造了可能。  

ABAQUS提供的材料损伤和失效模型如下：  

Shear、Ductile、JC Damage: 用于模拟金属材料剪切和延性损伤；  

FLD、FLSD、M-K、MSFLD Damage: 用于模拟金属薄壁件的缩颈破坏；  

Quade、Quads、Maxs、Maxe、Maxps、Maxpe Damage: 用于界面单元模拟的损伤破坏；  

Hashin Damage: 用于纤维增强复合材料的损伤模拟；  

Mullins Effect: 用于超弹性材料的损伤模拟；  

CDP：混凝土塑性损伤模型。  

同时用户还可以根据自己的需求，自定义材料子程序（UMAT/VUMAT），设计个性化的损伤材料模型。  

在研究生课题组中，有项目方向是做液体流动迁移扩散仿真，故考虑ABAQUS是否能够进行液体性质模拟。  

本作业主要应用Umeshmotion子程序模拟高温下冰块的热传递和融化过程。在此之前，我们在课堂上学习过支架的线性静力分析、压力容器内压静力分析、含切口板材单轴拉伸模拟、罐与接管的热分析、基体上薄膜脱粘分析等，结合这些基础，通过设定材料属性，使用ALE自适应网格控制，调用Umeshmotion子程序，来模拟高温下冰块的热传递和融化过程。ABAQUS的Umeshmotion利用自适应网格技术在计算过程中自动调整节点位置，由此可进行烧蚀、磨损等涉及节点移动的模型仿真。  

在ABAQUS中利用此可进行以下探究（本文仅进行案例复刻及一些改变）：  

热传递机制的模拟：在高温环境下，模拟冰块内部的热传递机制，包括传导、对流和辐射。ABAQUS提供了多种材料模型和边界条件来模拟这些热传递过程。例如，通过定义材料的热导率、比热容和热膨胀系数，以及设置对流换热系数和辐射参数，可以模拟冰块在高温环境中的热响应。  

融化动力学的探究：通过ABAQUS模拟冰块在高温条件下的融化速度和形态变化。ABAQUS的Umeshmotion子程序可以用来模拟冰块融化过程中体积的不断减少，这一仿真技巧也可以拓展应用到磨损、烧蚀、腐蚀等一系列涉及材料外形变化的仿真。  

温度分布的分析：利用ABAQUS模拟冰块在不同温度梯度下的内部温度分布。通过设置初始温度条件和对流换热系数，可以研究薄膜内的温度场，为后续的应力分析提供基础。  

物理性质变化的评估：在ABAQUS中模拟冰块融化过程中物理性质（如密度、热导率）的变化。这些性质的变化对热传递和融化过程有重要影响，可以通过定义温度-属性关系表来进行模拟。  

环境影响的考虑：研究环境因素（如压力、气流）对冰块热传递和融化过程的影响。ABAQUS允许设置复杂的边界条件和初始条件，以模拟实际环境中的热传递情况。  

热应力分析：由于温度场的变化会导致结构的热膨胀或收缩，从而产生热应力。ABAQUS可以进行热应力分析，考虑温度场如何随时间变化，并在分析中设置适当的边界条件。  

多物理场耦合分析：ABAQUS支持热-力耦合分析，可以同时考虑热传递和结构响应。可进行评估冰块在高温下的力学行为和稳定性。  

二、问题描述  

为冰块融化仿真，模拟高温下，冰块的热传递和融化过程。通过仿真可以得到冰块的体积不断减少的过程。主要用到Abaqus提供的Umeshmotion子程序。  

材料属性  

三、子程序Umeshmotion  

ABAQUS的子程序主要是使用Fortran语言编写的，IVF（Intel Visual Fortran）是常用的Fortran编译器之一。本文内容所用的程序如下：  

 SUBROUTINE Umeshmotion(UREF,ULOCAL,NODE,NNDOF,  

 *   LNODETYPE,ALOCAL,NDIM,TIME,DTIME,PNEWDT,  

 *   KSTEP,KINC,KMESHSWEEP,JMATYP,JGVBLOCK,LSMOOTH)  

C  

 INCLUDE 'ABA_PARAM.INC'  

C  

 DIMENSION ULOCAL(NDIM),JELEMLIST(100)  

 DIMENSION ALOCAL(NDIM,*),TIME(2)  

 DIMENSION JMATYP(*),JGVBLOCK(*)  

 DIMENSION ARRAY(1000)  

C LTRN指定输出坐标系，0为全局坐标系，1为局部坐标系  

 LTRN=0  

 CALL GETVRN(NODE,'NT',ARRAY,JRCD,JGVBLOCK,LTRN)  

 LSMOOTH = 1  

ULOCAL(NDIM)=ULOCAL(NDIM)-0.05*ARRAY(1)**2/10000.0  

 IF(LNODETYPE .EQ. 4)ULOCAL(NDIM) = ULOCAL(NDIM)*1.35  

 RETURN  

 END  

该模型采用的节点移动准则是节点随温度变化，在倒数第四行代码，节点随着正比于温度平方的速度沿着全局坐标系3方向进行移动。  

该程序虽短，但也具备了完备的Umeshmotion的功能，要点在于构造ULOCAL变量的变化关系。  

Umeshmotion子程序变量解释：  

1. ULOCAL自适应网格约束节点的网格位移或速度的组成部分，这些组成部分都在坐标系ALOCAL中已有描述。ULOCAL将作为由网格平滑算法确定的值传递到程序中。网格位移或速度的所有组件将被应用; 即您无法选择应用网格位移的方向。  

2. UREF作为自适应网格约束定义的一部分提供的、用户指定的位移或速度的值。该值基于任何幅度定义来更新，并用于与当前步骤相关联的自适应网格约束或默认斜坡幅度变量。  

3. NODE节点编号。  

4. NNDOF节点自由度数。  

5. LNODETYPE节点类型标志  

l LNODETYPE=1表明节点是位于自适应网格区域的内部；  

l LNODETYPE=2表明节点涉及绑定约束；  

l LNODETYPE=3表明节点位于自适应网格区域边界的拐角处；  

l LNODETYPE=4表明节点位于自适应网格区域边界的边缘；  

l LNODETYPE=5表明节点位于自适应网格区域边界的平面上；  

l LNODETYPE=6表明节点作为主节点参与约束（不是绑定约束）；  

l LNODETYPE=7表明节点作为从节点参与约束（不是绑定约束）；  

l LNODETYPE=10表明节点上作用了一个集中力。  

6. ALOCAL本地坐标系与节点处的自适应网格域的切线对齐。如果节点位于自适应网格域的内部，则ALOCAL将被设置为单位矩阵。在其他情况下，1方向是沿着平坦表面的边缘或平面。当NDIM=2时，2方向与表面正交。当NDIM=3时，2方向也位于平面的平面上，如果节点在边缘上，则该方向是任意的。当NDIM=3时，3方向垂直于表面，或者如果节点处于边缘，则为任意的。  

7. NDIM坐标维数。  

8. TIME(1) 当前分析步时间的值。  

9. TIME(2) 当前总时间的值。  

10. DTIME时间增量。  

11. KSTEP分析步编号。  

12. KINC增量步编号。  

13. KMESHSWEEP网格扫略数量。  

14. JMATYP必须传递到GETVRMAVGATNODE实用程序中以访问节点上的本地结果的变量。  

15. JGVBLOCK必须传递到GETVRN，GETNODETOELEMCONN和GETVRMAVGATNODE实用程序才能访问节点上的本地结果的变量。  

关于utility subroutine等详细信息可参考官方帮助文档。  

四、建模前置  

和以往课程里教学不同一点在于，本案例需要使用到ABAQUS的子程序调用，接下来简单介绍ABAQUS子程序调试。  

ABAQUS为用户提供了丰富而又灵活的用户子程序接口（USER SUBROUTINE），使得用户能够更灵活地解决一些问题，同时拓展了ABAQUS的功能。然而仅仅安装ABAQUS软件并不能直接使用到用户子程序的接口，需要关联两个软件Visual Studio（VS）和Intel visual Fortran（IVF）。  

1.安装前置软件。安装Visual Studio 2015 和 Intel visual Fortran 2016及ABAQUS软件。（需要先安装Visual Studio 2015 再安装 Intel visual Fortran2016 ，因为fortran运行器不能单独运行，需要安装在VS平台上，且这样IVF就能自动加载到VS的环境中，Intel visual Fortran的版本需要比Visual Studio 版本更高，ABAQUS的版本也需要适配，如下图，但安装顺序无要求，本内容的各应用版本为：Visual Studio 2015 和Intel visual Fortran 2016及ABAQUS 2024，仅供参考。在这篇文章中有更详细的介绍：https://mp.weixin.qq.com/s/QS70_OxhT-4C95um51RdRg）  

图常用的VS与IVF版本匹配  

PS：我在安装软件过程中遇到的问题：  

VS 2015安装包损坏或丢失，通过网上寻找解决方案，试过很多个并没有成功，最终在CSDN中这篇文章里第四个方法解决。（https://blog.csdn.net/hs_2017112123/article/details/122107693）。即“将你的安装包所在的路径复 制到安装失败的界面那个搜索包，例如我的packages包在此路径下：D:\vs2015\vs2015.pro_chs\packages，复 制然后贴到搜索包，丢包后可以重复复 制，继续复 制到安装完成即可。安装成功”。  

2.三软件关联。安装好软件之后，就是进行这三款软件的关联。使用Everything软件搜索Launcher.bat文件，注意需要是Abaqus下的这个文件，不是其他软件的。打开文件路径，使用记事本打开文件，进行修改Launcher.bat文件。  

首先打开是只有以下两行内容的：  

@echo off  

call "F:\Abaqus2024\SIMULIA\Commands\abq2024.bat" %*  

在这两行内容前添加vcvarsall.bat和ifortvars.bat文件的路径（使用Everything软件搜索，这两个文件在电脑中是唯一的），即添加以下两行内容，中间的路径依据个人修改，建议是复 制自己的路径。  

@call "F:\VS2016\VC\vcvarsall.bat" x64  

@call "F:\IVF\compilers_and_libraries_2016.1.146\windows\bin\ifortvars.bat" intel64 vs2015  

保存退出即可。  

3.验证子程序。在计算机“开始”界面找到“Dassault Systemes SIMULIA Established Products 2024”里的“Abaqus Verification”，右键找到文件所在位置，以管理员方式打开，即可进行验证。  

验证完毕后出现的日志文件中各项指标都是Pass即代表成功了。  

PS：我在此过程中遇到的问题：  

在搜寻vcvarsall.bat文件时，发现Unable to find vcvarsall.bat问题，这可能是下载vs时没有插件组选择上有缺漏，也可能是因为没有选C++桌面平台开发。网上的各种方法尝试并没有成功，找不到网上方法所对应的工具项，但我选择了对Visual Studio 2015 进行了重安装，选择修复，并勾选上VC++模块，才解决了这个问题。  

以上内容操作我是参考此视频讲解（https://www.bilibili.com/video/BV1tp4y187Aj/?spm_id_from=333.1387.favlist.content.click&vd_source=c8cf508bc4d5ce513c1f0cd4846a2bdd）和网上各种资料。  

五、建模步骤

1. 创建部件

图1 采用 Standard/Explicit模型

图2 创建部件Part-1

图3 绘制圆心(0，0)，半径50的圆

图4 设置深度100

2. 创建材料和截面属性

图5 编辑材料属性-传导率为1；密度为1e-9kg/mm3；杨氏模量1000；泊松比0.3；基本温度为20℃下膨胀系数1e-6；比热1e9

图6 创建截面Section-1，类别为“实体”、类型为“均质”，应用于材料Material-1

指派截面

图7 选择整个模块进行指派截面

图8 指派截面后

在“工具”-“set”-“管理器”中创建“集”.

图9 创建设置集

3. 定义装配件

图10 创建实例

4. 设置分析步

图11 创建分析步-选择温度-位移耦合

图12 编辑分析步-设置时间长度为500，几何非线性“开”

图13 设置增量步，最大增量步数设置为10000，增量步大小为10

图14 历程输出请求编辑-作用域为“集”，输出变量为VOLC，仅由自适应网格划分引起的面积或体积的改变

5. 调用子程序

图15 调用子程序，创建ALE自适应网格区域

图16 对分析步Step-1进行ALE自适应网格控制编辑

设置频率：1，对每个增量步重划扫掠网格：4

每4个增量步重划扫掠一次网格

图17 创建ALE自适应网格约束-选择“速度/角速度”

图18 选择Part-1-1，Set-all

图19 选择使用用户定义的子程序，即后续用到的Umeshmotion子程序

6. 设置载荷与边界条件

图20 创建边界条件-对位移固定

图21 选择对底面进行固定

图22 创建边界条件BC-2，选择类别为温度

图23 仍是对底面进行设置边界条件

图24 编辑好边界条件后

7. 划分网格

图25 设置网格控制属性，算法选用中性轴算法

图26 设置全局种子为10

图27 为部件划分网格

图28 指派单元类型，选择温度-位移耦合

8. 定义相互作用条件

图29 创建相互作用，选择表面热交换条件

图30 选择作用区域为除底面外区域

图31 设置膜层散热系数为0.2，环境温度为100

图32 已设置好的表面热交换条件相互作用

9. 提交作业

图 33 创建作业Job-1

图34 调用用户子程序 ice.for

六、结果分析

图35 作业运行监视

图36 S,Mises应力的0步长

图37 S,Mises应力的第10、100、200、300、400、500步长

图38 GRADT,Magnitude 最大处可达2.188e+01

图39 HFL,Magnitude 最大处可达2.188e+01

图40 HFL,在x,y,z方向

图41 LE Max,Principal，LE Mid,Principal，LE Min,Principal

图42 NT11 最大处达95.73

图43 RF,Magnitude 接触力变化

图44 RF,RF1,2,3 三个主方向上的接触力变化，说明在y轴方向上无

图45 U,Magnitude 位移变化

图46 U,U1,U2,U3 x,y,z三个主方向上的位移自由度

图47 冰块体积随时间的变化

七、总结

通过复刻案例，验证了ABAQUS是能够进行液体性质模拟的。本作业主要应用Umeshmotion子程序模拟高温下冰块的热传递和融化过程。使用ALE自适应网格控制，调用Umeshmotion子程序，来模拟高温下冰块的热传递和融化过程。得到应力、位移等性质变化云图，并且可通过可视化-动画-时间历程来动态的查看变化过程，后续将附以更加实际、复杂的条件，真实的模拟出冰块的融化过程情况，并进行分析。还有很多内容我还并不了解和理解，期待后续能够进一步学习使用有限元分析软件，用于项目课题中。