From 5012798576665071c1c865edd95ac1b9f68b59dd Mon Sep 17 00:00:00 2001 From: leiapple Date: Tue, 27 Feb 2024 14:20:49 +0100 Subject: [PATCH] Rice emit. --- matscipy/fracture_mechanics/crack.py | 151 ++++++++++++++++++++++++++- 1 file changed, 149 insertions(+), 2 deletions(-) diff --git a/matscipy/fracture_mechanics/crack.py b/matscipy/fracture_mechanics/crack.py index a60809c53..328a0cdb4 100644 --- a/matscipy/fracture_mechanics/crack.py +++ b/matscipy/fracture_mechanics/crack.py @@ -42,7 +42,8 @@ from matscipy.atomic_strain import atomic_strain from matscipy.elasticity import (rotate_elastic_constants, rotate_cubic_elastic_constants, - Voigt_6_to_full_3x3_stress) + Voigt_6_to_full_3x3_stress, + elastic_moduli) from matscipy.surface import MillerDirection, MillerPlane from matscipy.neighbours import neighbour_list from matscipy.optimize import ode12r @@ -341,8 +342,154 @@ def k1gsqG(self): return -2 / (self.a22 * ((self.mu1 + self.mu2) / (self.mu1 * self.mu2)).imag) -### +### Critical K_I for dislocation emission from (updated) Rice theory. + + def g1e(self, max_gamma): + """ + Critical energy release rate for dislocation emission. + J. R. Rice, JMPS, (1992). + The critical energy release rate is the unstable stacking + fault energy according to J. R. Rice. + + Parameters + ---------- + max_gamma: the unstable stacking fault energy of the dislocation slip system. + """ + + return max_gamma + + def g1e_fcc(self, 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 energyaccording to Ref. + (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. + """ + + if surface_energy > 3.45*max_gamma: + g1e_fcc = 0.145*surface_energy + 0.5*max_gamma + else: + g1e_fcc = max_gamma + + return g1e_fcc + + def k1e_iso(self, max_gamma, phi, theta): + """ + K1e, Rice criterion for dislocation emission from crack tip + in mode I fracture. J. R. Rice, JMPS, (1992). + The soluation is presented for isotropic elasticity. + + Parameters + ---------- + max_gamma: float + Maximum energy along the sliding path of the gamma line (maximum + stacking fault energy). + theta: float (units: rad.) + The angle between slip plane and crack plane. + 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. + + Returns + ------- + K_ie_iso : Critical K_I for dislocation emission under mode-I loading calculated + using isotropic elasticity. + """ + # calculate the isotropic elastic constant if not provided. + # mu: Shear modulus (GPa); B: bulk modulus (GPa); nu:Possion ratio + E, nu, Gm, B, K = elastic_moduli(self.C, l=np.array([1, 0, 0]), R=self.RotationMatrix, tol=1e-6) + E0 = E[0] + nu0 = nu[0,1] + G1e = 8 * max_gamma * ((1 + (1-nu0) * np.tan(phi)**2)/ + ((1 + np.cos(theta)) * np.sin(theta)**2)) + + K_ie_iso = np.sqrt(1.0e-3 * G1e / self.B) + + return K_ie_iso + + def k1e_aniso(self, 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). + + Parameters + ---------- + max_gamma: float + Maximum energy along the sliding path of the gamma line (maximum + stacking fault energy). + theta: float (units: rad.) + The angle between slip plane and crack plane. + 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. + + Returns + ------- + K_ie_aniso : Critical K_I for dislocation emission under mode-I loading calculated + using anisotropic elasticity. + """ + # assmeble the elastic matrix + Q = np.array([[self.C[0,0], self.C[0,5], self.C[0,4]], + [self.C[0,5], self.C[5,5], self.C[4,5]], + [self.C[0,4], self.C[4,5], self.C[4,4]]]) + + R = np.array([[self.C[0,5], self.C[0,1], self.C[0,3]], + [self.C[5,5], self.C[1,5], self.C[3,5]], + [self.C[4,5], self.C[1,4], self.C[3,4]]]) + + T = np.array([[self.C[5,5], self.C[1,5], self.C[3,5]], + [self.C[1,5], self.C[1,1], self.C[1,3]], + [self.C[3,5], self.C[1,3], self.C[3,3]]]) + + N1 = np.einsum('ik,kj', -np.linalg.inv(T), np.transpose(R)) + N2 = np.linalg.inv(T) + N3 = np.einsum('ik,kj,jl', R, np.linalg.inv(T), np.transpose(R)) - Q + + N = np.block([[N1, N2], [N3, np.transpose(N1)]]) + + # Solve the eigenvalue problem + Np, Nv = np.linalg.eig(N) + + # Check if the roots are conjugated with each other + for i in range(0,6,2): + if Np[i] != np.conjugate(Np[i+1]): + raise RuntimeError('Roots not in pairs.') + + A = Nv[0:3,0::2] + B = Nv[3:6,0::2] + L = 0.5 * (1.0j * np.einsum('ik,kj', A, np.linalg.inv(B))).real + S = np.array([[np.cos(theta), np.sin(theta), 0.], + [-np.sin(theta), np.cos(theta), 0.], + [0., 0., 1.0]]) + s0 = np.array([np.cos(phi), 0, np.sin(phi)]) + L_t = np.dot(S, np.dot(L, np.transpose(S))) + + coffe_elas = np.einsum('i,ij,j', s0, np.linalg.inv(L_t), np.transpose(s0)) + + ftheta_sigma_xx = ((self.mu2 / np.sqrt(np.cos(theta) + self.mu2 * np.sin(theta)) + - self.mu1 / np.sqrt(np.cos(theta) + self.mu1 * np.sin(theta))) + * self.mu1 * self.mu2 / (self.mu1 - self.mu2)).real + + ftheta_sigma_xy = ((1/np.sqrt(np.cos(theta) + self.mu1 * np.sin(theta)) + - 1/np.sqrt(np.cos(theta) + self.mu2 * np.sin(theta))) + * self.mu1 * self.mu2/(self.mu1 - self.mu2)).real + + ftheta_sigma_yy = ((self.mu1/np.sqrt(np.cos(theta) + self.mu2 * np.sin(theta)) + -self.mu2/np.sqrt(np.cos(theta)+self.mu1*np.sin(theta)))/(self.mu1-self.mu2)).real + + F12_theta = ((ftheta_sigma_yy - ftheta_sigma_xx) * np.sin(theta) * np.cos(theta) + + ftheta_sigma_xy * (np.cos(theta)**2 - np.sin(theta)**2)) + + K_ie_aniso = np.sqrt(1e-3) * np.sqrt(max_gamma * coffe_elas) / F12_theta + + return K_ie_aniso +### def displacement_residuals(r0, crack, x, y, ref_x, ref_y, kI, kII=0, power=1): """