HETVAL子程序？？

 

混凝土水化熱溫度場分析其實是相當于在混凝土的溫度場分析中加入了一個熱源，而這個熱源的放熱量是隨著時間變化的。由于在Abaqus中沒有直接功能來模擬隨著時間變化的熱源，所以需要借用HETVAL子程序來實現隨著時間變化的熱源功能，并將其耦合到混凝土溫度場的計算之中。

 

見圖1，HETVAL子程序用來提供傳熱分析模型的熱源，這個熱源是隨著時間變化的，這個與混凝土隨時間變化的水化放熱曲線是一致的（圖2）。另外，還可以將該熱源在HETVAL子程序中定義為受結果狀態變量影響而變化。這個特別重要，因為混凝土水化放熱的曲線隨著溫度的增加會放熱加快，而需要模擬這一現象就要考慮溫度場計算的溫度結果對混凝土水化放熱的影響。

 

圖1 abaqus幫助手冊關于HETVAL子程序的解釋（http://wufengyun.com:888/books/sub/default.htm）

 

圖2 混凝土水化放熱速率曲線 （圖片來自Fairbairn, Eduardo M. R. , and M. Azenha . "[RILEM State-of-the-Art Reports] Thermal Cracking of Massive Concrete Structures）

 

HETVAL的主體程序部分如圖3，其中子程序調用的形參解釋如下：

 

圖3 HETVAL主體子程序

 

---

 

 

01

 

 

子程序傳遞的變量

 

• CMNAME: 用戶自定義的材料名字。

• TEMP為含2個元素的數組，TEMP(1)為當前溫度，TEMP(2)為溫度增量。

• TIME為含2個元素的數組，TIME(1)為增量結束時的時間步長，TIME(2)為增量結束時的總時間。

• DTIME為時間步長。

• PREDEF(*)為包含所有用戶自定義場變量的數組。

• DPRED(*)為預定義場變量的增量數組。

 

---

 

 

02

 

子程序需要定義的變量

 

而FLUX(1)及其重要，間接表示定義的混凝土熱源，如圖4所示。FLUX(1)表示產生的熱流量，r，也就是單位時間單位體積的產熱量，單位為J/T/L^3，即為焦耳/(秒*立方米)。

 

圖4 FLUX(1)表示產生的熱流量

 

---

 

 

03

 

 

子程序需要更新的變量

 

STATEV(*)：用戶自定義的結果依賴的狀態變量數組。僅在熱傳遞分析中，STATEV在增量步開始時會傳入子程序，并且在增量步結束時傳遞回來，即更新狀態變量。

 

圖5 依賴于結果的STATEV狀態變量

 

---

 

 

04

 

 

水化放熱公式及實現程序

 

下面通過一個例子來探索下HETVAL子程序模擬混凝土水化熱形成的溫度場。

 

水化放熱曲線很復雜，工程中一般采用簡化的經驗擬合表達式，比如混凝土水化放熱的公式我們這里采用復合指數式，如式1所示，因為很多文獻指出用該水化放熱曲線模擬的溫度場與實測的溫度場模擬較好。Q0為混凝土水化放熱的最大值，取為364000J/kg（一般試驗或文獻里給的都是kJ/kg這種單位，但Abaqus里的標準單位為J，這里要統一單位），而b、c這兩個參數取值可在文獻中查到，不同水泥類型的參數取值不一樣，對于52.5普通硅酸鹽水泥，這里b取為0.36，c為0.74。注意到式1中的t的單位為天。

 

式1：

然后對式1求導，可得到單位立方米混凝土單位時間的水化放熱率，qv (J/m^3/hr)，見式2。這里將上式中的時間單位改為了hr，因為一般混凝土溫度場監測以小時為單位進行記錄。這里的qv就是HETVAL中的FLUX(1)。具體HETVAL子程序如下：

 

式2：

 

在該子程序中最重要為熱源放熱率FLUX(1)，并且及時用statev結果狀態變量來存儲，這樣方便在后處理中輸出查看，我稱之為實時調試，也就是可以實時監控自己定義或想要輸出變量在有限元計算中隨時間的變化規律。這里statev(1)記錄的是熱源放熱率flux(1)（時間單位為小時），statev(2)記錄的是時變的水化放熱總量（注意到時間單位也換為了小時，因為模型中的混凝土熱工性能參數，比如熱傳導系數用的時間單位也是小時，這里統一對應起來），而statev(3)、statev(4)和statev(5)分別為Q0，b，c，這樣可以驗證下子程序的正確性，看下輸出這些確定值是否在計算過程中發生了變化。

 

另外，特別需要注意的是，下列程序中的Cemass這個變量一定要出現在FLUX(1)的計算表達式中，否則計算出來的溫度變化極其小。這并不是說模型是有錯的，而是因為Cemass表示每立方米混凝土中的水泥質量（kg/m3），因為HETVAL子程序中出現的FLUX(1)表示單位體積單位時間的水化熱源放熱量，強調是單位體積，也就是單位立方米混凝土，所以這里的水泥放熱量應該是單位立方米混凝土中所有水泥的累加水泥放熱量，否則如果不乘以Cemass，則計算的單位質量水泥（kg）的水化放熱量，因此導致計算的溫度結果幾乎不變化。

      SUBROUTINE HETVAL(CMNAME,TEMP,TIME,DTIME,STATEV,FLUX,     1 PREDEF,DPRED)C      INCLUDE 'ABA_PARAM.INC'C      CHARACTER*80 CMNAMEC      DIMENSION TEMP(2),STATEV(20),PREDEF(*),TIME(2),FLUX(2),     1 DPRED(*)      PARAMETER(b=0.36D0,c=0.74D0,Q0=364000D0,Cemass=500)     C     qv(t)      FLUX(1)=Cemass*Q0*b*c/24*(Time(2)/24)**(C-1)*exp(-b*(Time(2)/24)**c)   c      FLUX(1)=10000000000 C     qv(t)      STATEV(1)=FLUX(1)C     Q(t)      STATEV(2)=Cemass*Q0*(1.0D0-exp(-b*(Time(2)/24)**c))C     Q0      STATEV(3)=Q0C     b      STATEV(4)=bC     c      STATEV(5)=c      RETURN      END

 

 

---

 

 

05

 

模型材料設置

這里我采用一個簡單的混凝土長方體來模擬水化熱溫度場，其中關于水化放熱的設置見圖6。heat generation這里一定要點選，否則模型不會考慮HETVAL子程序。同時，user material中需要定義結果狀態變量，這里我定義20個，事實上子程序中需要輸出的狀態變量僅為5個，但這樣設置一是不影響，而是方便后續額外增加新的狀態變量，利于觀察調試模型。

 

圖6 Heat Generation設置

 

---

 

 

06

 

水化熱溫度場計算結果

 

圖7為混凝土水化熱溫度場在11.6小時的溫度等值線圖，這里的NT11表示節點溫度（節點溫度也就一個自由度）?？梢钥闯?，等值線圖幾乎一樣均勻，這是由于整個混凝土塊與外界保持絕熱狀態，每個混凝土單元都在放熱，也就是每個混凝土單元之間不存在溫度梯度（圖8熱流量密度HFL等值線圖也可證明，可看出熱流量密度HFL很?。?。而圖9表示了混凝土水化熱溫度場計算結果隨時間的變化曲線，可以看出來混凝土一直在升溫，也就沒有熱傳導，相當于均勻升溫。圖10的動畫用等值線圖的方式表示了混凝土水化熱溫度場隨著時間的變化，但實際上在每個時刻的溫度場等值線圖是一致的。

 

圖7 NT11溫度等值線圖

 

圖8 HFL熱流量密度等值線圖

 

圖9 模型任一單元的溫升時程圖

 

 

圖10 模型溫度場的時變動畫

 

---

 

 

07

 

水化放熱速率及累計放熱量計算結果

圖11可發現SDV1（混凝土水化放熱速率，FLUX(1)）的等值線圖，發現等值線都是一致均勻的，這是因為每個混凝土單元都在采用同樣的水化放熱速率曲線，水化放熱速率都是一致的。另外，利用XYData和XYPlots，選擇圖12所示的某一混凝土單元，觀察SDV1隨時間的變化規律。圖13中也繪制了混凝土水化放速率隨著時間的變化曲線。同時圖14也繪制了混凝土累計水化放熱量隨著時間的變化曲線。因此，子程序在計算過程中的正確性得以保證?？梢钥闯觯艧崴俾视?突然增大到最大值，然后逐漸減少，在24小時后放熱速率變得較小且趨于恒定，但仍舊在放熱（之前已經釋放掉大部分的熱量）。

 

圖11 SDV1（FLUX1(1)）等值線圖

 

圖12 利用XYData繪制SDV1（FLUX1(1)）隨時間的變化曲線

 

圖13 SDV2（混凝土累計水化放熱量）等值線圖

 

而基于選中單元的中心（與SDV1一致），繪制SDV2(混凝土累計水化放熱量)隨時間的變化規律如圖14所示。

 

圖14 利用XYData繪制SDV2（混凝土累計水化放熱量）隨時間的變化曲線

 

最后，有需要歡迎通過微信公眾號聯系我們。

微信公眾號：320科技工作室。