matscipy.fracture_mechanics.crack

Functions

G_to_strain(G, E, nu, orig_height)

Convert from energy release rate G to strain for thin strip geometry

counted(f)

deformation_gradient_residual(r0, crack, x, ...)

deformation_gradient_residuals(r0, crack, x, ...)

Return actual displacement field minus ideal displacement field.

displacement_residual(r0, crack, x, y, ...)

displacement_residuals(r0, crack, x, y, ...)

Return actual displacement field minus ideal displacement field, divided by r**alpha.

find_tip_broken_bonds(atoms, cutoff[, ...])

Find position of the tip from the atom coordination, i.e. broken bonds.

find_tip_coordination(a[, bondlength, bulk_nn])

Find position of tip in crack cluster from coordination

find_tip_stress_field(atoms[, r_range, ...])

Find the position of crack tip by fitting to the isotropic K-field stress

fit_crack_stress_field(atoms[, r_range, ...])

Perform a least squares fit of near-tip stress field to isotropic solution

get_energy_release_rate(atoms)

Return the current energy release rate G for atoms

get_strain(atoms)

Return the current strain on thin strip configuration atoms

get_stress_intensity_factor(atoms[, ...])

Compute stress intensity factor K_I

isotropic_modeII_crack_tip_displacement_field(K, ...)

Compute Irwin singular crack tip displacement field for mode II fracture.

isotropic_modeI_crack_tip_displacement_field(K, ...)

Compute Irwin singular crack tip displacement field for mode I fracture.

isotropic_modeI_crack_tip_stress_field(K, r, t)

Compute Irwin singular crack tip stress field for mode I fracture.

plot_stress_fields(atoms[, r_range, ...])

Fit and plot atomistic and continuum stress fields

print_crack_system(directions)

Pretty printing of crack crystallographic coordinate system

strain_to_G(strain, E, nu, orig_height)

Convert from strain to energy release rate G for thin strip geometry

thin_strip_displacement_y(x, y, strain, a, b)

Return vertical displacement ramp used to apply initial strain to slab

Classes

ConstantStrainRate(orig_height, delta_strain)

Constraint which increments epsilon_yy at a constant strain rate

CubicCrystalCrack(crack_surface, crack_front)

Crack in a cubic crystal.

IsotropicStressField([K, x0, y0, sxx0, ...])

Calculator to return Irwin near-tip stress field at atomic sites

RectilinearAnisotropicCrack()

Near field solution for a crack in a rectilinear anisotropic elastic medium.

SinclairCrack(crk, cryst, calc, k[, kI, ...])

Flexible boundary conditions for a Mode I crack as described in
matscipy.fracture_mechanics.crack.counted(f)
class matscipy.fracture_mechanics.crack.RectilinearAnisotropicCrack

Bases: object

Near field solution for a crack in a rectilinear anisotropic elastic medium. See: G. C. Sih, P. C. Paris and G. R. Irwin, Int. J. Frac. Mech. 1, 189 (1965)

Methods

deformation_gradient(r, theta, kI[, kII])

Deformation gradient tensor in mode I/II fracture.

displacements(r, theta, kI[, kII])

Displacement field in mode I/II fracture.

k1g(surface_energy)

K1G, Griffith critical stress intensity in mode I fracture

rtheta(u, v, k)

Invert displacement field in mode I fracture, i.e. compute r and theta from displacements.

stresses(r, theta, kI[, kII])

Stress field in mode I/II fracture.

k1gsqG

set_plane_strain

set_plane_stress
__init__()
set_plane_stress(a11, a22, a12, a16, a26, a66)
set_plane_strain(b11, b22, b33, b12, b13, b23, b16, b26, b36, b66)
displacements(r, theta, kI, kII=0.0)

Displacement field in mode I/II fracture. Positions are passed in cylinder coordinates.

Parameters:

  • r (array_like) – Distances from the crack tip.

  • theta (array_like) – Angles with respect to the plane of the crack.

  • k (float) – Stress intensity factor.

Returns:

  • u (array) – Displacements parallel to the plane of the crack.

  • v (array) – Displacements normal to the plane of the crack.

deformation_gradient(r, theta, kI, kII=0.0)

Deformation gradient tensor in mode I/II fracture. Positions are passed in cyclinder coordinates.

Parameters:

  • r (array_like) – Distances from the crack tip.

  • theta (array_like) – Angles with respect to the plane of the crack.

  • k (float) – Stress intensity factor.

Returns:

  • du_dx (array) – Derivatives of displacements parallel to the plane within the plane.

  • du_dy (array) – Derivatives of displacements parallel to the plane perpendicular to the plane.

  • dv_dx (array) – Derivatives of displacements normal to the plane of the crack within the plane.

  • dv_dy (array) – Derivatives of displacements normal to the plane of the crack perpendicular to the plane.

stresses(r, theta, kI, kII=0.0)

Stress field in mode I/II fracture. Positions are passed in cylinder coordinates.

Parameters:

  • r (array_like) – Distances from the crack tip.

  • theta (array_like) – Angles with respect to the plane of the crack.

  • k (float) – Stress intensity factor.

Returns:

  • sig_x (array) – Diagonal component of stress tensor parallel to the plane of the crack.

  • sig_y (array) – Diagonal component of stress tensor normal to the plane of the crack.

  • sig_xy (array) – Off-diagonal component of the stress tensor.

rtheta(u, v, k)

Invert displacement field in mode I fracture, i.e. compute r and theta from displacements.

k1g(surface_energy)

K1G, Griffith critical stress intensity in mode I fracture

k1gsqG()
matscipy.fracture_mechanics.crack.displacement_residuals(r0, crack, x, y, ref_x, ref_y, kI, kII=0, power=1)

Return actual displacement field minus ideal displacement field, divided by r**alpha.

matscipy.fracture_mechanics.crack.displacement_residual(r0, crack, x, y, ref_x, ref_y, kI, kII=0, mask=None, power=1)
matscipy.fracture_mechanics.crack.deformation_gradient_residuals(r0, crack, x, y, cur, ref_x, ref_y, ref, kI, cutoff, kII=0)

Return actual displacement field minus ideal displacement field.

matscipy.fracture_mechanics.crack.deformation_gradient_residual(r0, crack, x, y, cur, ref_x, ref_y, ref, kI, cutoff, kII=0, mask=None)
class matscipy.fracture_mechanics.crack.CubicCrystalCrack(crack_surface, crack_front, C11=None, C12=None, C44=None, stress_state='plane strain', C=None, Crot=None, cauchy_born=None)

Bases: object

Crack in a cubic crystal.

Methods

crack_tip_position(x, y, ref_x, ref_y, x0, ...)

Return an estimate of the real crack tip position by minimizing the mean square error of the current configuration relative to the displacement field obtained for a stress intensity factor k from linear elastic fracture mechanics.

crack_tip_position_y(x, y, ref_x, ref_y, x0, ...)

Return an estimate of the y-coordinate of the real crack tip position assuming the stress intensity factor is k.

deformation_gradient(ref_x, ref_y, x0, y0, kI)

Deformation gradient for a list of cartesian positions.

deformation_gradient_from_cartesian_coordinates(dx, ...)

Displacement field in mode I fracture from cartesian coordinates.

deformation_gradient_from_cylinder_coordinates(r, ...)

Displacement field in mode I fracture from cylindrical coordinates.

displacements(ref_x, ref_y, x0, y0, kI[, kII])

Displacement field for a list of cartesian positions.

displacements_from_cartesian_coordinates(dx, ...)

Displacement field in mode I/II fracture from cartesian coordinates.

displacements_from_cylinder_coordinates(r, ...)

Displacement field in mode I/II fracture from cylindrical coordinates.

k1g(surface_energy)

Compute Griffith critical stress intensity in mode I fracture.

scale_displacements(x, y, ref_x, ref_y, ...)

Rescale atomic positions from stress intensity factor old_k to the new stress intensity factor new_k.

stresses(ref_x, ref_y, x0, y0, kI[, kII])

Stress field for a list of cartesian positions

stresses_from_cartesian_coordinates(dx, dy, kI)

Stress field in mode I fracture from cartesian coordinates

stresses_from_cylinder_coordinates(r, theta, kI)

Stress field in mode I fracture from cylindrical coordinates

k1gsqG
__init__(crack_surface, crack_front, C11=None, C12=None, C44=None, stress_state='plane strain', C=None, Crot=None, cauchy_born=None)

Initialize a crack in a cubic crystal with elastic constants C11, C12 and C44 (or optionally a full 6x6 elastic constant matrix C). The crack surface is given by crack_surface, the cracks runs in the plane given by crack_front.

k1g(surface_energy)

Compute Griffith critical stress intensity in mode I fracture.

Parameters:

surface_energy (float) – Surface energy of the respective crystal surface.

Returns:

k1g – Stress intensity factor.

Return type:

float

k1gsqG()
displacements_from_cylinder_coordinates(r, theta, kI, kII=0)

Displacement field in mode I/II fracture from cylindrical coordinates.

displacements_from_cartesian_coordinates(dx, dy, kI, kII=0)

Displacement field in mode I/II fracture from cartesian coordinates.

displacements(ref_x, ref_y, x0, y0, kI, kII=0)

Displacement field for a list of cartesian positions.

Parameters:

  • ref_x (array_like) – x-positions of the reference crystal.

  • ref_y (array_like) – y-positions of the reference crystal.

  • x0 (float) – x-coordinate of the crack tip.

  • y0 (float) – y-coordinate of the crack tip.

  • k (float) – Stress intensity factor.

Returns:

  • ux (array_like) – x-displacements.

  • uy (array_like) – y-displacements.

deformation_gradient_from_cylinder_coordinates(r, theta, kI, kII=0)

Displacement field in mode I fracture from cylindrical coordinates.

deformation_gradient_from_cartesian_coordinates(dx, dy, kI, kII=0)

Displacement field in mode I fracture from cartesian coordinates.

deformation_gradient(ref_x, ref_y, x0, y0, kI, kII=0)

Deformation gradient for a list of cartesian positions.

Parameters:

  • ref_x (array_like) – x-positions of the reference crystal.

  • ref_y (array_like) – y-positions of the reference crystal.

  • x0 (float) – x-coordinate of the crack tip.

  • y0 (float) – y-coordinate of the crack tip.

  • k (float) – Stress intensity factor.

Returns:

  • ux (array_like) – x-displacements.

  • uy (array_like) – y-displacements.

crack_tip_position(x, y, ref_x, ref_y, x0, y0, kI, kII=0, mask=None, residual_func=<function displacement_residual>, method='Powell', return_residuals=False)

Return an estimate of the real crack tip position by minimizing the mean square error of the current configuration relative to the displacement field obtained for a stress intensity factor k from linear elastic fracture mechanics.

Parameters:

  • x (array_like) – x-positions of the atomic system containing the crack.

  • y (array_like) – y-positions of the atomic system containing the crack.

  • ref_x (array_like) – x-positions of the reference crystal.

  • ref_y (array_like) – y-positions of the reference crystal.

  • x0 (float) – Initial guess for the x-coordinate of the crack tip.

  • y0 (float) – Initial guess for the y-coordinate of the crack tip.

  • k (float) – Stress intensity factor.

  • mask (array_like, optional) – Marks the atoms to use for this calculation.

  • residual_func (function) – Function returning per-atom residuals to be minimized.

  • method (str) – Optimization method. See method argument of scipy.optimize.minimize. Additionally, ‘leastsq’ invokes scipy.optimize.leastsq.

  • return_residuals (bool) – Function returns residuals if set to True.

Returns:

  • x0 (float) – x-coordinate of the crack tip.

  • y0 (float) – y-coordinate of the crack tip.

  • residuals (array, optional) – Per-atom residuals at end of optimization.

crack_tip_position_y(x, y, ref_x, ref_y, x0, y0, kI, kII=0, mask=None)

Return an estimate of the y-coordinate of the real crack tip position assuming the stress intensity factor is k.

Parameters:

  • x (array_like) – x-positions of the atomic system containing the crack.

  • y (array_like) – y-positions of the atomic system containing the crack.

  • ref_x (array_like) – x-positions of the reference crystal.

  • ref_y (array_like) – y-positions of the reference crystal.

  • x0 (float) – Initial guess for the x-coordinate of the crack tip.

  • y0 (float) – Initial guess for the y-coordinate of the crack tip.

  • k (float) – Stress intensity factor.

  • mask (array_like, optional) – Marks the atoms to use for this calculation.

Returns:

y0 – y-coordinate of the crack tip

Return type:

float

scale_displacements(x, y, ref_x, ref_y, old_k, new_k)

Rescale atomic positions from stress intensity factor old_k to the new stress intensity factor new_k. This is useful for extrapolation of relaxed positions. Note - this only works for pure mode I fracture

Parameters:

  • x (array_like) – x-positions of the atomic system containing the crack.

  • y (array_like) – y-positions of the atomic system containing the crack.

  • ref_x (array_like) – x-positions of the reference crystal.

  • ref_y (array_like) – y-positions of the reference crystal.

  • old_k (float) – Current stress intensity factor.

  • new_k (float) – Stress intensity factor for the output positions.

Returns:

  • x (array_like) – New x-positions of the atomic system containing the crack.

  • y (array_like) – New y-positions of the atomic system containing the crack.

stresses_from_cylinder_coordinates(r, theta, kI, kII=0)

Stress field in mode I fracture from cylindrical coordinates

stresses_from_cartesian_coordinates(dx, dy, kI, kII=0)

Stress field in mode I fracture from cartesian coordinates

stresses(ref_x, ref_y, x0, y0, kI, kII=0)

Stress field for a list of cartesian positions

Parameters:

  • ref_x (array_like) – x-positions of the reference crystal.

  • ref_y (array_like) – y-positions of the reference crystal.

  • x0 (float) – x-coordinate of the crack tip.

  • y0 (float) – y-coordinate of the crack tip.

  • k (float) – Stress intensity factor.

Returns:

  • sig_x (array_like) – xx-component of the stress tensor

  • sig_y (array_like) – yy-component of the stress tensor

  • sig_xy (array_like) – xy-component of the stress tensor

class matscipy.fracture_mechanics.crack.SinclairCrack(crk, cryst, calc, k, kI=None, kII=0, alpha=0.0, vacuum=6.0, variable_alpha=True, variable_k=False, alpha_scale=None, k_scale=None, extended_far_field=False, rI=0.0, rIII=0.0, cutoff=0.0, incl_rI_f_alpha=False, is_3D=False, theta_3D=0.0, cont_k='k1', precon_recompute_interval=10)

Bases: object

Flexible boundary conditions for a Mode I crack as described in

Sinclair, J. E. The Influence of the Interatomic Force Law and of Kinks on the

Propagation of Brittle Cracks. Philos. Mag. 31, 647–671 (1975)

Attributes:
k

Methods

strain_err(cutoff[, seperate_surface])

Function that returns the atomistic corrector strain error Dv for each atom using the norm of the difference between the corrector on each atom and all those around it within some cutoff.

u_cle([alpha, kI, kII])

Returns CLE displacement solution at current (alpha, K)

update_atoms([use_alpha_3D])

Update self.atoms from degrees of freedom (self.u, self.alpha, self.kI, self.kII)

animate

arc_length_continuation

convergence_line_plot

dfk_dk_approx

fit_cle

get_crack_tip_force

get_deformation_gradient

get_dofs

get_forces

get_k_force

get_potential_energy

get_precon

get_xdot

load_cb_model

optimize

pack

plot

rescale_k

save_cb_model

set_atoms

set_dofs

set_shiftmask

unpack

update_precon

write_atoms_to_file
__init__(crk, cryst, calc, k, kI=None, kII=0, alpha=0.0, vacuum=6.0, variable_alpha=True, variable_k=False, alpha_scale=None, k_scale=None, extended_far_field=False, rI=0.0, rIII=0.0, cutoff=0.0, incl_rI_f_alpha=False, is_3D=False, theta_3D=0.0, cont_k='k1', precon_recompute_interval=10)
Parameters:

  • CubicCrystalCrack (crk - instance of) –

  • crystal (cryst - Atoms object representing undeformed) –

  • calculator (calc - ASE-compatible) –

  • kI) (k - stress intensity factor in mode I (deprecated name for) –

  • I (rI - Size of region) –

  • II (kII - stress intensity factor in mode) –

  • position (alpha - crack tip) –

  • cell (vacuum - amount of vacuum to add to unit) –

  • True (extended_far_field - if) –

  • BCs) (crack tip position can vary (flexible) –

  • True – (needed for arc-length continuation)

  • vary (stress intensity factor can) – (needed for arc-length continuation)

  • True –

  • contrib (crack tip force includes region III) –

  • I –

  • III (rIII - Size of region) –

  • cutoff (cutoff - Potential) –

  • forces (incl_rI_f_alpha - Whether to include region I when calculating the f_alpha) –

  • preconditioner (precon_recompute_interval - How many preconditioned steps between re-computing the) –

property k
pack(u, alpha, k)
unpack(x, reshape=False, defaults=None)
get_dofs()
set_dofs(x)
u_cle(alpha=None, kI=None, kII=None)

Returns CLE displacement solution at current (alpha, K)

Note that this is now pre-multipled by the stress intensity factor k

fit_cle(r_fit=20.0, variable_alpha=True, variable_k=True, x0=None, grid=None)
get_deformation_gradient(x, y, kI=None, kII=None, de=0)
set_shiftmask(radial_dist)
update_atoms(use_alpha_3D=False)

Update self.atoms from degrees of freedom (self.u, self.alpha, self.kI, self.kII)

save_cb_model()
load_cb_model()
set_atoms(atoms)
get_crack_tip_force(forces=None, mask=None, full_array_output=False)
get_xdot(x1, x2, ds=None)
get_k_force(x1, xdot1, ds)
get_forces(x1=None, xdot1=None, ds=None, forces=None, mask=None)
dfk_dk_approx(xdot1)
update_precon(x, F=None)
get_precon(x, F)
optimize(ftol=0.001, steps=20, dump=False, args=None, precon=False, method='krylov', check_grad=True, dump_interval=10, verbose=0)
get_potential_energy()
rescale_k(new_kI)
arc_length_continuation(x0, x1, N=10, ds=0.01, ftol=0.01, direction=1, max_steps=10, continuation=False, traj_file='x_traj.h5', traj_interval=1, otf_traj=False, precon=False, ds_max=inf, ds_min=0, ds_aggressiveness=2, opt_method='krylov', cos_alpha_min=0.9, parallel=False, pipe_output=None, data_queue=None, kill_confirm_queue=None, allow_alpha_backtracking=False, follow_G_contour=False)
plot(ax=None, regions='1234', rzoom=inf, bonds=None, cutoff=2.8, tip=False, regions_styles=None, atoms_args=None, bonds_args=None, tip_args=None)
animate(x, k1g, regions='12', rzoom=inf, cutoff=2.8, frames=None, callback=None, plot_tip=True, regions_styles=None, atoms_args=None, bonds_args=None, tip_args=None)
write_atoms_to_file()
strain_err(cutoff, seperate_surface=False)

Function that returns the atomistic corrector strain error Dv for each atom using the norm of the difference between the corrector on each atom and all those around it within some cutoff. Also has an option to return states adjacent to the free surface of the crack seperately for comparison.

Parameters:

  • cutoff (float) – cutoff around each atom to use to find the strain error norm.

  • seperate_surface (bool) – whether or not to return the surface atoms in a seperate array.

Returns:

  • r (array of atom radial distance from the centre of the crack system)

  • dv (norm strain error for each atom in r.)

convergence_line_plot(num=0)
matscipy.fracture_mechanics.crack.isotropic_modeI_crack_tip_stress_field(K, r, t, xy_only=True, nu=0.5, stress_state='plane strain')

Compute Irwin singular crack tip stress field for mode I fracture.

Parameters:

  • K (float) – Mode I stress intensity factor. Units should match units of r.

  • r (array_like) – Radial distances from crack tip. Can be a multidimensional array to evaluate stress field on a grid.

  • t (array_like) – Angles from horzontal line y=0 ahead of crack tip, measured anticlockwise. Should have same shape as r.

  • xy_only (bool) – If True (default) only xx, yy, xy and yx components will be set.

  • nu (float) – Poisson ratio. Used only when xy_only=False, to determine zz stresses

  • stress_state (str) – One of”plane stress” or “plane strain”. Used if xyz_only=False to determine zz stresses.

Returns:

sigma

Return type:

array with shape r.shape + (3,3)

matscipy.fracture_mechanics.crack.isotropic_modeI_crack_tip_displacement_field(K, G, nu, r, t, stress_state='plane strain')

Compute Irwin singular crack tip displacement field for mode I fracture.

Parameters:

  • K (float) – Mode I stress intensity factor. Units should match units of G and r.

  • G (float) – Shear modulus. Units should match units of K and r.

  • nu (float) – Poisson ratio.

  • r (array_like) – Radial distances from crack tip. Can be a multidimensional array to evaluate stress field on a grid.

  • t (array_like) – Angles from horizontal line y=0 ahead of crack tip, measured anticlockwise. Should have same shape as r.

  • stress_state (str) – One of”plane stress” or “plane strain”. Used if xyz_only=False to determine zz stresses.

Returns:

  • u (array)

  • v (array) – Displacements. Same shape as r and t.

matscipy.fracture_mechanics.crack.isotropic_modeII_crack_tip_displacement_field(K, G, nu, r, t, stress_state='plane strain')

Compute Irwin singular crack tip displacement field for mode II fracture.

Parameters:

  • K (float) – Mode I stress intensity factor. Units should match units of G and r.

  • G (float) – Shear modulus. Units should match units of K and r.

  • nu (float) – Poisson ratio.

  • r (array_like) – Radial distances from crack tip. Can be a multidimensional array to evaluate stress field on a grid.

  • t (array_like) – Angles from horizontal line y=0 ahead of crack tip, measured anticlockwise. Should have same shape as r.

  • stress_state (str) – One of”plane stress” or “plane strain”. Used if xyz_only=False to determine zz stresses.

Returns:

  • u (array)

  • v (array) – Displacements. Same shape as r and t.

class matscipy.fracture_mechanics.crack.IsotropicStressField(K=None, x0=None, y0=None, sxx0=0.0, syy0=0.0, sxy0=0.0, nu=0.5, stress_state='plane strain')

Bases: object

Calculator to return Irwin near-tip stress field at atomic sites

Methods

get_stresses
__init__(K=None, x0=None, y0=None, sxx0=0.0, syy0=0.0, sxy0=0.0, nu=0.5, stress_state='plane strain')
get_stresses(atoms)
matscipy.fracture_mechanics.crack.strain_to_G(strain, E, nu, orig_height)

Convert from strain to energy release rate G for thin strip geometry

Parameters:

  • strain (float) – Dimensionless ratio (current_height - orig_height)/orig_height

  • E (float) – Young’s modulus relevant for a pull in y direction sigma_yy/eps_yy

  • nu (float) – Poission ratio -eps_yy/eps_xx

  • orig_height (float) – Unstrained height of slab

Returns:

G – Energy release rate in units consistent with input (i.e. in eV/A**2 if eV/A/fs units used)

Return type:

float

matscipy.fracture_mechanics.crack.G_to_strain(G, E, nu, orig_height)

Convert from energy release rate G to strain for thin strip geometry

Parameters:

  • G (float) – Energy release rate in units consistent with E and orig_height

  • E (float) – Young’s modulus relevant for a pull in y direction sigma_yy/eps_yy

  • nu (float) – Poission ratio -eps_yy/eps_xx

  • orig_height (float) – Unstrained height of slab

Returns:

strain – Dimensionless ratio (current_height - orig_height)/orig_height

Return type:

float

matscipy.fracture_mechanics.crack.get_strain(atoms)

Return the current strain on thin strip configuration atoms

Requires unstrained height of slab to be stored as OrigHeight key in atoms.info dictionary.

Also updates value stored in atoms.info.

matscipy.fracture_mechanics.crack.get_energy_release_rate(atoms)

Return the current energy release rate G for atoms

Result is computed assuming thin strip geometry, and using stored Young’s modulus and Poission ratio and original slab height from atoms.info dictionary.

Also updates G value stored in atoms.info dictionary.

matscipy.fracture_mechanics.crack.get_stress_intensity_factor(atoms, stress_state='plane strain')

Compute stress intensity factor K_I

Calls get_energy_release_rate() to compute G, then uses stored YoungsModulus and PoissionRatio_yz values from atoms.info dictionary to compute K_I.

Also updates value stored in atoms.info dictionary.

matscipy.fracture_mechanics.crack.fit_crack_stress_field(atoms, r_range=(0.0, 50.0), initial_params=None, fix_params=None, sigma=None, avg_sigma=None, avg_decay=0.005, calc=None, verbose=False)

Perform a least squares fit of near-tip stress field to isotropic solution

Stresses on the atoms are fit to the Irwin K-field singular crack tip solution, allowing the crack position, stress intensity factor and far-field stress components to vary during the fit.

Parameters:

  • atoms (Atoms object) –

    Crack system. For the initial fit, the following keys are used from the info dictionary:

    • YoungsModulus

    • PossionRatio_yx

    • G — current energy release rate

    • strain — current applied strain

    • CrackPos — initial guess for crack tip position

    The initial guesses for the stress intensity factor K and far-field stress sigma0 are computed from YoungsModulus, PoissonRatio_yx, G and strain, assuming plane strain in thin strip boundary conditions.

    On exit, new K, sigma0 and CrackPos entries are set in the info dictionary. These values are then used as starting guesses for subsequent fits.

  • r_range (sequence of two floats, optional) – If present, restrict the stress fit to an annular region r_range[0] <= r < r_range[1], centred on the previous crack position (from the CrackPos entry in atoms.info). If r_range is None, fit is carried out for all atoms.

  • initial_params (dict) – Names and initial values of parameters. Missing initial values are guessed from Atoms object.

  • fix_params (dict) – Names and values of parameters to fix during the fit, e.g. {y0: 0.0} to constrain the fit to the line y=0

  • sigma (None or array with shape (len(atoms), 3, 3)) – Explicitly provide the per-atom stresses. Avoids calling Atoms’ calculators get_stresses() method.

  • avg_sigma (None or array with shape (len(atoms), 3, 3)) – If present, use this array to accumulate the time-averaged stress field. Useful when processing a trajectory.

  • avg_decay (real) – Factor by which average stress is attenuated at each step. Should be set to dt/tau where dt is MD time-step and tau is a characteristic averaging time.

  • calc (Calculator object, optional) –

    If present, override the calculator used to compute stresses on the atoms. Default is atoms.get_calculator.

    To use the atom resolved stress tensor pass an instance of the AtomResolvedStressField class.

  • verbose (bool, optional) – If set to True, print additional information about the fit.

Returns:

params – Fitted parameters, in a form suitable for passing to IsotropicStressField constructor. These are the stress intensity factor K, the centre of the stress field (x0, y0), and the far field contribution to the stress (sxx0, syy0, sxy0).

Return type:

dict with keys [K, x0, y0, sxx0, syy0, sxy0]

matscipy.fracture_mechanics.crack.find_tip_coordination(a, bondlength=2.6, bulk_nn=4)

Find position of tip in crack cluster from coordination

matscipy.fracture_mechanics.crack.find_tip_broken_bonds(atoms, cutoff, bulk_nn=4, boundary_thickness=None)

Find position of the tip from the atom coordination, i.e. broken bonds. Using the C implementation of ‘neighbour_list’. Returns the tip’s position in cartesian coordinates.

Parameters:

  • atoms (ase.Atoms) – Atomic configuration.

  • cutoff (float) – Cutoff distance for neighbour search.

  • bulk_nn (integer) – Number of nearest neighbours for the standard bulk configuration.

  • boundary_buffer (float) – Thickness of the boundaries. Defaults to cutoff distance.

Returns:

tip_position – The x and y values are found. The z value is calculated as the midpoint of the depth.

Return type:

numpy array

matscipy.fracture_mechanics.crack.find_tip_stress_field(atoms, r_range=None, initial_params=None, fix_params=None, sigma=None, avg_sigma=None, avg_decay=0.005, calc=None)

Find the position of crack tip by fitting to the isotropic K-field stress

Fit is carried out using fit_crack_stress_field(), and parameters have the same meaning as there.

See also

fit_crack_stress_field

matscipy.fracture_mechanics.crack.plot_stress_fields(atoms, r_range=None, initial_params=None, fix_params=None, sigma=None, avg_sigma=None, avg_decay=0.005, calc=None)

Fit and plot atomistic and continuum stress fields

Firstly a fit to the Irwin K-field solution is carried out using fit_crack_stress_field(), and parameters have the same meaning as for that function. Then plots of the \(\sigma_{xx}\), \(\sigma_{yy}\), \(\sigma_{xy}\) fields are produced for atomistic and continuum cases, and for the residual error after fitting.

matscipy.fracture_mechanics.crack.thin_strip_displacement_y(x, y, strain, a, b)

Return vertical displacement ramp used to apply initial strain to slab

Strain is increased from 0 to strain over distance \(a <= x <= b\). Region \(x < a\) is rigidly shifted up/down by strain*height/2.

Here is an example of how to use this function on an artificial 2D square atomic lattice. The positions are plotted before (left) and after (right) applying the displacement, and the horizontal and vertical lines show the strain (red), a (green) and b (blue) parameters.

import matplotlib.pyplot as plt
import numpy as np

w = 1; h = 1; strain = 0.1; a = -0.5; b =  0.0
x = np.linspace(-w, w, 20)
y = np.linspace(-h, h, 20)
X, Y = np.meshgrid(x, y)
u_y = thin_strip_displacement_y(X, Y, strain, a, b)

for i, disp in enumerate([0, u_y]):
    plt.subplot(1,2,i+1)
    plt.scatter(X, Y + disp, c='k', s=5)
    for y in [-h, h]:
        plt.axhline(y, color='r', linewidth=2, linestyle='dashed')
        plt.axhline(y*(1+strain), color='r', linewidth=2)
    for x, c in zip([a, b], ['g', 'b']):
        plt.axvline(x, color=c, linewidth=2)
generated/thin-strip-displacement-y.png
Parameters:

  • x (array) –

  • y (array) – Atomic positions in unstrained slab, centered on origin x=0,y=0

  • strain (float) – Far field strain to apply

  • a (float) – x coordinate for beginning of strain ramp

  • b (float) – x coordinate for end of strain ramp

matscipy.fracture_mechanics.crack.print_crack_system(directions)

Pretty printing of crack crystallographic coordinate system

Specified by list of Miller indices for crack_direction (x), cleavage_plane (y) and crack_front (z), each of which should be a sequence of three floats

class matscipy.fracture_mechanics.crack.ConstantStrainRate(orig_height, delta_strain, mask=None)

Bases: object

Constraint which increments epsilon_yy at a constant strain rate

Rescaling is applied only to atoms where mask is True (default is all atoms)

Methods

apply_strain(atoms[, rigid_constraints])

Applies a constant strain to the system.

adjust_forces

adjust_positions

copy
__init__(orig_height, delta_strain, mask=None)
adjust_forces(atoms, forces)
adjust_positions(atoms, newpos)
copy()
apply_strain(atoms, rigid_constraints=False)

Applies a constant strain to the system.

Parameters:

  • atoms (ASE.atoms) – Atomic configuration.

  • rigid_constraints (boolean) – Apply (or not apply) strain to every atom. i.e. allow constrainted atoms to move during strain application