MP/Source/HomogenizedMat.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 homogenizedmat(Density, Em, PRm, Ep, PRp, cm, cd_fastest )

  IMPLICIT NONE


  INTEGER :: myid
  REAL*8 :: density
  REAL*8 :: em, gm, prm, ep,prp, gp
  REAL*8 :: cm, c2

  REAL*8 :: cd_fastest

  REAL*8 :: e_homog
  REAL*8 :: pr_homog

  INTEGER, parameter :: numatv = 2

  REAL*8, DIMENSION(1:6,1:6,1:numatv) :: l_tensor, m_tensor


  REAL*8, DIMENSION(1:6,1:6) :: l_bar, m_bar, lo

  REAL*8, parameter :: alpha1=1.0, alpha2= 0.0

!-----Create material constant matrices


  l_tensor = 0.d0
  m_tensor = 0.d0
  c2 = 1.d0 - cm
  gm = em/(2.d0*(1.d0+prm))
  gp = ep/(2.d0*(1.d0+prp))

  m_tensor(1,1,1) = 1.d0 / em
  m_tensor(1,2,1) = - prm / em
  m_tensor(1,3,1) = - prm / em
  m_tensor(2,1,1) = m_tensor(1,2,1) 
  m_tensor(2,2,1) = 1.d0 / em
  m_tensor(2,3,1) = - prm / em
  m_tensor(3,1,1) = m_tensor(1,3,1)
  m_tensor(3,2,1) = m_tensor(2,3,1)
  m_tensor(3,3,1) = 1.d0 /em
  m_tensor(4,4,1) = 1.d0/gm
  m_tensor(5,5,1) = 1.d0/gm
  m_tensor(6,6,1) = 1.d0/gm

  m_tensor(1,1,2) = 1.d0 / ep
  m_tensor(1,2,2) = - prp / ep
  m_tensor(1,3,2) = - prp / ep
  m_tensor(2,1,2) = m_tensor(1,2,2) 
  m_tensor(2,2,2) = 1.d0 / ep
  m_tensor(2,3,2) = - prp / ep
  m_tensor(3,1,2) = m_tensor(1,3,2)
  m_tensor(3,2,2) = m_tensor(2,3,2)
  m_tensor(3,3,2) = 1.d0 /ep
  m_tensor(4,4,2) = 1.d0/gp
  m_tensor(5,5,2) = 1.d0/gp
  m_tensor(6,6,2) = 1.d0/gp

  CALL invert2(m_tensor(:,:,1), l_tensor(:,:,1), 6)
  CALL invert2(m_tensor(:,:,2), l_tensor(:,:,2), 6)


  CALL compositestiffnes(l_tensor(1:6,1:6,1), l_tensor(1:6,1:6,2), m_tensor(1:6,1:6,1), m_tensor(1:6,1:6,2), cm, c2, &
       alpha1, alpha2, l_bar, m_bar, lo)


  e_homog = 1.d0/m_bar(1,1)
  pr_homog = -m_bar(1,2)*1./m_bar(1,1) ! nu = -M(1,2)*E

!  IF(myid.EQ.0)THEN
  print*,'Homogenized Stiffness =', e_homog
  print*,'Homogenized Poisson'//"'"//'s ratio',pr_homog
!  ENDIF

  cd_fastest = max( cd_fastest, &
             sqrt(e_homog*(1.d0-pr_homog)/density/(1.d0+pr_homog)/(1.d0-2.d0*pr_homog )) )

  RETURN
END SUBROUTINE homogenizedmat

SUBROUTINE compositestiffnes(Lm, L2, Mm, M2, cm, c2, alpha1, alpha2, L_bar, M_bar, Lo)

  IMPLICIT NONE


  REAL*8 :: alpha1, alpha2, cm, c2

  REAL*8, DIMENSION(1:6,1:6) :: lm, l2, mm, m2
  REAL*8, DIMENSION(1:6,1:6) :: l_bar, m_bar

  REAL*8, DIMENSION(1:6,1:6) :: dident = reshape( &
       (/1.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
         0.d0, 1.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
         0.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, &
         0.d0, 0.d0, 0.d0, 1.d0, 0.d0, 0.d0, &
         0.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, &
         0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 1.d0 /),(/6,6/) )

  REAL*8, DIMENSION(1:6,1:6) :: lo, mo 
  REAL*8, DIMENSION(1:6,1:6) :: p
  REAL*8, DIMENSION(1:6,1:6) :: s, sinv
  REAL*8, DIMENSION(1:6,1:6) :: am, a2, a, b, c, lst, a3
  
  INTEGER :: i,j

 

!  Comparison medium,  Lo = (a1*L1 + a2*L2)

  lo = alpha1*lm + alpha2*l2

  CALL invert2(lo, mo, 6)

  
!  /* Type of inclusion */

! Returns Eshelby tensor S

  CALL eshelby(0, mo, lo, p, s)

  
!          -1
! Returns S

  CALL invert2(s, sinv, 6)

!  nas_AB (Lo,SI,A,6,6,6); nas_AB (A,B,Lst,6,6,6);
! 
! I - S

 b = dident - s

! Lo * ( I - S)

 a = matmul(lo,sinv)
!
! Hill's constrant tensor
!
!  *    -1
! L  = S * Lo * ( I - S)

 lst = matmul(a, b)

!!$ PRINT*,'L*'
!!$
!!$ DO i = 1, 6
!!$    PRINT*,Lst(i,:)
!!$ ENDDO


! /* Compute L(bar) */

 a = lst + lm
 b = lst + l2

 CALL invert2(a,am,6)
 CALL invert2(b,a2,6)
 c = cm*am + c2*a2
 CALL invert2(c,a,6)
 l_bar = a - lst

 a = lst + l_bar

! compute Am

 a2 = lst + lm 
 CALL invert2(a2,c,6)
 am = matmul(c,a)
 
! compute A2

 a2 = lst + l2
 CALL invert2(a2,c,6)
 a2 = matmul(c,a)

! Check concentration factors

 DO i = 1, 6
    DO j = 1, 6
       c(i,j) = cm*am(i,j) + c2*a2(i,j)
       IF ( i .EQ. j .AND. (c(i,j) <= 0.9999 .OR. c(i,j) >=  1.0001))THEN
          print*,'a) Incorrect Concentration tensors Am, A2'
          stop
       ENDIF
       IF ( i.NE.j .AND.( c(i,j) >= 1.e-10 .OR. c(i,j) <= -1.e-10))THEN
          print*,'b) Incorrect Concentration tensors Am, A2'
          stop
       ENDIF
    ENDDO
 ENDDO

 CALL invert2(l_bar, m_bar, 6)

END SUBROUTINE compositestiffnes

SUBROUTINE eshelby(Shape, Mo, Lo, P, S)

  INTEGER :: shape
  REAL*8, DIMENSION(1:6,1:6) :: mo, lo
  REAL*8, DIMENSION(1:6,1:6) :: p

  REAL*8, DIMENSION(1:6,1:6) :: s
  

  REAL*8 :: e, g, nu, k, si, th


  s = 0.d0
  p = 0.d0

  ! Type of inclusion

  SELECT CASE (shape)

  CASE (0) ! /* Spherical inclusion in Isotropic medium */
     e = 1.d0/mo(1,1)
     g = 1.d0/mo(4,4)
     nu = -mo(1,2)*e
     k = e/(3.d0*(1.d0 - 2.d0*nu))
     si = (1.d0 + nu)/(3.d0*(1.d0 - nu)) 
     th = 2.d0*(4.d0 - 5.d0*nu)/(15.d0*(1.d0 - nu))
     
     p(1,1) = 1.d0/3.d0*(si/3.d0/k + th/g)
     p(1,2) = 1.d0/3.d0*(si/3.d0/k - th/2.d0/g)
     p(1,3) = 1.d0/3.d0*(si/3.d0/k - th/2.d0/g)
     p(2,1) = 1.d0/3.d0*(si/3.d0/k - th/2.d0/g)
     p(2,2) = 1.d0/3.d0*(si/3.d0/k + th/g)
     p(2,3) = 1.d0/3.d0*(si/3.d0/k - th/2.d0/g)
     p(3,1) = 1.d0/3.d0*(si/3.d0/k - th/2.d0/g)
     p(3,2) = 1.d0/3.d0*(si/3.d0/k - th/2.d0/g)
     p(3,3) = 1.d0/3.d0*(si/3.d0/k + th/g)
     p(4,4) = th/g
     p(5,5) = th/g
     p(6,6) = th/g
     s = matmul(lo,p)

  CASE default
     print*, "Non existing shape of inclusion"
     stop
  END SELECT

END SUBROUTINE eshelby

SUBROUTINE invert2( Inv_in, a, nrow)
                                                                                                          
! *--------------------------------------------------------------------*
! |                                                                    |
! |   Inverts a matrix by gauss jordan elimination and computes the    |
! |   the determinant. The matrix may be symmetric or nonsymmetric.    |
! |                                                                    |
! |    <a>    : a square, double precision matrix                      |
! |    <nrow> : dimensioned number of rows for 'a'. also the actual    |
! |             sizes for inversion computations                       |
! |    <det>  : determinant of the matrix ( zero if matrix is singular)|
! |                                                                    |
! |     NOTE  : the matrix is overwritten by the inverse.              |
! |                                                                    |
! *--------------------------------------------------------------------*
  IMPLICIT NONE
                                                                                                          
!  Argument variables
                                                                                                          
  INTEGER nrow
  REAL*8, DIMENSION(1:nrow,1:nrow) :: a, inv_in
  REAL*8 :: det
                                                                                 
!  Local variables
                                                                                 
  INTEGER ncol, k, j, i
                                                                                                          
  a = inv_in
                                                                                 
!  Invert matrix
                                                                                 
  det  = 1.d0
  ncol = nrow
                                                                                 
  DO  k = 1, nrow
                                                                                 
     IF( a(k,k) .EQ. 0.d0 ) THEN ! Check for singular matrix
        det = 0.d0
        print*,' Matrix is singular for inversion'
        stop
     END IF
                                                                                 
     det = det * a(k, k)
                                                                                 
     DO j = 1, ncol
        IF( j .NE. k ) a(k,j) = a(k,j) / a(k,k)
     END DO
                                                                                 
     a(k,k) = 1.d0 / a(k,k)
                                                                                 
     DO  i = 1, nrow
        IF( i .EQ. k) cycle
        DO j = 1, ncol
           IF( j .NE. k) a(i,j) = a(i,j) - a(i,k) * a(k,j)
        END DO
        a(i,k) = -a(i,k) * a(k,k)
     END DO
                                                                                 
  END DO
                                                                                                          
                                                                                 
  RETURN
END SUBROUTINE invert2