高温模拟下冰块的热传递和融化过程有限元分析

  问题描述  

为冰块融化仿真，模拟高温下，冰块的热传递和融化过程。通过仿真可以得到冰块的体积不断减少的过程。主要用到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节点类型标志

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

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

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

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

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

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

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

lLNODETYPE=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子程序，来模拟高温下冰块的热传递和融化过程。得到应力、位移等性质变化云图，并且可通过可视化-动画-时间历程来动态的查看变化过程，后续将附以更加实际、复杂的条件，真实的模拟出冰块的融化过程情况，并进行分析。还有很多内容我还并不了解和理解，期待后续能够进一步学习使用有限元分析软件，用于项目课题中。