如何在ABAQUS中计算热应变

11月前浏览872

为了探究在ABAQUS中热应变的计算原理，本文将首先建立单个单元的有限元模型，并利用单轴拉伸验证有限元模型的正确性。在单个单元有限元模型的基础上，考虑了平均热膨胀系数为常数以及随温度变化的两种情况，利用ABAQUS/CAE分别计算两种情况下的热应变变化，并与MATLAB的计算结果进行对比分析。随后，本文将根据前面的验证结果给出ABAQUS中热应变的计算方法，利用用户子程序UEXPAN编写热应变的计算程序，并与ABAQUS/CAE的计算结果进行对比，验证计算方法的正确性。

1. 利用单轴拉伸验证单个单元有限元模型

取杨氏模量E为2×10¹¹Pa，泊松比v为0.3，建立单个单元的有限元模型，并在三个单元面上分别施加对称约束，使该单元处于单轴拉伸状态，如图1.1所示。

图1.1 单个单元施加对称约束

在沿Y轴的方向上施加100MPa的拉伸载荷，可计算得到单元应力分布如图1.2所示。

图1.2 单个单元应力分布

由此可以计算得到在100MPa拉伸应力作用下，单元在Y方向产生的弹性应变为：

而在X方向和Z方向产生的弹性应变为：

这与ABAQUS的计算结果完全吻合，如1.3所示。

图1.3 单个单元在单轴拉伸下的应变E22

2. 利用ABAQUS/CAE计算热应变

2.1 平均热膨胀系数为常数值

在单轴拉伸的基础上，定义分析步thermo_expan，假定该单元初始温度为0ºC，利用Predefined Field将20ºC的温度值施加到整个单元上，则在该分析步中，单元温度将从0ºC线性增长至20ºC。假设单元的热膨胀是各向同性的，平均热膨胀系数α取为2×10^-5 ºC^-1，则在该分析步结束时在该单元上产生的热应变为：

此时单元上的应变为弹性应变和热应变的叠加，即：

因此在三个方向上的应变分量为：

在ABAQUS中分别输出三个主应变分量以及热应变（需要在场输出中勾选THE），得到其变化如图2.1所示。

图2.1 单元应变变化

图2.1中E为总的应变，THE为热应变，0到1s为施加拉伸载荷的阶段，此时热应变THE为0；1到2s时随着温度值从0ºC线性增长至20ºC，热应变值增加。最终计算得到的应变值与前面的计算结果是完全一致的。

2.2 平均热膨胀系数随温度变化

假定材料的平均热膨胀系数α随温度变化，不同温度下的取值如表2.1所示。

表1 不同温度下的平均热膨胀系数

在ABAQUS中，热应变通过下式进行计算¹：

其中：

α(θ,f_β)为平均热膨胀系数；

θ为当前温度；

θ^I为初始温度；

f_β为场变量当前的取值；

f_β^I为场变量的初始值；

θ⁰为热膨胀系数的参考温度；

上式中的第二项代表由于初始温度θ^I和参考温度θ⁰之间的温度差产生的热应变。这一项能够保证当参考温度不等于初始温度时初始热应变为0。

以同样的方式建立单个单元的有限元模型，并利用Predefined Field在单元上定义从0ºC到20ºC的温度变化，分别取参考温度θ⁰为0ºC和10ºC。在Matlab中利用上述公式编制热应变计算程序，程序如下所示。

%程序用于计算热应变Avg_alpha_data=[0 2e-5;10 3e-5;20 6e-5]; %平均热膨胀系数随温度变化Time_start=0;      %计算起始时间Time_end=1;        %计算终止时间Temp_start=0;      %起始温度Temp_end=20;       %终止温度Temp_ref=10;        %参考温度sample_num=100;    %增量步数量Time=linspace(Time_start,Time_end,sample_num)';Temp=linspace(Temp_start,Temp_end,sample_num)';thermal_strain=zeros(sample_num,1);for i=1:sample_num    %插值计算初始温度下的平均热膨胀系数    Avg_alpha_ini=interp1(Avg_alpha_data(:,1),...        Avg_alpha_data(:,2),Temp_start);    %插值计算当前温度下的平均热膨胀系数    Avg_alpha=interp1(Avg_alpha_data(:,1),...        Avg_alpha_data(:,2),Temp(i));    %计算热应变    thermal_strain(i)=Avg_alpha*(Temp(i)-Temp_ref)-...        Avg_alpha_ini*(Temp_start-Temp_ref);endplot(Time,thermal_strain)

将计算结果与ABAQUS进行对比，结果如图2.2所示。

图2.2 热应变随温度变化曲线

从图2.2中可以看出，ABAQUS与Matlab的计算结果完全一致，由于公式中的第二项总可以保证当参考温度θ⁰不等于初始温度θ^I时，因此单元的初始应变总是等于0。并且参考温度θ⁰的不同会对ABAQUS计算得到的应变值产生很大的影响，因此在定义随温度变化的平均热膨胀系数时，必须准确地指定与之对应的参考温度值。

3. 利用用户子程序UEXPAN计算热应变

ABAQUS用户子程序UEXPAN可以用于建立热应变与温度和/或预定义场变量和状态变量之间的复杂关系，通过指定热应变增量或热应变随温度的变化率实现对热应变的更新，该用户子程序将在单元的所有积分点上被调用，关于用户子程序UEXPAN的更多介绍参见如何在ABAQUS中定义热膨胀系数。

本文以表2.1中的平均线膨胀系数为例，通过UEXPAN来计算单元的热应变变化。计算思路为首先计算当前增量步i的起始时刻温度值θ^s和终止时刻温度值θ^f：

其中TEMP(1)为增量步终止时刻对应的温度值，TEMP(2)为当前增量步相对于前一增量步的温度增量。

随后根据当前增量步i的起始时刻温度值θ^s和终止时刻温度值θ^f计算前一增量步i-1和当前增量步i在终止时刻的热应变值：

其中θ⁰为参考温度，θ^I为初始温度，即材料热应变为0时的温度值。

因此从前一增量步i-1到当前增量步i产生的热应变增量为：

将热应变增量赋予给UEXPAN中的变量EXPAN即可实现热应变的更新。

利用上述算法编制得到的Fortran程序如下所示。

C     程序用于计算热应变      BLOCK DATA      IMPLICIT NONEC     该子程序用于存储全局变量数据C     定义随温度变化的热膨胀系数C     注意：由于FORTRAN中不能创建维度可变的数组,因此定义数组时只能使用常数值      REAL::MAT_TEMP(3)       !用于存放材料数据中的温度值      REAL::MAT_ALPHA(3)      !用于存在材料数据中与温度对应的平均热膨胀系数      INTEGER::DATA_NUM       !材料数据对个数            REAL::TEMP_INI          !分析开始时的初始温度(零应变时的温度)      REAL::TEMP_REF          !计算平热应变的参考温度C     通过关键词COMMON声明为全局变量      COMMON /THERMAL_STRAIN_DATA/ MAT_TEMP,MAT_ALPHA,DATA_NUM      COMMON /TEMP_PARA/ TEMP_INI,TEMP_REF            INTEGER::UEXPAN_ISREAD    !标识符,判断UEXPAN是否已经被调用      COMMON /SUBROUTINE_PARA/ UEXPAN_ISREAD      C     定义与热膨胀系数有关的数据      DATA MAT_TEMP /0,10,20/      DATA MAT_ALPHA /2E-5,3E-5,6E-5/      DATA DATA_NUM /3/      C     定义标识符为0,当UEXPAN首次被读取后,标识符将被置1      DATA UEXPAN_ISREAD /0/      C     定义初始温度和参考温度      DATA TEMP_INI /20/      DATA TEMP_REF /0/                  END BLOCK DATA                        SUBROUTINE UEXPAN(EXPAN,DEXPANDT,TEMP,TIME,DTIME,PREDEF,     1 DPRED,STATEV,CMNAME,NSTATV,NOEL)C      INCLUDE 'ABA_PARAM.INC'C      CHARACTER*80 CMNAMEC      DIMENSION EXPAN(*),DEXPANDT(*),TEMP(2),TIME(2),PREDEF(*),     1 DPRED(*),STATEV(NSTATV)C     TEMP_S为增量步起始时刻的温度值C     TEMP_F为增量步终止时刻的温度值C     TEMP_INI为分析初始温度C     TEMP_REF为计算热应变的参考温度C     MAT_TEMP和MAT_ALPHA用于定义随温度变化的热膨胀系数C     UEXPAN_ISREAD为标识符,用于输出定义的全局变量C     ALPHA_S为增量步起始时刻对应的平均热膨胀系数C     ALPHA_F为增量步终止时刻对应的平均热膨胀系数C     ALPHA_INI为分析初始温度对应的平均热膨胀系数C     THE_S为增量步起始时刻对应的热应变C     THE_F为增量步终止时刻对应的热应变C     变量声明以及读取全局变量C     注意使用COMMON读取全局变量时COMMON内部的变量定义顺序必须与BLOCK中的顺序一致      REAL::TEMP_S,TEMP_F,TEMP_INI,TEMP_REF      REAL::MAT_TEMP(3),MAT_ALPHA(3)      INTEGER::DATA_NUM      INTEGER::UEXPAN_ISREAD      REAL::ALPHA_S,ALPHA_F,ALPHA_INI      REAL::THE_S,THE_F            COMMON /TEMP_PARA/ TEMP_INI,TEMP_REF      COMMON /THERMAL_STRAIN_DATA/ MAT_TEMP,MAT_ALPHA,DATA_NUM      COMMON /SUBROUTINE_PARA/ UEXPAN_ISREAD      C*******************初始数据检查******************************C     检查全局变量是否已经成功读取到UEXPAN中C     将全局变量写入到ABAQUS的输出文件(.DAT)中C     在ABAQUS中通道号6为DAT文件      IF(UEXPAN_ISREAD.EQ.0)THEN          WRITE(6,*)'MATERIAL_DATA_TEMP',MAT_TEMP          WRITE(6,*)'MATERIAL_DATA_ALPHA',MAT_ALPHA          WRITE(6,*)'NUMBER OF DATA PAIR',DATA_NUM          WRITE(6,*)'INITIAL TEMPERATURE',TEMP_INI          WRITE(6,*)'REFERENCE TEMPERATURE',TEMP_REF          UEXPAN_ISREAD=1    !标识符置1,保证只会在第一次调用UEXPAN时输出数据      END IFC************************************************************C     **************计算热应变*********************      TEMP_F=TEMP(1)      TEMP_S=TEMP(1)-TEMP(2)      !TEMP(2)为相对于前一步的温度增量C     调用插值子程序计算平均热膨胀系数      CALL LINEAR_INTERPOLATION(TEMP_S,DATA_NUM,MAT_TEMP,     &MAT_ALPHA,ALPHA_S)   !增量步起始时刻的平均热膨胀系数            CALL LINEAR_INTERPOLATION(TEMP_F,DATA_NUM,MAT_TEMP,     &MAT_ALPHA,ALPHA_F)   !增量步终止时刻的平均热膨胀系数            CALL LINEAR_INTERPOLATION(TEMP_INI,DATA_NUM,MAT_TEMP,     &MAT_ALPHA,ALPHA_INI)  !分析初始温度对应的热膨胀系数      C     根据平均热膨胀系数计算热应变      THE_S=ALPHA_S*(TEMP_S-TEMP_REF)-     &ALPHA_INI*(TEMP_INI-TEMP_REF)            THE_F=ALPHA_F*(TEMP_F-TEMP_REF)-     &ALPHA_INI*(TEMP_INI-TEMP_REF)            EXPAN(1)=THE_F-THE_S     C     ********************************************      RETURN      END                        SUBROUTINE LINEAR_INTERPOLATION(INTER_X,NUM_DATA,XDATA,YDATA,     &INTER_Y)C     该子程序用于分段线性插值C     插值点必须位于数据点范围内,不能外插C     INTER_X(1)为需要进行插值的点(每次只能插值一个点)C     XDATA(N)和YDATA(N)为插值所用的数据点C     NUM_DATA(1)为所用数据点的个数NC     X1,X2和Y1,Y2为插值点所位于的曲线段C     注意:必须保证数据点个数与NUM_DATA一致C     INTER_Y(1)为插值点插值得到的结果          IMPLICIT NONEC     变量声明          REAL::INTER_X          INTEGER::NUM_DATA          REAL::XDATA(NUM_DATA),YDATA(NUM_DATA),INTER_Y          REAL::X1,X2,Y1,Y2C     索引指标声明          INTEGER::IC     判断插值点所位于的区间段      IF(INTER_X.EQ.XDATA(1))THEN    !插值点刚好为数据点的第一个起始点          INTER_Y=YDATA(1)          RETURN      END IF      DO I=1,NUM_DATA          IF(INTER_X.LE.XDATA(I))THEN              X2=XDATA(I)              Y2=YDATA(I)              X1=XDATA(I-1)              Y1=YDATA(I-1)              EXIT    !跳出循环并不再执行          END IF      END DO      C     根据曲线段进行插值      INTER_Y=(Y2-Y1)/(X2-X1)*(INTER_X-X1)+Y1                                   RETURN      END SUBROUTINE LINEAR_INTERPOLATION

上述程序中除了包含有用于计算热应变的UEXPAN之外，还定义有用于存储全局变量数据的子程序BLOCK DATA和用于进行线性插值计算的子程序LINEAR_INTERPOLATION。在通过子程序进行计算时，有时候可能需要在子程序中使用大量的常量，一般可以采用三种方式来读取常量：一种为将这些常量存放在文本文件中，然后在子程序计算时通过OPEN命令读取该文件以获取文件中的数据；另一种方式为将这些常量拟合为特殊的表达式，例如多项式，然后在子程序中定义多项式的系数，通过计算表达式来获取相关的数据；最简单的方式为直接在子程序中通过数组来存放这些数据，并在需要时直接调用数组。

第一种方式适用于数据量比较大的情况，但缺点是在子程序中大量地重复读写文件将会严重拖慢计算效率；第二种方式由于采用表达式进行计算，因此使用起来更加灵活，但缺点是并不是所有的数据都能拟合为特殊的表达式；第三种方式最为简单，但缺点是当数据量较大时编写程序较为繁琐。

本文由于使用的材料数据较少，因此采用第三种方式来存储相关的数据。并且，本文将这些材料数据定义为全局变量，并通过子程序BLOCK DATA进行存储，这样的优势有两点：首先，在整个分析过程中，这些材料数据只会在程序运行时被读取一次，而不会在程序运行时被重复读取，这样可以极大地提升计算效率；其次，当程序中存在多个子程序时，通过定义全局变量可以非常方便地在子程序之间进行数据传递，这与ABAQUS状态变量STATEV的功能非常类似。

为了便于程序调试，可以将需要调试的变量通过WRITE命令写入文件，在用户子程序中，可以写入文件的通道号（unitnumber）为15至18以及100以上的通道号，而其他通道号为ABAQUS保留的通道号，不能用于文件写入。对于ABAQUS/Standard，用户还可以使用通道号7将调试数据写入信息文件（.msg），或者使用通道号6将调试数据写入打印输出文件（.dat）。用户无需使用OPEN命令来打开该通道号，因为ABAQUS将自动打开通道号。