Numerical factorization of a matrix-function with exponential factors in an anti-plane problem for a crack with process zone.

In this paper, we consider an interface mode III crack with a process zone located in front of the fracture tip. The zone is described by imperfect transmission conditions. After application of the Fourier transform, the original problem is reduced to a vectorial Wiener-Hopf equation whose kernel contains oscillatory factors. We perform the factorization numerically using an iterative algorithm and discuss convergence of the method depending on the problem parameters. In the analysis of the solution, special attention is paid to its behaviour near both ends of the process zone. Qualitative analysis was performed to determine admissible values of the process zone's length for which equilibrium cracks exist. This article is part of the theme issue 'Modelling of dynamic phenomena and localization in structured media (part 1)'.

Introduction

The theory of brittle fracture, developed by Griffith using the energy balance approach [1] and later reformulated by Irwin in terms of stress fields, implies the existence of a stress singularity at the tip of a sharp crack. It is also widely recognized that, during the fracture process, phenomena occurring close to the crack tip are closely connected to the microstructure of the material. Here, we recall the models of Irwin [2] and Orowan [3] that account for plastic behaviour near the crack tip. In [4], Neuber investigates the mechanisms of physically nonlinear stress concentrations by replacing the sharp crack with a blunt notch. Barenblatt introduced the cohesive zone model in [5,6], describing the resisting forces that occur when material elements are being pulled apart. This approach allows the elimination of the stress singularity at the crack tip, which is useful for classic finite-element modelling [7]. Cohesive zone models have been widely used to study many physical phenomena, such as crack growth in viscoelastic materials [8], debonding in composite materials [9], the fracture of adhesive joints [10], among others.

All models of fracture belong to the class of mixed boundary value problems which, in turn, can be reduced to Wiener–Hopf functional equations [11-15]. This applies to modelling cracks through a formulation that is continuous, for both static and steady-state approaches [16,17] and for dynamic problems of fractures propagating within a discrete structure. The seminal contribution made by Slepyan to modelling numerous fracture problems should be mentioned. Among others, he has developed a unique technique for tackling dynamic problems in discrete structures [16,18]. This method is extremely effective for considering both lattice structures composed of masses and connecting springs [19], and structures made of masses and beams [20,21]. The developed technique allows one to determine the fundamental properties of the solution, in particular those related to the nature of crack propagation or phase transitions, and to apply this knowledge to a wide array of applications [22]. However, most of the problems solved so far only use a scalar Wiener–Hopf equation. This is because more advanced applications result in coupled systems, taking the form of a vectorial Wiener–Hopf problem with matrix kernels [23,24]. The same issue of complicated boundary value problems leading to vectorial Wiener–Hopf-type equations is also encountered during the study of certain classes of contact problems [25,26].

The issue of multiplicative decomposition, or factorization, is one of the main stages of the Wiener–Hopf method. In the case of a scalar problem, the solution can be found explicitly [27]. For higher-order problems, however, constructive solutions are only known for specific classes of matrices. Notable techniques include efficient approximation using rational functions [28-30], while numerical factorization of generalized Khrapkov–Daniele matrices is discussed in [24]. A comprehensive review of matrix factorization techniques can be found in [31].

In this paper, we analyse the static loading of a semi-infinite interface crack between two dissimilar elastic materials. We consider a formulation where the process zone reflects the bridging effect along a finite part of the interface in front of the crack. The contact in this zone is modelled by the so-called weak imperfect interface [32] which describes soft thin adhesive joints [33]. The Wiener–Hopf kernel corresponding to this problem is a matrix-valued function containing oscillating terms. There are several approaches to matrix factorization available in this case (see [34-37]); however, for this particular problem the most appropriate is that proposed in [38]. We discuss the numerical algorithm's peculiarities and demonstrate its efficiency. Finally, we determine all fracture mechanics parameters, evaluate a condition for the existence of an equilibrium state of the crack under remote loading and compute the corresponding length of the process zone.

Problem statement

Let us consider an infinite plane occupied by two different linearly elastic and isotropic materials with shear moduli μ, j = 1, 2. These materials are joined along a linear interface. We introduce the coordinate system (x, y), where the x-axis is directed along that interface, the y-axis is directed towards the first material and the origin is located at the point that separates the intact and damaged regions. On the interface y = 0, we place the process zone, −L < x < 0, which separates the cracked area, x < − L, from the region of ideal contact between the materials, x > 0, which is ahead of the crack tip. Here, L is the length of the process zone mentioned above. An external out-of-plane load p(x) is applied at the crack faces, between some predefined points x = − b and x = − a, with a, b > 0, such that it is self-balanced in terms of the principal force The problem formulation outlined above is shown in figure 1.

Figure 1.

Mode III crack with process zone. (Online version in colour.)

Governing equations and boundary conditions

The equilibrium condition takes the form of the Laplace equation where u = u(x, y) is the out-of-plane component of the displacement in the jth material. From [32], we have that the transmission conditions in the process zone y = 0,  − L < x < 0 are given by the following relations: where σ((x, y) = μ(∂u/∂y)(x, y) is the shear stress in the jth material and k is the interface parameter (which is similar to a spring constant [39]). Note that, in contrast to Barenblatt's model where the cohesive forces are predefined a priori while the length of the respective zone is computed to eliminate the stress tip singularity, our model defines only a relationship between the jump of displacements, [[u]], and traction, σ(1), along the interface. Meanwhile, the length of the process zone is determined from additional fracture conditions related to the strength of this zone.

Along the ideal contact region, y = 0, x > 0, there is a continuity of displacements and tractions. At the crack surface y = 0, x < − L, we have the applied external load σ((x, 0) = p(x). In summary, we have the following set of boundary conditions over the whole interface y = 0 It is convenient to introduce new functions W,  V,  Φ,  Ψ and g As a result, relations (2.4)1, (2.4)3 and (2.4)4 can be rewritten in the form where H(x) is the Heaviside function.

Following a similar approach to that outlined in [40], we look for a solution with the following asymptotic behaviour of the displacements and tractions where U, U, σ, τ, σ0 are some unknown constants. In the following, it will be shown that conditions above lead to a unique and physically justified solution. From (2.7), using the notation introduced in (2.5), we obtain The problem is now formulated in terms of the equilibrium equation (2.2) and boundary conditions (2.4)2, (2.6). For further analysis, it is also important to take into account the asymptotic relations (2.7) and (2.8). We note, according to (2.3), that U = kσ and U = kσ.

Reduction to a vectorial Wiener–Hopf problem

To simplify later analysis, it is convenient to introduce the following constants: We denote the one-sided Fourier transforms of the functions W, V, Φ, Ψ, g using the standard notation Applying the Fourier transform to (2.2) and (2.4)–(2.6), with respect to the variable x, we find that the image of displacement is given by the following expression: where is unknown.

We introduce new functions P−, G, F, related to the applied load p(x − L): alongside new unknown functions Z+, R−: In (2.10)–(2.13), the tildes (∼) correspond to Fourier transforms, while the ‘ ± ’ superscript denotes that the transform is regular when ± Im(t) > 0. Here is called a jump, and is the average. Note that the auxiliary parameter γ is introduced with the sole limitation that Re(γ) > 0, in order to eliminate the possible singularity of the auxiliary function Z+(t) at zero.

The function is regular in the half-plane ∓Im(t) > 0 and has a branch cut {t | Re(t) = 0,  ± Im(t) > 0}. It is also important to mention that the asymptotic behaviour of can be obtained from (2.7) and (2.8) by means of Abel-type theorems [11]. Indeed, functions decay as t−1 at infinity, while decays like t−1/2 (in the corresponding half-planes). At the zero point, t = 0, we find that the functions are bounded, while function has a singularity of order t−1/2.

In addition, the original problem transforms to a matrix Wiener–Hopf equation along the real axis: where and The total index of this matrix, κ = κ1 + κ2, is zero since det M(t) = const [41]. In general, partial indices κ1, κ2 may be non-zero, which affects the uniqueness of the solution and requires the fulfilment of additional conditions for the function f(t). The two simplest cases L = 0 (see [42]) and L → ∞ (for example, [32,40]) both yield a unique solution. Therefore, one can expect that matrix M(t) will admit a stable factorization for an arbitrary finite L≠0.

Solution of the matrix Wiener–Hopf equation

Let us introduce function D(t) defined on the real axis and perform both the additive and multiplicative splittings where and Here ln ± (t) is a branch of log t along the imaginary axis to be regular in the half-plane ± Im(t) > 0 and Li2(z) is the dilogarithm, defined by (see for example [43]) According to [43], the imaginary part of the dilogarithm Li2(1 − t) has a jump across the negative half of the real axis Im(t) = 0 while the real part is continuous in the entire complex plane.

In addition to this, we observe that Q ± (t)∼1 at infinity, while it decays as when t approaches zero from the corresponding half-plane. Similarly, functions Q ± (t) in (3.3)2 give a closed form factorization of the function Q(t) = |t|(d + |t|)−1, representing the Wiener–Hopf kernel of the equation (2.14) in the case L → ∞. It is worth noting that the case of an infinite plane with a semi-infinite interfacial crack was considered in [33], where a function similar to 1/Q(t) was factorized numerically.

Iterative procedure

Following the approach outlined in [38], we approximate Z+, R− (2.13) by functions Z+, R−, which are determined using the recurrence relations where and The iterations start with the initial condition and n* is the required number of iterations. Note that in (3.6) the functions inside the large brackets are continuous along the real axis, while they decay at infinity as t−1/2 and t−1 for (3.6)1 and (3.6)2, respectively. The convergence of this procedure has been discussed in [38]. As soon as Z+ and R− are evaluated, the calculation of and becomes straightforward. Indeed, from (2.14) we derive a scalar Wiener–Hopf equation with respect to the unknown functions , : which can be solved analytically. In general, for a given function h(t) that decays with the growth of |t|, the decomposition h(t) = h+(t) + h−(t) is accomplished by means of Cauchy-type integrals. Their limiting values along the real axis are defined by the Sokhotsky–Plemelj formulae [27]

Convergence

To provide an illustrative numerical example, we consider a symmetric uniform loading of unit magnitude, distributed along the line segment ( − b, − a) (figure 1). We normalize the damage zone length by a parameter l* = 10(b − a) such that is a dimensionless parameter.

For this formulation, the distributions of errors |R−(t) − R−(t)|, |Z+(t) − Z+(t)| along the real axis are shown in figure 2. It is clear that the numerical error decreases significantly as the iterative process continues. Note that the highest error occurs near t = 0, even though the functions R−, Z+ are bounded at the zero point. This is because the next term in their asymptotic expansion is of the order t1/2, and as such calculating the coefficients corresponding to these terms yields computational errors, which result in the noticeable inaccuracies at t = 0. It should be noted, however, that this error does decrease as the iterative process continues.

Figure 2.

Reduction of the error with growing number of iterations for L* = 0.5. (Online version in colour.)

Going further, in figures 3 and 4, we show that the iterative method's rate of convergence increases with the growth of the process zone L*. To demonstrate this we introduce the discrepancy δF, which represents the accuracy of the solution to the Wiener–Hopf equation (2.14) where F is the left-hand side of the second component of the vectorial Wiener–Hopf equation (2.14), while F is given in (2.12). In (3.12) and below, we denote the norm on by ∥ · ∥2.

Figure 3.

The discrepancy of the equation and the norms of errors versus the number of iterations for various numbers of integration points N. (Online version in colour.)

Figure 4.

The discrepancy of the equation and the norms of errors versus the number of iterations for a fixed number of integration points N = 800. (Online version in colour.)

Additionally, in these figures, we also show the distributions of three measures: the quantity δF mentioned above, and the norms of difference between the consecutive iterative steps ∥R− − R−∥2, ∥Z+ − Z+∥2, to provide a fuller picture of the numerical method's rate of convergence.

One can observe in figure 3 that the iterative method (3.6) converges faster (in terms of the number of iterations) and gives more accurate results for higher values of L*. However, it is important to account for the fact that functions Z+(t) and R−(t) oscillate with frequencies that grow as we increase L* and, consequently, a larger number of mesh nodes N are needed to account for this. This, in turn, makes each individual iteration more time-consuming.

To further clarify this point, and the dependence of the final algorithm on the number of integration points N, in figure 4 the same graphs are provided for fixed N = 800 with differing values of L*. It is clear that the number of iterations needed to reach the saturation limit decreases with increasing L*. On the other hand, for a fixed value of N, the numerical method gives more accurate results for smaller values of L*. The latter issue is eliminated by taking the number of mesh nodes N as a function of L*, as demonstrated by the results in figure 3.

It is worth noting that calculating the functions Z+, R− at each step of the iterative process (3.6) only involves evaluating the singular integrals within (3.10), while K ±  and T ±  remain the same for every iteration step. The integration in (3.6) was performed numerically on a finite part of the real axis, while taking into account the asymptotic behaviour of the integrand at infinity. Finally, to increase the algorithm's efficiency, the integration interval was partitioned non-uniformly.

In table 1, we show the time costs for different computational scenarios. The calculations were performed on a computer with an Intel Core i3-2310M CPU @ 2.10 GHz  ×  2 processor.

Table 1.

Time costs.

	N = 800			N = 1100	N = 1600	N = 5000
	L* = 0.1	L* = 0.5	L* = 2	L* = 0.1	L* = 0.5	L* = 2
number of iterations before saturation	11	8	6	12	10	8
average time per iteration (s)	29	29	29	50	70	320
total time (s)	330	249	195	650	784	2685
error δFn	7 × 10−5	7 × 10−5	10−4	10−5	5 × 10−6	4 × 10−6

Analysis of the numerical results

With the iterative scheme for solving the vectorial Wiener–Hopf problem in place, and the accuracy of the numerical algorithm demonstrated, we can begin the evaluation and analysis of the system behaviour. We begin with an examination of the fracture mechanics parameters, before moving on to the process zone itself, and finally determining the fracture's equilibrium conditions.

Evaluation of the fracture mechanics parameters

The known asymptotics of functions W, Φ, V, Ψ in the limits as x → 0− and x →  − L (2.7)–(2.8) allow us to define the behaviour of the stresses and displacements at both ends of the process zone. For example, the behaviour of the stress functions as x → 0+ follows immediately from (2.6) as: . Noting this, linear elastic fracture mechanics (LEFM) allows us to immediately introduce the stress intensity factor and analyse its dimensionless value K/K, where K is the stress intensity factor in the absence of the damage zone (L = 0), taken from [42]. Next, we evaluate the parameter σ, that is the stress level at the left edge of the process zone (2.7). Then, combining the relationship between the jump of displacements and traction (2.3) in this zone with the known asymptotics yields U = kσ, from which we can obtain the displacement of the crack opening at x = − L in its dimensionless form The rates at which the values of K/K and u* stabilize throughout the iterative process for various L are provided in figure 5.

Figure 5.

Convergence of the values of LEFM parameters with the iterations. (Online version in colour.)

Meanwhile, to demonstrate the numerical error in approximating the fracture parameters, we introduce the following measures The rate at which these errors decrease throughout the iterative process, alongside the dependence of the fracture parameters on the length of the process zone, are shown in figures 6 and 7. It should be noted that K/K = 1 and u* = 0 when L = 0.

Figure 6.

Normalized stress intensity factor. (Online version in colour.)

Figure 7.

Normalized crack opening at x = − L. (Online version in colour.)

It can easily be seen in figures 5–7 that the outlined method is able to achieve a high accuracy of approximation for the stress intensity factor and crack opening displacement. This is achieved after only a small number of iterations, and with a high level of stability of the numerical algorithm.

Note that figures 6 and 7 establish implicit relationships between the length of the process zone and two fracture parameters under consideration. They will be used later to identify admissible length of the process zone for a stable crack, see §4c.

Displacements inside the fracture zone

To conduct an examination of the process zone, y = 0, − L < x < 0, we recall from (2.6) that W(x) corresponds to the jump of displacements along the interface [[u]](x, 0) = W(x)(1 − H(x)). In fact, the displacements u(x, 0), (j = 1, 2) can be obtained in terms of this function by applying the inverse Fourier transform to (2.11), while for the traction we simply use the transmission conditions (2.3)

To demonstrate the solution convergence for the function W(x), we consider the case[1] when L* = 0.5. Similarly to the previous subsection, we introduce the new measure where W(x) is the function W(x) at the nth iteration. The values of this measure, alongside the corresponding dimensionless solution for W(x)/l* along the process zone, are provided in figure 8. It can be clearly seen that after the second iteration the curves W(x) almost coincide.

Figure 8.

Convergence of W(x) for L* = 0.5.

With the rapid convergence of W(x) established, we approximate the function W(x) with W(x), where n* is the number of required iterations to achieve a desired level of accuracy. Using this approach, the solution for function W(x), for different values of L*, are depicted in figure 9.

Figure 9.

Jump of displacements for various L*.

Therefore, as we can see from (4.4), an understanding of the system's behaviour within the process zone can be obtained simply from the analysis of W(x). In particular, one can observe, that both the displacement's discontinuity and the shear stress increase in the direction from the crack tip towards the opposite end of this zone.

Equilibrium conditions

In the previous analyses, we assumed that the length of the process zone L was predefined, but in reality, for equilibrium cracks, there are two critical conditions which ensure that the fracture will not propagate and that neither end of the process zone will move. These conditions can be expressed in terms of the crack opening displacement and stress intensity factor as follows: where U, K are material parameters. By normalizing (4.6), these equilibrium conditions can be placed in terms of u* and KIII/KIII where In this form, the equilibrium conditions can be related back to the normalized cohesive zone length L* using observations from §4a. Combining the results displayed in figures 6 and 7, the relationship between u* and KIII/KIII is immediately obtained, as provided in the top left of figure 10. This line represents the boundary between unstable and equilibrium cracks. Indeed, provided critical values of the crack opening displacement U* and fracture toughness K/KIII are taken from the white region (see the red marker in figure 10), conditions (4.7) are automatically satisfied. From this, alongside previously obtained relations, admissible values of L* for equilibrium cracks can be defined.

Figure 10.

Admissible values of L*. (Online version in colour.)

To better illustrate this point, consider taking two arbitrary points A,  B along the boundary between the unstable and equilibrium regimes (figure 10). Then, by definition, at the point A we have u* = U*, while the corresponding value of L* can obtained from figure 7, which we shall denote by L*. Furthermore, it is clear that condition (4.7)1 is satisfied for any L* < L*. Similarly, point B defines the critical process zone length L*, with condition (4.7)2 holding for any L* > L*. It should be noted, as seen in figure 10, that the values of (L*,  L*) obtained from each condition are identical

To summarize, the region (L*,  L*) obtained using the approach outlined above defines the range of the normalized length of the damage zone, L*, over which the equilibrium state (4.6) is possible. Furthermore, if the critical fracture parameters for a given problem fall within the grey region in the top left of figure 10, then for any L* at least one condition in (4.7) fails and the crack becomes unstable.

Impact of the damage zone's stiffness

To conclude the analysis of the given problem, we examine how the parameter of the process zone, k, affects other fracture parameters. This is achieved through an examination of the normalized constant d*, defined as where constant l*, which was introduced to ensure a proper scaling, was previously determined in (3.11), while the parameter d was introduced in (2.9). Note that all of the results presented in the previous subsections were obtained for a fixed value of d* = 0.5. The relationships between the normalized process zone length and the stress intensity factor, crack opening and critical process zone stress are provided in figure 11 for varying values of the parameter d*.

Figure 11.

Stress intensity factor, crack opening and critical stress in the process zone for various d*.

It is readily apparent from the figure above that as the constant k increases so do the value of KIII and u, while the critical stress σ decays. In other words, an increase in the parameter k leads to a weaker response to the crack within the process zone, while for smaller k this zone demonstrates tougher behaviour as reflected in the growth of σ (in the right part of figure 11). Note that, in the limit as k → 0, the solution is regular at x = 0 while the crack tip and stress singularity are located at the opposite edge of the damaged zone x = − L.

Discussions and conclusion

In this paper, we considered the static loading of a mode III crack at the interface of two elastic materials, with a process zone modelled by relations (2.3). The original problem, given in terms of equation (2.2) and boundary conditions (2.4)–(2.6), was transformed to a matrix Wiener–Hopf equation (2.14). The main difficulty in solving this equation was the presence of exponential factors inside matrix M.

It has been demonstrated that an iterative approach, of the kind first proposed in [38], is highly effective at resolving this Wiener–Hopf equation, thereby both solving the initial problem and determining important fracture mechanics parameters for the system. Namely, the stress intensity factor, K, at the crack tip and the critical opening, u, at the opposite end of the process zone. A thorough analysis of the numerical approach used showed the high level of stability and accuracy, alongside the rapid rate of convergence, for the entire solution. Of particular note was the fact that the values of both K and u stabilize after only a small number of iterations, although the required number of iterations is higher for smaller lengths of the process zone L.

The fact that the method outlined here is more efficient for larger L is not unexpected, as the initial condition (3.8) of the iterative process corresponds to the case L → ∞. As such, although the iterative method does converge for all positive L > 0, it is clear from the presented results that both the solution accuracy and convergence rate improve with the growth of L. There are however solutions to this issue. The computational performance, in this case, could be greatly increased through a rescaling of the x-axis, although this would introduce an additional challenge in the manipulation of the right-hand side of (2.14). Alternatively, for small L there exists an alternative asymptotic method for the matrix factorization, which was presented in [35].

In addition, it was demonstrated that the relationship between critical values of the crack opening displacement and the stress intensity factor represent a boundary between unstable and equilibrium cracks. This result was achieved through qualitative analysis, performed to determine admissible values of length of the process zone (4.9) for which both necessary conditions for equilibrium in (4.6) are satisfied.

Meanwhile, the impact of the stiffness parameter, 1/k, on the solution was shown to become more pronounced as k increases, with the crack opening growing while the traction inside the damage zone decays. Conversely, reducing k results in a higher resistance to the crack opening within the process zone. Indeed, for a significantly stiff process zone (k → 0), the stress singularity at the right edge vanishes and a high-stress concentration appears near the point x = − L.

Finally, the model considered in this paper can be used to analyse steady-state propagation of an interface crack with a developed damage zone. Similar models with a semi-infinite crack and various imperfect material interfaces were considered in [32]. In the former problem, the crack speed can influence the size of the damage zone and consequently affect the stability of the crack propagation.