ABAQUS中的应力应变描述

更多阅读：sppy.site

应力应变描述

文档《Abaqus Analysis User’s Guide》->《Introduction, Spatial Modeling, and Execution》->《Introduction》->《Abaqus syntax and conventions》->《Conventions》->《Stress and strain measures》

应力描述

在ABAQUS中采用Cauchy应力作为应力描述方式，Cauchy应力也被称为真实应力，其表示当前构型中单位面积上的力。

应变描述

对于几何非线性（有限变形）分析，存在多种的应变描述方式。不同于真实应力，并没有哪种应变描述可以被明确地称为真实应变。

对于小变形，只有一种应变描述方式

在大应变分析中，对于相同的物理变形过程，选择不同的应变描述方式会得到不同的应变数值。实际上，应变描述方式的选择取决于分析类型、材料行为以及（某种程度上的）用户偏好。

Abaqus/Standard和Abaqus/Explicit的应变描述存在差异：

Abaqus/Standard：默认情况（小变形）下，应变输出是总应变（变量符号E）；对于大应变（有限变形）下的壳单元、薄膜单元以及实体单元，可以选择两种应变描述方式：对数应变（变量符号LE）、名义应变（变量符号NE）。
Abaqus/Explicit：对数应变（变量符号LE）是默认应变输出量，也可以选择输出名义应变（变量符号NE），总应变（变量符号E）是不可用的。

总应变

总应变，即Total strain，变量符号E

在参考构型中，对应变率进行数值积分得到的总应变
εn+1=ΔR⋅εn⋅ΔRT+Δε\boldsymbol{\varepsilon}^{n+1}=\Delta\boldsymbol{R}\cdot\boldsymbol{\varepsilon}^n\cdot\Delta\boldsymbol{R}^\mathrm{T}+\Delta\boldsymbol{\varepsilon} εn+1=ΔR⋅εn⋅ΔRT+Δε
如果单元采用了共旋坐标系，则上式可简化为
εn+1=εn+Δε\boldsymbol{\varepsilon}^{n+1}=\boldsymbol{\varepsilon}^n+\Delta\boldsymbol{\varepsilon} εn+1=εn+Δε
应变增量 Δε\Delta\boldsymbol{\varepsilon}Δε 可通过对变形率 D\boldsymbol{D}D 在时间增量 Δt\Delta tΔt 上积分得到
Δε=∫tntn+1Ddt\Delta\boldsymbol{\varepsilon}=\int^{t^{n+1}}_{t^n} \boldsymbol{D}~\mathrm{d}t Δε=∫tntn+1D dt
这种应变测量适用于弹性-(粘)塑性或弹性蠕变材料，因为塑性应变和蠕变应变是通过相同的积分过程获得的。在这类材料中，弹性应变很小，总应变近似等于塑性应变或蠕变应变。

格林应变

格林应变，即Green’s strain，变量符号E

在Abaqus/Standard中，对于小应变的壳单元和梁单元，默认应变描述是格林应变（变量符号仍为E）
εG=12(FT⋅F−I)\boldsymbol{\varepsilon}^G=\frac{1}{2}(\boldsymbol{F}^\mathrm{T}\cdot\boldsymbol{F}-\boldsymbol{I}) εG=21(FT⋅F−I)
式中，F\boldsymbol{F}F 为变形梯度，I\boldsymbol{I}I 为单位张量。格林应变适用于小应变大旋转的情况。在有限变形分析中，小应变壳和梁不应使用弹塑性或超弹性本构关系，因为可能会得到不正确的分析结果或发生程序故障。

名义应变

名义应变，即Nominal strain，变量符号NE

εN=V−I=∑i=13(λi−1)niniT\boldsymbol{\varepsilon}^N=\boldsymbol{V}-\boldsymbol{I}=\sum^3_{i=1}(\lambda_i-1)\boldsymbol{n}_i\boldsymbol{n}_i^\mathrm{T} εN=V−I=i=1∑3(λi−1)niniT

式中，V=F⋅FT\boldsymbol{V}=\sqrt{\boldsymbol{F}\cdot\boldsymbol{F}^\mathrm{T}}V=F⋅FT 为左拉伸张量，λi\lambda_iλi 与 ni\boldsymbol{n}_ini 分别为其特征值与特征(列)向量。

对数应变

对数应变，即Logarithmic strain，变量符号LE

εL=ln⁡V=∑i=13ln⁡λininiT\boldsymbol{\varepsilon}^L=\ln\boldsymbol{V}=\sum^3_{i=1}\ln\lambda_i\boldsymbol{n}_i\boldsymbol{n}_i^\mathrm{T} εL=lnV=i=1∑3lnλininiT

超弹性材料的应变输出为对数应变。对于超粘弹性材料，对数弹性应变EE是由当前松弛应力状态计算得到，粘弹性应变CE则等于LE-EE。

UMAT中的输入量

应变量

在有限变形问题中，应变分量在UMAT调用之前被旋转，并且近似于对数应变，即 εL\boldsymbol{\varepsilon}^LεL 。

变形梯度

传递到UMAT中的并非变形梯度 F\boldsymbol{F}F ，而是修正变形梯度 F‾\overline{\boldsymbol{F}}F
F‾=F(J‾J)1n\overline{\boldsymbol{F}}=\boldsymbol{F}\bigg(\frac{\overline{J}}{J}\bigg)^{\frac{1}{n}} F=F(JJ)n1
式中，J=det⁡(F)J=\det(\boldsymbol{F})J=det(F) 是高斯点上的雅克比行列式，对于三维单元，n=3n=3n=3 ；对于二维单元，n=2n=2n=2 。

J‾\overline{J}J 则是整个单元上的平均雅克比行列式
J‾=1Vel∫JdVel\overline{J}=\frac{1}{V_{el}}\int J~\mathrm{d}V_{el} J=Vel1∫J dVel

根据单元积分点位置，计算离散加权平均得到 J‾\overline{J}J

旋转增量矩阵

对于可逆的变形梯度 F\boldsymbol{F}F ，由矩阵的极分解定理可得
F=R⋅U=V⋅R\boldsymbol{F}=\boldsymbol{R}\cdot\boldsymbol{U}=\boldsymbol{V}\cdot\boldsymbol{R} F=R⋅U=V⋅R

式中，U=FT⋅F\boldsymbol{U}=\sqrt{\boldsymbol{F}^\mathrm{T}\cdot\boldsymbol{F}}U=FT⋅F 、V=F⋅FT\boldsymbol{V}=\sqrt{\boldsymbol{F}\cdot\boldsymbol{F}^\mathrm{T}}V=F⋅FT 。于是可得对应的正交旋转矩阵 R\boldsymbol{R}R
R=F⋅(FT⋅F)−1=(F⋅FT)−1⋅F\boldsymbol{R}=\boldsymbol{F}\cdot\big(\sqrt{\boldsymbol{F}^\mathrm{T}\cdot\boldsymbol{F}}\big)^{-1}=\big(\sqrt{\boldsymbol{F}\cdot\boldsymbol{F}^\mathrm{T}}\big)^{-1}\cdot\boldsymbol{F} R=F⋅(FT⋅F)−1=(F⋅FT)−1⋅F
UMAT中的输入量中有两个变形梯度，即 F0\boldsymbol{F}_0F0 与 F1\boldsymbol{F}_1F1 ，分别对应该增量步变形前后的变形梯度。于是可得两个旋转矩阵，即 R0\boldsymbol{R}_0R0 与 R1\boldsymbol{R}_1R1 ，进一步可定义旋转增量矩阵 ΔR\Delta\boldsymbol{R}ΔR
R0⋅ΔR=R1⟹ΔR=R0−1⋅R1\boldsymbol{R}_0\cdot\Delta\boldsymbol{R}=\boldsymbol{R}_1\quad\Longrightarrow\quad\Delta\boldsymbol{R}=\boldsymbol{R}_0^{-1}\cdot\boldsymbol{R}_1 R0⋅ΔR=R1⟹ΔR=R0−1⋅R1

旋转增量矩阵 ΔR\Delta\boldsymbol{R}ΔR 作为UMAT输入量，是直接给出的

验证

在Abaqus/Standard中，采用几何非线性（有限变形）分析，选择输出对数应变。

在UMAT中输出：旋转增量矩阵DROT、变形梯度DFGRD0和DFGRD1、应变STRAN、应变增量DSTRAN、时间增量DTIME。Fortran程序：

在ABAQUS中，应变及应变增量为工程应变

      PRINT*,'DROT=...'WRITE(*,12),DROT(1,:)WRITE(*,12),DROT(2,:)WRITE(*,12),DROT(3,:)WRITE(*,*)PRINT*,'F0=...'WRITE(*,12),DFGRD0(1,:)WRITE(*,12),DFGRD0(2,:)WRITE(*,12),DFGRD0(3,:)WRITE(*,*)PRINT*,'F1=...'WRITE(*,12),DFGRD1(1,:)WRITE(*,12),DFGRD1(2,:)WRITE(*,12),DFGRD1(3,:)WRITE(*,*)PRINT*,'STRAIN=...'WRITE(*,12),(/STRAN(1),STRAN(4)/2.0D0,STRAN(5)/2.0D0/)WRITE(*,12),(/STRAN(4)/2.0D0,STRAN(2),STRAN(6)/2.0D0/)WRITE(*,12),(/STRAN(5)/2.0D0,STRAN(6)/2.0D0,STRAN(3)/)WRITE(*,*)PRINT*,'DSTRAIN=...'WRITE(*,12),(/DSTRAN(1),DSTRAN(4)/2.0D0,DSTRAN(5)/2.0D0/)WRITE(*,12),(/DSTRAN(4)/2.0D0,DSTRAN(2),DSTRAN(6)/2.0D0/)WRITE(*,12),(/DSTRAN(5)/2.0D0,DSTRAN(6)/2.0D0,DSTRAN(3)/)WRITE(*,*)PRINT*,'DTIME=',DTIMEWRITE(*,*)PRINT*,'######################################'WRITE(*,*)

输出数据保存在.log文件中，并从中提出一组数据进行验证。MATLAB程序：

% 从'.log'文件中提取数据
clc;clear;
format long;J_average=1;% 假设整个单元上的平均雅克比行列式为1
n=3;% 三维单元funEG=@(F) (F'*F-eye(3))/2;% 格林应变
funF=@(F) (J_average/det(F))^(-1/n)*F;% 由修正变形梯度得到变形梯度
funEL=@(F) funm(sqrtm(F*F'), @log);% 对数应变
funR=@(F) F/sqrtm(F'*F);% 旋转张量DROT=...[0.9999883445678  -0.0045470816350   0.00162319964590.0045463626195   0.9999895656127   0.0004463772085-0.0016252124224  -0.0004389923516   0.9999985829841];F0=...[0.9162657696006  -0.5340116556566   0.01929695098570.0000000000000   1.8334564757673   0.0000000000000-0.2250158864930  -0.1009421521995   0.5917174211528];F1=...[0.9128914668255  -0.5485035655311   0.01841876334290.0000000000000   1.8631389668749   0.0000000000000-0.2265230025927  -0.0988821595961   0.5845670848930];STRAIN0=...[-0.0579149200638  -0.2779025519040  -0.1523878429897-0.2779025519040   0.5481697880120  -0.1121959991532-0.1523878429897  -0.1121959991532  -0.4882665215795];DSTRAIN=...[-0.0040262445184  -0.0045467488471  -0.0029882636856-0.0045467488471   0.0160594161089  -0.0004426873816-0.0029882636856  -0.0004426873816  -0.0120095770493];STRAIN1=STRAIN0+DSTRAIN;DTIME=3.098242187499978E-002;

应变

由 F0\boldsymbol{F}_0F0 计算得到 ε0L\boldsymbol{\varepsilon}^L_0ε0L，并与STRAIN0对比

funEL(F0)=-0.060803391544017  -0.274684825559404  -0.155797481282906-0.274684825559404   0.552143951989882  -0.103910410215967-0.155797481282906  -0.103910410215967  -0.489335618314627STRAIN0=-0.057914920063800  -0.277902551904000  -0.152387842989700-0.277902551904000   0.548169788012000  -0.112195999153200-0.152387842989700  -0.112195999153200  -0.488266521579500

也可以考虑变形梯度的修正关系

funEL(funF(F0))=-0.060135077500271  -0.274684825559404  -0.155797481282906-0.274684825559404   0.552812266033628  -0.103910410215966-0.155797481282906  -0.103910410215966  -0.488667304270882

同样地，可以对比ε1L\boldsymbol{\varepsilon}^L_1ε1L 与STRAIN1。

旋转增量矩阵

由 F0\boldsymbol{F}_0F0 、F1\boldsymbol{F}_1F1 分别得到 R0\boldsymbol{R}_0R0 、R1\boldsymbol{R}_1R1 ，进一步计算出 ΔR\Delta\boldsymbol{R}ΔR ，并与DROT对比

funR(F0)\funR(F1)=0.999993311642542  -0.003361293648173   0.0014416570993980.003361951026421   0.999994245655608  -0.000453807191346-0.001440123424377   0.000458650936687   0.999998857841268DROT=0.999988344567800  -0.004547081635000   0.0016231996459000.004546362619500   0.999989565612700   0.000446377208500-0.001625212422400  -0.000438992351600   0.999998582984100

速度梯度

变形率 D\boldsymbol{D}D
D:=ΔεΔt\boldsymbol{D}:=\frac{\Delta\boldsymbol{\varepsilon}}{\Delta t} D:=ΔtΔε
自旋率张量 W\boldsymbol{W}W
W:=2Δt(ΔR−I)⋅(ΔR+I)−1\boldsymbol{W}:=\frac{2}{\Delta t}(\Delta\boldsymbol{R}-\boldsymbol{I})\cdot(\Delta\boldsymbol{R}+\boldsymbol{I})^{-1} W:=Δt2(ΔR−I)⋅(ΔR+I)−1
速度梯度 L\boldsymbol{L}L
L:=2Δt(F1−F0)⋅(F1+F0)−1\boldsymbol{L}:=\frac{2}{\Delta t}(\boldsymbol{F}_1-\boldsymbol{F}_0)\cdot(\boldsymbol{F}_1+\boldsymbol{F}_0)^{-1} L:=Δt2(F1−F0)⋅(F1+F0)−1
应有 D+W=L\boldsymbol{D}+\boldsymbol{W}=\boldsymbol{L}D+W=L

D=DSTRAIN/DTIME;
W=2/DTIME*(DROT-eye(3))/(DROT+eye(3));
L=2/DTIME*(F1-F0)/(F1+F0);D+W=-0.129952543241347  -0.293505063321040  -0.0440265164519600.000000000002498   0.518339598294507  -0.000000000000382-0.148874070703533  -0.028576680247064  -0.387625508999842L=-0.129950746561754  -0.293503500675591  -0.0440269566853450   0.518338007574153                   0-0.148875559349057  -0.028577393730105  -0.387626288891340