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],  The problem formulation outlined above is shown in figure 1.

(a) Governing equations and boundary conditions
The equilibrium condition takes the form of the Laplace equation where u j = u j (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: [[u]](x, 0) = u 1 (x, 0 + ) − u 2 (x, 0 − ) = kσ (1) yz (x, 0) and [[σ yz ]](x, 0) = σ (1) yz (x, 0 + ) − σ (2) yz (x, 0 − ) = 0, where σ (j) yz (x, y) = μ j (∂u j /∂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) yz , 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 σ (j) yz (x, 0) = p j (x). In summary, we have the following set of boundary conditions over the whole interface y = 0 [[u]](x, 0) = 0, x > 0, [[u]](x, 0) = kμ 1 ∂u 1 ∂y (x, 0), −L < x < 0, As a result, relations (2.4) 1 , (2.4) 3 and (2.4) 4 can be rewritten in the form and σ 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 l , U r , σ l , τ r , σ 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 l = kσ l and U r = kσ r .

(b) 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 j using the standard notationW 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: whereW − (t) is unknown. We introduce new functions P − , G, F, related to the applied load p(x − L): (2.12) 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 [[g]] − =g − 1 −g − 2 is called a jump, and g − = (g − 1 +g − 2 )/2 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.
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 L (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 L (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 L (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 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 Li 2 (z) is the dilogarithm, defined by (see for example [43]) According to [43], the imaginary part of the dilogarithm Li 2 (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 √ t 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.

(a) Iterative procedure
Following the approach outlined in [38], we approximate Z + , R − (2.13) by functions Z + n * , R − n * , which are determined using the recurrence relations where n = 1, n * and The iterations start with the initial condition 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] h 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 | 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 t 1/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.
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 n , which represents the accuracy of the solution to the Wiener-Hopf equation (2.14)  2 3 4 5 6 7 8 9 10 11 12 13 14 15  10  where F n 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 L 2 norm on R by · 2 . Additionally, in these figures, we also show the distributions of three measures: the quantity δF n mentioned above, and the norms of difference between the consecutive iterative steps R − n+1 − R − n 2 , Z + n+1 − Z + n 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 + n , R − n 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

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.
(a) Evaluation of the fracture mechanics parameters  and analyse its dimensionless value K L III /K III , where K III is the stress intensity factor in the absence of the damage zone (L = 0), taken from [42]. Next, we evaluate the parameter σ l , 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 l = kσ l , 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 L III /K III and u * l stabilize throughout the iterative process for various L are provided in figure 5.
Meanwhile, to demonstrate the numerical error in approximating the fracture parameters, we introduce the following measures δK (n) 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 L III /K III = 1 and u * l = 0 when L = 0.  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.

(b) 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 j (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 n (x) is the function W(x) at the nth iteration. The values of this measure, alongside the corresponding dimensionless solution for W n (x)/l * along the process zone, are provided in figure 8. It can be clearly seen that after the second iteration the curves W n (x) almost coincide.
With the rapid convergence of W n (x) established, we approximate the function W(x) with W n * (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.
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.

(c) 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 C , K C are material parameters. By normalizing (4.6), these equilibrium conditions can be placed in terms of u * l and K L III /K III u * l < U * C , where the relationship between u * l and K L III /K III 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 * C and fracture toughness K C /K III 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.
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 * l = U * C , while the corresponding value of L * can obtained from figure 7, which we shall denote by L * A . Furthermore, it is clear that condition (4.7) 1 is satisfied for any L * < L * A . Similarly, point B defines the critical process zone length L * B , with condition (4.7) 2 holding for any L * > L * B . It should be noted, as seen in figure 10, that the values of (L * B , L * A ) obtained from each condition are identical To summarize, the region (L * B , L * A ) 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.

(d) 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 Figure 11. Stress intensity factor, crack opening and critical stress in the process zone for various d * .
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 * .
It is readily apparent from the figure above that as the constant k increases so do the value of K L III and u l , while the critical stress σ l 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 σ l (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 ( 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 III , at the crack tip and the critical opening, u l , 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 III and u l 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.
Data accessibility. This article has no additional data. Competing interests. We declare we have no competing interests.