calculation of shear modulus for cubic crystals can be simplifieid, tests added for different crystal symmetries

examples/grid/material.yaml

@@ -9,7 +9,7 @@ phase:
     lattice: cF
     mechanical:
       output: [F, P, F_e, F_p, L_p, O]
-      elastic: {type: Hooke, C_11: 106.75e+9, C_12: 60.41e+9, C_44: 28.34e+9, modulus_type: 'Voigt'}
+      elastic: {type: Hooke, C_11: 106.75e+9, C_12: 60.41e+9, C_44: 28.34e+9}
       plastic:
         type: phenopowerlaw
         N_sl: [12]

src/homogenization_mechanical_RGC.f90

@@ -9,6 +9,7 @@
 submodule(homogenization:mechanical) RGC
   use rotations
   use lattice
+  use phase
 
   type :: tParameters
     integer, dimension(:), allocatable :: &
@@ -652,9 +653,9 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
 
     real(pReal), dimension(6,6) :: C
 
-
     C = phase_homogenizedC66(material_phaseID(co,ce),material_phaseEntry(co,ce))                    ! damage not included!
-    equivalentMu = lattice_equivalent_mu(C,'voigt')
+
+    equivalentMu = lattice_isotropic_mu(C,'voigt',phase_lattice_structure(co,ce))
 
   end function equivalentMu
 

src/lattice.f90

@@ -376,8 +376,8 @@ module lattice
 
   public :: &
     lattice_init, &
-    lattice_equivalent_nu, &
+    lattice_isotropic_nu, &
-    lattice_equivalent_mu, &
+    lattice_isotropic_mu, &
     lattice_symmetrize_33, &
     lattice_symmetrize_C66, &
     lattice_SchmidMatrix_slip, &
@@ -2149,10 +2149,11 @@ end function getlabels
 !> @brief Equivalent Poisson's ratio (ν)
 !> @details https://doi.org/10.1143/JPSJ.20.635
 !--------------------------------------------------------------------------------------------------
-pure function lattice_equivalent_nu(C,assumption) result(nu)
+pure function lattice_isotropic_nu(C,assumption,lattice) result(nu)
 
   real(pReal), dimension(6,6), intent(in) :: C                                                      !< Stiffness tensor (Voigt notation)
   character(len=5),            intent(in) :: assumption                                             !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
+  character(len=2),            intent(in) :: lattice
   real(pReal) :: nu
 
   real(pReal) :: K, mu
@@ -2172,20 +2173,22 @@ pure function lattice_equivalent_nu(C,assumption) result(nu)
     error stop 'invalid assumption'
   end if
 
-  mu = lattice_equivalent_mu(C,assumption)
+  mu = lattice_isotropic_mu(C,assumption,lattice)
   nu = (1.5_pReal*K-mu)/(3.0_pReal*K+mu)
 
-end function lattice_equivalent_nu
+end function lattice_isotropic_nu
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Equivalent shear modulus (μ)
 !> @details https://doi.org/10.1143/JPSJ.20.635
+!> @details Nonlinear Mechanics of Crystals 10.1007/978-94-007-0350-6, pp 563
 !--------------------------------------------------------------------------------------------------
-pure function lattice_equivalent_mu(C,assumption) result(mu)
+pure function lattice_isotropic_mu(C,assumption,lattice) result(mu)
 
   real(pReal), dimension(6,6), intent(in) :: C                                                      !< Stiffness tensor (Voigt notation)
   character(len=5),            intent(in) :: assumption                                             !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
+  character(len=2),            intent(in) :: lattice
   real(pReal) :: mu
 
   logical                     :: error
@@ -2193,18 +2196,32 @@ pure function lattice_equivalent_mu(C,assumption) result(mu)
 
 
   if     (IO_lc(assumption) == 'voigt') then
-    mu = (1.0_pReal*(C(1,1)+C(2,2)+C(3,3)) -1.0_pReal*(C(1,2)+C(2,3)+C(1,3)) +3.0_pReal*(C(4,4)+C(5,5)+C(6,6))) &
+      select case(lattice)
-       / 15.0_pReal
+        case('cF','cI')
+          mu = ( C(1,1) - C(1,2) + C(4,4)*3.0_pReal) / 5.0_pReal
+        case default
+          mu = (  C(1,1)+C(2,2)+C(3,3) &
+                -(C(1,2)+C(2,3)+C(1,3)) &
+                +(C(4,4)+C(5,5)+C(6,6)) * 3.0_pReal &
+               ) / 15.0_pReal
+      end select
+
   elseif (IO_lc(assumption) == 'reuss') then
-    call math_invert(S,error,C)
+      select case(lattice)
-    if (error) error stop 'matrix inversion failed'
+        case('cF','cI')
-    mu = 15.0_pReal &
+          mu = 1.0_pReal &
-       / (4.0_pReal*(S(1,1)+S(2,2)+S(3,3)) -4.0_pReal*(S(1,2)+S(2,3)+S(1,3)) +3.0_pReal*(S(4,4)+S(5,5)+S(6,6)))
+               / ((4.0_pReal/(5.0_pReal * (C(1,1)-C(1,2)))) + (3.0_pReal/(5.0_pReal*C(4,4))))
+        case default
+          call math_invert(S,error,C)
+          if (error) error stop 'matrix inversion failed'
+          mu = 15.0_pReal &
+              / (4.0_pReal*(S(1,1)+S(2,2)+S(3,3)) -4.0_pReal*(S(1,2)+S(2,3)+S(1,3)) +3.0_pReal*(S(4,4)+S(5,5)+S(6,6)))
+      end select
   else
     error stop 'invalid assumption'
   end if
 
-end function lattice_equivalent_mu
+end function lattice_isotropic_mu
 
 
 !--------------------------------------------------------------------------------------------------
@@ -2270,16 +2287,42 @@ subroutine selfTest
 
   call random_number(C)
   C(1,1) = C(1,1) + C(1,2) + 0.1_pReal
+  C(1,3) = C(1,2)
+  C(3,3) = C(1,1) 
   C(4,4) = 0.5_pReal * (C(1,1) - C(1,2))
-  C = lattice_symmetrize_C66(C,'cI')
+  C(6,6) = C(4,4)
-  if (dNeq(C(4,4),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt'
+
-  if (dNeq(C(4,4),lattice_equivalent_mu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss'
+  C_cI = lattice_symmetrize_C66(C,'cI')
+  if (dNeq(C_cI(4,4),lattice_isotropic_mu(C_cI,'voigt','cI'),1.0e-12_pReal)) error stop 'isotropic_mu/cI/voigt'
+  if (dNeq(C_cI(4,4),lattice_isotropic_mu(C_cI,'reuss','cI'),1.0e-12_pReal)) error stop 'isotropic_mu/cI/reuss'
+
+  lambda = C_cI(1,2)
+  if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_cI,'voigt','cI')), &
+          lattice_isotropic_nu(C_cI,'voigt','cI'),1.0e-12_pReal)) error stop 'isotropic_nu/cI/voigt'
+  if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_cI,'reuss','cI')), &
+          lattice_isotropic_nu(C_cI,'reuss','cI'),1.0e-12_pReal)) error stop 'isotropic_nu/cI/reuss'
+
+
+  C_hP = lattice_symmetrize_C66(C,'hP')
+  if (dNeq(C(4,4),lattice_isotropic_mu(C_hP,'voigt','hP'),1.0e-12_pReal)) error stop 'isotropic_mu/hP/voigt'
+  if (dNeq(C(4,4),lattice_isotropic_mu(C_hP,'reuss','hP'),1.0e-12_pReal)) error stop 'isotropic_mu/hP/reuss'
+ 
+  lambda = C_hP(1,2)
+  if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_hP,'voigt','hP')), &
+          lattice_isotropic_nu(C_hP,'voigt','hP'),1.0e-12_pReal)) error stop 'isotropic_nu/hP/voigt'
+  if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_hP,'reuss','hP')), &
+          lattice_isotropic_nu(C_hP,'reuss','hP'),1.0e-12_pReal)) error stop 'isotropic_nu/hP/reuss'
+
+  C_tI = lattice_symmetrize_C66(C,'tI')
+  if (dNeq(C(6,6),lattice_isotropic_mu(C_tI,'voigt','tI'),1.0e-12_pReal)) error stop 'isotropic_mu/tI/voigt'
+  if (dNeq(C(6,6),lattice_isotropic_mu(C_tI,'reuss','tI'),1.0e-12_pReal)) error stop 'isotropic_mu/tI/reuss'
+
+  lambda = C_tI(1,2)
+  if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_tI,'voigt','tI')), &
+          lattice_isotropic_nu(C_tI,'voigt','tI'),1.0e-12_pReal)) error stop 'isotropic_nu/tI/voigt'
+  if (dNeq(lambda*0.5_pReal/(lambda+lattice_isotropic_mu(C_tI,'reuss','tI')), &
+          lattice_isotropic_nu(C_tI,'reuss','tI'),1.0e-12_pReal)) error stop 'isotropic_nu/tI/reuss'
 
-  lambda = C(1,2)
-  if (dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'voigt')), &
-          lattice_equivalent_nu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_nu/voigt'
-  if (dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'reuss')), &
-          lattice_equivalent_nu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_nu/reuss'
 
 end subroutine selfTest
 

src/phase.f90

@@ -347,6 +347,7 @@ module phase
     phase_K_T, &
     phase_mu_phi, &
     phase_mu_T, &
+    phase_lattice_structure, &
     phase_results, &
     phase_allocateState, &
     phase_forward, &
@@ -434,6 +435,20 @@ subroutine phase_init
 end subroutine phase_init
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Get lattice structure for a given phase
+!--------------------------------------------------------------------------------------------------
+function phase_lattice_structure(co,ce) result(lattice)
+
+  integer, intent(in) :: co, ce
+
+  character(len=2)    :: lattice
+
+  lattice = phase_lattice(material_phaseID(co,ce))
+
+end function phase_lattice_structure
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief Allocate the components of the state structure for a given phase
 !--------------------------------------------------------------------------------------------------

src/phase_mechanical.f90

@@ -169,14 +169,16 @@ submodule(phase) mechanical
       integer,     intent(in) :: ph, en
     end function elastic_C66
 
-    pure module function elastic_mu(ph,en) result(mu)
+    pure module function elastic_mu(ph,en,isotropic_bound) result(mu)
       real(pReal) :: mu
       integer, intent(in) :: ph, en
+      character(len=5), intent(in) :: isotropic_bound
     end function elastic_mu
 
-    pure module function elastic_nu(ph,en) result(nu)
+    pure module function elastic_nu(ph,en,isotropic_bound) result(nu)
       real(pReal) :: nu
       integer, intent(in) :: ph, en
+      character(len=5), intent(in) :: isotropic_bound
     end function elastic_nu
 
   end interface

src/phase_mechanical_elastic.f90

@@ -8,7 +8,6 @@ submodule(phase:mechanical) elastic
       C_33, &
       C_44, &
       C_66
-    character(len=pStringLen) :: modulus_type
   end type tParameters
 
   type(tParameters), allocatable, dimension(:) :: param
@@ -58,8 +57,6 @@ module subroutine elastic_init(phases)
       if (phase_lattice(ph) == 'tI') &
         prm%C_66 = polynomial(elastic%asDict(),'C_66','T')
 
-      prm%modulus_type=elastic%get_asString('modulus_type',defaultVal='Voigt')
-
     end associate
   end do
 
@@ -105,17 +102,18 @@ end function elastic_C66
 !--------------------------------------------------------------------------------------------------
 !> @brief return shear modulus
 !--------------------------------------------------------------------------------------------------
-pure module function elastic_mu(ph,en) result(mu)
+pure module function elastic_mu(ph,en,isotropic_bound) result(mu)
 
   integer, intent(in) :: &
     ph, &
     en
+  character(len=5), intent(in) :: isotropic_bound
   real(pReal) :: &
     mu
 
   associate(prm => param(ph))
 
-  mu = lattice_equivalent_mu(elastic_C66(ph,en),prm%modulus_type)
+    mu = lattice_isotropic_mu(elastic_C66(ph,en),isotropic_bound,phase_lattice(ph))
 
   end associate
 
@@ -125,17 +123,18 @@ end function elastic_mu
 !--------------------------------------------------------------------------------------------------
 !> @brief return Poisson ratio
 !--------------------------------------------------------------------------------------------------
-pure module function elastic_nu(ph,en) result(nu)
+pure module function elastic_nu(ph,en,isotropic_bound) result(nu)
 
   integer, intent(in) :: &
     ph, &
     en
+  character(len=5), intent(in) :: isotropic_bound
   real(pReal) :: &
     nu
 
   associate(prm => param(ph))
 
-  nu = lattice_equivalent_nu(elastic_C66(ph,en),prm%modulus_type)
+    nu = lattice_isotropic_nu(elastic_C66(ph,en),isotropic_bound,phase_lattice(ph))
 
   end associate
 

src/phase_mechanical_plastic_dislotungsten.f90

@@ -35,6 +35,8 @@ submodule(phase:plastic) dislotungsten
       P_nS_neg
     integer :: &
       sum_N_sl                                                                                      !< total number of active slip system
+    character(len=5)                                     :: &
+      isotropic_bound
     character(len=pStringLen), allocatable, dimension(:) :: &
       output
     logical :: &
@@ -130,6 +132,8 @@ module function plastic_dislotungsten_init() result(myPlasticity)
 #else
     prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray)
 #endif
+    
+    prm%isotropic_bound = pl%get_asString('isotropic_bound',defaultVal='Voigt')
 
 !--------------------------------------------------------------------------------------------------
 ! slip related parameters
@@ -333,7 +337,7 @@ module function dislotungsten_dotState(Mp,ph,en) result(dotState)
             dot_rho_dip => dotState(indexDotState(ph)%rho_dip(1):indexDotState(ph)%rho_dip(2)), &
             dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)))
 
-    mu = elastic_mu(ph,en)
+    mu = elastic_mu(ph,en,prm%isotropic_bound)
     T = thermal_T(ph,en)
 
     call kinetics(Mp,T,ph,en,&
@@ -384,7 +388,7 @@ module subroutine dislotungsten_dependentState(ph,en)
 
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
-    dst%tau_pass(:,en) = elastic_mu(ph,en)*prm%b_sl &
+    dst%tau_pass(:,en) = elastic_mu(ph,en,prm%isotropic_bound)*prm%b_sl &
                        * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
 
     Lambda_sl_inv = 1.0_pReal/prm%D &

src/phase_mechanical_plastic_dislotwin.f90

@@ -74,6 +74,8 @@ submodule(phase:plastic) dislotwin
       fcc_twinNucleationSlipPair                                                                    ! ToDo: Better name? Is also use for trans
     character(len=:),          allocatable                 :: &
       lattice_tr
+    character(len=5)                                       :: &
+      isotropic_bound
     character(len=pStringLen), allocatable, dimension(:) :: &
       output
     logical :: &
@@ -186,6 +188,8 @@ module function plastic_dislotwin_init() result(myPlasticity)
     prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray)
 #endif
 
+   prm%isotropic_bound = pl%get_asString('isotropic_bound',defaultVal='Voigt')
+
 !--------------------------------------------------------------------------------------------------
 ! slip related parameters
     N_sl         = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
@@ -644,8 +648,8 @@ module function dislotwin_dotState(Mp,ph,en) result(dotState)
             dot_f_tw => dotState(indexDotState(ph)%f_tw(1):indexDotState(ph)%f_tw(2)), &
             dot_f_tr => dotState(indexDotState(ph)%f_tr(1):indexDotState(ph)%f_tr(2)))
 
-    mu = elastic_mu(ph,en)
+    mu = elastic_mu(ph,en,prm%isotropic_bound)
-    nu = elastic_nu(ph,en)
+    nu = elastic_nu(ph,en,prm%isotropic_bound)
     T = thermal_T(ph,en)
 
     f_matrix = 1.0_pReal &
@@ -732,7 +736,7 @@ module subroutine dislotwin_dependentState(ph,en)
 
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
-    mu = elastic_mu(ph,en)
+    mu = elastic_mu(ph,en,prm%isotropic_bound)
     sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,en))
     sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en))
 
@@ -930,8 +934,8 @@ pure subroutine kinetics_tw(Mp,T,abs_dot_gamma_sl,ph,en,&
 
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
-    mu = elastic_mu(ph,en)
+    mu = elastic_mu(ph,en,prm%isotropic_bound)
-    nu = elastic_nu(ph,en)
+    nu = elastic_nu(ph,en,prm%isotropic_bound)
     Gamma_sf = prm%Gamma_sf%at(T)
 
     tau_hat = 3.0_pReal*prm%b_tw(1)*mu/prm%L_tw &
@@ -1006,8 +1010,8 @@ pure subroutine kinetics_tr(Mp,T,abs_dot_gamma_sl,ph,en,&
 
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
-    mu = elastic_mu(ph,en)
+    mu = elastic_mu(ph,en,prm%isotropic_bound)
-    nu = elastic_nu(ph,en)
+    nu = elastic_nu(ph,en,prm%isotropic_bound)
     Gamma_sf = prm%Gamma_sf%at(T)
 
     tau_hat = 3.0_pReal*prm%b_tr(1)*mu/prm%L_tr &

src/phase_mechanical_plastic_nonlocal.f90

@@ -115,6 +115,8 @@ submodule(phase:plastic) nonlocal
       sum_N_sl = 0
     integer,     dimension(:),     allocatable :: &
       colinearSystem                                                                                !< colinear system to the active slip system (only valid for fcc!)
+    character(len=5)                           :: &
+      isotropic_bound          
     character(len=pStringLen), dimension(:), allocatable :: &
       output
     logical :: &
@@ -241,6 +243,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
     prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray)
 #endif
 
+    prm%isotropic_bound = pl%get_asString('isotropic_bound',defaultVal='Voigt')
     prm%atol_rho = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
 
     ini%N_sl     = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
@@ -609,8 +612,8 @@ module subroutine nonlocal_dependentState(ph, en)
 
   associate(prm => param(ph),dst => dependentState(ph), stt => state(ph))
 
-  mu = elastic_mu(ph,en)
+  mu = elastic_mu(ph,en,prm%isotropic_bound)
-  nu = elastic_nu(ph,en)
+  nu = elastic_nu(ph,en,prm%isotropic_bound)
   rho = getRho(ph,en)
 
   stt%rho_forest(:,en) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) &
@@ -880,8 +883,8 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
 
   associate(prm => param(ph),dst => dependentState(ph),del => deltaState(ph))
 
-  mu = elastic_mu(ph,en)
+  mu = elastic_mu(ph,en,prm%isotropic_bound)
-  nu = elastic_nu(ph,en)
+  nu = elastic_nu(ph,en,prm%isotropic_bound)
 
   !*** shortcut to state variables
   forall (s = 1:prm%sum_N_sl, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)
@@ -994,8 +997,8 @@ module subroutine nonlocal_dotState(Mp,timestep, &
 
   associate(prm => param(ph), dst => dependentState(ph), dot => dotState(ph), stt => state(ph))
 
-  mu = elastic_mu(ph,en)
+  mu = elastic_mu(ph,en,prm%isotropic_bound)
-  nu = elastic_nu(ph,en)
+  nu = elastic_nu(ph,en,prm%isotropic_bound)
   Temperature = thermal_T(ph,en)
 
   tau = 0.0_pReal