Displacement-based formulation of Koiter’s method: application to multi-modal post-buckling ﬁnite element analysis of plates

.

In particular, the multi-modal formulation of Koiter's approach provides a systematic and efficient procedure to assess the nonlinear behavior of the structure in cases where several buckling modes interact, such as in structures highly optimized for buckling [24,25,26] and imperfectionsensitive shell structures [27,28,29,30,31,32].In such designs, small imperfections due to variations in manufacturing parameters can induce different bifurcation paths [33,34], which can be studied by Koiter's perturbation analysis.
Plates may also show clustered buckling modes, usually when high aspect ratios are involved [16].Madeo et al. [16] evaluated variable stiffness plates where 4 modes were required in the multi-modal expansion to obtain a satisfactory approximation of the post-buckling behavior.In these studies, the authors used a bi-linear mixed 4-node plate element (MISS-4 [35,36,37]) with first-order shear deformation theory.
The aim of the present work is to derive a complete ORCID(s): 0000-0001-9711-0991 (S.G.P. Castro); 0000-0001-7511-2910 (E.L. Jansen) displacement-based formulation for the single-and multimodal Koiter's asymptotic analysis.The formulation and notation herein proposed attempts to keep a close correspondence between the theory and the implemented algorithms, and thereby helpful in addressing the issues experienced in the past with the finite element implementation of Koiter's asymptotic method, as highlighted by Casciaro [9].
The following discussion starts revisiting the double expansion of the total potential functional using Koiter's theory, as proposed by Budiansky [38] for a single-mode asymptotic analysis.The formulation is then expanded for the more general case of a multi-modal asymptotic analysis.General-purpose functional derivatives are presented and expressions for these functional derivatives are later obtained using von Kármán non-linear plate kinematics.These expressions can be used with any displacement-based analytical, semi-analytical or numerical approach; and in the present study it is proposed to approximate the displacement fields using the Bogner-Fox-Schmit (BFS) element [39], one of the most accurate rectangular finite elements for predicting the deflection of thin-walled shells, as stated by Zienkiewicz & Taylor [40,p. 153].The conventional BFS achieves third-order interpolation for the deflection using only 4 nodes per element and 6 degrees-of-freedom per node.An enrichment is proposed to the conventional BFS element to achieve third-order interpolation also for the inplane displacements, keeping only 4 nodes per element and increasing to 10 degrees-of-freedom per node.This enrichment of the in-plane displacements enable high accuracy while approximating the second-order displacement fields, required in Koiter's method.The developed formulation and implementation is verified against existing literature for single-and multi-modal asymptotic expansion of metallic and composite plates, under various loading and boundary conditions.Reference results for multi-modal analysis that are not readily available in the literature are reported, including the second-order displacement fields.

Expansion of the Total Potential Energy Functional
Koiter's asymptotic theory was first presented in 1945 [1].Budiansky [38] presents a complete review of Koiter's asymptotic method, proposing a general-purpose formulation based on the total potential energy functional of the system.The present work elaborates on this formulation and expands it to the complete multi-modal Koiter's analysis.Given a total potential energy functional that has the displacements and the a scalar load multiplier as unknowns [ , ].At a pre-buckling static equilibrium: where 0 0 0 is a known solution corresponding to the prebuckling load 0 .The formulation herein presented is compatible with either a linear or a nonlinear pre-buckling state.In Eq. 1, the notation ′ is used instead of to conveniently express the functional variation as a tensor product between the Fréchet derivative ′ and the variation of the vector containing all degrees-of-freedom [38].The first functional variation becomes a first-order tensor multiplying a vector; the second variation becomes a second-order tensor multiplying two vectors, and so forth, as demonstrated in Eq. 2.
Assume that there exists at least one point of equilibrium that intersects [ ( ), ] at a bifurcation point [ , ], such that = 0 0 0 .This bifurcation point is evaluated according to Appendix A. Hence, at the bifurcation point, the following static equilibrium ought to exist: Assuming that is continuous at the vicinity of , a Taylor expansion around the equilibrium point [ , ] can be applied on Eq. 3 to find approximations for new state of equilibrium [ ( ), ], ∀ > .The error of this approximation becomes zero when → .For this first Taylor expansion, a displacement perturbation is assumed such that: with → = 0 0 0 because crosses the bifurcation point.The notation ( ) = ( ) [ , ], is adopted to obtain the result of the first Taylor expansion shown in Eq. 5, that approximates the total potential energy functional at [ , ], with the new displacement state being = + .
It is worth emphasizing that each Frechét derivatives from the first expansion of Eq. 5 is calculated at [ , ], and not [ 0 0 0 , ]; where = 0 0 0 and 0 0 0 is a linear or nonlinear pre-buckling state.Terms ( ) represent ℎ -order Frechét derivatives of the total potential energy functional about the displacement degrees of freedom .Equation 2 demonstrates how these Frechét derivatives of the total potential energy functional ultimately generate ℎ -order tensors.Vectors and are 1 -order tensors, such that ( ) ( −1) should be read as ℎ -order tensor products.Defining the following notation for the derivatives of in terms of [38]: To approximate the total potential energy functional at [ , ] each Frechét derivate of Eq. 5 undergoes a second Taylor expansion that assumes ( ) continuous for any load increment − , resulting in: The expressions in Eq. 7 are readily applied into Eq. 5 to obtain the approximation for the total potential energy functional at the new equilibrium, i.e. ′ [ , ], which becomes:

Single-Mode Asymptotic Analysis
Koiter [1] proposed to express − and − using the following asymptotic expansion: where: (1) is a scalar parameter.
(2) is a first-order field, taken directly from one or a linear combination of multiple linear buckling modes.Vector is customarily rescaled dividing by the maximum normal displacement amplitude and multiplying by the plate or shell thickness.(3) is a second-order field that provides a correction to the first-order field.(4) The third-order field , and higher, are assumed to have a negligible contribution.(5) and are respectively first-and second-order load parameters to be determined.Equation 9 is a reduced-order model (ROM) relating the load and displacement around the equilibrium point [ , ].Note that this ROM could have been built around any equilibrium point, a property that is explored in the Koiter-Newton approach [41,42,43,44].Koiter's expansion given by Eq. 9 is directly used into Eq.8 to render: Note in Eq. 10 that only the terms up to 3 are shown.Equation 10 must be satisfied regardless the values of and , such that each term must vanish separately.The arbitrary value = can be used [38], with being the first eigenvector when = 1, or a composition of eigenvectors obtained at the bifurcation point [ , ].With = the orthogonality of the second-order field leads to Hence, Eq. 10 can be used to obtain the equations for and in the single-mode expansion: Equation 11 is in agreement with Casciaro [9].Note in Eq. 11 that can be calculated using only , whereas additionally needs the second-order field .Budiansky [38] suggests to compute using the terms for 2 in Eq. 10: which must hold for all arbitrary variations , such that: For problems in solid mechanics ′′ will generally be an invertible positive-definite square matrix [38], such that the solution of Eq. 12 can be: Note in Eqs. 12 and 13 that vector ̄ ̄ ̄ is appropriately used instead of because Eq. 13 allows multiple solutions of ̄ ̄ ̄ for different multipliers applied to .However, one must guarantee that any calculated ̄ ̄ ̄ is orthogonal to such that the second-order displacement vector becomes a valid vector basis for the reduced-order model of Eq. 9.
Gram-Schmidt orthogonalization [45] can be used to obtain the orthogonal component of ̄ ̄ ̄ , where first the projection of ̄ ̄ ̄ onto is obtained as: , the orthogonal can be calculated subtracting the projection from ̄ ̄ ̄ , using:

Multi-Modal Asymptotic Analysis
The possibility to perform post-buckling analysis in structures with clustered buckling modes is one of the apparent advantages of Koiter's method compared to other methods [9,16].The necessity of considering multiple modes in the asymptotic expansion has been demonstrated by many authors, e.g. when the vibration frequencies are closely spaced [46]; variable-angle tow plates with clustered buckling modes [16]; and imperfection sensitive structures [47,38,48,28,49,29,50,22]. The single-mode asymptotic expansion of Eq. 9 can be generalized to a multi-modal asymptotic expansion, as shown in Eqs. 15 and 16 [51,12]: where summation convention is applied for repeated indices , , = 1, 2, ⋯ is applied.Equation 16 is a reduced-order model consisting of a system of equations obtained for = 1, 2, ⋯ , , which has 1 , 2 , ⋯ , unknowns.The value correspond to the ℎ linear buckling eigenvalue , always re-scaled by dividing with the maximum out-of-plane displacement and multiplying by the plate thickness.Finding the right number of linear buckling modes in the multi-modal analysis is a common question [9], and an accepted criterion is to select a number of modes that lies within 10%-20% departing from the first critical load [9].
Note that Eq. 15 consists on a reduced-order model to calculate displacements based on a pre-buckled state with known linear buckling modes and known second-order displacement fields .As in the case of the single-mode expansion, for plates and shells it is customary to re-scale dividing by the maximum normal displacement amplitude of and multiplying by the plate or shell thickness.The coefficients for = 1, … , are found for each load after solving the system of equations given by Eq. 16.Solving Eq. 16 requires the calculation of all coefficients and .The expressions given by Eqs. 15 and 16 are applied to the expanded total potential energy functional of Eq. 8.The terms multiplying and are collected, analogously to the terms multiplying 2 and 3 for the single-mode expansion in Eq. 10.The collected terms for the multimodal expansion are shown in Eq. 17, where the following orthogonality property of the linear buckling modes is used: ⟨ , ⟩ = 0, ∀ ≠ ; leading to ′′ = 0, ∀ ≠ ; • ′′ = 0, ∀ ≠ ; and For the expanded equilibrium to be stationary, each term in Eq. 17 must vanish separately, similarly to the singlemode expansion.Assuming = in Eq. 17, the expressions for and can be obtained, as respectively given in Eqs.18 and 19.
Note in the calculation of that the second-order fields are needed.As in the single-mode expansion case, a nonorthogonal second-order field ̄ ̄ ̄ is calculated first.Note that the terms in brackets multiplying in Eq. 17 are obtained for any ℎ mode.Therefore, the contribution for all = 1, … , modes are added and the following equation for ̄ ̄ ̄ is obtained: The orthogonal second-order field vectors in the multimodal asymptotic expansion can be obtained after successive Gram-Schmidt orthogonalization [45] operations, used to remove the components of ̄ ̄ ̄ that are parallel to all linear buckling modes used in the multi-modal expansion , with = 1, 2, … , .This orthogonalization procedure is formulated in Eq. 21.

Functional Derivatives using Von Kármán Kinematics
This section demonstrates how to calculate the functional derivatives leading to the ℎ -order tensors ( ) , • ( ) and

Strains
Von Kármán proposed a kinematic relation for plates that neglect various non-linear terms that come from the full nonlinear Green-Lagrange strain-displacement relations [52, section 2.2.2].Von Kármán kinematics are also referred to in the literature as Donnell-type [53,54] or Kirchhoff-Love non-linear equations.Using classical equivalent single layer theory, the three-dimensional strains are expressed as ( , , ) = ( , )+ ( , ), such that the extensional and rotational strains can be defined in terms of the in-plane displacement field variables ( , ), ( , ), ( , ) as: Assuming the following approximation for the displacement field: where , , are known shape functions; each component , , can be written adopting summation convention for repeated indices as: with = 1, 2, ⋯ , .The strain variations can be represented as: where the ′ (prime) symbol is used to denote a Frechét's differentiation.Adopting index notation to represent the strains, such that 1 = , 2 = , 3 = , 1 = , 2 = , 3 = ; the first and second variations of the extensional and rotational strains become: With these definitions, the first Frechét's differentiation of the strains can be represented as: The second differentiation: Note in Eq. 28 that ′′ represents a symmetric secondorder tensor, which is an important property to be considered while implementing the method.
For the differentiations with respect to , we must recall that all functional expansions were calculated about the bifurcation point [ , ], such that the strains and stresses are those corresponding to the displacement = 0 0 0 , with = .Starting with and , using the notation If the pre-buckling state 0 0 0 is evaluated linearly for a plate with no bending-extension coupling and only in-plane prebuckling loads, 0 = 0 and the initial post-buckling analysis is greatly simplified.Nevertheless, the formulation presented next is valid for the more general case of 0 ≠ 0.
The second differentiation with respect to gives: For ′ and ′ , the first differentiation with respect to gives: For the second differentiations with respect to , •• ′ = 0 0 0 and •• ′ = 0 0 0. For ′′ and ′′ all derivatives with respect to are zero.

Stresses
Based on Eqs.22 -31 it is straightforward to compute the corresponding stresses.Using classical constitutive relations for laminated composites [55] and adopting the index notation: 1 = , 2 = and 3 = ; 1 = , 2 = and 3 = ; the stress-strain relations can be written as: where represents the plate membrane stiffness; the membrane-bending coupling stiffness; and the bending stiffness; all for , = 1, 2, 3.The first Frechét derivative of the stress terms are: Recalling from Eq. 28 that ′′ = 0, the second Frechét derivatives are: Note that ′′ , ′′ are symmetric second-order tensors.The first derivatives with respect to can be readily computed as: • ′′ = 0 • ′′ = 0 Finally, the second derivatives about are:

Functional Derivatives
Assuming a general loading vector with distributed forces ̂ ̂ ̂ at the plate boundaries Ω, the total potential energy can be written as: where Ω = and summation convention is adopted for terms with repeated index with = 1, 2, 3.The stationary total potential energy at [ , ] is defined as ′ = ′ [ , ], calculated using the first variation of = [ , ]: The variation is defined as = = {⋯ , , ⋯} ⊺ , such that the first Frechét derivative of the total potential energy becomes: Resuming with the second Frechét derivative, now replacing with = = {⋯ , , ⋯} ⊺ : Note that ′′ in Eq. 42 represents a second-order tensor with all geometrically non-linear terms present.If a linear pre-buckling state is assumed, ′′ becomes the linear constitutive stiffness matrix of the system [9].Continuing with the third Frechét derivative, now using = = {⋯ , , ⋯} ⊺ : Lastly, using = = {⋯ , , ⋯} ⊺ , the fourth Frechét derivative gives: It is now possible to compute the functional differentiations corresponding to the second expansion of the total potential energy functional, about .From Eq. 42, the first differentiation about becomes: The second differentiation of ′′ about becomes: The first differentiation of ′′′ can be calculated based on Eq. 43 as: The second differentiation of ′′′ about can be calculated from Eq. 47: The first and second differentiation of with respect to can be calculated based on Eq. 44:

Modified Bogner-Fox-Schmit Finite Element
The Bogner-Fox-Schmit (BFS) finite element [39] is a classical C1 contiguous confirming plate element obtained by taking tensor products of cubic Hermite splines.The BFS is one of the most accurate rectangular finite elements for thin-walled shells, as stated by Zienkiewicz & Taylor [40,p. 153], and therefore has been chosen to implement the formulation herein developed.With only 4 nodes per element, the standard BFS element approximates the out-of-plane displacements using 3 -order polynomials, which is still a reasonable low-order interpolation for plates and very simple to implement [56], in contrast with triangular elements which use higher order polynomials [56], such as the Argyris element [57].
A drawback of the standard BFS element when used in the context of Koiter's method is the linear interpolation of the in-place displacements [39], resulting in a poor convergence when trying to represent the second-order displacement fields and all dependent quantities such as the factors of Eq. 19.
A modified BFS element is proposed to enable third-order interpolation also for the in-plane displacements.Additional 4 degrees-of-freedom per node are included, being the first derivatives of the in-plane displacements: , , , , , and , .Hence, no nodes are added and the resulting element has still 4 nodes with 10 degrees-of-freedom per node.This enhanced BFS element is referred to in the following discussion as BFSC: Bogner-Fox-Schmit-Castro.The in-plane , and out-of-plane displacements are approximated using: where contains the 10 degrees-of-freedom of the ℎ node.The matrices containing the shape functions , , are defined as: with , , , calculated using natural coordinates [58,59,60]: where , are respectively the finite element dimensions along , .The values of , given in Eq. 54 were adopted for each of the four nodes.
In the present study, only rectangular elements were used, such that the natural coordinates can be defined simply as: = 2 ∕ − 1; and = 2 ∕ − 1.All derivatives of , , required in the strain equations can then be calculated in terms of the natural coordinates using Eq.55.All integrations over the finite element domains are performed numerically using standard Guass-quadrature and a mesh of 4 × 4 integration points per element, unless otherwise specified.

Implementation
The main challenge to implement Koiter's asymptotic approach is the calculation of the third-and fourth-order tensors ′′′ and .A naive implementation can easily blow the memory of any modern computer.For example, a mesh with only 20 nodes, or 200 degrees of freedom using the BFSC element, would generate a complete tensor with 200 4 = 1.6 × 10 9 elements to be stored, requiring 11.92 GB of memory if double precision is used.Moreover, calculating all 1.6 × 10 9 terms is computationally too expensive.For slightly larger problems found in any engineering application such naive approach becomes just impracticable.Many drawbacks experienced in the past with Koiter's asymptotic method are due to mistakes in implementation more than to intrinsic defects of the asymptotic approach itself, as stated by Casciaro [9].The implementation herein proposed carefully takes into account all non-zero terms derived in Section 5.3 while calculating the higher-order tensors.The strategy adopted consists on evaluating the products involving highorder tensors ′′′ and on-the-fly, already during the numerical integration of each element.For instance, defined in Eq. 44 is never fully calculated.Instead, the tensor product is evaluated and stored during the finite element assembly, using the following rearrangement of Eq. 44: where Ω represents the domain of one finite element.Note that only tensors up to second-order are fully calculated and stored, e.g.= = = , such that is calculated with Eq. 28 and with Eq. 34.The implementation carried out for the present study is based on Python [61], NumPy [62] and Cython [63], where all tensor products are efficiently evaluated using NumPy [62], as illustrated in the following pseudo-code: The authors could verify the benefit of the proposed formulation and notation while implementing the method, achieving a one-to-one correspondence between the method and implemented algorithms.

Results for single-mode Koiter's asymptotic expansion 8.1. Isotropic plates
Lanzo et al. [64] performed initial post-buckling analysis in isotropic plates with various aspect ratios, loading and boundary conditions, as illustrated in Figure 1.Table 1 Table 2 shows the convergence of , often referred to in the literature dealing single-mode asymptotic expansion as the " −factor".The results from Lanzo et al. [64] are compared against the BFSC implementation with progressively refined meshes.The mesh refinement is controlled by the number of elements along , named , calculating the number of elements along using = ∕ × .The numerical integration was done using 4x4 integration points, which provided a converged and stable integration scheme, i.e. we verified that adding more integration points did not change the results and using less integration points (3x3) resulted in unstable behaviour for models B1, B2, C1 and C2.
The first buckling mode 1 1 1 and corresponding secondorder field 11 11 11 for each all models with ∕ = 3 are given in Figures 3 -4.Those results correspond to a mesh of ( = 60) × ( = 20) BFSC elements.The plots are created using a visualization mesh of 900 × 300 points along and , respectively; calculating the consistent displacement according to Eq. 52 at each plotting point.Note the complexity of the second-order displacement, which are in-plane for this case.The fast convergence for shown in Table 2 is a direct consequence of the fast convergence of the second-order field .The enrichment on and enabled by the BFSC element is the main reason for this fast-converging behavior.

Composite plates
A composite plate presented by Phan & Reddy [65] and Bilotta et al. [66] will be used as a reference for verification of the single-mode expansion for composite plates.The following material properties are used: Using the load and boundary conditions of model 1 of Fig. 1, the plate geometry is defined based on = = 1 , and thickness of ℎ = 0.1 .The laminate stacking sequence is 0∕90∕90∕0 , with the principal direction being the plate's − .Table 3 shows the convergence of the critical linear buckling eigenvalue for the BFSC element.The results from Phan & Reddy [65] are exact solutions obtained with classical plate theory, whereas the results from [66] are based on bi-quadratic finite elements and a mesh of 25x25 elements.Note the good agreement between the BFSC element and the literature for all anisotropic ratios 1 ∕ 2 , already with a mesh of only 4x4 elements.
The convergence of the coefficient is shown in Table 4.The second-order displacement field 11 11 11 was solved with Eq. 13 and orthogonalized with Eq. 14.Note in Figure 6

Results for multi-mode Koiter's asymptotic expansion 9.1. Isotropic plates
Tiso [10] presented multi-modal results for model A1 using a geometry of: = 0.14 , = 0.10 and thickness of 1 ; with isotropic material properties = 70 and = 0.3.With this configuration, the first two buckling loads are close and respectively equal to 2828.15 and 2866.50 .Tiso's definition of coefficients can be translated to the notation herein presented as given in Eq. 57, which has a different index ordering when compared to Eq. 19.Despite the different index ordering, Eq. 19 and Eq.57 are equivalent in cases where = 0.
Table 5 show the verification of the formulation herein presented against the results from Tiso [10] for a multimodal expansion using 2 modes.
Jansen & Rahman [67] used the general purpose finite element code DIANA [68] to investigate the dynamic response of plates and shells using reduced-order models based on  Koiter's theory.The plate therein investigated has a geometry of = 0.6 , = 0.2 , total thickness of 1 .The same plate was used to verify the formulation and implementation herein presented for the multi-modal asymptotic expansion of a composite plate.Using the loading and boundary condition of model A1 (Figure 1, the analyses in DIANA are carried out using two element types: 1) The eight-node quadrilateral shell element CQ40L, based on first-order shear deformation kinematics (FSDT), mesh with 16x48 elements; 2) Allman-type triangular element [69,70], with kinematics based on the classical laminated plate the-ory (CLPT).The Allman-type element triangular mesh has 24x72x2 elements.Table 6 presents the linear buckling eigenvalues and Table 7 presents the results for the coefficients using 5 modes, calculated per Eq.19 for a plate with isotropic properties = 70 and = 0.3.The relative error is calculated about the results using DIANA's Allman-type element.The BFSC mesh has 16x48 elements, and the authors verified that all models are converged.There is clearly a close agreement between DIANA's implementation and the present formulation, with a maximum difference of 9%, attributed to differences verified for the linear buck-    6.For many coefficients, there is a sign switch because are symmetric eigenmodes about the plate mid-plane, providing interchangeably positive and negative solution for some of the coefficients.
Figure 7 illustrates the second-order field modes for the isotropic plate used in the multi-modal asymptotic ex-     Again, a maximum error of 9% is verified for the coefficients, attributed to the differences in values between the Allman-type and BFSC elements.

Composite plates
Figure 8 illustrates the second-order field modes for the composite plate used in the multi-modal asymptotic expansion.

Conclusion
A complete displacement-based formulation for the Koiter's multi-modal asymptotic expansion was presented.The formulation and notation adopted enabled a close correspondence between the formulation and implemented algorithms, facilitating further developments on the field of initial post- buckling analysis and reduced-order models based on Koiter's method.The present work used state-of-the-art collaborative tools to implement the developed methods: Python [61], NumPy [62] and Cython [63].
A modified Bogner-Fox-Schmit (BFS) finite element  was proposed and referred to as Bogner-Fox-Schimt-Castro (BFSC), which has enriched in-place displacements when compared to the conventional BFS element to enable an accurate representation of the second-order displacement fields.This element has only 4 nodes per element and 10 degrees-of-freedom per node, which are: , , , , , , , , , , , , , , and , .The results herein presented showed that the third-order interpolation of provided led to an ultra-fast convergence of the linear buckling eigenvalues, whereas the third-order interpolation of and led to a ultra-fast convergence of the second-order displacement fields and consequently of the factors obtained in Koiter's analysis.
Future studies should investigate the same strategy herein adopted to enrich the in-plane interpolation of the BFS element to create modified versions of the TUBA3 family [57].Other approaches for enriching the membrane displacement field could also be investigated, such as the interpolation covers suggested by Jun et al. [71] for triangular shell elements.Another alternative is to use fifth-order Hermite polynomials enhancing the BFS element to also interpolate the rotational strains , and , , keeping 1 continuity even with distorted meshes, thus keeping the high convergence rates herein observed for non-rectangular meshes [72].
In the present study, Von Kármán nonlinear kinematics were used to derive the ℎ -order tensors obtained from the total potential energy of the system.The effect of using more complete nonlinear kinematics such as Sanders [49,52] and Timoshenko & Gere [52] should also be investigated, especially for shells with single and double curvature.
The displacement-based notation herein presented simplifies formulating and implementing the compatibility be- for isotropic plate used in multi-modal asymptotic expansion analysis; columns from left to right are: , , tween different semi-analytical domains [73,74,75,76], even for Koiter's asymptotic analysis.Such possibilities would allow one to apply semi-analytical methods for stiffened panels and shells using efficient displacement representations such as those provided by Legendre hierarchical polynomials [77,78,79,80].Finally, the effect of initial imperfections should be included in the presented formulation to enable asymptotic perturbation analyses, especially important in imperfection sensitive shells.

A. Linear Buckling Analysis
In Eq. 5 it is known from the equilibrium condition that ′ [ , ] = 0. Dividing the remaining terms by ‖ ‖: Equation 61 represents the Neutral Equilibrium Criterion [29], which can be used to compute and , i.e. the linear buckling eigenvalues and eigenmodes, based on a prebuckling state [ 0 0 0 , 0 ].For a system with degrees-offreedom there are pairs , , so that one should select the mode that corresponds to the lowest positive , named in the previous discussions for single-mode asymptotic expansion.In most commonly used eigenvalue solvers, the eigenvalues can be computed and sorted in ascendant order, such that only a few first pairs of eigenvalues and eigenvectors need to be calculated.

Table 1
Convergence of BFSC element for rst buckling mode, isotropic case

Table 3
Convergence of BFSC element for rst buckling mode and model 1 with ∕ = 1, anisotropic case

Table 5 Multi
-modal expansion coecients Tiso for isotropic plate

Table 6
Eigenvalues for isotropic plate used in multi-modal analysis, units in[ ∕ ]

Table 8
presents the linear buckling eigenvalues and Table9present results for a composite plate with layup of +45∕ − 45∕90∕0 and material properties 11 =

Table 8
Eigenvalues for composite plate used in multi-modal analysis, units in[ ∕ ]