Skip to content

Commit

Permalink
TST: Added more mode II unit tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@Fraser-Birks
Fraser-Birks committed Oct 26, 2023
1 parent 1c00d45 commit 450e0f3
Show file tree
Hide file tree
Showing 2 changed files with 126 additions and 6 deletions.
12 changes: 6 additions & 6 deletions matscipy/fracture_mechanics/crack.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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:
Expand Down
120 changes: 120 additions & 0 deletions tests/test_cubic_crystal_crack.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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:
Expand Down

0 comments on commit 450e0f3

Please sign in to comment.