Nonlinear Static Stability of Imperfect Bio-Inspired Helicoidal Composite Beams

: The objective of this manuscript is to develop, for the ﬁrst time, a mathematical model for the prediction of buckling, postbuckling, and nonlinear bending of imperfect bio-inspired helicoidal composite beams with nonlinear rotation angle. The equilibrium nonlinear integrodifferential equations of imperfect (curved) helicoidal composite beams are derived from the Euler–Bernoulli kinematic assumption. The differential integral quadrature method (DIQM) and Newton-iterative method are employed to evaluate the response of imperfect helicoidal composite beams. Following the validation of the proposed model, numerical studies are performed to quantify the effect of rotation angle, imperfection amplitude, and foundation stiffness on postbuckling and bending behaviors of helicoidal composite beams. The perfect beam buckles through a pitchfork bifurcation. However, the imperfect beam snaps through the buckling type. The critical buckling load increases with the increasing value of elastic foundation constants. However, the nonlinear foundation constant has no effect in the case of perfect beams. The present model can be exploited in the analysis of bio-inspired structure, which has a failure similar to a metal and low interlaminar shear stress, and is used extensively in numerous engineering applications.


Introduction
The bodies of animals and plants are naturally designed in helicoidal arrangement to resist and protect them from enemies. This arrangement provides their bodies with excellent mechanical properties, which inspires researchers and scientists to design a manmade composite structure in the same way. Helicoidal structure, known as "Bouligand structure", is one of the exceptional and predominant arrangements noted in exoskeletons of the arthropod, crustaceans, sapidus, and insects. The helicoidal structures are portrayed by gradually changing the rotation angle of each lamina in the bulk unit [1]. In the case of helicoidal arrangement, the discontinuity of in-plane shear properties gradually decreases. Therefore, the debonding resistance, toughness, strength, and damage tolerance can be improved and designed [2]. Liu et al. [3] illustrated that the delamination resistance increases by decreasing the inter-ply angles in helicoidal laminates. Therefore, the helicoidal composite laminated (COL) structure has been exploited for numerous anti-impact applications, such as tanks, warcraft, warship, and blades of turbine [4]. Moreover, they may be used as alternatives for classical orthotropic laminated structures, which are used in many mechanical, military, civil, aerospace, and aeronautics industries [5]. Cheng et al. [1] and Shang et al. [6] examined the response of bio-inspired helicoidal composite beams with different orientation angles under the transverse point load. Grunenfelder et al. [7] and Jiang et al. [4] experimentally analyzed the impact resistance of helicoidal composite panels. Ginzburg et al. [8] explored the damage tolerance of helicoidal composite panels due to the low velocity impact and proved that they have the ability to sustain a transverse load up to 73%, which is more than the cross-ply scheme. Yang et al. [9] experimentally and numerically studied the multiscale finite element model, low-velocity impact response, and energy absorption capacity of bio-inspired CFRP laminates. Golewski [10] developed a special construction and material conditions to decrease the drawbacks of vibrations on concrete structures. Gul and Aydogdu [11] studied mechanical behaviors of Timoshenko nanobeams using the doublet mechanics theory. Yang et al. [12] and Yin et al. [13] presented a theoretical analysis to investigate the crack-driving force and toughening mechanism of bio-inspired helical structures. Lyratzakis et al. [14] reduced the induced vibrations due to the high-speed train on adjacent reinforced concrete buildings using single or double expanded polystyrene geofoam-filled trenches.
The geometries of beam, plate, and shell are not completely flat, but have small imperfections and curvatures [15]. Lacarbonara et al. [16] explored the two modes that are activated around a veering of frequencies, which result from the nonlinear stretching of imperfect beams. Emam [17], Emam and Nayfeh [18], Gupta et al. [19], and Gunda et al. [20] derived analytical solutions for postbuckling and vibration behaviors of COL beams with and without imperfections. Li and Qiao [21] studied nonlinear pre-and postbuckling behaviors of imperfect anisotropic composite beams by employing the Newton-iterative method and Galerkin method. Emam and Eltaher [22] introduced a mathematical model to predict the hygrothermal influence on buckling and postbuckling behaviors of composite imperfect beams. Emam et al. [23] studied postbuckling and vibration behaviors of imperfect multilayer nonlocal nanobeams under the pre-stress compressive load. Mohamed et al. [24] established a numerical model to analyze the nonlinear free and forced steady state vibrations of imperfect beams using DIQM. Guo et al. [25] exploited a domain decomposition to investigate the static and dynamic response of COL curved beams with elastic boundary conditions. Eltaher et al. [15] presented the influence of curved periodic and nonperiodic profiles of beams on the buckling, postbuckling, and dynamic response using DIQM. Wang et al. [26] and Yang et al. [27] introduced a mathematical model to predict buckling and postbuckling behaviors of composite beams reinforced with graphene platelets under electrical force. Mohamed et al. [28,29] exploited DIQM and the energy-equivalent method to study buckling and postbuckling behaviors of higher order carbon nanotubes. Eltaher et al. [30,31] studied nonlinear buckling, postbuckling, and vibration behaviors of imperfect carbon nanotubes (CNTs) modeled as a beam structure using analytical and numerical methods. Song et al. [32] presented a comparison between the matched-interface-boundary method and its interpolation for vibrations of stepped structures. Emam and Lacarbonara [33] developed a general formulation to investigate nonlinear buckling and postbuckling behaviors of curved beams, including the extensibility and shear deformability of the structure. Boutahar et al. [34] studied bending vibratory behaviors of FG beams using the refined beam theory. Zhang et al. [35] exploited a new method of curvature constrained interpolation in static and buckling behaviors of straight and strong curved thin beams. Almitani et al. [36] developed exact solutions for nonlinear bending, buckling, and postbuckling behaviors of imperfect helicoidal composite beams. Karamanli and Vo [37] studied vibrations of curved zigzag nanobeams using the doublet mechanics theory and finite element procedure.
Although many bio-inspired structures under an impact load have been studied experimentally, very limited research has been performed to present the mechanical bending, buckling, and postbuckling response of curved bio-inspired composite beams. Therefore, this paper aims to comprehensively cover this topic. The rest of the article is organized as follows. Section 2 presents the constitutive equations, material distribution, geometrical configuration, and problem formulation. Section 3 illustrates the numerical solution proce-dure for static bending and buckling response using the differential integral quadrature method. Section 4 discusses the numerical studies and parametric analyses, while Section 5 introduces the main conclusions and novel points derived from the parametric studies.

Problem Formulation
A laminated composite beam with N L uniform layers, total thickness h, length L, and width b are presented in the analysis. With respect to the classical beam theory, the axial (U) and lateral displacements (W) of a generic point at (x, 0,ẑ) can be written as follows [24]: whereû andŵ are the axial and lateral mid-plane displacements, andŵ 0 is the initial rise. The normal strain due to deformation is given by where ε 0 is the normal strain and κ 0 is the curvature of the mid-plane, which are defined as follows [29,38]: The force (N) and moment (M) resultants can be defined as follows [17]: The laminated axial, coupling, and bending stiffness are A ij , B ij , and D ij , respectively, which can be described by where the reduced-transformed stiffness of a single orthotropic lamina Q ij is defined as follows [17]: The transformed matrix T is defined as follows [39]: where θ is the angle of fibers at kth lamina. The plane reduced stiffness Q ij can be evaluated by [17] Q 11 = E 1 1 − ν 12 ν 21 , Q 12 = ν 12 E 2 1 − ν 12 ν 21 , where E 1 , E 2 , ν 12 , and G 12 are four independent material constants. The equations of motion of helicoidal COL beams can be represented by [15] m ∂ 2û ∂t 2 +μ 0 ∂û ∂t where m is the mass per unit length,F is the axial load,F w is the transverse load,μ 0 andμ 1 are the axial and transverse damping coefficients, respectively. As the in-plane inertia and damping are insignificant on the transverse deflection, they can be neglected. Equations (12) and (13) can be reduced into one equation as follows: whereP is the axial imposed force,k s is the elastic shear stiffness of the foundation,k L and k NL are the linear and nonlinear elastic foundation constants, respectively. Moreover,q and F are the distributed transverse and axial loads along the beam length, and Ω is the forced frequency defining the following quantities: The nondimensional governing equation is as follows: ..
where α = , Ω =Ω Herein, the buckling and bending problems are studied. The static equilibrium equation can be acquired by dropping the time dependent terms in Equation (16). The result is as follows [24]: During fabrication or due to heating and cooling processes, the structure may exhibit an initial curved shape as a form of imperfection. Therefore, the initial imperfection can supposedly have the following form, which is accompanied by the first buckling mode as follows: where A 0 is the amplitude of initial curvature. The boundary conditions in dimensionless form for S-S and C-C beams, respectively are as follows:

Numerical Technique
To define the mesh grid points, the shifted Chebyshev-Gauss-Lobatto grid points are used as follows [15]: where N is the number of grid points. According to the DIQM, the first-order derivative of a continuous function y(x) is as follows: Weighting coefficients can be evaluated as follows [40]: From Equation (24), the first-order derivative of a function can be written in a matrix form as follows: where The higher order matrices can be obtained using the matrix multiplication as follows [24]: The definite integral of a continuous function y(x) over the domain can be obtained as follows: where K is the pseudo-inverse of matrix C (1) and the row vector S = [S 1 , S 2 , . . . S N ]. Additional explanations regarding the DIQM and its conversions are comprehensively presented by Equations (19) and (30).

Buckling Problem
To compute the critical buckling load and postbuckling configuration, the external transverse load q should be set to zero. Therefore, one obtains Equation (31) can be written as follows: As known, Equation (32) is a nonlinear nonhomogeneous fourth-order ordinary differential equation, whose exact solution is not available. Therefore, DIQM is employed to discretize Equation (33) as follows: N1 is a row matrix whose elements are the difference between the last row and the first row of matrix C (1) , • stands for the matrix Hadamard product, and the column vector w is defined as follows: where w i = w(x i ). The initial shape of imperfection w 0 (x) is discretized as the known vector as . The corresponding boundary conditions in Equation (21) can be discretized in the same manner. The system of nonlinear algebraic Equation (36) can be written as follows: Notably, the discretized boundary conditions are appropriately substituted in the nonlinear algebraic system of Equation (37). The critical buckling load and postbuckling configuration are computed numerically using the Newton-iterative method. The solution of the linearized form of Equation (34) is used as the initial value for Newton's method.

Bending Problem
For the static bending problem, the external axial load P is set to zero and the following equilibrium equation is obtained: The applied transverse load q(x) can be expressed as follows: iformly distributed load (39) where q 0 is the intensity of the load, δ(.) is the Dirac delta function, and x p is the application position of the point load.
It is difficult to analytically solve the nonlinear equation described in Equation (38). Therefore, the numerical DIQM is used to solve it.
In the case of point load, Equation (38) is as follows: It is well-known that the Dirac delta function has the following properties: To use the properties of the Dirac delta function, Equation (40) is integrated over the domain [41] as follows: The left-hand side of Equation (43) is integrated numerically using the integral operator defined in Equation (30), as follows: Using Equation (41), Equation (44) can be simplified to the following: Discretizing Equation (45) by DIQM, the results in the matrix form are as follows: Mathematics 2022, 10, 1084 8 of 20 Using the uniform load rather than the point load in Equation (40), the function F q can be computed as follows:

Numerical Results
In this section, the numerical results for buckling and nonlinear bending behaviors of perfect and imperfect helicoidal composite beams embedded on elastic foundations are investigated. First, a comparison is carried out to show the validity of the present model. Second, parametric studies are presented to illustrate the significance of helicoidal rotation angle, imperfection amplitude, and elastic foundation constants on buckling, postbuckling, and bending behaviors of helicoidal composite beams.

Validation
To verify the accuracy of the present method, the numerical results for nonlinear buckling and postbuckling configurations of composite beams obtained by the present model are compared with those in the literature. Herein, the critical buckling load and postbuckling configuration of perfect and imperfect composite beams are compared with the analytical solutions presented by Emam and Nayfeh [18] and Emam [17]. For comparison purposes, the beam with six layers has the following bulk and material properties: The following layups have been considered [18] In Table 1, the first three buckling loads of S-S and C-C perfect composite beams with different layups are tabulated and compared with those reported [18]. Moreover, Table 2 presents a comparison of the first buckling load of C-C imperfect composite beams for different laminates. The load-deflection response associated with postbuckling behaviors of C-C perfect and imperfect composite beams for different layups are compared with those obtained by [17], as shown in Figure 1. The excellent agreement with Refs. [17,18] can be observed, which confirm the validity of the present beam model and solution methodology.

Parametric Studies
In this study, the material properties of Jiang et al. [4] are used, as presented in Table 3. Here, it is assumed that L = 100h. Table 4 presents the specifications of the selected layup configurations. Unidirectional and quasi-isotropic laminates are used as references. Table 3. Material and bulk properties.

Material Properties
Bulk Properties    Table 5 tabulates the critical buckling load of S-S and C-C perfect and imperfect composite beams with different layup specifications and various values of imperfection amplitude. As concluded in the case of imperfection amplitude A 0 ≤ 3, the smallest buckling load is observed in the case of HE3 configuration, and the highest buckling load is observed in the case of UD configuration, which has the largest buckling stiffness due to the ability of UD to resist an axial load that is aligned with its orientation (i.e., A and D stiffness have the largest values in the case of UD and smallest values in the case of HE3). The lamination schemes can be arranged in descending order, according to the buckling stiffness, as follows: UD, HR1, HE1, HR2, HR3, HE2, QI, and HE3, respectively. In the case of A 0 = 4, the optimum configuration that can sustain the largest buckling load is HR3. However, the UD has a null load solution and QI has a negative buckling load. Table 5. Critical buckling load (N ) of perfect and imperfect composite beams for different layups (k s = k L = k NL = 0, L = 100h). Variations of the critical buckling load with the dimensionless imperfection amplitude at different layup specifications are presented in Figure 2. As compared with the UD layup, the HR and HE layups with small rotation angles can improve the critical buckling load. Variations of the critical buckling load with the dimensionless imperfection amplitude at different layup specifications are presented in Figure 2. As compared with the UD layup, the HR and HE layups with small rotation angles can improve the critical buckling load.  Figure 3 displays the influence of elastic foundation constants on the critical buckling load of S-S and C-C perfect and imperfect helicoidal composite beams with layup sequences HR1 and HE1. It can be observed that an increase in both shear and linear foundation constants leads to an increase in the buckling load of perfect and imperfect helicoidal composite beams. On the other hand, an increase in the nonlinear foundation constants leads to an increase in the buckling load of imperfect beams and has no effect on the buckling load of perfect beams. Furthermore, it can be interpreted that the effect of shear parameter on the buckling load is more prominent than the effect of linear and nonlinear foundation constants.  Figure 3 displays the influence of elastic foundation constants on the critical buckling load of S-S and C-C perfect and imperfect helicoidal composite beams with layup sequences HR1 and HE1. It can be observed that an increase in both shear and linear foundation constants leads to an increase in the buckling load of perfect and imperfect helicoidal composite beams. On the other hand, an increase in the nonlinear foundation constants leads to an increase in the buckling load of imperfect beams and has no effect on the buckling load of perfect beams. Furthermore, it can be interpreted that the effect of shear parameter on the buckling load is more prominent than the effect of linear and nonlinear foundation constants.  Figure 3 displays the influence of elastic foundation constants on the critical buckling load of S-S and C-C perfect and imperfect helicoidal composite beams with layup sequences HR1 and HE1. It can be observed that an increase in both shear and linear foundation constants leads to an increase in the buckling load of perfect and imperfect helicoidal composite beams. On the other hand, an increase in the nonlinear foundation constants leads to an increase in the buckling load of imperfect beams and has no effect on the buckling load of perfect beams. Furthermore, it can be interpreted that the effect of shear parameter on the buckling load is more prominent than the effect of linear and nonlinear foundation constants.    Figure 4 demonstrates the nonlinear response in prebuckling and postbuckling states of helicoidal composite beams with layup sequence HR1 for various values of imperfection amplitude. The solid lines indicate stable paths, while the dotted lines indicate unstable ones. As opposed to the perfect beams, which show a pitchfork bifurcation, the buckling of imperfect beams exhibits a snap-through behavior. For the perfect beam, a zero equilibrium path is observed in the prebuckling region (P < P cr ). Herein, P cr is the bifurcation point. At critical buckling load P cr , the beam loses its stability, and two symmetrical nonzero stable paths emerge. In snap-through buckling, at each value of the applied axial load P, the composite beam has a nonzero solution. In addition to the primary solution, two secondary solutions appear at P = P cr . Herein, P cr represents the turning point connecting the secondary solutions.  In Figure 5, the load-deflection curves associated with the postbuckling response of S-S and C-C perfect helicoidal composite beams are displayed. It can be observed that the rotation angle of helicoidal composite beams has a great influence on the postbuckling response. For smaller values of axial load, the maximum deflection of HE1 and HE2 is smaller than the maximum deflection of HR1 and HR3, respectively. However, as the axial load increases, this trend is reversed. In Figure 5, the load-deflection curves associated with the postbuckling response of S-S and C-C perfect helicoidal composite beams are displayed. It can be observed that the rotation angle of helicoidal composite beams has a great influence on the postbuckling response. For smaller values of axial load, the maximum deflection of HE1 and HE2 is smaller than the maximum deflection of HR1 and HR3, respectively. However, as the axial load increases, this trend is reversed. Figures 6 and 7 present the load-deflection curves of imperfect helicoidal composite beams in the prebuckling and postbuckling domains of S-S and C-C boundary conditions, respectively. It can be observed that the rotation angle of helicoidal composite beams has a great influence on the buckling response.

UD
S-S supported beam; and (b) C-C Supported beam In Figure 5, the load-deflection curves associated with the postbuckling response of S-S and C-C perfect helicoidal composite beams are displayed. It can be observed that the rotation angle of helicoidal composite beams has a great influence on the postbuckling response. For smaller values of axial load, the maximum deflection of HE1 and HE2 is smaller than the maximum deflection of HR1 and HR3, respectively. However, as the axial load increases, this trend is reversed.

Bending Analysis
Herein, the dimensionless load parameter 0 is defined as follows: for uniform load and (49)

Bending Analysis
Herein, the dimensionless load parameter q 0 is defined as follows: for uniform load and (49) for point load Tables 6 and 7 present the value of imposed force amplitude in the case of uniform and point loads that provide the maximum definite deflection at different layup configurations without elastic foundation effects. It can be observed that to impose the definite deflection (w max = 0.1 mm) for all of the configurations, the UD and HR1 need higher loads than the other configurations. This indicates that the rigidity of UD and HR1 is higher as compared with the other configurations for both C-C and S-S supports. However, HE3 and HR1 are more flexible than the other configurations, which indicates that a small load can be applied to deflect the beam by 0.1 mm. For the same configuration and imposed deflection, the rigidity of C-C beam is higher than the S-S beam by around five times. Table 6. Uniform loadq 0 (N/m) applied to perfect composite beams for different layups and maximum deflection value (k L = k s = k NL = 0 ). In Figures 8 and 9, the load-deflection curves associated with the nonlinear bending response of S-S and C-C perfect and imperfect composite beams under the uniform lateral load with different layups are displayed. It can be observed that the rotation angle of helicoidal composite beams and the imperfection amplitude lead to a change in the trend of load-deflection response. In the case of perfect beams, the maximum deflection of HR1 is higher than the maximum deflection of HE1, and the maximum deflection of HR3 is higher than the maximum deflection of HE3. The reverse occurs for imperfect beams. These observations are valid for S-S and C-C boundary conditions. response of S-S and C-C perfect and imperfect composite beams under the uniform lateral load with different layups are displayed. It can be observed that the rotation angle of helicoidal composite beams and the imperfection amplitude lead to a change in the trend of load-deflection response. In the case of perfect beams, the maximum deflection of HR1 is higher than the maximum deflection of HE1, and the maximum deflection of HR3 is higher than the maximum deflection of HE3. The reverse occurs for imperfect beams. These observations are valid for S-S and C-C boundary conditions. The effect of elastic foundations on the nonlinear bending curves of helicoidal composite beams subjected to the point load with layup specification HR1 are studied in Figure 10. It can be observed that at a given applied lateral load, the maximum deflection becomes smaller as the elastic foundation coefficients increase. Comparing the effects of elastic foundations on different boundary conditions shows that the influence of elastic foundation coefficients is more effective on S-S helicoidal composite beams. The effect of elastic foundations on the nonlinear bending curves of helicoidal composite beams subjected to the point load with layup specification HR1 are studied in Figure 10. It can be observed that at a given applied lateral load, the maximum deflection becomes smaller as the elastic foundation coefficients increase. Comparing the effects of elastic foundations on different boundary conditions shows that the influence of elastic foundation coefficients is more effective on S-S helicoidal composite beams.

UD
The effect of elastic foundations on the nonlinear bending curves of helicoidal composite beams subjected to the point load with layup specification HR1 are studied in Figure 10. It can be observed that at a given applied lateral load, the maximum deflection becomes smaller as the elastic foundation coefficients increase. Comparing the effects of elastic foundations on different boundary conditions shows that the influence of elastic foundation coefficients is more effective on S-S helicoidal composite beams.

Concluding Remarks
This study presents a numerical analysis for the buckling, postbuckling, and nonlinear bending response of helicoidal composite perfect and imperfect beams. Herein, S-S and C-C boundary conditions are considered. Verification studies indicate that the DIQM Figure 10. Effect of elastic foundations on the nonlinear bending response of S-S and C-C perfect and imperfect helicoidal composite beams with layup sequence HR1 x p = 0.5, L = 100h . (a) k Nl = k s = 0; (b) k L = k s = 0; (c) k L = k NL = 0.

Concluding Remarks
This study presents a numerical analysis for the buckling, postbuckling, and nonlinear bending response of helicoidal composite perfect and imperfect beams. Herein, S-S and C-C boundary conditions are considered. Verification studies indicate that the DIQM is an accurate method for the analysis of buckling and postbuckling behaviors of imperfect beams. Several numerical results are presented to study the influence of helicoidal rotation angle, amplitude of initial imperfection, and elastic foundation coefficients on nonlinear bending, buckling, and postbuckling behaviors of composite beams. From the numerical results, the main conclusions are summarized as follows: • For large values of imperfection amplitude, HR and HE layups enhanced the critical buckling load.

•
The perfect beam buckles through a pitchfork bifurcation. However, the imperfect beam snaps through the buckling type. • The rotation angle of helicoidal composite beams reversed the trend of postbuckling response.

•
In the case of A 0 = 4, the optimum configuration that can sustain the largest buckling load is HR3. However, the UD has a null load solution and QI has a negative buckling load. • An increase in the nonlinear foundation constant leads to an increase in the buckling load of imperfect beams and has no effect on the buckling load of perfect beams.

•
The rigidity of UD and HR1 is higher than the other configurations for both C-C and S-S supports. However, HE3 and HR1 are more flexible than the other configurations. • For the same configuration and imposed deflection, the rigidity of C-C beam is higher than the S-S beam by around five times. • The proposed model can be used in the design and analysis of aerospace, naval, vehicles, and biomedical structures, which are manufactured from specific, highstrength, composite laminated perfect and curved beams.