ON DIFFERENT SHELL MODELING IN NUMERICAL ANALYSES OF A COMPOSITE STIFFENED PANEL UNDER UNIFORM PRESSURE

This work deals with the evaluation of the different modeling options based on the FEM to reproduce the structural performance of a composite stiffened panel in postbuckling regime under uniform pressure. Two fundamental modeling options are used for this purpose: monolithic and multi-part approaches. These strategies are based respectively on standard and 7-parameter finite element shell formulations as underlying mechanical theories. The numerical results obtained by means of these alternatives are compared with the experimental response of the panel with special emphasis in the postbuckling evolution of the behavior of the structure, and incorporating the initial geometric imperfections of the panel in the simulations.


INTRODUCTION
Composite materials have the potential to produce highly efficient structures.In the last two decades they have been extensively used in many industrial applications as is the case of stiffened panels, which are ideally suited for fuselage and wing skin panels of aircrafts.Postbuckling designs of such composite stiffened panels allow us to significantly improve the performance of these components, leading to notable cost reductions.The structural collapse of these panels is usually characterized by the appearance of distinct delaminations and skinstringer debonding processes at several locations.Therefore, in order to fully exploit this postbuckling load carrying capacity, high fidelity numerical models are especially required [1], [6], [16].
From the structural point of view, an essential problem in dealing with composite laminates is an accurate prediction of interlaminar stresses, which play a significant role in the failure analysis of these specimens.However, the classical Kirchhoff-Love and Reissner-Mindlin shell models do not include the entire set of relevant stresses.On the computational side, low-order elements are often preferred due to their relative simplicity in implementation and to the lower numerical costs when compared to higher-order interpolations.Following these arguments, different improved shell formulations using low-order elements have been proposed in the literature in order to develop accurate models of thin-walled composite structures.Mixed finite shell element formulations have attracted especial attention in the last years, that can incorporate the mechanical information in the thickness direction [3], [13], [17] among others.This research activity has led to the development of several mechanical theories for the numerical treatment of composite structures in many engineering applications (see [15], [21], [22]).
In this work, the response of a composite stiffened panel subjected to uniform pressure load is analyzed experimentally and numerically through FEM.Different modeling strategies of this specimen based on standard and 7-parameter shell formulations were carried out.This allow reliable and efficient predictions of the stress state at most critical areas in the structure during the loading process to be obtained.
The outline of the present paper is as follows: in Section 2 the main characteristics of the composite stiffened panel under study and the experimental programme performed are introduced.Sections 3 and 4 present the main features of the 7-parameter shell model and the Shell-solid-shell approaches respectively.The resulting FE models using both modeling options are introduced in Section 5.The buckling analysis of the structure is briefly presented in Section 6. Finally Section 7 demonstrates the applicability of the two computational approaches mentioned above for the postbuckling analysis of a composite stiffened panel by comparing the numerical results with the experimental data.In addition, Section 7 shows that the use of these model allow the most critical regions of the specimen to be determined during the loading procedure.

Specimen details
Figure 1 shows a sketch of the specimen analyzed.It is a curved panel, whose main dimensions are: radius 1520 mm, arc length 615 mm and width 594 mm.There are two circumferential stringers, with Ω-transverse section, which are separated by a distance of 265 mm between their corresponding central lines.The panel is manufactured of carbon-epoxy laminates with unidirectional and fabric plies.A complete description of the laminate disposition at each of the entities that compose the structural set can be found in [23].The stringers were co-bonded to the skin of the panel, including an adhesive layer of 0.184 mm in thickness, whose mechanical properties are E = 3.216 × 10 9 Pa, ν = 0.3.

Test description and experimental results
The panel was subjected to uniform pressure by means of a combination of gas and liquid using a pressure box specifically designed for this experiment, see Figures 2 and 3. To register the mechanical response of the specimen during the experimental procedure, the panel was monitored using 28 strain gages and a displacement transducer LVDT [23].
Prior to testing the specimen, and in order to identify the initial imperfections, the surface of the specimen was measured along several circumferential lines by means of a com- parison clock, once it was settled down in the experimental device [23].This geometrical data will be incorporated in subsequent numerical simulations, see Section 7. Afterwards, several tests were conducted for loads below 3.5 times the lowest critical buckling load numerically predicted, see [23].The pressure was removed at the end of each test, checking that no evidence of damage was detected.Figure 4 depicts the experimental evolution of the displacement transducer LVDT located at two different locations on the skin of the panel, which are denoted as T-1 and T-2.Note the quite repetitive response of the structure in the different tests, pointing out again the absence of significant damage on the specimen up to this load level.A final test until panel collapse was carried out, which occurred at 0.45 bars.

BASIC RELATIONS OF THE 7-PARAMETER SHELL MODEL
In this section, the basic features regarding the 7-parameter shell model and the corresponding FE formulation are briefly reviewed.This shell model was originally developed in [14], see also [7], [9] and [13] for more detailed descriptions and discussions concerning this structural model.

Geometry and displacements approximation
According to Figure 5, convective coordinates are introduced 1 .The covariant tangent vectors are obtained by partial derivatives of the position vectors with respect to the curvilinear coordinates: A α amd a α being where X and x are the position vectors of an arbitrary point in the shell body, R and r denote the position vectors of the corresponding points on the shell midsurface, A 3 and a 3 (also denoted as D and d respectively) are the so-called director vector of the shell in both configurations, and h corresponds to the initial shell thickness.Thus, θ 1 and θ 2 denote the in-plane shell coordinates and θ 3 is the thickness coordinate.In the present shell model, the parametrization of the kinematic field of the continuum body is approximated by assuming a linear variation of the displacements in the thickness direction.Hence, the position vectors in the shell body in both configurations are expressed as: By definition, the kinematic field is expressed as the difference of the position vectors of an arbitrary point of the continuum 3D body in the current and in the reference configurations: u = x − X.Therefore, The parametrization of the displacements is decomposed into: three translational degrees of freedom v corresponding to the shell midsurface, and three difference degrees of freedom w which denote the relative displacement between the mid and the upper surfaces of the body.This difference vector is used for updating of the director vector of the shell along the deformation process (playing a similar role as the rotational degrees of freedom in the Kirchhoff-Love and in the Reissner-Mindlin formulations).It is worth mentioning that, as was deeply discussed in [14], this parametrization of the displacements incurs the so-called Poisson thickness locking.This intrinsic defect is alleviated through a hybrid FE formulation utilizing the Enhanced Assumed Strain (EAS) Method ( [26]) , see Section 3.2.1.

Variational basis
The FE formulation of the present shell model is based on the three-field Hu-Washizu functional, where the displacement, strain and stress fields are the independent tensorial quantities.In the EAS Method, following the approach proposed in [7], the extra strain field Ẽ supplements the compatible Green-Lagrange strain tensor E u (derived from the kinematic equation).Hence, the strain field finally results: E = E u + Ẽ.In a Total Lagrangian formulation, the Hu-Washizu functional yields, The first term of ( 6) describes the internal potential, W int 3D is the internal strain energy stored and S denotes the Second Piola-Kirchhoff stress tensor.Finally, Π ext corresponds to the external potential.
Note that the "compatible" part of the strain tensor E u can be expressed as: where F denotes the Deformation Gradient tensor, the Green-Lagrange strain tensor being referred to the curvilinear basis on the shell midsurface.Hence the components of E u in matrix notation can be split into the constant (α ij ) and the linear (β ij ) parts of the strain distribution across the thickness direction: Imposing the orthogonality condition between the independent stress field S and the enhanced strain field Ẽ: Therefore, the stress field is completely removed from the present formulation.Thus, the first variation of the functional ( 6) is obtained via the directional derivative concept [12].The resulting nonlinear set of equations is solved iteratively by means of the corresponding linearization procedure [7].Moreover, once the discretization process of the kinematic and the enhanced strain fields is accomplished, it is possible to preserve the pure displacement FE formulation by eliminating the enhanced strains using a standard static condensation scheme.

Bilinear interpolation of the displacements
First-order elements with linear interpolation of the displacement field in-plane direction are now presented.The discretization procedure of the structure results in bi-dimensional FE meshes, where the nodal locations correspond to each of the four corners of the element of the midsurface.Within the standard isoparametric concept, both the geometry of the elements and displacements are interpolated using the same bilinear shape functions, denoted by N, through the nodal coordinates X h and the nodal displacements d u , where the subscript h indicates the finite element approximation, N refers to the number of element nodes (4 in case of linear elements).The two standard bilinear shape functions N k (ξ, η) in the natural space are computed as: where ξ k , η k are the nodal coordinate values in the natural space of the element.

Locking treatment
As known from standard low-order finite element approximations, different locking pathologies appear leading to unrealistic (over-stiffened) numerical predictions.In this regard, several numerical procedures have been developed during the last three decades in order to remove this deficiency, where the Enhanced Assumed Strain (EAS) [26] and the Assumed Natural Strain (ANS) [4] Methods can be considered as those most widely employed for this purpose.In the present FE formulation, a combination of both computational techniques is used to reduce locking effects following the scheme proposed in [2].Specifically, the EAS Method is employed for the treatment of the Membrane and Poisson Thickness locking.The ANS Method is used to avoid artificial thickness strains (especially encountered in curved structures, a phenomenon known as Curvature Thickness locking) according to the procedure proposed in [5].Finally, to avoid Transverse Shear locking, a comparison between the EAS and the ANS Methods is carried out in this research.
By enhancing the twelve components of the strain resultants (α ij , β ij ) by means of the EAS Method, see equation (8), this results in 26 enhancing parameters (denoted as EAS26).Considering the ANS Method to alleviate Transverse Shear locking acting on the constant part of the transverse shear strain resultants, α 13 , α 23 in equation ( 8), the number of enhancing parameters is reduced to 22 (denoted as EAS22 ANS).To illustrate the combination of both procedures, see [8].

Extension to laminated structures
The extension of the present FE formulation to laminated structures is accomplished by means of the Equivalent Single Layer approach [15].Models based on this strategy are generated by replacing the heterogeneous laminate by a weighted average of the physical properties of each ply across the thickness.Particularly, this scheme is performed by conducting the numerical integration of the constitutive tensor along thickness direction using a Gauss integration scheme with two quadrature points per layer.Thus, the integration of the element stiffness matrix reads where C is the matrix of the constitutive tensor, µ denotes the shift tensor (which maps the shell body quantities onto the shell midsurface [7], ngauss is the number of quadrature points and N L is the number of layers.
The implementation of this laminated 7-parameter shell element is performed into the commercial FEA code ABAQUS using the user subroutine UEL.This element will denoted as Three-dimensional 7-parameter shell element and TShell in abbreviated form in what follows.
For testing the implementation, a representative simulation under geometrically nonlinear regime is presented in this research (see [25] for a complete description of the computational aspects and for an exhaustive testing activity on this element).The selected problem consists of a cantilever beam under distributed end forces.This classical application allows the out-of-plane bending performance of the elements implemented to be evaluated.It has been extensively employed as a popular nonlinear benchmark problem by a large number of authors, see [5] and [20] among others for isotropic material law.One mesh density of 8 × 1 (identified by Mesh 8) elements is employed, where the first number refers to the in-plane density whereas the second one denotes one element in the thickness direction.
The main characteristics of the application are depicted in Figure 6.In particular, the laminate has an stacking sequence [45/-45] S and the zero-degree direction coincides with to the x-axis, Figure 6.The orthotropic material properties selected for this application are also reported in Figure 6. Figure 6.Cantilever composite beam subjected to end forces.Geometry and material data Figure 7 shows the tip displacements, in x and z directions, which are denoted by U and W respectively, of the points 1 and 2 (see Figure 6).In this plot, the mechanical response of the structure using the aforementioned Mesh 8 is compared with the solution provided by the ABAQUS' S4R element.It is worth mentioning that the most satisfactory agreement between the S4R and the TShell 7-parameter elements is achieved when the combined usage of the ANS and the EAS Methods is employed for the numerical alleviation of locking effects.In addition, note that the TShell 7-parameter behaved slightly stiffer than the reduced integration element S4R.

SHELL-SOLID-SHELL APPROACH
The use of standard three-dimensional "brick-solid" elements in the numerical simulations of slender structural components in bending dominated applications involve expensive computational procedures.This is a consequence of the high mesh density required, including several elements across the thickness direction and also throughout in-plane direction to preserve an adequate length side vs. thickness aspect ratio.An adequate modeling option including the complete relevant stress components can be carried out by means of combining standard shell and solid elements.Specifically, both typologies of elements are tied together with one or two solid elements across the thickness (see Figure 8), the approach being denoted as Shell-solid-shell.This alternative allows the complete three-dimensional stress field at the solid element to be obtained, while at the same time, the model behaves satisfactorily in bending applications due to the presence of the shell elements.
The location of the brick elements within the finite element set is generally determined by the presence of a specific mechanical feature at which a more accurate representation of the stress field is needed.This is the case of the adhesive layer in co-bonded and cocured skin-stringer joints in composite stiffened panels, where to include the corresponding offset parameter of the shell elements to avoid material overlapping, see Figure 8.b, is also required.The coupling between the shell and the solid elements in this approach is carried out through several numerical procedures, as for instance the * TIE constraint in ABAQUS [18].It is worth mentioning that this modeling option has been applied satisfactorily in some engineering applications concerning composite stiffened panels [24].In addition, note that the numerical costs involved can be further decreased by using reduced integration scheme elements for the displacement interpolation of both element types.

NUMERICAL ANALYSIS. FE MODELS
The numerical analysis of the structure was conducted using ABAQUS.Due to geometrical characteristics of the specimen described in Section 2, a conventional shell model was assumed to represent the most suitable option for this purpose.Particularly, first-order elements were selected for the discretization of the structure.Within the linear elements, for finite strain applications, it is possible to choose between fully integrated elements S4 (even including the Incompatible Modes Method, which can be seen as a particularization of the EAS Method) and reduced integrated elements S4R with hourglass stabilization [18].
In the case of composite structures, both elements are able to model laminated components including the possibility of performing the numerical integration of the constitutive law over the thickness by means of either the Simpson or the Gauss integration rules.The main difference between these typologies of elements relies on the computational costs involved in the numerical simulations, where the reduced integrated scheme elements allows the running times providing accurate results to be significantly reduced.
Accordingly, a monolithic FE model, composed of a single part, was developed from the nominal dimensions appearing in Section 2. The resulting mesh followed the same topological considerations exposed in [11], having 12200 linear quadrilateral elements of type S4R and 12221 degrees of freedom.However, in order to obtain reliable and efficient predictions of the stress state at the most critical locations of the structure during the loading process, two additional models based on the Shell-solid-shell approach and on the Three-dimensional 7-parameter shell element were performed.
First, using the Shell-solid-shell approach (called as multi-part model in the following), the skin and the stringers of the panel were modeled as independent entities using the S4R elements for their discretizations.The adhesive layer was modeled by a thin layer of C3D8R elements (standard first-order reduced integrated brick elements) with one element through the thickness, connecting both stringer and skin surfaces, see Figure 8.The resulting finite element mesh using this approach, with the same element size than the monolithic model, now has 15400 C3D8R elements for the adhesive layer, 3600 S4R for each of the stringers and 6600 S4R for the skin of the panel (see [11] for more details).
Second, a monolithic-like finite element model was also performed.In this modeling option, the Three-dimensional 7-parameter shell elements were included at the skin-stringer joints, substituting in this way the original S4R elements.Hence, the entire panel was modeled as an unique entity, only differing from the previous monolithic shell model in the shell element typology employed at each location of the specimen.The mesh density used was coincident with the monolithic model, therefore resulting in a combination of 10600 S4R elements and 1600 Three-dimensional 7-parameter elements, and with the same number of nodes.Note that as the seventh parameter in the Three-dimensional element was inserted via the EAS Method, the complete model has the same number of degrees of freedom than the monolithic model.
The load was simulated as a uniform pressure acting on all the elements of the skin.The clamped line generated by the fastener system would be located somewhere across the overlap region between the panel and the pressure box, which is about 5 mm in width, see [11] for more details concerning the numerical modeling of the clamping system.
Additionally, it is worth mentioning that with reference to the model which includes the Three-dimensional 7-parameter elements, some shortcomings are encountered in the application of the uniform pressure over these elements.As was already mentioned in Section 3, these elements were inserted into the commercial software ABAQUS by means of a userdefined subroutine.Specifically, all the types of external mechanical loads can be directly applied to these user-defined elements with the exception of the pressure loads on the element surface.This is due to the fact that the numerical solver does not have information with regard to the shape function employed for the displacement interpolation within the element, and therefore it is not able of transforming this pressure into equivalent nodal forces [18].
In view of this limitation and for visualization purposes (that is not supported for user-defined elements), a set of dummy elements overlap the elements implemented, which have in comparison with the Three-dimensional elements: (i) the same nodal locations and elements connectivity, (ii) the same order of displacement interpolations, and (iii) negligible material properties in order to not disturb the numerical results of the model.Thus, the external pressure was applied uniformly on the skin surface elements (S4R) at the central and at the extreme portions of the skin, and on the dummy elements at the skin-stringer joints.

BUCKLING ANALYSIS
Prior to conducting the postbuckling simulations of the panel, a buckling analysis of the structure was performed.This preliminary study serves as a validation procedure of the more sophisticated numerical models, as is the case of the multi-part one [11] (Shell-solidshell approach), and for other two fundamental reasons: 1. From the experimental point of view, this analysis was required to determine: (i) the most suitable locations on the specimen of the strain gages and the displacement transducer, and (ii) the final loads of each of the tests performed, in order to characterize the postbuckling evolution of the structure, see [23].
2. From the computational point of view, the buckling analysis provides the information concerning the displacement field associated to each buckling load, in order to perturb the nominal geometry of the panel according to the scheme of equation ( 14) given below.Therefore, the postbuckling analysis of the structure can be conducted based on the actual geometry of the panel, which would contribute to obtain more accurate numerical predictions, especially with regard to the final deformed shape of the panel.
With regard to the model that includes the Three-dimensional 7-parameter elements, no buckling analysis is presented in this investigation.This is due to the fact that the buckling loads and mode shapes obtained are associated to local instabilities of the dummy elements.This is a consequence of the very low mechanical properties assigned to these elements in comparison with the actual properties of the materials of the panel.However, due to the strong similarities with respect to the monolithic model, this model does not require a specific validation process, as was the case with the multi-part model.Thus, this monolithic-like model with 7-parameter shell elements can be directly used in the postbuckling analysis of the panel, see [25].

POSTBUCKLING ANALYSIS
The presence of initial geometric imperfections can affect notably the mechanical evolution of slender structures [1], [6], [19].Consequently, an experimental procedure that allow their determination and their subsequent incorporation into numerical simulations was developed and satisfactory applied in [11].This procedure approximates the actual geometry of the panel by a certain combination of the deformed shapes obtained in buckling analysis through a computational scheme based on the Least Square Method.
In particular, this approximation is performed by modifying the ideal position of the nodes by the addition of a combination of eigenvectors obtained from the buckling analysis of the problem.In this way, if x o i (n), i = 1, 2, 3, are the ideal coordinates of the node n, and u k i (n) are the components of the displacement of the node n corresponding to the k th buckling mode, a particular perturbed position of node n, x i (n), is obtained combining M buckling modes (usually the lowest M modes) in the following form: a k constituting any arbitrary coefficients set.Thus, in this work, from 2 to 30 buckling mode shapes were considered for the application of the former procedure, different perturbed geometries having been obtained.Hence, Table 1 shows the coefficients a i of equation ( 14) when 2 (label P02), 5 (label P05) and 10 (label P10) modes are used.Notice that these scaling buckling factors imply that initial geometric imperfections were about 50% of the skin thickness.The postbuckling analysis of the specimen was conducted including the geometrical nonlinearities existing in the commercial FE software ABAQUS/Standard and using a Newton-Raphson algorithm to solve the problem, and considering the perturbed geometry of the panel by scaling the buckling modes with the coefficients previously determined.

Postbuckling analysis by means of the Shell-solid-shell approach
Evolutions of T-1 and T-2 displacements for some combination of buckling modes are shown in Figure 9.Note the strong dependency of the postbuckling evolution on the specific combination of buckling modes.Thus, some cases (P02 and P05) reached the total load that appears in these plots (fixed in this case in 35000 Pa, around 3.5 times the lowest buckling load numerically predicted), whereas other combinations (P10, P15, P20, P25 and P30) did not.These difficulties in achieving equilibrium states along the resolution process can be minimized using the ABAQUS's * stabilize option, which introduces a pseudoviscous behaviour in the analysis which allows a certain amount of energy to be dissipated, see [11].By including this option, all cases reached 0.35 bars pressure and, comparing evolution in Figures 9 and 10, the solutions for most of the cases are only slightly altered by the * stabilize option.As was discussed in [11], by performing the comparison between the numerical predictions with experimental data, two main aspects should be commented.First, the discrepancies between the numerical and the experimental data with regards to the evolution of T-2 might be associated to the position of this point near the boundary, and therefore the actual clamping conditions could affect its evolution.Second, with reference to the evolution of T-1 with the current measurements of the skin of the panel, it was observed that the use of combinations P05 or P06 was the best option to simulate the mechanical response of the specimen.This is due to the fact that only four lines placed between the stringers were measured to approximate the actual geometry of the panel.Hence, if higher modes are considered, as the maximum geometrical imperfection abandons the central line, affecting the right and left portions of the skin of the panel (see [11]), the process does not have enough data to adjust properly initial geometry along axial direction.if better agreements are needed, more lines will have to be measured along the axial direction, then allowing the incorporation to the procedure of higher buckling modes.In analyzing the deformed shape of the panel predicted through the present modeling (considering the first 5 buckling modes, P05), Figure 11 shows the distributions of the radial displacements of the panel with respect to the frame represented in Figure 1.
To conclude the analysis using this approach, Figure 12 shows the stress components (peel and tangential stresses) distribution along the adhesive layers for a pressure of 3.5 bars.The locations of the most severe conditions, which are pointed out in this plot, coincide with those appearing in the experimental analysis reported in [25] (see also [11]).This section is devoted to the postbuckling analysis of the panel under study using the Three-dimensional 7-parameter element and considering the perturbation of the geometry of the panel based on the experimental measurements aforementioned in Section 2.
For the present application, this 7-parameter shell element allows an estimation of the peel and the transverse shear stress components to be made without introducing any additional numerical feature, as was performed in the case of the TIE constraint in the Shell-solid-shell approach (corresponding to the multi-part model).
In addition, due to the fact that special interest is focused on the stress distribution at the adhesive, only the stress and the strain components at the top and at the bottom of this layer are defined as State Dependent Variables (SDVs) with respect to the cylindrical frame shown in Figure 1 [25].In this way, the number of SDVs is reduced to 24 for each element.It is worth mentioning at this point that although the SDVs are computed at the 7-parameter Three-dimensional elements, they are transferred to the dummy elements for visualization purposes.
Regarding the treatment of the initial geometric imperfections by means of this approach, a specific explanation is required.There exists a strong similarity between the numerical buckling mode shapes corresponding to the monolithic and multi-part models, where the structural instabilities affect more significantly the central region of the skin located between the stringers [11].Since the Three-dimensional 7-parameter elements were only placed at the skin-stringer joints (a region not very influenced by the buckling modes shapes), it can be as-sumed that the coefficients a i of equation (14) for the present modeling option are equivalent to those presented in Table 1 (see also [25]).Therefore, these combinations of the buckling modes were employed in this case in order to obtain the perturbed geometries of the panel.
Figure 13 depicts the correlation between the experimental data corresponding to the LVDT at the location T-1 and the numerical results including the Three-dimensional 7parameter shell element.First of all, it is observed once again the strong dependency of the final deformed shape of the panel with respect to the initial geometric imperfections (depending on the number of buckling modes selected to perturb the nominal geometry of the panel).In this case, all combinations exhibit a satisfactory agreement with respect to the experimental data up to around the first buckling load predicted, see Figure 13.b [23].However, for loads higher than this, similar to the multi-part model, the best adjustments are obtained when the first 5 (P05) and 6 (P06) buckling modes are considered until a load magnitude of 0.2 bars, i.e., around 2P cr , which is similar to the case of the multi-part model by means of the Shell-solid-shell approach.In order to overcome severe converging difficulties in achieving equilibrium solutions, the complete set of computations included the ABAQUS's * stabilize option, using the default magnitude parameters.The simulations concluded at around 2P cr , as can be observed in Figure 13.It should be also commented that the artificial viscous energy inserted by this parameter into the simulations did not exceed the 5% of the total strain energy of the model.Note that in this specific application, these convergence difficulties are mainly associated with: (i) the 7-parameter shell element selected at the skin/stringer joints (see Section 3), and (ii) the structural instabilities as a consequence of the external loads.
Since no structural instabilities were encountered at this load level, the numerical difficulties in achieving equilibrium solutions are mainly attributed to the difference between the mechanical interpretation of the nodal degrees of freedom 4, 5 and 6, which are related to the difference vector w in this 7-parameter shell element, whereas in standard shell elements, these degrees of freedom usually correspond to the nodal rotations [25].
In analyzing the deformed shape of the panel predicted through the present modeling option (considering the first 5 buckling modes, P05), Figure 14 shows the distributions of the radial displacements of the panel with respect to the frame depicted in Figure 1.From this plot, it is observed that the final deformed shape of the panel is characterized by the appearance of two main half-waves along the circumferential direction, where one of them is subdivided into two sub-waves.Here, this subdivision is clearly noticeable for a load magnitude of 0.21 bars, at which the computation was terminated due to convergence issues.This fact agrees with the experimental observations [23], although slightly differing from the multi-part model (Figure 11), where the subdivision of one of the main half-waves was observed in this 7-parameter model at lower load levels, particularly at around 0.15 bars, see Figure 11. Figure 14.Evolution of radial displacement for P05 combination in the postbuckling analysis using the Three-dimensional 7-parameter shell element Finally, one of the main contributions concerning the use of Three-dimensional 7parameter elements is that these elements allow us to estimate the entire stress and strain distributions at specific locations in the specimen without any additional numerical "feature".Figure 15 depicts the distributions of the peel and the transverse shear stress components for a pressure of 2.1 bars, which correspond to the State Dependent Variables (SDVs) UVARM1, UVARM2 and UVARM3, respectively.The locations at which the most severe conditions achieved are highlighted in this diagram as in Figure 12.Particularly, in this plot it is observed again that the most critical areas are located principally at the circumferential borders of the specimen, being in this way coincident with those areas identified in the experimental study (see [23]) as well as in the multi-part approach (Section 7.1).

CONCLUDING REMARKS
This contribution analyzes the applicability of different numerical strategies based on FEM to predict the response of a composite stiffened panel subjected to uniform pressure.
Computational models were developed by means of several structural approaches: (i) a monolithic shell model, (ii) Shell-solid-shell model, and (iii) a monolithic-like model in-Figure 15.Distribution of stress components at the adhesive layers for an applied pressure of 2.1 bars when 5 modes were used to adjust the initial geometry of the panel using the Three-dimensional 7-parameter shell element at the skin/stringer joints cluding Three-dimensional shell elements at the skin/stringer joints.Particularly, in the Shellsolid-shell option, solid elements were used to model adhesive layers, and shell elements for the skin and for the stringers.The two latter modeling options lead to the prediction of peel and shear stresses fields around the adhesive layer along the loading process.
All models performed were able to incorporate the initial geometrical imperfections determined in the experimental part of this research using a computational procedure based on the Least Square Method.Simulations including the initial geometric imperfections by means of this procedure showed satisfactory agreements with the experimental data in the postbuckling range using the first five or six buckling modes.
In this regard, the Shell-solid-shell approach exhibited a satisfactory accuracy up to 3.5 times the first critical pressure numerically determined.
For the monolithic-like model including the Three-dimensional 7-parameter shell elements at the skin/stringer joints, in contrast to the other two approaches, strong difficulties in achieving equilibrium solutions along the simulation were encountered.The computations were conducted using the ABAQUS' * stabilize option.The maximum load reached was around 0.2 bars for all the cases considered.
The stress and the strain distributions at the skin/stringer joints were estimated through two different approaches: Shell-solid-shell technique and the monolithic-like model using 7parameter shell elements.Both models predict the most critical areas at these joints near the circumferential extremes of the panel.Consequently these approaches correlate satisfactorily with the experimental observations.

Figure 1 .Figure 2 .
Figure 1.Geometry and material data of the specimen

Figure 3 .Figure 4 .
Figure 3. Top view of the experimental set: pressure box and panel

Figure 5 .
Figure 5. Kinematics of the 7-parameter shell model

Figure 7 .
Figure 7. Tip cantilever composite beam subjected to end forces.Load deflection curves obtained with the implemented TShell7p elements in comparison with the ABAQUS' S4R elements

Figure 8 .
Figure 8.(a) Shell-solid-shell modeling approach, (b) Application of the Shell-solid-shell option to a skin-stringer joint in a composite stiffened panel

Figure 9 .
Figure 9. Evolutions of T1 and T2 displacement for some adjusted configurations

Figure 10 .
Figure 10.Evolution of T1 displacement for some adjusted configurations with ABAQUS's * stabilize option

Figure 11 . 1  2 Figure 12 .
Figure 11.Evolution of radial displacement for P05 combination in the postbuckling analysis

Figure 13 .
Figure13.Evolution of T1 displacement for some adjusted configurations including the Three-dimensional 7-parameter shell elements at the skin/stringer joints