1.簡介

　　熱膨脹效應：

　　(1)可以通過熱膨脹系數(thermal expansioncoefficient)來定義，ABAQUS可以通過熱膨脹系數來計算熱應變;

　　(2)熱膨脹效應可以是各向同性、正交異性或各向異性的;

　　(3)被定義為從參考溫度(reference temperature)到計算溫度產生的膨脹效應;

　　(4)可以被定義為一個依賴于溫度和/或場變量的函數;

　　(5)在ABAQUS/Standard中，如果是連續實體單元，可以被定義為一個分布函數;

　　(6)在ABAQUS/Standard中可以直接使用用戶子程序UEXPAN來定義(如果熱應變是一個關于場變量和狀態變量的復雜函數)。

　　2.計算熱應變

　　ABAQUS通過熱膨脹系數α來定義從參考溫度θ0到計算溫度θ產生的總的熱應變，如圖2.1所示。

　　圖2.1 熱膨脹系數的定義

　　產生的熱應變可以通過下式進行計算：

　　其中：

　　α(θ,fβ)為熱膨脹系數;

　　θ為當前溫度;

　　θI為初始溫度;

　　fβ為場變量當前的取值;

　　fβI為場變量的初始值;

　　θ0為熱膨脹系數的參考溫度;

　　上式中的第二項代表由于初始溫度θI和參考溫度θ0之間的溫度差產生的熱應變。這一項能夠保證當參考溫度不等于初始溫度時初始熱應變為0。

　　當熱膨脹系數不為關于溫度和場變量的函數時，可以不定義參考溫度θ0。

　　3.熱應變的轉換

　　在ABAQUS中通常是以表格的形式來輸入平均熱膨脹系數。但有的時候，提供給你的熱膨脹數據可能是以微分的形式來表示的：

　　即提供的是應變-溫度曲線(如圖2.1所示)的斜率。為了將這些數據轉換為ABAQUS支持的平均熱膨脹系數，需要進行一個從參考溫度θ0到當前溫度θ的積分運算：

　　例如，假設為一組常數值：dα1/dθ位于θ0和θ1之間;dα2/dθ位于θ1和θ2之間;dα3/dθ位于θ2和θ3之間;則有：

　　因此在ABAQUS中需要輸入的平均熱膨脹系數可表示為：

　　通過在ABAQUS/Standard中使用用戶子程序UEXPAN可以定義與溫度和/或場變量有關的熱應變增量;如果熱應變增量是與狀態變量有關的函數，則必須使用用戶子程序UEXPAN。

　　如果材料的熱膨脹行為是各向同性的，則在用戶子程序UEXPAN中只需要定義一個各向同性熱應變增量(Δε=Δε11=Δε22=Δε33);如果是正交異性的，則需要定義熱應變增量在三個主方向上的分量(Δε11,Δε22,Δε33);如果是各向異性的，則需要定義熱應變增量的六個分量(Δε11,Δε22,Δε33,Δε12,Δε13,Δε23)。

　　4.熱應力

　　如果一個結構不能自由膨脹，則溫度的改變將會在結構內部產生應力。例如，考慮一個長度為L的桿單元，桿的兩端被完全約束。橫截面面積，楊氏模量E和熱膨脹系數α均為常數。對于該一維問題，桿單元上產生的應力可以通過Hooke定律計算得到：

　　其中εx為總應變，εxth為熱應變，Δθ為溫度的該變量。由于桿單元是完全約束的，因此εx=0。如果桿單元兩端節點的溫度相同，則應力為：

　　受約束的熱膨脹可以產生很大的應力值。對于典型的結構鋼而言，150oC的溫度變化都將使得材料屈服。因此，在包含熱載荷的問題中，應該非常小心地處理邊界條件以避免對熱膨脹造成過約束。

　　5.用戶子程序UEXPAN

　　5.1簡介

　　用戶子程序UEXPAN：

　　(1)可以用于將熱應變增量定義為關于溫度、預定義場變量和狀態變量的函數;

　　(2)被用于建立熱應變與溫度和/或預定義場變量和狀態變量之間的復雜關系;可以在子程序中對熱應變進行更新;

　　(3)將會在單元的所有積分點上被調用;

　　(4)在熱力耦合分析中的每次迭代過程中，每個積分點會調用該子程序兩次。

　　5.2用戶子程序接口

　　SUBROUTINE UEXPAN(EXPAN,DEXPANDT,TEMP,TIME,DTIME,PREDEF,1 DPRED,STATEV,CMNAME,NSTATV,NOEL)CINCLUDE 'ABA_PARAM.INC'CCHARACTER*80 CMNAMECDIMENSION EXPAN(*),DEXPANDT(*),TEMP(2),TIME(2),PREDEF(*),1 DPRED(*),STATEV(NSTATV)user coding to define EXPAN, DEXPANDT and updateSTATEV if necessary.RETURNEND

　　5.3可以定義的變量

　　EXPAN (*)

　　熱應變的增量。熱應變分量的數量和順序取決于熱膨脹行為的類型。

　　對于各相同性熱膨脹，只需要定義一個分量;

　　對于正交異性熱膨脹，需要定義Δε11th,Δε22th,Δε33th三個分量;

　　對于各相異性熱膨脹，需要分別定義Δε11th,Δε22th,Δε33th, Δε12th,Δε13th,Δε23th六個分量;如果是平面應力，則只需要定義Δε11th,Δε22th,Δε33th三個分量。

　　DEXPANDT (*)

　　熱應變隨溫度的變化率，?εth/?θ。分量的數量和順序取決于熱膨脹行為的類型，與熱應變增量類似。

　　5.4可以更新的變量

　　STATEV (NSTATEV)

　　用戶自定義的狀態變量。除了熱力耦合分析，狀態變量將會在增量步的開始時間被傳遞到子程序中，并可以在增量步結束時被更新。對于熱力耦合分析，在每個積分點，每個增量步中子程序將被調用兩次。在第一次調用時，狀態變量的取值將會在增量步的初始時刻被傳遞到子程序中，并在增量步結束時被更新。在第二次調用時，第一次更新后的狀態變量值將被再次傳遞到子程序中，并在增量步的結束時刻可以再次更新。

　　用戶子程序UEXPAN允許熱應變增量和狀態變量是弱相關的。熱應變對狀態變量的導數不會包含在Jacobian矩陣中。

　　5.5用于獲取相關信息的變量

　　TEMP (1)

　　當前溫度(增量步的終止時刻)。

　　TEMP (2)

　　溫度增量。

　　TIME (1)

　　增量步終止時刻對應的分析步時間。

　　TIME (2)

　　增量步終止時刻對應的總時間。

　　DTIME

　　時間增量。

　　PREDEF (*)

　　包含有所有用戶自定義場變量的數組(分析開始時的初始值和當前的取值)。

　　CMNAME

　　用戶指定的材料名稱。

　　NSTATEV

　　與該材料有關的狀態變量的數量。

　　NOEL

　　用戶自定義單元號。

　　6.平均熱膨脹系數的計算

　　圖6.1為通過材料屬性模擬軟件JmatPro計算得到的SAE3140鋼在冷卻過程中奧氏體向馬氏體轉變的熱應變溫度曲線。

　　圖6.1 熱應變-溫度曲線

　　從圖6.1中可以看出隨著溫度的下降，材料產生熱應變也隨之降低，當溫度下降到約300oC時，熱應變出現突然增長的現象，這是由奧氏體向馬氏體轉變過程中的體積膨脹引起的。因此該突變點對應的溫度值即為馬氏體轉變開始溫度Ms。

　　采用節3中的方法計算平均熱膨脹系數：

　　其中ε0th為參考溫度θ0對應的熱應變，從圖6.1中可以看出取值為0，αi為溫度θi對應的平均熱膨脹系數。

　　在Matlab中編制平均熱膨脹系數的計算程序，如下所示。

　　%程序用于根據熱應變溫度曲線計算平均熱膨脹系數exp_data=importdata('thermal_temp_curve.txt'); %導入熱應變溫度數據temp_data=exp_data(:,1); %溫度數據strain_data=exp_data(:,2); %熱應變數據T_ref=1395; %參考溫度strain_ref=interp1(temp_data,strain_data,T_ref); %參考熱應變sample_num=500; %計算樣點T=linspace(1300,25,sample_num)';strain=interp1(temp_data,strain_data,T);alpha=zeros(sample_num,1);%計算平均熱膨脹系數for i=1:sample_num alpha(i)=(strain(i)-strain_ref)/(T(i)-T_ref);endplot(T,alpha)

　　通過程序計算得到不同參考溫度下以及JmatPro給出的平均熱膨脹系數如圖6.2所示。

　　圖6.2 平均熱膨脹系數隨溫度變化

　　從圖6.2中可以看出，當參考溫度θ0取為1395oC時，計算得到的平均熱膨脹系數隨溫度的變化與JmatPro給出的結果幾乎是重合的，這是因為JmatPro將應變為0的點(圖6.1)對應的溫度作為了參考溫度。隨著參考溫度的變化，平均熱膨脹系數也隨之改變。因此，為了在ABAQUS中正確地定義熱膨脹行為，必須準確地指定定義平均熱膨脹系數所采用的參考溫度，否則將無法獲得正確的熱應變值。

　　通過圖6.2可以看出，當參考溫度θ0取為某一值時，可以近似使得奧氏體和馬氏體的熱膨脹系數為一個不隨溫度變化的常數，如圖6.2中的綠色曲線所示，當奧氏體幾乎完全轉變為馬氏體之后，綠色曲線的末端近似為一條水平線，這代表著馬氏體的熱膨脹系數為一個常數值。而根據平均熱膨脹系數的定義式不難看出，該參考溫度的準確值可以通過如圖6.3所示的方法進行求解。

　　首先用直線分別擬合對應奧氏體和馬氏體體積分數為1時的熱應變-溫度曲線，兩條直線的交點即為參考溫度θ0。從圖6.3中可以看出參考溫度θ0約為881oC。

　　圖6.3 參考溫度的確定

　　兩條直線的斜率即為奧氏體和馬氏體的熱膨脹系數，其取值如表6.1所示。

　　表6.1 奧氏體和馬氏體的熱膨脹系數

　　利用上面計算得到的參考溫度重新計算平均熱膨脹系數隨溫度的變化，結果如圖6.4所示。

　　圖6.4 平均熱膨脹系數隨溫度的變化

　　圖6.4中兩條虛線分別代表奧氏體和馬氏體的體積分數為1時的熱膨脹系數，從圖中可以看出，在使用了圖6.3中確定的參考溫度后，奧氏體和馬氏體的熱膨脹系數基本不隨溫度變化。換句話說，通過合理地選取參考溫度，可以使得奧氏體和馬氏體的熱膨脹系數不隨溫度變化。

　　7.通過熱膨脹系數計算熱應變

　　通過上面分析過程可以看出，如果已知某一參考溫度θ0下奧氏體和馬氏體的平均熱膨脹系數為一個不隨溫度變化的常數，以及奧氏體向馬氏體轉變過程中奧氏體和馬氏體的體積分數(假定不存在擴散型相變)，則可以直接計算如圖6.1所示的熱應變-溫度曲線。

　　而奧氏體向馬氏體轉變的過程為非擴散型相變，可以通過Koistinen-Marburger模型來描述：

　　其中，vm和va分別為馬氏體和奧氏體的體積分數，b為與材料有關的常數，Ms為馬氏體轉變開始溫度。通過圖6.1中的拐點可以大致判斷出Ms約為300oC。

　　常數b可以通過馬氏體體積分數與溫度之間的關系來確定，這里直接給出通過JmatPro計算得到的結果，如表7.1所示。

　　由上式可得奧氏體體積分數和溫度之間滿足關系：

　　因此，將奧氏體體積分數的自然對數與溫度的數據進行線性擬合即可得到常數b的取值為0.019oC-1，Ms取為310.5oC，擬合曲線如圖7.1所示。

　　在已知任意溫度下奧氏體和馬氏體的體積分數后，則對應的熱膨脹系數可表示為：

　　熱應變可表示為：

　　其中θI為初始時刻的溫度，即熱應變為0對應的溫度值，在圖6.1中θI為1395oC。

　　圖7.1 擬合Koistinen-Marburger模型常數b

　　通過上述方法在Matlab中編寫計算程序，如下所示。

　　%程序用于計算奧氏體向馬氏體轉變過程中的熱應變隨溫度變化sample_num=200; %計算樣點數CTE_a=22.44e-6; %奧氏體熱膨脹系數CTE_m=15e-6; %馬氏體熱膨脹系數T_ini=1395; %初始溫度T_ref=881; %參考溫度T_Ms=310.5; %馬氏體轉變開始溫度KM_b=0.019; %Koistinen-Marburger模型系數bVa=ones(sample_num,1); %奧氏體體積分數Vm=zeros(sample_num,1); %馬氏體體積分數alpha=zeros(sample_num,1); %熱膨脹系數strain=zeros(sample_num,1); %熱應變T=linspace(1400,20,200)'; %溫度for i=1:sample_num %計算奧氏體和馬氏體體積分數 if T(i)<=T_Ms Vm(i)=1-exp(-1*KM_b*(T_Ms-T(i))); Va(i)=1-Vm(i); end %計算熱膨脹系數 alpha(i)=Vm(i)*CTE_m+Va(i)*CTE_a; %計算熱應變 strain(i)=alpha(i)*(T(i)-T_ref)-... alpha(i)*(T_ini-T_ref);endfigure;plot(T,Vm,T,Va);figure;plot(T,alpha);figure;plot(T,strain);

　　計算得到的熱應變隨溫度的變化如圖7.2所示。

　　圖7.2 熱應變-溫度變化曲線

　　圖7.2中黑色曲線為通過JmatPro計算得到的熱應變曲線，紅色曲線為馬氏體熱膨脹系數取值為10.46×10-6oC-1時計算得到的熱應變曲線，該熱膨脹系數值是通過擬合JmatPro給出的熱應變曲線得到的，但似乎與黑色曲線有較大的差異;反而是將馬氏體熱膨脹系數取值為15×10-6oC-1時與黑色曲線較為吻合，本文認為這主要是由于擬合直線斜率時溫度范圍選取不當引起的。在擬合馬氏體的熱膨脹系數時，本文選取了20oC~200oC的溫度范圍，而從表7.1中可以看出，當溫度位于200oC時，此時馬氏體的體積分數為90%，因此此時對應的熱膨脹系數并非馬氏體體積分數為100%時的熱膨脹系數。如果采用該溫度段進行擬合，將造成擬合得到的熱膨脹系數偏小。如果縮小擬合時選取的溫度范圍，預計能夠獲得較好的結果。皮特寵物博客

(素材來自網絡，侵刪)