Spectral formulation of the boundary integral equation method for antiplane problems
SSpectral formulation of the boundary integral equation method for antiplane problems
Kunnath Ranjith École Centrale School of Engineering Mahindra University Bahadurpally, Jeedimetla Telangana 500043, India E-mail: [email protected]
Abstract
A spectral formulation of the boundary integral equation method for antiplane problems is presented. The boundary integral equation method relates the displacement discontinuity and the shear stress at an interface between two half-planes. It involves evaluating a space-time convolution of the shear stress or the displacement discontinuity at the interface. In the spectral formulation, the convolution with respect to the spatial coordinate is performed in the spectral domain. The leads to greater numerical efficiency. Prior work on the spectral formulation of the boundary integral equation method has performed the elastodynamic convolution of the displacement discontinuity at the interface. In the present work, the convolution is performed of the shear stress at the interface. The formulation is validated by numerically calculating the response of the interface to harmonic and to impulsive disturbances, and comparing with known analytical solutions. To illustrate use of the method, dynamic rupture propagation with a slip-weakening friction law is simulated.
Keywords:
Boundary integral equation, elasticity, waves, fracture, spectral, modal ntroduction :
Failure of interfaces between materials due to dynamic crack or rupture propagation is often observed in nature and in technology. The most striking case of such failure is earthquake rupture. Technological materials such as composites also frequently fail by dynamic cracking or debonding at interfaces. Study of conditions leading to crack or rupture nucleation, propagation and arrest is of considerable theoretical and practical interest. The presence of strong non-linearities such as friction and plasticity, and the complicated dynamics due to elastic wave propagation gives rise to considerable richness in such problems. Several numerical studies of dynamic crack and rupture propagation have been reported in the literature. Andrews and Ben-Zion (1997) and Harris and Day (1997) have developed finite difference methods for the study of dynamic rupture propagation. Finite element simulations of interface rupture propagation have been reported by Baillet et al. (2005), Kammer et al. (2012), Kammer et al. (2014), Di Bartolomeo et al. (2012) and Tonazzi et al. (2013). The boundary integral equation method is also widely used to study dynamic fracture problems when the fracture is confined to a plane. The main advantage of the method is that field quantities on the fracture plane can be studied without the need to evaluate them at locations away from the fracture plane. Conventional computational methods such as the finite element method and the finite difference method require discretization of the complete domain, and are therefore less efficient. The boundary integral equation method was introduced by Kostrov (1966). In his formulation, the displacement discontinuities (i.e. slip or crack opening) on the fracture plane were expressed as a space-time convolution of the traction components of stress on the plane. An equivalent form was introduced by Budiansky and Rice (1979) in which the traction components on the fracture plane were written as a space-time convolution of the displacement continuities on the fracture plane. A spectral ormulation of the Budiansky-Rice approach was proposed by Morrissey and Geubelle (1997) for antiplane problems and by Geubelle and Rice (1995) for general 3D problems. In the spectral formulation, the convolution over the spatial coordinates is substituted with a multiplication operation in the spectral domain. The formulation relies on the numerical efficiency of the fast Fourier transform in comparison to numerical integration. In addition, due to the spectral nature of the method, it is suited to parallel computing. The spectral formulation for planar bi-material interfaces was introduced by Geubelle and Breitenfeld (1997) and Breitenfeld and Geubelle (1998). Hajarolasvadi and Elbanna (2017) have developed a hybrid finite difference – spectral boundary integral equation method scheme for rupture propagation while Ma et al. (2018) have proposed a hybrid finite element – spectral boundary integral equation method scheme. Ranjith (2015) developed a spectral formulation of the original equations of Kostrov (1966) for 2D plane strain. The method involves performing convolutions over the time-history of the traction components of stress at the interface. This is to be contrasted with the convolutions done over the time-history of the displacement discontinuity at the interface (or of the displacements at the interface) in previous work on the spectral formulation. In the present paper, the work of Ranjith (2015) is extended by developing the spectral formulation for 2D antiplane strain problems. The formulation is also validated and implemented. An illustration of the application of the method to simulating dynamic rupture propagation is given.
Governing equations :
Consider a planar interface between two elastic half-planes as in Fig. 1. It is assumed that both half-planes have the same shear modulus, , and mass density, . Let a Cartesian µ ρ oordinate system, , be located such that the interface between the half-planes is at and the - coordinate is normal to the half-planes. Further, let denote the time. We assume that all field quantities are independent of the - coordinate. Further, we assume that the displacement field only has the antiplane component, i.e., the displacement field is given by (1) Consequently, the non-zero components of the strain field are (2) The corresponding stress components are (3) The linear momentum balance equation is (4) Substituting Eqs. (3) and (2) into the above equation, we get the 2D scalar wave equation for the displacement field, , (5) where is the shear wave speed of the solid. x i , i = x = x t x u = u = u = u ( x , x , t ) ε = ε = ∂ u ∂ x , ε = ε = ∂ u ∂ x . σ = µ ε = σ , σ = µ ε = σ . ∂ σ ∂ x + ∂ σ ∂ x = ρ ∂ u ∂ t ∂ u ∂ x + ∂ u ∂ x = c s ∂ u ∂ t c s = µ / ρ e expand the displacement in a Fourier series with respect to the coordinate along the interface, (6) where is the spectral amplitude and is the wavenumber. Considering a single spectral component, and substituting into Eq. (5), we get (7) Further, taking the Laplace transform ,defined by (8) where is the Laplace variable, we get (9) The above is an ordinary differential equation for . It has two solutions (10) where (11) and is an arbitrary amplitude. In the upper half-plane, since the displacement field has to be bounded, we neglect the solution in Eq. (10) with the positive exponent. Thus, the single spectral component of the displacement (or more precisely, its Laplace transform) which satisfies the governing equation is u ( x , x , t ) = Ψ ( k , x , t ) e ikx k ∑ , Ψ k − k Ψ + ∂Ψ∂ x = c s ∂ Ψ∂ t f ( p ) = f ( t ) e − pt dt ∞ ∫ , p d Ψ dx − k + p c s ⎛⎝⎜⎜ ⎞⎠⎟⎟ Ψ = Ψ ( x , p ) Ψ ( k , x , p ) = U ∓ ( k , p ) e ± | k | α x , α = + p | k | c s U ∓ ( k , p ) (12) Clearly, can be interpreted as the amplitude of the spectral component of the displacement at the plane since (13) For the lower half-plane, a similar development is carried out and the displacement component in the lower half-plane can be shown to be (14) As before, has the interpretation as the amplitude of the spectral component of the displacement at the plane . The discontinuity in the displacement at the plane is (15) Thus, the spectral amplitude of the displacement discontinuity at is . The traction component of the stress tensor is continuous at the interface between the two half-planes. From Eq. (3), it is given by . (16) For a single spectral component, from Eqs. (12) and (14), we get (17) u ( x , x , p ) = U + ( k , p ) e − | k | α x e ikx . U + ( k , p ) x = + u ( x , x = + , p ) = U + ( k , p ) e ikx . u ( x , x , p ) = U − ( k , p ) e + | k | α x e ikx . U − ( k , p ) x = − x = δ ( x , p ) = u ( x , x = + , p ) − u ( x , x = − , p ) = ( U + ( k , p ) − U − ( k , p )) e ikx ≡ D ( k , p ) e ikx . x = D ( k , p ) x = τ ( x , p ) = σ ( x , x = ± , p ) τ ( x , p ) = − µ | k | α U + ( k , p ) e ikx = + µ | k | α U − ( k , p ) e ikx . hus, (18) where is the spectral amplitude of the shear stress on the interface. Thus, the spectral amplitudes of the shear stress at the interface and the displacement discontinuity are related by the equation (19) Multiplying throughout by and rearranging terms, the above equation can be rewritten as (20) In the limit , the coefficient of is clearly . This corresponds to the radiation damping effect, i.e., an instantaneous response in the slip velocity due to a shear stress disturbance. Extracting this response explicitly, Eq. (20) can be rewritten as (21) where , with . τ ( x , p ) = − µ | k | α U + − U − ) e ikx = − µ | k | α De ikx ≡ T ( k , p ) e ikx T ( k , p ) T ( k , p ) = − µ | k | α D ( k , p ). p µ p c s D ( k , p ) = − p | k | c s + p | k | c s ⎛⎝⎜⎜⎜⎜⎜⎜ ⎞⎠⎟⎟⎟⎟⎟⎟ T ( k , p ). p → ∞ T ( k , p ) − µ p c s D ( k , p ) + T ( k , p ) = F ( k , p ), F ( k , p ) = C ( k , p ) T ( k , p ) C ( k , p ) = − p | k | c s + p | k | c s ⎛⎝⎜⎜⎜⎜⎜⎜ ⎞⎠⎟⎟⎟⎟⎟⎟⎡⎣⎢⎢⎢⎢⎢⎢ ⎤⎦⎥⎥⎥⎥⎥⎥ ow, as , and hence is a well-behaved function of time. Taking inverse Laplace transform of the above equation, and noting that is the Laplace transform of the time derivative of , we get (22) where (23) The convolution kernel can be explicitly evaluated. Introducing the non-dimensional time variable, , it is easily shown that the convolution kernel , (24) where is the Bessel function of the first kind of order 1. This is to be contrasted with the result of Morrissey and Geubelle (1995) for the spectral formulation involving convolution over the displacement discontinuity at the interface. The convolution kernel they obtained is The present formulation can be adapted to incorporate initial stresses, , which may thought of as the stresses that would act at the plane in the absence of a displacement discontinuity there. In that case, we write the spectral amplitude of a single spectral component as (25) Modal Analysis:
In numerical simulations involving complex histories of the stress, and displacement discontinuity, , the time-convolution in Eq. (23) is carried out numerically. In order to C (., p ) → p → ∞ pD ( k , p ) D ( k , t ) µ ! D ( k , t )2 c s + T ( k , t ) = F ( k , t ), F ( k , t ) = C ( k , t − t ') T ( k , t ') dt ' t ∫ = | k | c s C (| k | c s ( t − t ')) T ( k , t ') dt ' t ∫ . C (.) Τ = | k | c s t C ( Τ ) = J ( Τ ) J (.) J ( Τ ) / Τ . τ o ( x , t ) x = T ( k , t ) τ ( x , t ) − τ o ( x , t ) = T ( k , t ) e ikx . τ ( x , t ), δ ( x , t ) erify the implementation of the convolution and to understand its numerical behavior, a modal analysis is carried out. A displacement discontinuity in a single Fourier mode of the form (26) where is the Heaviside step function, is applied and the time response of the resulting stress amplitude, , is determined both analytically and numerically. Taking Laplace transform of Eq. (26), it is clear that for the applied loading the spectral amplitude (27) Substituting in Eq. (19), we get (28) Taking inverse Laplace transform, it is easily seen that (29) where . The same quantity, , can be calculated numerically from Eq. (22). Substituting Eq. (26) in Eq. (22), we get . (30) where . The above equation is a Volterra equation of the second kind which can be solved numerically. The numerical solution is shown for different values of the ! δ ( x , t ) = H ( t ) e ikx , H ( t ) T ( t ) D ( k , p ) = p . T ( k , p ) = − µ | k |2 p + p | k | c s . r ( Τ ) ≡ − T ( Τ ) µ / 2 c s = + J ( Τ ') Τ ' ( Τ − Τ ') Τ ∫ d Τ ' Τ = | k | c s t r ( Τ ) + C ( Τ − Τ ' Τ ∫ ) r ( Τ ') d Τ ' = r ( Τ ) C ( Τ ) = J ( Τ ) iscretization parameter . Fig. 2 shows the comparison of the analytical solution given by Eq. (29) with the numerical solution of Eq. (30). We see good agreement between the analytical solution and the numerical solution when the discretization parameter Thus, the numerical implementation of the convolution in Eq. (22) is validated. Impulsive loading of two half-planes:
Consider a frictionless interface between two half-planes that is loaded by a pair of opposing impulsive line loads at the origin acting in the - direction. The impulsive load can be considered to give rise to the initial shear stress on the interface given by (31) where is the magnitude of the load and is the Dirac delta function (not to be confused with the same notation for the slip). The traction-free condition is (32) The interface between the two half-planes is discretized by a grid of elements. The size of each element is . Substituting Eq. (31) and Eq. (32) into Eq. (25), the spectral amplitude is determined by a fast Fourier transform (FFT) operation. It is then substituted into Eq. (22) and is determined numerically. An inverse FFT is performed on to determine the slip rate, , and then by numerical integration, the slip, The time step, , for the numerical integration of Eq.(22) is chosen so that it is a small fraction of the time take for the shear wave to traverse on grid spacing, , where is the Courant parameter. For the chosen discretization, the highest wavenumber is . We saw in the modal analysis that for an accurate numerical solution of the modal equation, ΔΤ = | k | c s Δ t ΔΤ = x τ o ( x , t ) = P δ ( x ) δ ( t ) P δ (.) τ ( x , t ) = N Δ x T ( k , t ) ! D ( k , t ) ! D ( k , t ) ! δ ( x , t ) δ ( x , t ). Δ t Δ t = β Δ x / c s β k max = π / Δ x he discretization parameter has to be at least 0.1. Applying this condition for the highest wavenumber , it is clear that we need . In other words, we need In practice, we are able to use much higher values of without introducing significant numerical errors. The slip at a location is plotted in Fig. 3 and compared with the analytical solution. The chosen parameters are and . We see good agreement between the analytical and numerical solutions. The oscillations in the numerical solution (which are not present in the analytical solution) can be ascribed to the limitation imposed on the spectral representation by a Gibbs-type phenomenon due to the Dirac delta function loading. Simulation of dynamic rupture propagation with a slip-weakening friction law:
To illustrate the use of the spectral numerical method proposed here, dynamic rupture propagation at an interface between two half-planes is simulated with a slip-weakening friction law acting at the interface between the half-planes. The frictional strength is given by (33) where and are the peak and residual frictional strength and is the critical slip required for attaining the residual frictional strength. An interface of length is considered with high-strength barriers of length at the left and right edges, to stop the rupture. It is discretized into elements with grid size . The normal stress is uniform throughout the length . The initial shear stress, , is shown in Fig. 4. It has a uniform value, , except in a central portion of length where the ΔΤ = | k | c s Δ t k max πβ ≤ β ≤ β x = X N = β = τ f ( δ ) = τ s − ( τ s − τ r ) δ / δ c , δ < δ c τ r , δ ≥ δ c ⎧⎨⎪⎩⎪ τ s τ r δ c L L s N Δ x = L / N L τ o τ bgo L nuc hear stress, , is above the peak residual frictional strength, . The frictional strength in the barriers at the left and right edges is set to a very high value. The simulation parameters are summarized in Table 1. The initial slip, slip velocity and shear stress are taken to be zero. The following algorithm is then implemented: 1) Perform FFT on the shear stress difference at previous time step, , to determine the spectral amplitudes, . 2)
Evaluate the convolution term Eq. (23) using the time-history of . It must be noted that the value of at the current time step, , does not contribute to the integral in Eq. (23) since . Thus, we know at the current time step the value of the combination (34) 3)
Perform inverse FFT to evaluate the value at the current time step of the combination (35) 4)
Update the slip at current time step using slip velocity distribution from the previous time step: (36) 5)
Use the friction law Eq. (33) to update the strength, . 6)
Determine the slip velocity at current time step, , and shear stress at current time step, , by solving Eq. (35) simultaneously with the friction law, Eq. (33). τ nuco τ s τ ( x , t − Δ t ) − τ o T ( k , t − Δ t ) T T T ( k , t ) C (0) = J (0) = µ ! D ( k , t )2 c s + T ( k , t ) = F ( k , t ) µ ! δ ( x , t )2 c s + τ ( x , t ) − τ o = f ( x , t ) δ ( x , t ) ! δ ( x , t − Δ t ) δ ( x , t ) = δ ( x , t − Δ t ) + ! δ ( x , t − Δ t ) Δ t τ s ! δ ( x , t ) τ ( x , t ) . If the material has attained the critical slip, , set and it follows that . b.
If and , set and it follows that . c.
If and , set and it follows that . The results of the simulation are shown in Figs. 5 to 7. The slip, slip velocity and shear stress along the interface at intervals of one second of time from to are shown in the figures. The time history of slip, slip velocity and shear stress history at one point on the fault are also shown in Figs. 8 to 10. With the chosen values of and , the numerical results appear to be well-resolved. The results obtained are also in qualitative agreement with previous studies of slip rupture propagation nucleated by an asperity, such as, Hajarolasvadi and Elbanna (2017).
Discussion:
The spectral formulation introduced here is an alternative to the one of Morrissey and Geubelle (1997). It appears that the present formulation is more easily extended to an interface between dissimilar elastic solids, than the one of Morrissey and Geubelle (1997). Taking advantage of the continuity of the traction component of stress at the interface, it is easily seen that Eq. (21) generalizes to (37) where the superscript “+” implies properties and functions relevant to the upper half-plane and the superscript “-“ implies properties and functions relevant to the lower half-plane. δ ≥ δ c τ = τ r µ ! δ / 2 c s = τ o + f ( x , t ) − τ r ( x , t ) δ < δ c τ f > f ! δ = τ ( x , t ) = τ o + f ( x , t ) δ < δ c τ f < f τ = τ f µ ! δ / 2 c s = τ o + f ( x , t ) − τ f ( x , t ) t = t = N β pD ( k , p ) + c s + µ + + c s − µ − ⎡⎣⎢⎢ ⎤⎦⎥⎥ T ( k , p ) = c s + µ + C + ( k , p ) + c s − µ − C − ( k , p ) ⎡⎣⎢⎢ ⎤⎦⎥⎥ T ( k , p ), aking inverse Laplace transform of the above equation, we get the spectral form of the boundary integral equation for a bi-material interface as (38) Thus, the convolution kernel has a simple form. The generalization of the approach of Morrissey and Geubelle (1997) to a bi-material interface was done by Geubelle and Breitenfeld (1997). It may be noted that in their formulation, the convolution kernel could not be expressed analytically and had to be found by a numerical inverse Laplace transformation. In contrast, we are able to find a closed-form expression for the convolution kernel for the bi-material interface, as shown above. Conclusions:
A new spectral formulation of the antiplane boundary integral equation method has been developed. It is based on performing elastodynamic convolution on the time-history of the shear stress on the interface. In contrast, prior approaches have performed convolution on the time-history of the displacement discontinuity at the interface. The formulation has been validated by performing a modal analysis and comparing the numerical and analytical solutions. The response of the half-planes to opposing impulse line loads on the interface has also been studied numerically and good agreement has been seen with the analytical solution. Dynamic rupture propagation with a slip-weakening friction law has also been simulated and the method is seen to give well-resolved numerical results. Further, it has been shown that the formulation is easily extended to incorporate material dissimilarity across the interface. ∂ D ∂ t + c s + µ + + c s − µ − ⎡⎣⎢⎢ ⎤⎦⎥⎥ T ( k , t ) = c s + µ + | k | c s + J (| k | c s + ( t − t ')) + c s − µ − | k | c s − J (| k | c s − ( t − t ')) ⎡⎣⎢⎢ ⎤⎦⎥⎥ t ∫ T ( k , t ') dt '. cknowledgement: This work is supported by Science and Engineering Research Board (SERB), India under grant no. CRG/2019/006080.
Data Availability Statement:
Data available on request.
References:
Andrews, D. J., Ben-Zion, Y., 1995. Wrinkle-like slip pulse on a fault between different materials.
J. Geophys. Res.
J. Tribol.
Int. J. Fract.,
93, 13-38. Budiansky, B. and Rice, J.R., 1979. An integral equation for dynamic elastic response of an isolated 3-D crack,
Wave Motion , 1, 187-92. Cochard, A. and Madariaga, R., 1994. Dynamic faulting under rate-dependent friction,
PAGEOPH , 142, 419-445. Di Bartolomeo, M., Meziane, A., Massi, F., Baillet, L., Fregolent, A., 2010. Dynamic rupture at a frictional interface between dissimilar materials with asperities.
Tribol. Int.
43, 1620–1630. Geubelle, P. H. and Breitenfeld, M. S., 1997. Numerical analysis of dynamic debonding under anti-plane shear loading,
Int. J. Fract. , 85, pp. 265-282. Geubelle, P. H. and Rice, J. R., 1995. A spectral method for three-dimensional elastodynamic fracture problems,
J. Mech. Phys. Solids , 43, 1791-1824. ajarolasvadi, S. and Elbanna, A. E, 2017. A new hybrid numerical scheme for modelling elastodynamics in unbounded media with near-source heterogeneities,
Geophys. J. Int.,
Bull. Seismol. Soc. Amer.
87, 1267-1280. Kammer, D. S., Yastrebov, V. A., Spijker, P., Molinari, J. F., 2012. On the propagation of slip fronts at frictional interfaces.
Tribol. Lett.
48, 27–32. Kammer, D. S., Yastrebov, V. A., Anciaux, G., Molinari, J.-P., 2014. The existence of a critical length scale in regularised friction.
J. Mech. Phys. Solids
63, 40-50. Kostrov, B.V., 1966. Unsteady propagation of longitudinal shear cracks,
J. Appl. Math. Mech. , 30, 1241-1248. Ma, X., Hajarolasvadi, S., Albertini, G., Kammer, D. S, Elbanna, A. E., 2018. A hybrid finite element-spectral boundary integral approach: Applications to dynamic rupture modeling in unbounded domains.
Int. J. Numer. Anal. Methods Geomech.
Int. J. Num. Methods Engg ., 40, 1181-1196. Ranjith, K., 2015. Spectral formulation of the elastodynamic boundary integral equations for bi-material interfaces,
Int. J. Solids Struct ., 59, 29-36. Tonazzi, D., Massi, F., Culla, A., Baillet, L., Fregolent, A., Berthier, Y., 2013. Instability scenarios between elastic media under frictional contact.
Mech. Syst. Signal Process.
40, 754–766.
Parameter Symbol Value
Shear wave speed (km/s) 3,464 Density (kg/m ) 2670.00 Background shear stress (MPa) 70.00 Nucleation shear stress (MPa) 81.60 Nucleation zone size (km) 3.00 Interface length (km) 100.00 High-strength barrier size (km) 35.00 Number of elements 1024 Courant parameter 0.30 Peak frictional strength (MPa) 81.24 Residual frictional strength (MPa) 63.00 Critical slip-weakening distance, (m) 0.40 Table 1: Parameters used in the simulation of slip rupture propagation c s ρ τ bgo τ nuco L nuc L L s N β τ s τ r δ c Figure 1: Interface between two half-planes is located at . Displacement is in the direction. Discontinuity in displacement can develop at the interface in response to applied stresses at the boundary. x = x x x Figure 2: Comparison of analytical and numerical solutions of modal analysis for different values of the discretization parameter r( T ) T = |k|c s tAnalytical ∆ T = 0.10 ∆ T = 0.25 ∆ T = 0.50
Figure 3: Comparison of analytical and numerical solutions for problem of slip response to pair of opposing impulsive line loads at the frictionless interface -2-1 0 1 2 3 4 5 6 7 8 0 0.5 1 1.5 2 µ | X | δ / P c s c s t/XAnalyticalNumerical Figure 4: Initial shear stress and peak frictional strength along the interface are shown. In a region of 3 km at the centre, the initial shear stress is greater than the peak frictional strength in order to initiate a rupture
50 55 60 65 70 75 80 85 90-15 -10 -5 0 5 10 15Distance along interface x (km)Initial stress (MPa)Peak frictional strength (MPa) Figure 5: Slip along the interface at time instants separated by 1 second S li p ( m ) Position x (km) t = 1 st = 2 st = 3 st = 4 st = 5 s Figure 6: Slip velocity along the interface at time instants separated by 1 second -2 0 2 4 6 8 10 12 14-15 -10 -5 0 5 10 15 S li p v e l o c it y ( m / s ) Position x (km)t = 1 st = 2 st = 3 st = 4 st = 5 s Figure 7: Shear stress along the interface at time instants separated by 1 second
62 64 66 68 70 72 74 76 78 80 82-15 -10 -5 0 5 10 15 I n t e rf ac i a l s h ea r s t r e ss τ ( M P a ) Position x (km) t = 1 st = 2 st = 3 st = 4 st = 5 s Figure 8: Time history of slip at a point 4.5 km to the right of the centre S li p ( m ) Time t (s)
Figure 9: Time history of slip velocity at a point 4.5 km to the right of the centre -1 0 1 2 3 4 5 6 7 0 2 4 6 8 10 12 S li p v e l o c it y ( m / s ) Time t (s)
Figure 10: Time history of shear stress at a point 4.5 km to the right of the centre.
62 64 66 68 70 72 74 76 78 80 82 0 2 4 6 8 10 12 I n t e rf ac i a l s h ea r s t r e ss ( M P a ))