Crack tip solution for Mode III cracks in spring interfaces

Engineering Fracture Mechanics 288 (2023) 109293 Available online 4 May 2023 0013-7944/© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Contents lists available at ScienceDirect Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech Crack tip solution for Mode III cracks in spring interfaces Sara Jiménez-Alfaroa,b, Vladislav Mantičb,∗ aInstitut Jean Le Rond d’Alembert, Sorbonne Université, Pl. Jussieu, 4, Paris, 75005, France bGrupo de Elasticidad y Resistencia de Materiales, Escuela Técnica Superior de Ingenieria, Universidad de Sevilla, Camino de los Descubrimientos, Sevilla, 41092, Spain ARTICLE INFO Keywords: Antiplane strain Adhesive interface Interface fracture Asymptotic solution Logarithmic stress singularity ABSTRACT Considering an infinite linear elastic isotropic solid in antiplane shear, the Mode III crack-tip solution for a semi-infinite crack located in a straight spring interface is systematically studied for the first time. A new analytic expression for this crack-tip solution is given in the form of a double asymptotic series of the main and the so-called associated shadow terms. It is shown that the series of the shadow terms associated with a main term is infinite, and all shadow terms include logarithmic terms. Thus, although the interface tractions are bounded, the linear elastic solution at this crack-tip has a logarithmic stress singularity which is comprehensively analysed. Noteworthy, the character of this stress singularity is very different from the well-known square root singularity at the crack tip in the classical fracture mechanics. A key advantage of the present approach is its simplicity, as only elementary mathematical tools are employed, and also its easy implementation in a computer algebra software. The latter fact is very relevant because the expressions of higher-order shadow terms become increasingly complicated, so their generation by a computer code becomes crucial. The present results allow the implementation of new enriched or singular crack-tip finite elements for such cracks, and the automatic generation of analytic solutions for benchmark problems for testing the finite-element codes using these special elements. Such codes can be applied to efficient numerical modelling of interface cracks, e.g., in adhesively bonded joints with a thin adhesive layer. 1. Introduction The classical Winkler spring interface model [1], which considers a continuous distribution of independent linear elastic springs, is widely used to solve many engineering problems due to simplicity of analytical studies and for the sake of computational efficiency in numerical models. In particular, thin adhesive layers are often modelled using the related shear-lag model [2] and its generalizations [3,4]. A recent review of the spring model can be found in [5], and comprehensive classifications of interface models including the spring model in [6,7]. Design and analysis of adhesively bonded joints working under different loading conditions requires suitable computational models for crack onset and propagation in adhesive layers [8]. In this context, although several models on different scales exist for the study of crack initiation and growth in adhesive layers, one of the most widely used in analytical and computational studies is the spring interface, also referred to as weak, imperfect, linear elastic or Robin interface, see [9–14]. It is well known that according to a general mathematical analysis of singularities in elliptic partial differential equations developed in [15,16], the displacements and stresses around a singular point can be described analytically through an asymptotic series expansion, including power-logarithmic terms in the radial variable. ∗Corresponding author. E-mail address: [email protected] (V. Mantič). https://doi.org/10.1016/j.engfracmech.2023.109293 Received 13 January 2023; Received in revised form 18 April 2023; Accepted 22 April 2023 Engineering Fracture Mechanics 288 (2023) 109293 2 S. Jiménez-Alfaro and V. Mantič Nomenclature 𝑎(𝑙) 𝑗,𝑘, 𝑏(𝑙) 𝑗,𝑘 Constant coefficients in the asymptotic series 𝐁𝑗,𝑘 Vector of coefficients in the system of equations 𝑒𝑗Relative error in the Robin boundary condition 𝑘Stiffness of the linear-elastic spring distribution 𝐾𝑗Generalized Stress Intensity Factors (GSIFs) 𝜇Shear modulus 𝐺Energy Release Rate (ERR) 𝐆𝑗,𝑘 Right-hand side vector of the system of equations 𝐌𝑗,𝑘 Matrix of the system of equations 𝐿Characteristic length in the studied problem 𝑁Number of main terms considered in the asymptotic series 𝑆𝑗Number of shadow terms considered for the 𝑗th main term 𝑢𝑒, 𝑢𝑜Even and odd part of the out-of-plane displacement solution 𝑢𝑗𝑗th element of the asymptotic series 𝑢𝑘 𝑗𝑘th shadow term associated to the 𝑗th main term 𝑥, 𝑦 Cartesian coordinates 𝛿Virtual crack length increment 𝛾Governing parameter of the problem 𝛤1,𝛤2Boundary parts 𝜆𝑗Singularity exponent 𝛺+,𝛺−Upper and lower half domain 𝑟, 𝜃 Polar coordinates 𝜎𝑟𝑧, 𝜎𝜃𝑧 Stress components in the cylindrical coordinates system BVP Boundary Value Problem ERR Energy Release Rate FEM Finite Element Method GSIF Generalized Stress Intensity Factor LEFM Linear Elastic Fracture Mechanics N–N Neumann–Neumann problem N–R Neumann–Robin problem R–N Robin–Neumann problem Considering a particular case of a crack in a thin adhesive layer modelled by a spring interface, the question about the character of the asymptotic crack-tip solution arises, which seems to be quite different from that of a crack in a homogeneous linear elastic material described by the classical Linear Elastic Fracture Mechanics (LEFM). As it can be seen in the asymptotic series for the crack tip solution in LEFM, deduced in [17,18] for plane problems (Mode I and II), and in [19] for antiplane shear problems (Mode III), stresses near the crack tip have singular solutions, with square root singularities, multiplied by Stress Intensity Factors (SIFs). In particular, the asymptotic series for the Mode III crack-tip solution can be easily obtained by separation of variables applied to the displacement solution in polar coordinates, as shown, e.g., in [20]. The behaviour of stresses and displacements around the tip of a crack situated in a spring interface has not been fully characterized so far, although some relevant results have been obtained by several authors [10–13,21–25], indicating that a logarithmic stress singularity takes place at such a crack-tip. Especially significant for the present work are the contributions [22,23], for the antiplane elastic problem, and [24], for the plane elastic problem, who deduced the first terms of the asymptotic expansions at the tip of a crack situated at a spring type interface. While an approach based on the well-known Wieghardt–Williams procedure was used in [22], more advanced mathematical tools such as singular integral equations and Mellin transformations were used in [23,24]. Regarding the related asymptotic series for corners (angular sectors) of any interior angle with spring type boundary condition, several relevant results were obtained by various authors. Truncated asymptotic series including also logarithmic terms were introduced in [22,26–28] to approximate the linear elastic solutions in corners of various specific angles with spring boundary condition at one or both faces, for antiplane and plane elasticity. However, the number of terms in these asymptotic series was quite limited. A general mathematical approach for the Laplace equation in corners with one or two spring (Robin) boundary conditions, based on the solution of suitably defined recursive boundary value problems (BVPs), was introduced in [29]. These developments motivated a more computationally oriented procedure presented in [30], for the singularity analysis of linear elastic corners under antiplane shear with a spring boundary condition on one corner face. Engineering Fracture Mechanics 288 (2023) 109293 3 S. Jiménez-Alfaro and V. Mantič Fig. 1. A semi-infinite crack in a straight spring (adhesive) interface between two half-planes. Notice that the spring distribution is acting in the out-of-plane direction. In the present paper, the antiplane solution for a semi-infinite crack located in a straight spring-interface between two linear elastic isotropic half-planes of the same material is deduced in the form of an asymptotic series, for the first time. This series represents the local solution in the neighbourhood of the tip of a finite crack in a linear elastic (spring type) interface. Due to symmetry arguments the solution is decomposed into a symmetric and a skew-symmetric part with respect to the interface line. Whereas the symmetric part of this solution is given by a series of smooth functions, the skew-symmetric part is given by the singular Mode III solution. An analytic expression for the displacements in this Mode III solution is given in the form of a double asymptotic series of the main and the so-called associated shadow terms. It is shown that the series of the shadow terms associated to a main term is infinite and all shadow terms include logarithmic terms. Therefore, the solution for stresses at the crack tip includes a logarithmic stress singularity, this singularity being comprehensively analysed. This Mode III solution, can be considered as the crack-tip solution in Mode III for a crack in a locally straight spring-interface between linear elastic adherents of the same isotropic material. A key advantage of the present novel procedure is its simplicity and the fact that it provides the complete asymptotic series for the crack-tip solution, which is an advantage in comparison to some previous related works where only the first terms were presented. Moreover, it can be easily implemented in a computer algebra software, like Mathematica [31], which has been shown to be a key feature as the higher-order shadow terms become progressively increasingly complicated. The automatic generation of the expressions for the terms of the asymptotic series can be very useful for their usage in the implementation of new enriched or singular crack-tip finite elements for cracks propagating in spring interfaces, and also for the generation of testing problems for FEM codes. In Section 2the antiplane shear problem for a semi-infinite crack at a straight spring-interface is formulated, and the symmetry arguments are used to reduce it to a corner-singularity problem for a half-plane. Then, in Section 3the analytic procedure based on the complex variable formulation is presented, a pseudocode implementing this procedure being shown in Section 4. The analytical expressions for the main terms and for the first two shadow terms are presented in Section 5. Graphical representations of the singular solution, given by a number of terms in the asymptotic series providing a sufficiently high accuracy of plots, are shown and discussed in Section 6. Finally, the Energy Release Rate (ERR) is computed for a crack in a spring interface by the Irwin crack closure integral in Section 7. 2. Description of the model As mentioned in Section 1, the aim of this paper is to study an antiplane elastic problem of a semi-infinite crack in a straight spring interface between two half-planes (two half-spaces in 3D view) of the same linear elastic isotropic materials. Therefore, linear elastic springs are acting in the out-of-plane 𝑧-direction. The solution of the problem is described in terms of the out-of-plane displacement 𝑢(𝑥, 𝑦) = 𝑢𝑧(𝑥, 𝑦), a harmonic function that can be approximated by an asymptotic series expansion in the vicinity of the crack-tip. The crack-tip is located at the origin 𝑂of the Cartesian coordinate system, see Fig. 1. The considered BVP (Pc), defined in the two-half planes, (𝛺+∪𝛺−)⊂R2, bonded by a semi-infinite spring interface, represented in Fig. 1, is formulated as follows 𝛥𝑢 = 0 in 𝛺+∪𝛺−,(1) (Pc)𝜎𝑦𝑧 =𝑘[𝑢(𝑥, 0+) − 𝑢(𝑥, 0−)]for 𝑥 > 0and 𝑦= 0,(2) Engineering Fracture Mechanics 288 (2023) 109293 4 S. Jiménez-Alfaro and V. Mantič 𝜎𝑦𝑧 = 0 for 𝑥 < 0and 𝑦= 0.(3) Along the positive part of the 𝑥-axis the adhesive interface between half-planes 𝛺+and 𝛺−is modelled by a continuous distribution of linear elastic springs, known in the applied-mathematics community as the Robin boundary condition, cf. [29,32], which is characterized by the stiffness 𝑘=𝑘𝑧>0, associated to the out-of-plane 𝑧-direction. The semi-infinite crack located along the negative part of the 𝑥-axis is described by the stress-free condition. The stress component 𝜎𝑦𝑧(𝑥, 𝑦)is related to 𝑢(𝑥, 𝑦)using the constitutive law 𝜎𝑦𝑧 =𝜇 𝑢,𝑦and 𝜎𝑥𝑧 =𝜇 𝑢,𝑥, with 𝜇 > 0being the shear elastic modulus of the bulk material. Therefore, (2) and (3) can be rewritten as 𝜕𝑢 𝜕𝑦 =𝑘 𝜇[𝑢(𝑥, 0+) − 𝑢(𝑥, 0−)]if 𝑥 > 0and 𝑦= 0,(4) 𝜕𝑢 𝜕𝑦 = 0 if 𝑥 < 0and 𝑦= 0.(5) Similarly as in [33], the solution 𝑢(𝑥, 𝑦)is decomposed as the sum of an even (symmetric) and an odd (skew-symmetric) part with respect to the coordinate 𝑦, as 𝑢e(𝑥, 𝑦) = 𝑢(𝑥, 𝑦) + 𝑢(𝑥, −𝑦) 2=𝑢e(𝑥, −𝑦),(6) 𝑢o(𝑥, 𝑦) = 𝑢(𝑥, 𝑦) − 𝑢(𝑥, −𝑦) 2= −𝑢o(𝑥, −𝑦).(7) The superscripts eand odenote the even and the odd part, respectively, the sum of them gives the original-problem solution, 𝑢(𝑥, 𝑦) = 𝑢e(𝑥, 𝑦) + 𝑢o(𝑥, 𝑦).(8) As we can see from (6) and (7), once we know 𝑢eor 𝑢oin a half plane, 𝛺+or 𝛺−, it is also given in the other half-plane. Applying the definitions for 𝑢e(𝑥, 𝑦)and 𝑢o(𝑥, 𝑦), and the equilibrium of tractions on the interface (i.e., the continuity of 𝜎𝑦𝑧 and of 𝑢,𝑦 across the interface), the BVP (Pc)previously presented is split into two problems defined, e.g., in the half-plane 𝛺+. The first one is for the even part of 𝑢(𝑥, 𝑦) 𝛥𝑢e= 0 in 𝛺+,(9) (Pe)𝜕𝑢e 𝜕𝑦 = 0 if 𝑥 > 0and 𝑦= 0,(10) 𝜕𝑢e 𝜕𝑦 = 0 if 𝑥 < 0and 𝑦= 0,(11) and the second for the odd part of 𝑢(𝑥, 𝑦) 𝛥𝑢o= 0 in 𝛺+,(12) (Po)𝜕𝑢o 𝜕𝑦 =𝛾𝑢o(𝑥, 0+)if 𝑥 > 0and 𝑦= 0,(13) 𝜕𝑢o 𝜕𝑦 = 0 if 𝑥 < 0and 𝑦= 0.(14) A governing parameter 𝛾of dimension 1∕length used in (13) is defined as 𝛾=2𝑘 𝜇.(15) The solutions to the BVP (Pe)in the Cartesian coordinates system are given by any finite linear combination or an infinite series of a subset of homogeneous harmonic polynomials, which include an even power of 𝑦, as it was obtained in [33]. These homogeneous harmonic polynomials can be written in the polar coordinates system as 𝑢𝑒 𝑗(𝑟, 𝜃) = 𝑟𝜆𝑗cos(𝜆𝑗𝜃),for 𝜆𝑗=𝑗− 1, 𝑗 ∈N.(16) For example, for 𝑗= 1,2and 3, these polynomial solutions to the BVP (Pe)written in Cartesian coordinates are: 𝑢𝑒 𝑗(𝑥, 𝑦)=1,𝑥and 𝑥2−𝑦2. Notice that in (16) the index 𝑗∈N1to avoid an infinite energy in the vicinity of the crack-tip, which corresponds to the case 𝑗≤0and 𝑗∈Z. The BVP named as (Po)is solved and studied in the following sections. This problem is also called the Robin–Neumann problem because of the boundary conditions imposed on the boundary of the half-plane 𝛺+, the 𝑥-axis. For the sake of notation simplicity, in the next sections the odd part of the complete problem solution, which corresponds to a spring interface with a crack under Mode III loading condition, will be also denoted as 𝑢=𝑢o. Noteworthy, the solution 𝑢𝑒and 𝑢𝑜are not unique, because no remote boundary conditions are prescribed at infinity. 1The following subsets of real numbers will be used in this work: Zthe set of integer numbers, Nthe set of natural numbers. Engineering Fracture Mechanics 288 (2023) 109293 5 S. Jiménez-Alfaro and V. Mantič Fig. 2. Scheme of the problem in the vicinity of the singularity. Notice that the spring distribution is acting in the out-of-plane direction. 3. Deduction of the asymptotic series expansion The solution to the Robin–Neumann problem (Po)in 𝛺+can be represented by a doubly infinite asymptotic series expansion, which can be considered as a particular case of the general asymptotic expansion deduced in [30] for corners, by setting the inner angle of corner (angular sector) as 𝜔=𝜋, i.e. defining a half-plane. The 𝑗th element of this series, denoted as 𝑢𝑗, is defined by the sum of a main term (main singularity) 𝑢(0) 𝑗and an infinite series of shadow terms (shadow singularities) 𝑢(𝑘) 𝑗. This doubly infinite asymptotic series expansion can be approximated by a truncated series of 𝑁terms 𝑢𝑗, each one with 𝑆𝑗shadow terms, 𝑢(𝑟, 𝜃) ≈ 𝑁 ∑ 𝑗=1 𝐾𝑗𝑢𝑗(𝑟, 𝜃) = 𝑁 ∑ 𝑗=1 𝐾𝑗⎛⎜⎜⎝ 𝑢(0) 𝑗(𝑟, 𝜃) + 𝑆𝑗 ∑ 𝑘=1 𝑢(𝑘) 𝑗(𝑟, 𝜃)⎞⎟⎟⎠ ,(17) For 𝑁and 𝑆𝑗going to infinity the complete doubly infinite asymptotic series is recovered. The coefficients 𝐾𝑗are called Generalized Stress Intensity Factors (GSIFs), and are defined by the overall configuration of the problem, i.e., the loads applied and the whole geometry of the specimen. The expression (17) is written in a polar coordinates system that is located at the origin 𝑂, the crack tip, as it can be seen in Fig. 2. For the sake of simplicity in the application of the boundary conditions, first we find the asymptotic series expansions of the solution for the so-called Neumann–Robin (N–R) problem, see Fig. 2(a). Then, by a change of variable 𝜃=𝜋−𝜃, (18) the solution to the Robin–Neumann (R–N) problem is obtained, see Fig. 2(b). This change is made to follow the tradition of representing a crack growing from left to right, so in Fig. 2(b),𝛤2corresponds to the upper face of the crack. The main terms 𝑢(0) 𝑗are calculated replacing the Robin boundary condition (13) by a homogeneous Neumann boundary condition, i.e., solving the Neumann–Neumann problem. Its analytical solution is given by an asymptotic series whose terms are in fact defined in (16). As follows from (16) and (17), the dimension of 𝐾𝑗is (length)1−𝜆𝑗= (length)2−𝑗. It is important to highlight that only 𝑗∈N are considered, since an infinite strain energy is obtained for 𝑗≤0, as explained above. The shadow terms 𝑢(𝑘) 𝑗in the N–R problem, with 𝑘= 1,…, 𝑆𝑗, are obtained by solving the following recursive BVPs 𝛥𝑢(𝑘) 𝑗= 0 in 𝛺+,(19) (Pk 𝑗) 𝜕𝑢(𝑘) 𝑗 𝜕𝜃 = 0 on 𝛤1,(20) 1 𝑟 𝜕𝑢(𝑘) 𝑗 𝜕𝜃 = −𝛾𝑢(𝑘−1) 𝑗on 𝛤2,(21) where 𝛤1and 𝛤2are the boundary parts defined by 𝑦= 0 and 𝑥 > 0and 𝑥 < 0, respectively, see Fig. 2(a). Following [30], terms 𝑢(𝑘) 𝑗(𝑟, 𝜃)are given by the complex variable expression written as 𝑢(𝑘) 𝑗(𝑧) = 𝑘 ∑ 𝑙=0 [𝑎(𝑙) 𝑗,𝑘 ⋅Im{𝑧𝜆𝑗+𝑘log𝑙𝑧}+𝑏(𝑙) 𝑗,𝑘 ⋅Re{𝑧𝜆𝑗+𝑘log𝑙𝑧}],(22) where the complex number associated to a point in the upper half-plane is 𝑧=𝑥+𝑖𝑦 =𝑟exp(𝑖𝜃), and 𝑎(𝑙) 𝑗,𝑘 and 𝑏(𝑙) 𝑗,𝑘 are real constant coefficients providing real valued displacements 𝑢(𝑘) 𝑗(𝑧). Since 𝑘is the index related to each shadow term in the series, we can observe that the maximum power of the logarithmic terms in the expression of 𝑢𝑗is increased each time we include an additional shadow term in the series. The expression in (22) can be rewritten using real variables and applying the binomial expansion as 𝑢(𝑘) 𝑗(𝑟, 𝜃) = 𝑟𝜆𝑗+𝑘 𝑘 ∑ 𝑙=0 𝑙 ∑ 𝑚=0 (𝑙 𝑚)log𝑚(𝑟)𝜃𝑙−𝑚[𝑎(𝑙) 𝑗,𝑘 sin 𝑡(𝑙−𝑚) 𝑗,𝑘 (𝜃) + 𝑏(𝑙) 𝑗,𝑘 cos 𝑡(𝑙−𝑚) 𝑗,𝑘 (𝜃)],(23) Engineering Fracture Mechanics 288 (2023) 109293 6 S. Jiménez-Alfaro and V. Mantič where (𝑙 𝑚)=𝑙! 𝑚!(𝑙−𝑚)! and 𝑡(𝑙−𝑚) 𝑗,𝑘 (𝜃)=(𝜆𝑗+𝑘)𝜃+𝜋(𝑙−𝑚) 2.(24) Regarding the homogeneous Neumann boundary condition in (20), we have 𝜕𝑢(𝑘) 𝑗 𝜕𝜃 ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐𝜃=0 =𝑟(𝜆𝑗+𝑘) 𝑘 ∑ 𝑙=0 𝑎(𝑙) 𝑗,𝑘 ((𝜆𝑗+𝑘) log𝑙𝑟+𝑙log𝑙−1 𝑟)= 0.(25) In order to fulfil (25) for all 𝑟 > 0, coefficients 𝑎(𝑙) 𝑗,𝑘 for all 𝑗, 𝑘 and 𝑙are zero. The Robin boundary condition in (21) is written using real variables and simplifying terms 𝑟𝜆𝑗+𝑘−1 in both sides as 𝑘 ∑ 𝑙=0 𝑏(𝑙) 𝑗,𝑘 𝑙 ∑ 𝑚=0 (𝑙 𝑚)log𝑚(𝑟)𝜋𝑙−𝑚[𝑙−𝑚 𝜋cos 𝑡(𝑙−𝑚) 𝑗,𝑘 (𝜋)−(𝜆𝑗+𝑘) sin 𝑡(𝑙−𝑚) 𝑗,𝑘 (𝜋)] = −𝛾 𝑘−1 ∑ 𝑙=0 𝑏(𝑙) 𝑗,𝑘−1 𝑙 ∑ 𝑚=0 (𝑙 𝑚)log𝑚(𝑟)𝜋𝑙−𝑚cos 𝑡(𝑙−𝑚) 𝑗,𝑘−1(𝜋).(26) Since (26) must be verified for all 𝑟 > 0, coefficients 𝑏(𝑙) 𝑗,𝑘 are determined by comparing terms with the same power of log𝑚𝑟on both sides of the equation, leading to a linear system of 𝑘+ 1 equations. For each 𝑘th shadow term, this system is expressed in matrix form as 𝐌𝑗,𝑘𝐁𝑗,𝑘 =𝐆𝑗,𝑘,(27) to determine the vector of coefficients 𝐁𝑗,𝑘 =[𝑏(0) 𝑗,𝑘,…, 𝑏(𝑚) 𝑗,𝑘 ,…, 𝑏(𝑘) 𝑗,𝑘]𝑇. The terms in the vector 𝐆𝑗,𝑘 =[𝐺(0) 𝑗,𝑘,…, 𝐺(𝑚),…, 𝐺(𝑘) 𝑗,𝑘]𝑇are written, after some simplifications, as 𝐺(𝑚) 𝑗,𝑘 =𝛾 𝑘−1 ∑ 𝑙=𝑚 𝑏(𝑙) 𝑗,𝑘−1(𝑙 𝑚)𝜋𝑙−𝑚cos [𝜋(𝑙−𝑚) 2],for 0≤𝑚 < 𝑘, (28) 𝐺(𝑘) 𝑗,𝑘 = 0.(29) This last element of the vector 𝐆𝑗,𝑘 is always zero because it is associated with the 𝑘th power of the logarithm, and in the (𝑘− 1)th shadow term (the shadow term used to define the vector 𝐆𝑗,𝑘) this power does not appear. The overall structure of the square matrix 𝐌𝑗,𝑘 is 𝐌𝑗,𝑘 =⎡⎢⎢⎢⎢⎢⎣ 0⋯𝑀(0,𝑙) 𝑗,𝑘 ⋯𝑀(0,𝑘) 𝑗,𝑘 0⋱ ⋮ ⋮ ⋱ ⋱ ⋮ ⋮0⋱ ⋮ 0⋯ ⋯ 0 0 ⎤⎥⎥⎥⎥⎥⎦ ,(30) where each term of the matrix is expressed, after some simplifications, as 𝑀(𝑚,𝑙) 𝑗,𝑘 =(𝑙 𝑚)𝜋𝑙−𝑚[𝑙−𝑚 𝜋cos [𝜋(𝑙−𝑚) 2]− (𝜆𝑗+𝑘) sin [𝜋(𝑙−𝑚) 2]] for 𝑚 < 𝑙, (31) 𝑀(𝑚,𝑙) 𝑗,𝑘 = 0 for 𝑚≥𝑙. (32) Since the first column of the matrix 𝐌𝑗,𝑘 is zero, the first coefficient 𝑏(0) 𝑗,𝑘 is undefined, and the terms 𝑢(𝑘) 𝑗(𝑟, 𝜃)are not unique. However, it can be verified that the singularity order related to the coefficient 𝑏(0) 𝑗,𝑘,𝜆𝑗+𝑘, is equal to the singularity order of other main term in the series, 𝑗′> 𝑗, which is 𝜆𝑗′=𝜆𝑗+𝑘. The difference between two particular solutions of Eqs. (19)–(21) is therefore the solution to the Neumann–Neumann problem, which means that it is included in the main term in the other 𝑗′th element of the series, as it was explained for a general case in [30]. Thus, for the sake of simplicity and without lost of generality we can assume 𝑏(0) 𝑗,𝑘 = 0,for 𝑘≥1.(33) Hence, the non-uniqueness of the solution would be solved if we considered an infinite series for the solution, which is, unfortunately, numerically impossible. On the other hand, coefficients 𝑏(𝑚) 𝑗,𝑘 for 𝑚≠0are expressed taking advantage of the fact that the matrix 𝐌𝑗,1has upper triangular form, and therefore they can be obtained by a procedure similar to the Gaussian back substitution, 𝑏(𝑘) 𝑗,𝑘 = 𝐺(𝑘−1) 𝑗,𝑘 𝑀(𝑘−1,𝑘) 𝑗,𝑘 ,(34) 𝑏(𝑙) 𝑗,𝑘 =1 𝑀(𝑙−1,𝑙) 𝑗,𝑘 [𝐺(𝑙−1) 𝑗,𝑘 − 𝑘 ∑ 𝑠=𝑙+1 𝑀(𝑙−1,𝑠) 𝑗,𝑘 𝑏(𝑠) 𝑗,𝑘],for 1≤𝑙≤𝑘− 1.(35) Engineering Fracture Mechanics 288 (2023) 109293 7 S. Jiménez-Alfaro and V. Mantič Notice that the initial vector of coefficients, for 𝑘= 1, is defined through the main term, 𝑢(0) 𝑗(𝑟, 𝜃). This is equivalent to set 𝑏(0) 𝑗,0= 1 in the iterative process we have described. Moreover, remember that in this section we have obtained the solution to the Neumann–Robin (N–R) problem in 𝛺+. Therefore, the change of variable given in (18) is applied to (23) to get the shadow terms in the so-called Robin–Neumann (R–N) problem. As it has been shown in this section, solving N–R has less computational complexity than solving R–N, because coefficients 𝑎(𝑙) 𝑗,𝑘 are reduced to zero by applying the Neumann boundary condition at 𝜃= 0. Finally, with a finite number of terms in the asymptotic series the Robin boundary condition is not exactly fulfilled. Nevertheless, the error decreases in the vicinity of the crack tip when more shadow terms are included in the series, i.e., 𝑆𝑗is increased. Therefore, a suitable choice of 𝑆𝑗depends on the desired accuracy of the solution, which can be defined through the relative error 𝑒𝑗(𝑟), 𝑒𝑗(𝑟) = 1 𝑟 𝜕𝑢𝑗 𝜕𝜃 (𝑟, 0) 𝛾𝑢𝑗(𝑟, 0) − 1 ? ≈ 0.(36) The minimum number of shadow terms 𝑆𝑗to be included in the truncated series depends on the range of 𝑟where an accuracy, defined by a maximum value of error, is requested. A study of the function 𝑒(𝑟)in this range is required to decide a suitable value of 𝑆𝑗. Notice that higher 𝑆𝑗lead to higher computational complexities. The order (in terms of the power of 𝑟) of the error of a truncated series increases with increasing 𝑆𝑗. This leads to more accurate approximations of the exact solution in a vicinity of the crack tip by the series with increasing 𝑆𝑗, but outside this region the higher-order terms diverge more rapidly than the lower order terms, so the truncated series with a higher 𝑆𝑗can provide worse approximations for the solution far away from the crack tip, see e.g. Fig. 9. Recall that a similar behaviour can be observed in approximations of smooth functions by the Taylor series. 4. Pseudocode computing asymptotic series expansion A key advantage of the above procedure to deduce the asymptotic series expansion for a semi-infinite crack in a straight spring interface is that its implementation in computer algebra software is straightforward and very efficient. Thus, this implementation allows to obtain a fully analytical solution in terms of truncated series, which depends on how many shadow terms are considered in the series associated to a 𝑗th main term. In the following lines we propose a pseudocode implementing this procedure for R–N problem. In this code, 𝐵[𝑘, 𝑙]represents the unknown coefficients. Define 𝛾, 𝐺, 𝑗, 𝑆𝑗 𝜆𝑗=𝑗− 1, 𝑢(0) 𝑗(𝑟, 𝜃) = 𝑟𝜆𝑗cos 𝜆𝑗𝜃, and 𝐵[0,0] = 1 𝑡[𝜃, 𝑘, 𝑙, 𝑚]=(𝜆𝑗+𝑘)𝜃+𝜋(𝑙−𝑚) 2 For 𝑘= 1 To 𝑆𝑗Do 𝑀[𝑚, 𝑙]=(𝑙 𝑚)𝜋(𝑙−𝑚)[𝑙−𝑚 𝜋cos [𝜋(𝑙−𝑚) 2]− (𝜆𝑗+𝑘) sin [𝜋(𝑙−𝑚) 2]] 𝐺[𝑚]=𝛾 𝑘−1 ∑ 𝑙=𝑚 𝐵[𝑘− 1, 𝑙](𝑙 𝑚)𝜋𝑙−𝑚cos [𝜋(𝑙−𝑚) 2] 𝐵[𝑘, 𝑘]=𝐺[𝑘− 1] 𝑀[𝑘− 1, 𝑘]and 𝐵[𝑘, 0]= 0 𝐅𝐨𝐫 𝑙=𝑘− 1 𝐓𝐨 1𝐃𝐨 𝐵[𝑘, 𝑙]=1 𝑀[𝑙− 1, 𝑙][𝐺[𝑙− 1] − 𝑘 ∑ 𝑠=𝑙+1 𝑀[𝑙− 1, 𝑠]𝐵[𝑘, 𝑠]] 𝐄𝐧𝐝𝐅𝐨𝐫 EndFor 𝑢𝑗(𝑟, 𝜃) = 𝑢(0) 𝑗(𝑟, 𝜃) + 𝑆𝑗 ∑ 𝑘=1 𝑟𝜆𝑗+𝑘 𝑘 ∑ 𝑙=0 𝐵[𝑘, 𝑙] 𝑙 ∑ 𝑚=0 (𝑙 𝑚)log𝑚(𝑟)𝜃𝑙−𝑚cos (𝑡[𝜃, 𝑘, 𝑙, 𝑚]) 𝑢𝑗(𝑟, 𝜃) = 𝑢𝑗(𝑟, 𝜋 −𝜃) 𝑒𝑗(𝑟) = 1 𝑟 𝜕𝑢𝑗 𝜕𝜃 (𝑟, 0) 𝛾𝑢𝑗(𝑟, 0) − 1 Engineering Fracture Mechanics 288 (2023) 109293 8 S. Jiménez-Alfaro and V. Mantič 𝜎𝑟𝑧𝑗 (𝑟, 𝜃) = 𝜇 𝜕𝑢𝑗(𝑟, 𝜃) 𝜕𝑟 𝜎𝜃𝑧𝑗 (𝑟, 𝜃) = 𝜇 𝑟 𝜕𝑢𝑗(𝑟, 𝜃) 𝜕𝜃 5. Explicit expressions of the first two terms of the asymptotic series solution In this section, we deduce explicit expressions of the first terms of the asymptotic series solution for a semi-infinite crack situated in a straight spring interface, following the theoretical procedure explained above. We first solve the Neumann–Robin and then the Robin–Neumann problem shown in Fig. 2, deducing the expressions of the first two shadow terms, i.e., 𝑆𝑗= 2, for the sake of simplicity, giving 𝑢𝑗(𝑟, 𝜃) = 𝑢(0) 𝑗(𝑟, 𝜃) + 𝑢(1) 𝑗(𝑟, 𝜃) + 𝑢(2) 𝑗(𝑟, 𝜃),(37) where the main term representing a solution to the homogeneous Neumann–Neumann problem for the upper half-plane is expressed as 𝑢(0) 𝑗(𝑟, 𝜃) = 𝑟𝑗−1 cos [(𝑗− 1)𝜃],(38) see (16) and the associated discussion. 5.1. Explicit expressions for the Neumann–Robin problem For a better understanding of the matrix procedure explained above for the Neumann–Robin problem, see (27), the coefficients 𝑏(𝑙) 𝑗,𝑘 in (22), are explicitly computed for 𝑘= 1 and 2. As it was previously mentioned, the recursive explicit expressions for 𝑏(𝑙) 𝑗,𝑘 can be obtained by a procedure similar to the Gaussian back substitution (34)–(35), which also has lower computational complexity. For the first shadow term, 𝑘= 1, the following system of equations 𝐌𝑗,1𝐁𝑗,1=𝐆𝑗,1is solved, applying (28)–(29) to compute vector 𝐆𝑗,1and (31)–(32) to compute matrix 𝐌𝑗,1. The resulting system has the form (0 −𝑗𝜋 0 0 )⎛⎜⎜⎝ 𝑏(0) 𝑗,1 𝑏(1) 𝑗,1⎞⎟⎟⎠ =(𝛾 0).(39) This gives 𝑏(1) 𝑗,1=−𝛾 𝑗𝜋 and 𝑏(0) 𝑗,1= 0, see (33). By substituting these coefficients in (23), the final expression of the first shadow term, 𝑘= 1, for the N–R problem, is obtained, 𝑢(1) 𝑗(𝑟, 𝜃) = 𝛾𝑟𝑗 𝑗𝜋 [𝜃sin (𝑗𝜃) − log 𝑟cos (𝑗𝜃)].(40) Similarly, the system of equations for the second shadow term, 𝑘= 2, is obtained, ⎛⎜⎜⎜⎝ 0 −𝜋(1 + 𝑗) −2𝜋 0 0 −2𝜋(1 + 𝑗) 0 0 0 ⎞⎟⎟⎟⎠⎛⎜⎜⎜⎜⎝ 𝑏(0) 𝑗,2 𝑏(1) 𝑗,2 𝑏(2) 𝑗,2 ⎞⎟⎟⎟⎟⎠ =⎛⎜⎜⎜⎜⎝ 0 −𝛾2 𝑗𝜋 0 ⎞⎟⎟⎟⎟⎠ .(41) The solution of this system gives the expressions for the coefficients 𝑏(0) 𝑗,2= 0,𝑏(1) 𝑗,2=−𝛾2 𝑗𝜋2(1+𝑗)2and 𝑏(2) 𝑗,2=𝛾2 2𝜋2𝑗(1+𝑗), see (33). The final expression for the second shadow term, 𝑘= 2, is 𝑢(2) 𝑗(𝑟, 𝜃) = −𝛾2𝑟1+𝑗 𝑗(𝑗+ 1)𝜋2[𝑈1(𝑟, 𝜃) + 𝑈2(𝑟, 𝜃)],(42) where 𝑈1(𝑟, 𝜃) = −1 (𝑗+ 1) (𝜃sin [(1 + 𝑗)𝜃]− cos [(1 + 𝑗)𝜃]log 𝑟)(43) 𝑈2(𝑟, 𝜃) = 1 2[cos [(1 + 𝑗)𝜃]𝜃2+ log 𝑟2𝜃sin [(1 + 𝑗)𝜃]− log2𝑟cos [(1 + 𝑗)𝜃]].(44) By combining the previous expressions (38),(40) and (42), the asymptotic series solution of the N–R problem in (37) for 𝑗= 1 with two shadow terms, 𝑆1= 2, is obtained 𝑢1(𝑟, 𝜃) = 1 + 𝛾𝑟 𝜋[𝜃sin (𝜃) − cos (𝜃) log 𝑟]+𝑟2𝛾2 4𝜋2{−𝜃2cos 2𝜃+𝜃sin 2𝜃−(cos 2𝜃+ 2𝜃sin 2𝜃)log 𝑟+ cos 2𝜃log2𝑟}.(45) Noteworthy, the main term 𝑢(0) 1is given by a unit rigid body translation in the 𝑧−direction. The stress components associated with this solution are 𝜎𝜃𝑧1(𝑟, 𝜃) = 𝛾𝜇 𝜋[𝜃cos 𝜃+ sin 𝜃+ sin 𝜃log 𝑟]+𝑟𝜇𝛾2 4𝜋2[sin 2𝜃(1 + 2𝜃2)−4𝜃cos 2𝜃log 𝑟− 2 sin 2𝜃log2𝑟],(46) Engineering Fracture Mechanics 288 (2023) 109293 9 S. Jiménez-Alfaro and V. Mantič 𝜎𝑟𝑧1(𝑟, 𝜃) = 𝛾𝜇 𝜋[𝜃sin 𝜃− cos 𝜃− cos 𝜃log 𝑟]−𝑟𝜇𝛾2 4𝜋2[cos 2𝜃(1 + 2𝜃2)+4𝜃sin 2𝜃log 𝑟− 2 cos 2𝜃log2𝑟],(47) As follows from expressions (46) and (47), both stress components have a logarithmic singularity at the singular point 𝑟= 0. Nevertheless, as it could be expected, 𝜎𝜃𝑧1is bounded along the negative part of the 𝑥-axis (𝜃=𝜋) in the neighbourhood of this singular point, due to the Robin boundary condition there, whereas it vanishes along the positive part of the 𝑥-axis (𝜃= 0). It should be mentioned that Sinclair [22] presented the singular solution for the displacement 𝑢1(𝑟, 𝜃)up to the first shadow term and for the stress 𝜎𝜃𝑧1(𝑟, 𝜃)up to the second shadow term, matching the corresponding parts of the expressions (45) and (46), respectively. 5.2. Explicit expressions for the Robin–Neumann problem As explained in Section 3, the asymptotic series solution of the R–N problem in (37) with two shadow terms, 𝑆1= 2, can be easily obtained by applying the change of variable (18) to the asymptotic series solution of the N–R problem in (45), 𝑢1(𝑟, 𝜃) = 1 + 𝛾𝑟 𝜋[(𝜋−𝜃) sin 𝜃+ cos 𝜃log 𝑟] +𝛾2𝑟2 4𝜋2[−(𝜋−𝜃) sin 2𝜃− (𝜋−𝜃)2cos 2𝜃−(−2(𝜋−𝜃) sin 2𝜃+ cos 2𝜃)log 𝑟+ cos 2𝜃log2𝑟].(48) The stresses corresponding to these displacements are expressed as 𝜎𝜃𝑧1(𝑟, 𝜃) = − 𝛾𝜇 𝜋[sin 𝜃− (𝜋−𝜃) cos 𝜃+ sin 𝜃log 𝑟] +𝑟𝜇𝛾2 4𝜋2[(2(𝜋−𝜃)2+ 1) sin 2𝜃+ 4(𝜋−𝜃) cos 2𝜃log 𝑟− 2 sin 2𝜃log2𝑟],(49) and 𝜎𝑟𝑧1(𝑟, 𝜃) = 𝛾𝜇 𝜋[(𝜋−𝜃) sin 𝜃+ cos 𝜃+ cos 𝜃log 𝑟] +𝑟𝜇𝛾2 4𝜋2[−(1 + 2(𝜋−𝜃)2)cos 2𝜃+ 4(𝜋−𝜃) sin 2𝜃log 𝑟+ 2 cos 2𝜃log2𝑟].(50) Similarly as in the N–R problem, these stresses have a logarithmic singularity at 𝑟= 0, and 𝜎𝜃𝑧1is bounded along the positive part of the 𝑥-axis (𝜃= 0) near the singular point. Remark. According to Mishuris [23], displacements and stresses in the upper half domain of this R–N problem are expressed as 𝑢𝑀≈1 𝜋𝜇 {𝜇∗𝜏∗𝜋 1 + 𝜇∗ 𝑎− 0+ [𝐺𝑀(−1) − 𝑐− 0−𝑎− 0(1 − log 𝑟)]𝑟cos 𝜃+𝑎− 0(𝜋−𝜃)𝑟sin 𝜃},(51) 𝜎𝜃𝑧𝑀 ≈ − 1 𝜋[𝑎− 0log 𝑟sin 𝜃+ [𝐺𝑀(−1) − 𝑐− 0] sin 𝜃−𝑎− 0(𝜋−𝜃) cos 𝜃],(52) 𝜎𝑟𝑧𝑀 ≈1 𝜋[𝑎− 0log 𝑟cos 𝜃+ [𝐺𝑀(−1) − 𝑐− 0] cos 𝜃+𝑎− 0(𝜋−𝜃) sin 𝜃],(53) where the original author’s notation, with a subscript ‘‘𝑀’’ indicating a solution or parameter by Mishuris [23], is employed. In the present case, the materials, in the upper and in the lower half-planes, are considered identical. For this reason, 𝜇∗= 1. Moreover, the function 𝐺𝑀(𝑠) = 𝜇𝜕𝑢 𝜕𝜃 (𝑠, 𝜋)evaluated in 𝑠= −1 is 𝐺𝑀(−1) = 0. The coefficient 𝜏∗=2 𝛾. Therefore, we get 𝑢𝑀≈ 𝑎− 0 𝛾𝜇 −(𝑐− 0+𝑎− 0) 𝜋𝜇 𝑟cos 𝜃+ 𝑟𝑎− 0 𝜋𝜇 cos 𝜃log 𝑟+ 𝑟𝑎− 0 𝜋𝜇 (𝜋−𝜃) sin 𝜃, (54) 𝜎𝜃𝑧𝑀 ≈−1 𝜋[𝑎− 0sin 𝜃log 𝑟−𝑐− 0sin 𝜃−𝑎− 0(𝜋−𝜃) cos 𝜃],(55) 𝜎𝑟𝑧𝑀 ≈1 𝜋[𝑎− 0cos 𝜃log 𝑟−𝑐− 0cos 𝜃+𝑎− 0(𝜋−𝜃) sin 𝜃],(56) which can be compared to our solution for 𝑗= 1 and 𝑆1= 1, 𝑢1(𝑟, 𝜃) = 1 + 𝑟𝛾 𝜋cos 𝜃log 𝑟+𝑟𝛾 𝜋(𝜋−𝜃) sin 𝜃, (57) 𝜎𝜃𝑧1(𝑟, 𝜃) = − 𝛾𝜇 𝜋[sin 𝜃log 𝑟+ sin 𝜃− (𝜋−𝜃) cos 𝜃](58) 𝜎𝑟𝑧1(𝑟, 𝜃) = 𝛾𝜇 𝜋[cos 𝜃log 𝑟+ cos 𝜃+ (𝜋−𝜃) sin 𝜃].(59) If we normalize Mishuri’s solution in (54) by considering 𝑎− 0=𝛾𝜇, we get the same solution as in (57) except of a term given by a constant times 𝑟cos 𝜃, which represents the following higher order main term for 𝑗= 2 in (38). It is easy to see that this difference is associated with the choice of 𝑏(0) 𝑗,𝑘 = 0 in (33). Actually, taking 𝑎− 0=𝛾𝜇 and 𝑐− 0= −𝛾𝜇, the solution obtained in [23] in (54)–(56) is exactly the same as the one deduced in the present paper given by expressions (57)–(58), for 𝑗= 1 and 𝑆1= 1. Engineering Fracture Mechanics 288 (2023) 109293 16 S. Jiménez-Alfaro and V. Mantič References [1] Winkler E. Die Lehre von der Elasticität und Festigkeit mit besondere Rücksicht auf ihre Anwendung in der Technik, für polytechnische Schulen, Bauakademien, Ingenieure, Maschienenbauer, Architecten, etc.. Prag: Verlag von H. Dominicus; 1867. [2] Volkersen O. Die Nietkraftverteilung in Zugbeanspruchten Nietverbindungen mit Konstanten Laschen-querschnitten. Luftfahrtforschung 1938;15:41–7. [3] Goland M, Reissner E. The stresses in cemented joints. J Appl Mech 1944;11:A17–27. [4] Hart-Smith LJ. Analysis and design of advanced composite bonded joints. Technical report, NASA CR-2218, 1974. [5] Dillard DA, Mukherjee B, Karnal P, Batra RC, Frechette J. A review of Winkler’s foundation and its profound influence on adhesion and soft matter applications. Soft Matter 2018;14:3669–83. [6] Benveniste Y, Miloh T. Imperfect soft and stiff interfaces in two-dimensional elasticity. Mech Mater 2001;33:309–23. [7] Javili A, Kaessmair S, Steinmann P. General imperfect interfaces. Comput Methods Appl Mech Engrg 2014;275:76–97. [8] Tserpes K, Barroso-Caro A, Carraro PA, Carrillo Beber V, Floros I, Gamon W, Kozłowski M, Santandrea F, Shahverdi M, Skejić D, Bedon C, Rajčić V. A review on failure theories and simulation models for adhesive joints. J Adhes 2022;98(12):1855–915. [9] Prandtl L. Ein Gedankenmodell für den Zerreißvorgang spröder Körper (A thought model for the fracture of brittle solids). Z Phys Angenw Math Mech 1933;13(2):129–33, Translation: Knauss W. G., Int J Fract (2011) 171:105–109. [10] Destuynder P, Michavila F, Santos A, Ousset Y. Some theoretical aspects in computational analysis of adhesive lap joints. Internat J Numer Methods Engrg 1992;35:1237–62. [11] Erdogan F. Fracture mechanics of interfaces. In: Rossmanith HP, editor. Damage and failure of interfaces. Rotterdam: Balkema Publishers; 1997, p. 3–36. [12] Geymonat G, Krasucki F, Lenci S. Mathematical analysis of a bonded joint with a soft thin adhesive. Math Mech Solids 1999;4:201–25. [13] Lenci S. Analysis of a crack at a weak interface. Int J Fract 2001;108:275–90. [14] Mantič V, Távara L, Blázquez A, Graciani E, París F. A linear elastic - brittle interface model: Application for the onset and propagation of a fibre-matrix interface crack under biaxial transverse loads. Int J Fract 2015;195:15–38. [15] Kondratiev VA. Boundary-value problems for elliptic equations in domains with conical or angular points. Trans Moscow Math Soc 1967;16:227–313. [16] Costabel M, Dauge M. Construction of corner singularities for Agmon-Douglis-Nirenberg elliptic systems. Math Nachr 1993;162(1):209–37. [17] Wieghardt K. Über das Spalten und Zerreissen elastischer Körper (On splitting and cracking of elastic bodies). Z Math Phys 1907;55(1–2):60–103, Translation: Rossmanith H.P., Fatigue Fract Engng Mater Struct (1995) 12:1371–1405. [18] Williams ML. On the stress distribution at the base of a stationary crack. J Appl Mech 1957;24(1):109–14. [19] Paris PC, Sih GC. Stress analysis of cracks, Vol. 381. ASTM STP26584S, Philadelphia: ASTM; 1965, p. 30–81. [20] Yosibash Z. Singularities in elliptic boundary value problems and elasticity and their connection with failure. New York: Springer; 2012. [21] Destuynder P, Michavila F, Ousset Y, Santos A. Utilisation du taux de restitution de l’energie dans l’analyse a deux échelles de l’endommagement d’un joint collé. C R Acad Sci Paris 1990;310 Serie I:161–5. [22] Sinclair GB. A note on the removal of further breakdowns in classical solutions of Laplace’s equation on sectorial regions. J Elasticity 1999;56:247–52. [23] Mishuris G. Interface crack and nonideal interface concept (Mode III). Int J Fract 2001;107:279–96. [24] Mishuris GS, Kuhn G. Asymptotic behaviour of the elastic solution near the tip of a crack situated at a nonideal interface. ZAMM Z Angew Math Mech 2001;81(12):811–26. [25] Antipov YA, Avila-Pozos O, Kolaczkowski ST, Movchan AB. Mathematical model of delamination cracks on imperfect interfaces. Int J Solids Struct 2001;38(36–37):6665–97. [26] Sinclair GB. On the influence of cohesive stress-separation laws on elastic stress singularities. J Elasticity 1996;44:203–21. [27] Sinclair GB. A note on the influence of cohesive stress-separation laws on elastic stress singularities in antiplane shear. J Elasticity 2009;94:87–93. [28] Sinclair GB. On the influence of adhesive stress-separation laws on elastic stress singularities. J Elasticity 2015;118:187–206. [29] Mghazli Z. Regularity of an elliptic problem with mixed Dirichlet-Robin boundary conditions in a polygonal domain. Calcolo 1992;29(3–4):241—267. [30] Jiménez-Alfaro S, Villalba V, Mantič V. Singular elastic solutions in corners with spring boundary conditions under anti-plane shear. Int J Fract 2020;223:197–220. [31] Wolfram S. Mathematica, A system for doing mathematics by computer. Redwood City: Addison-Wesley; 1991. [32] Medková D. The Laplace equation: Boundary value problems on bounded and unbounded Lipschitz domains. Cham: Springer; 2018. [33] Baranova S, Mogilevskaya SG, Mantič V, Jiménez-Alfaro S. Analysis of the antiplane problem with an embedded zero thickness layer described by the Gurtin-Murdoch model. J Elasticity 2020;140(2):171–95. [34] Entov VM, Salganik RL. On the Prandtl brittle fracture model. Mech Solids 1968;3(6):79–89, (Translated from Russian). [35] Fernlund G, Spelt JK. Analytical method for calculating adhesive joint fracture parameters. Eng Fract Mech 1991;40:119—132. [36] Krenk S. Energy release rate of symmetric adhesive joints. Eng Fract Mech 1992;43(4):549–59. [37] Bruno D, Greco F. Mixed mode delamination in plates: a refined approach. Int J Solids Struct 2001;38:9149–77. [38] Shahin K, Taheri F. The strain energy release rates in adhesively bonded balanced and unbalanced specimens and lap joints. Int J Solids Struct 2008;45:6284—6300. [39] Carpinteri A, Cornetti P, Pugno N. Edge debonding in FRP strengthened beams: Stress versus energy failure criteria. Eng Struct 2009;31:2436–47. [40] Irwin GR. Analysis of stresses and strains near the end of a crack traversing a plate. J Appl Mech 1957;24:361–4. [41] Jiménez-Alfaro S, Mantič V. FEM benchmark problems for cracks with spring boundary conditions under antiplane shear loadings. Aerotec Missili Spaz 2020;99(4):309–19. [42] Mantič V, Vázquez-Sánchez A, Romero-Laborda M, Muñoz-Reja M, Jiménez-Alfaro S, Távara L. New crack-tip element for the logarithmic stress-singularity of cracks on spring interfaces in mode III, Comput Methods Appl Mech Eng (to be submitted). [43] Alberty J, Carstensen C, Funken SA. Remarks around 50 lines of matlab: short finite element implementation. Numer Algorithms 1999;20(2):117–37. [44] Yosibash Z, Schapira Y. Edge stress intensity functions along elliptic and part-elliptic 3D cracks. Eng Fract Mech 2021;245:107477.