! ! CalculiX - A 3-dimensional finite element program ! Copyright (C) 1998-2015 Guido Dhondt ! ! This program is free software; you can redistribute it and/or ! modify it under the terms of the GNU General Public License as ! published by the Free Software Foundation(version 2); ! ! ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License ! along with this program; if not, write to the Free Software ! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. ! subroutine umat_xtal_wrapper(amat,iel,iint,kode,elconloc,emec, & emec0, & beta,xokl,voj,xkl,vj,ithermal,t1l,dtime,time,ttime, & icmd,ielas,mi,nstate_,xstateini,xstate,stre,stiff, & iorien,pgauss,orab,kstep,kinc,pnewdt) ! ! calculates stiffness and stresses for a nonlinear material ! defined by an ABAQUS umat routine ! ! icmd=3: calculates stress at mechanical strain ! else: calculates stress at mechanical strain and the stiffness ! matrix ! ! INPUT: ! ! amat material name ! iel element number ! iint integration point number ! ! kode material type (-100-#of constants entered ! under *USER MATERIAL): can be used for materials ! with varying number of constants ! ! elconloc(21) user defined constants defined by the keyword ! card *USER MATERIAL (max. 21, actual # = ! -kode-100), interpolated for the ! actual temperature t1l ! ! emec(6) Lagrange mechanical strain tensor (component order: ! 11,22,33,12,13,23) at the end of the increment ! (thermal strains are subtracted) ! emec0(6) Lagrange mechanical strain tensor at the start of the ! increment (thermal strains are subtracted) ! beta(6) residual stress tensor (the stress entered under ! the keyword *INITIAL CONDITIONS,TYPE=STRESS) ! ! xokl(3,3) deformation gradient at the start of the increment ! voj Jacobian at the start of the increment ! xkl(3,3) deformation gradient at the end of the increment ! vj Jacobian at the end of the increment ! ! ithermal 0: no thermal effects are taken into account ! >0: thermal effects are taken into account (triggered ! by the keyword *INITIAL CONDITIONS,TYPE=TEMPERATURE) ! t1l temperature at the end of the increment ! dtime time length of the increment ! time step time at the end of the current increment ! ttime total time at the start of the current step ! ! icmd not equal to 3: calculate stress and stiffness ! 3: calculate only stress ! ielas 0: no elastic iteration: irreversible effects ! are allowed ! 1: elastic iteration, i.e. no irreversible ! deformation allowed ! ! mi(1) max. # of integration points per element in the ! model ! nstate_ max. # of state variables in the model ! ! xstateini(nstate_,mi(1),# of elements) ! state variables at the start of the increment ! xstate(nstate_,mi(1),# of elements) ! state variables at the end of the increment ! ! stre(6) Piola-Kirchhoff stress of the second kind ! at the start of the increment ! ! iorien number of the local coordinate axis system ! in the integration point at stake (takes the value ! 0 if no local system applies) ! pgauss(3) global coordinates of the integration point ! orab(7,*) description of all local coordinate systems. ! If a local coordinate system applies the global ! tensors can be obtained by premultiplying the local ! tensors with skl(3,3). skl is determined by calling ! the subroutine transformatrix: ! call transformatrix(orab(1,iorien),pgauss,skl) ! ! ! OUTPUT: ! ! xstate(nstate_,mi(1),# of elements) ! updated state variables at the end of the increment ! stre(6) Piola-Kirchhoff stress of the second kind at the ! end of the increment ! stiff(21): consistent tangent stiffness matrix in the material ! frame of reference at the end of the increment. In ! other words: the derivative of the PK2 stress with ! respect to the Lagrangian strain tensor. The matrix ! is supposed to be symmetric, only the upper half is ! to be given in the same order as for a fully ! anisotropic elastic material (*ELASTIC,TYPE=ANISO). ! ! This routine allows for the use of an ABAQUS umat user subroutine ! in CalculiX. ! ! Note that the following fields are not supported ! so far: ! sse, - not used ! spd, - not used ! scd, - not used ! rpl, - not used ! ddsddt, - not used ! drplde, - not used ! drpldt, - not used ! predef, - not used ! dpred, - not used ! drot, - MUST CALCULATE ! pnewdt, - Not supported, no error, just no effect ! celent, - not used ! layer, - not used ! kspt - not used ! ! Furthermore, the following fields have a different meaning in ! ABAQUS and CalculiX: ! ! stran: in CalculiX: Lagrangian strain tensor (Converted before XTAL) ! in ABAQUS: logarithmic strain tensor ! dstran: in CalculiX: Lagrangian strain increment tensor (Converted) ! in ABAQUS: logarithmic strain increment tensor ! temp: in CalculiX: temperature at the end of the increment ! in ABAQUS: temperature at the start of the increment ! dtemp: in CalculiX: zero ! in ABAQUS: temperature increment ! ! Because of this, this routine should only be used for small ! deformations and small rotations (in that case all strain ! measures basically reduce to the infinitesimal strain). ! -- With the strain conversions, this can now be used for large ! strain problems as well. B.D. Huddleston (2017) ! implicit none ! character*80 amat ! integer ithermal,icmd,kode,ielas,iel,iint,nstate_,mi(*),i,j, & iorien, & ndi,nshr,ntens,nprops,layer,kspt,kstep,kinc,kal(2,6),kel(4,21), & j1,j2,j3,j4,j5,j6,j7,j8,jj ! real*8 elconloc(21),stiff(21),emec(6),emec0(6),beta(6), & stre(6), & vj,t1l,dtime,xkl(3,3),xokl(3,3),voj,pgauss(3),orab(7,*), & time,ttime,skl(3,3),xa(3,3),ya(3,3,3,3),xstate(nstate_,mi(1),*), & xstateini(nstate_,mi(1),*) ! real*8 ddsdde(6,6),sse,spd,scd,rpl,ddsddt(6),drplde(6), & drpldt,stran(6),dstran(6),abqtime(2),predef,temp,dtemp, & dpred,drot(3,3),celent,pnewdt, notdstran(6) ! real*8 utemp(3,3),vtemp(3,3),wtemp(3),ztemp(3,3),btemp(3),svderr, & rvtemp(3),roto(3,3),rot(3,3),fi(3,3),fit(3,3),xb(3,3) ! kal=reshape((/1,1,2,2,3,3,1,2,1,3,2,3/),(/2,6/)) ! kel=reshape((/1,1,1,1,1,1,2,2,2,2,2,2,1,1,3,3,2,2,3,3,3,3,3,3, & 1,1,1,2,2,2,1,2,3,3,1,2,1,2,1,2,1,1,1,3,2,2,1,3, & 3,3,1,3,1,2,1,3,1,3,1,3,1,1,2,3,2,2,2,3,3,3,2,3, & 1,2,2,3,1,3,2,3,2,3,2,3/),(/4,21/)) ! drot=reshape((/1.d0,0.d0,0.d0,0.d0,1.d0,0.d0,0.d0,0.d0,1.d0/), & (/3,3/)) ! ! calculating the mechanical strain ! do i=1,6 stran(i)=emec0(i) dstran(i)=emec(i)-emec0(i) enddo ! ntens=6 ! do i=1,nstate_ xstate(i,iint,iel)=xstateini(i,iint,iel) enddo ! abqtime(1)=time-dtime abqtime(2)=ttime+time-dtime ! temp=t1l dtemp=0.d0 ! ndi=3 nshr=3 ntens=ndi+nshr ! nprops=-kode-100 c nprops=21 ! ! taking local material orientations into account ! if(iorien.ne.0) then call transformatrix(orab(1,iorien),pgauss,skl) ! ! rotating the stress into the local system ! xa(1,1)=stre(1) xa(1,2)=stre(4) xa(1,3)=stre(5) xa(2,1)=stre(4) xa(2,2)=stre(2) xa(2,3)=stre(6) xa(3,1)=stre(5) xa(3,2)=stre(6) xa(3,3)=stre(3) ! do jj=1,6 stre(jj)=0.d0 j1=kal(1,jj) j2=kal(2,jj) do j3=1,3 do j4=1,3 stre(jj)=stre(jj)+ & xa(j3,j4)*skl(j3,j1)*skl(j4,j2) enddo enddo enddo ! ! rotating the strain into the local system ! xa(1,1)=stran(1) xa(1,2)=stran(4) xa(1,3)=stran(5) xa(2,1)=stran(4) xa(2,2)=stran(2) xa(2,3)=stran(6) xa(3,1)=stran(5) xa(3,2)=stran(6) xa(3,3)=stran(3) ! do jj=1,6 stran(jj)=0.d0 j1=kal(1,jj) j2=kal(2,jj) do j3=1,3 do j4=1,3 stran(jj)=stran(jj)+ & xa(j3,j4)*skl(j3,j1)*skl(j4,j2) enddo enddo enddo ! ! rotating the strain increment into the local system ! xa(1,1)=dstran(1) xa(1,2)=dstran(4) xa(1,3)=dstran(5) xa(2,1)=dstran(4) xa(2,2)=dstran(2) xa(2,3)=dstran(6) xa(3,1)=dstran(5) xa(3,2)=dstran(6) xa(3,3)=dstran(3) ! do jj=1,6 dstran(jj)=0.d0 j1=kal(1,jj) j2=kal(2,jj) do j3=1,3 do j4=1,3 dstran(jj)=dstran(jj)+ & xa(j3,j4)*skl(j3,j1)*skl(j4,j2) enddo enddo enddo endif ! ! changing physical strain into engineering strain ! do i=4,6 stran(i)=2.d0*stran(i) dstran(i)=2.d0*dstran(i) enddo ! ! Change strain from lagrangian to logarithmic ! call pushforward(xkl, stran, 3) call pushforward(xkl, dstran, 3) ! ! Calculate drot ! call get_drot(notdstran, drot, xokl, xkl, dtime, ntens) ! ! Call the CPFEM routine ! call umat_xtal(stre, xstate(1,iint,iel), ddsdde, sse, spd, & scd, rpl, ddsddt, drplde, drpldt, & stran, dstran, abqtime, dtime, temp, & dtemp, predef, dpred, amat, ndi, & nshr, ntens, nstate_, elconloc, nprops, & pgauss, drot, pnewdt, celent, xokl, & xkl, iel, iint, layer, kspt, & kstep, kinc) ! ! pull back cauchy stress to 2nd PK stress ! xa(1,1) = stre(1) xa(2,2) = stre(2) xa(3,3) = stre(3) xa(1,2) = stre(4) xa(2,1) = stre(4) xa(1,3) = stre(5) xa(3,1) = stre(5) xa(2,3) = stre(6) xa(3,2) = stre(6) do i=1,3 do j=1,3 fi(i,j) = xkl(i,j) enddo enddo call invert_matrix(fi,3) do i=1,3 do j=1,3 fit(i,j) = fi(j,i) enddo enddo call matrix_product(fi, xa, 3, xb) call matrix_product(xb, fit, 3, xa) stre(1) = xa(1,1)*vj stre(2) = xa(2,2)*vj stre(3) = xa(3,3)*vj stre(4) = xa(1,2)*vj stre(5) = xa(1,3)*vj stre(6) = xa(2,3)*vj ! ! taking local material orientations into account ! if(iorien.ne.0) then ! ! rotating the stress into the global system ! xa(1,1)=stre(1) xa(1,2)=stre(4) xa(1,3)=stre(5) xa(2,1)=stre(4) xa(2,2)=stre(2) xa(2,3)=stre(6) xa(3,1)=stre(5) xa(3,2)=stre(6) xa(3,3)=stre(3) ! do jj=1,6 stre(jj)=0.d0 j1=kal(1,jj) j2=kal(2,jj) do j3=1,3 do j4=1,3 stre(jj)=stre(jj)+ & xa(j3,j4)*skl(j1,j3)*skl(j2,j4) enddo enddo enddo endif ! ! calculate the stiffness matrix (the matrix is symmetrized) ! if(icmd.ne.3) then stiff(1)=ddsdde(1,1) stiff(2)=(ddsdde(1,2)+ddsdde(2,1))/2.d0 stiff(3)=ddsdde(2,2) stiff(4)=(ddsdde(1,3)+ddsdde(3,1))/2.d0 stiff(5)=(ddsdde(2,3)+ddsdde(3,2))/2.d0 stiff(6)=ddsdde(3,3) stiff(7)=(ddsdde(1,4)+ddsdde(4,1))/2.d0 stiff(8)=(ddsdde(2,4)+ddsdde(4,2))/2.d0 stiff(9)=(ddsdde(3,4)+ddsdde(4,3))/2.d0 stiff(10)=ddsdde(4,4) stiff(11)=(ddsdde(1,5)+ddsdde(5,1))/2.d0 stiff(12)=(ddsdde(2,5)+ddsdde(5,2))/2.d0 stiff(13)=(ddsdde(3,5)+ddsdde(5,3))/2.d0 stiff(14)=(ddsdde(4,5)+ddsdde(5,4))/2.d0 stiff(15)=ddsdde(5,5) stiff(16)=(ddsdde(1,6)+ddsdde(6,1))/2.d0 stiff(17)=(ddsdde(2,6)+ddsdde(6,2))/2.d0 stiff(18)=(ddsdde(3,6)+ddsdde(6,3))/2.d0 stiff(19)=(ddsdde(4,6)+ddsdde(6,4))/2.d0 stiff(20)=(ddsdde(5,6)+ddsdde(6,5))/2.d0 stiff(21)=ddsdde(6,6) ! if(iorien.ne.0) then ! ! rotating the stiffness coefficients into the global system ! call anisotropic(stiff,ya) ! do jj=1,21 j1=kel(1,jj) j2=kel(2,jj) j3=kel(3,jj) j4=kel(4,jj) stiff(jj)=0.d0 do j5=1,3 do j6=1,3 do j7=1,3 do j8=1,3 stiff(jj)=stiff(jj)+ya(j5,j6,j7,j8)* & skl(j1,j5)*skl(j2,j6)*skl(j3,j7)*skl(j4,j8) enddo enddo enddo enddo enddo endif endif return end ! SUBROUTINE GET_DROT(DSTRAN,DR,DFGRD0,DFGRD1,DTIME,NTENS) C======================================================================+ C----------------------------------------------------------------------+ C---------------- Calculation of rotation increment -----------------+ C-------------------- for the jaumann derivative --------------------+ C---------------------- method of hughes-winget ---------------------+ C-------------------- Courtesy of Dr. Youssef Hammi ------------------+ C----------------------------------------------------------------------+ C======================================================================+ C---------- C INPUT : C----------- C DFGRD0 : DEFORMATION GRADIENT AT TIME T C DFGRD1 : DEFORMATION GRADIENT AT TIME T+DT C DTIME : TIME INCREMENT C NTENS : NUMBER OF TENSOR COMPONENTS C----------- C OUTPUT : C----------- C DR : ROTATION INCEMENT MATRIX C DSTRAN : STRAIN INCREMENT TENSOR C----------------------------------------------------------------------+ C======================================================================= C----------------------------------------------------------------------+ C.1----- Precision implicit double precision (a-h,o-z) C.2----- Parameter parameter (n=3) C.3----- Dimension dimension dstran(ntens) dimension dfgrd0(3,3),dfgrd1(3,3),dr(3,3),df(3,3),f0i(3,3) dimension dqdq(3,3),dfm(3,3),dfp(3,3) dimension vl(3,3),dw(3,3),dm(3,3),dp(3,3) dimension tdr(3,3),ti(3,3) C.4----- Real,Integer,Complex,Double precision,Logical,Character C.5----- Common C.6----- Equivalence C.7----- Save,Data data half,one,two/0.5d0,1.d0,2.d0/ C.8----- Functions Definition C----------------------------------------------------------------------+ C======================================================================+ C----------------------------------------------------------------------+ C------- Deformation Gradient Increment DF C DF = F1 * F0^-1 do i=1,N do j=1,N F0I(i,j) = dfgrd0(i,j) end do end do call invert_matrix(F0I,N) call matrix_product(dfgrd1, F0I, N, dF) C----------------------------------------------------------------------+ C----------------------------------------------------------------------+ C------- Velocity Gradient L(t+Dt/2) C do i=1,N do j=1,N dfm(i,j) = df(i,j) dfp(i,j) = df(i,j) end do dfm(i,i) = dfm(i,i) - one dfp(i,i) = dfp(i,i) + one end do call invert_matrix(dfp,N) call matrix_product(dfm, dfp, N, vl) twodt = two/dtime do i=1,N do j=1,N vl(i,j) = twodt*vl(i,j) end do end do C----------------------------------------------------------------------+ C------- Strain Increment Tensor Dstran and Plastic Spin DW C dstran(1) = vl(1,1) dstran(2) = vl(2,2) dstran(3) = vl(3,3) dstran(4) = vl(1,2) + vl(2,1) if(ntens.gt.4) then dstran(5) = vl(1,3) + vl(3,1) dstran(6) = vl(2,3) + vl(3,2) end if do i=1,N do j=1,N dw(i,j) = half*(vl(i,j)-vl(j,i)) end do end do C----------------------------------------------------------------------+ C------- Increment Rotation Matrix DR (Jaumann) C halfdt = half*dtime do i=1,N do j=1,N dp(i,j) = halfdt*dw(i,j) dm(i,j) = -dp(i,j) end do dp(i,i) = dp(i,i) + one dm(i,i) = dm(i,i) + one end do call invert_matrix(dm,N) call matrix_product(dm, dp, N, dr) C----------------------------------------------------------------------+ C======================================================================+ C----------------------------------------------------------------------+ C--------- formats -------------------------------------------------+ C----------------------------------------------------------------------+ 101 format(1x,3(d18.12,2x)) c----------------------------------------------------------------------+ c======================================================================+ return end C subroutine invert_matrix(matrix, size) C implicit none C integer size, twosize, i, j, ii, jj C real*8 matrix(size, size) real*8 augmented(size, size*2) real*8 scale C twosize = size*2 C----- Augment the original matrix ------------------------------------+ do i=1,size do j=1,twosize if (j.le.size) then augmented(i,j) = matrix(i,j) elseif ((i+size).eq.j) then augmented(i,j) = 1 else augmented(i,j) = 0 end if end do end do C----- Forward elimination --------------------------------------------+ do i=1,size scale = augmented(i,i) do ii=1,twosize augmented(i,ii) = augmented(i,ii) / scale end do do j=i+1,size scale = augmented(j,i) do ii=1,twosize augmented(j,ii) = augmented(j,ii) - augmented(i,ii)*scale end do end do end do C----- Back propagation -----------------------------------------------+ do i=1,size-1 jj = size - i do j=jj+1,size scale = augmented(jj,j) do ii=1,twosize augmented(jj,ii) = augmented(jj,ii) - augmented(j,ii)*scale end do end do end do C----- Transfer inverted matrix back to matrix ------------------------+ do i=1,size do j=1,size matrix(i,j) = augmented(i,j+size) end do end do C return end C subroutine matrix_product(A, B, N, C) implicit none integer N, i, j, k real*8 A(N,N) real*8 B(N,N) real*8 C(N,N) do i=1,N do j=1,N C(i,j) = 0.0d0 do k=1,N C(i,j) = C(i,j) + A(i,k)*B(k,j) end do end do end do return end ! subroutine pushforward(F, x, n) implicit none integer i, j, k, n real*8 F(n,n), x(6) real*8 Fi(n,n), FiT(n,n), xfor(n,n), x1(n,n) do i=1,n do j=1,n Fi(i,j) = F(i,j) enddo enddo call invert_matrix(Fi, n) do i=1,n do j=1,n FiT(i,j) = Fi(j,i) enddo enddo xfor(1,1) = x(1) xfor(2,2) = x(2) xfor(3,3) = x(3) xfor(1,2) = x(4) xfor(2,1) = x(4) xfor(1,3) = x(5) xfor(3,1) = x(5) xfor(2,3) = x(6) xfor(3,2) = x(6) call matrix_product(FiT, xfor, n, x1) call matrix_product(x1, Fi, n, xfor) x(1) = xfor(1,1) x(2) = xfor(2,2) x(3) = xfor(3,3) x(4) = xfor(1,2) x(5) = xfor(1,3) x(6) = xfor(2,3) return end