Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014 Porto, Portugal, 30 June - 2 July 2014 A. Cunha, E. Caetano, P. Ribeiro, G. Müller (eds.) ISSN: 2311-9020; ISBN: 978-972-752-165-4 Reduced-order models for the analysis of non-linear vibrations of variable stiffness composite laminated plates Hamed Akhavan, Pedro Ribeiro Department of Mechanical Eng., Faculty of Engineering, University of Porto, R. Dr. Roberto Frias, 4200-465 Porto, Portugal email: [email protected], [email protected] ABSTRACT: This study focuses on reduced-order models to analyse non-linear oscillations of rectangular variable stiffness composite laminated (VSCL) plates. Each lamina of these VSCL plates has curvilinear fibres that are distributed by shifting along one of the Cartesian directions. Two model reduction approaches are applied: modal reduction alone and static condensation followed by modal reduction. In both cases, a p-version finite element that takes into account the effect of large amplitude displacements in the framework of third-order shear deformation theory (TSDT) is used to compute the linear modes, which are used as a basis. In the first approach, a reduced set that contains bending and in-plane modes is selected. The mode selection is discussed and it is confirmed that, also in this TSDT based model and in the non-linear regime, in-plane modes are essential for accuracy. In the second approach, the in-plane inertia is neglected in the full model and static condensation is applied, so that the influence of in-plane displacements is kept; only then is modal reduction implemented. This model has the advantage of leading to a smaller number of degrees of freedom and of automatically taking care of the in-plane displacements, but loses the effect of in-plane inertia. The accuracy and computational performance of both models in the definition of periodic responses of VSCL plates is compared. For that purpose, the laminates are subjected to transverse harmonic loads and the frequency responses are determined using the shooting method, with fifth-order Cash-Karp Runge-Kutta method and adaptive stepsize control. For excitation frequencies around the first natural frequency, the amplitudes of vibration, the time histories and the stability status of the diverse solutions are compared. KEY WORDS: Variable stiffness; Curvilinear fibres; Non-linear; Periodic vibration; Reduced-order model. 1 INTRODUCTION Variable stiffness composite laminated (VSCL) plates are a general type of composite plates, in comparison to constant stiffness composite laminated (CSCL) plates, see Ref. [1]. In the VSCL plates investigated here, fibres follow a curvilinear path (their path is straight in CSCL plates). Mechanical properties of composite plates with curvilinear fibres can change with fibre paths. Using curvilinear fibres instead of straight ones can improve these properties in specific structures (e.g. at the corners of aircraft fuselage windows). Some favourable factors of VSCL against CSCL plates in static non-linear analyses and free linear vibration are mentioned, for example, in Refs. [2-5]. Few works [6-8] addressed free and forced non-linear vibration of VSCL plates by using thin plate theory (CPT) and first-order shear deformation theory (FSDT). In references [2-8], the p-version finite element, with hierarchic shape functions, is used. In references [6, 7] the harmonic balance method (HBM) is applied, whereas in ref. [8] Newmark method is used. A comprehensive review on VSCL plates and methods used to analyse them can be found in Ref. [1]. In the current paper, we intend to apply a p-version finite element model, based on a third-order shear deformation theory (TSDT), to VSCL plates with shifted curvilinear fibres [9]. The finite element model may be condensed statically [10] and is reduced using modal coordinates [11]. The choice of modes is analysed in two model reduction procedures. Periodic solutions of the plate under harmonic loading, in the non-linear regime, are calculated using the shooting method. The monodromy matrix, which is a by-product of the shooting method, is employed to define the stability of the solutions. The response of VSCL plates to harmonic excitations with frequencies close to their fundamental frequencies is investigated. For that purpose, the frequency response curves are given, as well as phase plots and time histories. Fourier spectra are used to evaluate the harmonic content of the solution. 2 2.1 PERIODIC VIBRATION OF VSCL PLATES p-Version finite element model Rectangular, symmetric about the middle plane, clamped laminates are displayed in Figure 1, as well as Cartesian coordinates with origin at the plate centroid. Fibres in the ith lamina can be either straight ( Τ 0 = Τ 1 ), where Τ 0 and Τ 1 are the fibre angles at the centre and vertical edges, or curvilinear ( Τ 0 ≠ Τ 1 ). The fibre angles in the ith lamina change as θi x 2 Τ1 - Τ 0 x a Τ 0 . Then, a configuration like [〈 Τ 0 , Τ 1 〉1,…, 〈 Τ 0 , Τ 1 〉i,…, 〈 Τ 0 , Τ 1 〉k]sym represents a 2k-layered symmetric laminate with different fibre angles in each layer. Displacement field, ix, y, z, t , i u, v, w , following a third-order shear deformation theory (TSDT) in x, y, z coordinates, respectively, is (with ) u u 0 z x cz 3 x w0 x , 1803 Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014 v v 0 z y cz 3 y w0 y , (1) w w0 , in which i x, y, t , i u, v, w are displacements at the midplane and i x, y, t , i x, y are rotations around axis y and x. 0 a T0 T1 h z x z Figure 1. Up is a ply with a reference curvilinear fibre path and down is a four-layer symmetric laminate. T0 and T1 have negative values in this picture. A p-version finite element method (FEM) with hierarchical basis functions that can be applied to thin and thick plates, because it is based on a TSDT [2-4], is used. The displacement components are related to the generalized displacements by where 0 T Nu 0 0 0 0 0 T Nw 0 0 N i x, y , i u, w, x , y 0 0 0 N x T 0 q u t q v t q w t q t T q t N 0 0 0 0 yz 1 3cz y w 0 y , (3) 2 xz 1 3cz 2 x w 0 x . 0 0 Q44 0 0 0 0 0 Q55 0 0 1 0 2 0 23 0 13 Q66 12 (4) where subscripts 1, 2, and 3 represent, respectively, the fibre orientation, the in-plane orientation perpendicular to 1, and the transverse direction. Stresses and strains in plate coordinates x and y are found using transformation matrices (see Ref. [14]). Comparing constant stiffness to variable stiffness composites, it is worthy to say that due to curvilinear fibres in VSCL, directions 1 and 2 are not anymore constant. Defining virtual works of inertia, internal and external forces and then applying them in the principle of virtual work, leads to the equations of motion, as shown below (details on the matrices are given in [3]). y u 0 N u T 0 v 0 0 w 0 x 0 y 0 1 Q11 Q12 Q 2 21 Q22 0 23 0 0 0 13 12 0 0 x y w 0 x w 0 y 2cz 3 2 w 0 xy , Each ply with curvilinear fibre is a non-homogenous orthotropic fibrous composite layer, in which stresses and strains are related as [14] T1 b xy u 0 y v 0 x z cz 3 x y y x (2) x y y are vectors composed of in- plane, out-of-plane and rotational shape functions, found for instance in [12, 13]. The strain-displacement relations can be approximated using the von Kármán strains by x u 0 x z cz 3 x x cz 3 2 w0 x 2 w0 x 2 , 2 y v 0 y z cz 3 y y cz 3 2 w0 y 2 w0 y 2 , M 11 0 M 22 sym K 11 K 12 L L 22 K L sym 0 0 M 33 0 0 M 34 M 44 u t 0 q q 0 v t t M 35 q w t 0 q t M 55 q q t 0 u q t 0 v q t K 35 L w K L45 q t q t K 55 L x y 0 0 K 33 L 0 0 K 34 L K L44 x y 0 0 K 13 NL q w t 23 0 0 K NL q w t 31 32 33 K NL q w t K NL q w t K NL q w t 0 0 0 0 0 0 0 0 f w t . 0 0 0 0 0 0 0 0 q u t 0 q v t 0 q w t 0 q t 0 q t x y (5) 2 z 0, 1804 The vector on the right-hand side of equation (5) represents generalized external forces. For example, taking 7 shape functions in x and y directions – equal to a vector of forty nine 2-D shape functions – for each of five variables in equation (2), leads a full model in Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014 equation (5) with n=245 DOF (hereafter, n represents the number of degrees of freedom, or “order”, of the full model). 2.2 Static condensation Due to the clamped boundaries of the plate, absence of inplane external force and the fact that the transverse deflection is not bigger than two times the plate’s thickness, the in-plane displacements are very small and the corresponding in-plane inertia is not vital to the solution. Then the in-plane inertia can be eliminated from the formulation. This elimination process is known as static condensation [10] and reduces the order of the problem. Although it generally changes the bandwidth of the stiffness matrix, in the present TSDT formulation – and as occurs with Kirchhoff [15] or FSDT [12] approaches – the bandwidth of the linear matrices are not increased by static condensation. The equations of motion (5) are condensed as w t K 33 M 35 q K 34 L L t 0 q K L44 t sym M 55 q K NL q w t 0 0 q w t f w t 0 0 q t 0 . sym 0 q t 0 M 33 sym M 34 M 44 x y with T Shooting method To apply the shooting method [16], the system of m secondorder ordinary differential equations of motion (11) is transformed to a system of 2m first-order ordinary differential equations using q m y , 0 M 3 1 K K T K K 23 NL sym K K 11 L 12 L 22 L 13 NL 23 NL . (7) If the full model is of n=245 DOF, the condensed model would be of n=147 DOF. Both equation (5) and (6) may simply be written as t K Lqt K NL qt qt f t . Mq 2.3 2.4 y 0 0 M y - M . 0 q m 0 K L K NL q m p (12) stiffness terms. Apart from the improvements in the integration of first-order differential equations (12), the procedure here applied is similar to the one described in Refs. [17, 18]. x non-linear stiffness matrix K NL is not. In the above example, the full model was of order n=245 DOF, taking m=10 modes results in the reduced modal model of m=10 DOF. Shooting with fifth-order Runge-Kutta method with stepsize adjustment algorithm based on the embedded fourth-order q w t K 35 Runge-Kutta formulas [16] (with constant parameters found L by Cash and Karp) is applied, starting by using as first guess K L45 q x t 55 y K L q y t of the solution of equation (12) only with linear (6) q m y 13 K NL K 33 NL 2 K NL linear stiffness matrix K L are diagonal, but the generalized (8) Analysis in modal coordinates If we know that the excitation of the VSCL plate centres around a specific frequency, many of the modes probably will not be excited and we can assume the forced vibration to be the superposition of only a few modes [11], let’s say m modes. The truncated modal matrix nm assembled by m normal modes ( i as eigenvectors of linear matrix) can appear as i ... m , m n (9) Then the displacements can be written as qt q m t (10) Pre-multiplying equation (8) by the transpose T reduces the n-DOF system to one with a dimension equal to the number of modes used, m-DOF. Subsequently, instead of solving n coupled equations (8), we need to solve the m coupled equations represented by m t K L q m t K NL q m t q m t pt . Mq (11) We are using the modes of the linear system in a non-linear problem; therefore, the generalized mass matrix M and the 3.1 NUMERICAL TESTS AND VALIDATION How to select modes in two reduced-order models In the first approach, reduced-order model I (model I, in brief), the full model of equation (5) is transferred to a reduced set of modal coordinates. The modes to include in the truncated modal matrix , should be the more influential bending modes along with influential in-plane modes. Here, expression “bending mode” refers to normal modes of the linear system that actually contain transverse and rotational generalized displacements; on the other hand, “in-plane mode” identifies the normal modes with only membrane generalized displacements. It is obvious from the mass and linear stiffness matrices of equation (5), of model I, that membrane and bending modes are uncoupled. We note that shear is also present in the model, and we may have bending coupled with shear, and modes where shear is more important; for the sake of simplicity, we will designate modes where bending is important as “bending modes”. In the second approach, reduced-order model II (model II), the condensed model of equation (6), which has not in-plane inertia but is still affected by the influence of in-plane displacements, is transformed to modal coordinates. In this second approach, only bending modes are selected. The procedure of selecting influential modes is programmed and implemented as an in-house code in FORTRAN. In reduced-order model I, firstly all modes of vibration of the linear problem, including in-plane normal modes, are found. Then, equation (11) is solved in modal coordinates by the shooting method, using only one bending mode. A symmetric distributed transverse load with a fixed excitation frequency is applied, because this is the type of load considered in this work. The process is repeated for all non-in-plane modes, one by one. As an example, taking a vector of forty nine 2-D 1805 Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014 shape functions for each variable of equation (2), this step (Step 1) repeats for each of 3×49 modes, where transverse displacements are coupled with cross section rotations. In this step, those bending modes that result in bigger central deflection are candidates to enter the next step. In the next step, Step 2, equation (11) is again solved, but with two or three modes. One of these modes is the bending mode selected from the previous step, while the other one or two modes is/are selected among the in-plane modes of the linear problem. This stage may be repeated for all candidate bending mode selected in Step 1, combined with different inplane modes. Again, like in the end of step one, groups of bending and in-plane modes with bigger deflections are selected. The procedure ends when enough influential modes Table 1. Frequency ratio ±0.2 ±0.6 ±1 3.2 are found. Here, the word enough means a number of modes that gives convergence in deflection computation. To select normal modes for reduced-order model II, the same procedure is applied, but the first step is sufficient, due to the lack of in-plane modes in linear model II. In this model, the influence of in-plane displacements is considered in the reduced non-linear matrix. Then, the analysis to find forced responses of vibration of VSCL plates is performed by applying enough influential modes, which were selected by one of the aforementioned procedures, resulting in two types of reduced-order models. In the following section the two types of reduced order procedures are compared. against amplitude of vibration for an isotropic plate. HFEM and shooting Ref. [18] FEM without IDI Ref. [19] Elliptic function [20, 21] Newmark method Model I (8 modal coordinates) Model I (15 modal coordinates) Model II (3 modal coordinates) Model II (15 modal coordinates) 0.2707 1.4387 0.8950 1.2114 1.0777 - 0.1180 1.4195 0.8905 1.2083 1.0700 1.2429 0.1200 1.4195 0.8951 1.2117 1.0822 1.2540 0.2400 1.4394 0.8961 1.0825 - 0.131 1.418 0.914 1.225 1.126 1.298 0.2797 1.4418 0.9044 1.2674 1.1007 1.3332 0.1338 1.4202 0.8905 1.2106 1.0744 1.2494 0.1203 1.4189 0.8904 1.2102 1.0748 1.2498 Comparison of two reduced-order models First, a square isotropic plate with fully clamped boundaries is analysed. Table 1 compares the amplitude ratios of transverse central deflection of an undamped steel plate, as a function of excitation frequency, calculated in Refs. [18-21] – where models and methods of solution different from the ones here proposed were applied – with: the proposed reducedorder models I and II and the shooting method; the full model and Newmark method. Newmark method [22] is applied to the full model, equation (5), to find deflections without transformation to the modal coordinates, providing reference solutions of the present model. Two sets of results are taken from Ref. [19]. In one of them, an approach close to a harmonic balance procedure with one harmonic is applied to an h-version finite element model. In the other, following [20] and [21], a one mode approach (one linear mode, from a model with in-plane deformation but without inertia) is employed to obtain a Duffing type equation that is exactly solved using elliptic functions, while in Ref. [18] the shooting method was applied. The plate is square with length a=500 mm and thickness h=2.0833 mm, and is subjected to a harmonic distributed load , , N/m2. In Ref. [18], the p-version finite element method (HFEM) has just three out-of-plane and six in-plane unidimensional shape functions. In the table, IDI stands for in-plane displacement and inertia. In the full model plus Newmark method, 245 DOF are used; the same 245 DOF full model is employed to obtain the reduced models that lead to the four columns with results from the present paper, solved with the shooting method. Increasing the number of modal coordinates from 3 to 15 in model II, the frequencies change slightly. With either 3 or 15 modal coordinates, model II frequencies are close to the ones provided by the full model, Newmark solution, with the exception of the far away from resonance case Wmax/h=0.2. Also, convergence to the full model solutions was not yet achieved in model I using 15 modal coordinates; the values are generally close but marginally worse than the ones of model II. It should be mentioned that the run-time of model I, with the same number of modal coordinates, is larger than the run-time of model II. Table 2. Comparison of the non-linear frequency ratio in function of the amplitude of vibration angle-ply plate, by the full model plus Newmark method. 0.2016 0.5997 1.0005 1806 p-version FEM CPT [24] 0.7340 1.0085 1.1756 0.2 0.6 1 p-version FEM FSDT [23] 0.7268 1.005 1.171 h-version FEM CPT [25] 0.7219 1.0085 1.1787 , for a CSCL Newmark method 245 DOF 0.7275 1.0097 1.1814 Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014 In Table 2, a comparison is carried out with data available on Refs. [23] and [24] (computed using the p-version FEM, FSDT and thin plate theory, or CPT, theories, respectively) and with data computed by the h-version FEM using CPT in Ref. [25]. The full model (245 DOF) is solved in the timedomain by Newmark method [22], with a small stiffnessproportional damping parameter , [26]. The CSCL plate analysed is an angle-ply three-layer plate [45, -45, 45], with the following mechanical and geometric characteristics: , , , , , , The amplitude of deflection at centre of the plate with a sinusoidal wave with 3454.665 is given. The comparison of Table 2 allows us to verify the full TSDT model. Table 3. Comparison of the amplitude of vibration of a CSCL angle-ply plate, in function of frequency, different combinations of modes in modal coordinates and a full model by Newmark method. 0.7275 1.0097 1.1814 by Newmark method, 245 DOF 0.2 0.6 1.0 Bending modes 20 modes 30 modes 0.1935 0.1935 0.5089 0.5091 0.8244 0.8239 by Reduced-order model I Including bending/in-plane modes 40 modes 19 modes 0.1935 0.1979 0.5092 0.5925 0.8243 0.9706 Table 4. Influence of bending modes and bending/in-plane modes in the amplitude of vibration in function of excitation frequency . Frequency of external load, Deflection ratio by Newmak method, full model, 600 0.1175 650 0.1275 700 0.1404 750 0.1572 800 0.1799 850 0.2115 900 0.2570 950 0.3233 1000 0.4160 1050 0.5315 1100 0.6572 1150 0.7839 1200 0.9062 1250 1300 - , using , for an angle-ply plate , reduced-order model I 5 modes 0.1174 (-0.04%) 0.1273 (-0.15%) 0.1401 (-0.24%) 0.1566 (-0.40%) 0.1780 (-1.05%) 0.2074 (-1. 94%) 0.2516 (-2.09%) 0.3120 (-3.50%) 0.3929 (-5.55%) 0.4886 (-8.07%) 0.5952 (-9.44%) 0.7046 (-10.1%) 0.7920 (-12.6%) 0.9086 1.0041 Only bending modes are taken 10 modes 15 modes 20 modes 0.1172 0.1173 0.1173 (-0.22%) (-0.19%) (-0.16%) 0.1271 0.1272 0.1272 (-0.29%) (-0.23%) (-0.22%) 0.1398 0.1399 0.1399 (-0.46%) (-0.38%) (-0.36%) 0.1562 0.1564 0.1565 (-0.61%) (-0.50%) (-0.46%) 0.1783 0.1786 0.1787 (-0.89%) (-0.71%) (-0.66%) 0.2087 0.2091 0.2094 (-1.33%) (-1.11%) (-1.01%) 0.2513 0.2524 0.2528 (-2.20%) (-1.78%) (-1.65%) 0.3112 0.3137 0.3143 (-3.75%) (-2.97%) (-2.77%) 0.3927 0.3967 0.3979 (-5.61%) (-4.65%) (-4.35%) 0.4909 0.4977 0.4998 (-7.64%) (-6.36%) (-5.96%) 0.5976 0.6069 0.6107 (-9.07%) (-7.66%) (-7.08%) 0.7054 0.7179 0.7224 (-10.0%) (-8.42%) (-7.85%) 0.8103 0.8258 0.8312 (-10.6%) (-8.87%) (-8.27%) 0.9087 0.9298 0.9363 1.0053 1.0301 1.0375 25 modes 0.1174 (-0.11%) 0.1273 (-0.15%) 0.1400 (-0.25%) 0.1567 (-0.32%) 0.1791 (-0.46%) 0.2100 (-0.70%) 0.2540 (-1.15%) 0.3172 (-1.88%) 0.4033 (-3.06%) 0.5097 (-4.11%) 0.6246 (-4.96%) 0.7406 (-5.53%) 0.8532 (-5.84%) 0.9622 1.0670 19 bending and in-plane modes are taken 0.1178 (0.23%) 0.1280 (0.37%) 0.1457 (3.76%) 0.1571 (-0.04%) 0.1798 (-0.05%) 0.2114 (-0.04%) 0.2569 (-0.03%) 0.3231 (-0.08%) 0.4163 (0.07%) 0.5320 (0.10%) 0.6575 (0.05%) 0.7839 (0.00%) 0.9065 (0.03%) 1.0241 1.1371 1807 Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014 The full model, to some extent validated in Table 2, is used as reference in Table 3, in another test of the model reduction procedures. The plate is the one of Table 2. In Table 1, increasing modal coordinates from 8 to 15 in model I, so that more in-plane modes are considered, the relation between frequency and amplitude of vibration was more accurately computed. Table 3 also shows that, in reduced-order model I and in the non-linear regime, in-plane modes are essential for accuracy. Here, four analyses of non-linear frequencies of undamped plate [45, -45, 45] are performed using modal coordinates (of a model with 245 DOF): three of them include only the first 20, 30, and 40 bending modes of vibration and the last analysis involves a combination of 19 bending and inplane modes of vibration, selected by the technique written in section 3.1. In reduced-order model I, the results demonstrate that increasing the number of bending modes alone cannot guarantee the accuracy in the non-linear regime, but taking less bending modes along with in-plane modes can profoundly improve the precision of the model. Again, it is worthy to note that reduced-order model II (that excludes in-plane inertia) already considers the in-plane displacements, during static condensation. For the same plate and loading, a smaller model (using 9 bidimensional shape functions, i.e. 45 DOF) is applied in Table 4 to compare the amplitudes of vibration computed by Newmark method (full model) with the ones computed by the shooting method (reduced-order model I). Again, different reduced-order models with different number of bending normal modes in modal coordinates (5 to 25 modes) are considered and compared with a model that incorporates a combination of 19 bending/in-plane modes chosen as explained in section 3.1. In this table, the numbers in parenthesis represent the percentage of the difference with respect to Newmark results (full model). The table verifies that deflections using reduced-order model I and taking into account only bending normal modes are always smaller than reference deflections. Away from resonance, models only with bending modes give accurate results. However, when the non-linearity is more important, deflections with 19 bending and in-plane modes are in a better agreement with the reference data, than deflections computed with more bending modes, but without in-plane modes. a) 1,5 b) CSCL - T1=45 (Stable) VSCL - T1=45 (Unstable) VSCL - T1=0 (Stable) CSCL - T1=0 (Unstable) VSCL - T1=25 (Stable) VSCL - T1=25 (Unstable) 1 1,2 CSCL - T=45 VSCL - T=0 VSCL - T=25 W/h 0,8 0,4 W/h 0 0,5 -0,4 -0,8 0 -1,2 0 0,5 1 ω/ωl1 1,5 0 CSCL - T=45 d) 0,8 VSCL - T=25 10 / 0 0,5 / 0,6 0,4 0 -0,5 0,015 CSCL - T1=45 VSCL - T1=0 VSCL - T1=25 _ VSCL - T=0 -1 0,01 Time (s) c) 20 0,005 1 0,2 -10 _∕ 0 -20 0 1 2 3 4 5 6 7 8 9 10 Figure 2. Plate’s frequency response: a) maximum transverse displacement at centre of the plate with respect to frequency with stable and unstable solutions; b), c) and d) displacements along one cycle, phase plot, and Fourier series at excitation frequency equal to 0.378 , stable solution. 1808 Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014 3.3 Forced vibration of VSCL plates with reduced-order model A test case of a clamped three-layered VSCL plate with fibre configuration [〈45,T〉,〈-45,-T〉,〈45,T〉] is analysed, taking constant Τ 0 45 in all layers and different Τ 1 T as 0, 25, 45 (while Τ 1 45 exemplifies the angle-ply CSCL plate [45,45,45]). The plate characteristics are: , , , , , , The thickness of all layers is the same. Reduced-order models with 6 modal coordinates are used in the analysis; the corresponding full models had 245 DOF (49 bi-dimensional shape functions per displacement and cross section rotation component) and the statically condensed models had 147 DOF. The responses of two VSCL plates, T=0 and 25, and a CSCL plate, T=45, to a uniform distributed harmonic excitation of amplitude 3454.665 are represented in Figure 2. The plates are undamped. Instabilities when the Floquet multiplier crosses the unit circle through +1 are found and represented by a dashed line. The displacements along one cycle, phase plots and Fourier spectra show that, at some excitation frequencies, higher harmonics have a large influence in the response. The width of the resonant region shrinks slightly by changing T from 45 (i.e. CSCL plate) to 0 or 25 (i.e. VSCL plate). If the excitation frequency is equal to 0.378 , a similar pattern of displacements along a period of vibration and similar phase plots are observed for CSCL and VSCL plates, but the vibration period changes from plate to plate. 4 CONCLUSION A procedure based on the shooting method and on a Thirdorder Shear Deformation p-version finite element was proposed to investigate the geometrically non-linear forced vibration of variable stiffness composite plates. The method can find instabilities and periodic oscillations with higher harmonics, when they occur, for instance due to the appearance of superharmonics in the response. The technique takes advantage of two reduction methods, namely static condensation and modal coordinates, to decrease the order of matrices. If static condensation is not applied, then in-plane modes must be considered in the reduced order model, to achieve accurate results in the non-linear regime. Static condensation permits to take in-plane displacements into account indirectly, so that one only has to select bending modes. ACKNOWLEDGMENTS This investigation is carried out under supports of Portuguese Science and Technology Foundation (FCT) under PhD grant SFRH/BD/81707/2011 and project PTDC/EMEPME/098967/2008, which are gratefully acknowledged. [3] H. Akhavan, P. Ribeiro and M.F.S.F. de Moura, Large deflection and stresses in variable stiffness composite laminates with curvilinear fibres, International Journal of Mechanical Sciences, 73, 14-26, 2013. [4] H. Akhavan, P. Ribeiro and M.F.S.F.de Moura, Composite laminates with linearly varying fiber angles under static and dynamic loads, Proceeding of 54th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference and co-located Events, Boston, MA, 2013, AIAA 2013-1565. [5] S. Yazdani, P. Ribeiro and J.D. Rodrigues, A p-version layerwise model for large deflection of composite plates with curvilinear fibres, Composite Structures, 108, 181-190, 2014. [6] A. Houmat, Nonlinear free vibration of laminated composite rectangular plates with curvilinear fibers, Composite Structures, 106, 211-224, 2013. [7] P. Ribeiro, Non-linear free periodic vibrations of variable stiffness composite laminated plates, Nonlinear Dynamics, 70, 2, 1535-1548, 2012. [8] P. Ribeiro and H. Akhavan, Non-linear vibrations of variable stiffness composite laminated plates, Composite Structures, 94, 8, 2424-2432, 2012. [9] C. Waldhart, Analysis of tow-placed, variable-stiffness laminates, M.Sc. thesis, Blacksburg, Virginia Tech, 1996. [10] L. Meirovitch, Computational methods in structural dynamics, Sidthoff & Noordhoff, Rockville, MD, US, 1990. [11] W.T. Thomson and M.D. Dahleh, Eds., Theory of vibration with applications, Prentice-Hall, Inc., New Jersey, US, fifth edition, 1998. [12] P. Ribeiro, A hierarchical finite element for geometrically non-linear vibration of thick plates, Meccanica, 38, 1, 117-132, 2003. [13] P. Ribeiro and M. Petyt, Nonlinear vibration of plates by the hierarchical finite element and continuation methods, International Journal of Mechanical Sciences, 41, 4-5, 437–459, 1999. [14] J.N., Reddy, Ed., Mechanics of laminated composite plates: theory and analysis, CRC Press, Boca Raton, US, 1997. [15] W. Han and M. Petyt, Linear vibration analysis of laminated rectangular plates using the hierarchical finite element method– I. Free vibration analysis, Computers & Structures, 61, 4, 705-712, 1996. [16] W.H. Press, S.A. Teukolsky, B.P. Flannery and W.T. Vetterling, Eds., Numerical recipes in Fortran, Cambridge University Press, Cambridge, New York, second edition, 1992. [17] P. Ribeiro, Non-linear forced vibrations of thin/thick beams and plates by the finite element and shooting methods, Computers & Structures, 82, 17-19, 1413-1423, 2004. [18] P. Ribeiro, Periodic vibration of plates with large displacements, AIAA Journal, 40, 1, 185-188, 2002. [19] C. Mei and K. Decha-Umphai, Finite element method for nonlinear forced vibrations of rectangular plates, AIAA Journal, 23, 7, 1104-1110, 1985. [20] J.G. Eisley, Nonlinear vibration of beams and rectangular plates, Zeitschrift fur angewandte Mathematik and Physik, 15, 167-175, 1964. [21] C.S. Hsu, On the application of elliptic functions in nonlinear forced oscillations, Quarterly on Applied Mathematics, 17, 393-407, 1960. [22] K.J. Bathe and E.L. Wilson, Eds., Numerical methods in finite element analysis, Prentice-Hall, cop., Englewood Cliffs, US, 1976. [23] P. Ribeiro, Forced periodic vibrations of laminated composite plates by a p-version, first order shear deformation, finite element, Composites Science and Technology, 66, 11-12, 1844-1856, 2006. [24] P. Ribeiro, and M. Petyt, Non-linear vibration of composite laminated plates by the hierarchical finite element method, Composite Structures, 46, 3, 197-208, 1999. [25] C.K. Chiang, C. Mei, and C.E. Gray, Finite element large-amplitude free and forced vibrations of rectangular thin composite plates, Journal of Vibration and Acoustics-Transactions of the ASME, 113, 3, 309-315, 1991. [26] S. Adhikari, Damping models for structural vibration, Ph.D. thesis, Engineering Department, Cambridge University, Cambridge, UK, 2000. REFERENCES [1] P. Ribeiro, H. Akhavan, A. Teter and J. Warminski, A review on the mechanical behaviour of curvilinear fibre composite laminated panels, Journal of Composite Materials, DOI: 10.1177/0021998313502066. [2] H. Akhavan and P. Ribeiro, Natural modes of vibration of variable stiffness composite laminates with curvilinear fibers, Composite Structures, 93, 11, 3040-3047, 2011. 1809

© Copyright 2020