SciELO - Scientific Electronic Library Online

vol.16 número6Simplified expressions for dynamic behavior of cylindrical shells uncoupled and coupled with liquidsAnalysis of the contribution of the block-soil contact in piled foundations índice de autoresíndice de assuntospesquisa de artigos
Home Pagelista alfabética de periódicos  

Serviços Personalizados




Links relacionados


Latin American Journal of Solids and Structures

versão impressa ISSN 1679-7817versão On-line ISSN 1679-7825

Lat. Am. j. solids struct. vol.16 no.6 Rio de Janeiro  2019  Epub 22-Jul-2019 


A methodology to investigate the influence of the fillet radius on the buckling of integrally stiffened panels

Fernando Gonçalves Garciaa  b

Roberto Ramos Jr.b  *

aEmbraer S.A. - Av. Brig. Faria Lima, 2170, São José dos Campos, SP, Brazil. E-mail:

bUniversidade de São Paulo - Escola Politécnica - Departamento de Engenharia Mecânica, Av. Prof. Mello Moraes, 2231, São Paulo, SP, Brazil. E-mail:


The use of integrally stiffened panels (ISP's) in wings of small and medium-sized aircraft is frequent in aviation. These thin-walled structures are often subjected to compressive loads during their life-time and, although the finite element method has become an important tool for engineers to obtain the buckling load of structural members, some results based in charts published by NACA are still used in the preliminary stages of the wing design in many aircraft companies. However, these charts used for calculation of the critical compressive stress consider only idealized stiffened panels and neglect geometric details as the fillet radius used in the current design of ISP's. The objectives of this paper are twofold: (i) to show that the charts published by NACA provide good results for the critical buckling stresses for several geometries of ISP's, when compared to finite element results with proper boundary conditions, providing fillet radii are also neglected in the finite element models and (ii) to show that the values of the critical buckling stresses for local instability of ISP's may be significantly increased when one considers the effect of the fillet radius, meaning that this parameter should also be considered in the optimal design of such structures. Several numerical results obtained with finite element simulations based on different geometrical parameters of ISP's are presented in this study.

Keywords: Integrally stiffened panels; buckling; fillet radius; finite element method; wing; aircraft structures


Aircraft wing structures are composed of spars, ribs and stiffened panels that, working together, are able to withstand the aerodynamic loads. However, during operation, the wing upper panel is mainly subjected to compressive loads that may lead to different kinds of structural instability (local or global). In addition to that, considering the progressive improvement of wing geometry, use of new material alloys and new manufacturing techniques − which have enabled the production of even more light-weight structures − buckling has become increasingly important in the structural design of these thin-walled reinforced panels.

In particular, use of integrally stiffened panels (ISP's) in wings of small and medium-sized aircraft is quite frequent in aviation. Differently from riveted joints, ISP's are machined from a unique raw material plate, so that a single machined part can save many manufacturing costs due to the reduction of, e.g., riveting assembly, stocking parts, tooling expenses and manufacturing time ‎[1]. Besides, ISP's can be machined in many different configurations, some of them taking advantage of this machining flexibility to create configurations with an improved buckling capacity like isogrid and orthogrid stringers, sub-stiffeners, integrally T stiffeners, curvilinear stiffeners, skin strips and buckling containment features [2-5]. It is also important to mention the new manufacturing techniques that are being developed, like, in near-term, high-speed machining from near-net-shape plates or extrusions, and, in longer-term, 3D printing (wire-based additive layer manufacture or metal power deposition techniques), which improve manufacturing time and allow the creation of new stiffened panel shapes [6]. However, 3D-printed products have also limitations in their structural performance [7] and, thus, there is still a long way to go to use 3D-printing for manufacturing primary structures as ISP's in aviation.

An extensive review of the research conducted during the period 2000-2012 on the buckling and post-buckling of stiffened plates and shells was presented by Ni, Prusty and Hellier ‎[8]. Their review was concentrated in two of the main research methods: numerical analysis and experimental methods. With the advent and improvement of digital computers, numerical methods, as the Finite Element Method (FEM), have been intensively used in structural analysis, including problems related to buckling and post-buckling behavior of ISP's. Indeed, the capacity of FEM solvers has significantly increased in the last years, making it possible to represent almost any kind of structure, including geometrical details if any, for a buckling analysis using solid elements at an acceptable computational cost. Despite the undeniable importance of these numerical methods in the buckling analysis of structural members, it is not surprising that some results obtained in initial studies in this area are still used in preliminary stages of the wing design, as the charts published by NACA for calculation of the critical compressive stress of idealized web- and T-stiffened panels ‎(see, e.g., Boughan and Baab [9] and Seide and Stein [10]). A possible explanation for such a long use of these charts may be attributed to the ease and saving of time that can be obtained in the initial stages of the wing design when just a rough idea about the geometry of the panel is enough to proceed to the next steps of design. However, these charts are based on simple geometries and neglect geometric details as the fillet radius used in the current design of ISP's. In fact, most of the work published up to now, involving the design and mass optimization of an ISP, consider only the following five geometrical dimensions of the panel (see Fig. 1):

  • bs − reinforcement pitch;

  • bw − reinforcement height of the web, measured from the skin mid-plane;

  • ts − skin thickness;

  • tw − web thickness;

  • a − panel length (rib spacing).

Figure 1 Cross-section of an ISP and its main dimensions. 

The addition of a new design variable, as the fillet radius R , can lead to a new optimal panel, minimizing some geometrical limitation effects, since the fillet radius reinforces the panel junctions and, thus, increase the post-buckling strength. The choice of fillet radii also influences positively the durability and damage tolerance considerations, as does the inclusion of skin crenulations and stiffener pad-ups, as mentioned by Quinn, Murphy and Cervy ‎[11].

It is important to stress that the ISP's geometry considered in this study is the blade stiffened panel, shown in Fig.1, which is usually found in current aircraft wings. This kind of geometry, including the effect of the fillet radius on the buckling performance of an ISP, was also investigated by Fenner and Watson ‎[12] and Fenner ‎[13] who used the finite element code MSC NASTRAN with 3D solid elements in their analyses. They concluded that the local buckling load could be significantly improved by taking the fillet radius into account, whereas the global (overall) buckling load had only a small improvement. However, they investigated only a particular ISP geometry and, thus, it is not possible to infer whether the results are also valid for other ISP's geometries.

Yet, regarding the optimization design of wing upper panels, Barkanov et al. [14,15] developed a weight optimization methodology employing the method of experimental design, response surface technique and both linear [14] and non-linear [15] buckling analyses. ANSYS and NASTRAN finite element codes were used to investigate the buckling behavior of three rib bays laminated composite panels with T, I and HAT-stiffeners. Due to the large dimension of numerical problems and large number of optimization tasks to be solved, the authors used only shell elements to model the panel. In this sense, it must be added that the use of complete buckling charts, as a function of all significant geometrical dimensions of a stiffened panel, is a must to a correct wing upper panel optimization process and would certainly speed up the initial stages of the wing design.

Nowadays, the search for simple and valid formulas for the buckling coefficient, ks , is still important and useful in the stability analysis of thin-walled structures. Using the Rayleigh-Ritz method, Jahanpour and Roozbahani [16] evaluated the buckling loads of a simply supported rectangular plate under biaxial and shear loads. Then, after applying regression techniques and interpolation, the authors proposed a concise formula to obtain the buckling load coefficient for a variety of in-plane load combinations. The results were validated using finite element simulations and pointed to a maximum difference of 20%. Kadkhodayan and Maarefdoust [17] also determined the buckling coefficient of thin rectangular plates under in-plane loads for several boundary conditions. They considered an elastic-plastic behavior for the plate material, described by the Ramberg-Osgood model, and used the Generalized Differential Quadrature technique to solve the buckling problem. The analyses were carried out considering two theories of plasticity: the deformation theory and the incremental theory. The results were compared with data previously published in the literature.

Considering the present lack of studies about the buckling analyses of integrally stiffened panels that are able to capture the effect of the fillet radius, the objectives of this paper are:

  1. to show that the results provided by the charts published by NACA (‎[9], ‎[10]) give in general very good results for the critical compressive stress values for several geometries of ISP's: this will be demonstrated by comparing the results obtained with the curves presented by Boughan and Baab [9] and by Seide and Stein [10] with several finite element results for different ISP's geometries (proper boundary conditions and null values for the fillet radii are considered in these analyses). In particular, it will be shown that the agreement is not quite good for some values of tw/ts and bw/bs ratios, and designers must take care with the critical compressive stress values given by the NACA charts in these specific situations;

  2. to show that the values of the critical compressive stress for local instability of ISP's may be significantly increased when one considers the effect of the fillet radius (filleted junction) along the line junction between stiffener web and skin, meaning that this parameter should also be considered in an optimal design of such thin-walled structures. Several finite element models, for fillet radii values ranging from 0 up to 31 mm, have been created in order to study the influence of this parameter on the local buckling compressive stress. It is worth noting that, although Fenner and Watson ‎[12] have already presented a similar study, only a small range of fillet radius values (from 1 up to 5 mm) was considered by them, what limits their conclusions to a quite narrow range of the non-dimensional parameter related to the fillet radius.

Summing up, following the methodology which will be described in this work, it will be possible to create new charts for the critical compressive stress, which will now depend on another non-dimensional parameter related to the value of the fillet radius used in the panel. Although it is not our intention to present such general charts in this work, some preliminary results will be presented and discussed.


2.1 The classical analysis method

The classical plate theory shows that a rectangular plate buckles in a mode which is dependent on the aspect ratio ( a/bs ), on the applied loading form and on the boundary conditions ‎(see, e.g., Timoshenko and Gere [18]). Equation (1) defines the critical buckling stress for an axially compressed plate, where Ec is the Young's compression modulus; ν is the Poisson ratio and ks is the buckling coefficient which depends on the boundary conditions of the plate.

σcr=ksπ2Ec12(1ν2)(tsbs)2 (1)

For ISP's, Eq.(1) can be used for both local and global buckling critical stresses: results for the local buckling coefficient, ks,loc , are presented in Boughan and Baab [9] chart, based on the Principles of Moment Distribution, presented by Lundquist, Stowell and Schuette [19], whereas results for the global buckling coefficient, ks,glob , are presented in Seide and Stein [10] charts for analysis of the stability under compression of simply supported rectangular plates with one, two, three, or an infinite number of identical and equally spaced longitudinal stiffeners. Later on, Seide [20] proposed a simple expression for calculation of the effective moment of inertia of the stiffeners that could be used in conjunction with the buckling charts presented previously in [10].

Figure 2 shows how the web to skin junction is represented in different models: (a) represents an idealized junction which is used in 2D finite element models (using plate or shell elements) as well as in analytical models. In these two cases only the nominal distances ( a , bs , bw , …) are used to represent the geometry of the stiffened panel and no consideration is made to the fillet radius or to its effect on the buckling load; (b) represents a square junction (or fillet-less junction) which is used in finite element models with solid elements, but it keeps on neglecting the fillet radius when compared to other panel dimensions. It is not possible, thus, to gain insight on the benefits of using higher fillet radius with such simplified finite element models; and (c) represents the real junction geometry, which considers the fillet radius dimension, and is also modeled with solid finite elements.

Figure 2 Different representations of the stiffened web to skin junction. 

It is worth noting that, although Eq.(1) refers to nominal dimensions of the stiffened panel, as those depicted in Fig.1, the effective distances that should be used in real junctions to provide a correct value of the buckling stress are not exactly the nominal ones, since the stiffeners and the fillet radius have both an effect in the real boundary conditions of the associated idealized junction. To illustrate this point, let us take, for instance, a longitudinally compressed long plate (i.e., with a relatively high aspect ratio a/bs ), simply supported in all four edges, for which the buckling coefficient predicted by the classical theory is ks 4 (see, e.g., Young and Budynas [21], table 15.2). However, if the lateral unloaded edges are considered clamped, this value increases to ks 6.96, once kept the same aspect ratio (see, e.g., Young and Budynas [21], table 15.2). Thus, the classical plate theory shows that the presence of a sturdy stiffener in the lateral edges of a plate can increase the critical buckling stress about 74%. In the case of stiffened plates, the nominal dimensions shown in Fig.1 are used for ease, but the real boundary condition is not of a simply supported plate nor of a clamped plate, but something between both of them, meaning that the correct value of ks for the idealized junction is such that 4 < ks < 6.96. In this case, the correct value of the buckling coefficient must take into account all those aspects (as the fillet radius) that cannot be considered in an idealized junction model.

2.2 The finite element analysis method

Using the software Abaqus, a parametric finite element model was created, allowing the simulation of hundreds of similar geometric models just modifying some of the input data. Parabolic hexahedron solid elements (HEX20) were chosen to represent the ISP with its fillet radius. HEX20 elements have six quadrilateral faces with 20 nodes (8 nodes in the vertices and 12 mid-edge nodes), each node having three translational degrees of freedom. The accuracy of solutions using meshes based on HEX20 elements is the highest if compared with meshes based on other 3D solid element types. The following meshing rules were used to create accurate ISP finite element models for buckling analysis based on recommendations proposed by Fenner [13] and validated by a mesh convergence analysis (presented in Appendix 1):

  • A minimum of 2 elements was used along the thickness (in skin and web meshes);

  • A minimum of 4 elements was used along the fillet radius edge;

  • An average aspect ratio of 6.7 was used and the maximum aspect ratio was 27;

  • At least 50 elements were used along the panel length, with a minimum of 5 elements along a wavelength.

Fenner, Watson and Featherston [22] discussed three different stiffened panel models to represent the infinite length panel condition: single bay, double half bay, and quad-half bay. As 3D solid FE models have many degrees of freedom and a large computational cost is involved, the panel representation chosen in this study was the double half bay with one central stiffener, because it is the simplest way to reproduce the conditions existing in an infinite length and infinite width panel, whose results provide conservative values for the buckling stress. Figure 3 represents, in the highlighted area, the portion of the panel that was modeled.

Figure 3 ISP and the double half bay panel representation (highlighted). 

For determination of the local buckling coefficient, a hypothesis of cyclic symmetry is used in the boundary condition to create an infinite panel condition (i.e., infinite length and width), irrespective the number of stringers of the panel, what provides a conservative value for ks,loc . Cyclic symmetry means that the portion of the modeled panel will be repeated infinitely in both directions by using an appropriate symmetric boundary condition. Regarding the determination of the global buckling coefficient, Seide and Stein [10] assumed only a hypothesis of infinite length panel, and the number of stringers was an input data of the problem. In this work, however, a unique hypothesis of infinite panel, in both length and width directions (i.e., considering an infinite number of stringers), are used to determine both local and global buckling coefficients.

Figure 4 shows the boundary conditions used in all simulations in this study: the load was applied by imposing displacements ( u=δx and v=δy ) at two of the four faces of the model (faces 1 and 2, respectively), so an approximated stress state of unitary compressive stress in the longitudinal direction ( σx = -1MPa) and of null stress in the transversal direction ( σy = 0) is generated in the skin in the pre-buckling state, according to Eq.(2) and Eq.(3). Besides, displacements in the longitudinal direction are constrained for all nodes at face 3, and displacements in the transversal direction are also constrained for all nodes at face 4. It means that faces 1 and 3 are constrained to rotate around y-direction, whereas faces 2 and 4 are constrained to rotate around x-direction. At last, a vertical displacement constraint ( w = 0) is imposed for all nodes at the intersection of the mid-plane of the skin thickness and the transversal plane that represents the support provided by the rib in the skin attachment. The following paragraph explains how to find the position where this vertical displacement constraint must be applied. Using all these boundary conditions, it is possible to represent the conditions in a panel of infinite length and infinite width.

{εx=δxa=1Ec[σxνσy]σx1MPaσy0δx=aEc (2)

{εy=δybs=1Ec[σyνσx]σx1MPaσy0δy=+νbsEc (3)

Figure 4 Boundary conditions and forced displacements to create a unitary compressive stress in x-direction. 

For the double half bay model, Fenner [13] still suggests to concatenate the results of the FE model associated with an odd number of half-wavelengths with those associated with an even number of half-wavelengths. Instead, in the present work three models are simulated for each panel configuration and their results are compared in order to obtain the minimum eigenvalue. The first model is the one with a simple support in the middle of the panel length representing a rib support (Fig. 5.a), which will always be related to modes with an odd number of half-wavelengths, due to the boundary conditions imposed to the model. After discovering the critical odd mode, two additional models are generated forcing buckling in even modes with one less and one more half-wavelength than the critical odd mode previously found. To create these two models with even numbers of half-wavelengths, it is only necessary to offset the middle panel simple support by a distance of a/2n , as shown in Fig. 5.b, where n is the respective even number of half-wavelengths.

Figure 5 Examples of local buckling mode shapes for (a) n = 7 and (b) n = 6 half-wavelengths. 

For the purpose of modeling the post-buckling behavior, concatenation of the results obtained with these three models is not suitable because the separation of mode shapes prevents mode jumping between odd and even mode shapes. Thus, to analyze both buckling and post-buckling behavior with a single model, Fenner, Watson and Featherston [22] have used the quad-half bay model to represent the stiffened panel, instead using the double half bay model. However, the quad-half bay representation leads to an increased solution time as the number of elements in the model is doubled, when compared to that employed in the double-half bay representation.

2.3 Simulation settings

According to Boughan and Baab [9], the local buckling coefficient stress, σcr , as given in Eq.(1), depends basically on the following panel dimensions (see Fig.1 for reference): bs , bw , ts and tw . These four dimensions can generate only three independent non-dimensional parameters which are: ts/bs - this one being explicitly given in Eq.(1) -, tw/ts and bw/bs - these two last non-dimensional parameters being used to define the value of the local buckling coefficient ks,loc as presented in the charts published by Boughan and Baab [9]. When a new dimensional parameter is included, as is the case of the fillet radius ( R ), a corresponding non-dimensional parameter must be created in order to check the possible effects on the ks,loc values. The chosen ratio will be defined as R/bwf , where bwf=bwts/2 corresponds to the free height of the stiffener. Thus, in the range 0Rbwf , the associated range for the chosen non-dimensional ratio is 0R/bwf1 .

Regarding the determination of the global buckling stress, Seide and Stein [10] has shown that it also depends on the variable a (the panel length or rib spacing) besides the aforementioned dimensions ( bs , bw , ts and tw ). Thus, considering these five dimensions, four independent non-dimensional parameters may be created, one of them being the ts/bs ratio which already appears explicitly in Eq.(1). The three remaining non-dimensional parameters used by Seide and Stein [10] to obtain charts for the global buckling stress coefficient, ks,glob , are the ratios a/bs , EcIw/bsD , and Aw/bsts , where:

  • Iw is the modified moment of inertia of the stringer about the skin mid-plane (see also Seide [20]);

  • D=Ects312(1ν2) is the flexural rigidity of the plate;

  • Aw is the stringer area.

Thus, according to Seide and Stein [10], the global buckling stress coefficient, ks,glob , can be expressed as:

ks,glob=f(abs,EcIwbsD,Awbsts) (4)

With the introduction of the fillet radius ( R ) as a new variable to be considered in the global buckling analysis, and its corresponding non-dimensional parameter R/bwf , the coefficient ks,glob is expressed in a more general form as:

ks,glob=f(abs,EcIwbsD,Awbsts,Rbwf) (5)

In this case, the parameters Aw and Iw should be corrected in order to take into account the effect of the fillet radius, as shown in Fig.6.

Figure 6 Skin (blue) and stringer (green) areas used to calculate the global buckling critical stress. 

In order to obtain new local ( ks,loc ) and global ( ks,glob ) buckling coefficient charts according to the new sets of non-dimensional parameters, some routines were created using Python, Visual Basic Excel and Abaqus enabling: (i) to automatically generate several finite element models with different geometric configurations, (ii) to run linear buckling analyses for each of these models using Abaqus, (iii) to read the finite element results, (iv) to identify local and global modes and, finally, (v) to calculate the local and global buckling coefficients for each configuration.

To distinguish between local and global buckling modes and to find the corresponding number of half wave-lengths, n , the displacements of the panel edges are compared to a cosine curve. Besides, by knowing the eigenvalue λ , and that the enforced displacements yield to a uniform unitary stress, i.e., |σx|=1 (see Eq.(2) and Eq.(3)), it is possible to find the buckling coefficient using Eq.(1), resulting:

ks=λ|σx|12(1ν2)π2Ec(bsts)2=λtsbs2π2D (6)


Two different comparisons will be presented in this section considering the results obtained with the two analysis methods presented in sections 2.1 and 2.2:

  1. Numerical results obtained with finite element models for ISP's with R = 0 (no fillet radius) will be compared with results obtained for the local buckling coefficient of web-stiffened panels following Boughan and Baab [9] methodology. The purpose of this comparison is to determine how close (or how far) the results obtained with these two alternative methods are from each other. A total of eleven values for the ratio tw/ts , and nine values for the ratio bw/bs are considered in this comparison, summing up 99 different ISP's geometries. The results are presented in section 3.1.

  2. Numerical results obtained with finite element models for three different R/bwf ratios ( R/bwf=0 , R/bwf=0.18 and R/bwf=0.32 ) will be compared to results obtained for the global buckling coefficient following Seide and Stein [10] methodology (considering an infinite number of longitudinal stiffeners). The results are presented in section 3.2.

3.1 Local buckling coefficient results

Figure 7 reproduces the local buckling coefficient ( ks,loc ) charts presented by Boughan and Baab [9] for square-junctions (that is, for junctions with R = 0) showing the dependence of ks,loc on both non-dimensional ratios tw/ts and bw/bs . The dotted curved line in the middle of the chart divides it into two distinct regions: in the upper one the skin buckles restrained by the stiffeners, whereas in the lower one the stiffeners buckle restrained by the skin. Figure 7 is valid for materials with ν = 0.3, but similar charts can be created for different values of the Poisson ratio using the same theory. For each fixed value of tw/ts , among a total of eleven values in the range 0.5 ≤ tw/ts ≤ 2.0, finite element simulations were carried out considering R = 0 (square-junctions) and nine different values of bw/bs in order to obtain the associated linear buckling coefficient ( ks,loc ) using the methodology presented in the section 2.2. Thus, a total of 99 numerical values for ks,loc were obtained for different ISP's geometries and plotted in Fig.7.

Figure 7 Local buckling coefficient chart (adapted from Boughan and Baab [9]) and finite element results. 

The main difference between the two approaches is that the Principles of Moment Distribution connects skin and stiffener by a unique line in its junction (see Fig. 2a), whereas the proposed finite element modelling consider an extended region at the stiffener base (see Fig. 2b). In order to evaluate the numerical differences among the values obtained for the ks,loc using the two methods presented in sections 2.1 and 2.2, let us define the average absolute difference (AAD)j , for a given and fixed value of tw/ts=(tw/ts)j,1j11 , as:

(AAD)j=19i=19(AD)ij=19i=19|ks,loc,1,ijks,loc,2,ijks,loc,1,ij|,fortwts=(twts)j,1j11 (7)

where ks,loc,m,ij is the local buckling coefficient value obtained by the classical ( m = 1) or finite element ( m = 2) model, for bw/bs=(bw/bs)i , 1i9 , and for tw/ts=(tw/ts)j , 1j11 , and (AD)ij is the related absolute difference, for given i and j values.

Figure 8 shows that the calculated average absolute differences are, in general, smaller than 2.0% and are only greater than that at the extremities of the chart (i.e., for tw/ts = 0.5 and tw/ts = 2.0) where the thicknesses differences between skin and stiffener are bigger. The minimum absolute average difference occurs for tw/ts = 1.2, and its related average absolute difference is (AAD)0.27% . It is important to mention that some of the geometrical dimensions used for the simulations were randomly chosen to remove any features induced by a particular choice. It is also possible to notice, from Fig.7, that the results obtained by the Principles of Moment Distribution are conservative for local buckling of filletless ISP's for the majority of the cases (the maximum non-conservative error was about 4.5%).

Figure 8 Average absolute differences for ks,loc between classical and 3D solid finite element models. 

Figure 9 shows some of the (AD)ij values plotted against the non-dimensional ratio ts/bs . The (AD)ij values considered in Fig. 9 are those obtained for tw/ts=0.5 and tw/ts=2.0 and all values of bw/bs=(bw/bs)i , 1i9 . Despite the relatively high dispersion in the results, the two dotted lines shown in Fig.9, for each of the values of tw/ts , seem to indicate a tendency of smaller (AD)ij values as the ratio ts/bs diminishes.

Figure 9 Absolute differences for the local buckling coefficient for tw/ts = 0.5 and tw/ts = 2.0. 

Thus, in the case of very thin plates (i.e., for quite small values of ts/bs ), what must be the case considered in the charts presented by Boughan and Baab [9], the absolute differences in the ks,loc values, using the classical and finite element models, are not so great (the results obtained with the charts agree reasonably well with finite element results). The use of such very thin plates was probably observed in most of the ISP's used in aircraft construction at the time when Boughan and Baab [9] charts were created. However, nowadays, aircraft have larger wingspan and lower wing main box heights leading to thicker skins and higher ts/bs values. In these cases, of higher ts/bs values and bigger differences between stiffener and skin thicknesses (i.e., higher or smaller values of tw/ts ), greater differences between the classical theory and finite element results are observed and the use of Boughan and Baab [9] charts must be taken with care, since the values provided by them are not always on the safe side.

3.2 Global buckling coefficient results

Seide and Stein [10] presented 12 stiffened panel global buckling coefficient charts for plates with one, two, three and an infinite number of longitudinal stiffeners. To validate the methodology presented in Section 2.2 concerning the determination of the global buckling coefficient, several finite element results will be compared with the corresponding values presented by Seide and Stein [10] for a plate with an infinite number of longitudinal stiffeners and a ratio Aw/bsts=0.4 .

Figure 10 shows, in the background, the curves presented in the Seide and Stein [10] chart. In this same figure, the continuous blue lines represent the ks,glob values obtained by an Excel VBA routine that uses the Rayleigh-Ritz energy method (see Seide and Stein [10]) and the modified moment of inertia introduced by Seide [20]. The marks shown in Fig.10 represent the ks,glob results obtained with a total of 114 finite element models using the methodology described in section 2.2. The dimensions of the ISP's were chosen so that some pre-defined non-dimensional ratios could be obtained. The value of the ratio Aw/bsts was considered fixed ( Aw/bsts=0.4 ); values for the ratio EcIw/bsD varied in the interval [5, 300] and a total of 8 different values for this ratio were considered ( EcIw/bsD = 5, 20, 60, 100, 150, 200, 250, 300). Finally, values for the ratio a/bs varied in the interval [1, 8], and only integer values for this ratio were considered in this range. For each fixed set of values of the non-dimensional ratios, three different fillet radius ratios were considered in the analyses ( R/bwf = 0, R/bwf = 0.18 and R/bwf = 0.32).

Symbols:for R/bwf = 0,for R/bwf = 0.18,for R/bwf = 0.32. Adapted from Seide and Stein [10].

Figure 10 Chart validation for ks,glob (for an infinite number of longitudinal stiffeners and Aw/bsts = 0.4). 

It is worth noting that in Seide and Stein [10] work, all the presented curves were limited up to ks,glob = 4, due to skin local buckling between stiffeners. However, since the local buckling coefficients (shown in Fig. 7) and the global buckling coefficients (shown in Fig. 10) are obtained separately, it is not necessary to include this limitation in the global buckling chart. Thus, the continuous blue lines presented in Fig. 10 were not limited to ks = 4, what has been confirmed by the numerical results obtained with some of the finite element simulations, also shown in Fig.10.

It is possible to observe that the finite element results and the results presented by Seide and Stein [10] are very close in all cases, for ISP's with or without fillet radii. In fact, for each fixed set of values taken for the ratios a/bs and EcIw/bsD , the differences in the ks,glob values obtained for the three different values considered for the non-dimensional ratio R/bwf are minimal (the maximum difference between them was about 2.6%).

It can be concluded, thus, that the analytical formulas provided in [10] and [18] can be effectively used to provide good estimates of the global buckling critical stress for both fillet-less and filleted ISP's. Besides, it was possible to notice that, for 0 ≤ R/bwf ≤ 0.32, the differences in the calculated values of the global buckling stress are insignificant: the maximum non-conservative difference between finite element results and analytical results for all the 114 studied cases was about 2.1%. Regarding the non-dimensional parameters that control the global buckling behavior in filleted junctions, these results allow us to conclude that the ratios a/bs , EcIw/bsD and Aw/bsts are the most important ones, since the fillet radius effect is also taken in account in the stringer properties Iw and Aw .


A typical configuration was chosen to study the behavior of both local and global stress buckling values as the fillet radius dimension changes. Based on Figure 1, the ISP dimensions taken in this case study are: bs = 90mm, h = 35mm, ts = tw = 3.5mm, and a = 450mm, resulting in the following non-dimensional ratios: bw/bs = 0.369, tw/ts = 1, ts/bs = 0.039 and a/bs = 5. The material of the panel is aluminum with Young's compression modulus Ec = 73.1 GPa and Poisson's ratio ν = 0.33. Following the methodology presented in section 2.2, a total of 32 finite element models was created, each one with a different value of the fillet radius in the range 0 ≤ R ≤ 31mm, considering small increments of 1 mm for each fillet radius value.

Figure 11 shows the values of the global buckling coefficient, ks,glob , obtained for the 32 different values taken for the fillet radius, considering the two analysis methods: i) the classical formulation based on Seide and Stein [10] and on Seide [20] works (as discussed in section 2.1) and ii) the finite element method (as discussed in sections 2.2 and 2.3). The value of the panel cross-section area, as the fillet radius increases, is also shown in Fig. 11.

As it can be seen in Fig. 11, ks,glob has a smooth decrease as the fillet radius increases, for the studied configuration. However, the increase in the fillet radius leads also to a corresponding increase in the panel cross-section area, which surmounts the opposite effect in the ks,glob reduction. Considering the combination of these two effects, there is in fact an increase in the global buckling load as the fillet radius increases, as will be discussed later. Still concerning the global buckling stress, the values found by the finite element simulations are close to those predicted by the analytical theory presented by Seide and Stein [10], as it can be noted in Fig. 11, as well as in the results presented in section 3.2. However, it is important to register that the estimates given by the analytical theory may be non-conservative in some cases, as it can be seen in Fig. 11 in the range 0 ≤ R/bwf ≤ 0.5.

Figure 11 Effect of the fillet radius of a typical ISP geometry on the global buckling coefficient. 

Regarding the effect of the non-dimensional parameter R/bwf on the local buckling coefficient, Fig. 12 shows the behavior of ks,loc as the ratio R/bwf increases. In this graph, two non-dimensional ratios were fixed ( bw/bs = 0.369 and tw/ts = 1) and five different values for the ratio ts/bs were considered, namely ts/bs = 0.020, ts/bs = 0.033, ts/bs = 0.039, ts/bs = 0.058 and ts/bs = 0.120. For each set of fixed non-dimensional ratios ( bw/bs , tw/ts and ts/bs ), several linear buckling analyses using finite element models (as described in sections 2.2 and 2.3) were carried out for different values of the non-dimensional ratio R/bwf , providing the values of ks,loc that generated the five different curves depicted in Fig. 12 (each one for a given fixed value of the ratio ts/bs ) as functions of the ratio R/bwf . For the five studied cases, it can be noticed that ks,loc has an almost parabolic variation as the radius increases. Yet, it can be noticed that, as the fillet radius increases, its effect on the critical local buckling load becomes even more pronounced, mainly when small ts/bs ratios are considered.

Figure 12 Effect of the fillet radius in ISP local buckling coefficient (for bw/bs = 0.369 and tw/ts = 1). 

It can be noticed, also, that every curve starts at the same ks,loc value (obtained for R/bwf = 0), which can be obtained in Fig. 7 chart and provides ks,loc 3.971 (considering in this case ν = 0.33). For this specific case (i.e. for R/bwf = 0 and tw/ts = 1), the ks,loc practically does not depend on ts/bs ratio, as discussed in section 3.1. However, when the ratio R/bwf is considered in the modeling, there is a noticeable dependence of the local buckling coefficient on the ratio ts/bs . Yet, it can be seen in Fig. 12 that the ks,loc value may be even higher than 6.96 for some ISP filleted junction configurations, contradicting the classical theory results. In other words, this dependence of ks,loc on ts/bs ratio shall be seen as a way to correct the exponent 2 of the ratio ts/bs in Eq.(1), when ratios R/bwf > 0 are considered.

It is also possible to show that any two different configurations that have the same values for the four non-dimensional ratios presented in section 2.3 (i.e., ts/bs , tw/ts , bw/bs and R/bwf ) will present the same value of the local buckling coefficient, ks,loc , what shows that these four non-dimensional parameters are enough to represent the local buckling behavior of any ISP.

Regarding the increase in the capacity of the ISP, considered in this case study, to sustain compressive loads, Fig.13 shows the behavior of the total compressive load applied to the panel (leading either to global or local buckling failure modes) as the fillet radius increases. As mentioned before, it can be seen that, despite the decrease of the global buckling coefficient as the non-dimensional parameter R/bwf increases, the capacity of the panel to sustain higher compressive loads (considering global buckling) increases monotonically as R/bwf increases: from Pcr,glob = 118 kN (for R/bwf = 0) to Pcr,glob = 188 kN (for R/bwf 1), what represents an increase of 59% in the global buckling load. In the same way, a similar behavior is achieved in the local buckling behavior: the critical local buckling load increases from Pcr,loc = 170 kN (for R/bwf = 0) to Pcr,loc = 1,880 kN (for R/bwf 1), what gives an increase of 1006% in the local buckling load.

Figure 13 Effect of R/bwf on the global and local bucking loads in the studied ISP. 


The validity of the NACA buckling charts of the 1940's for integrally stiffened panels (ISP's) has been verified in this work by using 3D finite element models, with an infinite length and infinite width panel condition, for several different values of the non-dimensional parameters that influence the local and global buckling failures. For local buckling, the numerical results obtained with the finite element models for non-filleted ISP's presented a very good correlation with the results presented in the Boughan and Baab [9] chart. For the global buckling case, several finite element simulations were run for different values of the non-dimensional parameters, and the results showed a quite good agreement with those presented by Seide and Stein [10] and it was shown (for a specific case, i.e, for an infinite number of longitudinal stiffeners and a non-dimensional ratio Aw/bsts = 0.4) that the influence of the fillet radius value on the global buckling coefficient is not significant.

Additional results provided in this work showed that the influence of the fillet radius in the local buckling of ISP's is very important and cannot be neglected. It was also shown that the increase in the local buckling coefficient varies in an almost parabolic shape as the fillet radius increases, and it was possible to notice that the percentual increase in the local buckling stress was much higher than the percentual increase in the cross-sectional area for the studied case. The increase in the local buckling coefficient due to the increase in the fillet radius may be attributed to the consequent reductions in the skin and stiffener effective lengths ( bs and bw , respectively).

The present study was limited to ISP's of a blade stiffener type, but the idea can also be applied to other modern panel concepts as iso-grid, sub-stiffeners, and buckling containment features among others. Probably, the fillet radius inclusion could increase even more the local buckling load of these solutions and its benefit can only be noticed with the use of 3D solid FE models or with experimental tests. The effects of the fillet radius inclusion on the post-buckling behavior of ISP's were not considered in this work and are left to a further research.


[1] Munroe, J., Wilkins, K., Gruber, M. (2000). Integral Airframe Structures (IAS) - Validated feasibility study of integrally stiffened metallic fuselage panels for reducing manufacturing costs. NASA Contractor Report. NASA/CR-2000-209337. [ Links ]

[2] Akl, W., El-Sabbagh, A., Baz, A. (2008). Optimization of the static and dynamic characteristics of plates with isogrid stiffeners, Finite Elements in Analysis and Design, 44: 513-23. ]

[3] Quinn, D., Murphy, A., McEwan, W., Lemaitre, F. (2010). Non-prismatic sub-stiffening for stiffened panel plates-Stability behavior and performance gains, Thin-Walled Structures, 48: 401-13. ]

[4] Watson, A., Featherston, C.A., Kennedy, D. (2007). Optimization of postbuckled stiffened panels with multiple stiffener sizes. In: 48th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Honolulu, USA, vol. AIAA-2, 8 p. [ Links ]

[5] Özakça, M., Murphy, A., van der Veen, S., (2006). Buckling and post-buckling of sub-stiffened or locally tailored aluminium panels. In: 25th Int. Congress of the Aeronautical Sciences (ICAS-2006), Hamburgo, 16 p. [ Links ]

[6] Houston, G., Quinn, D., Murphy, A., Bron, F., (2016). Wing panel design with novel skin-buckling containment features. Journal of Aircraft, 53: 416-26. ]

[7] Joshi, S.C., Sheikh, A.A. (2015). 3D printing in aerospace and its long-term sustainability. Virtual and Physical Prototyping, 10: 175-85. ]

[8] Ni, X-Y, Prusty, B.G., Hellier, A.K. (2015). Buckling and post-buckling of isotropic and composite stiffened panels: a review on analysis and experiments (2000-2012). Trans RINA, Intl J Marit Eng, 157. ]

[9] Boughan, R.B., Baab, G.W. (1944). Charts for calculation of the critical compressive stress for local instability of idealized web- and T-stiffened panels. NACA-WR-L-204, NACA-ARR-L4H29. Hampton: Langley Research Center, 16p. [ Links ]

[10] Seide, P., Stein, M. (1949). Compressive buckling of simply supported plates with longitudinal stiffeners. NACA-TN-1825, Langley Aeronautical Laboratory, 24p. [ Links ]

[11] Quinn, D., Murphy, A., Cervi, L. (2011). Fatigue performance of aircraft panels with novel skin buckling containment features. Proc. Inst. Mech. Eng. - Part G: J Aerospace Eng, 225: 791-806. [ Links ]

[12] Fenner, P.E., Watson, A. (2012). Finite element buckling analysis of stiffened plates with filleted junctions. Thin-Walled Structures, 59: 171-80. ]

[13] Fenner, P.E. (2014). Comparing the accuracy of VICONOPT to FEM for analysing aircraft wing skin type panels. Master Thesis, Loughborough University, United Kingdom, 230p. [ Links ]

[14] Barkanov, E., Ozolins, O., Eglitis, E., Almeida, F., Bowering, M.C., Watson, G. (2014). Optimal design of composite lateral wing upper covers. Part I: Linear buckling analysis. Aerospace Science and Tech, 38: 1-8. ]

[15] Barkanov, E., Eglitis, E.,Almeida, F., Bowering, M.C., Watson, G. (2016). Optimal design of composite lateral wing upper covers. Part II: Nonlinear buckling analysis. Aerospace Science and Tech 51: 87-95. ]

[16] Jahanpour, A., Roozbahani, F. (2016). An applicable formula for elastic buckling of rectangular plates under biaxial and shear loads. Aerospace Science and Tech, 56: 100-11. ]

[17] Kadkhodayan, M.,Maarefdoust, M. (2014). Elastic/plastic buckling of isotropic thin plates subjected to uniform and linearly varying in-plane loading using incremental and deformation theories. Aerospace Science and Tech, 32: 66-83. ]

[18] Timoshenko,S.P., Gere, J.M. (1961). Theory of Elastic Stability, 2nd ed, McGraw Hill. [ Links ]

[19] Lundquist, E.E., Stowell, E.Z., Schuette, E.H. (1943). Principles of moment distribution applied to stability of structures composed of bars or plates. NACA-Tech. Report 809, Hampton: Langley Research Center, 14p. [ Links ]

[20] Seide, P. (1953). The effect of longitudinal stiffeners located on one side of a plate on the compressive buckling stress of the plate-stiffener combination. NACA-TN-2873, Langley Aeronautical Laboratory, 67p. [ Links ]

[21] Young, W.C., Budynas, R.G. (2002). Roark's formulas for stress and strain, 7th ed, McGraw-Hill. [ Links ]

[22] Fenner, P., Watson, A., Featherston, C. (2016). Modelling infinite length panels using the finite element method. Int J Structural Stability and Dynamics, 16. ]


Several finite element models were simulated in this work to evaluate the influence of the fillet radius on the buckling of ISP's. A description of the main features of the finite element modeling used to determine the buckling loads was described in section 2.2. Besides, a mesh sensitivity analysis was also performed to find the best element size at an acceptable computational cost; a summary of such analysis is provided in this appendix.

As shown in the preceding sections, the local buckling coefficient ks,loc depends on four non-dimensional ratios (namely ts/bs , tw/ts , bw/bs and R/bwf ), whereas the global buckling coefficient ks,glob depends mainly on three non-dimensional ratios (namely a/bs , EcIw/bsD and Aw/bsts ), since the fillet radius effect (related to R/bwf ) is already taken in account in the stringer properties Iw and Aw . Thus, in the finite element models used in this work, some dimensions were randomly chosen. Several tests were made involving hundreds of configurations and different values of the fillet radius ( R ), width ( bs ), length ( a ), stringer height ( bw ) and thicknesses ( tw and ts ). Studies about the mesh sensitivity for all these different configurations were also carried out, partially based on suggestions on meshing rules proposed by Fenner [13].

The mesh convergence analysis allowed the identification of the best element type to be used, as well its size in the panel cross-section and its length in the longitudinal direction. It was identified that the best element to be used in this case is the parabolic hexahedron solid element, because it is a high order element and, thus, less elements are needed in the longitudinal direction, when compared to meshes based in tetrahedral or pentahedral elements.

Several types of cross-section partitions were also tested, as shown in Fig. 14(a). However, it was noticed that the use of Abaqus automatic mesher using all sort of partitions shown in Fig.14(a) could create many distorted elements. To avoid that, only two partitions were kept: one in the skin thickness mid-plane and another in the panel mid-length. These two partitions are required to create the vertical displacement constraint ( w = 0) as described in section 2.2. The option for a specific rule for element sizing was chosen observing that: (i) a convergence process is achieved and (ii) distorted elements, as depicted in Fig. 14(b-c), are avoided. Figure 14(d) shows that a minimum of two elements along the thickness, in addition to an adequate element size, can prevent the creation of distorted elements.

Figure 14 ISP cross-section representation: studied partitions and meshing element sizes. 

Figure 15 shows that, for the specific case analised in section 4 (i.e., for that typical configuration), an element with 2 mm of size in the cross-section has a negligible deviation from its mesh convergence asymptotic result. Based on this mesh convergence analysis, the four meshing rules described in section 2.2 were considered for modeling all other ISP's geometries.

Figure 15 Mesh convergence analysis for a typical configuration. 

Received: March 13, 2019; Revised: June 05, 2019; Accepted: June 17, 2019; Published: June 19, 2019

* Corresponding author

Author’s Contributions:

Conceptualization, FG Garcia and R Ramos Jr; Methodology, FG Garcia and R Ramos Jr; Formal analysis and investigation, FG Garcia and R Ramos Jr; Validation: FG Garcia; Writing - original draft, FG Garcia; Writing - review & editing, FG Garcia and R Ramos Jr; Supervision, R Ramos Jr.


Pablo Andrés Muñoz Rojas.

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License