黏弹性：广义Maxwell模型的UMAT实现
这篇文章是广义Maxwell黏弹性模型的UMAT实现。之前已经推导过基础的黏弹性公式：
以及广义Maxwell的黏弹性公式：
这篇就把前面推送里的公式落实在Abaqus UMAT自定义本构模型代码里。其他求解器的自定义本构模型公式写法也都相似，可以参考。
1. UMAT公式-勘误与补充
前面两篇文章果然出现了疏漏和笔误。本节做一点勘误和补充。
勘误1 - 第一篇UMAT代码中的错误
在第一篇黏弹性推送中，
那段UMAT代码里，从第34行开始，计算瞬时模量的公式，Enu1都应该在分母而不是分子：
C 原始（错误）代码：
IF(JSTEP(2).NE.21) Then
ElamInst = Elambda1*Enu1
EmuInst = Emu1*Enu1
EBULK3 = 3*ElamInst+2*EmuInst
EG2 = 2*EmuInst
EG = EmuInst
EG3 = 3*EmuInst
ELAM = ElamInst
C 修正的代码应该是：
IF(JSTEP(2).NE.21) Then
ElamInst = Elambda1/Enu1
EmuInst = Emu1/Enu1
EBULK3 = 3*ElamInst+2*EmuInst
EG2 = 2*EmuInst
EG = EmuInst
EG3 = 3*EmuInst
ELAM = ElamInst
对应的公式为：
注意到这里 对应的是长期模量。想要获得瞬时模量，需要：
勘误2 - 从一般形式到黏壶应变形式的推导过程
顺便把上一篇文章里很“显然”的那个推导来详细推一遍。
一阶广义Maxwell黏弹性模型，从一般形式推导到使用瞬时模量和黏壶应变的形式。
在第一篇黏弹性文章里，应力更新公式写为：(上一篇文章里的公式 写错了，等号右侧第2,3项漏乘了2倍。第一篇里的(7-2),(7-3)是对的，下面的形式也是对的)
考虑 ，式(1)可以改写为：
考虑到：
以及
式(2)可以再改写为：
从而得到：
式(5)就是使用瞬时弹性常数和黏壶应变表示的一阶广义Maxwell黏弹性模型的应力更新公式。
改用K,G和归一化模量系数的形式
在Abaqus中，黏弹性模型的剪切和体积响应是独立定义的，并且使用了归一化的无量纲系数来描述第i阶黏弹性组件的刚度与瞬时刚度之比值。对于一阶广义Maxwell黏弹性模型，需要输入的参数就包括 。它们和Lamé常数的关系为：
式(5)可以改写为：
这个形式瞅着别扭，是因为一般在弹性力学里面，使用K和G描述本构关系时，都会把应变张量做极分解，分解成应变偏张量和应变球张量。（在 黏弹性：从一维模型到三维UMAT本构实现 文章里面式(4-1)介绍过相关推导）
在代码实现里面，其实还是用Lamé常数来得更直观、更便捷。
2. Abaqus-UMAT子程序调试方法
24年1月，仿真秀的 九千CAE 老师在直播中介绍了使用Visual Studio对Abaqus子程序进行调试的方法。直播回放链接在这里。
九千CAE：ABAQUS UMAT与VUMAT材料子程序入门方法
具体方法流程是这样的：
（当然首先要正确安装Abaqus和Fortran，以及Visual Studio）
第一步，修改.env环境文件。
在目录：
X:\SIMULIA\EstProducts\20xx\win_b64\SMA\site\win86_64.env # （旧版本为abaqus_v6.env）
这个环境文件中，compile_fortran= 命令下，
'/Od', '/Ob0', # <-- Disable Optimization for Debugging
'/Zi', # <-- Debugging Information
去掉前面的注释符号 #（此处已经去掉）；
在 link_sl 和 link_exe命令下，
'/debug', # <-- Debugging
去掉/debug前面的注释符号。
第二步，在子程序中合适的位置插入暂停代码。
一般可以是输入代码：
write(*,*) "Please input an integer"
read(*,*) tempRead
第三步，用命令行提交作业。
在inp文件和.for文件对应的目录下，右键-在此处打开Powershell窗口。（这种方式比直接使用命令行cd到指定目录来得更快）。abaqus提交任务的命令格式为：
abaqus job=xxx user=xxx int cpus=x ask=no
job和user后面的文件名都不必加扩展名。cpus指示多核，对于简单的验证案例来说也可以不加。具体语法可以参考帮助文档。
由于前面在子程序里添加了输入的代码，这里如果编译正确，程序会停在第一次调用UMAT的时候，在Powershell窗口上要求用户输入信息以继续。
第四步，用Visual Studio附加到进程
启动Visual Studio，打开.inp文件所在的文件夹，点开这个项目调用编译的.for文件。在Visual Studio菜单栏点击【调试】-【附加到进程】，搜索standard.exe（如果是显式分析就找到explicit.exe），附加。
有需要的话，在代码后面的位置增加断点。然后回到命令行窗口，随意输入一个数值。这时你的UMAT运行状态就可以实时显示在Visual Studio窗口上了。包括各类变量的数值都可以在界面上找到。
（这部分我就不放截图了，感兴趣的朋友如果看文字觉得不太清晰，可以去找视频。上面的直播回放是免费的）
3. 与Abaqus内置黏弹性本构的对比验证
关于广义Maxwell模型的阶次：
二阶广义Maxwell模型，其实是可以直接向下兼容一阶模型的，只需要把二阶对应的两个弹性常数和弛豫时间都设置为零即可。所以这段代码就直接写成二阶的形式了。
关于额外的黏弹性应变：
对于高阶广义Maxwell模型，前面的推导中引入了额外的变量（黏壶应变 ）。它是一个应变张量，在Abaqus中以向量形式存储，有6个分量。二阶广义Maxwell模型就有两个黏壶应变，所以一共需要12个标量来存储相应的信息。这就需要在材料设置里，添加12个状态变量（STATEV）。前6个分量存储第一个黏壶应变，后6个分量存储第二个黏壶应变。
读者如果要修改为更高阶的Maxwell黏弹性模型，还需要再增加更多的状态变量以匹配其阶数。
UMAT验证案例：
这次也使用Abaqus内置的线弹性-黏弹性模型与UMAT计算结果做对比验证。下面的参数都是随便写的，不具有实际物理意义。瞬时弹性参数取为：
UMAT输入的参数按顺序依次为：
注：本表格中的参数保留小数点后2位，进行了四舍五入。后续计算结果中的误差是由于此处舍入导致，并非UMAT列式问题。
根据前面式 (7)，可以计算出对应的Abaqus内置本构的黏弹性参数输入为：
| 参数 | 数值 |
|---|---|
| g_1 | 0.25 |
| k_1 | 0.28 |
| tau_1 | 0.5 |
| g_2 | 0.2 |
| k_2 | 0.17 |
| tau_2 | 1.5 |
验证模型仅包含2个边长为1的单元。三面设置对称边界条件，Z+面上其他三个点施加Z+向强制拉伸位移，边缘顶点处施加X,Y,Z方向强制位移，以验证拉伸和剪切混合加载条件下的黏弹性响应。2个单元分别赋予Abaqus内置黏弹性本构和UMAT本构。
设置两个分析步，分别是静态结构分析步（研究瞬时弹性响应）和3s的粘性分析步（研究黏性响应）。提交计算后，结果为：
两条曲线数值在第三位有效数字后面有细微差异，主要来自前面参数输入时的舍入误差。结果表明当前UMAT能够准确模拟线弹性-二阶广义Maxwell黏弹性本构模型的响应。
4. inp输入文件与UMAT源代码
免责声明：作者编程水平一般，对Fortran 77的语法也不是特别熟悉。这段UMAT代码仅实现了基本功能，未做更多扩展。本节的代码只用于展示UMAT编写思路。如果读者将这段代码用于自己的项目，请仔细检查其功能。作者不对代码的计算结果和可能造成的其他影响负责。
这段代码里也没有包含黏弹性应变能的部分。如果读者有需要，也可以参考前面推送里的公式，自己推导对应二阶广义Maxwell黏弹性本构的应变能密度表达式。
C=====================================================C
C 公 众号 《CAE知识地图》 C
C 作者 毕小喵 C
C 线弹性-二阶广义Maxwell黏弹性 C
C=====================================================C
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
+ RPL,DDSDDT,DRPLDE,DRPLDT,
+ STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
+ NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
+ CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATV),
+ DDSDDE(NTENS,NTENS),
+ DDSDDT(NTENS),DRPLDE(NTENS),
+ STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
+ PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
+ JSTEP(4),DSDEV1(NTENS,NTENS),DSDEV2(NTENS,NTENS),
+ STREND(NTENS),Ealpha1(NTENS),Dalpha1(NTENS),
+ Ealpha2(NTENS),Dalpha2(NTENS)
C
IF(NPROPS.EQ.5) THEN
Elam0 = PROPS(1)
Emu0 = PROPS(2)
Elam1 = PROPS(3)
Emu1 = PROPS(4)
Etau1 = PROPS(5)
Elam2 = 0
Emu2 = 0
Etau2 = 0
ELSE IF(NPROPS.EQ.8) THEN
Elam0 = PROPS(1)
Emu0 = PROPS(2)
Elam1 = PROPS(3)
Emu1 = PROPS(4)
Etau1 = PROPS(5)
Elam2 = PROPS(6)
Emu2 = PROPS(7)
Etau2 = PROPS(8)
ELSE
WRITE(6,*) 'This program only fit 5 or 8 input variables.
1 please read the source code and obey the rule.'
ENDIF
C
IF(TIME(2).EQ.0.0) THEN
Ealpha1 = 0.0
Ealpha2 = 0.0
STATEV = 0.0
ENDIF
C
C determine step
C
IF(JSTEP(2).NE.21) Then
C
DDSDDE=0.0
Do K1=1,NDI
Do K2=1,NDI
DDSDDE(K1,K2)=Elam0
ENDDO
DDSDDE(K1,K1)=Elam0+2*Emu0
ENDDO
Do K1=NDI+1,NTENS
DDSDDE(K1,K1)=Emu0
ENDDO
STREND = STRAN + DSTRAN
STRESS=MATMUL(DDSDDE,STREND)
C
ELSE
C
C Begin viscoelastic
C read visco strain alpha
C
DO K1=1,NTENS
Ealpha1(K1) = STATEV(K1)
ENDDO
DO K1=NTENS+1,NTENS+NTENS
Ealpha2(K1-NTENS) = STATEV(K1)
ENDDO
C
TermDAE1 = DTIME/(DTIME+2*Etau1)
Dalpha1 = TermDAE1*(DSTRAN+2*STRAN-2*Ealpha1)
Ealpha1 = Ealpha1+Dalpha1
TermDAE2 = DTIME/(DTIME+2*Etau2)
Dalpha2 = TermDAE2*(DSTRAN+2*STRAN-2*Ealpha2)
Ealpha2 = Ealpha2+Dalpha2
C
C save visco strain alpha
C
DO K1=1,NTENS
STATEV(K1) = Ealpha1(K1)
ENDDO
DO K1=NTENS+1,NTENS+NTENS
STATEV(K1) = Ealpha2(K1-NTENS)
ENDDO
C
C evaluate STREND
C
DDSDDE=0.0
DO K1=1,NDI
DO K2=1,NDI
DDSDDE(K1,K2)=Elam0
DSDEV1(K1,K2)=Elam1
DSDEV2(K1,K2)=Elam2
ENDDO
DDSDDE(K1,K1)=Elam0+2*Emu0
DSDEV1(K1,K1)=Elam1+2*Emu1
DSDEV2(K1,K1)=Elam2+2*Emu2
ENDDO
DO K1=NDI+1,NTENS
DDSDDE(K1,K1)=Emu0
DSDEV1(K1,K1)=Emu1
DSDEV2(K1,K1)=Emu2
ENDDO
STREND = STRAN + DSTRAN
STRESS = MATMUL(DDSDDE,STREND)-MATMUL(DSDEV1,Ealpha1)
+ -MATMUL(DSDEV2,Ealpha2)
C
C evaluate DDSDDE
C
TermA = Elam0+2*Emu0-TermDAE1*(Elam1+2*Emu1)
+ -TermDAE2*(Elam2+2*Emu2)
TermB = Elam0-TermDAE1*Elam1-TermDAE2*Elam2
TermC = Emu0-TermDAE1*Emu1-TermDAE2*Emu2
DO K1=1,NDI
DO K2=1,NDI
DDSDDE(K1,K2)=TermB
ENDDO
DDSDDE(K1,K1)=TermA
ENDDO
DO K1=NDI+1,NTENS
DDSDDE(K1,K1)=TermC
ENDDO
C
ENDIF
RETURN
END
下面是验证使用的inp文件。
*Heading
** Job name: shear_verify Model name: visco_5_8
** Generated by: Abaqus/CAE 2022
*Preprint, echo=NO, model=NO, history=NO, contact=NO
**
** PARTS
**
*Part, name=PART-1-1
*Node
1, 1., 1., 1.
2, 1., 0., 1.
3, 1., 1., 0.
4, 1., 0., 0.
5, 0., 1., 1.
6, 0., 0., 1.
7, 0., 1., 0.
8, 0., 0., 0.
*Element, type=C3D8
1, 5, 6, 8, 7, 1, 2, 4, 3
*Nset, nset=SET-1, generate
1, 8, 1
*Elset, elset=SET-1
1,
** Section: Section-1-SET-1
*Solid Section, elset=SET-1, material=ABQ_VISCO
,
*End Part
**
*Part, name=PART-1-UMAT-1
*Node
1, 1., 2.5, 1.
2, 1., 1.5, 1.
3, 1., 2.5, 0.
4, 1., 1.5, 0.
5, 0., 2.5, 1.
6, 0., 1.5, 1.
7, 0., 2.5, 0.
8, 0., 1.5, 0.
*Element, type=C3D8
1, 5, 6, 8, 7, 1, 2, 4, 3
*Nset, nset=SET-1, generate
1, 8, 1
*Elset, elset=SET-1
1,
** Section: Section-2-SET-1
*Solid Section, elset=SET-1, material=UMAT-VISCO8
,
*End Part
**
**
** ASSEMBLY
**
*Assembly, name=Assembly
**
*Instance, name=PART-1-1, part=PART-1-1
*End Instance
**
*Instance, name=PART-1-UMAT-1, part=PART-1-UMAT-1
*End Instance
**
*Nset, nset=FIX, instance=PART-1-1
8,
*Nset, nset=FIX, instance=PART-1-UMAT-1
8,
*Nset, nset=SET-6, instance=PART-1-1
1, 2, 5, 6
*Elset, elset=SET-6, instance=PART-1-1
1,
*Nset, nset=SET-7, instance=PART-1-UMAT-1
1, 2, 5, 6
*Elset, elset=SET-7, instance=PART-1-UMAT-1
1,
*Nset, nset=Set-9, instance=PART-1-UMAT-1
1,
*Nset, nset=Set-9, instance=PART-1-1
1,
*Nset, nset=XSYMM, instance=PART-1-1, generate
5, 8, 1
*Nset, nset=XSYMM, instance=PART-1-UMAT-1, generate
5, 8, 1
*Elset, elset=XSYMM, instance=PART-1-1
1,
*Elset, elset=XSYMM, instance=PART-1-UMAT-1
1,
*Nset, nset=YSYMM, instance=PART-1-1, generate
2, 8, 2
*Nset, nset=YSYMM, instance=PART-1-UMAT-1, generate
2, 8, 2
*Elset, elset=YSYMM, instance=PART-1-1
1,
*Elset, elset=YSYMM, instance=PART-1-UMAT-1
1,
*Nset, nset=ZSYMM, instance=PART-1-1
3, 4, 7, 8
*Nset, nset=ZSYMM, instance=PART-1-UMAT-1
3, 4, 7, 8
*Elset, elset=ZSYMM, instance=PART-1-1
1,
*Elset, elset=ZSYMM, instance=PART-1-UMAT-1
1,
*Elset, elset=_SURF-1_S3, internal, instance=PART-1-1
1,
*Elset, elset=_SURF-2_S3, internal, instance=PART-1-UMAT-1
1,
*Surface, type=ELEMENT, name=SURF-1
_SURF-1_S3, S3
*Surface, type=ELEMENT, name=SURF-2
_SURF-2_S3, S3
*End Assembly
**
** MATERIALS
**
*Material, name=ABQ_HYPER
*Hyperelastic, ogden
0.4, 6.4, 0.1
*Material, name=ABQ_VISCO
*Elastic, moduli=INSTANTANEOUS
100., 0.3
*Viscoelastic, time=PRONY
0.25, 0.28, 0.5
0.2, 0.17, 1.5
*Material, name=UMAT-VISCO5
*Depvar
2,
*User Material, constants=5
5., 3., 2.5, 1.5, 0.5
*Material, name=UMAT-VISCO8
*Depvar
12,
*User Material, constants=8
57.69, 38.46, 17.31, 9.62, 0.5, 8.65, 7.69, 1.5
** ----------------------------------------------------------------
**
** STEP: Step-static
**
*Step, name=Step-static, nlgeom=YES
*Static
0.1, 1., 1e-05, 0.1
**
** BOUNDARY CONDITIONS
**
** Name: Disp-BC-1 Type: 位移/转角
*Boundary
SET-6, 3, 3, 0.5
** Name: Disp-BC-2 Type: 位移/转角
*Boundary
SET-7, 3, 3, 0.5
** Name: Disp-BC-3 Type: 对称/反对称/完全固定
*Boundary
FIX, PINNED
** Name: Disp-BC-4 Type: 对称/反对称/完全固定
*Boundary
XSYMM, XSYMM
** Name: Disp-BC-5 Type: 对称/反对称/完全固定
*Boundary
YSYMM, YSYMM
** Name: Disp-BC-6 Type: 对称/反对称/完全固定
*Boundary
ZSYMM, ZSYMM
** Name: shear Type: 位移/转角
*Boundary
Set-9, 1, 1, 0.3
Set-9, 2, 2, 0.2
**
** OUTPUT REQUESTS
**
*Restart, write, frequency=0
**
** FIELD OUTPUT: F-Output-1
**
*Output, field
*Contact Output
CDISP, CSTRESS
**
** FIELD OUTPUT: F-Output-2
**
*Node Output
CF, RF, U
**
** FIELD OUTPUT: F-Output-3
**
*Element Output, directions=YES
LE, PE, PEEQ, PEMAG, S, SDV
**
** HISTORY OUTPUT: H-Output-1
**
*Output, history, variable=PRESELECT
*End Step
** ----------------------------------------------------------------
**
** STEP: Step-visco
**
*Step, name=Step-visco, nlgeom=YES, inc=1000
*Visco, cetol=0.05
0.2, 3., 3e-05, 0.2
**
** OUTPUT REQUESTS
**
*Restart, write, frequency=0
**
** FIELD OUTPUT: F-Output-4
**
*Output, field
*Contact Output
CDISP, CSTRESS
**
** FIELD OUTPUT: F-Output-5
**
*Node Output
CF, RF, U
**
** FIELD OUTPUT: F-Output-6
**
*Element Output, directions=YES
LE, PE, PEEQ, PEMAG, S, SDV
**
** HISTORY OUTPUT: H-Output-2
**
*Output, history, variable=PRESELECT
*End Step
如果这篇文章对你有启发或对你的研究工作有一点帮助，是我的荣幸。
