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

author https://doi.org/10.1590/1679-78255549 Abstract 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


INTRODUCTION
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][3][4][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): • s b − reinforcement pitch; • w b − reinforcement height of the web, measured from the skin mid-plane;  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 panels 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, s k , 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: 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 w s t t and w s b b ratios, and designers must take care with the critical compressive stress values given by the NACA charts in these specific situations; 1. 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 nondimensional 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.

The classical analysis method
The classical plate theory shows that a rectangular plate buckles in a mode which is dependent on the aspect ratio ( s a b ), on the applied loading form and on the boundary conditions (see, e.g., Timoshenko and Gere [18] [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 , s b , w b , …) 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. 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 s a b ), simply supported in all four edges, for which the buckling coefficient predicted by the classical theory is s k  4 (see, e.g., Young and Budynas [21], table 15.2).
However, if the lateral unloaded edges are considered clamped, this value increases to s k  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 s k for the idealized junction is such that 4 < s k < 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.

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 panels 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. 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 , s loc k . 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 (  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 2 a n , as shown in Fig. 5.b, where n is the respective even number of 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.

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  , the associated range for the chosen non-dimensional 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 ( s b , w b , s t and w t ). Thus, considering these five dimensions, four independent non-dimensional parameters may be created, one of them being the s s t b ratio which already appears explicitly in Eq.
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 In this case, the parameters w A and w I should be corrected in order to take into account the effect of the fillet radius, as shown in Fig.6. panels (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.,

VERIFICATION OF THE ANALYSIS METHOD
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.   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

Local buckling coefficient results
where , , ,   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 w s t t = 0.5 and w s t t = 2.0) where the thicknesses differences between skin and stiffener are bigger. The minimum absolute average difference occurs for w s t t = 1.2, and its related average absolute difference is ( . 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%).
Despite the relatively high dispersion in the results, the two dotted lines shown in Fig.9, for each of the values of w s t t , seem to indicate a tendency of smaller ( ) ij AD values as the ratio s s t b diminishes.  [9] charts were created. However, nowadays, aircraft have larger wingspan and lower wing main box heights leading to thicker skins and higher s s t b values. In these cases, of higher s s t b values and bigger differences between stiffener and skin thicknesses (i.e., higher or smaller values of w s t t ), 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.

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   [10].
It is worth noting that in Seide and Stein [10] work, all the presented curves were limited up to , s glob k = 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 s k = 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 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 ≤ wf R b ≤ 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

ADDITIONAL RESULTS AND DISCUSSION
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:  Figure 11 shows the values of the global buckling coefficient, , s glob k , 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, , s glob k 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 , s glob k 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 ≤ wf R b ≤ 0.5.  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

CONCLUSIONS
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 w s s A b t = 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 ( s b and w b , 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.

APPENDIX 1: MESH CONVERGENCE ANALYSIS
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 , 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 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.