RVE预测力学性能

ZengPeng · 发表于 2023-9-6 16:38:19

目前我所遇到的问题是，发生损伤后，应力依然不减！希望大佬能帮忙看看

SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATV),
1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
4 JSTEP(4)
C DDSDDE-雅可比矩阵 SSE-弹性应变能 SPD-塑性耗能 SCD-蠕变耗能
C 必须定义，应力stress;状态变量statev；雅可比矩阵DDSDDE
if (CMNAME .EQ. 'FIBER') THEN
CALL UMAT_FIBER(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
else if (CMNAME .EQ. 'MATRIX') THEN
CALL UMAT_MATRIX(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
end if
RETURN
END
C **********************************************************
C *subroutines are divided *
C **********************************************************
SUBROUTINE UMAT_FIBER(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATV),
1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
4 JSTEP(4)
C read matrial porperty
parameter(zero = 0.d0, one = 1.d0, two = 2.d0)
DIMENSION STRANT(6)
E1 = PROPS(1)
E2 = PROPS(2)
E3 = PROPS(3)
G12 = PROPS(4)
G13 = PROPS(5)
G23 = PROPS(6)
NU12 = PROPS(7)
NU13 = PROPS(8)
NU23 = PROPS(9)
Xt = PROPS(10)
Xc = PROPS(11)
G_IC_f = PROPS(12)
S = PROPS(13)
ft = STATEV(1)
fc = STATEV(2)
fs1t = STATEV(3)
fs2t = STATEV(4)
fs3t = STATEV(5)
DF = STATEV(6)
C update DDSDDE
NU21 = NU12
NU32 = NU23
NU31 = NU13
C DE1 = (one - DF)*E1
C DE2 = (one - DF)*E2
C DE3 = (one - DF)*E3
C DG12 = (one-DF)*G12
C DG13 = (one-DF)*G13
C DG23 = (one-DF)*G23
C DNU12 = (one-DF)*NU12
C DNU13 = (one-DF)*NU13
C DNU23 = (one-DF)*NU23
C DNU21 = DNU12
C DNU31 = DNU13
C DNU32 = DNU23
gg = one / (one - NU12*NU21 - NU23*NU32 - NU31*NU13
* - two*NU21*NU32*NU13)
C11 = E1 * ( one - NU23*NU32 ) * gg
C22 = E2 * ( one - NU13*NU31 ) * gg
C33 = E3 * ( one - NU12*NU21 ) * gg
C12 = E1 * ( NU21 + NU31*NU23 ) * gg
C13 = E1 * ( NU31 + NU21*NU32 ) * gg
C23 = E2 * ( NU32 + NU12*NU31 ) * gg
C21 = E1 * ( NU21 + NU31*NU23 ) * gg
C31 = E1 * ( NU31 + NU21*NU32 ) * gg
C32 = E2 * ( NU32 + NU12*NU31 ) * gg
C44 = G12
C55 = G13
C66 = G23
do i = 1, NTENS
do j = 1, NTENS
DDSDDE(i,j)= 0
end do
end do
DDSDDE(1,1) = (1-DF)*C11
DDSDDE(2,2) = (1-DF)*C22
DDSDDE(3,3) = (1-DF)*C33
DDSDDE(4,4) = (1-DFS1)*C44
DDSDDE(5,5) = (1-DFS2)*C55
DDSDDE(6,6) = (1-DFS3)*C66
DDSDDE(1,2) = (1-DF)*C12
DDSDDE(2,1) = (1-DF)*C21
DDSDDE(1,3) = (1-DF)*C13
DDSDDE(3,1) = (1-DF)*C31
DDSDDE(2,3) = (1-DF)*C23
DDSDDE(3,2) = (1-DF)*C32
do i=1,NTENS
STRANT(i) = STRAN(i) + DSTRAN(i)
end do
do i = 1, NTENS
do j = 1, NTENS
STRESS(i)=STRESS(i)+DDSDDE(i,j)*DSTRAN(j)
end do
end do
if (STRESS(1) .GE. zero) then
X1 = Xt
ft = STRESS(1)/X1
if (ft .gt. one) then
equivalent_strain_ft0 = Xt/E1
equivalent_strain_ft1 = (two* G_IC_f)/(Xt*celent)
D1 = (equivalent_strain_ft1*(STRANT(1)-
* equivalent_strain_ft0))/
* (STRANT(1)*(equivalent_strain_ft1-
* equivalent_strain_ft0))
else
D1 = zero
end if
else
X2 = -1*Xc
fc = STRESS(1)/X2
if (fc .gt. one) then
C take care of +-
equivalent_strain_fc0 = -one*STRESS(1)/E1
equivalent_strain_fc1 = two* G_IC_f/(Xc*celent)
D2 = (equivalent_strain_fc1*(STRANT(1)-
* equivalent_strain_fc0))/
* (STRANT(1)*(equivalent_strain_fc1-
* equivalent_strain_fc0))
else
D2 = zero
end if
end if
DF = one-(one-D1)*(one-D2)
if(STRESS(4) .GT. zero) then
fs1t = STRESS(4)/S
if (fs1t .GE. one) then
equivalent_strain_ft0 = Xt/E1
equivalent_strain_ft1 = two* G_IC_f/(Xt*celent)
DS1 = (equivalent_strain_ft1*(STRANT(1)-
* equivalent_strain_ft0))/
* (STRANT(1)*(equivalent_strain_ft1-
* equivalent_strain_ft0))
else
DS1 = zero
end if
else
S = -one*S
fs1c = STRESS(4)/S
if (fs1c .gt. one) then
equivalent_strain_fc0 = STRESS(1)/E1
equivalent_strain_fc1 = two* G_IC_f/(Xc*celent)
DS2 = (equivalent_strain_fc1*(STRANT(1)-
* equivalent_strain_fc0))/
* (STRANT(1)*(equivalent_strain_fc1-
* equivalent_strain_fc0))
else
DS2 = zero
end if
end if
DFS1 = one-(one-DS1)*(one-DS2)
if(STRESS(5) .GT. zero) then
fs2t = STRESS(5)/S
if (fs2t .GT. one) then
equivalent_strain_ft0 = Xt/E1
equivalent_strain_ft1 = two* G_IC_f/(Xt*celent)
DS3 = (equivalent_strain_ft1*(STRANT(1)-
* equivalent_strain_ft0))/
* (STRANT(1)*(equivalent_strain_ft1-
* equivalent_strain_ft0))
else
DS3 = zero
end if
else
S = -one*S
fs2c = STRESS(4)/S
if (fs2c .gt. one) then
equivalent_strain_fc0 = STRESS(1)/E1
equivalent_strain_fc1 = two* G_IC_f/(Xc*celent)
DS4 = (equivalent_strain_fc1*(STRANT(1)-
* equivalent_strain_fc0))/
* (STRANT(1)*(equivalent_strain_fc1-
* equivalent_strain_fc0))
else
DS4 = zero
end if
end if
DFS2 = one-(one-DS3)*(one-DS4)
if(STRESS(6) .GT. zero) then
fs3t = STRESS(6)/S
if (fs3t .GT. one) then
equivalent_strain_ft0 = Xt/E1
equivalent_strain_ft1 = two* G_IC_f/(Xt*celent)
DS5 = (equivalent_strain_ft1*(STRANT(1)-
* equivalent_strain_ft0))/
* (STRANT(1)*(equivalent_strain_ft1-
* equivalent_strain_ft0))
else
DS5 = zero
end if
else
S = -one*S
fs3c = STRESS(6)/S
if (fs3c .gt. one) then
equivalent_strain_fc0 = S11/E1
equivalent_strain_fc1 = two* G_IC_f/(Xc*celent)
DS6 = (equivalent_strain_fc1*(STRANT(1)-
* equivalent_strain_fc0))/
* (STRANT(1)*(equivalent_strain_fc1-
* equivalent_strain_fc0))
else
DS6 = zero
end if
end if
DFS3 = one-(one-DS5)*(one-DS6)
do i = 1, NTENS
do j = 1, NTENS
DDSDDE(i,j)= 0
end do
end do
if (DF .ge. zero) then
DF = max(DF,STATEV(6))
END IF
DDSDDE(1,1) = (1-DF)*C11
DDSDDE(2,2) = (1-DF)*C22
DDSDDE(3,3) = (1-DF)*C33
DDSDDE(4,4) = (1-DFS1)*C44
DDSDDE(5,5) = (1-DFS2)*C55
DDSDDE(6,6) = (1-DFS3)*C66
DDSDDE(1,2) = (1-DF)*C12
DDSDDE(2,1) = (1-DF)*C21
DDSDDE(1,3) = (1-DF)*C13
DDSDDE(3,1) = (1-DF)*C31
DDSDDE(2,3) = (1-DF)*C23
DDSDDE(3,2) = (1-DF)*C32
do i = 1, NTENS
do j = 1, NTENS
STRESS(i)=STRESS(i)+DDSDDE(i,j)*DSTRAN(j)
end do
end do
STATEV(1) = ft
STATEV(2) = fc
STATEV(3) = fs1t
STATEV(4) = fs2t
STATEV(5) = fs3t
STATEV(6) = DF
RETURN
END
C *************************************
C *************************************
SUBROUTINE UMAT_MATRIX(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATV),
1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
4 JSTEP(4)
parameter(zero = 0.d0, one = 1.d0, two = 2.d0, half = 0.5d0,
* three = 3.d0)
DIMENSION STRANT(6)
Em = PROPS(1)
Gm = PROPS(2)
NUm = PROPS(3)
Xmt = PROPS(4)
Xmc = PROPS(5)
Gms = PROPS(6)
G_IC_m = PROPS(7)
amt = STATEV(5)
amc = STATEV(6)
Dmt = STATEV(7)
Dmc = STATEV(8)
ADM = STATEV(9)
c DEm = (one-ADM)*Em
C DGm = (one-ADM)*Gm
C DNUm = (one-ADM)*NUm
C computer stiffiness matrix
EG=Em/(two*(one + NUm))
ELAM=(NUm*Em)/((one + NUm)*(one - two*NUm))
AC11 = TWO*EG+ELAM
AC22 = TWO*EG+ELAM
AC33 = TWO*EG+ELAM
AC44 = Gm/two
AC55 = Gm/two
AC66 = Gm/two
AC12 = ELAM
AC13 = ELAM
AC21 = ELAM
AC23 = ELAM
AC31 = ELAM
AC32 = ELAM
C ***************************
C reset DDSDDE
do i = 1, NTENS
do j = 1, NTENS
DDSDDE(i,j)= 0
end do
end do
DDSDDE(1,1) = (one - ADM)*AC11
DDSDDE(2,2) = (one - ADM)*AC22
DDSDDE(3,3) = (one - ADM)*AC33
DDSDDE(4,4) = (one - ADM)*AC44
DDSDDE(5,5) = (one - ADM)*AC55
DDSDDE(6,6) = (one - ADM)*AC66
DDSDDE(1,2) = (one - ADM)*AC12
DDSDDE(2,1) = (one - ADM)*AC21
DDSDDE(1,3) = (one - ADM)*AC13
DDSDDE(3,1) = (one - ADM)*AC31
DDSDDE(2,3) = (one - ADM)*AC23
DDSDDE(3,2) = (one - ADM)*AC32
C update strain
do i=1,NTENS
STRANT(i) = STRAN(i) + DSTRAN(i)
end do
C update stress
do i = 1, NTENS
do j = 1, NTENS
STRESS(i)=STRESS(i)+DDSDDE(i,j)*DSTRAN(j)
end do
end do
C caculate parametry in matrix damage
if (STRESS(1) .ge. 0) then
if (STRESS(2) .le. 0) then
a = 0
else
a = STRESS(2)
end if
if (STRANT(2) .le. 0) then
b = 0
else
b = STRANT(2)
end if
else
if ((-STRESS(2)) .le. 0) then
c = 0
else
c = -STRESS(2)
end if
if ((-STRANT(2)) .le. 0) then
d = 0
else
d = -STRANT(2)
end if
end if
C **************************
if (STRESS(1) .GE. zero) then
AX1 = Xmt
else
AX1 = -Xmt
end if
if (STRESS(3) .LT. zero) then
AX2 = -one*Xmc
else
AX2 = Xmc
end if
amt = STRESS(1)/AX1
amc = STRESS(3)/Xmc
if (amt .ge. one) then
equivalent_stress_mt = (celent*(a*b)+STRESS(4)*STRANT(4)+
* STRESS(6)*STRANT(6))/
* (celent*((b**two+STRANT(4)**two+
* STRANT(6)**two)**half))
equivalent_strain_mt1 = (two*G_IC_m*(amt**half))/
* equivalent_stress_mt
equivalent_strain_mt = celent*((b**two+STRANT(4)**two+
* STRANT(6)**two)**half)
equivalent_strain_mt0 = equivalent_strain_mt/(amt**half)
Dmt = (equivalent_strain_mt1*(equivalent_strain_mt-
* equivalent_strain_mt0))/
* (equivalent_strain_mt*(equivalent_strain_mt1-
* equivalent_strain_mt0))
else
Dmt = zero
end if
c write(*,*) 'Dmt=',Dmt
c write(*,*) 'mt1=',equivalent_strain_mt1
c write(*,*) 'mt=',equivalent_strain_mt
c write(*,*) 'mt0=',equivalent_strain_mt0
c write(*,*) 'fenmu=',equivalent_strain_mt*(equivalent_strain_mt1-
c * equivalent_strain_mt0)
c WRITE(*,*) 'FENZI2=',equivalent_strain_mt-equivalent_strain_mt0
c WRITE(*,*)'FENZI=',equivalent_strain_mt1*(equivalent_strain_mt-
c * equivalent_strain_mt0)
C WRITE(*,*)'FENMU=',equivalent_strain_mt1-
C * equivalent_strain_mt0
if (amc .ge. one) then
equivalent_strain_mc = celent*((d**two+STRANT(4)**two)**half)
1 /(amc**half)
ABC = (celent*c*d+STRESS(4)*STRANT(4))
equivalent_stress_mc = ABC/equivelent_strain_mc
equivalent_strain_mc0 = equivalent_strain_mc/(amc**two)
equivalent_strain_mc1 = (two*G_IC_m*(amc**half))/
* (equivalent_stress_mc)
Dmc = (equivalent_strain_mc1*(equivalent_strain_mc-
* equivalent_strain_mc0))/
* (equivalent_strain_mc*(equivalent_strain_mc1-
* equivalent_strain_mc0))
else
Dmc = zero
end if
c write(*,*) 'Dmc=',Dmc
C update stiffnesee matrix
C reset DDSDDE
ADM = one-(one-Dmt)*(one-Dmc)
c DEG=((one-ADF)*Em)/(two*(one + (one-ADF)*NUm))
c DELAM=((one-ADF)*NUm*(one-ADF)*Em)/((one + (one-ADF)*NUm)*
c * (one - two*(one-ADF)*NUm))
c AC11 = TWO*DEG+DELAM
c AC22 = TWO*DEG+DELAM
c AC33 = TWO*DEG+DELAM
c AC44 = DEG
c AC55 = DEG
c AC66 = DEG
c AC12 = DELAM
c AC13 = DELAM
c AC21 = DELAM
c AC23 = DELAM
c AC31 = DELAM
c AC32 = DELAM
C computer stiffiness matrix
do i = 1, NTENS
do j = 1, NTENS
DDSDDE(i,j)= 0
end do
end do
IF (ADM .gt. STATEV(9)) then
ADM = ADM
else
ADM = STATEV(9)
end if
DDSDDE(1,1) = (one - ADM)*AC11
DDSDDE(2,2) = (one - ADM)*AC22
DDSDDE(3,3) = (one - ADM)*AC33
DDSDDE(4,4) = (one - ADM)*AC44
DDSDDE(5,5) = (one - ADM)*AC55
DDSDDE(6,6) = (one - ADM)*AC66
DDSDDE(1,2) = (one - ADM)*AC12
DDSDDE(2,1) = (one - ADM)*AC21
DDSDDE(1,3) = (one - ADM)*AC13
DDSDDE(3,1) = (one - ADM)*AC31
DDSDDE(2,3) = (one - ADM)*AC23
DDSDDE(3,2) = (one - ADM)*AC32
C recompute stress
do i = 1, NTENS
do j = 1, NTENS
STRESS(i)=STRESS(i)+DDSDDE(i,j)*DSTRAN(j)
end do
end do
STATEV(5) = amt
STATEV(6) = amc
STATEV(7) = Dmt
STATEV(8) = Dmc
STATEV(9) = ADM
RETURN
END

复制代码

ZengPeng · 发表于 2023-9-6 16:42:45

后面，我在计算应力是加上STRESS(I) = 0后，能完全运算，但是得到的效果就是，所有的单元应力都是同时到达损伤，这明显和实际不符合。为什么不能出现那种某个地方先出现损伤，然后损伤逐渐扩大的情况呢，希望有大佬指点一二。

avatarr · 发表于 2023-12-7 12:38:39

请问楼主的问题解决了嘛

账号		自动登录	找回密码
密码			注册

账号		自动登录	找回密码
密码			立即注册

[子程序] RVE预测力学性能