Int J Fract (2015) 192:47–56 DOI 10.1007/s10704-014-9984-y
ORIGINAL PAPER
Numerical simulation of mode-III fracture incorporating interfacial mechanics Lauren A. Ferguson · Mallikarjunaiah Muddamallappa · Jay R. Walton
Received: 16 September 2014 / Accepted: 5 December 2014 / Published online: 19 December 2014 © Springer Science+Business Media Dordrecht 2014
Abstract Continuum surface methods, including the Sendova–Walton theory, offer a novel approach to fracture modeling in which boundary mechanics are used to augment the classical linear elastic fracture mechanics model for improved prediction of material behavior near fracture surfaces. These methods would be extremely useful in design simulations, but would require numerical implementation which to date has not been available. This has not been previously addressed due to the higher-order tangential derivatives appearing in the fracture surface boundary conditions which make standard implementation techniques, such as the finite element method, a challenge to implement. We propose a method for this implementation which involves reformulating the fracture boundary conditions to remove these higher-order derivatives in the case of mode-III fracture. We also present the initial results of our finite element implementation, which verify the improved stress and displacement field predictions near fracture surfaces. L. A. Ferguson (B) Air Force Research Laboratory, Wright-Patterson AFB, Dayton, OH, USA e-mail:
[email protected] L. A. Ferguson University of Dayton Research Institute, 300 College Park Avenue, Dayton, OH 45469-0060, USA M. Muddamallappa · J. R. Walton Department of Mathematics, Texas A&M University, College Station, TX, USA
Keywords Surface tension · Continuum surface methods · Sendova–Walton fracture theory · Mode-III fracture
1 Introduction Over the last few decades, system design and material qualification have become increasingly reliant on computational modeling tools. Simulation tests may be run significantly faster, cheaper, and under a larger set of conditions than experimental tests, which is particularly advantageous for the continued development and application of advanced materials. However, these computational tools are only as good as the underlying material and mechanical models they employ. We propose a numerical implementation for a promising new fracture model that improves the prediction of material behavior near fracture surfaces. One of the challenges of modeling fracture is that material behavior is dominated by different mechanisms depending on the proximity to a fracture surface. At the macroscale far from any cracks, continuum mechanics is sufficient and the standard linear elastic fracture mechanics (LEFM) model accurately predicts bulk displacement, stress, and strain fields. However, in the neighborhood of a crack, material behavior is governed by atomic forces which are not accounted for by LEFM. Consequently, LEFM is inaccurate in this neighborhood and erroneously predicts an infinite crack-tip stress/strain and an elliptical crack profile.
123
48
Various ideas to augment or replace LEFM with the required interatomic physics have been considered, including the quasicontinuum method (Miller and Tadmor 2007), variational models (Bourdin et al. 2008), and the strictly atomistic computational molecular dynamics (Abraham 2001). However, these models tend to have limited usage to specific length scales and certain mechanical or geometrical conditions, and often require parameter calibration for individual problems. An alternative approach which has garnered a lot of recent attention adds boundary mechanics to LEFM to reduce or eliminate the crack-tip stress singularity. In these so-called continuum surface methods, a twodimensional dividing surface (as defined by Gibbs, see Bumstead et al. 1928) is used to describe the interfacial region around the fracture edge and is endowed with the same physical properties attributed to this region. The choice of constitutive relation for the surface stress tensor has led to two distinct branches of work on this approach. The first, developed in Kim et al. (2010a, b), Antipov and Schiavone (2011), and Kim et al. (2011a, b, c), employs a Gurtin–Murdoch surface elasticity (Gurtin and Murdoch 1978). As shown in Walton (2012) and later in Kim et al. (2013), the use of surface elasticity results in a reduction of the cracktip stress from the square-root singularity of LEFM to a much weaker logarithmic one, but is insufficient to completely remove this singularity. The second branch, developed in Oh et al. (2006), Sendova and Walton (2010a, b), Rajagopal and Walton (2011), Zemlyanova and Walton (2012) and Walton (2014), assumes that the surface stress is a product of the surface tension and a projection tensor into the plane tangent to the deformed fracture surface [as first proposed in Slattery et al. (2004)]. In this Sendova– Walton theory, an appropriately chosen surface tension can reduce or eliminate the strong LEFM crack-tip stress singularity. For example, it is shown in Sendova and Walton (2010a) that for a Griffith crack under pure mode-I loading, a constant surface tension, like the surface elasticity models, results in a logarithmic stress singularity at the crack tips. However, a curvaturedependent surface tension is sufficient to remove this singularity entirely. This choice also predicts a more realistic crack profile with a cusp-like opening shape. Both branches of the continuum surface method could significantly improve system design simulations. However, the direct numerical implementation of this method has not been previously addressed. The finite
123
L. A. Ferguson et al.
element method (FEM), which is typically used to conduct such simulations, is not straightforward to apply here due to the higher-order tangential derivatives that appear in the fracture surface boundary conditions. We have developed a FEM program for the surface tension approach that accomplishes this implementation by transforming the crack surface boundary condition, using its corresponding Green’s function, into an expression that no longer contains higher-order derivatives. The result is a numerical model that can be used to accurately predict displacement, stress, and strain everywhere in the body. In this paper, we focus on quasistatic mode-III fracture under antiplane strain conditions. In this case, we show that the numerical model predicts finite crack-tip stresses and a cusped opening profile, in good agreement with the Sendova–Walton theory.
2 Mode-III fracture problem formulation We consider a two-dimensional (2D) infinite plate Ω ∞ with an under far-field antiplane shear loading σ23 included straight transverse crack Σ with half-length a, as shown in Fig. 1. Note that under antiplane strain conditions, there is no displacement in the x and y directions i.e., u 1 = u 2 = 0, u 3 = u 3 (X 1 , X 2 ).
(1)
We abuse notation throughout by interchanging X 1 , X 2 , X 3 with x, y, z, respectively, e.g., by u 3 we mean displacement in the z-direction. We also nondimensionalize by both the crack half-length (a) and the shear modulus. 2.1 Governing equations and boundary conditions Differential momentum balance (DMB) yields the standard equilibrium governing equation for both the Sendova–Walton and LEFM models, given by u 3 + β = 0, in Ω,
(2)
where β is a body force. Here we have used Hooke’s law as the standard constitutive equation for isotropic linear elasticity. In the case of antiplane shear, this is reduced to τ13 = u 3,1 , (3) τ23 = u 3,2 ,
Numerical simulation of mode-III fracture incorporating interfacial mechanics
49
∞ σ23 ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗
y
(b, b)
Q
L
R
z y
x −a
T
z
1 C (0, 0) ⊗⊗⊗⊗⊗
a
x
B
∞ −σ23
Fig. 2 Finite computational domain Q for 2D mode-III fracture with superposition
and u 3 = 0, on Γ B .
∞ σ23 Fig. 1 Infinite plate domain Ω for 2D mode-III fracture
after nondimensionalization, where τi j are components of the stress tensor. Since the DMB is a linear system, we employ the superposition principle to combine the solution to this system with another solution that has a fixed stress ∞ on the boundary. Effectively, this results τ = −σ23 in moving the far-field loading to the crack surface. We then make use of symmetry about the x and y axes to reduce the problem to the upper right quadrant. Consider the resulting finite computational domain Q = [0, b]2 shown in Fig. 2, where b is the body half-length which we may take as large as required to approximate the infinite problem. The undeformed upper-right crack surface is ΓC = [0, 1]×{0+ } after nondimensionalization. The DMB in (2) is now restricted to this domain. After superposition, the top and right edges are traction-free, i.e., ∇u 3 · n = 0, on ΓT ∪ Γ R ,
(4)
where n is the outward unit normal. We assume frictionless symmetry on the left and bottom edges, so that u 3,1 = 0, on Γ L ,
(5)
(6)
The significant distinction between LEFM and the Sendova–Walton model occurs at the crack surface, where LEFM assumes a traction-free boundary condition. For the Sendova–Walton model, we now ascribe a surface tension which is dependent on the in-plane stretching, i.e., γ˜ = γ0 + γ1 u 3,11 (x, 0).
(7)
The resulting jump momentum balance (JMB) together with the superposed far-field loading gives the following boundary condition over the (upper-right) crack surface ∞ , on ΓC . u 3,2 = −γ1 u 3,111 − σ23
(8)
2.2 Weak formulation For the weak form, we multiply the DMB in (2) by a scalar test function v and integrate over the domain. Integration by parts then yields ∇u 3 · ∇v − v∇u 3 · n = vβ, (9) Q
∂Q
Q
to which we apply the above boundary conditions. The only nonzero boundary integral which remains is that over the crack surface, so that the weak form becomes
123
50
L. A. Ferguson et al.
∇u 3 · ∇v + Q
ΓC
v uC 3,2 =
vβ,
(10)
Q
where we use u C 3,2 as a placeholder for the JMB condition in (8) that must be applied here. We cannot apply this condition directly since it results in a weak form with higher-order derivatives that cannot be eliminated via integration by parts. Instead, we will reformulate this condition using a Green’s function.
To compute the Green’s function, we look at the homogeneous equation 0 = L {y}(x) = −y ,
(17)
and choose two solutions satisfying the boundary conditions (13b) y1 (x) := x y2 (x) := 1 − x
⇒
y1 (0) = 0,
(18)
y2 (1) = 0.
2.3 Reformulation with Green’s function The corresponding Wronskian is For the JMB condition in (8), we rearrange terms to obtain 1 ∞ ], on ΓC . (11) − u 3,111 = [u 3,2 + σ23 γ1 This can be viewed as a second-order ordinary differential equation (ODE) for u 3,1 over the interval [0, 1]. We note that by the symmetry condition on Γ L and regularization on ΓC ∪ Γ B , we require that u 3,1 (0, 0) = u 3,1 (1, 0) = 0.
(12)
Thus the solution to the ODE (11) is the solution to the following boundary value problem (BVP). Definition 1 (Two-point BVP) Find y(x) := u 3,1 (x, 0) solving
W (y1 , y2 )(q) = y1 (q)y2 (q) − y1 (q)y2 (q) = −1. (19) Thus the Green’s function is given by G(x, q) =
(13a)
y(0) = y(1) = 0,
(13b)
where L {·} is the second-order, linear differential operator L {y}(x) = −y (x),
(14)
and the right-hand-side function is given by f (x) =
1 ∞ u 3,2 (x, 0) + σ23 . γ1
(15)
Let G(x, q) denote the Green’s function corresponding to this BVP. Then the solution to (13) is given by 1 u 3,1 (x, 0) = y(x) = G { f }(x) := G(x, q) f (q)dq. 0
(16)
123
(20)
Note that
1 0
x G(x, q) dq = − (1 − x). 2
(21)
Applying these results to the solution (16) yields u 3,1 (x, 0) =
L {y}(x) = f (x), 0 ≤ x ≤ 1,
−q(1 − x), 0 ≤ q ≤ x ≤ 1, −x(1 − q), 0 ≤ x ≤ q ≤ 1.
=
1 γ1 1 γ1
1 0 1
∞ G(x, q)[u 3,2 (q, 0) + σ23 ] dq
G(x, q) u 3,2 (q, 0) dq −
0
∞ σ23 x(1−x). 2γ1
(22) Now recall that the Hilbert transform gives the Dirichletto-Neumann map u 3,2 (x, 0+ ) = H {u 3,1 }
(23)
dq 1 ∞ u 3,1 (q, 0+ ) = − π −∞ q−x
(24)
1 1 2q = − u 3,1 (q, 0+ ) 2 dq, π 0 q − x2
(25)
where − indicates a Cauchy principal value and we have used the fact that u 3,1 (x, 0) is odd and u 3 vanishes outside the crack surface. Applying the BVP solution (22) to this equation yields
Numerical simulation of mode-III fracture incorporating interfacial mechanics
1 1 1 2q G(q, r ) u 3,2 (r, 0) dr 2 dq − π γ1 0 0 q − x2 1 σ∞ 2q dq. (26) − 23 − q(1 − q) 2 2π γ1 0 q − x2
u 3,2 (x, 0+ ) =
Since the Green’s function is smooth, we change the order of integration in the double integral term and swap dummy variables to obtain 1 h(x) := −
2r G(r, q) u 3,2 (q, 0) dq 2 dr r − x2 0 0 1 1 2r = u 3,2 (q, 0) dq − G(r, q) 2 dr r − x2 0 0 1 = k(x, q) u 3,2 (q, 0) dq, (27) 0
1 2r k(x, q) := − G(r, q) 2 dr r − x2 0 q 1 1 dr = −r (1 − q) + r+x r−x 0 1 1 1 dr + −q(1 − r ) + r + x r − x q q q q dr dr = −(1 − q) 2 dr − x +x 0 0 r +x 0 r −x 1 1 1 dr dr +q 2 dr − (x + 1) + (x − 1) q q r +x q r −x = (q + x) ln(q + x) + (q − x) ln |q − x| −q(1 + x) ln(1 + x) − q(1 − x) ln(1 − x).
(28)
1 1 ∇u 3 · ∇v + v(x, 0) π γ1 0 Q 1 · k(x, q) u 3,2 (q, 0) dq d x 0
q(1 − q) dq + q+x
=
vβ + Q
∞ σ23 2π γ1
1
v(x, 0) g(x) d x.
(31)
0
Note that this weak form has no higher-order derivatives, thus the standard FEM can now be applied. Specifically, we will solve the following.
A(u, v) = L(v), ∀v ∈ V,
(32)
where A(·, ·) and L(·) are the bilinear and linear forms, respectively, given by A(u, v) =
∇u · ∇v Q
+
1
v(x, 0)
0
0
L(v) =
1
vβ +
∂ K (x, q) u(q, y) dq d x, ∂y y=0
(33) 1
v(x, 0) G(x) d x,
(34)
0
Q
where the kernels K and G are given by
The single integral term in (26) simplifies as q(1 − q) dq g(x) := q−x 0 0 1 1 1 dq (1 − q + x) dq − x(1 + x) = q + x 0 0 1 1 1 dq (1 − q − x) dq + x(1 − x) + 0 0 q−x
1−x 1+x + x(1 − x) ln . = 1 − x(1 + x) ln x x 1
Definition 2 (Variational Problem) Find u ∈ V such that
where the inner integral simplifies to
where k(x, q) and g(x) are given in (28) and (29), respectively. Applying this result to (10) yields the final weak form
1
51
1
1 k(x, q), π γ1 σ∞ G(x) = 23 g(x), 2π γ1
K (x, q) =
(35) (36)
and V is the solution and test space given by
V = w(x) ∈ H (Q) : w
1
ΓB
=0 .
(37)
(29) Applying these results to (26) yields the reformulated JMB condition 1 1 C k(x, q) u 3,2 (q, 0) dq u 3,2 (x, 0) = π γ1 0 σ∞ (30) − 23 g(x), on ΓC , 2π γ1
3 FEM implementation We used two different methods for implementing this problem using finite elements. In the first, we solve the Fredholm integral equation which defines the reformulated JMB condition (30) to obtain a numerical
123
52
L. A. Ferguson et al.
approximation to the Neumann boundary condition for u 3,2 (x, 0) over the crack surface. The result can be applied to the original weak formulation (10) with standard finite elements. In the second method, we solve the variational problem given in Definition 2 with the reformulated JMB condition included directly. We call this a nonlocal boundary FEM approach since the reformulated JMB (30) is a nonlocal integral operator on the crack surface. We will describe these methods in more detail, and show in Sect. 5 that they are in good agreement. In each case, we used the deal. II library (Bangerth et al. 2007, 2013) to develop the code for FEM execution. For simplicity, we currently neglect body forces (i.e., β = 0). However, we note that they may be easily added if desired.
where u(qi ) is the solution to (39) at the quadrature point qi ∈ [0, 1]. The general solution to (39) at any point x ∈ [0, 1] can then be obtained using the Nyström method (Atkinson 1997)
3.1 Fredholm approximation approach
uC 3,2 (x, 0)
∇u 3 · ∇v + Q
ΓC
v uC 3,2 = 0,
0
(39) where K (x, q) and G(x) are given in (35) and (36), respectively. For the FEM implementation of this weak form, we assemble the system using Gaussian quadrature to approximate the integrals. In this case, all we require is the value of u C 3,2 at the quadrature points. To obtain these values, we have developed an algorithm for numerically solving the Fredholm integral equation above. We use composite Simpson’s rule to approximate the integral term, i.e., 1 n
K (x, q)u(q) dq ≈ w j K (x, q j )u(q j ), (40) j=0
where we have discretized the domain [0, 1] by a set of quadrature points {q j }nj=0 (for n even) with corresponding weights {w j }. Although the kernel K given in (35) and (28) is not smooth, all its singular points are removable. For example, via L’Hôpital’s rule, lim (q − x) ln |q − x| = 0.
x→q
123
n
w j K (qi , q j )u(q j )
j=0
= −G(qi ), for i = 0, . . . , n,
≈
n
w j K (x, q j )u(q j ) + G(x).
(42)
(43)
This formula is used to obtain values of u C 3,2 (x, 0) at the Gauss quadrature points required in the FEM system matrix assembly of the original weak form (38).
(38)
where u C 3,2 (x, 0) solves the Fredholm integral equation of the second kind [cf. (30)] 1 u(x)− K (x, q) u(q) dq = −G(x), for 0 ≤ x ≤ 1,
0
u(qi ) −
j=0
We consider the original weak form (10) in the absence of body forces
Since for a fixed x the number of singular points is finite, i.e., a set of measure zero in the interval [0, 1], we can equivalently replace the kernel by its continuous extension. We use collocation to obtain the approximate value of the solution function at the Simpson quadrature points, i.e., we solve the system
(41)
3.2 Nonlocal boundary approach One of the main advantages of FEM is that the assembly of the stiffness matrix can be split up over the elements, with the contribution over each element computed locally and then added into the global matrix. For the variational problem described in Definition 2, this is exactly what is done in the bulk, i.e., for the first term in the bilinear form (33). However, we must modify how we compute these contributions for the boundary term
1 0
1
v(x, 0) 0
∂ K (x, q) u(q, y) dq d x, ∂y y=0
(44)
resulting from the nonlocal integral operator in the reformulated JMB (30). In particular, after discretizing this double integral term over the triangulation, we must keep track of two elements: one over which the outer integral is taken and one for the inner integral. In this sense, we say that the FEM approach is nonlocal on the boundary. Fortunately, such nonlocal terms, like the standard local integrals, may still be computed via Gaussian quadrature. Recall that in FEM the test and solution
Numerical simulation of mode-III fracture incorporating interfacial mechanics
functions are first approximated by the set of global shape functions {ϕi } so that the only nonzero contributions of the double integral occur when there is shape function support over both elements. Let K 1 and K 2 be two elements from our triangulation that have nonzero intersection with the crack surface ΓC . Then we must assemble double-integral terms of the form D=
K 1 ∩ΓC
ϕ1 (x, 0)
K 2 ∩ΓC
∂ ϕ2 (q, y) K (x, q) dq d x, ∂y y=0
(45) where ϕ1 and ϕ2 have support over K 1 and K 2 , respectively. Although arbitrary quadrilaterals could be used, we assume a rectangular mesh for simplicity, so that these elements have the form K 1 = [a, b] × [0, c],
K 2 = [A, B] × [0, C].
(46)
We also structure the mesh so that the crack tip coincides with an element vertex. In this case, the standard bilinear mappings which transform the unit reference element Kˆ = [0, 1]2 to these real elements reduce to x(x) ˆ = a +(b − a)x, ˆ
q(q) ˆ = A + (B − A)q, ˆ
(47)
ˆ 0 = 0,
y( yˆ ) = C yˆ .
(48)
Applying these mappings converts the double integral in (45) to D=
1
1
ϕ1 (x(x), ˆ 0)
K (x(x), ˆ q(q)) ˆ ∂ yˆ ∂ ϕ2 (q(q), (B − A) d qˆ (b−a) d xˆ ˆ y( yˆ )) · ∂ yˆ yˆ =0 ∂ y (49) 1 1 = ϕˆ 1 (x, ˆ 0) K (x(x), ˆ q(q)) ˆ 0 0 1 ∂ ϕˆ2 (q, (b−a)(B − A) d qˆ d x, ˆ (50) ˆ yˆ ) · ∂ yˆ yˆ =0 C 0
0
where ϕˆ1 and ϕˆ2 are the corresponding reference shape functions. Finally, we replace both integrals with a 1D Gauss quadrature formula
1 0
f (x) ˆ d xˆ =
n
t=1
ωˆ t f (xˆt ),
(51)
53
where {xˆt } and {ωˆ t } are the quadrature points and weights, respectively. This yields n B−A ωˆ t ωˆ s ϕˆ1 (xˆt , 0) K (x(xˆt ), q(xˆs )) C t,s ∂ ϕˆ2 (xˆs , yˆ ) . · (52) ∂ yˆ
D = (b−a)
yˆ =0
This term is assembled in the stiffness matrix as usual, with an extra loop for the additional quadrature rule. However, it should be noted that this will typically result in a nonsymmetric matrix. In addition, note that the kernel function K is computed at the global quadrature points in the real elements K 1 and K 2 , not at the reference quadrature points. In effect, this approach is doing the same thing as the Fredholm approximation approach discussed in Sect. 3.1, in that the interior integral is approximated using a quadrature rule. The significant difference is that in this nonlocal boundary method, the collocation points over the surface are not solved for separately, but rather included and solved for in the global system.
4 Parameter determination The only undetermined parameter in the Sendova– Walton theory for mode-III fracture is the surface tension constant γ1 . This value is a material property that can be determined using calibration experiments. However, even without this calibration, we can still get an idea of the correct range for this value. In this section, we will prove that there is a lower bound on this range that will guarantee bounded crack-tip stresses and a cusped opening profile. We showed in Sects. 2.3 and 3.1 that the JMB condition may be reformulated as a Fredholm equation of the second kind (39). In standard form, this equation is given by γ1 u(x) − K [u](x) = −
∞ σ23 g(x), for 0 ≤ x ≤ 1, 2π (53)
where K is the integral operator 1 1 k(x, q) u(q) dq, (54) K [u](x) = π 0 and the kernel k and data function g are given by (28) and (29), respectively.
123
54
L. A. Ferguson et al. 0.0005
Theorem 1 The Fredholm integral equation in (53) has a unique, continuous solution u(x) ∈ [0, 1] for all but countably many values of γ1 .
Since the JMB condition requires that u 3,2 (x, 0) be a solution to this Fredholm equation, then by Hooke’s law (3), we have shown that the shear stress τ23 along the crack surface must be bounded for all but countably many values of γ1 . To find a lower bound on these values, we need to bound the spectral radius of the integral operator K , which may be obtained using the Hilbert-Schmidt norm. Thus, we require that |γ1 | > ||K || H S =
1 k(x, q) L 2 ≥ ρ(K ). π
(55)
We used Mathematica to approximate the value of this norm using numerical integration. The commands we used are shown in the code snippet below, which resulted in the value γ1min = 0.1023473774. Thus, we are guaranteed bounded stresses by taking |γ1 | > γ1min . Listing 1. Mode-III kernel L 2 -norm approximation
k[x_,q_] := (q+x)∗Log[Abs[q+x]] − q∗(1+x)∗Log[Abs[1+x]] + (q−x)∗Log[Abs[q−x]] − q∗(1−x)∗Log[Abs[1−x] ] ; val = NIntegrate[ (k[x,q])^2 , {x,0 ,1} , {q,0 ,1} , MaxRecursion→20, WorkingPrecision→20]; norm = NumberForm[Sqrt[ val ] /Pi, 10]
123
0.0004
Displacement (u3)
Proof We already noted in (41) that k(x, q) has a limit of zero when x → q, and this holds even for q = 1. Similarly, the other logarithmic singularities are removable by their coefficients. Thus the kernel is a continuous function everywhere in [0, 1]2 except for the set of removable discontinuity points S = {(x, q)} ∈ [0, 1]2 : x = {1, q}, which is a set of measure zero in R2 . Thus k is square integrable, i.e., k(x, q) ∈ L 2 ([0, 1]2 ). It follows by definition that K is a Hilbert-Schmidt operator. As a consequence, K is bounded, compact, and has a countable spectrum. Applying the Fredholm Alternative, we see that (53) has a unique solution in C[0, 1] if and only if γ1 is not an eigenvalue of K . Since the spectrum is countable, this concludes the proof.
Nonlocal Boundary Fredholm Approximation
0.0003
0.0002
0.0001
0 0
0.2
0.4
0.6
0.8
1
Position (x) on crack surface ∞ = 0.005, Fig. 3 Comparison of FEM implementations for σ23 γ1 = 1
5 Simulation results In this section, we show some displacement and stress results for our simulation. We take throughout a body half-length b = 8. In each case, we ran the code using an initially uniform coarse mesh which was adaptively refined twelve times so that the final mesh had more than 1,300,000 elements and degrees of freedom, with over 1,400 elements along the crack surface. The results of our two different FEM implementation approaches agree very closely, as shown in Fig. 3. We also computed the 2 -norm of the difference of the crack surface displacement vectors for these approaches as 1.9312E−04, which is reasonable given that they use different quadrature rules. For the remainder of this section, all results were computed using the nonlocal boundary implementation. The crack surface displacement in the above figure shows the cusp-like opening predicted by the Sendova– Walton theory for appropriate surface tension values. Figure 4 further demonstrates how the surface tension parameter γ1 affects the crack shape. For values above the threshold computed in Sect. 4 (γ1min ≈ 0.1), we see the expected cusped opening shape and significantly reduced displacement. Even for smaller nonzero values we see cusping, although the rest of the surface profile is closer to the standard elliptical LEFM solution, which is recovered for γ1 = 0. A similar comparison of the out-of-plane shear stress (τ23 ) near the crack tip is shown in Fig. 5. This stress is substantially reduced for nonzero surface tension values, even more so for those values above γ1min for which
Numerical simulation of mode-III fracture incorporating interfacial mechanics 0.006
0 0.01 0.05 0.5 1 4
∞= σ23 0.005 0.01 0.05
0.004
Displacement (u 3 )
Displacement (u 3 )
0.005
γ1 =
0.005
55
0.004 0.003 0.002
0.003
0.002
0.001
0.001
0
0 0
0.2
0.4
0.6
0.8
1
0
∞ = 0.005 Fig. 4 Surface displacement for varying γ1 with σ23
0.05
0.6
0.0035
0 0.01 0.05 0.5 1 4
0.8
1
∞= σ23 0.005 0.01 0.05
0.003 0.0025
Stress (τ23 )
Stress (τ23 )
0.4
∞ with γ = 1 Fig. 6 Surface displacement for varying σ23 1
γ1 =
0.04
0.2
Position (x) on crack surface
Position (x) on crack surface
0.03
0.02
0.002 0.0015 0.001
0.01
0
0.0005 0
1
1.02
1.04
1.06
1.08
1.1
Position (x) outside crack surface ∞ = 0.005 Fig. 5 Near-tip stress for varying γ1 with σ23
the theory predicts finite stresses. We also looked at the crack surface displacement and near-tip stress for various loading values, as shown in Figs. 6 and 7, respectively. As expected, we observe higher displacement and stress values for larger loading.
6 Conclusions We have successfully demonstrated an approach for the direct numerical implementation of the surface tension class of continuum surface methods (i.e., the Sendova– Walton theory) using FEM in the case of mode-III fracture. This approach reformulates the JMB condition using its corresponding Green’s function, which removes the inherent higher-order tangential derivatives. The actual FEM implementation was conducted
1
1.02
1.04
1.06
1.08
1.1
Position (x) outside crack surface ∞ with γ = 1 Fig. 7 Near-tip stress for varying σ23 1
using two different methods. In the first, the reformulated JMB condition was solved separately to obtain a numerical approximation to the Neumann condition on the crack surface, which was then included using a standard FEM procedure. In the second, a nonlocal boundary FEM computation was employed to include the Neumann condition directly during the assembly of the global system matrix. We showed in Sect. 5 that the FEM implementations agree well with each other. More importantly, they both confirm the predictions of the Sendova–Walton theory. In particular, a nonzero stretching-dependent surface tension produces bounded crack-tip stresses and a cusplike crack opening profile. The FEM implementation of the Sendova–Walton theory allows it to be applied almost immediately to design simulations, providing improved behavior prediction near fracture surfaces.
123
56
We are currently developing a corresponding implementation for mode-I fracture and plan to extend the ideas presented here to mixed-mode fracture and more complex geometries. In addition, we expect that this implementation could be extended to the general theory developed in Sendova and Walton (2010a) for nonlinear elasticity and nonlinear surface mechanics. Acknowledgments This work was supported by the Air Force Office of Scientific Research project 12RX13COR.
References Abraham FF (2001) The atomic dynamics of fracture. J Mech Phys Solids 49(9):2095–2111 Antipov YA, Schiavone P (2011) Integro-differential equation for a finite crack in a strip with surface effects. Q J Mech Appl Math 64(1):87–106 Atkinson KE (1997) The numerical solution of integral equations of the second kind. Cambridge University Press, Cambridge Bangerth W, Hartmann R, Kanschat G (2007) Deal. II— a general-purpose object-oriented finite element library. ACM Trans Math Softw 33(4):24/1–24/27 Bangerth W, Heister T, Heltai L, Kanschat G, Kronbichler M, Maier M, Turcksin B, Young T (2013) The deal. ii library, verion 8.1. arXiv preprint: arxiv.org/abs/1312.2266v4 Bourdin B, Francfort GA, Marigo J-J (2008) The variational approach to fracture. J Elast 91(1–3):5–148 Bumstead HA, Van Name RG, Longley WR (eds) (1928) The collected works of J. Willard Gibbs, vol 1. Longmans, Green, and Co., New York Gurtin ME, Murdoch AI (1978) Surface stress in solids. Int J Solids Struct 14(6):431–440 Kim CI, Schiavone P, Ru C-Q (2010a) The effects of surface elasticity on an elastic solid with mode-III crack: complete solution. J Appl Mech Trans ASME 77(2): 021011 Kim CI, Schiavone P, Ru C-Q (2010b) Analysis of a modeIII crack in the presence of surface elasticity and a prescribed non-uniform surface traction. Z Angew Math Phys 61(3):555–564
123
L. A. Ferguson et al. Kim CI, Schiavone P, Ru C-Q (2011a) Analysis of plane-strain crack problems (mode-I & mode-II) in the presence of surface elasticity. J Elast 104:397–420 Kim CI, Schiavone P, Ru C-Q (2011b) The effect of surface elasticity on a mode-III interface crack. Arch Mech 63(3):267– 286 Kim CI, Schiavone P, Ru C-Q (2011c) Effect of surface elasticity on an interface crack in plane deformations. Proc R Soc A 467(2136):3530–3549 Kim CI, Ru C-Q, Schiavone P (2013) A clarification of the role of crack-tip conditions in linear elasticity with surface effects. Math Mech Solids 18(1):59–66 Miller RE, Tadmor EB (2007) Hybrid continuum mechanics and atomistic methods for simulating materials deformation and failure. MRS Bull 32(11):920–926 Oh E-S, Walton JR, Slattery JC (2006) A theory of fracture based upon an extension of continuum mechanics to the nanoscale. J Appl Mech - T ASME 73(5):792–798 Rajagopal KR, Walton JR (2011) Modeling fracture in the context of a strain-limiting theory of elasticity: a single antiplane shear crack. Int J Fract 169(1):39–48 Sendova T, Walton JR (2010a) A new approach to the modeling and analysis of fracture through extension of continuum mechanics to the nanoscale. Math Mech Solids 15(3):368– 413 Sendova T, Walton JR (2010b) The effect of surface tension in modeling interfacial fracture. In: Todorov MD, Christov CI (eds) Application of mathematics in technical and natural sciences, AIP conference proceedings, vol 1301, pp 291– 300 Slattery JC, Oh E-S, Fu KB (2004) Extension of continuum mechanics to the nanoscale. Chem Eng Sci 59(21):4621– 4635 Walton JR (2012) A note on fracture models incorporating surface elasticity. J Elast 109(1):95–102 Walton JR (2014) Plane-strain fracture with curvature-dependent surface tension: mixed-mode loading. J Elast 114(1):127– 142 Zemlyanova AY, Walton JR (2012) Modeling of a curvilinear planar crack with a curvature-dependent surface tension. SIAM J Appl Math 72(5):1474–1492