diff --git a/PRIVATE b/PRIVATE index 76bb51348..bc6de828c 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 76bb51348de75207d483d369628670e5ae51dca9 +Subproject commit bc6de828cc4ee9c941b37113ca49fcf51abd3512 diff --git a/python/damask/_colormap.py b/python/damask/_colormap.py index 40aab652c..3060e99e0 100644 --- a/python/damask/_colormap.py +++ b/python/damask/_colormap.py @@ -5,7 +5,6 @@ import colorsys from pathlib import Path from typing import Sequence, Union, TextIO - import numpy as np import matplotlib as mpl if os.name == 'posix' and 'DISPLAY' not in os.environ: @@ -17,9 +16,9 @@ from PIL import Image from . import util from . import Table -_eps = 216./24389. -_kappa = 24389./27. -_ref_white = np.array([.95047, 1.00000, 1.08883]) # Observer = 2, Illuminant = D65 +_EPS = 216./24389. +_KAPPA = 24389./27. +_REF_WHITE = np.array([.95047, 1.00000, 1.08883]) # Observer = 2, Illuminant = D65 # ToDo (if needed) # - support alpha channel (paraview/ASCII/input) @@ -522,10 +521,10 @@ class Colormap(mpl.colors.ListedColormap): f_z = (lab[0]+16.)/116. - lab[2]/200. return np.array([ - f_x**3. if f_x**3. > _eps else (116.*f_x-16.)/_kappa, - ((lab[0]+16.)/116.)**3 if lab[0]>_kappa*_eps else lab[0]/_kappa, - f_z**3. if f_z**3. > _eps else (116.*f_z-16.)/_kappa - ])*(ref_white if ref_white is not None else _ref_white) + f_x**3. if f_x**3. > _EPS else (116.*f_x-16.)/_KAPPA, + ((lab[0]+16.)/116.)**3 if lab[0]>_KAPPA*_EPS else lab[0]/_KAPPA, + f_z**3. if f_z**3. > _EPS else (116.*f_z-16.)/_KAPPA + ])*(ref_white if ref_white is not None else _REF_WHITE) @staticmethod def _xyz2lab(xyz: np.ndarray, ref_white: np.ndarray = None) -> np.ndarray: @@ -537,8 +536,8 @@ class Colormap(mpl.colors.ListedColormap): http://www.brucelindbloom.com/index.html?Eqn_Lab_to_XYZ.html """ - ref_white = ref_white if ref_white is not None else _ref_white - f = np.where(xyz/ref_white > _eps,(xyz/ref_white)**(1./3.),(_kappa*xyz/ref_white+16.)/116.) + ref_white = ref_white if ref_white is not None else _REF_WHITE + f = np.where(xyz/ref_white > _EPS,(xyz/ref_white)**(1./3.),(_KAPPA*xyz/ref_white+16.)/116.) return np.array([ 116.0 * f[1] - 16.0, diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index e1d53ca83..e67149dea 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -4,6 +4,7 @@ !> @details List of files needed by MSC.Marc !-------------------------------------------------------------------------------------------------- #include "parallelization.f90" +#include "constants.f90" #include "IO.f90" #include "YAML_types.f90" #include "YAML_parse.f90" diff --git a/src/constants.f90 b/src/constants.f90 new file mode 100644 index 000000000..dd26ce78c --- /dev/null +++ b/src/constants.f90 @@ -0,0 +1,15 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, KU Leuven +!> @brief physical constants +!-------------------------------------------------------------------------------------------------- +module constants + use prec + + implicit none + public + + real(pReal), parameter :: & + T_ROOM = 300.0_pReal, & !< Room temperature in K + K_B = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin + +end module constants diff --git a/src/phase.f90 b/src/phase.f90 index 22c35416b..214b0f5fa 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -5,6 +5,7 @@ !-------------------------------------------------------------------------------------------------- module phase use prec + use constants use math use rotations use IO diff --git a/src/phase_mechanical_eigen_thermalexpansion.f90 b/src/phase_mechanical_eigen_thermalexpansion.f90 index 1c08140ae..a84dd7797 100644 --- a/src/phase_mechanical_eigen_thermalexpansion.f90 +++ b/src/phase_mechanical_eigen_thermalexpansion.f90 @@ -58,14 +58,14 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics) associate(prm => param(kinematics_thermal_expansion_instance(p))) kinematic_type => kinematics%get(k) - prm%T_ref = kinematic_type%get_asFloat('T_ref', defaultVal=0.0_pReal) + prm%T_ref = kinematic_type%get_asFloat('T_ref', defaultVal=T_ROOM) prm%A(1,1,1) = kinematic_type%get_asFloat('A_11') - prm%A(1,1,2) = kinematic_type%get_asFloat('A_11,T',defaultVal=0.0_pReal) + prm%A(1,1,2) = kinematic_type%get_asFloat('A_11,T', defaultVal=0.0_pReal) prm%A(1,1,3) = kinematic_type%get_asFloat('A_11,T^2',defaultVal=0.0_pReal) if (any(phase_lattice(p) == ['hP','tI'])) then prm%A(3,3,1) = kinematic_type%get_asFloat('A_33') - prm%A(3,3,2) = kinematic_type%get_asFloat('A_33,T',defaultVal=0.0_pReal) + prm%A(3,3,2) = kinematic_type%get_asFloat('A_33,T', defaultVal=0.0_pReal) prm%A(3,3,3) = kinematic_type%get_asFloat('A_33,T^2',defaultVal=0.0_pReal) end if do i=1, size(prm%A,3) @@ -98,14 +98,14 @@ module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) associate(prm => param(kinematics_thermal_expansion_instance(ph))) Li = dot_T * ( & - prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient - + prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient - + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient + prm%A(1:3,1:3,1) & ! constant coefficient + + prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient + + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient ) / & (1.0_pReal & - + prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. & - + prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. & - + prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. & + + prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1.0_pReal & + + prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2.0_pReal & + + prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3.0_pReal & ) end associate dLi_dTstar = 0.0_pReal diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 1b39f51ee..45abf48c1 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -1,13 +1,15 @@ submodule(phase:mechanical) elastic type :: tParameters - real(pReal) :: & + real(pReal),dimension(3) :: & C_11 = 0.0_pReal, & C_12 = 0.0_pReal, & C_13 = 0.0_pReal, & C_33 = 0.0_pReal, & C_44 = 0.0_pReal, & C_66 = 0.0_pReal + real(pReal) :: & + T_ref end type tParameters type(tParameters), allocatable, dimension(:) :: param @@ -28,7 +30,7 @@ module subroutine elastic_init(phases) phase, & mech, & elastic - + logical :: thermal_active print'(/,1x,a)', '<<<+- phase:mechanical:elastic init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:elastic:Hooke init -+>>>' @@ -45,15 +47,35 @@ module subroutine elastic_init(phases) associate(prm => param(ph)) - prm%C_11 = elastic%get_asFloat('C_11') - prm%C_12 = elastic%get_asFloat('C_12') - prm%C_44 = elastic%get_asFloat('C_44') + prm%T_ref = elastic%get_asFloat('T_ref', defaultVal=T_ROOM) + prm%C_11(1) = elastic%get_asFloat('C_11') + prm%C_11(2) = elastic%get_asFloat('C_11,T', defaultVal=0.0_pReal) + prm%C_11(3) = elastic%get_asFloat('C_11,T^2',defaultVal=0.0_pReal) + + prm%C_12(1) = elastic%get_asFloat('C_12') + prm%C_12(2) = elastic%get_asFloat('C_12,T', defaultVal=0.0_pReal) + prm%C_12(3) = elastic%get_asFloat('C_12,T^2',defaultVal=0.0_pReal) + + prm%C_44(1) = elastic%get_asFloat('C_44') + prm%C_44(2) = elastic%get_asFloat('C_44,T', defaultVal=0.0_pReal) + prm%C_44(3) = elastic%get_asFloat('C_44,T^2',defaultVal=0.0_pReal) + if (any(phase_lattice(ph) == ['hP','tI'])) then - prm%C_13 = elastic%get_asFloat('C_13') - prm%C_33 = elastic%get_asFloat('C_33') + prm%C_13(1) = elastic%get_asFloat('C_13') + prm%C_13(2) = elastic%get_asFloat('C_13,T', defaultVal=0.0_pReal) + prm%C_13(3) = elastic%get_asFloat('C_13,T^2',defaultVal=0.0_pReal) + + prm%C_33(1) = elastic%get_asFloat('C_33') + prm%C_33(2) = elastic%get_asFloat('C_33,T', defaultVal=0.0_pReal) + prm%C_33(3) = elastic%get_asFloat('C_33,T^2',defaultVal=0.0_pReal) + end if + + if (phase_lattice(ph) == 'tI') then + prm%C_66(1) = elastic%get_asFloat('C_66') + prm%C_66(2) = elastic%get_asFloat('C_66,T', defaultVal=0.0_pReal) + prm%C_66(3) = elastic%get_asFloat('C_66,T^2',defaultVal=0.0_pReal) end if - if (phase_lattice(ph) == 'tI') prm%C_66 = elastic%get_asFloat('C_66') end associate end do @@ -69,21 +91,44 @@ module function elastic_C66(ph,en) result(C66) integer, intent(in) :: & ph, & en + real(pReal), dimension(6,6) :: C66 + real(pReal) :: T associate(prm => param(ph)) C66 = 0.0_pReal - C66(1,1) = prm%C_11 - C66(1,2) = prm%C_12 - C66(4,4) = prm%C_44 + T = thermal_T(ph,en) + + C66(1,1) = prm%C_11(1) & + + prm%C_11(2)*(T - prm%T_ref)**1 & + + prm%C_11(3)*(T - prm%T_ref)**2 + + C66(1,2) = prm%C_12(1) & + + prm%C_12(2)*(T - prm%T_ref)**1 & + + prm%C_12(3)*(T - prm%T_ref)**2 + + C66(4,4) = prm%C_44(1) & + + prm%C_44(2)*(T - prm%T_ref)**1 & + + prm%C_44(3)*(T - prm%T_ref)**2 + if (any(phase_lattice(ph) == ['hP','tI'])) then - C66(1,3) = prm%C_13 - C66(3,3) = prm%C_33 + C66(1,3) = prm%C_13(1) & + + prm%C_13(2)*(T - prm%T_ref)**1 & + + prm%C_13(3)*(T - prm%T_ref)**2 + + C66(3,3) = prm%C_33(1) & + + prm%C_33(2)*(T - prm%T_ref)**1 & + + prm%C_33(3)*(T - prm%T_ref)**2 + end if - if (phase_lattice(ph) == 'tI') C66(6,6) = prm%C_66 + if (phase_lattice(ph) == 'tI') then + C66(6,6) = prm%C_66(1) & + + prm%C_66(2)*(T - prm%T_ref)**1 & + + prm%C_66(3)*(T - prm%T_ref)**2 + end if C66 = lattice_symmetrize_C66(C66,phase_lattice(ph)) diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 102e009fe..1dab5dfd4 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -7,9 +7,6 @@ !-------------------------------------------------------------------------------------------------- submodule(phase:plastic) dislotungsten - real(pReal), parameter :: & - kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin - type :: tParameters real(pReal) :: & D = 1.0_pReal, & !< grain size @@ -344,7 +341,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, & 0.0_pReal, & prm%dipoleformation) - v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(kB*T))*prm%f_at/(2.0_pReal*PI*kB*T)) & + v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(K_B*T))*prm%f_at/(2.0_pReal*PI*K_B*T)) & * (1.0_pReal/(d_hat+prm%d_caron)) dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron) ! ToDo: Discuss with Franz: Stress dependency? end where @@ -475,7 +472,7 @@ pure subroutine kinetics(Mp,T,ph,en, & if (present(tau_pos_out)) tau_pos_out = tau_pos if (present(tau_neg_out)) tau_neg_out = tau_neg - associate(BoltzmannRatio => prm%Q_s/(kB*T), & + associate(BoltzmannRatio => prm%Q_s/(K_B*T), & b_rho_half => stt%rho_mob(:,en) * prm%b_sl * 0.5_pReal, & effectiveLength => dst%Lambda_sl(:,en) - prm%w) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 10522737c..034b3f503 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -9,9 +9,6 @@ !-------------------------------------------------------------------------------------------------- submodule(phase:plastic) dislotwin - real(pReal), parameter :: & - kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin - type :: tParameters real(pReal) :: & Q_cl = 1.0_pReal, & !< activation energy for dislocation climb @@ -31,7 +28,7 @@ submodule(phase:plastic) dislotwin delta_G = 1.0_pReal, & !< Free energy difference between austensite and martensite i_tr = 1.0_pReal, & !< adjustment parameter to calculate MFP for transformation h = 1.0_pReal, & !< Stack height of hex nucleus - T_ref = 0.0_pReal, & + T_ref = T_ROOM, & a_cI = 1.0_pReal, & a_cF = 1.0_pReal real(pReal), dimension(2) :: & @@ -597,7 +594,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) shearBandingContribution: if (dNeq0(prm%v_sb)) then - E_kB_T = prm%E_sb/(kB*T) + E_kB_T = prm%E_sb/(K_B*T) call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? do i = 1,6 @@ -694,8 +691,8 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) * (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (mu*prm%b_sl(i)), & 1.0_pReal, & prm%ExtendedDislocations) - v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) & - * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal) + v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(K_B*T)) & + * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(K_B*T)) - 1.0_pReal) dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) & / (d_hat-prm%d_caron(i)) @@ -907,7 +904,7 @@ pure subroutine kinetics_sl(Mp,T,ph,en, & significantStress: where(tau_eff > tol_math_check) stressRatio = tau_eff/prm%tau_0 StressRatio_p = stressRatio** prm%p - Q_kB_T = prm%Q_sl/(kB*T) + Q_kB_T = prm%Q_sl/(K_B*T) v_wait_inverse = exp(Q_kB_T*(1.0_pReal-StressRatio_p)** prm%q) & / prm%v_0 v_run_inverse = prm%B/(tau_eff*prm%b_sl) @@ -980,7 +977,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& (prm%L_tw*prm%b_sl(i))*& - (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,en)-tau(i)))) + (1.0_pReal-exp(-prm%V_cs/(K_B*T)*(dst%tau_r_tw(i,en)-tau(i)))) else Ndot0=0.0_pReal end if @@ -1049,7 +1046,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& (prm%L_tr*prm%b_sl(i))*& - (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,en)-tau(i)))) + (1.0_pReal-exp(-prm%V_cs/(K_B*T)*(dst%tau_r_tr(i,en)-tau(i)))) else Ndot0=0.0_pReal end if diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index f53a16042..09c5f4d0f 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -19,9 +19,6 @@ submodule(phase:plastic) nonlocal type(tGeometry), dimension(:), allocatable :: geom - real(pReal), parameter :: & - kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin - ! storage order of dislocation types integer, dimension(*), parameter :: & sgl = [1,2,3,4,5,6,7,8] !< signed (single) @@ -1094,9 +1091,9 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & ! thermally activated annihilation of edge dipoles by climb rhoDotThermalAnnihilation = 0.0_pReal - D_SD = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature)) ! eq. 3.53 + D_SD = prm%D_0 * exp(-prm%Q_cl / (K_B * Temperature)) ! eq. 3.53 v_climb = D_SD * mu * prm%V_at & - / (PI * (1.0_pReal-nu) * (dUpper(:,1) + dLower(:,1)) * kB * Temperature) ! eq. 3.54 + / (PI * (1.0_pReal-nu) * (dUpper(:,1) + dLower(:,1)) * K_B * Temperature) ! eq. 3.54 forall (s = 1:prm%sum_N_sl, dUpper(s,1) > dLower(s,1)) & rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * v_climb(s) / (dUpper(s,1) - dLower(s,1)), & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & @@ -1671,9 +1668,9 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T, activationEnergy_P = criticalStress_P * activationVolume_P tauRel_P = min(1.0_pReal, tauEff / criticalStress_P) tPeierls = 1.0_pReal / prm%nu_a & - * exp(activationEnergy_P / (kB * T) & + * exp(activationEnergy_P / (K_B * T) & * (1.0_pReal - tauRel_P**prm%p)**prm%q) - dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (kB * T) & + dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (K_B * T) & * (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) * tauRel_P**(prm%p-1.0_pReal), & 0.0_pReal, & tauEff < criticalStress_P) @@ -1685,8 +1682,8 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T, criticalStress_S = prm%Q_sol / activationVolume_S tauRel_S = min(1.0_pReal, tauEff / criticalStress_S) tSolidSolution = 1.0_pReal / prm%nu_a & - * exp(prm%Q_sol / (kB * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q) - dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (kB * T) & + * exp(prm%Q_sol / (K_B * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q) + dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (K_B * T) & * (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal)* tauRel_S**(prm%p-1.0_pReal), & 0.0_pReal, & tauEff < criticalStress_S)