v3d4_NeoHookeanCompress.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_neohookeancompress(coor,matcstet,lmcstet,R_in,d,ci, &
     s11,s22,s33,s12,s23,s13,numnp,nstart,nend,numcstet,numat_vol, &
     xmu,xlamda)

!________________________________________________________________________
!
!  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
  REAL*8  :: xmu
!--   coordinate array
  REAL*8, DIMENSION(1:3,1:numnp) :: coor
  REAL*8  :: xlamda
!--   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 subtractions
  REAL*8 :: x14, x24, x34, y14, y24, y34, z14, z24, z34
!--  
!  REAL*8 :: c11, c21, c31
!--   6*volume, inverse(6*volume),  and the volume      
  REAL*8 :: vx6, vx6inv, 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,j2,j2inv
  REAL*8 :: c11, c12, c13, c21, c22, c23, c31, c32, c33

  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)
     
     
     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
     
     c11 =    y24*z34 - z24*y34
     c21 = -( x24*z34 - z24*x34 )
     c31 =    x24*y34 - y24*x34
     
     vx6 = -( x14*c11 + y14*c21 + z14*c31 )

! See the maple worksheet 'V3D4.mws' for the derivation of Vx6

!         Vx6 = x2*y3*z4 - x2*y4*z3 - y2*x3*z4 + y2*x4*z3 + z2*x3*y4 
!     $        - z2*x4*y3 - x1*y3*z4 + x1*y4*z3 + x1*y2*z4 - x1*y2*z3
!     $        - x1*z2*y4 + x1*z2*y3 + y1*x3*z4 - y1*x4*z3 - y1*x2*z4 
!     $        + y1*x2*z3 + y1*z2*x4 - y1*z2*x3 - z1*x3*y4 + z1*x4*y3
!     $        + z1*x2*y4 - z1*x2*y3 - z1*y2*x4 + z1*y2*x3

     vx6inv = 1.d0 / vx6

! See the maple worksheet 'V3D4.mws' for the derivation of [B]

     b1  = -(y3*z4 - y4*z3 - y2*z4 + y2*z3 + z2*y4 - z2*y3) * vx6inv
     b2  =  (x3*z4 - x4*z3 - x2*z4 + x2*z3 + z2*x4 - z2*x3) * vx6inv
     b3  = -(x3*y4 - x4*y3 - x2*y4 + x2*y3 + y2*x4 - y2*x3) * vx6inv
     b4  =  (y3*z4 - y4*z3 - y1*z4 + y1*z3 + z1*y4 - z1*y3) * vx6inv
     b5  = -(x3*z4 - x4*z3 - x1*z4 + x1*z3 + z1*x4 - z1*x3) * vx6inv
     b6  =  (x3*y4 - x4*y3 - x1*y4 + x1*y3 + y1*x4 - y1*x3) * vx6inv
     b7  = -(y2*z4 - z2*y4 - y1*z4 + y1*z2 + z1*y4 - z1*y2) * vx6inv
     b8  =  (x2*z4 - z2*x4 - x1*z4 + x1*z2 + z1*x4 - z1*x2) * vx6inv
     b9  = -(x2*y4 - y2*x4 - x1*y4 + x1*y2 + y1*x4 - y1*x2) * vx6inv
     b10 =  (y2*z3 - z2*y3 - y1*z3 + y1*z2 + z1*y3 - z1*y2) * vx6inv
     b11 = -(x2*z3 - z2*x3 - x1*z3 + x1*z2 + z1*x3 - z1*x2) * vx6inv
     b12 =  (x2*y3 - y2*x3 - x1*y3 + x1*y2 + y1*x3 - y1*x2) * vx6inv
!
! Calculate displacement gradient (H)
!
     dudx = b1*u1 + b4*u2 + b7*u3 + b10*u4
     dvdy = b2*v1 + b5*v2 + b8*v3 + b11*v4
     dwdz = b3*w1 + b6*w2 + b9*w3 + b12*w4
     dudy = b2*u1 + b5*u2 + b8*u3 + b11*u4
     dvdx = b1*v1 + b4*v2 + b7*v3 + b10*v4
     dvdz = b3*v1 + b6*v2 + b9*v3 + b12*v4
     dwdy = b2*w1 + b5*w2 + b8*w3 + b11*w4
     dudz = b3*u1 + b6*u2 + b9*u3 + b12*u4
     dwdx = b1*w1 + b4*w2 + b7*w3 + b10*w4
! 
! deformation gradients F
!
     f11 = 1.d0 + dudx
     f12 = dudy
     f13 = dudz
     f21 = dvdx
     f22 = 1.d0 + dvdy
     f23 = dvdz
     f31 = dwdx
     f32 = dwdy
     f33 = 1.d0 + dwdz

     j1 = f11*(f22*f33-f23*f32)+f12*(f31*f23-f21*f33)+f13*(-f31*f22+f21*f32)
     IF(j1.LE.0.d0) THEN
        WRITE(*,*)'area has become zero for element',i
        stop
     ENDIF
     j2 = j1*j1

!          T
!     C = F  F = right Cauchy-Green deformation tensor
!
!     Divided by third invariant squared

     j2inv = 1.d0/j2
     
     c11 = (f11*f11+f21*f21+f31*f31)*j2inv
     c12 = (f11*f12+f21*f22+f31*f32)*j2inv
     c13 = (f11*f13+f21*f23+f31*f33)*j2inv
     c21 = c12
     c22 = (f12*f12+f22*f22+f32*f32)*j2inv
     c23 = (f12*f13+f22*f23+f32*f33)*j2inv
     c31 = c13
     c32 = c23
     c33 = (f13*f13+f23*f23+f33*f33)*j2inv       

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

     s11(i) = xmu*(1.d0-(c22*c33-c23*c32)) + xlamda*log(j1)*(c22*c33-c23*c32)
     s22(i) = xmu*(1.d0-(c11*c33-c31*c13)) + xlamda*log(j1)*(c11*c33-c31*c13)
     s33(i) = xmu*(1.d0-(c11*c22-c12*c21)) + xlamda*log(j1)*(c11*c22-c12*c21)
     s12(i) = (-c12*c33+c13*c32)*(xmu - xlamda*log(j1))
     s23(i) = (-c11*c23+c13*c21)*(xmu - xlamda*log(j1))
     s13(i) = (c12*c23-c13*c22)*(xmu - xlamda*log(j1))
         

! calculate the volume

     vol = vx6 / 6.d0

! ASSEMBLE THE INTERNAL FORCE VECTOR
!
! local node 1
     r_in(k1n1) = r_in(k1n1) - vol* &
          ( s11(i)*b1*(1.d0+dudx) + s22(i)*b2*dudy + s33(i)*b3*dudz &
          + s12(i)*( b2*(1.d0+dudx) + b1*dudy ) &
          + s23(i)*( b3*dudy + b2*dudz ) &
          + s13(i)*( b3*(1.d0+dudx) + b1*dudz ) )
     r_in(k2n1) = r_in(k2n1) - vol* &
          ( s11(i)*b1*dvdx + s22(i)*b2*(1.d0+dvdy) + s33(i)*b3*dvdz &
          + s12(i)*( b1*(1.d0+dvdy) + b2*dvdx ) &
          + s23(i)*( b3*(1.d0+dvdy) + b2*dvdz ) &
          + s13(i)*( b3*dvdx + b1*dvdz ) )
     r_in(k3n1)   = r_in(k3n1)   - vol* &
          ( s11(i)*b1*dwdx + s22(i)*b2*dwdy + s33(i)*b3*(1.d0+dwdz) &
          + s12(i)*( b2*dwdx + b1*dwdy ) &
          + s23(i)*( b3*dwdy + b2*(1.d0 + dwdz) ) &
          + s13(i)*( b3*dwdx + b1*(1.d0 + dwdz) ) )
! local node 2 
     r_in(k1n2) = r_in(k1n2) - vol* &
          ( s11(i)*b4*(1.d0+dudx) + s22(i)*b5*dudy + s33(i)*b6*dudz &
          + s12(i)*( b5*(1.d0+dudx) + b4*dudy ) &
          + s23(i)*( b6*dudy + b5*dudz ) &
          + s13(i)*( b6*(1.d0+dudx) + b4*dudz ) )
     r_in(k2n2) = r_in(k2n2) - vol* &
          ( s11(i)*b4*dvdx + s22(i)*b5*(1.d0+dvdy) + s33(i)*b6*dvdz &
          + s12(i)*( b4*(1.d0+dvdy) + b5*dvdx ) &
          + s23(i)*( b6*(1.d0+dvdy) + b5*dvdz ) &
          + s13(i)*( b6*dvdx + b4*dvdz ) )
     r_in(k3n2)   = r_in(k3n2)   - vol* &
          ( s11(i)*b4*dwdx + s22(i)*b5*dwdy + s33(i)*b6*(1.d0+dwdz) &
          + s12(i)*( b5*dwdx + b4*dwdy ) &
          + s23(i)*( b6*dwdy + b5*(1.d0 + dwdz) ) &
          + s13(i)*( b6*dwdx + b4*(1.d0 + dwdz) ) )
! local node 3 
     r_in(k1n3) = r_in(k1n3) - vol* &
          ( s11(i)*b7*(1.d0+dudx) + s22(i)*b8*dudy + s33(i)*b9*dudz &
          + s12(i)*( b8*(1.d0+dudx) + b7*dudy ) &
          + s23(i)*( b9*dudy + b8*dudz ) &
          + s13(i)*( b9*(1.d0+dudx) + b7*dudz ) ) 
     r_in(k2n3) = r_in(k2n3) - vol* &
          ( s11(i)*b7*dvdx + s22(i)*b8*(1.d0+dvdy) + s33(i)*b9*dvdz &
          + s12(i)*( b7*(1.d0+dvdy) + b8*dvdx ) &
          + s23(i)*( b9*(1.d0+dvdy) + b8*dvdz ) &
          + s13(i)*( b9*dvdx + b7*dvdz ) )
     r_in(k3n3)   = r_in(k3n3)   - vol* &
          ( s11(i)*b7*dwdx + s22(i)*b8*dwdy + s33(i)*b9*(1.d0+dwdz) &
          + s12(i)*( b8*dwdx + b7*dwdy ) &
          + s23(i)*( b9*dwdy + b8*(1.d0 + dwdz) ) &
          + s13(i)*( b9*dwdx + b7*(1.d0 + dwdz) ) )
! local node 4         
     r_in(k1n4) = r_in(k1n4) - vol* &
          ( s11(i)*b10*(1.d0+dudx) + s22(i)*b11*dudy+s33(i)*b12*dudz &
          + s12(i)*( b11*(1.d0+dudx) + b10*dudy ) &
          + s23(i)*( b12*dudy + b11*dudz ) &
          + s13(i)*( b12*(1.d0+dudx) + b10*dudz ) )
     r_in(k2n4) = r_in(k2n4) - vol* &
          ( s11(i)*b10*dvdx + s22(i)*b11*(1.d0+dvdy)+s33(i)*b12*dvdz &
          + s12(i)*( b10*(1.d0+dvdy) + b11*dvdx ) &
          + s23(i)*( b12*(1.d0+dvdy) + b11*dvdz ) &
          + s13(i)*( b12*dvdx + b10*dvdz ) )
     r_in(k3n4)   = r_in(k3n4)   - vol* &
          ( s11(i)*b10*dwdx + s22(i)*b11*dwdy+s33(i)*b12*(1.d0+dwdz) &
          + s12(i)*( b11*dwdx + b10*dwdy ) &
          + s23(i)*( b12*dwdy + b11*(1.d0 + dwdz) ) &
          + s13(i)*( b12*dwdx + b10*(1.d0 + dwdz) ) )

          
  ENDDO
  
  RETURN
END SUBROUTINE v3d4_neohookeancompress