Skip to content

Latest commit

 

History

History
281 lines (266 loc) · 12.8 KB

example_neohooke-LSDYNA.md

File metadata and controls

281 lines (266 loc) · 12.8 KB
Raw

Example: Neo-Hookean Material (in LS-DYNA)

The following example discusses the implementation of the standard Neo-Hookean material model to show the usage of the tensor toolbox in LS-DYNA. The material model describes hyperelasticity and is formulated in the spatial (Eulerian) configuration for 3D and axisymmetric computations (@todo not yet adapted for 2D). Be aware that in LS-Dyna umat (stress and history update) and utan (tangent) are split up, as noted in the section "User-defined tangent routine "utan"".

Kinematics

In LS-DYNA we first have to set the option "IHYPER=1" in the material card keyword (e.g. *MAT_USER_DEFINED_MATERIAL_MODELS) to be able to access the deformation gradient (see section "Example material card"). The latter is then stored on top of the defined number of history variables. So, in case we don't require any history variables for our material model (NHV=0), as for this hyperelastic model, the deformation gradient is stored in the initial "hsv" entries 1 to 9 columnwise. As a consquence, we can extract it via the function

type(Tensor2) :: defoGrad_F
defoGrad_F  = defoGrad(hsv(1:9))

which takes the list entries 1 to 9 and stores them into the second order tensor defoGrad_F.

Equations

The Cauchy stress tensor for the Neo-Hookean material can be computed from

The Eulerian tangent modulus is computed as

With

  • Jacobian, determinant of the deformation gradient

  • First and second Lame parameters and
  • Second order unit/identity tensor
  • Fourth order unit/identity tensor

User-defined material routine "umat"

The commented program

In LS-Dyna we can expand the existing file dyn21umats.F (under Windows). At the beginning of the Fortran file include the tensor toolbox.

      include  'ttb/ttb_library.f'

Then scroll down to an unused user-defined material, for instance umat43, reading

      subroutine umat43 (cm,eps,sig,epsp,hsv,dt1,capa,etype,tt,
     1 temper,failel,crv,nnpcrv,cma,qmat,elsiz,idele,reject)

Firstly, we have to declare the usage of the tensor module and the LS-DYNA extension.

c Use the tensor toolbox
      use Tensor
      use TensorXLSDYNA

Next, the standard LS-Dyna declarations follow and don't need to be altered.

c Standard LS-Dyna declarations
      include 'nlqparm'
      include 'bk06.inc'
      include 'iounits.inc'
      real(kind=8), dimension(*) :: cm, eps, sig, hsv
      dimension crv(lq1,2,*),cma(*),qmat(3,3)
      logical failel,reject
      integer nnpcrv(*)
      character*5 etype
      INTEGER8 idele

Now we declare our variables.

c declaration
      ! Deformation gradient (unsymmetric second order tensor)
       type(Tensor2) :: defoGrad_F
      ! Jacobian, determinant of the deformation gradient
       double precision :: det_F
      ! Cauchy stress tensor; unit tensor
       type(Tensor2) :: cauchyStress_sig, Eye
      ! material parameters
      double precision YoungsMod_E, Poisson_nu, lame_lambda, shearMod_mu

Extract the material parameters from the 'cm'-array

      YoungsMod_E = cm(1)
      Poisson_nu =  cm(2)

Compute the Lame parameters lambda and mu (or G)

      lame_lambda = YoungsMod_E * Poisson_nu
     &              / ( (1.+Poisson_nu)*(1.-2.*Poisson_nu) )
      shearMod_mu = .5 * YoungsMod_E / (1.+Poisson_nu)

Get the unit tensor via the tensor toolbox

      Eye = identity2(Eye)

Extract the deformation gradient from the history 'hsv' and transform it into the unsymmetric second order tensor as described above via the function 'defoGrad(*)' from the module TensorXLSDYNA

      defoGrad_F = defoGrad( hsv(1:9) )

Compute the Jacobian as the determinant of the deformation gradient

      det_F = det(defoGrad_F)

Compute the Cauchy stress for the Neo-Hookean material

      cauchyStress_sig = 1./det_F * (
     &      shearMod_mu * ( (defoGrad_F * transpose(defoGrad_F)) - Eye )
     &      + lame_lambda * log(det_F) * Eye )

Transform the stress tensor into the 'sig' array

      sig(1:6) = asarray(voigt(cauchyStress_sig),6)
c
      return
      end

The plain program

      subroutine umat43 (cm,eps,sig,epsp,hsv,dt1,capa,etype,tt,
     1 temper,failel,crv,nnpcrv,cma,qmat,elsiz,idele,reject)
c
      use Tensor
      use TensorXLSDYNA
c
      include 'nlqparm'
      include 'bk06.inc'
      include 'iounits.inc'
      real(kind=8), dimension(*) :: cm, eps, sig, hsv
      dimension crv(lq1,2,*),cma(*),qmat(3,3)
      logical failel,reject
      integer nnpcrv(*)
      character*5 etype
      INTEGER8 idele
c
      type(Tensor2) :: defoGrad_F
      double precision :: det_F
      type(Tensor2) :: cauchyStress_sig, Eye
      double precision YoungsMod_E, Poisson_nu, lame_lambda, shearMod_mu
c
      YoungsMod_E = cm(1)
      Poisson_nu =  cm(2)
c
      lame_lambda = YoungsMod_E * Poisson_nu
     &              / ( (1.+Poisson_nu)*(1.-2.*Poisson_nu) )
      shearMod_mu = .5 * YoungsMod_E / (1.+Poisson_nu)
c
      Eye = identity2(Eye)
c
      defoGrad_F = defoGrad( hsv(1:9) )
c    
      det_F = det(defoGrad_F)
c
      cauchyStress_sig = 1./det_F * (
     &      shearMod_mu * ( (defoGrad_F * transpose(defoGrad_F)) - Eye )
     &      + lame_lambda * log(det_F) * Eye )
c
      sig(1:6) = asarray(voigt(cauchyStress_sig),6)
c
      return
      end

User-defined tangent routine "utan"

The commented program

In LS-Dyna the computation of the material (stress and history update) and the tangent are split up into two separate subroutines. The tangent is implemented in the file dyn21utan.F into the subroutines utanXX, where the material id XX has to correspond to the id of the umatXX subroutine it belongs to. We have already included the necessary files in the dyn21umats.F file, so for the tangent we can directly scroll down to the subroutine utan43.

      subroutine utan43(cm,eps,sig,epsp,hsv,dt1,unsym,capa,etype,tt,
     1 temper,es,crv,nnpcrv,failel,cma,qmat)

We again load the modules

c Use the tensor toolbox
      use Tensor
      use TensorXLSDYNA

Use the standard LS-Dyna declarations

c Standard LS-Dyna declarations
      include 'nlqparm'
      real(kind=8), dimension (*) :: cm, eps, sig, hsv
      dimension crv(lq1,2,*),cma(*)
      integer nnpcrv(*)
      dimension es(6,*),qmat(3,3)
      logical failel,unsym
      character*5 etype

Declare our variables now also containing the fourth order tensor tangent_E for our tangent modulus.

      ! Deformation gradient (unsymmetric second order tensor)
       type(Tensor2) :: defoGrad_F
      ! Jacobian, determinant of the deformation gradient
       double precision :: det_F
      ! unit tensor
       type(Tensor2) :: Eye
      ! Fourth order Eulerian tangent modulus
       type(Tensor4) :: tangent_E
      ! material parameters
      double precision YoungsMod_E, Poisson_nu, lame_lambda, shearMod_mu

As in the umat subroutine:

c Extract the material parameters from the 'cm'-array
      YoungsMod_E = cm(1)
      Poisson_nu = cm(2)
c Compute the Lame parameters lambda and mu (or G)
      lame_lambda = YoungsMod_E * Poisson_nu
     &              / ((1.+Poisson_nu)*(1.-2.*Poisson_nu))
      shearMod_mu = .5*YoungsMod_E / (1.+Poisson_nu)
c Get the unit tensor via the tensor toolbox
      Eye = identity2(Eye)
c Extract the deformation gradient from the history 'hsv'
      defoGrad_F = defoGrad( hsv(1:9) )
c Compute the Jacobian as the determinant of the deformation gradient      
      det_F = det(defoGrad_F)

In order to finally compute the Eulerian tangent for the Neo-Hookean material

      tangent_E = 1./det_F * (
     &       lame_lambda * (Eye.dya.Eye) ! Be aware of the required parentheses,
                                         ! Else error "An arithmetic or LOGICAL type is required in this context."
     &       + ( 2. * ( shearMod_mu - lame_lambda * log(det_F) )
     &              * identity4(Eye) )
     &       )

Eventually, we pack the fourth order tensor into the 6x6 matrix es

      es(1:6,1:6) = asarray(voigt(tangent_E),6,6)
c      
      return
      end

and end the subroutine.

The plain program

      subroutine utan43(cm,eps,sig,epsp,hsv,dt1,unsym,capa,etype,tt,
     1 temper,es,crv,nnpcrv,failel,cma,qmat)
c
      use Tensor
      use TensorXLSDYNA
c
      include 'nlqparm'
      real(kind=8), dimension (*) :: cm, eps, sig, hsv
      dimension crv(lq1,2,*),cma(*)
      integer nnpcrv(*)
      dimension es(6,*),qmat(3,3)
      logical failel,unsym
      character*5 etype
c
       type(Tensor2) :: defoGrad_F
       double precision :: det_F
       type(Tensor2) :: Eye
       type(Tensor4) :: tangent_E
      double precision YoungsMod_E, Poisson_nu, lame_lambda, shearMod_mu
c
      YoungsMod_E = cm(1)
      Poisson_nu = cm(2)
c
      lame_lambda = YoungsMod_E * Poisson_nu
     &              / ((1.+Poisson_nu)*(1.-2.*Poisson_nu))
      shearMod_mu = .5*YoungsMod_E / (1.+Poisson_nu)
c
      Eye = identity2(Eye)
c
      defoGrad_F = defoGrad( hsv(1:9) )
c     
      det_F = det(defoGrad_F)
c
      tangent_E = 1./det_F * (
     &       lame_lambda * (Eye.dya.Eye)
     &       + ( 2. * ( shearMod_mu - lame_lambda * log(det_F) )
     &              * identity4(Eye) )
     &       )
c
      es(1:6,1:6) = asarray(voigt(tangent_E),6,6)
c      
      return
      end

Example material card

The material id in the option "MT" must equal the id in umatXX and utanXX. The value of "NHV" sets the number of used history variables, here none are used in the material model. The option "IHYPER=1" stores the deformation gradient in the history "hsv" on top of the defined "NHV". The parameters "P1" and "P2" contain the Young's modulus in "cm(1)" and the Poisson ratio in "cm(2)", respectively.