diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 30c663545..6a26b703e 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -199,23 +199,12 @@ core: - release # Needs closer look -#Phenopowerlaw_singleSlip: -# stage: fortran -# script: Phenopowerlaw_singleSlip/test.py -# except: -# - master -# - release - -# Conversion to pytest ongoing -J2_plasticBehavior: - stage: fortran - script: - - module load $IntelMarc $HDF5Marc $MSC - - MSC_VERSION=2020 J2_plasticBehavior/test.py - except: - - master - - release - +# Phenopowerlaw_singleSlip: +# stage: fortran +# script: Phenopowerlaw_singleSlip/test.py +# except: +# - master +# - release ################################################################################################### SpectralRuntime: @@ -270,7 +259,9 @@ merge_into_master: - cd $DAMASKROOT - export TESTEDREV=$(git describe) # might be detached from development branch - echo $TESTEDREV > python/damask/VERSION - - git commit python/damask/VERSION -m "[skip ci] updated version information after successful test of $TESTEDREV" + - > + git diff-index --quiet HEAD || + git commit python/damask/VERSION -m "[skip ci] updated version information after successful test of $TESTEDREV" - export UPDATEDREV=$(git describe) # tested state + 1 commit - git checkout master - git merge $UPDATEDREV -s recursive -X ours # conflicts occur only for inconsistent state diff --git a/PRIVATE b/PRIVATE index 726a2c08e..d39905021 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 726a2c08efb67701531384ef9bc6352f4bacd75a +Subproject commit d399050216d814627edbe0ff1c05afe8f76c7b40 diff --git a/examples/config/Kinematics_Thermal_Expansion.config b/examples/config/Kinematics_Thermal_Expansion.config deleted file mode 100644 index 0cedff049..000000000 --- a/examples/config/Kinematics_Thermal_Expansion.config +++ /dev/null @@ -1,2 +0,0 @@ -(kinematics) thermal_expansion -thermal_expansion11 0.0000231 diff --git a/examples/config/Phase_Nonlocal_Aluminum.config b/examples/config/Phase_Nonlocal_Aluminum.config deleted file mode 100644 index 0d835d5aa..000000000 --- a/examples/config/Phase_Nonlocal_Aluminum.config +++ /dev/null @@ -1,63 +0,0 @@ -[Aluminum] -plasticity nonlocal -/nonlocal/ - -(output) rho_sgl_mob_edg_pos -(output) rho_sgl_imm_edg_pos -(output) rho_sgl_mob_edg_neg -(output) rho_sgl_imm_edg_neg -(output) rho_sgl_mob_scr_pos -(output) rho_sgl_imm_scr_pos -(output) rho_sgl_mob_scr_neg -(output) rho_sgl_imm_scr_neg -(output) rho_dip_edg -(output) rho_dip_scr -(output) rho_forest -(output) gamma -(output) tau_pass -(output) v_edg_pos -(output) v_edg_neg -(output) v_scr_pos -(output) v_scr_neg - -lattice_structure fcc -Nslip 12 # number of slip systems - -burgers 2.86e-10 # Burgers vector in m -rhoSglEdgePos0 0.25e10 # Initial positive edge single dislocation density in m/m**3 (per slip family) -rhoSglEdgeNeg0 0.25e10 # Initial negative edge single dislocation density in m/m**3 (per slip family) -rhoSglScrewPos0 0.25e10 # Initial positive screw single dislocation density in m/m**3 (per slip family) -rhoSglScrewNeg0 0.25e10 # Initial negative screw single dislocation density in m/m**3 (per slip family) -rhoDipEdge0 1e8 # Initial edge dipole dislocation density in m/m**3 (per slip family) -rhoDipScrew0 1e8 # Initial screw dipole dislocation density in m/m**3 (per slip family) -rhoSglScatter 0 # standard deviation of scatter in initial single dislocation density -#rhoSglRandom 1e12 # randomly distributed total dislocation density (sum over all slip systems and types) in m/m**3 -#rhoSglRandomBinning 1 # binning size of randomly distributed dislocations (number of dislocations per ip volume) -minimumDipoleHeightEdge 2e-9 # minimum distance for stable edge dipoles in m (per slip family) -minimumDipoleHeightScrew 2e-9 # minimum distance for stable screw dipoles in m (per slip family) -lambda0 80 # prefactor for mean free path -edgeMultiplication 0.1 # factor to which edges contribute to multiplication -atomicVolume 1.7e-29 # atomic volume in m**3 -selfdiffusionPrefactor 1e-4 # prefactor for self-diffusion coefficient in m**2/s -selfdiffusionEnergy 2.3e-19 # activation enthalpy for seld-diffusion in J -solidSolutionEnergy 2e-19 # activation energy of solid solution particles in J -solidSolutionConcentration 1e-5 # concentration of solid solution in parts per b^3 -solidSolutionSize 2 # size of solid solution obstacles in multiples of burgers vector length -peierlsStressEdge 1e5 # Peierls stress for edges in Pa (per slip family) -peierlsStressScrew 1e5 # Peierls stress for screws in Pa (per slip family) -doublekinkWidth 10 # width of double kinks in multiples of burgers vector length b -viscosity 1e-4 # viscosity for dislocation glide in Pa s -p 1 # exponent for thermal barrier profile -q 1 # exponent for thermal barrier profile -attackFrequency 50e9 # attack frequency in Hz -surfaceTransmissivity 1.0 # transmissivity of free surfaces for dislocation flux -grainboundaryTransmissivity 0.0 # transmissivity of grain boundaries for dislocation flux (grain bundaries are identified as interfaces with different textures on both sides); if not set or set to negative number, the subroutine automatically determines the transmissivity at the grain boundary -interaction_SlipSlip 0 0 0.625 0.07 0.137 0.137 0.122 # Dislocation interaction coefficient -linetension 0.8 # constant indicating the effect of the line tension on the hardening coefficients (0 to 1) -edgejog 1.0 # fraction of annihilated screw dipoles that forms edge jogs (0 to 1) -shortRangeStressCorrection 0 # switch for use of short range correction stress -cutoffRadius 1e-3 # cutoff radius for dislocation stress in m -CFLfactor 2.0 # safety factor for CFL flux check (numerical parameter) -significantRho 1e6 # minimum dislocation density considered relevant in m/m**3 -#significantN 0.1 # minimum dislocation number per ip considered relevant -randomMultiplication 0 # switch for probabilistic extension of multiplication rate diff --git a/examples/config/Phase_Nonlocal_Nickel.config b/examples/config/Phase_Nonlocal_Nickel.config deleted file mode 100644 index 7bc2e1581..000000000 --- a/examples/config/Phase_Nonlocal_Nickel.config +++ /dev/null @@ -1,62 +0,0 @@ -[Ni_nonlocal] - -elasticity hooke -plasticity nonlocal -/nonlocal/ -(output) rho_sgl_mob_edg_pos -(output) rho_sgl_imm_edg_pos -(output) rho_sgl_mob_edg_neg -(output) rho_sgl_imm_edg_neg -(output) rho_sgl_mob_scr_pos -(output) rho_sgl_imm_scr_pos -(output) rho_sgl_mob_scr_neg -(output) rho_sgl_imm_scr_neg -(output) rho_dip_edg -(output) rho_dip_scr -(output) rho_forest -(output) gamma -(output) tau_pass -(output) v_edg_pos -(output) v_edg_neg -(output) v_scr_pos -(output) v_scr_neg - - -lattice_structure fcc -Nslip 12 # number of slip systems per family -burgers 2.48e-10 # Burgers vector in m -rhoSglEdgePos0 6e10 # Initial positive edge single dislocation density in m/m**3 -rhoSglEdgeNeg0 6e10 # Initial negative edge single dislocation density in m/m**3 -rhoSglScrewPos0 6e10 # Initial positive screw single dislocation density in m/m**3 -rhoSglScrewNeg0 6e10 # Initial negative screw single dislocation density in m/m**3 -rhoDipEdge0 0 # Initial edge dipole dislocation density in m/m**3 -rhoDipScrew0 0 # Initial screw dipole dislocation density in m/m**3 -rhoSglScatter 0 -minimumDipoleHeightEdge 2.6e-9 # 3.0e-9 # minimum distance for stable edge dipoles in m -minimumDipoleHeightScrew 12.0e-9 # 50e-9 # minimum distance for stable screw dipoles in m -lambda0 45 # 33 # prefactor for mean free path -edgeMultiplication 0.1 -randomMultiplication 0 -atomicVolume 1.2e-29 -selfdiffusionPrefactor 1.9e-4 # Gottstein p.168 # prefactor for self-diffusion coefficient -selfdiffusionEnergy 5.1e-19 # Gottstein p.168 # activation energy self-diffusion -solidSolutionEnergy 1.8e-19 # activation energy of solid solution particles in J -solidSolutionConcentration 5e-7 # 1e-7 -solidSolutionSize 1.0 -peierlsStressEdge 1e5 # Peierls stress for edges in Pa (per slip family) -peierlsStressScrew 1e5 # Peierls stress for screws in Pa (per slip family) -doublekinkWidth 10 # width of double kinks in multiples of burgers vector length b -viscosity 1e-3 # viscosity for dislocation glide in Pa s -p 1 # exponent for thermal barrier profile -q 1 # exponent for thermal barrier profile -attackFrequency 50e9 # attack frequency in Hz -surfaceTransmissivity 1.0 # transmissivity of free surfaces for dislocation flux -grainBoundaryTransmissivity 0.0 -significantRho 1e8 # dislocation density considered relevant in m/m**3 -significantN 1 -shortRangeStressCorrection 0 -CFLfactor 1.1 # safety factor for CFL flux check (numerical parameter) -r 1 -interaction_SlipSlip 0 0 0.625 0.07 0.137 0.137 0.122 # Dislocation interaction coefficient -linetension 0.8 -edgejog 0.01 # 0.2 diff --git a/examples/config/Phase_Phenopowerlaw_Aluminum.yaml b/examples/config/Phase_Phenopowerlaw_Aluminum.yaml deleted file mode 100644 index 2a4afb219..000000000 --- a/examples/config/Phase_Phenopowerlaw_Aluminum.yaml +++ /dev/null @@ -1,16 +0,0 @@ -Aluminum: - lattice: cF - mechanical: - output: [F, P, F_e, F_p, L_p, O] - elastic: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: Hooke} - plastic: - N_sl: [12] - a_sl: 2.25 - dot_gamma_0_sl: 0.001 - h_0_sl-sl: 75e6 - h_sl-sl: [1, 1, 1.4, 1.4, 1.4, 1.4, 1.4] - n_sl: 20 - output: [xi_sl, gamma_sl] - type: phenopowerlaw - xi_0_sl: [31e6] - xi_inf_sl: [63e6] diff --git a/examples/config/Phase_Phenopowerlaw_BCC-Ferrite.yaml b/examples/config/Phase_Phenopowerlaw_BCC-Ferrite.yaml deleted file mode 100644 index 5b714af7c..000000000 --- a/examples/config/Phase_Phenopowerlaw_BCC-Ferrite.yaml +++ /dev/null @@ -1,17 +0,0 @@ -# Tasan et.al. 2015 Acta Materalia -# Tasan et.al. 2015 International Journal of Plasticity -# Diehl et.al. 2015 Meccanica -Ferrite: - lattice: cI - mechanical: - elastic: {C_11: 233.3e9, C_12: 135.5e9, C_44: 118.0e9, type: Hooke} - plastic: - N_sl: [12, 12] - a_sl: 2.0 - dot_gamma_0_sl: 0.001 - h_0_sl-sl: 1000.0e6 - h_sl-sl: [1, 1.4, 1, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4] - n_sl: 20 - type: phenopowerlaw - xi_0_sl: [95.e6, 96.e6] - xi_inf_sl: [222.e6, 412.7e6] diff --git a/examples/config/homogenization/mechanical/RGC_8grains.yaml b/examples/config/homogenization/mechanical/RGC_8grains.yaml index 9200c364d..9b70e3645 100644 --- a/examples/config/homogenization/mechanical/RGC_8grains.yaml +++ b/examples/config/homogenization/mechanical/RGC_8grains.yaml @@ -6,5 +6,5 @@ a_g: [0.0, 0.0, 0.0] c_alpha: 2.0 cluster_size: [2, 2, 2] - output: [M, Delta_V, avg_a_dot, max_a_dot] + output: [M, Delta_V, avg_dot_a, max_dot_a] xi_alpha: 10.0 diff --git a/examples/config/phase/mechanical/eigen/thermalexpansion_Al.yaml b/examples/config/phase/mechanical/eigen/thermalexpansion_Al.yaml index 1876141f6..aa132609e 100644 --- a/examples/config/phase/mechanical/eigen/thermalexpansion_Al.yaml +++ b/examples/config/phase/mechanical/eigen/thermalexpansion_Al.yaml @@ -1,5 +1,5 @@ type: thermalexpansion references: - en.wikipedia.org/wiki/Thermal_expansion -A_11: [23.1e-6] +A_11: 23.1e-6 T_ref: 293.15 diff --git a/examples/config/phase/mechanical/eigen/thermalexpansion_Au.yaml b/examples/config/phase/mechanical/eigen/thermalexpansion_Au.yaml index 5d8030e1e..34e71a1bc 100644 --- a/examples/config/phase/mechanical/eigen/thermalexpansion_Au.yaml +++ b/examples/config/phase/mechanical/eigen/thermalexpansion_Au.yaml @@ -1,5 +1,5 @@ type: thermalexpansion references: - en.wikipedia.org/wiki/Thermal_expansion -A_11: [14e-6] +A_11: 14e-6 T_ref: 293.15 diff --git a/examples/config/phase/mechanical/eigen/thermalexpansion_C35E.yaml b/examples/config/phase/mechanical/eigen/thermalexpansion_C35E.yaml index dea09aa43..64d772512 100644 --- a/examples/config/phase/mechanical/eigen/thermalexpansion_C35E.yaml +++ b/examples/config/phase/mechanical/eigen/thermalexpansion_C35E.yaml @@ -1,5 +1,7 @@ type: thermalexpansion references: - en.wikipedia.org/wiki/Thermal_expansion, fitted from image description -A_11: [12.70371e-6, 7.54e-9, -1.0e-11] +A_11: 12.70371e-6 +A_11,T: 7.54e-9 +A_11,T^2: -1.0e-11 T_ref: 273.0 diff --git a/examples/config/phase/mechanical/eigen/thermalexpansion_Cu.yaml b/examples/config/phase/mechanical/eigen/thermalexpansion_Cu.yaml index 3eb41aa90..ea09c7fb3 100644 --- a/examples/config/phase/mechanical/eigen/thermalexpansion_Cu.yaml +++ b/examples/config/phase/mechanical/eigen/thermalexpansion_Cu.yaml @@ -1,5 +1,5 @@ type: thermalexpansion references: - en.wikipedia.org/wiki/Thermal_expansion -A_11: [17e-6] +A_11: 17e-6 T_ref: 293.15 diff --git a/examples/config/phase/mechanical/eigen/thermalexpansion_Fe.yaml b/examples/config/phase/mechanical/eigen/thermalexpansion_Fe.yaml index f36250490..b104e36f5 100644 --- a/examples/config/phase/mechanical/eigen/thermalexpansion_Fe.yaml +++ b/examples/config/phase/mechanical/eigen/thermalexpansion_Fe.yaml @@ -1,5 +1,5 @@ type: thermalexpansion references: - en.wikipedia.org/wiki/Thermal_expansion -A_11: [11.8e-6] +A_11: 11.8e-6 T_ref: 293.15 diff --git a/examples/config/phase/mechanical/eigen/thermalexpansion_W.yaml b/examples/config/phase/mechanical/eigen/thermalexpansion_W.yaml index a6d069bd3..d223cdd7d 100644 --- a/examples/config/phase/mechanical/eigen/thermalexpansion_W.yaml +++ b/examples/config/phase/mechanical/eigen/thermalexpansion_W.yaml @@ -1,5 +1,5 @@ type: thermalexpansion references: - en.wikipedia.org/wiki/Thermal_expansion -A_11: [4.5e-6] +A_11: 4.5e-6 T_ref: 293.15 diff --git a/examples/config/phase/mechanical/eigen/thermalexpansion_X20Cr13.yaml b/examples/config/phase/mechanical/eigen/thermalexpansion_X20Cr13.yaml index 7842b9d6d..f7c55d485 100644 --- a/examples/config/phase/mechanical/eigen/thermalexpansion_X20Cr13.yaml +++ b/examples/config/phase/mechanical/eigen/thermalexpansion_X20Cr13.yaml @@ -1,5 +1,6 @@ type: thermalexpansion references: - en.wikipedia.org/wiki/Thermal_expansion, fitted from image description -A_11: [11.365e-6, 5.0e-9] +A_11: 11.365e-6 +A_11,T: 5.0e-9 T_ref: 273.0 diff --git a/examples/config/phase/mechanical/elastic/Hooke_Au.yaml b/examples/config/phase/mechanical/elastic/Hooke_Au.yaml index 9a4da774d..5683bfb24 100644 --- a/examples/config/phase/mechanical/elastic/Hooke_Au.yaml +++ b/examples/config/phase/mechanical/elastic/Hooke_Au.yaml @@ -2,7 +2,8 @@ type: Hooke references: - J.P. Hirth and J. Lothe, Theory of Dislocations, 1982, - John Wiley & Sons + John Wiley & Sons, + page 837 C_11: 186e9 C_12: 157e9 C_44: 42e9 diff --git a/examples/config/phase/mechanical/elastic/Hooke_Fe.yaml b/examples/config/phase/mechanical/elastic/Hooke_Fe.yaml index 0e5d7db5c..b2912035f 100644 --- a/examples/config/phase/mechanical/elastic/Hooke_Fe.yaml +++ b/examples/config/phase/mechanical/elastic/Hooke_Fe.yaml @@ -2,7 +2,8 @@ type: Hooke references: - J.P. Hirth and J. Lothe, Theory of Dislocations, 1982, - John Wiley & Sons + John Wiley & Sons, + page 837 C_11: 242e9 C_12: 146.5e9 -C_44: 11.2e9 +C_44: 112e9 diff --git a/examples/config/phase/mechanical/elastic/Hooke_Ni.yaml b/examples/config/phase/mechanical/elastic/Hooke_Ni.yaml index c76632cdd..15b684e4c 100644 --- a/examples/config/phase/mechanical/elastic/Hooke_Ni.yaml +++ b/examples/config/phase/mechanical/elastic/Hooke_Ni.yaml @@ -2,7 +2,8 @@ type: Hooke references: - J.P. Hirth and J. Lothe, Theory of Dislocations, 1982, - John Wiley & Sons + John Wiley & Sons, + page 837 C_11: 246.5e9 C_12: 147.3e9 C_44: 124.7e9 diff --git a/examples/config/phase/mechanical/plastic/nonlocal_Al.yaml b/examples/config/phase/mechanical/plastic/nonlocal_Al.yaml new file mode 100644 index 000000000..a4a4babb1 --- /dev/null +++ b/examples/config/phase/mechanical/plastic/nonlocal_Al.yaml @@ -0,0 +1,49 @@ +type: nonlocal +references: + C. Kords, + On the role of dislocation transport in the constitutive description of crystal plasticity, + RWTH Aachen 2013 +output: [rho_u_ed_pos, rho_b_ed_pos, rho_u_ed_neg, rho_b_ed_neg, rho_u_sc_pos, rho_b_sc_pos, rho_u_sc_neg, rho_b_sc_neg, rho_d_ed, rho_d_sc] +N_sl: [12] + +b_sl: [2.86e-10] +V_at: 0.017e-27 # Omega +d_ed: [1.6e-9] +d_sc: [10.e-9] +i_sl: [60] # k_2 (lambda_0 in Tab. 7.1) +f_ed_mult: 0.1 # k_1 + +rho_u_ed_neg_0: [1.25e9] # 6e10 / (12*4) +rho_u_ed_pos_0: [1.25e9] # 6e10 / (12*4) +rho_u_sc_neg_0: [1.25e9] # 6e10 / (12*4) +rho_u_sc_pos_0: [1.25e9] # 6e10 / (12*4) +rho_d_ed_0: [0] +rho_d_sc_0: [0] + +D_0: 7.e-29 +Q_cl: 0.0 # no temperature dependency +Q_sol: 2.00272e-19 # 1.25 eV +c_sol: 1.5e-6 # correct unit? +f_sol: 2.0 # d_obst in multiples of b +tau_Peierls_ed: [.1e6] +tau_Peierls_sc: [.1e6] +w: 10 # w_k in multiple of b + +p_sl: 1 +q_sl: 1 + +nu_a: 50.e9 +B: 1.e-2 +f_ed: 1.0 # k_3 + +h_sl-sl: [0, 0, 0.625, 0.07, 0.137, 0.137, 0.122] # Table 3.4 + +chi_GB: 0.0 # full blocking at GB +chi_surface: 1.0 # no blocking at surface + +f_F: 0.0 # no line tension correction +sigma_rho_u: 0 # no random distribution + + +short_range_stress_correction: false +rho_significant: 1e6 diff --git a/examples/config/phase/mechanical/plastic/nonlocal_Ni.yaml b/examples/config/phase/mechanical/plastic/nonlocal_Ni.yaml new file mode 100644 index 000000000..81b4becd1 --- /dev/null +++ b/examples/config/phase/mechanical/plastic/nonlocal_Ni.yaml @@ -0,0 +1,49 @@ +type: nonlocal +references: + C. Kords, + On the role of dislocation transport in the constitutive description of crystal plasticity, + RWTH Aachen 2013 +output: [rho_u_ed_pos, rho_b_ed_pos, rho_u_ed_neg, rho_b_ed_neg, rho_u_sc_pos, rho_b_sc_pos, rho_u_sc_neg, rho_b_sc_neg, rho_d_ed, rho_d_sc] +N_sl: [12] + +b_sl: [2.48e-10] +V_at: 0.012e-27 # Omega +d_ed: [2.6e-9] +d_sc: [12.e-9] +i_sl: [45] # k_2 +f_ed_mult: 0.1 # k_1 + +rho_u_ed_neg_0: [6.e10] # 2.88e12 / (12*4) +rho_u_ed_pos_0: [6.e10] # 2.88e12 / (12*4) +rho_u_sc_neg_0: [6.e10] # 2.88e12 / (12*4) +rho_u_sc_pos_0: [6.e10] # 2.88e12 / (12*4) +rho_d_ed_0: [0] +rho_d_sc_0: [0] + +D_0: 3.e-53 +Q_cl: 0.0 # no temperature dependency +Q_sol: 1.79444e-19 # 1.12 eV +c_sol: 5.e-7 # correct unit? +f_sol: 1. # d_obst +tau_Peierls_ed: [.1e6] +tau_Peierls_sc: [.1e6] +w: 10 # w_k + +p_sl: 1 +q_sl: 1 + +nu_a: 50.e9 +B: 1.e-3 +f_ed: 0.01 # k_3 + +h_sl-sl: [0, 0, 0.625, 0.07, 0.137, 0.137, 0.122] # Table 3.4 + +chi_GB: 0.0 # full blocking at GB +chi_surface: 1.0 # no blocking at surface + +f_F: 0.0 # no line tension correction +sigma_rho_u: 0 # no random distribution + + +short_range_stress_correction: false +rho_significant: 1e6 diff --git a/examples/config/phase/mechanical/plastic/phenopowerlaw_Al.yaml b/examples/config/phase/mechanical/plastic/phenopowerlaw_Al.yaml new file mode 100644 index 000000000..4e4bcf588 --- /dev/null +++ b/examples/config/phase/mechanical/plastic/phenopowerlaw_Al.yaml @@ -0,0 +1,15 @@ +type: phenopowerlaw +references: + - W.F. Hosford et al., + Acta Metallurgica, 8(3), 187-199, 1960, + 10.1016/0001-6160(60)90127-9, + fitted from Fig. 5 +output: [xi_sl, gamma_sl] +N_sl: [12] +n_sl: 20 +a_sl: 3.1 +h_0_sl-sl: 1.7e8 +xi_0_sl: [5.0e6] +xi_inf_sl: [37.5e6] +h_sl-sl: [1, 1, 1.4, 1.4, 1.4, 1.4, 1.4] +dot_gamma_0_sl: 4.5e-3 diff --git a/examples/config/phase/mechanical/plastic/phenopowerlaw_Au.yaml b/examples/config/phase/mechanical/plastic/phenopowerlaw_Au.yaml index c5e63e8a1..cb2e6fea1 100644 --- a/examples/config/phase/mechanical/plastic/phenopowerlaw_Au.yaml +++ b/examples/config/phase/mechanical/plastic/phenopowerlaw_Au.yaml @@ -9,9 +9,9 @@ references: output: [xi_sl, gamma_sl] N_sl: [12] n_sl: 83.3 -dot_gamma_0_sl: 0.001 -h_0_sl-sl: 75.0e6 -h_sl-sl: [1, 1, 1.4, 1.4, 1.4, 1.4, 1.4] a_sl: 1.0 +h_0_sl-sl: 75.0e6 xi_0_sl: [26.25e6] xi_inf_sl: [53.0e6] +h_sl-sl: [1, 1, 1.4, 1.4, 1.4, 1.4, 1.4] +dot_gamma_0_sl: 0.001 diff --git a/examples/config/phase/mechanical/plastic/phenopowerlaw_Cu.yaml b/examples/config/phase/mechanical/plastic/phenopowerlaw_Cu.yaml new file mode 100644 index 000000000..528f55763 --- /dev/null +++ b/examples/config/phase/mechanical/plastic/phenopowerlaw_Cu.yaml @@ -0,0 +1,15 @@ +type: phenopowerlaw +references: + - T Takeuchi, + Transactions of the Japan Institute of Metals 16(10), 629-640, 1975, + 10.2320/matertrans1960.16.629, + fitted from Fig. 3b +output: [xi_sl, gamma_sl] +N_sl: [12] +n_sl: 20 +a_sl: 1.0 +h_0_sl-sl: 2.4e8 +xi_0_sl: [1.5e6] +xi_inf_sl: [112.5e6] +h_sl-sl: [1, 1, 1.4, 1.4, 1.4, 1.4, 1.4] +dot_gamma_0_sl: 3e-3 diff --git a/examples/config/phase/mechanical/plastic/phenopowerlaw_DP-steel-ferrite.yaml b/examples/config/phase/mechanical/plastic/phenopowerlaw_DP-steel-ferrite.yaml new file mode 100644 index 000000000..84c0ecdc9 --- /dev/null +++ b/examples/config/phase/mechanical/plastic/phenopowerlaw_DP-steel-ferrite.yaml @@ -0,0 +1,14 @@ +type: phenopowerlaw +references: + - C.C. Tasan et al., + Acta Materialia, 81, 386-400, 2014, + 10.1016/j.actamat.2014.07.071 +output: [xi_sl, gamma_sl] +N_sl: [12, 12] +n_sl: 20 +a_sl: 2.25 +h_0_sl-sl: 1.0e9 +xi_0_sl: [95.e6, 96.e6] +xi_inf_sl: [222.e6, 412.e6] +h_sl-sl: [1, 1.4, 1, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4] +dot_gamma_0_sl: 0.001 diff --git a/python/damask/VERSION b/python/damask/VERSION index b1f60bf84..d9f554441 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha3-278-g23b9361bf +v3.0.0-alpha4-1-gb567194fe diff --git a/python/tests/reference/Result/4grains2x4x3.material.yaml b/python/tests/reference/Result/4grains2x4x3.material.yaml index 160db3ace..aebdc0338 100644 --- a/python/tests/reference/Result/4grains2x4x3.material.yaml +++ b/python/tests/reference/Result/4grains2x4x3.material.yaml @@ -633,7 +633,7 @@ homogenization: a_g: [0.0, 0.0, 0.0] c_alpha: 2.0 cluster_size: [2, 2, 2] - output: [M, Delta_V, avg_a_dot, max_a_dot] + output: [M, Delta_V, avg_dot_a, max_dot_a] xi_alpha: 10.0 phase: diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index b0bb96402..5fe754ce3 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -32,14 +32,12 @@ #include "phase_mechanical_plastic_nonlocal.f90" #include "phase_mechanical_eigen.f90" #include "phase_mechanical_eigen_cleavageopening.f90" -#include "phase_mechanical_eigen_slipplaneopening.f90" #include "phase_mechanical_eigen_thermalexpansion.f90" #include "phase_thermal.f90" #include "phase_thermal_dissipation.f90" #include "phase_thermal_externalheat.f90" #include "phase_damage.f90" #include "phase_damage_isobrittle.f90" -#include "phase_damage_isoductile.f90" #include "phase_damage_anisobrittle.f90" #include "homogenization.f90" #include "homogenization_mechanical.f90" diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index ff985eb6a..312c79aa7 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -258,7 +258,6 @@ subroutine readVTI(grid,geomSize,origin,material) character(len=:), allocatable :: temp real(pReal), dimension(:), allocatable :: delta - integer, dimension(:), allocatable :: stringPos integer :: i @@ -266,9 +265,9 @@ subroutine readVTI(grid,geomSize,origin,material) call IO_error(error_ID = 844, ext_msg = 'coordinate order') temp = getXMLValue(header,'WholeExtent') - if (any([(IO_floatValue(temp,IO_stringPos(temp),i),i=1,5,2)] /= 0)) & + if (any([(IO_intValue(temp,IO_stringPos(temp),i),i=1,5,2)] /= 0)) & call IO_error(error_ID = 844, ext_msg = 'coordinate start') - c = [(IO_floatValue(temp,IO_stringPos(temp),i),i=2,6,2)] + c = [(IO_intValue(temp,IO_stringPos(temp),i),i=2,6,2)] temp = getXMLValue(header,'Spacing') delta = [(IO_floatValue(temp,IO_stringPos(temp),i),i=1,3)] diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 967ddb73a..26e0a909e 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -261,6 +261,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) PetscErrorCode :: ierr integer :: i, j, k, ce + phi_current = x_scal !-------------------------------------------------------------------------------------------------- ! evaluate polarization field diff --git a/src/homogenization_mechanical_RGC.f90 b/src/homogenization_mechanical_RGC.f90 index a4b6adf13..ac7dc63ea 100644 --- a/src/homogenization_mechanical_RGC.f90 +++ b/src/homogenization_mechanical_RGC.f90 @@ -720,10 +720,10 @@ module subroutine RGC_results(ho,group) case('Delta_V') call results_writeDataset(dst%volumeDiscrepancy,group,trim(prm%output(o)), & 'volume discrepancy','m³') - case('max_a_dot') + case('max_dot_a') call results_writeDataset(dst%relaxationrate_max,group,trim(prm%output(o)), & 'maximum relaxation rate','m/s') - case('avg_a_dot') + case('avg_dot_a') call results_writeDataset(dst%relaxationrate_avg,group,trim(prm%output(o)), & 'average relaxation rate','m/s') end select diff --git a/src/phase.f90 b/src/phase.f90 index 47c7a06de..da1ff9d14 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -61,9 +61,6 @@ module phase grain end type tDebugOptions - logical, dimension(:), allocatable, public :: & ! ToDo: should be protected (bug in Intel Compiler) - phase_localPlasticity !< flags phases with local constitutive law - type(tPlasticState), allocatable, dimension(:), public :: & plasticState type(tState), allocatable, dimension(:), public :: & @@ -283,16 +280,6 @@ module phase dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) end subroutine damage_anisobrittle_LiAndItsTangent - module subroutine damage_isoductile_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) - integer, intent(in) :: ph, me - real(pReal), intent(in), dimension(3,3) :: & - S - real(pReal), intent(out), dimension(3,3) :: & - Ld !< damage velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - end subroutine damage_isoductile_LiAndItsTangent - end interface @@ -413,7 +400,7 @@ end subroutine phase_init subroutine phase_allocateState(state, & NEntries,sizeState,sizeDotState,sizeDeltaState) - class(tState), intent(out) :: & + class(tState), intent(inout) :: & state integer, intent(in) :: & NEntries, & diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index ca1855a00..e9b62b702 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -11,7 +11,6 @@ submodule(phase) damage enum, bind(c); enumerator :: & DAMAGE_UNDEFINED_ID, & DAMAGE_ISOBRITTLE_ID, & - DAMAGE_ISODUCTILE_ID, & DAMAGE_ANISOBRITTLE_ID end enum @@ -39,10 +38,6 @@ submodule(phase) damage logical, dimension(:), allocatable :: mySources end function isobrittle_init - module function isoductile_init() result(mySources) - logical, dimension(:), allocatable :: mySources - end function isoductile_init - module subroutine isobrittle_deltaState(C, Fe, ph, me) integer, intent(in) :: ph,me @@ -59,10 +54,6 @@ submodule(phase) damage S end subroutine anisobrittle_dotState - module subroutine isoductile_dotState(ph,me) - integer, intent(in) :: ph,me - end subroutine isoductile_dotState - module subroutine anisobrittle_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group @@ -73,11 +64,6 @@ submodule(phase) damage character(len=*), intent(in) :: group end subroutine isobrittle_results - module subroutine isoductile_results(phase,group) - integer, intent(in) :: phase - character(len=*), intent(in) :: group - end subroutine isoductile_results - end interface contains @@ -131,7 +117,6 @@ module subroutine damage_init if (damage_active) then where(isobrittle_init() ) phase_damage = DAMAGE_ISOBRITTLE_ID - where(isoductile_init() ) phase_damage = DAMAGE_ISODUCTILE_ID where(anisobrittle_init()) phase_damage = DAMAGE_ANISOBRITTLE_ID endif @@ -178,7 +163,7 @@ module function phase_f_phi(phi,co,ce) result(f) en = material_phaseEntry(co,ce) select case(phase_damage(ph)) - case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ISODUCTILE_ID,DAMAGE_ANISOBRITTLE_ID) + case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ANISOBRITTLE_ID) f = 1.0_pReal & - phi*damageState(ph)%state(1,en) case default @@ -242,7 +227,7 @@ module function integrateDamageState(dt,co,ce) result(broken) zeta = damper(damageState(ph)%dotState(:,me),source_dotState(1:size_so,1),source_dotState(1:size_so,2)) damageState(ph)%dotState(:,me) = damageState(ph)%dotState(:,me) * zeta & - + source_dotState(1:size_so,1)* (1.0_pReal - zeta) + + source_dotState(1:size_so,1)* (1.0_pReal - zeta) r(1:size_so) = damageState(ph)%state (1:size_so,me) & - damageState(ph)%subState0(1:size_so,me) & - damageState(ph)%dotState (1:size_so,me) * dt @@ -304,9 +289,6 @@ module subroutine damage_results(group,ph) case (DAMAGE_ISOBRITTLE_ID) sourceType call isobrittle_results(ph,group//'damage/') - case (DAMAGE_ISODUCTILE_ID) sourceType - call isoductile_results(ph,group//'damage/') - case (DAMAGE_ANISOBRITTLE_ID) sourceType call anisobrittle_results(ph,group//'damage/') @@ -332,9 +314,6 @@ function phase_damage_collectDotState(ph,me) result(broken) sourceType: select case (phase_damage(ph)) - case (DAMAGE_ISODUCTILE_ID) sourceType - call isoductile_dotState(ph,me) - case (DAMAGE_ANISOBRITTLE_ID) sourceType call anisobrittle_dotState(mechanical_S(ph,me), ph,me) ! correct stress? diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index 721c27ce2..9984e5514 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -94,8 +94,8 @@ module function anisobrittle_init() result(mySources) Nmembers = count(material_phaseID==p) call phase_allocateState(damageState(p),Nmembers,1,1,0) - damageState(p)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) - if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' + damageState(p)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal) + if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi' end associate @@ -110,7 +110,7 @@ end function anisobrittle_init !-------------------------------------------------------------------------------------------------- -!> @brief calculates derived quantities from state +!> @brief !-------------------------------------------------------------------------------------------------- module subroutine anisobrittle_dotState(S, ph,me) diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index 296f7d80e..aed6a91e9 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -66,8 +66,8 @@ module function isobrittle_init() result(mySources) Nmembers = count(material_phaseID==ph) call phase_allocateState(damageState(ph),Nmembers,1,1,1) - damageState(ph)%atol = src%get_asFloat('isobrittle_atol',defaultVal=1.0e-3_pReal) - if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' + damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal) + if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi' end associate diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 deleted file mode 100644 index ea6d8acdc..000000000 --- a/src/phase_damage_isoductile.f90 +++ /dev/null @@ -1,127 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine incorporating isotropic ductile damage source mechanism -!> @details to be done -!-------------------------------------------------------------------------------------------------- -submodule(phase:damage) isoductile - - type:: tParameters !< container type for internal constitutive parameters - real(pReal) :: & - gamma_crit, & !< critical plastic strain - q - character(len=pStringLen), allocatable, dimension(:) :: & - output - end type tParameters - - type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances) - - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -module function isoductile_init() result(mySources) - - logical, dimension(:), allocatable :: mySources - - class(tNode), pointer :: & - phases, & - phase, & - sources, & - src - integer :: Ninstances,Nmembers,ph - character(len=pStringLen) :: extmsg = '' - - - mySources = source_active('isoductile') - if(count(mySources) == 0) return - - print'(/,a)', ' <<<+- phase:damage:isoductile init -+>>>' - print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT) - - - phases => config_material%get('phase') - allocate(param(phases%length)) - - do ph = 1, phases%length - if(mySources(ph)) then - phase => phases%get(ph) - sources => phase%get('damage') - - associate(prm => param(ph)) - src => sources%get(1) - - prm%q = src%get_asFloat('q') - prm%gamma_crit = src%get_asFloat('gamma_crit') - -#if defined (__GFORTRAN__) - prm%output = output_as1dString(src) -#else - prm%output = src%get_as1dString('output',defaultVal=emptyStringArray) -#endif - - ! sanity checks - if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' - if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' - - Nmembers=count(material_phaseID==ph) - call phase_allocateState(damageState(ph),Nmembers,1,1,0) - damageState(ph)%atol = src%get_asFloat('isoductile_atol',defaultVal=1.0e-3_pReal) - if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' - - end associate - -!-------------------------------------------------------------------------------------------------- -! exit if any parameter is out of range - if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isoductile)') - endif - enddo - - -end function isoductile_init - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates derived quantities from state -!-------------------------------------------------------------------------------------------------- -module subroutine isoductile_dotState(ph, me) - - integer, intent(in) :: & - ph, & - me - - - associate(prm => param(ph)) - damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me)) & - / (prm%gamma_crit*damage_phi(ph,me)**prm%q) - end associate - -end subroutine isoductile_dotState - - -!-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file -!-------------------------------------------------------------------------------------------------- -module subroutine isoductile_results(phase,group) - - integer, intent(in) :: phase - character(len=*), intent(in) :: group - - integer :: o - - associate(prm => param(phase), stt => damageState(phase)%state) - outputsLoop: do o = 1,size(prm%output) - select case(trim(prm%output(o))) - case ('f_phi') - call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³') - end select - enddo outputsLoop - end associate - -end subroutine isoductile_results - -end submodule isoductile diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 0f4e039c7..104127b53 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -15,7 +15,6 @@ submodule(phase) mechanical PLASTICITY_NONLOCAL_ID, & KINEMATICS_UNDEFINED_ID, & KINEMATICS_CLEAVAGE_OPENING_ID, & - KINEMATICS_SLIPPLANE_OPENING_ID, & KINEMATICS_THERMAL_EXPANSION_ID end enum @@ -274,7 +273,6 @@ module subroutine mechanical_init(materials,phases) ! initialize plasticity allocate(plasticState(phases%length)) allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID) - allocate(phase_localPlasticity(phases%length), source=.true.) call plastic_init() diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index b34aee58b..019838689 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -13,10 +13,6 @@ submodule(phase:mechanical) eigen logical, dimension(:), allocatable :: myKinematics end function damage_anisobrittle_init - module function damage_isoductile_init() result(myKinematics) - logical, dimension(:), allocatable :: myKinematics - end function damage_isoductile_init - module function thermalexpansion_init(kinematics_length) result(myKinematics) integer, intent(in) :: kinematics_length logical, dimension(:,:), allocatable :: myKinematics @@ -70,7 +66,6 @@ module subroutine eigendeformation_init(phases) allocate(model_damage(phases%length), source = KINEMATICS_UNDEFINED_ID) where(damage_anisobrittle_init()) model_damage = KINEMATICS_cleavage_opening_ID - where(damage_isoductile_init()) model_damage = KINEMATICS_slipplane_opening_ID end subroutine eigendeformation_init @@ -201,11 +196,6 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS active = .true. - case (KINEMATICS_slipplane_opening_ID) - call damage_isoductile_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, en) - Li = Li + my_Li - dLi_dS = dLi_dS + my_dLi_dS - active = .true. end select if(.not. active) return diff --git a/src/phase_mechanical_eigen_slipplaneopening.f90 b/src/phase_mechanical_eigen_slipplaneopening.f90 deleted file mode 100644 index 2fd14700d..000000000 --- a/src/phase_mechanical_eigen_slipplaneopening.f90 +++ /dev/null @@ -1,184 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine incorporating kinematics resulting from opening of slip planes -!> @details to be done -!-------------------------------------------------------------------------------------------------- -submodule(phase:eigen) slipplaneopening - - integer, dimension(:), allocatable :: damage_isoductile_instance - - type :: tParameters !< container type for internal constitutive parameters - integer :: & - sum_N_sl !< total number of cleavage planes - real(pReal) :: & - dot_o, & !< opening rate of cleavage planes - q !< damage rate sensitivity - real(pReal), dimension(:), allocatable :: & - g_crit - real(pReal), dimension(:,:,:), allocatable :: & - P_d, & - P_t, & - P_n - end type tParameters - - type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances) - - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -module function damage_isoductile_init() result(myKinematics) - - logical, dimension(:), allocatable :: myKinematics - - integer :: p,i - character(len=pStringLen) :: extmsg = '' - integer, dimension(:), allocatable :: N_sl - real(pReal), dimension(:,:), allocatable :: d,n,t - class(tNode), pointer :: & - phases, & - phase, & - mech, & - pl, & - kinematics, & - kinematic_type - - - myKinematics = kinematics_active2('isoductile') - if(count(myKinematics) == 0) return - print'(/,a)', ' <<<+- phase:mechanical:eigen:slipplaneopening init -+>>>' - print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT) - - - phases => config_material%get('phase') - allocate(param(phases%length)) - - do p = 1, phases%length - if(myKinematics(p)) then - phase => phases%get(p) - mech => phase%get('mechanical') - pl => mech%get('plastic') - - kinematics => phase%get('damage') - - associate(prm => param(p)) - kinematic_type => kinematics%get(1) - - prm%dot_o = kinematic_type%get_asFloat('dot_o') - prm%q = kinematic_type%get_asFloat('q') - N_sl = pl%get_as1dInt('N_sl') - prm%sum_N_sl = sum(abs(N_sl)) - - d = lattice_slip_direction (N_sl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) - t = lattice_slip_transverse(N_sl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) - n = lattice_slip_normal (N_sl,phase%get_asString('lattice'),& - phase%get_asFloat('c/a',defaultVal=0.0_pReal)) - allocate(prm%P_d(3,3,size(d,2)),prm%P_t(3,3,size(t,2)),prm%P_n(3,3,size(n,2))) - - do i=1, size(n,2) - prm%P_d(1:3,1:3,i) = math_outer(d(1:3,i), n(1:3,i)) - prm%P_t(1:3,1:3,i) = math_outer(t(1:3,i), n(1:3,i)) - prm%P_n(1:3,1:3,i) = math_outer(n(1:3,i), n(1:3,i)) - enddo - - prm%g_crit = kinematic_type%get_as1dFloat('g_crit',requiredSize=size(N_sl)) - - ! expand: family => system - prm%g_crit = math_expand(prm%g_crit,N_sl) - - ! sanity checks - if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_n' - if (prm%dot_o <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_sdot0' - if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' anisoDuctile_critLoad' - -!-------------------------------------------------------------------------------------------------- -! exit if any parameter is out of range - if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(slipplane_opening)') - - end associate - endif - enddo - - -end function damage_isoductile_init - - -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient -!-------------------------------------------------------------------------------------------------- -module subroutine damage_isoductile_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) - - integer, intent(in) :: & - ph, me - real(pReal), intent(in), dimension(3,3) :: & - S - real(pReal), intent(out), dimension(3,3) :: & - Ld !< damage velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - - integer :: & - i, k, l, m, n - real(pReal) :: & - traction_d, traction_t, traction_n, traction_crit, & - udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt - - - associate(prm => param(ph)) - Ld = 0.0_pReal - dLd_dTstar = 0.0_pReal - do i = 1, prm%sum_N_sl - - traction_d = math_tensordot(S,prm%P_d(1:3,1:3,i)) - traction_t = math_tensordot(S,prm%P_t(1:3,1:3,i)) - traction_n = math_tensordot(S,prm%P_n(1:3,1:3,i)) - - traction_crit = prm%g_crit(i)* damage_phi(ph,me) - - udotd = sign(1.0_pReal,traction_d)* prm%dot_o* ( abs(traction_d)/traction_crit & - - abs(traction_d)/prm%g_crit(i))**prm%q - udott = sign(1.0_pReal,traction_t)* prm%dot_o* ( abs(traction_t)/traction_crit & - - abs(traction_t)/prm%g_crit(i))**prm%q - udotn = prm%dot_o* ( max(0.0_pReal,traction_n)/traction_crit & - - max(0.0_pReal,traction_n)/prm%g_crit(i))**prm%q - - if (dNeq0(traction_d)) then - dudotd_dt = udotd*prm%q/traction_d - else - dudotd_dt = 0.0_pReal - endif - if (dNeq0(traction_t)) then - dudott_dt = udott*prm%q/traction_t - else - dudott_dt = 0.0_pReal - endif - if (dNeq0(traction_n)) then - dudotn_dt = udotn*prm%q/traction_n - else - dudotn_dt = 0.0_pReal - endif - - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & - + dudotd_dt*prm%P_d(k,l,i)*prm%P_d(m,n,i) & - + dudott_dt*prm%P_t(k,l,i)*prm%P_t(m,n,i) & - + dudotn_dt*prm%P_n(k,l,i)*prm%P_n(m,n,i) - - Ld = Ld & - + udotd*prm%P_d(1:3,1:3,i) & - + udott*prm%P_t(1:3,1:3,i) & - + udotn*prm%P_n(1:3,1:3,i) - enddo - - end associate - -end subroutine damage_isoductile_LiAndItsTangent - -end submodule slipplaneopening diff --git a/src/phase_mechanical_eigen_thermalexpansion.f90 b/src/phase_mechanical_eigen_thermalexpansion.f90 index 3cfeb2f06..dba02d70e 100644 --- a/src/phase_mechanical_eigen_thermalexpansion.f90 +++ b/src/phase_mechanical_eigen_thermalexpansion.f90 @@ -29,7 +29,6 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics) logical, dimension(:,:), allocatable :: myKinematics integer :: Ninstances,p,i,k - real(pReal), dimension(:), allocatable :: temp class(tNode), pointer :: & phases, & phase, & @@ -57,25 +56,26 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics) do k = 1, kinematics%length if(myKinematics(k,p)) then associate(prm => param(kinematics_thermal_expansion_instance(p))) - kinematic_type => kinematics%get(k) + 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=0.0_pReal) - ! read up to three parameters (constant, linear, quadratic with T) - temp = kinematic_type%get_as1dFloat('A_11') - prm%A(1,1,1:size(temp)) = temp - temp = kinematic_type%get_as1dFloat('A_33',defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp)) - prm%A(3,3,1:size(temp)) = temp - do i=1, size(prm%A,3) - prm%A(1:3,1:3,i) = lattice_applyLatticeSymmetry33(prm%A(1:3,1:3,i),& - phase%get_asString('lattice')) - enddo + 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,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,3) = kinematic_type%get_asFloat('A_33,T^2',defaultVal=0.0_pReal) + endif + do i=1, size(prm%A,3) + prm%A(1:3,1:3,i) = lattice_applyLatticeSymmetry33(prm%A(1:3,1:3,i),phase_lattice(p)) + enddo end associate endif enddo enddo - end function thermalexpansion_init diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 7c9accc8d..0fa742b01 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -247,9 +247,8 @@ module function plastic_dislotungsten_init() result(myPlasticity) endIndex = endIndex + prm%sum_N_sl stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:) dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:) - plasticState(ph)%atol(startIndex:endIndex) = 1.0e-2_pReal - ! global alias - plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) + if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal) allocate(dst%threshold_stress(prm%sum_N_sl,Nmembers), source=0.0_pReal) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index f1f39c71f..c70d33f78 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -436,15 +436,14 @@ module function plastic_dislotwin_init() result(myPlasticity) endIndex = endIndex + prm%sum_N_sl stt%gamma_sl=>plasticState(ph)%state(startIndex:endIndex,:) dot%gamma_sl=>plasticState(ph)%dotState(startIndex:endIndex,:) - plasticState(ph)%atol(startIndex:endIndex) = 1.0e-2_pReal - ! global alias - plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) + if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw stt%f_tw=>plasticState(ph)%state(startIndex:endIndex,:) dot%f_tw=>plasticState(ph)%dotState(startIndex:endIndex,:) - plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tw',defaultVal=1.0e-7_pReal) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tw',defaultVal=1.0e-6_pReal) if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tw' startIndex = endIndex + 1 diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index aee9bc95d..6694d6548 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -139,8 +139,6 @@ module function plastic_isotropic_init() result(myPlasticity) dot%gamma => plasticState(ph)%dotState(2,:) plasticState(ph)%atol(2) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) if (plasticState(ph)%atol(2) < 0.0_pReal) extmsg = trim(extmsg)//' atol_gamma' - ! global alias - plasticState(ph)%slipRate => plasticState(ph)%dotState(2:2,:) end associate diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 461196d3f..69be33958 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -194,8 +194,6 @@ module function plastic_kinehardening_init() result(myPlasticity) dot%accshear => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' - ! global alias - plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:) o = plasticState(ph)%offsetDeltaState startIndex = endIndex + 1 diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index efd45e97b..fbb092d0e 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -81,11 +81,11 @@ submodule(phase:plastic) nonlocal c_sol, & !< concentration of solid solution in atomic parts p, & !< parameter for kinetic law (Kocks,Argon,Ashby) q, & !< parameter for kinetic law (Kocks,Argon,Ashby) - eta, & !< viscosity for dislocation glide in Pa s + B, & !< drag coefficient in Pa s nu_a, & !< attack frequency in Hz chi_surface, & !< transmissivity at free surface chi_GB, & !< transmissivity at grain boundary (identified by different texture) - f_c, & !< safety factor for CFL flux condition + C_CFL, & !< safety factor for CFL flux condition f_ed_mult, & !< factor that determines how much edge dislocations contribute to multiplication (0...1) f_F, & f_ed, & @@ -112,7 +112,7 @@ submodule(phase:plastic) nonlocal nonSchmid_pos, & nonSchmid_neg !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws) integer :: & - sum_N_sl + sum_N_sl = 0 integer, dimension(:), allocatable :: & colinearSystem !< colinear system to the active slip system (only valid for fcc!) character(len=pStringLen), dimension(:), allocatable :: & @@ -211,6 +211,7 @@ module function plastic_nonlocal_init() result(myPlasticity) phases => config_material%get('phase') + allocate(geom(phases%length)) allocate(param(phases%length)) @@ -230,15 +231,14 @@ module function plastic_nonlocal_init() result(myPlasticity) mech => phase%get('mechanical') pl => mech%get('plastic') - phase_localPlasticity(ph) = .not. pl%contains('nonlocal') - + plasticState(ph)%nonlocal = pl%get_asBool('flux',defaultVal=.True.) #if defined (__GFORTRAN__) prm%output = output_as1dString(pl) #else prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray) #endif - prm%atol_rho = pl%get_asFloat('atol_rho',defaultVal=1.0e4_pReal) + prm%atol_rho = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) ! This data is read in already in lattice prm%mu = lattice_mu(ph) @@ -317,13 +317,13 @@ module function plastic_nonlocal_init() result(myPlasticity) prm%rho_significant = pl%get_asFloat('rho_significant') prm%rho_min = pl%get_asFloat('rho_min', 0.0_pReal) - prm%f_c = pl%get_asFloat('f_c',defaultVal=2.0_pReal) + prm%C_CFL = pl%get_asFloat('C_CFL',defaultVal=2.0_pReal) prm%V_at = pl%get_asFloat('V_at') prm%D_0 = pl%get_asFloat('D_0') prm%Q_cl = pl%get_asFloat('Q_cl') prm%f_F = pl%get_asFloat('f_F') - prm%f_ed = pl%get_asFloat('f_ed') !,'edgejogs' + prm%f_ed = pl%get_asFloat('f_ed') prm%w = pl%get_asFloat('w') prm%Q_sol = pl%get_asFloat('Q_sol') prm%f_sol = pl%get_asFloat('f_sol') @@ -331,7 +331,7 @@ module function plastic_nonlocal_init() result(myPlasticity) prm%p = pl%get_asFloat('p_sl') prm%q = pl%get_asFloat('q_sl') - prm%eta = pl%get_asFloat('eta') + prm%B = pl%get_asFloat('B') prm%nu_a = pl%get_asFloat('nu_a') ! ToDo: discuss logic @@ -363,8 +363,8 @@ module function plastic_nonlocal_init() result(myPlasticity) if (any(prm%peierlsstress < 0.0_pReal)) extmsg = trim(extmsg)//' tau_peierls' if (any(prm%minDipoleHeight < 0.0_pReal)) extmsg = trim(extmsg)//' d_ed or d_sc' - if (prm%eta <= 0.0_pReal) extmsg = trim(extmsg)//' eta' - if (prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' + if (prm%B < 0.0_pReal) extmsg = trim(extmsg)//' B' + if (prm%Q_cl < 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' if (prm%nu_a <= 0.0_pReal) extmsg = trim(extmsg)//' nu_a' if (prm%w <= 0.0_pReal) extmsg = trim(extmsg)//' w' if (prm%D_0 < 0.0_pReal) extmsg = trim(extmsg)//' D_0' @@ -373,7 +373,7 @@ module function plastic_nonlocal_init() result(myPlasticity) if (prm%rho_min < 0.0_pReal) extmsg = trim(extmsg)//' rho_min' if (prm%rho_significant < 0.0_pReal) extmsg = trim(extmsg)//' rho_significant' if (prm%atol_rho < 0.0_pReal) extmsg = trim(extmsg)//' atol_rho' - if (prm%f_c < 0.0_pReal) extmsg = trim(extmsg)//' f_c' + if (prm%C_CFL < 0.0_pReal) extmsg = trim(extmsg)//' C_CFL' if (prm%p <= 0.0_pReal .or. prm%p > 1.0_pReal) extmsg = trim(extmsg)//' p_sl' if (prm%q < 1.0_pReal .or. prm%q > 2.0_pReal) extmsg = trim(extmsg)//' q_sl' @@ -417,7 +417,6 @@ module function plastic_nonlocal_init() result(myPlasticity) allocate(geom(ph)%V_0(Nmembers)) call storeGeometry(ph) - plasticState(ph)%nonlocal = pl%get_asBool('nonlocal') if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) & call IO_error(212,ext_msg='IPneighborhood does not exist') @@ -489,10 +488,9 @@ module function plastic_nonlocal_init() result(myPlasticity) stt%gamma => plasticState(ph)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers) dot%gamma => plasticState(ph)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers) del%gamma => plasticState(ph)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers) - plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-2_pReal) + plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-6_pReal) if(any(plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) & extmsg = trim(extmsg)//' atol_gamma' - plasticState(ph)%slipRate => plasticState(ph)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers) stt%rho_forest => plasticState(ph)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nmembers) stt%v => plasticState(ph)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers) @@ -643,7 +641,7 @@ module subroutine nonlocal_dependentState(ph, en, ip, el) !################################################################################################# rho0 = getRho0(ph,en) - if (.not. phase_localPlasticity(ph) .and. prm%shortRangeStressCorrection) then + if (plasticState(ph)%nonlocal .and. prm%shortRangeStressCorrection) then invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en)) invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,en)) @@ -739,17 +737,6 @@ module subroutine nonlocal_dependentState(ph, en, ip, el) enddo endif -#ifdef DEBUG - if (debugConstitutive%extensive & - .and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)& - .or. .not. debugConstitutive%selective)) then - print'(/,a,i8,1x,i2,1x,i1,/)', '<< CONST >> nonlocal_microstructure at el ip ',el,ip - print'(a,/,12x,12(e10.3,1x))', '<< CONST >> rhoForest', stt%rho_forest(:,en) - print'(a,/,12x,12(f10.5,1x))', '<< CONST >> tauThreshold / MPa', dst%tau_pass(:,en)*1e-6 - print'(a,/,12x,12(f10.5,1x),/)', '<< CONST >> tauBack / MPa', dst%tau_back(:,en)*1e-6 - endif -#endif - end associate end subroutine nonlocal_dependentState @@ -994,13 +981,13 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & gdot !< shear rates real(pReal), dimension(param(ph)%sum_N_sl) :: & tau, & !< current resolved shear stress - vClimb !< climb velocity of edge dipoles + v_climb !< climb velocity of edge dipoles real(pReal), dimension(param(ph)%sum_N_sl,2) :: & rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) dLower, & !< minimum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws real(pReal) :: & - selfDiffusion !< self diffusion + D_SD if (timestep <= 0.0_pReal) then plasticState(ph)%dotState = 0.0_pReal @@ -1025,17 +1012,9 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) -#ifdef DEBUG - if (debugConstitutive%basic & - .and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip) & - .or. .not. debugConstitutive%selective)) then - print'(a,/,10(12x,12(e12.5,1x),/))', '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip - print'(a,/,4(12x,12(e12.5,1x),/))', '<< CONST >> gdot / 1/s',gdot - endif -#endif - !**************************************************************************** - !*** limits for stable dipole height + + ! limits for stable dipole height do s = 1,ns tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,en) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal @@ -1052,8 +1031,8 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & dUpper = max(dUpper,dLower) - !**************************************************************************** - !*** dislocation multiplication + + ! dislocation multiplication rhoDotMultiplication = 0.0_pReal isBCC: if (phase_lattice(ph) == 'cI') then forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal) @@ -1068,16 +1047,16 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & else isBCC rhoDotMultiplication(:,1:4) = spread( & (sum(abs(gdot(:,1:2)),2) * prm%f_ed_mult + sum(abs(gdot(:,3:4)),2)) & - * sqrt(stt%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4) + * sqrt(stt%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4) ! eq. 3.26 endif isBCC forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) !**************************************************************************** - !*** calculate dipole formation and annihilation + ! dipole formation and annihilation - !*** formation by glide + ! formation by glide do c = 1,2 rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%b_sl & * ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile @@ -1102,7 +1081,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & enddo - !*** athermal annihilation + ! athermal annihilation rhoDotAthermalAnnihilation = 0.0_pReal forall (c=1:2) & rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%b_sl & @@ -1117,13 +1096,13 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & * 0.25_pReal * sqrt(stt%rho_forest(s,en)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed - !*** thermally activated annihilation of edge dipoles by climb + ! thermally activated annihilation of edge dipoles by climb rhoDotThermalAnnihilation = 0.0_pReal - selfDiffusion = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature)) - vClimb = prm%V_at * selfDiffusion * prm%mu & - / ( kB * Temperature * PI * (1.0_pReal-prm%nu) * (dUpper(:,1) + dLower(:,1))) + D_SD = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature)) ! eq. 3.53 + v_climb = D_SD * prm%mu * prm%V_at & + / (PI * (1.0_pReal-prm%nu) * (dUpper(:,1) + dLower(:,1)) * kB * Temperature) ! eq. 3.54 forall (s = 1:ns, dUpper(s,1) > dLower(s,1)) & - rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * vClimb(s) / (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) & - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have @@ -1231,24 +1210,23 @@ function rhoDotFlux(timestep,ph,en,ip,el) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) - forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) !**************************************************************************** !*** calculate dislocation fluxes (only for nonlocal plasticity) rhoDotFlux = 0.0_pReal - if (.not. phase_localPlasticity(ph)) then + if (plasticState(ph)%nonlocal) then !*** check CFL (Courant-Friedrichs-Lewy) condition for flux if (any( abs(gdot) > 0.0_pReal & ! any active slip system ... - .and. prm%f_c * abs(v0) * timestep & + .and. prm%C_CFL * abs(v0) * timestep & > IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) #ifdef DEBUG if (debugConstitutive%extensive) then print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', & maxval(abs(v0), abs(gdot) > 0.0_pReal & - .and. prm%f_c * abs(v0) * timestep & + .and. prm%C_CFL * abs(v0) * timestep & > IPvolume(ip,el) / maxval(IParea(:,ip,el))), & ' at a timestep of ',timestep print*, '<< CONST >> enforcing cutback !!!' @@ -1436,22 +1414,16 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) if (neighbor_e <= 0 .or. neighbor_i <= 0) then !* FREE SURFACE - !* Set surface transmissivity to the value specified in the material.config forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface) elseif (neighbor_phase /= ph) then !* PHASE BOUNDARY - !* If we encounter a different nonlocal phase at the neighbor, - !* we consider this to be a real "physical" phase boundary, so completely incompatible. - !* If one of the two phases has a local plasticity law, - !* we do not consider this to be a phase boundary, so completely compatible. - if (.not. phase_localPlasticity(neighbor_phase) .and. .not. phase_localPlasticity(ph)) & + if (plasticState(neighbor_phase)%nonlocal .and. plasticState(ph)%nonlocal) & forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pReal elseif (prm%chi_GB >= 0.0_pReal) then - !* GRAIN BOUNDARY ! - !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config) + !* GRAIN BOUNDARY if (any(dNeq(phase_orientation0(ph)%data(en)%asQuaternion(), & phase_orientation0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. & - (.not. phase_localPlasticity(neighbor_phase))) & + plasticState(neighbor_phase)%nonlocal) & forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB) else !* GRAIN BOUNDARY ? @@ -1645,13 +1617,13 @@ end subroutine stateInit !-------------------------------------------------------------------------------------------------- !> @brief calculates kinetics !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Temperature, ph) +pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T, ph) integer, intent(in) :: & c, & !< dislocation character (1:edge, 2:screw) ph real(pReal), intent(in) :: & - Temperature !< temperature + T !< T real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: & tau, & !< resolved external shear stress (without non Schmid effects) tauNS, & !< resolved external shear stress (including non Schmid effects) @@ -1669,22 +1641,15 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem tauEff, & !< effective shear stress tPeierls, & !< waiting time in front of a peierls barriers tSolidSolution, & !< waiting time in front of a solid solution obstacle - vViscous, & !< viscous glide velocity dtPeierls_dtau, & !< derivative with respect to resolved shear stress dtSolidSolution_dtau, & !< derivative with respect to resolved shear stress - meanfreepath_S, & !< mean free travel distance for dislocations between two solid solution obstacles - meanfreepath_P, & !< mean free travel distance for dislocations between two Peierls barriers - jumpWidth_P, & !< depth of activated area - jumpWidth_S, & !< depth of activated area - activationLength_P, & !< length of activated dislocation line - activationLength_S, & !< length of activated dislocation line + lambda_S, & !< mean free distance between two solid solution obstacles + lambda_P, & !< mean free distance between two Peierls barriers activationVolume_P, & !< volume that needs to be activated to overcome barrier activationVolume_S, & !< volume that needs to be activated to overcome barrier activationEnergy_P, & !< energy that is needed to overcome barrier - activationEnergy_S, & !< energy that is needed to overcome barrier criticalStress_P, & !< maximum obstacle strength - criticalStress_S, & !< maximum obstacle strength - mobility !< dislocation mobility + criticalStress_S !< maximum obstacle strength associate(prm => param(ph)) v = 0.0_pReal @@ -1695,57 +1660,42 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem if (abs(tau(s)) > tauThreshold(s)) then !* Peierls contribution - !* Effective stress includes non Schmid constributions - !* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity - tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) ! ensure that the effective stress is positive - meanfreepath_P = prm%b_sl(s) - jumpWidth_P = prm%b_sl(s) - activationLength_P = prm%w *prm%b_sl(s) - activationVolume_P = activationLength_P * jumpWidth_P * prm%b_sl(s) + tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) + lambda_P = prm%b_sl(s) + activationVolume_P = prm%w *prm%b_sl(s)**3.0_pReal criticalStress_P = prm%peierlsStress(s,c) activationEnergy_P = criticalStress_P * activationVolume_P - tauRel_P = min(1.0_pReal, tauEff / criticalStress_P) ! ensure that the activation probability cannot become greater than one + tauRel_P = min(1.0_pReal, tauEff / criticalStress_P) tPeierls = 1.0_pReal / prm%nu_a & - * exp(activationEnergy_P / (kB * Temperature) & + * exp(activationEnergy_P / (kB * T) & * (1.0_pReal - tauRel_P**prm%p)**prm%q) - if (tauEff < criticalStress_P) then - dtPeierls_dtau = tPeierls * prm%p * prm%q * activationVolume_P / (kB * Temperature) & - * (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) * tauRel_P**(prm%p-1.0_pReal) - else - dtPeierls_dtau = 0.0_pReal - endif + dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (kB * 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) - !* Contribution from solid solution strengthening - !* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity + ! Contribution from solid solution strengthening tauEff = abs(tau(s)) - tauThreshold(s) - meanfreepath_S = prm%b_sl(s) / sqrt(prm%c_sol) - jumpWidth_S = prm%f_sol * prm%b_sl(s) - activationLength_S = prm%b_sl(s) / sqrt(prm%c_sol) - activationVolume_S = activationLength_S * jumpWidth_S * prm%b_sl(s) - activationEnergy_S = prm%Q_sol - criticalStress_S = activationEnergy_S / activationVolume_S - tauRel_S = min(1.0_pReal, tauEff / criticalStress_S) ! ensure that the activation probability cannot become greater than one + lambda_S = prm%b_sl(s) / sqrt(prm%c_sol) + activationVolume_S = prm%f_sol * prm%b_sl(s)**3.0_pReal / sqrt(prm%c_sol) + criticalStress_S = prm%Q_sol / activationVolume_S + tauRel_S = min(1.0_pReal, tauEff / criticalStress_S) tSolidSolution = 1.0_pReal / prm%nu_a & - * exp(activationEnergy_S / (kB * Temperature)* (1.0_pReal - tauRel_S**prm%p)**prm%q) - if (tauEff < criticalStress_S) then - dtSolidSolution_dtau = tSolidSolution * prm%p * prm%q * activationVolume_S / (kB * Temperature) & - * (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal)* tauRel_S**(prm%p-1.0_pReal) - else - dtSolidSolution_dtau = 0.0_pReal - endif + * 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) & + * (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) !* viscous glide velocity tauEff = abs(tau(s)) - tauThreshold(s) - mobility = prm%b_sl(s) / prm%eta - vViscous = mobility * tauEff - !* Mean velocity results from waiting time at peierls barriers and solid solution obstacles with respective meanfreepath of - !* free flight at glide velocity in between. - !* adopt sign from resolved stress + v(s) = sign(1.0_pReal,tau(s)) & - / (tPeierls / meanfreepath_P + tSolidSolution / meanfreepath_S + 1.0_pReal / vViscous) - dv_dtau(s) = v(s)**2.0_pReal * (dtSolidSolution_dtau / meanfreepath_S + mobility /vViscous**2.0_pReal) - dv_dtauNS(s) = v(s)**2.0_pReal * dtPeierls_dtau / meanfreepath_P + / (tPeierls / lambda_P + tSolidSolution / lambda_S + prm%B /(prm%b_sl(s) * tauEff)) + dv_dtau(s) = v(s)**2.0_pReal * (dtSolidSolution_dtau / lambda_S + prm%B / (prm%b_sl(s) * tauEff**2.0_pReal)) + dv_dtauNS(s) = v(s)**2.0_pReal * dtPeierls_dtau / lambda_P + endif enddo diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index af6734add..45e7a5f14 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -247,7 +247,6 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) stt%xi_twin = spread(xi_0_tw, 2, Nmembers) dot%xi_twin => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) - if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl @@ -255,15 +254,12 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) dot%gamma_slip => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' - ! global alias - plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw stt%gamma_twin => plasticState(ph)%state (startIndex:endIndex,:) dot%gamma_twin => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) - if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' end associate diff --git a/src/prec.f90 b/src/prec.f90 index 1f884880f..7613cfe46 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -51,10 +51,7 @@ module prec end type type, extends(tState) :: tPlasticState - logical :: & - nonlocal = .false. - real(pReal), pointer, dimension(:,:) :: & - slipRate !< slip rate + logical :: nonlocal = .false. end type type :: tSourceState