Study of the unfolding failure of curved composite laminates.

Universidad de Sevilla Escuela T´ecnica Superior de Ingenier´ıa PhD Thesis: Study of the Unfolding Failure of Curved Composite Laminates A Thesis submitted by Juan Manuel Gonz´ alez Cantero for the Degree of Doctor of Philosophy in the Universidad de Sevilla Supervised by Enrique Graciani D´ ıaz Federico Par´ ıs Carballo 2017 Abstract The use of composite materials in structural applications has been considerably increased in the last years. This increment has introduced the use of composite laminates in complex components including highly curved zones. Highly curved laminates are prone to unfolding failure, a delamination caused by a bending moment which tries to flatten the curvature. Unfolding failure is classically associated to the Interlaminar Normal Stress (INS), which is characterized by the InterLaminar Tensile Strength (ILTS). The ILTS is typically obtained by a four-point bending test applied to an L-shaped specimen with all the plies oriented at 0o. If this test procedure is applied to composite laminates with differently oriented plies an apparent ILTS is obtained. The apparent ILTS exhibits a thicknessdependence which has not been physically explained yet. INS may be obtained by using FEM. Notwithstanding, analytical models are desirable for obtaining lower modelling and calculation times. Nowadays, analytical models have shown high errors in the calculation of the INS for some kinds of geometries and loading states. These errors are due to a free-edge effect associated with the curvature changes, e.g., in the joint of a curved part with a straight arm. Differences in INS values due to a change of curvature may be even of a 100% with respect to the actual value. The PhD project has studied, in a first stage, the four-point bending test, which is the one typically employed for ILTS determination, obtaining a non-linear model for the calculation of the load-displacement and loadbending moment distributions. The model has been successfully correlated with experimental results. The second stage of the PhD has been focused to the INS calculation, developing two different bi-dimensional models based on a series expansion of the displacement and a higher-order moments definition in the stresses. The first model is based in monomials, with very low computational times and a good accuracy, and the second model is based in Legendre polynomials, with still low computational times and a very high accuracy. Both models have allowed the change-of-curvature problem to be solved and to 4 analyse in detail the different parameters involving the stresses distribution near to free-edges. The third stage of the PhD has been focused on the expansion of the bidimensional models to obtain a three-dimensional model in order to consider also the effect of the torsion and the anticlastic effect combined with the curvature, as well as the effect of the finite width of the specimens. This model let us also to calculate accurately the out-of-plane stress state. Finally, in the fourth stage of the PhD, the stress calculation methodologies developed have been used for the data reduction of an existing test campaign. Results suggest that in laminates with a low apparent ILTS, intralaminar failures in the matrix direction may have taken place before delamination. Comparing the load causing the failure with the predictions made for the first intralaminar/interlaminar failure, and observing the crack locations in tested coupons, a feasible explanation for the thickness dependence has been found. This explanation consists in distinguishing two kinds of unfolding failures: first, the traditional unfolding initiated by the INS and, second, the newly defined induced unfolding, which is initiated by a failure associated to intralaminar stresses that propagates as a delamination due to the presence of the high INS. Based in the previous results, some optimization guidelines have been defined for the design and sizing of curved composite laminates based in the unfolding failure. Acknowledgements I wish to thank Federico Par´ıs and Enrique Graciani for their support which has made possible this project. Special thanks and recognition are also due to FIDAMC, the Spanish Foundation for the Development, Application and Research on Composite Materials, whose support has contributed to the development of the project. A special recognition goes to Airbus Group Innovations (AGI) for the opportunity of working in one of its main researching centres in Germany, the Airbus Group Chair of Seville for its support and Airbus Operations S.L. for its collaboration in the project. Finally, I wish to thank my family, my friends and my partners in FIDAMC for their support. 6 Contents List of Figures 11 List of Acronyms 15 1 Introduction 17 1.1 Motivations and aims of the project . . . . . . . . . . . . . . 19 1.2 State of the art in stresses calculation in curved laminates . 22 1.3 State of the art in failure criteria . . . . . . . . . . . . . . . 26 1.3.1 Failure initiation criteria . . . . . . . . . . . . . . . . 27 1.3.2 Failure propagation criteria . . . . . . . . . . . . . . 29 1.4 Thesisoverview......................... 30 2 Analysis of the four-point bending test 33 2.1 Standardized testing procedure for assessment of ILTS . . . 34 2.2 Non-linear model for obtaining the bending moment . . . . . 36 2.3 Simplification of the non-linear model . . . . . . . . . . . . . 48 2.4 Comparison of the non-linear model with experimental results 50 2.5 Comparison of the non-linear model with the ASTM procedure 53 3 Bi-dimensional models for evaluating interlaminar stresses 55 3.1 Theoretical basis of the regularized models . . . . . . . . . . 56 3.2 Approximations of the regularized models . . . . . . . . . . 61 3.3 Theoretical basis of the non-regularized models . . . . . . . 67 3.4 Monomials based model (MBM) . . . . . . . . . . . . . . . . 71 3.4.1 Development of the MBM for straight beams . . . . . 72 3.4.2 Development of the MBM for curved beams . . . . . 82 3.4.3 MBM numerical limitations . . . . . . . . . . . . . . 90 3.5 Legendre polynomials based model (LPBM) . . . . . . . . . 91 3.5.1 Characteristics of the Legendre polynomials . . . . . 92 3.5.2 Development of the LPBM for straight beams . . . . 95 3.5.3 Development of the LPBM for curved beams . . . . . 99 8 CONTENTS 3.6 Resolution of the regularized problem by using the non-regularized models.............................. 104 3.7 Resolution of the homogeneous problem . . . . . . . . . . . 107 3.8 Resolution of the joints between components . . . . . . . . . 109 3.8.1 Continuity of the circumferential stresses between components ......................... 110 3.8.2 Continuity of the shear stresses between components 111 3.8.3 Continuity of the circumferential displacements between components . . . . . . . . . . . . . . . . . . . 112 3.8.4 Continuity of the transverse displacements between components ....................... 114 3.8.5 Application of the boundary conditions and the continuity conditions . . . . . . . . . . . . . . . . . . . . 116 3.9 Results of the non-regularized models and validation by FEM 117 3.9.1 L-shaped beam under bending moment . . . . . . . . 117 3.9.2 L-shaped beam under a compressive load . . . . . . . 119 3.9.3 Joggle under a tensile load . . . . . . . . . . . . . . . 121 3.10 Numerical characteristics of the MBM and LPBM non regularizedmodels.......................... 123 3.11 Regularization distance . . . . . . . . . . . . . . . . . . . . . 125 3.12 Off-topic applications of the non regularized models . . . . . 129 3.12.1 Joint of beams with different material properties . . . 129 3.12.2 Non-regularized effects due to the boundary conditions 131 3.12.3 Non-regularized effects due to punctual loads . . . . . 133 4 Three-dimensional models for evaluating interlaminar stresses137 4.1 3D model with double regularization for singly-curved beams under bending moment and residual strains . . . . . . . . . 139 4.1.1 Theoretical development . . . . . . . . . . . . . . . . 141 4.1.2 Loads application and boundary conditions . . . . . . 149 4.1.3 Particular cases . . . . . . . . . . . . . . . . . . . . . 156 4.2 Analysis of the effect of residual stresses in the unfolding failure161 4.2.1 Homogeneous anisotropic materials . . . . . . . . . . 162 4.2.2 Composite laminates . . . . . . . . . . . . . . . . . . 164 4.2.3 Finite elements comparison . . . . . . . . . . . . . . 166 4.3 Analysis of the three-dimensional effects over the unfolding failure .............................. 166 4.3.1 Homogeneous anisotropic materials . . . . . . . . . . 167 4.3.2 Composite laminates . . . . . . . . . . . . . . . . . . 169 4.3.3 Finite elements comparison . . . . . . . . . . . . . . 170 4.4 Three-dimensional non-regularized models . . . . . . . . . . 173 CONTENTS 9 4.4.1 Development of the 3D non-regularized model in flat laminates ........................ 175 4.4.2 Development of the 3D non-regularized model in curved laminates ........................ 190 4.5 Applications of the 3D non-regularized models . . . . . . . . 209 4.5.1 Through-the-width distributions . . . . . . . . . . . . 209 4.5.2 Free-edge effects . . . . . . . . . . . . . . . . . . . . . 211 5 Failure mechanisms involved in the unfolding failure 215 5.1 The concept of the induced unfolding failure mechanism . . 216 5.2 Experimental results . . . . . . . . . . . . . . . . . . . . . . 218 5.2.1 Traditional unfolding . . . . . . . . . . . . . . . . . . 219 5.2.2 Induced unfolding . . . . . . . . . . . . . . . . . . . . 220 6 Design recommendations 227 6.1 Evaluation of the maximum bending moment with the combination of traditional and induced unfolding . . . . . . . . . 228 6.2 Optimization recommendations based on both failure mechanisms.............................. 234 7 Concluding remarks and future developments 237 A Matrices of the differential equation of the 2D non-regularized models 243 A.1 MBM matrices for the straight beam . . . . . . . . . . . . . 243 A.2 MBM matrices for the curved beam . . . . . . . . . . . . . . 244 A.3 LPBM matrices for the straight beam . . . . . . . . . . . . . 245 A.4 LPBM matrices for the curved beam . . . . . . . . . . . . . 246 B Matrices expressions of the 3D non-regularized models 247 B.1 Stiffness matrices for the flat laminate . . . . . . . . . . . . 247 B.2 Auxiliary matrices for the flat laminate . . . . . . . . . . . . 249 B.3 Stiffness matrices for the curved laminate . . . . . . . . . . . 250 B.4 Auxiliary matrices for the curved laminate . . . . . . . . . . 256 Bibliography 263 Curriculum Vitae 275 List of Acronyms 16 Chapter 1 Introduction Nowadays, composite materials are used in many different sectors and applications, from the aeronautical and automotive sectors to the rail and marine sectors. In their origins, composite materials were used mainly in lightweight structures involving thin laminates where the curvature radius is much higher than the thickness and the interlaminar stresses are typically neglected due to their extremely low values. Notwithstanding, in the last years the use of this kind of materials has been considerably increased, requiring to incorporate the use of composite laminates in more complex components, where thicker laminates and highly-curved laminates are included. The introduction of the composite materials in this kind of components and structures has caused the apparition of new failure mechanisms which were not previously considered. The unfolding failure is found within those failure mechanisms. The unfolding failure consists in a delamination given in highly-curved and moderately-curved oriented-fibre composite laminates which is produced preliminarily by the interlaminar stresses induced by the curvature. This kind of failure is generally associated with an opening bending moment. When a curved component is loaded under a bending moment which tries to open the curvature trying to give a flatter geometry to the curved laminate the plies try to separate like the pages of a book, see Figure 1.1. This effect is due to the different length between the plies located at the lowest and at the highest radius. Interlaminar normal stresses (INS) appear as the opposition of the plies to be separated, and, when they are high enough, the delamination involving the unfolding failure is produced. Some examples of composite laminates including the aforementioned highly-curved parts and, therefore, prone to the unfolding failure, are the L-shaped beams, the T-shaped beams, the Ω-shaped beams, the joggles or the corrugated laminates [1,2]. Introduction 18 Figure 1.1: When the laminate is opened by a bending moment the plies try to separate. The unfolding failure is a failure mechanism, mainly characteristic of the composite laminates. This fact is due to the low strength of these materials in the interlaminar direction (typically in the range of 40-80 MPa in CFRP with epoxy resins) respect to the intralaminar direction (which, for CFRP in the fibre direction, is typically in the range of 2000-3000 MPa). Thus, through-the-thickness failures are considerably more probable in this kind of materials than in isotropic materials in which strength is the same in all directions. The unfolding failure prediction is typically carried out in three different steps. The first step consists in obtaining the loads and boundary conditions of the curved laminate. This step is conditioned by the model used in the second step, which consists in the stress calculation. The stress calculation may be done by a Finite Element (FE) model or by analytical procedures. Thus, if a bi-dimensional analytical model is considered for the stress calculation, the loads obtained in the first step must be compatible with this model. In many cases non-linearities must be considered for the loads calculation, e.g., in the four-point bending test [3]. Finally, once the stresses have been calculated, a failure criterion is applied. The current PhD project is focused mainly in the second step. Different analytical stress calculation methodologies are developed in the context of the interlaminar stresses in curved laminates in order to reduce the computational times (in comparison with the FE models) and to improve the accuracy (in comparison with the nowadays analytical procedures which present high errors). The third step, relative to the failure criteria, is also studied with the detailed analysis of the failure mechanism involving the delamination of the unfolding failure, which is validated with some experimental results. The state of the art for the second and third steps in the unfolding failure prediction is analysed individually in the following sections, after introducing the motivations and aims of the project. The first step is not 19 Introduction commented in those sections while its analysis depends on each kind of problem. Notwithstanding, in the present document the particular case of the four-point bending problem is analysed in further chapters due to its importance for the experimental estimation of the InterLaminar Tensile Strength (ILTS). 1.1 Motivations and aims of the project The Thesis is oriented to analyse two main problems that have been found in the calculation of the unfolding failure reserve factor in typical components used in the aeronautical industry. The first problem is a thickness dependence of the ILTS, which is typically obtained with a four-point bending test over a curved composite laminate. This thickness dependence is widely observed in the unfolding calculations in the aeronautical companies and it is not demonstrated physically. The second problem is the lack of accuracy of the analytical models when applied to some kinds of loading states. These loading states are dominated by non-regularized effects which are not taken into account in the traditional analytical calculations, and nowadays are obtained typically by using a detailed FE model. However, an analytical or semi-analytical procedure is desired to improve the computational times allowing to reduce times in design problems. When a quasi-isotropic composite laminate is tested according to the four-point bending test standardized in the procedure ASTM D 6415/D 6415M [3], a curved beam strength (CBS) is obtained as the bending moment per unit of width applied to the curved zone of the L-sectioned beam at the instant of the failure. Considering this bending moment and by using a certain analytical or numerical calculation method for calculating the stresses (e.g., Ko and Jackson’s equations [4] or Kedward’s formula [1]), a maximum INS may be obtained at the instant of failure. The four-point bending test has neglectable interlaminar shear stresses (ISS) in the curved zone, an then, the traditional delamination failure criteria (e.g., the Kim and Soni criterion [5]), which depends generally on the INS and the ISS, may be applied as a maximum stress failure criterion. Therefore, considering that the failure is due to the INS as has been traditionally done, the maximum INS at the instant of failure may be considered as an apparent ILTS. The apparent ILTS, obtained as commented before from quasi-isotropic laminates, has typically shown a thickness dependence as shown in the example of Figure 1.2, where thinner specimens have lower strengths than thicker ones. This thickness dependence is shown by Edwards and Thompson [6]. In the Figure, the apparent ILTS has been expressed non-dimensionally re- 1.1. Motivations and aims of the project 20 spect to the strength in direction 2 (matrix in-plane direction of a ply), S22. The cause of this thickness dependence has not been physically explained yet although it is widely observed in the aeronautical sector. In this way, nowadays applications obtain the apparent ILTS from the four-point bending test depending on the thickness, and it is used when it is applied to other loading states and other kind of specimens (such as T-shaped beams or joggles). The application of this dependency to this other kind of specimens has not been demonstrated to be valid and it may cause errors in the estimation of the final failure. Hence, it requires a deeper study of the stresses and the failure mechanisms to determine the cause of the thickness dependence. Figure 1.2: Example of thickness dependence of the apparent ILTS obtained from quasi-isotropic laminates. The thickness dependence has not been shown in unidirectional laminates [0]n, as can be seen in the unidirectional specimens of Hoffmann et al. [7]. Furthermore, it is neither observed in some fabric composites, as can be seen in the results of the maximum INS obtained by Avalon and Donaldson [8]. The thickness dependence is against the size effects explained by Wisnom [9]. In accordance with the size effects, thicker laminates should have lower ILTS. Hoffmann et al. [7] tried to demonstrate the size effects in the ILTS. Notwithstanding, the ILTS obtained in the composite laminates was not that expected, obtaining the lowest strength in the thinnest specimens and the highest strength in the medium thickness specimens, which is also against the apparent ILTS distribution of Figure 1.2, where a decrement of the ILTS is obtained also in the medium thickness specimens. Edwards and Thompson [6] attributed the thickness dependence to the fact that, due to manufacturing defects, the curved specimens have small 21 Introduction interlaminar cracks or delaminations prior to the tests, and these delaminations propagate at a lower load than the corresponding to the ILTS value, causing the apparent ILTS to become lower than the real one. Makeev et al. [10] studied the role of the porosity in the delamination, showing that, with a small increment in the porosity and in the void size, the strength may decrease significantly. In this way, other possible explanation for the decrease in the strength in thinner laminates may be given by considering the presence of a more significant void content in those specimens. Fletcher et al. [11,12] and Kim et al. [13] assumed that the delamination is highly dominated by the free-edge effects, using a resin edge to reduce the free-edge effects and to improve the strength. Following their reasoning, other possible explanation consists in the free-edge effect, it requiring a deeper study on these stress concentrators respect to the thickness. Regarding the second motivation of the Thesis, consisting in the lack of accuracy in the stress state obtained by the analytical methodologies for some loading states, it may be due to several causes. Traditional analytical methodologies calculate only the regularized stress state, which consists in considering the curved laminate isolated, and generally loaded only under end loads, although some authors consider also distributed loads. However, there are some effects which may cause perturbations in the stresses obtaining a non-regularized stress state. This is the case, e.g., of the free-edge effects, where the free-edge boundary condition causes a stress concentration in a localized zone near to the edge. Other cause for the appearance of non-regularized stresses is the change of curvature of the laminate. Highly curved laminates are commonly linked with flat laminates. This is the case of the L-shaped beam, where the curved part of the L is joined to two straight arms, or the case of a joggle, where several curvature changes are given along the section. Considering the particular case of a L-shaped beam loaded under an end bending moment, if a traditional stress calculation method is used, such as the Ko and Jackson’s equations [4], the curved part of the section has a constant INS in the circumferential direction, while null INS are given in the straight arms. Therefore, a discontinuity is given at the joint of the curved part with the arms, see Figure 1.3(a), where the solution obtained with Ko and Jackson’s equations is depicted. This discontinuity is given in all the components of stresses and displacements, and not only in the INS. Hence, the stress tensor is modified by non-regularized effects causing a smooth variation of the stresses between the straight arm and the curved part in order to accomplish the equilibrium equations and the compatibility of the displacements, see Figure 1.3(b). This effect has been shown by Most et al. [14] and it is studied in the current Thesis. 1.2. State of the art in stresses calculation in curved laminates 22 (a) (b) Figure 1.3: INS colour plot in a L-shaped beam loaded under end bending moment. (a) Regularized stresses. (b) Non-regularized stresses. The non-regularized effects due to the change of curvature are localized near to the section where it takes place, modifying the stresses respect to the regularized ones only in this zone. Therefore, they are especially important in laminates with a short curved part respect to the thickness, or when the maximum regularized value is obtained near to the change of curvature. The aims of the Thesis are oriented to improve the calculation and prediction of the unfolding failure in curved composite laminates. The first aim of the project is oriented to develop a novel analytical procedure for obtaining the stresses in curved laminates, in order to predict the non-regularized effects without requiring a FE model. The use of an analytical procedure allows the stresses to be obtained with a satisfactory accuracy in smaller computational times than FE, which is especially useful for optimization problems where many evaluations of the stresses have to be done. Second, an explanation for the thickness dependence is looked for, trying to find the real failure mechanism involving the unfolding failure, in order to develop a new, more realistic, failure criterion. Moreover, the first objective helps to the resolution of the second one, as a satisfactory accuracy in the stress calculation is required for understanding the failure mechanism of the component. 1.2 State of the art in stresses calculation in curved laminates The calculation of composite laminates is based on the Classical Laminate Theory (CLT) [15, chap. 4] and on the First-order Shear Deformation Theory (FSDT) [16], which consist in a plate theory where the composite laminate is approximated by homogeneous equivalent material through the stiffness matrix A-B-D, and where the INS is neglected. Therefore, this 23 Introduction theory cannot be applied directly to the highly-curved laminates for obtaining the interlaminar stress. Notwithstanding, the CLT may be applied to curved laminates for obtaining the intralaminar stresses with a good accuracy (see [17], for instance), or for obtaining the deformation of curved laminates under different kind of loads (see [18], for instance). Lekhnitskii et al. [19, chap. 3] obtained the closed-form solution of the through-the-thickness stresses in curved bi-dimensional beams made by homogeneous anisotropic material under end load and bending moment. Hence, it cannot be applied to composite laminates while composite laminates are heterogeneous at a laminate level (non-homogeneities associated with the presence of fibre and matrix inside each ply are not considered in the present document). Notwithstanding, many authors have approximated the laminate by homogeneous equivalent material and have applied the Lekhnitskii’s equations [20–23] to determine the INS. The maximum INS for a curved beam under bending moment obtained by using the Lekhnitskii’s equations has been approximated by a closed form equation given by Kedward et al. [1]. This equation is independent of the stiffness properties of the material, depending only on the bending moment and on the geometrical parameters, and being only accurate when the circumferential stiffness is lesser than 20 times the radial stiffness. The INS caused by an axial force may be neglected while the dominant load for those stresses is the bending moment. With reference to the shear problem, the maximum of the ISS may be approximated by the maximum ISS in a flat laminate, given by Cui et al. [24]. The extension of the Lekhnitskii’s equations to a laminate problem was solved by Ko and Jackson [4], allowing the effect of the stacking sequence over the INS distribution to be obtained. Ko and Jackson’s equations constitute a linear system of equations where the number of equations depends on the number of plies Np, so that it requires to solve 3Npequations for the bending moment and 8Npequations for the end load (which may be decomposed in end axial force and end shear force). Smidt [25] and Shenoi and Wang [26] developed a similar model to the Ko and Jackson’s equations, applied to curved sandwich beams under bending load. Some authors have been working in the last years in the development of simpler models to solve the INS in curved laminates considering the stacking sequence. One of the main objectives of those models is to obtain approximations of the Ko and Jackson’s equations which may be implemented in a FE software to calculate the INS in curved laminates. Other objectives for obtaining approximated models are to calculate problems which cannot be solved with Ko and Jackson’s equations, such as double-curvature problems, three-dimensional cases or distributed loads (e.g., pressure), or to 1.2. State of the art in stresses calculation in curved laminates 24 reduce computational times, e.g., for solving optimization problems. The CLT may be extended for calculating the INS from the in-plane stresses by using the equilibrium equations. Notwithstanding, the CLT applied to curved laminates requires a modification of the A-B-D matrix so that considering the different lengths of the plies in the inner and outer radius. This modification implies logarithmic expressions in the terms of the stiffness matrix which has been presented by several authors as Qatu [27], Kress et al. [28] or Guedes and S´a [29]. Other approximated models begin with the Kirchoff assumption used in the CLT and consider different hypothesis. That is the case of Kress [30] and Kress et al. [28], who used the aforementioned assumption only in the circumferential displacements, calculating the expression of the through-the- thickness displacements. Some authors have considered the hypothesis of thin laminates tR (where tis the thickness and Ris the mean radius) in the development of an approximated model for curved beams, as Rodr´ıguez de la Cruz et al. [31] for sandwich curved beams, Lin and Hsieh [17] for variable curvature beams, Vargas et al. [32] for studying the thermal and mechanical stresses in unsymmetrical curved laminates or Bruno et al. [33] for studying the delamination propagation in curved laminates. Those models have a high accuracy when the tRhypothesis is accomplished. The series expansion in the displacements and/or the stresses is one of the most used approximations for calculating the stresses in thicker curved laminates where t∼R. Carrera et al. [34] presented a compilation of several methods based on series expansions in the displacements, not only for curved laminates. Those models usually start from the Timoshenko theory [35,36] expanding the displacements with additional terms. The most basic series expansion is done by using monomial functions in the thickness, as is the case of Matsunaga [37,38]. Other authors have preferred a Fourier series as Ren [39], who expanded the INS by using a Fourier series of sines, or the use of zig-zag functions as Carrera [40] and Icardi and Ferrero [41]. More advanced models including series expansions are based in more complex numerical techniques as the differential quadrature method [42,43] or the higher-order shear deformation theory [44,45]. Those analytical procedures have been compared with experimental results in several loading states and geometries. For example, Cui and Ruiz [46] and Wisnom et al. [47] tested C-specimens under combined shear, axial and bending loads, Cui et al. [24] tested curved beams under pure bending, and McRobbie et al. [48] tested open-rings of composite laminates and sandwiches under bending due to end loads. Many of the mentioned analytical approximations are oriented to develop 25 Introduction elements for a FE software. Finite elements allows the INS in curved laminates to be accurately obtained, but they require to refine the mesh in the interlaminar direction in a three-dimensional model [49]. The refinement of the mesh in the interlaminar direction requires also refining the mesh in the intralaminar plane in order to maintain the aspect ratio, so the final model has an excessive number of elements consuming too much computational time to obtain the solution. In that way, many authors have worked in the development of new elements for curved composites which do not require so much interlaminar refinement to obtain accurate INS. Rattanawangcharoen [50] developed a three-dimensional element with the main aim of analysing weakly bonded cylindrical panels. Kress and Winkler [51] developed a shell element for corrugated laminates. Cinefra and Carrera [52] and Vidal et al. [53] developed cylindrical shell elements using a layerwise theory. Finally, Fraternali and Bilotti [54] developed non-linear elements for curved laminates. The calculation of the INS may be also influenced by the residual stresses due to the manufacturing process and the cooling of the laminate after the curing process. The different mechanisms producing residual stresses are summarized by Wisnom et al. [55]. The residual stresses have been typically associated to the in-plane stresses and interlaminar residual stresses have been typically neglected, except in the free-edges. However, when the laminate is curved, residual INS appear due, as the residual in-plane stresses, to the different expansion coefficients between the plies. Those residual INS in curved laminates were shown experimentally and numerically by Kim and Lee [56] and by Takagaki et al. [57] in the case of cylindrical laminates, and analytically by Huang and Tauchert [58] in doubly-curved laminates. Notwithstanding, the aforementioned analytical and numerical procedures do not consider other effects which may influence in a high percentage the stress distributions. Those procedures consider typically the curved laminate isolated, calculating only the regularized stresses. However, some factors may cause a non-regularized distribution near to a particular focus. One of those focuses is the change of curvature, as shown by Most et al. [14]. The compatibility and equilibrium in a joint between two laminates with different curvature cause a perturbation over their stress distributions, which is especially important in the case of INS. Therefore, the change of curvature is an important factor to be considered in the unfolding calculation. The change of curvature is related also with the laminates of variable curvature or with the doubly-curved laminates. The double curvature in thick laminates has been widely studied in the literature, where the works of Qatu [59], Asadi et al. [60,61] and Tornabene et al. [62] may be mentioned. The most well-known non-regularized focus is the free-edge of the lami- 1.4. Thesis overview 32 folding. Finally, Chapter 6 summarizes the main results obtained in the previous Chapters in order to develop some design recommendations for curved composite laminates prone to unfolding failure. The PhD project has been developed in the facilities of FIDAMC, the Spanish Foundation for Research, Development and Application of Composite Materials, with the academic supervision of the University of Seville. Chapter 2 Analysis of the four-point bending test The traditional unfolding failure is associated with the interlaminar normal stresses (INS). The unfolding failure is generally caused by a bending moment which tries to open the curvature, causing that the different plies try to separate and inducing interlaminar tensile stresses. Hence, the traditional unfolding failure is mainly characterized by the interlaminar tensile strength (ILTS). There are several procedures to assess the ILTS property, which are summarized by Makeev et al. [10]. The simplest method to obtain this strength consists in a tensile test in the interlaminar direction, which is standardized by the test procedure ASTM D7291 [80]. Several authors, such as Hoffmann et al. [81], have characterized the ILTS using this procedure. Makeev et al. [10] showed that the ILTS obtained using this procedure is highly dominated by the size effects, it being different from the one obtained using different procedures and not being very representative for a curved specimen. The curved beam ILTS is obtained by using a four-point bending test, which is standardized by the test procedure ASTM D6415 [3]. This test procedure is complexer than the tensile test [80] but it is more widely used since it is less affected by the size effects. Notwithstanding, the curved beam is prone to manufacturing defects in form of porosity. Several authors such as Makeev et al. [10] and Jackson and Martin [101] have commented the effect of the porosity over the ILTS. Furthermore, this testing procedure has shown a thickness dependence of the ILTS, which have not been physically explained yet (see Edwards and Thompson [6]). The topic of the thickness dependence is dealt with in Chapter 5. The present Chapter studies the four-point bending test since it is the most common testing procedure for obtaining the ILTS, which is a key 2.1. Standardized testing procedure for assessment of ILTS 34 parameter in the analysis of unfolding failure. After analysing the testing procedure, a non-linear model for obtaining the bending moment in the curved part of the sample from the applied load is developed and finally it is compared with the methodology proposed by the ASTM procedure and compared with some experimental results. 2.1 Standardized testing procedure for assessment of ILTS The most widely used testing procedure for assessment of ILTS of composite materials is the four-point bending test, whose testing configuration is schematically depicted in Figure 2.1. The geometry of the specimen is determined by the mean radius of the curved part R, the thickness tand the angle in the undeformed configuration ϕi. The opening angle of the L-shaped section is obtained as Θi= 180o−2ϕi, so for a typical shape with Θi= 90othe angle ϕihas a value of 45o. Figure 2.1: Four-point bending test set-up. The testing set-up is defined by the radius of the rollers Rr, the distance between the upper rollers Lsand the horizontal distance between one upper roller and the nearest downer roller dx. The distance l0(indicated in Figure 2.1) is calculated from the previous values as follows: 35 Analysis of the four-point bending test l0=dx+ (2Rr+t) sin ϕi cos ϕi .(2.1) The total applied load per unit of width is denominated as Pand the relative displacements of the rollers is denominated as δ. This testing procedure includes the complexity of being highly non-linear due to the change of the position of the points where the rollers contact the sample when the load is increased, so a large displacements model is required. This non-linearity affects to the calculation of the bending moment in the curved part of the sample, which is required for calculating the INS and therefore for obtaining the ILTS from the failure load. The testing procedure ASTM D 6415/D 6415M [3] estimates this nonlinearity by using the hypothesis of considering that the arms of the L- shaped specimens remain straight during the loading process and turn to an angle ϕ(which for a null load has a value of ϕi) which is obtained from the experimental displacement δ. Therefore, the bending moment obtained from the testing procedure at the failure load does not depend only on this failure load, it depending also on the displacement measured at failure. In this way, the procedure does not need the stiffness property of the material as an input while both load and displacements are obtained experimentally. Considering the undeformed shape, the bending moment (per unit of width) is obtained from the applied load (per unit of width) as follows: Mus 0=−P 2 cos ϕi l0=−P 2 cos ϕi dx+ (2Rr+t) sin ϕi cos ϕi .(2.2) Considering that the straight arm remains straight and turn to a new angle ϕ, the corrected bending moment yield: M0=−P 2 cos ϕ dx+ (2Rr+t) sin ϕ cos ϕ,(2.3) where the new angle ϕis obtained from the displacement δas follows: dy=dxtan ϕi+2Rr+t cos ϕi−δ, (2.4a) ϕ= arcsin −dx(2Rr+t) + dypd2 x+d2 y−4R2 r−4Rrt−t2 d2 x+d2 y!.(2.4b) Once the bending moment is obtained the stresses are calculated typically by using the Lekhtniskii’s equations [19, chap. 3] or the Ko and Jackson 2.2. Non-linear model for obtaining the bending moment 36 equations [4]. The calculation of the stresses from the loads are analysed in further chapters. The testing procedure ASTM D 6415/D 6415M [3] recommends the use of the Kedward’s formula [1], an approximation of the Lekhnitskii’s INS maximum value: σr,max =−3M0 t√4R2−t2.(2.5) The approach of considering that the straight arms remain straight may be very conservative when the non-linearity is very pronounced, like in the thinner specimens. Therefore, for improving the calculation of the bending moment, a novel non-linear model is developed in the following sections. 2.2 Non-linear model for obtaining the bending moment In order to improve the calculation of the bending moment in a fourpoint bending test a novel model considering the non-linearities due to the large displacements is developed. The non-linearities are considered by calculating the contact point of the rollers with the specimen in the deformed configuration. The specimen is modelled by using a Timoshenko beam model [35, 36] and considering the symmetry of the problem according to Figure 2.2. Figure 2.2: Beam approximation of the specimen in a four-point bending test. 37 Analysis of the four-point bending test The whole beam is divided in 3 different beams. Beam 1 is defined as the part of the beam located between the two rollers, beam 2 is defined as the part of the beam located between the inner roller and the beginning of the curved part, and beam 3 is defined as the curved part. Furthermore, point 0 is defined as the point of contact with the outer roller, point 1 is defined as the point of contact with the inner roller, point 2 is defined as the point of joint between beam 2 and 3, and point 3 is defined as the symmetry point of the whole beam. Beginning with the analysis of beam 1, the deformed configuration is depicted in Figure 2.3. Figure 2.3: Deformed configuration of beam 1. Let us assume that when the upper roller is displaced a distance δin the vertical direction, the deformed shape has an angular change in the slope of the mean line at point 0 defined by ϕ0and an angular change in the slope at point 1 defined by ϕ1. The initial horizontal distance between the contact points translated to the mean line of the beam is called L0. This horizontal distance in the deformed shape, translated to the undeformed one according to Figure 2.3, is called L. The translation from the undeformed shape to the deformed shape is done by using the displacements at point 0, u0and w0, and the displacements at point 1, u1and w1. In accordance with Timoshenko beam theory and defining the coordinate system (s, z) shown in Figure 2.3, the displacements are approximated considering that plane sections remain plane: 2.2. Non-linear model for obtaining the bending moment 38 u(s, z) = u0(s) + zφ(s), w(s, z) = w0(s).(2.6) In this way, considering the displacements given in equation (2.6) and considering large displacements the angles ϕ0and ϕ1are calculated as: ϕ0= arctan φ0, ϕ1= arctan φ1,(2.7) where φ0=φ(0) and φ1=φ(L), with φ(s) being the slope change of the section. Considering a large displacements model, φ(s) is the tangent corresponding to ϕ(s). By using the definition of the displacements (2.6), the strains can be written as: εs=du0(s) ds +zdφ(s) ds , εz= 0,(2.8a) γsz =φ(s) + dw0(s) ds .(2.8b) In accordance with Figure 2.3, distance Lmay be calculated as: L=dx+ˆ Rr(sin(ϕi−ϕo) + sin(ϕi−ϕ1)) + (w1−wo) sin ϕi+ (uo−u1) cos ϕi, (2.9) where the parameter ˆ Rris defined as: ˆ Rr=Rr+t 2.(2.10) Substituting (2.7) into (2.9): L=dx+ˆ Rr sin ϕi−φocos ϕi p1 + φ2 o +sin ϕi−φ1cos ϕi p1 + φ2 1!+ (w1−wo) sin ϕi+ (uo−u1) cos ϕi,(2.11) The vertical force applied by the rollers is P/2, and then considering that the applied load in the outer roller has an angle ϕi−ϕoand applying the load translated to the undeformed configuration, it holds that: N(s) = −Psin ϕo 2 cos(ϕi−ϕo), Q(s) = −Pcos ϕo 2 cos(ϕi−ϕo),(2.12a) M(s) = −Pcos ϕo 2 cos(ϕi−ϕo)s, (2.12b) 39 Analysis of the four-point bending test where N(s) is the axial force, Q(s) is the shear force and M(s) is the bending moment. Substituting (2.7) into (2.12): N(s) = −Pφo 2(cos ϕi+φosin ϕi), Q(s) = −P 2(cos ϕi+φosin ϕi),(2.13a) M(s) = −P 2(cos ϕi+φosin ϕi)s. (2.13b) The constitutive equations of the material in a 2D approach yield for a ply p(with p= 1,2, ..., Np): σs=Cp 11εs+Cp 13εz, σz=Cp 13εs+Cp 33εz,(2.14a) τsz =Cp 55γsz,(2.14b) where Cp ij, for i, j = 1,3 and i=j= 5, are the stiffnesses of a ply pin the local coordinate system (s,z) considering a bi-dimensional assumption. Introducing the definition of the strains (2.8) and integrating the constitutive equations (2.14), the classical stiffness matrix A-B-D and the stiffness C are obtained (see Jones [15], Chapter 4) applied to a bi-dimensional problem: N(s) M(s)=A B B Ddu0(s)/ds dφ(s)/ds ,(2.15a) Q(s) = Cφ(s) + dw(s) ds ,(2.15b) where the components A,Band Dare given by: A= Np X p=1 Cp 11(zp−zp−1),(2.16a) B= Np X p=1 Cp 11 1 2(z2 p−z2 p−1),(2.16b) D= Np X p=1 Cp 11 1 3(z3 p−z3 p−1),(2.16c) and the stiffness Cis modified by considering a parabolic distribution of the shear stress: 2.2. Non-linear model for obtaining the bending moment 40 C= Np X p=1 Cp 55 3 2zp−zp−1−4 3t2(z3 p−z3 p−1).(2.17) For the sake of simplicity only symmetrical laminates are considered, and, as a consequence, B= 0. Hence, the displacement component φ(s) is obtained integrating equation (2.15a), considering (2.13b), as follows: dφ(s) ds =M(s) D−→ φ(s) = P 4D(cos ϕi+φosin ϕi)L2 cos2ϕi−s2+φ1, (2.18) where φ1is the value of φ(s) at point 1 (φ1=φ(L/ cos ϕi)). Considering that φ(0) = φ0the next equation for the force Pis obtained: P=4D(φo−φ1) cos2ϕi(cos ϕi+φosin ϕi) L2.(2.19) The displacement component u(s) is obtained integrating equation (2.15a), considering (2.13a), as follows: du(s) ds =N(s) A−→ u(s) = Pφo 2A(cos ϕi+φosin ϕi)L cos ϕi−s+u1,(2.20) where u1is the value of u(s) at point 1. Considering that u(0) = u0the following expression is obtained: (u0−u1) cos ϕi=PLφo 2A(cos ϕi+φosin ϕi).(2.21) Finally, integrating equation (2.15b), the transversal displacement yields: w(s) = P 2C(cos ϕi+φosin ϕi)L cos ϕi−s+φ1L cos ϕi−s+ +P 4D(cos ϕi+φosin ϕi)2L3 3 cos3ϕi−L2 cos2ϕi s+s3 3+w1,(2.22) where w1is the value of w(s) at point 1. Considering that w(0) = w0the following expression is obtained: 41 Analysis of the four-point bending test w1−w0=−P 2C(cos ϕi+φosin ϕi) L cos ϕi− P 4D(cos ϕi+φosin ϕi) 2L3 3 cos3ϕi−φ1 L cos ϕi .(2.23) Substituting (2.19) in (2.23), yields: (w1−w0) sin ϕi=−PL tan ϕi 2C(cos ϕi+φosin ϕi)−L 3(2φo+φ1) tan ϕi.(2.24) Therefore, the length Lis obtained by substituting (2.21) and (2.24) into (2.9), yielding: L= dx+ˆ Rrsin ϕi−φocos ϕi √1+φ2 o +sin ϕi−φ1cos ϕi √1+φ2 1 1 + 1 3(2φo+φ1) tan ϕi+P 2(cos ϕi+φosin ϕi)tan ϕi C−φo A.(2.25) Notice that considering a null load Pand, therefore, null values of φ0 and φ1, the initial length yield: L0=dx+ 2 ˆ Rrsin ϕi.(2.26) Beam 2 is defined between point 1, with a change in the slope of ϕ1as seen previously, and point 2, with a change in the slope of ϕ2. The deformed configuration is depicted in Figure 2.4. A local coordinate system (s2,z) is defined in accordance with Figure 2.4, where s2=s−(L/ cos ϕi). According to Figure 2.4, the length of beam 2, given by the length of the projection of the deformed beam 2 onto the undeformed configuration, is given by L2as follows: L2=Ls 2−Rsin ϕi−ˆ Rrsin(ϕi−ϕ1)−w1sin ϕi+u1cos ϕi.(2.27) Substituting (2.7) into (2.27): L2=Ls 2−Rsin ϕi−ˆ Rr(sin ϕi−φ1cos ϕi) p1 + φ2 1−w1sin ϕi+u1cos ϕi.(2.28) 2.3. Simplification of the non-linear model 48 2.3 Simplification of the non-linear model The model developed in the previous section is complex to solve due to the non-linearity of the equations. In the present section, several hypotheses are applied, which have been verified previously to be applicable in typical configurations using CFRP laminates, in order to obtain a closedform solution of the model. Considering the hypothesis ACD/t2, R2N2M2,EV REI, (φ0−φ1)L2Land φ2 11, the equations are simplified as follows: P=4D(φo−φ1) cos2ϕi(cos ϕi+φosin ϕi) L2,(2.53a) L= dx+ˆ Rrsin ϕi−φocos ϕi √1+φ2 o + sin ϕi−φ1cos ϕi 1 + 1 3(2φo+φ1) tan ϕi ,(2.53b) L2=Ls 2−Rsin ϕi−ˆ Rr(sin ϕi−φ1cos ϕi)−w1sin ϕi+u1cos ϕi,(2.53c) N2=−P(φo−φ1) cos ϕi 2(cos ϕi+φosin ϕi)(cos ϕi+φ1sin ϕi),(2.53d) M2=−PL 2(cos ϕi+φosin ϕi) cos ϕi ,(2.53e) φ1=φ2−M2 D L2 cos ϕi ,(2.53f) w1=w2−M2 D L2 2 2 cos2ϕi +φ2 L2 cos ϕi ,(2.53g) φ2=−R EI M2ϕi,(2.53h) u1=u2=M2 R2 EI (sin ϕi−ϕi),(2.53i) w2=M2 R2 EI (cos ϕi−1).(2.53j) Substituting (2.53h) and (2.53j) in (2.53g) and (2.53f): φ1=−M2R EI ϕi+1 D L2 cos ϕi,(2.54) w1=M2R2 EI (cos ϕi−1) −1 D L2 2 2 cos2ϕi−Rϕi EI L2 cos ϕi.(2.55) 49 Analysis of the four-point bending test Substituting (2.53a) in (2.53e): φ1=M2L 2Dcos ϕi +φo.(2.56) For typical test configurations, the distance L2is very small compared with the distance Land it may be approximated by being invariant with the value of the initial length: L2=Ls 2−Rsin ϕi−ˆ Rrsin ϕi.(2.57) Dividing (2.53d) between (2.53e) and substituting (2.56): N2=−M2 2cos ϕi 2D(cos ϕi+φ1sin ϕi).(2.58) Therefore, neglecting the angle φ1and substituting in (2.39a), the axial force in the central point may be approximated by: N3=−M2 2 2Dcos ϕi .(2.59) Substituting (2.56) in (2.54) and in (2.53b): φo=−M2R EI ϕi+1 Dcos ϕiL2+L 2,(2.60) L= dx+ˆ Rrsin ϕi−φocos ϕi √1+φ2 o + sin ϕi−φocos ϕi−M2L 2D 1 + φo+M2L 6Dcos ϕitan ϕi .(2.61) Using the Taylor expansion in φ0= 0 in the numerator of (2.61) and neglecting the powers higher or equal than two (considering φ01): L=dx+ˆ Rr2 sin ϕi−2φocos ϕi−M2L 2D 1 + φo+M2L 6Dcos ϕitan ϕi .(2.62) Substituting (2.60) in (2.62) and developing the equation: L=dx+R EI ϕi+1 Dcos ϕiL2+L 3M2Ltan ϕi+ ˆ Rr2 sin ϕi+ 2M2R EI ϕicos ϕi+1 DL2+L 4.(2.63) 2.4. Comparison of the non-linear model with experimental results 50 which constitutes a second-order equation in Lwhich may be solved analytically: L=−BL−pB2 L−4ALCL 2AL ,(2.64) AL=M2 3Dcos ϕi tan ϕi,(2.65) BL=L2M2 Dcos ϕi tan ϕi+RM2 EI ϕitan ϕi+ˆ RrM2 2D−1,(2.66) CL=dx+ 2 ˆ Rrsin ϕi+RM2 EI ϕicos ϕi+M2L2 D.(2.67) Therefore, the resolution of the model is reduced to consider a given value of the bending moment in the symmetry point M3∼ =M2<0 and calculating the length Lfrom (2.64), where the length L2is given by (2.57). Once Lhas been obtained the angle φ0is obtained from (2.60) and the angle φ1from (2.54). Finally, the load necessary to obtain the bending moment M3is calculated from (2.53a), and the displacement of the rollers is given by (2.52). If desired, the axial force at the central point is given by (2.59). If the input parameter is the load Pinstead of the bending moment M3, it is necessary to do an iterative process to obtain the desired load. 2.4 Comparison of the non-linear model with experimental results For the validation of the non-linear model, the results obtained with it are compared with experimental results. The experimental tests have been carried out by Ram´ırez [102]. The tests have been carried out over a standard CFRP which has been previously characterized, obtaining the material properties given in Table 2.1. E11 152 GPa E22 8.7 GPa ν12 0.37 G12 10 GPa Table 2.1: Ply properties Two kinds of specimens are considered: •Specimen SP1: [45,0,0,0,-45]S,t= 2 mm, R= 6 mm, W= 27 mm. 51 Analysis of the four-point bending test •Specimen SP2: [45,0,0,-45,90,45,0]S,t= 2.7 mm, R= 6.35 mm, W= 27 mm. The 0odirection is defined as the curved direction, while the 90odirection is defined as the through-the-width direction. The test configuration is defined by a roller radius of Rr= 5 mm, and the distances are Ls= 22 mm and dx= 7 mm. Three coupons of each kind of specimen are tested, which are numbered as SPx.y, where x indicates the type of the specimen (1 or 2) and y indicates the number of coupon (from 1 to 3). The numerical evaluation has been carried out by considering the hypotheses of Plane Strain and Plane Stress for obtaining the constitutive equations (2.14) from the 3D constitutive equations. Additionally, the Classic Laminate Theory (CLT) has been applied to obtain the classical 6x6 A-B-D matrix (see for instance [15], Chapter 4) and, considering that the out-of-plane forces and moments are null (the forces and moments associated to the transverse stress σyand the shear stress τsy, where yis the throughthe-width coordinate), so that an equivalent 2x2 A-B-D matrix for equation (2.15a) may be obtained by using the aforementioned consideration. Figures 2.5 and 2.6 show for SP1 and SP2 specimens respectively the numerical results obtained with the Plane Strain assumption, the Plane Stress assumption and with the 2x2 A-B-D matrix obtained from the 6x6 A-B-D matrix of the CLT. Those numerical results are compared with the experimental results obtained in the tests. Notice that the Plane Stress assumption is not physically realistic as it implies discontinuities in the through-the-width normal strain. Furthermore, in laminates with ±45o, Plane Strain assumption is also non realistic as implies discontinuities in the through-the-width interlaminar shear stress. Therefore, only the model based on the CLT A-B-D matrix is physically realistic. It may be observed that the results of the model evaluated with the stiffness properties of the CLT are always located between the results of the Plane Strain assumption and the Plane Stress assumption. The Plane Strain assumption is the stiffest assumption, and the Plane Stress assumption is the less stiff assumption. It can be obtained that the range between both Plane Strain and Plane Stress assumptions is higher when higher is the number of ±45oplies. Notwithstanding, the model based on the CLT is usually nearer to the Plane Stress model. The experimental results show a high agreement with the model based on the CLT approximation. It shows a high accuracy of this assumption compared with the Plane Strain and the Plane Stress assumptions, especially 2.4. Comparison of the non-linear model with experimental results 52 Figure 2.5: Numerical and experimental force-displacement relations in SP1. Figure 2.6: Numerical and experimental force-displacement relations in SP2. in laminates with a high percentage of ±45oplies. The non-linearity and the slope of the curves modified by the change of the point of contact of the specimen with the rollers are captured with a 53 Analysis of the four-point bending test high accuracy by the analytical model. 2.5 Comparison of the non-linear model with the ASTM procedure The model developed may be used to monitor the test in real time, validating that the experimental results are according to the analytical results and detecting if any kind of defects or errors are present during the test. However, the main utility of the model is not the monitoring of the test, but the calculation of the bending moment from the applied load. The calculation of the stresses involving the unfolding failure is done by means of the bending moment, so high errors in this calculation implies a poor estimation of the unfolding failure and, according to it, a poor estimation of the ILTS by using the four-point bending test, which is normalized by the procedure ASTM D6415 [3]. The ASTM procedure requires both load and displacement experimental parameters to estimate the ILTS. However, the analytical procedure developed in the present Chapter only requires one of the two parameters (force or displacement), and the other parameter can be calculated. This is useful for comparison with the experimental load-displacement distribution and validation of the accuracy of the model. In this way, both analytical and ASTM procedures are applied to the experimental specimens SP1 and SP2 commented in the previous section. For introducing the data in the ASTM procedure, specimens SP1 are considered to fail, based on the experimental results, at a mean load of P1= 790 N with a displacement δ1= 3.4 mm, and specimens SP2 are considered to fail, based on the experimental results, at a mean load of P2= 1270 N with a displacement δ2= 2.8 mm. The analytical procedure developed in the present document applied to the specimens SP1 using the load P1results in an slope φ0= 0.124, an slope φ1= 0.041 and a total bending moment MSP1= 11.7 Nm. These slopes, by using (2.7), imply an angle ϕ0= 0.123 rad and ϕ1= 0.041 rad. By using the ASTM procedure, the angle φobtained is 0.698 rad, which implies an increment ϕi−φ= 0.0982 rad, which is located between the angles ϕ0and ϕ1. Therefore, by using 2.3, the bending moment obtained according to the ASTM procedure yields MAST M SP 1= 12.1 Nm, which is a 3% higher than the one obtained from the non-linear model developed in the present document. The analytical procedure developed in the present document applied to the specimens SP2 using the load P2results in an slope φ0= 0.0965, an slope φ1= 0.0289 and a total bending moment MSP1=−20.6 Nm. These 2.5. Comparison of the non-linear model with the ASTM procedure 54 slopes, by using (2.7), imply an angle ϕ0= 0.0962 rad and ϕ1= 0.0289 rad. By using the ASTM procedure, the angle φobtained is 0.7221 rad, which implies an increment ϕi−φ= 0.0764 rad, which is located between the angles ϕ0and ϕ1. Therefore, by using (2.3), the bending moment obtained according to the ASTM procedure yields MAST M SP 1= 21.1 Nm, which is a 2% higher than the one obtained from the non-linear model developed in the present document. Therefore, for the present specimens, the ASTM procedure has very low errors in the calculation of the bending moment applied to the curved zone, presenting a high accuracy. These errors may be higher for specimens with a lower bending stiffness. Notwithstanding, the ASTM procedure requires to obtain the experimental displacement, which is a difficult task in view of the test set-up. An inaccurate experimental measurement of the failure displacement may imply errors in the bending moment calculation higher than the ones presented in the aforementioned experimental results. In contrast, the present analytical method requires only the experimental load, which is easily obtained with a high accuracy from the experimental tests. Chapter 3 Bi-dimensional models for evaluating interlaminar stresses The unfolding failure, which is traditionally associated to the interlaminar tensile stresses, is typically calculated using two different techniques. The fist technique consists in the use of a FE software, which requires a highly refined mesh in the thickness of the laminate. If the mesh in the thickness has to be refined for obtaining a good accuracy in the interlaminar stresses it implies also a refinement in the in-plane direction. This fact causes that a very high number of elements are needed in the numerical model and, hence, high calculation times are required for obtaining a good estimation of the interlaminar stresses. Furthermore, the high calculation times hinder the optimization of geometries and stacking sequences. The second technique consists in the use analytical solutions. The typical analytical solutions are based in bi-dimensional approximations, which reduce the scope of applications. Furthermore, traditional calculation methods for the interlaminar stresses are regularized models. Regularized stresses and strains are defined as the stresses and strains far enough from a perturbation. A typical example of a perturbation is given by the free edge effect, where a stress concentrator is given. However, this is not the only perturbation present. In the case of the traditional sections prone to unfolding failure, another kind of perturbation is given by curvature changes, e.g., in a L-sectioned beam where the curved zone is joint to the straight arms. This effect has been observed numerically by Most et al. [14]. The perturbation due to the change of curvature allows equalling the displacements and stresses between both parts, while the regularized solutions are not similar. Respect to the interlaminar normal stress (INS), the INS has a through-the-thickness distribution in the section of the change of curvature given approximately (for the case of a bending moment loading) 3.1. Theoretical basis of the regularized models 56 by the average distribution between the regularized one of one side and the regularized one of the adjacent side. When another kind of load is applied, the INS distribution is also reduced in the side where the highest values are given. It implies that when the maximum regularized value of the INS is given just at the joint, the non-regularized value is actually much lower due to the change of curvature perturbation and, therefore, regularized models may underestimate the failure. Accordingly, there is a necessity to introduce a non-regularized analytical model which allows calculating in a quicker and more accurate way the non-regularized stresses, getting a similar accuracy that of the detailed FE models, but with lower computational times. This chapter develops two non-regularized models based on a series expansion of the displacements and on the definition of higher-order moments in the stresses. Before that, the basis of the regularized models developed by Lekhnitskii [19, chap. 3] and Ko and Jackson [4] are presented, as well as an approximation of the regularized model based on the Timoshenko beam theory developed by the authors [94], which constitutes the first-order model of the non-regularized models. 3.1 Theoretical basis of the regularized models The models developed by Lekhnitskii [19, chap. 3] and Ko and Jackson [4] are developed in the present section to make the notation uniform with the notation used in the project. Regularized models developed by Lekhnitskii [19, chap. 3] and Ko and Jackson [4] constitute the solution of an isolated curved beam loaded under end axial and shear forces, N(0) = N0and Q(0) = Q0, and bending moment, M(0) = M0, see Fig. 3.1. These forces and moment are defined from the circumferential and shear stresses, σθ(r, θ) and τrθ(r, θ), as follows: N(θ) = ZR+ R− σθ(r, θ)dr, M(θ) = ZR+ R− (r−R)σθ(r, θ)dr, (3.1a) Q(θ) = ZR+ R− τrθ(r, θ)dr, (3.1b) where ris the radial coordinate, θthe circumferential coordinate, Ris the mean radius, R−is the inner radius, R+is the outer radius and tis the thickness. Notice that forces and the moment are integrated only across the 57 Bi-dimensional models for evaluating interlaminar stresses thickness and not along the width (defining the width as the out-of-plane dimension), so they are actually forces and moment per unit width. Figure 3.1: Geometry and loads considered in the regularized models of a curved beam. Boundary conditions at the ends of the curved beam assure only that the mean value and the bending moment of the stress distribution are equals to the applied loads. Therefore, the solution is only accurate at a particular distance from the ends. That distance is denominated Regularization Distance (RD) and is analysed in a further section. The forces and the moment have a circumferential distribution due to the equilibrium equations as follows: N(θ) = N0cos θ−Q0sin θ, (3.2a) Q(θ) = Q0cos θ+N0sin θ, (3.2b) M(θ) = M0+R(N0−N(θ)).(3.2c) Lekhnitskii’s equations were developed for anisotropic homogeneous materials, considering the stress-strain relations of the material given by the following equations: εθ(r, θ) = 1 Eθ σθ(r, θ)−νθr Eθ σr(r, θ),(3.3a) εr(r, θ) = −νθr Eθ σθ(r, θ) + 1 Er σr(r, θ),(3.3b) γrθ(r, θ) = 1 Grθ τrθ(r, θ),(3.3c) where Eris the radial stiffness, Eθthe circumferential stiffness, Grθ is the shear stiffness in the bi-dimensional plane and νθr is the Poisson ratio. These 3.2. Approximations of the regularized models 64 eN(θ) = 1 R duo(θ) dθ +w1(θ) R, eM(θ) = 1 R du1(θ) dθ ,(3.21b) where eN(θ) and eM(θ) are the Classical Beam Theory (CBT) axial and bending strains. The strains vector, eT(θ), is defined as follows: eT(θ) = eN(θ) eM(θ).(3.22) Therefore, equation (3.21a) can be expressed in matrix notation as follows: εθ(θ, z) = TT(r)eT(θ),TT(r) = R r R r(r−R).(3.23) Neglecting the radial strain (due to the hypothesis of the constant value of the radial displacement) the relation between the circumferential stress and the strains in a ply p, according to (3.5), depends on the stiffness of the ply, Cp θθ, as follows: σp θ(θ, z) = Cp θθεθ(θ, z) = Cp θθTT(r)eT(θ)=(Sp θ(r))TeT(θ).(3.24) Plies are enumerated beginning with p= 1 in the ply situated in the lower radius up to p=Nplocated in the outer radius. Couplings between shear and axial stresses and strains are not given, since the material is considered to be orthotropic in the θ-raxes. Therefore, the relation between the circumferential stress and the strain vector is given by Sp θ(r), defined according to (3.24) as follows: Sp θ(r) = Cp θθT(r) = R rCp θθ 1 r−R.(3.25) Integrating the circumferential stress in order to obtain the forces vector, next expression is obtained: MT(θ) = ZR+ R−1 r−Rσp θ(r, θ)dr =ZR+ R−1 r−R(Sp θ(r))TeT(θ)dr =ˆ KσeT(θ),(3.26) where the equivalent homogeneous material stiffness matrix is given by: ˆ Kσ=A0B0 B0D0,(3.27) 65 Bi-dimensional models for evaluating interlaminar stresses with: A0= Np X p=1 Cp θθRlog R+ p R− p,(3.28a) B0= Np X p=1 Cp θθRR+ p−R− p−Rlog R+ p R− p,(3.28b) D0= Np X p=1 Cp θθRR2log R+ p R− p+R+ p−R− p 2(R+ p+R− p−4R),(3.28c) where R+ pand R− pare the values of rat the outer and inner radius respectively for the ply pand Npis the number of plies. Notice that R+ p=R− p+1. The curvature radius causes a coupling effect between the axial and bending parameters (B06= 0) even when a symmetric laminate is given. Substituting (3.26) in (3.24), the circumferential stress can be calculated from the axial force and the bending moment as follows: σp θ(θ, r) = (Sp θ(r))Tˆ Kσ−1MT(θ).(3.29) Once the circumferential stress has been obtained, the other stresses can be calculated by equilibrium. The equilibrium equations of the elasticity in polar coordinates are given by: ∂σθ ∂θ +1 r ∂(r2τrθ) ∂r = 0,(3.30a) ∂(rσr) ∂r +∂τrθ ∂θ =σθ.(3.30b) Thus, substituting (3.29) in (3.30) and integrating, the shear and radial stresses are given by: τp rθ(θ, r) = τp rθ(R− p, θ)R− p r2 −(Sp rθ(r))Tˆ Kσ−1dMT(θ) dθ ,(3.31a) σp r(θ, r) = σp r(R− p, θ)R− p r+dτp rθ(R− p, θ) dθ R− p rR− p r−1 + (Sp r,1(r))Tˆ Kσ−1MT(θ)+(Sp r,2(r))Tˆ Kσ−1d2MT(θ) dθ2,(3.31b) 3.2. Approximations of the regularized models 66 where: Sp rθ(r) = 1 r2Zr R− p r0Sp θ(r0)dr0=Cp θθR r2"r−R− p r−R− p 2(r+R− p−2R)#,(3.32a) Sp r,1(r) = 1 rZr R− p Sp θ(r0)dr0=Cp θθR r"log r R− p r−R− p−Rlog r R− p#,(3.32b) Sp r,2(r) = 1 rZr R− p Sp rθ(r0)dr0 =Cp θθR r  log r R− p+R− p r−1 r 2−R− p−R−R− p 2R− p r+R1−log r R− p ,(3.32c) where Sp rθ(r) relates the shear stress with the derivative of the strains vector, Sp r,1(r) relates the radial stress with the strain vector and Sp r,2(r) relates the radial stress with the second derivative of the strain vector. The boundary conditions for τp rθ(R− p, θ) and σp r(R− p, θ) are obtained from the previous lamina p−1 so that τp rθ(R− p, θ) = τp−1 rθ (R+ p−1, θ) and σp r(R− p, θ) = σp−1 r(R+ p−1, θ), beginning with a boundary condition for the first ply p= 1 of σ1 r(R−, θ) = σ−(θ) and τ1 rθ(R−, θ) = τ−(θ). The boundary conditions of σNp r(R+, θ) = σ+(θ) and τNp rθ (R+, θ) = τ+(θ) are automatically satisfied. In conclusion, the equations of the method may be easily solved to calculate the stresses in two simple steps. The first step consists in the calculation of the axial force and the bending moment from equations (3.17). For example, in the case of constant distributed loads (q(θ) = q0,p(θ) = p0, m(θ) = m0) the forces and the bending moment yield: N(θ)=(N0−Rp0) cos θ−(Q0+Rq0) sin θ+Rp0,(3.33a) Q(θ) = (Q0+Rq0) cos θ+ (N0−Rp0) sin θ−Rq0,(3.33b) M(θ) = M0+R(N0−N(θ)) −R(m0+Rq0)θ. (3.33c) Once forces and moment have been obtained, the second step is to calculate stresses. Circumferential stress is directly obtained from equation (3.29). Shear stress is calculated recursively from (3.31a), beginning with a boundary condition of τrθ(R−, θ) = τ−(θ). Deriving equation (3.31a) respect to θ,dτrθ(r, θ)/dθ is obtained recursively too, beginning with a boundary condition of dτrθ(R−, θ)/dθ =dτ−(θ)/dθ. Finally, the radial stress 67 Bi-dimensional models for evaluating interlaminar stresses is obtained recursively from (3.31b) beginning with a boundary condition σr(R−, θ) = σ−(θ). In the resolution procedure it is necessary to define three important matrices and vectors: •MT(θ) vector, containing the axial force and the bending moment. •Stiffness matrix ˆ Kσ, relating the axial force and the bending moment with the mean section strains. •Sp(r) vectors, relating the stresses with the mean section strains and their derivatives. They are calculated from the stress-strain relations and the T(r) vector, which relates the strain with the mean section strains. If the matrix equations are developed in conventional form by multiplying the matrices, equations of the model given by [94] are obtained, as the same hypotheses have been considered. Therefore, results of the present model are the same than those presented in the aforementioned reference. Some comparisons with the Lekhnitskii’s equations, the Ko and Jackson’s equations and FE results can be found in [103]. 3.3 Theoretical basis of the non-regularized models The non-regularized models developed in the present project are based on a series expansion of the displacements. Several authors have tried to model several kinds of structural problems based also in series expansions. Carrera et al. [34] presented a compilation of several methods based on series expansions in the displacements applied to vibration analysis. Several kinds of those models are commented in Chapter 1. Among the models based on monomials, the one developed by Matsunaga [37] is of special interest in the present work. Matsunaga considers, for flat laminates, a series expansion with monomial functions in the in-plane displacements, including n+ 1 terms in these displacements, and another series expansion with monomial functions in the out-of-plane displacement, including nterms in this displacement. The approximation by using one order less in the out-of-plane displacement than in the in-plane ones is typical also in the first order models, such as the Timoshenko beam theory [35,36] or the CLT [15]. Later, Matsunaga applied his model to curved laminates [38] 3.3. Theoretical basis of the non-regularized models 68 by using also a series expansion with monomial functions in the two displacements. Matsunaga’s model uses a series expansion based on monomial functions in the displacements in the next way: us(s, z) = uo(s) + zu1(s) + ∞ X i=2 ziun(s),(3.34a) uz(s, z) = w1(s) + ∞ X i=2 zi−1wn(s),(3.34b) where a curvilinear coordinate system (s, z) is used (see Fig. 3.3), defined by the axial coordinate sand the through-thickness coordinate z. This curvilinear coordinate system is related with the previously defined polar coordinate system (θ, r) by the relations s=Rθ and z=r−R. Figure 3.3: Definition of the curvilinear coordinate system. Approximating the displacements by a n-order model instead of the infinite summation: us(s, z) = uo(s) + zu1(s) + n X i=2 ziui(s),(3.35a) uz(s, z) = w1(s) + n X i=2 zi−1wi(s).(3.35b) Notice that the through-thickness displacement, uz(s, z), is approximated with one order less than the axial displacement, us(s, z). The present model is based on a similar series expansion but by using different functions in the series expansion in the different displacements, which allows an easier and a more efficient numerical resolution of the problem to be performed. Furthermore, Matsunaga used a secondary series expansion 69 Bi-dimensional models for evaluating interlaminar stresses in the parameters that has been suppressed in the present model. In this way, axial displacement is expanded by using functions fs i(z), and thoughthickness displacements are expanded by using different functions fz i(z), so they yield: us(s, z) = uo(s) + zu1(s) + n X i=2 fs i(z)ui(s),(3.36a) uz(s, z) = w1(s) + n X i=2 fz i(z)wi(s).(3.36b) Different kinds of non-regularized models can be chosen depending on the functions used in the series expansion of the axial displacements, fs i(z). The present document develops two kind of models. The first one is based on monomial functions (fs i(z) = zi), and is called monomials based model (MBM). Due to numerical instabilities in the stiffness matrix of the MBM, a second model is developed based on Legendre polynomials (fs i(z) = pi(2z/t), where pi(2z/t) is the i-th Legendre polynomial), and is called Legendre polynomials based model (LPBM). The functions used for the through-thickness displacement, fz i(z), are determined by the shear strains. In the CBT, which constitutes the first order model, axial strain depends on two CBT 1D strain components given in (3.21b), shear strain depends only on one CBT 1D strain component and the through-thickness strain is neglected. If a n-th order model is chosen it is desirable to increment in n−1 components each strain, so that the axial strain is defined by n+1 components, the through-thickness strain is defined by n−1 components and the shear strain is defined by ncomponents. In that way, the condition to determine the functions for the series expansion in the through-thickness displacement is that shear strain depends only on ncomponents. The aim of developing these non-regularized models is to estimate the perturbation introduced by a change of curvature, so that a more accurate stresses distribution may be obtained in this kind of geometries. The geometries considered consist in a section of a curved composite laminate, treated as a chain of constant curvature beams. Each beam may be composed by different materials or different stacking sequences, but all beams must have the same thickness, which is constant along the axial coordinate. In that way, some examples of geometries in which the model can be applied are sections of L-shaped components, Ω-shaped components, joggles or corrugated laminates (see Fig. 3.4). Notice that typically the sections end in a long enough straight part in which stresses tend to their regularized values. Therefore, they can be 3.3. Theoretical basis of the non-regularized models 70 Figure 3.4: Examples of geometries where the non-regularized models can be applied. Global coordinate system definition. considered semi-infinite. Another kind of boundary conditions can be chosen at the ends of the chain, e.g., an embedded end. A global coordinate system (sG,zG) is considered in the whole section as shown in Fig. 3.4. The axial coordinate sGis a curvilinear coordinate along the mean line of the section. The origin is located in the first joint of the chain in the case of a semi-infinite end beam, and it is oriented towards the following joint. The through-thickness coordinate zGhas its origin in the mean line and the orientation of it can be chosen by the user. The section can be decomposed in several constant-curvature beams, and a local coordinate system (sk,zk) is considered for each component k(see Fig. 3.5), where skhas the same sense than sGbut with its origin in the joint with the previous beam, and zkcoincides with the global through-thickness coordinate zG. Therefore, three kinds of components can be considered: •Straight beams (Beams 1, 3 and 5 in Fig. 3.5). •Clockwise curved beams, which have a local through-thickness coordinate zkdefined in the same sense than the radial coordinate (Beam 2 in Fig. 3.5). •Counter-clockwise curved beams, which have a local through-thickness coordinate zkdefined positive towards the center of curvature (Beam 4 in Fig. 3.5). The problem is solved by considering each component isolated and considering arbitrary stresses boundary conditions applied at both ends of the components (or applied at the only one end for the case of the semi-infinite 71 Bi-dimensional models for evaluating interlaminar stresses Figure 3.5: Decomposition of the section in several constant-curvature beams. end components). Displacements and stresses are obtained in the component depending on the arbitrary boundary conditions. Finally, the arbitrary boundary conditions are calculated by imposing continuity of the displacements and the stresses at both sides of every joint of the chain of beams. 3.4 Monomials based model (MBM) The first model developed in this Thesis can be obtained by choosing the monomials functions as fs i(z) functions, similarly to what is done in the Matsunaga models [37,38]. The difference with these models is given in the definition of the fz i(s) functions, which are not necessarily monomials. Therefore, the MBM is based on an approximation of the displacements as follows: us(s, z) = uo(s) + zu1(s) + n X i=2 ziui(s),(3.37a) uz(s, z) = w1(s) + n X i=2 fz i(z)wi(s).(3.37b) However, as will be seen later, for a better numerical resolution it is recommendable to use non-dimensional parameters, at least non-dimensional displacements and lengths and optionally non-dimensional stresses and pressures. Numerical errors are minimized when displacements and lengths are adimensionalized with the half of the thickness. Non-dimensional parameters are marked with an upper bar in the magnitude, so if uis a generic 3.4. Monomials based model (MBM) 72 displacement/length and σis a generic stress/pressure, the non-dimensional parameters ¯uand ¯σrespectively are defined as: ¯u=2u t,¯σ=σ Eref ,(3.38) where tis the thickness and Eref is a reference stiffness, e.g., the stiffness in the fibre direction. Hence, the series expansion using non-dimensional parameters can be developed alternatively to 3.37 as follows: ¯us(¯s, ¯z) = ¯uo(¯s) + ¯z¯u1(¯s) + n X i=2 ¯zi¯ui(¯s),(3.39a) ¯uz(¯s, ¯z) = ¯w1(¯s) + n X i=2 ¯ fz i(¯z) ¯wi(¯s).(3.39b) For the sake of simplicity, the bar indicating the non-dimensional magnitudes is obviated in the present section while all the magnitudes used are non-dimensional unless it is said otherwise. Nevertheless, the model can be implemented also with the dimensional magnitudes, obtaining a worse numerical behaviour. The model is developed for the isolated beam under arbitrary boundary conditions, and in the next sections the joint to other beams will be analysed. The model is developed first for the case of straight beams, and after that it is developed for curved beams, jointly for both clockwise and counterclockwise curved beams. This model has been published by the author in [95] considering only end loads applied. However, the model developed here considers also distributed loads in the transverse boundaries. 3.4.1 Development of the MBM for straight beams A straight beam with thickness t(considered non-dimensional with value t= 2) and length Lkis studied (see Figure 3.6). A local non-dimensional coordinate system skand zkis used as defined in Figure 3.5. The beam is loaded under distributed transverse and shear loads in the positive boundary of zk,σ+(sk) and τ+(sk), and under distributed transverse and shear loads in the negative boundary of zk,σ−(sk) and τ−(sk). Arbitrary stresses boundary conditions are applied in sk= 0, given by σ0 s(zk) and τ0 s(zk), and in sk=Lk, given by σL s(zk) and τL s(zk). For the sake of simplicity, the index kindicating the component is suppressed in the following. 73 Bi-dimensional models for evaluating interlaminar stresses Figure 3.6: Straight beam. Displacements The longitudinal and through-thickness displacements approximation of the straight beam is given by equations (3.39), where the axial displacements have been approximated by a series expansion by using monomial functions. The functions fz i(z) used for the through-thickness displacements have to be determined. The series expansion can be expressed in matrix notation by the definition of the following vectors of length n−1: uT(s) = u2(s)u3(s)... un(s),(3.40a) wT(s) = w2(s)w3(s)... wn(s),(3.40b) fsT(z) = z2z3... zn,(3.40c) fzT(z) = fz 2(z)fz 3(z)... fz n(z).(3.40d) Therefore, equation (3.39) can be expressed as follows: us(s, z) = uo(s) + zu1(s) + fsT(z)u(s),(3.41a) uz(s, z) = w1(s) + fzT(z)w(s).(3.41b) Strains Strains equations in the local Cartesian coordinate system yield: εs(s, z) = ∂us(s, z) ∂s , εz(s, z) = ∂uz(s, z) ∂z ,(3.42a) γsz(s, z) = ∂us(s, z) ∂z +∂uz(s, z) ∂s .(3.42b) 3.4. Monomials based model (MBM) 80 d2x(s) ds2=Gx(s) + gNN(s) + gMM(s) +gσ+σ+(s) + gσ−σ−(s) + gτ+ dτ+(s) ds +gτ− dτ−(s) ds ,(3.66a) x(s) = Ms(s) w(s),(3.66b) where the auxiliary matrices appearing in the previous equation (G,gN, gM,...) are given in Appendix A.1. The resolution of the previous equation is divided in two parts: obtaining a particular solution and solving the homogeneous equation. Notice that the axial force, the bending moment and the distributed loads affect only to the particular solution, which constitutes the regularized solution of the problem. This regularized solution coincides with the CLT solution. The homogeneous equation let us to calculate the perturbations over the regularized values so that the continuity conditions of stresses and displacements at the joint are satisfied. Constant distributed loads case Considering the straight beam under constant distributed loads σ+(s) = σ+ 0,σ−(s) = σ− 0,τ+(s) = τ+ 0and τ−(s) = τ− 0as a particular case, using equilibrium equations (3.53) and (3.54) and using the boundary conditions (3.55), the axial force, shear force and bending moment yield: N(s) = N0−(τ+ 0−τ− 0)s, (3.67a) Q(s) = Q0−(σ+ 0−σ− 0)s, (3.67b) M(s) = M0+Q0−t 2(τ+ 0+τ− 0)s−(σ+ 0−σ− 0)s2 2.(3.67c) Therefore, substituting (3.67) in (3.66a) the particular/regularized solution yields: xreg(s) = qNN(s) + qMM(s) + qσ+σ+ 0+qσ−σ− 0,(3.68) where auxiliary vectors qN,qM,qσ+and qσ−are defined by: qN=−G−1gN,qM=−G−1gM,(3.69a) qσ+=G−1(G−1gM−gσ+),qσ−=−G−1(G−1gM+gσ−).(3.69b) 81 Bi-dimensional models for evaluating interlaminar stresses The derivative of xreg(s), substituting (3.53a) and (3.54), yields: dxreg(s) ds =qQQ(s) + qτ+τ+ 0+qτ−τ− 0,(3.70) where the auxiliary vectors are given by: qQ=qM,qτ+=−qN−qM t 2,qτ−=qN−qM t 2.(3.71) The resolution of the homogeneous equation is analysed in section 3.7. If a not constant distributed loads is applied, a similar procedure may be developed for obtaining the solution of the homogeneous equation. Stresses calculation The stresses calculation is developed in detail in [95]. In summary, the stresses calculation in the MBM is similar than the calculation in the model of order 1 developed in section 3.2. Once the vector x(s), given in (3.66b), has been calculated, the axial higher order moments Ms(s) are directly obtained as the first n−1 components of x(s). Shear higher-order moments and transverse higher-order moments may be obtained from the equilibrium equations (3.57) and (3.59), respectively. Strains may be obtained by using the compliance equations (3.63). Once strains are obtained, stresses can be calculated in each ply by using the 2D constitutive equations (3.60). Notwithstanding, only axial stress calculated with the aforementioned procedure results accurate. Therefore, similarly than in the CBT, the other two stresses (the shear stress and the transverse stress) are obtained from the axial stress by using the equilibrium equations (3.51). In that way, stresses can be written as: σp θ(s, z) = Sp s(z)Tˇ KσMT(s),(3.72a) τp sz(s, z) = τp sz(s, z− p)−Sp sz(z)Tˇ Kσ dMT(s) ds ,(3.72b) σp z(s, z) = σp z(s, z− p)−dτp sz(s, z− p) ds z−z− p+Sp z,2(z)Tˇ Kσ d2MT(s) ds2,(3.72c) where MT(s) is a vector containing the axial and transverse moments: MT(s) =     N(s) M(s) Ms(s) Mz(s)     ,(3.73) 3.4. Monomials based model (MBM) 82 and the shape vectors Sp s(z) (defined by its components Sp si(z) for i= 0,1, ..., 2n−1), Sp sz(z) and Sp z,2(z) are given by: Sp si(z) = Qp 11zi, i = 0,1, ..., n, (3.74a) Sp si+n−1(z) = Qp 13zi−2, i = 2,3, ..., n, (3.74b) Sp sz(z) = Zz z− p Sp s(z0)dz0,(3.74c) Sp z,2(z) = Zz z− p Sp sz(z0)dz0.(3.74d) Notice that the shear stress and the transverse stress are obtained by integration of the equilibrium equations, which implies that boundary conditions for these stresses are required. In the case of the shear stress, τp sz(s, z− p) is required at the bottom of each ply p, which is obtained by the corresponding value at the top of the previous ply τp sz(s, z− p) = τp−1 sz (s, z+ p−1), beginning with a boundary condition for the first ply of τ1 sz(s, −t/2)=τ−(s). The boundary condition τNp sz (s, t/2) = τ+(s) is automatically accomplished. The same case is given in the transverse stress, where the boundary condition is obtained from the previous ply σp z(s, z− p) = σp−1 z(s, z+ p−1) beginning with a boundary condition σ1 z(s, −t/2) = σ−(s). The boundary condition σNp z(s, t/2) = σ+(s) is automatically accomplished. 3.4.2 Development of the MBM for curved beams A curved beam with thickness t(considered non-dimensional with value t= 2) and length Lk=RkΘk(where Rkis the mean radius and Θkis the angular length of the beam) is studied (see Figure 3.7). A local nondimensional coordinate system skand zkis used, similarly to that used for the straight beam. Notice that the only difference between the clockwise and the counter-clockwise curved beam lies in the direction of the zkaxis. The beam is loaded under distributed transverse and shear loads in the positive boundary of zk,σ+(sk) and τ+(sk), and under distributed transverse and shear loads in the negative boundary of zk,σ−(sk) and τ−(sk). Arbitrary stresses boundary conditions are applied in sk= 0, given by σ0 s(zk) and τ0 s(zk), and in sk=Lk, given by σL s(zk) and τL s(zk). These stresses must be in equilibrium. For the sake of simplicity, the index kindicating the component is suppressed in this section. 83 Bi-dimensional models for evaluating interlaminar stresses (a) (b) Figure 3.7: (a) Clockwise curved beam. (b) Counter-clockwise curved beam. Displacements The displacements approximation for this case is analogous to that made in the straight beam: us(s, z) = uo(s) + zu1(s) + fsT(z)u(s),(3.75a) uz(s, z) = w1(s) + fzT(z)w(s),(3.75b) where vector fs(z) contains the monomials with orders from 2 to nas given in (3.40c). Vector fz(z) has to be obtained from the shear strain similarly as done in the straight beam model. The calculation of this vector can be found in [95], and its components fz i(z), with i= 2,3, ..., n, are given by: fz i(z) = zi−1±i−1 iR zi, i = 2,3, ..., n. (3.76) When a symbol ±or ∓is used in the present section, the upper symbol is used for the clockwise beam, and the lower symbol is used for the counterclockwise beam. 3.4. Monomials based model (MBM) 84 Notice that the derivative of vector fz(z) yields: dfz(z) dz =R±z RHzˆ ˆ fs(z),(3.77) where matrix Hzand vector ˆ ˆ fs(z) are the same than in the straight beam, given by equations (3.49). Furthermore, according to (3.76), vector fz(z) can be expressed depending on vector fs(z) as follows: fz(z) = ±R(Hsfs(z) + hMz),(3.78) where matrix Hsis defined by its components Hs ij, with i, j = 2,3, ..., n, and vector hMis defined by its components hM i, with i= 2,3, ..., n: Hs ij =(i−1) δij i R2±δ(i−1)j R, hM i=±δi2 R.(3.79) Strains Strains definitions in the local polar coordinate system are given by: εs(s, z) = R R±z∂us(s, z) ∂s ±uz(s, z) R,(3.80a) γsz(s, z) = R R±zR±z R ∂us(s, z) ∂z ∓us(s, z) R+∂uz(s, z) ∂s ,(3.80b) εz(s, z) = ∂uz(s, z) ∂z .(3.80c) Substituting (3.75) in (3.80) and considering (3.77) and (3.78), strains can be written as: εs(s, z) = R R±zheN(s) + zeM(s) + fsT(z)es(s)i,(3.81a) γsz(s, z) = R R±zheQ(s) + fzT(z)eγ(s)i,(3.81b) εz(s, z) = R±z Rhˆ ˆ fs T(z)ez(s)i,(3.81c) where eN(s), eM(s) and eQ(s) are the CBT 1D strain components, and es(s), eγ(s) and ez(s) are the higher-order 1D strain components, which are defined as: 85 Bi-dimensional models for evaluating interlaminar stresses eN(s) = du0(s) ds ±w1(s) R,(3.82a) eM(s) = du1(s) ds +hMTw(s),(3.82b) eQ(s) = u1(s)∓u0(s) R+dw1(s) ds ,(3.82c) es(s) = du(s) ds +HsTw(s),(3.82d) eγ(s) = dw(s) ds +HγTu(s),(3.82e) ez(s) = HzTw(s).(3.82f) Equilibrium As in the straight beam case, the distributed loads appear only in the equilibrium equations, being the main difference with respect to the model developed in [95]. The 2D equilibrium equations expressed in the local polar coordinate system yield: R∂σs(s, z) ∂s +1 R±z ∂ ∂z (R±z)2τsz(s, z)= 0,(3.83a) ∂ ∂z [(R±z)σz(s, z)] + R∂τsz(s, z) ∂s =±σs(s, z),(3.83b) The CBT forces and moment and the higher-order moments are defined, as in the straight beam, by equations (3.52). It is demonstrated in [95] that functions fτ(z) and fσ(z) yield for the present case: fτ(z) = fz(z),fσ(z) = R±z R2ˆ ˆ fs(z).(3.84) In that way, integrating the equilibrium equations (3.83) and also the first one of them multiplied by z, the CBT forces and moment equilibrium equations yield: RdN(s) ds ±Q(s) + R+τ+(s)−R−τ−(s)=0,(3.85a) RdQ(s) ds +R+σ+(s)−R−σ−(s) = ±N(s),(3.85b) dM(s) ds +t 2RR+τ+(s) + R−τ−(s)=Q(s).(3.85c) 3.4. Monomials based model (MBM) 86 where R+is the radius in the positive boundary of zand R−is the radius in the negative boundary of z, given by: R+=R±t 2, R−=R∓t 2.(3.86) The CBT forces, N(s) and Q(s), and moment, M(s), are obtained directly by integrating equations (3.85) and introducing the boundary conditions given by equations (3.55). Prescribed values of these CBT forces and moment in s=Lmust be in equilibrium with the prescribed values in s= 0 and, therefore, boundary conditions in s=Lare automatically accomplished. The higher-order moments equilibrium equation in direction scan be obtained by integrating equation (3.83a) multiplied by fs(z), which yields: dMs(s) ds +R+ Rfst 2τ+(s)−R− Rfs−t 2τ−(s) = HγMτ(s).(3.87) Finally, the higher-order moments equilibrium equation in direction z can be obtained by integrating equation (3.83b) multiplied by fz(z) and using properties (3.77) and (3.78), which yields: dMτ(s) ds +R+ Rfzt 2σ+(s)−R− Rfz−t 2σ−(s) =HzMz(s) + HsMs(s) + hMM(s).(3.88) Therefore, higher-order moments equilibrium is given by equations (3.87) and (3.88), which cannot be directly solved as in the case of the CBT forces and moment. The boundary conditions at both ends for these higher order moments may be obtained from the arbitrary boundary conditions, by integrating them according to the higher-order moments definition given in (3.52). Constitutive law As in the straight beam case, the 2D constitutive law for an individual ply is given by equations (3.60). Substituting in these equations the strains, given by (3.81), and integrating the equations in order to obtain the forces and moments (including the higher-order one), relations (3.61) are obtained. However, in this case the components of matrices ˆ Kσand ˆ Kτare influenced by the curvature. Expressions of these matrices can be found in the Appendices of [95]. 87 Bi-dimensional models for evaluating interlaminar stresses The compliance matrices can be defined by inverting the stiffness matrices, ˇ Kσ=ˆ Kσ−1and ˇ Kτ=ˆ Kτ−1, so that decomposing the compliance matrices in several components as seen in [95], compliance equations yield (3.63). Solution procedure The solution procedure is analogous to the one developed in [95] but including the distributed loads. First, substituting (3.87) in (3.88) the transverse higher-order moments yield: Mz(s) = Hz−1R+ Rfzt 2σ+(s)−R− Rfz−t 2σ−(s) +Hγ−1d2Ms(s) ds2+R+ Rfst 2dτ+(s) ds −R− Rfs−t 2dτ−(s) ds  −HsMs(s)−hMM(s)) .(3.89) Introducing (3.63b) into (3.82e) longitudinal higher-order displacements yield: HγTu(s) = ˇ Cττ Mτ(s) + ˇcQτ Q(s)−dw(s) ds ,(3.90) Substituting (3.63a) into (3.82d) and (3.82f), using (3.90) to eliminate the derivative of the longitudinal displacement u(s), using (3.87) to eliminate the higher-order shear moments Mτ(s), and using equation (3.85b) to eliminate the derivative of the shear force Q(s), the following differential equation may be obtained (see [95] and [96] for more details): d2x(s) ds2=Gx(s) + gNN(s) + gMM(s) +gσ+σ+(s) + gσ−σ−(s) + gτ+ dτ+(s) ds +gτ− dτ−(s) ds ,(3.91) where the matrices of the previous equation are given in Appendix A.2 and vector x(s) is defined in (3.66b). The resolution of the previous equation is divided in two parts: obtaining a particular solution and solving the homogeneous equation. Notice that the axial force, the bending moment and the distributed loads affect only to the particular solution, which constitutes the regularized solution of the problem. This regularized solution converges to Ko & Jackson solution [4] 3.4. Monomials based model (MBM) 88 when the order increases. The homogeneous equation let us to calculate the perturbations over the regularized values so that the axial boundary conditions (the arbitrary boundary conditions) are satisfied. Constant distributed loads case Considering the curved beam loaded under constant distributed loads σ+(s) = σ+ 0,σ−(s) = σ− 0,τ+(s) = τ+ 0and τ−(s) = τ− 0as a particular case, using equilibrium equations (3.85) and using the boundary conditions (3.55), the axial force, shear force and bending moment can be written as: N(s) = (N0∓R+σ+ 0±R−σ− 0) cos s R−(±Q0+R+τ+ 0−R−τ− 0) sin s R ±R+σ+ 0∓R−σ− 0,(3.92a) Q(s) = ±(N0∓R+σ+ 0±R−σ− 0) sin s R±(±Q0+R+τ+ 0−R−τ− 0) cos s R ∓R+τ+ 0±R−τ− 0,(3.92b) M(s) = M0±R(N0∓R+σ+ 0±R−σ− 0)1−cos s R ±R(±Q0+R+τ+ 0−R−τ− 0) sin s R∓(R+)2τ+ 0−(R−)2τ− 0s R.(3.92c) Therefore, substituting (3.92) in (3.91), the particular/regularized solution yields: xreg(s) = qNN(s) + qMM(s) + qσ+σ+ 0+qσ−σ− 0,(3.93) where auxiliary vectors qN,qM,qσ+and qσ−are defined by: qN=−(R2G+I)−1(R2gN±RG−1gM),(3.94a) qM=−G−1gM,(3.94b) qσ+=−G−1gσ++ (R2G+I)−1G−1(∓gN+RgM)R+,(3.94c) qσ−=−G−1gσ−+ (R2G+I)−1G−1(±gN−RgM)R−.(3.94d) The derivative of xreg(s), substituting (3.85a) and (3.85c), yields: dxreg(s) ds =qQQ(s) + qτ+τ+ 0+qτ−τ− 0,(3.95) 89 Bi-dimensional models for evaluating interlaminar stresses where the auxiliary vectors are given by: qQ=qM∓qN R,(3.96a) qτ+=R+ R−qN−qM t 2,qτ−=R− RqN−qM t 2.(3.96b) The resolution of the homogeneous equation is analysed in section 3.7. If a non constant distributed loads is applied, a similar procedure may be developed for obtaining the solution of the homogeneous equation. Stresses calculation The stresses calculation is developed in detail in [95]. With a procedure analogous to the one employed in the straight beam case, once the vector x(s) given in (3.66b) has been calculated, the axial higher order moments Ms(s) are directly obtained as the first n−1 components of x(s). Shear higher-order moments and transverse higher-order moments may be obtained from the equilibrium equations (3.87) and (3.88) respectively. Strains may be obtained by using the compliance equations (3.63). Once strains are obtained, circumferential stress can be calculated in each ply by using the 2D constitutive equations (3.60). Finally, the shear and through-thickness stresses are obtained by integration of the equilibrium equations (3.83). In that way, stresses can be written as: σp s(s, z) = Sp s(z)Tˇ KσMT(s),(3.97a) τp sz(s, z) = τp sz(s, z− p)R±z− p R±z2 −Sp sz(z)Tˇ Kσ dMT(s) ds ,(3.97b) σp z(s, z) = R±z− p R±zσp z(s, z− p)±Rdτp sz(s, z− p) ds R±z− p R±z−1 ±Sp z,1(z)Tˇ KσMT(s) + Sp z,2(z)Tˇ Kσ d2MT(s) ds2,(3.97c) where the vector MT(s) is given in (3.73) and the shape vectors Sp s(z) (given by its components Sp si(z) for i= 0,1, ..., 2n−1), Sp sz(z), Sp z,1(z) and Sp z,2(z) are given by: Sp si(z) = R R±zQp 11zi, i = 0,1, ..., n, (3.98a) Sp si+n−1(z) = R±z RQp 13zi−2, i = 2,3, ..., n, (3.98b) 3.5. Legendre polynomials based model (LPBM) 96 where fs(z) is given in equation (3.107a) and it can be demonstrated, following the same procedure employed in the straight beam case of the MBM, that fz(z) = ˆ fs(z), given in equation (3.107b). Strains By using the strain-displacement relations in Cartesian coordinates (3.42) and the Legendre polynomials properties (3.110), strains can be written as: εs(s, z) = eN(s) + zeM(s) + fsT(z)es(s),(3.115a) γsz(s, z) = eQ(s) + ˆ fs T(z)eγ(s),(3.115b) εz(s, z) = ˆ ˆ fs T(z)ez(s),(3.115c) where eN(s), eM(s) and eQ(s) are the non-dimensional CBT 1D strain components, and es(s), eγ(s) and ez(s) are the non-dimensional higher-order 1D strain components, which are defined by: eN(s) = du0(s) ds ,(3.116a) eM(s) = du1(s) ds ,(3.116b) eQ(s) = u1(s) + dw1(s) ds +hQTu(s),(3.116c) es(s) = du(s) ds ,(3.116d) eγ(s) = dw(s) ds +HγTu(s),(3.116e) ez(s) = HzTw(s),(3.116f) where for the present case of a straight beam hQ=ρ. Notice that expressions (3.115) are the same than in the MBM but considering different matrices Hγand Hzand considering a null vector ρin the MBM case. Equilibrium The non-dimensional higher-order moments in the LPBM are defined as follows: N(s) = Z1 −1 σs(s, z)dz , M(s) = Z1 −1 σs(s, z)z dz , (3.117a) 97 Bi-dimensional models for evaluating interlaminar stresses Q(s) = Z1 −1 τsz(s, z)dz , (3.117b) Ms(s) = Z1 −1 σs(s, z)fs(z)dz , (3.117c) Mτ(s) = Z1 −1 τsz(s, z)fτ(z)dz , (3.117d) Mz(s) = Z1 −1 σz(s, z)fσ(z)dz . (3.117e) where it can be demonstrated with a similar procedure to that employed in the MBM that fτ(z) = ˆ fs(z) and fσ(z) = ˆ ˆ fs(z). Therefore, integrating the equilibrium equations (3.51), the forces and bending moment equilibrium equations yield as given in (3.53) and (3.54) considering the non-dimensional thickness t= 2, those equations being independent of the kind of model. The higher-order moments equilibrium equations yield: dMs(s) ds +fs(1) τ+(s)−fs(−1) τ−(s) = HγMτ(s) + hQQ(s),(3.118a) dMτ(s) ds +ˆ fs(1) σ+(s)−ˆ fs(−1) σ−(s) = HzMz(s).(3.118b) Constitutive law As in the straight beam case of the MBM, the 2D constitutive law for an individual ply is given by equations (3.60). Substituting in these equations the strains, given by (3.115), and integrating the equations in order to obtain the forces and moments (including the higher-order ones), relations (3.61) are obtained. However, in this case matrices ˆ Kσand ˆ Kτare given by equations (3.62) using the functions fs(z), ˆ fs(z) and ˆ ˆ fs(z) corresponding to the LPBM. Some integrals of these functions are solved in the Appendices of [96]. The compliance matrices can be defined by inverting the stiffness matrices, ˇ Kσ=ˆ Kσ−1and ˇ Kτ=ˆ Kτ−1, so that decomposing the compliance matrices in several components as seen in [96], compliance equations yield (3.63). 3.5. Legendre polynomials based model (LPBM) 98 Solution procedure The solution procedure is similar to the one developed in the MBM. Introducing (3.63b) into (3.116e), longitudinal higher-order displacements yield: HγTu(s) = ˇ Cττ Mτ(s) + ˇcQτ Q(s)−dw(s) ds ,(3.119) Substituting (3.63a) into (3.116d) and (3.116f), using (3.119) to eliminate the derivative of the longitudinal displacement u(s), using (3.118a) to eliminate the higher-order shear moments Mτ(s), and using equation (3.53b) to eliminate the derivative of the shear force Q(s), the following differential equation may be obtained (see [96] for more details): d2x(s) ds2=Gx(s) + gNN(s) + gMM(s) +gσ+σ+(s) + gσ−σ−(s) + gτ+ dτ+(s) ds +gτ− dτ−(s) ds ,(3.120) where the auxiliary matrices employed in the previous equation are given in Appendix A.3 and vector x(s) is defined in (3.66b). If a non constant distributed loads is applied, a similar procedure may be developed for obtaining the solution of the homogeneous equation. Constant distributed loads case Considering the straight beam loaded under constant distributed loads σ+(s) = σ+ 0,σ−(s) = σ− 0,τ+(s) = τ+ 0and τ−(s) = τ− 0as a particular case, using equilibrium equations (3.53) and (3.54) and using the boundary conditions (3.55), the axial force, shear force and bending moment yield (3.67). Therefore, as the forces and the bending moment are the same as in the straight beam of the MBM, and the differential equation (3.120) is the same as in the MBM, using appropriate values of the auxiliary matrices (Hγ, Hz,. . .), the solution of this problem is the same than in the MBM, given by equations (3.68)-(3.71). The resolution of the homogeneous equation is analysed in section 3.7. Stresses calculation With an analogous development to that carried out in the straight beam case of the MBM, once the vector x(s) given in (3.66b) has been calculated, 99 Bi-dimensional models for evaluating interlaminar stresses the axial higher order moments Ms(s) are directly obtained as the first n−1 components of x(s). Shear higher-order moments and transverse higherorder moments may be obtained from the equilibrium equations (3.118a) and (3.118b) respectively. Strains may be obtained by using the compliance equations (3.63). Once strains are obtained, circumferential stress can be calculated in each ply by using the 2D constitutive equations (3.60). Finally, the shear and through-thickness stresses are obtained by integration of the equilibrium equations (3.51). In that way, stresses yield: σp θ(s, z) = Sp s(z)Tˇ KσMT(s),(3.121a) τp sz(s, z) = τp sz(s, z− p)−Sp sz(z)Tˇ Kσ dMT(s) ds ,(3.121b) σp z(s, z) = σp z(s, z− p)−dτp sz(s, z− p) ds z−z− p+Sp z,2(z)Tˇ Kσ d2MT(s) ds2, (3.121c) where the vector MT(s) is given in (3.73) and the shape vectors Sp s(z) (given by its components Sp si(z) for i= 0,1, ..., 2n−1), Sp sz(z) and Sp z,2(z) are given by: Sp si(z) = Qp 11pi(z), i = 0,1, ..., n, (3.122a) Sp si+n−1(z) = Qp 13pi−2(z), i = 2,3, ..., n, (3.122b) Sp sz(z) = Zz z− p Sp s(z0)dz0,(3.122c) Sp z,2(z) = Zz z− p Sp sz(z0)dz0.(3.122d) Shear and through-thickness stresses have to be calculated iteratively as explained in the MBM, beginning with the boundary conditions τ1 sz(s, −1) = τ−(s) and σ1 z(s, −1) = σ−(s). The boundary conditions τNp sz (s, 1) = τ+(s) and σNp z(s, 1) = σ+(s) are automatically accomplished. 3.5.3 Development of the LPBM for curved beams A curved beam with thickness t(considered non-dimensional with value t= 2) and length Lk=RkΘk(where Rkis the mean radius and Θkis the angular length of the beam) is studied (see Figure 3.7). A local nondimensional coordinate system skand zkis used, defined as in the MBM. Notice that the difference of the clockwise and the counter-clockwise curved beam lies in the direction of the zkaxis. 3.5. Legendre polynomials based model (LPBM) 100 The beam is loaded under distributed transverse and shear loads in the positive boundary of zk,σ+(sk) and τ+(sk), and under distributed transverse and shear loads in the negative boundary of zk,σ−(sk) and τ−(sk). Arbitrary stresses boundary conditions are applied at sk= 0, given by σ0 s(zk) and τ0 s(zk), and at sk=Lk, given by σL s(zk) and τL s(zk). For the sake of simplicity, the index kindicating the component is not used in the model development. The model development for a curved beam in the LPBM is detailed in [96]. Consequently, only the main equations of the model are shown below. Displacements The displacements in matrix notation are expressed by a series expansion according to equations (3.114). In the LPBM for a curved beam, functions fz(z) yield: fz(z) = ˆ fs(z)±1 RΞ2fs(z) + Ξ1ˆ ˆ fs(z),(3.123) which accomplish property (3.77) and can be expressed also as follows: fz(z) = ±R(Hsfs(z) + hMz+hN),(3.124) where matrix Hsand vectors hNand hMare given by: Hs=±1 RΥ+1 R2(Ξ2+Ξ1ΥΥ),(3.125a) hM=±1 Rυ+1 R2Ξ1Υυ,(3.125b) hN=1 R2Ξ1υ.(3.125c) Strains Strains are given by: εs(s, z) = R R±zeN(s) + zeM(s) + fs(z)Tes(s),(3.126a) γsz(s, z) = R R±zeQ(s) + fz(z)Teγ(s),(3.126b) εz(s, z) = R±z R ˆ ˆ fs(z)Tez(s),(3.126c) 101 Bi-dimensional models for evaluating interlaminar stresses the components of the strains being defined as: eN(s) = du0(s) ds ±w1(s) R+hNTw(s),(3.127a) eM(s) = du1(s) ds +hMTw(s),(3.127b) eQ(s) = u1(s)∓u0(s) R+dw1(s) ds +hQTu(s),(3.127c) es(s) = du(s) ds +HsTw(s),(3.127d) eγ(s) = dw(s) ds +HγTu(s),(3.127e) ez(s) = HzTw(s),(3.127f) where auxiliary vector hQis given by: hQ=ρ∓1 RHγΞ2υ.(3.128) Equilibrium The forces, bending moment and higher-order moments are defined, as in the straight beam, by equations (3.117), by using functions fτ(z) and fσ(z) given by equations (3.84) with the corresponding functions fz(z) and ˆ ˆ fs(z) to the LPBM. The equilibrium equations for the forces and the bending moment remain the same than in the MBM, given by (3.85) with a non-dimensional thickness of t= 2. The equilibrium equations for the higher-order moments yield: dMs(s) ds +R+ Rτ+(s)fs(1) −R− Rτ−(s)fs(−1) = HγMτ(s) + hQQ(s), (3.129a) dMτ(s) ds +R+ Rσ+(s)fz(1) −R− Rσ−(s)fz(−1) =HzMz(s) + HsMs(s) + hMM(s) + hNN(s).(3.129b) Constitutive law As in the straight beam case of the MBM, the 2D constitutive law for an individual ply is given by equations (3.60). Substituting in those equations 3.5. Legendre polynomials based model (LPBM) 102 the strains, given by (3.115), and integrating the equations in order to obtain the forces and moments (including the higher-order ones), relations (3.61) are obtained. However, in this case components of matrices ˆ Kσand ˆ Kτ are different to those of the straight beam case due to the influence of the curvature. Expressions of these matrices can be found in [96]. The compliance matrices can be defined by inverting the stiffness matrices, ˇ Kσ=ˆ Kσ−1and ˇ Kτ=ˆ Kτ−1, so that decomposing the compliance matrices in several components as seen in [96], compliance equations yield (3.63). Solution procedure Introducing (3.63b) into (3.127e) longitudinal higher-order displacements yield: HγTu(s) = ˇ Cττ Mτ(s) + ˇcQτ Q(s)−dw(s) ds ,(3.130) Substituting (3.63a) into (3.127d) and (3.127f), using (3.130) to eliminate the derivative of the longitudinal displacement u(s), using (3.129a) to eliminate the higher-order shear moments Mτ(s), and using equation (3.85b) to eliminate the derivative of the shear force Q(s), the following differential equation is obtained: d2x(s) ds2=Gx(s) + gNN(s) + gMM(s) +gσ+σ+(s) + gσ−σ−(s) + gτ+ dτ+(s) ds +gτ− dτ−(s) ds ,(3.131) where the matrices of the previous equation are given in Appendix A.4 and vector x(s) is defined in (3.66b). Constant distributed loads case Considering the curved beam loaded under constant distributed loads σ+(s) = σ+ 0,σ−(s) = σ− 0,τ+(s) = τ+ 0and τ−(s) = τ− 0as a particular case, using equilibrium equations (3.85) and using the boundary conditions (3.55), the axial force, shear force and bending moment yield (3.92). Therefore, as the forces and the bending moment are the same than in the curved beam of the MBM, and the differential equation (3.131) is analogous to that in the MBM, the solution of this problem can be written as in the MBM, given by equations (3.93)-(3.96), but using the appropriate values of the auxiliary matrices (Hγ,Hz,. . .). 103 Bi-dimensional models for evaluating interlaminar stresses The resolution of the homogeneous equation is analysed in section 3.7. If a not constant distributed loads is applied, a similar procedure may be developed for obtaining the solution of the homogeneous equation. Stresses calculation Stresses are obtained from the higher-order moments as follows: σp s(s, z) = Sp s(z)Tˇ KσMT(s),(3.132a) τp sz(s, z) = τp sz(s, z− p)R±z− p R±z2 −Sp sz(z)Tˇ Kσ dMT(s) ds ,(3.132b) σp z(s, z) = R±z− p R±zσp z(s, z− p)±Rdτp sz(s, z− p) ds R±z− p R±z−1 ±Sp z,1(z)Tˇ KσMT(s) + Sp z,2(z)Tˇ Kσ d2MT(s) ds2,(3.132c) where the vector MT(s) is given in (3.73) and the shape vectors Sp s(z) (given by its components Sp si(z) for i= 0,1, ..., 2n−1), Sp sz(z), Sp z,1(z) and Sp z,2(z) are given by: Sp si(z) = R R±zQp 11pi(z), i = 0,1, ..., n, (3.133a) Sp si+n−1(z) = R±z RQp 13pi−2(z), i = 2,3, ..., n, (3.133b) Sp sz(z) = R (R±z)2Zz z− p (R±z0)Sp s(z0)dz0,(3.133c) Sp z,1(z) = 1 R±zZz z− p Sp s(z0)dz0,(3.133d) Sp z,2(z) = R R±zZz z− p Sp sz(z0)dz0.(3.133e) Shear and through-thickness stresses have to be calculated iteratively as explained in the MBM, beginning with the boundary conditions τ1 sz(s, −1) = τ−(s) and σ1 z(s, −1) = σ−(s). The boundary conditions τNp sz (s, 1) = τ+(s) and σNp z(s, 1) = σ+(s) are automatically accomplished. 3.6. Resolution of the regularized problem by using the non-regularized models 104 3.6 Resolution of the regularized problem by using the non-regularized models Both non-regularized models, MBM and LPBM, have been developed following the same procedure. Therefore, both models are reduced to a set of differential equations that can be written in the same manner, using different constant vectors and matrices. In this way, e.g., the regularized values of vector x(s), defined in (3.66b), for constant distributed loads case are given by equation (3.68) for the straight beam in both MBM and LPBM, and (3.93) for the curved beam in both MBM and LPBM. If another kind of distributed load is given the particular solution have to be obtained in a similar way. As described previously, once the regularized solution xreg(s) has been obtained, regularized higher order moments, displacements, strains and stresses can be also obtained. In the present section, results of the aforementioned regularized solutions are compared for several orders with the regularized solution given by Ko & Jackson [4] in the case of null distributed loads, showing the convergence of the models when the order increases. Further comparisons may be found in [95,108,109]. The problem with constant distributed loads cannot be solved with Ko & Jackson’s solution, then, it will be compared with finite elements solutions (for the non-regularized case) in a following section. Both models, MBM and LPBM, are developed by a polynomial series expansion until an order n, and consequently, both models have almost the same accuracy for the same order, except when the MBM is used with any order above its limitation, given in Table 3.1. However, as it will be shown later, LPBM requires higher computational times, therefore for low orders (under the limit seen in Table 3.1) it is preferable to use the MBM and for higher orders it is preferable to use the LPBM. Hence, the results presented in this and the following sections are valid for both models while the order is under the limit shown in Table 3.1, and only for the LPBM for higher orders. The model is applied to a L section forming 90obetween the straight arms with a radius to thickness ratio of R/t = 1.5 at the corner. The material selected is a conventional UD-CFRP with a ply thickness of 0.2 mm and the ply properties summarized in Table 3.2. The stacking sequence chosen is a quasi-isotropic lay-up of 16 plies, defined by [45,0,-45,90]2S. The regularized solution depends only on the forces, bending moment and distributed loads. In the present section, stresses due to the forces and the bending moment will be analysed and compared with Ko & Jackson’s 105 Bi-dimensional models for evaluating interlaminar stresses E11 150 GPa E22 10 GPa E33 10 GPa ν12 0.3ν13 0.3ν23 0.3 G12 6 GPa G13 6 GPa G23 3.85 GPa Table 3.2: Ply properties for the numerical examples solution. First, stresses due to an unitary bending moment, M0= 1 N, are depicted in Figure 3.10 by applying the non-regularized model for n= 1, 2, 4 and 8. Figure 3.10: Regularized solution of the non-regularized models for the bending moment loading. Notice that the shear stresses due to the bending moment are null, and they have not been depicted in Figure 3.10. Results of the model are very accurate for the regularized stresses due to the bending moment even for low orders, showing a fast convergence towards Ko and Jackson’s solution. Stresses due to an unitary axial force, N0= 1 N/mm, are depicted in Figure 3.11 by applying the non-regularized model for n= 1, 3, 6 and 12. 3.8. Resolution of the joints between components 112 The matrix KMis calculated by using least squares in both sides of equation (3.148): ∂ ∂KM Z1 −1 1 fz(z)−KM1 ˆ fs(z)2 dz!= 0,(3.149) where the derivative yields: 2Z1 −1 1 fz(z)−KM1 ˆ fs(z)1ˆ fs(z)Tdz = 0.(3.150) Considering that KMdoes not depend on z, it can be expressed as follows: KM=BMAM−1,(3.151a) AM=Z1 −11 ˆ fs(z)1ˆ fs(z)dz, (3.151b) BM=Z1 −11 fz(z)1ˆ fs(z)dz. (3.151c) Therefore, continuity of shear stresses is imposed in a weak sense equalling the higher-order moments expressed in the common basis: ¯ Mk−1 τ(L) = ¯ Mk τ(0).(3.152) where ¯ Mk−1 τ(L) and ¯ Mk τ(0) are obtained using (3.148). 3.8.3 Continuity of the circumferential displacements between components To impose continuity, circumferential higher-order displacement vector have to be equalled between both components: uk−1(L) = uk(0). Using equation (3.65) the equality can be written as: −dwk−1 ds (L) + ˇ Ck−1 ττ Mk−1 τ(L) + ˇck−1 Qτ Qk−1(L) =−dwk ds (0) + ˇ Ck ττ Mk τ(0) + ˇck Qτ Qk(0).(3.153) Notice that matrices Hγdo not depend on the typology of the component, so it has been cancelled at both sides of the previous equation. 113 Bi-dimensional models for evaluating interlaminar stresses Considering, for the sake of simplicity, the constant distributed loads case, the derivative of w(s) can be expressed depending on Mτ(s) and Ms(s) by using the derivative of equation (3.140), equation (3.144), equation (3.95) and considering that the non-regularized value is the regularized value plus the perturbation, what results in: dwk(sk) ds =hdXk z0(sk) ds dXk zL(sk) ds iΓk 0 Γk L       Mk s(0) Mk τ(0) Mk s(L) Mk τ(L)    −    Mk s,reg(0) Mk τ,reg(0) Mk s,reg(L) Mk τ,reg(L)        +qk zQQk(sk) + qk zτ+τk 0+(sk) + qk zτ−τk 0−(sk),(3.154) where qk zQ,qk zτ+and qk zτ−are respectively the last n−1 components of vectors qQ,qτ+and qτ−in the component k. Considering equations (3.129a), (3.93) and (3.95), vectors Mk s,reg(sk) and Mk τ,reg(sk) can be written as: Mk s,reg(sk) Mk τ,reg(sk)=KF  Nk(sk) Mk(sk) Qk(sk) +Kσ    σk 0+ σk 0− τk 0+ τk 0−     ,(3.155a) KF=qk sN qk sM o o o Hγ−1(qk sQ −hk Q),(3.155b) Kσ=qk sσ+qk sσ−o o o o Hγ−1(qk sτ++R+ Rfs(1)) Hγ−1(qk sτ−−R− Rfs(−1)), (3.155c) where the constants σk 0+,σk 0−,τk 0+ and τk 0−are the constant distributed loads σ+ 0,σ− 0,τ+ 0and τ− 0respectively in the component k. Substituting (3.154), (3.155a) and (3.148) in (3.153) the continuity condition of the circumferential displacements is obtained: 3.8. Resolution of the joints between components 114 Kk−1 sL     Mk−1 s(0) ¯ Mk−1 τ(0) Mk−1 s(L) ¯ Mk−1 τ(L)     +ˆ Kk−1 sL         Nk−1(0) Mk−1(0) Qk−1(0) Nk−1(L) Mk−1(L) Qk−1(L)         +ˆ ˆ Kk−1 sL     σk−1 0+ σk−1 0− τk−1 0+ τk−1 0−     =Kk s0    Mk s(0) ¯ Mk τ(0) Mk s(L) ¯ Mk τ(L)     +ˆ Kk s0         Nk(0) Mk(0) Qk(0) Nk(L) Mk(L) Qk(L)         +ˆ ˆ Kk s0    σk 0+ σk 0− τk 0+ τk 0−     .(3.156) If one of the two adjacent beams is considered semi-infinite it affects to the corresponding side of the equal sign in equation (3.156). This is due to the dependence of equation (3.156) with the value of the higher-order moments at both ends sk= 0 and sk=Lk. 3.8.4 Continuity of the transverse displacements between components As in the case of the shear moments, it is recommended to express the transverse displacement vector in the straight beam basis for imposing continuity of transverse displacements. The transverse displacements expressed in the straight beam basis and in the component own basis are given by: uz(s, z) = ¯w1(s) + ˆ fs T(z)¯w(s), uz(s, z) = w1(s) + fzT(z)w(s).(3.157) The continuity of displacements is imposed in a weak manner by using least squares. In the theoretical model w1(s) has been calculated independently of the rest of parameters, w1(s) being calculated once the rest of the problem has been solved, and, therefore, it is not considered in the least squares problem. ∂ ∂¯wZ1 −1ˆ fs T(z)¯w(s)−fzT(z)w(s)2dz= 0.(3.158) The derivative yield: 115 Bi-dimensional models for evaluating interlaminar stresses Z1 −1 ˆ fs(z)ˆ fs T(z)¯w(s)−fzT(z)w(s)dz = 0.(3.159) Therefore, Kwis obtained as following the same procedure used in the shear higher-order moments case: ¯w(s) = Kww(s),Kw=Aw−1Bw(3.160a) Aw=Z1 −1 ˆ fs(z)ˆ fs T(z)dz, Bw=Z1 −1 ˆ fs(z)fzT(z)dz. (3.160b) Aw,ij =ˆ AL,i+1,j+1, i, j = 1,2, ..., n −1.(3.160c) Bw,ij =ˆ BL,j+1,i+1, i, j = 1,2, ..., n −1.(3.160d) Equalling ¯wk−1(L) = ¯wk(0), substituting equations (3.140), (3.144) and (3.93) and knowing that Xk z0(L) = 0 and Xk zL(0) = 0, the continuity of the transversal displacements is imposed with: Kk−1 w    Xk−1 zL (L)Γk−1 L    ˆ Mk−1 s(0) ˆ Mk−1 τ(0) ˆ Mk−1 s(L) ˆ Mk−1 τ(L)     +qk−1 zN qk−1 zM Nk−1(L) Mk−1(L)    =Kk w    Xk z0(0)Γk 0    ˆ Mk s(0) ˆ Mk τ(0) ˆ Mk s(L) ˆ Mk τ(L)     +qk zN qk zM Nk(0) Mk(0)    .(3.161) Separating the perturbations as the non-regularized value minus the regularized one, and using equations (3.155a) for the regularized values, the previous equation can be expressed as follows: 3.9. Results of the non-regularized models and validation by FEM 116 Kk−1 zL     Mk−1 s(0) ¯ Mk−1 τ(0) Mk−1 s(L) ¯ Mk−1 τ(L)     +ˆ Kk−1 zL         Nk−1(0) Mk−1(0) Qk−1(0) Nk−1(L) Mk−1(L) Qk−1(L)         +ˆ ˆ Kk−1 zL     σk−1 0+ σk−1 0− τk−1 0+ τk−1 0−     =Kk z0    Mk s(0) ¯ Mk τ(0) Mk s(L) ¯ Mk τ(L)     +ˆ Kk z0         Nk(0) Mk(0) Qk(0) Nk(L) Mk(L) Qk(L)         +ˆ ˆ Kk z0    σk 0+ σk 0− τk 0+ τk 0−     .(3.162) If one of the two adjacent beams is considered semi-infinite it affects to the corresponding side of the equal sign in equation (3.162). 3.8.5 Application of the boundary conditions and the continuity conditions Each component has 4n−4 unknown variables given by the higher-order moments (radial and circumferential) at each end of the component, so (4n− 4)Ncunknown variables are given. Continuity conditions of stresses, (3.146) and (3.152), allows the unknown variables to be reduced to (2n−2)(Nc+1). Boundary conditions at both ends of the chain of beams, which can be given in displacements or in stresses, provide 4n−4 equations, so (2n− 2)(Nc−1) additional equations are required. If one end component is considered semi-infinite, the number of unknowns in the component is reduced to 2n−2 unknowns. Therefore, if semi-infinite components are assumed at both ends (2n−2)(Nc−1) additional equations are also required. These (2n−2)(Nc−1) equations are given by the continuity conditions of the displacements, (3.156) and (3.162). Notice that all the boundary conditions and continuity conditions are linear with the unknown variables, and, consequently, they can be solved easily as a problem of the form of Ax =b, where bdepends on the applied forces, bending moment and distributed loads. In the present case, vector xincludes all the circumferential higher-order moments Mk s(0) and/or Mk s(L), and all the shear higher-order moments expressed in the functions base of the straight beam ¯ Mk τ(0) and/or ¯ Mk τ(L). Notice also that the forces and the bending moment at the end of one component are the same than these at the beginning of the following component. 117 Bi-dimensional models for evaluating interlaminar stresses 3.9 Results of the non-regularized models and validation by FEM In the present section several geometries and loading states are analysed and compared with appropriate FE models to validate the semi-analytical model. First, a L-sectioned beam is analysed with two different loading states: a bending moment loading, similarly to a four-point bending test and a compression loading in one arm of the section. Finally, a joggle is analysed under tensile loading. The models can be also applied to other geometries including changes of curvature with constant thickness, such as Ω-shaped beams, C-shaped beams or corrugated laminates. 3.9.1 L-shaped beam under bending moment The typical loading state associated to the unfolding failure is the bending moment of a corner. A L-shaped section loaded under a pure bending moment calculated with a regularized model has a constant stress distribution in the circumferential direction due to the constant value of the bending moment in that direction. However, the change of curvature due to the straight arms cause a non-constant stress distribution in the zone where the change of curvature takes place (see [14]). For the numerical implementation, a material with ply properties of Table 3.2 is considered. The laminate considered is a quasi-isotropic laminate with stacking sequence [45 0 -45 90]S. The L-shaped section is determined by the thickness t= 1.6 mm and by the mean radius of the curved part R= 2.4 mm, see Figure 3.13. The applied load is M0= 1 Nm/m. Figure 3.13: L-shaped beam under bending moment. In the FE model, the straight arms are taken of length L= 10 mm, which has been considered long enough to avoid that the non-regularized effects of the end of the arms do no affect to the stresses distribution in the zone near 3.9. Results of the non-regularized models and validation by FEM 118 the change of curvature. Linear shell elements with four nodes have been used in the FE model with five elements per ply in the thickness direction. The ply properties are introduced considering a plane strain state. The hypotheses used for modelling the FE model are similar to the hypotheses employed in the semi-analytical model, and, consequently, similar solutions must be obtained. The results of the maximum INS depending on the axial coordinate s by using the LPBM with an order of n= 50, Ko & Jackson’s solution and the FEM results are depicted in Figure 3.14. The same solutions for the maximum/minimum ISS are depicted in Figure 3.15, where the Ko & Jackson’s model solution has been suppressed as the regularized shear stress in the problem is null in the whole section. Figure 3.14: Maximum INS in the L-shaped beam under bending moment. Notice that the FE model agrees well with the results of the LPBM, clearly showing the existence of the non-regularized effects due to the change of curvature and validating the non-regularized models developed in the present project. The shear stresses appear even without the presence of a shear force due to the non-regularized effects in order to accomplish the equilibrium equations. These shear stresses present their maximum value in the change of curvature point, with a discontinuity in their derivative in the aforementioned section. Notice that in the view of the obtained results, non-regularized models are not necessary to predict the unfolding failure in this case, as the maxi- 119 Bi-dimensional models for evaluating interlaminar stresses Figure 3.15: Maximum ISS in the L-shaped beam under bending moment. mum value of the INS is predicted accurately by the regularized model given by Ko & Jackson. The reason is that the maximum regularized INS is constant in the whole curved beam. Notwithstanding, in a problem where the maximum INS is given near to a change of curvature non-regularized effects are more relevant, as in the case of a L-shaped beam under a compressive load shown in the next section. 3.9.2 L-shaped beam under a compressive load A compressive loading in an arm of a L-shaped beam induces that the maximum bending moment is given at the change of curvature of the opposite arm. Therefore, the maximum INS is affected by the non-regularized effects and the use of a regularized model is not accurate in this kind of loading states. For the numerical implementation, the same material, configuration and FE model characteristics are used, since the loading state is changed as depicted in Figure 3.16. The results of the maximum INS depending on the axial coordinate s by using the LPBM with an order of n= 50, Ko & Jackson’s solution and the FE model predictions are depicted in Figure 3.17. The results of the maximum/minimum ISS are depicted in Figure 3.18, where the maximum value of the Ko & Jackson’s solution has been suppressed as the maximum value of the regularized shear stress in the problem is null in the whole beam. 3.9. Results of the non-regularized models and validation by FEM 120 Figure 3.16: L-shaped beam under a compressive load. Figure 3.17: Maximum INS in the L-shaped beam under compressive load. Notice that, as the bending problem analysed in the previous subsection, the FE model agrees well with the results of the LPBM. The maximum value of the INS in the present loading case is a 64% of the maximum value given by a regularized model. Consequently, regularized models are very conservative when the maximum INS is given near to a change of curvature. These results show the importance of considering the non-regularized effects in the calculation of the unfolding failure. The regularized value of the ISS is very close to the non-regularized value, with small discrepancies near to the changes of curvature. In this case, the non-regularized effects for the ISS are not so important due to the 121 Bi-dimensional models for evaluating interlaminar stresses Figure 3.18: Maximum ISS in the L-shaped beam under compressive load. presence of a dominant shear force. 3.9.3 Joggle under a tensile load The joggle is a very common element in riveted or bolted joints between panels where the continuity of the surface is important from an aerodynamic point of view. This element has two curved parts, commonly having the same mean radius R. One of this curved parts is clockwise and the other one is counter-clockwise. Consequently, the joint between them constitutes a change of curvature even if they both have the same mean radius. However, a short straight part is usually located between both curved parts, see Figure 3.19. Figure 3.19: Joggle under a tensile load. The importance of the non-regularized model in a joggle resides in the shortness of the curved elements and the high curvature, which implies that 3.12. Off-topic applications of the non regularized models 128 Figure 3.24: Complex plane representation of √λmin. Figure 3.25: Regularization distances for a single-ply laminate depending on the orientation. between a value of 0.65tand 2.25t. In the case of a composite laminate these values may change and it is necessary to evaluate each stacking sequence independently. 129 Bi-dimensional models for evaluating interlaminar stresses 3.12 Off-topic applications of the non regularized models The MBM and LPBM have been applied in the previous sections to composite sections loaded by end forces and bending moments including changes of curvature, in order to predict the effect of these changes of curvature on the stresses distribution. However, these models may be applied to many others applications [110]. In particular, in the present section, the models are applied to evaluate stresses at joints between beams with different material properties, to determine non-regularized effects at the vicinity of an embedment and to analyse the three point bending test. 3.12.1 Joint of beams with different material properties In a joint of two materials a discontinuity in the material properties takes place, which, due to the different Poisson contractions of the materials, causes the appearance of a singular stress field. This kind of problem may be analysed using the models developed in section 3.8 by considering two components and introducing different material properties in each component. As a numerical example, two homogeneous isotropic materials are considered with different Young modulus and the same Poisson coefficients, with the material properties shown in Table 3.3. Material 1 E1100 GPa ν10.3 Material 2 E210 GPa ν20.3 Table 3.3: Material properties for the numerical examples The two beams made with different materials are considered perfectly bonded at a flat surface as shown in Figure 3.26. Both beams are of constant thickness t. Furthermore, a tensile load N0is acting perpendicular to the joint surface. Figure 3.26: Joint of two beams with different material properties, axially loaded. 3.12. Off-topic applications of the non regularized models 130 For comparison reasons, the problem has been implemented in a FE software. Quadratic bi-dimensional elements with a plane strain assumption have been used, where a variable mesh size has been used concentrating the elements near to the free edges and near to the joint section. Figure 3.27 shows the distribution of the axial stress σsat the joint section obtained by using the LPBM with n= 500, which has required a computational time of 50 seconds, and compared with the FE results. Figure 3.27: Axial stress at the joint of two different materials. The LPBM, as the FEM or any model based on series expansion with polynomials, cannot predict an infinite value of the stresses at the singular point. However, the tendency is clearly captured, observing both tendencies to an infinite value in the free edges. The LPBM obtains a similar distribution to the FEM. Far away from the joint section the distribution of stresses tends to a constant value, it being the regularized value. Figure 3.28(a) shows the shear stress τsz and Figure 3.28(b) shows the through-thickness stress σzin the joint section obtained by using the LPBM with n= 500 and compared with the FE results. The through-thickness stress is different at both sides, as this stress is not continuous in the axial direction s, so the stress at both sides, corresponding to both materials, is depicted. The stiffer material is tensioned in the through-thickness direction, while the softer material is compressed. For the present material properties, the through-thickness stress is in the same order than the mean value of the axial stress. The normal through-thickness stress is not singular in the free edges. Notwithstanding, it has an infinite slope changing sharply from 0 to a finite value. The shear stress is singular similarly than the axial stress, 131 Bi-dimensional models for evaluating interlaminar stresses (a) (b) Figure 3.28: Stresses in the joint of two different materials. (a) Shear stress. (b) Through-thickness stress in both sides. however, it is neither captured by the LPBM nor by the FE model. A high agreement is observed between the LPBM and the FEM results. 3.12.2 Non-regularized effects due to the boundary conditions The way in which the load is applied may introduce also non-regularized effects. According to the Saint-Venant principle, at sufficiently large distances from load the way of application of two equivalent loads do not affect to the stress distribution. However, the way in which the load is applied affects to the nearest zone in the form of non-regularized effects. If boundary conditions are prescribed stresses, this kind of non-regularized effects may be calculated. Analogously, prescribed displacements may introduce non-regularized effects over the stress distribution. As a numerical example, a straight semi-infinite beam embedded, that is, with null transverse and longitudinal displacements at the end, is considered. It is loaded under a tensile load N0, and it is constituted by a homogeneous isotropic material with E= 100 GPa and ν= 0.3. For comparison reasons, the problem has been implemented in a FEM software. Quadratic bi-dimensional elements with a plane strain assumption have been used, where a variable mesh size has been used concentrating the elements near to the free edge in the thickness direction and near to the constrained section. Figure 3.29 shows the axial stress σsin the embedment section obtained by using the LPBM with n= 500, which has required a computational time 3.12. Off-topic applications of the non regularized models 132 of 25 seconds, and by using a FE model. Figure 3.29: Axial stress in an embedment under tensile load. The stress distribution obtained in Figure 3.29 is very similar to the stress distribution obtained in Figure 3.27 as the embedment is analogous to a joint with a material with an infinite stiffness. The singular behaviour of the stress is clearly captured with the model, and results have a high agreement with the FE results. Figure 3.30(a) shows the shear stress τsz and Figure 3.30(b) shows the through-thickness stress σzin the embedment section obtained by using the LPBM with n= 500 and the FEM. (a) (b) Figure 3.30: Stresses in an embedment under tensile load. (a) Shear stress. (b) Through-thickness stress. 133 Bi-dimensional models for evaluating interlaminar stresses Shear and through-thickness stresses have more pronounced distributions than in the case of the joint with a different material. Shear stress and through-thickness normal stress have a singular value at the free edge, and both have a high agreement with FE results. Another numerical example is analysed by using a pure bending moment load M0instead of the tensile load. Figure 3.31 shows the axial stresses σs in the embedment section obtained by using the LPBM with n= 500 with the mentioned bending moment, which has required a computational time of 25 seconds, and by using a FE model. Figure 3.31: Axial stress in an embedment under bending moment. A singular behaviour in the axial stresses distribution is again observed, with axial stresses tending to infinite in the free edges in the embedment section. A high agreement with the FE results is again observed. Figure 3.32(a) shows the shear stresses τsz and Figure 3.32(b) shows the through-thickness stresses σzin the embedment section obtained by using the LPBM with n= 500 and the FEM. In the present case the shear stresses tend to infinite in the free edge. Respect to the through-thickness stress, a tensile zone and a compressive zone is observed in difference with the tensile load case. A high agreement with the FE results is again observed. 3.12.3 Non-regularized effects due to punctual loads The last off-topic application analysed is the effect of punctual loads over the stress distribution. Theotokoglou and Sideridis [111] showed that near to a punctual load in flat laminates the stress distribution is distorted and the 3.12. Off-topic applications of the non regularized models 134 (a) (b) Figure 3.32: Stresses in an embedment under bending moment. (a) Shear stress. (b) Through-thickness stress. typical lineal distribution of axial stresses is not obtained. In particular, a three-point bending test is analysed. The stress distribution in a three-point bending test has been analysed by Makeev et al. [78]. The model developed can predict the effect of a punctual load by considering a external pressure modelled as a Dirac delta and obtaining the regularized expression from equation (3.131). However, the resolution of the regularized value in that case is not easy. As an alternative, by using the developed regularized solution, the punctual load is approximated by a very short straight component with a very high constant pressure whose resultant is the applied load. For the numerical example, the three-point bending test depicted in Figure 3.33 is considered with P= 1 N/mm, L= 2.5 mm and an eight-plies laminate with stacking sequence [45,0,-45,90]2S, a thickness per ply of 0.2 mm and the ply properties of table 3.2. Figure 3.33: Three-point bending test configuration. All loads are approximated as applied in a length δ= 0.05 mm as a pressure with the resultant force given in Figure 3.33. 135 Bi-dimensional models for evaluating interlaminar stresses The numerical example in the non-regularized case is evaluated by using the LPBM with an order of n= 100, which has required a computational time of 30 seconds. The regularized solution of the present problem consists in a linear distribution of axial stresses in the thickness, which increases or decreases linearly with the scoordinate according to the bending moment. The regularized shear stress has a quadratic distribution in the thickness. The axial stresses in the top face of the top ply of Figure 3.33 and in the bottom face of the bottom ply of Figure 3.33 depending on the scoordinate are depicted in Figure 3.34 in parallel with their regularized values. Figure 3.34: Axial stress in the top and bottom of a three-point bending test sample. The maximum of the axial stress in the bottom ply is a 15% higher than the one predicted by a regularized model, and similarly for the top face in the compression zone. This is due to the non-regularized effect of the punctual load, which changes the distribution of the axial stress near to the section of the load application. Far from the load application points, the axial stress tends to its regularized value. The maximum and minimum values of the shear stress depending on the scoordinate is depicted in Figure 3.35. The shear stress presents a singularity at the load application points, tendency which is captured by the non-regularized model. Far away from the load application points, the shear stress tend to its regularized value, which is typically used for the calculation of the interlaminar shear strength 3.12. Off-topic applications of the non regularized models 136 Figure 3.35: Shear stress in a three-point bending test. when this kind of test is used for characterization of the material. Chapter 4 Three-dimensional models for evaluating interlaminar stresses The unfolding failure is typically calculated, from an analytical point of view, by using bi-dimensional (2D) approximations as those presented in Chapter 3, considering in most cases a plane stress or a plane strain approximation, although other generalized plane assumptions may give better results. When a more complex geometry is given, or loads are not constant in the width direction, finite elements models are typically used for calculating the stresses in the curved laminates. By using finite elements models some three-dimensional (3D) effects appearing in curved laminates may be obtained, such as the free-edge effects, the torsion due to the ±45oplies or the anticlastic effect. In the present Chapter, two kinds of 3D models are developed. These 3D models are developed for a singly curved composite laminate loaded under pure bending moment. The pure bending moment case has been chosen because the failure is generally given in zones of the laminate where the bending moment is the predominant load. The first model developed consists in considering the regularized zone of the laminate, far away from the free-edges in the width direction and far away from a change of curvature in the curved direction. Thus, stresses and strains are considered to depend only on the interlaminar coordinate, which reduces considerably the complexity of the problem expressing the partial derivatives as absolute derivatives with respect to a single coordinate. Spencer et al. [97] developed a similar model, only applied to determining the residual stresses due to a homogeneous temperature change in the laminate. The present regularized model is an extension of this model to consider also the loads acting over the laminate. Several authors have developed approximated models in a 3D stresses state for estimating the stresses due 4.1. Regularized 3D model for singly-curved beams under bending moment and residual strains 144 where ¯ Cp ij (with i, j = 1,2, ..., 6) are the stiffnesses of the ply pexpressed in the cylindrical coordinate system, ∆εp θis the residual strain of the ply pin the circumferential direction, ∆εp yis the residual strain of the ply pin the axial direction, ∆εp ris the residual strain of the ply pin the radial direction and ∆γp θy is the residual shear in-plane strain of the ply p. Notice that plies may be differently oriented and, therefore, stresses have been marked with a superscript pindicating in which ply they are calculated. Resolution of the interlaminar shear strains The interlaminar shear strain functions, γry(r) and γrθ(r), can be easily determined by combining the corresponding equilibrium equations with the constitutive law. Hence, substituting the interlaminar shear constitutive equations of (4.10) in the equilibrium equations (4.8a) and (4.8b) they yield: ¯ Cp 44 d(rγry(r)) dr +¯ Cp 45 d(rγrθ(r)) dr = 0,(4.11a) ¯ Cp 45 d(r2γry(r)) dr +¯ Cp 55 d(r2γrθ(r)) dr = 0.(4.11b) The solution of the differential equations system (4.11) yield: γry(r) = R− p ζp 45r¯ Cp 55ηp 1−¯ Cp 45 R− p rηp 2,(4.12a) γrθ(r) = R− p ζp 45r−¯ Cp 45ηp 1+¯ Cp 44 R− p rηp 2,(4.12b) where R− pis the lowest radius of ply p(which constitutes the interface between ply pand ply p−1) and, equivalently, R+ pis the highest radius of ply p.ηp 1and ηp 2are the integration constants in ply p, and ζp 45 is defined by: ζp 45 =¯ Cp 44 ¯ Cp 55 −(¯ Cp 45)2.(4.13) Substituting equations (4.12) into the corresponding constitutive equations of (4.10), interlaminar shear stresses yield: τp ry(r) = R− p rηp 1,(4.14a) τp rθ(r) = R− p r2 ηp 2.(4.14b) Equations (4.14) may also be obtained by integrating directly equations (4.8a) and (4.8b). In view of (4.14), there are a total of 2Npintegration 145 Three-dimensional models for evaluating interlaminar stresses constants, ηp 1and ηp 2, with p= 1,2, . . . , Np, where Npis the total number of plies. An amount of (2Np−2) equations can be obtained by using the equilibrium between plies, which establish the continuity of the interlaminar stresses in the interfaces. The other two equations can be obtained by applying a boundary condition of the interlaminar stresses in the inner or in the outer radius of the laminate. If a boundary condition for one of the through-thickness shear stresses is applied on the inner or outer radius, the other boundary condition in the opposite radius is automatically satisfied due to the equilibrium equations. In the present model, a boundary condition of a null value of both shear interlaminar stresses in the inner radius is considered, which implies a null value of both interlaminar stresses in the outer radius. With these conditions, the integration constants yield ηp 1=ηp 2= 0 for every ply p, and, therefore, interlaminar shear stresses and strains are null in the whole thickness. Resolution of the radial displacements The radial displacements can be determined by the radial equilibrium equation. Substituting the constitutive law (4.10) into the equilibrium equation (4.8c) it yields: ¯ Cp 13 d(rεθ(r)) dr +¯ Cp 23 d(rεy(r)) dr +¯ Cp 33 d(rεr(r)) dr +¯ Cp 36 d(rγθy(r)) dr = ¯ Cp 11εθ(r) + ¯ Cp 12εy(r) + ¯ Cp 13εr(r) + ¯ Cp 16γθy(r) + δp 3−δp 1,(4.15) where: δp i=¯ Cp 1i∆εp θ+¯ Cp 2i∆εp y+¯ Cp 3i∆εp r+¯ Cp 6i∆γp θy, i = 1,2,3,6.(4.16) Notice that ¯ Cp ij =¯ Cp ji. Substituting the strains from equations (4.4), (4.5) and (4.6) into equation (4.15): ¯ Cp 33rd dr rdwo(r) dr −¯ Cp 11wo(r) = hp 1γ0 θy r2 R+hp 2ε0 θ+hp 3ε0 y+δp 3−δp 1r+¯ Cp 16Rγ0 yθ,(4.17) where the auxiliary parameters hp 1,hp 2and hp 3are defined as: 4.1. Regularized 3D model for singly-curved beams under bending moment and residual strains 146 hp 1=¯ Cp 16 −2¯ Cp 36, hp 2=¯ Cp 11 −¯ Cp 13, hp 3=¯ Cp 12 −¯ Cp 23.(4.18) Equation (4.17) can be integrated by, first, obtaining a particular solution and, second, solving the homogeneous equation. Therefore, integrating equation (4.17) the radial displacement yields: w0(r) = −hp 4Rγ0 yθ +rhp 5ε0 θ+hp 6ε0 y+hp 7+ r2 Rhp 8γ0 θy +ξp 1rκp ¯ Cp 13 +κp¯ Cp 33 +ξp 2r−κp ¯ Cp 13 −κp¯ Cp 33 ,(4.19) where ξp 1and ξp 2are two integration constants introduced by the homogeneous equation, which depend on the ply p. In (4.19), the parameter κp, establishing the radial dependence of the radial displacement, is defined as: κp=s¯ Cp 11 ¯ Cp 33 ,(4.20) and the auxiliary parameters are defined as: hp 4=¯ Cp 16 ¯ Cp 11 , hp 5=¯ Cp 11 −¯ Cp 13 ¯ Cp 33 −¯ Cp 11 , hp 6=¯ Cp 12 −¯ Cp 23 ¯ Cp 33 −¯ Cp 11 ,(4.21a) hp 7=δp 3−δp 1 ¯ Cp 33 −¯ Cp 11 , hp 8=¯ Cp 16 −2¯ Cp 36 4¯ Cp 33 −¯ Cp 11 .(4.21b) The particular solution of (4.19) is not valid in plies with κp= 1 or with κp= 2. These cases are considered in a following subsection. The particular solution is not valid when κp= 0 either, but this is not a realistic case, and it is not studied. Notice that κpreminds the Lekhnitskii’s parameter given in equation (3.10a) relating the circumferential and the radial stiffnesses. However, it is not exactly the same, coinciding only when the Lekhnitskii’s equations are applied considering a plane strain assumption. The ply integration constants ξp 1and ξp 2constitute a set of 2Npunknown variables. These variables are determined by imposing the continuity of w1(r) and σr(r) between plies, which implies a total of 2Np−2 equations, and by imposing the boundary conditions of σr(r) in both radial boundaries, R−and R+. These boundary conditions are usually applied with a null 147 Three-dimensional models for evaluating interlaminar stresses value, but they can also have another value, which will cause membrane forces, Nθ, which are defined as: Nθ=ZR+ R− σθ(r)dr. (4.22) Substituting the equilibrium equation (4.8c) into (4.22), the membrane forces yield: Nθ=ZR+ R− d(rσr(r)) dr dr =R+p+−R−p−.(4.23) Once the continuity between plies and the boundary conditions have been imposed, the unknown variables ξp 1and ξp 2are determined depending on the strain constant parameters ε0 θ,ε0 y,γ0 θy and γ0 yθ. These strain parameters define the circumferential strain, the axial strain and torsional strains, which can be fixed as inputs in the problem, or they can be calculated from the applied loads to the laminate. The calculation from applied loads is studied in a following section. In this last case, the unknown variables ξp 1and ξp 2are determined jointly with the strain constant parameters ε0 θ,ε0 y,γ0 θy and γ0 yθ. Strains and stresses calculation In order to apply the equilibrium equations of the radial stress, σr, between plies, it is necessary to obtain its closed-form equation. In this way, stresses expressions can be obtained from the strains expressions. Strains expressions can be obtained by substituting the radial displacement, w0(r), given in equation (4.19), into equations (4.6), which yield: εθ(r) = −R rhp 4γ0 yθ + (1 + hp 5)ε0 θ+hp 6ε0 y+hp 7+r Rhp 8γ0 θy+ ξp 1rκp−1 ¯ Cp 13 +κp¯ Cp 33 +ξp 2r−κp−1 ¯ Cp 13 −κp¯ Cp 33 ,(4.24a) εr(r) = hp 5ε0 θ+hp 6ε0 y+hp 7+ 2 r Rhp 8γ0 θy +κpξp 1rκp−1 ¯ Cp 13 +κp¯ Cp 33 −κpξp 2r−κp−1 ¯ Cp 13 −κp¯ Cp 33 .(4.24b) Substituting the strains given in equations (4.24), (4.4) and (4.5) into the constitutive law (4.10), the stresses can be written as: σp θ(r) = Hp 11ε0 θ+Hp 12ε0 y+R rHp 13γ0 yθ+Hp 14γ0 θy r R+Hp 15ξp 1rκp−1+Hp 16ξp 2r−κp−1+Hp 17, (4.25a) 4.1. Regularized 3D model for singly-curved beams under bending moment and residual strains 148 σp y(r) = Hp 21ε0 θ+Hp 22ε0 y+R rHp 23γ0 yθ+Hp 24γ0 θy r R+Hp 25ξp 1rκp−1+Hp 26ξp 2r−κp−1+Hp 27, (4.25b) σp r(r) = Hp 31ε0 θ+Hp 32ε0 y+R rHp 33γ0 yθ+Hp 34γ0 θy r R+Hp 35ξp 1rκp−1+Hp 36ξp 2r−κp−1+Hp 37, (4.25c) τp θy(r) = Hp 61ε0 θ+Hp 62ε0 y+R rHp 63γ0 yθ+Hp 64γ0 θy r R+Hp 65ξp 1rκp−1+Hp 66ξp 2r−κp−1+Hp 67, (4.25d) where the auxiliary ply constants Hp ij, with i= 1,2,3,6 and j= 1,2..., 7, are defined as follows: Hp i1=¯ Cp i1+ ( ¯ Cp i1+¯ Cp i3)hp 5, i = 1,2,3,6,(4.26a) Hp i2=¯ Cp i2+ ( ¯ Cp i1+¯ Cp i3)hp 6, i = 1,2,3,6,(4.26b) Hp i3=¯ Cp i6−¯ Cp i1hp 4, i = 1,2,3,6,(4.26c) Hp i4=¯ Cp i6+ ( ¯ Cp i1+ 2 ¯ Cp i3)hp 8, i = 1,2,3,6,(4.26d) Hp i5=¯ Cp i1+κp¯ Cp i3 ¯ Cp 13 +κp¯ Cp 33 , Hp i6=¯ Cp i1−κp¯ Cp i3 ¯ Cp 13 −κp¯ Cp 33 , i = 1,2,3,6,(4.26e) Hp i7= ( ¯ Cp i1+¯ Cp i3)hp 7−δp i, i = 1,2,3,6,(4.26f) Therefore, equation (4.25c) is used to apply the continuity of the transversal stresses between plies, in the same way than equation (4.19) is used to apply the continuity of the radial displacements. Displacements calculation The displacement field is defined by its components: the circumferential displacement u(r, θ, y), the axial displacement v(r, θ, y) and the radial displacement w(r, θ, y). Displacements are obtained by integrating the displacement/strain relations, which can be expressed in cylindrical coordinates as: εθ(r) = 1 r ∂u(r, θ, y) ∂θ +w(r, θ, y) r, εy(r) = ∂v(r, θ, y) ∂y , εr(r) = ∂w(r, θ, y) ∂r , (4.27a) γθy(r) = ∂u(r, θ, y) ∂y +1 r ∂v(r, θ, y) ∂θ , γry(r) = ∂w(r, θ, y) ∂y +∂v(r, θ, y) ∂r , (4.27b) γrθ(r) = ∂u(r, θ, y) ∂r −u(r, θ, y) r+1 r ∂w(r, θ, y) ∂θ .(4.27c) 149 Three-dimensional models for evaluating interlaminar stresses Therefore, by introducing (4.4), (4.5), (4.6) and the null values of the interlaminar shear strains (γrθ =γry = 0) into equations (4.27), and integrating the equations, displacements (cancelling the rigid solid motion that appears) are given by: u(r, θ, y) = yr Rγ0 θy +θrε0 θ,(4.28a) v(r, θ, y) = yε0 y+Rθγ0 yθ,(4.28b) w(r, θ, y) = w0(r).(4.28c) Hence, ε0 θis associated to a normal circumferential deformation, ε0 yis associated to a normal axial deformation, γ0 θy is associated to a torsional deformation and γ0 yθ is associated to an in-plane shear deformation in the width direction, as depicted in Figure 4.2. Figure 4.2: Deformations associated to the strain constant parameters. 4.1.2 Loads application and boundary conditions The forces and moments in the laminate are firstly defined for a general loading state depending on the two in-plane coordinates, θand y. After that, the non-dependence of these forces and moments with both in-plane coordinates is considered, in order to obtain their values from the applied loads. Once the values of the forces and moments have been obtained, a system of linear equations is developed in order to determine the strain constant parameters, ε0 θ,ε0 y,γ0 θy and γ0 yθ, and the integration constants of the radial displacement, ξp 1and ξp 2, depending on these loads. 4.1. Regularized 3D model for singly-curved beams under bending moment and residual strains 150 Equilibrium of forces and moments The forces and moments in a cylindrical shell for a general loading state are defined by the following expressions [116, chap. 1]: Nθ(θ, y) = ZR+ R− σθ(r, θ, y)dr, Mθ(θ, y) = ZR+ R− (r−R)σθ(r, θ, y)dr, (4.29a) Nθy(θ, y) = ZR+ R− τθy(r, θ, y)dr, Mθy(θ, y) = ZR+ R− (r−R)τθy(r, θ, y)dr, (4.29b) Ny(θ, y) = ZR+ R− r Rσy(r, θ, y)dr, My(θ, y) = ZR+ R− r R(r−R)σy(r, θ, y)dr, (4.29c) Nyθ(θ, y) = ZR+ R− r Rτθy(r, θ, y)dr, Myθ(θ, y) = ZR+ R− r R(r−R)τθy(r, θ, y)dr, (4.29d) Qθ(θ, y) = ZR+ R− τrθ(r, θ, y)dr, Qy(θ, y) = ZR+ R− r Rτry(r, θ, y)dr, (4.29e) where Nθ, as defined previously, is the membrane force or the circumferential force, Mθis the circumferential moment, Nyis the axial force, Myis the axial moment, Qθis the circumferential shear force, Qyis the axial shear force, Nθy is the in-plane shear force associated to the circumferential direction, Nyθ is the in-plane shear force associated to the axial direction, Mθy is the torsional moment associated to the circumferential direction and Myθ is the torsional moment associated to the axial direction. Notice that forces and moments are defined per unit width due to the shell assumption. Notice that torsional moment associated to the circumferential direction and the in-plane shear forces are related by the following equation: Mθy(θ, y) = R(Nyθ(θ, y)−Nθy(θ, y)).(4.30) The equilibrium equations of forces can be obtained by integrating equations (4.7) in the thickness and considering the value of normal interlaminar stresses and the null value of shear interlaminar stresses in the radial boundaries, yielding: ∂Nθ(θ, y) ∂θ +R∂Nyθ(θ, y) ∂y +Qθ(θ, y) = 0,(4.31a) ∂Nθy(θ, y) ∂θ +R∂Ny(θ, y) ∂y = 0,(4.31b) 151 Three-dimensional models for evaluating interlaminar stresses ∂Qθ(θ, y) ∂θ +R∂Qy(θ, y) ∂y +R+p+−R−p−=Nθ(θ, y).(4.31c) The equilibrium equations of moments can be obtained by integrating equations (4.7a) and (4.7b) in the thickness multiplied by (r−R) and considering a null value of shear interlaminar stresses in the radial boundaries. These equilibrium equations yield: ∂Mθ(θ, y) ∂θ +R∂Myθ(θ, y) ∂y =RQθ(θ, y),(4.32a) ∂Mθy(θ, y) ∂θ +R∂My(θ, y) ∂y =RQy(θ, y).(4.32b) Suppressing Mθy(θ, y) by substituting (4.30) in (4.32b): ∂Nyθ(θ, y) ∂θ −∂Nθy(θ, y) ∂θ +∂My(θ, y) ∂y =Qy(θ, y).(4.33) Considering that the zone under study is far enough from the free edges and the changes of curvature and under a loading state for which stresses do not depend on θand y, then the forces and moments neither depend on θnor on y. Substituting these conditions in equations (4.31c) and (4.32) the circumferential and axial shear forces, Qθand Qy, are null and the circumferential force, Nθ, is given by the membrane forces (4.23). The absence of the shear forces was previously demonstrated with the absence of the interlaminar shear stresses, τrθ and τry. However, these conditions are not sufficient to solve the problem. In that way, to get the conditions to solve the problem it is necessary to consider the non-dependence with each coordinate independently of the other. Hence, considering only the non-dependence of the forces and moments with the y-direction, equations (4.31) and (4.32a) yield: dNθ(θ) dθ +Qθ(θ) = 0,dQθ(θ) dθ +R+p+−R−p−=Nθ(θ),(4.34a) dNθy(θ) dθ = 0,(4.34b) dMθ(θ) dθ =RQθ(θ).(4.34c) Notice that Qθ(θ) is null and Nθ(θ) is given by the membrane forces in the zone far enough from the free edges and the changes of curvature. However, equation (4.34a) implies a sinusoidal variation of these forces plus the membrane force, so taking into account that far enough from the free 4.1. Regularized 3D model for singly-curved beams under bending moment and residual strains 152 edges and the changes of curvatures in the θ-direction they have to be equal to the membrane value, implies that Qθ(θ) = 0 and Nθ(θ) is given by the membrane forces (4.23). Substituting the null value of Qθ(θ) in equation (4.34c) it implies that Mθ(θ) is constant, so that Mθ(θ) = M0with M0being the applied bending moment at the θboundary. Equation (4.34b) implies also a constant value of Nθy(θ) = P0, with P0being the applied in-plane shear force in the circumferential direction at the θboundary. Notice that M0and P0are not respectively the total applied bending moment and in-plane shear force divided between the width of the specimen in the axial direction. This is due to the consideration of being far enough from the free edges in the y-direction. Near to the free edges the distribution may be different so Mθ(θ) and Nθy(θ) are not constant in the whole width. However, when the width is much higher than the thickness Mθ(θ) and Nθy(θ) can be approximated by the total applied bending moment and inplane shear force divided by the width. If a closed cylinder is calculated (a pipe), the values of M0and P0may not be easy to obtain. However, in this case the displacements have to take the same value at θ= 0 and at θ= 2π. Introducing this condition into (4.28) the two equations of M0and P0can be substituted by null values of ε0 θand γ0 yθ. In the same way, considering only the non-dependence of the forces and moments with the θ-direction, equations (4.31a), (4.31b) and (4.32a) yield: dNy(y) dy = 0,(4.35a) RdNyθ(y) dy +Qθ(y) = 0,dMyθ(y) dy =Qθ(y).(4.35b) Equation (4.35a) implies a constant value of the axial force, Ny(y) = N0, where N0is the applied axial load. On the other side, equations (4.35b) imply, suppressing Qθ(y) by substituting one equation in the other, a constant value of Myθ(y) + RNyθ(y) = T0, where T0is the applied torsional moment. Therefore, in the studied zone where stresses and strains depend only on the radial coordinate, the four remaining equations to close the problem and determining the strain constant parameters, ε0 θ,ε0 y,γ0 θy and γ0 yθ, are given by the applied loads Mθ=M0,Ny=N0,Nθy =P0and Myθ +RNyθ =T0. The forces and moments not included in the previous conditions, such as the membrane force (Nθ) or the axial moment (My), are not inputs of the problems as they are results of other loads or effects, such as the pressure or the anticlastic effect combined with the curvature, and they are calculated from the input loads. 153 Three-dimensional models for evaluating interlaminar stresses Application of the external loads To apply the external loads equation it is necessary to integrate equations (4.25) accordingly to the forces and moments definitions. Therefore, substituting the circumferential stress (4.25a) into the definition of the circumferential bending moment (4.29a) and considering that Mθ=M0, the applied bending moment per unit of width can be written as: M0= Np X p=1 ZR+ p R− p (r−R)σp θ(r)dr = Np X p=1 KM 1pε0 θ+KM 2pε0 y+KM 3pγ0 yθ +KM 4pγ0 θy +KM 5pξp 1+KM 6pξp 2+KM 7p,(4.36) where the auxiliary ply constants KM jp (with j= 1,2, ..., 7 and p= 1,2, ..., Np) are defined as follows: KM jp =Hp 1j R+ p−R− p 2R+ p+R− p−2R, j = 1,2,7,(4.37a) KM 3p=Hp 13RR+ p−R− p−Rlog R+ p R− p,(4.37b) KM 4p=Hp 14 (R+ p)3−(R− p)3 3R−(R+ p)2−(R− p)2 2,(4.37c) KM 5p=Hp 15 (R+ p)κp+1 −(R− p)κp+1 κp+ 1 −R(R+ p)κp−(R− p)κp κp,(4.37d) KM 6p=Hp 16 (R+ p)1−κp−(R− p)1−κp 1−κp+R(R+ p)−κp−(R− p)−κp κp.(4.37e) Substituting the axial stress (4.25b) into the definition of the axial force (4.29c) and considering that Ny=N0, the applied axial force per unit of length can be written as: N0= Np X p=1 ZR+ p R− p r Rσp y(r)dr = Np X p=1 KN 1pε0 θ+KN 2pε0 y+KN 3pγ0 yθ +KN 4pγ0 θy +KN 5pξp 1+KN 6pξp 2+KN 7p,(4.38) B.4. Auxiliary matrices for the curved laminate 256 Bε 16 = Np X p=1 ¯ Cp 16 Zz+ p z− p z+R Rzdz!, Bγ 66 = Np X p=1 ¯ Cp 66 Zz+ p z− p zdz!,(B.19c) kε 11 = Np X p=1 ¯ Cp 11 Zz+ p z− p fs(z)dz!,(B.20a) kγ 16 = Np X p=1 ¯ Cp 16 Zz+ p z− p R z+Rfs(z)dz!,(B.20b) kε 12 = Np X p=1 ¯ Cp 12 Zz+ p z− p z+R Rfy(z)dz!,(B.20c) kγ 26 = Np X p=1 ¯ Cp 26 Zz+ p z− p fy(z)dz!,(B.20d) kε 16 = Np X p=1 ¯ Cp 16 Zz+ p z− p z+R Rfs(z)dz!,(B.20e) kγ 66 = Np X p=1 ¯ Cp 66 Zz+ p z− p fs(z)dz!,(B.20f) kε 13 = Np X p=1 ¯ Cp 13 Zz+ p z− pz+R R2ˆ ˆ fs(z)dz!,(B.20g) kγ 36 = Np X p=1 ¯ Cp 36 Zz+ p z− p z+R R ˆ ˆ fs(z)dz!.(B.20h) B.4 Auxiliary matrices for the curved laminate The auxiliary matrices employed in Chapter 4 in the development of the non-regularized model for the curved laminates are defined as follows: HQ 00 =ˇ A44 −(ˇ kC 44)THγ−1hQ−hQTHQ 10,(B.21a) HQ 01 = (ˇ kC 44)THγ−1−hQTHQ 11,(B.21b) HQ 02 =ˇ A45 −(ˇ kC 45)THγ−1hQ−hQTHQ 12,(B.21c) HQ 03 = (ˇ kC 45)THγ−1−hQTHQ 13,(B.21d) 257 Matrices expressions of the 3D non-regularized models HQ 10 = (HγT)−1ˇ kC 44 −ˇ Kτ 44Hγ−1hQ,(B.21e) HQ 11 = (HγT)−1ˇ Kτ 44Hγ−1,(B.21f) HQ 12 = (HγT)−1ˇ kC 45 −ˇ Kτ 45Hγ−1hQ,(B.21g) HQ 13 = (HγT)−1ˇ Kτ 45Hγ−1,(B.21h) HQ 20 =ˇ A45 −(ˇ kC 45)THγ−1hQ−hQTHQ 30,(B.21i) HQ 21 = (ˇ kC 45)THγ−1−hQTHQ 31,(B.21j) HQ 22 =ˇ A55 −(ˇ kC 55)THγ−1hQ−hQTHQ 32,(B.21k) HQ 23 = (ˇ kC 55)THγ−1−hQTHQ 33,(B.21l) HQ 30 = (HγT)−1ˇ kC 45 −ˇ Kτ 45Hγ−1hQ,(B.21m) HQ 31 = (HγT)−1ˇ Kτ 45Hγ−1,(B.21n) HQ 32 = (HγT)−1ˇ kC 55 −ˇ Kτ 55Hγ−1hQ,(B.21o) HQ 33 = (HγT)−1ˇ Kτ 55Hγ−1.(B.21p) The matrix ˇ KAyield: ˇ KA=                 ˇ KA 11 ˇ KA 12 ˇ KA 13 ˇ KA 14 ˇ KA 15 ˇ KA 16 ˇ KA 17 ˇ KA 18 ˇ KA 21 ˇ KA 22 ˇ KA 23 ˇ KA 24 ˇ KA 25 ˇ KA 26 ˇ KA 27 ˇ KA 28 ˇ KA 31 ˇ KA 32 ˇ KA 33 ˇ KA 34 ˇ KA 35 ˇ KA 36 ˇ KA 37 ˇ KA 38 ˇ KA 41 ˇ KA 42 ˇ KA 43 ˇ KA 44 ˇ KA 45 ˇ KA 46 ˇ KA 47 ˇ KA 48 ˇ KA 51 ˇ KA 52 ˇ KA 53 ˇ KA 54 ˇ KA 55 ˇ KA 56 ˇ KA 57 ˇ KA 58 ˇ KA 61 ˇ KA 62 ˇ KA 63 ˇ KA 64 ˇ KA 65 ˇ KA 66 ˇ KA 67 ˇ KA 68 ˇ KA 71 ˇ KA 72 ˇ KA 73 ˇ KA 74 ˇ KA 75 ˇ KA 76 ˇ KA 77 ˇ KA 78 ˇ KA 81 ˇ KA 82 ˇ KA 83 ˇ KA 84 ˇ KA 85 ˇ KA 86 ˇ KA 87 ˇ KA 88 ˇ KA 91 ˇ KA 92 ˇ KA 93 ˇ KA 94 ˇ KA 95 ˇ KA 96 ˇ KA 97 ˇ KA 98 ˇ KA z1ˇ KA z2ˇ KA z3ˇ KA z4ˇ KA z5ˇ KA z6ˇ KA z7ˇ KA z8                 ,(B.22) where the columns of ˇ KAare defined from the columns of Kσas follows:      ˇ KA 11 ˇ KA 21 . . . ˇ KA z1      =     ˇ B11 ˇ D11 . . . ˇ kB 13     −     (ˇ kA 13)T (ˇ kB 13)T . . . ˇ Kσ 33      Hz−1hM,(B.23a) B.4. Auxiliary matrices for the curved laminate 258      ˇ KA 12 ˇ KA 22 . . . ˇ KA z2      =     (ˇ kA 11)T (ˇ kB 11)T . . . (ˇ Kσ 13)T     −     (ˇ kA 13)T (ˇ kB 13)T . . . ˇ Kσ 33      Hz−1Hs,(B.23b)      ˇ KA 13 ˇ KA 23 . . . ˇ KA z3      =R     ˇ A11 ˇ B11 . . . ˇ kA 13     −     (ˇ kA 13)T (ˇ kB 13)T . . . ˇ Kσ 33      Hz−1(Hγ−1hQ+RhN),(B.23c)      ˇ KA 14 ˇ KA 24 . . . ˇ KA z4      =     ˇ B12 ˇ D12 . . . ˇ kB 23      ,     ˇ KA 15 ˇ KA 25 . . . ˇ KA z5      =     (ˇ kA 13)T (ˇ kB 13)T . . . ˇ Kσ 33      Hz−1Hγ−1,(B.23d)      ˇ KA 16 ˇ KA 26 . . . ˇ KA z6      =     (ˇ kA 12)T (ˇ kB 12)T . . . (ˇ Kσ 23)T      ,     ˇ KA 17 ˇ KA 27 . . . ˇ KA z7      =     ˇ B16 ˇ D16 . . . ˇ kB 36     −1 R     ˇ A16 ˇ B16 . . . ˇ kA 36      ,(B.23e)      ˇ KA 18 ˇ KA 28 . . . ˇ KA z8      =     (ˇ kA 16)T (ˇ kB 16)T . . . (ˇ Kσ 36)T      ,(B.23f) The matrices ˇ KBand ˇ KCand the vector ˇ KDare obtained by the following equations: ˇ KB=ˇ KB 11 ˇ KB 12 ˇ KB 13 ˇ KB 14 ˇ KB 15 ˇ KB 16 ˇ KB 17 ˇ KB 21 ˇ KB 22 ˇ KB 23 ˇ KB 24 ˇ KB 25 ˇ KB 26 ˇ KB 27= −ˇ KA 21 ˇ KA 22 ˇ KA 31 ˇ KA 32−1ˇ KA 23 ˇ KA 24 ˇ KA 25 ˇ KA 26 ˇ KA 27 ˇ KA 28 −hMT ˇ KA 33 ˇ KA 34 ˇ KA 35 ˇ KA 36 ˇ KA 37 ˇ KA 38 −HsT, (B.24a) ˇ KC=KC 11 KC 12 KC 21 KC 22=−ˇ KA 21 ˇ KA 22 ˇ KA 31 ˇ KA 32−1Bε 11 Bγ 16 kε 11 kγ 16,(B.24b) ˇ KD=ˇ KD 1 ˇ KD 2=ˇ KA 21 ˇ KA 22 ˇ KA 31 ˇ KA 32−1ˇ kε 2 ˇ kε 3,(B.24c) where the components ˇ kε i,i= 1,2, ..., 9, z, are defined in equation (4.149). 259 Matrices expressions of the 3D non-regularized models Finally, matrices Gand Gcand vector gεare defined by the following expressions: G=GA−1GB,Gc=GA−1GC,gε=GA−1Hε,(B.25) where the auxiliary matrices GA,GBand GCand the auxiliary vector Hε are defined by the following expressions: GA=         GA 11 GA 12 GA 13 GA 14 GA 15 GA 16 GA 21 GA 22 GA 23 GA 24 GA 25 GA 26 GA 31 GA 32 GA 33 GA 34 GA 35 GA 36 GA 41 GA 42 GA 43 GA 44 GA 45 GA 46 GA 51 GA 52 GA 53 GA 54 GA 55 GA 56 GA 61 GA 62 GA 63 GA 64 GA 65 GA 66         ,(B.26a) GB=         GB 11 GB 12 GB 13 GB 14 GB 15 GB 16 GB 21 GB 22 GB 23 GB 24 GB 25 GB 26 GB 31 GB 32 GB 33 GB 34 GB 35 GB 36 GB 41 GB 42 GB 43 GB 44 GB 45 GB 46 GB 51 GB 52 GB 53 GB 54 GB 55 GB 56 GB 61 GB 62 GB 63 GB 64 GB 65 GB 66         ,(B.26b) GC=         GC 11 GC 12 GC 21 GC 22 GC 31 GC 32 GC 41 GC 42 GC 51 GC 52 GC 61 GC 62         ,Hε=         Hε 1 Hε 2 Hε 3 Hε 4 Hε 5 Hε 6         ,(B.26c) the components of matrix GAyielding: GA 11 =HQ 00 −ˇ KA 53 −ˇ KA 51 ˇ KB 11 −ˇ KA 52 ˇ KB 21,(B.27a) GA 12 =HQ 01 −ˇ KA 55 −ˇ KA 51 ˇ KB 13 −ˇ KA 52 ˇ KB 23,(B.27b) GA 13 =HQ 02,GA 14 =HQ 03, GA 15 =−1,GA 16 =hQT(HγT)−1,(B.27c) GA 21 =HQ 10 −ˇ KA 63 −ˇ KA 61 ˇ KB 11 −ˇ KA 62 ˇ KB 21,(B.27d) GA 22 =HQ 11 −ˇ KA 65 −ˇ KA 61 ˇ KB 13 −ˇ KA 62 ˇ KB 23,(B.27e) GA 23 =HQ 12,GA 24 =HQ 13,GA 25 =0,GA 26 =−(HγT)−1,(B.27f) GA 31 =HQ 20−ˇ KA 83+ˇ KA 73 R−ˇ KA 81 −ˇ KA 71 Rˇ KB 11−ˇ KA 82 −ˇ KA 72 Rˇ KB 21,(B.27g) B.4. Auxiliary matrices for the curved laminate 260 GA 32 =HQ 21 −ˇ KA 85 +ˇ KA 75 R−ˇ KA 81 −ˇ KA 71 Rˇ KB 13 −ˇ KA 82 −ˇ KA 72 Rˇ KB 23, (B.27h) GA 33 =HQ 22,GA 34 =HQ 23, GA 35 = 0,GA 36 =0,(B.27i) GA 41 =HQ 30 −ˇ KA 93 −ˇ KA 91 ˇ KB 11 −ˇ KA 92 ˇ KB 21,(B.27j) GA 42 =HQ 31 −ˇ KA 95 −ˇ KA 91 ˇ KB 13 −ˇ KA 92 ˇ KB 23,(B.27k) GA 43 =HQ 32,GA 44 =HQ 33,GA 45 =0,GA 46 =0,(B.27l) GA 51 =−ˇ KA 13 −ˇ KA 11 ˇ KB 11 −ˇ KA 12 ˇ KB 21,(B.27m) GA 52 =−ˇ KA 15 −ˇ KA 11 ˇ KB 13 −ˇ KA 12 ˇ KB 23,(B.27n) GA 53 = 0,GA 54 =0, GA 55 = 0,GA 56 =0,(B.27o) GA 61 =−ˇ KA z3−ˇ KA z1ˇ KB 11 −ˇ KA z2ˇ KB 21,(B.27p) GA 62 =−ˇ KA z5−ˇ KA z1ˇ KB 13 −ˇ KA z2ˇ KB 23,(B.27q) GA 63 =0,GA 64 =0,GA 65 =0,GA 66 =0,(B.27r) the components of matrix GByielding: GB 11 =ˇ KA 54 +ˇ KA 51 ˇ KB 12 +ˇ KA 52 ˇ KB 22,GB 12 =ˇ KA 56 +ˇ KA 51 ˇ KB 14 +ˇ KA 52 ˇ KB 24, (B.28a) GB 13 =ˇ KA 57 +ˇ KA 51 ˇ KB 15 +ˇ KA 52 ˇ KB 25,GB 14 =ˇ KA 58 +ˇ KA 51 ˇ KB 16 +ˇ KA 52 ˇ KB 26, (B.28b) GB 15 = 0,GB 16 =ˇ KA 51 ˇ KB 17 +ˇ KA 52 ˇ KB 27,(B.28c) GB 21 =ˇ KA 64 +ˇ KA 61 ˇ KB 12 +ˇ KA 62 ˇ KB 22,GB 22 =ˇ KA 66 +ˇ KA 61 ˇ KB 14 +ˇ KA 62 ˇ KB 24, (B.28d) GB 23 =ˇ KA 67 +ˇ KA 61 ˇ KB 15 +ˇ KA 62 ˇ KB 25,GB 24 =ˇ KA 68 +ˇ KA 61 ˇ KB 16 +ˇ KA 62 ˇ KB 26, (B.28e) GB 25 =0,GB 26 =ˇ KA 61 ˇ KB 17 +ˇ KA 62 ˇ KB 27,(B.28f) GB 31 =ˇ KA 84 −ˇ KA 74 R+ˇ KA 81 −ˇ KA 71 Rˇ KB 12 +ˇ KA 82 −ˇ KA 72 Rˇ KB 22,(B.28g) 261 Matrices expressions of the 3D non-regularized models GB 32 =ˇ KA 86 −ˇ KA 76 R+ˇ KA 81 −ˇ KA 71 Rˇ KB 14 +ˇ KA 82 −ˇ KA 72 Rˇ KB 24,(B.28h) GB 33 =ˇ KA 87 −ˇ KA 77 R+ˇ KA 81 −ˇ KA 71 Rˇ KB 15 +ˇ KA 82 −ˇ KA 72 Rˇ KB 25,(B.28i) GB 34 =ˇ KA 88 −ˇ KA 78 R+ˇ KA 81 −ˇ KA 71 Rˇ KB 16 +ˇ KA 82 −ˇ KA 72 Rˇ KB 26,(B.28j) GB 35 = 0,GB 36 =ˇ KA 81 −ˇ KA 71 Rˇ KB 17 +ˇ KA 82 −ˇ KA 72 Rˇ KB 27,(B.28k) GB 41 =ˇ KA 94 +ˇ KA 91 ˇ KB 12 +ˇ KA 92 ˇ KB 22,GB 42 =ˇ KA 96 +ˇ KA 91 ˇ KB 14 +ˇ KA 92 ˇ KB 24, (B.28l) GB 43 =ˇ KA 97 +ˇ KA 91 ˇ KB 15 +ˇ KA 92 ˇ KB 25,GB 44 =ˇ KA 98 +ˇ KA 91 ˇ KB 16 +ˇ KA 92 ˇ KB 26, (B.28m) GB 45 =0,GB 46 =ˇ KA 91 ˇ KB 17 +ˇ KA 92 ˇ KB 27,(B.28n) GB 51 =ˇ KA 14 +ˇ KA 11 ˇ KB 12 +ˇ KA 12 ˇ KB 22,GB 52 =ˇ KA 16 +ˇ KA 11 ˇ KB 14 +ˇ KA 12 ˇ KB 24, (B.28o) GB 53 =ˇ KA 17 +ˇ KA 11 ˇ KB 15 +ˇ KA 12 ˇ KB 25,GB 54 =ˇ KA 18 +ˇ KA 11 ˇ KB 16 +ˇ KA 12 ˇ KB 26, (B.28p) GB 55 =−1 R,GB 56 =ˇ KA 11 ˇ KB 17 +ˇ KA 12 ˇ KB 27 −hNT,(B.28q) GB 61 =ˇ KA z4+ˇ KA z1ˇ KB 12 +ˇ KA z2ˇ KB 22,GB 62 =ˇ KA z6+ˇ KA z1ˇ KB 14 +ˇ KA z2ˇ KB 24, (B.28r) GB 63 =ˇ KA z7+ˇ KA z1ˇ KB 15 +ˇ KA z2ˇ KB 25,GB 64 =ˇ KA z8+ˇ KA z1ˇ KB 16 +ˇ KA z2ˇ KB 26, (B.28s) GB 65 =0,GB 66 =ˇ KA 10,1ˇ KB 17 +ˇ KA 10,2ˇ KB 27 −HzT,(B.28t) the components of matrix GCyielding: GC 11 =Bε 12 +ˇ KA 51KC 11 +ˇ KA 52KC 21, GC 12 =Bγ 26 +ˇ KA 51KC 12 +ˇ KA 52KC 22,(B.29a) GC 21 =kε 12 +ˇ KA 61KC 11 +ˇ KA 62KC 21,GC 22 =kγ 26 +ˇ KA 61KC 12 +ˇ KA 62KC 22, (B.29b) GC 31 =Bε 16 −Aε 16 R+ˇ KA 81 −ˇ KA 71 RKC 11 +ˇ KA 82 −ˇ KA 72 RKC 21,(B.29c) B.4. Auxiliary matrices for the curved laminate 262 GC 32 =Bγ 66 −Aγ 66 R+ˇ KA 81 −ˇ KA 71 RKC 12 +ˇ KA 82 −ˇ KA 72 RKC 22,(B.29d) GC 41 =kε 16+ˇ KA 91KC 11+ˇ KA 92KC 21,GC 42 =kγ 66+ˇ KA 91KC 12+ˇ KA 92KC 22,(B.29e) GC 51 =Aε 11 +ˇ KA 11KC 11 +ˇ KA 12KC 21, GC 52 =Aγ 16 +ˇ KA 11KC 12 +ˇ KA 12KC 22,(B.29f) GC 61 =kε 13 +ˇ KA 10,1KC 11 +ˇ KA 10,2KC 21,(B.29g) GC 62 =kγ 36 +ˇ KA 10,1KC 12 +ˇ KA 10,2KC 22,(B.29h) and the components of matrix Hεyielding: Hε 1=ˇ KA 51 ˇ KD 1+ˇ KA 52 ˇ KD 2−ˇ kε 5,Hε 2=ˇ KA 61 ˇ KD 1+ˇ KA 62 ˇ KD 2−ˇ kε 6,(B.30a) Hε 3=ˇ KA 81 −ˇ KA 71 Rˇ KD 1+ˇ KA 82 −ˇ KA 72 Rˇ KD 2−ˇ kε 8+ˇ kε 7 R,(B.30b) Hε 4=ˇ KA 91 ˇ KD 1+ˇ KA 92 ˇ KD 2−ˇ kε 9, Hε 5=ˇ KA 11 ˇ KD 1+ˇ KA 12 ˇ KD 2−ˇ kε 1,(B.30c) Hε 6=ˇ KA 10,1ˇ KD 1+ˇ KA 10,2ˇ KD 2−ˇ kε z.(B.30d) Bibliography [1] K.T. Kedward, R.S. Wilson, and S. K. Mclean. Flexure of simply curved composite shapes. Composites, 20:527–536, 1989. [2] P.C. Paul, C.R. Saff, K.B. S¨anger, M.A. Mahler, H.P. Kan, and E.F. Kautz. Out of plane analysis for composite structures. In Eighth DoD/NASA/FAA Conference on Fibrous Composites in Structural Design, pages 263–279, 1992. Norfolk, Virginia. [3] ASTM D 6415/D 6415M - 06a. Standard test method for measuring the curved beam strength of a fibre-reinforced polymer-matrix composite. ASTM Standards, 2007. [4] W.L. Ko and R.H. Jackson. Multilayer theory for delamination analysis of a composite curved bar subjected to end forces and end moments. NASA Technical Memorandum 4139, 1989. [5] R.Y. Kim and S.R. Soni. Failure of composite laminates due to combined interlaminar normal and shear stresses. In Composites’86: recent advances in Japan and the United States, Proceedings of Japan-U.S. CCM-III, pages 341–350, 1986. [6] T. Edwards and J. Thompson. Spar corner radius integrity for the A400M wing. Applied Mechanics and Materials, 3-4:197–202, 2005. [7] M. Hoffmann, K. Zimmermann, B. Bautz, and P. Middendorf. Size effect on through-thickness strength properties of 3D loaded composite laminates. In Proceedings of the 17th European Conference on Composite Materials ECCM17, 2016. 26-30 of June, Munich, Germany. [8] S.C. Avalon and S.L. Donaldson. Strength of composite angle brackets with multiple geometries and nanofiber-enhanced resins. Journal of Composite Materials, 45(9):1017–1030, 2010. [9] M.R. Wisnom. Size effects in the testing of fibre-composite materials. Composites Science and Technology, 59:1937–1957, 1999. Bibliography 264 [10] A. Makeev, G. Seon, Y. Nikishkov, and E. Lee. Methods for assessment of interlaminar tensile strength of composite materials. Journal of Composite Materials, 49(7):783–794, 2015. [11] T.A. Fletcher, T. Kim, T.J. Dodwell, R. Butler, R. Scheichl, and R. Newley. The influence of free edges on curved beam strength. In Proceedings of the 21th International Conference on Composite Materials ICCM20, 2015. 19-24 of July, Copenhaguen, Denmark. [12] T.A. Fletcher, A.K. Reinarz, T.J. Dodwell, R. Butler, R. Scheichl, and R Newley. Efficient modelling and accurate certification of curved aerospace laminates. In Proceedings of the 17th European Conference on Composite Materials ECCM17, 2016. 26-30 of June, Munich, Germany. [13] T. Kim, T.A. Fletcher, T.J. Dodwell, R. Butler, R. Scheichl, J. Ankersen, and R. Newley. The effect of free edges on inter-laminar performance of curved laminates. In 56th AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 2015. [14] J. Most, D. Stegmair, and D. Petry. Error estimation between simple, closed-form analytical formulae and full-scale FEM for interlaminar stress prediction in curved laminates. Composite Structures, 131:72– 81, 2015. [15] R.M. Jones. Mechanics of Composite Materials. Mc Graw-Hill, 1975. [16] J.M. Whitney and N.J. Pagano. Shear deformation in heterogeneous anisotropic plates. Journal of Applied Mechanics, 37:1031–1036, 1970. [17] K.C. Lin and C.M. Hsieh. The closed form general solutions of 2-D curved laminated beams of variable curvatures. Composite Structures, 79(4):606–618, 2007. [18] K.C. Lin and C.W. Lin. Finite deformation of 2-D laminated curved beams with variable curvatures. International Journal of Non-Linear Mechanics, 46(10):1293–1304, 2011. [19] S.G. Lekhnitskii, S.W. Tsai, and T. Cheron. Anisotropic Plates. Gordon and Breach Science Publishers, 1968. [20] W.L. Ko. Delamination stresses in semicircular laminated composite curved bars. NASA Technical Memorandum 4026, 1988. 265 Bibliography [21] A. Sharma and C.E. Bakis. Analysis of elastic stresses in thick, polarorthotropic, C-shaped rings. Journal of Composite Materials, 38:1619– 1638, 2004. [22] A. Sharma and C.E. Bakis. C-shape specimen for tensile radial strength of thick, filament-wound rings. Journal of Composite Materials, 40(2):97–116, 2006. [23] A. Schmitz and P. Horst. Bending deformation limits of corrugated unidirectionally reinforced composites. Composite Structures, 107:103–111, 2014. [24] W. Cui, T. Liu, J. Len, and R. Ruo. Interlaminar tensile strength (ILTS) measurement of woven glass/polyester laminates specimen using four-point curved beam specimen. Composites Part A, 27A:1097– 1105, 1996. [25] S. Smidt. Bending of curved sandwich beams. Composite Structures, 33:211–225, 1995. [26] R.A. Shenoi and W. Wang. Through-thickness stresses in curved composite laminates and sandwich beams. Composites Science and Technology, 61(11):1501–1512, 2001. [27] M.S. Qatu. Theories and analyses of thin and moderately thick laminated composite curved beams. International Journal of Solids and Structures, 30(20):2743–2756, 1993. [28] G. Kress, R. Roos, M. Barbezat, C. Dransfeld, and P. Ermanni. Model for interlaminar normal stress in singly curved laminates. Composite Structures, 69(4):458–469, 2005. [29] R.M. Guedes and A. S´a. Numerical analysis of singly curved shallow composite panels under three-point bend load. Composite Structures, 83(2):212–220, 2008. [30] G. Kress. Minimized computational effort for the thick-walled composite tube problem. Computers and Structures, 54(4):633–639, 1995. [31] V. Rodr´ıguez de la Cruz, D. Fern´andez Caballero, F. Mujika, and J.M. Mu˜noz Guijosa. Analysis of out-of-plane stresses in sandwich beams subjected to pure bending with large deflections. Journal of Composite Materials, 47(22):2809–2822, 2012. Bibliography 272 [99] J.M. Gonz´alez-Cantero, E. Graciani, F. Par´ıs, and B. L´opez-Romano. Stacking sequence optimization of curved UD-CFRP laminates for improving unfolding strength considering thermal residual stresses. In Proceedings of the 17th European Conference on Composite Materials ECCM17, 2016. 26-30 of June, Munich, Germany. [100] J.M. Gonz´alez-Cantero, E. Graciani, B. L´opez-Romano, and F. Par´ıs. Induced unfolding: a new mechanism of unfolding failure in composite laminates. Submitted to Composites Science and Technology, 2017. [101] W.C. Jackson and R.H. Martin. An interlaminar tensile strength specimen. Composite Materials: Testing and Design (11th Volume), ASTM STP 1206, pages 333–354, 1993. [102] L. Ram´ırez L´opez. Fabricaci´on y ensayo de probetas ILTS. Degree Final Project, Superior Engineering School, University of Seville, 2015. [103] J.M. Gonz´alez-Cantero, E. Graciani, A. Bl´azquez, and F. Par´ıs. Analytic evaluation of radial stresses in unfolding failure of composite materials. Comparison with numerical solutions. In Proceedings of the 16th European Conference on Composite Materials ECCM16, 2014. 22-26 of June, Seville, Spain. [104] G. Arfken. Orthogonal Polynomials. Mathematical Methods for Physicists, 3rd edition, 1985. [105] M. Abramowitz and I.A. Stegun. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, chapter 8. Dover Books on Mathematics, 1965. [106] D. Zwillinger. Handbook of Differential Equations, chapter 22. Boston, MA: Academic Press, 3rd edition, 1997. [107] M.E.A. El-Mikkawy and G.-S. Cheon. Combinatorial and hypergeometric identities via the Legendre polynomials - A computational approach. Applied Mathematics and Computation, 166(1):181–195, 2005. [108] J.M. Gonz´alez-Cantero, E. Graciani, F. Par´ıs, B. L´opez-Romano, and D. Meizoso-Latova. Semi-analytic solution of non-regularized unfolding stresses in composite beams employing a series approximation based on Legendre polynomials. In Proceedings of the 21th International Conference on Composite Materials ICCM20, 2015. 19-24 of July, Copenhaguen, Denmark. 273 Bibliography [109] J.M. Gonz´alez-Cantero, E. Graciani, F. Par´ıs, D. Meizoso-Latova, and J.A. Rodr´ıguez-S´anchez. Resultados semi-anal´ıticos de tensiones no regularizadas en barras curvas de material compuesto usando un desarrollo en serie de funciones monomiales. In Proceedings of the XI Congreso Nacional de Materiales Compuestos, Matcomp’15, 2015. 6- 8 of July, M´ostoles, Spain. [110] J.M. Gonz´alez-Cantero, E. Graciani, F. Par´ıs, and B. L´opez-Romano. Semi-analytical stress calculation in composite materials using a 2D series expansion of the displacements. In Proceedings of the 19th International Conference on Composite Structures ICCS19, 2016. 5-9 of September, Porto, Portugal. [111] E.E. Theotokoglou and E. Sideridis. Study of composite beams in asymmetric four-point bending. Journal of Reinforced Plastics and Composites, 30(13):1125–1137, 2011. [112] Y.C. Chiang. Curved laminate analysis. Structural Engineering and Mechanics, 39(2):169–186, 2011. [113] C.P. Wu and Y.W. Chi. Three-dimensional nonlinear analysis of laminated cylindrical shells under cylindrical bending. European Journal of Mechanics A/Solids, 24:837–856, 2005. [114] O.G. Kravchenko, C. Li, A. Strachanand, S.G. Kravchenko, and R.B. Pipes. Prediction of the chemical and thermal shrinkage in a thermoset polymer. Composites: Part A, 66:35–43, 2014. [115] O.G. Kravchenko, S.G. Kravchenko, and R.B. Pipes. Chemical and thermal shrinkage in thermosetting prepreg. Composites Part A, 80:72–81, 2016. [116] Wilhelm Fl¨ugge. Stresses in Shells. Springer-Verlag Berlin, second edition, 1973. [117] R. Roos. Model for Interlaminar Normal Stresses in Doubly Curved Laminates. PhD thesis, Swiss Federal Institute of Technology, Zurich, 2008. [118] A. Bl´azquez, V. Mantiˇc, F. Par´ıs, and N.L. McCartney. BEM analysis of damage progress in 0/90 laminates. Engineering Analysis with Boundary Elements, 33:762–769, 2009. Bibliography 274 [119] F. Par´ıs, A. Bl´azquez, N.L. McCartney, and V. Mantiˇc. Characterization and evolution of matrix and interface related damage in [0/90]S laminates under tension. Part I: Numerical predictions. Composites Science and Technology, 70:1168–1175, 2010. [120] F. Par´ıs, A. Bl´azquez, N.L. McCartney, and A. Barroso. Characterization and evolution of matrix and interface related damage in [0/90]S laminates under tension. Part II: Experimental evidence. Composites Science and Technology, 70:1176–1183, 2010. Curriculum Vitae Name: Juan Manuel Gonz´alez Cantero Date of Birth: December 30, 1990 Nationality: Spanish E-mail: Juan.Man[email protected] Education: •PhD, Universidad de Sevilla, Escuela T´ecnica Superior de Ingenieros (Sept 2014 - Jul 2017). •63 credits in Advanced Design in Mechanical Engineering Master at the University of Seville (2014). •Aeronautical Engineer at the Engineering School of Seville University with an average mark of 8.93/10 (Sept 2008 - Oct 2013). Awards: –Award for the Outstanding Graduate of the Year. –Best final year project on an aeronautical topic (C´atedra EADS). •Erasmus FORMAT-EO (FORmation of Multi-disciplinary Approaches to Training in Earth Observation) intensive course in the University of Leicester (July 2013). Working experience: •PhD Student and stress engineer in FIDAMC (Spain) since February 2014, working on national and European R&D projects. •Internship at Airbus Group Innovations, within the department ”Lightweight Design & Optimization” at the Airbus Group site Ottobrunn (Germany) since the 01/04/2016 until 31/07/2016. •Internship at the Armament and Experimentation Logistic Center of the Spanish Air Force at August 2013. Curriculum Vitae 276 Publications: •J.M. Gonz´alez-Cantero, E. Graciani, A. Bl´azquez, and F. Par´ıs. A new analytical model for evaluating interlaminar stresses in the unfolding failure of composite laminates. Composite Structures, 147:260–273, 2016. •J.M. Gonz´alez-Cantero, E. Graciani, F. Par´ıs, and B. L´opez-Romano. Semi-analytic model to evaluate non-regularized stresses causing unfolding failure in composites. Composite Structures, 171:77–91, 2017. •J.M. Gonz´alez-Cantero, E. Graciani, B. L´opez-Romano, and F. Par´ıs. Induced unfolding: a new mechanism of unfolding failure in composite laminates. Submitted to Composites Science and Technology, 2017. •J.M. Gonz´alez-Cantero, E. Graciani, B. L´opez-Romano, and F. Par´ıs. Enhanced semi-analytic model for the evaluation of non-regularized stresses in the unfolding failure by using Legendre polynomials. In preparation. •J.M. Gonz´alez-Cantero. 3D analytical model for evaluating stresses in curved aerostructural composite laminates. ECSE2016-1, Airbus Group Chair of Aeronautical Studies, School of Engineering, University of Seville, 2016. Conferences: •J.M. Gonz´alez-Cantero, E. Graciani, A. Bl´azquez, and F. Par´ıs. Analytic evaluation of radial stresses in unfolding failure of composite materials. Comparison with numerical solutions. In Proceedings of the 16th European Conference on Composite Materials ECCM16, 2014. 22-26 of June, Seville, Spain. •J.M. Gonz´alez-Cantero, E. Graciani, F. Par´ıs, B. L´opez-Romano, and J.A. Rodr´ıguez-S´anchez. Study of the composite ILTS test in CFRP angles by means of analytical and FEM models. In 7th International Conference on Composites Testing and Model Identification, 2015. IMDEA, Getafe, Spain. •J.M. Gonz´alez-Cantero, E. Graciani, F. Par´ıs, D. Meizoso-Latova, and J.A. Rodr´ıguez-S´anchez. Resultados semi-anal´ıticos de tensiones no regularizadas en barras curvas de material compuesto usando un desarrollo en serie de funciones monomiales. In Proceedings of the XI Congreso Nacional de Materiales Compuestos, Matcomp’15, 2015. 6-8 of July, M´ostoles, Spain. 277 Curriculum Vitae •J.M. Gonz´alez-Cantero, E. Graciani, F. Par´ıs, B. L´opez-Romano, and D. Meizoso-Latova. Semi-analytic solution of non-regularized unfolding stresses in composite beams employing a series approximation based on Legendre polynomials. In Proceedings of the 21th International Conference on Composite Materials ICCM20, 2015. 19-24 of July, Copenhaguen, Denmark. •J.M. Gonz´alez-Cantero, E. Graciani, F. Par´ıs, and B. L´opez-Romano. Stacking sequence optimization of curved UD-CFRP laminates for improving unfolding strength considering thermal residual stresses. In Proceedings of the 17th European Conference on Composite Materials ECCM17, 2016. 26-30 of June, Munich, Germany. •J.M. Gonz´alez-Cantero, E. Graciani, F. Par´ıs, and B. L´opez-Romano. Semi-analytical stress calculation in composite materials using a 2D series expansion of the displacements. In Proceedings of the 19th International Conference on Composite Structures ICCS19, 2016. 5-9 of September, Porto, Portugal.