abaqus umat介电弹性体Mooney-Rivlin model + electrostatic

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)CINCLUDE 'ABA_PARAM.INC'CCHARACTER*8 CMNAMEDIMENSION 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)CDIMENSION BBAR(6),DISTGR(3,3),BBARBBAR(6)CPARAMETER(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)CC ----------------------------------------------------------------C UMAT FOR COMPRESSIBLE MOONEY-RIVLIN MODELC W=C10(I-3)+C01(I-3)+1/D1 (J-1)*(J-1)C ----------------------------------------------------------------C PROPS(1) - C10 C1C PROPS(2) - C01 C2C PROPS(3) - D1 Bulk modulusC PROPS(4) - E1 E-field in X directionC PROPS(5) - E2 E-field in Y directionC PROPS(6) - E3 E-field in Z directionC PROPS(7) – EPSILON Vaccum permittivity*relative permittivityC ----------------------------------------------------------------CC Read Material PROPERTIESCC10=PROPS(1)C01=PROPS(2)D1 =PROPS(3)EPSILON=PROPS(7)CC JACOBIAN AND DISTORTION TENSORCDET=DFGRD1(1, 1)*DFGRD1(2, 2)*DFGRD1(3, 3)1 -DFGRD1(1, 2)*DFGRD1(2, 1)*DFGRD1(3, 3)IF(NSHR.EQ.3) THENDET=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 IFSCALE=DET**(-ONE/THREE)DO K1=1, 3DO K2=1, 3DISTGR(K2, K1)=SCALE*DFGRD1(K2, K1)END DOEND DOCC CALCULATE LEFT CAUCHY-GREEN TENSOR (B=F.Ft)CBBAR(1)=DISTGR(1, 1)**2+DISTGR(1, 2)**2+DISTGR(1, 3)**2BBAR(2)=DISTGR(2, 1)**2+DISTGR(2, 2)**2+DISTGR(2, 3)**2BBAR(3)=DISTGR(3, 3)**2+DISTGR(3, 1)**2+DISTGR(3, 2)**2BBAR(4)=DISTGR(1, 1)*DISTGR(2, 1)+DISTGR(1, 2)*DISTGR(2, 2)1 +DISTGR(1, 3)*DISTGR(2, 3)IF(NSHR.EQ.3) THENBBAR(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 IFCC CALCULATAE (B.Bt)CBBARBBAR(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)CC CALCULATE THE STRESSCTRBBAR=(BBAR(1)+BBAR(2)+BBAR(3))TRBBAR1=(BBAR(1)+BBAR(2)+BBAR(3))/THREETRBBAR2=(BBAR(1)**2+BBAR(2)**2+BBAR(3)**2)/THREETRBBARBBAR=(BBARBBAR(1)+BBARBBAR(2)+BBARBBAR(3))/THREEEG1=TWO*C10/DETEG2=TWO*C01/DETEK=TWO/D1*(TWO*DET-ONE)PR=TWO/D1*(DET-ONE)DO K1=1,NDISTRESS(K1)=EG1*(BBAR(K1)-TRBBAR1)+PR1 +EG2*(TRBBAR*BBAR(K1)-TRBBAR2-BBARBBAR(K1)2 +TRBBARBBAR)END DODO K1=NDI+1,NDI+NSHRSTRESS(K1)=EG1*BBAR(K1)1 +EG2*(TRBBAR*BBAR(K1)-BBARBBAR(K1))END DOcC CALCULATE THE STIFFNESSCEG23=EG1*TWO/THREEEG24=EG2*TWO/THREEDDSDDE(1, 1)= EG23*(BBAR(1)+TRBBAR1)+EK1 +EG24*(-(BBAR(1)*(3*BBAR(1)+5*(BBAR(2)+BBAR(3))))/22 +BBAR(4)**2+BBAR(5)**2)DDSDDE(2, 2)= EG23*(BBAR(2)+TRBBAR1)+EK1 +EG24*(-(BBAR(2)*(5*BBAR(1)+3*BBAR(2)+5*BBAR(3)))/22 +BBAR(4)**2+BBAR(6)**2)DDSDDE(3, 3)= EG23*(BBAR(3)+TRBBAR1)+EK1 +EG24*(-(BBAR(3)*(5*(BBAR(1)+BBAR(2))+3*BBAR(3)))/22 +BBAR(5)**2+BBAR(6)**2)DDSDDE(1, 2)=-EG23*(BBAR(1)+BBAR(2)-TRBBAR1)+EK1 +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)+EK1 +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)+EK1 +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)/TWO1 +EG24*(-((3*(BBAR(1)+BBAR(2))+5*BBAR(3))*BBAR(4))/42 +(BBAR(5)*BBAR(6))/2)DDSDDE(2, 4)= EG23*BBAR(4)/TWO1 +EG24*(-((3*(BBAR(1)+BBAR(2))+5*BBAR(3))*BBAR(4))/42 +(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)**22 +3*(BBAR(1)+BBAR(2))*BBAR(3)+2*(16*BBAR(4)**23 +5*(BBAR(5)**2+BBAR(6)**2)))/8)IF(NSHR.EQ.3) THENDDSDDE(1, 5)= EG23*BBAR(5)/TWO1 +EG24*(-((3*BBAR(1)+5*BBAR(2)+3*BBAR(3))*BBAR(5))/42 +(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)/TWO1 +EG24*(-((3*BBAR(1)+5*BBAR(2)+3*BBAR(3))*BBAR(5))/42 +(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)/TWO1 +EG24*((2*BBAR(4)*BBAR(5)-(5*BBAR(1)+3*(BBAR(2)2 +BBAR(3)))*BBAR(6))/4)DDSDDE(3, 6)= EG23*BBAR(6)/TWO1 +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))/TWO1 +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)**23 +32*BBAR(5)**2+10*BBAR(6)**2)/8)DDSDDE(6, 6)= EG1*(BBAR(2)+BBAR(3))/TWO1 +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)/TWO1 +EG24*((22*BBAR(4)*BBAR(5)+(-9*BBAR(1)+13*(BBAR(2)2 +BBAR(3)))*BBAR(6))/8)DDSDDE(4,6)= EG1*BBAR(5)/TWO1 +EG24*((-3*(BBAR(1)+3*BBAR(2)+BBAR(3))*BBAR(5))/82 +(3*BBAR(4)*BBAR(6))/4)DDSDDE(5,6)= EG1*BBAR(4)/TWO1 +EG24*(((13*(BBAR(1)+BBAR(2))-9*BBAR(3))*BBAR(4)2 +22*BBAR(5)*BBAR(6))/8)END IFDO K1=1, NTENSDO K2=1, K1-1DDSDDE(K1, K2)=DDSDDE(K2, K1)END DOEND DOCC Calculate the inverse of deformation gradientCDFGRDM1_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,3DO J=1,3DFGRDM1_INV(I,J)=DFGRDM1_INV(I,J)/DETEND DOEND DOCC Calculate the current true electric fieldCEn(1)=PROPS(14)En(2)=PROPS(15)En(3)=PROPS(16)DO I=1,3Et(I)=0END DODO I=1,3DO J=1,3Et(J)=Et(J)+En(I)*DFGRDM1_INV(I,J)END DOEND DOEkEk=0.0D0DO I=1, 3EkEk=EkEk+Et(I)*Et(I)ENDDODO I=1,3DO J=1,3Stress_Max(I,J)=0.0D0ENDDOENDDOCC Update the true stress due to polarizationCDO I=1,3DO J=1,3IF (I==J) THENStress_Max(I,J)=EPSILON*Et(I)*Et(J)-0.50D0*EPSILON*EkEk(1)ELSEStress_Max(I,J)=EPSILON*Et(I)*Et(J)ENDIFENDDOENDDOSTRESS(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)RETURNEND