v3d4_nl_arruda_boyce.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_nl_arruda_boyce(coor,matcstet,lmcstet,R_in,d, &
      s11,s22,s33,s12,s23,s13, &
      numnp,nstart,nend,numcstet,numat_vol, &
      mu,kapp)

!!****f* Rocfrac/Rocfrac/Source/v3d4_NL_ARRUDA_BOYCE.f90
!!
!!  NAME
!!    v3d4_NL_ARRUDA_BOYCE
!!
!!  FUNCTION
!!
!!    Computes the internal force for a large deformation, 
!!    hyperelastic arruda-boyce 4-node tetrahedral.
!!
!!  INPUTS
!!
!!   NumNP    -- Number of nodes
!!   numcstet -- Number of elements
!!   Coor     -- number of coordinates
!!   Matcstet  -- Material id
!!   d -- nodal displacement
!!   mu -- shear modulus
!!   kappa -- bulk modulus
!!
!!   lmcstet  -- Nodal connectivity
!!   nstart, nend -- element beginning and end loop counter
!!   numat_vol -- number of materials
!!
!!  OUTPUT
!!
!!    Rin -- internal force vector
!!    S11,S22,S33,S12,S23,S13 -- Stresses at gauss points
!!
!!****


      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
!--   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
!--   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
!--  Coordinate subtractions
      REAL*8 :: x14, x24, x34, y14, y24, y34, z14, z24, z34
!--  Coordinate subtractions: These are to speed up B calculation
      REAL*8 :: x12, x13, y12, y13, z12, z13
      REAL*8 :: val11, val21, val31
! -- 
      REAL*8 :: f11, f12, f13, f21, f22, f23, f31, f32, f33
      REAL*8, DIMENSION(1:numat_vol) :: mu, kapp
      REAL*8, DIMENSION(1:3,1:3) :: fij, cij

      INTEGER :: kk,ll,mm

      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)


         x12 = x1 - x2          ! not used in vol. calc
         x13 = x1 - x3          ! not used in vol. calc
         x14 = x1 - x4
         x24 = x2 - x4
         x34 = x3 - x4
         y12 = y1 - y2          ! not used in vol. calc
         y13 = y1 - y3          ! not used in vol. calc
         y14 = y1 - y4
         y24 = y2 - y4
         y34 = y3 - y4
         z12 = z1 - z2          ! not used in vol. calc
         z13 = z1 - z3          ! not used in vol. calc
         z14 = z1 - z4
         z24 = z2 - z4
         z34 = z3 - z4

         val11 =    y24*z34 - z24*y34
         val21 = -( x24*z34 - z24*x34 )
         val31 =    x24*y34 - y24*x34

         vx6 = -( x14*val11 + y14*val21 + z14*val31 )

         vx6inv = 1.d0 / vx6

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

         b1  = (y34*z24 - y24*z34) * vx6inv
         b2  = (z34*x24 - z24*x34) * vx6inv
         b3  = (x34*y24 - x24*y34) * vx6inv
         b4  = (y13*z14 - y14*z13) * vx6inv
         b5  = (z13*x14 - z14*x13) * vx6inv
         b6  = (x13*y14 - x14*y13) * vx6inv
         b7  = (y14*z12 - y12*z14) * vx6inv
         b8  = (z14*x12 - z12*x14) * vx6inv
         b9  = (x14*y12 - x12*y14) * vx6inv
         b10 = (y12*z13 - y13*z12) * vx6inv
         b11 = (z12*x13 - z13*x12) * vx6inv
         b12 = (x12*y13 - x13*y12) * vx6inv
! 
! 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
!
! Arruda-Boyce Nonlinear Elasticity Model 
!             
! -- NOTE: ci(7,j) : shear modulus 


!--   (1) Deformation Gradient Fij 
  
         fij(1,1) = f11
         fij(1,2) = f12
         fij(1,3) = f13 
         fij(2,1) = f21 
         fij(2,2) = f22
         fij(2,3) = f23 
         fij(3,1) = f31
         fij(3,2) = f32
         fij(3,3) = f33

!--  (2) right Cauchy-Green deformation Tensor Cij = Fmi Fmj 
!             T
!        C = F  F
 
         DO kk=1,3
            DO ll=1,3
               cij(kk,ll) = 0.d0
               DO mm=1,3
                  cij(kk,ll)=cij(kk,ll)+fij(mm,kk)*fij(mm,ll)
               ENDDO
            ENDDO
         ENDDO

         CALL arruda_boyce(cij, &
         s11(i),s22(i),s33(i),s12(i),s23(i),s13(i),i, &
         mu(j),kapp(j))
         
! 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*f11 + s22(i)*b2*f12 + s33(i)*b3*f13 &
             + s12(i)*( b2*f11 + b1*f12 ) &
             + s23(i)*( b3*f12 + b2*f13 ) &
             + s13(i)*( b3*f11 + b1*f13 ) )
         r_in(k2n1) = r_in(k2n1) - vol* &
             ( s11(i)*b1*f21 + s22(i)*b2*f22 + s33(i)*b3*f23 &
             + s12(i)*( b1*f22 + b2*f21 ) &
             + s23(i)*( b3*f22 + b2*f23 ) &
             + s13(i)*( b3*f21 + b1*f23 ) )
         r_in(k3n1)   = r_in(k3n1)   - vol* &
             ( s11(i)*b1*f31 + s22(i)*b2*f32 + s33(i)*b3*f33 &
             + s12(i)*( b2*f31 + b1*f32 ) &
             + s23(i)*( b3*f32 + b2*f33 ) &
             + s13(i)*( b3*f31 + b1*f33 ) )
! local node 2 
         r_in(k1n2) = r_in(k1n2) - vol* &
             ( s11(i)*b4*f11 + s22(i)*b5*f12 + s33(i)*b6*f13 &
             + s12(i)*( b5*f11 + b4*f12 ) &
             + s23(i)*( b6*f12 + b5*f13 ) &
             + s13(i)*( b6*f11 + b4*f13 ) )
         r_in(k2n2) = r_in(k2n2) - vol* &
             ( s11(i)*b4*f21 + s22(i)*b5*f22 + s33(i)*b6*f23 &
             + s12(i)*( b4*f22 + b5*f21 ) &
             + s23(i)*( b6*f22 + b5*f23 ) &
             + s13(i)*( b6*f21 + b4*f23 ) )
         r_in(k3n2)   = r_in(k3n2)   - vol* &
             ( s11(i)*b4*f31 + s22(i)*b5*f32 + s33(i)*b6*f33 &
             + s12(i)*( b5*f31 + b4*f32 ) &
             + s23(i)*( b6*f32 + b5*f33 ) &
             + s13(i)*( b6*f31 + b4*f33 ) )
! local node 3 
         r_in(k1n3) = r_in(k1n3) - vol* &
             ( s11(i)*b7*f11 + s22(i)*b8*f12 + s33(i)*b9*f13 &
             + s12(i)*( b8*f11 + b7*f12 ) &
             + s23(i)*( b9*f12 + b8*f13 ) &
             + s13(i)*( b9*f11 + b7*f13 ) )
         r_in(k2n3) = r_in(k2n3) - vol* &
             ( s11(i)*b7*f21 + s22(i)*b8*f22 + s33(i)*b9*f23 &
             + s12(i)*( b7*f22 + b8*f21 ) &
             + s23(i)*( b9*f22 + b8*f23 ) &
             + s13(i)*( b9*f21 + b7*f23 ) )
         r_in(k3n3)   = r_in(k3n3)   - vol* &
             ( s11(i)*b7*f31 + s22(i)*b8*f32 + s33(i)*b9*f33 &
             + s12(i)*( b8*f31 + b7*f32 ) &
             + s23(i)*( b9*f32 + b8*f33 ) &
             + s13(i)*( b9*f31 + b7*f33 ) )
! local node 4         
         r_in(k1n4) = r_in(k1n4) - vol* &
             ( s11(i)*b10*f11 + s22(i)*b11*f12+s33(i)*b12*f13 &
             + s12(i)*( b11*f11 + b10*f12 ) &
             + s23(i)*( b12*f12 + b11*f13 ) &
             + s13(i)*( b12*f11 + b10*f13 ) )
         r_in(k2n4) = r_in(k2n4) - vol* &
             ( s11(i)*b10*f21 + s22(i)*b11*f22+s33(i)*b12*f23 &
             + s12(i)*( b10*f22 + b11*f21 ) &
             + s23(i)*( b12*f22 + b11*f23 ) &
             + s13(i)*( b12*f21 + b10*f23 ) )
         r_in(k3n4)   = r_in(k3n4)   - vol* &
             ( s11(i)*b10*f31 + s22(i)*b11*f32+s33(i)*b12*f33 &
             + s12(i)*( b11*f31 + b10*f32 ) &
             + s23(i)*( b12*f32 + b11*f33 ) &
             + s13(i)*( b12*f31 + b10*f33 ) )

      x1 = coor(1,n1) + u1
      x2 = coor(1,n2) + u2
      x3 = coor(1,n3) + u3
      x4 = coor(1,n4) + u4
      y1 = coor(2,n1) + v1
      y2 = coor(2,n2) + v2
      y3 = coor(2,n3) + v3
      y4 = coor(2,n4) + v4
      z1 = coor(3,n1) + w1
      z2 = coor(3,n2) + w2
      z3 = coor(3,n3) + w3
      z4 = coor(3,n4) + w4
     
      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
     
      vx6 = -( x14*val11 + y14*val21 + z14*val31 )

      IF(vx6.LT.0.d0) THEN
         print*,'ROCFRAC ::        ERROR * ERROR * ERROR'
         print*,'ROCFRAC :: DEFORMED ELEMENT VOLUME TURNED NEGATIVE'
         stop
      ENDIF
          
      ENDDO

      RETURN
      END