On FEM analysis of Cosserat-type stiffened shells: static and stability linear analysis

The present research investigates the theory and numerical analysis of shells stiffened with beams in the framework based on the geometrically exact theories of shells and beams. Shell’s and beam’s kinematics are described by the Cosserat surface and the Cosserat rod, respectively, which are consistent including deformation and strain measures. A FEM approximation of the virtual work principle leads to the conforming shell and beam FE with 6 DoFs (including the drilling rotation for shells) in each node. Examples of static and stability linear analyses are included. Novel design formulas for the stability of stiffened shells are included.


Introduction
Shells are commonly used as a part of various structures. They are considered as an economical, efficient, and aesthetic choice since are able to bear high loads with low weight kept. Additional reinforcement could be achieved by using locally placed stiffeners. Stiffened shells have various applications, like structural elements in vessels, aircraft, and aerospace vehicles. They are also used as bridge decks and as floor slabs in buildings. Inspiration for this type of structures might be found in nature: leaves and insects' wings are often membranes enforced with veins.
Over the years, great effort has been done to investigate and understand stiffened plates and shell behaviour in various circumstances: under static and dynamic loading, deformation in thermal environment, free vibration, etc. Nevertheless, still novel contributions are provided to better understand their behaviour and capabilities in up-to-date applications.
The current study is focused on the employment geometrically exact shell and beam theories to stiffened shells analysis as the original contribution. This approach results in compatible kinematics' formulations of a shell as the Cosserat surface and beam as the Cosserat rod. The virtual work principle, applicable to the stiffened shell as an unseparated structure, is used to develop the finite element method (FEM) approximate solution. Consistency of shell and beam finite elements (FE) is achieved. In the present paper, the static linear analysis and stability analysis of various stiffened homogenous shells are conducted, including curved and branched shells. Since stiffened shells are widely used in engineering design, current parametric studies consist not only of rare data (e.g. displacements and critical load factors), but also a proposition of approximate design formulas given.
The present paper is prepared on the occasion of the 65th birthday of Prof. Altenbach, in recognition of his research in the field of the shell structures. It should be noted that the models of complex models of beams, plates, and shells are in the focus of interests of Prof. Altenbach, see e.g. papers [1][2][3][4][5] and books [6,7].

Literature overview
Many papers, devoted to the problem of static response and buckling of plates and shells stiffened with beams, have been published over the years. The present overview will be focused on those papers, in which reference solutions (mainly which employs FEM) for isotropic stiffened plates and shells could be found.
In 1940s, three reports have been published, which provide reference results for the buckling of unstiffened curved panels [8] and stiffened curved panels [9,10]. The results are obtained by the analytical approach, and even nowadays many papers refer to those as reference solutions. Buckling of stiffened panels and shells was also investigated by Timoshenko, by analytical methods. His results could be found e.g. in [11].
Pioneer FEM analysis of stiffened shells has been carried out in McBeam's thesis in 1968 [12]. Since 1980s, FEM is widely employed for the analysis of stiffened plates and shells. Constant development of FEM theory and codes, together with the rapid growth of computer resources, indicates the necessity of re-considering even some basic problems, keeping in mind that FEM is an approximate method and needs careful verification.
In 1970s and 1980s, several results based on numerical calculations were published. Rossow and Ibrahimkhail [13] used FEM enriched with a constraint method for the linear study of stiffened plates. Sobel and Agrwal [14] adopted the numerical code based on Donell's theory to obtain buckling loads of curved panels. Mizusawa and co-workers [15] used the Rayleigh-Ritz method with B-spline functions as coordinates to investigate buckling of skew plates in a wide range parametric analysis. Bathe and Bolourchi [16] presented a novel finite element for different aims, including linear and nonlinear static analyses. Those were shell elements with max. 16 nodes and beam elements with max. 4 nodes, both with the Lagrangian interpolation functions. The integral equation technique was used by Srinivasan and Thiruvenkatachari [17] to analyse the static and dynamic responses of an annular eccentrically stiffened sector. Deb and Booton [18] proposed an 8-noded Serendipity plate FE and a 3-noded stiffener element dedicated to the linear analysis.
In the last decade of the twentieth century, further development of computational methods was observed. Bhimaraddi and co-workers [19] developed the special-purpose FE for annular plates. Palani and coworkers [20] proposed an isoparametric FE with quadratic shape function applied for the analysis of stiffened plates and shells. Nevertheless, numerical examples consist only of plate tasks. Kolli and Chandrashekharat [21] provided a formulation of a FE for stiffeners which incorporates the effect of lateral strain, with verification calculation conducted for isotropic shells. Bedair [22] original contribution was the usage of sequential quadratic programming in stiffened plates analysis. Work by Satish Kumar and Mukhopadhyay [23] is devoted to composite structures, but it should be mentioned as the first FEM implementation (to author's best knowledge), where the placement of stiffener is independent from plate nodal lines. Sadek and Tawfik [24] followed this concept, developing a 9-noded plate FE with stiffener placement parallel to nodal lines.
Here the review reaches the twenty-first century, with the paper by Wen and co-workers [25] who used the boundary element method to analyse stiffened plates. Interaction forces between a plate and stiffeners are considered. Peng and co-workers [26] employed the element-free Galerkin method to perform the static analysis of stiffened plates. This approach allows independency of discretization of a plate and a stiffener. The differential quadrature element method was used by Jiang and co-workers [27] to investigate the buckling of stiffened cylindrical panels under axial compression. Voros [28] applied the novel stiffener FE with 7 DOFs in each node for the plate buckling analysis. Torsion warping and flexion-torsion interaction are included in stiffener FE.
Ojeda's thesis [29] consists of geometrically linear and nonlinear analyses of stiffened plates with FEM employed. The Bedair's paper [30] is focused on the design of stiffened box girders (multibranched shells) in buckling. A numerical procedure is presented there, which uses a combination of energy formulations and mathematical programming techniques. Jafarpour Hamedani and Rahbar Ranji [31] adopted conventional and super-elements to a FEM analysis of buckling of stiffened plates. Again, stiffener placement is arbitrary. A parametric study of stiffeners' stiffness, various boundary conditions, and loads is carried out. Nguyen-Thoi and co-workers [32] proposed the usage of triangular elements with the cell-based smoothed discrete shear gap method to the static linear, frequency, and buckling analyses of stiffened plates. Tran and co-workers [33] performed a FEM analysis of buckling and ultimate strength of curved panels in uniaxial compression. A design procedure for such structures is proposed. Shi and co-workers [34] analysed a static response and buckling loads of plates stiffened with arbitrary placed curved stiffeners. Bedair [35] committed the paper that was devoted to the problem of the design of stiffened box girders, in case of web shear buckling. Panda and Barik [36] developed a 4-noded stiffened plate FE, based on an isoparametric approach with the shearlocking problem solved. The approach was verified by a buckling analysis of plates (square, skew, rectangular, circular). Hosseini and Soltani [37] adopted the meshless collocation method to the linear static bending analysis of stiffened plates. Zhang and Xu [38] provided results of static analysis of stiffened plates with partially composite action. Investigations are based on an analytical approach.
Summarizing, many investigations were dedicated to stiffened plates and shells. The special-purpose and universal FEM applications were developed. Other methods were also used, proving their accuracy and efficiency. Many reference solutions are available, which is useful in verifying novel theories and codes. Nevertheless, curved and branched shell examples are insufficiently represented, even though such shells are widely used in the design of real structures.

Theory
The Cosserat approach to modelling shells and beams is widely used in the literature, see reviews for shells [43,44] and the references which are given in [45][46][47][48]. Below the derivation of the model using a 3D-to-2D reduction is considered.

3D-to-2D reduction
What distinguishes the 6-parameter nonlinear shell theory (see e.g. fundamental works by Zhilin [49] and Chróścielewski et. al. [50]) and the exact nonlinear rod theory [51] from other well-established theories is the exact reduction of 3D equilibrium equations to 2D equilibrium equations, without any kinematic assumptions (Fig. 1). It results in exact virtual work principles with the material law excluded from them. Integration of the 3D balance law could be done by an application of a standard through the thickness integration, for further details see [50,51]. Both theories are capable of describing structures that experience finite translations, finite rotations, and small strains.
Integrations performed on the shell and rod bodies separately are described in previously mentioned papers. The present paper employs results obtained there, and the following derivations start from the reduced 2D shell and 1D rod. Nevertheless, some comments on the integration over the body of stiffened shell should be enlisted: • a division between shell and rod bodies is arbitrary, but must be chosen to fulfil some requirements on their regularity; • in some cases, the problem of doubly integrated volume may arise. A similar problem is observed in junctions of multibranched shells (see e.g. [52,53]); • shell and rod are perfectly bonded in initial and any following configuration, without possible separation or slipping; • no limitations on the rod's cross section are introduced. A centre of mass and shear centre do not have to lay directly on a curve that defines the reduced beam. Formal derivations will start from the configuration as in Fig. 2. There is an arbitrarily chosen part Π of the surface M, with the edge divided into 3 parts: ∂Π (arbitrary cut), ∂ M f (free edge), and Γ (free edge where stiffener is placed). The stiffing rod is drawn outwards shell in Fig. 2 for clarity. The current study omits multibranched shells, where Γ could be also a singular curve on which branches meet. This extension could be straightforwardly adapted from previous studies. Every point of the configuration is defined by the position vector y.
Various loads are acting on the body: surface force and couple loads (f and c); stress resultant and stress couple vectors (n n and m n ) on the edge ∂Π; external edge loads (n * and m * ). There are also loads acting on the rod: distributed resultant force and couple (f Γ and c Γ ); stress resultant and couple (n + , n − , m + and m − ). Stress resultant and couples are equal to external end loads (n + = n * Γ , m + = m * Γ ) at the free end of the rod. The novelty of the current approach is in taking into account the interaction between shell and rod, which is represented by distributed resultant forces f Γ and couples c Γ . They act on both sides of the imaginary cross section, with opposite signs, and are treated as resultant of normal and traction stresses acting between the 3D shell and the 3D rod before the dimension reduction.
The total force acting on the shell is given as and the total moment is The total force acting on the rod is given as and the total moment is Thus, the next step of the derivation is to transform those formulas to surface integrals over Π and curvilinear integrals over Γ with appropriate boundary conditions. This is performed with the usage of surface Cauchy theorem and surface divergence operator Div. As the result of Eqs. (1-4) conversion (for details see e.g. [52]) the following formulas are obtained: where ad: E 3 → so (3) is an isomorphism, where so (3) is a vector space of skew-symmetric tensors (see [50] for details). The symbol (.) denotes differentiation with respect to an arc coordinate of the rod's reference curve.
The global equilibrium is achieved when Meanwhile, the local equilibrium conditions, which result from Eqs. (5)- (8), are as follows: -for a regular point of the shell: Div where F = ∇ y. Those are supplemented with boundary conditions for the shell: at the edge ∂ M f : The next step towards FEM implementation is to build a virtual work principle based on local equilibrium conditions and appropriate boundary conditions. Details of derivations could be found, e.g. in [54] (for shells) and in [55] (for rods). It is important to mention that further derivations will be made without any simplifications nor additional assumptions. So, the following formula could be written: where v and w denote any vector field over the area Π . After transformations the following formula is obtained: where W = ad w, ∇v, and ∇w are surface gradients of vector fields. Similar transformations could be done with local equilibrium conditions for rod, obtaining and after transformations Equations (10) and (12) are added to obtain the virtual work principle for a stiffened shell. After this summation is made, the term Γ (n Γ · v + m Γ · w) dl, which includes contact forces and couples, vanishes. Since this is possible at this point, the rest of the derivations which lead to the finite element formulation (e.g. interpolation of quantities) are carried out in a known manner and will be omitted here.
Equations (10) and (12) consist of terms which can be interpreted as the virtual strain and curvature tensors for the shell and for the rod

Kinematics of stiffened shell
It is assumed that the position of each point of undeformed shell reference configuration is given by vector x (Fig. 3). At each point tangent space T x M exists, defined by two vectors t 0 β = x, β . There is also an orthogonal vector The structure tensor T 0 ∈ SO(3) is defined by the orthogonal rigid triad t 0 i (x). Stiffeners are placed along some curves, chosen so that one of the vectors t 0 1 or t 0 2 remains tangent to them. Stiffener and shell share the same position vector and microstructure tensor and remain tied during motion. The current position and orientation are given by where Strains and curvatures at the reference surfaces are given by: If the index of the axis parallel to the beam is chosen (β = 1 or β = 2), strains in the beam could be easily recovered from (17). Virtual strains and curvatures given by (13) now can be denoted as: Again, it is enough to choose the proper index β to get virtual strains of the beam.

Finite elements
In the FEM approximation, shell and beam-compatible Lagrangian FE are developed and used to conduct the linear static and stability analyses. Appropriate linear material laws for shell and beam are applied in the FEM analysis. Numerical examples are limited to the cases where stiffeners are placed exactly at the shell FE edges and they share nodes, so that problem of nonconforming shell and beam discretizations is not addressed in the current study.
In numerical studies in the present paper, conforming finite elements used for shell and beam domain are used (Fig. 4). The same number of degrees of freedom (3 translations and 3 rotation parameters) is present in each shared node, and conforming definition of shape functions (a 3rd-order polynomial) is chosen for shell and beam FE. Both types of elements are derived by the appropriate application of the theory described above. What's important, shell elements are equipped with so-called drilling rotation, which acts with respect to shell normal. This rotation corresponds with beam rotation related to in-plane bending.

Shell and beam 6-parameter FE
In the current study, fully integrated CAMe16 shell finite elements are used to simulate shell structure. They are 16-node C 0 Lagrangian elements with 6 parameters (3 translations u, v, w, and 3 rotation parameters ψ 1 , ψ 2 , ψ 3 ) in each node. They were widely used in previous studies of unstiffened shells, e.g. [56][57][58]. They are characterized by negligible hourglassing and locking effects without additional techniques implemented. The constitutive relation for elastic, homogenous, isotropic Cosserat shells is discussed, e.g. in [57]. In the current study, parameters E (Young's modulus) and v (Poisson's ratio) are various in numerical examples, while other parameters are kept the same, namely N = √ 2/2 and α t = 0.01. In beam's discretization, B6e4 finite elements are used. These are 4-node C 0 Lagrangian elements. Their application is shown, e.g. in [55], and similar FE was also developed by Smoleński in [51]. The constitutive relation for the beam is given as where σ = N 1 T 2 T 3 M s M 2 M 3 T and ε = ε 11 ε 12 ε 13 κ 11 κ 2 1 κ 31 T denotes forces and strains, respectively, while C denotes a constitutive matrix. Since in applications eccentric beams are often used (see Fig. 5), relation (19) could be rearranged to form where In the present paper, the constitutive matrix for a beam with respect to its centre of mass is given in standard form asĈ where inertia moments are calculated with respect to axest 0 2 ,t 0 3 (Fig. 5).

Solution techniques
A linearized equation of equilibrium (summed Eqs. (10) and (12)) after discretization for a typical finite element has a form: In Eq. (23) symbols not mentioned before stand for: K (e) . Matrices and vectors are calculated for each finite element and are aggregated to obtain global matrices and vectors, denoted without (e) index.
To compute results for the linear analysis, a set of linear equations is solved to obtain translations and rotation parameters for the prescribed load. In (24) K M denotes the tangent material stiffness matrix obtained in the reference (undeformed) configuration. The theory itself was not linearized; thus, it is not the static linear analysis in the rigorous form.
In the stability analysis an eigenproblem in the form is solved, giving eigenvalues λ (critical load multipliers) and eigenvectors v (buckling modes). Matrices K M and K G are calculated at a basis of the deformed configuration, resulting from the appropriate linear analysis.

Numerical examples
All calculations are carried out in the author's code written in Fortran. PARDISO sparse linear solver is employed in the linear analysis, and FEAST sparse eigenproblem solver is used in the stability analysis. Commonly, Young's modulus is denoted as E, Poisson's ratio as v, and a shell thickness as h.

Slab supported by a single column
The first example is not exactly a stiffened shell, but a slab with a tied perpendicular column. The task was previously employed in research [59][60][61][62]. This analysis's goal is to demonstrate that the connection between finite elements of two types is properly modelled. Geometry is shown in Fig. 6. The slab's thickness is constant and equal to h = 1, and the column's cross section is circular with radius r = 0.25. The whole structure is made of a linear elastic material with parameters E = 10 6 and v = 0.3. Two load cases are considered: A with concentrated loads and B with concentrated moments, both with equal resultant moment on the structure. The discretization of the slab is regular and consists of 4 × 4, 8 × 8, or 16 × 16 CAMe16 elements, while the column is always divided into 4 B6e4 elements. Expected deformation in both cases is achieved, namely uniform twisting of the column and almost rigid rotation of the plate. The first results compared with reference are the drilling rotation value at the slab's vertex, and they are collected in Table 1. Values of translation of vertex in the direction of the acting force in case A (or corresponding values in case B) are enlisted in Table 2. Obtained results are identical in both cases. Mesh density variation does not cause a meaningful change in the absolute value of displacement, but, what is worth noticing, the absolute difference between meshes 4 × 4 and 8 × 8 is smaller than between 8 × 8 and 16 × 16. It could be explained by analysing drilling rotations values in the close neighbourhood of the connection point of column and slab (Fig. 7). The value of the rotation of the middle node is dependent only on the column's properties and its FE discretization (the same in every analysis). A significant jump of rotation value between the middle and next node is observed; it rises along with mesh refinement. Nonmonotonous distribution of rotation angle values is observed approximately in the range of single FE and nearly constant value of the angle in the far field.

Stiffened simply supported plate
In this example, simply supported plates with various dimensions ratios (square and rectangular), stiffeners placement, and parameters are analysed (Fig. 8). The common attribute of these analyses is taking advantage  . 7 The supported slab, drilling rotation values Fig. 8 The simply supported plate geometry, loads, stiffener's cross section of symmetry, so that torsional and lateral bending stiffnesses of beams, lying on symmetry planes, can be neglected. This example tests the capability of the current theory and its implementation to predict the linear behaviour in bending under distributed or concentrated transverse loads. Material and geometrical data are listed in Table 3. In the square plate, only the x-stiffener is used and this plate is analysed in two cases: with a concentric or an eccentric stiffener. In the rectangular plate, the x-stiffener and y-stiffener are different from each other, but common for two analyses (with concentrated or distributed load). The square plate with a single stiffener placed along the x-axis was analysed in papers [13,[20][21][22][24][25][26]32,37]. The second case, the rectangular plate with two orthogonal stiffeners, was the topic of research in [ [24] 4.632* 1.424* Wen et al. [25] 4.553* 1.335* Peng et al. [26] 4.87** 1.36** Nguyen-Thoi et al. [32] 4.498* 1.431* Hosseini and Soltani [37] 4.8363* 1.3647* *Explicitly given in reference **Read from a graph in reference  [20] 14.95** -Kolli et al. [21] 8.702* 1.240* Bedair [22] 8.00* 1.12* Peng et al. [26] 8.70** 1.30** Tamijani et al. [63] 9.75** -Hosseini and Soltani [37] 8.875** -Zhang and Xu [38] 8.91265* 1.23191* Cunha et al. [64] 11.2** -*Explicitly given in a reference **Read from a graph in a reference The results of the displacement of the centre point are collected in Tables 4 and 5. It is worth noticing that some of the reference results are obtained from graphs by their digitalization. Present results are given for 3 different discretizations. The obtained results are comparable with those from literature.
Additionally, graphs of displacements along selected paths are collected in Figs. 9 and 10. When the path is placed along stiffener, every used discretization (including discretization with the single shell and single beam elements) gives almost the same result. When paths located on the unstiffened area are analysed, the simplest discretization is revealed as not sufficient. This is more visible in cases with distributed loads. Summarizing,  when exactly these or similar tasks are analysed, not only middle point deflection should be compared with reference values, but also deflections along paths placed on unstiffened areas, to check both stiffened and unstiffened shell bending action.

Buckling of a plate with a single stiffener
In the third example, a wide range of parametric analysis of plate buckling is carried out. Rectangular, skew, unstiffened, and stiffened plates with various stiffener parameters are taken into account. The obtained results are compared to those taken from [11] (analytical solutions), [15] (Rayleigh-Ritz method with B-spline functions), and [65] (FEM method, based on Mindlin theory). For unstiffened skew plates, reference values of critical load are taken from [66] (analytical results for thin shells, span to thickness ratio equal to 1000). Although results will be presented as nondimensional values, it is worth to enlist explicit values of dimensions and loads, boundary conditions, and material parameters, to make further comparisons possible. Therefore, as shown in Fig. 11, dimensions are a = b = 1, t = 0.01, and linear elastic behaviour is described by Hooke's law with parameters E = 2e11, v = 0.3. The skew angle is denoted as ϕ (equal to zero for the rectangular plate).
Motion in the z-direction is prohibited on the whole edge, at two corners displacement in y-directions is equal to zero, and finally at one of the corners displacements along the x-axis are prohibited. Two edges of the plate are loaded in-plane with line loads heading in opposite directions. In order to induce a uniform   Additionally, stiffener bending stiffness around local axis 2 is denoted as E J 2 , and its torsional stiffness is given as G J s (G = E 2(1+v) is Kirchhoff modulus). Stiffness for bending in the second direction is neglected, just like eccentricity. It is assumed that the total cross-sectional area of the stiffened plate is a sum of bt and A s . For clarity, a constitutive matrix for beam could be given in explicit form as FE discretizations are denoted as M × N , where M and N are the numbers of shell elements along the x-and y-axis, while M denotes also the number of beam elements used. Tables 6,7,8,9,10,11, and 12 present current and reference values of buckle load multiplier λ crit = q crit /q ref . In the first test, the buckling load for the unstiffened plate with various skew angles was analysed ( Table 6). Present results reveal better correspondence with those obtained by the analytical method [66] than those from reference FEM solutions [31,36,63,65]. For present results, 1st, 2nd, and 3rd load multipliers are presented in a single table's cell, while reference results consist of 1st load multipliers.
Significant differences have been noticed for plates with skew angle 30 • and 45 • , both in unstiffened ( Table 6) and stiffened cases (Table 10). In other cases, good agreement between the present and reference solutions has been obtained. It is worth to note that enlisting 2nd and 3rd critical loads is an original contribution, deformation modes paired with those are usually more complicated than 1st one, and then it is more demanding for a used method to calculate them.

Buckling of an axially loaded cylindrical panel
In this example, the buckling of axially loaded of panel curved in one direction is analysed. The aim and scope of the simulation are similar to the previous example investigation of critical loads and corresponding modes in various parametric studies. Geometry, FEM discretization and boundary conditions are shown in Fig. 12. Two types of boundary conditions are used: in type A, all outer edges are supported in the z-direction, in type B in the radial direction. Additional boundary conditions are used to prevent rigid motions. Axial compression is applied as distributed loads placed along curved edges. The reference value of the load is q ref = π 2 E D/B 2 , and in latter considerations buckling load multiplier λ crit = q crit /q ref is given as a result. The first step of the study is the stability analysis of the unstiffened cylindrical shell. This type of problem was analysed, e.g. in studies [8][9][10][11]14,27,33,[67][68][69][70][71][72][73]. Attention must be paid to the definition of nondimensional parameters that represent the curvature and dimension of the panel. Two definitions are used among papers: Explicit shells dimensions used in the present study are a = 1000, a/b = 1.0, h = 10.0, E = 210,000, v = 0.3, while R is calculated from relation (27).
Results collected in Table 13 show good agreement between the present and reference results from paper [72]. The current investigation is extended to the first 3 values of critical load multiplier comparing to reference. Additionally, the influence of boundary conditions is studied for the most curved panel, showing minor differences. For further studies, shell discretization 8 × 8e16 is chosen. For stiffened shells, for both longitudinal and circumferential stiffeners, 8 4-node B6 elements are used. Among the years, various formulas were developed to predict the buckling of a cylindrical panel. Timoshenko and Gere proposed the following formula with the assumption of deformation as a product of trigonometric functions in both directions. Stowell in the report [8] cited equation given by Redshaw in 1934 and proposed new formula where λ ∞ is buckling coefficient for the flat shell (λ ∞ = 4 for rectangular, simply supported plate). Stowell assumed deformation given by multiplication of trigonometric (longitudinal direction) and exponential functions (chordwise direction). The same formula was obtained by Tran [73], by analysis of limit states (plate when Z → 0 and cylinder when Z → ∞). Furthermore, Domb and Leigh [71] proposed an empirical curve based on nonlinear FEM analysis of aircraft fuselage panels:  Finally, Martins and co-workers ( [67,72] proposed few formulas for buckling coefficient, e.g. for short panels under uniform compression. In [72] formula is given, where and In paper [67] simple formula was presented. The graph in Fig. 13 shows data generated by formulas (28)- (35), compared to present results obtained by FEM and linear buckling analysis. The present results are in good agreement with reference curves. From Z = 1 up to Z = 10 best convergence is with Martins' formula (10), then the transition zone begins, and for Z ≥ 50 Stowell's (30) and Domb's and Leigh's (31) formulas are in best agreement.
From this point, attention is paid to stiffened curved cylindrical panels. Figures 14 and 15 consist of results obtained for stiffened panels, with different stiffener types and panel's curvatures. Buckling modes are shown for curvature parameters in range Z ∈ 10 1.0 ; 10 2.0 . As in the unstiffened panel case, 1st mode shape changes gradually. For very stiff longitudinal stiffener this change is less pronounced, while in other cases it is clearly visible.
The obtained data points were used to determine empirical design curves obtained with data fitting. A least-square optimization algorithm was used. Those curves are based on the formula and values of coefficients are enlisted in Table 14 along with basic fitting statistics.

Buckling of a box section
In the last numerical example, the buckling of a box section is analysed. A literature survey and description of references where the selected design problems of box sections are addressed could be found in the paper     Figure 16 also describes boundary conditions which prevent rigid motions. Two separate load cases are considered, namely compression and torsion. In both cases point loads are defined by P ref = 10 4 and load multiplier λ. To avoid problems with stress concentrations under point loads on free upper and lower edges of the box is stiffened by quasi-rigid rims, made of beams with cross section 40x40.
The box is stiffened in two independent ways: longitudinal (two stiffeners on opposite sides of the box) or circumferential (four stiffeners which form a rim around the box). Both stiffeners are placed eccentrically, on the inner side of the box, with cross-sectional dimensions h s × b s (various in the following analysis). In the present analysis, the box is unstiffened or stiffened with only one type of stiffener. The whole structure is built of a linear elastic material with constants E = 1.07e5, v = 0.34.
Shell discretization consists of 2x4x8 CAMe16 finite elements. Stiffeners are placed exactly on node lines, discretized with an appropriate number of 4-node B6 finite elements.  The first analysed case is the compression test. Results (1st, 2nd, and 3rd critical load multipliers and modes) for the unstiffened box section, and two cases of stiffeners are collected in Tables 15 and 16. Longitudinal stiffeners have a significant influence on buckling load, with over 3 times increase of 1st buckling load with the addition of 5.3% of total structure volume. For circumferential stiffener, an 18% increase of buckling load is achieved for the biggest stiffener.
The second analysis is devoted to the case of torsion. Critical load multipliers and modes (1st, 2nd, and 3rd mode) are collected in Tables 17 and 18. As previously, longitudinal stiffeners have a bigger influence on critical load multipliers, up to 232% increase of first critical load.
The last study is a case of combined loading in which both compression and torsion loads are applied simultaneously. Both of them are proportional to the load multiplier λ. Stiffeners with dimensions h s = 8, b s = 2 are used. Results are presented in the interaction curves in Fig. 17, where loads are normalized with respect to buckling loads in pure compression and pure torsion.

Conclusions
In the present paper theoretical background and FEM application of stiffened shell are presented. The theory should be considered as exact theory, since it takes advantage of the exact reduction of 3D equilibrium equations to 2D equilibrium equations, without any kinematic assumptions. The virtual work principle is derived since it is necessary for FEM formulation.  The present paper also consists of a wide range of FEM analyses, including static linear response and buckling problems. Results are compared to those obtained from previous studies with excellent agreement achieved. Novel results are presented, especially including those from the analysis of branched shells, which were seldom investigated so far. Some design formulas for stability analysis of stiffened curved panels and box sections have been proposed at the basis of obtained results. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Funding This work was supported by the National Science Centre, Poland, with the Grant DEC-2012/05/D/ST8/02298.