! ====================================================================== ! THE UMAT TESTING PROGRAM ! 2012 (c) tomasz.garbowski@put.poznan.pl, aleksander.marek.pl@gmail.com ! ====================================================================== ! ! UNCOMMENT THE FOLLOWING LINES FOR UMAT TESTING ! program test ! ! implicit none ! integer, parameter :: nprops=20, ntens=4, nstatv=9, ndi=3, nshr=1 !! real*8 :: stress(ntens) = (/ -4.39484d-02, -1.37232d0, !! 2 -1.75092d0 /) ! real*8 :: stress(ntens) = (/ 0d-02, 0d0, ! 2 -0d0, 0d0 /) ! real*8 :: statev(nstatv) ! real*8 :: ddsdde(ntens,ntens) ! real*8 :: stran(ntens) ! real*8 :: dstran(ntens) = (/ 9.6401d-03, -6.1055d-03, 0d0, ! 2 -2.7975d-2 /) ! real*8 :: props(nprops) = (/ 1.d2, 2.d2, 3.d2, 2.d-1, 25.d-2, ! 1 3.d-1, 8.d1, 6.5d1, 4.d1, 1.3d0, 0.d1, 0.d1, 0.d1, 0.d1, 0.d1, ! 2 0.d1, 0.d1, 0.d1, ! 3 1.d2, 1.d0 /) ! real*8 :: pnewdt = 1 ! ! call umat(STRESS,STATEV,DDSDDE,STRAN,DSTRAN, ! & NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,PNEWDT) ! ! end program ! ====================================================================== ! if test program is ON the following lines should be uncomment ! ! subroutine umat(stress,statev,ddsdde,stran,dstran, ! & ndi,nshr,ntens,nstatv,props,nprops,pnewdt) ! ! implicit real*8 (a-h, o-z) ! ---------------------------------------------------------------------- ! if test program is OFF the following lines should be uncomment 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,KSTEP,KINC) INCLUDE 'ABA_PARAM.INC' 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) ! START DIMENSION Cun(6,6), Cut(NTENS,NTENS), Pun(6,6), Qun(6), C(9), 1 D(NTENS,NTENS), P(NTENS,NTENS), Q(NTENS), 2 s_star(NTENS), x(NTENS+1), x2(NTENS), G(NTENS), 3 R1(NTENS), R2(NTENS), R3(NTENS,NTENS), Z(NTENS,NTENS), 4 F(NTENS+1), dFdX(NTENS+1,NTENS+1), dXdF(NTENS+1,NTENS+1), 5 dX(NTENS+1), dFdE(NTENS+1,NTENS), dXdE(NTENS+1,NTENS) real*8, parameter :: ONE=1.D0, TWO=2.D0, ZERO=0.D0, tol=1.5e-08 ! ********************************************************************** ! USER-DEFINED PROPERTIES (NPROPS=20): ! ! DESCRIPTIONS: ! ! Elastic properties ! PROPS(1) : E1 (Elastic modulus in MD /MPa) ! PROPS(2) : E2 (Elastic modulus in CD /MPa) ! PROPS(3) : E3 (Elastic modulus in ZD /MPa) ! PROPS(4) : mu12 (In-plane Poisson's ratio) ! PROPS(5) : mu13 (In-plane Poisson's ratio) ! PROPS(6) : mu23 (In-plane Poisson's ratio) ! PROPS(7) : G12 (In-plane shear modulus /MPa) ! PROPS(8) : G13 ! PROPS(9) : G23 ! Inelastic parameters ! PROPS(10) : sigma0 (reference strength) ! PROPS(10) : ft1 (Hoffman's potential R11+) ! PROPS(11) : ft2 (Hoffman's potential R22+) ! PROPS(12) : ft3 (Hoffman's potential R33+) ! PROPS(13) : fc1 (Hoffman's potential R11-) ! PROPS(14) : fc2 (Hoffman's potential R22-) ! PROPS(15) : fc3 (Hoffman's potential R33-) ! PROPS(16) : fs12 (Hoffman's potential R12) ! PROPS(17) : fs13 (Hoffman's potential R13) ! PROPS(18) : fs23 (Hoffman's potential R22) ! Hardening parameters ! PROPS(19): rQ (Hardening parameter /MPa) ! PROPS(20): rN (exponent) ! ! if props(12) ==0, routine will put ft3=ft2 ! if props(15) ==0, routine will calculate this value to satisfy C7+C8+C9=0) ! if props(19) ==0, rQ = s0/1e-8^rN ! !*********************************************************************** ! USER-DEFINED STATE VARIABLES (NSTATEV=9): ! ! DESCRIPTIONS: ! ! STATEV(1:NTENS): epsilon^el (del = dstran - dep) ! STATEV(NTENS+1:2*NTENS): epsilon^pl (dpl = dlambda*normal) ! STATEV(2*NTENS+1): epsilon^pl_eff (Effective plastic strain [sqrt(2/3*dep'*dep)] ! !*********************************************************************** ! READ IN-PARAMETERS ! if model is simplified to isotropy if (props(2)==0) props(2)=props(1) !e2 if (props(3)==0) props(3)=props(1) !e3 if (props(5)==0) props(5)=props(4) !nu13 if (props(6)==0) props(6)=props(4) !nu23 if (props(7)==0) props(7)=props(1)/2/(1+props(4)) !g12 if (props(8)==0) props(8)=props(1)/2/(1+props(4)) !g13 if (props(9)==0) props(9)=props(1)/2/(1+props(4)) !g23 ! if model is simplified to mises plasticity or hill if (props(11)==0) props(11)=props(10) !ft2 if (props(12)==0) props(12)=props(10) !ft3 if (props(13)==0) props(13)=props(10) !fc1 if (props(14)==0) props(14)=props(11) !fc2 if (props(15)==0) props(15)=props(12) !fc3 if (props(16)==0) props(16)=props(10)/sqrt(3.d0) !fs12 if (props(17)==0) props(17)=props(10)/sqrt(3.d0) !fs13 if (props(18)==0) props(18)=props(10)/sqrt(3.d0) !fs23 print *,props ! Parameters E1=PROPS(1) E2=PROPS(2) E3=PROPS(3) V12=PROPS(4) V13=PROPS(5) V23=PROPS(6) G12=PROPS(7) G13=PROPS(8) G23=PROPS(9) V21=V12*E2/E1 V31=V13*E3/E1 V32=V23*E3/E2 s0=PROPS(10) ft1=PROPS(10)/s0 ft2=PROPS(11)/s0 ft3=PROPS(12)/s0 fc1=PROPS(13)/s0 fc2=PROPS(14)/s0 fc3=PROPS(15)/s0 fs12=PROPS(16)/s0 fs13=PROPS(17)/s0 fs23=PROPS(18)/s0 rQ=PROPS(19) rN=PROPS(20) if (rN>1) rN=ONE/rN if (rQ==ZERO) rQ=s0/(tol**rN) ! ====================================================================== ! material compliance matrix /3D/ Cun(1,:) = (/ ONE/E1, -V21/E2, -V31/E3, ZERO, ZERO, ZERO /) Cun(2,:) = (/ -V12/E1, ONE/E2, -V32/E3, ZERO, ZERO, ZERO /) Cun(3,:) = (/ -V13/E1, -V23/E2, ONE/E3, ZERO, ZERO, ZERO /) Cun(4,:) = (/ ZERO, ZERO, ZERO, ONE/G12, ZERO, ZERO /) Cun(5,:) = (/ ZERO, ZERO, ZERO, ZERO, ONE/G13, ZERO /) Cun(6,:) = (/ ZERO, ZERO, ZERO, ZERO, ZERO, ONE/G23 /) ! Hoffman model /extended Hill model/ C(1) = ONE/TWO*( ONE/(ft1*fc1)+ONE/(ft2*fc2)-ONE/(ft3*fc3)) C(2) = ONE/TWO*(-ONE/(ft1*fc1)+ONE/(ft2*fc2)+ONE/(ft3*fc3)) C(3) = ONE/TWO*( ONE/(ft1*fc1)-ONE/(ft2*fc2)+ONE/(ft3*fc3)) C(4) = ONE/fs12**TWO C(5) = ONE/fs13**TWO C(6) = ONE/fs23**TWO C(7) = (fc1 - ft1)/(fc1*ft1) C(8) = (fc2 - ft2)/(fc2*ft2) if (ft3==0) ft3=ft2 if (ft3==0) fc3 = ft3/(1-ft3*(C(7)+C(8))) C(9) = (fc3 - ft3)/(fc3*ft3) ! P-MATRIX Pun(1,:)=2.D0*(/ C(1)+C(3), -C(1), -C(3), ZERO, ZERO, ZERO /) Pun(2,:)=2.D0*(/ -C(1), C(2)+C(1), -C(2), ZERO, ZERO, ZERO /) Pun(3,:)=2.D0*(/ -C(3), -C(2), C(3)+C(2), ZERO, ZERO, ZERO /) Pun(4,:)=2.D0*(/ ZERO, ZERO, ZERO, C(4), ZERO, ZERO /) Pun(5,:)=2.D0*(/ ZERO, ZERO, ZERO, ZERO, C(5), ZERO /) Pun(6,:)=2.D0*(/ ZERO, ZERO, ZERO, ZERO, ZERO, C(6) /) ! Q-MATRIX Qun = (/ C(7), C(8), C(9), ZERO, ZERO, ZERO /) ! Equivalent plastic strain /ep_eff = sqrt(2/3*epl'*epl) eps_eff = STATEV(2*NTENS+1) if (eps_eff==ZERO) eps_eff = tol**two ! Reduction to PS or PE conditions ! ====================================================================== DO I=1,NTENS DO J=1,NTENS Cut(I,J)=ZERO P(I,J)=ZERO END DO Q(I)=ZERO END DO do i=1,NDI do j=1,NDI Cut(i,j)=Cun(i,j) P(i,j)=Pun(i,j) end do Q(i)=Qun(i) end do do i=1,nshr do j=1,nshr Cut(NDI+i,NDI+i) = Cun(3+i,3+i) P(ndi+i,ndi+i) = Pun(3+i,3+i) end do Q(NDI+i) = Qun(3+i) end do ! ====================================================================== ! material stiffness matrix /3D/ call matinv(ntens,Cut,D) ! Eye-matrix do i=1,ntens do j=1,ntens R3(i,j)=0.d0 if (i==j) R3(i,j)=1.d0 end do end do Z = R3 do j=1,nshr Z(ndi+j,ndi+j) = 5.d-1 end do ! TRIAL stress s_star = STRESS + matmul(D,DSTRAN) ! Yield function sy = s0 + rQ*(eps_eff)**rN ! sigma_y ds = rN*rQ*(eps_eff)**(rN-ONE) ! d.sigma_y/d.ep_eff R1 = matmul(P,s_star) seff = ONE/TWO*dot_product(s_star,R1) + sy*dot_product(Q,s_star) flyt = seff - sy**TWO ! gradient : d.flyt/d.sigma G = R1 + Q*sy ! d.ep_eff/d.lambda R2 = matmul(Z,G) de = sqrt(2.D0/3.D0*dot_product(G,R2)) ! flyt correction if (abs(flyt)100) then PNEWDT=0.25 return end if end do ! Specific plastic dissipation ! SPD = x1*de*sy ! Volumetric heat generation per unit time ! RPL = PROPS(21)*SPD/DTIME ! stress update STRESS = x2 ! =========================================================================== ! tangent material stiffnes matrix dFdE(1:NTENS,:) = D dXdE = matmul(dxdF,dFdE) DDSDDE = dXdE(2:NTENS+1,:) ! UPDATES STATEV(NTENS+1:2*NTENS) = STATEV(NTENS+1:2*NTENS) + x1*G(1:NTENS) STATEV(1:NTENS) = STATEV(1:NTENS) +DSTRAN-STATEV(NTENS+1:2*NTENS) STATEV(2*NTENS+1) = ep_eff endif return end subroutine matinv(n,a,ai) ! This subroutine inverts a matrix "a" and returns the inverse in "ai" ! n - Input by user, an integer specifying the size of the matrix to be inverted. ! a - Input by user, an n by n real array containing the matrix to be inverted. ! ai - Returned by subroutine, an n by n real array containing the inverted matrix. ! d - Work array, an n by 2n real array used by the subroutine. ! io - Work array, a 1-dimensional integer array of length n used by the subroutine. ! From http://www.mae.usu.edu/faculty/wphillips/MAE6510.html implicit real*8 (a-h, o-z) character*80 CMNAME dimension a(n,n), ai(n,n), io(n), d(n,2*n) ! Fill in the "io" and "d" matrix. do i=1,n io(i)=i end do do i=1,n do j=1,n d(i,j)=a(i,j) if(i.eq.j)then d(i,n+j)=1.0D0 else d(i,n+j)=0.0D0 endif end do end do ! Scaling do i=1,n m=1 do k=2,n if(abs(d(i,k)).gt.abs(d(i,m)))m=k end do tmp=d(i,m) do k=1,2*n d(i,k)=d(i,k)/tmp end do end do ! Lower Elimination do i=1,n-1 ! Pivoting m=i do j=i+1,n if(abs(d(io(j),i)).gt.abs(d(io(m),i)))m=j end do itmp=io(m) io(m)=io(i) io(i)=itmp ! Scale the Pivot element to unity r=d(io(i),i) do k=1,2*n d(io(i),k)=d(io(i),k)/r end do do j=i+1,n r=d(io(j),i) do k=1,2*n d(io(j),k)=d(io(j),k)-r*d(io(i),k) end do end do end do ! Upper Elimination r=d(io(n),n) do k=1,2*n d(io(n),k)=d(io(n),k)/r end do do i=n-1,1,-1 do j=i+1,n r=d(io(i),j) do k=1,2*n d(io(i),k)=d(io(i),k)-r*d(io(j),k) end do end do end do ! Fill Out "ai" matrix do i=1,n do j=1,n ai(i,j)=d(io(i),n+j) end do end do return end subroutine matinv