Buckling is a key sizing driver on the design process of composite panels used in the aerospace industry. Therefore, a wide variety of methods have been developed to determine the critical buckling load of structures, from simple analytical solutions to more complex finite element methods.
Typically, during the structural sizing process the methodology to be used for panel stability is gradually evolving from simplified closed formulae using basic panel idealization thru mathematical energy based methods with stiffened panel idealization and concluding with complex finite analysis models including detailed component idealization called virtual testing.
The state of the art of mathematical approach is limited to rectangular solid stiffened panels. In this paper an evolution from solid to holed panels, including nonrectangular panels, submitted to any inplane load combination is proposed.
Several hole shapes and panels with anisotropic layups are considered. Panel thickness evolution is not considered in the theory development, however the capability of doing panel assemblies using this methodology override this limitation.
A brief energy method basis is described below.
The RayleighRitz procedure is a well know energy method based on the principle of virtual work, which states that a continuous body is in equilibrium if and only if the virtual work of all forces, internal and external , acting on the body is zero under a virtual displacement (Ref [1]):

(1) 
The virtual work done by actual body forces and surface forces moving through the virtual displacement is:

(2) 
And the total internal virtual work stored in the body is:

(3) 
When the equation is applied to elastic bodies, the strain density function exists such that and the principle of virtual work can be transformed into a minimization problem known as Principle of Minimum Total Potential Energy.

(4) 
The sum is the total potential energy and the body is in equilibrium when the variation of this magnitude is null, reaching a minimum.
The RayleighRitz procedure consist on expressing the potential energy functional as a function of the plate displacements. In the literature, this deflection shape is approximated by a trigonometric series which it must comply with the boundary conditions. For instance for simply supported conditions trigonometric series of sins are recommended.
But other formulations can be used like Chebyshev polynomials (Ref [2]). It has been considered due to their superior rate of convergence and the accumulated experience with them in previous works.

(4) 
For the sake of simplicity, only the terms depending on the perpendicular displacements are retained on this paper, but the method is applicable to all the degrees of freedom of the plate including inplane displacements.
Introducing this deflection shape function into the potential energy function defined in the previous section a function depending only on the solution unknowns is obtained. The equilibrium of the structure is reached, according to the principle of minimum total potential energy, by obtaining the minimum of the function with respect to the unknown coefficients.
As it was explained in the introduction, the energy methods are typically limited to basic geometries due to the difficulty of finding a shape function that defines the panel deflection along the integration domain while fulfilling the boundary conditions. To solve this constraint, it is necessary to carry out a transformation of the domain of integration, switching from to , using a variable change and a panel division that will allow to convert the holed panel into an assembly of 4 basic rectangular panels, which constitutes one of the innovations of this work.
The panel will be divided in four subpanels so that the holed edge is split in four edges, each one belonging to one of the subpanels. The division is made joining the center of the hole with the four corners of the panel, as shown in next figure (1).
Each of these subpanels is transformed into an isoparametric square with sides of length 2, with the new axes passing through the midpoint of each edge:Once the change of variable has been defined to simplify the problem, the RayleighRitz procedure defined in the previous chapter is now applied to the assembly of the panels. Therefore, the Total Potential energy result from the internal strain energy and the work performed by the external loads:

(5) 
For a laminated plate of anisotropic material, the internal strain energy including the coupling effect is:

(6) 
Taking into account the classical laminate theory (Ref [3]), the strain at any ply of the laminate can be expressed as the addition of the midplane strain and the curvature contribution:

(7) 

(8)

Where the midplane strain and the curvature can be expressed in terms of the displacements field as follows:
Joining the last two equations into equation (9), leads to:

(10)

Where the index k denotes the subpanel identification and [A], [B] and [D] are the laminate stiffness matrices according to the classical laminate theory.
The same way, the work done by the external loads for each subpanel is defined as follows:

(11)

As previously indicated, the minimum potential energy principle leads to an optimization problem to obtain the solution. The boundary conditions are imposed by means of a set of linear equations that the coefficients of the shape function have to fulfill:

(12) 
Where {S} is a compact form of all the unknown coefficients, d represents the independent term and [C] is the matrix of coefficients derived from the boundary equations imposition.
The minimization of function together with the set of constraints in the unknown coefficients lead to a constrained optimization that can be transformed into an unconstrained one by means of the Lagrange multipliers method. This method reduces the problem with n variables and I constraints to a new system with variables without any constraints.
A new objective function is defined adding to the original function the constraints equations multiplied by arbitrary parameters (Lagrange multipliers):

(13) 
In this case, the structure is composed by 4 interconnected panels, therefore, the modified total potential energy has to be computed for the complete structure as the sum of the total potential energy of each individual panel in addition to the terms coming from the Lagrange multipliers method:

(14) 
Being the expression of the strain energy and the external forces work quadratic in the coefficients, once they bare derived with respect to these coefficients a linear system of equations is obtained, that can be expressed in matrix form as follows:

(14)  

(15) 
Some of the unknown coefficients will not be independent but can be derived from the relations imposed by the boundary conditions. To identify the dependent coefficients, the rectangular matrix [C] used to define the restrictions can be transformed by reordering the unknown coefficients rewriting the set of equations as shown below:

(16)  

(17) 
is a square nonsingular matrix and are the dependent coefficients that can be derived from the independent coefficients :

(18) 
Substituting and removing the following system is obtained to determine :

(19) 
Where:

(20)  

(21)  

(22) 
The buckling onset of the structure is obtained from the indeterminate solution of this system, giving rise to an eigenvalue problem providing the buckling load level corresponding to each buckling mode (eigenvalues) and the corresponding buckling deflection shape (eigenvectors).

(23) 
In order to verify the correctness of the method, the results for a single panel subjected to different loading combinations and imposing different boundary conditions are compared against FE results (NASTRAN).
The following table (1) summarizes the geometry and properties used:
Parameter  Units  Value 
Ply thickness  mm  0.184 
Layup    [45,45,0,90]_{s} 
Loading  N/mm  1 
The panel is always simply supported for the following scenarios, except for the free edge section where one of the sides has been left free.
The method has been implemented in a software program that has been called BUHO which stands for “BUckling of composite panels with HOles”.
The next table (2) show the first two buckling modes for a 400x200 panel with longitudinal compression loading ( = 1 N/mm) and a hole located in the center of the panel, with a radius varying from 10 mm to 70 mm.
Radius  10  30  
Mode 1  Mode 2  Mode 1  Mode 2  
BUHO  
NASTRAN  
Error  1.5%  0.6%  0.9%  0.7% 
Radius  50  70  
Mode 1  Mode 2  Mode 1  Mode 2  
BUHO  
NASTRAN  
Error  1.2%  0.8%

1.4%  0.8% 
Next table (3) shows a trapezoidal panel loaded with N_{x}, N_{y} and N_{xy }and the same radius as the previous scenario. In this case it can be seen that, when the hole is too close to the edge, the division of the subpanels is too deformed and could potentially lead to a higher error between Nastran and the method.
Radius  10  30  
Mode 1  Mode 2  Mode 1  Mode 2  
BUHO  
NASTRAN  
Error  2.5%  1.6%  3.1%  1.8% 
Radius  50  70  
Mode 1  Mode 2  Mode 1  Mode 2  
BUHO  
NASTRAN  
Error  4.3%  3.3%  6.4%  4.8% 
One of the advantages of using Chebyshev polynomials is that the boundary conditions are not imposed, as it was explained in the previous section. In the next table (4) it will be shown an example where one of the edges has been left free to illustrate that the results are equally good. The panel is loaded with N_{x}.
Radius  10  30  
Mode 1  Mode 2  Mode 1  Mode 2  
BUHO  
NASTRAN  
Error  3.0%  1.3%  2.1%  2.8%  
Radius  50  70  
Mode 1  Mode 2  Mode 1  Mode 2  
BUHO  
NASTRAN  
Error  3.7%  4.5%  4.4%  4.6% 
The most significant advantages of using the program instead of Nastran and conventional methods are:
In order to account for all real scenarios, the number of different shapes of holes should be increased. Accessibility, associated structure and structural requirements are some of the variables when deciding the shape of the hole to be applied on the panel. The methodology to include these shapes will be the same, applying a different change of variable for each case.
In some cases, the panel needs to be reinforced to increase its stiffness, leading to the addition of longitudinal or transverse stiffeners, as it is shown in the next figure (3):
This scenario would be approached by placing the stiffener in one of the edges, so when the panel division is performed, the stringer is not deformed. In addition, the strain energy and the external forces work for each stiffener should be added to the overall solution.
The developed method in this work considers the Classical laminate theory in which the transverse shear effect is neglected. It happens when a high stiffness is supposed, which means that the normal of the midplane remains straight according to the curvature radius of the midplane after deflection. The hypothesis is valid for thin panels, but implies that results are nonconservative. Therefore, buckling could occur sooner than expected when the relation between the width of the panel and its thickness is high enough. This effect is not significant for the cases studied in the project; however, the addition of this step constitutes a natural next step.
[1]  J. Reddy, Energy Principles and variational methods in applied mechanics, Texas: Texas A&M University, 1984. 
[2]  D. H. John.C. Mason, Chebysehv Polynomials, 2002. 
[3]  S. G. Kollár LP, Mechanics of Composite Structures, Cambridge: Cambridge University Press, 2003. 
[4]  M. H. Collado, Buckling Insability of Composite Material Panels with Holes, Master Thesis, 2018. 
Published on 20/06/24
Accepted on 27/01/24
Submitted on 31/05/23
Volume 08  COMUNICACIONES MATCOMP21 (2022) Y MATCOMP23 (2023), Issue Núm. 5  Materiales y Estructuras, 2024
Licence: Other
Are you one of the authors of this document?