diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index bb9083efb..0de483d70 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -174,8 +174,12 @@ pure function kinematics_thermal_expansion_initialStrain(ipc, ip, el) offset = thermalMapping(homog)%p(ip,el) kinematics_thermal_expansion_initialStrain = & - (temperature(homog)%p(offset) - lattice_referenceTemperature(phase)) * & - lattice_thermalExpansion33(1:3,1:3,phase) + (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**1 / 1. * & + lattice_thermalExpansion33(1:3,1:3,1,phase) + & ! constant coefficient + (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**2 / 2. * & + lattice_thermalExpansion33(1:3,1:3,2,phase) + & ! linear coefficient + (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**3 / 3. * & + lattice_thermalExpansion33(1:3,1:3,3,phase) ! quadratic coefficient end function kinematics_thermal_expansion_initialStrain @@ -215,9 +219,16 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc, TDot = temperatureRate(homog)%p(offset) TRef = lattice_referenceTemperature(phase) - Li = TDot* & - lattice_thermalExpansion33(1:3,1:3,phase)/ & - (1.0_pReal + lattice_thermalExpansion33(1:3,1:3,phase)*(T - TRef)) + Li = TDot * ( & + lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**0 & ! constant coefficient + + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 & ! linear coefficient + + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**2 & ! quadratic coefficient + ) / & + (1.0_pReal \ + + lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**1 / 1. & + + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. & + + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**3 / 3. & + ) dLi_dTstar3333 = 0.0_pReal end subroutine kinematics_thermal_expansion_LiAndItsTangent diff --git a/src/lattice.f90 b/src/lattice.f90 index f46d0ef93..37393b82e 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1082,9 +1082,10 @@ module lattice lattice_nu, & lattice_trans_mu, & lattice_trans_nu + real(pReal), dimension(:,:,:,:), allocatable, public, protected :: & ! with higher-order parameters (e.g. temperature-dependent) + lattice_thermalExpansion33 real(pReal), dimension(:,:,:), allocatable, public, protected :: & lattice_thermalConductivity33, & - lattice_thermalExpansion33, & lattice_damageDiffusion33, & lattice_vacancyfluxDiffusion33, & lattice_vacancyfluxMobility33, & @@ -1380,8 +1381,8 @@ subroutine lattice_init allocate(lattice_C3333(3,3,3,3,Nphases), source=0.0_pReal) allocate(lattice_trans_C66(6,6,Nphases), source=0.0_pReal) allocate(lattice_trans_C3333(3,3,3,3,Nphases), source=0.0_pReal) + allocate(lattice_thermalExpansion33 (3,3,3,Nphases), source=0.0_pReal) ! constant, linear, quadratic coefficients allocate(lattice_thermalConductivity33 (3,3,Nphases), source=0.0_pReal) - allocate(lattice_thermalExpansion33 (3,3,Nphases), source=0.0_pReal) allocate(lattice_damageDiffusion33 (3,3,Nphases), source=0.0_pReal) allocate(lattice_vacancyfluxDiffusion33 (3,3,Nphases), source=0.0_pReal) allocate(lattice_vacancyfluxMobility33 (3,3,Nphases), source=0.0_pReal) @@ -1545,11 +1546,17 @@ subroutine lattice_init case ('thermal_conductivity33') lattice_thermalConductivity33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt) case ('thermal_expansion11') - lattice_thermalExpansion33(1,1,section) = IO_floatValue(line,chunkPos,2_pInt) + do i = 2_pInt, min(4,chunkPos(1)) ! read up to three parameters (constant, linear, quadratic with T) + lattice_thermalExpansion33(1,1,i-1_pInt,section) = IO_floatValue(line,chunkPos,i) + enddo case ('thermal_expansion22') - lattice_thermalExpansion33(2,2,section) = IO_floatValue(line,chunkPos,2_pInt) + do i = 2_pInt, min(4,chunkPos(1)) ! read up to three parameters (constant, linear, quadratic with T) + lattice_thermalExpansion33(2,2,i-1_pInt,section) = IO_floatValue(line,chunkPos,i) + enddo case ('thermal_expansion33') - lattice_thermalExpansion33(3,3,section) = IO_floatValue(line,chunkPos,2_pInt) + do i = 2_pInt, min(4,chunkPos(1)) ! read up to three parameters (constant, linear, quadratic with T) + lattice_thermalExpansion33(3,3,i-1_pInt,section) = IO_floatValue(line,chunkPos,i) + enddo case ('specific_heat') lattice_specificHeat(section) = IO_floatValue(line,chunkPos,2_pInt) case ('vacancyformationenergy') @@ -1756,22 +1763,24 @@ subroutine lattice_initializeStructure(myPhase,CoverA,CoverA_trans,a_fcc,a_bcc) end select end select - lattice_thermalConductivity33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_thermalConductivity33(1:3,1:3,myPhase)) - lattice_thermalExpansion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_thermalExpansion33(1:3,1:3,myPhase)) - lattice_DamageDiffusion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_DamageDiffusion33(1:3,1:3,myPhase)) - lattice_vacancyfluxDiffusion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_vacancyfluxDiffusion33(1:3,1:3,myPhase)) - lattice_vacancyfluxMobility33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_vacancyfluxMobility33(1:3,1:3,myPhase)) - lattice_PorosityDiffusion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_PorosityDiffusion33(1:3,1:3,myPhase)) + forall (i = 1_pInt:3_pInt) & + lattice_thermalExpansion33 (1:3,1:3,i,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& + lattice_thermalExpansion33 (1:3,1:3,i,myPhase)) + + lattice_thermalConductivity33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& + lattice_thermalConductivity33 (1:3,1:3,myPhase)) + lattice_DamageDiffusion33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& + lattice_DamageDiffusion33 (1:3,1:3,myPhase)) + lattice_vacancyfluxDiffusion33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& + lattice_vacancyfluxDiffusion33 (1:3,1:3,myPhase)) + lattice_vacancyfluxMobility33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& + lattice_vacancyfluxMobility33 (1:3,1:3,myPhase)) + lattice_PorosityDiffusion33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& + lattice_PorosityDiffusion33 (1:3,1:3,myPhase)) lattice_hydrogenfluxDiffusion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_hydrogenfluxDiffusion33(1:3,1:3,myPhase)) - lattice_hydrogenfluxMobility33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& - lattice_hydrogenfluxMobility33(1:3,1:3,myPhase)) + lattice_hydrogenfluxDiffusion33(1:3,1:3,myPhase)) + lattice_hydrogenfluxMobility33 (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& + lattice_hydrogenfluxMobility33 (1:3,1:3,myPhase)) select case(lattice_structure(myPhase)) !--------------------------------------------------------------------------------------------------