Closed-form methodology for the structural analysis of stiffened composite plates with cutouts and non-uniform lay-up

Composite Structures 343 (2024) 118284 Available online 14 June 2024 0263-8223/© 2024 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Composite Structures journal homepage: www.elsevier.com/locate/compstruct Closed-form methodology for the structural analysis of stiffened composite plates with cutouts and non-uniform lay-up A. Blázquez a,∗, D. Pastorino b, B. López-Romano c, F. París a aGrupo de Elasticidad y Resistencia de Materiales, Escuela Técnica Superior de Ingeniería, Universidad de Sevilla, Camino de los Descubrimientos s/n, Sevilla, 41092, Spain bAirbus Operations, S.L., Getafe, 28906, Spain cFIDAMC, Foundation for the Research, Development and Application of Composite Materials, Avda. Rita Levi Montalcini 29, 28906 Getafe, Madrid, Spain ARTICLE INFO Keywords: Anisotropy Composites Stress analysis Stiffened plate Cutout Closed-form Lekhnitskii formalism ABSTRACT Building upon the methodology developed in prior studies, which utilized the Lekhnitskii formalism to address composite plates featuring cutouts, varying thicknesses, and different stacking sequences, the present work introduces the incorporation of stiffeners. The structural configuration comprises plates and stringers, with each stringer positioned between two plate regions and subjected to both stretching and bending loads. It is important to note that no coupling between the bending and stretching responses is accounted for; therefore, lay-ups must be symmetric, and stringers must be embedded and bisymmetric. The stringers are modeled using the Euler–Bernoulli formulation with the free torsion hypothesis, while the plates are modeled using the Kirchhoff–Love formulation. Several benchmark problems are analyzed, and the results are compared with those obtained using finite element analysis (utilizing Abaqus software), demonstrating a satisfactory agreement while also showcasing competitive computational efficiency. Thus, the present methodology provides the industry with a novel tool that enables efficient parametric analysis and facilitates the most promising configurations during the initial phases of the design to be selected by engineers. 1. Introduction In many engineering applications, it is usual to find structural elements formed by a combination of beams and shells. They are called stiffened panels resulting in a significant weight reduction of structures. Shell component is usually identified as the skin, and the beams as the stringers (see Fig. 1). The external structures of the wings and fuselage of aircrafts or the structure of a ship’s hull are significant examples of structures constituted by this type of elements. Although the final design is typically verified using sophisticated models, often based on Finite Element Analysis (FEA), numerous decisions need to be made in the early design stages. In such situations, where certain data may be unknown, it becomes attractive to utilize alternative tools that enable quick estimations of the effects of several alternatives, or the development of some parametric analysis. It is during these initial stages that analytical or pseudo-analytical formulations prove to be highly beneficial. In this context, the authors developed a Lekhnitskii-based closedform methodology for analyzing composite plates with elliptical cutouts in Pastorino et al. [1,2,3,4] for stretching, and in Pastorino et al. [5] for bending problems. Building upon these previous works, this ∗Corresponding author. E-mail address: [email protected] (A. Blázquez). paper extends the methodology to the analysis of stiffened flat panels. Specifically, it focuses on the case of embedded stiffeners and symmetric laminates, thus maintaining the decoupling between membrane and bending effects. The inclusion of one-dimensional stringers into this methodology is a clear enhancement that increases its range of applicability. The plate model described in Pastorino et al. [5] adopts the Kirchhoff–Love plate theory. Coherently, the beam model considered here will follow the Euler–Bernoulli beam theory, which is the simplest model for a beam. This model is based on the assumption that the transversal section does not deform in its own plane and remains plane and perpendicular to the deformed axis of the beam. Despite its relative simplicity, the application of this model to anisotropic beams involves some difficulties in the estimation of the equivalent elastic modulus and the bending stiffness. In this regard, Vinson and Sierakowski [6], Omri and Vladimir [7], Bauchau and Craig [8], Kassapoglou [9] and Librescu and Song [10] should be considered. In the frame of analytical or semi-analytical methodologies, very few investigations link one-dimensional formulations together with two-dimensional formulations for plates with or without cutouts. https://doi.org/10.1016/j.compstruct.2024.118284 Received 18 March 2024; Received in revised form 27 May 2024; Accepted 12 June 2024 Composite Structures 343 (2024) 118284 2 A. Blázquez et al. Fig. 1. Thermoplastic composite panel stiffened by five T-section longitudinal stringers. In Oden and Ripperger [11] it is assumed a shear lag mechanism, the skin developing mainly shear stress and the stringer axial stress. For stiffened flat panels, Lekhnitskii [12] considers a homogeneous equivalent plate modifying the stiffness properties of the structure. Similarly, Kassapoglou [9] provides a method for considering stiffened skins by smearing the stiffness properties of both the skin and the stringers in order to consider an equivalent section stiffness. In previous work, Kassapoglou [13] develops an energy-based closed-form methodology for the estimation of the stresses at the skin–stringer interface. However, these methodologies are rough approximations based on very simplified assumptions that provide a general and imprecise estimate of the solution to the problem, suitable mainly for obtaining qualitative and overall results, but not for fine details or exact precision. In contrast, the methodology proposed here identifies each stiffener and each area of the plate individually, providing very satisfactory results of the displacements and the section force distributions in each of them. 2. Definition of the generic stiffened plate under study To the best of the authors’ knowledge, there is no closed-form methodologies in the peer-reviewed literature that couples the formulation of anisotropic plates with cutouts and a one-dimensional beam under axial, bending, and torsional loads. Hence, the contents of this article will result in an innovative contribution to the current state of the art. The formulation here presented considers the following hypotheses: 1. Geometry is comprised of a set of plate regions and stringers so that each stringer is always placed between two plate regions. Each plate region can have a cutout or none. 2. The structures might be subjected to stretching (or membrane), and bending loads, but no transverse load (pressure). 3. One-dimensional beam is considered based upon the Euler– Bernoulli formulation and free torsion hypothesis, warping being negligible; and two-dimensional plate structures are based upon the Kirchhoff–Love formulation. 4. There is no bending–stretching coupled response. Hence, the stringers are straight, the shear center and the barycenter of the beams coincide each other (bi-symmetrical section), and they are placed in the middle plane of the plates which is a symmetric material plane. Fig. 2 illustrates a generic example of the problem under consideration. The case shown in the figure comprises five regions and two stringers. Each component has its own set of axes, with all of them featuring axes 𝐱3perpendicular to the mid-plane of the structure and pointing in the same direction. Laminated sequences may vary between regions, and different sections of the beam can be considered for the stringers. Fig. 2. Generic example of a plate with cutouts and embedded stringers. The following section is dedicated to recalling the Euler–Bernoulli beam model, which is considered for each stringer, and the Kirchhoff– Love plate theory, which is considered for each subplate (with or without cutouts). This is with the main purpose of showing the notation used by Bauchau and Craig [8] which is employed here. Next, the closeform methodology presented for plates with cutouts in the authors’ previous works, Pastorino et al. [4],Pastorino [14], are extended to stiffened plates, paying special attention to the coupling between beam and plate models. Some benchmark problems are considered for validation. Finally, a section is dedicated to the conclusions and subsequent developments. 3. Beam formulation Let us define a local coordinate system located at the barycenter of each section. The 𝐱1axis aligns with the undeformed beam axis, while the 𝐱3axis is perpendicular to the plane of the stiffened plate. The 𝐱2axis is oriented to form a rectangular right-hand system. The combination of the Euler–Bernoulli and the free torsion hypotheses implies that the cross-section is rigid in its own plane and it remains plane and normal to the beam axis after deformation. Thus four degrees of freedom are used to describe the kinematic of the beam: 𝑢1,𝑢2,𝑢3, 𝜙1. Overline marks indicate that the variable is defined at the section level, which means it is a function of 𝑥1coordinate only (dependence is omitted for concision). In the beam model, assuming a homogeneous cross-section, 𝜎11 = 𝐸1𝜀11, where 𝐸1is the elastic modulus of the material in the axial direction of the beam; and 𝜏1𝑠=𝐺1𝑠𝛾1𝑠, with 𝑠= 2,3denoting the direction of tangential stresses, and 𝐺1𝑠representing the corresponding shear modulus of the material. The relationships between stresses and strains become more intricate for a composite laminate. When considering thin-walled sections, the direction 𝑠follows the mid-line of the wall, let 𝐸(𝑘) 1and 𝐺(𝑘) 1𝑠represent, respectively, the elastic modulus in the axial direction and the shear modulus in the plane 1𝑠of the 𝑘th layer, as proposed by Kassapoglou [9]: 𝐸(𝑘) 1=1 ℎ(𝑘)𝑎(𝑘) 11 𝐺(𝑘) 1𝑠=1 ℎ(𝑘)𝑎(𝑘) 66 (1) where ℎ(𝑘)is the thickness of the ply, and 𝑎(𝑘) 11 and 𝑎(𝑘) 66 are, respectively, the first and sixth diagonal terms of the inverse of the extensional stiffness matrix. Considering this, it is appropriate to redefine the barycenter based on a definition weighted by the elastic modulus. In other words, the origin of the local coordinate system is that point for which the following conditions hold: 𝐾 ∑ 𝑘=1 (∫𝑘 𝑥2𝐸(𝑘) 1d𝑎)= 0 (2a) 𝐾 ∑ 𝑘=1 (∫𝑘 𝑥3𝐸(𝑘) 1d𝑎)= 0 (2b) Composite Structures 343 (2024) 118284 3 A. Blázquez et al. where 𝐾is the number of layers of the laminate, and 𝑘is the cross-sectional area of the 𝑘ply. Following the nomenclature proposed by Bauchau and Craig [8], and considering that there is no bending coupling, the section forces become:  𝑁1= 𝐾 ∑ 𝑘=1 ∫𝑘 𝜎11d𝑎=𝑆0𝜀1(3a)  𝑀2= 𝐾 ∑ 𝑘=1 ∫𝑘 𝑥3𝜎11d𝑎=𝐻𝑐 22𝜅2(3b)  𝑀3= 𝐾 ∑ 𝑘=1 ∫𝑘 (−𝑥2)𝜎11d𝑎=𝐻𝑐 33𝜅3(3c)  𝑀1= 𝐾 ∑ 𝑘=1 ∫𝑘(𝜏13𝑥2−𝜏12𝑥3)d𝑎=𝐻𝑐 11𝜅1(3d) where 𝑆0= 𝐾 ∑ 𝑘=1 (∫𝑘 𝐸(𝑘) 1d𝑎)(4a) 𝐻𝑐 22 = 𝐾 ∑ 𝑘=1 (∫𝑘 𝑥2 3𝐸(𝑘) 1d𝑎)(4b) 𝐻𝑐 33 = 𝐾 ∑ 𝑘=1 (∫𝑘 𝑥2 2𝐸(𝑘) 1d𝑎)(4c) 𝐻𝑐 11 =𝐽 𝐴 𝐾 ∑ 𝑘=1 (∫𝑘 𝐺(𝑘) 1𝑠d𝑎)(4d) The parameter 𝑆0represents the axial stiffness, 𝐻𝑐 22 and 𝐻𝑐 33 denote the bending stiffnesses about the 𝑥2and 𝑥3axes, respectively, and 𝐻𝑐 11 represents the torsion stiffness. Notably, in accordance with hypothesis 4 detailed in Section 2, the term corresponding to cross-bending stiffness (𝐻𝑐 23 in the notation used) is null, as there is no coupling between bending in the plane of the plate and in the plane perpendicular to the plate. The superscript 𝑐is employed to indicate that the bending stiffnesses are calculated at the barycenter. Expressions for torsion constants, 𝐽, corresponding to different section shapes can be found in [15], or for an open thin-walled section in [8]. The system is completed with the equilibrium equations where the external applied distributed beam loads: 𝑝1,𝑝2,𝑝3,𝑔1,𝑔2and 𝑔3appear. Loads 𝑝1,𝑝2and 𝑝3are distributed forces along each axes, and 𝑔1,𝑔2 and 𝑔3are distributed moments around each axes. After some basic substitutions, the here-considered beam model is governed by the following equations: 𝑆0 d2𝑢1 d𝑥2 1 +𝑝1= 0 (5a) 𝐻𝑐 33 d4𝑢2 d𝑥4 1 +𝑞3= 0 (5b) 𝐻𝑐 22 d4𝑢3 d𝑥4 1 +𝑞2= 0 (5c) 𝐻𝑐 11 d2 𝜙1 d𝑥2 1 +𝑔1= 0 (5d) where the new external loads, 𝑞2and 𝑞3, are given by: 𝑞2= − 𝑝3−d𝑔2 d𝑥1 (6a) 𝑞3= − 𝑝2+d𝑔3 d𝑥1 (6b) Eqs. (5a)a and (5a)b belong to the stretching problems, and (5a)c and (5a)d belong to the bending problem, and they are decoupled each other in the case considered here. Fig. 3. Component geometry and cutout parameters definition. In order to approximate the kinematics of the beam, although there are many possibilities, polynomial series are selected in this work: 𝑢1(𝑠) = 𝑓10 + 𝑁𝑏 ∑ 𝑛=1 (𝑓1𝑛𝑠𝑛)(7a) 𝑢2(𝑠) = 𝑓20 + 𝑁𝑏 ∑ 𝑛=1 (𝑓2𝑛𝑠𝑛)(7b) 𝑢3(𝑠) = 𝑓30 + 𝑁𝑏 ∑ 𝑛=1 (𝑓3𝑛𝑠𝑛)(7c)  𝜙1(𝑠) = 𝑓40 + 𝑁𝑏 ∑ 𝑛=1 (𝑓4𝑛𝑠𝑛)(7d) Here, 𝑠represents a coordinate along the bar axis, while 𝑁𝑏denotes the order of the approximation, and the coefficients of the approximation function are denoted as 𝑓𝑖𝑗 with 𝑖= 1,…,4and 𝑗= 0,1,…, 𝑁𝑏. In Eqs. (7), the same order is applied to all movements; however, it is not mandatory to use the same order for all kinematic variables. 4. Kirchhoff–Love plate theory The formulation of the plates has been presented in previous articles addressing membrane [4] and bending [5] problems. In essence, the approach involves dividing the plate into regions, each featuring either one elliptical cutout or none. Each subplate must be composed of the same material (laminate setup), although the material may differ from one region to another. The incorporation of the Lekhnitskii formalism for infinite plates with cutouts [12], along with the utilization of Ogonowski’s formulation [16] for considering finite sizes on each subplate, and the incorporation of corresponding boundary conditions along the interfaces, completes the problem. While the specific details are available in [4,5], a summary is provided below for completeness. Each plate is defined by an arbitrary external shape (𝛤𝐸) and an elliptically shaped internal cutout (𝛤𝐶) as shown in Fig. 3. The ellipse major and minor semi-axes are denoted by 𝑎and 𝑏. A local coordinate system (𝐱1,𝐱2)is defined by the ellipse axes. Both 𝛤𝐸and 𝛤𝐶contours are covered by the arc-length 𝑠, which is directed counter-clockwise at the internal and clockwise at the external boundary. Tangential (𝐬) and normal (𝐧) directions comprise a set of coordinates that follow the boundary, 𝐧being always directed outward from the component. According to the Classical Laminate Theory (CLT) hypotheses, the segment perpendicular to the mid-plane of the plate remains straight, undeformed, and perpendicular to the deformed mid-surface. Thus, the kinematics of the plate can be defined by 3 degrees of freedom: 𝑢1,𝑢2, and 𝑢3. Composite Structures 343 (2024) 118284 4 A. Blázquez et al. Hat marks indicate that the variable is defined at the thickness segment level, which means it is a function of 𝑥1and 𝑥2coordinates (dependence is omitted for concision). Strains can be obtained directly from the small strains 𝜀−𝑢equations, and stresses from the constitutive law for each layer. Section forces are obtained by integrating the stresses:  𝑁1=∫ℎ∕2 −ℎ∕2 𝜎11d𝑥3=𝐴11 𝜀1+𝐴12 𝜀2+𝐴13 𝛾12 (8a)  𝑁2=∫ℎ∕2 −ℎ∕2 𝜎22d𝑥3=𝐴21 𝜀1+𝐴22 𝜀2+𝐴23 𝛾12 (8b)  𝑁12 =∫ℎ∕2 −ℎ∕2 𝜎12d𝑥3=𝐴31 𝜀1+𝐴32 𝜀2+𝐴33 𝛾12 (8c)  𝑀1=∫ℎ∕2 −ℎ∕2 𝑥3𝜎11d𝑥3=𝐷11 𝜅1+𝐷12 𝜅2+𝐷13 𝜅12 (8d)  𝑀2=∫ℎ∕2 −ℎ∕2 𝑥3𝜎22d𝑥3=𝐷21 𝜅1+𝐷22 𝜅2+𝐷23 𝜅12 (8e)  𝑀12 =∫ℎ∕2 −ℎ∕2 𝑥3𝜎12d𝑥3=𝐷31 𝜅1+𝐷32 𝜅2+𝐷33 𝜅12 (8f) The variables 𝜀1and 𝜀2are the unit elongations along 𝑥1and 𝑥2 axis respectively, and with 𝛾12 are variables of the stretching problem. The variables 𝜅1and 𝜅2correspond to the curvatures of the deformed surface of the plate in 13 and 23 planes respectively, and with the cross curvature 𝜅12 are variables of the bending problem. Here, 𝐴𝑖𝑗 and 𝐷𝑖𝑗 are components of the stiffness matrix of the plate, [17]. It is noteworthy that no coupling terms (𝐵𝑖𝑗 ) are considered in this work, assuming hypothesis 4 outlined in Section 2. The system is complemented by the Equilibrium equations that for plates under membrane loads can be rewritten in terms of an Airy stress function 𝜑, Lekhnitskii [12]. 𝑎22 𝜕4𝜑 𝜕𝑥4 1 − 2𝑎26 𝜕4𝜑 𝜕𝑥3 1𝜕𝑥2 +(2𝑎12 +𝑎66)𝜕4𝜑 𝜕𝑥2 1𝜕𝑥2 2 + −2𝑎16 𝜕4𝜑 𝜕𝑥1𝜕𝑥3 2 +𝑎11 𝜕4𝜑 𝜕𝑥4 2 = 0 (9) Section forces  𝑁1, 𝑁2and  𝑁12 are obtained from 𝜑:  𝑁1=𝜕2𝜑 𝜕𝑥2 2  𝑁2=𝜕2𝜑 𝜕𝑥2 1  𝑁12 = − 𝜕2𝜑 𝜕𝑥1𝜕𝑥2 (10) While the equilibrium equation in the transverse direction was developed in Pastorino et al. [5] in terms of the deflection 𝑤(notice that 𝑤=𝑢3) as follows, 𝐷11 𝜕4𝑤 𝜕𝑥4 1 + 4𝐷16 𝜕4𝑤 𝜕𝑥3 1𝜕𝑥2 + 2 (𝐷12 + 2𝐷66)𝜕4𝑤 𝜕𝑥2 1𝜕𝑥2 2 + +4𝐷26 𝜕4𝑤 𝜕𝑥1𝜕𝑥3 2 +𝐷22 𝜕4𝑤 𝜕𝑥4 2 = 0 (11) Bending problem section forces  𝑀1, 𝑀2and  𝑀12 are obtained from Eqs. (8d) to (8f), and  𝑄1and  𝑄2from the equations of equilibrium of bending moments. Both (9) and (11) are fourth-order homogeneous partial derivative equations, that can be expressed in the form: 𝐴1 𝜕4𝑊 𝜕𝑥4 1 +𝐴2 𝜕4𝑊 𝜕𝑥3 1𝜕𝑥2 +𝐴3 𝜕4𝑊 𝜕𝑥2 1𝜕𝑥2 2 + +𝐴4 𝜕4𝑊 𝜕𝑥1𝜕𝑥3 2 +𝐴5 𝜕4𝑊 𝜕𝑥4 2 = 0.(12) 𝑊being 𝜑for the stretching problem or 𝑤for the bending problem, with the proper 𝐴𝑖,𝑖= 1,…,5, that can be easily deduced. This differential equation can be integrated using four first-order operators (see [4,5,12]): 𝐹1(𝐹2(𝐹3(𝐹4(𝑊))))= 0,(13) 𝐹𝑗being: 𝐹𝑗=𝜕 𝜕𝑥2 −𝜇𝑗 𝜕 𝜕𝑥1 ,(14) where 𝜇𝑗are the roots of the characteristic equation: 𝐴5𝜇4+𝐴4𝜇3+𝐴3𝜇2+𝐴2𝜇+𝐴1= 0.(15) Four roots, 𝜇𝑗with 𝑗= 1,…,4, are obtained from (15), which are complex and conjugated in pairs [12]. Let us ordering them in the way that 𝜇3=𝜇1and 𝜇4=𝜇2. Following Lekhnitskii [12], integrating with the operators one at a time on the basis of Eq. (13), 𝐹4(𝜑) = 𝐺3, 𝐹3(𝐺3) = 𝐺2, 𝐹2(𝐺2) = 𝐺1, 𝐹1=𝜕𝐺1 𝜕𝑥2 −𝜇𝑗 𝜕𝐺1 𝜕𝑥1 ,(16) From the last equation of (16) it follows that: 𝐺1=𝑓1(𝑥1+𝜇1𝑥2)(17) where 𝑓1is any arbitrary function dependent on 𝑥1+𝜇1𝑥2. The stepby-step integration of Eqs. (16) results in the general solution of (12), see [12,18], as the sum of four complex function 𝑊𝑗(𝑧𝑗), being 𝑧𝑗= 𝑥1+𝜇𝑗𝑥2for 𝑗= 1,…,4and being 𝑧3=𝑧1and 𝑧4=𝑧2: 𝑊(𝑥1, 𝑥2) = 4 ∑ 𝑗=1 𝑊𝑗(𝑧𝑗)(18) Following the procedure developed by Lekhnitskii [12], the following conformal mapping: 𝜁𝑗(𝑧𝑗) = 𝑧𝑗±√𝑧2 𝑗−𝑎2−𝑏2𝜇2 𝑗 𝑎−𝑖𝑏𝜇𝑗 .(19) transforms an elliptic hole of semiaxis lengths 𝑎and 𝑏in the real plane, (𝑥1, 𝑥2), into a unitary circular one in the transformed complex plane, 𝜁𝑗. It has to be mentioned that the selection of the sign ahead of the square root is not straight forward and the algorithm developed by Koussios [19] is used. Now, based on the Lekhnitskii formalism, the formulation is redefined using the derivatives of the functions 𝑊𝑗(𝑧𝑗). Let us call them 𝛺𝑗(𝑧𝑗) = d𝑊𝑗∕d𝑧𝑗,𝑗= 1,…,4, and define 𝛺(𝑥1, 𝑥2)as the sum of the 𝛺𝑗(𝑧𝑗). Taken into account that 𝛺3(𝑧3) =  𝛺1(𝑧1)and 𝛺4(𝑧4) =  𝛺2(𝑧2), it results that 𝛺(𝑥1, 𝑥2)is a real function: 𝛺(𝑥1, 𝑥2) = 4 ∑ 𝑗=1 𝛺𝑗(𝑧𝑗)=2Re [2 ∑ 𝑗=1 𝛺𝑗(𝑧𝑗)](20) This formulation allows a straightforward definition of the boundary conditions. However, it precludes using the boundary conditions in terms of 𝑊. This is unimportant for the stretching problem, but not for the bending problem, where the boundary conditions in deflection must be replaced by another one more or less equivalent condition affecting the slope at the edges in both normal and tangential directions. This is a conventional approach to the problem in the existing literature [12, 20,21]. To consider that the plate is finite the work of Ogonowski [16] is followed, using Laurent series to approximate 𝛺𝑗, [4,5,22,23]: 𝛺𝑗(𝑧𝑗) = 𝐶0𝑗+(𝐶′ 0𝑗+𝐶′′ 0𝑗𝑧𝑗)ln(𝜁𝑗)+ + 𝑁 ∑ 𝑛=1 (𝐶𝑛𝑗 𝜁−𝑛 𝑗+𝐶∗ 𝑛𝑗 𝜁𝑛 𝑗)(21) where 𝐶0𝑗,𝐶′ 0𝑗,𝐶′′ 0𝑗,𝐶𝑛𝑗 , and 𝐶∗ 𝑛𝑗 are unknown coefficients determined by imposing the boundary conditions, and 𝑁represents the number of approximation functions to be used. It is important to note that the term 𝐶0𝑗solely influences displacements, not forces. Additionally, the logarithmic term is associated with unbalanced forces (𝐶′ 0𝑗) and moments (𝐶′′ 0𝑗) at the cutout [12,24], (notice that 𝐶′′ 0𝑗= 0 for the stretching problem). The logarithmic terms and terms with negative exponent are omitted in plates without a cutout. Composite Structures 343 (2024) 118284 5 A. Blázquez et al. Fig. 4. Force and displacement diagrams for the stretching problem. 5. Boundary conditions Boundary conditions along the cutout boundaries, external boundaries, and shared boundaries between two regions are extensively discussed in [4,5]. They are imposed on sets of collocation points distributed along the boundaries of each subplate, producing a system of equations to obtain the coefficients 𝐶0𝑗,𝐶′ 0𝑗,𝐶′′ 0𝑗,𝐶𝑛𝑗 , and 𝐶∗ 𝑛𝑗 for each area. Here, we provide a concise overview, thus focusing on the conditions established between plate regions and the intervening stringer. Typically, the conditions along the cutout boundaries establish an analytical relationship between 𝛺1and 𝛺2. As a result, no collocation points are needed, which directly enables the calculation of the coefficients of the logarithmic terms and reduces the number of unknown coefficients for the exponential terms by half. Specifically, for selfequilibrated loads on 𝛤𝐶, it implies 𝐶′ 0𝑗= 0 with 𝑗= 1,…,4for the stretching problem (𝛺=d𝜑∕d𝑧𝑗), and 𝐶′ 0𝑗=𝐶′′ 0𝑗= 0 with 𝑗= 1,…,4 for the bending problem (𝛺=d𝑤∕d𝑧𝑗). Along the external boundaries, conditions involve imposing a proper combination of resultant forces or translations, and moments or rotations of the thickness line at any point (𝑠) along 𝛤𝐸. This is achieved using collocation points. Shared boundaries between regions require the imposition of the principle of action–reaction between section forces and the continuity of displacements and rotations. Collocation points are also utilized for this purpose. The number of collocation points placed along external and shared boundaries has to be selected to ensure that the number of equations in each boundary exceeds the number of unknowns, the resulting system of equations being solved using the least squared method. The number of collocation points is a parameter of the proposed procedure. Figs. 4 and 5illustrate the force and moment diagrams on one side, and the displacement diagram on the other side, for establishing the conditions between plate regions and the intervening stringer in stretching and bending problems, respectively. The reference system considered is oriented parallel to the beam system, i.e., 𝐬=𝐱1,𝐧=𝐱2, and 𝐳=𝐱3, with subplate 1 located at the positive 𝑥2coordinate and subplate 2 at the negative 𝑥2coordinate. It is worth noting that, considering the beam as a one-dimensional element, the distributed moments 𝑔𝑛and 𝑔𝑧are null. The distributed loads transferred from the subplate regions to the beam are as follows (coordinate dependencies are omitted for conciseness): 𝑝𝑠= 𝑁(1) 𝑠𝑛 − 𝑁(2) 𝑠𝑛 (22a) 𝑝𝑛= 𝑁(1) 𝑛− 𝑁(2) 𝑛(22b) 𝑝𝑧= 𝑉(1) 𝑛− 𝑉(2) 𝑛(22c) Fig. 5. Force and displacement diagrams for the bending problem. 𝑔𝑠= 𝑀(1) 𝑛− 𝑀(2) 𝑛(22d) As usual in plate problems,  𝑉𝑛is the effective transverse shear force:  𝑉𝑛= 𝑄𝑛+𝜕 𝑀𝑠𝑛∕𝜕𝑠. The compatibility condition implies the identification of all the movements of the plates and the beam along the beam axis. It should be noted that the equations involving the displacement 𝑢𝑧, which represents the deflection of the plate and the beam, need to be reconsidered due to the bending problem being formulated in terms of the derivative of the deflection. To address this, a manipulation, analogous to [5], based on the derivative with respect to 𝑠is employed. While some information may be lost, it can be easily proven that it does not affect the stress and resultants, nor the load transfer between the plates and the beam. Additionally, assuming negligible shear deformations, the conditions between rotations around the beam axis should be expressed in terms of the derivative of the deflection. 𝑢(1) 𝑠=𝑢(2) 𝑠=𝑢𝑠(23a) 𝑢(1) 𝑛=𝑢(2) 𝑛=𝑢𝑛(23b) 𝑢(1) 𝑧=𝑢(2) 𝑧=𝑢𝑧⟹𝜕 𝑢(1) 𝑧 𝜕𝑠 =𝜕 𝑢(2) 𝑧 𝜕𝑠 =𝜕 𝑢𝑧 𝜕𝑠 (23c)  𝜙(1) 𝑠= 𝜙(2) 𝑠= 𝜙𝑠⟹𝜕 𝑢(1) 𝑧 𝜕𝑛 =𝜕 𝑢(2) 𝑧 𝜕𝑛 = 𝜙𝑠(23d) Therefore, it is necessary to enforce twelve equations along the stringer. This is accomplished by employing a set of collocation points distributed along the stringer, and the resultant system is solved using the least squares method. Preliminary tests have demonstrated that a logarithmic distribution of collocation points, with increased density near the corners, is preferable to a uniform one (as is used for other boundaries). Special attention should be directed to the endpoints of the beam, where two collocation points must be considered, one at each end, to apply correctly the conditions at the extremes of the beam. 6. Model verification - benchmark The closed-form formulation developed for the analysis of anisotropic stiffened plates has been implemented by using the mathematical software Matlab 2018a, under the aforementioned hypotheses. The methodology developed is validated in this section by comparing its results with those obtained by finite element (FE) analysis. To eliminate other possible sources of discrepancies, FE models incorporate beam elements that adhere to the Euler–Bernoulli hypothesis and plate elements that adhere to the Kirchhoff–Love hypothesis in the problems of plate bending. Specifically, in the Abaqus 6.14 software, which is the Composite Structures 343 (2024) 118284 6 A. Blázquez et al. Table 1 Mechanical properties of a lamina. 𝐸11 𝐸22 𝐺12 𝜈12 Thickness 140 GPa 9 GPa 4.65 GPa 0.3 0.184 mm Table 2 Lay-ups used in the benchmark exercises. Ref. Lay-up LU1 [45,−45,0,90,0]2S LU2 [45,−45,0,90,0,0,−45,45]2S LU3 [45,−45,0,90,0,0,−45,45,0]2S LU4 [45,−45,0,90,0,0,−45,45,45,−45,0,−45,45]2S LU5 [−45,45,0,90,0]S LU6 [45,−45,−45,45,0,90]2S Fig. 6. Geometry, loads and lay-up for exercise 1. one used in this study, these elements are referred to as B33 (linear line, 2 nodes, and 6 degrees of freedom per node) for beams, and STRI3 (linear triangular, 3 nodes and 5 degrees of freedom per node) or STRI65 (quadratic quadrilateral, 6 nodes, and 5 degrees of freedom per node) for plates [25]. Nevertheless, for stretching problems, any element of the standard library can be considered, for example, S8R (quadratic quadrilateral, 8 nodes and 6 degrees of freedom per node). Table 1 presents the mechanical properties of the material considered in the examples, and Table 2 collects the lay-ups. Initially, a series of tests should be conducted to identify the optimal combination of parameters, ensuring both accurate solutions and efficiency. These parameters include the number of functions used to approximate the solution within each region of the plate (𝑁), the equation-unknown ratio (𝑅𝑢), or equivalently, the number of collocation points, and the number of functions utilized to approximate the solution on the stringer (𝑁𝑏). Previous studies dedicated to stretching [4] and bending [5] problems have established that a suitable balance is achieved with 𝑁= 20 and 𝑅𝑢= 2. The first benchmark exercise presented below is employed to select the number of functions to be used to describe the behavior of the stringer. The same combination is applied to the remaining exercises. 6.1. Benchmark exercise 1 The first exercise to be presented involves the structure depicted in Fig. 6. It comprises a plate with two circular cutouts and a stringer with I-section. The plate is divided into four regions, each featuring distinct lay-ups. The figure also illustrates the section and layup of the stringer. The structure is subjected to a combination of boundary conditions in the plane, with kinematic conditions in red, distributed loads in blue, and a concentrated load at one end of the stringer in green, resulting in a stretching problem. The combination of parameters 𝑁,𝑅𝑢, and 𝑁𝑏characterizes the solution process. According to Pastorino et al. [4,5], a satisfactory Fig. 7. Relative difference (𝑅𝑑) and computational time plotted against the truncation order for the beam (𝑁𝑏). balance between precision and time is achieved with 𝑁= 20 and 𝑅𝑢= 2. Consequently, for stiffened plates, these values are adopted. However, another parameter arises, namely, the number of functions to describe the behavior of the beam, denoted as 𝑁𝑏. To determine a suitable value, a series of tests are conducted by varying 𝑁𝑏from 2 to 40. It is important to note that this variation will also impact the number of collocation points along the beam axis, as determined by 𝑅𝑢. To evaluate the accuracy of the results, we use the parameter 𝑅𝑑, which is defined as the relative difference in axial force in the beam compared to the value obtained from a reference solution with 𝑁𝑏= 100: 𝑅𝑑= Max (|||𝑁1(𝑠) − 𝑁(ref) 1(𝑠)|||) |||𝑁(ref) 1(𝑠)||| (24) The evolution of the values of this parameter is shown in Fig. 7 as a function of the truncation order. The time of execution is also plotted to provide all the information needed for the selection of the optimum 𝑁𝑏. In conclusion, it can be inferred that a satisfactory level of accuracy (𝑅𝑑<1%) is achieved for 𝑁𝑏>10. Given that the increase in execution time is only modest, 𝑁𝑏= 40 (𝑅𝑑≈ 0.1%) is chosen and will be utilized in subsequent exercises. Fig. 8 illustrates the comparison of section force distribution in the plate regions obtained through the proposed approach with 𝑁= 20 and 𝑁𝑏= 40, and a FE model with 523626 degrees of freedom utilizing S8R and B33 elements. Additionally, Figs. 9 and 10 present the axial force and in-plane bending moment, respectively, along the stringer also compared with the results from the FE model. A high degree of agreement is observed, affirming the robustness of the proposed procedure. Concerning CPU time, the closed-form solution under the mentioned conditions required 8.6 s, whereas the FEM solution took 29 s when executed on an Intel Core i5-8250U CPU 1.60 GHz with 8 GB DDR3 RAM. 6.2. Benchmark exercise 2 The second benchmark exercise is illustrated in Fig. 11. It comprises a stiffened plate with three regions featuring distinct composite layups, each incorporating a circular cutout, and two stringers with a circular tubular section. The applied loads on the plate are depicted in blue, and the support conditions are highlighted in red. Additionally, concentrated moments at the ends of the stringers are indicated in green. It is important to emphasize that the simple support conditions (𝑢3= 0) have been replaced by the condition affecting the derivative of deflection along the boundary, implying that the deflection is Composite Structures 343 (2024) 118284 7 A. Blázquez et al. Fig. 8. Section forces for the example problem 1. Fig. 9. Axial force along the stringer for example 1. constant along the boundary but not necessarily zero, as explained further by Lekhnitskii [12] and Pastorino et al. [5]. These combined actions result in a bending state within the plate regions and induce out-of-plane bending and torsion in the stringers. The problem is solved using 𝑁= 20,𝑅𝑢= 2 and 𝑁𝑏= 40, while the FE model used for comparison used a STRI3 and B33 element mesh with 304,818 degrees of freedom. Bending result distributions on the plate regions are compared in Fig. 12, and the bending and torsion moments in the stringers in Figs. 13 and 14. Remarkably, an excellent agreement is achieved between the presented methodology and the finite element solution. However, a slight discrepancy is observed in the bending and torsional moments along Fig. 10. In-plane bending moment along the stringer for example 1. Fig. 11. Geometry, loads and lay-up for exercise 2. the stringers. In light of this difference, additional tests have been conducted, revealing that an improvement in result accuracy is more effectively attained by increasing the truncation order of the plate regions (𝑁), as opposed to increasing 𝑁𝑏. The rationale behind this observation is that the transfer of load at the joints between the plates and the beams is influenced by both 𝑁and 𝑁𝑏. Simply increasing one of them does not constitute the optimal strategy for achieving more satisfactory results. Due to space limitations, detailed results are not presented here. The CPU time for the closed-form solution was 7.43 s, while the FEM solution required 27 s on an Intel Core i5-8250U CPU 1.60 GHz with 8 GB DDR3 RAM. 6.3. Benchmark exercise 3 The geometry of the third and latest benchmark exercise presented here is depicted in Fig. 15. It features a trapezoidal stiffened plate comprising six plate regions with two cutouts and two stringers with an I-section. Boundary conditions are defined in terms of loads (colored in blue), displacements (colored in red), and two concentrated forces at the stringers (colored in green). As in previous examples, the problem is solved using the methodology presented here with 𝑁= 20,𝑅𝑢= 2, and 𝑁𝑑= 40. The obtained solution is compared with that from a finite element model based on S8R and B33 elements, featuring 517,666 degrees of freedom. Results for the plate regions are shown in Fig. 16, and along the stringers in Figs. 17 and 18. Once again, the agreement between both solutions is excellent, though some slight discrepancies are observed. Similar to benchmark example 2, the solution obtained with the proposed methodology will become closer to the finite element solution by increasing the number of functions used for the approximation in both the plate regions (𝑁) and the beams (𝑁𝑏). In this exercise, the close-form solution required 14.86 s and the FEM consumed 70 s. Composite Structures 343 (2024) 118284 8 A. Blázquez et al. Fig. 12. Section forces for the example problem 2. Fig. 13. Bending moment along the stringers for example 2. 7. Discussion and conclusions The present paper introduces an innovative closed-form approach to solve the structural problem of composite stiffened plates with cutouts. The lay-ups in the plate regions are assumed to be symmetrical and the stringers embedded, thereby decoupling the stretching and bending problems. This development broadens the application range of closed-form methodologies based on the Lekhnitskii formalism previously presented by the authors in [4,5]. Although the applicability of the proposed procedure may be considered limited, it is not uncommon to use designs with symmetrical Fig. 14. Torsion moment along the stringers for example 2. Fig. 15. Geometry, loads and lay-up for exercise 3. Composite Structures 343 (2024) 118284 9 A. Blázquez et al. Fig. 16. Section forces for the example problem 3. Fig. 17. Axial force along the stringers for example 3. planar configurations featuring approximately elliptical cutouts, where the methodology proposed here could be directly applicable. In cases where the panel configuration is not symmetrical, and even for curved panels, the aerospace industry often uses simplified methods in the initial design stages to help engineering decisions. This is in line with the primary objective of the methodology here presented, enabling efficient parametric analysis, facilitating the selection of the most promising configurations. The methodology has been successfully validated by comparing the obtained solutions in plates and beams with results from FE models Fig. 18. In-plane bending moment along the stringers for example 3. using standard elements of Abaqus, which adhere to the Kirchhoff– Love hypothesis for plate bending and the Euler–Bernoulli hypothesis for beams, respectively. The computational time is quite competitive and the accuracy is outstanding even for low approximation orders. The most innovative findings of the methodology are: •Frame. The methodology is developed to analyze the structural behavior of stiffened plates with non-uniform lay-ups, including elliptical cutouts and embedded stringers. •Accuracy. The tool has been demonstrated to be able to be used with an outstanding level of accuracy for the exercises presented.