223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'preparing-alpha4' into 'development' 

polishing and fixing bugs

See merge request damask/DAMASK!405
This commit is contained in:
Franz Roters 2021-06-29 14:12:00 +00:00
parent b9f273099d 7745640393
commit b69bae272f
45 changed files with 261 additions and 706 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
.gitlab-ci.yml
View File

@ -206,17 +206,6 @@ core:
# - master # - master
# - release # - 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
################################################################################################### ###################################################################################################
SpectralRuntime: SpectralRuntime:
stage: performance stage: performance

2
PRIVATE

@ -1 +1 @@
Subproject commit 79929319c0756bdff58dd69cf4969774a880c1bf Subproject commit 698b255cc0a6e496f1917992dd4dafc6d50b66b3

2
examples/config/Kinematics_Thermal_Expansion.config
View File

@ -1,2 +0,0 @@
(kinematics) thermal_expansion
thermal_expansion11 0.0000231

63
examples/config/Phase_Nonlocal_Aluminum.config
View File

@ -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

62
examples/config/Phase_Nonlocal_Nickel.config
View File

@ -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

16
examples/config/Phase_Phenopowerlaw_Aluminum.yaml
View File

@ -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]

17
examples/config/Phase_Phenopowerlaw_BCC-Ferrite.yaml
View File

@ -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]

2
examples/config/homogenization/mechanical/RGC_8grains.yaml
View File

@ -6,5 +6,5 @@
a_g: [0.0, 0.0, 0.0] a_g: [0.0, 0.0, 0.0]
c_alpha: 2.0 c_alpha: 2.0
cluster_size: [2, 2, 2] 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 xi_alpha: 10.0

2
examples/config/phase/mechanical/eigen/thermalexpansion_Al.yaml
View File

@ -1,5 +1,5 @@
type: thermalexpansion type: thermalexpansion
references: references:
- en.wikipedia.org/wiki/Thermal_expansion - en.wikipedia.org/wiki/Thermal_expansion
A_11: [23.1e-6] A_11: 23.1e-6
T_ref: 293.15 T_ref: 293.15

2
examples/config/phase/mechanical/eigen/thermalexpansion_Au.yaml
View File

@ -1,5 +1,5 @@
type: thermalexpansion type: thermalexpansion
references: references:
- en.wikipedia.org/wiki/Thermal_expansion - en.wikipedia.org/wiki/Thermal_expansion
A_11: [14e-6] A_11: 14e-6
T_ref: 293.15 T_ref: 293.15

4
examples/config/phase/mechanical/eigen/thermalexpansion_C35E.yaml
View File

@ -1,5 +1,7 @@
type: thermalexpansion type: thermalexpansion
references: references:
- en.wikipedia.org/wiki/Thermal_expansion, fitted from image description - 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 T_ref: 273.0

2
examples/config/phase/mechanical/eigen/thermalexpansion_Cu.yaml
View File

@ -1,5 +1,5 @@
type: thermalexpansion type: thermalexpansion
references: references:
- en.wikipedia.org/wiki/Thermal_expansion - en.wikipedia.org/wiki/Thermal_expansion
A_11: [17e-6] A_11: 17e-6
T_ref: 293.15 T_ref: 293.15

2
examples/config/phase/mechanical/eigen/thermalexpansion_Fe.yaml
View File

@ -1,5 +1,5 @@
type: thermalexpansion type: thermalexpansion
references: references:
- en.wikipedia.org/wiki/Thermal_expansion - en.wikipedia.org/wiki/Thermal_expansion
A_11: [11.8e-6] A_11: 11.8e-6
T_ref: 293.15 T_ref: 293.15

2
examples/config/phase/mechanical/eigen/thermalexpansion_W.yaml
View File

@ -1,5 +1,5 @@
type: thermalexpansion type: thermalexpansion
references: references:
- en.wikipedia.org/wiki/Thermal_expansion - en.wikipedia.org/wiki/Thermal_expansion
A_11: [4.5e-6] A_11: 4.5e-6
T_ref: 293.15 T_ref: 293.15

3
examples/config/phase/mechanical/eigen/thermalexpansion_X20Cr13.yaml
View File

@ -1,5 +1,6 @@
type: thermalexpansion type: thermalexpansion
references: references:
- en.wikipedia.org/wiki/Thermal_expansion, fitted from image description - 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 T_ref: 273.0

3
examples/config/phase/mechanical/elastic/Hooke_Au.yaml
View File

@ -2,7 +2,8 @@ type: Hooke
references: references:
- J.P. Hirth and J. Lothe, - J.P. Hirth and J. Lothe,
Theory of Dislocations, 1982, Theory of Dislocations, 1982,
John Wiley & Sons John Wiley & Sons,
page 837
C_11: 186e9 C_11: 186e9
C_12: 157e9 C_12: 157e9
C_44: 42e9 C_44: 42e9

5
examples/config/phase/mechanical/elastic/Hooke_Fe.yaml
View File

@ -2,7 +2,8 @@ type: Hooke
references: references:
- J.P. Hirth and J. Lothe, - J.P. Hirth and J. Lothe,
Theory of Dislocations, 1982, Theory of Dislocations, 1982,
John Wiley & Sons John Wiley & Sons,
page 837
C_11: 242e9 C_11: 242e9
C_12: 146.5e9 C_12: 146.5e9
C_44: 11.2e9 C_44: 112e9

3
examples/config/phase/mechanical/elastic/Hooke_Ni.yaml
View File

@ -2,7 +2,8 @@ type: Hooke
references: references:
- J.P. Hirth and J. Lothe, - J.P. Hirth and J. Lothe,
Theory of Dislocations, 1982, Theory of Dislocations, 1982,
John Wiley & Sons John Wiley & Sons,
page 837
C_11: 246.5e9 C_11: 246.5e9
C_12: 147.3e9 C_12: 147.3e9
C_44: 124.7e9 C_44: 124.7e9

49
examples/config/phase/mechanical/plastic/nonlocal_Al.yaml Normal file
View File

@ -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

49
examples/config/phase/mechanical/plastic/nonlocal_Ni.yaml Normal file
View File

@ -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

15
examples/config/phase/mechanical/plastic/phenopowerlaw_Al.yaml Normal file
View File

@ -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

6
examples/config/phase/mechanical/plastic/phenopowerlaw_Au.yaml
View File

@ -9,9 +9,9 @@ references:
output: [xi_sl, gamma_sl] output: [xi_sl, gamma_sl]
N_sl: [12] N_sl: [12]
n_sl: 83.3 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 a_sl: 1.0
h_0_sl-sl: 75.0e6
xi_0_sl: [26.25e6] xi_0_sl: [26.25e6]
xi_inf_sl: [53.0e6] 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

15
examples/config/phase/mechanical/plastic/phenopowerlaw_Cu.yaml Normal file
View File

@ -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

14
examples/config/phase/mechanical/plastic/phenopowerlaw_DP-steel-ferrite.yaml Normal file
View File

@ -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

2
python/tests/reference/Result/4grains2x4x3.material.yaml vendored
View File

@ -633,7 +633,7 @@ homogenization:
a_g: [0.0, 0.0, 0.0] a_g: [0.0, 0.0, 0.0]
c_alpha: 2.0 c_alpha: 2.0
cluster_size: [2, 2, 2] 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 xi_alpha: 10.0
phase: phase:

2
src/commercialFEM_fileList.f90
View File

@ -32,14 +32,12 @@
#include "phase_mechanical_plastic_nonlocal.f90" #include "phase_mechanical_plastic_nonlocal.f90"
#include "phase_mechanical_eigen.f90" #include "phase_mechanical_eigen.f90"
#include "phase_mechanical_eigen_cleavageopening.f90" #include "phase_mechanical_eigen_cleavageopening.f90"
#include "phase_mechanical_eigen_slipplaneopening.f90"
#include "phase_mechanical_eigen_thermalexpansion.f90" #include "phase_mechanical_eigen_thermalexpansion.f90"
#include "phase_thermal.f90" #include "phase_thermal.f90"
#include "phase_thermal_dissipation.f90" #include "phase_thermal_dissipation.f90"
#include "phase_thermal_externalheat.f90" #include "phase_thermal_externalheat.f90"
#include "phase_damage.f90" #include "phase_damage.f90"
#include "phase_damage_isobrittle.f90" #include "phase_damage_isobrittle.f90"
#include "phase_damage_isoductile.f90"
#include "phase_damage_anisobrittle.f90" #include "phase_damage_anisobrittle.f90"
#include "homogenization.f90" #include "homogenization.f90"
#include "homogenization_mechanical.f90" #include "homogenization_mechanical.f90"

5
src/grid/discretization_grid.f90
View File

@ -258,7 +258,6 @@ subroutine readVTI(grid,geomSize,origin,material)
character(len=:), allocatable :: temp character(len=:), allocatable :: temp
real(pReal), dimension(:), allocatable :: delta real(pReal), dimension(:), allocatable :: delta
integer, dimension(:), allocatable :: stringPos
integer :: i integer :: i
@ -266,9 +265,9 @@ subroutine readVTI(grid,geomSize,origin,material)
call IO_error(error_ID = 844, ext_msg = 'coordinate order') call IO_error(error_ID = 844, ext_msg = 'coordinate order')
temp = getXMLValue(header,'WholeExtent') 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') 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') temp = getXMLValue(header,'Spacing')
delta = [(IO_floatValue(temp,IO_stringPos(temp),i),i=1,3)] delta = [(IO_floatValue(temp,IO_stringPos(temp),i),i=1,3)]

1
src/grid/grid_damage_spectral.f90
View File

@ -261,6 +261,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer :: i, j, k, ce integer :: i, j, k, ce
phi_current = x_scal phi_current = x_scal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate polarization field ! evaluate polarization field

4
src/homogenization_mechanical_RGC.f90
View File

@ -720,10 +720,10 @@ module subroutine RGC_results(ho,group)
case('Delta_V') case('Delta_V')
call results_writeDataset(dst%volumeDiscrepancy,group,trim(prm%output(o)), & call results_writeDataset(dst%volumeDiscrepancy,group,trim(prm%output(o)), &
'volume discrepancy','m³') 'volume discrepancy','m³')
case('max_a_dot') case('max_dot_a')
call results_writeDataset(dst%relaxationrate_max,group,trim(prm%output(o)), & call results_writeDataset(dst%relaxationrate_max,group,trim(prm%output(o)), &
'maximum relaxation rate','m/s') 'maximum relaxation rate','m/s')
case('avg_a_dot') case('avg_dot_a')
call results_writeDataset(dst%relaxationrate_avg,group,trim(prm%output(o)), & call results_writeDataset(dst%relaxationrate_avg,group,trim(prm%output(o)), &
'average relaxation rate','m/s') 'average relaxation rate','m/s')
end select end select

15
src/phase.f90
View File

@ -61,9 +61,6 @@ module phase
grain grain
end type tDebugOptions 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 :: & type(tPlasticState), allocatable, dimension(:), public :: &
plasticState plasticState
type(tState), allocatable, dimension(:), public :: & type(tState), allocatable, dimension(:), public :: &
@ -283,16 +280,6 @@ module phase
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
end subroutine damage_anisobrittle_LiAndItsTangent 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 end interface
@ -413,7 +400,7 @@ end subroutine phase_init
subroutine phase_allocateState(state, & subroutine phase_allocateState(state, &
NEntries,sizeState,sizeDotState,sizeDeltaState) NEntries,sizeState,sizeDotState,sizeDeltaState)
class(tState), intent(out) :: & class(tState), intent(inout) :: &
state state
integer, intent(in) :: & integer, intent(in) :: &
NEntries, & NEntries, &

23
src/phase_damage.f90
View File

@ -11,7 +11,6 @@ submodule(phase) damage
enum, bind(c); enumerator :: & enum, bind(c); enumerator :: &
DAMAGE_UNDEFINED_ID, & DAMAGE_UNDEFINED_ID, &
DAMAGE_ISOBRITTLE_ID, & DAMAGE_ISOBRITTLE_ID, &
DAMAGE_ISODUCTILE_ID, &
DAMAGE_ANISOBRITTLE_ID DAMAGE_ANISOBRITTLE_ID
end enum end enum
@ -39,10 +38,6 @@ submodule(phase) damage
logical, dimension(:), allocatable :: mySources logical, dimension(:), allocatable :: mySources
end function isobrittle_init 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) module subroutine isobrittle_deltaState(C, Fe, ph, me)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
@ -59,10 +54,6 @@ submodule(phase) damage
S S
end subroutine anisobrittle_dotState 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) module subroutine anisobrittle_results(phase,group)
integer, intent(in) :: phase integer, intent(in) :: phase
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
@ -73,11 +64,6 @@ submodule(phase) damage
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
end subroutine isobrittle_results 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 end interface
contains contains
@ -131,7 +117,6 @@ module subroutine damage_init
if (damage_active) then if (damage_active) then
where(isobrittle_init() ) phase_damage = DAMAGE_ISOBRITTLE_ID where(isobrittle_init() ) phase_damage = DAMAGE_ISOBRITTLE_ID
where(isoductile_init() ) phase_damage = DAMAGE_ISODUCTILE_ID
where(anisobrittle_init()) phase_damage = DAMAGE_ANISOBRITTLE_ID where(anisobrittle_init()) phase_damage = DAMAGE_ANISOBRITTLE_ID
endif endif
@ -178,7 +163,7 @@ module function phase_f_phi(phi,co,ce) result(f)
en = material_phaseEntry(co,ce) en = material_phaseEntry(co,ce)
select case(phase_damage(ph)) 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 & f = 1.0_pReal &
- phi*damageState(ph)%state(1,en) - phi*damageState(ph)%state(1,en)
case default case default
@ -304,9 +289,6 @@ module subroutine damage_results(group,ph)
case (DAMAGE_ISOBRITTLE_ID) sourceType case (DAMAGE_ISOBRITTLE_ID) sourceType
call isobrittle_results(ph,group//'damage/') call isobrittle_results(ph,group//'damage/')
case (DAMAGE_ISODUCTILE_ID) sourceType
call isoductile_results(ph,group//'damage/')
case (DAMAGE_ANISOBRITTLE_ID) sourceType case (DAMAGE_ANISOBRITTLE_ID) sourceType
call anisobrittle_results(ph,group//'damage/') call anisobrittle_results(ph,group//'damage/')
@ -332,9 +314,6 @@ function phase_damage_collectDotState(ph,me) result(broken)
sourceType: select case (phase_damage(ph)) sourceType: select case (phase_damage(ph))
case (DAMAGE_ISODUCTILE_ID) sourceType
call isoductile_dotState(ph,me)
case (DAMAGE_ANISOBRITTLE_ID) sourceType case (DAMAGE_ANISOBRITTLE_ID) sourceType
call anisobrittle_dotState(mechanical_S(ph,me), ph,me) ! correct stress? call anisobrittle_dotState(mechanical_S(ph,me), ph,me) ! correct stress?

6
src/phase_damage_anisobrittle.f90
View File

@ -94,8 +94,8 @@ module function anisobrittle_init() result(mySources)
Nmembers = count(material_phaseID==p) Nmembers = count(material_phaseID==p)
call phase_allocateState(damageState(p),Nmembers,1,1,0) call phase_allocateState(damageState(p),Nmembers,1,1,0)
damageState(p)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) damageState(p)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal)
if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
end associate end associate
@ -110,7 +110,7 @@ end function anisobrittle_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state !> @brief
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine anisobrittle_dotState(S, ph,me) module subroutine anisobrittle_dotState(S, ph,me)

4
src/phase_damage_isobrittle.f90
View File

@ -66,8 +66,8 @@ module function isobrittle_init() result(mySources)
Nmembers = count(material_phaseID==ph) Nmembers = count(material_phaseID==ph)
call phase_allocateState(damageState(ph),Nmembers,1,1,1) call phase_allocateState(damageState(ph),Nmembers,1,1,1)
damageState(ph)%atol = src%get_asFloat('isobrittle_atol',defaultVal=1.0e-3_pReal) damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal)
if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
end associate end associate

127
src/phase_damage_isoductile.f90
View File

@ -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

2
src/phase_mechanical.f90
View File

@ -15,7 +15,6 @@ submodule(phase) mechanical
PLASTICITY_NONLOCAL_ID, & PLASTICITY_NONLOCAL_ID, &
KINEMATICS_UNDEFINED_ID, & KINEMATICS_UNDEFINED_ID, &
KINEMATICS_CLEAVAGE_OPENING_ID, & KINEMATICS_CLEAVAGE_OPENING_ID, &
KINEMATICS_SLIPPLANE_OPENING_ID, &
KINEMATICS_THERMAL_EXPANSION_ID KINEMATICS_THERMAL_EXPANSION_ID
end enum end enum
@ -274,7 +273,6 @@ module subroutine mechanical_init(materials,phases)
! initialize plasticity ! initialize plasticity
allocate(plasticState(phases%length)) allocate(plasticState(phases%length))
allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID) allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID)
allocate(phase_localPlasticity(phases%length), source=.true.)
call plastic_init() call plastic_init()

10
src/phase_mechanical_eigen.f90
View File

@ -13,10 +13,6 @@ submodule(phase:mechanical) eigen
logical, dimension(:), allocatable :: myKinematics logical, dimension(:), allocatable :: myKinematics
end function damage_anisobrittle_init 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) module function thermalexpansion_init(kinematics_length) result(myKinematics)
integer, intent(in) :: kinematics_length integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: myKinematics logical, dimension(:,:), allocatable :: myKinematics
@ -70,7 +66,6 @@ module subroutine eigendeformation_init(phases)
allocate(model_damage(phases%length), source = KINEMATICS_UNDEFINED_ID) allocate(model_damage(phases%length), source = KINEMATICS_UNDEFINED_ID)
where(damage_anisobrittle_init()) model_damage = KINEMATICS_cleavage_opening_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 end subroutine eigendeformation_init
@ -201,11 +196,6 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
Li = Li + my_Li Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS dLi_dS = dLi_dS + my_dLi_dS
active = .true. 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 end select
if(.not. active) return if(.not. active) return

184
src/phase_mechanical_eigen_slipplaneopening.f90
View File

@ -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

18
src/phase_mechanical_eigen_thermalexpansion.f90
View File

@ -29,7 +29,6 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
logical, dimension(:,:), allocatable :: myKinematics logical, dimension(:,:), allocatable :: myKinematics
integer :: Ninstances,p,i,k integer :: Ninstances,p,i,k
real(pReal), dimension(:), allocatable :: temp
class(tNode), pointer :: & class(tNode), pointer :: &
phases, & phases, &
phase, & phase, &
@ -61,21 +60,22 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
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) prm%A(1,1,1) = kinematic_type%get_asFloat('A_11')
temp = kinematic_type%get_as1dFloat('A_11') prm%A(1,1,2) = kinematic_type%get_asFloat('A_11,T',defaultVal=0.0_pReal)
prm%A(1,1,1:size(temp)) = temp prm%A(1,1,3) = kinematic_type%get_asFloat('A_11,T^2',defaultVal=0.0_pReal)
temp = kinematic_type%get_as1dFloat('A_33',defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp)) if (any(phase_lattice(p) == ['hP','tI'])) then
prm%A(3,3,1:size(temp)) = temp 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) do i=1, size(prm%A,3)
prm%A(1:3,1:3,i) = lattice_applyLatticeSymmetry33(prm%A(1:3,1:3,i),& prm%A(1:3,1:3,i) = lattice_applyLatticeSymmetry33(prm%A(1:3,1:3,i),phase_lattice(p))
phase%get_asString('lattice'))
enddo enddo
end associate end associate
endif endif
enddo enddo
enddo enddo
end function thermalexpansion_init end function thermalexpansion_init

5
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -247,9 +247,8 @@ module function plastic_dislotungsten_init() result(myPlasticity)
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:) stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:)
dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:) dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = 1.0e-2_pReal plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
! global alias if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:)
allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal) 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) allocate(dst%threshold_stress(prm%sum_N_sl,Nmembers), source=0.0_pReal)

7
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -436,15 +436,14 @@ module function plastic_dislotwin_init() result(myPlasticity)
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
stt%gamma_sl=>plasticState(ph)%state(startIndex:endIndex,:) stt%gamma_sl=>plasticState(ph)%state(startIndex:endIndex,:)
dot%gamma_sl=>plasticState(ph)%dotState(startIndex:endIndex,:) dot%gamma_sl=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = 1.0e-2_pReal plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
! global alias if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw endIndex = endIndex + prm%sum_N_tw
stt%f_tw=>plasticState(ph)%state(startIndex:endIndex,:) stt%f_tw=>plasticState(ph)%state(startIndex:endIndex,:)
dot%f_tw=>plasticState(ph)%dotState(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' if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tw'
startIndex = endIndex + 1 startIndex = endIndex + 1

2
src/phase_mechanical_plastic_isotropic.f90
View File

@ -139,8 +139,6 @@ module function plastic_isotropic_init() result(myPlasticity)
dot%gamma => plasticState(ph)%dotState(2,:) dot%gamma => plasticState(ph)%dotState(2,:)
plasticState(ph)%atol(2) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) 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' 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 end associate

2
src/phase_mechanical_plastic_kinehardening.f90
View File

@ -194,8 +194,6 @@ module function plastic_kinehardening_init() result(myPlasticity)
dot%accshear => plasticState(ph)%dotState(startIndex:endIndex,:) dot%accshear => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) 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' 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 o = plasticState(ph)%offsetDeltaState
startIndex = endIndex + 1 startIndex = endIndex + 1

178
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -81,11 +81,11 @@ submodule(phase:plastic) nonlocal
c_sol, & !< concentration of solid solution in atomic parts c_sol, & !< concentration of solid solution in atomic parts
p, & !< parameter for kinetic law (Kocks,Argon,Ashby) p, & !< parameter for kinetic law (Kocks,Argon,Ashby)
q, & !< 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 nu_a, & !< attack frequency in Hz
chi_surface, & !< transmissivity at free surface chi_surface, & !< transmissivity at free surface
chi_GB, & !< transmissivity at grain boundary (identified by different texture) 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_ed_mult, & !< factor that determines how much edge dislocations contribute to multiplication (0...1)
f_F, & f_F, &
f_ed, & f_ed, &
@ -112,7 +112,7 @@ submodule(phase:plastic) nonlocal
nonSchmid_pos, & nonSchmid_pos, &
nonSchmid_neg !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws) nonSchmid_neg !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws)
integer :: & integer :: &
sum_N_sl sum_N_sl = 0
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
colinearSystem !< colinear system to the active slip system (only valid for fcc!) colinearSystem !< colinear system to the active slip system (only valid for fcc!)
character(len=pStringLen), dimension(:), allocatable :: & character(len=pStringLen), dimension(:), allocatable :: &
@ -211,6 +211,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(geom(phases%length)) allocate(geom(phases%length))
allocate(param(phases%length)) allocate(param(phases%length))
@ -230,15 +231,14 @@ module function plastic_nonlocal_init() result(myPlasticity)
mech => phase%get('mechanical') mech => phase%get('mechanical')
pl => mech%get('plastic') pl => mech%get('plastic')
phase_localPlasticity(ph) = .not. pl%contains('nonlocal') plasticState(ph)%nonlocal = pl%get_asBool('flux',defaultVal=.True.)
#if defined (__GFORTRAN__) #if defined (__GFORTRAN__)
prm%output = output_as1dString(pl) prm%output = output_as1dString(pl)
#else #else
prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray) prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray)
#endif #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 ! This data is read in already in lattice
prm%mu = lattice_mu(ph) 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_significant = pl%get_asFloat('rho_significant')
prm%rho_min = pl%get_asFloat('rho_min', 0.0_pReal) 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%V_at = pl%get_asFloat('V_at')
prm%D_0 = pl%get_asFloat('D_0') prm%D_0 = pl%get_asFloat('D_0')
prm%Q_cl = pl%get_asFloat('Q_cl') prm%Q_cl = pl%get_asFloat('Q_cl')
prm%f_F = pl%get_asFloat('f_F') 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%w = pl%get_asFloat('w')
prm%Q_sol = pl%get_asFloat('Q_sol') prm%Q_sol = pl%get_asFloat('Q_sol')
prm%f_sol = pl%get_asFloat('f_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%p = pl%get_asFloat('p_sl')
prm%q = pl%get_asFloat('q_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') prm%nu_a = pl%get_asFloat('nu_a')
! ToDo: discuss logic ! 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%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 (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%B < 0.0_pReal) extmsg = trim(extmsg)//' B'
if (prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' 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%nu_a <= 0.0_pReal) extmsg = trim(extmsg)//' nu_a'
if (prm%w <= 0.0_pReal) extmsg = trim(extmsg)//' w' if (prm%w <= 0.0_pReal) extmsg = trim(extmsg)//' w'
if (prm%D_0 < 0.0_pReal) extmsg = trim(extmsg)//' D_0' 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_min < 0.0_pReal) extmsg = trim(extmsg)//' rho_min'
if (prm%rho_significant < 0.0_pReal) extmsg = trim(extmsg)//' rho_significant' 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%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%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' 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)) allocate(geom(ph)%V_0(Nmembers))
call storeGeometry(ph) call storeGeometry(ph)
plasticState(ph)%nonlocal = pl%get_asBool('nonlocal')
if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) & if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) &
call IO_error(212,ext_msg='IPneighborhood does not exist') 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) 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) 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) 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)) & if(any(plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) &
extmsg = trim(extmsg)//' atol_gamma' 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%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) 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) 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)) 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)) 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 enddo
endif 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 associate
end subroutine nonlocal_dependentState end subroutine nonlocal_dependentState
@ -994,13 +981,13 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
gdot !< shear rates gdot !< shear rates
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau, & !< current resolved shear stress 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) :: & real(pReal), dimension(param(ph)%sum_N_sl,2) :: &
rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) rhoDip, & !< current dipole dislocation densities (screw and edge dipoles)
dLower, & !< minimum stable dipole distance for edges and screws dLower, & !< minimum stable dipole distance for edges and screws
dUpper !< current maximum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws
real(pReal) :: & real(pReal) :: &
selfDiffusion !< self diffusion D_SD
if (timestep <= 0.0_pReal) then if (timestep <= 0.0_pReal) then
plasticState(ph)%dotState = 0.0_pReal 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) 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) 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 do s = 1,ns
tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,en) 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 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) dUpper = max(dUpper,dLower)
!****************************************************************************
!*** dislocation multiplication ! dislocation multiplication
rhoDotMultiplication = 0.0_pReal rhoDotMultiplication = 0.0_pReal
isBCC: if (phase_lattice(ph) == 'cI') then isBCC: if (phase_lattice(ph) == 'cI') then
forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal) 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 else isBCC
rhoDotMultiplication(:,1:4) = spread( & rhoDotMultiplication(:,1:4) = spread( &
(sum(abs(gdot(:,1:2)),2) * prm%f_ed_mult + sum(abs(gdot(:,3:4)),2)) & (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 endif isBCC
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) 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 do c = 1,2
rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%b_sl & rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
* ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile * ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile
@ -1102,7 +1081,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
enddo enddo
!*** athermal annihilation ! athermal annihilation
rhoDotAthermalAnnihilation = 0.0_pReal rhoDotAthermalAnnihilation = 0.0_pReal
forall (c=1:2) & forall (c=1:2) &
rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%b_sl & 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 * 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 rhoDotThermalAnnihilation = 0.0_pReal
selfDiffusion = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature)) D_SD = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature)) ! eq. 3.53
vClimb = prm%V_at * selfDiffusion * prm%mu & v_climb = D_SD * prm%mu * prm%V_at &
/ ( kB * Temperature * PI * (1.0_pReal-prm%nu) * (dUpper(:,1) + dLower(:,1))) / (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)) & 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) & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
- rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - 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 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) 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) 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) !*** calculate dislocation fluxes (only for nonlocal plasticity)
rhoDotFlux = 0.0_pReal rhoDotFlux = 0.0_pReal
if (.not. phase_localPlasticity(ph)) then if (plasticState(ph)%nonlocal) then
!*** check CFL (Courant-Friedrichs-Lewy) condition for flux !*** check CFL (Courant-Friedrichs-Lewy) condition for flux
if (any( abs(gdot) > 0.0_pReal & ! any active slip system ... 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) > 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 #ifdef DEBUG
if (debugConstitutive%extensive) then if (debugConstitutive%extensive) then
print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip 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 ', & print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', &
maxval(abs(v0), abs(gdot) > 0.0_pReal & 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))), & > IPvolume(ip,el) / maxval(IParea(:,ip,el))), &
' at a timestep of ',timestep ' at a timestep of ',timestep
print*, '<< CONST >> enforcing cutback !!!' 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 if (neighbor_e <= 0 .or. neighbor_i <= 0) then
!* FREE SURFACE !* 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) forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface)
elseif (neighbor_phase /= ph) then elseif (neighbor_phase /= ph) then
!* PHASE BOUNDARY !* PHASE BOUNDARY
!* If we encounter a different nonlocal phase at the neighbor, if (plasticState(neighbor_phase)%nonlocal .and. plasticState(ph)%nonlocal) &
!* 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)) &
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pReal forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pReal
elseif (prm%chi_GB >= 0.0_pReal) then elseif (prm%chi_GB >= 0.0_pReal) then
!* GRAIN BOUNDARY ! !* GRAIN BOUNDARY
!* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config)
if (any(dNeq(phase_orientation0(ph)%data(en)%asQuaternion(), & if (any(dNeq(phase_orientation0(ph)%data(en)%asQuaternion(), &
phase_orientation0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. & 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) forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB)
else else
!* GRAIN BOUNDARY ? !* GRAIN BOUNDARY ?
@ -1645,13 +1617,13 @@ end subroutine stateInit
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates kinetics !> @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) :: & integer, intent(in) :: &
c, & !< dislocation character (1:edge, 2:screw) c, & !< dislocation character (1:edge, 2:screw)
ph ph
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
Temperature !< temperature T !< T
real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: & real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: &
tau, & !< resolved external shear stress (without non Schmid effects) tau, & !< resolved external shear stress (without non Schmid effects)
tauNS, & !< resolved external shear stress (including 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 tauEff, & !< effective shear stress
tPeierls, & !< waiting time in front of a peierls barriers tPeierls, & !< waiting time in front of a peierls barriers
tSolidSolution, & !< waiting time in front of a solid solution obstacle tSolidSolution, & !< waiting time in front of a solid solution obstacle
vViscous, & !< viscous glide velocity
dtPeierls_dtau, & !< derivative with respect to resolved shear stress dtPeierls_dtau, & !< derivative with respect to resolved shear stress
dtSolidSolution_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 lambda_S, & !< mean free distance between two solid solution obstacles
meanfreepath_P, & !< mean free travel distance for dislocations between two Peierls barriers lambda_P, & !< mean free distance 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
activationVolume_P, & !< volume that needs to be activated to overcome barrier activationVolume_P, & !< volume that needs to be activated to overcome barrier
activationVolume_S, & !< 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_P, & !< energy that is needed to overcome barrier
activationEnergy_S, & !< energy that is needed to overcome barrier
criticalStress_P, & !< maximum obstacle strength criticalStress_P, & !< maximum obstacle strength
criticalStress_S, & !< maximum obstacle strength criticalStress_S !< maximum obstacle strength
mobility !< dislocation mobility
associate(prm => param(ph)) associate(prm => param(ph))
v = 0.0_pReal 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 if (abs(tau(s)) > tauThreshold(s)) then
!* Peierls contribution !* Peierls contribution
!* Effective stress includes non Schmid constributions tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s))
!* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity lambda_P = prm%b_sl(s)
tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) ! ensure that the effective stress is positive activationVolume_P = prm%w *prm%b_sl(s)**3.0_pReal
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)
criticalStress_P = prm%peierlsStress(s,c) criticalStress_P = prm%peierlsStress(s,c)
activationEnergy_P = criticalStress_P * activationVolume_P 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 & 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) * (1.0_pReal - tauRel_P**prm%p)**prm%q)
if (tauEff < criticalStress_P) then dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (kB * T) &
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), &
* (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) * tauRel_P**(prm%p-1.0_pReal) 0.0_pReal, &
else tauEff < criticalStress_P)
dtPeierls_dtau = 0.0_pReal
endif
!* Contribution from solid solution strengthening ! 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
tauEff = abs(tau(s)) - tauThreshold(s) tauEff = abs(tau(s)) - tauThreshold(s)
meanfreepath_S = prm%b_sl(s) / sqrt(prm%c_sol) lambda_S = prm%b_sl(s) / sqrt(prm%c_sol)
jumpWidth_S = prm%f_sol * prm%b_sl(s) activationVolume_S = prm%f_sol * prm%b_sl(s)**3.0_pReal / sqrt(prm%c_sol)
activationLength_S = prm%b_sl(s) / sqrt(prm%c_sol) criticalStress_S = prm%Q_sol / activationVolume_S
activationVolume_S = activationLength_S * jumpWidth_S * prm%b_sl(s) tauRel_S = min(1.0_pReal, tauEff / criticalStress_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
tSolidSolution = 1.0_pReal / prm%nu_a & tSolidSolution = 1.0_pReal / prm%nu_a &
* exp(activationEnergy_S / (kB * Temperature)* (1.0_pReal - tauRel_S**prm%p)**prm%q) * exp(prm%Q_sol / (kB * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q)
if (tauEff < criticalStress_S) then dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (kB * T) &
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), &
* (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal)* tauRel_S**(prm%p-1.0_pReal) 0.0_pReal, &
else tauEff < criticalStress_S)
dtSolidSolution_dtau = 0.0_pReal
endif
!* viscous glide velocity !* viscous glide velocity
tauEff = abs(tau(s)) - tauThreshold(s) 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)) & v(s) = sign(1.0_pReal,tau(s)) &
/ (tPeierls / meanfreepath_P + tSolidSolution / meanfreepath_S + 1.0_pReal / vViscous) / (tPeierls / lambda_P + tSolidSolution / lambda_S + prm%B /(prm%b_sl(s) * tauEff))
dv_dtau(s) = v(s)**2.0_pReal * (dtSolidSolution_dtau / meanfreepath_S + mobility /vViscous**2.0_pReal) 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 / meanfreepath_P dv_dtauNS(s) = v(s)**2.0_pReal * dtPeierls_dtau / lambda_P
endif endif
enddo enddo

4
src/phase_mechanical_plastic_phenopowerlaw.f90
View File

@ -247,7 +247,6 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
stt%xi_twin = spread(xi_0_tw, 2, Nmembers) stt%xi_twin = spread(xi_0_tw, 2, Nmembers)
dot%xi_twin => plasticState(ph)%dotState(startIndex:endIndex,:) dot%xi_twin => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) 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 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl 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,:) dot%gamma_slip => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) 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' 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 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw endIndex = endIndex + prm%sum_N_tw
stt%gamma_twin => plasticState(ph)%state (startIndex:endIndex,:) stt%gamma_twin => plasticState(ph)%state (startIndex:endIndex,:)
dot%gamma_twin => plasticState(ph)%dotState(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) 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 end associate

5
src/prec.f90
View File

@ -51,10 +51,7 @@ module prec
end type end type
type, extends(tState) :: tPlasticState type, extends(tState) :: tPlasticState
logical :: & logical :: nonlocal = .false.
nonlocal = .false.
real(pReal), pointer, dimension(:,:) :: &
slipRate !< slip rate
end type end type
type :: tSourceState type :: tSourceState