v3d4_NeoHookeanInCompressDef.f90

Go to the documentation of this file.

!*********************************************************************
!* Illinois Open Source License                                      *
!*                                                                   *
!* University of Illinois/NCSA                                       * 
!* Open Source License                                               *
!*                                                                   *
!* Copyright@2008, University of Illinois.  All rights reserved.     *
!*                                                                   *
!*  Developed by:                                                    *
!*                                                                   *
!*     Center for Simulation of Advanced Rockets                     *
!*                                                                   *
!*     University of Illinois                                        *
!*                                                                   *
!*     www.csar.uiuc.edu                                             *
!*                                                                   *
!* Permission is hereby granted, free of charge, to any person       *
!* obtaining a copy of this software and associated documentation    *
!* files (the "Software"), to deal with the Software without         *
!* restriction, including without limitation the rights to use,      *
!* copy, modify, merge, publish, distribute, sublicense, and/or      *
!* sell copies of the Software, and to permit persons to whom the    *
!* Software is furnished to do so, subject to the following          *
!* conditions:                                                       *
!*                                                                   *
!*                                                                   *
!* @ Redistributions of source code must retain the above copyright  * 
!*   notice, this list of conditions and the following disclaimers.  *
!*                                                                   * 
!* @ Redistributions in binary form must reproduce the above         *
!*   copyright notice, this list of conditions and the following     *
!*   disclaimers in the documentation and/or other materials         *
!*   provided with the distribution.                                 *
!*                                                                   *
!* @ Neither the names of the Center for Simulation of Advanced      *
!*   Rockets, the University of Illinois, nor the names of its       *
!*   contributors may be used to endorse or promote products derived * 
!*   from this Software without specific prior written permission.   *
!*                                                                   *
!* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
!* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
!* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
!* NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
!* COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       * 
!* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
!* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
!* USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
!*********************************************************************
!* Please acknowledge The University of Illinois Center for          *
!* Simulation of Advanced Rockets in works and publications          *
!* resulting from this software or its derivatives.                  *
!*********************************************************************
SUBROUTINE v3d4_neohookeanincompressdef(coor,matcstet,lmcstet,R_in,d,ci, &
     s11,s22,s33,s12,s23,s13,numnp,nstart,nend,numcstet,numat_vol, &
     xmu,xkappa)
  
!________________________________________________________________________
!
!  V3D4 - Performs displacement based computations for Volumetric 3D 
!         4-node tetrahedra large deformation (i.e. nonlinear) elastic 
!         elements with linear interpolation functions. (constant strain 
!         tetrahedra). Returns the internal force vector R_in.
!
!  DATE: 04.2000                  AUTHOR: SCOT BREITENFELD
!
! Source:
!    "Nonlinear continuum mechanics for finite element analysis"
!     Javier Bonet and Richard D. Wood
!     ISBN 0-521-57272-X

!________________________________________________________________________

  IMPLICIT NONE
!---- Global variables
  INTEGER :: numnp          ! number of nodal points
  INTEGER :: numcstet       ! number of CSTet elements
  INTEGER :: numat_vol      ! number of volumetric materials
!--   coordinate array
  REAL*8, DIMENSION(1:3,1:numnp) :: coor
!--   elastic stiffness consts
  REAL*8, DIMENSION(1:9,1:numat_vol) :: ci
!--   internal force
  REAL*8, DIMENSION(1:3*numnp) :: r_in
!--  displacement vector
  REAL*8, DIMENSION(1:3*numnp) :: d
!--   CSTet stress
  REAL*8, DIMENSION(1:numcstet) :: s11, s22, s33, s12, s23, s13
!--  x, y and z displacements of nodes
  REAL*8 :: u1,u2,u3,u4,v1,v2,v3,v4,w1,w2,w3,w4
!--   coordinate holding variable
  REAL*8 :: x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4
!--   coordinate holding variable
  REAL*8 :: x1d,x2d,x3d,x4d,y1d,y2d,y3d,y4d,z1d,z2d,z3d,z4d
!--  Coordinate subtractions
  REAL*8 :: x14, x24, x34, y14, y24, y34, z14, z24, z34
!--   6*volume, inverse(6*volume),  and the volume      
!      REAL*8 :: Vx6, vol
!--   spacial derivatives
  REAL*8 :: b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12
!--   partial derivatives of the displacement 
  REAL*8 :: dudx,dvdy,dwdz,dudy,dvdx,dvdz,dwdy,dudz,dwdx
!--   strains
  REAL*8 :: e11,e22,e33,e12,e23,e13
! -- connectivity table for CSTet
  INTEGER, DIMENSION(1:4,1:numcstet) :: lmcstet
! -- mat number for CST element
  INTEGER, DIMENSION(1:numcstet) :: matcstet
! -- node numbers
  INTEGER :: n1,n2,n3,n4
! -- dummy and counters
  INTEGER :: i,j,nstart,nend
  INTEGER :: k1n1,k1n2,k1n3,k1n4,k2n1,k2n2,k2n3,k2n4
  INTEGER :: k3n1,k3n2,k3n3,k3n4
! -- 
  REAL*8 :: f11, f12, f13, f21, f22, f23, f31, f32, f33
  REAL*8 :: j1,jsqrd,jsqrdinv
  REAL*8 :: c11, c12, c13, c21, c22, c23, c31, c32, c33
  REAL*8, DIMENSION(1:numat_vol) :: xmu, xkappa
  REAL*8, DIMENSION(1:4,1:4) :: eledb
  REAL*8,DIMENSION(1:3,1:3) :: xj
  REAL*8,DIMENSION(1:3,1:3) :: xjv0
  REAL*8 :: weight = 1.d0/6.d0
  INTEGER :: id, jd, in, ip
  REAL*8 :: xcoor
  REAL*8 :: detjb,detjbv0,evol,theta,press
  REAL*8 :: ic, iiic,vx6old0,vx6old
  REAL*8 :: onethird = 1.d0/3.d0
  REAL*8 :: val11, val21, val31
  
  REAL*8 :: v0x6, v0x6inv
  REAL*8 :: vx6, vx6inv
  REAL*8 :: vol, vol0
  REAL*8 :: term1, term2, cinv
  REAL*8 :: sh1, sh2, sh3, sh4, sh5, sh6, sh7, sh8, sh9, sh10, sh11, sh12
!--  Coordinate subtractions
  REAL*8 :: x14d, x24d, x34d, y14d, y24d, y34d, z14d, z24d, z34d
!--  Coordinate subtractions: These are to speed up B calculation
  REAL*8 :: x12d, x13d, y12d, y13d, z12d, z13d
  REAL*8 :: sigma(3,3),f(3,3),ft(3,3),spk2(3,3), btens(3,3)
  REAL*8 :: rmat3_det, detf,jinv
  
  DO i = nstart, nend
     j = matcstet(i)
     
     n1 = lmcstet(1,i)
     n2 = lmcstet(2,i)
     n3 = lmcstet(3,i)
     n4 = lmcstet(4,i)
     
     k3n1 = 3*n1
     k3n2 = 3*n2
     k3n3 = 3*n3
     k3n4 = 3*n4
     
     k2n1 = k3n1 - 1
     k2n2 = k3n2 - 1
     k2n3 = k3n3 - 1
     k2n4 = k3n4 - 1
     
     k1n1 = k3n1 - 2
     k1n2 = k3n2 - 2
     k1n3 = k3n3 - 2
     k1n4 = k3n4 - 2     
     ! k#n# dummy variables replaces:
     u1 = d(k1n1)           ! (3*n1-2)
     u2 = d(k1n2)           ! (3*n2-2)
     u3 = d(k1n3)           ! (3*n3-2)
     u4 = d(k1n4)           ! (3*n4-2)
     v1 = d(k2n1)           ! (3*n1-1)
     v2 = d(k2n2)           ! (3*n2-1)
     v3 = d(k2n3)           ! (3*n3-1)
     v4 = d(k2n4)           ! (3*n4-1)
     w1 = d(k3n1)           ! (3*n1)
     w2 = d(k3n2)           ! (3*n2)
     w3 = d(k3n3)           ! (3*n3)
     w4 = d(k3n4)           ! (3*n4)
     
     x1 = coor(1,n1)
     x2 = coor(1,n2)
     x3 = coor(1,n3)
     x4 = coor(1,n4)
     y1 = coor(2,n1)
     y2 = coor(2,n2)
     y3 = coor(2,n3)
     y4 = coor(2,n4)
     z1 = coor(3,n1)
     z2 = coor(3,n2)
     z3 = coor(3,n3)
     z4 = coor(3,n4)
     
     x1d = x1 + u1
     x2d = x2 + u2
     x3d = x3 + u3
     x4d = x4 + u4
     y1d = y1 + v1
     y2d = y2 + v2
     y3d = y3 + v3
     y4d = y4 + v4
     z1d = z1 + w1
     z2d = z2 + w2
     z3d = z3 + w3
     z4d = z4 + w4
!
!     creates the Jacobian matrix using equation ???
!     Using undeformed coordinates
     
     x14 = x1 - x4
     x24 = x2 - x4
     x34 = x3 - x4
     y14 = y1 - y4
     y24 = y2 - y4
     y34 = y3 - y4
     z14 = z1 - z4
     z24 = z2 - z4
     z34 = z3 - z4
     
     val11 =    y24*z34 - z24*y34
     val21 = -( x24*z34 - z24*x34 )
     val31 =    x24*y34 - y24*x34
     
     v0x6 = -( x14*val11 + y14*val21 + z14*val31 )
!
!     creates the Jacobian matrix using equation ???
!     Determines the determinant of xj and checks for zero values
!
!     Using deformed coordinates

     x12d = x1d - x2d          ! not used in vol. calc
     x13d = x1d - x3d          ! not used in vol. calc
     x14d = x1d - x4d
     x24d = x2d - x4d
     x34d = x3d - x4d
     y12d = y1d - y2d          ! not used in vol. calc
     y13d = y1d - y3d          ! not used in vol. calc
     y14d = y1d - y4d
     y24d = y2d - y4d
     y34d = y3d - y4d
     z12d = z1d - z2d          ! not used in vol. calc
     z13d = z1d - z3d          ! not used in vol. calc
     z14d = z1d - z4d
     z24d = z2d - z4d
     z34d = z3d - z4d
     
     val11 =    y24d*z34d - z24d*y34d
     val21 = -( x24d*z34d - z24d*x34d )
     val31 =    x24d*y34d - y24d*x34d
     
     vx6 = -( x14d*val11 + y14d*val21 + z14d*val31 )
     
     IF(vx6.LE.0.d0) THEN
        WRITE(*,100) i
        stop
     ENDIF
     

! calculate the volume
      
     vol0 = v0x6/6.d0
     vol  = vx6/6.d0
     
     v0x6inv = 1.d0/v0x6
     vx6inv = 1.d0/vx6
     
! See the maple worksheet 'V3D4.mws' for the derivation of Vx6

! See the maple worksheet 'V3D4.mws' for the derivation of [B]
! NOTE: Factored for a more equivalent/compact form then maple's

     b1  = ( (y3-y4)*(z2-z4) - (y2-y4)*(z3-z4) ) * v0x6inv
     b2  = ( (z3-z4)*(x2-x4) - (z2-z4)*(x3-x4) ) * v0x6inv
     b3  = ( (x3-x4)*(y2-y4) - (x2-x4)*(y3-y4) ) * v0x6inv
     b4  = ( (y1-y3)*(z1-z4) - (y1-y4)*(z1-z3) ) * v0x6inv
     b5  = ( (z1-z3)*(x1-x4) - (z1-z4)*(x1-x3) ) * v0x6inv
     b6  = ( (x1-x3)*(y1-y4) - (x1-x4)*(y1-y3) ) * v0x6inv
     b7  = ( (y1-y4)*(z1-z2) - (y1-y2)*(z1-z4) ) * v0x6inv
     b8  = ( (z1-z4)*(x1-x2) - (z1-z2)*(x1-x4) ) * v0x6inv
     b9  = ( (x1-x4)*(y1-y2) - (x1-x2)*(y1-y4) ) * v0x6inv
     b10 = ( (y1-y2)*(z1-z3) - (y1-y3)*(z1-z2) ) * v0x6inv
     b11 = ( (z1-z2)*(x1-x3) - (z1-z3)*(x1-x2) ) * v0x6inv
     b12 = ( (x1-x2)*(y1-y3) - (x1-x3)*(y1-y2) ) * v0x6inv
! 
! deformation gradients F
!
     f11 = 1.d0 + ( b1*u1 + b4*u2 + b7*u3 + b10*u4 ) ! 1 + ( dudx )
     f22 = 1.d0 + ( b2*v1 + b5*v2 + b8*v3 + b11*v4 ) ! 1 + ( dvdy )
     f33 = 1.d0 + ( b3*w1 + b6*w2 + b9*w3 + b12*w4 ) ! 1 + ( dwdz )
     f12 = b2*u1 + b5*u2 + b8*u3 + b11*u4 ! dudy
     f21 = b1*v1 + b4*v2 + b7*v3 + b10*v4 ! dvdx
     f23 = b3*v1 + b6*v2 + b9*v3 + b12*v4 ! dvdz
     f32 = b2*w1 + b5*w2 + b8*w3 + b11*w4 ! dwdy
     f13 = b3*u1 + b6*u2 + b9*u3 + b12*u4 ! dudz
     f31 = b1*w1 + b4*w2 + b7*w3 + b10*w4 ! dwdx
     
     theta = vol/vol0
     
     press = xkappa(j)*(theta-1.d0)
     
     j1 = f11*(f22*f33-f23*f32)+f12*(f31*f23-f21*f33)+f13*(-f31*f22+f21*f32)
     IF(j1.LE.0.d0) THEN
        WRITE(*,*)'Volume has become zero for element',i
        stop
     ENDIF
     jsqrd = j1*j1
     
! first invariant of C
     
     ic = f11**2+f21**2+f31**2+f12**2+f22**2+f32**2+f13**2+f23**2+f33**2
     
! Third invariant of C
     
     iiic = jsqrd        
     
!          T
!     C = F  F = right Cauchy-Green deformation tensor
!
     c11 = (f11*f11+f21*f21+f31*f31)
     c12 = (f11*f12+f21*f22+f31*f32)
     c13 = (f11*f13+f21*f23+f31*f33)
     c21 = c12 
     c22 = (f12*f12+f22*f22+f32*f32)
     c23 = (f12*f13+f22*f23+f32*f33) 
     c31 = c13
     c32 = c23
     c33 = (f13*f13+f23*f23+f33*f33)
!
!           T
!    b = F F  = left Cauchy-Green deformation tensor
!

     btens(1,1) = f11**2+f12**2+f13**2
     btens(1,2) = f11*f21+f12*f22+f13*f23
     btens(1,3) = f11*f31+f12*f32+f13*f33
     btens(2,1) = btens(1,2)
     btens(2,2) = f21**2+f22**2+f23**2
     btens(2,3) = f21*f31+f22*f32+f23*f33
     btens(3,1) = btens(1,3)
     btens(3,2) = btens(2,3)
     btens(3,3) = f31**2+f32**2+f33**2
!
!    Cauchy stress

     call neoinccauchystress(3,xmu(j),j1,btens,sigma)
     call addpres(3,press,sigma)

!
! Second Piola-Kirchoff tensor
! Eq. (5.28), pg. 124 

     jsqrdinv = 1.d0/jsqrd
     
     term1 = xmu(j)*iiic**(-onethird)
     term2 = press*j1
     
     cinv = (c22*c33-c23*c32)*jsqrdinv
     s11(i) = term1*(1.d0-onethird*ic*cinv) + term2*cinv
     cinv = (c11*c33-c31*c13)*jsqrdinv
     s22(i) = term1*(1.d0-onethird*ic*cinv) + term2*cinv
     cinv = (c11*c22-c12*c21)*jsqrdinv
     s33(i) = term1*(1.d0-onethird*ic*cinv) + term2*cinv
     cinv = (-c12*c33+c13*c32)*jsqrdinv
     s12(i) = cinv*(-term1*onethird*ic + term2)
     cinv = (-c11*c23+c13*c21)*jsqrdinv
     s23(i) = cinv*(-term1*onethird*ic + term2)
     cinv = ( c12*c23-c13*c22)*jsqrdinv
     s13(i) = cinv*(-term1*onethird*ic + term2)

     spk2(1,1) = s11(i)
     spk2(1,2) = s12(i)
     spk2(1,3) = s13(i)
     spk2(2,1) = spk2(1,2)
     spk2(2,2) = s22(i)
     spk2(2,3) = s23(i)
     spk2(3,1) = spk2(1,3)
     spk2(3,2) = spk2(2,3)
     spk2(3,3) = s33(i)

     f(1,1) = f11
     f(1,2) = f12
     f(1,3) = f13
     f(2,1) = f21
     f(2,2) = f22
     f(2,3) = f23
     f(3,1) = f31
     f(3,2) = f32
     f(3,3) = f33

     ft(:,:) = transpose(f(:,:))

     detf = rmat3_det(f)

     jinv=1.d0/detf

!     print*,'lkjk',sigma(1,1),sigma(2,2), sigma(3,3)
!     SIGMA(:,:) = Jinv*F(:,:)*SPK2(:,:)*FT(:,:)
!     print*,sigma(1,1),sigma(2,2), sigma(3,3)

     sh1  = (y34d*z24d - y24d*z34d) * vx6inv
     sh2  = (z34d*x24d - z24d*x34d) * vx6inv
     sh3  = (x34d*y24d - x24d*y34d) * vx6inv
     sh4  = (y13d*z14d - y14d*z13d) * vx6inv
     sh5  = (z13d*x14d - z14d*x13d) * vx6inv
     sh6  = (x13d*y14d - x14d*y13d) * vx6inv
     sh7  = (y14d*z12d - y12d*z14d) * vx6inv
     sh8  = (z14d*x12d - z12d*x14d) * vx6inv
     sh9  = (x14d*y12d - x12d*y14d) * vx6inv
     sh10 = (y12d*z13d - y13d*z12d) * vx6inv
     sh11 = (z12d*x13d - z13d*x12d) * vx6inv
     sh12 = (x12d*y13d - x13d*y12d) * vx6inv

! ASSEMBLE THE INTERNAL FORCE VECTOR
!
! local node 1
     r_in(k1n1) = r_in(k1n1) - vol* &
          ( sigma(1,1)*sh1 + sigma(1,2)*sh2 + sigma(1,3)*sh3 )
     r_in(k2n1) = r_in(k2n1) - vol* &
          ( sigma(1,2)*sh1 + sigma(2,2)*sh2 + sigma(2,3)*sh3 )
     r_in(k3n1) = r_in(k3n1) - vol* &
          ( sigma(1,3)*sh1 + sigma(2,3)*sh2 + sigma(3,3)*sh3 )
! local node 2 
     r_in(k1n2) = r_in(k1n2) - vol* &
          ( sigma(1,1)*sh4 + sigma(1,2)*sh5 + sigma(1,3)*sh6 )
     r_in(k2n2) = r_in(k2n2) - vol* &
          ( sigma(1,2)*sh4 + sigma(2,2)*sh5 + sigma(2,3)*sh6 )
     r_in(k3n2) = r_in(k3n2) - vol* &
          ( sigma(1,3)*sh4 + sigma(2,3)*sh5 + sigma(3,3)*sh6 )
! local node 3 
     r_in(k1n3) = r_in(k1n3) - vol* &
          ( sigma(1,1)*sh7 + sigma(1,2)*sh8 + sigma(1,3)*sh9 )
     r_in(k2n3) = r_in(k2n3) - vol* &
          ( sigma(1,2)*sh7 + sigma(2,2)*sh8 + sigma(2,3)*sh9 )
     r_in(k3n3) = r_in(k3n3) - vol* &
          ( sigma(1,3)*sh7 + sigma(2,3)*sh8 + sigma(3,3)*sh9 )
! local node 4  
     r_in(k1n4) = r_in(k1n4) - vol* &
          ( sigma(1,1)*sh10 + sigma(1,2)*sh11 + sigma(1,3)*sh12 )
     r_in(k2n4) = r_in(k2n4) - vol* &
          ( sigma(1,2)*sh10 + sigma(2,2)*sh11 + sigma(2,3)*sh12 )
     r_in(k3n4) = r_in(k3n4) - vol* &
          ( sigma(1,3)*sh10 + sigma(2,3)*sh11 + sigma(3,3)*sh12 )

     
  ENDDO
  
  RETURN
      
100 FORMAT(' Negative Jacobian for element: ',i10)
END SUBROUTINE v3d4_neohookeanincompressdef


FUNCTION rmat3_det ( a )
!
!*******************************************************************************
!
!! RMAT3_DET computes the determinant of a 3 by 3 matrix.
!
!
!  The determinant is the volume of the shape spanned by the vectors
!  making up the rows or columns of the matrix.
!
!  Formula:
!
!    det = a11 * a22 * a33 - a11 * a23 * a32
!        + a12 * a23 * a31 - a12 * a21 * a33
!        + a13 * a21 * a32 - a13 * a22 * a31
!
!  Parameters:
!
!    Input, real A(3,3), the matrix whose determinant is desired.
!
!    Output, real RMAT3_DET, the determinant of the matrix.
!
        IMPLICIT NONE
!
        REAL*8 a(3,3)
        REAL*8 rmat3_det
!
        rmat3_det =   a(1,1) * ( a(2,2) * a(3,3) - a(2,3) * a(3,2) ) &
             + a(1,2) * ( a(2,3) * a(3,1) - a(2,1) * a(3,3) ) &
             + a(1,3) * ( a(2,1) * a(3,2) - a(2,2) * a(3,1) )

        RETURN
      END FUNCTION rmat3_det

      SUBROUTINE neoinccauchystress(ndime,xmu,detf,btens,sigma)
!
!-----------------------------------------------------------------------
!
!     Determines the Cauchy stress tensor for nearly incompressible 
!     neo-Hookean materials
!     
!
!     ndime  -->  number of dimensions
!     xmu    -->  mu coefficient
!     detf   -->  determinant of F
!     btens  -->  right Cauchy-Green tensor
!     sigma  -->  Cauchy stress tensor
!
!-----------------------------------------------------------------------
!
        IMPLICIT NONE
        
        REAL*8, DIMENSION(1:3,1:3) :: sigma, btens
        REAL*8 :: detf, xmu
        INTEGER :: ndime
        
        REAL*8 :: trace, xm
        INTEGER :: id, jd
!
!     Obtains the Cauchy stress using equation ???
!     First finds the trace
!
      trace = 0.d0
      DO id = 1, 3
         trace = trace + btens(id,id)
      ENDDO

      xm = xmu/(detf**(5./3.))
      DO id = 1, ndime      
         DO jd = 1, ndime
            sigma(id,jd) = xm*btens(id,jd)
         ENDDO
         sigma(id,id)=sigma(id,id)-xm*trace/3.
      ENDDO

      RETURN
    END SUBROUTINE neoinccauchystress

    subroutine addpres(ndime,press,sigma)

!-----------------------------------------------------------------------
!
!     Adds the pressure to the deviatoric stresses
!
!     ndime  -->  number of dimensions
!     press  -->  pressure
!     sigma  -->  Cauchy stress tensor
!
!-----------------------------------------------------------------------
!

      implicit none

      integer :: ndime
      real*8, DIMENSION(1:3,1:3) :: sigma
      real*8 :: press

      integer :: id
!
!     Obtains the Cauchy stress from deviatoric and pressure
!     components using equation ???
!
      do id = 1,ndime      
         sigma(id,id) = sigma(id,id) + press
      enddo
      return
      end