文章详细介绍了如何在Abaqus中使用Fortran开发子程序，包括DLOAD、DFLUX、DISP等子程序的开发，以及位移边界条件的设定，同时提到了VGETVRM、USDFLD和VUSDFLD在不同分析类型中的应用，并简述了弹塑性UMAT的理论基础和开发。
子程序开发
vs2012、fortran2013、abaqus2016的安装与关联教程：BV1eZ4y1M76M
abaqus卸载教程：BV1co4y1f7oU
vs使用
新建项目
项目结构
默认生成的是fortran90格式的文件，即后缀.f90，但abaqus CAE中支持的是老的格式fortran77，即.for后缀
abaqus支持自由格式、固定格式，默认是识别固定格式，帮助文档中给出的例子也都是固定格式
但可在win86_64.env文件(新版本abaqus)中修改，添加’/free’，使变成支持自由格式
若在编写代码时更喜欢使用自由格式，可在编写完后通过转换工具Forquil(http://quill.fcode.cn/)进行在线转换成固定格式
运行与调试
创建项目时选择将解决方案与项目放在统一目录的效果，与勾选没啥区别
创建新的源程序文件
fortran语法
固定格式使用C或*注释，自由格式使用!注释
固定格式前六列不能使用，但在换行时第六列可以用数字(其他字符也行)表示换行，固定格式只支持每行70列
program Console1
implicit none
! Variables
! Body of Console1
print *, 'Hello World'
end program Console1
implicit none表示不使用fortran默认的变量精度设置，所有变量的类型自己声明
fortran90默认将所有变量定义为双精度，i、j、k、l、m、n默认是整型，
! 基本语法练习
program main
!进入主程序
implicit none! 用户自行控制变量的数据类型，不使用fortran默认的设置
! 整型数据
!integer :: x, y, add
!x = 1
!y = 2
!add = x + y
!write (*, '(i10)') add ! 使用write打印变量结果 *指代后面的变量add 括号中的i代表输出整型，10代表输出10位有效数字 若不够10位 前面用空格代替占位
! 实型 相当于float
real :: x
x = 1.628
print *, x
! 默认字节数 huge()返回该类型可取到的最大值
integer :: a
real :: b
real(16) :: c
print *, huge(a), huge(b), huge(c)
! 逻辑型
logical :: flag1, flag2
flag1 = .true.
flag2 = .false.
print *, flag1, flag2
! 字符型
character (5) :: ch1
ch1 = "abcde"
write(*, *) ch1(1:3)
! 数组
real, dimension(5) :: arr1
data arr1 /1.0, 2.0, 3.0, 4.0, 5.0/
print *, arr1
! 点乘
integer, dimension(3) :: arr2
data arr2 /1, 2, 3/
integer, dimension(3) :: arr3
data arr3 /3,2,1/
integer :: dotpro
dotpro = dot_product(arr2, arr3)
print *, dotpro
! 矩阵
integer,dimension(2,2) :: matr1
data matr1 /1,2,3,4/
print *, matr1(2,2)!putput is 4
! 矩阵乘法
integer, dimension(2,2) :: matr2, matr3, matr4
data matr2 /1,2,3,4/
data matr3 /1,2,3,4/
matr4 = matmul(matr2, matr3)
print *, matr4
! 循环语句：if-else if-else
integer :: x = 10
if(x .EQ. 10) then
print *, 'x=10'
else if(x .GT. 10) then
print *, 'x>10'
else
print *, 'x<10'
end if
! 循环语句：do循环
integer :: n
do n = 1, 6
print *, n
end do
! 循环语句：do while循环
integer :: n =1
do while (n <= 6)
print *, n
n = n + 1
end do
! 跳出控制exit---直接退出循环 相当于break
integer :: n
do n = 1, 100
print *,n
if (n=10) then
exit
end if
end do
! 跳出控制cycle---跳过本次循环后面的语句，直接进入下一次循环 相当于continue
integer :: n
do n = 1, 20
if(n==10) then
cycle
end if
print *, n
end do
! 文件操作--新建并写入数据
integer :: i
open(unit=2021, file=myfile.txt', status='new')
do i = 1, 5
write(2021,*) i
end do
close(2021)
! 文件操作---读出数据
integer :: i
integer, dimension(5) :: p
open(unit-2021,file='myfile.txt', status='old')
do i = 1,5
read(2021,*) p(i)
end do
do i = 1,5
print *, p(i)
end do
close(2021)
pause ! 控制运行结果窗口不自动退出
end program main ! 结束主程序
DLOAD
SUBROUTINE DLOAD(F,KSTEP,KINC,TIME,NOEL,NPT,LAYER,KSPT,
1 COORDS,JLTYP,SNAME)
C
INCLUDE 'ABA_PARAM.INC'
C
DIMENSION TIME(2), COORDS (3)
CHARACTER*80 SNAME
C user coding to define F
x = COORDS(1)
y = COORDS(2)
z = COORDS(3)
vmove = 31.4 * TIME(1)
xori = (1+2.5*vmove)*cos(vmove)
yori = (1+2.5*vmove)*sin(vmove)
IF ((abs(x-xori) .le. 1) .and. (abs(y-yori) .le. 1)) THEN
F = 400.0
ELSE
F = 0.0
ENDIF
RETURN
END
DFLUX
本节案例：将两块板焊接在一起
DFLUX介绍：
- 在传热或质量扩散分析中，定义非均匀分布通量作为位置、时间、温度、单元号、积分点号等的函数
- 分析中的每个基于单元或基于表面的（仅传热）的非均匀分布通量在每个通量积分点被调用
- 忽略与关联的非均匀分布通量定义可能出现的任何振幅引用
- 将节点用作通量积分点，用于一阶传热、一阶耦合温度-位移、一阶耦合热电结构和质量扩散单元
SUBROUTINE DFLUX(FLUX,SOL,KSTEP,KINC,TIME,NOEL,NPT,COORDS,
1 JLTYP,TEMP,PRESS,SNAME)
C
INCLUDE 'ABA_PARAM.INC'
C
DIMENSION FLUX(2), TIME(2), COORDS(3)
CHARACTER*80 SNAME
user coding to define FLUX(1) and FLUX(2)
RETURN
END
材料参数：传导、密度、比热容，如果是温度位移耦合还需定义膨胀、潜热
单元类型：普通热传导分析使用DC3D8、温度位移耦合使用C3D8T
划分网格时使用偏置加密网格
SUBROUTINE DFLUX(FLUX,SOL,KSTEP,KINC,TIME,NOEL,NPT,COORDS,
1 JLTYP,TEMP,PRESS,SNAME)
C
INCLUDE 'ABA_PARAM.INC'
DIMENSION FLUX(2), TIME(2), COORDS(3)
CHARACTER*80 SNAME
C user coding to define FLUX(1) and FLUX(2)
Q = 1000
v = 5
0 Rh = 0.005
d = v*TIME(2)
x = COORDS(1)
y = COORDS(2)
x0 = 0
y0 = 0
z0 = 0
PI = 3.1415
R=sqrt((x-x0-d)**2+(y-y0)**2)
FLUX(1)=3*Q/(PI*Rh**2)*exp(-3*R**2/Rh**2)
RETURN
END
位移边界条件子程序DISP
DISP：主要用于定义位移（有时也用于定义温度）边界条件或connector单元运动，会覆盖CAE界面上设置的相应边界条件或connector设置类型的具体数值
子程序接口
SUBROUTINE DISP(U,KSTEP,KINC,TIME,NODE,NOEL,JDOF,COORDS)
C
INCLUDE 'ABA_PARAM.INC'
C
DIMENSION U(3),TIME(3),COORDS(3)
C
user coding to define U
RETURN
END
接口变量说明：
- U(1)：位移、孔隙压强、温度等通过幅值传递给程序，一般用作位移
- U(2)：速度，即位移的一阶导
- U(3)：加速度，即位移的二阶导
- KSTEP：分析步编号
- KINC：增量步编号
- TIME(1)：分析步时间
- TIME(2)：总时间
- TIME(3)：增量步时间
- NODE：节点编号
- NOEL：单元编号
- JDOF：自由度
- COORS：位置坐标，一般为原始坐标(即变形前)，几何非线性计算时定义connector单元时为前一增量步的坐标，不可用
一个例子
因为子程序DISP中U(1)是位移表达式，而位移的位置和方向需要在CAE界面操作设定
位移大小虽然也能通过CAE操作给定，但当位移的大小与时间、空间等其他变量相关的一个函数时需要通过子程序来写入位移大小表达式(即把U(1)写成参数的函数）
定义位移的位置和方向，选择使用子程序定义位移，暂定值等于1，之后会被子程序覆盖掉
tips：使用子程序时，先划分网格稀疏一点，避免子程序编写错误造成计算时间太长等错误
在提交作业时选择子程序文件(注意导入前要在vs保存源文件)：
除了在CAE提交任务时打开子程序文件，还可以通过Abaqus Command命令行输入abaqus job=Jobname user=filename
但这种使用方式要求事先准备好inp文件和子程序文件(后缀为.for或.obj），并在命令行使用cd命令将工作目录调整至inp文件和子程序文件所在的文件夹
SUBROUTINE DISP(U,KSTEP,KINC,TIME,NODE,NOEL,JDOF,COORDS)
INCLUDE 'ABA_PARAM.INC'
DIMENSION U(3),TIME(3),COORDS(3)
C user coding to define U
C U(1) = 1
U(1) = 0.5 * sin(TIME(1)*6.28)
RETURN
END
应用实例
与时间相关的位移表达式
魔性的摇摇椅：三个实体方块5 * 5 * 5，一个梁，两个方块和一个板设置tie绑定
basic_DISP2
SUBROUTINE DISP(U,KSTEP,KINC,TIME,NODE,NOEL,JDOF,COORDS)
INCLUDE 'ABA_PARAM.INC'
DIMENSION U(3),TIME(3),COORDS(3)
C user coding to define U
IF (COORDS(1) <= 5) THEN
U(1) = 0.5 * sin(TIME(1) * 12.56)
ELSE
U(1) = -0.5 * sin(TIME(1) * 12.56)
END IF
RETURN
END
955
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
