C Mooney-Rivlin model + electrostatic energy function (E_MRmodel.f)
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,NSTATEV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*8 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATEV),
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 Et(3),Stress_Max(3,3),EkEk(1),DFGRDM1_INV(3,3),En(3),
5 TEMP(1),DTEMP(1),DFGRDP(3),DFGRDM1(3, 3),DFGRDM0(3, 3),
6 CBAR(3,3),DDSDDE_Max(6,6)
C
DIMENSION BBAR(6),DISTGR(3,3),BBARBBAR(6)
C
PARAMETER(ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0, FOUR=4.D0,
1 FIVE=5.D0, SIX=6.D0, SEVEN=7.D0, EIGHT=8.D0)
C
C ----------------------------------------------------------------
C UMAT FOR COMPRESSIBLE MOONEY-RIVLIN MODEL
C W=C10(I-3)+C01(I-3)+1/D1 (J-1)*(J-1)
C ----------------------------------------------------------------
C PROPS(1) - C10 C1
C PROPS(2) - C01 C2
C PROPS(3) - D1 Bulk modulus
C PROPS(4) - E1 E-field in X direction
C PROPS(5) - E2 E-field in Y direction
C PROPS(6) - E3 E-field in Z direction
C PROPS(7) – EPSILON Vaccum permittivity*relative permittivity
C ----------------------------------------------------------------
C
C Read Material PROPERTIES
C
C10=PROPS(1)
C01=PROPS(2)
D1 =PROPS(3)
EPSILON=PROPS(7)
C
C JACOBIAN AND DISTORTION TENSOR
C
DET=DFGRD1(1, 1)*DFGRD1(2, 2)*DFGRD1(3, 3)
1 -DFGRD1(1, 2)*DFGRD1(2, 1)*DFGRD1(3, 3)
IF(NSHR.EQ.3) THEN
DET=DET+DFGRD1(1, 2)*DFGRD1(2, 3)*DFGRD1(3, 1)
1 +DFGRD1(1, 3)*DFGRD1(3, 2)*DFGRD1(2, 1)
2 -DFGRD1(1, 3)*DFGRD1(3, 1)*DFGRD1(2, 2)
3 -DFGRD1(2, 3)*DFGRD1(3, 2)*DFGRD1(1, 1)
END IF
SCALE=DET**(-ONE/THREE)
DO K1=1, 3
DO K2=1, 3
DISTGR(K2, K1)=SCALE*DFGRD1(K2, K1)
END DO
END DO
C
C CALCULATE LEFT CAUCHY-GREEN TENSOR (B=F.Ft)
C
BBAR(1)=DISTGR(1, 1)**2+DISTGR(1, 2)**2+DISTGR(1, 3)**2
BBAR(2)=DISTGR(2, 1)**2+DISTGR(2, 2)**2+DISTGR(2, 3)**2
BBAR(3)=DISTGR(3, 3)**2+DISTGR(3, 1)**2+DISTGR(3, 2)**2
BBAR(4)=DISTGR(1, 1)*DISTGR(2, 1)+DISTGR(1, 2)*DISTGR(2, 2)
1 +DISTGR(1, 3)*DISTGR(2, 3)
IF(NSHR.EQ.3) THEN
BBAR(5)=DISTGR(1, 1)*DISTGR(3, 1)+DISTGR(1, 2)*DISTGR(3, 2)
1 +DISTGR(1, 3)*DISTGR(3, 3)
BBAR(6)=DISTGR(2, 1)*DISTGR(3, 1)+DISTGR(2, 2)*DISTGR(3, 2)
1 +DISTGR(2, 3)*DISTGR(3, 3)
END IF
C
C CALCULATAE (B.Bt)
C
BBARBBAR(1)=BBAR(1)*BBAR(1)+BBAR(4)*BBAR(4)+BBAR(5)*BBAR(5)
BBARBBAR(2)=BBAR(4)*BBAR(4)+BBAR(2)*BBAR(2)+BBAR(6)*BBAR(6)
BBARBBAR(3)=BBAR(5)*BBAR(5)+BBAR(6)*BBAR(6)+BBAR(3)*BBAR(3)
BBARBBAR(4)=BBAR(1)*BBAR(4)+BBAR(4)*BBAR(2)+BBAR(5)*BBAR(6)
BBARBBAR(5)=BBAR(1)*BBAR(5)+BBAR(4)*BBAR(6)+BBAR(5)*BBAR(3)
BBARBBAR(6)=BBAR(4)*BBAR(5)+BBAR(2)*BBAR(6)+BBAR(6)*BBAR(3)
C
C CALCULATE THE STRESS
C
TRBBAR=(BBAR(1)+BBAR(2)+BBAR(3))
TRBBAR1=(BBAR(1)+BBAR(2)+BBAR(3))/THREE
TRBBAR2=(BBAR(1)**2+BBAR(2)**2+BBAR(3)**2)/THREE
TRBBARBBAR=(BBARBBAR(1)+BBARBBAR(2)+BBARBBAR(3))/THREE
EG1=TWO*C10/DET
EG2=TWO*C01/DET
EK=TWO/D1*(TWO*DET-ONE)
PR=TWO/D1*(DET-ONE)
DO K1=1,NDI
STRESS(K1)=EG1*(BBAR(K1)-TRBBAR1)+PR
1 +EG2*(TRBBAR*BBAR(K1)-TRBBAR2-BBARBBAR(K1)
2 +TRBBARBBAR)
END DO
DO K1=NDI+1,NDI+NSHR
STRESS(K1)=EG1*BBAR(K1)
1 +EG2*(TRBBAR*BBAR(K1)-BBARBBAR(K1))
END DO
c
C CALCULATE THE STIFFNESS
C
EG23=EG1*TWO/THREE
EG24=EG2*TWO/THREE
DDSDDE(1, 1)= EG23*(BBAR(1)+TRBBAR1)+EK
1 +EG24*(-(BBAR(1)*(3*BBAR(1)+5*(BBAR(2)+BBAR(3))))/2
2 +BBAR(4)**2+BBAR(5)**2)
DDSDDE(2, 2)= EG23*(BBAR(2)+TRBBAR1)+EK
1 +EG24*(-(BBAR(2)*(5*BBAR(1)+3*BBAR(2)+5*BBAR(3)))/2
2 +BBAR(4)**2+BBAR(6)**2)
DDSDDE(3, 3)= EG23*(BBAR(3)+TRBBAR1)+EK
1 +EG24*(-(BBAR(3)*(5*(BBAR(1)+BBAR(2))+3*BBAR(3)))/2
2 +BBAR(5)**2+BBAR(6)**2)
DDSDDE(1, 2)=-EG23*(BBAR(1)+BBAR(2)-TRBBAR1)+EK
1 +EG24*((3*(2*BBAR(1)*BBAR(2)-(4*(BBAR(1)+BBAR(2))
2 *(BBAR(1)+BBAR(2)+BBAR(3)))/3-2*BBAR(4)**2))/2)
DDSDDE(1, 3)=-EG23*(BBAR(1)+BBAR(3)-TRBBAR1)+EK
1 +EG24*((3*(2*BBAR(1)*BBAR(3)-(4*(BBAR(1)+BBAR(3))
2 *(BBAR(1)+BBAR(2)+BBAR(3)))/3-2*BBAR(5)**2))/2)
DDSDDE(2, 3)=-EG23*(BBAR(2)+BBAR(3)-TRBBAR1)+EK
1 +EG24*((3*(2*BBAR(2)*BBAR(3)-(4*(BBAR(2)+BBAR(3))
2 *(BBAR(1)+BBAR(2)+BBAR(3)))/3-2*BBAR(6)**2))/2)
DDSDDE(1, 4)= EG23*BBAR(4)/TWO
1 +EG24*(-((3*(BBAR(1)+BBAR(2))+5*BBAR(3))*BBAR(4))/4
2 +(BBAR(5)*BBAR(6))/2)
DDSDDE(2, 4)= EG23*BBAR(4)/TWO
1 +EG24*(-((3*(BBAR(1)+BBAR(2))+5*BBAR(3))*BBAR(4))/4
2 +(BBAR(5)*BBAR(6))/2)
DDSDDE(3, 4)=-EG23*BBAR(4)
1 +EG24*((-2*(BBAR(1)+BBAR(2))+BBAR(3))*BBAR(4)
2 -3*BBAR(5)*BBAR(6))
DDSDDE(4, 4)= EG1*(BBAR(1)+BBAR(2))
1 +EG24*((13*BBAR(1)**2-6*BBAR(1)*BBAR(2)+13*BBAR(2)**2
2 +3*(BBAR(1)+BBAR(2))*BBAR(3)+2*(16*BBAR(4)**2
3 +5*(BBAR(5)**2+BBAR(6)**2)))/8)
IF(NSHR.EQ.3) THEN
DDSDDE(1, 5)= EG23*BBAR(5)/TWO
1 +EG24*(-((3*BBAR(1)+5*BBAR(2)+3*BBAR(3))*BBAR(5))/4
2 +(BBAR(4)*BBAR(6))/2)
DDSDDE(2, 5)=-EG23*BBAR(5)
1 +EG24*((-2*BBAR(1)+BBAR(2)-2*BBAR(3))*BBAR(5)
2 -3*BBAR(4)*BBAR(6))
DDSDDE(3, 5)= EG23*BBAR(5)/TWO
1 +EG24*(-((3*BBAR(1)+5*BBAR(2)+3*BBAR(3))*BBAR(5))/4
2 +(BBAR(4)*BBAR(6))/2)
DDSDDE(1, 6)=-EG23*BBAR(6)
1 +EG24*(-3*BBAR(4)*BBAR(5)+(BBAR(1)-2*(BBAR(2)
2 +BBAR(3)))*BBAR(6))
DDSDDE(2, 6)= EG23*BBAR(6)/TWO
1 +EG24*((2*BBAR(4)*BBAR(5)-(5*BBAR(1)+3*(BBAR(2)
2 +BBAR(3)))*BBAR(6))/4)
DDSDDE(3, 6)= EG23*BBAR(6)/TWO
1 +EG24*((2*BBAR(4)*BBAR(5)-(5*BBAR(1)+3*(BBAR(2)
2 +BBAR(3)))*BBAR(6))/4)
DDSDDE(5, 5)= EG1*(BBAR(1)+BBAR(3))/TWO
1 +EG24*((13*BBAR(1)**2+3*BBAR(1)*(BBAR(2)-2*BBAR(3))
2 +3*BBAR(2)*BBAR(3)+13*BBAR(3)**2+10*BBAR(4)**2
3 +32*BBAR(5)**2+10*BBAR(6)**2)/8)
DDSDDE(6, 6)= EG1*(BBAR(2)+BBAR(3))/TWO
1 +EG24*((13*BBAR(2)**2-6*BBAR(2)*BBAR(3)
2 +13*BBAR(3)**2+3*BBAR(1)*(BBAR(2)+BBAR(3))
3 +10*(BBAR(4)**2+BBAR(5)**2)+32*BBAR(6)**2)/8)
DDSDDE(4,5)= EG1*BBAR(6)/TWO
1 +EG24*((22*BBAR(4)*BBAR(5)+(-9*BBAR(1)+13*(BBAR(2)
2 +BBAR(3)))*BBAR(6))/8)
DDSDDE(4,6)= EG1*BBAR(5)/TWO
1 +EG24*((-3*(BBAR(1)+3*BBAR(2)+BBAR(3))*BBAR(5))/8
2 +(3*BBAR(4)*BBAR(6))/4)
DDSDDE(5,6)= EG1*BBAR(4)/TWO
1 +EG24*(((13*(BBAR(1)+BBAR(2))-9*BBAR(3))*BBAR(4)
2 +22*BBAR(5)*BBAR(6))/8)
END IF
DO K1=1, NTENS
DO K2=1, K1-1
DDSDDE(K1, K2)=DDSDDE(K2, K1)
END DO
END DO
C
C Calculate the inverse of deformation gradient
C
DFGRDM1_INV(1,1)=DFGRD1(2,2)*DFGRD1(3,3)
1 -DFGRD1(2,3)*DFGRD1(3,2)
DFGRDM1_INV(1,2)=-DFGRD1(1,2)*DFGRD1(3,3)
1 +DFGRD1(1,3)*DFGRD1(3,2)
DFGRDM1_INV(1,3)=DFGRD1(1,2)*DFGRD1(2,3)
1 -DFGRD1(1,3)*DFGRD1(2,2)
DFGRDM1_INV(2,1)=-DFGRD1(2,1)*DFGRD1(3,3)
1 +DFGRD1(2,3)*DFGRD1(3,1)
DFGRDM1_INV(2,2)=DFGRD1(1,1)*DFGRD1(3,3)
1 -DFGRD1(1,3)*DFGRD1(3,1)
DFGRDM1_INV(2,3)=-DFGRD1(1,1)*DFGRD1(2,3)
1 +DFGRD1(1,3)*DFGRD1(2,1)
DFGRDM1_INV(3,1)=DFGRD1(2,1)*DFGRD1(3,2)
1 -DFGRD1(2,2)*DFGRD1(3,1)
DFGRDM1_INV(3,2)=-DFGRD1(1,1)*DFGRD1(3,2)
1 +DFGRD1(1,2)*DFGRD1(3,1)
DFGRDM1_INV(3,3)=DFGRD1(1,1)*DFGRD1(2,2)
1 -DFGRD1(1,2)*DFGRD1(2,1)
DO I=1,3
DO J=1,3
DFGRDM1_INV(I,J)=DFGRDM1_INV(I,J)/DET
END DO
END DO
C
C Calculate the current true electric field
C
En(1)=PROPS(14)
En(2)=PROPS(15)
En(3)=PROPS(16)
DO I=1,3
Et(I)=0
END DO
DO I=1,3
DO J=1,3
Et(J)=Et(J)+En(I)*DFGRDM1_INV(I,J)
END DO
END DO
EkEk=0.0D0
DO I=1, 3
EkEk=EkEk+Et(I)*Et(I)
ENDDO
DO I=1,3
DO J=1,3
Stress_Max(I,J)=0.0D0
ENDDO
ENDDO
C
C Update the true stress due to polarization
C
DO I=1,3
DO J=1,3
IF (I==J) THEN
Stress_Max(I,J)=EPSILON*Et(I)*Et(J)-0.50D0*EPSILON*EkEk(1)
ELSE
Stress_Max(I,J)=EPSILON*Et(I)*Et(J)
ENDIF
ENDDO
ENDDO
STRESS(1)=STRESS(1)+Stress_Max(1,1)
STRESS(2)=STRESS(2)+Stress_Max(2,2)
STRESS(3)=STRESS(3)+Stress_Max(3,3)
STRESS(4)=STRESS(4)+Stress_Max(1,2)
STRESS(5)=STRESS(5)+Stress_Max(1,3)
STRESS(6)=STRESS(6)+Stress_Max(2,3)
RETURN
END
{var%20f='http://v.t.sina.com.cn/share/share.php?appkey=1515056452',u=z||d.location,p=['&url=',e(u),'&title=',e(t||d.title),'&source=',e(r),'&sourceUrl=',e(l),'&content=',c||'gb2312','&pic=',e(p||'')].join('');function%20a(){if(!window.open([f,p].join(''),'mb',['toolbar=0,status=0,resizable=1,width=440,height=430,left=',(s.width-440)/2,',top=',(s.height-430)/2].join('')))u.href=[f,p].join('');};if(/Firefox/.test(navigator.userAgent))setTimeout(a,0);else%20a();})(screen,document,encodeURIComponent,'','','https://www.xiaopingtou.cn//data/attach/topic/topicKPo7gB.jpg', '推荐 mhkeil51 的文章《abaqus umat介电弹性体Mooney-Rivlin model + electrostati》','https://www.xiaopingtou.net/article-56513.html','页面编码gb2312|utf-8默认gb2312'));){kind=link}