Merge remote-tracking branch 'origin/development' into python-improvements
This commit is contained in:
Martin Diehl
commit
61204dcd93
2
PRIVATE
2
PRIVATE
|
|
@ -1 +1 @@
|
|||
Subproject commit 45ef93dbfa3e0e6fa830914b3632e188c308a099
|
||||
Subproject commit 7846c71126705cc5d41dd79f2d595f4864434068
|
||||
|
|
@ -74,9 +74,23 @@ end subroutine CPFEM_initAll
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine CPFEM_init
|
||||
|
||||
integer(HID_T) :: fileHandle
|
||||
character(len=pStringLen) :: fileName
|
||||
|
||||
|
||||
print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(IO_STDOUT)
|
||||
|
||||
if (interface_restartInc > 0) call crystallite_restartRead
|
||||
|
||||
if (interface_restartInc > 0) then
|
||||
print'(/,a,i0,a)', ' reading restart information of increment from file'; flush(IO_STDOUT)
|
||||
write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
|
||||
fileHandle = HDF5_openFile(fileName)
|
||||
|
||||
call homogenization_restartRead(fileHandle)
|
||||
call constitutive_restartRead(fileHandle)
|
||||
|
||||
call HDF5_closeFile(fileHandle)
|
||||
endif
|
||||
|
||||
end subroutine CPFEM_init
|
||||
|
||||
|
|
@ -86,7 +100,19 @@ end subroutine CPFEM_init
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine CPFEM_restartWrite
|
||||
|
||||
call crystallite_restartWrite
|
||||
integer(HID_T) :: fileHandle
|
||||
character(len=pStringLen) :: fileName
|
||||
|
||||
|
||||
print*, ' writing field and constitutive data required for restart to file';flush(IO_STDOUT)
|
||||
|
||||
write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
|
||||
fileHandle = HDF5_openFile(fileName,'a')
|
||||
|
||||
call homogenization_restartWrite(fileHandle)
|
||||
call constitutive_restartWrite(fileHandle)
|
||||
|
||||
call HDF5_closeFile(fileHandle)
|
||||
|
||||
end subroutine CPFEM_restartWrite
|
||||
|
||||
|
|
|
|||
|
|
@ -43,7 +43,7 @@ void gethostname_c(char hostname[], int *stat){
|
|||
|
||||
|
||||
void getusername_c(char username[], int *stat){
|
||||
struct passwd *pw = getpwuid(geteuid());
|
||||
struct passwd *pw = getpwuid(getuid());
|
||||
if(pw && strlen(pw->pw_name) <= STRLEN){
|
||||
strncpy(username,pw->pw_name,STRLEN+1);
|
||||
*stat = 0;
|
||||
|
|
|
|||
|
|
@ -364,7 +364,8 @@ subroutine flux(f,ts,n,time)
|
|||
real(pReal), dimension(2), intent(out) :: &
|
||||
f
|
||||
|
||||
call thermal_conduction_getSourceAndItsTangent(f(1), f(2), ts(3), n(3),mesh_FEM2DAMASK_elem(n(1)))
|
||||
f(2) = 0.0_pReal
|
||||
call thermal_conduction_getSource(f(1), ts(3), n(3),mesh_FEM2DAMASK_elem(n(1)))
|
||||
|
||||
end subroutine flux
|
||||
|
||||
|
|
|
|||
|
|
@ -32,8 +32,8 @@
|
|||
#include "constitutive_plastic_disloTungsten.f90"
|
||||
#include "constitutive_plastic_nonlocal.f90"
|
||||
#include "constitutive_thermal.f90"
|
||||
#include "source_thermal_dissipation.f90"
|
||||
#include "source_thermal_externalheat.f90"
|
||||
#include "constitutive_thermal_dissipation.f90"
|
||||
#include "constitutive_thermal_externalheat.f90"
|
||||
#include "kinematics_thermal_expansion.f90"
|
||||
#include "constitutive_damage.f90"
|
||||
#include "source_damage_isoBrittle.f90"
|
||||
|
|
@ -51,4 +51,5 @@
|
|||
#include "homogenization_mech_none.f90"
|
||||
#include "homogenization_mech_isostrain.f90"
|
||||
#include "homogenization_mech_RGC.f90"
|
||||
#include "homogenization_thermal.f90"
|
||||
#include "CPFEM.f90"
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
|
|
@ -2,6 +2,16 @@
|
|||
!> @brief internal microstructure state for all damage sources and kinematics constitutive models
|
||||
!----------------------------------------------------------------------------------------------------
|
||||
submodule(constitutive) constitutive_damage
|
||||
enum, bind(c); enumerator :: &
|
||||
DAMAGE_UNDEFINED_ID, &
|
||||
DAMAGE_ISOBRITTLE_ID, &
|
||||
DAMAGE_ISODUCTILE_ID, &
|
||||
DAMAGE_ANISOBRITTLE_ID, &
|
||||
DAMAGE_ANISODUCTILE_ID
|
||||
end enum
|
||||
|
||||
integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:,:), allocatable :: &
|
||||
phase_source !< active sources mechanisms of each phase
|
||||
|
||||
interface
|
||||
|
||||
|
|
@ -119,24 +129,24 @@ module subroutine damage_init
|
|||
|
||||
phases => config_material%get('phase')
|
||||
|
||||
allocate(sourceState (phases%length))
|
||||
allocate(damageState (phases%length))
|
||||
allocate(phase_Nsources(phases%length),source = 0) ! same for kinematics
|
||||
|
||||
do ph = 1,phases%length
|
||||
phase => phases%get(ph)
|
||||
sources => phase%get('source',defaultVal=emptyList)
|
||||
phase_Nsources(ph) = sources%length
|
||||
allocate(sourceState(ph)%p(phase_Nsources(ph)))
|
||||
allocate(damageState(ph)%p(phase_Nsources(ph)))
|
||||
enddo
|
||||
|
||||
allocate(phase_source(maxval(phase_Nsources),phases%length), source = SOURCE_undefined_ID)
|
||||
allocate(phase_source(maxval(phase_Nsources),phases%length), source = DAMAGE_UNDEFINED_ID)
|
||||
|
||||
! initialize source mechanisms
|
||||
if(maxval(phase_Nsources) /= 0) then
|
||||
where(source_damage_isoBrittle_init (maxval(phase_Nsources))) phase_source = SOURCE_damage_isoBrittle_ID
|
||||
where(source_damage_isoDuctile_init (maxval(phase_Nsources))) phase_source = SOURCE_damage_isoDuctile_ID
|
||||
where(source_damage_anisoBrittle_init (maxval(phase_Nsources))) phase_source = SOURCE_damage_anisoBrittle_ID
|
||||
where(source_damage_anisoDuctile_init (maxval(phase_Nsources))) phase_source = SOURCE_damage_anisoDuctile_ID
|
||||
where(source_damage_isoBrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISOBRITTLE_ID
|
||||
where(source_damage_isoDuctile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISODUCTILE_ID
|
||||
where(source_damage_anisoBrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISOBRITTLE_ID
|
||||
where(source_damage_anisoDuctile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISODUCTILE_ID
|
||||
endif
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -189,16 +199,16 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi
|
|||
constituent = material_phasememberAt(grain,ip,el)
|
||||
do source = 1, phase_Nsources(phase)
|
||||
select case(phase_source(source,phase))
|
||||
case (SOURCE_damage_isoBrittle_ID)
|
||||
case (DAMAGE_ISOBRITTLE_ID)
|
||||
call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
|
||||
|
||||
case (SOURCE_damage_isoDuctile_ID)
|
||||
case (DAMAGE_ISODUCTILE_ID)
|
||||
call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
|
||||
|
||||
case (SOURCE_damage_anisoBrittle_ID)
|
||||
case (DAMAGE_ANISOBRITTLE_ID)
|
||||
call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
|
||||
|
||||
case (SOURCE_damage_anisoDuctile_ID)
|
||||
case (DAMAGE_ANISODUCTILE_ID)
|
||||
call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
|
||||
|
||||
case default
|
||||
|
|
@ -214,6 +224,111 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi
|
|||
end subroutine constitutive_damage_getRateAndItsTangents
|
||||
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
|
||||
!> using Fixed Point Iteration to adapt the stepsize
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module function integrateDamageState(dt,co,ip,el) result(broken)
|
||||
|
||||
real(pReal), intent(in) :: dt
|
||||
integer, intent(in) :: &
|
||||
el, & !< element index in element loop
|
||||
ip, & !< integration point index in ip loop
|
||||
co !< grain index in grain loop
|
||||
logical :: broken
|
||||
|
||||
integer :: &
|
||||
NiterationState, & !< number of iterations in state loop
|
||||
ph, &
|
||||
me, &
|
||||
so
|
||||
integer, dimension(maxval(phase_Nsources)) :: &
|
||||
size_so
|
||||
real(pReal) :: &
|
||||
zeta
|
||||
real(pReal), dimension(constitutive_source_maxSizeDotState) :: &
|
||||
r ! state residuum
|
||||
real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState
|
||||
logical :: &
|
||||
converged_
|
||||
|
||||
|
||||
ph = material_phaseAt(co,el)
|
||||
me = material_phaseMemberAt(co,ip,el)
|
||||
|
||||
converged_ = .true.
|
||||
broken = constitutive_damage_collectDotState(co,ip,el,ph,me)
|
||||
if(broken) return
|
||||
|
||||
do so = 1, phase_Nsources(ph)
|
||||
size_so(so) = damageState(ph)%p(so)%sizeDotState
|
||||
damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%subState0(1:size_so(so),me) &
|
||||
+ damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt
|
||||
source_dotState(1:size_so(so),2,so) = 0.0_pReal
|
||||
enddo
|
||||
|
||||
iteration: do NiterationState = 1, num%nState
|
||||
|
||||
do so = 1, phase_Nsources(ph)
|
||||
if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so)
|
||||
source_dotState(1:size_so(so),1,so) = damageState(ph)%p(so)%dotState(:,me)
|
||||
enddo
|
||||
|
||||
broken = constitutive_damage_collectDotState(co,ip,el,ph,me)
|
||||
if(broken) exit iteration
|
||||
|
||||
do so = 1, phase_Nsources(ph)
|
||||
zeta = damper(damageState(ph)%p(so)%dotState(:,me), &
|
||||
source_dotState(1:size_so(so),1,so),&
|
||||
source_dotState(1:size_so(so),2,so))
|
||||
damageState(ph)%p(so)%dotState(:,me) = damageState(ph)%p(so)%dotState(:,me) * zeta &
|
||||
+ source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta)
|
||||
r(1:size_so(so)) = damageState(ph)%p(so)%state (1:size_so(so),me) &
|
||||
- damageState(ph)%p(so)%subState0(1:size_so(so),me) &
|
||||
- damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt
|
||||
damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%state(1:size_so(so),me) &
|
||||
- r(1:size_so(so))
|
||||
converged_ = converged_ .and. converged(r(1:size_so(so)), &
|
||||
damageState(ph)%p(so)%state(1:size_so(so),me), &
|
||||
damageState(ph)%p(so)%atol(1:size_so(so)))
|
||||
enddo
|
||||
|
||||
if(converged_) then
|
||||
broken = constitutive_damage_deltaState(mech_F_e(ph,me),co,ip,el,ph,me)
|
||||
exit iteration
|
||||
endif
|
||||
|
||||
enddo iteration
|
||||
|
||||
broken = broken .or. .not. converged_
|
||||
|
||||
|
||||
contains
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief calculate the damping for correction of state and dot state
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
real(pReal) pure function damper(current,previous,previous2)
|
||||
|
||||
real(pReal), dimension(:), intent(in) ::&
|
||||
current, previous, previous2
|
||||
|
||||
real(pReal) :: dot_prod12, dot_prod22
|
||||
|
||||
dot_prod12 = dot_product(current - previous, previous - previous2)
|
||||
dot_prod22 = dot_product(previous - previous2, previous - previous2)
|
||||
if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then
|
||||
damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22)
|
||||
else
|
||||
damper = 1.0_pReal
|
||||
endif
|
||||
|
||||
end function damper
|
||||
|
||||
end function integrateDamageState
|
||||
|
||||
|
||||
!----------------------------------------------------------------------------------------------
|
||||
!< @brief writes damage sources results to HDF5 output file
|
||||
!----------------------------------------------------------------------------------------------
|
||||
|
|
@ -226,23 +341,23 @@ module subroutine damage_results(group,ph)
|
|||
|
||||
sourceLoop: do so = 1, phase_Nsources(ph)
|
||||
|
||||
if (phase_source(so,ph) /= SOURCE_UNDEFINED_ID) &
|
||||
if (phase_source(so,ph) /= DAMAGE_UNDEFINED_ID) &
|
||||
call results_closeGroup(results_addGroup(group//'sources/')) ! should be 'damage'
|
||||
|
||||
sourceType: select case (phase_source(so,ph))
|
||||
|
||||
case (SOURCE_damage_anisoBrittle_ID) sourceType
|
||||
call source_damage_anisoBrittle_results(ph,group//'sources/')
|
||||
|
||||
case (SOURCE_damage_anisoDuctile_ID) sourceType
|
||||
call source_damage_anisoDuctile_results(ph,group//'sources/')
|
||||
|
||||
case (SOURCE_damage_isoBrittle_ID) sourceType
|
||||
case (DAMAGE_ISOBRITTLE_ID) sourceType
|
||||
call source_damage_isoBrittle_results(ph,group//'sources/')
|
||||
|
||||
case (SOURCE_damage_isoDuctile_ID) sourceType
|
||||
case (DAMAGE_ISODUCTILE_ID) sourceType
|
||||
call source_damage_isoDuctile_results(ph,group//'sources/')
|
||||
|
||||
case (DAMAGE_ANISOBRITTLE_ID) sourceType
|
||||
call source_damage_anisoBrittle_results(ph,group//'sources/')
|
||||
|
||||
case (DAMAGE_ANISODUCTILE_ID) sourceType
|
||||
call source_damage_anisoDuctile_results(ph,group//'sources/')
|
||||
|
||||
end select sourceType
|
||||
|
||||
enddo SourceLoop
|
||||
|
|
@ -250,4 +365,123 @@ module subroutine damage_results(group,ph)
|
|||
end subroutine damage_results
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken)
|
||||
|
||||
integer, intent(in) :: &
|
||||
co, & !< component-ID of integration point
|
||||
ip, & !< integration point
|
||||
el, & !< element
|
||||
ph, &
|
||||
of
|
||||
integer :: &
|
||||
so !< counter in source loop
|
||||
logical :: broken
|
||||
|
||||
|
||||
broken = .false.
|
||||
|
||||
SourceLoop: do so = 1, phase_Nsources(ph)
|
||||
|
||||
sourceType: select case (phase_source(so,ph))
|
||||
|
||||
case (DAMAGE_ISODUCTILE_ID) sourceType
|
||||
call source_damage_isoDuctile_dotState(co, ip, el)
|
||||
|
||||
case (DAMAGE_ANISODUCTILE_ID) sourceType
|
||||
call source_damage_anisoDuctile_dotState(co, ip, el)
|
||||
|
||||
case (DAMAGE_ANISOBRITTLE_ID) sourceType
|
||||
call source_damage_anisoBrittle_dotState(mech_S(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)),&
|
||||
co, ip, el) ! correct stress?
|
||||
|
||||
end select sourceType
|
||||
|
||||
broken = broken .or. any(IEEE_is_NaN(damageState(ph)%p(so)%dotState(:,of)))
|
||||
|
||||
enddo SourceLoop
|
||||
|
||||
end function constitutive_damage_collectDotState
|
||||
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief for constitutive models having an instantaneous change of state
|
||||
!> will return false if delta state is not needed/supported by the constitutive model
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken)
|
||||
|
||||
integer, intent(in) :: &
|
||||
co, & !< component-ID of integration point
|
||||
ip, & !< integration point
|
||||
el, & !< element
|
||||
ph, &
|
||||
of
|
||||
real(pReal), intent(in), dimension(3,3) :: &
|
||||
Fe !< elastic deformation gradient
|
||||
integer :: &
|
||||
so, &
|
||||
myOffset, &
|
||||
mySize
|
||||
logical :: &
|
||||
broken
|
||||
|
||||
|
||||
broken = .false.
|
||||
|
||||
sourceLoop: do so = 1, phase_Nsources(ph)
|
||||
|
||||
sourceType: select case (phase_source(so,ph))
|
||||
|
||||
case (DAMAGE_ISOBRITTLE_ID) sourceType
|
||||
call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(co,ip,el), Fe, &
|
||||
co, ip, el)
|
||||
broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,of)))
|
||||
if(.not. broken) then
|
||||
myOffset = damageState(ph)%p(so)%offsetDeltaState
|
||||
mySize = damageState(ph)%p(so)%sizeDeltaState
|
||||
damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) = &
|
||||
damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) + damageState(ph)%p(so)%deltaState(1:mySize,of)
|
||||
endif
|
||||
|
||||
end select sourceType
|
||||
|
||||
enddo SourceLoop
|
||||
|
||||
end function constitutive_damage_deltaState
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief checks if a source mechanism is active or not
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function source_active(source_label,src_length) result(active_source)
|
||||
|
||||
character(len=*), intent(in) :: source_label !< name of source mechanism
|
||||
integer, intent(in) :: src_length !< max. number of sources in system
|
||||
logical, dimension(:,:), allocatable :: active_source
|
||||
|
||||
class(tNode), pointer :: &
|
||||
phases, &
|
||||
phase, &
|
||||
sources, &
|
||||
src
|
||||
integer :: p,s
|
||||
|
||||
phases => config_material%get('phase')
|
||||
allocate(active_source(src_length,phases%length), source = .false. )
|
||||
do p = 1, phases%length
|
||||
phase => phases%get(p)
|
||||
sources => phase%get('source',defaultVal=emptyList)
|
||||
do s = 1, sources%length
|
||||
src => sources%get(s)
|
||||
if(src%get_asString('type') == source_label) active_source(s,p) = .true.
|
||||
enddo
|
||||
enddo
|
||||
|
||||
|
||||
end function source_active
|
||||
|
||||
|
||||
end submodule constitutive_damage
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
|
|
@ -552,10 +552,8 @@ end function plastic_nonlocal_init
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief calculates quantities characterizing the microstructure
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module subroutine plastic_nonlocal_dependentState(F, instance, of, ip, el)
|
||||
module subroutine plastic_nonlocal_dependentState(instance, of, ip, el)
|
||||
|
||||
real(pReal), dimension(3,3), intent(in) :: &
|
||||
F
|
||||
integer, intent(in) :: &
|
||||
instance, &
|
||||
of, &
|
||||
|
|
@ -647,7 +645,7 @@ module subroutine plastic_nonlocal_dependentState(F, instance, of, ip, el)
|
|||
ph = material_phaseAt(1,el)
|
||||
me = material_phaseMemberAt(1,ip,el)
|
||||
invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me))
|
||||
invFe = matmul(constitutive_mech_Fp(ph)%data(1:3,1:3,me),math_inv33(F))
|
||||
invFe = math_inv33(constitutive_mech_Fe(ph)%data(1:3,1:3,me))
|
||||
|
||||
rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg)
|
||||
rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg)
|
||||
|
|
@ -976,13 +974,11 @@ end subroutine plastic_nonlocal_deltaState
|
|||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief calculates the rate of change of microstructure
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
module subroutine plastic_nonlocal_dotState(Mp, F, Temperature,timestep, &
|
||||
module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, &
|
||||
instance,of,ip,el)
|
||||
|
||||
real(pReal), dimension(3,3), intent(in) :: &
|
||||
Mp !< MandelStress
|
||||
real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
|
||||
F !< Deformation gradient
|
||||
real(pReal), intent(in) :: &
|
||||
Temperature, & !< temperature
|
||||
timestep !< substepped crystallite time increment
|
||||
|
|
@ -1149,7 +1145,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Temperature,timestep, &
|
|||
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
|
||||
- rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have
|
||||
|
||||
rhoDot = rhoDotFlux(F,timestep, instance,of,ip,el) &
|
||||
rhoDot = rhoDotFlux(timestep, instance,of,ip,el) &
|
||||
+ rhoDotMultiplication &
|
||||
+ rhoDotSingle2DipoleGlide &
|
||||
+ rhoDotAthermalAnnihilation &
|
||||
|
|
@ -1178,10 +1174,8 @@ end subroutine plastic_nonlocal_dotState
|
|||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief calculates the rate of change of microstructure
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
function rhoDotFlux(F,timestep, instance,of,ip,el)
|
||||
function rhoDotFlux(timestep,instance,of,ip,el)
|
||||
|
||||
real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
|
||||
F !< Deformation gradient
|
||||
real(pReal), intent(in) :: &
|
||||
timestep !< substepped crystallite time increment
|
||||
integer, intent(in) :: &
|
||||
|
|
@ -1293,7 +1287,7 @@ function rhoDotFlux(F,timestep, instance,of,ip,el)
|
|||
m(1:3,:,3) = -prm%slip_transverse
|
||||
m(1:3,:,4) = prm%slip_transverse
|
||||
|
||||
my_F = F(1:3,1:3,1,ip,el)
|
||||
my_F = constitutive_mech_F(ph)%data(1:3,1:3,of)
|
||||
my_Fe = matmul(my_F, math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,of)))
|
||||
|
||||
neighbors: do n = 1,nIPneighbors
|
||||
|
|
@ -1311,7 +1305,7 @@ function rhoDotFlux(F,timestep, instance,of,ip,el)
|
|||
|
||||
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
|
||||
neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el))
|
||||
neighbor_F = F(1:3,1:3,1,neighbor_ip,neighbor_el)
|
||||
neighbor_F = constitutive_mech_F(np)%data(1:3,1:3,no)
|
||||
neighbor_Fe = matmul(neighbor_F, math_inv33(constitutive_mech_Fp(np)%data(1:3,1:3,no)))
|
||||
Favg = 0.5_pReal * (my_F + neighbor_F)
|
||||
else ! if no neighbor, take my value as average
|
||||
|
|
|
|||
|
|
@ -3,13 +3,28 @@
|
|||
!----------------------------------------------------------------------------------------------------
|
||||
submodule(constitutive) constitutive_thermal
|
||||
|
||||
enum, bind(c); enumerator :: &
|
||||
THERMAL_UNDEFINED_ID ,&
|
||||
THERMAL_DISSIPATION_ID, &
|
||||
THERMAL_EXTERNALHEAT_ID
|
||||
end enum
|
||||
|
||||
type :: tDataContainer
|
||||
real(pReal), dimension(:), allocatable :: T
|
||||
end type tDataContainer
|
||||
integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: &
|
||||
thermal_source
|
||||
|
||||
type(tDataContainer), dimension(:), allocatable :: current
|
||||
|
||||
integer :: thermal_source_maxSizeDotState
|
||||
interface
|
||||
|
||||
module function source_thermal_dissipation_init(source_length) result(mySources)
|
||||
integer, intent(in) :: source_length
|
||||
logical, dimension(:,:), allocatable :: mySources
|
||||
end function source_thermal_dissipation_init
|
||||
|
||||
|
||||
module function source_thermal_externalheat_init(source_length) result(mySources)
|
||||
integer, intent(in) :: source_length
|
||||
logical, dimension(:,:), allocatable :: mySources
|
||||
|
|
@ -21,7 +36,7 @@ submodule(constitutive) constitutive_thermal
|
|||
end function kinematics_thermal_expansion_init
|
||||
|
||||
|
||||
module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase)
|
||||
module subroutine thermal_dissipation_getRate(TDot, Tstar,Lp,phase)
|
||||
integer, intent(in) :: &
|
||||
phase !< phase ID of element
|
||||
real(pReal), intent(in), dimension(3,3) :: &
|
||||
|
|
@ -29,18 +44,16 @@ submodule(constitutive) constitutive_thermal
|
|||
real(pReal), intent(in), dimension(3,3) :: &
|
||||
Lp !< plastic velocuty gradient for a given element
|
||||
real(pReal), intent(out) :: &
|
||||
TDot, &
|
||||
dTDot_dT
|
||||
end subroutine source_thermal_dissipation_getRateAndItsTangent
|
||||
TDot
|
||||
end subroutine thermal_dissipation_getRate
|
||||
|
||||
module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of)
|
||||
module subroutine thermal_externalheat_getRate(TDot, phase,of)
|
||||
integer, intent(in) :: &
|
||||
phase, &
|
||||
of
|
||||
real(pReal), intent(out) :: &
|
||||
TDot, &
|
||||
dTDot_dT
|
||||
end subroutine source_thermal_externalheat_getRateAndItsTangent
|
||||
TDot
|
||||
end subroutine thermal_externalheat_getRate
|
||||
|
||||
end interface
|
||||
|
||||
|
|
@ -49,14 +62,60 @@ contains
|
|||
!----------------------------------------------------------------------------------------------
|
||||
!< @brief initializes thermal sources and kinematics mechanism
|
||||
!----------------------------------------------------------------------------------------------
|
||||
module subroutine thermal_init
|
||||
module subroutine thermal_init(phases)
|
||||
|
||||
class(tNode), pointer :: &
|
||||
phases
|
||||
|
||||
class(tNode), pointer :: &
|
||||
phase, thermal, sources
|
||||
|
||||
integer :: &
|
||||
ph, so, &
|
||||
Nconstituents
|
||||
|
||||
|
||||
print'(/,a)', ' <<<+- constitutive_thermal init -+>>>'
|
||||
|
||||
allocate(current(phases%length))
|
||||
|
||||
allocate(thermalState (phases%length))
|
||||
allocate(thermal_Nsources(phases%length),source = 0)
|
||||
|
||||
do ph = 1, phases%length
|
||||
|
||||
Nconstituents = count(material_phaseAt == ph) * discretization_nIPs
|
||||
|
||||
allocate(current(ph)%T(Nconstituents))
|
||||
phase => phases%get(ph)
|
||||
if(phase%contains('thermal')) then
|
||||
thermal => phase%get('thermal')
|
||||
sources => thermal%get('source',defaultVal=emptyList)
|
||||
|
||||
thermal_Nsources(ph) = sources%length
|
||||
endif
|
||||
allocate(thermalstate(ph)%p(thermal_Nsources(ph)))
|
||||
enddo
|
||||
|
||||
allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID)
|
||||
|
||||
if(maxval(thermal_Nsources) /= 0) then
|
||||
where(source_thermal_dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID
|
||||
where(source_thermal_externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID
|
||||
endif
|
||||
|
||||
thermal_source_maxSizeDotState = 0
|
||||
PhaseLoop2:do ph = 1,phases%length
|
||||
|
||||
do so = 1,thermal_Nsources(ph)
|
||||
thermalState(ph)%p(so)%partitionedState0 = thermalState(ph)%p(so)%state0
|
||||
thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%partitionedState0
|
||||
enddo
|
||||
|
||||
thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, &
|
||||
maxval(thermalState(ph)%p%sizeDotState))
|
||||
enddo PhaseLoop2
|
||||
|
||||
! initialize source mechanisms
|
||||
if(maxval(phase_Nsources) /= 0) then
|
||||
where(source_thermal_dissipation_init (maxval(phase_Nsources))) phase_source = SOURCE_thermal_dissipation_ID
|
||||
where(source_thermal_externalheat_init(maxval(phase_Nsources))) phase_source = SOURCE_thermal_externalheat_ID
|
||||
endif
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!initialize kinematic mechanisms
|
||||
if(maxval(phase_Nkinematics) /= 0) where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) &
|
||||
|
|
@ -68,58 +127,236 @@ end subroutine thermal_init
|
|||
!----------------------------------------------------------------------------------------------
|
||||
!< @brief calculates thermal dissipation rate
|
||||
!----------------------------------------------------------------------------------------------
|
||||
module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, S, Lp, ip, el)
|
||||
module subroutine constitutive_thermal_getRate(TDot, T, ip, el)
|
||||
|
||||
integer, intent(in) :: &
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
real(pReal), intent(in) :: &
|
||||
T
|
||||
real(pReal), intent(in), dimension(:,:,:,:,:) :: &
|
||||
S, & !< current 2nd Piola Kirchhoff stress
|
||||
Lp !< plastic velocity gradient
|
||||
real(pReal), intent(inout) :: &
|
||||
TDot, &
|
||||
dTDot_dT
|
||||
T !< plastic velocity gradient
|
||||
real(pReal), intent(out) :: &
|
||||
TDot
|
||||
|
||||
real(pReal) :: &
|
||||
my_Tdot, &
|
||||
my_dTdot_dT
|
||||
my_Tdot
|
||||
integer :: &
|
||||
phase, &
|
||||
ph, &
|
||||
homog, &
|
||||
instance, &
|
||||
grain, &
|
||||
source, &
|
||||
constituent
|
||||
me, &
|
||||
so, &
|
||||
co
|
||||
|
||||
homog = material_homogenizationAt(el)
|
||||
instance = thermal_typeInstance(homog)
|
||||
|
||||
do grain = 1, homogenization_Nconstituents(homog)
|
||||
phase = material_phaseAt(grain,el)
|
||||
constituent = material_phasememberAt(grain,ip,el)
|
||||
do source = 1, phase_Nsources(phase)
|
||||
select case(phase_source(source,phase))
|
||||
case (SOURCE_thermal_dissipation_ID)
|
||||
call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
|
||||
S(1:3,1:3,grain,ip,el), &
|
||||
Lp(1:3,1:3,grain,ip,el), &
|
||||
phase)
|
||||
TDot = 0.0_pReal
|
||||
do co = 1, homogenization_Nconstituents(homog)
|
||||
ph = material_phaseAt(co,el)
|
||||
me = material_phasememberAt(co,ip,el)
|
||||
do so = 1, thermal_Nsources(ph)
|
||||
select case(thermal_source(so,ph))
|
||||
case (THERMAL_DISSIPATION_ID)
|
||||
call thermal_dissipation_getRate(my_Tdot, mech_S(ph,me),mech_L_p(ph,me),ph)
|
||||
|
||||
case (SOURCE_thermal_externalheat_ID)
|
||||
call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
|
||||
phase, constituent)
|
||||
case (THERMAL_EXTERNALHEAT_ID)
|
||||
call thermal_externalheat_getRate(my_Tdot, ph,me)
|
||||
|
||||
case default
|
||||
my_Tdot = 0.0_pReal
|
||||
my_dTdot_dT = 0.0_pReal
|
||||
end select
|
||||
Tdot = Tdot + my_Tdot
|
||||
dTdot_dT = dTdot_dT + my_dTdot_dT
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end subroutine constitutive_thermal_getRateAndItsTangents
|
||||
end subroutine constitutive_thermal_getRate
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function constitutive_thermal_collectDotState(ph,me) result(broken)
|
||||
|
||||
integer, intent(in) :: ph, me
|
||||
logical :: broken
|
||||
|
||||
integer :: i
|
||||
|
||||
|
||||
broken = .false.
|
||||
|
||||
SourceLoop: do i = 1, thermal_Nsources(ph)
|
||||
|
||||
if (thermal_source(i,ph) == THERMAL_EXTERNALHEAT_ID) &
|
||||
call source_thermal_externalheat_dotState(ph,me)
|
||||
|
||||
broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,me)))
|
||||
|
||||
enddo SourceLoop
|
||||
|
||||
end function constitutive_thermal_collectDotState
|
||||
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief integrate state with 1st order explicit Euler method
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module function integrateThermalState(Delta_t,co,ip,el) result(broken)
|
||||
|
||||
real(pReal), intent(in) :: Delta_t
|
||||
integer, intent(in) :: &
|
||||
el, & !< element index in element loop
|
||||
ip, & !< integration point index in ip loop
|
||||
co !< grain index in grain loop
|
||||
logical :: &
|
||||
broken
|
||||
|
||||
integer :: &
|
||||
ph, &
|
||||
me, &
|
||||
so, &
|
||||
sizeDotState
|
||||
|
||||
|
||||
ph = material_phaseAt(co,el)
|
||||
me = material_phaseMemberAt(co,ip,el)
|
||||
|
||||
broken = constitutive_thermal_collectDotState(ph,me)
|
||||
if(broken) return
|
||||
|
||||
do so = 1, thermal_Nsources(ph)
|
||||
sizeDotState = thermalState(ph)%p(so)%sizeDotState
|
||||
thermalState(ph)%p(so)%state(1:sizeDotState,me) = thermalState(ph)%p(so)%subState0(1:sizeDotState,me) &
|
||||
+ thermalState(ph)%p(so)%dotState(1:sizeDotState,me) * Delta_t
|
||||
enddo
|
||||
|
||||
end function integrateThermalState
|
||||
|
||||
|
||||
module subroutine thermal_initializeRestorationPoints(ph,me)
|
||||
|
||||
integer, intent(in) :: ph, me
|
||||
|
||||
integer :: so
|
||||
|
||||
|
||||
do so = 1, size(thermalState(ph)%p)
|
||||
thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state0(:,me)
|
||||
enddo
|
||||
|
||||
end subroutine thermal_initializeRestorationPoints
|
||||
|
||||
|
||||
|
||||
module subroutine thermal_windForward(ph,me)
|
||||
|
||||
integer, intent(in) :: ph, me
|
||||
|
||||
integer :: so
|
||||
|
||||
|
||||
do so = 1, size(thermalState(ph)%p)
|
||||
thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state(:,me)
|
||||
enddo
|
||||
|
||||
end subroutine thermal_windForward
|
||||
|
||||
|
||||
module subroutine thermal_forward()
|
||||
|
||||
integer :: ph, so
|
||||
|
||||
|
||||
do ph = 1, size(thermalState)
|
||||
do so = 1, size(thermalState(ph)%p)
|
||||
thermalState(ph)%p(so)%state0 = thermalState(ph)%p(so)%state
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end subroutine thermal_forward
|
||||
|
||||
|
||||
module subroutine thermal_restore(ip,el)
|
||||
|
||||
integer, intent(in) :: ip, el
|
||||
|
||||
integer :: co, ph, me, so
|
||||
|
||||
|
||||
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
ph = material_phaseAt(co,el)
|
||||
me = material_phaseMemberAt(co,ip,el)
|
||||
|
||||
do so = 1, size(thermalState(ph)%p)
|
||||
thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%partitionedState0(:,me)
|
||||
enddo
|
||||
|
||||
enddo
|
||||
|
||||
end subroutine thermal_restore
|
||||
|
||||
|
||||
!----------------------------------------------------------------------------------------------
|
||||
!< @brief Get temperature (for use by non-thermal physics)
|
||||
!----------------------------------------------------------------------------------------------
|
||||
module function thermal_T(ph,me) result(T)
|
||||
|
||||
integer, intent(in) :: ph, me
|
||||
real(pReal) :: T
|
||||
|
||||
|
||||
T = current(ph)%T(me)
|
||||
|
||||
end function thermal_T
|
||||
|
||||
|
||||
!----------------------------------------------------------------------------------------------
|
||||
!< @brief Set temperature
|
||||
!----------------------------------------------------------------------------------------------
|
||||
module subroutine constitutive_thermal_setT(T,co,ip,el)
|
||||
|
||||
real(pReal), intent(in) :: T
|
||||
integer, intent(in) :: co, ip, el
|
||||
|
||||
|
||||
current(material_phaseAt(co,el))%T(material_phaseMemberAt(co,ip,el)) = T
|
||||
|
||||
end subroutine constitutive_thermal_setT
|
||||
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief checks if a source mechanism is active or not
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function thermal_active(source_label,src_length) result(active_source)
|
||||
|
||||
character(len=*), intent(in) :: source_label !< name of source mechanism
|
||||
integer, intent(in) :: src_length !< max. number of sources in system
|
||||
logical, dimension(:,:), allocatable :: active_source
|
||||
|
||||
class(tNode), pointer :: &
|
||||
phases, &
|
||||
phase, &
|
||||
sources, thermal, &
|
||||
src
|
||||
integer :: p,s
|
||||
|
||||
phases => config_material%get('phase')
|
||||
allocate(active_source(src_length,phases%length), source = .false. )
|
||||
do p = 1, phases%length
|
||||
phase => phases%get(p)
|
||||
if (phase%contains('thermal')) then
|
||||
thermal => phase%get('thermal',defaultVal=emptyList)
|
||||
sources => thermal%get('source',defaultVal=emptyList)
|
||||
do s = 1, sources%length
|
||||
src => sources%get(s)
|
||||
if(src%get_asString('type') == source_label) active_source(s,p) = .true.
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
|
||||
|
||||
end function thermal_active
|
||||
|
||||
|
||||
end submodule constitutive_thermal
|
||||
|
|
|
|||
|
|
@ -4,7 +4,7 @@
|
|||
!> @brief material subroutine for thermal source due to plastic dissipation
|
||||
!> @details to be done
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
submodule(constitutive:constitutive_thermal) source_thermal_dissipation
|
||||
submodule(constitutive:constitutive_thermal) source_dissipation
|
||||
|
||||
integer, dimension(:), allocatable :: &
|
||||
source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism?
|
||||
|
|
@ -27,19 +27,20 @@ contains
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
module function source_thermal_dissipation_init(source_length) result(mySources)
|
||||
|
||||
integer, intent(in) :: source_length
|
||||
integer, intent(in) :: source_length
|
||||
logical, dimension(:,:), allocatable :: mySources
|
||||
|
||||
class(tNode), pointer :: &
|
||||
phases, &
|
||||
phase, &
|
||||
sources, &
|
||||
src
|
||||
sources, thermal, &
|
||||
src
|
||||
integer :: Ninstances,sourceOffset,Nconstituents,p
|
||||
|
||||
print'(/,a)', ' <<<+- source_thermal_dissipation init -+>>>'
|
||||
print'(/,a)', ' <<<+- thermal_dissipation init -+>>>'
|
||||
|
||||
mySources = thermal_active('dissipation',source_length)
|
||||
|
||||
mySources = source_active('thermal_dissipation',source_length)
|
||||
Ninstances = count(mySources)
|
||||
print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
|
||||
if(Ninstances == 0) return
|
||||
|
|
@ -50,19 +51,20 @@ module function source_thermal_dissipation_init(source_length) result(mySources)
|
|||
allocate(source_thermal_dissipation_instance(phases%length), source=0)
|
||||
|
||||
do p = 1, phases%length
|
||||
phase => phases%get(p)
|
||||
if(count(mySources(:,p)) == 0) cycle
|
||||
phase => phases%get(p)
|
||||
if(any(mySources(:,p))) source_thermal_dissipation_instance(p) = count(mySources(:,1:p))
|
||||
sources => phase%get('source')
|
||||
if(count(mySources(:,p)) == 0) cycle
|
||||
thermal => phase%get('thermal')
|
||||
sources => thermal%get('source')
|
||||
do sourceOffset = 1, sources%length
|
||||
if(mySources(sourceOffset,p)) then
|
||||
source_thermal_dissipation_offset(p) = sourceOffset
|
||||
associate(prm => param(source_thermal_dissipation_instance(p)))
|
||||
src => sources%get(sourceOffset)
|
||||
|
||||
src => sources%get(sourceOffset)
|
||||
prm%kappa = src%get_asFloat('kappa')
|
||||
Nconstituents = count(material_phaseAt==p) * discretization_nIPs
|
||||
call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,0,0,0)
|
||||
call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,0,0,0)
|
||||
|
||||
end associate
|
||||
endif
|
||||
|
|
@ -76,7 +78,7 @@ end function source_thermal_dissipation_init
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Ninstancess dissipation rate
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase)
|
||||
module subroutine thermal_dissipation_getRate(TDot, Tstar, Lp, phase)
|
||||
|
||||
integer, intent(in) :: &
|
||||
phase
|
||||
|
|
@ -86,14 +88,12 @@ module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT
|
|||
Lp
|
||||
|
||||
real(pReal), intent(out) :: &
|
||||
TDot, &
|
||||
dTDot_dT
|
||||
TDot
|
||||
|
||||
associate(prm => param(source_thermal_dissipation_instance(phase)))
|
||||
TDot = prm%kappa*sum(abs(Tstar*Lp))
|
||||
dTDot_dT = 0.0_pReal
|
||||
end associate
|
||||
|
||||
end subroutine source_thermal_dissipation_getRateAndItsTangent
|
||||
end subroutine thermal_dissipation_getRate
|
||||
|
||||
end submodule source_thermal_dissipation
|
||||
end submodule source_dissipation
|
||||
|
|
@ -4,7 +4,7 @@
|
|||
!> @author Philip Eisenlohr, Michigan State University
|
||||
!> @brief material subroutine for variable heat source
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
submodule(constitutive:constitutive_thermal) source_thermal_externalheat
|
||||
submodule(constitutive:constitutive_thermal) source_externalheat
|
||||
|
||||
|
||||
integer, dimension(:), allocatable :: &
|
||||
|
|
@ -13,7 +13,7 @@ submodule(constitutive:constitutive_thermal) source_thermal_externalheat
|
|||
|
||||
type :: tParameters !< container type for internal constitutive parameters
|
||||
real(pReal), dimension(:), allocatable :: &
|
||||
t_n, &
|
||||
t_n, &
|
||||
f_T
|
||||
integer :: &
|
||||
nIntervals
|
||||
|
|
@ -31,19 +31,20 @@ contains
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
module function source_thermal_externalheat_init(source_length) result(mySources)
|
||||
|
||||
integer, intent(in) :: source_length
|
||||
integer, intent(in) :: source_length
|
||||
logical, dimension(:,:), allocatable :: mySources
|
||||
|
||||
class(tNode), pointer :: &
|
||||
phases, &
|
||||
phase, &
|
||||
sources, &
|
||||
src
|
||||
sources, thermal, &
|
||||
src
|
||||
integer :: Ninstances,sourceOffset,Nconstituents,p
|
||||
|
||||
print'(/,a)', ' <<<+- source_thermal_externalHeat init -+>>>'
|
||||
print'(/,a)', ' <<<+- thermal_externalheat init -+>>>'
|
||||
|
||||
mySources = thermal_active('externalheat',source_length)
|
||||
|
||||
mySources = source_active('thermal_externalheat',source_length)
|
||||
Ninstances = count(mySources)
|
||||
print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
|
||||
if(Ninstances == 0) return
|
||||
|
|
@ -54,15 +55,16 @@ module function source_thermal_externalheat_init(source_length) result(mySources
|
|||
allocate(source_thermal_externalheat_instance(phases%length), source=0)
|
||||
|
||||
do p = 1, phases%length
|
||||
phase => phases%get(p)
|
||||
phase => phases%get(p)
|
||||
if(any(mySources(:,p))) source_thermal_externalheat_instance(p) = count(mySources(:,1:p))
|
||||
if(count(mySources(:,p)) == 0) cycle
|
||||
sources => phase%get('source')
|
||||
thermal => phase%get('thermal')
|
||||
sources => thermal%get('source')
|
||||
do sourceOffset = 1, sources%length
|
||||
if(mySources(sourceOffset,p)) then
|
||||
source_thermal_externalheat_offset(p) = sourceOffset
|
||||
associate(prm => param(source_thermal_externalheat_instance(p)))
|
||||
src => sources%get(sourceOffset)
|
||||
src => sources%get(sourceOffset)
|
||||
|
||||
prm%t_n = src%get_asFloats('t_n')
|
||||
prm%nIntervals = size(prm%t_n) - 1
|
||||
|
|
@ -70,9 +72,8 @@ module function source_thermal_externalheat_init(source_length) result(mySources
|
|||
prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n))
|
||||
|
||||
Nconstituents = count(material_phaseAt==p) * discretization_nIPs
|
||||
call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0)
|
||||
call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,1,1,0)
|
||||
end associate
|
||||
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
|
@ -95,7 +96,7 @@ module subroutine source_thermal_externalheat_dotState(phase, of)
|
|||
|
||||
sourceOffset = source_thermal_externalheat_offset(phase)
|
||||
|
||||
sourceState(phase)%p(sourceOffset)%dotState(1,of) = 1.0_pReal ! state is current time
|
||||
thermalState(phase)%p(sourceOffset)%dotState(1,of) = 1.0_pReal ! state is current time
|
||||
|
||||
end subroutine source_thermal_externalheat_dotState
|
||||
|
||||
|
|
@ -103,14 +104,13 @@ end subroutine source_thermal_externalheat_dotState
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns local heat generation rate
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of)
|
||||
module subroutine thermal_externalheat_getRate(TDot, phase, of)
|
||||
|
||||
integer, intent(in) :: &
|
||||
phase, &
|
||||
of
|
||||
real(pReal), intent(out) :: &
|
||||
TDot, &
|
||||
dTDot_dT
|
||||
TDot
|
||||
|
||||
integer :: &
|
||||
sourceOffset, interval
|
||||
|
|
@ -121,7 +121,7 @@ module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_d
|
|||
|
||||
associate(prm => param(source_thermal_externalheat_instance(phase)))
|
||||
do interval = 1, prm%nIntervals ! scan through all rate segments
|
||||
frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) &
|
||||
frac_time = (thermalState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) &
|
||||
/ (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment
|
||||
if ( (frac_time < 0.0_pReal .and. interval == 1) &
|
||||
.or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) &
|
||||
|
|
@ -130,9 +130,8 @@ module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_d
|
|||
prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries...
|
||||
! ...or extrapolate if outside of bounds
|
||||
enddo
|
||||
dTDot_dT = 0.0
|
||||
end associate
|
||||
|
||||
end subroutine source_thermal_externalheat_getRateAndItsTangent
|
||||
end subroutine thermal_externalheat_getRate
|
||||
|
||||
end submodule source_thermal_externalheat
|
||||
end submodule source_externalheat
|
||||
|
|
@ -25,10 +25,10 @@ subroutine damage_none_init
|
|||
if (damage_type(h) /= DAMAGE_NONE_ID) cycle
|
||||
|
||||
Nmaterialpoints = count(material_homogenizationAt == h)
|
||||
damageState(h)%sizeState = 0
|
||||
allocate(damageState(h)%state0 (0,Nmaterialpoints))
|
||||
allocate(damageState(h)%subState0(0,Nmaterialpoints))
|
||||
allocate(damageState(h)%state (0,Nmaterialpoints))
|
||||
damageState_h(h)%sizeState = 0
|
||||
allocate(damageState_h(h)%state0 (0,Nmaterialpoints))
|
||||
allocate(damageState_h(h)%subState0(0,Nmaterialpoints))
|
||||
allocate(damageState_h(h)%state (0,Nmaterialpoints))
|
||||
|
||||
allocate (damage(h)%p(Nmaterialpoints), source=1.0_pReal)
|
||||
|
||||
|
|
|
|||
|
|
@ -76,12 +76,12 @@ subroutine damage_nonlocal_init
|
|||
#endif
|
||||
|
||||
Nmaterialpoints = count(material_homogenizationAt == h)
|
||||
damageState(h)%sizeState = 1
|
||||
allocate(damageState(h)%state0 (1,Nmaterialpoints), source=1.0_pReal)
|
||||
allocate(damageState(h)%subState0(1,Nmaterialpoints), source=1.0_pReal)
|
||||
allocate(damageState(h)%state (1,Nmaterialpoints), source=1.0_pReal)
|
||||
damageState_h(h)%sizeState = 1
|
||||
allocate(damageState_h(h)%state0 (1,Nmaterialpoints), source=1.0_pReal)
|
||||
allocate(damageState_h(h)%subState0(1,Nmaterialpoints), source=1.0_pReal)
|
||||
allocate(damageState_h(h)%state (1,Nmaterialpoints), source=1.0_pReal)
|
||||
|
||||
damage(h)%p => damageState(h)%state(1,:)
|
||||
damage(h)%p => damageState_h(h)%state(1,:)
|
||||
|
||||
end associate
|
||||
enddo
|
||||
|
|
|
|||
|
|
@ -256,7 +256,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
|
|||
PetscObject :: dummy
|
||||
PetscErrorCode :: ierr
|
||||
integer :: i, j, k, cell
|
||||
real(pReal) :: Tdot, dTdot_dT
|
||||
real(pReal) :: Tdot
|
||||
|
||||
T_current = x_scal
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -278,7 +278,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
|
|||
cell = 0
|
||||
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
|
||||
cell = cell + 1
|
||||
call thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T_current(i,j,k), 1, cell)
|
||||
call thermal_conduction_getSource(Tdot, T_current(i,j,k), 1, cell)
|
||||
scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) &
|
||||
+ thermal_conduction_getMassDensity (1,cell)* &
|
||||
thermal_conduction_getSpecificHeat(1,cell)*(T_lastInc(i,j,k) - &
|
||||
|
|
|
|||
|
|
@ -16,6 +16,7 @@ module homogenization
|
|||
use thermal_conduction
|
||||
use damage_none
|
||||
use damage_nonlocal
|
||||
use HDF5_utilities
|
||||
use results
|
||||
|
||||
implicit none
|
||||
|
|
@ -26,6 +27,8 @@ module homogenization
|
|||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! General variables for the homogenization at a material point
|
||||
real(pReal), dimension(:), allocatable, public :: &
|
||||
homogenization_T
|
||||
real(pReal), dimension(:,:,:), allocatable, public :: &
|
||||
homogenization_F0, & !< def grad of IP at start of FE increment
|
||||
homogenization_F !< def grad of IP to be reached at end of FE increment
|
||||
|
|
@ -55,6 +58,9 @@ module homogenization
|
|||
num_homog !< pointer to mechanical homogenization numerics data
|
||||
end subroutine mech_init
|
||||
|
||||
module subroutine thermal_init
|
||||
end subroutine thermal_init
|
||||
|
||||
module subroutine mech_partition(subF,ip,el)
|
||||
real(pReal), intent(in), dimension(3,3) :: &
|
||||
subF
|
||||
|
|
@ -63,7 +69,15 @@ module homogenization
|
|||
el !< element number
|
||||
end subroutine mech_partition
|
||||
|
||||
module subroutine mech_homogenize(ip,el)
|
||||
module subroutine thermal_partition(T,ip,el)
|
||||
real(pReal), intent(in) :: T
|
||||
integer, intent(in) :: &
|
||||
ip, & !< integration point
|
||||
el !< element number
|
||||
end subroutine thermal_partition
|
||||
|
||||
module subroutine mech_homogenize(dt,ip,el)
|
||||
real(pReal), intent(in) :: dt
|
||||
integer, intent(in) :: &
|
||||
ip, & !< integration point
|
||||
el !< element number
|
||||
|
|
@ -91,7 +105,9 @@ module homogenization
|
|||
homogenization_init, &
|
||||
materialpoint_stressAndItsTangent, &
|
||||
homogenization_forward, &
|
||||
homogenization_results
|
||||
homogenization_results, &
|
||||
homogenization_restartRead, &
|
||||
homogenization_restartWrite
|
||||
|
||||
contains
|
||||
|
||||
|
|
@ -122,9 +138,10 @@ subroutine homogenization_init
|
|||
|
||||
|
||||
call mech_init(num_homog)
|
||||
call thermal_init()
|
||||
|
||||
if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init
|
||||
if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init
|
||||
if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init(homogenization_T)
|
||||
if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init(homogenization_T)
|
||||
|
||||
if (any(damage_type == DAMAGE_none_ID)) call damage_none_init
|
||||
if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init
|
||||
|
|
@ -168,10 +185,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
|
|||
converged = .false. ! pretend failed step ...
|
||||
subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation
|
||||
|
||||
if (homogState(ho)%sizeState > 0) &
|
||||
homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me)
|
||||
if (damageState(ho)%sizeState > 0) &
|
||||
damageState(ho)%subState0(:,me) = damageState(ho)%State0(:,me)
|
||||
if (homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me)
|
||||
if (damageState_h(ho)%sizeState > 0) damageState_h(ho)%subState0(:,me) = damageState_h(ho)%State0(:,me)
|
||||
|
||||
cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog)
|
||||
|
||||
|
|
@ -184,10 +199,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
|
|||
! wind forward grain starting point
|
||||
call constitutive_windForward(ip,el)
|
||||
|
||||
if(homogState(ho)%sizeState > 0) &
|
||||
homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me)
|
||||
if(damageState(ho)%sizeState > 0) &
|
||||
damageState(ho)%subState0(:,me) = damageState(ho)%State(:,me)
|
||||
if(homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me)
|
||||
if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%subState0(:,me) = damageState_h(ho)%State(:,me)
|
||||
|
||||
endif steppingNeeded
|
||||
elseif ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite
|
||||
|
|
@ -201,10 +214,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
|
|||
|
||||
call constitutive_restore(ip,el,subStep < 1.0_pReal)
|
||||
|
||||
if(homogState(ho)%sizeState > 0) &
|
||||
homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me)
|
||||
if(damageState(ho)%sizeState > 0) &
|
||||
damageState(ho)%State(:,me) = damageState(ho)%subState0(:,me)
|
||||
if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me)
|
||||
if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%subState0(:,me)
|
||||
endif
|
||||
|
||||
if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.]
|
||||
|
|
@ -257,7 +268,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
|
|||
do co = 1, myNgrains
|
||||
call crystallite_orientations(co,ip,el)
|
||||
enddo
|
||||
call mech_homogenize(ip,el)
|
||||
call mech_homogenize(dt,ip,el)
|
||||
enddo IpLooping3
|
||||
enddo elementLooping3
|
||||
!$OMP END PARALLEL DO
|
||||
|
|
@ -315,9 +326,64 @@ subroutine homogenization_forward
|
|||
|
||||
do ho = 1, size(material_name_homogenization)
|
||||
homogState (ho)%state0 = homogState (ho)%state
|
||||
damageState(ho)%state0 = damageState(ho)%state
|
||||
damageState_h(ho)%state0 = damageState_h(ho)%state
|
||||
enddo
|
||||
|
||||
end subroutine homogenization_forward
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine homogenization_restartWrite(fileHandle)
|
||||
|
||||
integer(HID_T), intent(in) :: fileHandle
|
||||
|
||||
integer(HID_T), dimension(2) :: groupHandle
|
||||
integer :: ho
|
||||
|
||||
|
||||
groupHandle(1) = HDF5_addGroup(fileHandle,'homogenization')
|
||||
|
||||
do ho = 1, size(material_name_homogenization)
|
||||
|
||||
groupHandle(2) = HDF5_addGroup(groupHandle(1),material_name_homogenization(ho))
|
||||
|
||||
call HDF5_read(groupHandle(2),homogState(ho)%state,'omega') ! ToDo: should be done by mech
|
||||
|
||||
call HDF5_closeGroup(groupHandle(2))
|
||||
|
||||
enddo
|
||||
|
||||
call HDF5_closeGroup(groupHandle(1))
|
||||
|
||||
end subroutine homogenization_restartWrite
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine homogenization_restartRead(fileHandle)
|
||||
|
||||
integer(HID_T), intent(in) :: fileHandle
|
||||
|
||||
integer(HID_T), dimension(2) :: groupHandle
|
||||
integer :: ho
|
||||
|
||||
|
||||
groupHandle(1) = HDF5_openGroup(fileHandle,'homogenization')
|
||||
|
||||
do ho = 1, size(material_name_homogenization)
|
||||
|
||||
groupHandle(2) = HDF5_openGroup(groupHandle(1),material_name_homogenization(ho))
|
||||
|
||||
call HDF5_write(groupHandle(2),homogState(ho)%state,'omega') ! ToDo: should be done by mech
|
||||
|
||||
call HDF5_closeGroup(groupHandle(2))
|
||||
|
||||
enddo
|
||||
|
||||
call HDF5_closeGroup(groupHandle(1))
|
||||
|
||||
end subroutine homogenization_restartRead
|
||||
|
||||
|
||||
end module homogenization
|
||||
|
|
|
|||
|
|
@ -52,12 +52,11 @@ submodule(homogenization) homogenization_mech
|
|||
end subroutine mech_RGC_averageStressAndItsTangent
|
||||
|
||||
|
||||
module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHappy)
|
||||
module function mech_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy)
|
||||
logical, dimension(2) :: doneAndHappy
|
||||
real(pReal), dimension(:,:,:), intent(in) :: &
|
||||
P,& !< partitioned stresses
|
||||
F,& !< partitioned deformation gradients
|
||||
F0 !< partitioned initial deformation gradients
|
||||
F !< partitioned deformation gradients
|
||||
real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
|
||||
real(pReal), dimension(3,3), intent(in) :: avgF !< average F
|
||||
real(pReal), intent(in) :: dt !< time increment
|
||||
|
|
@ -113,67 +112,73 @@ module subroutine mech_partition(subF,ip,el)
|
|||
ip, & !< integration point
|
||||
el !< element number
|
||||
|
||||
integer :: co
|
||||
real(pReal), dimension (3,3,homogenization_Nconstituents(material_homogenizationAt(el))) :: Fs
|
||||
|
||||
|
||||
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
|
||||
|
||||
case (HOMOGENIZATION_NONE_ID) chosenHomogenization
|
||||
crystallite_F(1:3,1:3,1,ip,el) = subF
|
||||
Fs(1:3,1:3,1) = subF
|
||||
|
||||
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
|
||||
call mech_isostrain_partitionDeformation(&
|
||||
crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
subF)
|
||||
call mech_isostrain_partitionDeformation(Fs,subF)
|
||||
|
||||
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
|
||||
call mech_RGC_partitionDeformation(&
|
||||
crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
subF,&
|
||||
ip, &
|
||||
el)
|
||||
call mech_RGC_partitionDeformation(Fs,subF,ip,el)
|
||||
|
||||
end select chosenHomogenization
|
||||
|
||||
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
call constitutive_mech_setF(Fs(1:3,1:3,co),co,ip,el)
|
||||
enddo
|
||||
|
||||
|
||||
end subroutine mech_partition
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Average P and dPdF from the individual constituents.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module subroutine mech_homogenize(ip,el)
|
||||
module subroutine mech_homogenize(dt,ip,el)
|
||||
|
||||
real(pReal), intent(in) :: dt
|
||||
integer, intent(in) :: &
|
||||
ip, & !< integration point
|
||||
el !< element number
|
||||
|
||||
integer :: co,ce
|
||||
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
|
||||
real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
|
||||
|
||||
|
||||
ce = (el-1)* discretization_nIPs + ip
|
||||
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
|
||||
|
||||
case (HOMOGENIZATION_NONE_ID) chosenHomogenization
|
||||
homogenization_P(1:3,1:3,ce) = crystallite_P(1:3,1:3,1,ip,el)
|
||||
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = crystallite_stressTangent(1,ip,el)
|
||||
homogenization_P(1:3,1:3,ce) = constitutive_mech_getP(1,ip,el)
|
||||
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = constitutive_mech_dPdF(dt,1,ip,el)
|
||||
|
||||
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
|
||||
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el)
|
||||
dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el)
|
||||
Ps(:,:,co) = constitutive_mech_getP(co,ip,el)
|
||||
enddo
|
||||
call mech_isostrain_averageStressAndItsTangent(&
|
||||
homogenization_P(1:3,1:3,ce), &
|
||||
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
|
||||
crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
dPdFs, &
|
||||
Ps,dPdFs, &
|
||||
homogenization_typeInstance(material_homogenizationAt(el)))
|
||||
|
||||
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
|
||||
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el)
|
||||
dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el)
|
||||
Ps(:,:,co) = constitutive_mech_getP(co,ip,el)
|
||||
enddo
|
||||
call mech_RGC_averageStressAndItsTangent(&
|
||||
homogenization_P(1:3,1:3,ce), &
|
||||
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
|
||||
crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
dPdFs, &
|
||||
Ps,dPdFs, &
|
||||
homogenization_typeInstance(material_homogenizationAt(el)))
|
||||
|
||||
end select chosenHomogenization
|
||||
|
|
@ -198,21 +203,17 @@ module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy)
|
|||
|
||||
integer :: co
|
||||
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
|
||||
real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
|
||||
real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
|
||||
|
||||
|
||||
if (homogenization_type(material_homogenizationAt(el)) == HOMOGENIZATION_RGC_ID) then
|
||||
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el)
|
||||
dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(subdt,co,ip,el)
|
||||
Fs(:,:,co) = constitutive_mech_getF(co,ip,el)
|
||||
Ps(:,:,co) = constitutive_mech_getP(co,ip,el)
|
||||
enddo
|
||||
doneAndHappy = &
|
||||
mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
crystallite_partitionedF0(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el),&
|
||||
subF,&
|
||||
subdt, &
|
||||
dPdFs, &
|
||||
ip, &
|
||||
el)
|
||||
doneAndHappy = mech_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ip,el)
|
||||
else
|
||||
doneAndHappy = .true.
|
||||
endif
|
||||
|
|
|
|||
|
|
@ -24,9 +24,6 @@ submodule(homogenization:homogenization_mech) homogenization_mech_RGC
|
|||
end type tParameters
|
||||
|
||||
type :: tRGCstate
|
||||
real(pReal), pointer, dimension(:) :: &
|
||||
work, &
|
||||
penaltyEnergy
|
||||
real(pReal), pointer, dimension(:,:) :: &
|
||||
relaxationVector
|
||||
end type tRGCstate
|
||||
|
|
@ -170,8 +167,7 @@ module subroutine mech_RGC_init(num_homogMech)
|
|||
nIntFaceTot = 3*( (prm%N_constituents(1)-1)*prm%N_constituents(2)*prm%N_constituents(3) &
|
||||
+ prm%N_constituents(1)*(prm%N_constituents(2)-1)*prm%N_constituents(3) &
|
||||
+ prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1))
|
||||
sizeState = nIntFaceTot &
|
||||
+ size(['avg constitutive work ','average penalty energy'])
|
||||
sizeState = nIntFaceTot
|
||||
|
||||
homogState(h)%sizeState = sizeState
|
||||
allocate(homogState(h)%state0 (sizeState,Nmaterialpoints), source=0.0_pReal)
|
||||
|
|
@ -180,8 +176,6 @@ module subroutine mech_RGC_init(num_homogMech)
|
|||
|
||||
stt%relaxationVector => homogState(h)%state(1:nIntFaceTot,:)
|
||||
st0%relaxationVector => homogState(h)%state0(1:nIntFaceTot,:)
|
||||
stt%work => homogState(h)%state(nIntFaceTot+1,:)
|
||||
stt%penaltyEnergy => homogState(h)%state(nIntFaceTot+2,:)
|
||||
|
||||
allocate(dst%volumeDiscrepancy( Nmaterialpoints), source=0.0_pReal)
|
||||
allocate(dst%relaxationRate_avg( Nmaterialpoints), source=0.0_pReal)
|
||||
|
|
@ -243,12 +237,11 @@ end subroutine mech_RGC_partitionDeformation
|
|||
!> @brief update the internal state of the homogenization scheme and tell whether "done" and
|
||||
! "happy" with result
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHappy)
|
||||
module function mech_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy)
|
||||
logical, dimension(2) :: doneAndHappy
|
||||
real(pReal), dimension(:,:,:), intent(in) :: &
|
||||
P,& !< partitioned stresses
|
||||
F,& !< partitioned deformation gradients
|
||||
F0 !< partitioned initial deformation gradients
|
||||
F !< partitioned deformation gradients
|
||||
real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
|
||||
real(pReal), dimension(3,3), intent(in) :: avgF !< average F
|
||||
real(pReal), intent(in) :: dt !< time increment
|
||||
|
|
@ -287,8 +280,8 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHa
|
|||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster
|
||||
allocate(resid(3*nIntFaceTot), source=0.0_pReal)
|
||||
allocate(tract(nIntFaceTot,3), source=0.0_pReal)
|
||||
allocate(resid(3*nIntFaceTot), source=0.0_pReal)
|
||||
allocate(tract(nIntFaceTot,3), source=0.0_pReal)
|
||||
relax = stt%relaxationVector(:,of)
|
||||
drelax = stt%relaxationVector(:,of) - st0%relaxationVector(:,of)
|
||||
|
||||
|
|
@ -346,17 +339,6 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHa
|
|||
if (residMax < num%rtol*stresMax .or. residMax < num%atol) then
|
||||
doneAndHappy = .true.
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! compute/update the state for postResult, i.e., all energy densities computed by time-integration
|
||||
do iGrain = 1,product(prm%N_constituents)
|
||||
do i = 1,3;do j = 1,3
|
||||
stt%work(of) = stt%work(of) &
|
||||
+ P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal)
|
||||
stt%penaltyEnergy(of) = stt%penaltyEnergy(of) &
|
||||
+ R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal)
|
||||
enddo; enddo
|
||||
enddo
|
||||
|
||||
dst%mismatch(1:3,of) = sum(NN,2)/real(nGrain,pReal)
|
||||
dst%relaxationRate_avg(of) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal)
|
||||
dst%relaxationRate_max(of) = maxval(abs(drelax))/dt
|
||||
|
|
@ -523,7 +505,6 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHa
|
|||
integer, dimension (3) :: iGrain3,iGNghb3,nGDim
|
||||
real(pReal), dimension (3,3) :: gDef,nDef
|
||||
real(pReal), dimension (3) :: nVect,surfCorr
|
||||
real(pReal), dimension (2) :: Gmoduli
|
||||
integer :: iGrain,iGNghb,iFace,i,j,k,l
|
||||
real(pReal) :: muGrain,muGNghb,nDefNorm
|
||||
real(pReal), parameter :: &
|
||||
|
|
@ -755,15 +736,9 @@ module subroutine mech_RGC_results(instance,group)
|
|||
associate(stt => state(instance), dst => dependentState(instance), prm => param(instance))
|
||||
outputsLoop: do o = 1,size(prm%output)
|
||||
select case(trim(prm%output(o)))
|
||||
case('W')
|
||||
call results_writeDataset(group,stt%work,trim(prm%output(o)), &
|
||||
'work density','J/m³')
|
||||
case('M')
|
||||
call results_writeDataset(group,dst%mismatch,trim(prm%output(o)), &
|
||||
'average mismatch tensor','1')
|
||||
case('R')
|
||||
call results_writeDataset(group,stt%penaltyEnergy,trim(prm%output(o)), &
|
||||
'mismatch penalty density','J/m³')
|
||||
case('Delta_V')
|
||||
call results_writeDataset(group,dst%volumeDiscrepancy,trim(prm%output(o)), &
|
||||
'volume discrepancy','m³')
|
||||
|
|
|
|||
|
|
@ -0,0 +1,39 @@
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
!> @author Martin Diehl, KU Leuven
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
submodule(homogenization) homogenization_thermal
|
||||
|
||||
|
||||
contains
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Allocate variables and set parameters.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module subroutine thermal_init()
|
||||
|
||||
|
||||
print'(/,a)', ' <<<+- homogenization_thermal init -+>>>'
|
||||
|
||||
allocate(homogenization_T(discretization_nIPs*discretization_Nelems))
|
||||
|
||||
|
||||
end subroutine thermal_init
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Partition T onto the individual constituents.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module subroutine thermal_partition(T,ip,el)
|
||||
|
||||
real(pReal), intent(in) :: T
|
||||
integer, intent(in) :: &
|
||||
ip, & !< integration point
|
||||
el !< element number
|
||||
|
||||
|
||||
call constitutive_thermal_setT(T,1,ip,el)
|
||||
|
||||
end subroutine thermal_partition
|
||||
|
||||
|
||||
end submodule homogenization_thermal
|
||||
|
|
@ -453,12 +453,13 @@ contains
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine lattice_init
|
||||
|
||||
integer :: Nphases, p,i
|
||||
integer :: Nphases, ph,i
|
||||
class(tNode), pointer :: &
|
||||
phases, &
|
||||
phase, &
|
||||
mech, &
|
||||
elasticity
|
||||
elasticity, &
|
||||
thermal
|
||||
|
||||
print'(/,a)', ' <<<+- lattice init -+>>>'; flush(IO_STDOUT)
|
||||
|
||||
|
|
@ -476,67 +477,71 @@ subroutine lattice_init
|
|||
lattice_mu, lattice_nu,&
|
||||
source=[(0.0_pReal,i=1,Nphases)])
|
||||
|
||||
do p = 1, phases%length
|
||||
phase => phases%get(p)
|
||||
do ph = 1, phases%length
|
||||
phase => phases%get(ph)
|
||||
mech => phase%get('mechanics')
|
||||
elasticity => mech%get('elasticity')
|
||||
lattice_C66(1,1,p) = elasticity%get_asFloat('C_11')
|
||||
lattice_C66(1,2,p) = elasticity%get_asFloat('C_12')
|
||||
lattice_C66(1,1,ph) = elasticity%get_asFloat('C_11')
|
||||
lattice_C66(1,2,ph) = elasticity%get_asFloat('C_12')
|
||||
|
||||
lattice_C66(1,3,p) = elasticity%get_asFloat('C_13',defaultVal=0.0_pReal)
|
||||
lattice_C66(2,2,p) = elasticity%get_asFloat('C_22',defaultVal=0.0_pReal)
|
||||
lattice_C66(2,3,p) = elasticity%get_asFloat('C_23',defaultVal=0.0_pReal)
|
||||
lattice_C66(3,3,p) = elasticity%get_asFloat('C_33',defaultVal=0.0_pReal)
|
||||
lattice_C66(4,4,p) = elasticity%get_asFloat('C_44',defaultVal=0.0_pReal)
|
||||
lattice_C66(5,5,p) = elasticity%get_asFloat('C_55',defaultVal=0.0_pReal)
|
||||
lattice_C66(6,6,p) = elasticity%get_asFloat('C_66',defaultVal=0.0_pReal)
|
||||
lattice_C66(1,3,ph) = elasticity%get_asFloat('C_13',defaultVal=0.0_pReal)
|
||||
lattice_C66(2,2,ph) = elasticity%get_asFloat('C_22',defaultVal=0.0_pReal)
|
||||
lattice_C66(2,3,ph) = elasticity%get_asFloat('C_23',defaultVal=0.0_pReal)
|
||||
lattice_C66(3,3,ph) = elasticity%get_asFloat('C_33',defaultVal=0.0_pReal)
|
||||
lattice_C66(4,4,ph) = elasticity%get_asFloat('C_44',defaultVal=0.0_pReal)
|
||||
lattice_C66(5,5,ph) = elasticity%get_asFloat('C_55',defaultVal=0.0_pReal)
|
||||
lattice_C66(6,6,ph) = elasticity%get_asFloat('C_66',defaultVal=0.0_pReal)
|
||||
|
||||
select case(phase%get_asString('lattice'))
|
||||
case('cF')
|
||||
lattice_structure(p) = lattice_FCC_ID
|
||||
lattice_structure(ph) = lattice_FCC_ID
|
||||
case('cI')
|
||||
lattice_structure(p) = lattice_BCC_ID
|
||||
lattice_structure(ph) = lattice_BCC_ID
|
||||
case('hP')
|
||||
lattice_structure(p) = lattice_HEX_ID
|
||||
lattice_structure(ph) = lattice_HEX_ID
|
||||
case('tI')
|
||||
lattice_structure(p) = lattice_BCT_ID
|
||||
lattice_structure(ph) = lattice_BCT_ID
|
||||
case('oP')
|
||||
lattice_structure(p) = lattice_ORT_ID
|
||||
lattice_structure(ph) = lattice_ORT_ID
|
||||
case('aP')
|
||||
lattice_structure(p) = lattice_ISO_ID
|
||||
lattice_structure(ph) = lattice_ISO_ID
|
||||
case default
|
||||
call IO_error(130,ext_msg='lattice_init: '//phase%get_asString('lattice'))
|
||||
end select
|
||||
|
||||
lattice_C66(1:6,1:6,p) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,p),phase%get_asString('lattice'))
|
||||
lattice_C66(1:6,1:6,ph) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,ph),phase%get_asString('lattice'))
|
||||
|
||||
lattice_nu(p) = lattice_equivalent_nu(lattice_C66(1:6,1:6,p),'voigt')
|
||||
lattice_mu(p) = lattice_equivalent_mu(lattice_C66(1:6,1:6,p),'voigt')
|
||||
lattice_nu(ph) = lattice_equivalent_nu(lattice_C66(1:6,1:6,ph),'voigt')
|
||||
lattice_mu(ph) = lattice_equivalent_mu(lattice_C66(1:6,1:6,ph),'voigt')
|
||||
|
||||
lattice_C66(1:6,1:6,p) = math_sym3333to66(math_Voigt66to3333(lattice_C66(1:6,1:6,p))) ! Literature data is in Voigt notation
|
||||
lattice_C66(1:6,1:6,ph) = math_sym3333to66(math_Voigt66to3333(lattice_C66(1:6,1:6,ph))) ! Literature data is in Voigt notation
|
||||
do i = 1, 6
|
||||
if (abs(lattice_C66(i,i,p))<tol_math_check) &
|
||||
call IO_error(135,el=i,ip=p,ext_msg='matrix diagonal "el"ement of phase "ip"')
|
||||
if (abs(lattice_C66(i,i,ph))<tol_math_check) &
|
||||
call IO_error(135,el=i,ip=ph,ext_msg='matrix diagonal "el"ement of phase "ip"')
|
||||
enddo
|
||||
|
||||
lattice_rho(ph) = phase%get_asFloat('rho', defaultVal=0.0_pReal)
|
||||
|
||||
! SHOULD NOT BE PART OF LATTICE BEGIN
|
||||
lattice_K(1,1,p) = phase%get_asFloat('K_11',defaultVal=0.0_pReal)
|
||||
lattice_K(2,2,p) = phase%get_asFloat('K_22',defaultVal=0.0_pReal)
|
||||
lattice_K(3,3,p) = phase%get_asFloat('K_33',defaultVal=0.0_pReal)
|
||||
lattice_K(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_K(1:3,1:3,p), &
|
||||
|
||||
if (phase%contains('thermal')) then
|
||||
thermal => phase%get('thermal')
|
||||
lattice_K(1,1,ph) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal)
|
||||
lattice_K(2,2,ph) = thermal%get_asFloat('K_22',defaultVal=0.0_pReal)
|
||||
lattice_K(3,3,ph) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal)
|
||||
lattice_K(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_K(1:3,1:3,ph), &
|
||||
phase%get_asString('lattice'))
|
||||
lattice_c_p(ph) = thermal%get_asFloat('c_p', defaultVal=0.0_pReal)
|
||||
endif
|
||||
|
||||
|
||||
lattice_D(1,1,ph) = phase%get_asFloat('D_11',defaultVal=0.0_pReal)
|
||||
lattice_D(2,2,ph) = phase%get_asFloat('D_22',defaultVal=0.0_pReal)
|
||||
lattice_D(3,3,ph) = phase%get_asFloat('D_33',defaultVal=0.0_pReal)
|
||||
lattice_D(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,ph), &
|
||||
phase%get_asString('lattice'))
|
||||
|
||||
lattice_c_p(p) = phase%get_asFloat('c_p', defaultVal=0.0_pReal)
|
||||
lattice_rho(p) = phase%get_asFloat('rho', defaultVal=0.0_pReal)
|
||||
|
||||
lattice_D(1,1,p) = phase%get_asFloat('D_11',defaultVal=0.0_pReal)
|
||||
lattice_D(2,2,p) = phase%get_asFloat('D_22',defaultVal=0.0_pReal)
|
||||
lattice_D(3,3,p) = phase%get_asFloat('D_33',defaultVal=0.0_pReal)
|
||||
lattice_D(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,p), &
|
||||
phase%get_asString('lattice'))
|
||||
|
||||
lattice_M(p) = phase%get_asFloat('M',defaultVal=0.0_pReal)
|
||||
lattice_M(ph) = phase%get_asFloat('M',defaultVal=0.0_pReal)
|
||||
! SHOULD NOT BE PART OF LATTICE END
|
||||
|
||||
call selfTest
|
||||
|
|
|
|||
|
|
@ -61,7 +61,7 @@ module material
|
|||
|
||||
type(tState), allocatable, dimension(:), public :: &
|
||||
homogState, &
|
||||
damageState
|
||||
damageState_h
|
||||
|
||||
type(Rotation), dimension(:,:,:), allocatable, public, protected :: &
|
||||
material_orientation0 !< initial orientation of each grain,IP,element
|
||||
|
|
@ -101,7 +101,7 @@ subroutine material_init(restart)
|
|||
|
||||
|
||||
allocate(homogState (size(material_name_homogenization)))
|
||||
allocate(damageState (size(material_name_homogenization)))
|
||||
allocate(damageState_h (size(material_name_homogenization)))
|
||||
|
||||
allocate(temperature (size(material_name_homogenization)))
|
||||
allocate(damage (size(material_name_homogenization)))
|
||||
|
|
|
|||
|
|
@ -101,9 +101,9 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
|
|||
if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
|
||||
|
||||
Nconstituents = count(material_phaseAt==p) * discretization_nIPs
|
||||
call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0)
|
||||
sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
|
||||
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'
|
||||
call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0)
|
||||
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
|
||||
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'
|
||||
|
||||
end associate
|
||||
|
||||
|
|
@ -146,7 +146,7 @@ module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el)
|
|||
damageOffset = material_homogenizationMemberAt(ip,el)
|
||||
|
||||
associate(prm => param(source_damage_anisoBrittle_instance(phase)))
|
||||
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal
|
||||
damageState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal
|
||||
do i = 1, prm%sum_N_cl
|
||||
traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i))
|
||||
traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i))
|
||||
|
|
@ -154,8 +154,8 @@ module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el)
|
|||
|
||||
traction_crit = prm%g_crit(i)*damage(homog)%p(damageOffset)**2.0_pReal
|
||||
|
||||
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
|
||||
= sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
|
||||
damageState(phase)%p(sourceOffset)%dotState(1,constituent) &
|
||||
= damageState(phase)%p(sourceOffset)%dotState(1,constituent) &
|
||||
+ prm%dot_o / prm%s_crit(i) &
|
||||
* ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + &
|
||||
(max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + &
|
||||
|
|
@ -185,7 +185,7 @@ module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, d
|
|||
|
||||
sourceOffset = source_damage_anisoBrittle_offset(phase)
|
||||
|
||||
dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
|
||||
dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent)
|
||||
|
||||
localphiDot = 1.0_pReal &
|
||||
+ dLocalphiDot_dPhi*phi
|
||||
|
|
@ -204,7 +204,7 @@ module subroutine source_damage_anisoBrittle_results(phase,group)
|
|||
integer :: o
|
||||
|
||||
associate(prm => param(source_damage_anisoBrittle_instance(phase)), &
|
||||
stt => sourceState(phase)%p(source_damage_anisoBrittle_offset(phase))%state)
|
||||
stt => damageState(phase)%p(source_damage_anisoBrittle_offset(phase))%state)
|
||||
outputsLoop: do o = 1,size(prm%output)
|
||||
select case(trim(prm%output(o)))
|
||||
case ('f_phi')
|
||||
|
|
|
|||
|
|
@ -87,9 +87,9 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
|
|||
if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
|
||||
|
||||
Nconstituents=count(material_phaseAt==p) * discretization_nIPs
|
||||
call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0)
|
||||
sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal)
|
||||
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'
|
||||
call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0)
|
||||
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal)
|
||||
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'
|
||||
|
||||
end associate
|
||||
|
||||
|
|
@ -128,7 +128,7 @@ module subroutine source_damage_anisoDuctile_dotState(co, ip, el)
|
|||
damageOffset = material_homogenizationMemberAt(ip,el)
|
||||
|
||||
associate(prm => param(source_damage_anisoDuctile_instance(phase)))
|
||||
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
|
||||
damageState(phase)%p(sourceOffset)%dotState(1,constituent) &
|
||||
= sum(plasticState(phase)%slipRate(:,constituent)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit)
|
||||
end associate
|
||||
|
||||
|
|
@ -154,7 +154,7 @@ module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, d
|
|||
|
||||
sourceOffset = source_damage_anisoDuctile_offset(phase)
|
||||
|
||||
dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
|
||||
dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent)
|
||||
|
||||
localphiDot = 1.0_pReal &
|
||||
+ dLocalphiDot_dPhi*phi
|
||||
|
|
@ -173,7 +173,7 @@ module subroutine source_damage_anisoDuctile_results(phase,group)
|
|||
integer :: o
|
||||
|
||||
associate(prm => param(source_damage_anisoDuctile_instance(phase)), &
|
||||
stt => sourceState(phase)%p(source_damage_anisoDuctile_offset(phase))%state)
|
||||
stt => damageState(phase)%p(source_damage_anisoDuctile_offset(phase))%state)
|
||||
outputsLoop: do o = 1,size(prm%output)
|
||||
select case(trim(prm%output(o)))
|
||||
case ('f_phi')
|
||||
|
|
|
|||
|
|
@ -74,9 +74,9 @@ module function source_damage_isoBrittle_init(source_length) result(mySources)
|
|||
if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit'
|
||||
|
||||
Nconstituents = count(material_phaseAt==p) * discretization_nIPs
|
||||
call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,1)
|
||||
sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
|
||||
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'
|
||||
call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,1)
|
||||
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
|
||||
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'
|
||||
|
||||
end associate
|
||||
|
||||
|
|
@ -124,13 +124,13 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, co, ip, el)
|
|||
strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%W_crit
|
||||
! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/prm%W_crit
|
||||
|
||||
if (strainenergy > sourceState(phase)%p(sourceOffset)%subState0(1,constituent)) then
|
||||
sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
|
||||
strainenergy - sourceState(phase)%p(sourceOffset)%state(1,constituent)
|
||||
if (strainenergy > damageState(phase)%p(sourceOffset)%subState0(1,constituent)) then
|
||||
damageState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
|
||||
strainenergy - damageState(phase)%p(sourceOffset)%state(1,constituent)
|
||||
else
|
||||
sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
|
||||
sourceState(phase)%p(sourceOffset)%subState0(1,constituent) - &
|
||||
sourceState(phase)%p(sourceOffset)%state(1,constituent)
|
||||
damageState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
|
||||
damageState(phase)%p(sourceOffset)%subState0(1,constituent) - &
|
||||
damageState(phase)%p(sourceOffset)%state(1,constituent)
|
||||
endif
|
||||
end associate
|
||||
|
||||
|
|
@ -158,8 +158,8 @@ module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLo
|
|||
|
||||
associate(prm => param(source_damage_isoBrittle_instance(phase)))
|
||||
localphiDot = 1.0_pReal &
|
||||
- phi*sourceState(phase)%p(sourceOffset)%state(1,constituent)
|
||||
dLocalphiDot_dPhi = - sourceState(phase)%p(sourceOffset)%state(1,constituent)
|
||||
- phi*damageState(phase)%p(sourceOffset)%state(1,constituent)
|
||||
dLocalphiDot_dPhi = - damageState(phase)%p(sourceOffset)%state(1,constituent)
|
||||
end associate
|
||||
|
||||
end subroutine source_damage_isoBrittle_getRateAndItsTangent
|
||||
|
|
@ -176,7 +176,7 @@ module subroutine source_damage_isoBrittle_results(phase,group)
|
|||
integer :: o
|
||||
|
||||
associate(prm => param(source_damage_isoBrittle_instance(phase)), &
|
||||
stt => sourceState(phase)%p(source_damage_isoBrittle_offset(phase))%state)
|
||||
stt => damageState(phase)%p(source_damage_isoBrittle_offset(phase))%state)
|
||||
outputsLoop: do o = 1,size(prm%output)
|
||||
select case(trim(prm%output(o)))
|
||||
case ('f_phi')
|
||||
|
|
|
|||
|
|
@ -78,9 +78,9 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
|
|||
if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
|
||||
|
||||
Nconstituents=count(material_phaseAt==p) * discretization_nIPs
|
||||
call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0)
|
||||
sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
|
||||
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'
|
||||
call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0)
|
||||
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
|
||||
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'
|
||||
|
||||
end associate
|
||||
|
||||
|
|
@ -119,7 +119,7 @@ module subroutine source_damage_isoDuctile_dotState(co, ip, el)
|
|||
damageOffset = material_homogenizationMemberAt(ip,el)
|
||||
|
||||
associate(prm => param(source_damage_isoDuctile_instance(phase)))
|
||||
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
|
||||
damageState(phase)%p(sourceOffset)%dotState(1,constituent) = &
|
||||
sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit
|
||||
end associate
|
||||
|
||||
|
|
@ -145,7 +145,7 @@ module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLo
|
|||
|
||||
sourceOffset = source_damage_isoDuctile_offset(phase)
|
||||
|
||||
dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
|
||||
dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent)
|
||||
|
||||
localphiDot = 1.0_pReal &
|
||||
+ dLocalphiDot_dPhi*phi
|
||||
|
|
@ -164,7 +164,7 @@ module subroutine source_damage_isoDuctile_results(phase,group)
|
|||
integer :: o
|
||||
|
||||
associate(prm => param(source_damage_isoDuctile_instance(phase)), &
|
||||
stt => sourceState(phase)%p(source_damage_isoDuctile_offset(phase))%state)
|
||||
stt => damageState(phase)%p(source_damage_isoDuctile_offset(phase))%state)
|
||||
outputsLoop: do o = 1,size(prm%output)
|
||||
select case(trim(prm%output(o)))
|
||||
case ('f_phi')
|
||||
|
|
|
|||
|
|
@ -10,6 +10,7 @@ module thermal_conduction
|
|||
use results
|
||||
use constitutive
|
||||
use YAML_types
|
||||
use discretization
|
||||
|
||||
implicit none
|
||||
private
|
||||
|
|
@ -24,7 +25,7 @@ module thermal_conduction
|
|||
|
||||
public :: &
|
||||
thermal_conduction_init, &
|
||||
thermal_conduction_getSourceAndItsTangent, &
|
||||
thermal_conduction_getSource, &
|
||||
thermal_conduction_getConductivity, &
|
||||
thermal_conduction_getSpecificHeat, &
|
||||
thermal_conduction_getMassDensity, &
|
||||
|
|
@ -38,25 +39,28 @@ contains
|
|||
!> @brief module initialization
|
||||
!> @details reads in material parameters, allocates arrays, and does sanity checks
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine thermal_conduction_init
|
||||
subroutine thermal_conduction_init(T)
|
||||
|
||||
integer :: Ninstances,Nmaterialpoints,h
|
||||
real(pReal), dimension(:), intent(inout) :: T
|
||||
|
||||
integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce
|
||||
class(tNode), pointer :: &
|
||||
material_homogenization, &
|
||||
homog, &
|
||||
homogThermal
|
||||
|
||||
|
||||
|
||||
print'(/,a)', ' <<<+- thermal_conduction init -+>>>'; flush(6)
|
||||
|
||||
Ninstances = count(thermal_type == THERMAL_conduction_ID)
|
||||
allocate(param(Ninstances))
|
||||
|
||||
material_homogenization => config_material%get('homogenization')
|
||||
do h = 1, size(material_name_homogenization)
|
||||
if (thermal_type(h) /= THERMAL_conduction_ID) cycle
|
||||
homog => material_homogenization%get(h)
|
||||
do ho = 1, size(material_name_homogenization)
|
||||
if (thermal_type(ho) /= THERMAL_conduction_ID) cycle
|
||||
homog => material_homogenization%get(ho)
|
||||
homogThermal => homog%get('thermal')
|
||||
associate(prm => param(thermal_typeInstance(h)))
|
||||
associate(prm => param(thermal_typeInstance(ho)))
|
||||
|
||||
#if defined (__GFORTRAN__)
|
||||
prm%output = output_asStrings(homogThermal)
|
||||
|
|
@ -64,21 +68,30 @@ subroutine thermal_conduction_init
|
|||
prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray)
|
||||
#endif
|
||||
|
||||
Nmaterialpoints=count(material_homogenizationAt==h)
|
||||
Nmaterialpoints=count(material_homogenizationAt==ho)
|
||||
|
||||
allocate (temperature (h)%p(Nmaterialpoints), source=thermal_initialT(h))
|
||||
allocate (temperatureRate(h)%p(Nmaterialpoints), source=0.0_pReal)
|
||||
allocate (temperature (ho)%p(Nmaterialpoints), source=thermal_initialT(ho))
|
||||
allocate (temperatureRate(ho)%p(Nmaterialpoints), source=0.0_pReal)
|
||||
|
||||
end associate
|
||||
enddo
|
||||
|
||||
ce = 0
|
||||
do el = 1, discretization_Nelems
|
||||
do ip = 1, discretization_nIPs
|
||||
ce = ce + 1
|
||||
ho = material_homogenizationAt(el)
|
||||
if (thermal_type(ho) == THERMAL_conduction_ID) T(ce) = thermal_initialT(ho)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end subroutine thermal_conduction_init
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief return heat generation rate
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
|
||||
subroutine thermal_conduction_getSource(Tdot, T,ip,el)
|
||||
|
||||
integer, intent(in) :: &
|
||||
ip, & !< integration point number
|
||||
|
|
@ -86,20 +99,17 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
|
|||
real(pReal), intent(in) :: &
|
||||
T
|
||||
real(pReal), intent(out) :: &
|
||||
Tdot, dTdot_dT
|
||||
integer :: &
|
||||
homog
|
||||
|
||||
Tdot = 0.0_pReal
|
||||
dTdot_dT = 0.0_pReal
|
||||
Tdot
|
||||
|
||||
homog = material_homogenizationAt(el)
|
||||
call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S,crystallite_Lp ,ip, el)
|
||||
integer :: &
|
||||
homog
|
||||
|
||||
homog = material_homogenizationAt(el)
|
||||
call constitutive_thermal_getRate(TDot, T,ip,el)
|
||||
|
||||
Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal)
|
||||
dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal)
|
||||
|
||||
end subroutine thermal_conduction_getSourceAndItsTangent
|
||||
end subroutine thermal_conduction_getSource
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -112,14 +122,16 @@ function thermal_conduction_getConductivity(ip,el)
|
|||
el !< element number
|
||||
real(pReal), dimension(3,3) :: &
|
||||
thermal_conduction_getConductivity
|
||||
|
||||
integer :: &
|
||||
grain
|
||||
co
|
||||
|
||||
|
||||
thermal_conduction_getConductivity = 0.0_pReal
|
||||
do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
|
||||
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
thermal_conduction_getConductivity = thermal_conduction_getConductivity + &
|
||||
crystallite_push33ToRef(grain,ip,el,lattice_K(:,:,material_phaseAt(grain,el)))
|
||||
crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el)))
|
||||
enddo
|
||||
|
||||
thermal_conduction_getConductivity = thermal_conduction_getConductivity &
|
||||
|
|
@ -138,14 +150,16 @@ function thermal_conduction_getSpecificHeat(ip,el)
|
|||
el !< element number
|
||||
real(pReal) :: &
|
||||
thermal_conduction_getSpecificHeat
|
||||
|
||||
integer :: &
|
||||
grain
|
||||
co
|
||||
|
||||
|
||||
thermal_conduction_getSpecificHeat = 0.0_pReal
|
||||
|
||||
do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
|
||||
+ lattice_c_p(material_phaseAt(grain,el))
|
||||
+ lattice_c_p(material_phaseAt(co,el))
|
||||
enddo
|
||||
|
||||
thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
|
||||
|
|
@ -164,15 +178,16 @@ function thermal_conduction_getMassDensity(ip,el)
|
|||
el !< element number
|
||||
real(pReal) :: &
|
||||
thermal_conduction_getMassDensity
|
||||
|
||||
integer :: &
|
||||
grain
|
||||
co
|
||||
|
||||
|
||||
thermal_conduction_getMassDensity = 0.0_pReal
|
||||
|
||||
|
||||
do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &
|
||||
+ lattice_rho(material_phaseAt(grain,el))
|
||||
+ lattice_rho(material_phaseAt(co,el))
|
||||
enddo
|
||||
|
||||
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &
|
||||
|
|
|
|||
|
|
@ -6,6 +6,7 @@ module thermal_isothermal
|
|||
use prec
|
||||
use config
|
||||
use material
|
||||
use discretization
|
||||
|
||||
implicit none
|
||||
public
|
||||
|
|
@ -15,22 +16,33 @@ contains
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief allocates fields, reads information from material configuration file
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine thermal_isothermal_init
|
||||
subroutine thermal_isothermal_init(T)
|
||||
|
||||
integer :: h,Nmaterialpoints
|
||||
real(pReal), dimension(:), intent(inout) :: T
|
||||
|
||||
integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce
|
||||
|
||||
print'(/,a)', ' <<<+- thermal_isothermal init -+>>>'; flush(6)
|
||||
|
||||
do h = 1, size(material_name_homogenization)
|
||||
if (thermal_type(h) /= THERMAL_isothermal_ID) cycle
|
||||
do ho = 1, size(thermal_type)
|
||||
if (thermal_type(ho) /= THERMAL_isothermal_ID) cycle
|
||||
|
||||
Nmaterialpoints = count(material_homogenizationAt == h)
|
||||
Nmaterialpoints = count(material_homogenizationAt == ho)
|
||||
|
||||
allocate(temperature (h)%p(Nmaterialpoints),source=thermal_initialT(h))
|
||||
allocate(temperatureRate(h)%p(Nmaterialpoints),source = 0.0_pReal)
|
||||
allocate(temperature (ho)%p(Nmaterialpoints),source=thermal_initialT(ho))
|
||||
allocate(temperatureRate(ho)%p(Nmaterialpoints),source = 0.0_pReal)
|
||||
|
||||
enddo
|
||||
|
||||
ce = 0
|
||||
do el = 1, discretization_Nelems
|
||||
do ip = 1, discretization_nIPs
|
||||
ce = ce + 1
|
||||
ho = material_homogenizationAt(el)
|
||||
if (thermal_type(ho) == THERMAL_isothermal_ID) T(ce) = thermal_initialT(ho)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end subroutine thermal_isothermal_init
|
||||
|
||||
end module thermal_isothermal
|
||||
|
|
|
|||
Loading…
Reference in New Issue