From dc406a01c02b8a6c7ec0a24330022c9428f1062c Mon Sep 17 00:00:00 2001 From: Luv Sharma Date: Wed, 10 Sep 2014 08:37:12 +0000 Subject: [PATCH] added funtions to get averged properties at integration points. --- code/crystallite.f90 | 29 ++++++ code/homogenization.f90 | 215 ++++++++++++++++++++++++++++++++++++++++ code/lattice.f90 | 41 ++++++++ 3 files changed, 285 insertions(+) diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 76c85d229..89a907781 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -115,6 +115,9 @@ module crystallite crystallite_init, & crystallite_stressAndItsTangent, & crystallite_orientations, & +#ifdef NEWSTATE + crystallite_push33ToRef, & +#endif crystallite_postResults private :: & crystallite_integrateStateFPI, & @@ -3413,6 +3416,32 @@ logical function crystallite_stateJump(g,i,e) end function crystallite_stateJump +#ifdef NEWSTATE +!-------------------------------------------------------------------------------------------------- +!> @brief Map 2nd order tensor to reference config +!-------------------------------------------------------------------------------------------------- +function crystallite_push33ToRef(g,i,e, tensor33) + use prec, only: & + pInt, & + pReal + use math, only: & + math_inv33 + + implicit none + real(pReal), dimension(3,3) :: crystallite_push33ToRef + real(pReal), dimension(3,3), intent(in) :: tensor33 + real(pReal), dimension(3,3) :: invFp + integer(pInt), intent(in):: & + e, & ! element index + i, & ! integration point index + g ! grain index + + invFp = math_inv33(crystallite_Fp(1:3,1:3,g,i,e)) + crystallite_push33ToRef = matmul(invFp,matmul(tensor33,transpose(invFp))) + +end function crystallite_push33ToRef + +#endif !-------------------------------------------------------------------------------------------------- !> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> intermediate acceleration of the Newton-Raphson correction diff --git a/code/homogenization.f90 b/code/homogenization.f90 index 7eea45481..7e5352ff7 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -65,6 +65,11 @@ module homogenization #ifdef NEWSTATE field_getDAMAGE, & field_putDAMAGE, & + field_getDamageMobility, & + field_getDamageDiffusion33, & + field_getThermalConductivity33, & + field_getMassDensity, & + field_getSpecificHeat, & #endif materialpoint_postResults private :: & @@ -943,6 +948,216 @@ end function homogenization_averageHeat #ifdef NEWSTATE !-------------------------------------------------------------------------------------------------- +!> @brief Returns average specific heat at each integration point +!-------------------------------------------------------------------------------------------------- +function field_getSpecificHeat(ip,el) + use mesh, only: & + mesh_element + use lattice, only: & + lattice_specificHeat + use material, only: & + material_phase, & + material_homog, & + field_thermal_type, & + FIELD_THERMAL_ADIABATIC_ID, & + FIELD_THERMAL_CONDUCTION_ID, & + homogenization_Ngrains + + implicit none + real(pReal) :: field_getSpecificHeat + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + integer(pInt) :: & + Ngrains, ipc + + field_getSpecificHeat =0.0_pReal + + select case(field_thermal_type(material_homog(ip,el))) + + case (FIELD_THERMAL_ADIABATIC_ID) + field_getSpecificHeat = 0.0_pReal + + case (FIELD_THERMAL_CONDUCTION_ID) + do ipc = 1, homogenization_Ngrains(mesh_element(3,el)) + field_getSpecificHeat = field_getSpecificHeat + lattice_specificHeat(material_phase(ipc,ip,el)) + enddo + + end select + + field_getSpecificHeat = field_getSpecificHeat /homogenization_Ngrains(mesh_element(3,el)) + +end function field_getSpecificHeat + +!-------------------------------------------------------------------------------------------------- +!> @brief Returns average mass density at each integration point +!-------------------------------------------------------------------------------------------------- +function field_getMassDensity(ip,el) + use mesh, only: & + mesh_element + use lattice, only: & + lattice_massDensity + use material, only: & + material_phase, & + material_homog, & + field_thermal_type, & + FIELD_THERMAL_ADIABATIC_ID, & + FIELD_THERMAL_CONDUCTION_ID, & + homogenization_Ngrains + + + implicit none + real(pReal) :: field_getMassDensity + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + integer(pInt) :: & + Ngrains, ipc + + field_getMassDensity =0.0_pReal + + select case(field_thermal_type(material_homog(ip,el))) + + case (FIELD_THERMAL_ADIABATIC_ID) + field_getMassDensity = 0.0_pReal + + case (FIELD_THERMAL_CONDUCTION_ID) + do ipc = 1, homogenization_Ngrains(mesh_element(3,el)) + field_getMassDensity = field_getMassDensity + lattice_massDensity(material_phase(ipc,ip,el)) + enddo + + end select + + field_getMassDensity = field_getMassDensity /homogenization_Ngrains(mesh_element(3,el)) + +end function field_getMassDensity +!------------------------------------------------------------------------------------------- +!> @brief Returns average conductivity tensor for thermal field at each integration point +!------------------------------------------------------------------------------------------- +function field_getThermalConductivity33(ip,el) + use mesh, only: & + mesh_element + use lattice, only: & + lattice_thermalConductivity33 + use material, only: & + material_phase, & + material_homog, & + field_thermal_type, & + FIELD_THERMAL_ADIABATIC_ID, & + FIELD_THERMAL_CONDUCTION_ID, & + homogenization_Ngrains + use crystallite, only: & + crystallite_push33ToRef + + + implicit none + real(pReal), dimension(3,3) :: field_getThermalConductivity33 + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + integer(pInt) :: & + Ngrains, ipc + + field_getThermalConductivity33 =0.0_pReal + + select case(field_thermal_type(material_homog(ip,el))) + + case (FIELD_THERMAL_ADIABATIC_ID) + field_getThermalConductivity33 = 0.0_pReal + + case (FIELD_THERMAL_CONDUCTION_ID) + do ipc = 1, homogenization_Ngrains(mesh_element(3,el)) + field_getThermalConductivity33 = field_getThermalConductivity33 + & + crystallite_push33ToRef(ipc,ip,el,lattice_thermalConductivity33(:,:,material_phase(ipc,ip,el))) + enddo + + end select + + field_getThermalConductivity33 = field_getThermalConductivity33 /homogenization_Ngrains(mesh_element(3,el)) + +end function field_getThermalConductivity33 +!-------------------------------------------------------------------------------------------------- +!> @brief Returns average diffusion tensor for damage field at each integration point +!-------------------------------------------------------------------------------------------------- +function field_getDamageDiffusion33(ip,el) + use mesh, only: & + mesh_element + use lattice, only: & + lattice_DamageDiffusion33 + use material, only: & + material_phase, & + material_homog, & + field_damage_type, & + FIELD_DAMAGE_LOCAL_ID, & + FIELD_DAMAGE_NONLOCAL_ID, & + homogenization_Ngrains + + implicit none + real(pReal), dimension(3,3) :: field_getDamageDiffusion33 + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + integer(pInt) :: & + Ngrains, ipc + + field_getDamageDiffusion33 =0.0_pReal + + select case(field_damage_type(material_homog(ip,el))) + + case (FIELD_DAMAGE_LOCAL_ID) + field_getDamageDiffusion33 = 0.0_pReal + + case (FIELD_DAMAGE_NONLOCAL_ID) + do ipc = 1, homogenization_Ngrains(mesh_element(3,el)) + field_getDamageDiffusion33 = field_getDamageDiffusion33 + lattice_DamageDiffusion33(:,:,material_phase(ipc,ip,el)) + enddo + + end select + + field_getDamageDiffusion33 = field_getDamageDiffusion33 /homogenization_Ngrains(mesh_element(3,el)) + +end function field_getDamageDiffusion33 +!-------------------------------------------------------------------------------------------------- +!> @brief Returns average mobility for damage field at each integration point +!-------------------------------------------------------------------------------------------------- +real(pReal) function field_getDamageMobility(ip,el) + use mesh, only: & + mesh_element + use lattice, only: & + lattice_damageMobility + use material, only: & + material_phase, & + material_homog, & + field_damage_type, & + FIELD_DAMAGE_LOCAL_ID, & + FIELD_DAMAGE_NONLOCAL_ID, & + homogenization_Ngrains + + implicit none + integer(pInt), intent(in) :: & + ip, & !< integration point number + el !< element number + integer(pInt) :: & + Ngrains, ipc + + field_getDamageMobility =0.0_pReal + + select case(field_damage_type(material_homog(ip,el))) + + case (FIELD_DAMAGE_LOCAL_ID) + field_getDamageMobility = 0.0_pReal + + case (FIELD_DAMAGE_NONLOCAL_ID) + do ipc = 1, homogenization_Ngrains(mesh_element(3,el)) + field_getDamageMobility = field_getDamageMobility + lattice_DamageMobility(material_phase(ipc,ip,el)) + enddo + + end select + + field_getDamageMobility = field_getDamageMobility /homogenization_Ngrains(mesh_element(3,el)) + +end function field_getDamageMobility +!-------------------------------------------------------------------------------------------------- !> @brief ToDo !-------------------------------------------------------------------------------------------------- real(pReal) function field_getDAMAGE(ip,el) diff --git a/code/lattice.f90 b/code/lattice.f90 index 0244aa2c9..a49009934 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -708,8 +708,16 @@ module lattice real(pReal), dimension(:,:,:), allocatable, public, protected :: & lattice_thermalConductivity33, & lattice_thermalExpansion33, & +#ifdef NEWSTATE + lattice_damageDiffusion33, & +#endif lattice_surfaceEnergy33 real(pReal), dimension(:), allocatable, public, protected :: & +#ifdef NEWSTATE + lattice_damageMobility, & + lattice_massDensity, & + lattice_specificHeat, & +#endif lattice_referenceTemperature enum, bind(c) enumerator :: LATTICE_undefined_ID, & @@ -963,6 +971,12 @@ subroutine lattice_init allocate(lattice_C3333(3,3,3,3,Nphases), source=0.0_pReal) allocate(lattice_thermalConductivity33(3,3,Nphases), source=0.0_pReal) allocate(lattice_thermalExpansion33 (3,3,Nphases), source=0.0_pReal) +#ifdef NEWSTATE + allocate(lattice_damageDiffusion33 (3,3,Nphases), source=0.0_pReal) + allocate(lattice_damageMobility ( Nphases), source=0.0_pReal) + allocate(lattice_massDensity ( Nphases), source=0.0_pReal) + allocate(lattice_specificHeat ( Nphases), source=0.0_pReal) +#endif allocate(lattice_surfaceEnergy33 (3,3,Nphases), source=0.0_pReal) allocate(lattice_referenceTemperature (Nphases), source=0.0_pReal) @@ -1088,6 +1102,33 @@ subroutine lattice_init lattice_surfaceEnergy33(3,3,section) = IO_floatValue(line,positions,2_pInt) case ('reference_temperature') lattice_referenceTemperature(section) = IO_floatValue(line,positions,2_pInt) +#ifdef NEWSTATE + case ('k_d11') + lattice_DamageDiffusion33(1,1,section) = IO_floatValue(line,positions,2_pInt) + case ('k_d12') + lattice_DamageDiffusion33(1,2,section) = IO_floatValue(line,positions,2_pInt) + case ('k_d13') + lattice_DamageDiffusion33(1,3,section) = IO_floatValue(line,positions,2_pInt) + case ('k_d21') + lattice_DamageDiffusion33(2,1,section) = IO_floatValue(line,positions,2_pInt) + case ('k_d22') + lattice_DamageDiffusion33(2,3,section) = IO_floatValue(line,positions,2_pInt) + case ('k_d23') + lattice_DamageDiffusion33(2,3,section) = IO_floatValue(line,positions,2_pInt) + case ('k_d31') + lattice_DamageDiffusion33(3,1,section) = IO_floatValue(line,positions,2_pInt) + case ('k_d32') + lattice_DamageDiffusion33(3,2,section) = IO_floatValue(line,positions,2_pInt) + case ('k_d33') + lattice_DamageDiffusion33(3,3,section) = IO_floatValue(line,positions,2_pInt) + case ('mobility_d') + lattice_DamageMobility(section) = IO_floatValue(line,positions,2_pInt) + case ('specific_heat') + lattice_specificHeat(section) = IO_floatValue(line,positions,2_pInt) + case ('mass_density') + lattice_massDensity(section) = IO_floatValue(line,positions,2_pInt) + +#endif end select endif enddo