umat_ciarlet_el.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	umat_ciarlet_el (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)
subroutine	voigt_33mat_to_6vec (mat, vec)
subroutine	voigt_tetrad_to_matrix (C, stiff)
subroutine	s2_ciarlet (Cb, lamda, mu, S2, Cbi, detC)
subroutine	ds2_ciarlet_de (lamda, mu, detC, Cbi, A)
subroutine	m33inv_det (A, AINV, DET, OK_FLAG)

Function/Subroutine Documentation

ds2_ciarlet_de()

subroutine ds2_ciarlet_de	(	double precision	lamda,
		double precision	mu,
		double precision	detC,
		double precision, dimension(3,3)	Cbi,
		double precision, dimension(3,3,3,3)	A
	)

C     computes linearization of 2PK stresses with E
C     for Ciarlet-model
      implicit none
C     input:
      double precision lamda,mu,detc,cbi(3,3)
C     output:
      double precision a(3,3,3,3)
C     local:
      double precision f1,f2
      integer i,j,k,l
      
      f1 = lamda*detc
      f2 = 2.d0*mu - lamda*(detc-1.d0)
      do i=1,3
         do j=1,3
            do k=1,3
               do l=1,3
                  a(i,j,k,l) =  
     &                 f1 * cbi(k,l)*cbi(i,j) +
     &                 f2 * cbi(i,k)*cbi(l,j)     
               enddo
            enddo
         enddo
      enddo
      return

m33inv_det()

subroutine m33inv_det	(	double precision, dimension(3,3), intent(in)	A,
		double precision, dimension(3,3), intent(out)	AINV,
		double precision	DET,
		logical, intent(out)	OK_FLAG
	)

C     
C
C     !******************************************************************
C     !  M33INV_DET  -  Compute the inverse of a 3x3 matrix.
C     !
C     !  A       = input 3x3 matrix to be inverted
C     !  AINV    = output 3x3 inverse of matrix A
C     !  OK_FLAG = (output) .TRUE. if the input matrix could be inverted, 
C     !            and .FALSE. if the input matrix is singular.
C     !******************************************************************

      IMPLICIT NONE

      DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN)  :: a
      DOUBLE PRECISION, DIMENSION(3,3), INTENT(OUT) :: ainv
      LOGICAL, INTENT(OUT) :: ok_flag

      DOUBLE PRECISION, PARAMETER :: eps = 1.0d-10
      DOUBLE PRECISION :: det
      DOUBLE PRECISION, DIMENSION(3,3) :: cofactor


      det =   a(1,1)*a(2,2)*a(3,3) - a(1,1)*a(2,3)*a(3,2)  
     &      - a(1,2)*a(2,1)*a(3,3) + a(1,2)*a(2,3)*a(3,1) 
     &      + a(1,3)*a(2,1)*a(3,2) - a(1,3)*a(2,2)*a(3,1)

      IF (abs(det) .LE. eps) THEN
         ainv = 0.0d0
         ok_flag = .false.
         RETURN
      END IF

      cofactor(1,1) = +(a(2,2)*a(3,3)-a(2,3)*a(3,2))
      cofactor(1,2) = -(a(2,1)*a(3,3)-a(2,3)*a(3,1))
      cofactor(1,3) = +(a(2,1)*a(3,2)-a(2,2)*a(3,1))
      cofactor(2,1) = -(a(1,2)*a(3,3)-a(1,3)*a(3,2))
      cofactor(2,2) = +(a(1,1)*a(3,3)-a(1,3)*a(3,1))
      cofactor(2,3) = -(a(1,1)*a(3,2)-a(1,2)*a(3,1))
      cofactor(3,1) = +(a(1,2)*a(2,3)-a(1,3)*a(2,2))
      cofactor(3,2) = -(a(1,1)*a(2,3)-a(1,3)*a(2,1))
      cofactor(3,3) = +(a(1,1)*a(2,2)-a(1,2)*a(2,1))

      ainv = transpose(cofactor) / det

      ok_flag = .true.

      RETURN

s2_ciarlet()

subroutine s2_ciarlet	(	double precision, dimension(3,3)	Cb,
		double precision	lamda,
		double precision	mu,
		double precision, dimension(3,3)	S2,
		double precision, dimension(3,3)	Cbi,
		double precision	detC
	)

C     computes 2PK stresses for Ciarlet-model
C     
      implicit none
C     input:
      double precision cb(3,3),lamda,mu
C     output:
      double precision s2(3,3), cbi(3,3), detc
C     local:
      double precision id(3,3)
      double precision f1
      integer i,j
      logical ok_flag

C     Set id:
      do i=1,3
         do j=1,3
            id(i,j)=0.d0
         enddo
         id(i,i)=1.d0
      enddo

C     Inverse of Cb and det(C):
      call m33inv_det(cb, cbi, detc, ok_flag)
C
C
      f1 = 0.5d0*lamda*(detc-1.d0)
      do i=1,3
         do j=1,3
            s2(i,j)=f1*cbi(i,j) + mu * ( id(i,j)-cbi(i,j) )
         enddo
      enddo
      return

umat_ciarlet_el()

subroutine umat_ciarlet_el	(	character*80	amat,
		integer	iel,
		integer	iint,
		integer	kode,
		real*8, dimension(21)	elconloc,
		real*8, dimension(6)	emec,
		real*8, dimension(6)	emec0,
		real*8, dimension(6)	beta,
		real*8, dimension(3,3)	xokl,
		real*8	voj,
		real*8, dimension(3,3)	xkl,
		real*8	vj,
		integer	ithermal,
		real*8	t1l,
		real*8	dtime,
		real*8	time,
		real*8	ttime,
		integer	icmd,
		integer	ielas,
		integer, dimension(*)	mi,
		integer	nstate_,
		real*8, dimension(nstate_,mi(1),*)	xstateini,
		real*8, dimension(nstate_,mi(1),*)	xstate,
		real*8, dimension(6)	stre,
		real*8, dimension(21)	stiff,
		integer	iorien,
		real*8, dimension(3)	pgauss,
		real*8, dimension(7,*)	orab
	)

!
!     calculates stiffness and stresses for a user defined material
!     law
!
!     icmd=3: calcutates 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 increment
!
!     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).
!                        Notice that the matrix is an integral part of the 
!                        fourth order material tensor, i.e. the Voigt notation
!                        is not used.
!
      implicit none
!
      character*80 amat
      integer ithermal,icmd,kode,ielas,iel,iint,nstate_,mi(*),iorien
      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
      real*8 xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*)
!
!     
!     Ela-Pla-modifications:
      real*8 c1(3,3), c1i(3,3), detc1, ds2de(3,3,3,3), s2pk1(3,3)
      real*8 lamda, mu

      lamda=elconloc(1)
      mu=elconloc(2)

C     ciarlet_doc
C
C     Compute pullback C1(3,3) of metric tensor g.
      c1  = matmul(transpose(xkl),xkl)
C
C     2PK-Stress at C:
      call s2_ciarlet(c1,lamda,mu, s2pk1,c1i,detc1)
C     Tangent/linearization:
      call ds2_ciarlet_de(lamda,mu,detc1,c1i,ds2de)

C     Dhondt-Voigt-Order, see umat_lin_iso_el.f
      call voigt_33mat_to_6vec(s2pk1,stre)
      call voigt_tetrad_to_matrix(ds2de,stiff)

C     ciarlet_cod

      return

voigt_33mat_to_6vec()

subroutine voigt_33mat_to_6vec	(	double precision, dimension(3,3)	mat,
		double precision, dimension(6)	vec
	)

      implicit none
C     in:
      double precision mat(3,3)
C     out:
      double precision vec(6)
C     11,22,33,12,13,23
      vec(1)=mat(1,1)
      vec(2)=mat(2,2)
      vec(3)=mat(3,3)
      vec(4)=mat(1,2)
      vec(5)=mat(1,3)
      vec(6)=mat(2,3)

      return

voigt_tetrad_to_matrix()

subroutine voigt_tetrad_to_matrix	(	double precision, dimension(3,3,3,3)	C,
		double precision, dimension(21)	stiff
	)

      implicit none
C     in:
      double precision c(3,3,3,3)
C     out:
      double precision stiff(21)
      
C     indices according to 
C     1) file anisotropic.f and 
C     2) html-documentation
C     -> Making C symmetric in the **last two** indices:
C     
C     stiff stores:
C     1111, 1122, 2222, 1133
C     2233, 3333, 1112, 2212
C
C     1111
      stiff(1)  = c(1,1,1,1)
C     1122 = 2211
      stiff(2)  = c(1,1,2,2)
C     2222
      stiff(3)  = c(2,2,2,2)
C     1133 = 3311
      stiff(4)  = c(1,1,3,3)
C     2233 = 3322
      stiff(5)  = c(2,2,3,3)
C     3333
      stiff(6)  = c(3,3,3,3)
C     1112 = 1121 = 1211 = 2111
      stiff(7)  = 0.5d0*(c(1,1,1,2)+c(1,1,2,1))
C     2212 = 1222 = 2122 = 2221
      stiff(8)  = 0.5d0*(c(2,2,1,2)+c(2,2,2,1))
C
C     stiff stores:
C     3312, 1212, 1113, 2213
C     3313, 1213, 1313, 1123
C     
C     3312 = 1233 = 2133 = 3321
      stiff(9)  = 0.5d0*(c(3,3,1,2)+c(3,3,2,1))
C     1212 = 1221 = 2112 = 2121
      stiff(10) = 0.5d0*(c(1,2,1,2)+c(1,2,2,1))
C     1113 = 1131 = 1311 = 3111
      stiff(11) = 0.5d0*(c(1,1,1,3)+c(1,1,3,1))
C     2213 = 1322 = 2231 = 3122
      stiff(12) = 0.5d0*(c(2,2,1,3)+c(2,2,3,1))
C     
C     3313 = 1333 = 3133 = 3331
      stiff(13) = 0.5d0*(c(3,3,1,3)+c(3,3,3,1))     
C     1213 = 1231 = 1312 = 1321 = 2113 = 2131 = 3112 = 3121:
      stiff(14) = 0.5d0*(c(1,2,1,3)+c(1,2,3,1))
C     1313 = 1331 = 3113 = 3131
      stiff(15) = 0.5d0*(c(1,3,1,3)+c(1,3,3,1))
C     1123 = 1132 = 2311 = 3211
      stiff(16) = 0.5d0*(c(1,1,2,3)+c(1,1,3,2))
C
C     stiff stores:
C     2223, 3323, 1223, 1323
C     2323
C
C     2223 = 2232 = 2322 = 3222
      stiff(17) = 0.5d0*(c(2,2,2,3)+c(2,2,3,2)) 
C     3323 = 2333 = 3233 = 3332
      stiff(18) = 0.5d0*(c(3,3,2,3)+c(3,3,3,2))
C     1223 = 1232 = 2123 = 2132 = 2312 = 2321 = 3212 = 3221
      stiff(19) = 0.5d0*(c(1,2,2,3)+c(1,2,3,2))
C     1323 = 1332 = 2313 = 2331 = 3123 = 3132 = 3213 = 3231
      stiff(20) = 0.5d0*(c(1,3,2,3)+c(1,3,3,2))
C     2323 = 2332
      stiff(21) = 0.5d0*(c(2,3,2,3)+c(2,3,3,2))

      return