高溫模擬下冰塊的熱傳遞和融化過程

一、科學問題提出

材料失效，這種現象在物理世界中普遍存在，力學范濤中主要代表物質失去承載能力的現象。進一步解釋材料失效，是指材料或結構在外部作用下（如應力、溫度、化學腐蝕等）無法繼續履行其設計功能的現象。

失效包括但不限于強度失效、疲勞失效、蠕變失效、腐蝕失效、磨損失效和應力腐蝕失效等。失效與破壞的區別在于失效是功能上的損失，而破壞是物理上的斷裂。失效的形式可以是彈性變形失效、塑性變形失效、韌性斷裂失效、脆性斷裂失效、疲勞斷裂失效和腐蝕失效。

材料失效的原因可能包括設計不當、材料缺陷、鑄造缺陷、焊接缺陷和熱處理缺陷等。失效機理通常分為斷裂、變形、磨損和腐蝕四種類型。影響材料失效的因素包括微觀結構和宏觀形態的變化，以及材料的力學性能、化學性能和環境因素。

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可以采用Latent heat方法模擬相變材料，定義相變溫度、相變潛熱和相變區間。

多物理場耦合分析：ABAQUS支持熱-力耦合分析，可以同時考慮熱傳遞和結構響應。可進行評估冰塊在高溫下的力學行為和穩定性。

二、問題描述

為冰塊融化仿真，模擬高溫下，冰塊的熱傳遞和融化過程。通過仿真可以得到冰塊的體積不斷減少的過程。主要用到Abaqus提供的Umeshmotion子程序。

材料屬性

屬性名稱	大小
傳導率	1
密度	1e-9kg/mm3
楊氏模量	1000
泊松比	0.3
基本溫度	20℃
膨脹系數	1e-6
比熱	1e9
膜層散熱系數	0.2
環境溫度	100℃

三、子程序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節點類型標志

? LNODETYPE=1表明節點是位于自適應網格區域的內部；

? LNODETYPE=2表明節點涉及綁定約束；

? LNODETYPE=3表明節點位于自適應網格區域邊界的拐角處；

? LNODETYPE=4表明節點位于自適應網格區域邊界的邊緣；

? LNODETYPE=5表明節點位于自適應網格區域邊界的平面上；

? LNODETYPE=6表明節點作為主節點參與約束（不是綁定約束）；

? LNODETYPE=7表明節點作為從節點參與約束（不是綁定約束）；

? 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子程序，來模擬高溫下冰塊的熱傳遞和融化過程。得到應力、位移等性質變化云圖，并且可通過可視化-動畫-時間歷程來動態的查看變化過程，后續將附以更加實際、復雜的條件，真實的模擬出冰塊的融化過程情況，并進行分析。還有很多內容我還并不了解和理解，期待后續能夠進一步學習使用有限元分析軟件，用于項目課題中。