Composite plate buckling. Best practices in Airbus Commercial

RESUMEN: En los últimos 20 años el poder computacional se ha incrementado exponencialmente lo que ha dado lugar a que se puedan construir modelos de elementos finitos muy detallados que representan de manera fiable la estructura que es objeto de dimensionado. Sin embargo, esta potencia de cálculo también se puede usar para resolver ecuaciones complejas que hasta este momento sólo podían resolverse basándose en simplificaciones.

Calcular el pandeo de una estructura es un factor clave si se quiere que esta sea óptima. Es esencial determinar si este pandeo es inestable, conduciendo a un colapso estructural, o si por el contrario existe capacidad de redistribución y post-pandeo. Diseñar una estructura con capacidad de redistribución de carga optimiza el peso incrementando la seguridad y, por tanto es el principal objetivo del dimensionado de algunos componentes estructurales tales como cajones de torsión o revestimiento de fuselajes.

Durante años Airbus ha desarrollado soluciones de pandeo, desde las típicas soluciones cerradas a complejos desarrollos energéticos. Este artículo resume la evolución de dichos métodos con especial énfasis en los desarrollos energéticos y en las prácticas seguidas por Airbus para su desarrollo. Finalmente se incluye un resumen de potenciales mejoras y desarrollos posibles para el futuro.

Palabras clave: Métodos analíticos, Rayleigh-Ritz, pandeo, métodos de elementos finitos, MEF

ABSTRACT: During the last 20 years computational progress has significantly increased the capacity to determine the behavior of structures under complex loading conditions through finite element models analysis. This computational power can be used to build complex detailed finite element models but also to resolve more complex equations which were not possible to be coded before.

Accurate buckling prediction is still a key factor on optimized structural designs. The accurate determination of the onset of buckling is absolutely essential to determine the initial point of instability and the subsequent post-buckling capability. Determining if the buckling will lead to immediate failure or the loading will be able to be redistributed afterwards is a major point to be considered. Designing a structural component with loading redistribution capability optimize the weight and increase the accuracy, so is a key objective of designers on structural components like torque boxes or fuselage skin.

During years Airbus has been focused on developing accurate buckling predictions from close form solutions to complex energy methods development. This paper aims to summarize the evolution of such methods with special emphasis in energy methods and the best practices followed in Airbus for development. On top a summary of potential room for improvement and further evolution will be proposed.

Keywords: Analytical method, Rayleigh-Ritz, buckling, finite element method, FEM

1. Introduction

There are typically two ways to approach the problem of the panel stability analysis. One solution will be the use of Finite Element (FE) solution. There are several well-known commercial FE solvers on use in the industry, NASTRAN, ABAQUS, etc. Second solution is the semi-analytical approach development based on energy methods (EM).

Both approaches are complementary and can be used at different maturity status of the design based on following assessments:

• Energy methods, (EM), are more CPU efficient than FE solutions keeping same level of accuracy.
• Massive calculations can be launched with EM but not with FE. On top on EM design criteria discriminating local buckling versus global ones can be fixed and also it can be integrated in complex analytical processes which it is not possible with FE solutions.
• Any geometry or boundary condition can be analyzed with FE, however EM is limited to rectangular geometries and defined boundary conditions, like simply supported or clamping. We will see later that this is in evolution now.

In conclusion, Airbus have considered that the development of both approaches leads to the most efficient solution. In this paper we will be limited to the internal developments for energy method approach (analytical solutions).

2. Computational power and method development

During the last 40 years the computational power has been developed overriding all the different limits law pronounce over such period of time [1], see figure 1:

Figura 1. Computational power evolution

This has permitted that methodologies such as described in [2] have become a reality as a computer program, evolving from tabulated graphs, such as can be seen in figure 2, to current computer programs that can be executed as many times as the user requires, calculating each time the specific equations linked to the exact geometry and loads given. The effort, interpretation, interpolations that an engineer should do are reduced to the minimum, increasing the engineering efficiency.

Figure 2. Typical charts used to predict onset of buckling for different geometries and load combinations

3. Energy Method

Energy method is based on the general theorem of the equilibrium of mechanical systems: the potential energy of the system reaches an extreme value at its equilibrium position.

Two states of the panel are considered: the equilibrium state that is immediately previous to the buckling one, (at which the middle surface of the panel remains straight under the action of the applied loads); and the equilibrium state that is immediately subsequent to the buckling one (at which the middle surface of the panels is bent). Calling Π_I the potential energy of the first panel state and Π_ΙΙ the potential energy of the second one, the relationship between both of them is given by:

ΠI = ΠII + δΠ

(1)

According to the potential energy theorem, the state I shows a stable equilibrium if the next condition is fulfilled:

ΠI < ΠII (δΠ > 0)

(2)

The equilibrium will be unstable if:

ΠI > ΠII (δΠ < 0)

(3)

Critical loads are those for which the state I is no longer stable, i.e:

δΠ = 0

(4)

The total energy of the system can be split in: the strain energy (U) and the potential energy of the applied loads (V).

To illustrate, the strain energy (U) of a simple flat panel can be expressed as follow:

$\delta U={\frac {1}{2}}\int _{0}^{a}\int _{0}^{b}{\left\{{\epsilon }_{x}^{0}\right\}}^{T}\left\{A\right\}\left\{{\epsilon }_{x}^{0}\right\}+2{\left\{{\epsilon }_{x}^{0}\right\}}^{T}\left\{B\right\}\left\{{\epsilon }_{x}^{0}\right\}+{\left\{\kappa \right\}}^{T}\left\{D\right\}\left\{\kappa \right\}dxdy$

(5)

Where [A], [B] and [D] are the laminate stiffness matrices according to the classical laminate theory.

The potential energy of the applied loads (V) will be:

$\delta V={\frac {1}{2}}\int _{0}^{a}\int _{0}^{b}\left[\lambda {F}_{x}{\left({\frac {\partial w}{\partial x}}\right)}^{2}+\lambda {F}_{y}{\left({\frac {\partial w}{\partial y}}\right)}^{2}+2\lambda {F}_{xy}{\frac {\partial w}{\partial x}}{\frac {\partial w}{\partial y}}\right]dxdy$

(6)

Where ʎ indicates the factor applied to the loads which makes the panel instable.

Integrating and equating equations (5) and (6) we can determine the deformed shape function (w(x,y)), and we can obtain the critical buckling load. However this solution cannot be solved mathematically and the energy method propose an approximated solution assuming a deformed shape function as follows:

1. A displacement set is assumed, such as complies automatically with boundary conditions (see 3.1). Typically Fourier series function are used using several coefficients A_ij. The more coefficients used the closer the solution is:

$w(x,y)=\sum _{i}^{}\sum _{j}^{}{A}_{ij}\cdot {w}_{ij}(x,y)$

(7)

2. The previous governing equations are developed making the derivatives on (5) and (6)

3. The displacement functions are those fulfilling (4), equating (5) and (6) and minimizing ${\textstyle \delta \Pi }$

${\frac {\partial \left(\delta \Pi \right)}{\partial {A}_{ij}}}=0$

(8)

4. This operation leads to a set of linear homogeneous equations of the unknown coefficients and the problem is converted in an eigenvalue problem, where the eigenvalue ( ${\textstyle \lambda }$ ) will be the onset of buckling and the eigenvector is the deformed shape:

$[[\delta U]+\lambda [\delta V]]\lbrace A\rbrace =\lbrace 0\rbrace ~$

(9)

3.1. Rayleigh-Ritz approach vs any boundary conditions

The Rayleigh-Ritz approach implies defining an approximate solution representing a displacement field that automatically complies with the boundary conditions.

For example, for a simply supported edge a simple Fourier series development can be proposed:

$w(x,y)=\sum _{i=1}^{{n}_{x}}\sum _{j=1}^{{n}_{y}}{A}_{ij}\cdot sin\,\left({\frac {i\pi x}{a}}\right)\,sin\,\left({\frac {j\pi y}{b}}\right)\,$

(10)

The greater the number of terms is (n_x and n_y), the more accurate the solution will be. These functions satisfy the geometrical boundary conditions (w = M = 0 null displacements and moments at the four edges).

Nevertheless, different boundary conditions are sometimes needed, one free edge for example. For those cases typical Fourier series are no more adequate [10]. To overcome this problem, and define various boundary conditions, these have to be imposed as a set of nc linear equations that must be subsequently satisfied by the solution:

${\textstyle {e}_{i}\,\left(c1,c2,\ldots ,cN\right)=\,{d}_{i}\,for\,i=}$$1\ldots {n}_{c}$

(11)

The minimization problem becomes a conditional optimization that will be solved by means of the well-known Lagrange multipliers method.

The Lagrange multiplier method consists in multiplying each of these constraint equations with arbitrary parameters ʎi and adding the result to the minimum potential energy equation, obtaining a new objective function Π₂

${\Pi }_{2}=\Pi +\sum _{}^{}{\lambda }_{i}\cdot ei\,(c1,c2,\ldots ,cN)\,$

(12)

The problem is now more complex, new equations are added but the flexibility and accuracy to represent different problems is significantly improved.

With this mathematical approach it is possible to find analytical solutions to buckling problems as shown in figure 3. An important remark is that with this solution only local buckling problems can be solved. For global stability another solution must be used, FE solvers.

Figure 3. Some of stability problems solved with Lagrange multipliers

3.2. Stiffener idealization

In order to represent properly the stiffeners attached to the panels, its contribution need to be added to the strain and potential energy of the system:

$\delta {U}_{stiff}={\frac {1}{2}}\sum _{i=1}^{{n}_{r}}\int _{0}^{a}\left[{E}_{stiff}{I}_{stiff}{\left({\frac {{\partial }^{2}w}{\partial {x}^{2}}}\right)}_{}^{2}\right]dx$ ${\delta V}_{stiff}={\frac {1}{2}}\sum _{i=1}^{{n}_{r}}{A}_{stiff}\int _{0}^{a}{\sigma }_{stiff}{\left({\frac {\partial w}{\partial x}}\right)}_{}^{2}dx$

(13)

The problem of this idealization is that an additional correction is needed on the flexural inertia respect the panel mid-line following [4]

${I}_{le}={I}_{o}+{\frac {1}{2}}\left[{\frac {{A}_{0}\cdot {\left({e}_{0}+{t}_{s}-{\frac {{t}_{1}}{2}}\right)}^{2}}{1+{Z}_{n}{\frac {{A}_{0}}{{b}_{1}{t}_{1}}}}}+\right.$$\left.{\frac {{A}_{0}\cdot {\left({e}_{0}+{t}_{s}-{\frac {{t}_{2}}{2}}\right)}^{2}}{1+{Z}_{n}{\frac {{A}_{0}}{{b}_{2}{t}_{2}}}}}\right]$

(14)

Where Z_n is a factor depending on the number of waves (figure 4):

Figure 4. Correcting factor for bending inertia

Bending inertia is one of the main parameters to guarantee that the buckling waves are contained in the panel in case of buckling and no catastrophic global buckling failure. Therefore, it is important to correctly capture its idealization. A formulation was developed to consider the correct position of the center of gravity without the need of correcting factors.

3.3. Curved panels

Many structures are not flat, typically an aircraft fuselage structure is a curved stiffened panel. Curvature effect increases the critical buckling load so an equivalent flat panel analysis will provide conservative results. Nevertheless, this approach leads to a weight increase as it is shown in figure 5. In this figure is shown that a 20% critical buckling onset increase can be expected considering the curvature of the panel for a typical fuselage of 2000 mm radius submitted to pure compression. The effect grows exponentially when the curvature increases.

Therefore, the problem should be adapted to take into account the curvature effect. This increases the complexity of the problem adding two new set of unknowns to the problem. Typically for flat panels it was enough to know the out of plane displacements (w(x,y)). Within a curve panel, the axial (u) and transverse displacements (v) enter now into the problem

definition.

Figure 5. Buckling onset load depending on panel curvature radius

The strain at any layer of the laminate (as for flat panels) can be expressed as the addition of the mid-plane strain plus the curvature contribution:

${\left\{\epsilon \right\}}_{k}=\left\{{\epsilon }^{o}\right\}+{\underline {z}}_{k}\left\{\kappa \right\}$

(15)

Where:

${\textstyle \left\{{\begin{matrix}{\epsilon }_{x}^{o}\\{\epsilon }_{y}^{o}\\{\gamma }_{xy}^{o}\end{matrix}}\right\}=}$$\left\{{\begin{matrix}{\frac {\partial u}{\partial x}}\\{\frac {\partial v}{\partial y}}+{\frac {w}{R}}\\{\frac {\partial u}{\partial y}}+{\frac {\partial v}{\partial x}}\end{matrix}}\right\}$ ${\textstyle \left\{{\begin{matrix}{\kappa }_{x}\\{\kappa }_{y}\\{\kappa }_{xy}\end{matrix}}\right\}=}$$\left\{{\begin{matrix}-{\frac {{\partial }^{2}w}{\partial {x}^{2}}}\\-{\frac {{\partial }^{2}w}{\partial {y}^{2}}}-f\left(R\right)\\-2{\frac {{\partial }^{2}w}{\partial x\partial y}}-f(R)\end{matrix}}\right\}$

(16)

In comparison with flat panels, the existence of panel curvature couples the mid-plane strains {ε}^o with the curvatures {κ}, or in other words, the in-plane displacements (u, v) with the out-of-plane displacements w. This coupling is the source of all the differences with the flat panel formulation and forces having to retain the three displacements (u, v, w) throughout the whole calculation process.

Additionally, the displacement functions cannot be only defined for w, but also for u and v. The full development must be done considering now these new equations and the number of new unknowns.

3.4. Transverse shear

The transverse shear stiffness is a resin dominated property in composite panels. This means that this transverse stiffness is an order of magnitude lower than in-plane properties, leading to the inapplicability of Kirchhoff’s hypothesis in bending for thick composite plates. In figure 6 it is shown how transverse shear stiffness influence grows when the t/b ratio increases.

Figure 6 Transverse shear effect on buckling load as function of thickness over panel width

Kirchhoff’s assumption establishes that a cross section of a plate submitted to bending remains on a plane surface normal to the neutral axis. This could be applicable to isotropic material but it leads to non-conservative results in composite, as it is shown in figure 6. To solve it Mindlin’s plate theory [5] is applied, which considers that plane sections remain plane in bending but non perpendicular to the neutral axis. This hypothesis introduces a rotation due to transverse shear to be added to the usual rotation due to classical plate bending theory. This statement implies that the transverse strains ɣ_xz and ɣ_yz are independent of z.

In reality, the transverse shear stresses are not constant in each ply of the laminate, therefore shear correction factors k₁ and k₂ are used in deriving the transverse shear constitutive relations for the plate. Following the procedures of [5] and [6] ɣ_xz and ɣ_yz are replaced by k₁·ɣ_xz and k₂·ɣ_yz

respectively in the strain energy density expression of an orthotropic ply due to transverse shear:

$\left\{{\frac {{Q}_{y}}{{Q}_{x}}}\right\}=\left[{\begin{matrix}{k}_{2}^{2}{A}_{44}&{k}_{1}{k}_{2}{A}_{45}\\{k}_{1}{k}_{2}{A}_{45}&{k}_{1}^{2}{A}_{55}\end{matrix}}\right]\left\{{\frac {{\gamma }_{yz}}{{\gamma }_{xz}}}\right\}$

(20)

where [Q]k is stiffness matrix of lamina k in local axes. The value of k1 , k2 is determined such that the coefficient of Qx in the two shear deformation of strain energy expressions is identical, but it will not be developed here, see reference [11].

This is adding again, 2 new unknowns to the problem that can be approximated by the same shape functions described previously. Therefore 5 unknowns are now needed to be solved following procedure indicated in section 3.

3.5. Validation

Once the development is completed an extensive and intensive validation is done in three different bases:

• Comparison with current close form solutions and widely extensive design curves in literature and standards.
• Comparison with finite element methods
• Comparison with real test data

Obviously non exhaustive validation can be presented in this brief paper but different illustrative examples are enclosed in the following figures, from figure 7 to figure 10, given the adequate confidence that the problem solution is well captured and therefore applicable for sizing and certification:

Figure 7 Solution [8] for combine buckling (z=f(width, thickness)) reproduced by Airbus method

Figure 8. Solution [9] for curved buckling plates repeated through Airbus method

Figure 9 Stringer sections under longitudinal load (T, I and C) and pure shear (Omega) from [10]

Figure 10 Stiffened panel test results vs Airbus (ISAMI) tool results

3.6. Implementation

The mathematical development of the solutions presented here is really complex and its implementation is not straight forward for an efficient and optimum CPU management (see figure 11 for examples of different performances on computational time vs operations)

For example, a cubic running time function will be the result of a naive multiplication of two nxn matrices. Matrix chain multiplication with brute force can be time exponential [5].

If A is a 10 × 30 matrix, B is a 30 × 5 matrix, and C is a 5 × 60 matrix, then:

• computing (AB)C needs (10×30×5) + (10×5×60) = 1500 + 3000 = 4500 operations, while
• computing A(BC) needs (30×5×60) + (10×30×60) = 9000 + 18000 = 27000 operations.

A quicker solution to this problem can be achieved by breaking up the problem into a set of related sub-problems, something that will be extensively used in the energy methods if the adequate performances would like to be achieved.

Figure 11. Number of operations vs time

High complexity problems bring to face quickly the limitation of the technology and therefore adequate technologies need to be understood and managed. Software development skills are important. A poor implementation could ruin a master math development.

4. Conclusions and future work

During the present paper, the different evolutions of energy methods inside Airbus has been presented. These solutions have the advantage of a high performance when compared with more general FE methods (~10 times faster), which makes them suitable to be used, for example, together with optimization engines. Moreover, they are suitable to be linked to sizing processes that chain different failure modes that can interact between them making prediction reliable, fast and safe as described in the introduction

Nevertheless, not all the effects and geometries are today covered by these solutions and therefore finite element models are used when applicability is superseded. Notwithstanding, the energy methods can be evolved to cover:

• Non-rectangular panels.
• Holes of different shapes and sizes.
• Evolutive and/or tapered sections
• Cylindrical, spherical (or domes) and conical shells
• Double curvature shells
• Global buckling of isogrid shells
• Post-buckling behavior

The development of these solutions, mainly the post-buckling behavior, will make a step forward on the buckling treatment reducing the need of finite element models and making available to the engineers a wide variety of solutions to support them in the design phase definitions of new composite structures.

Bibliografía

[1] Hennessy and Patterson, « Morgan-Grampian Publishers », Computer Architecture, 2019

[2] Handbook of structural stability, Buckling of flat panels, NACA-TN-3781, July 1957

[3] J. N. Reddy, Energy Principles and Variational Methods in Applied Mechanics, Wiley, 2nd edition, 2002.

[4] 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, 1953

[5] Mindlin, M. D, Journal of Applied Mechanics, Influence of Rotatory Inertia and Shear on Flexural Motions of Isotropic, Elastic Plates p.31-38. , 1951

[6] Whitney, J. M., Pagano, N. J., Journal of Applied Mechanics, p.1031-1036. Shear Deformation in Heterogeneous Anisotropic Plates, 1970

[7] https://en.wikipedia.org/wiki/Matrix_chain_multiplication

[8] Critical combinations of Shear and Direct Axial Stress for Curved Rectangular Panels, NACA TN 1928, 1949

[9] Flat panels in shear. Buckling of long panels with transverse stiffeners, ESDU 02.03.02, 1983

[10] Martin-Esteban, José Antonio. Energy method for buckling of CFRP interconnected plates with arbitrary boundary conditions, 20th International Conference on Composite Materials, July 2015

[11] Carrasco-Fernández, José. Anisotropic Stiffened Panel Buckling and Bending Analyses using Rayleigh-Ritz Method. International Conference on Advanced Computational Engineering and Experimenting. ACE-X 2011