Asymptotic solutions in anisotropic elastic multi-material corners with frictional contact

PhD Thesis Asymptotic solutions in anisotropic elastic multi-material corners with frictional contact María de los Ángeles Herrera Garrido Sevilla, 2024 Universidad de Sevilla Grupo de Elasticidad y Resistencia de Materiales Doctoral Program in Mechanical Engineering and Industrial Management Doctoral Thesis: Asymptotic solutions in anisotropic elastic multi-material corners with frictional contact Author: Mar´ ıa de los ´ Angeles Herrera Garrido Advisor: Prof. Vladislav Mantiˇ c 2024 Acknowledgements I would like to begin by expressing my heartfelt gratitude to everyone who has been involved, directly or indirectly, in the progress of my PhD thesis over the past nearly four years. I am deeply grateful to Prof. Vladislav Mantiˇ c, who has served as my advisor and my mentor throughout this thesis journey. It has been an extraordinary privilege to benefit from the time he has dedicated to me, the knowledge he has imparted, and the supportive environment he has created around him. Thank you for the opportunity to learn and grow, both academically and personally. I wish to express my gratitude to the entire Group of Elasticity and Strength of Materials (GERM). I extend my thanks to each and every member for their collective support. Without the combined efforts of professors, technical experts, PhD students, and administrative staff, the research group wouldn’t have gained the resources, knowledge, and recognition crucial for the completion of this thesis. Following this collective appreciation, I must acknowledge the individual efforts that have significantly impacted my path. Prof. Alberto Barroso, thank you for all the support you provided at the beginning of this journey. Dr. Israel Garc´ ıa, your questions and comments have helped to advance this work. Isabel del Principe, your assistance, always with a warm smile, has been deeply appreciated. Lastly, I am particularly grateful to Prof. Luis T´ avara and Dr. Mar´ ıa del Mar Mu˜ noz-Reja, for your mentorship, guidance, and the countless shared moments. Thank you for providing me with the opportunity to return to the academic world. How could I forget, even though it dates back to 2011, the place where the seed of research was planted in me. It was in the Departmento de Medios Continuos y Teor´ ıa de Estructuras at the Escuela Polit´ ecnica Superior de la Universidad de Sevilla, and for that, I will always be grateful for the trust and closeness shown by Dr. Beatriz Hortig´ on. During this thesis, I have had the chance to collaborate with Professor Roman Vodiˇ cka, from the Technical University of Koˇ sice. I owe him my gratitude for his ideas and the numerical results obtained through his boundary element code for the completion of Chapter 5. I will also never be able to adequately express the immense gratitude I feel for the invaluable opportunity I had to collaborate with Prof. Dominique Leguillon during my three-month stay at Sorbonne University for the development of the work corresponding to Chapter 9 of this thesis. I really appreciate how friendly and approachable he was, despite his high status in academia. His passing, not long after sharing this time together, makes all the shared moments, teachings, and dedication even more priceless. I would like to express my thanks to the CertBond – COST Action CA18120 (European Cooperation in Science and Technology) that supported part of my stay in Paris. Furthemore, I would like to thank the European Regional Development Fund as this thesis was supported by the Ministerio Espa˜ nol de Ciencia e Innovaci´ on and the Consejer´ ıa de Transformaci´ on Econ´ omica, Industria, Conocimiento y Universidades, Junta de Andaluc´ ıa [Projects: 1 2 M.A. Herrera-Garrido PGC2018-099197-B-I00, P18-FR-1928, US-1266016, PID2021-1233250B-I00]. Thanks to this funding, I have been able to thoroughly enjoy participating in multiple conferences over the past four years, all made possible by the generous support of the PIs of these projects, Prof. Vladislav Mantiˇ c, Prof. Alberto Barroso and Prof. Luis T´ avara. Returning to the University of Seville, I want to thank my friends in the department: Simone S., Ana M., Miguel M., Angel V. and Pablo O. I am especially thankful to Angel A., Carlos S., Mar´ ıa Luisa V., Teresa A. and Juanan C., who have been my constant support network and my endless source of encouragement and motivation. This achievement would not have been possible without their invaluable friendship. Juanca C., I first want to thank you for all your work, which forms the basis of much of Chapter 8 of this thesis, but above all, for your joy and your willingness to always be there and participate in everything. Sara J., thank you for welcoming me to the Sorbonne and making me feel closer to Seville during those months. Thank you for your optimism and for showing me that believing in people is possible. Seraf´ ın S., thank you for always being there, for getting me to talk; grabbing a coffee or ‘something’ with you is worth more than any session with a psychologist. Thank you, my friends Juanlu and Mar´ ıa. Starting each day with our walks has helped to refresh my mind. I want to express my gratitude to the most important people in my life. To my family, my mother, Socorro, for never doubting me and always supporting me. My brothers, Juanlu, Javi and Ra´ ul, and to my sisters-in-law, who are also like sisters, Eli and Bel´ en, thank you for giving me the encouragement and self-confidence needed to embark on this journey. And specially to my father, wherever you are, I hope you continue to be proud of me. This work is dedicated to you. Finally, thank you, Jos´ e. Thank you for supporting me and encouraging me to turn my life around and embark on this journey that is now coming to an end. Thank you for being there, for cheering me on, and for helping me through difficult times. Abstract The main objective of this thesis is to propose a methodology for characterizing singular elastic solutions near singular points in anisotropic single or multi-material corners through an asymptotic series expansion of power-law terms. The first part of this thesis conducts a state-of-the-art study, identifying various aspects yet to be covered, leading to the definition of five specific objectives. In the second part, a general matrix formalism for stress singularities analysis under generalized plane strain is reviewed and verified, along with its implementation in MATLAB code. Results obtained with its implementation are verified and compared with existing literature, revealing certain discrepancies, particularly in cases involving frictional contact when the plane x3=0 is not a plane of elastic symmetry. Special attention is paid to frictional solutions, including complex solutions and unexpected solutions showing coupling between in-plane and anti-plane modes despite elastic symmetries in the corner. The latter are validated using a numerical boundary element method. The third part of the thesis focuses on developing specific semi-analytical methodologies for two relevant problems: interface crack with frictional contact and frictional contact between two anisotropic wedges. These methodologies are based on the Stroh formalism and Coulomb friction law, allowing verification of solutions obtained through the general matrix formalism that could not be verified through comparison with other authors. In the fourth part, the SingSol web application is introduced. This web application, based on the Python language, represents an enhanced version of the MATLAB implementation of the general matrix formalism. It includes an optimization module that enables the computation of solutions in a fully autonomous way. Lastly, with the objective of analyzing the importance of considering frictional contact in the analysis of real engineering applications, a numerical analysis of the three-point end notched flexure test, is conducted, with a focus on the interface crack with frictional contact. A numerical methodology is developed to study crack growth under bending load, examining how the friction coefficient and stiffness discrepancies between materials affect the critical load and crack length increment. In conclusion, this thesis provides a comprehensive framework for analyzing stress singularities in multi-material corners with frictional contact, developing and validating specific semi-analytical methodologies and practical engineering applications. 3 Contents Contents 5 List of Tables 11 List of Figures 17 I Introduction 1 1 Introduction 3 1.1 Motivation.................................... 3 1.2 StateofArt ................................... 8 1.2.1 Singular solution with homogeneous boundary condition . . . . . . . 8 1.2.2 Singular solution with frictional contact condition . . . . . . . . . . . 9 1.2.2.1 Two-wedge frictional contact . . . . . . . . . . . . . . . . 11 1.2.2.2 The interface crack . . . . . . . . . . . . . . . . . . . . . 13 1.3 Objectivesofthethesis............................. 14 1.4 Outlineofthethesis............................... 16 II Matrix Formalism 19 2 Matrix formalism for homogeneous BC and IC 21 2.1 Introduction................................... 21 2.2 Stroh Formalism applied to stress singularity analysis of linear elastic anisotropic materials..................................... 22 2.3 Singular elastic solution in a single-material wedge. Transfer matrix . . . . . 25 2.4 Matrix formalism for boundary and interface conditions . . . . . . . . . . . . 27 2.4.1 Reference frame attached to a wedge face . . . . . . . . . . . . . . . 27 2.4.2 Boundary condition matrices . . . . . . . . . . . . . . . . . . . . . . 27 2.4.3 Interface conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.5 Characteristic system for the analysis of singularities in a multi-material elasticcorner .................................... 31 2.5.1 Multi-material-wedge transfer matrix . . . . . . . . . . . . . . . . . 32 2.5.2 Characteristic system assembly for a multi-material corner . . . . . . 33 2.5.2.1 Opencorner ......................... 34 2.5.2.2 Closed corner (periodic corner) . . . . . . . . . . . . . . . 36 2.5.3 Solution of the characteristic system. Singular elastic solution . . . . 38 5 12 M.A. Herrera-Garrido 3.5 Comparison of results for 0 ≤λ≤1 and ωobtained by the present code and those found in the literature, for open single-material orthotropic corners with at least one frictional boundary condition, Section 3.5.1.1. (*) Fulfils the frictional dissipation of energy condition. (o) No frictional tangential stress on the contact face, including anti-plane mode. (±)One contact face under compression and the other contact face under tension. . . . . . . . . . . . . . . . . . . . . . . . . 73 3.6 Comparison of results for 0 ≤λ≤1 and ωobtained by the present code and those found in the literature, for open single-material orthotropic corners with at least one frictional boundary condition, Section 3.5.1.1. (*) Fulfils the frictional dissipation of energy condition. (o) No frictional tangential stress on the contact face, including anti-plane mode. (±)One contact face under compression and the other contact face under tension. . . . . . . . . . . . . . . . . . . . . . . . . 74 3.7 Example 3.2.1. Comparison of results for 0 ≤λ≤1 and ωobtained by the present code and those found in the literature for a wedge sliding with friction over a semi-plane. (*) Fulfils the frictional dissipation of energy condition. . . . . 76 3.8 Numerical values for 0 <λ<1 and ωshown in Fig. 3.5, for open isotropic bi-material corners with frictional interface, see Section 3.5.1.2. (*) Fulfils the frictional dissipation of energy condition. . . . . . . . . . . . . . . . . . . . . . 77 3.9 Example 3.2.4. Numerical values of the 0 <λ<1 and ωshown in Fig. 3.6 for open orthotropic bi-material corners with friction interface described in the figure, see Section 3.5.1.2. (*) Fulfils the frictional dissipation of energy condition. (o) No frictional tangential stress on the contact face, including anti-plane mode. 79 3.10 Example 3.2.5. Numerical values for the singularity exponent 0 <λ<1 and direction of the tangential stress vector angle ωfor two single-material wedges sliding with friction in both interfaces, see Section 3.5.1.2. (*) Fulfils the frictional dissipation of energy condition. Shaded cells contain results obtained with the present code that match the results obtained in [111]. . . . . . . . . . . . . . 81 3.11 Example 3.3.1. Numerical values for the singularity exponent 0 <λ<1 and direction of the tangential stress vector angle ωfor an interface crack in an isotropic corner, see Section 3.5.2.1. (*) Fulfils the frictional dissipation of energy condition. (o) No frictional tangential stress in the interface, including antiplanemode...................................... 83 3.12 Example 3.3.2. Numerical values for the singularity exponent 0 <λ<1 and direction of the tangential stress vector angle ωfor an interface crack between anisotropic materials, see Section 3.5.2.1. (*) Fulfils the frictional dissipation of energy condition. (o) No frictional tangential stress in the interface, including anti-planemode. .................................. 84 3.13 Example 3.3.3. Numerical values for the singularity exponent 0 <λ<1 and direction of the tangential stress vector angle ωfor a crack terminating at a perfectly bonded interface in an isotropic corner with a friction coefficient in the crack µ=1, see Section 3.5.2.1. (*) Fulfils the frictional dissipation of energy condition. (o) No frictional tangential stress in the interface, including anti-plane mode. ........................................ 86 CHAPTER List of Tables 13 3.14 Example 3.4. Numerical values for the singularity exponent 0 <λ<1 and direction of the tangential stress vector angle ωfor two isotropic single-material wedges sliding with friction in both interfaces, see Section 3.5.2.2. In [1], a friction coefficient on one face is denoted as f, and a specific ratio of the friction coefficients on both faces as w. (*) Fulfils the frictional dissipation of energy condition. (o) No frictional tangential stress in the interface, including anti-plane mode. (±)One contact interface under compression and the other contact interfaceundertension. ................................. 87 4.1 Singularity exponent λand tangential stress angle ωsolutions obtained with both code versions for an anisotropic bi-material with friction contact interface crack. Elastic constants for the upper material and lower material are of material G with φ2=60◦and φ3=0◦, and material E with φ2=30◦and φ3=−2◦in Table 3.3, respectively. Friction coefficient µ=1..................101 4.2 Real part of the stress and displacement values obtained with both code versions at the friction contact interface of a bi-material. Elastic constants for the upper material and lower material are of material G with φ2=60◦and φ3=0◦, and material E with φ2=30◦and φ3=−2◦in Table 3.3, respectively. Friction coefficient µ=1...................................101 4.3 Imaginary part of the stress and displacement values obtained with both code versions at the friction contact interface of a bi-material. Elastic constants for the upper material and lower material are of material G with φ2=60◦and φ3=0◦, and material E with φ2=30◦and φ3=−2◦in Table 3.3, respectively. Friction coefficient µ=1...................................102 5.1 Engineering constants for the linear elastic materials used in the studied counterexamples. Elastic moduli and shear moduli in GPa. . . . . . . . . . . . . . . . 110 5.2 Eigensolutions found with the general matrix formalism for stress singularity problems Chapters 2 and 3 in an isotropic single-material wedge of aluminium with friction contact at θ=0◦and stress-free boundary condition at θ=300◦. Friction coefficient µ=1. .............................114 5.3 Eigensolutions found by the general matrix formalism for stress singularity problems (Chapters 2 and 3) in an orthotropic single-material wedge Fig. 5.1, made of material C in Table 5.1, with friction contact at θ=0◦and the stress-free boundary condition at θ=300◦. Friction coefficient µ=1.............120 5.4 Eigensolutions found with the general matrix formalism for stress singularity problems (Chapters 2 and 3) in an isotropic bi-material corner (Material A in Table 5.1), with friction contact interface at θ=0◦and stress-free boundary condition at θ=−60◦and at θ=300◦. Friction coefficient µ=1.........124 5.5 Eigensolutions found with the general matrix formalism for stress singularity problems (Chapters 2 and 3) in a closed bi-material corner schematized in Fig. 5.25, Mat. 1 is an isotropic adhesive and Mat. 2 an unidirectional fibre reinforced plastic ply (AS4/8552), material B and E in Table 5.1, respectively. Friction coefficient µ=1. (∗) for unexpected solutions . . . . . . . . . . . . . . . . . . . . . . 129 14 M.A. Herrera-Garrido 5.6 Summary of unexpected eigensolutions found by the general matrix formalism for stress singularity problems with different geometries and boundary conditions presented in Chapter 3. Material A and D elastic properties in Table 5.1. . . 133 6.1 Singularity exponent values λdepending on the PES in the bi-material, considering 0<Re(λ)<1. .................................146 6.2 Uniform stress solution modes for λ=1. .....................147 6.3 Analysis of the existence of a non-trivial solution with µ>0 and λ=1 2+n, n∈Zbased on the possible symmetries of the materials that form the interface crack.........................................151 6.4 Engineering constants for the materials used for the studied examples. Material A and B are graphite-epoxy fibre-reinforced composite materials. Shear and elasticmoduliinGPa. ...............................159 6.5 Numerically computed values of λand ˜ ωfor an interface crack between two dissimilar orthotropic materials as a function of the frictional coefficient µ. . . . 160 6.6 Numerically computed values of λand ˜ ωfor an interface crack between two laminaes of material A with φ(1) 2=45◦and φ(2) 2varying from −60◦to 90◦as a function of the frictional coefficient µ........................160 6.7 Numerically computed values of singularity exponent λand sliding angle ˜ ωfor an interface crack between an orthotropic material and an apparently anisotropic material with a friction coefficient µ=1 as a function of fibre angle of the anisotropic material φ(1) 3. The real values of λare computed by solving the system (6.48) and (7.44), and the complex values of λby a generalisation of the matrix formalism presented in Chapter 4. . . . . . . . . . . . . . . . . . . . . . 167 7.1 Engineering constants for the materials used for the studied examples in total precision. Shear and elastic moduli in GPa. . . . . . . . . . . . . . . . . . . . . 180 7.2 Parametrization of λand ˜ ωvs θ2in the stress free-stress free case. Friction coefficient µ=1 and θ0=−180◦..........................181 7.3 Parametrization of λand ˜ ωvs φ2in the stress free-stress free case. . . . . . . . . 183 7.4 Parametrization of λand ˜ ωvs µin the clamped-stress free case. . . . . . . . . . 184 8.1 Comparison of some of the most comprehensive computational methodologies for analyzing stress singularities identified in existing literature. . . . . . . . . . 190 8.2 Comparison of results for Example 3.1.1 obtained through optimization using (3.60) (first three columns) and (3.61) (last three columns). The results in the last three columns are derived using those from the first column as the initial startingpoint. ....................................198 8.3 Influence of the grid refinement (mesh refinement) on the number of solutions found and on computation time. The target number of solutions to find is 11. . . . 200 8.4 Comparison of results obtained with different minimization methods after the first minimization in Example 3.1.1. . . . . . . . . . . . . . . . . . . . . . . . . 202 8.5 Comparison of roots found using the Nelder-Mead minimization method with their corresponding minimum singular value σmin and the time taken, against the solutions presented in Table 3.4. . . . . . . . . . . . . . . . . . . . . . . . . . . 203 CHAPTER List of Tables 15 8.6 Comparison of roots found using the Powell minimization method with their corresponding minimum singular value σmin and the time taken, against the solutions presented in Table 3.4. The ’-’ symbol denotes solutions not found. . . . . 203 8.7 Comparison of roots found using the BFGS minimization method with their corresponding minimum singular value σmin and the time taken, against the solutions presented in Table 3.4. The ’-’ symbol denotes solutions not found. . . . . . 203 8.8 Comparison of roots found using the SLSQP minimization method with their corresponding minimum singular value σmin and the time taken, against the solutions presented in Table 3.4.The ’-’ symbol denotes solutions not found. . . . . 204 9.1 Dundurs mismatch parameter for isotropic bi-materials with E1=1 GPa and ν1=ν2=0.3.....................................227 10.1 Classification of the differences found between the results obtained with the present computational code based on the general matrix formalism and the results found in the literature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 234 List of Figures 1.1 Multi-material corner notation (2D view) . . . . . . . . . . . . . . . . . . . . . 4 1.2 The laminate saturated with transverse cracks and delaminations induced at their tips serves as an example of multiple multi-material corners in microscale (from [117])......................................... 5 1.3 Multi-material corners in an adhesively bonded double-lap joint between an aluminium plate and two carbon-fibre-reinforced polymer (CFRP) laminates (from Singular elastic solutions in anisotropic multi-material corners. Applications to composite[96]).................................... 5 1.4 F-18 wing to fuselage attachment. Example of metal-composite stepped adhesive joint (from Adhesively Bonded Joints [73]). . . . . . . . . . . . . . . . . . . . . 7 2.1 Schema of a single-material wedge. . . . . . . . . . . . . . . . . . . . . . . . . 26 2.2 Orthonormal basis attached to each face . . . . . . . . . . . . . . . . . . . . . . 27 2.3 Reference frame for an inclined plane . . . . . . . . . . . . . . . . . . . . . . . 30 2.4 Graphic representation of the real part of the determinant of the characteristic matrix of the multi-material corner. Example 2.1.12 . . . . . . . . . . . . . . . . 44 2.5 Singular stress and displacement fields for orthotropic material ‘G’ under fracture mode I. Results obtained from analytic expressions are represented by continuous lines, results obtained by the present code are represented by marks. . . . 49 2.6 Variation of the singularity exponent λvs the angle γ, which defines the vector ˜ nin the boundary condition ‘only u˜nrestricted’, increases on a single-material wedge of the orthotropic material ‘D’ with φ=0 (see Table 2.3). . . . . . . . . . 50 2.7 Singularity exponent λfor the tri-material junction with perfectly bonded interfaces and E3/E1=100, overlayed with [23, Fig. 9a] with perfect match of all results. Additionally, the computational code reveals a solution corresponding to the anti-plane mode (red-circles), not addressed in [23]. . . . . . . . . . . . . . . 50 3.1 Orthonormal bases attached to each face of a multi-material corner for frictional sliding contact. a) 3D view of an open corner with orthonormal basis on its outer faces, b) vectors in the plane defined by kand n, c) orthogonal projection on the x1x2-plane, d) orthogonal projection on the x1x3-plane. .............. 55 3.2 Schema of a corner with frictionless and frictional boundary and interface conditions. ....................................... 62 3.3 Fibre position for apparently anisotropic material . . . . . . . . . . . . . . . . . 69 17 18 M.A. Herrera-Garrido 3.4 Contour plot of σmin values of the SVD of the matrix Kcorner in Example 3.1.1, for 0 ≤λ≤1 and −π≤ω≤π. Red arrows are solutions founds also solving eigenequation from [121], black arrows are solutions that match with the antiplane mode and green arrows are unexpected solutions. . . . . . . . . . . . . . . 71 3.5 0 <λ<1 and ωresults in Table 3.8 presented together with [29, Fig. 3]. λBwith the eigenvalues plotted by Churchman and Hills [29] for the case of the bonded interface. See Ex 3.2.2 and Ex 3.2.3 in Table 3.8. . . . . . . . . . . . . . . . . . 77 3.6 Example 3.2.4. δ=λ−1 values versus upper wedge angle ϑ2for the case of a bi-material open corner with frictional contact interface presented in Table 3.9, superimposed with those obtained by Chen et al. [28, Fig. 5]. Mat. 1 and Mat. 2 are material E in Table 3.3, see Section 3.5.1.2. . . . . . . . . . . . . . . . . . . 80 3.7 Example 3.2.5. Schema of broken graphite/epoxy laminate as an open bi-material corner with frictional contact interface. Mat. 1 and Mat. 2 are material D in Table 3.3, and friction coefficient µ=0.5. Results presented in Table 3.10, see Section3.5.1.2. ................................... 80 3.8 Configuration used for interface crack examples. . . . . . . . . . . . . . . . . . 82 3.9 Corner configuration for examples in Table 3.13, crack reaching a perfectly bondedinterface. .................................. 85 4.1 Re(det(Kcorner)) as a function of λfor an anisotropic bi-material with a frictionless interface crack. Elastic constants for the upper material and lower material are of material G with φ2=60◦and φ3=0◦, and material E with φ2=30◦ and φ3=−2 in Table 3.3, respectively. This graph corresponds to the representation in the version of the code presented in Chapters 2 and 3, zoom not included. 96 4.2 Contour plot of the minimum singular values σmin of the SVD of the matrix Kcorner as a function of Re(λ)and Im(λ)for an anisotropic bi-material with a frictionless interface crack. Elastic constants for the upper material and lower material are of material G with φ2=60◦and φ3=0◦, and material E with φ2=30◦and φ3=−2 in Table 3.3, respectively. This graph corresponds to the representation in the new version of the code. . . . . . . . . . . . . . . . . . 96 4.3 Angular variations of the singular displacement and stress fields at r=1 for a frictionless delamination crack in a laminate with lower material G in Table 3.3 with φ2=60◦and φ3=0◦, and material E in Table 3.3 with φ2=30◦and φ3= −2. Figs. 4.3a and 4.3b depict the real and imaginary parts of the displacement field, respectively. Similarly, Figs. 4.3c and 4.3d illustrate the real and imaginary components of the stress field, respectively. The singularity exponent λ=0.5+ i0.000851280..................................... 98 4.4 In Fig. 4.4a, a contour plot illustrates the minimum singular values σmin resulting from the Singular Value Decomposition (SVD) of the matrix Kcorner(λ,ω). These values vary with Re(λ)and ωfor an anisotropic bi-material with a frictional interface crack, using elastic constants of materials G and E (see Table 3.3). The friction coefficient is set to µ=0. Additionally, in Fig. 4.4b, a plot shows the minimum singular values σmin of Kcorner with λ=0.5+i0.000851280 as a function of ωfor the same problem setup. These graphs correspond to the representation in the previous version of the code. . . . . . . . . . . . . . . . . . 98 CHAPTER List of Figures 19 4.5 3D contour plot of the minimum singular values σmin of the SVD of the matrix Kexpanded corner as a function of Re(λ),ω′and ω′′ with Im(λ) = 0 in Fig. 4.5a and Im(λ) = 0.000851280 for in Figs. 4.5b to 4.5d. The solved corner is an anisotropic bi-material interface crack with friction contact. Elastic constants for the upper material and lower material are of material G with φ2=60◦and φ3=0◦, and material E with φ2=30◦and φ3=−2◦in Table 3.3, respectively. Friction coefficient µ=0. This graph corresponds to the representation in the newversionofthecode. .............................. 99 4.6 2D contour plot of the minimum singular values σmin of the SVD of the matrix Kexpanded corner (0.5+i0.000851280,ω′,ω′′)as a function of ω′and ω′′. The solved corner is an anisotropic bi-material with friction contact interface crack. Elastic constants for the upper material and lower material are of material G with φ2=60◦and φ3=0◦, and material E with φ2=30◦and φ3=−2◦in Table 3.3, respectively. Friction coefficient µ=0. ......................100 4.7 Re(σθθ ), for a λ∈C, as a function of the crack tip distance rconsidering either λor Re(λ) in (2.98), for an anisotropic bi-material with a frictional interface crack. Elastic constants for the upper material and lower material are of material G with φ2=60◦and φ3=0◦, and material E with φ2=30◦and φ3=−2 in Table 3.3, respectively. Friction coefficient µ=1. Singularity exponent λ= 0.505290314 +i0.002086948. ...........................102 4.8 Illustrative case, Re(σθθ )for a λ∈C, as a function of the crack tip distance rconsidering either λor Re(λ) in (2.98), for an anisotropic bi-material with a frictional interface crack. Elastic constants for the upper material and lower material are of material G with φ2=60◦and φ3=0◦, and material E with φ2=30◦ and φ3=−2 in Table 3.3, respectively. Friction coefficient µ=1. Singularity exponent assumed for illustrative purposes λ=0.505290314 +i0.3. . . . . . . . 103 5.1 Uni-material problem schema. . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 5.2 Singularity exponent λvs friction coefficient µin a 300◦isotropic wedge, made of material A in Table 5.1, with frictional contact at θ0=0◦and stress free boundary condition at θ1=300◦...........................112 5.3 Tangential stress direction, angle ω, vs frictional coefficient in a 300◦isotropic material wedge, made of material A in Table 5.1, with frictional contact at θ0= 0◦and stress free boundary condition at θ1=300◦. ................112 5.4 Singularity exponent λvs wedge solid angle ϑin an isotropic wedge, made of material A in Table 5.1, with frictional contact at θ0=0◦with a friction coefficient µ=1 and stress free boundary condition at θ1=ϑ. Values of λ9,λ10,λ13,λ14 for wedge solid angle 235◦<ϑ<243◦obtained with the developed code in Chapter4.......................................113 5.5 Tangential stress direction angle ωvs wedge solid angle ϑin an isotropic material with contact friction at θ0=0 with µ=1 and stress free boundary condition at θ1=ϑ◦. Values of ω′ iand ω′′ iwith i=9,10,13,14 for wedge solid angle 235◦<ϑ<243◦obtained with the developed code in Chapter 4. . . . . . . . . . 114 5.6 Displacement of the nodes at the perimeter of the isotropic wedge with r=1, material A in Table 5.1, with the stress free boundary condition at θ=300◦and frictional contact at θ=0◦with friction coefficient µ=1. ............116 20 M.A. Herrera-Garrido 5.7 Comparison between the analytical and the numerical stresses on the frictional contact boundary at θ=0◦of an isotropic wedge, material A in Table 5.1, with stress free boundary condition at θ=300◦. Friction coefficient µ=1. . . . . . . 117 5.8 Singularity exponent λvs friction coefficient µin a 300◦orthotropic wedge, made of material C in Table 5.1, with frictional contact at θ0=0◦and stress free boundary condition at θ1=300◦...........................118 5.9 Direction of the tangential stress angle ωvs frictional coefficient in a 300◦orthotropic material with contact friction at θ0=0 and stress free boundary condition at θ1=300◦...................................118 5.10 Singularity exponent λvs solid angle ϑin deg. in an orthotropic wedge, made of material C in Table 5.1, with frictional contact at θ0=0◦with µ=1 and the stress free boundary condition at θ1=300◦. Values of λ13,λ14,λ17 and λ18 for wedge solid angle 316◦<ϑ<343◦obtained with the developed code in Chapter 4.119 5.11 Direction of the tangential stress angle ωvs solid angle ϑin deg. in an orthotropic material with frictional contact at θ0=0◦with µ=1 and the stress free boundary condition at θ1=300◦.Values of ω′ iand ω′′ iwith i=13,14,17,18 for wedge solid angle 316◦<ϑ<343◦obtained with the developed code in Chapter4.......................................120 5.12 Displacements corresponding to λ9in Table 5.3 for nodes on the perimeter of an orthotropic single-material wedge (Material C in Table 5.1) with friction contact at θ0=0◦with µ=1 and the stress-free boundary at θ1=300◦Fig. 5.1. . . . . . 121 5.13 Comparison between the analytical and the numerical stress fields corresponding to λ9in Table 5.3 on the frictional contact boundary at θ=0◦of an orthotropic single-material wedge (Material C in Table 5.1) with friction contact at θ0=0◦ with µ=1 and stress-free boundary at θ1=300◦Fig. 5.1. . . . . . . . . . . . . 122 5.14 Bi-material problem schema. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5.15 Singularity exponent λvs frictional coefficient µin an isotropic bi-material corner with stress free boundary conditions on both outer faces and frictional contact in the interface. Wedge 1 (material A in Table 5.1) goes from ϑ0=−60◦to ϑ1=0◦and wedge 2 (material A in Table 5.1) goes from ϑ1=0◦to ϑ2=300◦. . 124 5.16 The direction of the tangential stress angle ωvs frictional coefficient in an isotropic bi-material corner with the stress free boundary conditions on both outer faces and frictional contact in the interface. Wedge 1 goes from ϑ0=−60◦to ϑ1=0◦ and wedge 2 goes from ϑ1=0◦to ϑ2=300◦....................125 5.17 Singularity exponent λvs second wedge solid angle θ2=ϑin an isotropic bimaterial corner with the stress free boundary conditions on both outer faces and frictional contact in the interface with friction coefficient µ=1. First wedge goes from θ0=−60◦to θ1=0◦and second wedge goes from θ1=0◦to θ2=ϑ. Material A in Table 5.1 used for both materials in the corner. . . . . . . . . . . . 126 5.18 Direction of the tangential stress angle ωvs second wedge solid angle θ2=ϑ in an isotropic bi-material corner with the stress free boundary conditions on both outer faces and frictional contact in the interface with friction coefficient µ=1. First wedge goes from θ0=−60◦to θ1=0◦and second wedge goes from θ1=0◦to θ2=ϑ. Material A in Table 5.1 used for both materials in the corner.........................................127 CHAPTER List of Figures 21 5.19 Displacements of the nodes along the perimeter of the isotropic bi-material corner with r=1.....................................128 5.20 Comparison between the analytical and the numerical stress fields at the frictional contact boundary with µ=1 at θ=0◦of an isotropic bi-material corner (material A in Table 5.1) with the stress free boundary condition at θ=−60◦ and θ=300◦.....................................129 5.21 Singularity exponent λvs friction coefficient µin an orthotropic bi-material corner (Material C in Table 5.1) with the stress free boundary conditions on both outer faces and frictional contact in the interface. The first wedge goes from ϑ0=−60◦to ϑ1=0◦and second wedge goes from ϑ1=0◦to ϑ2=300◦. . . . 130 5.22 Direction of the tangential stress angle ωvs friction coefficient µin an orthotropic bi-material corner (Material C in Table 5.1) with the stress free boundary conditions on both outer faces and frictional contact in the interface. The first wedge goes from ϑ0=−60◦to ϑ1=0◦and second wedge goes from ϑ1=0◦ to ϑ2=300◦. ....................................131 5.23 Singularity exponent λvs second wedge solid angle ϑin an orthotropic bimaterial corner with the stress free boundary conditions on both outer faces and frictional contact in the interface with friction coefficient µ=1. The first wedge goes from θ0=−60◦to θ1=0◦and the second wedge goes from θ1=0◦to θ2=ϑ. .......................................132 5.24 Direction of the tangential stress angle ωvs second wedge solid angle ϑin an orthotropic bi-material corner with the stress free boundary conditions on both outer faces and frictional contact in the interface with friction coefficient µ=1. The first wedge goes from θ0=−60◦to θ1=0◦and the second wedge goes from θ1=0◦to θ2=ϑ. ..............................132 5.25 Closed bi-material corner schema. Mat. 1 is an isotropic adhesive and Mat. 2 an unidirectional fibre reinforced plastic ply (AS4/8552), material B and E in Table5.1,respectively ...............................133 6.1 Scheme of a semi-infinite interface crack with frictional contact with interface conditions and coordinate system definition. . . . . . . . . . . . . . . . . . . . . 140 6.2 Direction of the tangential stress in contact zone and the relative displacement betweenwedges. ..................................142 6.3 Limit states on the cone of the 3D Coulomb friction model. . . . . . . . . . . . . 148 6.4 (a) Graphical representation of (6.48) and (7.44) for 0 ≤λ≤1 and −π≤˜ ω≤π, for an interface crack between two apparently monoclinic materials. Material 1 and 2 are A in Table 6.4 and φ(1) 2=45◦and φ(2) 2=−60◦. Friction coefficient µ=0.5. Numerical values are shown in Table 6.6. For checking, (b) represents the contour plot of the minimum singular value σmin of the matrix Kcorner(λ,˜ ω) obtained by SVD, for 0 ≤λ≤1 and −π≤˜ ω≤π. Zero values of σmin represent the solutions obtained with the general code Chapter 3. Red circles correspond to solutions of the system (6.48) and (7.44). Green circles are solutions with λ=0.5.158 6.5 Possible orientation of the fibre in the apparently monoclinic composite material reinforced with long fibres (LFRC). From left to right, fibre lying in the planes x1−x3and x1−x2,respectively. ..........................159 CHAPTER 1 Introduction 1.1 Motivation Linear elastic fracture mechanics (LEFM) predicts that stresses near the tip of a crack or a v-notch or at junctions between different materials become unbounded [77]. Stress singularities can occur due to discontinuities such as changes in material properties (e.g., dissimilar material joints, interface cracks), jumps in boundary or interface conditions, non-smooth geometry (e.g., V-shaped cracks or notches). These stresses, known as singular stresses or stress singularities, can usually be represented through an asymptotic series expansion, given by products of power-logarithmic terms in radial coordinate and angular functions: σi j(r,θ) = N ∑ n=1 L ∑ l=0 Knlrλn−1loglr f (nl) i j (θ)+..., (1.1) where Knare the coefficients of the series expansion called Generalized Stress Intensity Factors (GSIFs), ris the distance to the singular point, λnare the characteristic exponents or singular exponents (eigenvalues) that depend on the elastic properties and the local geometry of the problem, f(nl) i jk (θ)are characteristic or singular angular shape functions (eigenfunctions), N≥1 can be finite or infinite number and L≥0 is a finite number. For the sake of simplicity in this Thesis, L=0 is assumed, as L≥1 is associated with special cases involving repeated roots of the characteristic equation of the corner, whose algebraic multiplicity is greater than the geometric multiplicity [37, 120, 121, 134]. With this assumption, the term l is omitted and expression (1.1) reduces to σi j(r,θ) = N ∑ n=1 Knrλn−1f(n) i j (θ)+.... (1.2) The terms within these series, known as power-type singularities, provide solutions for the elliptic system consisting of three partial differential equations for generalized plane strain. They also meet the boundary and interface conditions in the neighbourhood of the corner tip. It is then common to say that singular stresses follow the power law σi j(r,θ)≈rλ−1. As can be seen when analyzing the previous expression, for distances to the singular point very small, r→0, and values of λin the range 0 <Re λ<1, the stress tends to infinity. This is why they are called unbounded stresses. 3 4 M.A. Herrera-Garrido M 1 M 2 M 3 M k M m M M ... ... θ 0 θ 1 θ 2 θ 3 θ k-1 θ k θ m θ M-1 θ M ... ... bound. cond. bound. cond. θ 0 ,..., θ M materials wedges frictionless contact interface x 1 x 2 θ r (W≤M) θ m-1 W0 ,..., ϑϑ 0 ϑ 1 ϑ 2 ϑ 1-w ϑ w ϑ W ϑ 1-W ϑ W ϑ 0 ϑ 1 w perfectly bonded interface Figure 1.1: Multi-material corner notation (2D view) Points (or edges in 3D view) where such discontinuities originate singular stresses are called singular points (or singular edges, e.g. crack fronts). These points are prone to failure initiation due to high stress values in their neighbourhood [88, 159], which is generally referred to as a corner. The singular point itself is called a corner tip. A general configuration where several materials meet at a singular point, see Fig. 1.1 is called a multi-material corner. For the sake of brevity, we will often refer to single-material wedge as a material and to multi-material wedge with perfectly bonded materials simply as a wedge. Multi-material corners exist across different scales, from the nanoscale to macroscale, in a wide range of heterogeneous materials including composites, multilayers, and polycrystals (e.g., metals, ceramics, rocks, and photovoltaic cells), as well as in structural elements. For instance, as depicted in Fig. 1.2, numerous multi-material corners can be observed at the microscale, exemplified by a 90° ply block saturated with transverse cracks and delaminations induced at their tips, resulting in the separation of the 0°/90° interface. Each of these crack tips can be modeled as a multi-material corner in a singularity analysis. At meso-scale, Fig. 1.3 depicts a structural element containing multiple multi-material corners. This figure is a schematic representation of a double-lap joint, specifically depicting the junction between an aluminum plate and a composite laminate with an epoxy adhesive layer. Notably, this figure highlights several multi-material corners, each of which is subjected to a detailed analysis as presented in [11]. Another illustrative example of a structural element at macro-scale with multiple multi-material corners can be seen in Fig. 1.4, where the wing-to-fuselage attachment exemplifies a metal-composite stepped adhesive joint. Each contact zone between the graphite composite and the titanium structure leads to the formation of a multi-material corner. As mentioned earlier, stress singularities are not only produced by material discontinuities as those observed in adhesive joints, but also at jumps of interface or boundary conditions. Among these conditions, friction has numerous implications in many areas of science CHAPTER 15 Figure 1.2: The laminate saturated with transverse cracks and delaminations induced at their tips serves as an example of multiple multi-material corners in microscale (from [117]). Figure 1.3: Multi-material corners in an adhesively bonded double-lap joint between an aluminium plate and two carbon-fibre-reinforced polymer (CFRP) laminates (from Singular elastic solutions in anisotropic multi-material corners. Applications to composite [96]). and technology, it is an inherent phenomenon in the operation of mechanical systems where structural elements are in contact with each other, and is responsible for significant energy loss and wear. Frictional contact is often associated with stress concentrations or stress singularities with unbounded stresses, e.g., in material forming with frictional contact between the work-piece and the tool, punch indentation, fretting fatigue cracks, interface cracks between dissimilar materials, and in joints [5]. Frictional contact occurs in many configurations of multi-material corners in deformable structures at different scales, as might be expected, see Table 1.1, and it is therefore necessary to consider frictional contact in the singularity analysis of such corners for a large number of engineering applications. In structures where a singular point appears, when analyzing the possibility of failure or the growth of a pre-existing crack, the approach followed in the literature is through a local stress analysis. This is because it has been shown [88, 121, 159] that the initiation and propagation of the crack can be predicted if factors such as the GSIF, the energy release rate (G), the crack opening displacement (COD), or the J integral are known. The GSIF quantifies the effect of stress singularity at the crack tip or a discontinuity. It measures how much stresses intensify at the crack tip compared to the nominal stresses applied to the material. There are various methods for calculating the GSIFs, as described by 6 M.A. Herrera-Garrido Schematic corner Application of frictional contact • Sharp wedge in contact with a semiplane [28, 30, 51]. • Tool sliding over a semi-plane workpiece • Shear-induced material detachment • Closed interface crack [32]. • Flat-punch indentation reaching a stick zone [83] • End-point of the contact zone in an interface crack [27, 136] • Interface crack [96, 97] • 8Crack kink [1] • The crack reaches the interface and produces debonding • Stiff cutting tool penetrating into compliant material • The crack reaches the interface and open [28] • Fig. 1.2 • The crack reaches a perfectly bonded interface [33] • The crack reaches the interface and produces debonding [26] • Typical failure type in laminated materials Table 1.1: Representative cases of practical applications of multi-material corners with frictional contact CHAPTER 17 Figure 1.4: F-18 wing to fuselage attachment. Example of metal-composite stepped adhesive joint (from Adhesively Bonded Joints [73]). Giner et al. [53] or by Barroso [8, Ch.3] in the dedicated chapter. Since GSIFs depend not only on local geometry and material characteristics but also on the overall stresses to which the structure is subjected, except in very simple cases, numerical modeling of the global structure is necessary for their calculation. In the literature review conducted by Barroso [8, Ch.3] on the calculation of GSIFs, a classification into local and global techniques is provided. Local techniques for calculating GSIFs involve the use of asymptotic solutions (see Equation 1.2) in the vicinity of the corner tip. As noted by Barroso, there are two main approaches: those that propose numerically solving the model and applying a fitting to the asymptotic solution in a postprocess [9], and those that propose generating special singular elements containing the asymptotic solution and using them to solve the numerical model. Regardless of the method used, if a local approach is employed for obtaining GSIFs Kk, it is first necessary to determine the singularity exponents λi, as well as the characteristic functions fi jk associated with each exponent. Knowing the GSIFs Kkand the characteristic singularity exponent λkis not only useful for predicting crack propagation through LEFM or even crack initiation through various failure criteria developed within this framework, such as the Coupled Criterion of Finite Fracture Mechanic (CCFFM) [21, 85], the Theory of Critical Distances (TCD) [36] or the Cohesive Zone Model (CZM) [6, 43], but it can also assist in improving the design of our structure. By considering various design options, such as possible modifications in geometry or variations in material properties, having prior knowledge of the singular characteristic exponent and GSIFs allows for comparison among the different design available options. While it may not be possible to entirely eliminate the singularity, comparing these values enables the selection of an option where the singularity is less severe [10, 13, 78]. Moreover, in numerical analyses, such as those performed using the finite element method (FEM), it is advantageous to conduct prior stress analysis in regions where stress singularities are expected. This enables the incorporation of specialized singular elements into the mesh around these areas, facilitating model convergence [76, 95, 103]. On the other hand, ignoring singularities in numerical calculations significantly diminishes the accuracy of the results. The importance of accounting for these singularities becomes more pronounced with the 8 M.A. Herrera-Garrido increasing interest in expanding non-local fracture criteria [36, 85]. The analytical solution can also be utilized to assess corresponding FEM models, especially to ascertain whether the mesh is adequately refined near the corner to capture the expected asymptotic behavior predicted analytically. 1.2 State of Art Due to the relevance of the analysis of stress singularities in general, clarified previously, and the related importance gained by including friction conditions in the model, numerous works focused on obtaining the possible singularity exponents λkthat would complete the power-law asymptotic series expansion in (1.2) to define singular stresses and displacements around the singular point or stress singularity have been found in the literature. Below is a review of the main ideas and methodologies found in the literature. First, we will review the general case of homogeneous boundary conditions or interfaces, and then focus on cases that include frictional contact, specifically delving into two of the most interesting applications, the frictional sliding contact between two wedges and the interface crack. 1.2.1 Singular solution with homogeneous boundary condition The most relevant early studies of stress singularities in isotropic and single-material corners, with free or fixed boundaries, by Wieghardt [143] and Williams [145], see Vasilopoulos [139] for a careful review, were generalised in a large number of later works. In particular, many analytical, semi-analytical and numerical approaches have been proposed to study stress singularities in multi-material corners. For an overview of the main works from the early studies to more recent developments, refer to Carpinteri and Paggi [24] and Mantiˇ c et al. [96, 97]. Numerical approaches, such as those developed in [88, 159], are more general, as they can solve some stress singularity problems that cannot be solved by analytical or semi-analytical approaches. However, analytical or semi-analytical approaches, such as those developed in [38] for isotropic materials and in [35, 96, 97] for anisotropic materials, which provide an explicit closed-form expression of the characteristic corner equation (eigenequation), usually lead to higher accuracy of results. For isotropic bi-material and tri-material corners, basic eigenequations were deduced, e.g., in [18, 39] and [118], respectively. A semi-analytic procedure for isotropic multimaterial corners was introduced by Dempsey and Sinclair [38]. Carpinteri and Paggi [22] extended their study to include non-homogeneous isotropic materials, examining the asymptotic stress field at the vertex of multi-material junctions composed of angularly non-homogeneous material wedges. One of the earliest studies of corner singularities in anisotropic materials were those by Bogy [16, 17], analysing singular stresses in anisotropic single and bi-material corners, and by Ting and Chou [137], paying attention to configurations with multiple (repeated) roots of the characteristic equation of an anisotropic material (so-called mathematically degenerate materials) and also with multiple roots of the corner eigenequation. In the last 30 years, many authors have focused their studies on anisotropic multi-material corners. Costabel et al.[35] developed a semi-analytic code for the computation of stress singularities for elastic multi-material corner problems with perfectly bonded interfaces and with single roots of the material characteristic equation (so-called mathematically non - degenerate materials). Wu [153] used the Stroh formalism to deduce eigenequations for CHAPTER 19 multi-material anisotropic corners with a wide range of boundary conditions, but not including frictional contact. Using the Stroh formalism [123, 124], several related approaches for the singularity analysis of anisotropic multi-material corners with perfectly bonded interfaces were developed in 2003. Hwu [66] focused on materials that are mathematically non-degenerate, while Yin [156] and Barroso et al. [11] addressed both non-degenerate and degenerate materials. In particular, the code implemented in [11] by using the computer algebra software Mathematica [148] was able to calculate characteristic exponents for any multi-material open or closed (periodic) corner in generalized plane strain, considering perfectly bonded interfaces, with free or clamped boundary faces. The great novelty of this code was that for the first time it was possible to analyse corners with mathematically non-degenerate, degenerate (e.g., isotropic) and extraordinarily degenerate materials. The procedure presented in [11] was further developed by Mantiˇ c et al. [96, 97], who introduced a general matrix formalism able to analyse multi-material corners problems with many different boundary conditions, including also frictional contact for boundary and interface conditions, and with any kind of linear elastic materials. This matrix formalism was established, and various ad-hoc tests were conducted in [96, 97] for its validation. However, there is still a need to implement this formalism into a general computational code and conduct subsequent verification. 1.2.2 Singular solution with frictional contact condition Although several contact models are available at present, the most commonly used, especially in computational structural analysis, is the Signorini-Coulomb contact model based on the Signorini unilateral contact law and the Coulomb friction law, as it proposes a relatively simple and widely accepted formulation for frictional contact problems, [5, 79, 112, 149, 155]. This is the contact model that will be considered in the following. Due to the inherently nonlinear nature of frictional contact, it is more difficult to include such contact in the analysis of singular elastic solutions in multi-material corner problems, than considering only standard linear boundary or interface conditions such as free, clamped, symmetry, etc. Incorporating a true 3D formulation of the Coulomb friction law described by the Coulomb friction cone is particularly difficult. Thus, although frictional contact has been considered in corner singularity analysis by many authors, they have usually made some simplifying assumptions. In particular, previous works have usually considered linear elastic isotropic materials in contact under plane strain (or plane stress) conditions [1, 5, 29, 30, 32, 33, 63, 121, 144, 154]. Few attempts have been made to analyse corner singularities with frictional contact between anisotropic materials, and even fewer to consider generalised plane strain (GPS) conditions, with the elastic fields independent of the x3-coordinate, i.e. ui=ui(x1,x2)with i=1,2,3. In the following, we will focus on the stress singularity analysis of linear elastic anisotropic multi-material corners in GPS with frictional sliding contact. Several authors simplified the 3D Coulomb friction law, represented by the Coulomb cone, by adopting some simplifying hypotheses. As a result, the collinearity of the frictional tangential stress and relative tangential displacement vectors may not be satisfied in some GPS cases, e.g., with a general orientation of the planes of elastic symmetry (PES) of contacting materials. For instance, some authors arbitrarily restricted the relative displacement vector (sliding 10 M.A. Herrera-Garrido motion vector) in the contact plane to a predetermined direction. In particular, Poonsawat et al. [111] and Magnier et al. [94] simplified the Coulomb frictional contact problem by assuming continuity of displacement in the x3-direction between the materials in contact. The differences with respect to the true 3D Coulomb frictional contact model that such a simplification implies in the results of the singularity analysis are discussed in Section 3.5. A step forward was made by Leguillon [84], who generalised this approach considering any given direction of the relative tangential displacement vector in the contact plane. In another series of articles, Ting and Chou [136] and Sung and Chung [125] simplified the Coulomb frictional contact law, with the friction coefficient µ, basically assuming that the values of the frictional tangential stresses σ12 and σ23 are either ±µσ22 or zero (considering the x2-direction normal to the contact x1x3-plane). This assumption may not be valid in configurations where there is no PES coinciding with the x3=0 coordinate plane, and where the relative displacements in the x1- and x3-directions cannot generally be decoupled, leading to incorrect results. This difficult problem of singularity analysis of linear elastic anisotropic multi-material corners with frictional sliding contact under GPS conditions was solved independently and almost simultaneously in 2012 by Chen et al. [28] and Mantiˇ c et al. [96]1, see also [97], who introduced two new but distinct semi-analytic approaches, both based on the Stroh formalism of anisotropic elasticity [123, 124]. The former approach [28] was formulated for a particular case of two anisotropic wedges in a frictional sliding contact, deducing a system of two nonlinear eigenequations for the stress singularity exponent λand the angle ωof the frictional tangential stress vector measured in the contact plane. The latter approach [96, 97] formulated for multi-material corners with any finite number of faces with isotropic and anisotropic frictional contact, developed an original and general matrix formalism leading to the search for the minimum singular value σmin of the rectangular matrix of a nonlinear eigenvalue system. σmin was considered as a function of the stress singularity exponent λand the angles of the frictional tangential stress vectors ωwmeasured in different contact planes. However, as mentioned in the previous section, a general and fully operative computational implementation of the latter approach was not developed at that time, only a few preliminary results were presented, so neither any verification by comparison with some previously published results was carried out. In [96, 97], the singularity analysis of a closed corner with the configuration shown in Fig. 1.3-B was conducted. Material 1 and material 2 in Fig. 1.3 correspond to material B (the adhesive) and G (the CFRP) in Table 2.3, respectively. The friction coefficient at the vertical interface, ϑ=90◦, is µ=1, and the horizontal interface, ϑ=0◦, is perfectly bonded. While the expected values of the characteristic tangential stress angles, ω=0◦,±90◦, and ±180◦, were identified, the analysis revealed additional characteristic values with ω=0◦,±90◦, and ±180◦. Notably, these unexpected solutions represent singular elastic solutions with coupled in-plane and out-of-plane displacements. It is worth mentioning that these solutions appear despite the symmetry of the analyzed corner configuration with respect to the x3=0 coordinate plane, thus these solutions were referred to as strange. It is essential to note that although these unexpected solutions were highlighted, no demonstration of their possible existence was provided. Moreover, it is crucial to acknowledge that, to the best of the author’s knowledge, such solutions have not been indicated by any other author in the existing literature. Thus, further exploration is necessary to identify corner singularity 1Book manuscript submitted to the publisher in 2012. CHAPTER 111 problems where similar strange singular solutions may arise and to validate the viability of these unconventional results. It is crucial to ascertain whether these strange solutions are genuine and not merely an outcome of an error in the formalism presented in [96, 97] or the employed ”ad-hoc” methodology. This involves seeking more cases where such solutions, where in-plane and anti-plane modes cannot be decoupled, may appear. If the existence of these solutions, in which the in-plane and anti-plane modes cannot be decoupled due to the presence of friction, is verified, it would demonstrate that in any analysis of singular stresses where friction is involved, separate analyses of the in-plane and anti-plane modes are not sufficient. Instead, a 3D analysis or at least a generalized plane strain analysis is necessary. This is because conducting either an in-plane or anti-plane mode analysis could result in the overlooking of these singular solutions, thus making the Wieghardt-Williams series (1.1), which characterizes the stress state in the multi-material corner, incomplete. As mentioned earlier, many authors have focused their efforts on solving specific cases presented in Table 1.1, rather than developing a generic tool capable of obtaining singular stress and displacement solutions for multi-material corners without limitations on the number of materials, their geometry, or boundary and interface conditions, including frictional contact. Among the cases discussed in the Table 1.1, this thesis focus the analysis on two of those cases that have produced significant interest in the scientific community due to their relevance: specifically, the case of two wedges sliding against each other with friction or without friction and the case of the interface crack with contact between crack faces. 1.2.2.1 Two-wedge frictional contact The first application mentioned in Table 1.1, the problem of a sharp wedge in sliding contact with a semi-plane problem has been widely analyze, mostly for isotropic materials. In this problem a sharp wedge slide with or without friction over a semi-plane. Most authors make use of Dundurs parallelograms [44] to effectively display results for all possible combinations of materials, utilizing the Dundurs’ mismatch parameters for isotropic bi-materials, denoted as αand β. Dundurs [70], analyzes the stress concentration at a sharp wedge with frictionless contact by applying the Mellin transformation. The main result presented is that if sharp wedges between contacting bodies are unavoidable, the analysis indicates that undercutting these edges can effectively reduce locally high stresses. Later, Gdoutos and Theocaris [51], employing a different approach, succeeded in deducing an expression for eigenequations for frictional contact between a half-plane and a sharp wedge . Their findings were corroborated when they particularized their expression for the frictionless case and obtained the same expression as presented in [70]. In this scenario, for a rigid indenter, the geometry of the indenter does not impact the singularity, while for an elastic indenter, the singularity becomes more severe as the indenter angle increases. Additionally, if the friction coefficient is positive and increases, so does the singularity, whereas if the friction coefficient is negative and increases, the singularity decreases. A negative friction coefficient, as commonly understood in the literature, implies sliding in the opposite direction as defined for a positive friction coefficient. Comninou [30] continued the work done by Dundurs [70], applying the Mellin transformation to the sharp edge frictional sliding contact, and obtained the same results as in [51], while being able to indicate when the logarithmic singularities appear. The expressions provided by Comninou [30] and Gdoutos and Theocaris [51] were particularized for the case of a square block sliding over an elastically similar half-plane by Karuppanan et al.[74]. The Part II Matrix Formalism 19 CHAPTER 2 General matrix formalism for anisotropic multi-material corners with homogeneous boundary and interface conditions 2.1 Introduction The present chapter, based on [59], applies and further develops the general matrix formalism for singularity analysis of multi-material corners introduced in [96, 97]. The description of this formalism in [96, 97] was very concise, which in view of its generality and complexity could make difficult its comprehension by readers. Furthermore this matrix formalism has not been verified so far by its full computational implementation, except for by ad hoc implementations for some specific cases studied in [96, 97]. Therefore, the aim of the present chapter is to present this matrix formalism in a comprehensive way for the case of homogeneous and orthogonal boundary and interface conditions (i.e. omitting friction contact, which makes its presentation much simpler), revise this matrix formalism and provide complementary and clarifying explanations. However, the most relevant contribution of this chapter is the first fully general implementation of the formalism introduced in [96, 97], and an exhausting checking of this implementation by comparing its results with the numerous analytical and numerical results of corner singularity studies available in the literature. In Chapter 3, both boundary and interface conditions with frictional sliding contact will be included in the present matrix formalism, and its computational implementation will be described together with some examples of its successful testing. Although some corner problems lead to eigenequations with multiple (repeated) roots, and such configurations were studied by Dempsey and Sinclair [38] for isotropic materials, and by Ting and Chou [137], Wu [153] and Steigemann [122] for anisotropic materials, the present chapter is limited to single roots of the corner eigenequations, for the sake of simplicity. 21 22 M.A. Herrera-Garrido This chapter presents a new code developed in MATLAB [100, 101], which is based on the previous one [11]. The possibility to consider interfaces between the materials in the corner with frictionless sliding, and several boundary conditions with allowed or restricted displacement directions have been added to this code, see Fig. 1.1 for the notation used to describe a multi-material corner. The possibility of sliding between materials implies that the size of the matrix to be solved is not always 3 ×3 or 6 ×6, but 6W×6W, where Wis the number of wedges formed by angular sectors of materials perfectly bonded together. The following ingredients are key in the development of the present matrix formalism: - Stroh formalism for anisotropic elasticity [123, 124] - Transfer matrix for a single-material wedge proposed by Ting [135] - Matrix formalism for homogeneous orthogonal boundary and interface conditions proposed by Mantiˇ c et al. [96, 97]. We speak about homogeneous orthogonal boundary conditions when the stress and displacement vectors are orthogonal to each other. This happens when there is a coordinate system in which always one of the components of one of these vectors is null. This means, when the displacement in one direction is restricted (equal to zero), the stress is allowed, and when the displacement is allowed in one direction, the stress will vanish in that direction. For example, in the symmetry boundary condition the displacement in the normal direction is zero un=0, but the normal stress component is in general non zero σn=0. The content of this chapter is divided into sections according to the steps that the implemented computer code follows: Material characterization by applying first the Stroh formalism in Section 2.2, and then computing the transfer matrix in Section 2.3. In the next step we obtain the boundary and interface matrices employing the matrix formalism for homogeneous orthogonal boundary and interface conditions presented in Section 2.4. In Section 2.5, all these matrices are assembled to form the characteristic matrix of the corner defining a nonlinear eigensystem. The numerical solution of this eigensystem gives the characteristic (singular) exponents. For some special cases such as single-material and bi-material corners, reduced eigensystems are deduced. The characteristic angular functions of stresses and displacements are calculated individually for each characteristic exponent in Section 2.6. The entire implementation is summarised in Section 2.7. Finally, in Section 2.8, many numerical examples are shown in order to check our own code and also the closed-form eigenequations presented by other authors. These examples include several relevant cases such as V-notches, interface cracks, cracks terminating at an interface, multi-material closed corners. Noteworthy, the present code can be used to test new analytic formulas for the characteristic exponents or eigenequations developed for relevant engineering problems, which are very suitable for practical applications and parametric studies. 2.2 Stroh Formalism applied to stress singularity analysis of linear elastic anisotropic materials Stroh formalism is a powerful tool to solve linear anisotropic elastic problems. This formalism is based on the equilibrium equations in terms of displacements (in the absence of body CHAPTER 223 forces) obtained by using the constitutive law of linear anisotropic elasticity, σi j =Ci jkl εkl =Ci jkluk,l,(2.1) Ci jkluk,l j =0,(2.2) where Ci jkl is the positive definite (Ci jklεi jεkl >0 for non zero εi j =εji)and symmetric fourthorder tensor of elastic stiffnesses (Ci jkl =Cjikl =Ckli j),σi j is the Cauchy stress tensor, uithe displacement field and εi j the small strain tensor. Assuming a generalized plane strain state, where the displacements uidepend exclusively on the coordinates x1and x2,ui(x1,x2), we consider a solution of the system of equations (2.2) in the following form uk(x1,x2) = akf(z),(2.3) z=x1+px2,(2.4) where f(z)is an analytic function of a variable zgiven by a linear combination of x1and x2, and pand aiare constants to be determined. Taking into account the geometry of the present corner problem (Fig. 1.1), it will be useful to employ zwritten in polar coordinates z=r(cos(θ)+ psin(θ)).(2.5) Differentiating twice the expression (2.3) with respect to xland xjand substituting in (2.2), we obtain the quadratic eigenvalue problem for the number pand the vector a [Ci1k1+p(Ci1k2+Ci2k1)+ p2Ci2k2]ak=0.(2.6) To write this eigenvalue problem in matrix form, we define the following matrices (secondorder tensors) Qik =Ci1k1,Rik =Ci1k2,Tik =Ci2k2.(2.7) Then, the quadratic eigenvalue problem (2.6) is written as [Q+p(R+RT)+ p2T]a=0.(2.8) A non trivial solution for amay exist only if the determinant of the matrix in this linear system is zero, detQ+p(R+RT)+ p2T=0.(2.9) This equation is known as the Lekhnitskii-Stroh sextic polynomial equation in pfor anisotropic materials under generalized plane strain [7]. Lekhnitskii [90] showed that the six dimensionless roots of this equation pα(α=1,...,6), called also eigenvalues, are complex, i.e. pα∈C. In the present work, we will assume the following ordering of these six eigenvalues according to the sign of their imaginary part Im pα>0,pα+3=¯pα,α=1,2,3,(2.10) where Im denotes the imaginary part and the overbar denotes the complex conjugate value. For an eigenvalues pα, the associated eigenvectors aαcan be obtained from the linear system (2.8). The stress tensor components σi1and σi2can be obtained by the following expressions σi1= (Qik +pRik)akf′(z),(2.11) 24 M.A. Herrera-Garrido σi2= (Rki +pTik)akf′(z),(2.12) while σ33 is obtained from the constitutive law by taking into account that ε33 =0. Expressions (2.11) and (2.12) can be rewritten as σi1=−pbif′(z),σi2=bif′(z),(2.13) where the eigenvector bis obtained from the eigenvector aas b= (RT+pT)a=−1 p(Q+pR)a.(2.14) The eigenvectors aand bare ordered according to the ordering of pαin (2.10) aα+3=¯ aαand bα+3=¯ bα,α=1,2,3.(2.15) By introducing the stress function vector φas ϕi=bif(z),or in matrix notation φ=bf(z),(2.16) the expressions (2.13) take the following form σi1=−ϕi,2,σi2=ϕi,1.(2.17) The stress function vector φrepresents, together with the displacement vector u, the main variables of the Stroh formalism. For mathematically non-degenerate materials, i.e. those with the six different eigenvalues pα, we can superimpose the six solutions for displacement vector from (2.3), and the corresponding six solutions for the stress function vector from (2.16), giving, in view of (2.15), u(x1,x2) = 3 ∑ α=1 aαfα(zα)+ ¯ aαfα(¯ zα),(2.18) φ(x1,x2) = 3 ∑ α=1 bαfα(zα)+ ¯ bαfα(¯ zα).(2.19) Analogous fundamental representations of the Stroh formalism of the displacement and stress function vectors for degenerate materials, with multiple eigenvalues pα, can be found in [11, 96, 97, 134, 156]. For corner singularity analysis, the following expressions will be considered in the neighbourhood of the corner vertex fα(zα) = zλ αqα,fα+3(zα) = ¯zλ α¯qα,(2.20) with no summation over α, and where λis the dimensionless characteristic exponent, related with the stress singularity order δby δ=1−λ. To keep the strain energy bounded, and the resultant forces along the radial lines to be finite in the neighbourhood of the corner tip, the real part of λmust be non-negative, Re(λ)≥0. Rigid body translations correspond to λ=0. Although, according to the mathematical theory of linear elliptic systems [34, 77, 122], the general expression of fαmay involve also positive integer powers of a logarithmic term fα(zα) = zλ αlogpzαqα,with p=0,1,2... (2.21) CHAPTER 225 in the present work only power-law singularities (2.20) will be considered for the sake of simplicity. Substituting (2.20) in the representations (2.18) and (2.19) gives u=A⟨zλ ∗⟩q+¯ A⟨¯zλ ∗⟩˜ q,(2.22) φ=B⟨zλ ∗⟩q+¯ B⟨¯zλ ∗⟩˜ q,(2.23) where the column vectors of matrices Aand Bare the eigenvectors aαand bα, respectively, A=a1,a2,a3,B=b1,b2,b3,(2.24) and ⟨zλ ∗⟩denotes a diagonal matrix ⟨zλ ∗⟩=diag(zλ 1,zλ 2,zλ 3).(2.25) By using (2.5) we get zα=r(cos(θ)+ pαsin(θ)) = rζα(θ),(2.26) and substituting in (2.25) ⟨zλ ∗⟩=rλdiag(ζλ 1(θ),ζλ 2(θ),ζλ 3(θ)).(2.27) By substituting this expression in (2.22) and (2.23) we obtain u=rλ{A⟨ζλ ∗⟩q+¯ A⟨¯ ζλ ∗⟩˜ q},(2.28) φ=rλ{B⟨ζλ ∗⟩q+¯ B⟨¯ ζλ ∗⟩˜ q}.(2.29) These representations can be written in a compact way as w(r,θ) = rλXZλt,(2.30) where: w(r,θ) = u(r,θ) φ(r,θ),X=A¯ A B¯ B,Zλ(θ) = ⟨ζλ α(θ)⟩0 0⟨¯ ζλ α(θ)⟩,t=q ˜ q. (2.31) For expressions analogous to (2.30) and (2.31) for mathematically degenerate materials, see [11, 96, 97, 156]. 2.3 Singular elastic solution in a single-material wedge. Transfer matrix This section briefly presents the very useful transfer matrix concept introduced by Ting [135], see [11, 96, 97] for a more detailed analysis. Considering the special geometry of singularity problems for multi-material corners, the transfer matrix concept proposed by Ting [135] is very helpful to simplify the problem analysis. In the case of a single-material wedge number m, this matrix allows us to get a relation between displacements and the stress function vector, uand φ, on both faces of the wedge. 26 M.A. Herrera-Garrido Figure 2.1: Schema of a single-material wedge. That is, the transfer matrix Em(λ,θm,θm−1), gives a relation between the power-law solution values on the two faces of the wedge, wm(r,θm−1)and wm(r,θm), see Fig. 2.1, wm(r,θm) = Em(λ,θm,θm−1)wm(r,θm−1).(2.32) The transfer matrix Em, for mathematically non-degenerate materials, and for a given λ, depends on the properties of the wedge material and on the wedge geometry as Em(λ,θm,θm−1) = XZλ(θm,θm−1)X−1,(2.33) where Xis defined in (2.31), and Zλ(θm,θm−1) = ⟨ζλ α(θm,θm−1)⟩0 0⟨¯ ζλ α(θm,θm−1)⟩,(2.34) and ζα(θm,θm−1) = cos(θm−θm−1)+ pα(θm−1)sin(θm−θm−1),(2.35) with pα(θm−1) = pαcos(θm−1)−sin(θm−1) pαsin(θm−1)+cos(θm−1).(2.36) Analogous expressions for degenerate and extraordinary degenerate materials can be found in [11, 96, 97] Continuity problems can arise when evaluating the complex power function ζλ α(θm,θm−1) for which a proper choice of branch cut is necessary. Taking into account that 0 <θm− θm−1≤2π, and Im(pα(θ)) >0 for α=1,2,3, we use in our code the following argument function ψand radius ρ ψ=arg(ζα(θm,θm−1)) ∈(0,2π⟩,ρ=|ζα(θm−θm−1)|(2.37) in the (real variable) expansion of the complex power function ζλ α(θm,θm−1) = ρRe(λ)e−ψIm(λ)[cos(Re(λ)ψ+Im(λ)lnρ)+isin(Re(λ)ψ+Im(λ)lnρ)] (2.38) valid for any programming language. An equivalent possibility to avoid the continuity problem when coding in MATLAB is presented in [65, Section 4.4.1]. For the sake of notation simplicity, hereafter Em(λ,θm,θm−1)will be denoted as Em(λ). CHAPTER 227 2.4 Matrix formalism for boundary and interface conditions 2.4.1 Reference frame attached to a wedge face To define the boundary and interface conditions, it is first necessary to choose a suitable reference frame for displacement and traction vector, which facilitates the further development of the present matrix formalism. For this purpose, an orthonormal basis of vectors attached to the corner faces has been chosen as defined in [96, 97], see Fig. 2.2 (sr(ϑ),s3,n(ϑ)),(2.39) with Cartesian components sr(ϑ) =  −cosϑ −sinϑ 0 ,s3=  0 0 1 ,n(ϑ) =  −sinϑ cosϑ 0 .(2.40) x x 2 x 3 1 n ϑ ϑ0 s r n s (ϑ ) (ϑ ) (ϑ ) W ϑ 1 ϑ (ϑ ) s n(ϑ ) (ϑ ) 0 0 3 3 1 1 s r W-1 ss3 n(ϑ ) W-1 s3 s r (ϑ ) W W W r Figure 2.2: Orthonormal basis attached to each face These unique right-handed reference frames at each boundary face or interface (sr(ϑ), s3(ϑ),n(ϑ)) are defined for the imposition of boundary or interface conditions, for both displacements and stress function vectors. The vector sr(ϑ)should be used with caution in (2.41), where the vector sis considered always as counterclockwise on a solid boundary, whereas the direction of sron a wedge face can be clockwise or counterclockwise. 2.4.2 Boundary condition matrices As follows from the introduction section, the homogeneous boundary conditions considered in the present chapter fulfill the orthogonality condition between the boundary traction and displacement vectors. Therefore, there is a Cartesian coordinate system in which if a component of the stress function vector φis not zero then that component of the displacement vector uis zero, and vice versa. As the traction vector tnat a point of a wedge face, associated to the normal vector n defined in (2.40), can be calculated as tn(x1,x2) = −∂φ ∂s(x1,x2),(2.41) 34 M.A. Herrera-Garrido DBC(ϑw)and DI(ϑw)fulfill the orthogonality relations in (2.46) and (2.52), we can prove that the matrix Dcorner ext(ϑ)fulfills the following orthogonality relation Dcorner ext(ϑ)DT corner ext(ϑ) = I12W×12W,(2.65) where ϑ= (ϑ0,ϑ1,···,ϑW).(2.66) 2.5.2.1 Open corner Extended boundary and interface condition matrix. In the case of an open corner, the boundary and interface condition matrices, defined in Sections 2.4.2 and 2.4.3, are assembled in a real-valued 12W×12W extended matrix of boundary and interface conditions for an open multi-material corner Dcorner ext.(ϑ) = blocked diag[DBC(ϑ0),DI(ϑ1),...,DBC(ϑW)] =     DBC(ϑ0)06×12 06×12 ··· 06×6 012×6DI(ϑ1)012×12 ··· 012×6 . . .. . .. . ..... . . 06×606×12 06×12 ··· DBC(ϑW)     .(2.67) Recall that DBC and DIare 6 ×6 and 12 ×12 matrices, respectively. It is easy to prove that this matrix satisfies the orthogonality relation (2.65), as it is a block diagonal matrix whose main-diagonal blocks matrices (DBC(ϑ0),DI(ϑ1),...,DBC(ϑW)) are orthogonal. Open wedge characteristic matrix. As previously done for an elastic variable vector ww(r,ϑw)for a wedge win expression (2.47), when multiplying the matrix Dcorner ext. by the wcorner ext. vector, we get the wcorner PU which can be partitioned into sub-vectors of the prescribed and unknown variables Dcorner ext.(ϑ)wcorner ext. =wcorner PU,(2.68) where the 12W×1 vector wcorner PU =             wP(r,ϑ0) wU(r,ϑ0) wP(r,ϑ1) wU(r,ϑ1) . . . wP(r,ϑW) wU(r,ϑW)             ,(2.69) can be reduced to the 6W×1 vector of unknown variables, considering that the prescribed variables vanish, wcorner U =      wU(r,ϑ0) wU(r,ϑ1) . . . wU(r,ϑW)      .(2.70) CHAPTER 235 Note that the vectors wU(r,ϑw)for ϑ0and ϑWare 3×1 vectors, as they are the unknown variables for displacement and stress function vectors on the boundaries, while for the interfaces at ϑ1,...,ϑW−1they are 6×1 vectors. In view of the orthogonality property (2.65), relation (2.68) can be a rewritten as wcorner ext. =DT corner ext.(ϑ)wcorner PU,(2.71) which substituted into (2.64) leads to Kcorner ext.(λ)DT corner ext.(ϑ)wcorner PU =06W×1.(2.72) In this expression, some columns of the matrix obtained by multiplying Kcorner ext.(λ)by DT corner ext.(ϑ)can be removed as they would be multiplied by the prescribed zero values of wP(r,ϑw). This procedure leads to the final complex-valued square matrix with the reduced dimensions (6W×6W) Kcorner(λ) =      K1˜ DT BC(ϑ0)−˜ DT 1(ϑ1)06×6··· 06×606×3 06×3K2˜ DT 2(ϑ1)−˜ DT 1(ϑ2)··· 06×606×3 . . .. . .. . ..... . .. . . 06×306×606×6··· KW˜ DT 2(ϑW−1)−˜ DT BC(ϑW)     , (2.73) which introduced into (2.72) together with (2.70) gives the characteristic system for the singularity analysis of an open multi-material corner, which represents a nonlinear eigenvalue problem for this corner, Kcorner(λ)wcorner U =06W×1.(2.74) Examples of open corners composed by three wedges, i.e. W=3, are studied in Table 2.9. As it may be not so easily inferred from the general structure of Kcorner(λ)in (2.73) which is the structure of this matrix for a particular case of a corner with a small number of wedges W=1 or 2, these cases will be described in detail in the following. One-wedge corner. In the case of one wedge, W=1, Kcorner(λ)is a 6 ×6 matrix defined as Kcorner(λ) = K1(λ)˜ DT BC(ϑ0)−˜ DT BC(ϑ1),(2.75) where ˜ DT BC(ϑw), for w=0 and 1, are real 6×3 matrices. Actually, in the cases studied in the present work, where only orthogonal boundary conditions can be prescribed on both boundary surfaces, this can be even further reduced, cf. [11, 99, 153], Kcorner reduced(λ) =Du(ϑ1)K(1) 1(λ)˜ DT u(ϑ0)+ Du(ϑ1)K(2) 1(λ)˜ DT ϕ(ϑ0) +Dϕ(ϑ1)K(3) 1(λ)˜ DT u(ϑ0)+ Dϕ(ϑ1)K(4) 1(λ)˜ DT ϕ(ϑ0), (2.76) where K(i) 1are the 3 ×3 sub-matrices that form the 6 ×6wedge transfer matrix of the first and only wedge in this kind of corners K1(λ) = "K(1) 1(λ)K(2) 1(λ) K(3) 1(λ)K(4) 1(λ)#.(2.77) 36 M.A. Herrera-Garrido The characteristic system in this case reduces to Kcorner reduced(λ)wcorner U(r,ϑ0) = 03×1.(2.78) Results of singularity analysis for several examples of this kind of corners can be found in Tables 2.4 to 2.7. Two-wedge corner. For the corners composed by two wedges, W=2, the structure of the 12 ×12 corner characteristic matrix is Kcorner(λ) = K1˜ DT BC(ϑ0)−˜ DT 1(ϑ1)06×3 06×3K2˜ DT 2(ϑ1)−˜ DT BC(ϑ2).(2.79) Results of singularity analysis for several examples of two wedge corners can be found in Tables 2.8 and 2.9. 2.5.2.2 Closed corner (periodic corner) Extended interface condition matrix. The procedure to obtain the characteristic matrix of the corner is similar to the one presented in Section 2.5.2.1, but in contrast with open corners, in closed corners there are no boundary conditions and only interface conditions are prescribed. For this reason the extended matrix of interface conditions for a closed multimaterial corner is different from that for an open corner. The structure of this matrix for closed corners, defined by all the interface condition matrices, is somewhat more involved than for open corners, Dcorner ext.(ϑ) =         D2(ϑ0)06×12 06×12 ··· 06×12 D1(ϑW) ˜ D2(ϑ0)06×12 06×12 ··· 06×12 ˜ D1(ϑW) 012×6DI(ϑ1)012×12 ··· 012×12 012×6 . . .. . .. . ..... . .. . . 012×6012×12 012×12 ··· DI(ϑW−1)012×6         .(2.80) This is a 12W×12Wreal-valued matrix constructed in such a way that it fulfills the orthogonality condition (2.65), as shown in Appendix C. Remember that, as it is a closed corner, its first and last angle must meet the following relation: ϑW=ϑ0+360o. Closed corner characteristic matrix As for an open corner, when multiplying the matrix Dcorner ext. by the wcorner ext. vector, we get the 12W×1 vector wcorner PU wcorner PU =Dcorner ext.(ϑ)wcorner ext..(2.81) This vector can be partitioned in sub-vectors of the prescribed and unknown variables on boundary faces and interfaces wcorner PU =              wP(r,ϑ0) wU(r,ϑ0) wP(r,ϑ1) wU(r,ϑ1) . . . wP(r,ϑW−1) wU(r,ϑW−1)              .(2.82) CHAPTER 237 By omitting the prescribed variables in wcorner PU, the following 6W×1 vector wcorner U of unknown variables is obtained: wcorner U =      wU(r,ϑ0) wU(r,ϑ1) . . . wU(r,ϑW−1)      .(2.83) For closed corners, all the subvectors wU(r,ϑw)are 6×1 vectors. In view of the orthogonality relation (2.65), equation (2.81) can be rewritten as wcorner ext.=DT corner ext.(ϑ)wcorner PU.(2.84) Then, by substituting this equation in (2.64) we obtain Kcorner ext.(λ)DT corner ext.(ϑ)wcorner PU =06W×1.(2.85) In this expression some columns of the matrix given by multiplying Kcorner ext.(λ)by DT corner ext.(ϑ) can be removed as they would be multiplied by the prescribed zero values of wP(r,ϑw). This leads to a reduced 6W×6Wcomplex-valued matrix Kcorner(λ) =         K1˜ DT 2(ϑ0)−˜ DT 1(ϑ1)06×6··· 06×606×6 06×6K2˜ DT 2(ϑ1)−˜ DT 1(ϑ2)··· 06×606×6 . . .. . .. . ..... . .. . . 06×606×606×6··· KW−1˜ DT 2(ϑW−2)−˜ DT 1(ϑW−1) −˜ DT 1(ϑ0)06×606×6··· 06×6KW˜ DT 2(ϑW−1)         , (2.86) whose application in (2.85) together with (2.83) leads to the characteristic system for the singularity analysis of a closed multi-material corner, which represents a nonlinear eigenvalue problem for this corner, Kcorner(λ)wcorner U=06W×1.(2.87) Examples of closed corners composed by more than one wedge, i.e. W≥2, are studied in Table 2.11. Similarly as done for open corners, it may be not an easy task to infer from the general structure described in (2.86) the structure of Kcorner(λ)for the cases with only one wedge, i.e. W=1. Therefore, these special cases will be studied in detail in the following. One-wedge corner. Two cases can be considered in closed corners with only one wedge, i.e. W=1. The first case with all the materials perfectly bonded, and the second case with frictionless sliding on one of the interfaces. In the first case, as there is no boundary condition, Kcorner can be calculated as Kcorner(λ) = K1−I6×6.(2.88) This kind of corners is studied in Examples 2.3.1, 2.3.2 and 2.3.4 in Table 2.10. It is easy to see that in this case, singularity exponents in the range 0 <Re(λ)<1 may only exist when there are two or more dissimilar materials perfectly bonded, i.e. M≥2. In case that the corner is made by a single material, then λ=n, where nis an integer number. 38 M.A. Herrera-Garrido In the second case, when the corner is made by one single-material or multi-material wedge and one of its interface conditions is frictionless sliding contact, Kcorner can be calculated as Kcorner ext(λ)DT corner ext(ϑ0)wcorner PU(r,ϑ0) = 06×1,(2.89) K1(λ)−I6×6DT 2(ϑ0)˜ DT 2(ϑ0) DT 1(ϑ1)˜ DT 1(ϑ1) 0 wcorner U(r,ϑ0)=0,(2.90) K1(λ)DT 2(ϑ0)−DT 1(ϑ1)K1(λ)˜ DT 2(ϑ0)−˜ DT 1(ϑ1)0 wcorner U(r,ϑ0)=0,(2.91) K1(λ)˜ DT 2(ϑ0)−˜ DT 1(ϑ1)wcorner U(r,ϑ0) = 0,(2.92) Kcorner(λ) = K1˜ DT 2(ϑ0)−˜ DT 1(ϑ1).(2.93) Recall that ϑ1=ϑ0+360◦. These kind of corners are studied in Examples 2.3.3, 2.3.5 and 2.3.6 in Table 2.10. In the special case with the wedge given by a single material, λ=0.5+n, where nis an integer number. 2.5.3 Solution of the characteristic system. Singular elastic solution Avoiding the trivial solution for wcorner U, where wcorner U =0, any other solution of the corner eigenequation, (2.74) for open corners and (2.87) for closed corners, is a (right) null vector of the corner characteristic matrix Kcorner(λ). To get to this solution, we need first to find the characteristic (singular) values of λ, for which the complex-valued matrix Kcorner(λ) is singular with null determinant. Therefore, a straightforward method to find the singular values is finding the roots of the matrix determinant detKcorner(λ) = 0.(2.94) Up to this point, all the steps of the present procedure have been analytical except for the procedure for computing the roots of the Lekhnitskii-Stroh sextic characteristic polynomial of an anisotropic linear elastic material in those cases where its roots cannot be expressed in terms of radicals [19, 58]. In some specific cases, this sextic equation can be solved analytically, e.g., for all transversely isotropic materials with any spatial orientation [12, 129], and also for some classes of orthotropic materials [19]. Noteworthy, as follows from the procedure presented, all elements of the matrix Kcorner(λ) are complex analytic (holomorphic) functions of λ, thus, also the determinant of this matrix is a complex analytic function of λ. Nevertheless, an analytic solution of (2.94) is not possible in general case, except for very specific simple cases. Even in the case of a single isotropic material, in plane strain, with free boundary faces, (2.94) is a transcendental equation requiring a numerical solution [139]. For this reason, we use a numerical method called Muller method [104] to solve it. This is a standard iterative procedure suitable for searching for complex roots of complex analytic functions. Argument principle. Incorporated within the code is the Argument Principle [72], which aids in locating all the roots of the complex-valued analytic function detKcorner(λ)within a specified region of the complex plane. The procedure indicates the number of roots (real and complex) in a region bounded by a closed contour C. This method counts the number of times that the complex value of the analytic function detKcorner(λ)rotates around the origin CHAPTER 239 of the co-ordinate system, i.e. what is the multiple of 360◦by which the argument of this function increases or decreases along the contour C, J=1 2π[arg(detKcorner(λ))].(2.95) 2.6 Stresses and displacements Once that the characteristic exponents λof a corner have been found, we can calculate the characteristic stress and displacement fields. First, we have to substitute the obtained λ value into expression (2.74) or (2.87) considering a fixed value of the radial coordinate r, e.g. r=1. The solution for wcorner U, for this fixed value of r=1, is a (right) null vector of the matrix Kcorner(λ). As the obtained value for λis substituted in Kcorner(λ)we get a numerically defined matrix. To compute the null eigenvector of the numerically defined matrix Kcorner(λ)we use the Singular Value Decomposition (SVD) giving the right singular vector of Kcorner(λ)associated to the minimum singular value σmin ⪆0, which corresponds to wcorner U for r=1. The reason for using SVD is that sometimes we have found stability problems when computing the null eigenvector of Kcorner(λ), whereas SVD has shown to be a very robust procedure for this purpose. Next step is to complete wcorner U with zero values of wcorner P to get wcorner PU. Care must be taken with the size of the subvectors to be added, since as said before, for closed corners all subvectors wUare 6 ×1 vectors and so are the corresponding subvectors wP, but for open corners, the first and last subvectors wUare 3 ×1 vectors and the rest 6×1 vectors, the corresponding subvectors wPhave the same size. Introducing wcorner PU in expression (2.71) or (2.84), we obtain wcorner ext. This vector contains the displacement and stress function for both faces of each wedge, ww(r,ϑw−1)and ww(r,ϑw). The w1(r,ϑ0)vector for the first face of the first wedge, corresponds directly to the w1(r,θ0)vector of the first face of the first material of the first wedge. Using now the transfer matrix in (2.32), we get w1(r,θ1)of the second face of the first material of the first wedge. As within a wedge we only have perfectly bonded interfaces, we know that the w2(r,θ1)vector has the same values as w1(r,θ1). Continuing with this process we can compute wm(r,θm−1) for the first face of each material in the corner. Now, starting with each wm(r,θm−1)and using an analogous expression to (2.32), it is easy to get wm(r,θ)for each θwithin each material, wm(r,θ) = Em(λ,θ,θm−1)wm(r,θm−1)for θm−1≤θ≤θm.(2.96) Following the steps proposed by Ting ([134], Section 7.3) the displacement vectors and stress tensor in cylindrical coordinates are: ur=−sT r(θ)u(r,θ),uθ=nT(θ)u(r,θ),u3=sT 3(θ)u(r,θ),(2.97) and σrr =sT r(θ)φ,θ(r,θ) r,σθθ =nT(θ)φ,r(r,θ), σrθ=−nT(θ)φ,θ(r,θ) r=−sT r(θ)φ,r(r,θ), σr3=−sT 3(θ)φ,θ(r,θ) r,σθ3=sT 3(θ)φ,r(r,θ),(2.98) 40 M.A. Herrera-Garrido where the comma denotes differentiation. We can get φ,rand φ,θfrom (2.30). On the one hand, as E,X,Zand tdo not depend on rwe have ∂ ∂rw(r,θ) = λrλ−1Em(λ,θ,θm−1)XZλ(λ,θ,θm−1)t(2.99) for φ,r. On the other hand, as Edepends on θwe have, for non-degenerate materials, ∂ ∂θ w(r,θ) = ∂ ∂θ Em(λ,θ,θm−1)wm(r,θm−1),(2.100) where ∂ ∂θ Em(λ,θ,θm−1) = X∂ ∂θ Zmλ(λ,θ,θm−1)X−1,(2.101) ∂ ∂θ Zmλ(λ,θ,θm−1) = λZmλ−1(λ,θ,θm−1)∂ ∂θ Zm(λ,θ,θm−1)(2.102) ∂ ∂θ Zm(λ,θ,θm−1) = ⟨ζ′ α(θ,θm−1)⟩0 0⟨¯ ζ′ α(θ,θm−1)⟩,(2.103) ζ′ α(θ,θm−1) = ∂ ∂θ ζα(θ,θm−1) = −sin(θ−θm−1)+ pα(θm−1)cos(θ−θm−1).(2.104) These expressions allow us to compute stresses and displacements as functions of θfor a given value of r. An example to illustrate this procedure can be found in Section 2.8.3. 2.7 MATLAB implementation The semianalytic procedure described above has been implemented in MATLAB [101] using Symbolic Math Toolbox. This code calculates the singularity exponents and plots displacement and stress singular fields associated to a corner problem. This code is organized in 6 modules: 1. Data input 2. Definition of single-material wedges 3. Boundary and interface condition matrices 4. Characteristic system assembly 5. Solution of the characteristic system 6. Displacement and stress singular fields 2.7.1 Data input The first module reads the data entered by the user and performs the necessary calculations. This module gives the user two options, to enter all the values that define the problem through a text file or interactively. CHAPTER 241 2.7.2 Material definition This module is subdivided into three different functions, each function is for a type of material currently considered in the code. The options are isotropic, transversely isotropic or orthotropic. The code can be easily generalized to any other class of anisotropic materials, covering both, mathematically non-degenerated and degenerated materials following [11, 96, 97]. From the elastic constants and the initial and final angles for a single-material wedge number m,θm−1and θmrespectively, the corresponding function will apply the Stroh formalism and store the transfer matrix Em(λ)(2.33) as the definition of the single-material wedge. This module is based on the theory described in Sections 2.2 and 2.3. Special attention was paid to the manipulation of angles and complex numbers to ensure the continuity of analytic functions as discussed in Section 2.3. In MATLAB the function used to obtain the argument of a complex number is ‘angle()’ that returns the phase angle in the interval ⟨−π,π⟩. For the correct execution of (2.38), the following definition of ‘arg()’ function is used arg(z)def =(angle(z)+ 2πif angle(z)≤0 angle(z)if angle(z)>0 (2.105) 2.7.3 Boundary and interface condition matrices The first step of this module is the creation of the right-handed reference frames (sr,s3,n) attached to each face of the wedges that form the corner, see Fig. 2.2, and also (˜ m,˜ l,˜ n)for the cases of inclined planes shown in Fig. 2.3. Once we have the reference frames defined for each angle at which a boundary or interface condition is prescribed, the boundary and interface condition matrices are generated depending on the type of boundary or interface condition, see Tables 2.1 and 2.2. 2.7.4 Characteristic system assembly In this part, the Kwmatrix for each multi-material wedge (2.60) will be formed by the Em matrices (2.33) of the materials that form it. These Kwmatrices are combined with matrices DBC in (2.45) or DIin (2.51) generating the characteristic matrix Kcorner(λ)for open corners (2.73) or for closed corners (2.86). Special attention is paid to the size of the rectangular and square matrices used in assembling of the characteristic system. 2.7.5 Solution of the characteristic system This module is based on Section 2.5.3. Once we have the characteristic matrix of the system Kcorner(λ), (2.73) for open corners or (2.86) for closed corners, we initially define the interval of real values of λwhere the roots of detKcorner(λ)are searched for. This is the first part of the code where solution is numerical instead of analytical, since different values of λ will be substituted in the analytical expression of Kcorner(λ)matrix to solve numerically the determinant. If the code detects that there may be a root between two of the λvalues for which the determinant has been calculated, it will try to find a root in that interval using the Muller method [104]. After showing the roots automatically found, the software will show a graphic representation of Re(detKcorner(λ)) that could help the user to check the roots 42 M.A. Herrera-Garrido automatically found and indicate, in case that some root is missed, if new starting points for searching more roots are necessary. To illustarte this, the points where the curve intersects the abscissa axis in Fig. 2.4 are the points where the user should look for roots in case that the automatic procedure had not taken them into account. The code also gives the option to run the submodule ‘Argument Principle’ to detect possible roots (real or complex) in a region of complex plane that have not been detected automatically. 2.7.6 Displacement and stress singular fields At this point the code will perform the steps described in Section 2.6 to produce a graph like the one shown in Fig. 2.5. As can be seen from the formulas described in (2.98), for the calculation of stresses it is necessary to differentiate the stress function vector φwith respect to θor r. To speed up the calculation, the analytical expression of each derivative has been used instead of leaving it to the software to derive it in each case. In the case of the derivative with respect to r, it is very simple, since according to the form of the expression of φshown in (2.29), the differentiation with respect to rproduces the same function φmultiplied by the root λ, considering r=1, see (2.99). Slightly more complicated is the differentiation with respect to θ, for which we will have to differentiate with respect to θthe ‘transfer matrix’ Em, which is the only part of the expression dependent on θ. The expression for the derivative of Emis computed directly from (2.101)-(2.104). Special care must be taken regarding the compatibility of dimensions of many vectors and matrices used in the procedure. In the case a parameterization with respect to a corner parameter is requested, as in the example in Section 2.8.4, the code will repeat the modules 2 to 5 as many times as necessary. 2.8 Examples To verify the correct implementation of the present procedure for corner singularity analysis in the developed code, many examples are solved and the results are compared with the numeric solution of the closed-form eigenequations available in the literature. In Table 2.3 we list, in the full precision considered, the characteristics of the used materials. For orthotropic materials, the fibers in the x1x3-plane are oriented by angle φwith respect to the x2-axis. This angle is specified for each case of orthotropic material. Table 2.3: Engineering constants for the materials used for the studied examples. Shear and elastic moduli in GPa. Material E1E2E3G12 G13 G23 ν12 ν13 ν23 φ A 68.67 0.33 B 3 0.35 C 137.9 14.48 1 5.86 1 1 0.21 0 0 0 D 137.9 14.4795 14.4795 5.86 5.86 5.86 0.21 0.21 0.21 φ E 137.9 14.48 14.48 4.98 4.98 4.98 0.21 0.21 0.21 φ F 5.85 0.25 G 141.3 9.58 9.58 5 5 3.5 0.3 0.3 0.32 φ The studied problems are divided into open and closed corners and further subdivided into single-wedge and multi-wedge corners. The following acronyms are used: ‘BC’ Boundary CHAPTER 243 Condition, ‘IC’- Interface Condition, ‘B’ perfectly Bonded, ‘FL’ FrictionLess, ‘F’ stress Free, ‘C’ Clamped, ‘S’ Symmetry (frictionless contact), ‘A’ Antisymmetry, ‘(SY)’ SYmmetric loading, ‘(SK)’ SKew-symmetric loading, and ‘(A’) Antiplane shear. In the tables in this section, Example is abbreviated as ’Ex.’ and material as ’mat.’ and the last column is CPU time in a workstation (DELL Precision 5550). In most of the studied problems, the results found in the literature are usually for plane strain or plane stress only, whereas the presented formalism works in generalized plane strain. This means that this formalism in some cases finds solutions corresponding to the anti-plane shear that are not provided by some methods found in the literature. In those cases, in the corresponding table, the result obtained by the presented formalism is included beside a void cell corresponding to a result not covered by the eigenequation or closed-form expression for λfound in the literature 2.8.1 Solutions for open corners 2.8.1.1 Solutions for open corners with only one single-material or multi-material wedge The expression for Kcorner(λ)for this special case was deduced in Section 2.5.2.1. Single-material wedge. If the wedge is made of only one material, then K1(λ) = E1(λ) in (2.75). For isotropic materials, with stress free-stress free boundary conditions, results for λare compared with those by Vasilopoulus [139] with perfect match. For isotropic materials with different homogeneous boundary conditions, we compare with the results obtained by the closed-form expressions collected by Sinclair [121]. In Table 2.4 we show some of these results. Table 2.4: Comparison with the results found in the literature for open single-material isotropic corners with different boundary conditions,Section 2.8.1.1. Corner configuration Ex. BC1 BC2 Mat Results in literature [121] Present results Time (sec.) x 2 1 x ϑ0 ϑ 1 BC 1 BC 2 x 3 2.1.1 F / 0º F / 270º A 0.5444837368 (SY) 0.5444837368 2.5 0.6666666667 (A) 0.6666666667 0.9085291898 (SK) 0.9085291898 2.1.2 F / 0º C / 270º 0.3333333333 (A) 0.3333333333 3.1 0.8607568402 0.8607568402 2.1.3 C / 0º C / 270º 0.5904563986 (SK) 0.5904563986 2.9 0.6666666667 (A) 0.6666666667 0.7673218225 (SY) 0.7673218225 2.1.4 FL / 0º FL / 270º 0.3333333333 (SY) (SK) 0.3333333333 3.1 0.6666666667 (A) 0.6666666667 2.1.5 FL / 0º F / 270º 0.3333333333 0.3333333333 2.9 0.6666666667 0.6666666667 2.1.6 FL / 0º C / 270º 0.3333333333 0.3333333333 3.0 0.6666666667 0.6666666667 50 M.A. Herrera-Garrido sr n mñ γ γ 0 20 40 60 80 100 120 140 160 180 -0.2 0 0.2 0.4 0.6 0.8 1 Figure 2.6: Variation of the singularity exponent λvs the angle γ, which defines the vector ˜ nin the boundary condition ‘only u˜nrestricted’, increases on a single-material wedge of the orthotropic material ‘D’ with φ=0 (see Table 2.3). Figure 2.7: Singularity exponent λfor the tri-material junction with perfectly bonded interfaces and E3/E1=100, overlayed with [23, Fig. 9a] with perfect match of all results. Additionally, the computational code reveals a solution corresponding to the anti-plane mode (red-circles), not addressed in [23]. 2.9 Concluding remarks This chapter revisits and enhances the matrix formalism and computational procedure initially introduced in a concise manner in [96, 97], aiming for improved clarity for readers. However, the primary contribution of this chapter lies in its presentation of the first comprehensive implementation of the proposed procedure, extensively validated through numerous numerical examples. This implementation serves to validate the performance of the proposed general matrix formalism as well. The code underwent rigorous validation via multiple tests, comparing the obtained results with those presented by other authors. Thus, although frictional contact has not been addressed in this chapter, the Objective 1 of this thesis is partially achieved. The subsequent step involves replicating the methodology outlined in this chapter, with a specific CHAPTER 251 focus on frictional contact. Regarding Objective 1.3, which entails verifying both the matrix formalism and its implementation in a computational code, and Objective 1.4, involving the verification of methodologies developed by other authors in the literature, it is noteworthy that significant discrepancies were identified in one methodology. Specifically, this occurred in the context of frictionless contact between two anisotropic materials [28]. These discrepancies will be readdressed in the following chapter, which deals with cases involving frictional contact. Nevertheless, as evidenced by the tables of results, excellent agreement has been attained in the majority of cases. Only minor typographical errors in some previous publications required attention. It is observed that certain results obtained through the general matrix formalism are absent when employing the compared previous and specific methodologies. This discrepancy often arises because the solutions correspond to an antiplane mode, whereas the majority of methodologies found in the literature only yield results for the in-plane mode. CHAPTER 3 General matrix formalism for anisotropic multi-material corners with frictional boundary and interface conditions 3.1 Introduction The aim of this chapter, based on [60], is to add to the previously introduced matrix formalism and its MATLAB implementation in Chapter 2 the ability to study cases involving frictional sliding contact, either at the boundary of the multi-material corner or between any of its wedges. To achieve this, we will thoroughly revise and complete the matrix formalism for the singularity analysis of linear elastic anisotropic multi-material corners with frictional sliding contact under GPS, as succinctly introduced in [96, 97]. Secondly, we aim to develop its general computational implementation and test its performance and accuracy by comparing the obtained results with those available in the relevant literature. This should allow us, in the future, to solve virtually any singularity problem for anisotropic multi-material corners including frictional sliding contact and covering any number and type of homogeneous linear elastic materials in each wedge, any type of homogeneous boundary and interface conditions, any geometry of the corner with radial faces. This generalisation requires a substantial modification of the previous matrix formalism developed and implemented in Chapter 2, by including additional unknown variables giving angles ωwof the frictional tangential component of the stress vector and the relative tangential displacement vectors at the faces or interfaces with frictional contact, in view of their collinearity in the isotropic frictional contact. This modification leads to the solution of a non-linear eigensystem with a rectangular matrix for the singularity exponent λand angles ωw. Several numerical procedures are proposed to solve this non-linear eigensystem and their accuracy is tested. The present semi-analytical matrix formalism is a general mathematical tool capable of solving stress singularity problems under GPS for anisotropic multi-material corners with 53 54 M.A. Herrera-Garrido frictional sliding contact, and with the local geometry near the corner apex and boundary and interface conditions as described in Fig. 1.1. The study of open or closed (periodic) corners, from the most simple case of single-material and single-wedge case to the multi-material wedge and multi-wedge case of multi-material corners, is covered by this formalism. As in Chapter 2, only power singularities are considered, for the sake of simplicity. Thus, special cases of logarithmic or power-logarithmic singularities [38, 121, 122] are not covered by this formalism, although its semi-analytic nature would allow the study of such singularities if necessary. Theoretically, the geometry of each wedge in the corner is defined by a semiinfinite angular sector. This chapter is organised as follows: in Section 3.2, the matrices for frictional sliding boundary and interface conditions, one of the main novelties of this work, are introduced. The characteristic matrix of the non-linear eigensystem is assembled using the single-material wedge transfer matrices and the boundary and interface condition matrices in Section 3.3. The implementation of the code is summarised in Section 3.4, which develops some points for which the reader may need further clarification. The numerical results in Section 3.5 show perfect agreement with the results for isotropic corners and orthotropic corners, with one PES parallel to the contact plane and another PES perpendicular to the x3-axis, found in the literature, but significant differences are sometimes obtained with the results obtained by other authors for orthotropic corners with a general orientation of PES, the reasons for these differences are discussed. Finally, in Section 3.6 some concluding remarks are given. 3.2 Matrix formalism for boundary and interface conditions 3.2.1 Reference frame attached to a corner face To describe the isotropic frictional sliding contact within the framework of the present matrix formalism, it is first necessary to define two orthonormal bases for the sliding and frictional shear stress directions, making use of the orthonormal basis of vectors attached to the corner faces defined in (2.39) and (2.40), which make it easy to impose the necessary conditions, see Fig. 3.1, (k(ϑ,ω),m(ϑ,ω),n(ϑ)) (3.1) (nµ(ϑ,ω,µ),sµ(ϑ,ω,µ),m(ϑ,ω)),(3.2) with k(ϑ,ω) = cos(ω)sr(ϑ)+sin(ω)s3,(3.3) m(ϑ,ω) = −sin(ω)sr(ϑ)+cos(ω)s3,(3.4) nµ(ϑ,ω,µ) = n(ϑ)+ µk(ϑ,ω) p1+µ2,(3.5) sµ(ϑ,ω,µ) = k(ϑ,ω)−µn(ϑ) p1+µ2,(3.6) where µ≥0 is a given constant value of the kinetic (or dynamic) friction coefficient on the face of the wedge, and the unknown value of the angle ωgives the direction of the tangential CHAPTER 355 component τof the stress vector1tassociated with the normal vector n(ϑ), on the face at the angle ϑ, in the sense that τ=−|τ|k(ϑ,ω),(3.7) Therefore, the unit vector k(ϑ,ω)determines the direction of the tangential stress vector τ. Vectors nµand sµare used in Sections 3.2.3 and 3.2.5 and the reason for choosing these vectors is explained in Appendices D and E. The subscript whas been omitted in the above vector definitions, for the sake of simplicity, but each wedge face in a multi-material corner under friction can have a different value of µand ω. (a) sμ s3 μ s nμ s r s3 n mk n m sk n w μ r 0 0 ω w ω ω ω (b) nμ sμ n k (c) s nμ mks r n n n s s m μ μ μ r k (d) m k n s r s r mk sμ sμ nμ nμ s 3 ω s 3 ω w w ω 0 ω 0 Figure 3.1: Orthonormal bases attached to each face of a multi-material corner for frictional sliding contact. a) 3D view of an open corner with orthonormal basis on its outer faces, b) vectors in the plane defined by kand n, c) orthogonal projection on the x1x2-plane, d) orthogonal projection on the x1x3-plane. 3.2.2 Signorini-Coulomb conditions for frictional sliding contact on a boundary The frictional contact model adopted in the present work is the well-known Amontons- Coulomb friction model (or simply the Coulomb model). For frictional sliding contact as 1Usually referred to as the tangential stress vector. 56 M.A. Herrera-Garrido a boundary condition for the face number w=0 or W, a null displacement must be imposed in the normal direction to the face. uθ(r,ϑw) = 0.(3.8) In this work, only µ≥0 is considered, since a negative value of the friction coefficient considered by some other authors [30, 111, 144, 154] corresponds in the present formulation with a sliding direction in the opposite sense to that obtained with a positive friction coefficient, that is, with angles ω±180◦and ω, respectively. The shear stress vector τ(r,ϑw)and the normal stress σθθ (r,ϑw)are related by the Coulomb friction coefficient µ |τ(r,ϑw)|=−µσθθ (r,ϑw),(3.9) where the shear stress vector on the faces with the number w=0 or Wis τ(r,ϑw) = −σθrsr(ϑw)+σθ3(r,ϑw)s3.(3.10) The last condition for frictional sliding contact is that the (tangential) displacement vector, ui(r,ϑw), and the tangential stress vector τi(r,ϑw)must be parallel, idenoting the wedge number and wdenoting the face number, τi(r,ϑw) = −κui(r,ϑw)with (κ≤0 for i=1 and w=0 κ≥0 for i=Wand w=W (3.11) choosing the sign of κso that the direction of the physical tangential stress is opposite to the direction of sliding to satisfy the condition of frictional dissipation of energy. In Appendix D, these contact conditions are related to the reference frames presented in Section 3.2.1 for a better comprehension of the matrices presented in Section 3.2.3. The non-linear frictionless or frictional sliding contact conditions are first imposed as linear boundary or interface conditions, which is one of the main ideas behind this formalism. Then, after solving the corner eigenproblem, each obtained eigensolution is checked for the fulfilment of the compression condition on all contact surfaces, and for the frictional energy dissipation condition. Only eigensolutions verifying these non-linear conditions are valid. 3.2.3 Boundary condition matrices for frictional contact The frictional contact boundary conditions described in Section 3.2.2 are imposed by using suitable boundary condition matrices defined in the following. In contrast to the boundary conditions in Table 2.1, the matrices Duand Dϕfor the frictional contact boundary condition introduced in Table 3.1 are non-square real matrices and do not fulfil the same orthogonality relation. These matrices are rectangular, having one row more than the matrices for orthogonal boundary conditions. This extra row leads to four boundary conditions instead of three. This additional boundary condition is associated with the fact that a certain direction of shear component of the traction vector ωis assumed here, although, in general, the value of ωis unknown a priori and is computed by solving the corner eigenvalue problem. To express the boundary conditions of frictional contact, for boundary faces w=0 and W, as a linear relation in a way similar to that in (2.42) Du(ϑw,ωw)u(r,ϑw)+Dϕ(ϑw,ωw,µw)φ(r,ϑw) = 0,(3.12) CHAPTER 357 Table 3.1: Matrices Duand Dϕfor Coulomb sliding friction. Du(ϑ,ω)=[n(ϑ),0,m(ϑ,ω),0]TDϕ(ϑ,ω,µ)=[0,sµ(ϑ,ω,µ),0,m(ϑ,ω)]T ˜ Du(ϑ,ω)=[0,k(ϑ,ω)]T˜ Dϕ(ϑ,ω,µ)= [nµ(ϑ,ω,µ),0]T the following hypothesis must be assumed: the whole face is sliding in the same direction, with ωand µconstant on the whole face, and a monotonic loading is applied from the unloaded state. In this way, the non-linearities associated with the unilateral Signorini conditions of non-penetration and the Coulomb friction law are avoided. The main boundary-condition-matrix,DBC introduced in (2.45) for the homogeneous boundary condition, is expressed for the frictional sliding condition matrices in (3.12) as, see Table 3.1, DBC(ϑw,ωw,µw) = Du(ϑw,ωw)Dϕ(ϑw,ωw,µw) ˜ Du(ϑw,ωw)˜ Dϕ(ϑw,ωw,µw).(3.13) Similarly, as discussed in Section 2.4.2, in the frictional contact case, the matrices in Table 3.1 are not uniquely defined. It is easy to show that the matrix DBC is a real square matrix 6 ×6 that fulfils the orthogonality relation, cf. [11, 96, 97, 99], DBC(ϑw,ωw,µw)DT BC(ϑw,ωw,µw) = DT BC(ϑw,ωw,µw)DBC(ϑw,ωw,µw) = I6×6.(3.14) The orthogonality relation in (3.14) is equivalent to the following relations for the matrices in Table 3.1: DuDT u+DϕDT ϕ=I4×4,(3.15) Du˜ DT u+Dϕ˜ DT ϕ=04×2,(3.16) ˜ DuDT u+˜ DϕDT ϕ=02×4,(3.17) ˜ Du˜ DT u+˜ Dϕ˜ DT ϕ=I2×2.(3.18) For the sake of notation simplicity, the argument parameters ϑw,ωw, and µware omitted. When multiplying from the left the vector w(r,ϑw)by the matrix DBC in (2.45) or (3.13), the prescribed and unknown components of w(r,ϑw)appear grouped in two separate blocks, wP(r,ϑw)and wU(r,ϑw), respectively, DBC(ϑw,ωw,µw)w(r,ϑw) = wP(r,ϑw) wU(r,ϑw).(3.19) The vectors wPand wUare 3×1 vectors for the homogeneous (orthogonal) boundary conditions in Table 2.1, but, in case of frictional sliding contact, wPis a 4 ×1 vector and wUis a 2×1 vector. Taking into account the orthogonality relation in (3.14) and the fact that wP=0, w(r,ϑw) = DT BC(ϑw,ωw,µw)wP(r,ϑw) wU(r,ϑw)= ˜ DT u(ϑw,ωw) ˜ DT ϕ(ϑw,ωw,µw)wU(r,ϑw) = ˜ DT BC(ϑw,ωw,µw)wU(r,ϑw) (3.20) for w=0 and W, similarly to (2.48). 58 M.A. Herrera-Garrido 3.2.4 Signorini-Coulomb conditions for frictional sliding contact on an interface For frictional sliding contact as an interface condition for faces at angle ϑw, with 1 ≤w≤ W−1, for open corners, and 0 ≤w≤Wfor closed (periodic) corners, a null relative normal displacement must be imposed between the two wedges in contact; ∆uθ(r,ϑw) = uw θ(r,ϑw)−uw+1 θ(r,ϑw) = 0.(3.21) The tangential stress vector τ(r,ϑw)(3.10) and the normal stress are related by the Coulomb friction coefficient µas in (3.9). Additionally, a continuity condition of tractions (in fact an equilibrium of tractions) must be imposed for the materials in contact, tw(ϑw)−tw+1(ϑw) = 0.(3.22) At the interface, the relative tangential displacement vector is defined as ∆u(r,ϑw) = −∆ur(r,ϑw)sr(r,ϑw) + ∆u3(r,ϑw)s3(r,ϑw),(3.23) where ∆ur(r,ϑw) = uw r(r,ϑw)−uw+1 r(r,ϑw), ∆u3(r,ϑw) = uw 3(r,ϑw)−uw+1 3(r,ϑw).(3.24) To satisfy the condition of frictional dissipation of energy, the tangential stress vector τi, acting on the face of the wedge i=w,w+1, must be parallel to ∆u(r,ϑw)but with the opposite direction. Therefore, τi=−κ∆u(r,ϑw)with κ≥0,(3.25) or component-wise τi r(r,ϑw) = σi θr(r,ϑw) = −κuw r(r,ϑw)−uw+1 r(r,ϑw), τi 3(r,ϑw) = σi θ3(r,ϑw) = −κuw 3(r,ϑw)−uw+1 3(r,ϑw).(3.26) In Appendix E, these contact conditions are related to the reference frames presented in Section 3.2.1 for a better comprehension of the matrices presented in Section 3.2.5. 3.2.5 Interface condition matrices for frictional contact The frictional contact interface conditions described in Section 3.2.4 are imposed by using suitable interface condition matrices defined in the following. Similarly to the frictional sliding boundary condition, in the case of the frictional sliding contact interface condition, a certain friction shear direction (or sliding direction) given by the angle ωis assumed here, which means that there are seven interface conditions instead of the expected six conditions. Therefore, the matrices D1and D2, presented in Table 3.2, are rectangular 7 ×6 real matrices, while ˜ D1and ˜ D2are rectangular 5 ×6 real matrices. As in the previous section, the expressions of these matrices in Table 3.2 are just a convenient example of possible options to impose the frictional interface contact. CHAPTER 359 Table 3.2: Matrices D1,D2,˜ D1and ˜ D2for frictional interface contact conditions. D1(ϑ,ω,µ) = 1 √2−n(ϑ)03x103x3−m(ϑ,ω)03x1 03x1sµ(ϑ,ω,µ)−I3x303x1m(ϑ,ω)T Frictional D2(ϑ,ω,µ) = 1 √2n(ϑ)03x103x3m(ϑ,ω)03x1 03x1sµ(ϑ,ω,µ)I3x303x1m(ϑ,ω)T interface ˜ D1(ϑ,ω,µ) = 1 √2n(ϑ)03x1√2k(ϑ,ω)03×1m(ϑ,ω) 03×1nµ(ϑ,ω,µ)03×103×103×1T ˜ D2(ϑ,ω,µ) = 1 √2n(ϑ)03x103×1√2k(ϑ,ω)m(ϑ,ω) 03×1nµ(ϑ,ω,µ)03×103×103×1T To express the frictional contact interface condition as a linear relation for vectors of elastic variables, the same hypothesis as adopted in the previous section for the friction boundary condition must be applied, leading to D1(ϑw,ωw,µw)ww(r,ϑw)+D2(ϑw,ωw,µw)ww+1(r,ϑw) = 0.(3.27) The main interface-condition matrix defined for the interface frictional contact condition is a 12 ×12 real matrix DI, analogous to the matrix DBC for boundary conditions. For perfectly bonded interface and frictionless sliding condition, this matrix is defined as (2.51), whereas for the frictional contact interface as DI(ϑw,ωw,µw) = D1(ϑw,ωw,µw)D2(ϑw,ωw,µw) ˜ D1(ϑw,ωw,µw)˜ D2(ϑw,ωw,µw).(3.28) Again, it can easily be checked that the matrix DI(ϑw,ωw,µw)in (3.28) fulfils the orthogonality relation DI(ϑw,ωw,µw)DT I(ϑw,ωw,µw) = DT I(ϑw,ωw,µw)DI(ϑw,ωw,µw) = I12×12.(3.29) This orthogonality relation is analogous to the orthogonality relations for a perfectly bonded interface and a frictionless sliding interface deduced in (2.52), and is equivalent to the following relations for the matrices in Table 3.2 for the frictional sliding interface D1DT 1+D2DT 2=I7×7, D1˜ DT 1+D2˜ DT 2=07×5, ˜ D1DT 1+˜ D2DT 2=05×7, ˜ D1˜ DT 1+˜ D2˜ DT 2=I5×5. (3.30) For the sake of notation simplicity, the argument parameters ϑw,ωw, and µware omitted. Multiplying a 12 ×1 vector given by two 6 ×1 vectors w(r,ϑ)defined on the face wof the two consecutive wedges, wand w+1, from the left by the matrix DI, the vectors wP(r,ϑ) and wU(r,ϑ)of the prescribed and unknown elastic variables, respectively, are obtained DI(ϑw,ωw,µw)ww(r,ϑw) ww+1(r,ϑw)=wP(r,ϑw) wU(r,ϑw),(3.31) where wPis a 7×1 vector and wUis a 5×1 vector. Then, by multiplying (3.31) from the left by DT I, using the orthogonality relation (3.29) and taking into account that wP=0we get a very relevant expression of the interface elastic variables in terms of wU, ww(r,ϑw) ww+1(r,ϑw)=DT I(ϑw,ωw,µw)wP(r,ϑw) wU(r,ϑw)=˜ DT 1(ϑw,ωw,µw) ˜ DT 2(ϑw,ωw,µw)wU(r,ϑw).(3.32) 66 M.A. Herrera-Garrido 3.3.3 Solution of the characteristic system. Singular elastic solution The sought solution of the non-linear eigensystems in (3.49) or (3.57) is a (right) null vector of the matrix Kcorner(λ,ω), ignoring the trivial solutions with wcorner U =0. This means that the next step is to find the characteristic (singular) values of λand ω. As detailed in Chapter 2, in the simplest case, for corners where F=0, the matrix Kcorner is a square matrix and λis the only unknown. To obtain the λvalues, solutions of the characteristic equation of the corner, the roots of the determinant of the matrix Kcorner can be found detKcorner(λ) = 0.(3.58) In general the roots of (3.58) can be real or complex. If λis a complex root of (3.58) its complex conjugate is also a root of (3.58). For corners where F≥1, a pair of values of λand ω, non-trivial solution of the overdetermined homogeneous systems (3.49) or (3.57), for open or closed corners, respectively, should be found. In the following, three different ways to get the values of λand ωare proposed. In the first method the characteristic equation of the corner is written as a condition for a vanishing determinant det(¯ KT corner(λ,ω)Kcorner(λ,ω)) = 0,(3.59) which ensures that the columns of Kcorner are linearly dependent. This procedure has the advantage that the solution can be found in a similar way as in (3.58), but a disadvantage is that the algebraic multiplicity of each root is doubled. The second method is to evaluate the singular value decomposition (SVD) [54] of the matrix Kcorner(λ,ω)in (3.49) or (3.57), and to search for the vanishing minimum singular value by solving σmin(Kcorner(λ,ω)) = 0.(3.60) The last proposed method is to solve the following non-linear system of 6Wequations f(X) = Kcorner(λ,ω)wcorner U =0(3.61) for (6W+1)×1 vector of unknowns X=  wcorner U λ ω .(3.62) This non-linear system could be augmented by a normalisation condition for wcorner U leading to a system with the same number of equations and unknowns. Usually, the characteristic values of λare searched in the range 0 ≤Re(λ)≤1, as those singularity exponents λcorrespond to singular elastic solutions in the corner with unbounded stresses and strains at the corner tip, but with a finite elastic strain energy there. The characteristic values of ωare searched in the range −180o≤ωki≤180o, for each component (i=1,...,F). Recall that one ωkivalue is found per each interface or boundary face where a frictional contact is considered. Once the values of λand ω, the solution to the corner problem, are obtained, the singular stress and displacement fields can be computed using the same procedure as in Section 2.6. CHAPTER 367 Once displacements and stresses are obtained for a pair (λ,ω), it is possible to check the feasibility of the eigensolution by verifying the compression condition on the contact surfaces and the condition of frictional energy dissipation. If the computed relative sliding direction and the shear stress vector direction are parallel but with opposite directions, the frictional energy dissipation condition is fulfilled. This does not necessarily mean that eigensolutions which do not verify this condition cannot be useful, given that the friction energy dissipation condition is strictly speaking written in terms of the variation of relative displacement, not its value. Therefore, in some configurations such eigensolutions could be superposed with other solutions leading to fulfilment of the friction energy dissipation condition. Significant observations on the possibility of obtaining a complex singularity exponent λ∈C, and its implications in configurations involving frictional contact are presented in Section 3.6 and in Chapter 4 the necessary changes to the formalism and its implementation for obtaining valid solutions with λ∈Cwill be detailed. 3.4 MATLAB implementation New functions required to address corners with frictional sliding contact on outer boundaries and interfaces have been incorporated into the code developed in Section 2.7 for corners with homogeneous boundary and interface conditions. No modifications to the number of modules are needed; therefore, only the pertinent changes in the code are discussed in the following - Data input Friction coefficient can be introduced for some corner surfaces. - Definition of single-material wedges - Boundary and interface condition matrices The matrices described in Tables 3.1 and 3.2 are added for some corner surfaces. - Characteristic system assembly The main difference is the possibility that the matrix Kcorner is a rectangular matrix, which depends not only on λ, but also on as many angles ωkias faces with frictional sliding contact exist in the corner. - Solution of the characteristic system This module involves the most substantial changes. A new sub-module is developed for corners with frictional contacts. The non-linear eigensystem to be solved depends on the characteristic pair (λ,ω). In the most common and simplest case, when there is only one frictional contact, to help find the initial point for the search for the solution of (3.59), (3.60) or (3.61) a contour plot of values of σmin(λ,ω)for 0 ≤λ≤1 and −π≤ω≤πis computed, as shown in Fig. 3.4. In this plot the user can search manually for pairs of values of λand ωthat make the value of σmin(λ,ω)close to zero. These values are used as a starting point to find the solution of (3.61) with the Levenberg-Marquardt algorithm. With this algorithm, it is possible to find real and complex values of λ. The reason to use the third method (3.61) from those presented in Section 3.3.3, is that after several tests done with MATLAB, this method tends to converge to a more accurate solution than the methods based on (3.59) and (3.60). In some cases where frictional contact is imposed, the presented algorithm converges to λ∈C. In these cases, when ω∈R, typically ω=0,180◦, the results obtained with the code are valid, but when ω=0,180◦, the code cannot separate the real part of the 68 M.A. Herrera-Garrido solution from the imaginary part and ω∈Cis obtained. This solution is not valid, and the corner singularity problem needs to be solved by a different procedure presented in Chapter 4. - Displacement and stress singular fields This module requires to additionally substitute the value of ωin (3.49) or (3.57) and take into account the modified size of the boundary and interface matrices in the case of frictional contact. 3.5 Numerical examples Following the tests carried out in Chapter 2, the matrix formalism and the code developed in the present chapter are checked by comparing its results with those available in the literature. The materials employed in Section 2.8, featuring engineering constants listed in Table 2.3, are also utilized in the examples presented in this section, with minor alterations in the material ’D’. The updated engineering constants are provided in Table 3.3 with full precision Table 3.3: Engineering constants for the linear elastic materials used in the studied examples. Elastic moduli and shear moduli in GPa and fibre angles φiin degree. Material E1E2E3G12 G13 G23 ν12 ν13 ν23 φi A 68.67 0.33 B 3 0.35 C 137.89 14.48 1 5.86 1 1 0.21 0 0 0 D 137.9 14.48 14.48 5.86 5.86 5.86 0.21 0.21 0.21 φi E 137.9 14.48 14.48 4.98 4.98 4.98 0.21 0.21 0.21 φi F 5.85 0.25 G 141.3 9.58 9.58 5 5 3.5 0.3 0.3 0.32 φi In the numerical examples the boundary condition (BC) used are stress-free (F), clamped (C), symmetry (frictionless contact) (S), antisymmetry (A) and frictional contact (FC). The interface condition (IC) are perfectly bonded (B), frictionless interface (FL) and frictional contact (FC). The abbreviation “Ex” stands for the number of each example and is indicated in the corresponding column of each table when it contains more than one example or in the title of the table or figure when it contains only one example. For cases with boundary faces or interfaces with frictional contact, f≶0 is used if the authors do so in their articles and in that case µ=|f|. The symbol φicorresponds to the rotation angle around the xi−axis of the composite material fibre, see Fig. 3.3. If the composite material has the fibres lying on the x1x3-plane and parallel to x1or x3-axis the material behaves as an orthotropic material, if the fibre is lying on the x1x3or x1x2-plane rotated around the x2or x3-axis respectively an angle φi, black lines in Fig. 3.3, the composite material behaves as an apparently monoclinic material, and if the composite material has the fibre as the red line in Fig. 3.3, it behaves as an apparently anisotropic material. As remarked in the previous sections, once the λand ωsolutions in the studied problems are computed, it is necessary to check if the relative tangential displacement obtained and the friction shear stress satisfy the friction dissipation condition, having both the same direction but opposite orientation, expressions (3.11) and (3.25). Eigensolutions that meet the frictional dissipation of energy requirement are marked CHAPTER 369 Figure 3.3: Fibre position for apparently anisotropic material with (∗), and with (o)for solutions with no frictional tangential stress at the interface, including anti-plane mode. It must also be checked if the compression condition σθθ ≤0 is satisfied on the contact face. When this condition is not satisfied, the solution is marked with the symbol (±), which means that one contact face is under compression and the other face is under tension. In the examples shown in this section, the number of faces with frictional contact is F=1, except in examples 3.1.9, 3.1.14 and 3.4 where F=2. Noteworthy, problems with F>2 can be solved with the present matrix formalism and the developed computational code, but the authors have not found any previously solved problems with more than two friction contact surfaces in the literature, so a comparison is not possible. 3.5.1 Singular solutions for open corners 3.5.1.1 Singular solutions for open corners with only one single-material wedge Examples 3.1.1 - 3.1.4 For corners made of a single material, K1=E1, Table 3.4 presents the results for an isotropic wedge with different values of angles θ1and different boundary conditions, one of them being frictional sliding contact. These results are compared with those obtained using the closed form eigenequations by Sinclair [121], who observed that both positive and negative friction coefficients f≶0 should be tried to ensure that the shear stress opposes any sliding, i.e., ensure the frictional dissipation of energy condition. This can be achieved by selecting the appropriate sign of f. With the present formalism, with µ≥0, we can find the results for both signs of fcomputed with the closed-form expressions found in [121] considering −π≤ω≤π. Solutions in [121] correspond to the plane mode and coincide up to more than 9 decimal places with the results obtained with the present code. Furthermore, this present methodology allows us to find the solutions corresponding to the anti-plane mode, those with ω=±90◦, as well as some unexpected solutions, those with ω=0◦,±90◦,180◦, that involve modes in which the plane and anti-plane modes cannot be decoupled. To the best of the authors’ knowledge, these cases have not been previously reported by any other author for either isotropic or orthotropic materials, and their existence and properties will be studied in depth in a forthcoming article. The case with both boundaries of the corner having frictional sliding boundary condition with µ=1 is computed for a wedge angle of 300◦. The results obtained are the same as those for a wedge angle of 150◦with a frictional sliding boundary condition on one face and symmetry or antisymmetry boundary condition on the other face, Examples 3.1.3 and 3.1.4, respectively, in Table 3.4. The solution of Example 3.1.4 with f=−1 obtained by searching for roots of 70 M.A. Herrera-Garrido the closed-form corner-eigenequation given in [121] and also by the present formalism, with the angle ω=180◦, gives a complex singularity exponent λ∈C. This solution leads to an oscillatory stress field similar to the well-known oscillatory solution in the open model of an interface crack between dissimilar materials with the Dundurs parameter β=0 as analysed in [46, 147]. In Fig. 3.4, the contour plot of σmin values of the singular value decomposition (SVD) of the matrix Kcorner(λ,ω), matrix of the characteristic system for the singularity analysis, for the case in Example 3.1.1 is shown to illustrate how the solutions for λand ωare found. The darker blue zones correspond to values of σmin close to 0, so the singularity exponent λand the associated angle ωcan be found there. The red arrows correspond to solutions found also solving eigenequation from [121], the black arrows are the solutions corresponding to the anti-plane mode and the green arrows are the unexpected solutions mentioned above. These latter values and their existence will be studied in Chapter 5. Examples 3.1.5 - 3.1.14 For a single-material orthotropic corner, the results are compared with the results obtained by the closed-form expression of corner eigenequations deduced in [99]. These examples are shown in Table 3.5 for half-planes and in Table 3.6 for cracks with frictional-sliding contact boundary conditions. The values of ϑ0and ϑ1used for the computation of the solution are indicated, but the results were checked to be valid for any combination of ϑ1−ϑ0=180◦in Table 3.5 and ϑ1−ϑ0=360◦in Table 3.6. Examples 3.1.9 and 3.1.14 represent the cases with frictional sliding contact on the two boundary faces. In these cases, the sliding (or shear) direction angles are shown for both faces, ω0and ω1. In [99] the closed-form expressions were deduced only for the plane mode, in Tables 3.5 and 3.6 the results obtained by the present code corresponding to the anti-plane mode (ω=±90◦)and modes in which the plane and anti-plane modes cannot be decoupled (ω=0◦,±90◦,180◦) are also presented. Any combination of λ,ω0,ω1values that satisfies the conditions λ∈Z, i.e. λis an integer number, and ω1=−ω0±nπ, with n∈Z, is also a solution of Example 3.1.9. Example 3.1.8 corresponds to an orthotropic semi-plane with antisymmetric boundary condition at ϑ0and frictional sliding contact at ϑ1, the same boundary condition as in Example 3.1.4, resulting in a solution with complex singularity exponent λ∈Calso. In all examples, the compression condition at the frictional contact boundary is checked. As Examples 3.1.9 and 3.1.14 have frictional contact on both boundaries, both outer faces must simultaneously fulfil the compression condition. However, this requirement is not satisfied by the solutions marked with (±)in Tables 3.5 and 3.6. In Table 3.6 some λvalues are solution of the problem regardless of the value of ω, this is indicated as ∀ωin the table. CHAPTER 371 Figure 3.4: Contour plot of σmin values of the SVD of the matrix Kcorner in Example 3.1.1, for 0 ≤λ≤1 and −π≤ω≤π. Red arrows are solutions founds also solving eigenequation from [121], black arrows are solutions that match with the anti-plane mode and green arrows are unexpected solutions. Table 3.4: Comparison of results for 0 ≤λ≤1 and ωobtained by the present code and those found in the literature, for open single-material isotropic corners with one frictional boundary condition, Section 3.5.1.1. (*) Fulfils the frictional dissipation of energy condition. (o) No frictional tangential stress on the contact face, including anti-plane mode. Corner configuration Ex. BC1 BC2 Mat Results in literature [121] Present results λ,ω x 2 1 x ϑ0 ϑ 1 BC 1 BC 2 x 3 3.1.1 FC µ=1 ϑ0=0◦ F ϑ1=300◦ A 0.3916876291 f=1 0.3916876291 0.0◦ 0.5238206585 f=1 0.5238206585 0.0◦ - 0.3168546234 65.645◦ - 0.3168546234 - 65.645◦ - 0.6, 90.0◦(o) - 0.6, -90.0◦(o) - 0.7694649581 98.480◦(∗) - 0.7694649581 - 98.480◦(∗) - 0.209797619 145.101◦ - 0.209797619 −145.101◦ 0.1968137883 f=-1 0.1968137883 180.000◦(∗) 3.1.2 C ϑ1=300◦ 0.329046166 f=1 0.329046166 0.0◦(∗) 0.649215698 f=1 0.649215698 0.0◦(∗) 0.867301578 f=1 0.867301578 0.0◦ 72 M.A. Herrera-Garrido Corner configuration Ex. BC1 BC2 Mat Results in literature [121] Present results λ,ω x 2 1 x ϑ0 ϑ 1 BC 1 BC 2 x 3 3.1.2 FC µ=1 ϑ0=0◦ C ϑ1=300◦ A - 0.3, 90.0◦(o) - 0.3, -90.0◦(o) - 0.9, 90.0◦(o) - 0.9, -90.0◦(o) - 1 107.303◦(∗) - 1,-107.303◦(∗) 0.297128929 f=-1 0.297128929 180.0◦ 0.512105442 f=-1 0.512105442 180.0◦(∗) 3.1.3 S ϑ1=150◦ 0.484725982 f=1 0.484725982 0◦(∗) 0.760133083 0.760133083 f=1 0◦ –0.954171230 79.766◦ –0.954171230 -79.766◦ 3.1.4 A ϑ1=150◦ 0.245553923 0.245553923 f=1 0◦ – 0.285960827 44.056◦ – 0.285960827 - 44.056◦ – 0.6, 90.0◦(o) – 0.6, -90.0◦(o) 0.533683783 ±0.20979646i f=-1 0.533683783 ±0.20979646i 180.0◦ CHAPTER 373 Table 3.5: Comparison of results for 0 ≤λ≤1 and ωobtained by the present code and those found in the literature, for open single-material orthotropic corners with at least one frictional boundary condition, Section 3.5.1.1. (*) Fulfils the frictional dissipation of energy condition. (o) No frictional tangential stress on the contact face, including anti-plane mode. (±)One contact face under compression and the other contact face under tension. Corner configuration Ex. BC1 BC2 Mat Results in literature [99] λ,f Present results λ ω x x 2 x 3 1 ϑ0 180º BC1 BC2 3.1.5 F ϑ0=20◦ FC µ=1 ϑ1=200◦ C 0.423597837 f=1 0.423597837 0◦(∗) - 0.479710607 ±71.58632◦(∗) - 0.520289393 ±111.78927◦ 0.547802109 f=-1 0.547802109 ±180◦ 3.1.6 C ϑ0=20◦ 0.423597837 f=1 0.423597837 0◦ - 0.5 ±90◦(o) 0.547802109 f=-1 0.547802109 ±180◦(∗) 3.1.7 S ϑ0=20◦ 0.923597837 f=1 0.923597837 0◦ 0.047802109 f=-1 0.047802109 180◦(∗) 3.1.8 A ϑ0=20◦ 0.461798918± 0.1098658i f=1 0.461798918 ±0.1098658i 0◦ - 0.5 ±90◦(o) - 0.618471327 ±91.0944◦(∗) - 0.3815286737 ±92.0756◦ 0.707032932 f=-1 0.707032932 180◦(∗) 0.340769177 f=-1 0.340769177 180◦ 3.1.9 FC µ=1 ϑ0=20◦ 1f1=1 f2=1 1 ω0=0◦ ω1=0◦ 0.875795728 f1=-1, f2=1 0.875795728 ω0=180◦, ω1=0◦ 0.124204272 f1=1, f2=-1 0.124204272 ω0=0◦, ω1=180◦(±) 1, f1=−1, f2=−1 1, ω0=180◦, ω1=180◦ 74 M.A. Herrera-Garrido Table 3.6: Comparison of results for 0 ≤λ≤1 and ωobtained by the present code and those found in the literature, for open single-material orthotropic corners with at least one frictional boundary condition, Section 3.5.1.1. (*) Fulfils the frictional dissipation of energy condition. (o) No frictional tangential stress on the contact face, including anti-plane mode. (±)One contact face under compression and the other contact face under tension. Corner configuration Ex. BC1 BC2 Mat Results in literature [99] λ,f Present results λ,ω x x 2 x 3 1 ϑ 0 BC1 BC2 360º 3.1.10 F ϑ0=0◦ FC µ=1 ϑ1=360◦ C 0.5, f=±1 0.5, ∀ω(o) 0.222856347 f=1 0.222856347 0◦(∗) 0.722856347 f=1 0.722856347 0◦(∗) - 0.240280639 ±69.204◦(∗) - 0.740280639 ±69.204◦(∗) - 0.259719361 ±110.800◦ - 0.759719361 ±110.800◦ 0.277143653 f=-1 0.277143653 180◦ 0.777143653 f=-1 0.777143653 180◦ 3.1.11 C ϑ0=0◦ 0.222856347 f=1 0.222856347 0◦ 0.5, f=±1 0.5, ∀ω(o) 0.722856347 f=1 0.722856347 0◦ - 0.25, 90◦(o) - 0.75, 90◦(o) 0.277143653 f=-1 0.277143653 180◦(∗) 0.777143653 f=-1 0.777143653 180◦(∗) 3.1.12 S ϑ0=0◦ 0.472856347 f=1 0.472856347 0◦ 0.5, f=±1 0.5 ∀ω(o) 0.972856347 f=1 0.972856347 0◦ 0.027143653 f=-1 0.027143653 180◦(∗) 0.527143653 f=-1 0.527143653 180◦(∗) 3.1.13 A ϑ0=0◦ 0.222856347 f=1 0.222856347 0◦ 0.25, f=±1 0.25, ∀ω(o) 0.722856347 f=1 0.722856347 0◦ 0.75, f=±1 0.75, ∀ω(o) 0.277143653 f=-1 0.277143653 180◦(∗) 0.777143653 f=-1 0.777143653 180◦(∗) CHAPTER 375 Corner configuration Ex. BC1 BC2 Mat Results in literature [99] λ,f Present results λ,ω x x 2 x 3 1 ϑ 0 BC1 BC2 360º 3.1.14 FC µ=1 ϑ0=0◦ FC µ=1 ϑ1=360◦ C 0.5, f1=1, f2=1 0.5, ω0=0◦, ω1=0◦(±) 0.5, f1=-1, f2=1 0.5, ω0=180◦, ω1=0◦(o) 0.445712694 f1=-1, f2=1 0.445712694 ω0=180◦, ω1=0◦(±) 0.945712694 f1=-1, f2=1 0.945712694 ω0=180◦, ω1=0◦ 0.5, f1=1, f2=-1 0.5, ω0=0◦, ω1=180◦(o) 0.054287306 f1=1, f2=-1 0.054287306 ω0=0◦, ω1=180◦(∗) 0.554287306 f1=1, f2=-1 0.554287306 ω0=0◦, ω1=180◦(±) 0.5 f1=-1, f2=-1 0.5, ω0=180◦, ω1=180◦(±) 3.5.1.2 Singular solution for open corners with two single-material wedges Example 3.2.1 Gdoutos and Theocaris [51], employing the Williams [145] technique, deduced a corner eigenequation to study the existence of a power singularity at the apex of a plane indenter acting on an elastic half-plane, applying the Dundurs mismatch parameter β for bi-material isotropic corners for different values of the friction coefficient and the sharp edge angle. This eigenequation [51, Eq. (3)] is implemented to obtain the singularity exponent for the friction coefficients f=−1, −0.5, 0, 0.5, and 1, and the solid angles of the indenter, ϑ=30◦, 45◦, 60◦, 90◦and 120◦. The material used for the half plane is B and for the indenter is A in Table 3.3. All results obtained by solving the corner eigenequation [51, Eq. (3)] are compared with the results obtained for each case by the present general matrix formalism with perfect match at least up to 12 digits. Later, Comninou [30] studied the same case, but employing the Mellin transform that essentially leads to the same expression of eigenequation as in [51] and therefore to the same results. In Table 3.7 some of the results obtained are shown. The results shown for a friction coefficient µ=0 have been obtained in two different ways with the present formalism for verification purposes. First, using the interface condition matrices considering specifically the frictionless case as presented in Chapter 2, and then with the interface condition matrices in Table 3.2 considering µ=0 obtaining identical results up to 12 digits. Examples 3.2.2 and 3.2.3 Churchman and Hills [29] present a simpler corner eigenequation for the same case as in Example 3.2.1 but in isotropic corners with the same materials and ϑ2=90◦angle of indenter. In this kind of problem, as both wedges are made of the same material, the results are independent of the properties of the material. In Fig. 3.5, the results λand ωobtained with the present formalism are superimposed on the Fig. 3 in [29], and the numerical results are presented in Table 3.8. The positive values of −fin Fig. 3.5 correspond 82 M.A. Herrera-Garrido 3.5.2 Singular solutions for closed (periodic) corners 3.5.2.1 Singular solutions for closed (periodic) corners with one multi-material wedge Example 3.3.1. Frictional interface crack in an isotropic bi-material Interface cracks with frictionless and frictional contact in isotropic bi-materials were first studied by Comninou [31, 32]. The present results in Table 3.11 for the frictional interface crack are compared with the values of the singularity exponent λgiven by the closed-form expression deduced in [32]. The configuration used is shown in Fig. 3.8, where Mat1and Mat2are material B and A in Table 3.3, respectively, and IC1 in ϑ=−180◦is frictional contact and IC2 in ϑ=0◦ is a perfectly bonded interface condition. The friction coefficients for IC1 are indicated in Table 3.11 for each case in the first column2. None of the solutions obtained with f>0 with the eigenequation found in [32] satisfies the frictional dissipation of energy condition. This is because the compliant material is the upper material in the interface crack corner. Therefore, to obtain results that fulfil the frictional dissipation of energy condition, one possibility is to interchange the material position obtaining the same result as if we change the sign of f, which implies changing the direction of the friction shear stress considering compressions in the contact zone. The present matrix formalism does not need to use negative friction coefficients, as it has the possibility of having a direction of the tangential stress vector angle in the opposite direction. It has also been checked that for a friction coefficient µ=0 or in the homogeneous case (Mat1=Mat2), the combination of materials does not matter, λ=0.5, as explained in [31, 32]. The singular solutions associated with λ1,ω1and λ3,ω3represent a symmetric set of nonsymmetric solutions, each one can be obtained from the other simply by multiplying it by −1. Figure 3.8: Configuration used for interface crack examples. Example 3.3.2. Frictional interface crack in an anisotropic bi-material The results by Comninou [31, 32] for the contact model of interface cracks in isotropic bi-materials were later generalised to anisotropic bi-materials by many authors, see Section 1.2.2.2 for most relevant references. In Table 3.12, the present results are compared with those obtained by solving a closed-form eigenequation for λderived by Sung and Chung [125]. The configuration used is as shown in Fig. 3.8, with Mat1and Mat2indicated in Table 3.12 and the interface conditions IC1 and IC2 are frictional contact and perfectly bonded, respectively. Orthotropic 2In the expression proposed by Comninou [32], it should be taken into account the following relation between λComninou and λused in the present article λ=1−λComninou (3.63) CHAPTER 383 Table 3.11: Example 3.3.1. Numerical values for the singularity exponent 0 <λ<1 and direction of the tangential stress vector angle ωfor an interface crack in an isotropic corner, see Section 3.5.2.1. (*) Fulfils the frictional dissipation of energy condition. (o) No frictional tangential stress in the interface, including anti-plane mode. µResults in [32] f=µ Results in [32] f=−µλ1,ω1λ2,ω2λ3,ω3λ4,ω4 0 0.5 0.5 0.5 0.5 0.5 0.5 -90◦(o) 0◦90◦(o) 180◦ 0.2 0.486626904 0.513373096 0.5 0.513373096 0.5 0.486626904 -90◦(o) 0◦(∗)90◦(o) 180◦ 0.4 0.473300850 0.52669915 0.5 0.52669915 0.5 0.473300850 -90◦(o) 0◦(∗)90◦(o) 180◦ 0.5 0.466669922 0.533330078 0.5 0.533330078 0.5 0.466669922 -90◦(o) 0◦(∗)90◦(o) 180◦ 0.6 0.460067898 0.539932102 0.5 0.539932102 0.5 0.460067898 -90◦(o) 0◦(∗)90◦(o) 180◦ 0.8 0.446972184 0.553027816 0.5 0.553027816 0.5 0.446972184 -90◦(o) 0◦(∗)90◦(o) 180◦ 1 0.43405509 0.56594491 0.5 0.56594491 0.5 0.43405509 -90◦(o) 0◦(∗)90◦(o) 180◦ materials are represented with fibres that match the x1-axis, apparently monoclinic materials are orthotropic materials with the fibre in x1x3-plane but with a φ2angle measured from the x1-axis or with the fibre in x1x2-plane but with a φ3angle measured from the x1-axis and apparently anisotropic materials are orthotropic materials whose fibre is rotated an angle φ2 around the x2-axis and an angle φ3around the x3-axis, as represented in Fig. 3.3. Excellent agreement is found between [125] and the present results for orthotropic and apparently monoclinic materials with fibres in the x1x2-plane, but no coincidence is found either in the case of monoclinic materials with fibres in the x1x3-plane or in general anisotropic materials. This mismatch can be explained by analysing the frictional sliding contact model used by Sung and Chung in [125], where the angle ω, direction of the tangential stress, is not defined, assuming that the values of the frictional tangential stresses σi2with i=1,3 are either zero or ±µσ22. A procedure similar to that previously developed for open corners with frictional contact in Chapter 7 is applied in Chapter 6 to the frictional interface crack problem, whose results agree with those obtained by the present matrix formalism. As was checked for the isotropic case, it was also checked for any orthotropic or monoclinic material with a PES parallel to the coordinate planes x1x3or x1x2, in the frictionless case with µ=0 and in the case of similar or identical materials (Mat1=Mat2with φ1 i=φ2 i), that λ=0.5 as well; however, in the case of dissimilar general anisotropic materials with µ=0, λ=0.5±iε. Apart from the special cases (frictionless contact, similar or identical materials, antiplane case with zero compression and zero friction shear), if we only consider the singular solutions that meet the frictional dissipation of energy criterion, then Re(λ)>0.5. This observation is important as it shows that under friction and with dissimilar materials the singularity will be weaker than in the frictionless case or the case of similar or identical materials. When the last example in Table 3.12, corresponding to a bi-material formed by a semiplane of material G with ϕ2=60◦and a semi-plane of material E with ϕ2=30◦and ϕ3=−2◦ (an apparently anisotropic material), is solved with the present code, complex values are obtained for both λand ω. Therefore, this problem is also analysed with the code currently 84 M.A. Herrera-Garrido under development for the special case of complex λwhere two angles ω,ω′and ω′′, are obtained for the frictional shear vector. The treatment of complex singularity exponents λis discussed in more detail in Section 3.6. For completeness, the result obtained with this code are added to this table as well. The crack tip singularity in interface cracks with frictional contact was also studied by Leguillon [84] for isotropic and anisotropic bi-materials. In the model proposed in [84], the sliding direction is assumed as a kinematic condition and can be chosen in any direction between the plane mode and the anti-plane mode, this differs from the model considered in the present work. Thus, good agreement can only be achieved when comparing the results in [84] to those of the present study for isotropic materials and orthotropic materials with a stacking sequence of [0,90], when prescribing a direction of the tangential stress vector angle of 0◦or 90◦. Table 3.12: Example 3.3.2. Numerical values for the singularity exponent 0 <λ<1 and direction of the tangential stress vector angle ωfor an interface crack between anisotropic materials, see Section 3.5.2.1. (*) Fulfils the frictional dissipation of energy condition. (o) No frictional tangential stress in the interface, including anti-plane mode. µMat1Mat2 Results in [125] λ1,ω1λ2,ω2λ3,ω3λ4,ω4 0.5 G φi=0◦ E φi=0◦ 0.501651679 0.5 (o) 0.501651679 0.5 (o) 0.498348321 f=0.5 -90◦0◦(∗)90◦180◦ 10.503303268 0.5 (o) 0.503303268 0.5 (o) 0.496696732 f=1 -90◦0◦(∗)90◦180◦ 0.5 G φ3=60◦ E φ3=30◦ 0.502122213 0.5 (o) 0.502122213 0.5 (o) 0.496523204 f=0.5 -90◦-0◦(∗)90◦180◦ 10.503552288 0.5 (o) 0.503552288 0.5 (o) 0.489789778 f=1 -90◦-0◦(∗)90◦180◦ 0.5 G φ2=60◦ E φ2=30◦ 0.509595142 0.5 (o) 0.511519758 0.5 (o) 0.488480242 f=0.5 237.377◦−45.007◦(∗)57.377◦134.992◦ 10.519172879 0.5 (o) 0.523009419 0.5 (o) 0.476990581 f=1 237.377◦−45.007◦(∗)57.377◦134.992◦ 1G φ2=60 Eφ2=30◦ φ3=−2◦ 0.516405370 f=1 0.503192599 240.639◦ 0.523049834 −47.067◦(∗) 0.477013012 137.057◦ µMat1Mat2λ5,ω5λ6,ω6 0.5 G φ2=60◦ E φ2=30◦ 0.493212664 0.506787336 -135.007◦44.992◦(∗) 10.486431495 0.513568505 -135.007◦44.992◦(∗) 1G φ2=60◦ Eφ2=30◦ φ3=−2◦ 0.483420518 0.506309057±i0.0018419970 -137.876◦ω′=54.256◦,ω′′ =50.543◦ Example 3.3.3. Closed crack tip terminating at a perfectly bonded interface Fig. 3.9 represents the configuration considered in the examples presented in Table 3.13, where the crack is located in one of the bonded materials and not on the interface. This case can be studied as a closed corner consisting of a three-material wedge, and was studied by Comninou and Dundurs [33] among many other authors. Some tests have been carried out to compare the results of the present procedure with the results obtained from eigenequation [33, Eq. (5)], all with perfect agreement. In Table 3.13 some of these numerical results are shown. Mat1 and Mat3are material B and Mat2material A in Table 3.3. In Table 3.13, the angle θ1−θ0 CHAPTER 385 is indicated for each case. In this example, the inclination of the crack that terminates at the interface is varied by changing the angle 0 ≤θ0≤180◦, while keeping the angles θ1=0◦ and θ2=180◦fixed. The θ3=θ0+360◦is varying together with θ0. The friction coefficient used at the interface IC1 is µ=1. IC2 and IC3 are perfectly bonded interfaces. In Table 3.13, crack inclinations for which no singularity exponent was found in the range 0 <λ<1 have been marked with the symbol −. x x x 1 3 2 Figure 3.9: Corner configuration for examples in Table 3.13, crack reaching a perfectly bonded interface. It has also been checked that: • When θ0=90◦in Fig. 3.9, the results do not depend on the friction coefficient, this is evident due to the symmetry generated in the corner geometry. • When θ0approaches the values 0◦or 180◦, in Fig. 3.9, the results for λtend to the λ value in the case of a frictional interface crack. • The singular solutions associated with λ1,ω1and λ3,ω3represent a symmetric set of non-symmetric solutions, each one can be obtained from the other just by multiplying it by −1 as they correspond to the anti-plane mode. These solutions are also symmetric with respect to the solutions associated with λ1,ω1and λ3,ω3for a problem with a crack inclination θ′ 1−θ′ 0= (θ1−θ0)−θ2, which could be easily explained by symmetry arguments in view of the isotropy of the materials. 3.5.2.2 Singular solutions for closed corners with two single-material wedges Example 3.4 Arias et al. [1] studied the case of two isotropic materials sliding one over the other with two frictional contact interfaces for the case of tectonic plates. Table 3.14 shows no considerable difference between the results obtained with the present matrix formalism and the results obtained with the eigenequation presented by Arias et al. [1] for isotropic homogeneous corners. The material used in the calculation is material A from Table 3.3, but the same results are obtained for any isotropic material, as long as the same material is chosen for both wedges. As this example has frictional contact on IC1 and IC2, both interfaces must simultaneously fulfil the compression condition. This requirement is checked and in most cases presented in Table 3.14 it is not met. Solutions not fulfilling the compression condition simultaneously at both interfaces are marked with (±). 86 M.A. Herrera-Garrido Table 3.13: Example 3.3.3. Numerical values for the singularity exponent 0 <λ<1 and direction of the tangential stress vector angle ωfor a crack terminating at a perfectly bonded interface in an isotropic corner with a friction coefficient in the crack µ=1, see Section 3.5.2.1. (*) Fulfils the frictional dissipation of energy condition. (o) No frictional tangential stress in the interface, including anti-plane mode. θ1−θ0 Results in [33] λ1,ω1λ2,ω2λ3,ω3λ4,ω4 0◦0.43405509 0.5 0.43405509 0.5 0.56594491 -90◦(o) 0◦90◦(o) 180◦(∗) 15◦0.508410685 0.543314279 0.508410685 0.543314279 0.639143816 -90◦(o) 0◦90◦(o) 180◦(∗) 30◦0.633740388 0.594708947 0.633740388 0.594708947 0.683933616 -90◦(o) 0◦90◦(o) 180◦(∗) 45◦-0.656319424 - 0.656319424 0.669427348 -90◦(o) - 90◦(o) 180◦ 60◦-0.730292305 - 0.730292305 0.6347699 -90◦(o) - 90◦(o) 180◦ 75◦0.822212612 0.814682253 0.822212612 0.814682253 0.62589152 -90◦(o) 0◦(∗)90◦(o) 180◦ 90◦0.6633562 0.869773695 0.663356243 0.869773695 0.663356243 -90◦(o) 0◦(o) 90◦(o) 180◦(o) 105◦0.62589152 0.814682253 0.62589152 0.814682253 0.822212612 -90◦(o) 0◦90◦(o) 180◦(∗) 120◦0.6347699 0.730292305 0.6347699 0.730292305 - -90◦(o) 0◦90◦(o) 135◦0.669427348 0.656319424 0.669427348 0.656319424 - -90◦(o) 0◦90◦(o) 150◦0.683933616 0.594708947 0.683933616 0.594708947 0.633740388 -90◦(o) 0◦(∗)90◦(o) 180◦ 165◦0.639143816 0.543314279 0.639143816 0.543314279 0.508410685 -90◦(o) 0◦(∗)90◦(o) 180◦ 180◦0.56594491 0.5 0.56594491 0.5 0.43405509 -90◦(o) 0◦(∗)90◦(o) 180◦ CHAPTER 387 Table 3.14: Example 3.4. Numerical values for the singularity exponent 0 <λ<1 and direction of the tangential stress vector angle ωfor two isotropic single-material wedges sliding with friction in both interfaces, see Section 3.5.2.2. In [1], a friction coefficient on one face is denoted as f, and a specific ratio of the friction coefficients on both faces as w. (*) Fulfils the frictional dissipation of energy condition. (o) No frictional tangential stress in the interface, including anti-plane mode. (±)One contact interface under compression and the other contact interface under tension. Corner configuration ICI IC2 Results λin [1] Present results λ,ω1,ω2 µ=0.4 ϑ0=30◦ µ=0.4 ϑ1=330◦ 0.349057987 0.349057987 f=0.4, w=1ω0=ω1=0◦ 0.726220223 0.726220223 (±) f=0.4, w=1ω0=ω1=0◦ 0.734403364 0.734403364 (±) f=-0.4, w=- 1 ω0=180◦,ω1= 0◦ 0.147909364 0.147909364 f=-0.4, w=- 1 ω0=180◦,ω1= 0◦ 0.410676802 0.410676802 (∗) f=0.4, w=-1 ω0=0◦,ω1= 180◦ 0.704192317 0.704192317 (±) f=0.4, w=-1 ω0=0◦,ω1= 180◦ 0.349057987 0.349057987 f=-0.4, w=1ω0=ω1=180◦ 0.726220223 0.726220223 (±) f=-0.4, w=1ω0=ω1=180◦ -0.6 (o) ω0=ω1=±90◦ µ=0.5 ϑ0=30◦ µ=1 ϑ1=330◦ 0.351115240 0.351115240 f=-1, w=0.5 ω0=ω1=0◦ 0.731708957 0.731708955 (±) f=-1, w=0.5 ω0=ω1=0◦ 0.737966428 0.737966428 (±) f=-1, w=- 0.5 ω0=180◦,ω1= 0◦ 0.513337305 0.513337355 (±) f=1, w=-0.5 ω0=0◦,ω1= 180◦ 0.600884868 0.600884868 (±) f=1, w=-0.5 ω0=0◦,ω1= 180◦ 0.890144750 0.890144750 (±) f=1, w=-0.5 ω0=0◦,ω1= 180◦ 0.511522494 0.511522501 (±) f=1, w=0.5 ω0=ω1=π 0.712545420 0.712545420 (±) f=1, w=0.5 ω0=ω1=π -0.6 (o) ω0=ω1=±90◦ 88 M.A. Herrera-Garrido 3.6 Concluding remarks This chapter builds upon the review and improvement of the matrix formalism initiated in the preceding chapter, specifically targeting scenarios involving any finite number of frictional sliding contact surfaces. Leveraging the modular design introduced in the code developed in the previous chapter, additional modules necessary for addressing frictional contact are incorporated. This chapter represents a seamless extension of the preceding work, ensuring comprehensive coverage of frictional contact scenarios within the developed framework. With these steps, Objectives 1.1 and 1.2 are completed, as the matrix formalism presented in [96, 97] has been fully reviewed and implemented in the computational code. After performing this extension of the code, both the matrix formalism and the developed code were once again verified by comparing the results obtained with those obtained using methodologies developed by other authors found in the literature. In the examples presented in this chapter, all containing at least one contact surface with friction, more differences have been found than in the previous chapter, all justified. In the cases solved for isotropic materials, in addition to the solutions corresponding to the anti-plane mode not obtained using the methods proposed by other authors, unexpected elastic solutions have been found, in which the plane and anti-plane modes cannot be decoupled despite the elastic symmetry of the corner, Examples 3.1.1-3.1.4. In the case of orthotropic materials that maintain elastic symmetries with the coordinate planes, this type of unexpected elastic solution has also been found, Examples 3.1.5, 3.1.8 and 3.1.10. Apart from these comments, all results that can be compared show complete coincidence. On the other hand, once the corner loses elastic symmetry with respect to the coordinate planes, that is, once anisotropic materials are involved in the problem, large discrepancies start to emerge with the results obtained by the methodologies proposed by other authors. In the example 3.2.5 where comparisons were made with the numeric results presented by Poonsawat [111] and by Magnier [94], the differences found were justified by the use of a different model of the Coulomb law for friction, not as general as the one proposed in these chapter since it does not allow relative sliding in any direction. This is supported by the cases in which the corner does present certain symmetries in which the angle of the tangential stress direction ω=0◦or 90◦, the results coincide entirely. The same happens with the example of the interface crack with friction in anisotropic materials, Example 3.3.2, since Sung and Chung [125] do not consider the angle of tangential stress as an unknown variable neither. On the other hand, as it was previously announced in the conclusions of the previous chapter for the frictionless case, in example 3.2.4, despite Chen et al. [28] used a similar friction model, the results obtained using the code presented in this chapter and those shown in the graphs of the cited article do not match. Therefore, in addition to fulfilling Objective 3 of this thesis, methodologies capable of verifying these results will be developed in Chapters 6 and 7. In addition to the differences mentioned in the previous paragraph, the case of complex solutions needs to be highlighted. As shown in the numerical results, in cases of frictional contact between isotropic materials, or orthotropic materials with a PES x3=0, potentially valid solutions with λ∈Cand ω∈Rmay appear. In the case of frictionless contact on one face of an anisotropic material and complex λ∈C, it was observed by studying the solution in displacements that there are two different relative directions of the tangential stress vector given by angles ω′ uand ω′′ uon the surface in contact. Assuming a continuous variation of the friction coefficient mu from zero to positive values, it is considered that the CHAPTER 389 values of ω′ u,ω′′ uand λalso vary continuously. Therefore, another procedure is needed for modelling the frictional contact, which considers a complex λand two values of ω′ uand ω′′ uas unknowns on each frictional surface. These cases of complex λare not covered in the present chapter for the sake of brevity. However, by taking advantage of one of the benefits of code, which is its ability to be expanded, the code can be generalised to allow solutions of problems with frictional contact and complex singularity exponent λ∈C. This new procedure, which considers two angles, ω′and ω′′, of frictional shear vector for a λ∈C is presented in Chapter 4. With all this, Objective 1 is mostly achieved, as a review and verification of the general matrix formalism has been conducted, pending only the verification of results that could not be confirmed, which will be addressed in Chapters 6 and 7, along with a new semianalytic procedure to obtain valid results involving friction and complex singularity exponent λ∈C, which will be addressed in the next chapter of this thesis. Additionally, as unexpected solutions have been obtained, progress has been made towards achieving Objective 2 of this thesis. CHAPTER 4 Generalization of the matrix formalism for complex solutions in presence of frictional contact 4.1 Introduction The numerical results highlighted in Chapter 3 illustrate that in scenarios involving frictional contact between isotropic materials, or orthotropic materials with a Plane of Elastic Symmetry (PES) x3=0, potentially valid solutions may arise where λ∈Cand ω∈R. When frictionless contact occurs on one face of an anisotropic material with a complex λ∈C, an analysis of displacement solutions, by the code developed in Chapters 2 and 3, reveals the presence of two distinct relative sliding angles, ω′ uand ω′′ u, on the contacting surface. Assuming a smooth transition of the friction coefficient from zero to positive values, it is inferred that the parameters ω′ u,ω′′ u, and λexhibit continuous variation. Hence, an alternative methodology is required to model frictional contact, taking into account a complex λand considering two values of ω′ uand ω′′ uas variables on each frictional surface. These cases of complex λwere not covered in the previous chapter for the sake of brevity and simplicity of explanations. However, by taking advantage of one of the benefits of the Matlab computational code, which is its ability to be expanded, the code is generalised to allow solutions of problems with frictional contact and complex singularity exponent λ∈C. This new procedure, which considers two angles, ω′and ω′′, of frictional shear vector for a λ∈C has already been implemented and tested, and it is presented in this chapter, although it was already used in Example 3.3.2 in Table 3.12. The existence of λ∈Csolution entails oscillatory stress fields, which may involve an infinite number of compression and tension zones that interchange with each other in the contact zone near the corner apex, thus not meeting the compression condition. To give a physical meaning and utility to a complex singular eigensolution with λ∈Ceither for frictional or frictionless contact, we propose the concept of small scale tension zone, where the distance from the corner apex to the far end of the tension zone that is located farthest 91 98 M.A. Herrera-Garrido (a) (b) (c) (d) Figure 4.3: Angular variations of the singular displacement and stress fields at r=1 for a frictionless delamination crack in a laminate with lower material G in Table 3.3 with φ2=60◦and φ3=0◦, and material E in Table 3.3 with φ2=30◦and φ3=−2. Figs. 4.3a and 4.3b depict the real and imaginary parts of the displacement field, respectively. Similarly, Figs. 4.3c and 4.3d illustrate the real and imaginary components of the stress field, respectively. The singularity exponent λ=0.5+i0.000851280. (a) (b) Figure 4.4: In Fig. 4.4a, a contour plot illustrates the minimum singular values σmin resulting from the Singular Value Decomposition (SVD) of the matrix Kcorner(λ,ω). These values vary with Re(λ)and ωfor an anisotropic bi-material with a frictional interface crack, using elastic constants of materials G and E (see Table 3.3). The friction coefficient is set to µ=0. Additionally, in Fig. 4.4b, a plot shows the minimum singular values σmin of Kcorner with λ=0.5+i0.000851280 as a function of ωfor the same problem setup. These graphs correspond to the representation in the previous version of the code. CHAPTER 499 (a) (b) (c) (d) -4 -2 0 2 4 -4 -2 0 2 4 -0.5 0 0.5 1 1.5 Figure 4.5: 3D contour plot of the minimum singular values σmin of the SVD of the matrix Kexpanded corner as a function of Re(λ),ω′and ω′′ with Im(λ) = 0 in Fig. 4.5a and Im(λ) = 0.000851280 for in Figs. 4.5b to 4.5d. The solved corner is an anisotropic bi-material interface crack with friction contact. Elastic constants for the upper material and lower material are of material G with φ2=60◦and φ3=0◦, and material E with φ2=30◦and φ3=−2◦in Table 3.3, respectively. Friction coefficient µ=0. This graph corresponds to the representation in the new version of the code. ω′=ωRe(u) and ω′′ =ωIm(u) in Fig. 4.3b and the stress and displacement fields obtained are the same as shown in Fig. 4.3. 4.4.2 Frictional interface crack in an anisotropic bi-material After increasing the friction coefficient between the surfaces of the interface crack to µ=1, the outcomes for the same problem are now computed with both code versions, and the results are presented in Table 4.1. Focusing on the solutions obtained with the current version of the code presented in Table 4.1 with Im(λ) = 0, that is, the second and fourth solutions, the outcomes closely align, coinciding beyond the 10th decimal place with the results obtained with the previous version of the code. In these cases, ω′or ω′′ assumes the value of the tangential stress angle ωobtained from the previous code version, while the other one remains unconstrained. If from all the possible solution, it is chosen ω′′ =ω′, the solution exhibits multiplicity of 2. 100 M.A. Herrera-Garrido Figure 4.6: 2D contour plot of the minimum singular values σmin of the SVD of the matrix Kexpanded corner (0.5+ i0.000851280,ω′,ω′′)as a function of ω′and ω′′. The solved corner is an anisotropic bi-material with friction contact interface crack. Elastic constants for the upper material and lower material are of material G with φ2=60◦ and φ3=0◦, and material E with φ2=30◦and φ3=−2◦in Table 3.3, respectively. Friction coefficient µ=0. For the results with low Im(λ)in Table 4.1 obtained with the current code version, that is, the first and third solutions shown, Re(λ)coincide up to 8 decimal places with the value of the singularity exponent λ∈Robtained with the previous version of the code. The tangential stress angles ω′obtained with the current version of the code associate with those λvalues, match up to 6 decimal places with the tangential stress angle ωobtained with the previous version of the code. However, the solution with significant Im(λ), solution number six in Table 4.1, the outcome obtained from each code version diverge substantially. The accurate result aligns with the current code version. This is evident when comparing the displacement and stress results in the crack faces of the interface crack, as shown in Tables 4.2 and 4.3. In particular, the normal and tangential stress obtained with the previous code version do not satisfy the Coulomb law outlined in (3.9) |τ(r=1,ϑw)|=qσrθ(r=1,ϑ0)2+σ3θ(r=1,ϑ0)2 =q0.023328082+(−0.02708534)2=0.03574654 =−µσθθ (4.12) for the real part and |τ(r=1,ϑ0)|=qσrθ(r=1,ϑ0)2+σ3θ(r=1,ϑ0)2 =q0.015547452+(−0.02155939)2=0.02658064 =−µσθθ (4.13) for the imaginary part. Moreover, the tangential stress vector does not align with the tangential displacement vector, as specified by the condition in (3.25) σrθ(r=1,ϑ0) = 0.02332808 =−κu2 r(r=1,ϑ0)−u1 r(r=1,ϑ0)=−0.50624298κ σ3θ(r=1,ϑ0) = −0.02708534 =−κu2 3(r=1,ϑ0)−u1 3(r=1,ϑ0)=−0.54961245κ. (4.14) CHAPTER 4101 Table 4.1: Singularity exponent λand tangential stress angle ωsolutions obtained with both code versions for an anisotropic bi-material with friction contact interface crack. Elastic constants for the upper material and lower material are of material G with φ2=60◦and φ3=0◦, and material E with φ2=30◦and φ3=−2◦in Table 3.3, respectively. Friction coefficient µ=1. Previous code version (Chapter 3) λ0.483420518 0.503192599 0.523049834 0.477013012 0.506885327 ±i0.002658000 ω-137.88◦-119.36◦-47.07◦137.06◦51.02◦±i2.37◦ Current code version (Chapter 4) λ0.483420514 ±i0.000004831 0.503192599 0.523049832 ±i0.000008869 0.477013012 0.505290314 ±i0.002086948 ω′-137.88◦-119.36◦-47.07◦137.06◦51.80◦ ω′′ -172.33◦ω∈R-225.18◦ω∈R57.82◦ Table 4.2: Real part of the stress and displacement values obtained with both code versions at the friction contact interface of a bi-material. Elastic constants for the upper material and lower material are of material G with φ2=60◦ and φ3=0◦, and material E with φ2=30◦and φ3=−2◦in Table 3.3, respectively. Friction coefficient µ=1. Material uθuru3σθθ σrθσ3θ Previous version Mat1-0.00219808 0.31427963 -0.25153692 -0.03569923 0.02332808 -0.02708534 Mat2-0.19196335 0.29807553 Current version Mat1-0.00058743 0.48926522 -0.46465591 -0.03471296 0.02146826 -0.02727826 Mat2-0.3197165 0.56326195 as there is no κvalue that fulfill simultaneously both expressions. The same conclusion is obtained with the imaginary part of the stress and displacement results. In contrast, the real and imaginary part of the stress and displacement outcomes obtained with the current version of the code comply with (3.9) |τ(r=1,ϑw)|=qσrθ(r=1,ϑ0)2+σ3θ(r=1,ϑ0)2 =q0.021468262+(−0.02727826)2=0.03471296 =−µσθθ , (4.15) |τ(r=1,ϑ0)|=qσrθ(r=1,ϑ0)2+σ3θ(r=1,ϑ0)2 =q0.000490102+(−0.00077901))2=0.00092036 =−µσθθ , (4.16) and with (3.25) σrθ(r=1,ϑ0) = 0.02146826 =−κu2 r(r=1,ϑ0)−u1 r(r=1,ϑ0)=0.808981723κ σ3θ(r=1,ϑ0) = −0.02727826 =−κu2 3(r=1,ϑ0)−u1 3(r=1,ϑ0)=−1.027917861κ. (4.17) with κ=0.02653739 fulfilling both the collinearity between the relative tangential displacements vector and the tangential stress vector and the frictional dissipation of energy condition. Now that the existence of complex solutions in the presence of frictional contact has been established, it is possible to analyze the size of the proposed small-scale tension zone. This zone refers to the region where the distance from the corner apex to the far end of the tension zone, located farthest from the corner vertex, is much smaller than the relevant characteristic length of the boundary value problem under study. This phenomenon arises due to the oscillatory behavior of the normal stress field. In Fig. 4.7 it is represented the normal stress field along the interface crack with λ=0.505290314 +i0.002086948. In this case, it can be observed that with Im(λ)being so small, compression is maintained across 102 M.A. Herrera-Garrido Table 4.3: Imaginary part of the stress and displacement values obtained with both code versions at the friction contact interface of a bi-material. Elastic constants for the upper material and lower material are of material G with φ2=60◦and φ3=0◦, and material E with φ2=30◦and φ3=−2◦in Table 3.3, respectively. Friction coefficient µ=1. Material uθuru3σθθ σrθσ3θ Previous version Mat100.42363203 -0.41456792 -0.02651698 0.01554745 -0.021559391 Mat2-0.28032459 0.50432975 Current version Mat1-0.00126697 -0.14066748 0.17254811 -0.00092036 0.00049010 -0.00077901 Mat20.10304949 -0.21484174 Figure 4.7: Re(σθθ ), for a λ∈C, as a function of the crack tip distance rconsidering either λor Re(λ) in (2.98), for an anisotropic bi-material with a frictional interface crack. Elastic constants for the upper material and lower material are of material G with φ2=60◦and φ3=0◦, and material E with φ2=30◦and φ3=−2 in Table 3.3, respectively. Friction coefficient µ=1. Singularity exponent λ=0.505290314 +i0.002086948. the entire interval of rshown in Fig. 4.7, 10−10 ≤r≤10. Thus, the solution maintains the prescribed contact conditions at the corner. In Fig. 4.8, just for illustrate the effect of the imaginary part of the singularity exponent on the oscillatory behavior of the normal stress, the same solution in Fig. 4.7 is shown but with Im(λ) = 0.3. In Fig. 4.8, we can observe that while we have compression for r≥10−10, we have tension for r≤10−10. This zone is what we can define as the small-scale tension zone. There may be different factors that increase the complex part of the singularity exponent, in the interface crack, for example, when the material asymmetry at the corner increases, indicating greater anisotropy, solutions with singular exponents feature a larger imaginary part. In other cases, increasing the frictional coefficient or modifying the geometry may increase also the imaginary part of the singularity exponent. The higher the imaginary part of the singularity exponent, the more pronounced the oscillations in the singular stress fields, as evident from the following expressions: CHAPTER 4103 Figure 4.8: Illustrative case, Re(σθθ )for a λ∈C, as a function of the crack tip distance rconsidering either λor Re(λ) in (2.98), for an anisotropic bi-material with a frictional interface crack. Elastic constants for the upper material and lower material are of material G with φ2=60◦and φ3=0◦, and material E with φ2=30◦and φ3=−2 in Table 3.3, respectively. Friction coefficient µ=1. Singularity exponent assumed for illustrative purposes λ=0.505290314 +i0.3. ϕ,r=λrλ−1wm(r,θ) = λr(Re(λ)−1+iIm(λ))wm(r,θ) = λr(Re(λ)−1)r(iIm(λ))wm(r,θ) = λr(Re(λ)−1)wm(r,θ)(cos(Im(λ)ln(r))+isin(Im(λ)ln(r))) = λr(Re(λ)−1)Re[cos(Im(λ)ln(r))Re(wm(r,θ))−sin(Im(λ)ln(r))Im(wm(r,θ))] (4.18) 4.5 Concluding remarks The matrix formalism for stress singularity problems and its implementation in Matlab presented in Chapters 2 and 3 is capable of obtaining valid solutions with λ∈Cas long as there are no frictional contacts in the multi-material corner, or in cases where there is frictional contact and the in-plane and anti-plane modes are decoupled, that is, for values of ω=0◦,±90◦, or ±180◦. In cases involving friction but where ω=0◦,±90◦or ±180◦, the formalism fails, returning solutions that do not satisfy the necessary collinearity between relative displacement and the tangential stress. In this chapter, that deficiency in the general matrix formalism is remedied, as a methodology has been developed that, by separating the real and imaginary parts of the characteristic system, allows for the definition of two tangential stress directions, ω′and ω′′, for each interface or material face in contact with friction, corresponding to the real and imaginary parts of the solution, respectively. Subsequently, thanks to the modularity of the Matlab code, it was straightforward to implement these changes in the implementation. The code has been validated with case studies of cracks at interfaces between anisotropic materials. It has been verified that the code is able to reproduce the results obtained with 104 M.A. Herrera-Garrido previous versions of the code for cases of contact without friction or with zero friction. In addition, it has been shown that the code is able to find valid solutions for cases of contact with positive friction, which were not possible with the previous version of the matrix formalism and its computational code implementation. It has been checked that the solutions satisfy the Coulomb contact conditions and the collinearity between the tangential displacement vector and the tangential stress vector. The concept of small-scale tension zone has been proposed to give physical meaning and practical utility to the complex solutions. It has been observed that the existence of solutions with λ∈Cimplies oscillations in the stress field, which generate zones of tension and compression alternating in the contact zone near the corner vertex. It has been suggested that these tension zones can be negligible if the distance from the corner vertex to the end of the farthest tension zone is much smaller than the characteristic length of the problem. The effect of the imaginary part of λon the size of the small-scale stress zone has been analyzed. The procedure developed in this chapter concludes the review and enhancement of the matrix formalism presented in [96, 97], as well as its implementation into a computational code. The results obtained with this code could not be verified through comparison with results from other authors, as this decomposition into real and imaginary parts allowing for different angles of the tangential direction is entirely novel. Therefore, with the work presented in this chapter, the Objective 1 of this thesis is considered completed. CHAPTER 5 Unexpected singular elastic solutions coupling in-plane and out-of-plane displacements in presence of frictional contact 5.1 Introduction Traditionally, in the context of isotropic materials subjected to Generalized Plane Strain (GPS), where the solution in displacement (and stresses) remains independent of the x3- coordinate, denoted as ui=ui(x1,x2) (i=1,2,3), it is considered that the elastic solution can be decomposed into in-plane and anti-plane strain problems [4, 116]. Moreover, this decomposition applies not only to isotropic materials, but generally, if the plane x3=0 is plane of elastic symmetry (PES) (e.g. in fibre-reinforce materials with fibres parallel or perpendicular to the plane x3=0), any solution can be usually decomposed into in-plane and an anti-plane (out-of-plane) part [134, Section 3.1] and [14, 65]. Such decomposition allows for the independent study of two simpler problems, an in-plane and an anti-plane problem, rather than addressing a more complex problem. Nevertheless, as will be demonstrated in the present work, there are configurations of corner problems whose singular eigensolutions cannot be generally decomposed in this manner when frictional contact is included by considering the Coulomb law of friction. In [96, 97], it was observed for the first time that in the studied case of a closed corner formed by an isotropic material and an orthotropic material, with one interface condition being perfectly bonded and the other one with frictional contact, eigensolutions were obtained where the in-plane and anti-plane modes could not be decoupled. In such strange singular eigensolutions (eigenmodes) a specific tangential stress angle ω=0,±π/2,±π(defined in Chapter 3) appears in the frictional contact interface. Initially, these solutions could not be verified. However, by validating the general matrix formalism and the computational code in Chapter 3, in addition to observing such a solution, more instances where these strange (unexpected) solutions appear have been noticed. This type of solution appears also in Chap- 105 106 M.A. Herrera-Garrido ter 7, where a semi-analytical model is presented for studying stress singularities in the case of frictional contact between two wedges. Therefore, the present chapter aims to justify the existence of such unexpected singular elastic eigensolutions. Firstly, in Section 5.2.1, a review will be provided to justify why it is expected to be able to decouple the in-plane mode from the anti-plane mode when certain elastic symmetries exist in the corner. Section 5.2.2 focuses on frictional contact, which is the reason why this type of solutions may appear. Section 5.2.3 focuses on singularities analysis, and finally, Section 5.3 presents several counterexamples with eigensolutions where such decomposition is not possible due to the presence of a frictional contact zone in the problem. 5.2 GPS elastic problem with x3=0being a PES. Is decomposition into plane and anti-plane strain problems always possible? 5.2.1 Decomposition into plane and anti-plane strain problems for monoclinic materials with x3=0being a PES This way of proceeding, i.e., decomposing such corner singularity analysis into plane and anti-plane strain analysis, is based on the fact that the elastic solution of such problems can be decomposed into in-plane and out-plane parts in the bulk. This can be easily shown as follows. Consider a displacement field ui(x1,x2) (i=1,2,3)as a solution of a linear elastic Boundary Value Problem (BVP) in GPS defined in a domain Ω∈R2defining the transversal section of a cylindrical body. Let us split the components of the displacement field uito in-plane (blue) and out-of- plane (red) parts ui=  u1(x1,x2) u2(x1,x2) u3(x1,x2) .(5.1) Then, let us identify the in-plane and out-of-plane components of the tensor of small strains εi j(x1,x2) = 1 2(ui,j(x1,x2)+ uj,i(x1,x2)), εi j =   ε11 ε12 ε13 ε12 ε22 ε23 ε12 ε22 0   =   u1,11 2(u1,2+u2,1)1 2u3,1 1 2(u1,2+u2,1)u2,21 2u3,2 1 2u3,11 2u3,20   .(5.2) Let us assume that x3=0 is a PES of the considered linear elastic material. Therefore, the stiffness matrix of such material takes the following form Ci j =          C11 C12 C13 0 0 C16 C12 C22 C23 0 0 C26 C13 C23 C33 0 0 C36 000C44 C45 0 000C45 C55 0 C16 C26 C36 0 0 C66          .(5.3) CHAPTER 5107 Thus, the in-plane components of the stress tensor are defined as follows σ11 =C11ε11 +C12ε22 +C162ε12 σ22 =C12ε11 +C22ε22 +C262ε12 σ33 =C13ε11 +C23ε22 +C362ε12 σ12 =C16ε11 +C26ε22 +C662ε12. (5.4) Whereas, the out-of-plane components of σi j are defined as σ23 =C442ε23 +C452ε13 σ13 =C452ε23 +C552ε13.(5.5) Thus, the partition of the stress tensor σi j(x1,x2)into in-plane (blue) and out-of-plane (red) components is the following σi j =   σ11 σ12 σ13 σ12 σ22 σ23 σ13 σ23 σ33   .(5.6) Then, the equations of equilibrium for σi j adopt the following form in Ωwriting separately the in-plane part σ11,1+σ12,2+f1=0, σ12,1+σ22,2+f2=0,(5.7) and the out-of-plane part σ13,1+σ23,2+f3=0,(5.8) where fi(i=1,2,3)is the vector field of body forces. Combining (5.7), (5.4) and (5.2) it is easy to see that equations (5.7) include only second order derivatives of u1and u2. Whereas, combining (5.8), (5.5) and (5.2) it is easy to see that equations (5.8) include only second order derivatives of u3. Thus, we have just shown that assuming a PES x3=0 of the material, the solution of a BVP in Ωcan always be split into in-plane and out-of-plane parts. This decomposition can also be derived from the Lekhnitshki-Stroh complex variable representation of displacement ui(x1,x2)and stress function vector φi(x1,x2)[64, 65, 134], as the 3 ×3 matrices Aand Bin these representation can be decomposed into in-plane submatrix 2 (blue) and 1 anti-plane (red) sub-matrix A=   a11 a12 0 a21 a22 0 0 0 a33   ,(5.9) B=   b11 b12 0 b21 b22 0 0 0 b33   .(5.10) The expressions for these matrices can be found in [134, 6.3.7] and [134, 6.3.8] respectively. 114 M.A. Herrera-Garrido Figure 5.5: Tangential stress direction angle ωvs wedge solid angle ϑin an isotropic material with contact friction at θ0=0 with µ=1 and stress free boundary condition at θ1=ϑ◦. Values of ω′ iand ω′′ iwith i=9,10,13,14 for wedge solid angle 235◦<ϑ<243◦obtained with the developed code in Chapter 4. 5.3.1.3 Results compared with numerical methods Choosing ϑ=300◦and µ=1 it can be seen in both analysis that there are 13 pairs of λi−ωi solutions shown in Table 5.2. With the eigenequation 2.18 in [121], if f=1 we obtain the solutions λ2and λ3in Table 5.2 and with f=−1λ1,λ7and λ8, as Sinclair [121] only analyses the in-plane loading. With this configuration, 6 of the solution shown in Table 5.2 are unexpected, λ9−λ−14, apart from the out-of-plane solution, λ4−λ5. The strange singular solutions are marked with an asterisk in Table 5.2. Table 5.2: Eigensolutions found with the general matrix formalism for stress singularity problems Chapters 2 and 3 in an isotropic single-material wedge of aluminium with friction contact at θ=0◦and stress-free boundary condition at θ=300◦. Friction coefficient µ=1. i1 2 3 5 6 λi0.196813788 0.391687631 0.523820659 0.600000000 0.600000000 ωi180.000000000◦0.000000000◦0.000000000◦90.000000000◦−90.000000000◦ i9* 10* 11* 12* 13* λi0.769464958 0.769464958 0.316854623 0.316854623 0.209797619 ωi98.479800278◦−98.479800278◦65.645189499◦−65.645189499◦145.1005708073◦ i14* 7 8 λi0.209797619 0.701951807 + 0.147326550i 0.701951807 - 0.147326550i ωi−145.1005708073◦180.000000000◦180.000000000◦ To prove that the unexpected pairs of λi−ωiare in fact solutions of the problem, the analytical stress and displacement fields were computed and compared with the result of a numerical analysis. All the solution have been checked, but for sake of brevity only the analysis of the solution λ11 −ω11 is presented. The numerical method applied is BEM following the underneath steps: CHAPTER 5115 1. Define a wedge boundary given by the contact zone, an arc at the unit distance from the tip and the stress free zone. From the tip until a distance equal to 0.1 units, a bias is applied, the size of the smallest element is 10−6units, and the largest at 0.1 is 0.00346 units. From 0.1 to 1 all the elements have the same size, 0.00346 units. The elements arlong the arc are at a distance of 2◦ 2. Compute the displacements ux,uyand uzand the stresses σθθ ,σθr,σθ3,σrr and σr3us- ing the analytical solution at the defined mesh nodes. The displacement field obtained analytically with λ11 =0.316854623 is shown in Fig. 5.6. To improve the visibility, the results have been scaled. 3. Impose the computed displacements in the BEM nodes along the whole boundary and compute the stress field. 4. Verify that at the contact zone, the value of σθθ implies compression. Indeed, in this case as σθθ (θ0=0,r)<0, the compression is verified. 5. Check the fulfillment of the Coulomb law and compare the direction of the resultant tangential stress with the ωvalue obtained as the solution by the matrix general formalism. In this case, it has been verified that all the stress numerical values computed verified the following Coulomb relation σθθ (θ0=0◦,r) = qσ2 θr(θ0=0◦,r)+ σ2 θ3(θ0=0◦,r)(5.23) and that arccos σθr(θ0=0◦,r) µσθθ (θ0=0◦,r)≈65.64◦(5.24) and arcsin σθ3(θ0=0◦,r) µσθθ (θ0=0◦,r)≈65.64◦(5.25) corresponding with the ω11 value shown in Table 5.2. 6. Verify the fulfillment of the energy dissipation condition, the resultant tangential stress at the contact zone must oppose the tangential displacement. This step for this problem is not mandatory as the local fields related to the analysed strange solution we are studying might be part of a global solution that combines other local solutions and help to met the energy dissipation condition on the global problem. 7. In a bi-logarithmic scale, fit a linear function to the numerical result and compare the slope with the singular exponent solution of the problem. Fig. 5.7, where the stresses σθθ ,σθrand σθ3in bi-logarithmic scale are presented, shows an excellent agreement of the analytical and numerical eigensolution. The numerical solution present high scattering in the neartip zone with r<10−5due to the stress singularity, if these points are excluded when fitting the linear curve, the slopes for the fitted lines are close to singularity order −δ11 =λ11 −1=−0.683178812 • Linear slope for σθθ :−0.681742028 • Linear slope for σθr:−0.681819112 • Linear slope for σθ3:−0.681722431 116 M.A. Herrera-Garrido Figure 5.6: Displacement of the nodes at the perimeter of the isotropic wedge with r=1, material A in Table 5.1, with the stress free boundary condition at θ=300◦and frictional contact at θ=0◦with friction coefficient µ=1. 5.3.2 Orthotropic single-material wedge In this section, the same procedure applied in the previous section is repeated but the wedge material in Fig. 5.1 is the orthotropic material C in Table 5.1. 5.3.2.1 Influence of the friction coefficient As in the previous section, an attempt is made to find out whether the friction coefficient affects the number of unexpected singularities and in particular their presence. In Fig. 5.8 the variation of the singular exponent λvs the friction coefficient µis represented. The following conclusion can be observed: • No strange solutions are found in the frictionless case. • As in the case of the single isotropic material, when increasing the friction coefficient from the frictionless case, a solution with unexpected sliding angle appears, point A in Figs. 5.8 and 5.9. If the friction coefficient continues increasing, the singularity exponent remains constant but the sliding angle changes. If the friction coefficient continues increasing, more strange singularities appear, point B, C, D and E in in Figs. 5.8 and 5.9. • As in the isotropic case, due to the elastic symmetric with respect the coordinate plane x3=0, the asymmetric strange solutions appear in pairs of eigensolutions, which are reflection of each other with respect to the plane of material symmetry x3=0, i.e. λ15 =λ16 and with ω15 =−ω16. • At point F in Fig. 5.8, µ≈0.28, λ(F) 1=λ(F) 2and two solution with λ∈Rbecome conjugate complex for all µ⪆0.28. CHAPTER 5117 Figure 5.7: Comparison between the analytical and the numerical stresses on the frictional contact boundary at θ=0◦of an isotropic wedge, material A in Table 5.1, with stress free boundary condition at θ=300◦. Friction coefficient µ=1. • From all the solutions shown in Figs. 5.2 and 5.3, after checking the stress and displacement fields, it is possible to say that only λ6−ω6,λ7−ω7,λ9−ω9and λ10 −ω10 fulfill the energy dissipation condition. These solution fulfill the energy dissipation condition for all the values shown in Figs. 5.2 and 5.3. • The solutions corresponding to λ4−ω4and λ5−ω5are the antiplane solutions and have no compression in the contact zone, hence these solutions represent the frictionless case even with µ>0 and are independent of µ. • Looking to Fig. 5.9, it can be seen that, for higher friction coefficients, non-expected angles ωtend to ±90◦. This aspect was not observed in the isotropic case. 5.3.2.2 Influence of the material solid angle Also for orthotropic material it is interesting to analyse how the wedge solid angle affects the existence of these non-expected singular solutions with sliding angles different to 0◦,±90◦ and ±180◦. In Figs. 5.10 and 5.11 are represented the results for the singular exponent λand the direction of the tangential stress angle ωrespectively, as a function of the wedge solid angle ϑdefined in Fig. 5.1, starting in the semi-plane case with ϑ=180◦and ending with the crack case ϑ=360◦. The schema of the solved problem is presented in Fig. 5.1 with a friction coefficient µ=1 at ϑ=0◦. The engineering elastic properties of the orthotropic material are the same as used in the previous section, material C in Table 5.1. From the results shown in Figs. 5.10 and 5.11 similar remarks as in Section 5.3.1.2 can be drawn 118 M.A. Herrera-Garrido Figure 5.8: Singularity exponent λvs friction coefficient µin a 300◦orthotropic wedge, made of material C in Table 5.1, with frictional contact at θ0=0◦and stress free boundary condition at θ1=300◦. Figure 5.9: Direction of the tangential stress angle ωvs frictional coefficient in a 300◦orthotropic material with contact friction at θ0=0 and stress free boundary condition at θ1=300◦. • Strange solutions can appear for a ϑangle lower than 180◦. As saw in the previous section, with a lower frictional coefficient, the non-expected singularity exponents λ may not appear. • Increasing the wedge solid angle ϑchanges the singularity exponent λ, and in cases of non-expected singular solution, also implies changes in the direction of the tangential stress angle ω. CHAPTER 5119 Figure 5.10: Singularity exponent λvs solid angle ϑin deg. in an orthotropic wedge, made of material C in Table 5.1, with frictional contact at θ0=0◦with µ=1 and the stress free boundary condition at θ1=300◦. Values of λ13,λ14,λ17 and λ18 for wedge solid angle 316◦<ϑ<343◦obtained with the developed code in Chapter 4. • Complex singularity solutions λ∈Cappear when two different solution have the same singularity exponent λi=λjand direction of the tangential stress ωi=ωjfor i= j. If the solutions have an expected direction of the tangential stress angle ω, the complex solution have also one expected direction of the tangential stress. For example in Figs. 5.10 and 5.11 solution λ1−ω1and λ2−ω2are overlapped in the interval 232◦<ϑ<345 with ω=180◦. In the other hand, when the directions of the tangential stress of both solution that are overlapped are unexpected, the real and imaginary part of the solution have different directions of the tangential stress angles ω′and ω′′. For example, when the unexpected solution λ13 −ω13 (and λ14 −ω14) is overlapped by the solution λ17 −ω17 (λ18 −ω18) from 316◦<ϑ<343◦, both solution become conjugate complex and two sliding angle ω′and ω′′ appear. These solutions are computed with the generalization for complex solution in presence of frictional contact of the general matrix formalism presented in Chapter 4. 5.3.2.3 Results compared with numerical methods In Table 5.3 are presented all the results found for the singularity problem represented in Fig. 5.1 for an orthotropic wedge, material C in Table 5.1 with ϑ=300◦and µ=1. With the orthotropic material and this configuration, 5 pairs of unexpected solutions appear, marked with an asterisk (*) in Table 5.3. Only 4 out of the 18 presented solutions fulfill the energy dissipation condition in the frictional contact boundary, these solutions are λ6,λ7,λ9and λ10. To verify the non-expected solutions in Table 5.3 the same methodology to compare the analytical solution with a numerical solution presented in Section 5.3.1.3 is followed. In this study, for sake of brevity only the verification made for λ9=0.316545653 is presented. Following step 1 and 2, the displacement field obtained analytically with λ9= 0.316545653 is shown in Fig. 5.12, scaled to improve visibility. 120 M.A. Herrera-Garrido Figure 5.11: Direction of the tangential stress angle ωvs solid angle ϑin deg. in an orthotropic material with frictional contact at θ0=0◦with µ=1 and the stress free boundary condition at θ1=300◦.Values of ω′ iand ω′′ i with i=13,14,17,18 for wedge solid angle 316◦<ϑ<343◦obtained with the developed code in Chapter 4. Table 5.3: Eigensolutions found by the general matrix formalism for stress singularity problems (Chapters 2 and 3) in an orthotropic single-material wedge Fig. 5.1, made of material C in Table 5.1, with friction contact at θ=0◦and the stress-free boundary condition at θ=300◦. Friction coefficient µ=1. i 1 2 3 λi0.720418063+0.127401596i0.720418063-0.127401596i0.228460059 ωi180.000000000◦180.000000000◦180.000000000◦ i 4 5 6 7 9∗ λi0.600000000 0.600000000 0.570130765 0.351600130 0.316545653 ωi90.000000000◦−90.000000000◦0.000000000◦0.000000000◦58.4378279◦ i 10∗11∗12∗13∗14∗ λi0.316545653 0.258154147 0.258154147 0.688166081 0.688166081 ωi−58.4378279◦115.826327◦−115.826327◦98.31451952◦−98.31451952◦ i 15∗16∗17∗18∗ λi0.884911910 0.884911910 0.832422905 0.832422905 ωi80.41987908◦−80.41987908◦93.58249824◦−93.58249824◦ Those displacements are imposed in the BEM nodes to compute the stresses as indicated in step 3. The existence of the solution is verified by checking the compression in the friction contact zone and the fulfilment of the Coulomb law, step 4 and 5. In this case, all the numerical values of σθθ ,σθrand σθ3verifies the following relations σθθ (θ0=0◦,r)<0,(5.26) σθθ (θ0=0◦,r) = qσ2 θr(θ0=0◦,r)+ σ2 θ3(θ0=0◦,r),(5.27) CHAPTER 5121 Figure 5.12: Displacements corresponding to λ9in Table 5.3 for nodes on the perimeter of an orthotropic singlematerial wedge (Material C in Table 5.1) with friction contact at θ0=0◦with µ=1 and the stress-free boundary at θ1=300◦Fig. 5.1. and that arccos σθr(θ0=0◦,r) µσθθ (θ0=0◦,r)≈58.44◦,(5.28) and arcsin σθ3(θ0=0◦,r) µσθθ (θ0=0◦,r)≈58.44◦.(5.29) corresponding with the ω9value shown in Table 5.3. In this case, as it was pointed out for the corresponding analytical solution in 5.3.2.1, the numerical solution fulfills the energy dissipation condition. Nevertheless, as in the previous example, the local fields related to λ9 solution might be part of a global solution that combines other local solutions and help to met the energy dissipation condition on the global problem. Following step 7, the numerical and analytical stress fields at the frictional contact are presented in a bi-logarithmic scale in Fig. 5.13, showing excellent agreement. Similar to the previous example the near-tip solution, first 3 points in the curves, exhibits high scattering. If these points are excluded when fitting the linear curve, the slopes for the fitted lines are close to singularity order −δ9=λ9−1=−0.68348779 • Linear slope for σθθ :−0.6821795327, • Linear slope for σθr:−0.6821259201, • Linear slope for σθ3:−0.6794944953. 122 M.A. Herrera-Garrido Figure 5.13: Comparison between the analytical and the numerical stress fields corresponding to λ9in Table 5.3 on the frictional contact boundary at θ=0◦of an orthotropic single-material wedge (Material C in Table 5.1) with friction contact at θ0=0◦with µ=1 and stress-free boundary at θ1=300◦Fig. 5.1. 5.3.3 Isotropic bi-material corner In the previous sections it has been prove that unexpected solutions can appear in the most simple cases, isotropic or orthotropic single-material corners in frictional contact with a rigid obstacle. In this section it is checked that these unexpected solution also appear in presence of frictional contact between two isotropic materials. Fig. 5.14 outlines the problem configuration. The material used in the analyses performed below is material A from Table 5.1 for both sectors in the corner. Figure 5.14: Bi-material problem schema. 5.3.3.1 Influence of the frictional coefficient Following the same structure as in the previously studied cases, the initial analysis enables us to verify whether the friction coefficient at the interface between both materials has a sim- CHAPTER 5123 ilar effect regarding the occurrence of unexpected solutions to that observed for the friction coefficient at the frictional boundary. Figs. 5.15 and 5.16 illustrate the variation of the singularity exponent λand the direction of the tangential stress angle ω, respectively, in relation to the friction coefficient µ. It can be seen a notable correlation between the friction coefficient and the occurrence of unexpected solutions. The following remarks can be done: • No strange solutions are found in the frictionless case. • The first pair of unexpected solutions λ6and λ7need a higher value of the frictional coefficient that in the single-material case. Also in the bi-material case, the singularity exponent λvalue remains constant when increasing the frictional coefficient, but the direction of the tangential stress varies. • As both wedges are made of the same material, the symmetry at the corner persist. Therefore, as in the cases of single-material scenarios, unexpected asymmetric solutions should emerge symmetric pairs of asymmetric solutions, with one being the reflection of the other with respect to the symmetry plane x3=0, i.e ω6=−ω7. • The values of ω8and ω9are not represented in Fig. 5.16 as λ8=λ9=0.5 is a valid solution combined with any ωfor any friction coefficient µ. These solutions correspond with the stick case where no relative displacement appears between both materials. • Solution λ1and λ2correspond with the antiplane solution with no compression on the frictional contact interface. •λ3,λ6and λ7fulfill the energy dissipation criterion in the interval 0 <µ≤1, while λ4does not fulfill this criterion for any value in the interval 0 <µ≤1. λ5fulfill the energy dissipation criterion in the interval 0 <µ≤1 only when λ5<1. 5.3.3.2 Influence of the second material solid angle In the case of the bi-material corner Fig. 5.14, we have the option to vary two solid angles, θ0and θ2. In order to maintain a similar analysis to that conducted in the single-material case, only the angle θ2will be varied. Maintaining θ0=−60◦and θ1=0◦and starting from the semi-plane scenario with θ2=180◦, we gradually cover the entire corner, reaching θ2=θ0+360◦=300◦. In Figs. 5.17 and 5.18 are represented the results of the singularity exponent λand the direction of the tangential stress ωas a function of the second material solid angle θ2=ϑin Fig. 5.14. In this analysis the friction coefficient is set to µ=1. Based on the data analysis, the following observation can be listed • Only one pair of unexpected solutions, λ6and λ7, is found. With ϑ<292◦these solutions become λ6=λ7=1, solutions that are already present for any ϑsolid angle and any direction of the tangential stress. For that reason there are not ω6,ω7values for ϑ<292◦in Fig. 5.18. • Increasing the solid angle implies changes in the singularity exponent, and in the cases of strange singularities, changes in the tangential stress direction too. 226 M.A. Herrera-Garrido Figure 9.9: Crack length increment in ENF test under critical load in mm vs the elastic modulus of the stiffer material E2for different interface crack coefficient of friction µ escalates rapidly with an increase in the elastic modulus of the stiffer material. This behavior remains consistently similar, with a slight variation in slope, across all analyzed values of the coefficient of friction. 9.4.4 Relation of stress singularity with previous results In this section, the singularity exponent λis computed using the developed code presented in Chapters 2 and 3 across the range of coefficient of friction and elastic modulus values employed in the previous sections. The objective is to investigate whether any discernible relationships exist with the observations outlined in the preceding sections and the singularity exponent too. Figure 9.10: Singularity exponent λvs coefficient of friction µfor different elastic moduli of the stiffer material In Fig. 9.10, the relationship between the singularity exponent λand the coefficient of CHAPTER 9227 Elastic modulus E2(GPa) 1 2 3 5 7 10 15 Dundurs parameter β0 0.09523 0.14286 0.19048 0.21428 0.23377 0.25 Table 9.1: Dundurs mismatch parameter for isotropic bi-materials with E1=1 GPa and ν1=ν2=0.3. friction µis depicted for each value of the elastic modulus of the stiffer material in the bimaterial, E2. An apparently linear relationship is observed, with lower slopes obtained for smaller values of E2. This behavior is akin to that described in Fig. 9.6. As discussed in previous chapters, in the case of isotropic bi-materials, Comninou [32] derived the following analytical expression to obtain the singularity exponent for a frictional interface crack: λ=1−1 πacot(β µ).(9.16) This trigonometric function, for small values of µ, behaves similarly to a linear relationship with a slope of β π. Here, βis the Dundurs second mismatch parameter for isotropic bimaterials [44], that with E1=1 and ν1=ν2=0.3 takes the values specified in Table 9.1 for each E2used in the examples. This can explain the variation between the slopes in Figs. 9.6 and 9.10. Ultimately, it can be concluded that, based on Figs. 9.6 and 9.10, the critical load is apparently linearly related to the coefficient of friction µ, and this relationship may be governed by the Dundurs mismatch parameter β. 9.5 Concluding remarks In this chapter, it has been demonstrated that it is feasible to apply the coupled criterion of finite fracture mechanics to the three-point end notched flexure test considering friction contact between interface crack faces. A numerical methodology has been developed and implemented for various values of the coefficient of friction and the elastic modulus of one of the materials. Consequently, it has been numerically verified that both parameters, the coefficient of friction and the elastic moduli mismatch directly influence on both the critical load Pcand the crack length increment lc. Furthermore, similar to the singularity exponent λ, it has been observed that the critical load Pcmaintains a approximately linear relationship with the coefficient of friction for µ≤1. In the case of the singularity exponent λ, this approximately linear relationship exhibits a slope associated with the Dundurs second mismatch parameter β[44], which is directly dependent on the discrepancies of the bi-material elastic moduli. On the other hand, once the critical load Pcis obtained, the crack length increment lcfor that load can be analyzed. Both an increase in the coefficient of friction and an increase in the discrepancy between the elastic moduli have been found to correlate with an increase in the crack length. Nevertheless, it has been noted that the growth of the crack is more significantly affected by changes in the elastic moduli mismatch than by variations in the coefficient of friction µ. With this chapter, the Objective 5 of this thesis is achieved. Thus, all the objectives set at the beginning of this research work have been accomplished. Part VI From Reflection to Projection 229 CHAPTER 10 Summary and conclusions 10.1 Summary The main objective of this thesis was to propose a methodology for characterising the singular elastic solutions in the neighborhood of singular points by an asymptotic series expansion of power-law terms. Once the characteristic exponents of this asymptotic series are obtained, it would only be necessary to accurately calculate the GSIFs using one of several methods proposed in the literature. Once the GSIFs are determined, and thus the singular solution around the singular point, it can be leveraged for various purposes, but primarily to predict damage initiation. To provide an analysis of singular stresses, after conducting a state-of- the-art study, various objectives were proposed to address the observed open questions, all of which have been successfully achieved during the completion of this thesis. During the literature review, it was observed that among all the methodologies reviewed, the one that seemed most versatile and accurate, offering the possibility of solving a greater number of cases, was the semianalytical one proposed by Mantiˇ c et al. [96, 97], which in turn was an extension of the proposal in [11]. This methodology was neither verified nor computationally implemented previously, so the first specific objective was its review, implementation, and verification. As it is such a versatile tool, for its proper analysis and verification, this objective was divided into the Chapters 2 and 3. First conducting an analysis, implementation, and verification for all cases that did not involve frictional contact in Chapter 2, and in Chapter 3, an exhaustive analysis and verification of the general matrix formalism part involving friction was performed. In Chapter 3, a deficiency of the matrix formalism was observed when obtaining solutions in cases in which the singularity exponent λ∈Cin presence of friction contact. Therefore, in Chapter 4, a proposal for generalizing this matrix formalism is made to obtain valid and useful solutions with λ∈Cin presence of friction contact. Throughout the verification of the general matrix formalism and its implementation, unexpected elastic solutions were observed. These solutions in corners with elastic symmetries are solutions in which the in-plane and anti-plane modes cannot be decoupled. In the literature, only mention of such solutions was made in [96, 97], but no verification of these was presented. Therefore, in Chapter 5, to test the existence of such solutions, the analytical so- 231 232 M.A. Herrera-Garrido lutions obtained with the singularity analysis code are compared with the boundary element method. Furthermore, it was observed in the literature that most authors, when solving stress singularity problems involving frictional contact, except for Chen et al. [27, 28], despite solving the problems under generalized plane strain, apply the Coulomb law with certain restrictions that limit the freedom of the sliding direction. Therefore, it was deemed necessary to also develop simpler semi-analytical methodologies to solve specific cases. Among all the possible cases to analyze, the most studied ones in the literature were chosen: interface crack with frictional contact and frictional contact between two wedges. These objectives were developed in Chapters 6 and 7 respectively. Thanks to the methodologies developed in Chapters 6 and 7, the general matrix formalism presented in Chapters 2 and 3, was fully verified by being able to compare and obtain identical solutions in cases where the direction of the tangential stress vector angle can be different from 0◦or 90◦. However, this could not be achieved by comparing with the methodologies or numerical results found in the literature since, apart from Chen [27, 28], no other author studied solutions in which, in view of the lack of elastic symmetries of the anisotropic materials included in the corner, due to the elastic non-symmetries of the corner in presence of anisotropic materials, the in-plane and anti-plane modes cannot be decoupled. Once the code is fully verified and its great versatility is observed, it is proposed to offer this tool to the scientific community through a web application. The structure and performance of this application are presented in Chapter 8. Finally, due to the significant challenges posed by including friction in the developed analysis, emphasis is placed on the need to consider frictional contact when solving stress singularity problems. To observe the effect of friction in a real case study, the analysis presented investigates how the friction coefficient affects the maximum load applicable in the three-point bending test with a pre-crack using the coupled criterion. 10.2 Conclusions In this section, the intention is to delve deeper into the conclusions already mentioned in each chapter for each particular objective. Therefore, this section will be subdivided according to the objectives set in Section 1.3 10.2.1 Review and verification of the General Matrix Formalism and its implementation in a Matlab computational code This thesis reviews and complements the general matrix formalism and the computational procedure introduced in a very concise way in [96, 97], in order to be easily understood by the readers and users. However, the most relevant contribution of the present thesis regarding this objective is that it presents the first general implementation of the proposed procedure extensively and successfully tested by many numerical examples. The main features of the general matrix formalism for stress singularity problems and of the developed code are: - Based on the following ingredients * Stroh formalism for anisotropic elasticity [123, 124] CHAPTER 10 233 * Transfer matrix for a single-material wedge proposed by Ting [135] * Matrix formalism for homogeneous orthogonal boundary and interface conditions proposed by Mantiˇ c et al. [96, 97]. - Compute the singularity exponent λ, the angle ωof the friction tangential stress vector on each frictional contact surface (in cases including frictional contact), and the associated singular stress and displacement fields. - Applicable to corners with any finite number of single-material wedges - Applicable to any linear elastic and homogeneous material in each wedge - Able to analyse both open and closed (periodic) corners - Considers Generalised Plane Strain (GPS) ui(x,y),i=x,y,z, a kind of 2.5D elastic problem formulation - Considers homogeneous boundary conditions covering stress free and fixed faces, or faces with some restricted or allowed direction of displacements, defined either in the reference frame aligned with the cylindrical coordinate system or in an inclined reference frame - Any finite number of frictional surfaces, outer boundary surfaces or interfaces, in the corner, based on the 3D Coulomb law for isotropic friction represented by the Coulomb cone. - Collinearity of friction shear vector and the relative displacement vector is satisfied. - The frictional dissipation of energy condition and compression condition are checked. This implementation allows us to verify the correct performance of the proposed matrix formalism too. The code underwent validation through extensive testing, comparing its results with those presented by other authors or derived from analytical expressions proposed in literature. Consequently, we were able to validate methodologies and analytical expressions proposed by other authors in parallel.Throughout this validation process, discrepancies emerged between our computational results and those documented in the literature, alongside discovering results not previously reported. Rather than casting doubt on the code-derived results, these discrepancies served to highlight its efficacy by producing valid outcomes not attainable through alternative methodologies. The underlying reasons for these discrepancies are categorized in Table 10.1. After its verification, the following advantages of the Matlab code can be enumerated: - Versatility, since it can solve a multitude of problems with essentially any useful boundary and interface conditions, covering cases with any finite number of frictional sliding contact surfaces, and any linear elastic and homogeneous material in each wedge. Alongside the boundary conditions proposed in [96, 97], in Chapter 2, the boundary condition of symmetry and antisymmetry with inclined planes is also proposed as a novelty. This demonstrates how straightforward it is to expand the current formalism to provide greater versatility. This makes it a highly useful tool for addressing a multitude of stress singularity problems. 234 M.A. Herrera-Garrido No frictional contact Frictional contact Isotropic or orthotropic materials (All coordinate planes are PES) The literature only presents inplane mode results. In the literature a simplified Coulomb law assumes that the solution can always be decomposed into in-plane and antiplane modes . Other cases There is a disagreement found with the numerical results presented in [28] for the frictionless sliding contact case, but there is an agreement with those obtained by applying their methodology presented in [28]. In the literature the authors apply\propose simplifications to the 3D Coulomb law. There is a disagreement with the numerical results presented in [28] for the frictional sliding contact case, but there is an agreement with those obtained by applying the methodology in [28]. Table 10.1: Classification of the differences found between the results obtained with the present computational code based on the general matrix formalism and the results found in the literature. - Accuracy and reliability. It can be seen that all the numerical results computed by the code developed in the present work are essentially identical, at least up to 10 digits shown, to the results obtained by the numerical solution of closed-form eigenequations found in the literature. This excellent accuracy is a consequence of the semianalytic character of the present procedure, where numerical solution is used to find roots of the sextic Lekhnitskii-Stroh characteristic polynomial of an anisotropic material if its solution in radicals is not possible, and to find roots of the transcendental eigenequation of the corner problem, where we are looking for roots of a complex-valued analytic (holomorphic) function of one complex variable. The rest of the calculations are fully analytical, providing the maximum accuracy to the present procedure. In fact, an arbitrary accuracy can be achieved using modern computer algebra software when working with numbers stored with an arbitrarily high precision. - Efficiency, despite its complexity, the code maintains efficiency in computation, making it suitable for practical engineering applications. - Novelty, the code has discovered new results not previously documented in the literature. - Ease of use and extension for possible consideration of cases not considered in the current code version, since the code is programmed by easily modifiable modules. This is demonstrated in this document, since starting from the developed code presented in Chapter 2, first in Section 3.4 and later in Section 4.3 it is detailed the changes made in the code to add the possibility of solving problems with frictional contact, and to obtain valid singularity solutions with λ∈Cwith frictional contact. - Robustness, several self-checks can be performed: * When rotating the corner, all ϑwill change, and with them also all relevant matrices. Despite this, the final results for λ, and w(r,ϑ)(also ωin case of CHAPTER 10 235 considering frictional contact) have proven to be the same. * For a multi-material wedge, where there are one or more perfectly bonded interfaces, the formalism can be applied by using either the ‘perfectly bonded interface condition’ or exploiting the ‘wedge transfer matrix’, and the result must be the same in both cases, even having substantially different sizes of the corner characteristic matrix. The frictionless interface case can be studied using the frictionless interface condition or applying a friction coefficient µ=0. The size of the matrices and the solution method of the characteristic system are different, but the results obtained for λand w(r,ϑ)are the same. * When applying a friction coefficient µ=0 to a boundary and when applying the symmetry boundary condition (only uθrestricted), despite having different sizes of the characteristic matrix, the results are the same in both cases. * In the present formalism, the use of µ<0 is not intended. However, it has been used for verification purposes. The problem is first solved using the friction coefficient µ>0. If λ1,ω1is a solution to the problem, then when solving it again with the friction coefficient −µ, the solution should be λ1,ω2, where ω2=ω1+π. Regarding the complex singularity exponent λ∈C, leading to oscillatory singular solutions, in presence of frictional contact, as shown in the numerical results in Section 3.5, in cases of frictional contact between isotropic materials, or orthotropic materials with a PES x3=0, potentially valid solutions for in-plane mode with ω=0◦or 180◦and λ∈Cmay appear. In the case of frictionless contact on one face of an anisotropic material and complex λ∈C, it was observed by studying the solution in displacements that there are two different relative direction of the displacement vector angles ω′ uand ω′′ uon the surface in contact. Assuming a continuous variation of the friction coefficient from zero to positive values, it is considered that the values of ω′ u,ω′′ uand λvary continuously. Therefore, another procedure is needed for modelling the frictional contact, which considers a complex λand two values of ω′ uand ω′′ uas unknowns on each frictional surface. These cases of complex λare covered in Chapter 4. However, these results could only be verified through self-checks, as no procedure has been found in the literature that allows for an analysis of these characteristics, in which not one, but two angles of the direction of the tangential stress, corresponding to the real and imaginary solution, are considered as unknowns. Surprisingly, these solutions with complex λ∈Cand two values of ω′ uand ω′′ uhave also been observed in Chapter 5 for isotropic and orthotropic materials with PES x3=0. The concept of a small-scale tension zone has been proposed to give a physical meaning and utility to complex solutions involving oscillatory stress fields for λ∈C. It has been observed that the size of this zone depends on the imaginary value of the singularity exponent λ∈C, which in turn can be affected by the anisotropy of the materials, the friction coefficient, or the problem geometry. The developed code is a general tool useful for researchers who need to know the singularity exponents and the singular stress fields to improve or check their numerical results by FEM, or to those researchers who need to verify their analytic formulas or eigenequations developed for specific stress singularity problems. 242 M.A. Herrera-Garrido shown in [27] does not match with the results obtained with the present methodology and consequently neither with the results obtained with the matrix formalism presented in Chapters 2 and 3. Chen [27] presents a sliding angle ˜ ω=−4.186◦and λ=0.5 as the solution at the tip of the contact zone in a bi-material composed by two laminates with stacking sequence 45◦/−45◦. If we solve this same problem we can check that both obtained solutions with the methodology presented in Chapter 7 and with the matrix formalism in Chapters 2 and 3, the sliding angle ˜ ω=±90◦with λ=0.5 and ˜ ω=0◦or 180◦also with λ=0.5. The same example was solved also with the methodology proposed by Ting and Chou [136], obtaining a solution with ˜ ω=±90◦and λ=0.5 as well. 6) Apart from the results presented by Chen et al. [27, 28], no other methodology has been found in the literature that does not impose the sliding angle between the two materials in frictional contact. In the case of free-free boundary conditions, both the results obtained with the fixed geometry of two semi-planes presented by Chen [27] and the methodology valid for any geometry presented by Chen et al. [28] do not coincide with the results obtained through the general matrix formalism presented in Chapters 2 and 3 or with the results obtained through the methodology presented in Chapter 7. However, when comparing the results obtained through the matrix formalism with those obtained through the specific methodology for analyzing stress singularities in the frictional contact problem between two materials, perfect agreement is observed in all decimal places. The same agreement holds for the rest of the boundary conditions available in both methodologies. Thus, both the methodology presented in Chapter 7 and some features of the general matrix formalism (Chapters 2 and 3) that could not been compared before have been verified. 10.2.5 Integration of the developed general Matrix Formalism into Web Platform At this point, the general matrix formalism has been thoroughly verified through comparison with numerical results or expressions presented by other authors, as well as with results obtained through the methodologies presented in Chapters 6 and 7, and even through comparison with numerical methods in Chapter 5, demonstrating its significant utility. Therefore, the idea of offering this tool to the scientific community through a web application, called SingSol is proposed here. To achieve this, the Matlab code developed and tested in Chapters 2 and 3 has been completely translated into Python. This makes it easier to adapt it to an online tool without relying on third-party licenses. The Python implementation of the matrix formalism contains all the advantages mentioned in Section 10.2.1, but it also includes the following enhancements: -User-friendly interface: SingSol application stands out for its user-friendly and highly intuitive interface. The platform guides users through a step-by-step process, enabling them to input problem details efficiently. From choosing between an open corner or closed corner to the detailed specification of materials, boundary conditions, and interfaces, SingSol dynamically adapts the form based on user selections. These characteristics ensure that the data entry process is clear and straightforward. The user-friendly interface makes SingSol an accessible tool even for those not experts in the field, providing an efficient user experience. CHAPTER 10 243 -Object-oriented programming (OOP): The Python code has been programmed using OOP instead of procedural programming like the Matlab code. This improves code organization, modularity, and reusability by encapsulating data and functions into objects. Additionally, OOP promotes easier maintenance, scalability, and collaboration among developers. -Optimization module: Thanks to the modular structure of the code inherited from its Matlab version, no modifications to the base code were necessary to add the optimization module. This module enables obtaining singular solutions in corners with one or more frictional contacts without requiring the user to specify initial search values for λand ωin these solutions. For the development of this new optimization module, various parameters affecting the automatic acquisition of solutions to the problem were analyzed: • Objective function: Among the three methods proposed in Chapter 3 to solve the system of nonlinear equations, it is concluded that to efficiently obtain all results, it is necessary to minimize twice first using (3.60) and then perform a least squares adjustment using (3.61). • Mesh size: The mesh size must be made fine enough to obtain all the roots λand ω, but not too fine to reduce the unnecessary computation times. The preset value is a mesh refinement of 32 ×32 in case of F=1, and in general 32 ×32 ×···×32 | {z } Ftimes . • Filters definition: Different filters are defined to eliminate combinations of values λ and ωthat will not lead to optimal solutions of the problem. These values must be restrictive enough to eliminate combinations that are far from being optimal solutions, as analyzing them would result in a loss of computation time. However, they should not be too restrictive to eliminate points that could potentially converge to an optimal value. • Optimization methods: When comparing the minimization methods Nelder-Mead, Powell, BFGS, and SLSQP available in the Python Scipy library for minimizing the objective function (3.60), it is observed that, using a 32x32 mesh, only the Nelder- Mead method successfully obtains all solutions to the given problem, Example 3.1.1. While the Powell method may offer slightly better results in σmin, the Nelder-Mead method provides practically identical λand ωvalues in significantly less time. Therefore, for low-dimensional problems, the Nelder-Mead method is preferred due to its efficiency in obtaining all solutions promptly. However, for higher-dimensional problems, methods like Powell or BFGS with denser meshes may be more suitable. All these parameters have been configured with default values, which, based on the results verification conducted in Section 8.5, are deemed optimal. This verification involved comparing the automatically generated results from the Python code with those presented in Chapter 3. However, given the broad spectrum of problems this tool can tackle, these default settings may not be optimal for every scenario. Therefore, experienced users might need to fine-tune some or all of these parameters to better align with the specific requirements of their problem. 244 M.A. Herrera-Garrido 10.2.6 Analysis of the significance of considering frictional contact in real engineering scenarios Based on the extensive analysis presented in Chapter 9 regarding the application of the Coupled Criterion (CC) of Finite Fracture Mechanics (FFM) to predict the growth of a crack along an adhesively bonded interface, several key conclusions can be drawn: -Feasibility of CC Application: The study has successfully demonstrated the feasibility of applying the CC of FFM to predict crack initiation and propagation along adhesively bonded interfaces, specifically in the context of the three-point end-notched flexure (3P ENF) test configuration. This marks a significant step forward in understanding the behavior of frictional interface cracks and their influence on the reliability of adhesive joints. -Influence of Friction Coefficient: The analysis has revealed that the coefficient of friction between crack surfaces plays a substantial role in determining the critical load for crack propagation. Notably, a linear correlation has been observed between the critical load and the coefficient of friction. -Impact of Material Mismatch: The mismatch in elastic moduli between the materials forming the specimen also significantly affects the critical load and crack length increment. Particularly, when there is a substantial difference in elastic moduli, the critical load exhibits more pronounced variations. This highlights the importance of considering material properties when assessing the reliability of adhesive joints. -Crack Length Increment: The increment in crack length under critical load has been found to scale rapidly with an increase in the elastic modulus of the stiffer material. This observation underscores the critical role played by material stiffness in determining crack propagation behavior. -Singularity Exponent Relation: The analysis of the singularity exponent λhas revealed an approximately linear relationship with the coefficient of friction µ, with lower slopes obtained for smaller mismatch between the values of the elastic modulus of the bi-material. This relationship is governed by the Dundurs mismatch second parameter β. Given that the behavior of the singularity exponent with respect to the coefficient of friction and the difference in elastic moduli in the bi-material closely resembles the behavior of the critical load under these same two factors, it could be inferred that the nearly linear relationship between the maximum load and the coefficient of friction is also governed by the second Dundurs parameter. However, to confirm this hypothesis, further detailed analysis and specific experiments evaluating the relationship between these parameters in the particular context of the studied system would be necessary In summary, the comprehensive analysis presented in this chapter provides valuable insights into the complex mechanics of frictional interface cracks and their implications for the reliability of adhesive joints. By elucidating the relationship between critical load, material properties, and friction coefficient, this study contributes to advancing the understanding of fracture mechanics in adhesive joints applications. CHAPTER 11 Future developments 11.1 Introduction As was done in the previous chapter, this chapter is organized according to the objectives of this thesis. Relating each future development to its corresponding objective, this section aims to provide a comprehensive overview of the anticipated advancements within the scope of this thesis. 11.2 Review and verification of the General Matrix Formalism 11.2.1 Anisotropic frictional contact The theoretical formulation of the matrix formalism introduced in [96, 97] also considers anisotropic friction, which means that the sliding direction and the frictional shear direction can be different, and then as a result the friction coefficient µvaries as a function of these directions [102, 115, 149]. The computational implementation of this theoretical formulation in the present code is quite straightforward, and it would allow characterising stress singularities in the configurations where this direction-dependent friction anisotropy is relevant. 11.2.2 Power-logarithmic stress singularities It can be expected that the semi-analytic nature of the present procedure will allow a generalisation to also cover power-logarithmic stress singularities following the strategy in complex variable developed in [122], but the need to include the angle of the frictional shear vector among the unknowns will require some modifications of the strategy in [122]. 11.2.3 Non-homogeneous materials In the state of the art of this thesis, it was mentioned that authors such as Carpinteri and Paggi [22] include the possibility of non-homogeneous materials in their analysis of stress singularities in multi-material corners. This type of material refers to those that, when represented 245 246 M.A. Herrera-Garrido as wedges in a polar system, exhibit a variation in elastic properties with respect to the angle θ, but these elastic properties remain independent of the distance to the vertex. Although the matrix formalism presented by Mantiˇ c et al. [96, 97] does not directly address this possibility, its ability to analyze multi-material corners formed by any number of materials suggests the possibility of discretizing the evolution of the elastic properties of these non-homogeneous materials through multiple materials, each with elastic properties corresponding to its position in the non-homogeneous sector. Although it would be possible to use this multi-material sector approach without affecting the size of the matrix, the complexity of the transfer matrix would be increased due to multiple matrix products. Therefore, it is proposed to study this type of material as another special type, without resorting to discretization into multiple materials with different elastic properties by using the formalism proposed by Kirchner [75]. 11.3 Development of Specific Semi-Analytical Methods 11.3.1 Singularity exponent λvalues not accounted for by the developed semi-analytical methodology In Chapter 7, certain corner singularity cases were identified where values of the singularity exponent λ, such as λ∈Z, need to be studied as special cases. This is because, for instance, such values may lead to the indefiniteness of certain matrices within the system. While these cases were not included in Chapter 7, it would be useful to carry out an analysis of these cases in a manner similar to that performed in Chapter 6. 11.3.2 Complex solutions Currently, the semi-analytical methodologies presented in Chapters 6 and 7 are tailored to singular solutions where λ∈R. By adapting the approach of generalization for complex solutions developed in Chapter 4, we have the potential to expand the applicability of these methodologies to cover the solutions with complex singularity exponents λ∈Cas well. This expansion would equip us with two robust methodologies, providing a comprehensive framework for validating the outcomes obtained in Chapter 4. 11.3.3 Community sharing Following the example of the general matrix formalism introduced in Chapters 2 and 3, which has been made accessible to the scientific community via the tool presented in Chapter 8, the methodologies devised in Chapters 6 and 7—already implemented in MATLAB code for verification purposes—could likewise be shared with the scientific community. To achieve this, while the MATLAB files could be made available on platforms such as GitLab or GitHub, there’s a proposal to convert them into Python code. This would enable users to avoid dependencies on third-party licenses. CHAPTER 11 247 11.4 Integration of the developed General Matrix Formalism into Web Platform SingSol 11.4.1 Complex solution in Web App The complex singular solutions obtained in Chapter 4 were verified by obtaining the displacement and stress solutions that met the requirements of the Coulomb law and isotropic friction. However, as no other existing methodology offering solutions of this type has been found in the literature, it has not been possible to verify them through result comparisons, and therefore, they have not yet been included in the tool. Once the methodologies proposed in Chapters 6 and 7 are generalized to cover valid complex singular solutions, and it becomes possible to verify the results through comparison of singular solutions with friction contact and λ∈Cobtained via the generalization presented in Chapter 4, the methodology proposed in Chapter 4 could be added to the web application. 11.4.2 File loading capability Amplifying user functionality by enabling direct import of .JSON file into the application. This feature enables users to make modifications to previously generated .JSON file from similar cases, eliminating the need to input parameters manually. Furthermore, the direct import functionality allows users to explore solutions for the same problem with variations in singularity exponent interval or higher resolutions in stress and displacement parameters. This flexibility enhances the versatility of the application, addressing a broader range of user needs. 11.4.3 Parameterisation of variables Incorporating advanced parameterisation capabilities to explore the effect of varying some variables defining the corner problem. These ongoing improvements and planned enhancements underline SingSol’s commitment to providing a robust and user-friendly platform for the analysis of singular stress problems, addressing current challenges and anticipating future needs in the field of structural mechanics and material science. 11.5 Significance of considering frictional contact 11.5.1 Orthotropic and anisotropic bi-material The objective of the presented study of 3P ENF test was to analyze the importance of considering friction in cases involving contact between two materials. By validating this premise using the simple case of an isotropic bi-material, there was no immediate need to incorporate orthotropic, monoclinic, or generally anisotropic materials into the study. However, in scenarios involving such materials, it becomes necessary to account for mixed-mode fracture behavior between mode II and mode III. As demonstrated in Chapter 6, the plane strain and antiplane strain solutions cannot be decoupled in these cases. Consequently, there is a need to generalize the finite element model to accommodate generalized plane strain deformation. 248 M.A. Herrera-Garrido This broader approach could also facilitate understanding how the singular solutions obtained in Chapter 6, featuring perpendicular shear vectors, can be combined so that the linear combination fulfills the Coulomb friction mode and co-linearity of relative tangential displacements vector and frictional shear vector in the contact zone adjacent to the crack tip. 11.5.2 Dundurs mismatch parameter In the present study of 3P ENF test, it was observed that the critical load has a nearly linear relationship with the friction coefficient, and the slope of this relationship was closely related to the discrepancy between the elastic moduli of the materials forming the bi-material. A similar relationship exists between the singularity exponent and the friction coefficient, but in this case, it is known that the slope is approximately β/π, where βis the second Dundurs mismatch parameter of elastic bi-materials. At this point, it is raised whether the relationship between the friction coefficient and the critical load is also somehow influenced by the parameter β. APPENDIX A Appendix A Example of a boundary condition Taking as example the case of only urrestricted for a corner face at ϑ, where ϑ=ϑ0or ϑW, the Cartesian components of the vectors in (2.40) read sr(ϑ) =  −cos(ϑ) −sin(ϑ) 0 ,s3=  0 0 1 ,n(ϑ) =  −sin(ϑ) cos(ϑ) 0 .(A.1) For this case (only urrestricted) Du(ϑ) =   sT r(ϑ) 01×3 01×3  ,and Dϕ(ϑ) =   01×3 nT(ϑ) s3  .(A.2) Now, by substituting (A.1) and (A.2) in expression (2.42):  −cos(ϑ)−sin(ϑ)0 0 0 0 0 0 0   u1(r,ϑ) u2(r,ϑ) u3(r,ϑ) +  0 0 0 −sin(ϑ)cos(ϑ)0 0 0 1   ϕ1(r,ϑ) ϕ2(r,ϑ) ϕ3(r,ϑ) =0, (A.3) we obtain −cos(ϑ)u1(r,ϑ)−sin(ϑ)u2(r,ϑ) = 0, −sin(ϑ)ϕ1(r,ϑ)+ cos(ϑ)ϕ2(r,ϑ) = 0, ϕ3(r,ϑ) = 0. (A.4) By writing the left hand side of these equations in polar coordinates we get ur(r,ϑ) = 0, ϕθ(r,θ) = 0⇒σθθ (r,ϑ) = 0, ϕ3(r,θ) = 0⇒σθ3(r,θ) = 0. (A.5) These results match with the expected boundary condition for only urrestricted conditions. 249 250 M.A. Herrera-Garrido To see how expression (2.47) works, we can analyze it using the above example. First, we assemble the matrix DBC(ϑ), DBC(ϑ) =            sT r(ϑ) 01×3 01×3    01×3 nT(ϑ) sT 3     01×3 nT(ϑ) sT 3    sT r(ϑ) 01×3 01×3            =           −cos(ϑ)−sin(ϑ)0 0 0 0 0 0 0   0 0 0 −sin(ϑ)cos(ϑ)0 0 0 1    0 0 0 −sin(ϑ)cos(ϑ)0 0 0 1  −cos(ϑ)−sin(ϑ)0 0 0 0 0 0 0           , (A.6) then we substitute it in (2.47) leading to          −cos(ϑ)−sin(ϑ)0 0 0 0 0 0 0 −sin(ϑ)cos(ϑ)0 0 0 0 0 0 1 0 0 0 −cos(ϑ)−sin(ϑ)0 −sin(ϑ)cos(ϑ)0 0 0 0 0 0 1 0 0 0                   u1(r,ϑ) u2(r,ϑ) u3(r,ϑ) ϕ1(r,ϑ) ϕ2(r,ϑ) ϕ3(r,ϑ)          =          −cos(ϑ)u1(r,ϑ)−sin(ϑ)u2(r,ϑ) −sin(ϑ)ϕ1(r,ϑ) +cos(ϑ)ϕ2(r,ϑ) ϕ3(r,ϑ) −cos(ϑ)ϕ1(r,ϑ)−sin(ϑ)ϕ2(r,ϑ) −sin(ϑ)u1(r,ϑ) +cos(ϑ)u2(r,ϑ) u3(r,ϑ)          (A.7) and with (A.4)          −cos(ϑ)u1(r,ϑ)−sin(ϑ)u2(r,ϑ) −sin(ϑ)ϕ1(r,ϑ)+cos(ϑ)ϕ2(r,ϑ) ϕ3(r,ϑ) −cos(ϑ)ϕ1(r,ϑ)−sin(ϑ)ϕ2(r,ϑ) −sin(ϑ)u1(r,ϑ)+cos(ϑ)u2(r,ϑ u3(r,ϑ)          =          0 0 0 −cos(ϑ)ϕ1(r,ϑ)−sin(ϑ)ϕ2(r,ϑ) −sin(ϑ)u1(r,ϑ)+cos(ϑ)u2(r,ϑ) u3(r,ϑ)          =wP(r,ϑ) wU(r,ϑ), (A.8) we get the separation of the components of w(r,ϑ)into the prescribed (ur,ϕθ,ϕ3)and unknown (ϕr,uθ,ur)components. B Example of an interface condition The most difficult case of interface conditions dealt in this work is the case of frictionless sliding interface. Considering the same orthonormal basis as in (A.1), the matrices in Table 2.2 for this interface condition take the form D1(ϑ) = 1 √2          sin(ϑ)−cos(ϑ)0 0 0 0 0 0 0 −cos(ϑ)−sin(ϑ)0 0 0 0 −1 0 0 0 0 0 0 −1 0 0 0 0 0 0 −1 0 0 0 0 0 1          ,(A.9) D2(ϑ) = 1 √2          −sin(ϑ)cos(ϑ)0 0 0 0 0 0 0 −cos(ϑ)−sin(ϑ)0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1          ,(A.10) CHAPTER 11 251 ˜ D1(ϑ) = 1 √2          −sin(ϑ)cos(ϑ)0 0 0 0 0 0 0 −sin(ϑ)cos(ϑ)0 −√2cos(ϑ)−√2sin(ϑ)0 0 0 0 0 0 0 0 0 0 0 0 √2 0 0 0 0 0 0 0 0 0          ,(A.11) ˜ D2(ϑ) = 1 √2          −sin(ϑ)cos(ϑ)0 0 0 0 0 0 0 −sin(ϑ)cos(ϑ)0 0 0 0 0 0 0 −√2cos(ϑ)−√2sin(ϑ)0 0 0 0 0 0 0 0 0 0 0 0 √2 0 0 0          .(A.12) By substituting now (A.9) and (A.10) into (2.49) we obtain 1 √2          sin(ϑ)−cos(ϑ)0 0 0 0 0 0 0 −cos(ϑ)−sin(ϑ)0 0 0 0 −1 0 0 0 0 0 0 −1 0 0 0 0 0 0 −1 0 0 0 0 0 1                   uw 1(r,ϑ) uw 2(r,ϑ) uw 3(r,ϑ) ϕw 1(r,ϑ) ϕw 2(r,ϑ) ϕw 3(r,ϑ)          + 1 √2          −sin(ϑ)cos(ϑ)0 0 0 0 0 0 0 −cos(ϑ)−sin(ϑ)0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1                   uw+1 1(r,ϑ) uw+1 2(r,ϑ) uw+1 3(r,ϑ) ϕw+1 1(r,ϑ) ϕw+1 2(r,ϑ) ϕw+1 3(r,ϑ)          =0, (A.13) where superindices wand w+1 refer to two consecutive wedges, and ϑ=ϑw. By multiplying vectors and matrices in (A.13) we get the interface conditions in explicit form 1 √2(sin(ϑ)uw 1(r,ϑ)−cos(ϑ)uw 2(r,ϑ)−sin(ϑ)uw+1 1(r,ϑ)+cos(ϑ)uw+1 2(r,ϑ)) = 0, 1 √2(−cos(ϑ)ϕw 1(r,ϑ)−sin(ϑ)ϕw 2(r,ϑ)−cos(ϑ)ϕw+1 1(r,ϑ)−sin(ϑ)ϕw+1 2(r,ϑ)) = 0, 1 √2(−ϕw 1(r,ϑ)+ϕw+1 1(r,ϑ)) = 0, 1 √2(−ϕw 2(r,ϑ)+ϕw+1 2(r,ϑ)) = 0, 1 √2(−ϕw 3(r,ϑ)+ϕw+1 3(r,ϑ)) = 0, 1 √2(ϕw 3(r,ϑ)+ϕw+1 3(r,ϑ)) = 0. (A.14) By rewriting these conditions in polar coordinates it is easy to see what these equations imply uw θ(r,ϑ) = uw+1 θ(r,ϑ),(A.15) 258 M.A. Herrera-Garrido and in matrix form nT(ϑw)01×3nT(ϑw)01×3    uw(r,ϑw) φw(r,ϑw) uw+1(r,ϑw) φw+1(r,ϑw)    =wcorner U1.(A.56) It is known that the stress field must fulfil the Coulomb law, but the value of each component and, therefore, the value of the stress function vector is unknown. To determine the value of the stress function vector it is necessary to know its value corresponding to the nµ(ϑw) direction. As the stress function value is the same in the wedges wand w+1, expression (A.54)2is rewritten as h01×3nT µ(ϑw)01×3nT µ(ϑw)i    uw(r,ϑw) φw(r,ϑw) uw+1(r,ϑw) φw+1(r,ϑw)    =wcorner U2.(A.57) As both the difference and the sum of the tangential resultant displacement in the k(ϑw,ωk)direction are unknowns, expressions (A.54)3and (A.54)4respectively, the value of the displacement in the direction k(ϑw,ωk)for wedge wis kT(ϑw,ωk)uw(r,ϑw) = wcorner U3,(A.58) and in matrix form kT(ϑw,ωk)01×301×301×3    uw(r,ϑw) φw(r,ϑw) uw+1(r,ϑw) φw+1(r,ϑw)    =wcorner U3,(A.59) and similar for wedge w+1, 01×301×3kT(ϑw,ωk)01×3    uw(r,ϑw) φw(r,ϑw) uw+1(r,ϑw) φw+1(r,ϑw)    =wcorner U4.(A.60) The last expression defines the displacement vector as the sum of the displacements in the m(r,ϑw)direction (A.54)5 m(ϑw,ωk)01×3m(ϑw,ωk)01×3    uw(r,ϑw) φw(r,ϑw) uw+1(r,ϑw) φw+1(r,ϑw)    =wcorner U5,(A.61) By collecting the previous expressions of unknown variables in (A.54) in matrix form, the matrices ˜ D1and ˜ D2in Table 3.2 are constructed by means of the expressions (A.56), (A.57) and (A.59)-(A.61), and once again adding a factor 1 √2or √2 to ensure their orthonormality. Bibliography [1] Arias, R., R. Madariaga, and M. Adda-Bedia (2011). Singular elasto-static field near a fault king. Pure Appl Geophys 168, 2167–2179. [2] Audoly, B. (2000). Asymptotic study of the interfacial crack with friction. J. Mech. Phys. Solids 48, 1851–1964. [3] Banks-Sills, L. (2018). Interface Fracture and Delaminations in Composite Materials. Cham, Switzerland: Springer. [4] Barber, J. R. (2010). Elasticity. Dordrecht: Springer. [5] Barber, J. R. (2018). Contact mechanics. Springer, Cham. [6] Barenblatt, G. (1962). The mathematical theory of equilibrium cracks in brittle fracture. Volume 7 of Advances in Applied Mechanics, pp. 55–129. Elsevier. [7] Barnet, D. and H. Kirchner (1997). A proof of the equivalence of the Stroh and Lekhnitskii sextic equations for plane anisotropic elastottatics. Philos Mag A 76(1), 231–239. [8] Barroso, A. (2007). Caracterizaci´ on de estados singulares de tensi´ on en esquinas multiamteriales. Aplicaci´ on a uniones adhesivas con materiales compuestos. Ph. D. thesis, Departamento de Mec´ anica de Medios Continuos, Teor´ ıa de Estructuras e Ingenier´ ıa del Terreno. Universidad de Sevilla. [9] Barroso, A., E. Graciani, V. Mantiˇ c, and F. Par´ ıs (2012). A least squares procedure for the evaluation of multiple generalized stress intensity factors at 2D multimaterial corners by BEM. Eng Anal Bound Elem 36, 458–470. [10] Barroso, A., B. Lauke, V. Mantiˇ c, and F. Par´ ıs (2016). Tensile and shear strength of bimaterial interfaces within composite materials. Compos Sci Technol 124, 81–88. [11] Barroso, A., V. Mantiˇ c, and F. Par´ ıs (2003). Singularity analysis of anisotropic multimaterial corners. Int J Fract 119(1), 1–23. [12] Barroso, A., V. Mantiˇ c, and F. Par´ ıs (2009). Computing stress singularities in transversely isotropic multimaterial corners by means of explicit expressions of the orthonormalized Stroh-eigenvectors. Eng Fract Mech 79, 250–268. [13] Barroso, A., J. Mar´ ın, V. Mantiˇ c, and F. Paris (2022). Removing stress singularities in standard tests of composite materials. Revista de Materiales Compuestos. [14] Bl´ azquez, A., V. Mantiˇ c, and F. Par´ ıs (2006). Application of BEM to generalized plane problems for anisotropic elastic materials in presence of contact. Eng. Anal. Bound. Elem 30, 489–502. 259 260 M.A. Herrera-Garrido [15] Bogy, D. B. (1971a). On the plane elastostatic problem of a loaded crack terminating at a material interface. J Appl Mech 38(4), 911–918. [16] Bogy, D. B. (1971b). Two edge-bonded elastic wedges of different materials and wedge angles under surface tractions. J Appl Mech 32(2), 377–386. [17] Bogy, D. B. (1972). The plane solution of anisotropic elastic wedges under normal and shear loading. Transactions of ASME, Journal of Applied 39, 1103–1109. [18] Bogy, D. B. and K. C. Wang (1971). Stress singularities at interface corners in bonded dissimilar isotropic elastic materials. International Journal of Solids and Structures 7(8), 993–1005. [19] Bonnet, G. (2009). Orthotropic elastic media having a closed form expression of the Green tensor. Int J Solids Struct 46(5), 1240–1250. [20] Carlsson, L. A., J. W. Gillespie, and R. B. Pipes (1985). On the analysis and design of the end notched flexure (ENF) specimen for mode II testing. J Compos Mater 20, 594–604. [21] Carpinteri, A., P. Cornetti, N. Pugno, A. Sapora, and D. Taylor (2008). A finite fracture mechanics approach to structures with sharp v-notches. Eng Fract Mech 75, 1736–1752. [22] Carpinteri, A. and M. Paggi (2005). On the asymptotic stress field in angularly nonhomogeneous materials. Int J Fract 135, 267–283. [23] Carpinteri, A. and M. Paggi (2007). Analytical study of the singularities arising at multi-material interfaces in 2D linear elastic problems. Eng Fract Mech 74, 59–74. [24] Carpinteri, A. and M. Paggi (2009). Asymptotic analysis in linear elasticity: From the pioneering studies by Wieghardt and Irwin until today. Eng Fract Mech 76(12), 1771– 1784. [25] Chen, H. P. (1998). Stress singularities in anisotropic multimaterial wedges and junctions. Int. J. Solids Struct. 35(11), 1057–1073. [26] Chen, H. P. (2008). Singular stresses at the tip of a crack terminating at frictional interface of fibre reinforced composites. 16th. UK Conference on Computational Mechanics, pp. 1–4. [27] Chen, H. P. (2013). Frictional interface crack-tip singular stress field in anisotropic composites. XII International Conference on Computational Plasticity. Fundamentals and Applications, pp. 1087–1096. [28] Chen, H. P., Z. Guo, and X. Zhou (2012). Stress singularities of contact problems with a frictional interface in anisotropic bimaterials. Fatigue Fract. Eng. Mater. Struct. 35, 718–731. [29] Churchman, C. and D. Hills (2006). General results for complete contacts subject to oscillatory shear. J Mech Phys Solids 54, 1186–1205. [30] Comninou, M. (1976, 07). Stress singularity at a sharp edge in contact problems with friction. J Appl Math Phys 27, 493–499. CHAPTER Bibliography 261 [31] Comninou, M. (1977a). The interface crack. J Appl Mech 44, 631–636. [32] Comninou, M. (1977b). Interface crack with friction in the contact zone. J Appl Mech 44(4), 780–781. [33] Comninou, M. and J. Dundurs (1979). A closed crack tip terminating at an interface. J Appl Mech 46, 97–100. [34] Costabel, M. and M. Dauge (1993). Construction of corner singularities for Agmon- Douglis-Nirenberg elliptic systems. Mathematische Nachrichten 162(1), 209–237. [35] Costabel, M., M. Dauge, and Y. Lafranche (2001). Fast semi-analytic computation of elastic edge singularities. Comput. Meth. Appl. Mech. Eng. 190(24), 2111–2134. [36] D., T. (2007). The Theory of Critical Distance: A New Perspective in Fracture Mechanics. Oxford: Elsevier. [37] Dempsey, J. (1995). Power-logarithmic stress singularities at bi-material corners and interface cracks. J. Adhes. Sci. Technol. 9, 253–265. [38] Dempsey, J. P. and G. B. Sinclair (1979). On the stress singularities in the plane elasticity of the composite wedge. J Elast 9(4), 317–327. [39] Dempsey, J. P. and G. B. Sinclair (1981). On the singular behavior at the vertex of a bi-material wedge. J Elast 11(3), 373–391. [40] Deng, X. (1993). On stationary and moving interface cracks with frictionless contact in anisotropic bimaterials. Proc., Math. Phys. Sci. 443((1919)), 563–572. [41] Deng, X. (1994). An asymptotic analysis of stationary and moving cracks with frictional contact along bimaterial interfaces and in homogeneous solids. Int. J. Solids Structures 31((17)), 2407–2429. [42] Doitrand, A. and D. Leguillon (2018). 3D applicantion of the coupled criterion to crack initiation prediction in epoxy/aluminum specimens under four point bending. Int J Solids Strct 143, 175–182. [43] Dugdale, D. (1960). Yielding of steel sheets containing slits. J Mech Phys Solids 8(2), 100–104. [44] Dundurs, J. (1969). Elastic interaction of dislocations with inhomogeneities. In T. Mura (Ed.), Mathematical Theory of Dislocations, pp. 70–115. ASME, New York. [45] Dundurs, J. (1970). Some properties of elastic stresses in a composite. In A. Eringen (Ed.), Recent Advances in Engineering, Volume 5, pp. 203–216. Gordon and Breach, New York. [46] England, A. (1965). A crack between dissimilar media. J. Appl Mech 32, 400–402. [47] Gao, H. (1992). A closed interface crack in anisotropic bimaterials. Int. J. Fract. 55, 33–39. [48] Gao, H., M. Abbudi, and D. M. Barnett (1992). Interfacial crack-tip field in anisotropic elastic solids. J. Mech. Phys. Solids 40(2), 393–416. 262 M.A. Herrera-Garrido [49] Garc´ ıa, I., M. Paggi, and V. Mantiˇ c (2014). Fiber-size effects on the onset of fibermatrix debonding under transverse tension: A comparison between cohesive zone and finite fracture mechanics models. Engineering Fracture Mechanics 115, 96–110. [50] Gdoutos, E. (2005). Fracture Mechanics. Dordretch: Springer. [51] Gdoutos, E. and P. Theocaris (1975). Stress concentrations at the apex of a plane indenter acting on an elstic half plane. J Appl Mech 42(2), 688–692. [52] Gharpuray, V. M., J. Dundurs, and L. M. Keer (1991). A crack termianting at a slipping interface between two materials. J Appl Mech 58, 960–963. [53] Giner, E., N. Sakumar, F. Fuenmayor, and A. Vercher (2008). Singularity enrichment for complete sliding contact using the partition of unity finite element method. Int. J. Numer. Meth. Engng 76, 1402–1418. [54] Golub, G. and C. Loan (1985). Matrix Computations. Baltimore: The Johns Hopkins University Press. [55] Graciani, E. (2006). Formulation and implementation of the boundary element method for axisymmetric problems with contact. Application to the characterization of the fibermatrix interface in composite materials. (In Spanish). Ph. D. thesis, Universidad de Sevilla. [56] Graciani, E., V. Mantiˇ c, and F. Par´ ıs (2007). On the estimation of the first interpenetration point in the open model of interface cracks. Int. J. Fract. 143(3), 287–290. [57] Griffith, A. (1920). The phenomena of rupture and flow in solids. Philos. Trans. R. Soc. A 221, 163–198. [58] Head, A. K. (1979). The Galois unsolvability of the sextic equation of anisotropic elasticity. J Elast 9, 9–20. [59] Herrera-Garrido, M., V. Mantiˇ c, and A. Barroso (2021). A powerful matrix formalism for stress singularities in anisotropic multimaterial corners. homogeneous (orthogonal) boundary and interface conditions. Theor Appl Fract Mech 119, 103271. [60] Herrera-Garrido, M., V. Mantiˇ c, and A. Barroso (2023). A semi-analytical matrix formalism for stress singularities in anisotropic multi-material corners with frictional boundary and interface conditions. Theor. Appl. Fract. Mech. 127, 104160. [61] Hills, D. A. and J. R. Barber (1993). Interface cracks. Int. J. Mech. Sci. 35, 27–37. [62] Hills, D. A., P. A. Kelly, D. N. Dai, and A. M. Korsunsky (1996). Interface cracks. In Solution of Crack Problems: The Distributed Dislocation Technique, pp. 107–136. Dordrecht: Springer Netherlands. [63] Hong, K., M. Thouless, W. Lu, and J. Barber (2018). Asymptotic stress fields for complete contact between dissimilar materials. J Appl Mech 86(011009). [64] Hwu, C. (2010). Anisotropic Elastic Plates. New York: Springer. [65] Hwu, C. (2021). Anisotropic Elasticity with Matlab. Cham: Springer. CHAPTER Bibliography 263 [66] Hwu, C., M. Omiya, and K. Kishimoto (2003). A key matrix N for the stress singularity of the anisotropic elastic composite wedges. JSME Int. J. Series A 46, 40–50. [67] Ingebrigtsen, K. A. and A. Tonning (1969). Elastic surface waves in crystal. Phys. Rev. 184, 942–951. [68] Inglis, C. (1913). Stresses in a plate due to the presence of cracks and sharp corners. Trans. R. Inst. Nav. Archit. 60, 219. [69] J. Dundurs, M. C. (1979). Some consequences of the inequality conditions in contact and crack problems. J Elast 9(1), 71–82. [70] J. Dundurs, M. L. (1972). Stress concentration at a sharp edge in contact problems. J Elast 2(2), 109–112. [71] Jimenez-Alfaro, S. and D. Leguillon (2023). Finite fracture mechanics at the microscale. Application to bending test of micro cantilever beams. Eng Fract Mech 258, 108012. [72] Johnson, T. and W. Tucker (2009). Enclosing all zeros of an analytic function — a rigorous approach. Journal of Computational and Applied Mathematics 228(1), 418–423. [73] Johnson, W. (1988). Adhesively Bonded Joints: Testing, Analysis, and Design. ASTM International. [74] Karuppanan, S., C. Churchman, D. Hills, and E. Giner (2008). A sliding frictional contact between a square block and an elastically similar half-plane. Eur. J. Mech. A/Solids 27, 443–459. [75] Kirchner, H. O. K. (1989). Elastically anisotropic angularly inhomogeneous media I. A new formalism. Philos Mag B 60(3), 423–432. [76] Knees, D. (2005). Regularity results for quasilinear elliptic systems of power-law growth in nonsmooth domains- Boundary, transmission and crack problems-. Ph. D. thesis, Fakult¨ at Mathematik and Physik der Universit¨ at Stuttgart. [77] Kondratiev, V. A. (1967). Boundary-value problems for elliptic equations in domains with conical or angular points. Trans. Moscow Math. Soc. 16, 227–313. [78] Lauke, B. and A. Barroso (2011). Notched-butt test for the determination of adhesion strength at bimaterial interfaces. Compos Interfaces 18, 661–669. [79] Laursen, T. A. (2003). Computational Contact and Impact Mechanics: Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis. London: Berlin/Heidelberg, Germany. [80] Leblond, J. B. and J. Frelat (2004). Crack kinking from an initially closed, ordinary or interface crack, in the presence of friction. Eng. Fract. Mech. 71, 289–307. [81] Lee, J. and J. Barber (2006). An automated procedure for determining asymptotic elastic stress fields at singular points. J. Strain Anal. 41, 287–295. [82] Lee, J. and H. Gao (1994). A generalized Comninou contact model for interface cracks in anisotropic elastic solids. Int. J. Fract. 67, 53–68. 264 M.A. Herrera-Garrido [83] Lee, J. H., Y. Gao, A. F. Bower, H. Xu, and G. M. Pharr (2018). Stiffness of frictional contact of dissimilar elastic solids. J. Mech. and Phys. of Solids 112, 318–333. [84] Leguillon, D. (1999). Interface crack tip singularity with contact and friction. C.R. Acad. Sci. Paris 327, 437–442. [85] Leguillon, D. (2002). Strength or toughness? a criterion for crack onset at notch. Eur J Mech A Solids 21, 61–72. [86] Leguillon, D., J. Laurencin, and M. Dupeux (2003). Failure initiation in an epoxy joint betweeen two steel plates. Eur J Mech A Solids 22, 509–524. [87] Leguillon, D., E. Martin, O. , and R. Bermejo (2015). Application of the coupled stress-energy criterion to predict the fracture behaviour of layered ceramics designed with internal compressive stresses. Eur J Mech A Solids 54, 94–104. [88] Leguillon, D. and E. Sanchez-Palencia (1987). Computation of singular solutions in elliptic problems and elasticity. Paris: Masson. [89] Leguillon, D. and Z. Yosibash (2003). Crack onset at a v-notch. influence of the notch tip radius. Int J Fract 122, 1–21. [90] Lekhitskii, S. (1938). Some cases of the elastic equilibrium of a homogeneous cylinder with arbitrary anisotropy (in Russian). Appl. Math. Mech. (Prikladnaya Matematrikai Mekhanika) 2, 345–367. [91] Linkov, A. M. and V. F. Koshelev (2006). Multi-wedge points and multi-wedge elements in computational mechanics: evaluation of exponents and angular distribution. Int J Solids Struct 43, 5909–5930. [92] Liu, Y. J., S. Mukherjee, N. Nishimura, M. Schanz, W. Ye, , A. Sutradhar, E. Pan, A. Dumont, A. Frangi, and A. Saez (2011). Recent advances and emerging applications of the boundary element method. Appl Mech Rev 64, 030802. [93] Lothe, J. and D. M. Barnett (1976). On the existence of surface wave solutions for anisotropic half-spaces with free surface. Int. J. Fract. 47, 428–433. [94] Magnier, V., S. G´ ery, and D. Suzanne (2010). Asymptotic approach to analyze singular stress state in anisotropic multi-material: Application to the rivets. Int J Solids Struct 47(16), 2070–2080. [95] Maiti, S. K. (1992). A multicorner variable order singularity triangle to model neighbouring singularities. Int J Numer Meth Eng 35(2), 391–408. [96] Mantiˇ c, V., A. Barroso, and F. Par´ ıs (2014). Singular elastic solutions in anisotropic multimaterial corners. Applications to composites. In V. Mantiˇ c (Ed.), Mathematical Methods and Models in Composites, pp. 425–495. Imperial College Press. [97] Mantiˇ c, V., A. Barroso, and F. Par´ ıs (2023). Computational procedure for singularity analysis of anisotropic elastic multimaterial corners-Applications to composites and their joints. In V. Mantiˇ c (Ed.), Mathematical Methods and Models in Composites (Second Edition), pp. 613–695. World Scientific Publishing Company. CHAPTER Bibliography 265 [98] Mantiˇ c, V., A. Bl´ azquez, E. Correa, and F. Par´ ıs (2006). Analysis of interface cracks with contact in composites by 2D BEM. In M. A. M. Guagliano (Ed.), Fracture and Damage of Composites, pp. 189–248. Southampton: WIT Press. [99] Mantiˇ c, V., F. Par´ ıs, and J. Ca˜ nas (1997). Stress singulrities in 2D orthotropic corners. Int J Fract 83, 67–90. [100] MathWorks (2020). MATLAB User’s Guide. Natick, MA: MathWorks. [101] Moler, C. (2004). Numerical Computing with MATLAB. Society for Industrial and Applied Mathematics. [102] Mr´ oz, Z. and S. Stupkiewicz (1994). An anisotropic friction and wear model. Int J Solids Struct 31(8), 1113–1131. [103] Mukhopadhyay, N., S. Maiti, and A. Kakodkar (2000). Variable singularity boundary element and its applications in computation of SIFs. Comput Struct 77(2), 141–154. [104] Muller, D. E. (1956). A method for solving algebraic equations using an automatic computer. Mathematical Tables and Other Aids to Computation 10(56), 208–215. [105] Ni, L. and S. Nemat-Nasser (1991). Interface cracks in anisotropic dissimilar materials: an analytical solution. J. Mech. Phys. Solids 39, 113–144. [106] Ni, L. and S. Nemat-Nasser (1992). Interface cracks in anisotropic dissimilar materials: General case. Q Appl Math L(2), 305–322. [107] Ortiz, J. E., V. Mantiˇ c, and F. Paris (2006). A domain-independent integral for computation of stress intensity factors along three-dimensional crack fronts and edges by BEM. Int J Solids Struct 43, 5593–5512. [108] Papadakis, P. and I. Babuˇ ska (1995). A numerical procedure for the determination of certain quantities related to the stress intensity factors in two-dimensional elasticity. Comput. Methods Appl. Mech. Eng. 122, 69–92. [109] Picard, D., D. Leguillon, and C. Putot (2006). A method to estimate the influence of the notch-root radius on the fracture toughness measurement of ceramics. J Eur Ceram Soc 26, 1421–1427. [110] Picu, C. R. and V. Gupta (1996). Stress singularities at triple junctions with freely sliding grains. Int J Solids Struct 33(11), 1535–1541. [111] Poonsawat, P., A. Wijeyewickrema, and P. Karasudhi (2001). Singular stress fields of angle-ply and monoclinic bimaterial wedges. Int J Solids Struct 38, 91–113. [112] Popov, V. L. (2017). Contact Mechanics and Friction. Physical Principles and Applications. 2nd Edition. Berlin: Springer. [113] Qu, J. and J. Bassani (1898). Crack on bimaterial and bicrystal interfaces. J Mech Phys Solids 37, 417–433. [114] Rice, J. R. (1988). Elastic fracture concepts for interfacial cracks. J. Appl Mech 55, 98–103. 266 M.A. Herrera-Garrido [115] Rodriguez-Tembleque, L., R. Abascal, and M. H. F. Aliabadi (2012). Anisotropic wear framework for 3D contact and rolling problems. Int J Solids Struct s 241-244, 1–19. [116] Sadd, M. H. (2020). Elasticity. Theory, Applications, and Numerics. 4th Edition. Oxford: Elsevier. [117] S´ anchez-Carmona, S., A. Barroso, E. Correa, and M. Velasco (2023). Numerical and experimental analysis of ultra-thin plies in carbon/epoxy cross-ply laminates subjected to bending after tensile cyclic loading. Theoretical and Applied Fracture Mechanics 128, 104108. [118] Sator, C. and W. Becker (2012). Closed-form solutions for stress singularities at plane bi- and trimaterial junctions. Arch Appl Mech 82(5), 643–658. [119] Sih, G. C., P. Paris, and G. Irwin (1965). On cracks in rectilinearly anisotropic bodies. International Journal of Fracture Mechanics 1, 189–203. [120] Sinclair, G. B. (1999). A note on the removal of further breakdowns in classical solutions of Laplace’s equation on sectorial regions. J Elast 56, 247–252. [121] Sinclair, G. B. (2004). Stress singularities in classical elasticity–I: Removal, interpretation and analysis, –II: Asymptotic identification. Appl Mech Rev 57(4), 251–297 and 385–439. [122] Steigemann, M. (2015). Power-law solutions of anisotropic multi-material elasticity problems. J Elast 118, 63–87. [123] Stroh, A. N. (1958). Dislocations and cracks in anisotropic elasticity. Philos Mag 3, 625–646. [124] Stroh, A. N. (1962). Steady state problems in anisotropic elasticity. J Math Phys 41, 77–103. [125] Sung, J. and W. Chung (2003). Frictional interface crack in anisotropic bimaterial under combined shear and compression. Int J Solids Struct 40, 6839–6857. [126] Sung, K. and S. Purushothaman (1999). Development of conducting adhesive materials for microelectronic applications. J. Electron. Mater. 28(11), 1314–1318. [127] Suo, Z. (1990). Singularities, interfaces and cracks in dissimilar anisotropic media. Proc. R. Soc. Lond, A 427(1873), 331–358. [128] Syst` emes, D. (2023). Abaqus Analysis User’s Manual. Providence, RI: Dassault Syst` emes. [129] Tanuma, K. (1996). Surface-impedance tensor of transversely isotropic elastic materials. Q J Mech Appl Math 49, 29–48. [130] Theocaris, P. and E. Gdoutos (1977). Stress singularities in cracked composite fullplanes. Int J Fract 13(6), 763–773. [131] Ting, T. C. T. (1986). Explicit solution and invariance of the singularities at an interface crack in anisotropic composites. International Journal Solids Struct. 22(9), 965–983. CHAPTER Bibliography 267 [132] Ting, T. C. T. (1990). Interface cracks in anisotropic bimaterials. J. Mech. Phys. Solids 30(4), 505–513. [133] Ting, T. C. T. (1995). Generalized Dundurs constant for anisotropic bimaterials. Int. J. Solids Struct. 32(3/4), 483–500. [134] Ting, T. C. T. (1996). Anisotropic Elasticity, Theory and Applications. New York: Oxford Science Publications. [135] Ting, T. C. T. (1997). Stress singularities at the tip of interfaces in polycrystals. In H. P. Rossmanith (Ed.), Damage and Failure of Interfaces, pp. 75–82. Balkema Publishers, Rotterdam. [136] Ting, T. C. T. and S. Chou (1985). Stress singularities at tip of contact zone in anisotropic interface crack. In Developments in the Science and Technology of Composite Materials, First European Conference on Composite Materials (Second Edition), Bordeaux, pp. 164–169. European Association for Composite Materials. [137] Ting, T. C. T. and S. C. Chou (1981). Edge singularities in anisotropic composites. International Journal of Solids and Structures 17(11), 1057–1068. [138] Ting, T. T. C. and M. Z. Wang (1992). Generalized Stroh formalism for anisotropic elasticity for general boundary conditions. Acta Mechanica Sinica 8, 193–207. [139] Vasilopoulos, D. (1988). On the determination of higher order terms of singular elastic stress fields near corners. Numerische Mathematik 53(1–2), 51–95. [140] Walker, S. and R. Leine (2019). Set-valued anisotropic dry friction laws: formulation, experimental verification and instability phenomenon. Nonlinear Dyn. 96, 885–920. [141] Wang, X. and P. Schiavone (2022). An anisotropic elastic solid with an elliptical inhomogeneity under a non-uniform in-plane and anti-plane remote loading. Q J Mech Appl Math 75(4), 301–313. [142] Westergaard, H. (1933). Stresses at a crack, size of the crack, and the bending of reinforced concrete. ACI J. Proc. 30. [143] Wieghardt, K. (1907). ¨ Uber das Spalten und Zerreissen elastischer K¨ orper (On splitting and cracking of elastic bodies). Zeitschrift f¨ ur Mathematik und Physik 55(1-2), 60– 103. Translation: Rossmanith H.P. Fatigue Fract Engng Mater Struct (1995) 12:1371– 1405. [144] Wijeyewickrema, A. C., J. Dundurs, and L. M. Keer (1995). The Singular Stress Field of a Crack Terminating at a Frictional Interface Between Two Materials. J Appl Mech 62(2), 289–293. [145] Williams, M. L. (1952). Stress singularities resulting from various boundary conditions in angular corners of plates in extension. J Appl Mech 19, 526–528. [146] Williams, M. L. (1957). On the stress distribution at the base of a stationary crack. J. Appl Mech 24(1), 109–114.