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

Merge remote-tracking branch 'origin/development' into stress-ramp-loadcase

This commit is contained in:
Martin Diehl 2020-09-25 15:05:45 +02:00
parent 79d672f4a7 20f5278f51
commit 8081ed3693
28 changed files with 791 additions and 622 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitlab-ci.yml
View File

@ -344,6 +344,7 @@ Phenopowerlaw_singleSlip:
Pytest_grid: Pytest_grid:
stage: grid stage: grid
script: script:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
- cd pytest - cd pytest
- pytest - pytest
except: except:

2
PRIVATE

@ -1 +1 @@
Subproject commit 2550bf655aebe61ce459718323ba38b1ac205a7f Subproject commit dc568df60e36b659d9a1f84ac93fd4abb1b8fe3c

2
VERSION
View File

@ -1 +1 @@
v3.0.0-alpha-321-g0f6495430 v3.0.0-alpha-335-gdb8f6400f

6
python/damask/__init__.py
View File

@ -20,6 +20,12 @@ from ._result import Result # noqa
from ._geom import Geom # noqa from ._geom import Geom # noqa
from ._material import Material # noqa from ._material import Material # noqa
from . import solver # noqa from . import solver # noqa
from . import util # noqa
from . import seeds # noqa
from . import grid_filters # noqa
from . import mechanics # noqa
# deprecated # deprecated
Environment = _ Environment = _

24
python/damask/_geom.py
View File

@ -212,7 +212,7 @@ class Geom:
return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights) return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights)
@staticmethod @staticmethod
def from_Laguerre_tessellation(grid,size,seeds,weights,periodic=True): def from_Laguerre_tessellation(grid,size,seeds,weights,material=None,periodic=True):
""" """
Generate geometry from Laguerre tessellation. Generate geometry from Laguerre tessellation.
@ -226,6 +226,9 @@ class Geom:
Position of the seed points in meter. All points need to lay within the box. Position of the seed points in meter. All points need to lay within the box.
weights : numpy.ndarray of shape (seeds.shape[0]) weights : numpy.ndarray of shape (seeds.shape[0])
Weights of the seeds. Setting all weights to 1.0 gives a standard Voronoi tessellation. Weights of the seeds. Setting all weights to 1.0 gives a standard Voronoi tessellation.
material : numpy.ndarray of shape (seeds.shape[0]), optional
Material ID of the seeds. Defaults to None, in which case materials are
consecutively numbered.
periodic : Boolean, optional periodic : Boolean, optional
Perform a periodic tessellation. Defaults to True. Perform a periodic tessellation. Defaults to True.
@ -245,22 +248,22 @@ class Geom:
result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords]) result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords])
pool.close() pool.close()
pool.join() pool.join()
material = np.array(result.get()) material_ = np.array(result.get())
if periodic: if periodic:
material = material.reshape(grid*3) material_ = material_.reshape(grid*3)
material = material[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0] material_ = material_[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0]
else: else:
material = material.reshape(grid) material_ = material_.reshape(grid)
return Geom(material = material+1, return Geom(material = material_+1 if material is None else material[material_],
size = size, size = size,
comments = util.execution_stamp('Geom','from_Laguerre_tessellation'), comments = util.execution_stamp('Geom','from_Laguerre_tessellation'),
) )
@staticmethod @staticmethod
def from_Voronoi_tessellation(grid,size,seeds,periodic=True): def from_Voronoi_tessellation(grid,size,seeds,material=None,periodic=True):
""" """
Generate geometry from Voronoi tessellation. Generate geometry from Voronoi tessellation.
@ -272,15 +275,18 @@ class Geom:
Physical size of the geometry in meter. Physical size of the geometry in meter.
seeds : numpy.ndarray of shape (:,3) seeds : numpy.ndarray of shape (:,3)
Position of the seed points in meter. All points need to lay within the box. Position of the seed points in meter. All points need to lay within the box.
material : numpy.ndarray of shape (seeds.shape[0]), optional
Material ID of the seeds. Defaults to None, in which case materials are
consecutively numbered.
periodic : Boolean, optional periodic : Boolean, optional
Perform a periodic tessellation. Defaults to True. Perform a periodic tessellation. Defaults to True.
""" """
coords = grid_filters.cell_coord0(grid,size).reshape(-1,3) coords = grid_filters.cell_coord0(grid,size).reshape(-1,3)
KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds) KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
devNull,material = KDTree.query(coords) devNull,material_ = KDTree.query(coords)
return Geom(material = material.reshape(grid)+1, return Geom(material = (material_+1 if material is None else material[material_]).reshape(grid),
size = size, size = size,
comments = util.execution_stamp('Geom','from_Voronoi_tessellation'), comments = util.execution_stamp('Geom','from_Voronoi_tessellation'),
) )

113
python/damask/seeds.py Normal file
View File

@ -0,0 +1,113 @@
from scipy import spatial as _spatial
import numpy as _np
from . import util
from . import grid_filters
def from_random(size,N_seeds,grid=None,seed=None):
"""
Random seeding in space.
Parameters
----------
size : numpy.ndarray of shape (3)
Physical size of the seeding domain.
N_seeds : int
Number of seeds.
grid : numpy.ndarray of shape (3), optional.
If given, ensures that all seeds initiate one grain if using a
standard Voronoi tessellation.
seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
A seed to initialize the BitGenerator. Defaults to None.
If None, then fresh, unpredictable entropy will be pulled from the OS.
"""
rng = _np.random.default_rng(seed)
if grid is None:
coords = rng.random((N_seeds,3)) * size
else:
grid_coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
coords = grid_coords[rng.choice(_np.prod(grid),N_seeds, replace=False)] \
+ _np.broadcast_to(size/grid,(N_seeds,3))*(rng.random((N_seeds,3))*.5-.25) # wobble without leaving grid
return coords
def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,seed=None):
"""
Seeding in space according to a Poisson disc distribution.
Parameters
----------
size : numpy.ndarray of shape (3)
Physical size of the seeding domain.
N_seeds : int
Number of seeds.
N_candidates : int
Number of candidates to consider for finding best candidate.
distance : float
Minimum acceptable distance to other seeds.
periodic : boolean, optional
Calculate minimum distance for periodically repeated grid.
seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
A seed to initialize the BitGenerator. Defaults to None.
If None, then fresh, unpredictable entropy will be pulled from the OS.
"""
rng = _np.random.default_rng(seed)
coords = _np.empty((N_seeds,3))
coords[0] = rng.random(3) * size
i = 1
progress = util._ProgressBar(N_seeds+1,'',50)
while i < N_seeds:
candidates = rng.random((N_candidates,3))*_np.broadcast_to(size,(N_candidates,3))
tree = _spatial.cKDTree(coords[:i],boxsize=size) if periodic else \
_spatial.cKDTree(coords[:i])
distances, dev_null = tree.query(candidates)
best = distances.argmax()
if distances[best] > distance: # require minimum separation
coords[i] = candidates[best] # maximum separation to existing point cloud
i += 1
progress.update(i)
return coords
def from_geom(geom,selection=None,invert=False,average=False,periodic=True):
"""
Create seed from existing geometry description.
Parameters
----------
geom : damask.Geom
Geometry, from which the material IDs are used as seeds.
selection : iterable of integers, optional
Material IDs to consider.
invert : boolean, false
Do not consider the material IDs given in selection. Defaults to False.
average : boolean, optional
Seed corresponds to center of gravity of material ID cloud.
periodic : boolean, optional
Center of gravity with periodic boundaries.
"""
material = geom.material.reshape((-1,1),order='F')
mask = _np.full(geom.grid.prod(),True,dtype=bool) if selection is None else \
_np.isin(material,selection,invert=invert)
coords = grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3,order='F')
if not average:
return (coords[mask],material[mask])
else:
materials = _np.unique(material[mask])
coords_ = _np.zeros((materials.size,3),dtype=float)
for i,mat in enumerate(materials):
pc = (2*_np.pi*coords[material[:,0]==mat,:]-geom.origin)/geom.size
coords_[i] = geom.origin + geom.size / 2 / _np.pi * (_np.pi +
_np.arctan2(-_np.average(_np.sin(pc),axis=0),
-_np.average(_np.cos(pc),axis=0))) \
if periodic else \
_np.average(coords[material[:,0]==mat,:],axis=0)
return (coords_,materials)

17
python/tests/test_Geom.py
View File

@ -83,6 +83,13 @@ class TestGeom:
with pytest.raises(ValueError): with pytest.raises(ValueError):
Geom.load(tmpdir/'no_materialpoint.vtr') Geom.load(tmpdir/'no_materialpoint.vtr')
def test_invalid_material(self):
with pytest.raises(TypeError):
Geom(np.zeros((3,3,3),dtype='complex'),np.ones(3))
def test_cast_to_int(self):
g = Geom(np.zeros((3,3,3)),np.ones(3))
assert g.material.dtype in np.sctypes['int']
@pytest.mark.parametrize('compress',[True,False]) @pytest.mark.parametrize('compress',[True,False])
def test_compress(self,default,tmpdir,compress): def test_compress(self,default,tmpdir,compress):
@ -324,8 +331,8 @@ class TestGeom:
size = np.random.random(3) + 1.0 size = np.random.random(3) + 1.0
N_seeds= np.random.randint(10,30) N_seeds= np.random.randint(10,30)
seeds = np.random.rand(N_seeds,3) * np.broadcast_to(size,(N_seeds,3)) seeds = np.random.rand(N_seeds,3) * np.broadcast_to(size,(N_seeds,3))
Voronoi = Geom.from_Voronoi_tessellation( grid,size,seeds, periodic) Voronoi = Geom.from_Voronoi_tessellation( grid,size,seeds, np.arange(N_seeds)+5,periodic)
Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(N_seeds),periodic) Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(N_seeds),np.arange(N_seeds)+5,periodic)
assert geom_equal(Laguerre,Voronoi) assert geom_equal(Laguerre,Voronoi)
@ -337,7 +344,7 @@ class TestGeom:
weights= np.full((N_seeds),-np.inf) weights= np.full((N_seeds),-np.inf)
ms = np.random.randint(1, N_seeds+1) ms = np.random.randint(1, N_seeds+1)
weights[ms-1] = np.random.random() weights[ms-1] = np.random.random()
Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,weights,np.random.random()>0.5) Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,weights,periodic=np.random.random()>0.5)
assert np.all(Laguerre.material == ms) assert np.all(Laguerre.material == ms)
@ -349,7 +356,7 @@ class TestGeom:
material = np.ones(grid) material = np.ones(grid)
material[:,grid[1]//2:,:] = 2 material[:,grid[1]//2:,:] = 2
if approach == 'Laguerre': if approach == 'Laguerre':
geom = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(2),np.random.random()>0.5) geom = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(2),periodic=np.random.random()>0.5)
elif approach == 'Voronoi': elif approach == 'Voronoi':
geom = Geom.from_Voronoi_tessellation(grid,size,seeds, np.random.random()>0.5) geom = Geom.from_Voronoi_tessellation(grid,size,seeds, periodic=np.random.random()>0.5)
assert np.all(geom.material == material) assert np.all(geom.material == material)

36
python/tests/test_seeds.py Normal file
View File

@ -0,0 +1,36 @@
import pytest
import numpy as np
from scipy.spatial import cKDTree
from damask import seeds
from damask import Geom
class TestSeeds:
@pytest.mark.parametrize('grid',[None,np.ones(3,dtype='i')*10])
def test_from_random(self,grid):
N_seeds = np.random.randint(30,300)
size = np.ones(3) + np.random.random(3)
coords = seeds.from_random(size,N_seeds,grid)
assert (0<=coords).all() and (coords<size).all()
@pytest.mark.parametrize('periodic',[True,False])
def test_from_Poisson_disc(self,periodic):
N_seeds = np.random.randint(30,300)
N_candidates = N_seeds//15
distance = np.random.random()
size = np.ones(3)*distance*N_seeds
coords = seeds.from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=periodic)
min_dists, _ = cKDTree(coords,boxsize=size).query(coords, 2) if periodic else \
cKDTree(coords).query(coords, 2)
assert (0<= coords).all() and (coords<size).all() and np.min(min_dists[:,1])>=distance
def test_from_geom(self):
grid = np.random.randint(10,20,3)
N_seeds = np.random.randint(30,300)
size = np.ones(3) + np.random.random(3)
coords = seeds.from_random(size,N_seeds,grid)
geom_1 = Geom.from_Voronoi_tessellation(grid,size,coords)
coords,material = seeds.from_geom(geom_1)
geom_2 = Geom.from_Voronoi_tessellation(grid,size,coords,material)
assert (geom_2.material==geom_1.material).all()

68
src/constitutive_plastic_disloTungsten.f90
View File

@ -17,18 +17,18 @@ submodule(constitutive:constitutive_plastic) plastic_disloTungsten
D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient
Q_cl = 1.0_pReal !< activation energy for dislocation climb Q_cl = 1.0_pReal !< activation energy for dislocation climb
real(pReal), allocatable, dimension(:) :: & real(pReal), allocatable, dimension(:) :: &
b_sl, & !< magnitude of burgers vector [m] b_sl, & !< magnitude of Burgers vector [m]
D_a, & D_a, &
i_sl, & !< Adj. parameter for distance between 2 forest dislocations i_sl, & !< Adj. parameter for distance between 2 forest dislocations
atomicVolume, & !< factor to calculate atomic volume f_at, & !< factor to calculate atomic volume
tau_0, & !< Peierls stress tau_Peierls, & !< Peierls stress
!* mobility law parameters !* mobility law parameters
delta_F, & !< activation energy for glide [J] Q_s, & !< activation energy for glide [J]
v0, & !< dislocation velocity prefactor [m/s] v_0, & !< dislocation velocity prefactor [m/s]
p, & !< p-exponent in glide velocity p, & !< p-exponent in glide velocity
q, & !< q-exponent in glide velocity q, & !< q-exponent in glide velocity
B, & !< friction coefficient B, & !< friction coefficient
kink_height, & !< height of the kink pair h, & !< height of the kink pair
w, & !< width of the kink pair w, & !< width of the kink pair
omega !< attempt frequency for kink pair nucleation omega !< attempt frequency for kink pair nucleation
real(pReal), allocatable, dimension(:,:) :: & real(pReal), allocatable, dimension(:,:) :: &
@ -142,7 +142,7 @@ module function plastic_disloTungsten_init() result(myPlasticity)
phase%get_asFloat('c/a',defaultVal=0.0_pReal)) phase%get_asFloat('c/a',defaultVal=0.0_pReal))
if(trim(phase%get_asString('lattice')) == 'bcc') then if(trim(phase%get_asString('lattice')) == 'bcc') then
a = pl%get_asFloats('nonSchmid_coefficients',defaultVal = emptyRealArray) a = pl%get_asFloats('a_nonSchmid',defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1) prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1) prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
else else
@ -158,17 +158,17 @@ module function plastic_disloTungsten_init() result(myPlasticity)
rho_mob_0 = pl%get_asFloats('rho_mob_0', requiredSize=size(N_sl)) rho_mob_0 = pl%get_asFloats('rho_mob_0', requiredSize=size(N_sl))
rho_dip_0 = pl%get_asFloats('rho_dip_0', requiredSize=size(N_sl)) rho_dip_0 = pl%get_asFloats('rho_dip_0', requiredSize=size(N_sl))
prm%v0 = pl%get_asFloats('v_0', requiredSize=size(N_sl)) prm%v_0 = pl%get_asFloats('v_0', requiredSize=size(N_sl))
prm%b_sl = pl%get_asFloats('b_sl', requiredSize=size(N_sl)) prm%b_sl = pl%get_asFloats('b_sl', requiredSize=size(N_sl))
prm%delta_F = pl%get_asFloats('Q_s', requiredSize=size(N_sl)) prm%Q_s = pl%get_asFloats('Q_s', requiredSize=size(N_sl))
prm%i_sl = pl%get_asFloats('i_sl', requiredSize=size(N_sl)) prm%i_sl = pl%get_asFloats('i_sl', requiredSize=size(N_sl))
prm%tau_0 = pl%get_asFloats('tau_peierls', requiredSize=size(N_sl)) prm%tau_Peierls = pl%get_asFloats('tau_Peierls', requiredSize=size(N_sl))
prm%p = pl%get_asFloats('p_sl', requiredSize=size(N_sl), & prm%p = pl%get_asFloats('p_sl', requiredSize=size(N_sl), &
defaultVal=[(1.0_pReal,i=1,size(N_sl))]) defaultVal=[(1.0_pReal,i=1,size(N_sl))])
prm%q = pl%get_asFloats('q_sl', requiredSize=size(N_sl), & prm%q = pl%get_asFloats('q_sl', requiredSize=size(N_sl), &
defaultVal=[(1.0_pReal,i=1,size(N_sl))]) defaultVal=[(1.0_pReal,i=1,size(N_sl))])
prm%kink_height = pl%get_asFloats('h', requiredSize=size(N_sl)) prm%h = pl%get_asFloats('h', requiredSize=size(N_sl))
prm%w = pl%get_asFloats('w', requiredSize=size(N_sl)) prm%w = pl%get_asFloats('w', requiredSize=size(N_sl))
prm%omega = pl%get_asFloats('omega', requiredSize=size(N_sl)) prm%omega = pl%get_asFloats('omega', requiredSize=size(N_sl))
prm%B = pl%get_asFloats('B', requiredSize=size(N_sl)) prm%B = pl%get_asFloats('B', requiredSize=size(N_sl))
@ -176,7 +176,7 @@ module function plastic_disloTungsten_init() result(myPlasticity)
prm%D = pl%get_asFloat('D') prm%D = pl%get_asFloat('D')
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%atomicVolume = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal prm%f_at = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal
prm%D_a = pl%get_asFloat('D_a') * prm%b_sl prm%D_a = pl%get_asFloat('D_a') * prm%b_sl
prm%dipoleformation = pl%get_asBool('dipole_formation_factor', defaultVal = .true.) prm%dipoleformation = pl%get_asBool('dipole_formation_factor', defaultVal = .true.)
@ -186,16 +186,16 @@ module function plastic_disloTungsten_init() result(myPlasticity)
rho_dip_0 = math_expand(rho_dip_0, N_sl) rho_dip_0 = math_expand(rho_dip_0, N_sl)
prm%q = math_expand(prm%q, N_sl) prm%q = math_expand(prm%q, N_sl)
prm%p = math_expand(prm%p, N_sl) prm%p = math_expand(prm%p, N_sl)
prm%delta_F = math_expand(prm%delta_F, N_sl) prm%Q_s = math_expand(prm%Q_s, N_sl)
prm%b_sl = math_expand(prm%b_sl, N_sl) prm%b_sl = math_expand(prm%b_sl, N_sl)
prm%kink_height = math_expand(prm%kink_height, N_sl) prm%h = math_expand(prm%h, N_sl)
prm%w = math_expand(prm%w, N_sl) prm%w = math_expand(prm%w, N_sl)
prm%omega = math_expand(prm%omega, N_sl) prm%omega = math_expand(prm%omega, N_sl)
prm%tau_0 = math_expand(prm%tau_0, N_sl) prm%tau_Peierls = math_expand(prm%tau_Peierls, N_sl)
prm%v0 = math_expand(prm%v0, N_sl) prm%v_0 = math_expand(prm%v_0, N_sl)
prm%B = math_expand(prm%B, N_sl) prm%B = math_expand(prm%B, N_sl)
prm%i_sl = math_expand(prm%i_sl, N_sl) prm%i_sl = math_expand(prm%i_sl, N_sl)
prm%atomicVolume = math_expand(prm%atomicVolume, N_sl) prm%f_at = math_expand(prm%f_at, N_sl)
prm%D_a = math_expand(prm%D_a, N_sl) prm%D_a = math_expand(prm%D_a, N_sl)
! sanity checks ! sanity checks
@ -203,17 +203,17 @@ module function plastic_disloTungsten_init() result(myPlasticity)
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 (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0' if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0'
if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0' if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0'
if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0' if (any(prm%v_0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0'
if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl'
if (any(prm%delta_F <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s' if (any(prm%Q_s <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s'
if (any(prm%tau_0 < 0.0_pReal)) extmsg = trim(extmsg)//' tau_peierls' if (any(prm%tau_Peierls < 0.0_pReal)) extmsg = trim(extmsg)//' tau_Peierls'
if (any(prm%D_a <= 0.0_pReal)) extmsg = trim(extmsg)//' D_a or b_sl' if (any(prm%D_a <= 0.0_pReal)) extmsg = trim(extmsg)//' D_a or b_sl'
if (any(prm%atomicVolume <= 0.0_pReal)) extmsg = trim(extmsg)//' f_at or b_sl' if (any(prm%f_at <= 0.0_pReal)) extmsg = trim(extmsg)//' f_at or b_sl'
else slipActive else slipActive
rho_mob_0= emptyRealArray; rho_dip_0 = emptyRealArray rho_mob_0= emptyRealArray; rho_dip_0 = emptyRealArray
allocate(prm%b_sl,prm%D_a,prm%i_sl,prm%atomicVolume,prm%tau_0, & allocate(prm%b_sl,prm%D_a,prm%i_sl,prm%f_at,prm%tau_Peierls, &
prm%delta_F,prm%v0,prm%p,prm%q,prm%B,prm%kink_height,prm%w,prm%omega, & prm%Q_s,prm%v_0,prm%p,prm%q,prm%B,prm%h,prm%w,prm%omega, &
source = emptyRealArray) source = emptyRealArray)
allocate(prm%forestProjection(0,0)) allocate(prm%forestProjection(0,0))
allocate(prm%h_sl_sl (0,0)) allocate(prm%h_sl_sl (0,0))
@ -354,7 +354,7 @@ module subroutine plastic_disloTungsten_dotState(Mp,T,instance,of)
dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation
0.0_pReal, & 0.0_pReal, &
prm%dipoleformation) prm%dipoleformation)
v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%atomicVolume/(2.0_pReal*pi*kB*T)) & v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%f_at/(2.0_pReal*pi*kB*T)) &
* (1.0_pReal/(dip_distance+prm%D_a)) * (1.0_pReal/(dip_distance+prm%D_a))
dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,of))/(dip_distance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency? dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,of))/(dip_distance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency?
end where end where
@ -477,12 +477,12 @@ pure subroutine kinetics(Mp,T,instance,of, &
if (present(tau_pos_out)) tau_pos_out = tau_pos if (present(tau_pos_out)) tau_pos_out = tau_pos
if (present(tau_neg_out)) tau_neg_out = tau_neg if (present(tau_neg_out)) tau_neg_out = tau_neg
associate(BoltzmannRatio => prm%delta_F/(kB*T), & associate(BoltzmannRatio => prm%Q_s/(kB*T), &
dot_gamma_0 => stt%rho_mob(:,of)*prm%b_sl*prm%v0, & dot_gamma_0 => stt%rho_mob(:,of)*prm%b_sl*prm%v_0, &
effectiveLength => dst%Lambda_sl(:,of) - prm%w) effectiveLength => dst%Lambda_sl(:,of) - prm%w)
significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check)
StressRatio = (abs(tau_pos)-dst%threshold_stress(:,of))/prm%tau_0 StressRatio = (abs(tau_pos)-dst%threshold_stress(:,of))/prm%tau_Peierls
StressRatio_p = StressRatio** prm%p StressRatio_p = StressRatio** prm%p
StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
@ -490,7 +490,7 @@ pure subroutine kinetics(Mp,T,instance,of, &
t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength)
t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos)
vel = prm%kink_height/(t_n + t_k) vel = prm%h/(t_n + t_k)
dot_gamma_pos = dot_gamma_0 * sign(vel,tau_pos) * 0.5_pReal dot_gamma_pos = dot_gamma_0 * sign(vel,tau_pos) * 0.5_pReal
else where significantPositiveTau else where significantPositiveTau
@ -500,10 +500,10 @@ pure subroutine kinetics(Mp,T,instance,of, &
if (present(ddot_gamma_dtau_pos)) then if (present(ddot_gamma_dtau_pos)) then
significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check) significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check)
dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
* (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0 * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls
dtk = -1.0_pReal * t_k / tau_pos dtk = -1.0_pReal * t_k / tau_pos
dvel = -1.0_pReal * prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal
ddot_gamma_dtau_pos = dot_gamma_0 * dvel* 0.5_pReal ddot_gamma_dtau_pos = dot_gamma_0 * dvel* 0.5_pReal
else where significantPositiveTau2 else where significantPositiveTau2
@ -512,7 +512,7 @@ pure subroutine kinetics(Mp,T,instance,of, &
endif endif
significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check) significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check)
StressRatio = (abs(tau_neg)-dst%threshold_stress(:,of))/prm%tau_0 StressRatio = (abs(tau_neg)-dst%threshold_stress(:,of))/prm%tau_Peierls
StressRatio_p = StressRatio** prm%p StressRatio_p = StressRatio** prm%p
StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
@ -520,7 +520,7 @@ pure subroutine kinetics(Mp,T,instance,of, &
t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength)
t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos)
vel = prm%kink_height/(t_n + t_k) vel = prm%h/(t_n + t_k)
dot_gamma_neg = dot_gamma_0 * sign(vel,tau_neg) * 0.5_pReal dot_gamma_neg = dot_gamma_0 * sign(vel,tau_neg) * 0.5_pReal
else where significantNegativeTau else where significantNegativeTau
@ -530,10 +530,10 @@ pure subroutine kinetics(Mp,T,instance,of, &
if (present(ddot_gamma_dtau_neg)) then if (present(ddot_gamma_dtau_neg)) then
significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check) significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check)
dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
* (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0 * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls
dtk = -1.0_pReal * t_k / tau_neg dtk = -1.0_pReal * t_k / tau_neg
dvel = -1.0_pReal * prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal
ddot_gamma_dtau_neg = dot_gamma_0 * dvel * 0.5_pReal ddot_gamma_dtau_neg = dot_gamma_0 * dvel * 0.5_pReal
else where significantNegativeTau2 else where significantNegativeTau2

128
src/constitutive_plastic_dislotwin.f90
View File

@ -16,38 +16,38 @@ submodule(constitutive:constitutive_plastic) plastic_dislotwin
real(pReal) :: & real(pReal) :: &
mu = 1.0_pReal, & !< equivalent shear modulus mu = 1.0_pReal, & !< equivalent shear modulus
nu = 1.0_pReal, & !< equivalent shear Poisson's ratio nu = 1.0_pReal, & !< equivalent shear Poisson's ratio
D0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient
Qsd = 1.0_pReal, & !< activation energy for dislocation climb Q_cl = 1.0_pReal, & !< activation energy for dislocation climb
omega = 1.0_pReal, & !< frequency factor for dislocation climb omega = 1.0_pReal, & !< frequency factor for dislocation climb
D = 1.0_pReal, & !< grain size D = 1.0_pReal, & !< grain size
p_sb = 1.0_pReal, & !< p-exponent in shear band velocity p_sb = 1.0_pReal, & !< p-exponent in shear band velocity
q_sb = 1.0_pReal, & !< q-exponent in shear band velocity q_sb = 1.0_pReal, & !< q-exponent in shear band velocity
CEdgeDipMinDistance = 1.0_pReal, & !< adjustment parameter to calculate minimum dipole distance D_a = 1.0_pReal, & !< adjustment parameter to calculate minimum dipole distance
i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning
tau_0 = 1.0_pReal, & !< strength due to elements in solid solution tau_0 = 1.0_pReal, & !< strength due to elements in solid solution
L_tw = 1.0_pReal, & !< Length of twin nuclei in Burgers vectors L_tw = 1.0_pReal, & !< Length of twin nuclei in Burgers vectors
L_tr = 1.0_pReal, & !< Length of trans nuclei in Burgers vectors L_tr = 1.0_pReal, & !< Length of trans nuclei in Burgers vectors
xc_twin = 1.0_pReal, & !< critical distance for formation of twin nucleus x_c_tw = 1.0_pReal, & !< critical distance for formation of twin nucleus
xc_trans = 1.0_pReal, & !< critical distance for formation of trans nucleus x_c_tr = 1.0_pReal, & !< critical distance for formation of trans nucleus
V_cs = 1.0_pReal, & !< cross slip volume V_cs = 1.0_pReal, & !< cross slip volume
sbResistance = 1.0_pReal, & !< value for shearband resistance xi_sb = 1.0_pReal, & !< value for shearband resistance
sbVelocity = 1.0_pReal, & !< value for shearband velocity_0 v_sb = 1.0_pReal, & !< value for shearband velocity_0
E_sb = 1.0_pReal, & !< activation energy for shear bands E_sb = 1.0_pReal, & !< activation energy for shear bands
SFE_0K = 1.0_pReal, & !< stacking fault energy at zero K Gamma_sf_0K = 1.0_pReal, & !< stacking fault energy at zero K
dSFE_dT = 1.0_pReal, & !< temperature dependence of stacking fault energy dGamma_sf_dT = 1.0_pReal, & !< temperature dependence of stacking fault energy
gamma_fcc_hex = 1.0_pReal, & !< Free energy difference between austensite and martensite delta_G = 1.0_pReal, & !< Free energy difference between austensite and martensite
i_tr = 1.0_pReal, & !< adjustment parameter to calculate MFP for transformation i_tr = 1.0_pReal, & !< adjustment parameter to calculate MFP for transformation
h = 1.0_pReal !< Stack height of hex nucleus h = 1.0_pReal !< Stack height of hex nucleus
real(pReal), allocatable, dimension(:) :: & real(pReal), allocatable, dimension(:) :: &
b_sl, & !< absolute length of burgers vector [m] for each slip system b_sl, & !< absolute length of Burgers vector [m] for each slip system
b_tw, & !< absolute length of burgers vector [m] for each twin system b_tw, & !< absolute length of Burgers vector [m] for each twin system
b_tr, & !< absolute length of burgers vector [m] for each transformation system b_tr, & !< absolute length of Burgers vector [m] for each transformation system
Delta_F,& !< activation energy for glide [J] for each slip system Q_s,& !< activation energy for glide [J] for each slip system
v0, & !< dislocation velocity prefactor [m/s] for each slip system v_0, & !< dislocation velocity prefactor [m/s] for each slip system
dot_N_0_tw, & !< twin nucleation rate [1/m³s] for each twin system dot_N_0_tw, & !< twin nucleation rate [1/m³s] for each twin system
dot_N_0_tr, & !< trans nucleation rate [1/m³s] for each trans system dot_N_0_tr, & !< trans nucleation rate [1/m³s] for each trans system
t_tw, & !< twin thickness [m] for each twin system t_tw, & !< twin thickness [m] for each twin system
CLambdaSlip, & !< Adj. parameter for distance between 2 forest dislocations for each slip system i_sl, & !< Adj. parameter for distance between 2 forest dislocations for each slip system
t_tr, & !< martensite lamellar thickness [m] for each trans system and instance t_tr, & !< martensite lamellar thickness [m] for each trans system and instance
p, & !< p-exponent in glide velocity p, & !< p-exponent in glide velocity
q, & !< q-exponent in glide velocity q, & !< q-exponent in glide velocity
@ -209,23 +209,23 @@ module function plastic_dislotwin_init() result(myPlasticity)
rho_mob_0 = pl%get_asFloats('rho_mob_0', requiredSize=size(N_sl)) rho_mob_0 = pl%get_asFloats('rho_mob_0', requiredSize=size(N_sl))
rho_dip_0 = pl%get_asFloats('rho_dip_0', requiredSize=size(N_sl)) rho_dip_0 = pl%get_asFloats('rho_dip_0', requiredSize=size(N_sl))
prm%v0 = pl%get_asFloats('v_0', requiredSize=size(N_sl)) prm%v_0 = pl%get_asFloats('v_0', requiredSize=size(N_sl))
prm%b_sl = pl%get_asFloats('b_sl', requiredSize=size(N_sl)) prm%b_sl = pl%get_asFloats('b_sl', requiredSize=size(N_sl))
prm%Delta_F = pl%get_asFloats('Q_s', requiredSize=size(N_sl)) prm%Q_s = pl%get_asFloats('Q_s', requiredSize=size(N_sl))
prm%CLambdaSlip = pl%get_asFloats('i_sl', requiredSize=size(N_sl)) prm%i_sl = pl%get_asFloats('i_sl', requiredSize=size(N_sl))
prm%p = pl%get_asFloats('p_sl', requiredSize=size(N_sl)) prm%p = pl%get_asFloats('p_sl', requiredSize=size(N_sl))
prm%q = pl%get_asFloats('q_sl', requiredSize=size(N_sl)) prm%q = pl%get_asFloats('q_sl', requiredSize=size(N_sl))
prm%B = pl%get_asFloats('B', requiredSize=size(N_sl), & prm%B = pl%get_asFloats('B', requiredSize=size(N_sl), &
defaultVal=[(0.0_pReal, i=1,size(N_sl))]) defaultVal=[(0.0_pReal, i=1,size(N_sl))])
prm%tau_0 = pl%get_asFloat('tau_0') prm%tau_0 = pl%get_asFloat('tau_0')
prm%CEdgeDipMinDistance = pl%get_asFloat('D_a') prm%D_a = pl%get_asFloat('D_a')
prm%D0 = pl%get_asFloat('D_0') prm%D_0 = pl%get_asFloat('D_0')
prm%Qsd = pl%get_asFloat('Q_cl') prm%Q_cl = pl%get_asFloat('Q_cl')
prm%ExtendedDislocations = pl%get_asBool('extend_dislocations',defaultVal = .false.) prm%ExtendedDislocations = pl%get_asBool('extend_dislocations',defaultVal = .false.)
if (prm%ExtendedDislocations) then if (prm%ExtendedDislocations) then
prm%SFE_0K = pl%get_asFloat('Gamma_sf_0K') prm%Gamma_sf_0K = pl%get_asFloat('Gamma_sf_0K')
prm%dSFE_dT = pl%get_asFloat('dGamma_sf_dT') prm%dGamma_sf_dT = pl%get_asFloat('dGamma_sf_dT')
endif endif
prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation',defaultVal = .false.) prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation',defaultVal = .false.)
@ -238,29 +238,29 @@ module function plastic_dislotwin_init() result(myPlasticity)
! expand: family => system ! expand: family => system
rho_mob_0 = math_expand(rho_mob_0, N_sl) rho_mob_0 = math_expand(rho_mob_0, N_sl)
rho_dip_0 = math_expand(rho_dip_0, N_sl) rho_dip_0 = math_expand(rho_dip_0, N_sl)
prm%v0 = math_expand(prm%v0, N_sl) prm%v_0 = math_expand(prm%v_0, N_sl)
prm%b_sl = math_expand(prm%b_sl, N_sl) prm%b_sl = math_expand(prm%b_sl, N_sl)
prm%Delta_F = math_expand(prm%Delta_F, N_sl) prm%Q_s = math_expand(prm%Q_s, N_sl)
prm%CLambdaSlip = math_expand(prm%CLambdaSlip, N_sl) prm%i_sl = math_expand(prm%i_sl, N_sl)
prm%p = math_expand(prm%p, N_sl) prm%p = math_expand(prm%p, N_sl)
prm%q = math_expand(prm%q, N_sl) prm%q = math_expand(prm%q, N_sl)
prm%B = math_expand(prm%B, N_sl) prm%B = math_expand(prm%B, N_sl)
! sanity checks ! sanity checks
if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0'
if ( prm%Qsd <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl'
if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0' if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0'
if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0' if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0'
if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0' if (any(prm%v_0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0'
if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl'
if (any(prm%Delta_F <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s' if (any(prm%Q_s <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s'
if (any(prm%CLambdaSlip <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl' if (any(prm%i_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl'
if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B' if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B'
if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p_sl' if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p_sl'
if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q_sl' if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q_sl'
else slipActive else slipActive
rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray
allocate(prm%b_sl,prm%Delta_F,prm%v0,prm%CLambdaSlip,prm%p,prm%q,prm%B,source=emptyRealArray) allocate(prm%b_sl,prm%Q_s,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray)
allocate(prm%forestProjection(0,0),prm%h_sl_sl(0,0)) allocate(prm%forestProjection(0,0),prm%h_sl_sl(0,0))
endif slipActive endif slipActive
@ -279,7 +279,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%t_tw = pl%get_asFloats('t_tw', requiredSize=size(N_tw)) prm%t_tw = pl%get_asFloats('t_tw', requiredSize=size(N_tw))
prm%r = pl%get_asFloats('p_tw', requiredSize=size(N_tw)) prm%r = pl%get_asFloats('p_tw', requiredSize=size(N_tw))
prm%xc_twin = pl%get_asFloat('x_c_tw') prm%x_c_tw = pl%get_asFloat('x_c_tw')
prm%L_tw = pl%get_asFloat('L_tw') prm%L_tw = pl%get_asFloat('L_tw')
prm%i_tw = pl%get_asFloat('i_tw') prm%i_tw = pl%get_asFloat('i_tw')
@ -300,7 +300,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%r = math_expand(prm%r,N_tw) prm%r = math_expand(prm%r,N_tw)
! sanity checks ! sanity checks
if ( prm%xc_twin < 0.0_pReal) extmsg = trim(extmsg)//' x_c_twin' if ( prm%x_c_tw < 0.0_pReal) extmsg = trim(extmsg)//' x_c_twin'
if ( prm%L_tw < 0.0_pReal) extmsg = trim(extmsg)//' L_tw' if ( prm%L_tw < 0.0_pReal) extmsg = trim(extmsg)//' L_tw'
if ( prm%i_tw < 0.0_pReal) extmsg = trim(extmsg)//' i_tw' if ( prm%i_tw < 0.0_pReal) extmsg = trim(extmsg)//' i_tw'
if (any(prm%b_tw < 0.0_pReal)) extmsg = trim(extmsg)//' b_tw' if (any(prm%b_tw < 0.0_pReal)) extmsg = trim(extmsg)//' b_tw'
@ -324,8 +324,8 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%h = pl%get_asFloat('h', defaultVal=0.0_pReal) ! ToDo: How to handle that??? prm%h = pl%get_asFloat('h', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%i_tr = pl%get_asFloat('i_tr', defaultVal=0.0_pReal) ! ToDo: How to handle that??? prm%i_tr = pl%get_asFloat('i_tr', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%gamma_fcc_hex = pl%get_asFloat('delta_G') prm%delta_G = pl%get_asFloat('delta_G')
prm%xc_trans = pl%get_asFloat('x_c_tr', defaultVal=0.0_pReal) ! ToDo: How to handle that??? prm%x_c_tr = pl%get_asFloat('x_c_tr', defaultVal=0.0_pReal) ! ToDo: How to handle that???
prm%L_tr = pl%get_asFloat('L_tr') prm%L_tr = pl%get_asFloat('L_tr')
prm%h_tr_tr = lattice_interaction_TransByTrans(N_tr,pl%get_asFloats('h_tr_tr'), & prm%h_tr_tr = lattice_interaction_TransByTrans(N_tr,pl%get_asFloats('h_tr_tr'), &
@ -351,7 +351,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%s = math_expand(prm%s,N_tr) prm%s = math_expand(prm%s,N_tr)
! sanity checks ! sanity checks
if ( prm%xc_trans < 0.0_pReal) extmsg = trim(extmsg)//' x_c_trans' if ( prm%x_c_tr < 0.0_pReal) extmsg = trim(extmsg)//' x_c_trans'
if ( prm%L_tr < 0.0_pReal) extmsg = trim(extmsg)//' L_tr' if ( prm%L_tr < 0.0_pReal) extmsg = trim(extmsg)//' L_tr'
if ( prm%i_tr < 0.0_pReal) extmsg = trim(extmsg)//' i_tr' if ( prm%i_tr < 0.0_pReal) extmsg = trim(extmsg)//' i_tr'
if (any(prm%t_tr < 0.0_pReal)) extmsg = trim(extmsg)//' t_tr' if (any(prm%t_tr < 0.0_pReal)) extmsg = trim(extmsg)//' t_tr'
@ -366,15 +366,15 @@ module function plastic_dislotwin_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! shearband related parameters ! shearband related parameters
prm%sbVelocity = pl%get_asFloat('v_sb',defaultVal=0.0_pReal) prm%v_sb = pl%get_asFloat('v_sb',defaultVal=0.0_pReal)
if (prm%sbVelocity > 0.0_pReal) then if (prm%v_sb > 0.0_pReal) then
prm%sbResistance = pl%get_asFloat('xi_sb') prm%xi_sb = pl%get_asFloat('xi_sb')
prm%E_sb = pl%get_asFloat('Q_sb') prm%E_sb = pl%get_asFloat('Q_sb')
prm%p_sb = pl%get_asFloat('p_sb') prm%p_sb = pl%get_asFloat('p_sb')
prm%q_sb = pl%get_asFloat('q_sb') prm%q_sb = pl%get_asFloat('q_sb')
! sanity checks ! sanity checks
if (prm%sbResistance < 0.0_pReal) extmsg = trim(extmsg)//' xi_sb' if (prm%xi_sb < 0.0_pReal) extmsg = trim(extmsg)//' xi_sb'
if (prm%E_sb < 0.0_pReal) extmsg = trim(extmsg)//' Q_sb' if (prm%E_sb < 0.0_pReal) extmsg = trim(extmsg)//' Q_sb'
if (prm%p_sb <= 0.0_pReal) extmsg = trim(extmsg)//' p_sb' if (prm%p_sb <= 0.0_pReal) extmsg = trim(extmsg)//' p_sb'
if (prm%q_sb <= 0.0_pReal) extmsg = trim(extmsg)//' q_sb' if (prm%q_sb <= 0.0_pReal) extmsg = trim(extmsg)//' q_sb'
@ -386,8 +386,8 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%D = pl%get_asFloat('D') prm%D = pl%get_asFloat('D')
twinOrSlipActive: if (prm%sum_N_tw + prm%sum_N_tr > 0) then twinOrSlipActive: if (prm%sum_N_tw + prm%sum_N_tr > 0) then
prm%SFE_0K = pl%get_asFloat('Gamma_sf_0K') prm%Gamma_sf_0K = pl%get_asFloat('Gamma_sf_0K')
prm%dSFE_dT = pl%get_asFloat('dGamma_sf_dT') prm%dGamma_sf_dT = pl%get_asFloat('dGamma_sf_dT')
prm%V_cs = pl%get_asFloat('V_cs') prm%V_cs = pl%get_asFloat('V_cs')
endif twinOrSlipActive endif twinOrSlipActive
@ -602,7 +602,7 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
Lp = Lp * f_unrotated Lp = Lp * f_unrotated
dLp_dMp = dLp_dMp * f_unrotated dLp_dMp = dLp_dMp * f_unrotated
shearBandingContribution: if(dNeq0(prm%sbVelocity)) then shearBandingContribution: if(dNeq0(prm%v_sb)) then
BoltzmannRatio = prm%E_sb/(kB*T) BoltzmannRatio = prm%E_sb/(kB*T)
call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design?
@ -613,10 +613,10 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
tau = math_tensordot(Mp,P_sb) tau = math_tensordot(Mp,P_sb)
significantShearBandStress: if (abs(tau) > tol_math_check) then significantShearBandStress: if (abs(tau) > tol_math_check) then
StressRatio_p = (abs(tau)/prm%sbResistance)**prm%p_sb StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb
dot_gamma_sb = sign(prm%sbVelocity*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau) dot_gamma_sb = sign(prm%v_sb*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau)
ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%sbResistance & ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%xi_sb &
* (abs(tau)/prm%sbResistance)**(prm%p_sb-1.0_pReal) & * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) &
* (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal)
Lp = Lp + dot_gamma_sb * P_sb Lp = Lp + dot_gamma_sb * P_sb
@ -654,7 +654,7 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
Gamma, & !< stacking fault energy Gamma, & !< stacking fault energy
tau, & tau, &
sigma_cl, & !< climb stress sigma_cl, & !< climb stress
b_d !< ratio of burgers vector to stacking fault width b_d !< ratio of Burgers vector to stacking fault width
real(pReal), dimension(param(instance)%sum_N_sl) :: & real(pReal), dimension(param(instance)%sum_N_sl) :: &
dot_rho_dip_formation, & dot_rho_dip_formation, &
dot_rho_dip_climb, & dot_rho_dip_climb, &
@ -675,7 +675,7 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
call kinetics_slip(Mp,T,instance,of,dot_gamma_sl) call kinetics_slip(Mp,T,instance,of,dot_gamma_sl)
dot%gamma_sl(:,of) = abs(dot_gamma_sl) dot%gamma_sl(:,of) = abs(dot_gamma_sl)
rho_dip_distance_min = prm%CEdgeDipMinDistance*prm%b_sl rho_dip_distance_min = prm%D_a*prm%b_sl
slipState: do i = 1, prm%sum_N_sl slipState: do i = 1, prm%sum_N_sl
tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i))
@ -701,12 +701,12 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
!@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 !@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i))) sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i)))
if (prm%ExtendedDislocations) then if (prm%ExtendedDislocations) then
Gamma = prm%SFE_0K + prm%dSFE_dT * T Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T
b_d = 24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu)* Gamma/(prm%mu*prm%b_sl(i)) b_d = 24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu)* Gamma/(prm%mu*prm%b_sl(i))
else else
b_d = 1.0_pReal b_d = 1.0_pReal
endif endif
v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Qsd/(kB*T)) & v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) &
* (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal) * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal)
dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,of) & dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,of) &
@ -768,7 +768,7 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of)
sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,of)) sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,of))
sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,of)) sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,of))
Gamma = prm%SFE_0K + prm%dSFE_dT * T Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T
!* rescaled volume fraction for topology !* rescaled volume fraction for topology
f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,of)/prm%t_tw ! this is per system ... f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,of)/prm%t_tw ! this is per system ...
@ -776,7 +776,7 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of)
! ToDo ...Physically correct, but naming could be adjusted ! ToDo ...Physically correct, but naming could be adjusted
inv_lambda_sl_sl = sqrt(matmul(prm%forestProjection, & inv_lambda_sl_sl = sqrt(matmul(prm%forestProjection, &
stt%rho_mob(:,of)+stt%rho_dip(:,of)))/prm%CLambdaSlip stt%rho_mob(:,of)+stt%rho_dip(:,of)))/prm%i_sl
if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) & if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) &
inv_lambda_sl_tw = matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_twin) inv_lambda_sl_tw = matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_twin)
@ -805,21 +805,21 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of)
!* threshold stress for growing twin/martensite !* threshold stress for growing twin/martensite
if(prm%sum_N_tw == prm%sum_N_sl) & if(prm%sum_N_tw == prm%sum_N_sl) &
dst%tau_hat_tw(:,of) = Gamma/(3.0_pReal*prm%b_tw) & dst%tau_hat_tw(:,of) = Gamma/(3.0_pReal*prm%b_tw) &
+ 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip burgers here correct? + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct?
if(prm%sum_N_tr == prm%sum_N_sl) & if(prm%sum_N_tr == prm%sum_N_sl) &
dst%tau_hat_tr(:,of) = Gamma/(3.0_pReal*prm%b_tr) & dst%tau_hat_tr(:,of) = Gamma/(3.0_pReal*prm%b_tr) &
+ 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip burgers here correct? + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct?
+ prm%h*prm%gamma_fcc_hex/ (3.0_pReal*prm%b_tr) + prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr)
dst%V_tw(:,of) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,of)**2.0_pReal dst%V_tw(:,of) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,of)**2.0_pReal
dst%V_tr(:,of) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,of)**2.0_pReal dst%V_tr(:,of) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,of)**2.0_pReal
x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the burgers vector for slip and is the same for twin and trans x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans
dst%tau_r_tw(:,of) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%xc_twin)+cos(pi/3.0_pReal)/x0) dst%tau_r_tw(:,of) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the burgers vector for slip x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip
dst%tau_r_tr(:,of) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%xc_trans)+cos(pi/3.0_pReal)/x0) dst%tau_r_tr(:,of) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
end associate end associate
@ -927,8 +927,8 @@ pure subroutine kinetics_slip(Mp,T,instance,of, &
significantStress: where(tau_eff > tol_math_check) significantStress: where(tau_eff > tol_math_check)
stressRatio = tau_eff/prm%tau_0 stressRatio = tau_eff/prm%tau_0
StressRatio_p = stressRatio** prm%p StressRatio_p = stressRatio** prm%p
BoltzmannRatio = prm%Delta_F/(kB*T) BoltzmannRatio = prm%Q_s/(kB*T)
v_wait_inverse = prm%v0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q) v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q)
v_run_inverse = prm%B/(tau_eff*prm%b_sl) v_run_inverse = prm%B/(tau_eff*prm%b_sl)
dot_gamma_sl = sign(stt%rho_mob(:,of)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau) dot_gamma_sl = sign(stt%rho_mob(:,of)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau)

6
src/constitutive_plastic_isotropic.f90
View File

@ -14,7 +14,7 @@ submodule(constitutive:constitutive_plastic) plastic_isotropic
M, & !< Taylor factor M, & !< Taylor factor
dot_gamma_0, & !< reference strain rate dot_gamma_0, & !< reference strain rate
n, & !< stress exponent n, & !< stress exponent
h0, & h_0, &
h_ln, & h_ln, &
xi_inf, & !< maximum critical stress xi_inf, & !< maximum critical stress
a, & a, &
@ -109,7 +109,7 @@ module function plastic_isotropic_init() result(myPlasticity)
prm%xi_inf = pl%get_asFloat('xi_inf') prm%xi_inf = pl%get_asFloat('xi_inf')
prm%dot_gamma_0 = pl%get_asFloat('dot_gamma_0') prm%dot_gamma_0 = pl%get_asFloat('dot_gamma_0')
prm%n = pl%get_asFloat('n') prm%n = pl%get_asFloat('n')
prm%h0 = pl%get_asFloat('h_0') prm%h_0 = pl%get_asFloat('h_0')
prm%M = pl%get_asFloat('M') prm%M = pl%get_asFloat('M')
prm%h_ln = pl%get_asFloat('h_ln', defaultVal=0.0_pReal) prm%h_ln = pl%get_asFloat('h_ln', defaultVal=0.0_pReal)
prm%c_1 = pl%get_asFloat('c_1', defaultVal=0.0_pReal) prm%c_1 = pl%get_asFloat('c_1', defaultVal=0.0_pReal)
@ -310,7 +310,7 @@ module subroutine plastic_isotropic_dotState(Mp,instance,of)
/ prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n) / prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n)
endif endif
dot%xi(of) = dot_gamma & dot%xi(of) = dot_gamma &
* ( prm%h0 + prm%h_ln * log(dot_gamma) ) & * ( prm%h_0 + prm%h_ln * log(dot_gamma) ) &
* abs( 1.0_pReal - stt%xi(of)/xi_inf_star )**prm%a & * abs( 1.0_pReal - stt%xi(of)/xi_inf_star )**prm%a &
* sign(1.0_pReal, 1.0_pReal - stt%xi(of)/xi_inf_star) * sign(1.0_pReal, 1.0_pReal - stt%xi(of)/xi_inf_star)
else else

70
src/constitutive_plastic_kinehardening.f90
View File

@ -9,15 +9,15 @@ submodule(constitutive:constitutive_plastic) plastic_kinehardening
type :: tParameters type :: tParameters
real(pReal) :: & real(pReal) :: &
gdot0 = 1.0_pReal, & !< reference shear strain rate for slip n = 1.0_pReal, & !< stress exponent for slip
n = 1.0_pReal !< stress exponent for slip dot_gamma_0 = 1.0_pReal !< reference shear strain rate for slip
real(pReal), allocatable, dimension(:) :: & real(pReal), allocatable, dimension(:) :: &
theta0, & !< initial hardening rate of forward stress for each slip h_0_f, & !< initial hardening rate of forward stress for each slip
theta1, & !< asymptotic hardening rate of forward stress for each slip h_inf_f, & !< asymptotic hardening rate of forward stress for each slip
theta0_b, & !< initial hardening rate of back stress for each slip h_0_b, & !< initial hardening rate of back stress for each slip
theta1_b, & !< asymptotic hardening rate of back stress for each slip h_inf_b, & !< asymptotic hardening rate of back stress for each slip
tau1, & xi_inf_f, &
tau1_b xi_inf_b
real(pReal), allocatable, dimension(:,:) :: & real(pReal), allocatable, dimension(:,:) :: &
interaction_slipslip !< slip resistance from slip activity interaction_slipslip !< slip resistance from slip activity
real(pReal), allocatable, dimension(:,:,:) :: & real(pReal), allocatable, dimension(:,:,:) :: &
@ -125,7 +125,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
phase%get_asFloat('c/a',defaultVal=0.0_pReal)) phase%get_asFloat('c/a',defaultVal=0.0_pReal))
if(trim(phase%get_asString('lattice')) == 'bcc') then if(trim(phase%get_asString('lattice')) == 'bcc') then
a = pl%get_asFloats('nonSchmid_coefficients',defaultVal = emptyRealArray) a = pl%get_asFloats('a_nonSchmid',defaultVal = emptyRealArray)
if(size(a) > 0) prm%nonSchmidActive = .true. if(size(a) > 0) prm%nonSchmidActive = .true.
prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1) prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1) prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
@ -138,37 +138,37 @@ module function plastic_kinehardening_init() result(myPlasticity)
phase%get_asString('lattice')) phase%get_asString('lattice'))
xi_0 = pl%get_asFloats('xi_0', requiredSize=size(N_sl)) xi_0 = pl%get_asFloats('xi_0', requiredSize=size(N_sl))
prm%tau1 = pl%get_asFloats('xi_inf_f', requiredSize=size(N_sl)) prm%xi_inf_f = pl%get_asFloats('xi_inf_f', requiredSize=size(N_sl))
prm%tau1_b = pl%get_asFloats('xi_inf_b', requiredSize=size(N_sl)) prm%xi_inf_b = pl%get_asFloats('xi_inf_b', requiredSize=size(N_sl))
prm%theta0 = pl%get_asFloats('h_0_f', requiredSize=size(N_sl)) prm%h_0_f = pl%get_asFloats('h_0_f', requiredSize=size(N_sl))
prm%theta1 = pl%get_asFloats('h_inf_f', requiredSize=size(N_sl)) prm%h_inf_f = pl%get_asFloats('h_inf_f', requiredSize=size(N_sl))
prm%theta0_b = pl%get_asFloats('h_0_b', requiredSize=size(N_sl)) prm%h_0_b = pl%get_asFloats('h_0_b', requiredSize=size(N_sl))
prm%theta1_b = pl%get_asFloats('h_inf_b', requiredSize=size(N_sl)) prm%h_inf_b = pl%get_asFloats('h_inf_b', requiredSize=size(N_sl))
prm%gdot0 = pl%get_asFloat('dot_gamma_0') prm%dot_gamma_0 = pl%get_asFloat('dot_gamma_0')
prm%n = pl%get_asFloat('n') prm%n = pl%get_asFloat('n')
! expand: family => system ! expand: family => system
xi_0 = math_expand(xi_0, N_sl) xi_0 = math_expand(xi_0, N_sl)
prm%tau1 = math_expand(prm%tau1, N_sl) prm%xi_inf_f = math_expand(prm%xi_inf_f, N_sl)
prm%tau1_b = math_expand(prm%tau1_b, N_sl) prm%xi_inf_b = math_expand(prm%xi_inf_b, N_sl)
prm%theta0 = math_expand(prm%theta0, N_sl) prm%h_0_f = math_expand(prm%h_0_f, N_sl)
prm%theta1 = math_expand(prm%theta1, N_sl) prm%h_inf_f = math_expand(prm%h_inf_f, N_sl)
prm%theta0_b = math_expand(prm%theta0_b,N_sl) prm%h_0_b = math_expand(prm%h_0_b, N_sl)
prm%theta1_b = math_expand(prm%theta1_b,N_sl) prm%h_inf_b = math_expand(prm%h_inf_b, N_sl)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! sanity checks ! sanity checks
if ( prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0' if ( prm%dot_gamma_0 <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0'
if ( prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' n' if ( prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' n'
if (any(xi_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_0' if (any(xi_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_0'
if (any(prm%tau1 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_f' if (any(prm%xi_inf_f <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_f'
if (any(prm%tau1_b <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_b' if (any(prm%xi_inf_b <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_b'
!ToDo: Any sensible checks for theta? !ToDo: Any sensible checks for theta?
else slipActive else slipActive
xi_0 = emptyRealArray xi_0 = emptyRealArray
allocate(prm%tau1,prm%tau1_b,prm%theta0,prm%theta1,prm%theta0_b,prm%theta1_b,source=emptyRealArray) allocate(prm%xi_inf_f,prm%xi_inf_b,prm%h_0_f,prm%h_inf_f,prm%h_0_b,prm%h_inf_b,source=emptyRealArray)
allocate(prm%interaction_SlipSlip(0,0)) allocate(prm%interaction_SlipSlip(0,0))
endif slipActive endif slipActive
@ -303,16 +303,16 @@ module subroutine plastic_kinehardening_dotState(Mp,instance,of)
dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) & dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) &
* ( prm%theta1 & * ( prm%h_inf_f &
+ (prm%theta0 - prm%theta1 + prm%theta0*prm%theta1*sumGamma/prm%tau1) & + (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) &
* exp(-sumGamma*prm%theta0/prm%tau1) & * exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) &
) )
dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * & dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * &
( prm%theta1_b + & ( prm%h_inf_b + &
(prm%theta0_b - prm%theta1_b & (prm%h_0_b - prm%h_inf_b &
+ prm%theta0_b*prm%theta1_b/(prm%tau1_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))& + prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))&
) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%theta0_b/(prm%tau1_b+stt%chi0(:,of))) & ) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%h_0_b/(prm%xi_inf_b+stt%chi0(:,of))) &
) )
end associate end associate
@ -442,14 +442,14 @@ pure subroutine kinetics(Mp,instance,of, &
enddo enddo
where(dNeq0(tau_pos)) where(dNeq0(tau_pos))
gdot_pos = prm%gdot0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active gdot_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos) * sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos)
else where else where
gdot_pos = 0.0_pReal gdot_pos = 0.0_pReal
end where end where
where(dNeq0(tau_neg)) where(dNeq0(tau_neg))
gdot_neg = prm%gdot0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 gdot_neg = prm%dot_gamma_0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg) * sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg)
else where else where
gdot_neg = 0.0_pReal gdot_neg = 0.0_pReal

330
src/constitutive_plastic_nonlocal.f90
View File

@ -46,58 +46,58 @@ submodule(constitutive:constitutive_plastic) plastic_nonlocal
type :: tInitialParameters !< container type for internal constitutive parameters type :: tInitialParameters !< container type for internal constitutive parameters
real(pReal) :: & real(pReal) :: &
rhoSglScatter, & !< standard deviation of scatter in initial dislocation density sigma_rho_u, & !< standard deviation of scatter in initial dislocation density
rhoSglRandom, & random_rho_u, &
rhoSglRandomBinning random_rho_u_binning
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
rhoSglEdgePos0, & !< initial edge_pos dislocation density rho_u_ed_pos_0, & !< initial edge_pos dislocation density
rhoSglEdgeNeg0, & !< initial edge_neg dislocation density rho_u_ed_neg_0, & !< initial edge_neg dislocation density
rhoSglScrewPos0, & !< initial screw_pos dislocation density rho_u_sc_pos_0, & !< initial screw_pos dislocation density
rhoSglScrewNeg0, & !< initial screw_neg dislocation density rho_u_sc_neg_0, & !< initial screw_neg dislocation density
rhoDipEdge0, & !< initial edge dipole dislocation density rho_d_ed_0, & !< initial edge dipole dislocation density
rhoDipScrew0 !< initial screw dipole dislocation density rho_d_sc_0 !< initial screw dipole dislocation density
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
N_sl N_sl
end type tInitialParameters end type tInitialParameters
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: & real(pReal) :: &
atomicVolume, & !< atomic volume V_at, & !< atomic volume
Dsd0, & !< prefactor for self-diffusion coefficient D_0, & !< prefactor for self-diffusion coefficient
selfDiffusionEnergy, & !< activation enthalpy for diffusion Q_cl, & !< activation enthalpy for diffusion
atol_rho, & !< absolute tolerance for dislocation density in state integration atol_rho, & !< absolute tolerance for dislocation density in state integration
significantRho, & !< density considered significant rho_significant, & !< density considered significant
significantN, & !< number of dislocations considered significant rho_min, & !< number of dislocations considered significant
doublekinkwidth, & !< width of a doubkle kink in multiples of the burgers vector length b w, & !< width of a doubkle kink in multiples of the Burgers vector length b
solidSolutionEnergy, & !< activation energy for solid solution in J Q_sol, & !< activation energy for solid solution in J
solidSolutionSize, & !< solid solution obstacle size in multiples of the burgers vector length f_sol, & !< solid solution obstacle size in multiples of the Burgers vector length
solidSolutionConcentration, & !< 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)
viscosity, & !< viscosity for dislocation glide in Pa s eta, & !< viscosity for dislocation glide in Pa s
fattack, & !< attack frequency in Hz nu_a, & !< attack frequency in Hz
surfaceTransmissivity, & !< transmissivity at free surface chi_surface, & !< transmissivity at free surface
grainboundaryTransmissivity, & !< transmissivity at grain boundary (identified by different texture) chi_GB, & !< transmissivity at grain boundary (identified by different texture)
CFLfactor, & !< safety factor for CFL flux condition f_c, & !< safety factor for CFL flux condition
fEdgeMultiplication, & !< 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)
linetensionEffect, & f_F, &
edgeJogFactor, & f_ed, &
mu, & mu, &
nu nu
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
minDipoleHeight_edge, & !< minimum stable edge dipole height d_ed, & !< minimum stable edge dipole height
minDipoleHeight_screw, & !< minimum stable screw dipole height d_sc, & !< minimum stable screw dipole height
peierlsstress_edge, & tau_Peierls_ed, &
peierlsstress_screw, & tau_Peierls_sc, &
lambda0, & !< mean free path prefactor for each i_sl, & !< mean free path prefactor for each
burgers !< absolute length of burgers vector [m] b_sl !< absolute length of Burgers vector [m]
real(pReal), dimension(:,:), allocatable :: & real(pReal), dimension(:,:), allocatable :: &
slip_normal, & slip_normal, &
slip_direction, & slip_direction, &
slip_transverse, & slip_transverse, &
minDipoleHeight, & ! edge and screw minDipoleHeight, & ! edge and screw
peierlsstress, & ! edge and screw peierlsstress, & ! edge and screw
interactionSlipSlip ,& !< coefficients for slip-slip interaction h_sl_sl ,& !< coefficients for slip-slip interaction
forestProjection_Edge, & !< matrix of forest projections of edge dislocations forestProjection_Edge, & !< matrix of forest projections of edge dislocations
forestProjection_Screw !< matrix of forest projections of screw dislocations forestProjection_Screw !< matrix of forest projections of screw dislocations
real(pReal), dimension(:,:,:), allocatable :: & real(pReal), dimension(:,:,:), allocatable :: &
@ -244,7 +244,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
phase%get_asFloat('c/a',defaultVal=0.0_pReal)) phase%get_asFloat('c/a',defaultVal=0.0_pReal))
if(trim(phase%get_asString('lattice')) == 'bcc') then if(trim(phase%get_asString('lattice')) == 'bcc') then
a = pl%get_asFloats('nonSchmid_coefficients',defaultVal = emptyRealArray) a = pl%get_asFloats('a_nonSchmid',defaultVal = emptyRealArray)
if(size(a) > 0) prm%nonSchmidActive = .true. if(size(a) > 0) prm%nonSchmidActive = .true.
prm%nonSchmid_pos = lattice_nonSchmidMatrix(ini%N_sl,a,+1) prm%nonSchmid_pos = lattice_nonSchmidMatrix(ini%N_sl,a,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(ini%N_sl,a,-1) prm%nonSchmid_neg = lattice_nonSchmidMatrix(ini%N_sl,a,-1)
@ -253,7 +253,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
prm%nonSchmid_neg = prm%Schmid prm%nonSchmid_neg = prm%Schmid
endif endif
prm%interactionSlipSlip = lattice_interaction_SlipBySlip(ini%N_sl, & prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl, &
pl%get_asFloats('h_sl_sl'), & pl%get_asFloats('h_sl_sl'), &
phase%get_asString('lattice')) phase%get_asString('lattice'))
@ -279,112 +279,112 @@ module function plastic_nonlocal_init() result(myPlasticity)
enddo enddo
enddo enddo
ini%rhoSglEdgePos0 = pl%get_asFloats('rho_u_ed_pos_0', requiredSize=size(ini%N_sl)) ini%rho_u_ed_pos_0 = pl%get_asFloats('rho_u_ed_pos_0', requiredSize=size(ini%N_sl))
ini%rhoSglEdgeNeg0 = pl%get_asFloats('rho_u_ed_neg_0', requiredSize=size(ini%N_sl)) ini%rho_u_ed_neg_0 = pl%get_asFloats('rho_u_ed_neg_0', requiredSize=size(ini%N_sl))
ini%rhoSglScrewPos0 = pl%get_asFloats('rho_u_sc_pos_0', requiredSize=size(ini%N_sl)) ini%rho_u_sc_pos_0 = pl%get_asFloats('rho_u_sc_pos_0', requiredSize=size(ini%N_sl))
ini%rhoSglScrewNeg0 = pl%get_asFloats('rho_u_sc_neg_0', requiredSize=size(ini%N_sl)) ini%rho_u_sc_neg_0 = pl%get_asFloats('rho_u_sc_neg_0', requiredSize=size(ini%N_sl))
ini%rhoDipEdge0 = pl%get_asFloats('rho_d_ed_0', requiredSize=size(ini%N_sl)) ini%rho_d_ed_0 = pl%get_asFloats('rho_d_ed_0', requiredSize=size(ini%N_sl))
ini%rhoDipScrew0 = pl%get_asFloats('rho_d_sc_0', requiredSize=size(ini%N_sl)) ini%rho_d_sc_0 = pl%get_asFloats('rho_d_sc_0', requiredSize=size(ini%N_sl))
prm%lambda0 = pl%get_asFloats('i_sl', requiredSize=size(ini%N_sl)) prm%i_sl = pl%get_asFloats('i_sl', requiredSize=size(ini%N_sl))
prm%burgers = pl%get_asFloats('b_sl', requiredSize=size(ini%N_sl)) prm%b_sl = pl%get_asFloats('b_sl', requiredSize=size(ini%N_sl))
prm%lambda0 = math_expand(prm%lambda0,ini%N_sl) prm%i_sl = math_expand(prm%i_sl,ini%N_sl)
prm%burgers = math_expand(prm%burgers,ini%N_sl) prm%b_sl = math_expand(prm%b_sl,ini%N_sl)
prm%minDipoleHeight_edge = pl%get_asFloats('d_ed', requiredSize=size(ini%N_sl)) prm%d_ed = pl%get_asFloats('d_ed', requiredSize=size(ini%N_sl))
prm%minDipoleHeight_screw = pl%get_asFloats('d_sc', requiredSize=size(ini%N_sl)) prm%d_sc = pl%get_asFloats('d_sc', requiredSize=size(ini%N_sl))
prm%minDipoleHeight_edge = math_expand(prm%minDipoleHeight_edge, ini%N_sl) prm%d_ed = math_expand(prm%d_ed,ini%N_sl)
prm%minDipoleHeight_screw = math_expand(prm%minDipoleHeight_screw,ini%N_sl) prm%d_sc = math_expand(prm%d_sc,ini%N_sl)
allocate(prm%minDipoleHeight(prm%sum_N_sl,2)) allocate(prm%minDipoleHeight(prm%sum_N_sl,2))
prm%minDipoleHeight(:,1) = prm%minDipoleHeight_edge prm%minDipoleHeight(:,1) = prm%d_ed
prm%minDipoleHeight(:,2) = prm%minDipoleHeight_screw prm%minDipoleHeight(:,2) = prm%d_sc
prm%peierlsstress_edge = pl%get_asFloats('tau_peierls_ed', requiredSize=size(ini%N_sl)) prm%tau_Peierls_ed = pl%get_asFloats('tau_Peierls_ed', requiredSize=size(ini%N_sl))
prm%peierlsstress_screw = pl%get_asFloats('tau_peierls_sc', requiredSize=size(ini%N_sl)) prm%tau_Peierls_sc = pl%get_asFloats('tau_Peierls_sc', requiredSize=size(ini%N_sl))
prm%peierlsstress_edge = math_expand(prm%peierlsstress_edge, ini%N_sl) prm%tau_Peierls_ed = math_expand(prm%tau_Peierls_ed,ini%N_sl)
prm%peierlsstress_screw = math_expand(prm%peierlsstress_screw,ini%N_sl) prm%tau_Peierls_sc = math_expand(prm%tau_Peierls_sc,ini%N_sl)
allocate(prm%peierlsstress(prm%sum_N_sl,2)) allocate(prm%peierlsstress(prm%sum_N_sl,2))
prm%peierlsstress(:,1) = prm%peierlsstress_edge prm%peierlsstress(:,1) = prm%tau_Peierls_ed
prm%peierlsstress(:,2) = prm%peierlsstress_screw prm%peierlsstress(:,2) = prm%tau_Peierls_sc
prm%significantRho = pl%get_asFloat('rho_significant') prm%rho_significant = pl%get_asFloat('rho_significant')
prm%significantN = pl%get_asFloat('rho_num_significant', 0.0_pReal) prm%rho_min = pl%get_asFloat('rho_min', 0.0_pReal)
prm%CFLfactor = pl%get_asFloat('f_c',defaultVal=2.0_pReal) prm%f_c = pl%get_asFloat('f_c',defaultVal=2.0_pReal)
prm%atomicVolume = pl%get_asFloat('V_at') prm%V_at = pl%get_asFloat('V_at')
prm%Dsd0 = pl%get_asFloat('D_0') !,'dsd0' prm%D_0 = pl%get_asFloat('D_0')
prm%selfDiffusionEnergy = pl%get_asFloat('Q_cl') !,'qsd' prm%Q_cl = pl%get_asFloat('Q_cl')
prm%linetensionEffect = pl%get_asFloat('f_F') prm%f_F = pl%get_asFloat('f_F')
prm%edgeJogFactor = pl%get_asFloat('f_ed') !,'edgejogs' prm%f_ed = pl%get_asFloat('f_ed') !,'edgejogs'
prm%doublekinkwidth = pl%get_asFloat('w') prm%w = pl%get_asFloat('w')
prm%solidSolutionEnergy = pl%get_asFloat('Q_sol') prm%Q_sol = pl%get_asFloat('Q_sol')
prm%solidSolutionSize = pl%get_asFloat('f_sol') prm%f_sol = pl%get_asFloat('f_sol')
prm%solidSolutionConcentration = pl%get_asFloat('c_sol') prm%c_sol = pl%get_asFloat('c_sol')
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%viscosity = pl%get_asFloat('eta') prm%eta = pl%get_asFloat('eta')
prm%fattack = pl%get_asFloat('nu_a') prm%nu_a = pl%get_asFloat('nu_a')
! ToDo: discuss logic ! ToDo: discuss logic
ini%rhoSglScatter = pl%get_asFloat('sigma_rho_u') ini%sigma_rho_u = pl%get_asFloat('sigma_rho_u')
ini%rhoSglRandom = pl%get_asFloat('random_rho_u',defaultVal= 0.0_pReal) ini%random_rho_u = pl%get_asFloat('random_rho_u',defaultVal= 0.0_pReal)
if (pl%contains('random_rho_u')) & if (pl%contains('random_rho_u')) &
ini%rhoSglRandomBinning = pl%get_asFloat('random_rho_u_binning',defaultVal=0.0_pReal) !ToDo: useful default? ini%random_rho_u_binning = pl%get_asFloat('random_rho_u_binning',defaultVal=0.0_pReal) !ToDo: useful default?
! if (rhoSglRandom(instance) < 0.0_pReal) & ! if (rhoSglRandom(instance) < 0.0_pReal) &
! if (rhoSglRandomBinning(instance) <= 0.0_pReal) & ! if (rhoSglRandomBinning(instance) <= 0.0_pReal) &
prm%surfaceTransmissivity = pl%get_asFloat('chi_surface',defaultVal=1.0_pReal) prm%chi_surface = pl%get_asFloat('chi_surface',defaultVal=1.0_pReal)
prm%grainboundaryTransmissivity = pl%get_asFloat('chi_GB', defaultVal=-1.0_pReal) prm%chi_GB = pl%get_asFloat('chi_GB', defaultVal=-1.0_pReal)
prm%fEdgeMultiplication = pl%get_asFloat('f_ed_mult') prm%f_ed_mult = pl%get_asFloat('f_ed_mult')
prm%shortRangeStressCorrection = pl%get_asBool('short_range_stress_correction', defaultVal = .false.) prm%shortRangeStressCorrection = pl%get_asBool('short_range_stress_correction', defaultVal = .false.)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! sanity checks ! sanity checks
if (any(prm%burgers < 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' if (any(prm%b_sl < 0.0_pReal)) extmsg = trim(extmsg)//' b_sl'
if (any(prm%lambda0 <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl' if (any(prm%i_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl'
if (any(ini%rhoSglEdgePos0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_ed_pos_0' if (any(ini%rho_u_ed_pos_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_ed_pos_0'
if (any(ini%rhoSglEdgeNeg0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_ed_neg_0' if (any(ini%rho_u_ed_neg_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_ed_neg_0'
if (any(ini%rhoSglScrewPos0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_sc_pos_0' if (any(ini%rho_u_sc_pos_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_sc_pos_0'
if (any(ini%rhoSglScrewNeg0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_sc_neg_0' if (any(ini%rho_u_sc_neg_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_u_sc_neg_0'
if (any(ini%rhoDipEdge0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_d_ed_0' if (any(ini%rho_d_ed_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_d_ed_0'
if (any(ini%rhoDipScrew0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_d_sc_0' if (any(ini%rho_d_sc_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_d_sc_0'
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%viscosity <= 0.0_pReal) extmsg = trim(extmsg)//' eta' if (prm%eta <= 0.0_pReal) extmsg = trim(extmsg)//' eta'
if (prm%selfDiffusionEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' if (prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl'
if (prm%fattack <= 0.0_pReal) extmsg = trim(extmsg)//' nu_a' if (prm%nu_a <= 0.0_pReal) extmsg = trim(extmsg)//' nu_a'
if (prm%doublekinkwidth <= 0.0_pReal) extmsg = trim(extmsg)//' w' if (prm%w <= 0.0_pReal) extmsg = trim(extmsg)//' w'
if (prm%Dsd0 < 0.0_pReal) extmsg = trim(extmsg)//' D_0' if (prm%D_0 < 0.0_pReal) extmsg = trim(extmsg)//' D_0'
if (prm%atomicVolume <= 0.0_pReal) extmsg = trim(extmsg)//' V_at' ! ToDo: in disloTungsten, the atomic volume is given as a factor if (prm%V_at <= 0.0_pReal) extmsg = trim(extmsg)//' V_at' ! ToDo: in disloTungsten, the atomic volume is given as a factor
if (prm%significantN < 0.0_pReal) extmsg = trim(extmsg)//' rho_num_significant' if (prm%rho_min < 0.0_pReal) extmsg = trim(extmsg)//' rho_min'
if (prm%significantrho < 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%CFLfactor < 0.0_pReal) extmsg = trim(extmsg)//' f_c' if (prm%f_c < 0.0_pReal) extmsg = trim(extmsg)//' f_c'
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'
if (prm%linetensionEffect < 0.0_pReal .or. prm%linetensionEffect > 1.0_pReal) & if (prm%f_F < 0.0_pReal .or. prm%f_F > 1.0_pReal) &
extmsg = trim(extmsg)//' f_F' extmsg = trim(extmsg)//' f_F'
if (prm%edgeJogFactor < 0.0_pReal .or. prm%edgeJogFactor > 1.0_pReal) & if (prm%f_ed < 0.0_pReal .or. prm%f_ed > 1.0_pReal) &
extmsg = trim(extmsg)//' f_ed' extmsg = trim(extmsg)//' f_ed'
if (prm%solidSolutionEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' Q_sol' if (prm%Q_sol <= 0.0_pReal) extmsg = trim(extmsg)//' Q_sol'
if (prm%solidSolutionSize <= 0.0_pReal) extmsg = trim(extmsg)//' f_sol' if (prm%f_sol <= 0.0_pReal) extmsg = trim(extmsg)//' f_sol'
if (prm%solidSolutionConcentration <= 0.0_pReal) extmsg = trim(extmsg)//' c_sol' if (prm%c_sol <= 0.0_pReal) extmsg = trim(extmsg)//' c_sol'
if (prm%grainboundaryTransmissivity > 1.0_pReal) extmsg = trim(extmsg)//' chi_GB' if (prm%chi_GB > 1.0_pReal) extmsg = trim(extmsg)//' chi_GB'
if (prm%surfaceTransmissivity < 0.0_pReal .or. prm%surfaceTransmissivity > 1.0_pReal) & if (prm%chi_surface < 0.0_pReal .or. prm%chi_surface > 1.0_pReal) &
extmsg = trim(extmsg)//' chi_surface' extmsg = trim(extmsg)//' chi_surface'
if (prm%fEdgeMultiplication < 0.0_pReal .or. prm%fEdgeMultiplication > 1.0_pReal) & if (prm%f_ed_mult < 0.0_pReal .or. prm%f_ed_mult > 1.0_pReal) &
extmsg = trim(extmsg)//' f_ed_mult' extmsg = trim(extmsg)//' f_ed_mult'
endif slipActive endif slipActive
@ -620,16 +620,16 @@ module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el)
! coefficients are corrected for the line tension effect ! coefficients are corrected for the line tension effect
! (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals) ! (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals)
if (any(lattice_structure(material_phaseAt(1,el)) == [LATTICE_bcc_ID,LATTICE_fcc_ID])) then if (any(lattice_structure(material_phaseAt(1,el)) == [LATTICE_bcc_ID,LATTICE_fcc_ID])) then
myInteractionMatrix = prm%interactionSlipSlip & myInteractionMatrix = prm%h_sl_sl &
* spread(( 1.0_pReal - prm%linetensionEffect & * spread(( 1.0_pReal - prm%f_F &
+ prm%linetensionEffect & + prm%f_F &
* log(0.35_pReal * prm%burgers * sqrt(max(stt%rho_forest(:,of),prm%significantRho))) & * log(0.35_pReal * prm%b_sl * sqrt(max(stt%rho_forest(:,of),prm%rho_significant))) &
/ log(0.35_pReal * prm%burgers * 1e6_pReal))** 2.0_pReal,2,prm%sum_N_sl) / log(0.35_pReal * prm%b_sl * 1e6_pReal))** 2.0_pReal,2,prm%sum_N_sl)
else else
myInteractionMatrix = prm%interactionSlipSlip myInteractionMatrix = prm%h_sl_sl
endif endif
dst%tau_pass(:,of) = prm%mu * prm%burgers & dst%tau_pass(:,of) = prm%mu * prm%b_sl &
* sqrt(matmul(myInteractionMatrix,sum(abs(rho),2))) * sqrt(matmul(myInteractionMatrix,sum(abs(rho),2)))
!*** calculate the dislocation stress of the neighboring excess dislocation densities !*** calculate the dislocation stress of the neighboring excess dislocation densities
@ -731,7 +731,7 @@ module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el)
where(rhoTotal > 0.0_pReal) rhoExcessGradient_over_rho = rhoExcessGradient / rhoTotal where(rhoTotal > 0.0_pReal) rhoExcessGradient_over_rho = rhoExcessGradient / rhoTotal
! ... gives the local stress correction when multiplied with a factor ! ... gives the local stress correction when multiplied with a factor
dst%tau_back(s,of) = - prm%mu * prm%burgers(s) / (2.0_pReal * PI) & dst%tau_back(s,of) = - prm%mu * prm%b_sl(s) / (2.0_pReal * PI) &
* ( rhoExcessGradient_over_rho(1) / (1.0_pReal - prm%nu) & * ( rhoExcessGradient_over_rho(1) / (1.0_pReal - prm%nu) &
+ rhoExcessGradient_over_rho(2)) + rhoExcessGradient_over_rho(2))
enddo enddo
@ -841,7 +841,7 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) & forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) &
rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t))
gdotTotal = sum(rhoSgl(:,1:4) * v, 2) * prm%burgers gdotTotal = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
@ -850,10 +850,10 @@ module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
forall (i=1:3,j=1:3,k=1:3,l=1:3) & forall (i=1:3,j=1:3,k=1:3,l=1:3) &
dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) & dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) &
+ prm%Schmid(i,j,s) * prm%Schmid(k,l,s) & + prm%Schmid(i,j,s) * prm%Schmid(k,l,s) &
* sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%burgers(s) & * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) &
+ prm%Schmid(i,j,s) & + prm%Schmid(i,j,s) &
* ( prm%nonSchmid_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & * ( prm%nonSchmid_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) &
- prm%nonSchmid_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%burgers(s) - prm%nonSchmid_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s)
enddo enddo
end associate end associate
@ -929,8 +929,8 @@ module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el)
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
enddo enddo
dUpper(:,1) = prm%mu * prm%burgers/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau))
dUpper(:,2) = prm%mu * prm%burgers/(4.0_pReal * PI * abs(tau)) dUpper(:,2) = prm%mu * prm%b_sl/(4.0_pReal * PI * abs(tau))
where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) & where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) &
dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1)) dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1))
@ -1041,7 +1041,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
my_rhoSgl0 = rho0(:,sgl) my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of)
gdot = rhoSgl(:,1:4) * v * spread(prm%burgers,2,4) gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
#ifdef DEBUG #ifdef DEBUG
if (debugConstitutive%basic & if (debugConstitutive%basic &
@ -1060,8 +1060,8 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
enddo enddo
dLower = prm%minDipoleHeight dLower = prm%minDipoleHeight
dUpper(:,1) = prm%mu * prm%burgers/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau))
dUpper(:,2) = prm%mu * prm%burgers/(4.0_pReal * PI * abs(tau)) dUpper(:,2) = prm%mu * prm%b_sl/(4.0_pReal * PI * abs(tau))
where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) & where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) &
dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1)) dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1))
@ -1075,18 +1075,18 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
rhoDotMultiplication = 0.0_pReal rhoDotMultiplication = 0.0_pReal
isBCC: if (lattice_structure(ph) == LATTICE_bcc_ID) then isBCC: if (lattice_structure(ph) == LATTICE_bcc_ID) 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)
rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%burgers(s) & ! assuming double-cross-slip of screws to be decisive for multiplication rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication
* sqrt(stt%rho_forest(s,of)) / prm%lambda0(s) ! & ! mean free path * sqrt(stt%rho_forest(s,of)) / prm%i_sl(s) ! & ! mean free path
! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation ! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation
rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%burgers(s) & ! assuming double-cross-slip of screws to be decisive for multiplication rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication
* sqrt(stt%rho_forest(s,of)) / prm%lambda0(s) ! & ! mean free path * sqrt(stt%rho_forest(s,of)) / prm%i_sl(s) ! & ! mean free path
! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation ! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation
endforall endforall
else isBCC else isBCC
rhoDotMultiplication(:,1:4) = spread( & rhoDotMultiplication(:,1:4) = spread( &
(sum(abs(gdot(:,1:2)),2) * prm%fEdgeMultiplication + 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(:,of)) / prm%lambda0 / prm%burgers, 2, 4) * sqrt(stt%rho_forest(:,of)) / prm%i_sl / prm%b_sl, 2, 4)
endif isBCC endif isBCC
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),of) forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),of)
@ -1097,20 +1097,20 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
!*** formation by glide !*** formation by glide
do c = 1,2 do c = 1,2
rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%burgers & 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
+ rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile + rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile
+ abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) ! positive mobile --> negative immobile + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) ! positive mobile --> negative immobile
rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pReal * dUpper(:,c) / prm%burgers & rhoDotSingle2DipoleGlide(:,2*c) = -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
+ rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile + rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile
+ abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c))) ! negative mobile --> positive immobile + abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c))) ! negative mobile --> positive immobile
rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%burgers & rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
* rhoSgl(:,2*c+3) * abs(gdot(:,2*c)) ! negative mobile --> positive immobile * rhoSgl(:,2*c+3) * abs(gdot(:,2*c)) ! negative mobile --> positive immobile
rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%burgers & rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
* rhoSgl(:,2*c+4) * abs(gdot(:,2*c-1)) ! positive mobile --> negative immobile * rhoSgl(:,2*c+4) * abs(gdot(:,2*c-1)) ! positive mobile --> negative immobile
rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) & rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) &
@ -1123,7 +1123,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
!*** 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%burgers & rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%b_sl &
* ( 2.0_pReal * (rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) + rhoSgl(:,2*c) * abs(gdot(:,2*c-1))) & ! was single hitting single * ( 2.0_pReal * (rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) + rhoSgl(:,2*c) * abs(gdot(:,2*c-1))) & ! was single hitting single
+ 2.0_pReal * (abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single + 2.0_pReal * (abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single
+ rhoDip(:,c) * (abs(gdot(:,2*c-1)) + abs(gdot(:,2*c)))) ! single knocks dipole constituent + rhoDip(:,c) * (abs(gdot(:,2*c-1)) + abs(gdot(:,2*c)))) ! single knocks dipole constituent
@ -1132,13 +1132,13 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
if (lattice_structure(ph) == LATTICE_fcc_ID) & if (lattice_structure(ph) == LATTICE_fcc_ID) &
forall (s = 1:ns, prm%colinearSystem(s) > 0) & forall (s = 1:ns, prm%colinearSystem(s) > 0) &
rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) & rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) &
* 0.25_pReal * sqrt(stt%rho_forest(s,of)) * (dUpper(s,2) + dLower(s,2)) * prm%edgeJogFactor * 0.25_pReal * sqrt(stt%rho_forest(s,of)) * (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%Dsd0 * exp(-prm%selfDiffusionEnergy / (kB * Temperature)) selfDiffusion = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature))
vClimb = prm%atomicVolume * selfDiffusion * prm%mu & vClimb = prm%V_at * selfDiffusion * prm%mu &
/ ( kB * Temperature * PI * (1.0_pReal-prm%nu) * (dUpper(:,1) + dLower(:,1))) / ( kB * Temperature * PI * (1.0_pReal-prm%nu) * (dUpper(:,1) + dLower(:,1)))
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) * vClimb(s) / (dUpper(s,1) - dLower(s,1)), &
@ -1252,7 +1252,7 @@ function rhoDotFlux(F,Fp,timestep, instance,of,ip,el)
my_rhoSgl0 = rho0(:,sgl) my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) !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,instance),of) !ToDo: MD: I think we should use state0 here
gdot = rhoSgl(:,1:4) * v * spread(prm%burgers,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,instance),of) forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),of)
@ -1264,14 +1264,14 @@ function rhoDotFlux(F,Fp,timestep, instance,of,ip,el)
!*** 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%CFLfactor * abs(v0) * timestep & .and. prm%f_c * 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%CFLfactor * abs(v0) * timestep & .and. prm%f_c * 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 !!!'
@ -1334,8 +1334,8 @@ function rhoDotFlux(F,Fp,timestep, instance,of,ip,el)
neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_instance),no),0.0_pReal) neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_instance),no),0.0_pReal)
endforall endforall
where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN & where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%rho_min &
.or. neighbor_rhoSgl0 < prm%significantRho) & .or. neighbor_rhoSgl0 < prm%rho_significant) &
neighbor_rhoSgl0 = 0.0_pReal neighbor_rhoSgl0 = 0.0_pReal
normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), & normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), &
IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! normal of the interface in (average) deformed configuration (pointing neighbor => me) IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! normal of the interface in (average) deformed configuration (pointing neighbor => me)
@ -1460,7 +1460,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,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 !* Set surface transmissivity to the value specified in the material.config
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%surfaceTransmissivity) 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 we encounter a different nonlocal phase at the neighbor,
@ -1469,13 +1469,13 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e)
!* we do not consider this to be a phase boundary, so completely compatible. !* 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 (.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%grainboundaryTransmissivity >= 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) !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config)
if (any(dNeq(material_orientation0(1,i,e)%asQuaternion(), & if (any(dNeq(material_orientation0(1,i,e)%asQuaternion(), &
material_orientation0(1,neighbor_i,neighbor_e)%asQuaternion())) .and. & material_orientation0(1,neighbor_i,neighbor_e)%asQuaternion())) .and. &
(.not. phase_localPlasticity(neighbor_phase))) & (.not. phase_localPlasticity(neighbor_phase))) &
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%grainboundaryTransmissivity) forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB)
else else
!* GRAIN BOUNDARY ? !* GRAIN BOUNDARY ?
!* Compatibility defined by relative orientation of slip systems: !* Compatibility defined by relative orientation of slip systems:
@ -1631,7 +1631,7 @@ subroutine stateInit(ini,phase,NipcMyPhase,instance)
associate(stt => state(instance)) associate(stt => state(instance))
if (ini%rhoSglRandom > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume if (ini%random_rho_u > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume
do e = 1,discretization_nElem do e = 1,discretization_nElem
do i = 1,discretization_nIP do i = 1,discretization_nIP
if (material_phaseAt(1,e) == phase) volume(material_phasememberAt(1,i,e)) = IPvolume(i,e) if (material_phaseAt(1,e) == phase) volume(material_phasememberAt(1,i,e)) = IPvolume(i,e)
@ -1639,11 +1639,11 @@ subroutine stateInit(ini,phase,NipcMyPhase,instance)
enddo enddo
totalVolume = sum(volume) totalVolume = sum(volume)
minimumIPVolume = minval(volume) minimumIPVolume = minval(volume)
densityBinning = ini%rhoSglRandomBinning / minimumIpVolume ** (2.0_pReal / 3.0_pReal) densityBinning = ini%random_rho_u_binning / minimumIpVolume ** (2.0_pReal / 3.0_pReal)
! fill random material points with dislocation segments until the desired overall density is reached ! fill random material points with dislocation segments until the desired overall density is reached
meanDensity = 0.0_pReal meanDensity = 0.0_pReal
do while(meanDensity < ini%rhoSglRandom) do while(meanDensity < ini%random_rho_u)
call random_number(rnd) call random_number(rnd)
phasemember = nint(rnd(1)*real(NipcMyPhase,pReal) + 0.5_pReal) phasemember = nint(rnd(1)*real(NipcMyPhase,pReal) + 0.5_pReal)
s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal) s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal)
@ -1656,15 +1656,15 @@ subroutine stateInit(ini,phase,NipcMyPhase,instance)
from = 1 + sum(ini%N_sl(1:f-1)) from = 1 + sum(ini%N_sl(1:f-1))
upto = sum(ini%N_sl(1:f)) upto = sum(ini%N_sl(1:f))
do s = from,upto do s = from,upto
noise = [math_sampleGaussVar(0.0_pReal, ini%rhoSglScatter), & noise = [math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u), &
math_sampleGaussVar(0.0_pReal, ini%rhoSglScatter)] math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u)]
stt%rho_sgl_mob_edg_pos(s,e) = ini%rhoSglEdgePos0(f) + noise(1) stt%rho_sgl_mob_edg_pos(s,e) = ini%rho_u_ed_pos_0(f) + noise(1)
stt%rho_sgl_mob_edg_neg(s,e) = ini%rhoSglEdgeNeg0(f) + noise(1) stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1)
stt%rho_sgl_mob_scr_pos(s,e) = ini%rhoSglScrewPos0(f) + noise(2) stt%rho_sgl_mob_scr_pos(s,e) = ini%rho_u_sc_pos_0(f) + noise(2)
stt%rho_sgl_mob_scr_neg(s,e) = ini%rhoSglScrewNeg0(f) + noise(2) stt%rho_sgl_mob_scr_neg(s,e) = ini%rho_u_sc_neg_0(f) + noise(2)
enddo enddo
stt%rho_dip_edg(from:upto,e) = ini%rhoDipEdge0(f) stt%rho_dip_edg(from:upto,e) = ini%rho_d_ed_0(f)
stt%rho_dip_scr(from:upto,e) = ini%rhoDipScrew0(f) stt%rho_dip_scr(from:upto,e) = ini%rho_d_sc_0(f)
enddo enddo
enddo enddo
endif endif
@ -1732,14 +1732,14 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem
!* Effective stress includes non Schmid constributions !* 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 !* 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 tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) ! ensure that the effective stress is positive
meanfreepath_P = prm%burgers(s) meanfreepath_P = prm%b_sl(s)
jumpWidth_P = prm%burgers(s) jumpWidth_P = prm%b_sl(s)
activationLength_P = prm%doublekinkwidth *prm%burgers(s) activationLength_P = prm%w *prm%b_sl(s)
activationVolume_P = activationLength_P * jumpWidth_P * prm%burgers(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) ! ensure that the activation probability cannot become greater than one
tPeierls = 1.0_pReal / prm%fattack & tPeierls = 1.0_pReal / prm%nu_a &
* exp(activationEnergy_P / (kB * Temperature) & * exp(activationEnergy_P / (kB * Temperature) &
* (1.0_pReal - tauRel_P**prm%p)**prm%q) * (1.0_pReal - tauRel_P**prm%p)**prm%q)
if (tauEff < criticalStress_P) then if (tauEff < criticalStress_P) then
@ -1752,14 +1752,14 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem
!* 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 !* 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%burgers(s) / sqrt(prm%solidSolutionConcentration) meanfreepath_S = prm%b_sl(s) / sqrt(prm%c_sol)
jumpWidth_S = prm%solidSolutionSize * prm%burgers(s) jumpWidth_S = prm%f_sol * prm%b_sl(s)
activationLength_S = prm%burgers(s) / sqrt(prm%solidSolutionConcentration) activationLength_S = prm%b_sl(s) / sqrt(prm%c_sol)
activationVolume_S = activationLength_S * jumpWidth_S * prm%burgers(s) activationVolume_S = activationLength_S * jumpWidth_S * prm%b_sl(s)
activationEnergy_S = prm%solidSolutionEnergy activationEnergy_S = prm%Q_sol
criticalStress_S = activationEnergy_S / activationVolume_S 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 tauRel_S = min(1.0_pReal, tauEff / criticalStress_S) ! ensure that the activation probability cannot become greater than one
tSolidSolution = 1.0_pReal / prm%fattack & tSolidSolution = 1.0_pReal / prm%nu_a &
* exp(activationEnergy_S / (kB * Temperature)* (1.0_pReal - tauRel_S**prm%p)**prm%q) * exp(activationEnergy_S / (kB * Temperature)* (1.0_pReal - tauRel_S**prm%p)**prm%q)
if (tauEff < criticalStress_S) then if (tauEff < criticalStress_S) then
dtSolidSolution_dtau = tSolidSolution * prm%p * prm%q * activationVolume_S / (kB * Temperature) & dtSolidSolution_dtau = tSolidSolution * prm%p * prm%q * activationVolume_S / (kB * Temperature) &
@ -1770,7 +1770,7 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem
!* viscous glide velocity !* viscous glide velocity
tauEff = abs(tau(s)) - tauThreshold(s) tauEff = abs(tau(s)) - tauThreshold(s)
mobility = prm%burgers(s) / prm%viscosity mobility = prm%b_sl(s) / prm%eta
vViscous = mobility * tauEff vViscous = mobility * tauEff
!* Mean velocity results from waiting time at peierls barriers and solid solution obstacles with respective meanfreepath of !* Mean velocity results from waiting time at peierls barriers and solid solution obstacles with respective meanfreepath of
@ -1805,7 +1805,7 @@ pure function getRho(instance,of,ip,el)
getRho(:,mob) = max(getRho(:,mob),0.0_pReal) getRho(:,mob) = max(getRho(:,mob),0.0_pReal)
getRho(:,dip) = max(getRho(:,dip),0.0_pReal) getRho(:,dip) = max(getRho(:,dip),0.0_pReal)
where(abs(getRho) < max(prm%significantN/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%significantRho)) & where(abs(getRho) < max(prm%rho_min/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
getRho = 0.0_pReal getRho = 0.0_pReal
end associate end associate
@ -1830,7 +1830,7 @@ pure function getRho0(instance,of,ip,el)
getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal) getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal)
getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal) getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal)
where(abs(getRho0) < max(prm%significantN/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%significantRho)) & where(abs(getRho0) < max(prm%rho_min/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
getRho0 = 0.0_pReal getRho0 = 0.0_pReal
end associate end associate

160
src/constitutive_plastic_phenopowerlaw.f90
View File

@ -8,28 +8,28 @@ submodule(constitutive:constitutive_plastic) plastic_phenopowerlaw
type :: tParameters type :: tParameters
real(pReal) :: & real(pReal) :: &
gdot0_slip = 1.0_pReal, & !< reference shear strain rate for slip dot_gamma_0_sl = 1.0_pReal, & !< reference shear strain rate for slip
gdot0_twin = 1.0_pReal, & !< reference shear strain rate for twin dot_gamma_0_tw = 1.0_pReal, & !< reference shear strain rate for twin
n_slip = 1.0_pReal, & !< stress exponent for slip n_sl = 1.0_pReal, & !< stress exponent for slip
n_twin = 1.0_pReal, & !< stress exponent for twin n_tw = 1.0_pReal, & !< stress exponent for twin
spr = 1.0_pReal, & !< push-up factor for slip saturation due to twinning f_sl_sat_tw = 1.0_pReal, & !< push-up factor for slip saturation due to twinning
c_1 = 1.0_pReal, & c_1 = 1.0_pReal, &
c_2 = 1.0_pReal, & c_2 = 1.0_pReal, &
c_3 = 1.0_pReal, & c_3 = 1.0_pReal, &
c_4 = 1.0_pReal, & c_4 = 1.0_pReal, &
h0_SlipSlip = 1.0_pReal, & !< reference hardening slip - slip h_0_sl_sl = 1.0_pReal, & !< reference hardening slip - slip
h0_TwinSlip = 1.0_pReal, & !< reference hardening twin - slip h_0_tw_sl = 1.0_pReal, & !< reference hardening twin - slip
h0_TwinTwin = 1.0_pReal, & !< reference hardening twin - twin h_0_tw_tw = 1.0_pReal, & !< reference hardening twin - twin
a_slip = 1.0_pReal a_sl = 1.0_pReal
real(pReal), allocatable, dimension(:) :: & real(pReal), allocatable, dimension(:) :: &
xi_slip_sat, & !< maximum critical shear stress for slip xi_inf_sl, & !< maximum critical shear stress for slip
H_int, & !< per family hardening activity (optional) h_int, & !< per family hardening activity (optional)
gamma_twin_char !< characteristic shear for twins gamma_tw_char !< characteristic shear for twins
real(pReal), allocatable, dimension(:,:) :: & real(pReal), allocatable, dimension(:,:) :: &
interaction_SlipSlip, & !< slip resistance from slip activity h_sl_sl, & !< slip resistance from slip activity
interaction_SlipTwin, & !< slip resistance from twin activity h_sl_tw, & !< slip resistance from twin activity
interaction_TwinSlip, & !< twin resistance from slip activity h_tw_sl, & !< twin resistance from slip activity
interaction_TwinTwin !< twin resistance from twin activity h_tw_tw !< twin resistance from twin activity
real(pReal), allocatable, dimension(:,:,:) :: & real(pReal), allocatable, dimension(:,:,:) :: &
P_sl, & P_sl, &
P_tw, & P_tw, &
@ -78,8 +78,8 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
N_sl, N_tw N_sl, N_tw
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
xi_slip_0, & !< initial critical shear stress for slip xi_0_sl, & !< initial critical shear stress for slip
xi_twin_0, & !< initial critical shear stress for twin xi_0_tw, & !< initial critical shear stress for twin
a !< non-Schmid coefficients a !< non-Schmid coefficients
character(len=pStringLen) :: & character(len=pStringLen) :: &
extmsg = '' extmsg = ''
@ -120,7 +120,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
phase%get_asFloat('c/a',defaultVal=0.0_pReal)) phase%get_asFloat('c/a',defaultVal=0.0_pReal))
if(phase%get_asString('lattice') == 'bcc') then if(phase%get_asString('lattice') == 'bcc') then
a = pl%get_asFloats('nonSchmid_coefficients',defaultVal=emptyRealArray) a = pl%get_asFloats('a_nonSchmid',defaultVal=emptyRealArray)
if(size(a) > 0) prm%nonSchmidActive = .true. if(size(a) > 0) prm%nonSchmidActive = .true.
prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1) prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1) prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
@ -128,36 +128,36 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
prm%nonSchmid_pos = prm%P_sl prm%nonSchmid_pos = prm%P_sl
prm%nonSchmid_neg = prm%P_sl prm%nonSchmid_neg = prm%P_sl
endif endif
prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(N_sl, & prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl, &
pl%get_asFloats('h_sl_sl'), & pl%get_asFloats('h_sl_sl'), &
phase%get_asString('lattice')) phase%get_asString('lattice'))
xi_slip_0 = pl%get_asFloats('xi_0_sl', requiredSize=size(N_sl)) xi_0_sl = pl%get_asFloats('xi_0_sl', requiredSize=size(N_sl))
prm%xi_slip_sat = pl%get_asFloats('xi_inf_sl', requiredSize=size(N_sl)) prm%xi_inf_sl = pl%get_asFloats('xi_inf_sl', requiredSize=size(N_sl))
prm%H_int = pl%get_asFloats('h_int', requiredSize=size(N_sl), & prm%h_int = pl%get_asFloats('h_int', requiredSize=size(N_sl), &
defaultVal=[(0.0_pReal,i=1,size(N_sl))]) defaultVal=[(0.0_pReal,i=1,size(N_sl))])
prm%gdot0_slip = pl%get_asFloat('dot_gamma_0_sl') prm%dot_gamma_0_sl = pl%get_asFloat('dot_gamma_0_sl')
prm%n_slip = pl%get_asFloat('n_sl') prm%n_sl = pl%get_asFloat('n_sl')
prm%a_slip = pl%get_asFloat('a_sl') prm%a_sl = pl%get_asFloat('a_sl')
prm%h0_SlipSlip = pl%get_asFloat('h_0_sl_sl') prm%h_0_sl_sl = pl%get_asFloat('h_0_sl_sl')
! expand: family => system ! expand: family => system
xi_slip_0 = math_expand(xi_slip_0, N_sl) xi_0_sl = math_expand(xi_0_sl, N_sl)
prm%xi_slip_sat = math_expand(prm%xi_slip_sat,N_sl) prm%xi_inf_sl = math_expand(prm%xi_inf_sl,N_sl)
prm%H_int = math_expand(prm%H_int, N_sl) prm%h_int = math_expand(prm%h_int, N_sl)
! sanity checks ! sanity checks
if ( prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0_sl' if ( prm%dot_gamma_0_sl <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0_sl'
if ( prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_sl' if ( prm%a_sl <= 0.0_pReal) extmsg = trim(extmsg)//' a_sl'
if ( prm%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//' n_sl' if ( prm%n_sl <= 0.0_pReal) extmsg = trim(extmsg)//' n_sl'
if (any(xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_0_sl' if (any(xi_0_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_0_sl'
if (any(prm%xi_slip_sat <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_sl' if (any(prm%xi_inf_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_inf_sl'
else slipActive else slipActive
xi_slip_0 = emptyRealArray xi_0_sl = emptyRealArray
allocate(prm%xi_slip_sat,prm%H_int,source=emptyRealArray) allocate(prm%xi_inf_sl,prm%h_int,source=emptyRealArray)
allocate(prm%interaction_SlipSlip(0,0)) allocate(prm%h_sl_sl(0,0))
endif slipActive endif slipActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -167,50 +167,50 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
twinActive: if (prm%sum_N_tw > 0) then twinActive: if (prm%sum_N_tw > 0) then
prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase%get_asString('lattice'),& prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal)) phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%interaction_TwinTwin = lattice_interaction_TwinByTwin(N_tw,& prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,&
pl%get_asFloats('h_tw_tw'), & pl%get_asFloats('h_tw_tw'), &
phase%get_asString('lattice')) phase%get_asString('lattice'))
prm%gamma_twin_char = lattice_characteristicShear_twin(N_tw,phase%get_asString('lattice'),& prm%gamma_tw_char = lattice_characteristicShear_twin(N_tw,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal)) phase%get_asFloat('c/a',defaultVal=0.0_pReal))
xi_twin_0 = pl%get_asFloats('xi_0_tw',requiredSize=size(N_tw)) xi_0_tw = pl%get_asFloats('xi_0_tw',requiredSize=size(N_tw))
prm%c_1 = pl%get_asFloat('c_1',defaultVal=0.0_pReal) prm%c_1 = pl%get_asFloat('c_1',defaultVal=0.0_pReal)
prm%c_2 = pl%get_asFloat('c_2',defaultVal=1.0_pReal) prm%c_2 = pl%get_asFloat('c_2',defaultVal=1.0_pReal)
prm%c_3 = pl%get_asFloat('c_3',defaultVal=0.0_pReal) prm%c_3 = pl%get_asFloat('c_3',defaultVal=0.0_pReal)
prm%c_4 = pl%get_asFloat('c_4',defaultVal=0.0_pReal) prm%c_4 = pl%get_asFloat('c_4',defaultVal=0.0_pReal)
prm%gdot0_twin = pl%get_asFloat('dot_gamma_0_tw') prm%dot_gamma_0_tw = pl%get_asFloat('dot_gamma_0_tw')
prm%n_twin = pl%get_asFloat('n_tw') prm%n_tw = pl%get_asFloat('n_tw')
prm%spr = pl%get_asFloat('f_sl_sat_tw') prm%f_sl_sat_tw = pl%get_asFloat('f_sl_sat_tw')
prm%h0_TwinTwin = pl%get_asFloat('h_0_tw_tw') prm%h_0_tw_tw = pl%get_asFloat('h_0_tw_tw')
! expand: family => system ! expand: family => system
xi_twin_0 = math_expand(xi_twin_0,N_tw) xi_0_tw = math_expand(xi_0_tw,N_tw)
! sanity checks ! sanity checks
if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0_tw' if (prm%dot_gamma_0_tw <= 0.0_pReal) extmsg = trim(extmsg)//' dot_gamma_0_tw'
if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//' n_tw' if (prm%n_tw <= 0.0_pReal) extmsg = trim(extmsg)//' n_tw'
else twinActive else twinActive
xi_twin_0 = emptyRealArray xi_0_tw = emptyRealArray
allocate(prm%gamma_twin_char,source=emptyRealArray) allocate(prm%gamma_tw_char,source=emptyRealArray)
allocate(prm%interaction_TwinTwin(0,0)) allocate(prm%h_tw_tw(0,0))
endif twinActive endif twinActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! slip-twin related parameters ! slip-twin related parameters
slipAndTwinActive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then slipAndTwinActive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then
prm%h0_TwinSlip = pl%get_asFloat('h_0_tw_sl') prm%h_0_tw_sl = pl%get_asFloat('h_0_tw_sl')
prm%interaction_SlipTwin = lattice_interaction_SlipByTwin(N_sl,N_tw,& prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,&
pl%get_asFloats('h_sl_tw'), & pl%get_asFloats('h_sl_tw'), &
phase%get_asString('lattice')) phase%get_asString('lattice'))
prm%interaction_TwinSlip = lattice_interaction_TwinBySlip(N_tw,N_sl,& prm%h_tw_sl = lattice_interaction_TwinBySlip(N_tw,N_sl,&
pl%get_asFloats('h_tw_sl'), & pl%get_asFloats('h_tw_sl'), &
phase%get_asString('lattice')) phase%get_asString('lattice'))
else slipAndTwinActive else slipAndTwinActive
allocate(prm%interaction_SlipTwin(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0 allocate(prm%h_sl_tw(prm%sum_N_sl,prm%sum_N_tw)) ! at least one dimension is 0
allocate(prm%interaction_TwinSlip(prm%sum_N_tw,prm%sum_N_sl)) ! at least one dimension is 0 allocate(prm%h_tw_sl(prm%sum_N_tw,prm%sum_N_sl)) ! at least one dimension is 0
prm%h0_TwinSlip = 0.0_pReal prm%h_0_tw_sl = 0.0_pReal
endif slipAndTwinActive endif slipAndTwinActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -237,7 +237,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:) stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_slip = spread(xi_slip_0, 2, NipcMyPhase) stt%xi_slip = spread(xi_0_sl, 2, NipcMyPhase)
dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:) dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
@ -245,7 +245,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw endIndex = endIndex + prm%sum_N_tw
stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:) stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_twin = spread(xi_twin_0, 2, NipcMyPhase) stt%xi_twin = spread(xi_0_tw, 2, NipcMyPhase)
dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
@ -354,20 +354,20 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
sumGamma = sum(stt%gamma_slip(:,of)) sumGamma = sum(stt%gamma_slip(:,of))
sumF = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char) sumF = sum(stt%gamma_twin(:,of)/prm%gamma_tw_char)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices ! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices
c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%c_1*sumF** prm%c_2) c_SlipSlip = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2)
c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%c_3 c_TwinSlip = prm%h_0_tw_sl * sumGamma**prm%c_3
c_TwinTwin = prm%h0_TwinTwin * sumF**prm%c_4 c_TwinTwin = prm%h_0_tw_tw * sumF**prm%c_4
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate left and right vectors ! calculate left and right vectors
left_SlipSlip = 1.0_pReal + prm%H_int left_SlipSlip = 1.0_pReal + prm%h_int
xi_slip_sat_offset = prm%spr*sqrt(sumF) xi_slip_sat_offset = prm%f_sl_sat_tw*sqrt(sumF)
right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip & right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_inf_sl+xi_slip_sat_offset)) **prm%a_sl &
* sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_inf_sl+xi_slip_sat_offset))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! shear rates ! shear rates
@ -378,11 +378,11 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! hardening ! hardening
dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * & dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * &
matmul(prm%interaction_SlipSlip,dot%gamma_slip(:,of)*right_SlipSlip) & matmul(prm%h_sl_sl,dot%gamma_slip(:,of)*right_SlipSlip) &
+ matmul(prm%interaction_SlipTwin,dot%gamma_twin(:,of)) + matmul(prm%h_sl_tw,dot%gamma_twin(:,of))
dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%interaction_TwinSlip,dot%gamma_slip(:,of)) & dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%h_tw_sl,dot%gamma_slip(:,of)) &
+ c_TwinTwin * matmul(prm%interaction_TwinTwin,dot%gamma_twin(:,of)) + c_TwinTwin * matmul(prm%h_tw_tw,dot%gamma_twin(:,of))
end associate end associate
end subroutine plastic_phenopowerlaw_dotState end subroutine plastic_phenopowerlaw_dotState
@ -460,29 +460,29 @@ pure subroutine kinetics_slip(Mp,instance,of, &
enddo enddo
where(dNeq0(tau_slip_pos)) where(dNeq0(tau_slip_pos))
gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active gdot_slip_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active
* sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos) * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_sl, tau_slip_pos)
else where else where
gdot_slip_pos = 0.0_pReal gdot_slip_pos = 0.0_pReal
end where end where
where(dNeq0(tau_slip_neg)) where(dNeq0(tau_slip_neg))
gdot_slip_neg = prm%gdot0_slip * 0.5_pReal & ! only used if non-Schmid active, always 1/2 gdot_slip_neg = prm%dot_gamma_0_sl * 0.5_pReal & ! only used if non-Schmid active, always 1/2
* sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip, tau_slip_neg) * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_sl, tau_slip_neg)
else where else where
gdot_slip_neg = 0.0_pReal gdot_slip_neg = 0.0_pReal
end where end where
if (present(dgdot_dtau_slip_pos)) then if (present(dgdot_dtau_slip_pos)) then
where(dNeq0(gdot_slip_pos)) where(dNeq0(gdot_slip_pos))
dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_sl/tau_slip_pos
else where else where
dgdot_dtau_slip_pos = 0.0_pReal dgdot_dtau_slip_pos = 0.0_pReal
end where end where
endif endif
if (present(dgdot_dtau_slip_neg)) then if (present(dgdot_dtau_slip_neg)) then
where(dNeq0(gdot_slip_neg)) where(dNeq0(gdot_slip_neg))
dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_slip/tau_slip_neg dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_sl/tau_slip_neg
else where else where
dgdot_dtau_slip_neg = 0.0_pReal dgdot_dtau_slip_neg = 0.0_pReal
end where end where
@ -524,15 +524,15 @@ pure subroutine kinetics_twin(Mp,instance,of,&
enddo enddo
where(tau_twin > 0.0_pReal) where(tau_twin > 0.0_pReal)
gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_tw_char)) & ! only twin in untwinned volume fraction
* prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin * prm%dot_gamma_0_tw*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_tw
else where else where
gdot_twin = 0.0_pReal gdot_twin = 0.0_pReal
end where end where
if (present(dgdot_dtau_twin)) then if (present(dgdot_dtau_twin)) then
where(dNeq0(gdot_twin)) where(dNeq0(gdot_twin))
dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin dgdot_dtau_twin = gdot_twin*prm%n_tw/tau_twin
else where else where
dgdot_dtau_twin = 0.0_pReal dgdot_dtau_twin = 0.0_pReal
end where end where

4
src/damage_nonlocal.f90
View File

@ -134,7 +134,7 @@ function damage_nonlocal_getDiffusion(ip,el)
damage_nonlocal_getDiffusion = 0.0_pReal damage_nonlocal_getDiffusion = 0.0_pReal
do grain = 1, homogenization_Ngrains(homog) do grain = 1, homogenization_Ngrains(homog)
damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + & damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + &
crystallite_push33ToRef(grain,ip,el,lattice_DamageDiffusion(1:3,1:3,material_phaseAt(grain,el))) crystallite_push33ToRef(grain,ip,el,lattice_D(1:3,1:3,material_phaseAt(grain,el)))
enddo enddo
damage_nonlocal_getDiffusion = & damage_nonlocal_getDiffusion = &
@ -157,7 +157,7 @@ real(pReal) function damage_nonlocal_getMobility(ip,el)
damage_nonlocal_getMobility = 0.0_pReal damage_nonlocal_getMobility = 0.0_pReal
do ipc = 1, homogenization_Ngrains(material_homogenizationAt(el)) do ipc = 1, homogenization_Ngrains(material_homogenizationAt(el))
damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_DamageMobility(material_phaseAt(ipc,el)) damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_M(material_phaseAt(ipc,el))
enddo enddo
damage_nonlocal_getMobility = damage_nonlocal_getMobility/& damage_nonlocal_getMobility = damage_nonlocal_getMobility/&

6
src/grid/discretization_grid.f90
View File

@ -107,9 +107,9 @@ subroutine discretization_grid_init(restart)
call discretization_init(materialAt, & call discretization_init(materialAt, &
IPcoordinates0(myGrid,mySize,grid3Offset), & IPcoordinates0(myGrid,mySize,grid3Offset), &
Nodes0(myGrid,mySize,grid3Offset),& Nodes0(myGrid,mySize,grid3Offset),&
merge((grid(1)+1) * (grid(2)+1) * (grid3+1),& ! write bottom layer merge((grid(1)+1) * (grid(2)+1) * (grid3+1),& ! write top layer...
(grid(1)+1) * (grid(2)+1) * grid3,& ! do not write bottom layer (is top of rank-1) (grid(1)+1) * (grid(2)+1) * grid3,& ! ...unless not last process
worldrank<1)) worldrank+1==worldsize))
FEsolving_execElem = [1,product(myGrid)] ! parallel loop bounds set to comprise all elements FEsolving_execElem = [1,product(myGrid)] ! parallel loop bounds set to comprise all elements
FEsolving_execIP = [1,1] ! parallel loop bounds set to comprise the only IP FEsolving_execIP = [1,1] ! parallel loop bounds set to comprise the only IP

88
src/homogenization_mech_RGC.f90
View File

@ -4,20 +4,20 @@
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Relaxed grain cluster (RGC) homogenization scheme !> @brief Relaxed grain cluster (RGC) homogenization scheme
!> Nconstituents is defined as p x q x r (cluster) !> N_constituents is defined as p x q x r (cluster)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(homogenization) homogenization_mech_RGC submodule(homogenization) homogenization_mech_RGC
use rotations use rotations
type :: tParameters type :: tParameters
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
Nconstituents N_constituents
real(pReal) :: & real(pReal) :: &
xiAlpha, & xi_alpha, &
ciAlpha c_Alpha
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
dAlpha, & D_alpha, &
angles a_g
integer :: & integer :: &
of_debug = 0 of_debug = 0
character(len=pStringLen), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
@ -163,20 +163,20 @@ module subroutine mech_RGC_init(num_homogMech)
prm%output = homogMech%get_asStrings('output',defaultVal=emptyStringArray) prm%output = homogMech%get_asStrings('output',defaultVal=emptyStringArray)
#endif #endif
prm%Nconstituents = homogMech%get_asInts('cluster_size',requiredSize=3) prm%N_constituents = homogMech%get_asInts('cluster_size',requiredSize=3)
if (homogenization_Ngrains(h) /= product(prm%Nconstituents)) & if (homogenization_Ngrains(h) /= product(prm%N_constituents)) &
call IO_error(211,ext_msg='clustersize (mech_rgc)') call IO_error(211,ext_msg='clustersize (mech_rgc)')
prm%xiAlpha = homogMech%get_asFloat('xi_alpha') prm%xi_alpha = homogMech%get_asFloat('xi_alpha')
prm%ciAlpha = homogMech%get_asFloat('c_alpha') prm%c_alpha = homogMech%get_asFloat('c_alpha')
prm%dAlpha = homogMech%get_asFloats('D_alpha', requiredSize=3) prm%D_alpha = homogMech%get_asFloats('D_alpha', requiredSize=3)
prm%angles = homogMech%get_asFloats('a_g', requiredSize=3) prm%a_g = homogMech%get_asFloats('a_g', requiredSize=3)
NofMyHomog = count(material_homogenizationAt == h) NofMyHomog = count(material_homogenizationAt == h)
nIntFaceTot = 3*( (prm%Nconstituents(1)-1)*prm%Nconstituents(2)*prm%Nconstituents(3) & nIntFaceTot = 3*( (prm%N_constituents(1)-1)*prm%N_constituents(2)*prm%N_constituents(3) &
+ prm%Nconstituents(1)*(prm%Nconstituents(2)-1)*prm%Nconstituents(3) & + prm%N_constituents(1)*(prm%N_constituents(2)-1)*prm%N_constituents(3) &
+ prm%Nconstituents(1)*prm%Nconstituents(2)*(prm%Nconstituents(3)-1)) + prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1))
sizeState = nIntFaceTot & sizeState = nIntFaceTot &
+ size(['avg constitutive work ','average penalty energy']) + size(['avg constitutive work ','average penalty energy'])
@ -197,8 +197,8 @@ module subroutine mech_RGC_init(num_homogMech)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! assigning cluster orientations ! assigning cluster orientations
dependentState(homogenization_typeInstance(h))%orientation = spread(eu2om(prm%angles*inRad),3,NofMyHomog) dependentState(homogenization_typeInstance(h))%orientation = spread(eu2om(prm%a_g*inRad),3,NofMyHomog)
!dst%orientation = spread(eu2om(prm%angles*inRad),3,NofMyHomog) ifort version 18.0.1 crashes (for whatever reason) !dst%orientation = spread(eu2om(prm%a_g*inRad),3,NofMyHomog) ifort version 18.0.1 crashes (for whatever reason)
end associate end associate
@ -229,8 +229,8 @@ module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! compute the deformation gradient of individual grains due to relaxations ! compute the deformation gradient of individual grains due to relaxations
F = 0.0_pReal F = 0.0_pReal
do iGrain = 1,product(prm%Nconstituents) do iGrain = 1,product(prm%N_constituents)
iGrain3 = grain1to3(iGrain,prm%Nconstituents) iGrain3 = grain1to3(iGrain,prm%N_constituents)
do iFace = 1,6 do iFace = 1,6
intFace = getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain intFace = getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain
aVect = relaxationVector(intFace,instance,of) ! get the relaxation vectors for each interface from global relaxation vector array aVect = relaxationVector(intFace,instance,of) ! get the relaxation vectors for each interface from global relaxation vector array
@ -290,7 +290,7 @@ module procedure mech_RGC_updateState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! get the dimension of the cluster (grains and interfaces) ! get the dimension of the cluster (grains and interfaces)
nGDim = prm%Nconstituents nGDim = prm%N_constituents
nGrain = product(nGDim) nGrain = product(nGDim)
nIntFaceTot = (nGDim(1)-1)*nGDim(2)*nGDim(3) & nIntFaceTot = (nGDim(1)-1)*nGDim(2)*nGDim(3) &
+ nGDim(1)*(nGDim(2)-1)*nGDim(3) & + nGDim(1)*(nGDim(2)-1)*nGDim(3) &
@ -324,12 +324,12 @@ module procedure mech_RGC_updateState
!------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------
! computing the residual stress from the balance of traction at all (interior) interfaces ! computing the residual stress from the balance of traction at all (interior) interfaces
do iNum = 1,nIntFaceTot do iNum = 1,nIntFaceTot
faceID = interface1to4(iNum,param(instance)%Nconstituents) ! identifying the interface ID in local coordinate system (4-dimensional index) faceID = interface1to4(iNum,param(instance)%N_constituents) ! identifying the interface ID in local coordinate system (4-dimensional index)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! identify the left/bottom/back grain (-|N) ! identify the left/bottom/back grain (-|N)
iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index) iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index)
iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) iGrN = grain3to1(iGr3N,param(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceN = getInterface(2*faceID(1),iGr3N) intFaceN = getInterface(2*faceID(1),iGr3N)
normN = interfaceNormal(intFaceN,instance,of) normN = interfaceNormal(intFaceN,instance,of)
@ -337,7 +337,7 @@ module procedure mech_RGC_updateState
! identify the right/up/front grain (+|P) ! identify the right/up/front grain (+|P)
iGr3P = iGr3N iGr3P = iGr3N
iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate system (3-dimensional index) iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate system (3-dimensional index)
iGrP = grain3to1(iGr3P,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) iGrP = grain3to1(iGr3P,param(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceP = getInterface(2*faceID(1)-1,iGr3P) intFaceP = getInterface(2*faceID(1)-1,iGr3P)
normP = interfaceNormal(intFaceP,instance,of) normP = interfaceNormal(intFaceP,instance,of)
@ -393,7 +393,7 @@ module procedure mech_RGC_updateState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! compute/update the state for postResult, i.e., all energy densities computed by time-integration ! compute/update the state for postResult, i.e., all energy densities computed by time-integration
do iGrain = 1,product(prm%Nconstituents) do iGrain = 1,product(prm%N_constituents)
do i = 1,3;do j = 1,3 do i = 1,3;do j = 1,3
stt%work(of) = stt%work(of) & stt%work(of) = stt%work(of) &
+ P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal) + P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal)
@ -450,18 +450,18 @@ module procedure mech_RGC_updateState
! ... of the constitutive stress tangent, assembled from dPdF or material constitutive model "smatrix" ! ... of the constitutive stress tangent, assembled from dPdF or material constitutive model "smatrix"
allocate(smatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal) allocate(smatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal)
do iNum = 1,nIntFaceTot do iNum = 1,nIntFaceTot
faceID = interface1to4(iNum,param(instance)%Nconstituents) ! assembling of local dPdF into global Jacobian matrix faceID = interface1to4(iNum,param(instance)%N_constituents) ! assembling of local dPdF into global Jacobian matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! identify the left/bottom/back grain (-|N) ! identify the left/bottom/back grain (-|N)
iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate sytem iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate sytem
iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate into global grain ID iGrN = grain3to1(iGr3N,param(instance)%N_constituents) ! translate into global grain ID
intFaceN = getInterface(2*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system intFaceN = getInterface(2*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system
normN = interfaceNormal(intFaceN,instance,of) normN = interfaceNormal(intFaceN,instance,of)
do iFace = 1,6 do iFace = 1,6
intFaceN = getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface intFaceN = getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface
mornN = interfaceNormal(intFaceN,instance,of) mornN = interfaceNormal(intFaceN,instance,of)
iMun = interface4to1(intFaceN,param(instance)%Nconstituents) ! translate the interfaces ID into local 4-dimensional index iMun = interface4to1(intFaceN,param(instance)%N_constituents) ! translate the interfaces ID into local 4-dimensional index
if (iMun > 0) then ! get the corresponding tangent if (iMun > 0) then ! get the corresponding tangent
do i=1,3; do j=1,3; do k=1,3; do l=1,3 do i=1,3; do j=1,3; do k=1,3; do l=1,3
smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) & smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) &
@ -476,13 +476,13 @@ module procedure mech_RGC_updateState
! identify the right/up/front grain (+|P) ! identify the right/up/front grain (+|P)
iGr3P = iGr3N iGr3P = iGr3N
iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate sytem iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate sytem
iGrP = grain3to1(iGr3P,param(instance)%Nconstituents) ! translate into global grain ID iGrP = grain3to1(iGr3P,param(instance)%N_constituents) ! translate into global grain ID
intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identifying the connecting interface in local coordinate system intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identifying the connecting interface in local coordinate system
normP = interfaceNormal(intFaceP,instance,of) normP = interfaceNormal(intFaceP,instance,of)
do iFace = 1,6 do iFace = 1,6
intFaceP = getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface intFaceP = getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface
mornP = interfaceNormal(intFaceP,instance,of) mornP = interfaceNormal(intFaceP,instance,of)
iMun = interface4to1(intFaceP,param(instance)%Nconstituents) ! translate the interfaces ID into local 4-dimensional index iMun = interface4to1(intFaceP,param(instance)%N_constituents) ! translate the interfaces ID into local 4-dimensional index
if (iMun > 0) then ! get the corresponding tangent if (iMun > 0) then ! get the corresponding tangent
do i=1,3; do j=1,3; do k=1,3; do l=1,3 do i=1,3; do j=1,3; do k=1,3; do l=1,3
smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) & smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) &
@ -522,12 +522,12 @@ module procedure mech_RGC_updateState
! computing the global stress residual array from the perturbed state ! computing the global stress residual array from the perturbed state
p_resid = 0.0_pReal p_resid = 0.0_pReal
do iNum = 1,nIntFaceTot do iNum = 1,nIntFaceTot
faceID = interface1to4(iNum,param(instance)%Nconstituents) ! identifying the interface ID in local coordinate system (4-dimensional index) faceID = interface1to4(iNum,param(instance)%N_constituents) ! identifying the interface ID in local coordinate system (4-dimensional index)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! identify the left/bottom/back grain (-|N) ! identify the left/bottom/back grain (-|N)
iGr3N = faceID(2:4) ! identify the grain ID in local coordinate system (3-dimensional index) iGr3N = faceID(2:4) ! identify the grain ID in local coordinate system (3-dimensional index)
iGrN = grain3to1(iGr3N,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) iGrN = grain3to1(iGr3N,param(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceN = getInterface(2*faceID(1),iGr3N) ! identify the interface ID of the grain intFaceN = getInterface(2*faceID(1),iGr3N) ! identify the interface ID of the grain
normN = interfaceNormal(intFaceN,instance,of) normN = interfaceNormal(intFaceN,instance,of)
@ -535,7 +535,7 @@ module procedure mech_RGC_updateState
! identify the right/up/front grain (+|P) ! identify the right/up/front grain (+|P)
iGr3P = iGr3N iGr3P = iGr3N
iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identify the grain ID in local coordinate system (3-dimensional index) iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identify the grain ID in local coordinate system (3-dimensional index)
iGrP = grain3to1(iGr3P,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index) iGrP = grain3to1(iGr3P,param(instance)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identify the interface ID of the grain intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identify the interface ID of the grain
normP = interfaceNormal(intFaceP,instance,of) normP = interfaceNormal(intFaceP,instance,of)
@ -664,7 +664,7 @@ module procedure mech_RGC_updateState
real(pReal) :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb real(pReal) :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb
real(pReal), parameter :: nDefToler = 1.0e-10_pReal real(pReal), parameter :: nDefToler = 1.0e-10_pReal
nGDim = param(instance)%Nconstituents nGDim = param(instance)%N_constituents
rPen = 0.0_pReal rPen = 0.0_pReal
nMis = 0.0_pReal nMis = 0.0_pReal
@ -685,11 +685,11 @@ module procedure mech_RGC_updateState
!----------------------------------------------------------------------------------------------- !-----------------------------------------------------------------------------------------------
! computing the mismatch and penalty stress tensor of all grains ! computing the mismatch and penalty stress tensor of all grains
grainLoop: do iGrain = 1,product(prm%Nconstituents) grainLoop: do iGrain = 1,product(prm%N_constituents)
Gmoduli = equivalentModuli(iGrain,ip,el) Gmoduli = equivalentModuli(iGrain,ip,el)
muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain
bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector
iGrain3 = grain1to3(iGrain,prm%Nconstituents) ! get the grain ID in local 3-dimensional index (x,y,z)-position iGrain3 = grain1to3(iGrain,prm%N_constituents) ! get the grain ID in local 3-dimensional index (x,y,z)-position
interfaceLoop: do iFace = 1,6 interfaceLoop: do iFace = 1,6
intFace = getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain intFace = getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain
@ -699,7 +699,7 @@ module procedure mech_RGC_updateState
+ int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal)) + int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal))
where(iGNghb3 < 1) iGNghb3 = nGDim where(iGNghb3 < 1) iGNghb3 = nGDim
where(iGNghb3 >nGDim) iGNghb3 = 1 where(iGNghb3 >nGDim) iGNghb3 = 1
iGNghb = grain3to1(iGNghb3,prm%Nconstituents) ! get the ID of the neighboring grain iGNghb = grain3to1(iGNghb3,prm%N_constituents) ! get the ID of the neighboring grain
Gmoduli = equivalentModuli(iGNghb,ip,el) ! collect the shear modulus and Burgers vector of the neighbor Gmoduli = equivalentModuli(iGNghb,ip,el) ! collect the shear modulus and Burgers vector of the neighbor
muGNghb = Gmoduli(1) muGNghb = Gmoduli(1)
bgGNghb = Gmoduli(2) bgGNghb = Gmoduli(2)
@ -728,9 +728,9 @@ module procedure mech_RGC_updateState
!------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------
! compute the stress penalty of all interfaces ! compute the stress penalty of all interfaces
do i = 1,3; do j = 1,3; do k = 1,3; do l = 1,3 do i = 1,3; do j = 1,3; do k = 1,3; do l = 1,3
rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xiAlpha & rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xi_alpha &
*surfCorr(abs(intFace(1)))/prm%dAlpha(abs(intFace(1))) & *surfCorr(abs(intFace(1)))/prm%D_alpha(abs(intFace(1))) &
*cosh(prm%ciAlpha*nDefNorm) & *cosh(prm%c_alpha*nDefNorm) &
*0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) & *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) &
*tanh(nDefNorm/num%xSmoo) *tanh(nDefNorm/num%xSmoo)
enddo; enddo;enddo; enddo enddo; enddo;enddo; enddo
@ -885,8 +885,8 @@ module procedure mech_RGC_updateState
associate(prm => param(instance)) associate(prm => param(instance))
F = 0.0_pReal F = 0.0_pReal
do iGrain = 1,product(prm%Nconstituents) do iGrain = 1,product(prm%N_constituents)
iGrain3 = grain1to3(iGrain,prm%Nconstituents) iGrain3 = grain1to3(iGrain,prm%N_constituents)
do iFace = 1,6 do iFace = 1,6
intFace = getInterface(iFace,iGrain3) intFace = getInterface(iFace,iGrain3)
aVect = relaxationVector(intFace,instance,of) aVect = relaxationVector(intFace,instance,of)
@ -916,8 +916,8 @@ module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ins
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: instance integer, intent(in) :: instance
avgP = sum(P,3) /real(product(param(instance)%Nconstituents),pReal) avgP = sum(P,3) /real(product(param(instance)%N_constituents),pReal)
dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%Nconstituents),pReal) dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%N_constituents),pReal)
end subroutine mech_RGC_averageStressAndItsTangent end subroutine mech_RGC_averageStressAndItsTangent
@ -975,7 +975,7 @@ pure function relaxationVector(intFace,instance,of)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! collect the interface relaxation vector from the global state array ! collect the interface relaxation vector from the global state array
iNum = interface4to1(intFace,param(instance)%Nconstituents) ! identify the position of the interface in global state array iNum = interface4to1(intFace,param(instance)%N_constituents) ! identify the position of the interface in global state array
if (iNum > 0) then if (iNum > 0) then
relaxationVector = state(instance)%relaxationVector((3*iNum-2):(3*iNum),of) relaxationVector = state(instance)%relaxationVector((3*iNum-2):(3*iNum),of)
else else

8
src/homogenization_mech_isostrain.f90
View File

@ -13,7 +13,7 @@ submodule(homogenization) homogenization_mech_isostrain
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
integer :: & integer :: &
Nconstituents N_constituents
integer(kind(average_ID)) :: & integer(kind(average_ID)) :: &
mapping mapping
end type end type
@ -51,7 +51,7 @@ module subroutine mech_isostrain_init
homogMech => homog%get('mech') homogMech => homog%get('mech')
associate(prm => param(homogenization_typeInstance(h))) associate(prm => param(homogenization_typeInstance(h)))
prm%Nconstituents = homogMech%get_asInt('N_constituents') prm%N_constituents = homogMech%get_asInt('N_constituents')
select case(homogMech%get_asString('mapping',defaultVal = 'sum')) select case(homogMech%get_asString('mapping',defaultVal = 'sum'))
case ('sum') case ('sum')
prm%mapping = parallel_ID prm%mapping = parallel_ID
@ -107,8 +107,8 @@ module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dP
avgP = sum(P,3) avgP = sum(P,3)
dAvgPdAvgF = sum(dPdF,5) dAvgPdAvgF = sum(dPdF,5)
case (average_ID) case (average_ID)
avgP = sum(P,3) /real(prm%Nconstituents,pReal) avgP = sum(P,3) /real(prm%N_constituents,pReal)
dAvgPdAvgF = sum(dPdF,5)/real(prm%Nconstituents,pReal) dAvgPdAvgF = sum(dPdF,5)/real(prm%N_constituents,pReal)
end select end select
end associate end associate

34
src/kinematics_cleavage_opening.f90
View File

@ -12,10 +12,10 @@ submodule(constitutive:constitutive_damage) kinematics_cleavage_opening
integer :: & integer :: &
sum_N_cl !< total number of cleavage planes sum_N_cl !< total number of cleavage planes
real(pReal) :: & real(pReal) :: &
sdot0, & !< opening rate of cleavage planes dot_o, & !< opening rate of cleavage planes
n !< damage rate sensitivity q !< damage rate sensitivity
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
critLoad g_crit
real(pReal), dimension(:,:,:,:), allocatable :: & real(pReal), dimension(:,:,:,:), allocatable :: &
cleavage_systems cleavage_systems
end type tParameters end type tParameters
@ -70,21 +70,21 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin
N_cl = kinematic_type%get_asInts('N_cl') N_cl = kinematic_type%get_asInts('N_cl')
prm%sum_N_cl = sum(abs(N_cl)) prm%sum_N_cl = sum(abs(N_cl))
prm%n = kinematic_type%get_asFloat('q') prm%q = kinematic_type%get_asFloat('q')
prm%sdot0 = kinematic_type%get_asFloat('dot_o') prm%dot_o = kinematic_type%get_asFloat('dot_o')
prm%critLoad = kinematic_type%get_asFloats('g_crit',requiredSize=size(N_cl)) prm%g_crit = kinematic_type%get_asFloats('g_crit',requiredSize=size(N_cl))
prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),& prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal)) phase%get_asFloat('c/a',defaultVal=0.0_pReal))
! expand: family => system ! expand: family => system
prm%critLoad = math_expand(prm%critLoad,N_cl) prm%g_crit = math_expand(prm%g_crit,N_cl)
! sanity checks ! sanity checks
if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%sdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o' if (prm%dot_o <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o'
if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit' if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range ! exit if any parameter is out of range
@ -128,13 +128,13 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S,
dLd_dTstar = 0.0_pReal dLd_dTstar = 0.0_pReal
associate(prm => param(kinematics_cleavage_opening_instance(material_phaseAt(ipc,el)))) associate(prm => param(kinematics_cleavage_opening_instance(material_phaseAt(ipc,el))))
do i = 1,prm%sum_N_cl do i = 1,prm%sum_N_cl
traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset)**2.0_pReal traction_crit = prm%g_crit(i)* damage(homog)%p(damageOffset)**2.0_pReal
traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i))
if (abs(traction_d) > traction_crit + tol_math_check) then if (abs(traction_d) > traction_crit + tol_math_check) then
udotd = sign(1.0_pReal,traction_d)* prm%sdot0 * ((abs(traction_d) - traction_crit)/traction_crit)**prm%n udotd = sign(1.0_pReal,traction_d)* prm%dot_o * ((abs(traction_d) - traction_crit)/traction_crit)**prm%q
Ld = Ld + udotd*prm%cleavage_systems(1:3,1:3,1,i) Ld = Ld + udotd*prm%cleavage_systems(1:3,1:3,1,i)
dudotd_dt = sign(1.0_pReal,traction_d)*udotd*prm%n / (abs(traction_d) - traction_crit) dudotd_dt = sign(1.0_pReal,traction_d)*udotd*prm%q / (abs(traction_d) - traction_crit)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & 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) & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
+ dudotd_dt*prm%cleavage_systems(k,l,1,i) * prm%cleavage_systems(m,n,1,i) + dudotd_dt*prm%cleavage_systems(k,l,1,i) * prm%cleavage_systems(m,n,1,i)
@ -142,9 +142,9 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S,
traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i))
if (abs(traction_t) > traction_crit + tol_math_check) then if (abs(traction_t) > traction_crit + tol_math_check) then
udott = sign(1.0_pReal,traction_t)* prm%sdot0 * ((abs(traction_t) - traction_crit)/traction_crit)**prm%n udott = sign(1.0_pReal,traction_t)* prm%dot_o * ((abs(traction_t) - traction_crit)/traction_crit)**prm%q
Ld = Ld + udott*prm%cleavage_systems(1:3,1:3,2,i) Ld = Ld + udott*prm%cleavage_systems(1:3,1:3,2,i)
dudott_dt = sign(1.0_pReal,traction_t)*udott*prm%n / (abs(traction_t) - traction_crit) dudott_dt = sign(1.0_pReal,traction_t)*udott*prm%q / (abs(traction_t) - traction_crit)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & 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) & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
+ dudott_dt*prm%cleavage_systems(k,l,2,i) * prm%cleavage_systems(m,n,2,i) + dudott_dt*prm%cleavage_systems(k,l,2,i) * prm%cleavage_systems(m,n,2,i)
@ -152,9 +152,9 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S,
traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i))
if (abs(traction_n) > traction_crit + tol_math_check) then if (abs(traction_n) > traction_crit + tol_math_check) then
udotn = sign(1.0_pReal,traction_n)* prm%sdot0 * ((abs(traction_n) - traction_crit)/traction_crit)**prm%n udotn = sign(1.0_pReal,traction_n)* prm%dot_o * ((abs(traction_n) - traction_crit)/traction_crit)**prm%q
Ld = Ld + udotn*prm%cleavage_systems(1:3,1:3,3,i) Ld = Ld + udotn*prm%cleavage_systems(1:3,1:3,3,i)
dudotn_dt = sign(1.0_pReal,traction_n)*udotn*prm%n / (abs(traction_n) - traction_crit) dudotn_dt = sign(1.0_pReal,traction_n)*udotn*prm%q / (abs(traction_n) - traction_crit)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & 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) & dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
+ dudotn_dt*prm%cleavage_systems(k,l,3,i) * prm%cleavage_systems(m,n,3,i) + dudotn_dt*prm%cleavage_systems(k,l,3,i) * prm%cleavage_systems(m,n,3,i)

40
src/kinematics_slipplane_opening.f90
View File

@ -12,10 +12,10 @@ submodule(constitutive:constitutive_damage) kinematics_slipplane_opening
integer :: & integer :: &
sum_N_sl !< total number of cleavage planes sum_N_sl !< total number of cleavage planes
real(pReal) :: & real(pReal) :: &
sdot0, & !< opening rate of cleavage planes dot_o, & !< opening rate of cleavage planes
n !< damage rate sensitivity q !< damage rate sensitivity
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
critLoad g_crit
real(pReal), dimension(:,:,:), allocatable :: & real(pReal), dimension(:,:,:), allocatable :: &
P_d, & P_d, &
P_t, & P_t, &
@ -70,8 +70,8 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi
associate(prm => param(kinematics_slipplane_opening_instance(p))) associate(prm => param(kinematics_slipplane_opening_instance(p)))
kinematic_type => kinematics%get(k) kinematic_type => kinematics%get(k)
prm%sdot0 = kinematic_type%get_asFloat('dot_o') prm%dot_o = kinematic_type%get_asFloat('dot_o')
prm%n = kinematic_type%get_asFloat('q') prm%q = kinematic_type%get_asFloat('q')
N_sl = pl%get_asInts('N_sl') N_sl = pl%get_asInts('N_sl')
prm%sum_N_sl = sum(abs(N_sl)) prm%sum_N_sl = sum(abs(N_sl))
@ -89,15 +89,15 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi
prm%P_n(1:3,1:3,i) = math_outer(n(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 enddo
prm%critLoad = kinematic_type%get_asFloats('g_crit',requiredSize=size(N_sl)) prm%g_crit = kinematic_type%get_asFloats('g_crit',requiredSize=size(N_sl))
! expand: family => system ! expand: family => system
prm%critLoad = math_expand(prm%critLoad,N_sl) prm%g_crit = math_expand(prm%g_crit,N_sl)
! sanity checks ! sanity checks
if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_n' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_n'
if (prm%sdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_sdot0' if (prm%dot_o <= 0.0_pReal) extmsg = trim(extmsg)//' anisoDuctile_sdot0'
if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisoDuctile_critLoad' if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' anisoDuctile_critLoad'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range ! exit if any parameter is out of range
@ -150,27 +150,27 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S
traction_t = math_tensordot(S,prm%P_t(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_n = math_tensordot(S,prm%P_n(1:3,1:3,i))
traction_crit = prm%critLoad(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage traction_crit = prm%g_crit(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage
udotd = sign(1.0_pReal,traction_d)* prm%sdot0* ( abs(traction_d)/traction_crit & udotd = sign(1.0_pReal,traction_d)* prm%dot_o* ( abs(traction_d)/traction_crit &
- abs(traction_d)/prm%critLoad(i))**prm%n - abs(traction_d)/prm%g_crit(i))**prm%q
udott = sign(1.0_pReal,traction_t)* prm%sdot0* ( abs(traction_t)/traction_crit & udott = sign(1.0_pReal,traction_t)* prm%dot_o* ( abs(traction_t)/traction_crit &
- abs(traction_t)/prm%critLoad(i))**prm%n - abs(traction_t)/prm%g_crit(i))**prm%q
udotn = prm%sdot0* ( max(0.0_pReal,traction_n)/traction_crit & udotn = prm%dot_o* ( max(0.0_pReal,traction_n)/traction_crit &
- max(0.0_pReal,traction_n)/prm%critLoad(i))**prm%n - max(0.0_pReal,traction_n)/prm%g_crit(i))**prm%q
if (dNeq0(traction_d)) then if (dNeq0(traction_d)) then
dudotd_dt = udotd*prm%n/traction_d dudotd_dt = udotd*prm%q/traction_d
else else
dudotd_dt = 0.0_pReal dudotd_dt = 0.0_pReal
endif endif
if (dNeq0(traction_t)) then if (dNeq0(traction_t)) then
dudott_dt = udott*prm%n/traction_t dudott_dt = udott*prm%q/traction_t
else else
dudott_dt = 0.0_pReal dudott_dt = 0.0_pReal
endif endif
if (dNeq0(traction_n)) then if (dNeq0(traction_n)) then
dudotn_dt = udotn*prm%n/traction_n dudotn_dt = udotn*prm%q/traction_n
else else
dudotn_dt = 0.0_pReal dudotn_dt = 0.0_pReal
endif endif

30
src/kinematics_thermal_expansion.f90
View File

@ -11,7 +11,7 @@ submodule(constitutive:constitutive_thermal) kinematics_thermal_expansion
real(pReal) :: & real(pReal) :: &
T_ref T_ref
real(pReal), dimension(3,3,3) :: & real(pReal), dimension(3,3,3) :: &
expansion = 0.0_pReal A = 0.0_pReal
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param type(tParameters), dimension(:), allocatable :: param
@ -64,13 +64,13 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi
! read up to three parameters (constant, linear, quadratic with T) ! read up to three parameters (constant, linear, quadratic with T)
temp = kinematic_type%get_asFloats('A_11') temp = kinematic_type%get_asFloats('A_11')
prm%expansion(1,1,1:size(temp)) = temp prm%A(1,1,1:size(temp)) = temp
temp = kinematic_type%get_asFloats('A_22',defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp)) temp = kinematic_type%get_asFloats('A_22',defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp))
prm%expansion(2,2,1:size(temp)) = temp prm%A(2,2,1:size(temp)) = temp
temp = kinematic_type%get_asFloats('A_33',defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp)) temp = kinematic_type%get_asFloats('A_33',defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp))
prm%expansion(3,3,1:size(temp)) = temp prm%A(3,3,1:size(temp)) = temp
do i=1, size(prm%expansion,3) do i=1, size(prm%A,3)
prm%expansion(1:3,1:3,i) = lattice_applyLatticeSymmetry33(prm%expansion(1:3,1:3,i),& prm%A(1:3,1:3,i) = lattice_applyLatticeSymmetry33(prm%A(1:3,1:3,i),&
phase%get_asString('lattice')) phase%get_asString('lattice'))
enddo enddo
@ -98,9 +98,9 @@ pure module function kinematics_thermal_expansion_initialStrain(homog,phase,offs
associate(prm => param(kinematics_thermal_expansion_instance(phase))) associate(prm => param(kinematics_thermal_expansion_instance(phase)))
initialStrain = & initialStrain = &
(temperature(homog)%p(offset) - prm%T_ref)**1 / 1. * prm%expansion(1:3,1:3,1) + & ! constant coefficient (temperature(homog)%p(offset) - prm%T_ref)**1 / 1. * prm%A(1:3,1:3,1) + & ! constant coefficient
(temperature(homog)%p(offset) - prm%T_ref)**2 / 2. * prm%expansion(1:3,1:3,2) + & ! linear coefficient (temperature(homog)%p(offset) - prm%T_ref)**2 / 2. * prm%A(1:3,1:3,2) + & ! linear coefficient
(temperature(homog)%p(offset) - prm%T_ref)**3 / 3. * prm%expansion(1:3,1:3,3) ! quadratic coefficient (temperature(homog)%p(offset) - prm%T_ref)**3 / 3. * prm%A(1:3,1:3,3) ! quadratic coefficient
end associate end associate
end function kinematics_thermal_expansion_initialStrain end function kinematics_thermal_expansion_initialStrain
@ -133,14 +133,14 @@ module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, i
associate(prm => param(kinematics_thermal_expansion_instance(phase))) associate(prm => param(kinematics_thermal_expansion_instance(phase)))
Li = TDot * ( & Li = TDot * ( &
prm%expansion(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient
+ prm%expansion(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient + prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient
+ prm%expansion(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient
) / & ) / &
(1.0_pReal & (1.0_pReal &
+ prm%expansion(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. & + prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. &
+ prm%expansion(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. & + prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. &
+ prm%expansion(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. & + prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. &
) )
end associate end associate
dLi_dTstar = 0.0_pReal dLi_dTstar = 0.0_pReal

40
src/lattice.f90
View File

@ -394,13 +394,13 @@ module lattice
! SHOULD NOT BE PART OF LATTICE BEGIN ! SHOULD NOT BE PART OF LATTICE BEGIN
real(pReal), dimension(:), allocatable, public, protected :: & real(pReal), dimension(:), allocatable, public, protected :: &
lattice_mu, lattice_nu, & lattice_mu, lattice_nu, &
lattice_damageMobility, & lattice_M, &
lattice_massDensity, & lattice_rho, &
lattice_specificHeat lattice_c_p
real(pReal), dimension(:,:,:), allocatable, public, protected :: & real(pReal), dimension(:,:,:), allocatable, public, protected :: &
lattice_C66, & lattice_C66, &
lattice_thermalConductivity, & lattice_K, &
lattice_damageDiffusion lattice_D
integer(kind(lattice_UNDEFINED_ID)), dimension(:), allocatable, public, protected :: & integer(kind(lattice_UNDEFINED_ID)), dimension(:), allocatable, public, protected :: &
lattice_structure lattice_structure
! SHOULD NOT BE PART OF LATTICE END ! SHOULD NOT BE PART OF LATTICE END
@ -465,11 +465,11 @@ subroutine lattice_init
allocate(lattice_structure(Nphases),source = lattice_UNDEFINED_ID) allocate(lattice_structure(Nphases),source = lattice_UNDEFINED_ID)
allocate(lattice_C66(6,6,Nphases), source=0.0_pReal) allocate(lattice_C66(6,6,Nphases), source=0.0_pReal)
allocate(lattice_thermalConductivity (3,3,Nphases), source=0.0_pReal) allocate(lattice_K (3,3,Nphases), source=0.0_pReal)
allocate(lattice_damageDiffusion (3,3,Nphases), source=0.0_pReal) allocate(lattice_D (3,3,Nphases), source=0.0_pReal)
allocate(lattice_damageMobility,& allocate(lattice_M,&
lattice_massDensity,lattice_specificHeat, & lattice_rho,lattice_c_p, &
lattice_mu, lattice_nu,& lattice_mu, lattice_nu,&
source=[(0.0_pReal,i=1,Nphases)]) source=[(0.0_pReal,i=1,Nphases)])
@ -517,22 +517,22 @@ subroutine lattice_init
! SHOULD NOT BE PART OF LATTICE BEGIN ! SHOULD NOT BE PART OF LATTICE BEGIN
lattice_thermalConductivity(1,1,p) = phase%get_asFloat('K_11',defaultVal=0.0_pReal) lattice_K(1,1,p) = phase%get_asFloat('K_11',defaultVal=0.0_pReal)
lattice_thermalConductivity(2,2,p) = phase%get_asFloat('K_22',defaultVal=0.0_pReal) lattice_K(2,2,p) = phase%get_asFloat('K_22',defaultVal=0.0_pReal)
lattice_thermalConductivity(3,3,p) = phase%get_asFloat('K_33',defaultVal=0.0_pReal) lattice_K(3,3,p) = phase%get_asFloat('K_33',defaultVal=0.0_pReal)
lattice_thermalConductivity(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_thermalConductivity(1:3,1:3,p), & lattice_K(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_K(1:3,1:3,p), &
phase%get_asString('lattice')) phase%get_asString('lattice'))
lattice_specificHeat(p) = phase%get_asFloat('c_p',defaultVal=0.0_pReal) lattice_c_p(p) = phase%get_asFloat('c_p', defaultVal=0.0_pReal)
lattice_massDensity(p) = phase%get_asFloat('rho', defaultVal=0.0_pReal) lattice_rho(p) = phase%get_asFloat('rho', defaultVal=0.0_pReal)
lattice_DamageDiffusion(1,1,p) = phase%get_asFloat('D_11',defaultVal=0.0_pReal) lattice_D(1,1,p) = phase%get_asFloat('D_11',defaultVal=0.0_pReal)
lattice_DamageDiffusion(2,2,p) = phase%get_asFloat('D_22',defaultVal=0.0_pReal) lattice_D(2,2,p) = phase%get_asFloat('D_22',defaultVal=0.0_pReal)
lattice_DamageDiffusion(3,3,p) = phase%get_asFloat('D_33',defaultVal=0.0_pReal) lattice_D(3,3,p) = phase%get_asFloat('D_33',defaultVal=0.0_pReal)
lattice_DamageDiffusion(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_DamageDiffusion(1:3,1:3,p), & lattice_D(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,p), &
phase%get_asString('lattice')) phase%get_asString('lattice'))
lattice_DamageMobility(p) = phase%get_asFloat('M',defaultVal=0.0_pReal) lattice_M(p) = phase%get_asFloat('M',defaultVal=0.0_pReal)
! SHOULD NOT BE PART OF LATTICE END ! SHOULD NOT BE PART OF LATTICE END
call selfTest call selfTest

38
src/source_damage_anisoBrittle.f90
View File

@ -12,11 +12,11 @@ submodule (constitutive:constitutive_damage) source_damage_anisoBrittle
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: & real(pReal) :: &
sdot_0, & !< opening rate of cleavage planes dot_o, & !< opening rate of cleavage planes
n !< damage rate sensitivity q !< damage rate sensitivity
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
critDisp, & !< critical displacement s_crit, & !< critical displacement
critLoad !< critical load g_crit !< critical load
real(pReal), dimension(:,:,:,:), allocatable :: & real(pReal), dimension(:,:,:,:), allocatable :: &
cleavage_systems cleavage_systems
integer :: & integer :: &
@ -75,18 +75,18 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
N_cl = src%get_asInts('N_cl',defaultVal=emptyIntArray) N_cl = src%get_asInts('N_cl',defaultVal=emptyIntArray)
prm%sum_N_cl = sum(abs(N_cl)) prm%sum_N_cl = sum(abs(N_cl))
prm%n = src%get_asFloat('q') prm%q = src%get_asFloat('q')
prm%sdot_0 = src%get_asFloat('dot_o') prm%dot_o = src%get_asFloat('dot_o')
prm%critDisp = src%get_asFloats('s_crit', requiredSize=size(N_cl)) prm%s_crit = src%get_asFloats('s_crit', requiredSize=size(N_cl))
prm%critLoad = src%get_asFloats('g_crit', requiredSize=size(N_cl)) prm%g_crit = src%get_asFloats('g_crit', requiredSize=size(N_cl))
prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),& prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal)) phase%get_asFloat('c/a',defaultVal=0.0_pReal))
! expand: family => system ! expand: family => system
prm%critDisp = math_expand(prm%critDisp,N_cl) prm%s_crit = math_expand(prm%s_crit,N_cl)
prm%critLoad = math_expand(prm%critLoad,N_cl) prm%g_crit = math_expand(prm%g_crit,N_cl)
#if defined (__GFORTRAN__) #if defined (__GFORTRAN__)
prm%output = output_asStrings(src) prm%output = output_asStrings(src)
@ -95,10 +95,10 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
#endif #endif
! sanity checks ! sanity checks
if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%sdot_0 <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o' if (prm%dot_o <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o'
if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit' if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit'
if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit' if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
@ -152,14 +152,14 @@ module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i))
traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i))
traction_crit = prm%critLoad(i)*damage(homog)%p(damageOffset)**2.0_pReal traction_crit = prm%g_crit(i)*damage(homog)%p(damageOffset)**2.0_pReal
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
= sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & = sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
+ prm%sdot_0 / prm%critDisp(i) & + prm%dot_o / prm%s_crit(i) &
* ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%n + & * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + &
(max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%n + & (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + &
(max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%n) (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%q)
enddo enddo
end associate end associate

16
src/source_damage_anisoDuctile.f90
View File

@ -12,9 +12,9 @@ submodule(constitutive:constitutive_damage) source_damage_anisoDuctile
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: & real(pReal) :: &
n !< damage rate sensitivity q !< damage rate sensitivity
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
critPlasticStrain !< critical plastic strain per slip system gamma_crit !< critical plastic strain per slip system
character(len=pStringLen), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
output output
end type tParameters end type tParameters
@ -68,11 +68,11 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
src => sources%get(sourceOffset) src => sources%get(sourceOffset)
N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray)
prm%n = src%get_asFloat('q') prm%q = src%get_asFloat('q')
prm%critPlasticStrain = src%get_asFloats('gamma_crit',requiredSize=size(N_sl)) prm%gamma_crit = src%get_asFloats('gamma_crit',requiredSize=size(N_sl))
! expand: family => system ! expand: family => system
prm%critPlasticStrain = math_expand(prm%critPlasticStrain,N_sl) prm%gamma_crit = math_expand(prm%gamma_crit,N_sl)
#if defined (__GFORTRAN__) #if defined (__GFORTRAN__)
prm%output = output_asStrings(src) prm%output = output_asStrings(src)
@ -81,8 +81,8 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
#endif #endif
! sanity checks ! sanity checks
if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
NipcMyPhase=count(material_phaseAt==p) * discretization_nIP NipcMyPhase=count(material_phaseAt==p) * discretization_nIP
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
@ -127,7 +127,7 @@ module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
associate(prm => param(source_damage_anisoDuctile_instance(phase))) associate(prm => param(source_damage_anisoDuctile_instance(phase)))
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
= sum(plasticState(phase)%slipRate(:,constituent)/(damage(homog)%p(damageOffset)**prm%n)/prm%critPlasticStrain) = sum(plasticState(phase)%slipRate(:,constituent)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit)
end associate end associate
end subroutine source_damage_anisoDuctile_dotState end subroutine source_damage_anisoDuctile_dotState

14
src/source_damage_isoDuctile.f90
View File

@ -12,8 +12,8 @@ submodule (constitutive:constitutive_damage) source_damage_isoDuctile
type:: tParameters !< container type for internal constitutive parameters type:: tParameters !< container type for internal constitutive parameters
real(pReal) :: & real(pReal) :: &
critPlasticStrain, & !< critical plastic strain gamma_crit, & !< critical plastic strain
N q
character(len=pStringLen), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
output output
end type tParameters end type tParameters
@ -64,8 +64,8 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
associate(prm => param(source_damage_isoDuctile_instance(p))) associate(prm => param(source_damage_isoDuctile_instance(p)))
src => sources%get(sourceOffset) src => sources%get(sourceOffset)
prm%N = src%get_asFloat('q') prm%q = src%get_asFloat('q')
prm%critPlasticStrain = src%get_asFloat('gamma_crit') prm%gamma_crit = src%get_asFloat('gamma_crit')
#if defined (__GFORTRAN__) #if defined (__GFORTRAN__)
prm%output = output_asStrings(src) prm%output = output_asStrings(src)
@ -74,8 +74,8 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
#endif #endif
! sanity checks ! sanity checks
if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%critPlasticStrain <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
NipcMyPhase=count(material_phaseAt==p) * discretization_nIP NipcMyPhase=count(material_phaseAt==p) * discretization_nIP
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
@ -120,7 +120,7 @@ module subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
associate(prm => param(source_damage_isoDuctile_instance(phase))) associate(prm => param(source_damage_isoDuctile_instance(phase)))
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%N)/prm%critPlasticStrain sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit
end associate end associate
end subroutine source_damage_isoDuctile_dotState end subroutine source_damage_isoDuctile_dotState

18
src/source_thermal_externalheat.f90
View File

@ -13,8 +13,8 @@ submodule(constitutive:constitutive_thermal) source_thermal_externalheat
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
time, & t_n, &
heat_rate f_T
integer :: & integer :: &
nIntervals nIntervals
end type tParameters end type tParameters
@ -64,10 +64,10 @@ module function source_thermal_externalheat_init(source_length) result(mySources
associate(prm => param(source_thermal_externalheat_instance(p))) associate(prm => param(source_thermal_externalheat_instance(p)))
src => sources%get(sourceOffset) src => sources%get(sourceOffset)
prm%time = src%get_asFloats('t_n') prm%t_n = src%get_asFloats('t_n')
prm%nIntervals = size(prm%time) - 1 prm%nIntervals = size(prm%t_n) - 1
prm%heat_rate = src%get_asFloats('f_T',requiredSize = size(prm%time)) prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n))
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP NipcMyPhase = count(material_phaseAt==p) * discretization_nIP
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0)
@ -121,13 +121,13 @@ module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_d
associate(prm => param(source_thermal_externalheat_instance(phase))) associate(prm => param(source_thermal_externalheat_instance(phase)))
do interval = 1, prm%nIntervals ! scan through all rate segments do interval = 1, prm%nIntervals ! scan through all rate segments
frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - prm%time(interval)) & frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) &
/ (prm%time(interval+1) - prm%time(interval)) ! fractional time within segment / (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment
if ( (frac_time < 0.0_pReal .and. interval == 1) & if ( (frac_time < 0.0_pReal .and. interval == 1) &
.or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) & .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) &
.or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) & .or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) &
TDot = prm%heat_rate(interval ) * (1.0_pReal - frac_time) + & TDot = prm%f_T(interval ) * (1.0_pReal - frac_time) + &
prm%heat_rate(interval+1) * frac_time ! interpolate heat rate between segment boundaries... prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries...
! ...or extrapolate if outside of bounds ! ...or extrapolate if outside of bounds
enddo enddo
dTDot_dT = 0.0 dTDot_dT = 0.0

4
src/thermal_adiabatic.f90
View File

@ -169,7 +169,7 @@ function thermal_adiabatic_getSpecificHeat(ip,el)
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat & thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat &
+ lattice_specificHeat(material_phaseAt(grain,el)) + lattice_c_p(material_phaseAt(grain,el))
enddo enddo
thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat & thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat &
@ -195,7 +195,7 @@ function thermal_adiabatic_getMassDensity(ip,el)
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity & thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity &
+ lattice_massDensity(material_phaseAt(grain,el)) + lattice_rho(material_phaseAt(grain,el))
enddo enddo
thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity & thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity &

6
src/thermal_conduction.f90
View File

@ -127,7 +127,7 @@ function thermal_conduction_getConductivity(ip,el)
thermal_conduction_getConductivity = 0.0_pReal thermal_conduction_getConductivity = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
thermal_conduction_getConductivity = thermal_conduction_getConductivity + & thermal_conduction_getConductivity = thermal_conduction_getConductivity + &
crystallite_push33ToRef(grain,ip,el,lattice_thermalConductivity(:,:,material_phaseAt(grain,el))) crystallite_push33ToRef(grain,ip,el,lattice_K(:,:,material_phaseAt(grain,el)))
enddo enddo
thermal_conduction_getConductivity = thermal_conduction_getConductivity & thermal_conduction_getConductivity = thermal_conduction_getConductivity &
@ -153,7 +153,7 @@ function thermal_conduction_getSpecificHeat(ip,el)
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
+ lattice_specificHeat(material_phaseAt(grain,el)) + lattice_c_p(material_phaseAt(grain,el))
enddo enddo
thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
@ -180,7 +180,7 @@ function thermal_conduction_getMassDensity(ip,el)
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &
+ lattice_massDensity(material_phaseAt(grain,el)) + lattice_rho(material_phaseAt(grain,el))
enddo enddo
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &