基于Abaqus用户子程序UEXTERNALDB调用外部文件数据
Abaqus用户子程序UEXTERNALDB常用于在Abaqus/Standard求解器运行时，实现Abaqus与外部程序之间的通信，如交换数据等。本文将通过一个简单案例，介绍如何基于Abaqus用户子程序UEXTERNALDB实现对外部文件数据的调用。
计算案例
考虑线弹性断裂力学中的一个典型案例：一个中心裂纹板的裂纹面上承受有非均布载荷作用，这个非均布载荷通过一个分段三次样条函数来描述。在第i个分段上，该非均布载荷可表示为：
式中：αi，βi，γi和δi为分段样条函数的系数，si-1和si为第i个分段的区间下限和上限。
为了求解中心裂纹板在该非均布载荷作用下的应力分布，我们可以通过Abaqus用户子程序DLOAD或解析场来定义该非均布载荷。因此，本文建立了如图1所示的1/4对称有限元模型。图中，紫色箭头代表了DLOAD子程序的载荷定义域，也就是非均布载荷的作用区域。
图1 中心裂纹板1/4对称有限元模型
事实上，当这个非均布载荷非常简单时，我们可以直接在DLOAD用户子程序中定义该非均布载荷的表达式。然而，这种方法的缺点是，每当我们需要更改该非均布载荷时，则需要重新修改并编译子程序。并且，当需要修改的系数较多时，直接修改子程序将非常繁琐。因此，一种更为简便的方式是将非均布载荷的系数数据存储在外部文件中，以供DLOAD子程序调用。这样，当我们需要修改非均布载荷的系数时，只需要更改外部文件的数据即可，而无需修改并重新编译子程序。
为此，我们构建了一个如下所示的子程序。该子程序由DLOAD子程序和UEXTERNALDB子程序构成。UEXTERNALDB子程序用于读取定义非均布载荷系数的外部文件，并将数据传递给出DLOAD子程序。而DLOAD子程序则根据传递得到的系数实现非均布载荷的施加。
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C 程序用于在裂纹面上施加非均布载荷
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C-------------------------------------------------------------------------
C 读取存放非均布载荷分段样条曲线系数的外部文本文件
C-------------------------------------------------------------------------
SUBROUTINE UEXTERNALDB(LOP,LRESTART,TIME,DTIME,KSTEP,KINC)
C
INCLUDE 'ABA_PARAM.INC'
C
DIMENSION TIME(2)
DIMENSION INTV(1),REALV(1)
CHARACTER*256 CHARV(1)
CHARACTER*256 OUTDIR !Abaqus结果文件输出目录
CHARACTER*14 EDS_FILE_NAME !非均布载荷系数外部文件名
INTEGER::I !循环计数指标
INTEGER::FILE_PT_INDEX=1 !外部文件读取索引
INTEGER,PARAMETER::EDS_COEFF_NUM=6 !非均布载荷系数数量
REAL::EDS_DATA_ALPHA(EDS_COEFF_NUM) !分段样条曲线系数Alpha
REAL::EDS_DATA_BETA(EDS_COEFF_NUM) !分段样条曲线系数Beta
REAL::EDS_DATA_GAMMA(EDS_COEFF_NUM) !分段样条曲线系数Gamma
REAL::EDS_DATA_DELTA(EDS_COEFF_NUM) !分段样条曲线系数Delta
REAL::EDS_DATA_LOWER_RANGE(EDS_COEFF_NUM) !分段样条曲线区间下限
REAL::EDS_DATA_UPPER_RANGE(EDS_COEFF_NUM) !分段样条曲线区间上限
COMMON /EDS_DATA/ EDS_DATA_ALPHA, !定义为全局变量
& EDS_DATA_BETA,EDS_DATA_GAMMA,
& EDS_DATA_DELTA,EDS_DATA_LOWER_RANGE,
& EDS_DATA_UPPER_RANGE
C 读取Abaqus结果文件输出目录(工作目录)
CALL GETOUTDIR(OUTDIR,LENOUTDIR)
EDS_FILE_NAME="EDS_Coeff.txt"
C 读取非均布载荷分段样条曲线系数
IF(LOP.EQ.0)THEN
OPEN(110,FILE=TRIM(OUTDIR)//"\"//EDS_FILE_NAME)
DO WHILE(.NOT.EOF(110))
READ(110,*) EDS_DATA_ALPHA(FILE_PT_INDEX),
& EDS_DATA_BETA(FILE_PT_INDEX),
& EDS_DATA_GAMMA(FILE_PT_INDEX),
& EDS_DATA_DELTA(FILE_PT_INDEX),
& EDS_DATA_LOWER_RANGE(FILE_PT_INDEX),
& EDS_DATA_UPPER_RANGE(FILE_PT_INDEX)
FILE_PT_INDEX=FILE_PT_INDEX+1
END DO
CLOSE(110)
END IF
C 输出非均布载荷分段样条曲线系数至MSG文件
CALL STDB_ABQERR(1,"---------------------------",INTV,REALV,CHARV)
CALL STDB_ABQERR(1,"EDS Coefficient Output ",INTV,REALV,CHARV)
C 输出非均布载荷分段样条曲线系数Alpha
CALL STDB_ABQERR(1,"---------------------------",INTV,REALV,CHARV)
CALL STDB_ABQERR(1,"Coefficient Alpha Table: ",INTV,REALV,CHARV)
DO I=1,EDS_COEFF_NUM
INTV(1)=I
REALV(1)=EDS_DATA_ALPHA(I)
CALL STDB_ABQERR(1,"Alpha(%I)=%R ",INTV,REALV,CHARV)
END DO
CALL STDB_ABQERR(1,"---------------------------",INTV,REALV,CHARV)
C 输出非均布载荷分段样条曲线系数Beta
CALL STDB_ABQERR(1,"Coefficient Beta Table: ",INTV,REALV,CHARV)
DO I=1,EDS_COEFF_NUM
INTV(1)=I
REALV(1)=EDS_DATA_BETA(I)
CALL STDB_ABQERR(1,"Beta(%I)=%R ",INTV,REALV,CHARV)
END DO
CALL STDB_ABQERR(1,"---------------------------",INTV,REALV,CHARV)
C 输出非均布载荷分段样条曲线系数Gamma
CALL STDB_ABQERR(1,"Coefficient Gamma Table: ",INTV,REALV,CHARV)
DO I=1,EDS_COEFF_NUM
INTV(1)=I
REALV(1)=EDS_DATA_GAMMA(I)
CALL STDB_ABQERR(1,"Gamma(%I)=%R ",INTV,REALV,CHARV)
END DO
CALL STDB_ABQERR(1,"---------------------------",INTV,REALV,CHARV)
C 输出非均布载荷分段样条曲线系数Delta
CALL STDB_ABQERR(1,"Coefficient Delta Table: ",INTV,REALV,CHARV)
DO I=1,EDS_COEFF_NUM
INTV(1)=I
REALV(1)=EDS_DATA_DELTA(I)
CALL STDB_ABQERR(1,"Delta(%I)=%R ",INTV,REALV,CHARV)
END DO
CALL STDB_ABQERR(1,"---------------------------",INTV,REALV,CHARV)
C 输出非均布载荷分段样条曲线区间下限
CALL STDB_ABQERR(1,"EDS Lower Limit Table: ",INTV,REALV,CHARV)
DO I=1,EDS_COEFF_NUM
INTV(1)=I
REALV(1)=EDS_DATA_LOWER_RANGE(I)
CALL STDB_ABQERR(1,"Lower Limit(%I)=%R ",INTV,REALV,CHARV)
END DO
CALL STDB_ABQERR(1,"---------------------------",INTV,REALV,CHARV)
C 输出非均布载荷分段样条曲线区间上限
CALL STDB_ABQERR(1,"EDS Upper Limit Table: ",INTV,REALV,CHARV)
DO I=1,EDS_COEFF_NUM
INTV(1)=I
REALV(1)=EDS_DATA_UPPER_RANGE(I)
CALL STDB_ABQERR(1,"Upper Limit(%I)=%R ",INTV,REALV,CHARV)
END DO
CALL STDB_ABQERR(1,"---------------------------",INTV,REALV,CHARV)
RETURN
END
C-------------------------------------------------------------------------
C 施加非均布载荷分布至裂纹面
C-------------------------------------------------------------------------
SUBROUTINE DLOAD(F,KSTEP,KINC,TIME,NOEL,NPT,LAYER,KSPT,
& COORDS,JLTYP,SNAME)
C
INCLUDE 'ABA_PARAM.INC'
C
DIMENSION TIME(2),COORDS(3)
CHARACTER*80 SNAME
CHARACTER*256 CHARV(1)
CHARACTER*256 OUTDIR
INTEGER::I !循环计数索引
INTEGER::INTERVAL_INDEX !分段样条区间号
C 读取全局变量数据
INTEGER,PARAMETER::EDS_COEFF_NUM=6 !非均布载荷系数数量
REAL::EDS_DATA_ALPHA(EDS_COEFF_NUM) !分段样条曲线系数Alpha
REAL::EDS_DATA_BETA(EDS_COEFF_NUM) !分段样条曲线系数Beta
REAL::EDS_DATA_GAMMA(EDS_COEFF_NUM) !分段样条曲线系数Gamma
REAL::EDS_DATA_DELTA(EDS_COEFF_NUM) !分段样条曲线系数Delta
REAL::EDS_DATA_LOWER_RANGE(EDS_COEFF_NUM) !分段样条曲线区间下限
REAL::EDS_DATA_UPPER_RANGE(EDS_COEFF_NUM) !分段样条曲线区间上限
COMMON /EDS_DATA/ EDS_DATA_ALPHA, !定义为全局变量
& EDS_DATA_BETA,EDS_DATA_GAMMA,
& EDS_DATA_DELTA,EDS_DATA_LOWER_RANGE,
& EDS_DATA_UPPER_RANGE
C 判断当前材料点是否位于裂纹面
IF ((COORDS(1).LT.EDS_DATA_LOWER_RANGE(1)).OR.
& (COORDS(1).GT.EDS_DATA_UPPER_RANGE(EDS_COEFF_NUM)))THEN
CALL STDB_ABQERR(-3,"Loading region is beyond crack surface",
& INTV,REALV,CHARV)
END IF
C 判断当前材料点所在分段样条区间号
DO I=1,EDS_COEFF_NUM
IF ((COORDS(1).GE.EDS_DATA_LOWER_RANGE(I)).AND.
& (COORDS(1).LE.EDS_DATA_UPPER_RANGE(I)))THEN
INTERVAL_INDEX=I
EXIT
END IF
END DO
C 计算当前材料点上的非均布载荷值
F=EDS_DATA_ALPHA(INTERVAL_INDEX)*COORDS(1)**3+
& EDS_DATA_BETA(INTERVAL_INDEX)*COORDS(1)**2+
& EDS_DATA_GAMMA(INTERVAL_INDEX)*COORDS(1)+
& EDS_DATA_DELTA(INTERVAL_INDEX)
RETURN
END
下面，本文将对该子程序的实现流程做简要介绍。在UEXTERNALDB子程序中，字符串变量EDS_FILE_NAME用于定义存放非均布载荷外部数据的文件名。EDS_COEFF_NUM变量为非均布载荷的分段数，为了保证数据能被正常读取，必须填写正确的分段数。EDS_DATA_ALPHA等变量则用于定义非均布载荷的系数和区间上、下限。这里，COMMON指令用于将这些变量定义为全局变量，使其能够被传递到DLOAD子程序中。然后，UEXTERNALDB子程序将根据指定的文件名找到对应的外部数据文件，并通过READ函数循环读取外部文件数据，并将其存放到相应的变量中。注意，这里我们使用了UEXTERNALDB子程序提供的LOP变量来判断读取数据文件的时机。这里，我们设置为当LOP变量取值为0时开始读取外部文件数据，这代表数据文件是在计算开始时被调用的。我们也可以通过判断LOP变量的取值来实现在计算过程的其他阶段来读取外部文件。事实上，这正是使用UEXTERNALDB子程序读取外部文件的优势。它允许我们可以控制在计算过程的特定阶段来读写外部文件数据。当然，我们也可以在DLOAD子程序中通过READ函数来读取外部文件。然而，由于DLOAD子程序将在每次进入积分点时均被调用一次，这意味在每个积分点上，都将执行一次文件读写，这将严重降低计算效率。虽然我们可以通过指定增量步或分析步号来判断外部文件的读写时机，但这并不是一种较好的解决方法。在DLOAD子程序中，我们将同样通过COMMON指令来接收全局变量的取值，并基于外部文件定义的数据来实现非均布载荷的施加，这里不再赘述。需要该案例完整计算文件的读者可在公众 号回复ARB_CFT获取。
UEXTERNALDB子程序简介
用户子程序UEXTERNALDB：
l 可以在每次分析的开始时刻、每个增量步的开始时刻、每个增量步的结束时刻以及每次分析的结束时刻被分别调用一次（此外，也会在重启动分析的开始时刻被调用一次）；
l 可以用于与其他软件或Abaqus/Standard中的其他子程序进行通信；
l 可以用于在分析的开始时刻打开用于其他用户子程序的外部文件，或者在分析的结束时刻关闭这些外部文件；
l 可以用于在每个增量步的开始时刻计算或读取历程信息。这些信息能够被写入用户自定义COMMON块或外部文件，以便在分析中被其他子程序调用；
l 还可以用于将用户计算的历程信息的当前值写入到外部文件。
子程序接口
需要定义的变量
用户子程序UEXTERNALDB在被调用后无需用户更改子程序中变量的取值。
提供参考信息的变量
当用户子程序UEXTERNALDB被调用后，以下变量的取值可以被用户读取，以便为需要执行的操作提供参考信息。
LOP
LOP=0表明子程序正在分析的开始时刻被调用。
LOP=1表明子程序正在当前分析增量步的开始时刻被调用。如果在当前分析增量步下计算无法收敛，需要更小的时间增量，则该子程序会在分析增量步的开始时刻被多次调用。
LOP=2表明子程序正在当前分析增量步的结束时刻被调用。此时，所有在重启动分析中用户需要用到信息必须被写入到外部文件。
LOP=3表明子程序正在分析的结束时刻被调用。
LOP=4表明子程序正在重启动分析的开始时刻被调用。此时，所有必须的外部文件都应该被打开并合理放置，所有重启动分析用到信息也应该从外部文件中被读入。
LOP=5表明子程序正在当前分析步的开始时刻被调用。变量KSTEP中包含了当前的分析步号。
LOP=6表明子程序正在当前分析步的结束时刻被调用。变量KSTEP中包含了当前的分析步号。
TIME(1)
当前分析步时间
TIME(2)
总的分析步时间
DTIME
时间增量
KSTEP
当前分析步号。如果LOP=4，则KSTEP为重启动分析步号。
KINC
当前增量步号。如果LOP=4，则KINC为重启动增量步号。
