Free Vibration Analysis of Curved Laminated Composite Beams with Different Shapes, Lamination Schemes, and Boundary Conditions

A general formulation is considered for the free vibration of curved laminated composite beams (CLCBs) with alterable curvatures and diverse boundary restraints. In accordance with higher-order shear deformation theory (HSDT), an improved variational approach is introduced for the numerical modeling. Besides, the multi-segment partitioning strategy is exploited for the derivation of motion equations, where the CLCBs are separated into several segments. Penalty parameters are considered to handle the arbitrary boundary conditions. The admissible functions of each separated beam segment are expanded in terms of Jacobi polynomials. The solutions are achieved through the variational approach. The proposed methodology can deal with arbitrary boundary restraints in a unified way by conveniently changing correlated parameters without interfering with the solution procedure.


Introduction
The moderately tall curved laminated composite beams (CLCBs) are widely applied in engineering fields. Besides, the laminated composite material has excellent mechanical properties including high compression-resistance capacity, high strength-to-weight and stiffness-to-weight ratios, preeminent corrosion-resistance, and powerful customizable capacity [1,2]. In practical use, the CLCBs may commonly undergo various kinds of complicated dynamic loads and other complex work environments, which lead to excessive vibration and fatigue damage of the structure. Dynamic modelling is the precondition for understanding the vibration characteristics of CLCBs. Thus, this paper aims to evaluate the vibration features of CLCBs with alterable curvatures by a modified variational approach in the framework of higher-order shear deformation theory (HSDT).
Numerous works have been carried out on the vibration problems of CLCBs. A group of equations were constructed by Qatu [3] for vibration analysis of simply supported CLCBs, which

Description of the CLCB Model
The schematic plot and geometric parameters of a CLCB are presented in Figure 1. The curved beam has curvature radius R ϕ , thickness h, and width b. The orthogonal coordinates (ϕ, y, z) are located at the middle surface of the CLCB. The ϕ-, y-, and zcoordinates are along the central line (indicated as red line), width, and thickness orientations, respectively. The curved beam comprises several orthotropic layers, where Z k represents the distance from the top surface of the k'th layer to the middle surface. The fiber orientation angle of the k'th lamina with respect to the ϕ-axis of the CLCBs is defined as α k f iber , and the case of α k f iber = 0 is presented for a better understanding. Vibration of the CLCB occurs in the ϕ-z plane and can be characterized by the central line. In engineering applications, different types of CLCBs with diverse curvatures may be encountered, e.g., hyperbolic, parabolic, elliptic, and circular ones. For each type, their radii of curvatures R ϕ (ϕ) can be described in terms of different geometric parameters. Different types of CLCBs are shown in Figure 2; note that the x-axis, which is perpendicular to the z-axis in the ϕ-z plane, is added for a better description.
For an elliptic curved beam (Figure 2a), R ϕ (ϕ) is defined in terms of semimajor axis length a e and semiminor axis length b e : a 2 e sin 2 ϕ + b 2 e cos 2 ϕ 3 . (1) For a parabolic curved beam (Figure 2b), R ϕ (ϕ) can be expressed as where R 0 and R 1 indicate the horizontal radius.
where R0 and R1 indicate the horizontal radius. The hyperbolic CLCB (Figure 2c) is represented by the following equations: where Rs denotes the distance from ox to o x ′ ′ (revolution axis); ah signifies the length of the semitransverse axis, and bh indicates that of the semiconjugate axis. Figure 2d shows the circular curved beam, where R indicates the mean radius.

Energy Expression of the CLCB
The energy expressions of the CLCB are established in the framework of HSDT. Hence, the displacement components of the CLCB for an arbitrary point can be written as in which u and w indicate the displacements along the φ and z orientations at the middle surface, The hyperbolic CLCB (Figure 2c) is represented by the following equations: where R s denotes the distance from ox to o x (revolution axis); a h signifies the length of the semitransverse axis, and b h indicates that of the semiconjugate axis. Figure 2d shows the circular curved beam, where R indicates the mean radius.

Energy Expression of the CLCB
The energy expressions of the CLCB are established in the framework of HSDT. Hence, the displacement components of the CLCB for an arbitrary point can be written as W(ϕ, y, z, t) = w(ϕ, y, t), in which u and w indicate the displacements along the ϕ and z orientations at the middle surface, respectively. In addition, ψ ϕ denotes the rotation of the transverse normal with respect to the ϕ direction. φ ϕ and λ ϕ signify the higher-order terms connected to the Taylor series. Mathematically, The linear strain-displacement relationships, which consider z/R ϕ terms, can be given as follows: where ε 0 ϕ and γ 0 ϕz signify the middle surface normal and shear strains, respectively, and k i(i=0,1,2) ϕ and k i(i=0,1,2) ϕz indicate curvature changes. They are defined as Materials 2020, 13, 1010 6 of 22 The force (N) and moment resultants (M) of the CLCB can be represented by middle surface strains and curvature changes: where N ϕ , N y , N ϕy , N 1 ϕ , N 1 y , and N 1 ϕy signify the in-plane forces; M ϕ , M y , M ϕy , M 2 ϕ , M 2 y , and M 2 ϕy denote the moment resultants; and Q ϕ , Q y , P ϕ , P y , Q 1 ϕ , Q 1 y , P 2 ϕ , and P 2 y indicate the transverse shear force resultants. The corresponding extensional stiffnesses A ij (i, j = 1, 2, 6), coupling stiffnesses B ij , bending stiffnesses C ij , and other stiffness coefficients (D ij , E ij , F ij , G ij and H ij ) can be expressed as where S k ij indicates transformed reduced stiffness coefficients for the k'th layer. In addition, they are defined as where α k f iber denotes the fiber orientation angle of the k'th lamina with respect to the ϕ-axis of the CLCBs and S k e f (e, f = 1, 2, 4, 5, and 6) is the material properties. With regard to CLCBs, the following parameters (N y , N ϕy , N 1 y , and N 1 ϕy ; Q y , P y , Q 1 y , and P 2 y ; M y , M ϕy , M 2 y , and M 2 ϕy ) are presumed to be zero. Hence, Equation (8) can be rewritten as where Materials 2020, 13, 1010 8 of 22 On the foundation of previous discussions, the kinetic energy of the CLCB may be represented as where (I 0 , I 1 , I 2 , I 3 , I 4 , I 5 , Besides, the strain energy may be expressed as

Variational Formulation for Curved Beam
In general, by solving the governing equations associated with boundary conditions, the eigenvalue solutions can be gained. However, it is not an easy job to conduct this process, as choosing suitable admissible functions for arbitrary boundary conditions is difficult. Alternatively, a new procedure by which the problem is expressed in a modified variational form may be developed to simplify the solution.
In this research, a highly efficient and accurate variational approach is adopted and vibration analysis of moderately tall CLCBs with diverse curvatures is conducted. First, the total variational functional total is obtained, the terms of which include the kinetic energy T, strain energy U, and strain energy function connected with the interface and boundary restraint p f . Then, the displacement field is expanded in terms of Jacobi polynomials, which contain unknown Jacobi expanded coefficients. By conducting the variational operation (i.e., δ total = 0) with respect to unknown Jacobi expanded coefficients, an equation in matrix form can be acquired. Through solving this equation, the frequencies and corresponding mode shapes can be easily determined.
Besides, the formulation is derived by employing the multi-segment partitioning technique. The curved beam is separated into N ϕ identical segments along the ϕ orientation. As each segment becomes a free-free substructure, the displacement discretization with respect to admissible functions becomes more convenient. This is related to the fact that the continuity conditions of the CLCB need not be imposed, as their satisfaction is implemented in a variational statement. As such, the issue is simplified to simulating the interaction of each beam segment with common boundaries. This will make the computer implementation of CLCBs much simpler: (a) The selections of admissible functions become flexible as both boundary and continuity conditions are relaxed; (b) by flexibly choosing the appropriate admissible functions, the convergence performance of the present methodology can be quite fast and stable; and (c) the accuracy of the modified variational method can be substantially improved by segmentation, especially for higher modes. Then, the vibration problem may be distinguished by an improved variational principle [24][25][26], which turns into finding the minimum value of a variational functional as where T i is the maximum kinetic energy of the i'th beam segment and U i is the corresponding maximum strain energy. To construct a uniform model to handle arbitrary boundary conditions, as well as ensure numerical stability, the technique of penalty function is exploited. The strain energy function connected with the boundary elastic restraint can be shown as where k τ (τ = u, w, ϕ, φ, ν) indicates the penalty terms expressing the elastic stiffness at both ends of the curved beam, and η τ (τ = u, w, ϕ, φ, ν) represents the continuity or boundary coefficients. By properly defining η τ and k τ , both continuity conditions for the interface and the arbitrary boundary conditions can be conveniently obtained. For the continuity conditions of two neighbored curved beam domains, the continuity coefficients η τ (τ = u, w, ϕ, φ, ν) = 1, while for boundary conditions, different combinations of boundary coefficients η τ and k τ are imposed to achieve various boundary conditions (Table 1). Meanwhile, it is worth noting that both classical and elastic boundary restraints can be achieved by properly selecting the values of k τ (τ = u, w, ϕ, φ, ν); note that k τ0 and k τ1 represent the penalty terms at edges ϕ = ϕ 0 and ϕ = ϕ 1 , respectively. The elastic boundary restraints represent a boundary condition between simply supported and clamped boundary conditions, which can be modeled by springs at the edges. For instance, by setting one or some of the k τ (τ = u, w, ϕ, φ, ν) at certain values, the elastic boundary conditions can be conveniently obtained. Table 1. Boundary coefficients and penalty parameters for various boundary conditions. More information about the penalty terms to handle the boundary conditions can be found in previous investigations [27][28][29][30][31]. Consequently, the variational form for a CLCB subjected to arbitrary boundary conditions is

Solution Procedure
By introducing the modified variational approach in conjunction with the multi-segment partitioning strategy, the choice of admissible displacements for each beam segment can be flexible. This is due to the fact that continuity and boundary conditions of the curved beam are relaxed by the functional total . The total allows the use of identical displacement functions for every curved beam segment. In the present analysis, the displacement components for each divided segment of the CLCB are congruously represented by means of Jacobi polynomials.
where U m , W m , ψ ϕ m , Φ ϕ m , and Υ ϕ m indicate the relevant Jacobi expanded coefficients; P (α,β) m (ϕ) represents the Jacobi polynomial of order m, which is related to the displacement components along the central line orientation; and ω and t signify angular frequency and time, respectively. Besides, the maximum value of m or the truncation terms are represented by M. By choosing different combinations of Jacobi parameters α and β, various orthogonal polynomials can be obtained such as Chebyshev, Gegenbauer, and Legendre polynomials. More details about these can be discovered in [32][33][34]. Substituting Equations (13) and (15)-(17) into Equation (18) and incorporating Equation (19), the following can be achieved: Conducting the variation operation for total , i.e., δ total = 0, with regard to Jacobi expanded coefficients (U m , V m , ψ ϕ m , Φ ϕ m , and Υ ϕ m ), the following characteristic equations can be obtained: in which K signifies the stiffness matrix, M indicates the mass matrix, and E denotes the undetermined coefficients. Apparently, through the solution of Equation (21), the eigenvalues and eigenvectors can be acquired.

Results and Discussions
In this part, the convergence performance, reliability, efficiency, and accuracy of the current methodology are studied by a number of numerical cases. For the sake of brevity, the free, clamped, simply supported, slided, and elastic boundary conditions are represented by F, C, SS, SD, and E, respectively, as shown in Table 1. Besides, three kinds of elastic boundary restraints (E1, E2, and E3) are taken into consideration. Then, a two-letter string is utilized to represent the boundary conditions of two ends, e.g., SS-E indicates the curved beam subjected to the simply supported boundary condition at edge ϕ = ϕ 0 and the elastic one at edge ϕ = ϕ 1 .
First, it is essential to investigate the convergence performance of the present modified variational approach for vibration analysis of the curved beam. To begin with, the comparison of dimensionless frequencies Ω for a C-C circular CLCB with respect to different N ϕ is presented in Table 2. The non-dimensional frequency Ω is defined as Ω = ω R 2 (ϕ 1 − ϕ 0 ) 2 12ρ/(E 1 h 2 ). Two thicknesses h = 0.1 and 0.2 m are considered. The material properties are ρ = 1580 kg/m 3 , E 2 = 10 GPa, E 1 = 15E 2 , G 12 = G 13 = G 23 = 0.55E 2 , µ 12 = 0.27, and [α 1 f iber /α 2 f iber ] =[0 • /90 • ] (i.e., the fiber orientation angles of the first and second lamina are 0 • and 90 • , respectively). In addition, the geometric parameters are ϕ 0 = −π/3 and ϕ 1 = π/3. The Jacobi parameters are α = β = 0 and the truncation number is M = 8. The number of segments is selected as N ϕ = 2, 4, 6, 8, and 10. In addition, the results from exact solutions exploiting CBT [3] and exact serious results based on FSDT [14] are shown for comparison. As N ϕ is increased, fast convergence of Ω can be observed. In addition, the present results coincide well with the ones by FSDT and CBT. As a consequence, with small N ϕ , accurate results of Ω can be achieved. For subsequent cases, unless otherwise specified, N ϕ = 8 is chosen by default. Secondly, the influence of penalty terms k τ (τ = u, w, ϕ, φ, ν) (see Equation (19)) on Ω is investigated. The dimensionless frequency Ω versus penalty parameters k t for a C-C curved beam is shown in Figure 3. The free vibration of CLCBs with arbitrary boundary conditions including both classical (clamped, free, and simply supported boundary conditions) and elastic boundary conditions is investigated. Here, for the convergence performance of penalty terms, the typical C-C boundary conditions are selected. Note that for other boundary conditions, similar results can be achieved. The material properties, geometric parameters, Jacobi parameters (α and β), and number of truncation terms M are the same as those of Table 2. For each case, one type of k τ is varied from 10 0 to 10 16 while the others are unchanged (= 10 14 ). When penalty parameters are absent or small, pseudo-rigid modes might emerge, indicating that the continuity conditions of the interface may not be imposed properly. By augmenting k τ , the continuity condition can be satisfied. It is observed that when k τ is in the range of 10 10 -10 16 , the solution becomes very stable, with Ω remaining the same. Therefore, k τ(u,v,ϕ) = 10 14 and k τ(φ,ν) = 10 12 are employed as the coupling parameters in the following discussions. It should be noted that for different types of curved beams (e.g., elliptical, paraboloidal, and hyperbolical ones), when k τ(u,v,ϕ) > 10 10 and k τ(φ,ν) > 10 8 , similarly, Ω converges to certain values (not shown). Thus, the determinations of penalty and boundary parameters are consistent for different curved beams.
It has been mentioned before that the present approach can be applicable to the prediction of the vibration behavior of a curved laminated beam subjected to elastic boundary conditions. Hence, the impact of boundary parameters on the vibration characteristics should be analyzed. The dimensionless frequency Ω versus boundary parameters k τ for a C-E curved laminated beam is shown in Figure 4, with the other parameters being the same as those of Figure 3. The edge ϕ = ϕ 0 is clamped while the border ϕ = ϕ 1 is elastically supported. For the edge ϕ = ϕ 1 , one set of boundary springs vary between 10 0 and 10 16 while the others are presumed to be infinity (= 10 14 ). The variations of Ω with k τ are similar to those in Figure 3, and the results can reflect that the determination of k τ in Table 1 for different boundary conditions should be appropriate. It should be noted that for different types of curved beams (e.g., elliptical, paraboloidal, and hyperbolical ones), when ( ) , , u v k τ ϕ > 10 10 and ( ) , k τ φ ν > 10 8 , similarly, Ω converges to certain values (not shown). Thus, the determinations of penalty and boundary parameters are consistent for different curved beams. It has been mentioned before that the present approach can be applicable to the prediction of the vibration behavior of a curved laminated beam subjected to elastic boundary conditions. Hence, the impact of boundary parameters on the vibration characteristics should be analyzed. The dimensionless frequency Ω versus boundary parameters k τ for a C-E curved laminated beam is shown in Figure 4, with the other parameters being the same as those of Figure 3. The edge φ = φ0 is clamped while the border φ = φ1 is elastically supported. For the edge φ = φ1, one set of boundary springs vary between 10 0 and 10 16 while the others are presumed to be infinity (= 10 14 ). The variations of Ω with k τ are similar to those in Figure 3, and the results can reflect that the determination of k τ in Table 1 for different boundary conditions should be appropriate.  Thirdly, the Jacobi orthogonal polynomials including the number of truncation terms and Jacobi parameters (α and β) should be investigated. The dimensionless frequency Ω versus truncation terms M for a circular CLCB and elliptical CLCB with the C-C boundary condition is displayed in Figure  5a,b, respectively, the other parameters remaining the same as those of  Thirdly, the Jacobi orthogonal polynomials including the number of truncation terms and Jacobi parameters (α and β) should be investigated. The dimensionless frequency Ω versus truncation terms M for a circular CLCB and elliptical CLCB with the C-C boundary condition is displayed in Figure 5a,b, respectively, the other parameters remaining the same as those of Figure 3. Obviously, for both cases, Ω converges very fast with the increase in M. As M becomes larger than 8, Ω (first five modes) converges to certain values. Subsequently, M = 8 is selected for all the numerical examples. To demonstrate the effect of Jacobi parameters (α and β) on the vibration behavior of the curved beam, the percentage error of Ω for various combinations of α and β is illustrated in Figure 6. Several combinations of (α, β) including (α, β) = (−0.5, −0.5), (−0.5, 0), (0.5, 0), (0, 0.5), (0.5, 0.5), and (1, 1) are chosen. The results of (α, β) = (0, 0) are used as the base results, with the percentage error defined as Ω α,β − Ω α=0,β=0 /Ω α=0,β=0 . In addition, three types of boundary restraints are considered: F-F, F-C, and C-C boundary conditions. It is shown that the maximum value of the percentage error is no more than 2 × 10 −2 , which may indicate that the Jacobi parameters may not affect the value of Ω. Thus, different admissible displacement functions with various Jacobi parameters can be used to establish the formulation, leading to flexible selections of admissible displacements. For the following calculations, unless otherwise stated, (α, β) = (0, 0) is utilized.  After comprehending the characteristics of several parameters, it is necessary to show the reliability and precision of the present approach. Due to lack of experimental data, the present results are compared to those from results by other numerical studies in the literature. To further validate the present method, the finite element method (FEM) is also utilized for comparison. Table 3 shows the comparison of Ω (first eight modes) for the C-C elliptical, paraboloidal, and hyperbolical CLCBs.   After comprehending the characteristics of several parameters, it is necessary to show the reliability and precision of the present approach. Due to lack of experimental data, the present results are compared to those from results by other numerical studies in the literature. To further validate the present method, the finite element method (FEM) is also utilized for comparison. Table 3 shows the comparison of Ω (first eight modes) for the C-C elliptical, paraboloidal, and hyperbolical CLCBs.  After comprehending the characteristics of several parameters, it is necessary to show the reliability and precision of the present approach. Due to lack of experimental data, the present results are compared to those from results by other numerical studies in the literature. To further validate the present method, the finite element method (FEM) is also utilized for comparison. Table 3 shows the comparison of Ω (first eight modes) for the C-C elliptical, paraboloidal, and hyperbolical CLCBs. The material parameters are ρ = 1500 kg/m 3 , E 2 = 10 GPa, E 1 = 15E 2 , G 12 = G 13 = 0.5E 2 , G 23 = 0. for all the cases. The numbers of elements for the elliptical, paraboloidal, and hyperbolical beams are 5808, 1310, and 1330, respectively. As shown in Table 3, the errors of the results by the present method and those by the FEM are not larger than 2.92%, which is small and indicates the accuracy of the present methodology. The natural frequencies depend on the different radii of curvature, lamination schemes, and boundary conditions. The differences of frequencies at the C-C boundary condition with respect to different types of CLCB (i.e., elliptical, paraboloidal, and hyperbolical) should be related to the radii of curvature changes. Then, the comparison of Ω for circular CLCBs subjected to various boundary conditions is displayed in Table 4. Results from the decomposition approach [17] and wave solution method [35] are given as a contrast. Three boundary restraints F-F, F-C, and C-SS are taken into account. For all the cases, the frequencies are in complete agreement. To further validate the present method, several mode shapes for curved composite laminated beams including elliptical, paraboloidal, and hyperbolical ones in the x-z plane are presented in Figure 7. The mode shapes obtained by the FEM through ABAQUS are presented for comparison. Apparently, all the mode shapes from the present methodology and FEM match well with each other. On the whole, the present methodology is capable of solving the vibration problem of the CLCB subjected to arbitrary boundary conditions. Finally, several novel results are presented, which can be served as benchmark solutions. The frequencies (first five modes) for the elliptical, paraboloidal, and hyperbolical CLCB with diverse lamination schemes and boundary conditions are shown in Tables 5-7   Several mode shapes of elliptical ( Figure A1), paraboloidal ( Figure A2), and hyperbolical ( Figure A3 Tables 5-7, respectively. To demonstrate the material properties, the mode shapes of isotropic beams (elliptical, paraboloidal, and hyperbolical) are also given in Figures A1-A3. The material properties of the isotropic beams are ρ = 1500 kg/m 3 , E = 100 GPa, and µ = 0.25. The following aspects should be noted. First, for all the cases, the mode shapes for laminated composite beams and isotropic beams are different. With regard to the F-C boundary condition, obvious differences can be observed for the first-order mode, while for the C-C boundary condition, the mode shapes are similar for the first-and second-order modes, although distinct differences can be observed for the third-order mode. Secondly, with various radii of curvature changes and boundary conditions, the corresponding mode shapes are also quite different. Thirdly, by choosing different lamination schemes, the material properties change, leading to variations in mode shapes of the CLCBs. By properly selecting the lamination schemes, high strength ratio and corrosion resistance can be obtained, which may make the laminated composite materials perform better than others.
validate the present method, several mode shapes for curved composite laminated beams including elliptical, paraboloidal, and hyperbolical ones in the x-z plane are presented in Figure 7. The mode shapes obtained by the FEM through ABAQUS are presented for comparison. Apparently, all the mode shapes from the present methodology and FEM match well with each other. On the whole, the present methodology is capable of solving the vibration problem of the CLCB subjected to arbitrary boundary conditions.

Conclusions
This paper proposes a unified formulation for dynamic analysis of different types of CLCBs subjected to general boundary conditions. In the framework of HSDT, an improved variational approach along with a multi-segment partitioning technique is exploited to construct the theoretical model. The displacement components of each curved beam segment are presented in terms of Jacobi orthogonal polynomials. Through a series of numerical cases, the fast convergence performance and stable, highly efficient, and precise features of the current methodology have been demonstrated. It has been proved that the present methodology can be applicable for vibration cases with respect to both classic and elastic boundary conditions. Furthermore, it should be noted that the current approach leads to flexible choices of admissible displacement functions, which is one prominent advantage compared to other methods. This paper presents the free vibration results, while the dynamic analysis of CLCBs under external excitation by the present approach has not been studied. Besides, linear geometric analysis of CLCBs by a unified formulation has been conducted. The present approach should also be appropriate for non-linear geometric analysis of CLCBs by considering the non-linear geometric effect (e.g., introducing the damping force). These will be considered in our future study.
Author Contributions: Conceptualization, X.Z.; validation, H.L. and Y.Y.; formal analysis, B.Q.; writing-original draft preparation, B.Q.; writing-review and editing, X.Z. and Q.W.; funding acquisition, B.Q., X.Z. and Q.W. All authors have read and agree to the published version of the manuscript.