From 450e0f3432ea8b32f8b25d1afcbf8342dd42ffc7 Mon Sep 17 00:00:00 2001 From: Fraser-Birks Date: Thu, 26 Oct 2023 15:37:25 +0100 Subject: [PATCH] TST: Added more mode II unit tests --- matscipy/fracture_mechanics/crack.py | 12 +-- tests/test_cubic_crystal_crack.py | 120 +++++++++++++++++++++++++++ 2 files changed, 126 insertions(+), 6 deletions(-) diff --git a/matscipy/fracture_mechanics/crack.py b/matscipy/fracture_mechanics/crack.py index 64a60990..0ccbd4f7 100644 --- a/matscipy/fracture_mechanics/crack.py +++ b/matscipy/fracture_mechanics/crack.py @@ -222,11 +222,11 @@ def deformation_gradient(self, r, theta, kI, kII=0.0): f_II = kII / np.sqrt(2*math.pi*r) h4 = self.inv_mu1_mu2 - du_dx_II = f_I*( self.inv_mu1_mu2*( self.p2/h2 - self.p1/h3 ) ).real - du_dy_II = f_I*( h4*( (self.mu2*self.p2)/h2 - (self.mu1*self.p1)/h3 ) ).real + du_dx_II = f_II*( self.inv_mu1_mu2*( self.p2/h2 - self.p1/h3 ) ).real + du_dy_II = f_II*( h4*( (self.mu2*self.p2)/h2 - (self.mu1*self.p1)/h3 ) ).real - dv_dx_II = f_I*( self.inv_mu1_mu2*( self.q2/h2 - self.q1/h3 ) ).real - dv_dy_II = f_I*( h4*( (self.mu2*self.q2)/h2 - (self.mu1*self.q1)/h3 ) ).real + dv_dx_II = f_II*( self.inv_mu1_mu2*( self.q2/h2 - self.q1/h3 ) ).real + dv_dy_II = f_II*( h4*( (self.mu2*self.q2)/h2 - (self.mu1*self.q1)/h3 ) ).real du_dx = du_dx_I + du_dx_II du_dy = du_dy_I + du_dy_II @@ -493,7 +493,7 @@ def displacements_from_cartesian_coordinates(self, dx, dy, kI, kII=0): """ abs_dr = np.sqrt(dx*dx+dy*dy) theta = np.arctan2(dy, dx) - return self.displacements_from_cylinder_coordinates(abs_dr, theta, kI, kII=0) + return self.displacements_from_cylinder_coordinates(abs_dr, theta, kI, kII=kII) def displacements(self, ref_x, ref_y, x0, y0, kI, kII=0): @@ -1689,7 +1689,7 @@ def arc_length_continuation(self, x0, x1, N=10, ds=0.01, ftol=1e-2, with h5py.File(traj_file, 'a') as hf: x_traj = hf['x'] x_traj.resize((x_traj.shape[0] + 1, x_traj.shape[1])) - x_traj[-1, :] = x2 + x_traj[-1, :] = np.append(x2[:-1],[self.kI,self.kII]) row += 1 break except: diff --git a/tests/test_cubic_crystal_crack.py b/tests/test_cubic_crystal_crack.py index 08adfcd1..da45d877 100644 --- a/tests/test_cubic_crystal_crack.py +++ b/tests/test_cubic_crystal_crack.py @@ -201,6 +201,7 @@ def test_anisotropic_near_field_solution(self): self.assertTrue(np.max(residual[mask]) < 0.2) def test_consistency_of_deformation_gradient_and_stress(self): + #test pure mode 1 first for C11, C12, C44, surface_energy, k1 in self.materials: crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44) k = crack.k1g(surface_energy)*k1 @@ -224,8 +225,60 @@ def test_consistency_of_deformation_gradient_and_stress(self): self.assertAlmostEqual(eps[0, 0], eps_xx) self.assertAlmostEqual(eps[1, 1], eps_yy) self.assertAlmostEqual(eps[0, 1], eps_xy) + + #now test mode 2 + for C11, C12, C44, surface_energy, k1 in self.materials: + crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44) + k = crack.k1g(surface_energy)*k1 + for i in range(10): + x = np.random.uniform(-10, 10) + y = np.random.uniform(-10, 10) + F = crack.deformation_gradient(x, y, 0, 0, 0, kII=k) + eps = (F+F.T)/2-np.eye(2) + sig_x, sig_y, sig_xy = crack.stresses(x, y, 0, 0, 0,kII=k) + eps_xx = crack.crack.a11*sig_x + \ + crack.crack.a12*sig_y + \ + crack.crack.a16*sig_xy + eps_yy = crack.crack.a12*sig_x + \ + crack.crack.a22*sig_y + \ + crack.crack.a26*sig_xy + # Factor 1/2 because elastic constants and matrix product are + # expressed in Voigt notation. + eps_xy = (crack.crack.a16*sig_x + \ + crack.crack.a26*sig_y + \ + crack.crack.a66*sig_xy)/2 + self.assertAlmostEqual(eps[0, 0], eps_xx) + self.assertAlmostEqual(eps[1, 1], eps_yy) + self.assertAlmostEqual(eps[0, 1], eps_xy) + + #now test mixed mode + for C11, C12, C44, surface_energy, k1 in self.materials: + crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44) + k = crack.k1g(surface_energy)*k1 + for i in range(10): + x = np.random.uniform(-10, 10) + y = np.random.uniform(-10, 10) + F = crack.deformation_gradient(x, y, 0, 0, k/2, k/2) + eps = (F+F.T)/2-np.eye(2) + sig_x, sig_y, sig_xy = crack.stresses(x, y, 0, 0, k/2, k/2) + eps_xx = crack.crack.a11*sig_x + \ + crack.crack.a12*sig_y + \ + crack.crack.a16*sig_xy + eps_yy = crack.crack.a12*sig_x + \ + crack.crack.a22*sig_y + \ + crack.crack.a26*sig_xy + # Factor 1/2 because elastic constants and matrix product are + # expressed in Voigt notation. + eps_xy = (crack.crack.a16*sig_x + \ + crack.crack.a26*sig_y + \ + crack.crack.a66*sig_xy)/2 + self.assertAlmostEqual(eps[0, 0], eps_xx) + self.assertAlmostEqual(eps[1, 1], eps_yy) + self.assertAlmostEqual(eps[0, 1], eps_xy) + def test_consistency_of_deformation_gradient_and_displacement(self): + #test mode 1 first eps = 1e-3 for C11, C12, C44, surface_energy, k1 in self.materials: crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44) @@ -249,7 +302,54 @@ def test_consistency_of_deformation_gradient_and_displacement(self): F_num = np.transpose([[du_dx, du_dy], [dv_dx, dv_dy]]) self.assertArrayAlmostEqual(F, F_num, tol=1e-4) + #now test mode 2 + for C11, C12, C44, surface_energy, k1 in self.materials: + crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44) + k = crack.k1g(surface_energy)*k1 + for i in range(10): + x = i+1 + y = i+1 + F = crack.deformation_gradient(x, y, 0, 0, 0, k) + # Finite difference approximation of deformation gradient + u, v = crack.displacements(x, y, 0, 0, 0, k) + ux0, vx0 = crack.displacements(x-eps, y, 0, 0, 0, k) + uy0, vy0 = crack.displacements(x, y-eps, 0, 0, 0, k) + ux1, vx1 = crack.displacements(x+eps, y, 0, 0, 0, k) + uy1, vy1 = crack.displacements(x, y+eps, 0, 0, 0, k) + du_dx = (ux1-ux0)/(2*eps) + du_dy = (uy1-uy0)/(2*eps) + dv_dx = (vx1-vx0)/(2*eps) + dv_dy = (vy1-vy0)/(2*eps) + du_dx += np.ones_like(du_dx) + dv_dy += np.ones_like(dv_dy) + F_num = np.transpose([[du_dx, du_dy], [dv_dx, dv_dy]]) + self.assertArrayAlmostEqual(F, F_num, tol=1e-4) + + #now test mixed mode + for C11, C12, C44, surface_energy, k1 in self.materials: + crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44) + k = crack.k1g(surface_energy)*k1 + for i in range(10): + x = i+1 + y = i+1 + F = crack.deformation_gradient(x, y, 0, 0, k/2, k/2) + # Finite difference approximation of deformation gradient + u, v = crack.displacements(x, y, 0, 0, k/2, k/2) + ux0, vx0 = crack.displacements(x-eps, y, 0, 0, k/2, k/2) + uy0, vy0 = crack.displacements(x, y-eps, 0, 0, k/2, k/2) + ux1, vx1 = crack.displacements(x+eps, y, 0, 0, k/2, k/2) + uy1, vy1 = crack.displacements(x, y+eps, 0, 0, k/2, k/2) + du_dx = (ux1-ux0)/(2*eps) + du_dy = (uy1-uy0)/(2*eps) + dv_dx = (vx1-vx0)/(2*eps) + dv_dy = (vy1-vy0)/(2*eps) + du_dx += np.ones_like(du_dx) + dv_dy += np.ones_like(dv_dy) + F_num = np.transpose([[du_dx, du_dy], [dv_dx, dv_dy]]) + self.assertArrayAlmostEqual(F, F_num, tol=1e-4) + def test_elastostatics(self): + #check for pure mode I eps = 1e-3 for C11, C12, C44, surface_energy, k1 in self.materials: crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44) @@ -268,6 +368,26 @@ def test_elastostatics(self): # equilibrium) self.assertAlmostEqual(divsx, 0.0, places=3) self.assertAlmostEqual(divsy, 0.0, places=3) + eps = 1e-3 + + #also check for mixed mode + for C11, C12, C44, surface_energy, k1 in self.materials: + crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44) + k = crack.k1g(surface_energy)*k1 + for i in range(10): + x = i+1 + y = i+1 + # Finite difference approximation of the stress divergence + sxx0x, syy0x, sxy0x = crack.stresses(x-eps, y, 0, 0, k/3,(2*k)/3) + sxx0y, syy0y, sxy0y = crack.stresses(x, y-eps, 0, 0, k/3,(2*k)/3) + sxx1x, syy1x, sxy1x = crack.stresses(x+eps, y, 0, 0, k/3,(2*k)/3) + sxx1y, syy1y, sxy1y = crack.stresses(x, y+eps, 0, 0, k/3,(2*k)/3) + divsx = (sxx1x-sxx0x)/(2*eps) + (sxy1y-sxy0y)/(2*eps) + divsy = (sxy1x-sxy0x)/(2*eps) + (syy1y-syy0y)/(2*eps) + # Check that divergence of stress is zero (elastostatic + # equilibrium) + self.assertAlmostEqual(divsx, 0.0, places=3) + self.assertAlmostEqual(divsy, 0.0, places=3) def test_J_integral(self): if quadrature is None: