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

barnett_lothe_L(C[, n_omega])

Barnett-Lothe tensor L for a straight dislocation along z, from the 6x6 elastic constants C, via the integral formalism

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

hill_poisson_ratio(C)

Isotropic Poisson ratio of the Voigt-Reuss-Hill average of the 6x6 elastic constant matrix C (independent of the orientation of C).

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)
matscipy.fracture_mechanics.crack.hill_poisson_ratio(C)

Isotropic Poisson ratio of the Voigt-Reuss-Hill average of the 6x6 elastic constant matrix C (independent of the orientation of C).

matscipy.fracture_mechanics.crack.barnett_lothe_L(C, n_omega=64)

Barnett-Lothe tensor L for a straight dislocation along z, from the 6x6 elastic constants C, via the integral formalism

L = -(1/pi) int_0^pi [(mn)(nn)^-1(nm) - (mm)] domega,

with m = (cos w, sin w, 0), n = (-sin w, cos w, 0) and (ab)_ik = C_ijkl a_j b_l. Unlike the Stroh eigenvector route this is well conditioned for (nearly) isotropic materials, where the Stroh eigenvalues are degenerate. The integrand is smooth and pi-periodic, so the uniform rule converges exponentially (machine precision by n_omega ~ 64). For an isotropic material L = diag(mu/(1-nu), mu/(1-nu), mu).

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.

g1e(max_gamma)

Critical energy release rate for dislocation emission on a slip plane that coincides with the crack plane (theta = 0) under mode II loading along the slip direction, G_IIe = gamma_us (J.

g1e_fcc(surface_energy, max_gamma)

Critical energy release rate for dislocation emission in fcc crystals.

k1e_aniso(max_gamma, phi, theta)

K1e, Rice criterion for dislocation emission from crack tip in mode I fracture according to anisotropic elasticity Ref.

k1e_iso(max_gamma, phi, theta[, nu])

K1e, Rice criterion for dislocation emission from crack tip in mode I fracture.

k1g(surface_energy)

K1G, Griffith critical stress intensity in mode I fracture

k_from_g(G)

Mode I stress intensity factor corresponding to energy release rate G, using the same anisotropic relation as k1g(), so that k1g(surface_energy) == k_from_g(2 * surface_energy).

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()
g1e(max_gamma)

Critical energy release rate for dislocation emission on a slip plane that coincides with the crack plane (theta = 0) under mode II loading along the slip direction, G_IIe = gamma_us (J. R. Rice, JMPS, 1992).

For emission on an inclined slip plane under mode I loading use k1e_iso() or k1e_aniso().

Parameters:

max_gamma (float) – Unstable stacking fault energy of the slip system.

g1e_fcc(surface_energy, max_gamma)

Critical energy release rate for dislocation emission in fcc crystals. The critical energy release rate depends both on the unstable stacking fault energy and surface energy according to P. Andric & W. A. Curtin, JMPS, 2017.

Parameters:

  • surface_energy (surface energy of the crack plane)

  • max_gamma (the unstable stacking fault energy of the dislocation slip system.)

k_from_g(G)

Mode I stress intensity factor corresponding to energy release rate G, using the same anisotropic relation as k1g(), so that k1g(surface_energy) == k_from_g(2 * surface_energy).

k1e_iso(max_gamma, phi, theta, nu=None)

K1e, Rice criterion for dislocation emission from crack tip in mode I fracture. J. R. Rice, JMPS, (1992).

The critical energy release rate is Rice’s isotropic result

G_Ie = 8 gamma_us (1 + (1 - nu) tan^2 phi) / ((1 + cos theta) sin^2 theta)

converted to K with k_from_g(), so the result is in the same units as k1g() and the two can be compared directly.

Parameters:

  • max_gamma (float) – Maximum energy along the sliding path of the gamma line (unstable stacking fault energy).

  • phi (float (units: rad.)) – The angle between the slip direction and a vector lying on the slip plane and perpendicular to the crack-front direction.

  • theta (float (units: rad.)) – The angle between slip plane and crack plane.

  • nu (float, optional) – Poisson ratio. Defaults to the Voigt-Reuss-Hill average of the elastic constants C (exact for an isotropic material).

Returns:

  • K_ie_iso (Critical K_I for dislocation emission under mode-I loading calculated)

  • using isotropic elasticity.

k1e_aniso(max_gamma, phi, theta)

K1e, Rice criterion for dislocation emission from crack tip in mode I fracture according to anisotropic elasticity Ref. Beltz and Rice, JMPS, (1994).

Requires the 6x6 elastic constants in the crack frame, self.C (set by CubicCrystalCrack). The result is in the same units as k1g(); for an isotropic material it equals k1e_iso().

Parameters:

  • max_gamma (float) – Maximum energy along the sliding path of the gamma line (unstable stacking fault energy).

  • phi (float (units: rad.)) – The angle between the slip direction and a vector lying on the slip plane and perpendicular to the crack-front direction.

  • theta (float (units: rad.)) – The angle between slip plane and crack plane.

Returns:

  • K_ie_aniso (Critical K_I for dislocation emission under mode-I loading calculated)

  • using anisotropic elasticity.

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.

k1e_aniso(max_gamma, phi, theta)

Critical mode I stress intensity for dislocation emission, Beltz-Rice anisotropic criterion; same units as k1g().

k1e_iso(max_gamma, phi, theta[, nu])

Critical mode I stress intensity for dislocation emission, Rice's isotropic criterion; same units as k1g().

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

k1e_iso(max_gamma, phi, theta, nu=None)

Critical mode I stress intensity for dislocation emission, Rice’s isotropic criterion; same units as k1g(). See RectilinearAnisotropicCrack.k1e_iso().

k1e_aniso(max_gamma, phi, theta)

Critical mode I stress intensity for dislocation emission, Beltz-Rice anisotropic criterion; same units as k1g(). See RectilinearAnisotropicCrack.k1e_aniso().

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.

  • kI (float) – Mode I stress intensity factor.

  • kII (float) – Mode II stress intensity factor.

Returns:

F – (N, 2, 2) in-plane deformation gradient (identity included) in derivative-first order, [[du/dx, dv/dx], [du/dy, dv/dy]]; i.e. the transpose of the component-first F_ab = dx_a/dX_b.

Return type:

array_like

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