Nonlinear Dynamics of the High-Speed Rotating Plate

High speed rotating blades are crucial components of modern large aircraft engines. The rotating blades under working condition frequently suffer from the aerodynamic, elastic and inertia loads, which may lead to large amplitude nonlinear oscillations. This paper investigates nonlinear dynamic responses of the blade with varying rotating speed in supersonic airflow. The blade is simplified as a pre-twist and presetting cantilever composite plate. Warping effect of the rectangular cross-section of the plate is considered. Based on the first-order shear deformation theory and von-Karman nonlinear geometric relationship, nonlinear partial differential dynamic equations of motion for the plate are derived by using Hamilton’s principle. Galerkin approach is applied to discretize the partial differential governing equations of motion to ordinary differential equations. Asymptotic perturbation method is exploited to derive four-degree-of-freedom averaged equation for the case of 1 : 3 internal resonance-1/2 sub-harmonic resonance. Based on the averaged equation, numerical simulation is used to analyze the influence of the perturbation rotating speed on nonlinear dynamic responses of the blade. Bifurcation diagram, phase portraits, waveforms and power spectrum prove that periodic motion and chaotic motion exist in nonlinear vibration of the rotating cantilever composite plate.


Introduction
High-speed rotating blades are essential components of modern large aircraft engines.Blades are designed with large aspect ratios and thin-wall structure in order to raise operational efficiency.In the actual working condition, the rotating blades are frequently subjected to the aerodynamic, elastic, and inertia loads.Various types of excitation lead to large amplitude nonlinear parametric vibrations of the blades, which can result in the resonance phenomena and undesirable disasters, especially when the rotating blades operate with high speed and huge centrifugal force.According to the investigation, the vibration failure of the aircraft engine is more than 60% of the total failure, while the vibration failure of the blade accounts for more than 70% of the total vibration failure.Resonance and flutter produced by the forced vibration and self-excited vibration are the main reason leading to the blade failure.So, it is very important that a reasonable model is established to accurately predict nonlinear vibration characteristics and other complex dynamics of the blade.
In recent years, considerable attention has been given to the studies on the vibration characteristics of rotating blades.Lin et al. [1] deduced governing differential equations and general elastic boundary conditions of a nonuniform pretwist Timoshenko beam by using Hamilton's principle.Yao et al. [2] investigated the nonlinear dynamic responses of the thin-walled rotating cantilever beam.Sarkar and Ganguli [3] assumed the modal function as polynomials satisfying all four boundary conditions and discussed free vibration of the nonhomogeneous Timoshenko beam.Georgiades et al. [4] derived equations of motion of a rotating composite Timoshenko beam by utilizing Hamilton's principle.Xie et al. [5,6] numerically investigated the effect of symmetric and asymmetric shroud gaps, rotational speeds, and the aerodynamic force amplitude on dynamic characteristics of the rotating Euler-Bernoulli beam with a mass point at the free end.Huang and Kuang [7] investigated the effect of a near root local blade crack on the stability of a bladed disk, in which an individual blade is modeled simply as a cantilever Euler-Bernoulli beam.
Currently, plenty of theoretical analyses are focused on the beam model of the blade, such as the simple Euler-Bernoulli beam and the pretwist Timoshenko beam.But, theoretical studies on the plate model of the blade are still few.Young and Liou [8] established equations of a rotating cantilever plate with a time-varying speed and numerically investigated the effect of the Coriolis force on boundaries of the unstable regions.Yoo and Kim [9] derived linear equations of motion and analyzed free vibration of rotating cantilever plates.Sinha and Turner [10] derived governing partial differential equations of motion for a rotating pretwist plate to investigate the static and dynamic frequencies of the blade.Li and Zhang [11] presented a dynamic model of a functionally graded rectangular plate undergoing large overall motions.Kouchakzadeh et al. [12] analyzed the nonlinear aeroelasticity of a laminated composite plate in supersonic airflow by using the classical plate theory.Malgaca et al. [13] tried to control the vibration of the rotating blade at different speeds by utilizing root-embedded piezoelectric materials.Kam et al. [14] proposed a procedure to analyze the structural failure of a composite wind blade subjected to quasi-static loads.Tang and Chen [15] presented the solvability conditions of nonlinear partial differential equations for in-plane moving plates of two cases, in which the one case is without internal resonance and another case is under internal resonances, respectively.Sun et al. [16] used a quadratic layerwise theory and a new dynamic model to study dynamic behaviors for a multilayer pretwist rotating blade.Wang et al. [17] established a time-dependent nonlinear model of a flexible bladerotor-bearing system by using the Lagrange method.Shariyat et al. [18] performed the nonlinear dynamic analysis of rectangular composite plates.Zhang and Li [19] adopted the Lagrangian method to acquire dynamic equations of motion for the pretwist and predeformed rotating cantilever plate subjected to the harmonic aerodynamic force.Mendonça et al. [20] considered internal damping in the shaft to study dynamic behaviors of the rotors mounted on composite shafts.Banichuk et al. [21,22] investigated the stability and bifurcation of the rotating blade under different conditions.
Based on the shallow shell theory and the Ritz method, Leissa et al. [23] explored frequencies and mode shapes of turbomachinery blades with the coupling of bending and twist.Kee and Kim [24] assumed blades as the moderately thick open cylindrical shell models.Sun et al. [25] applied the general shell theory to investigate the influence of parameters on natural frequency and damping characteristics of the shell model blade.Sinha and Zylka [26] simplified the rotating pretwisted turbomachinery cantilevered airfoil as an anisotropic shell and derived the free vibration equations of motion for the transverse deflection of the shell including the warping effect.Volker and Joachim [27] analyzed the aeroelastic phenomena of twenty compressor blades simplified as the spring-damper models mounted on the hub.Ekici et al. [28] proposed a nonlinear harmonic balance method to compute the unsteady selfexcited aerodynamic of asymmetry turbomachinery blades.Farhadi and Hosseini-Hashemi [29] studied aeroelastic behaviors of a rotating thick plate in the supersonic airflow.Lachenal et al. [30] presented the design, analysis, and realization of a zero stiffness twist morphing wind turbine blade subjected to the gust loads.
Dynamic modeling of the rotating blade requires an accurate expression of the aerodynamic force.But, the physical mechanism of aerodynamic interaction has been mainly investigated by experimental methods.Models have been proposed to deal with aerodynamic interaction problems, which is done with by the piston theory widely.Ashley and Zartarian [31] firstly proposed the quasi-steady piston theory to deal with aerodynamic interaction problems.Navazi and Haddadpour [32] studied the thermal stability of the functionally graded plate subjected to the aerodynamic load obtained from the first-order piston theory.According to large deformation geometric relationship, the piston theory and the quasi-static thermal stress theory, Yuan and Qiu [33] established the aerodynamic model of a composite stiffened panel and used Hamilton's principle to derive the equations of motion for the system.Yang et al. [34] applied a modified local piston theory to analyze aeroelastic behaviors of curved panels.
Although extensive studies have been carried out on rotating cantilever beams, studies on nonlinear dynamic responses of the pretwist, presetting rotating cantilever plates are still few.In this paper, nonlinear dynamic behaviors of the blade with varying rotating speeds under the supersonic airflow are investigated.Considering the shear deformation and the warping effect, equations of motion for the cantilever plate are derived by using Hamilton's principle.The Galerkin approach is applied to discretize the partial differential governing equations of motion to ordinary differential equations.The asymptotic perturbation method is exploited to obtain averaged equations of the system in the case of 1 : 3 internal resonance-1/2 subharmonic resonance.Based on the averaged equations, numerical simulation is applied to investigate the bifurcation and chaotic dynamics of the rotating cantilever composite plate.In order to analyze the internal resonance, we choose the perturbation rotating speed as the controlling parameter to investigate nonlinear behaviors of the pretwist and presetting rotating cantilever plate.From the results of the numerical simulation, it is found that the system performs periodic and chaotic motions under specific conditions.It is observed that the perturbation rotating speed has a significant influence on the nonlinear dynamic behaviors of the rotating plate.Since we can control the responses of the system from the chaotic motions to the periodic motions by changing the perturbation rotating speed, we can control the large amplitude nonlinear vibrations of the blade.

Equations of Motion for the Rotating Cantilevered Blade
The schematic diagram of the rotating cantilever blade is shown in Figure 1.The blade is simplified as a pretwist and 2 International Journal of Aerospace Engineering presetting rotating cantilever plate in the following dynamical analysis.Rotating speed of the blade is Ω t = Ω 0 + Ω 1 cos ωt , where Ω 0 is the steady-state rotating speed and Ω 1 cos ωt is the periodic perturbation.The shape of the blade is a rectangular plate, which is characterized by the span length L, the chord length C, and the thickness h.The cross-section of the cantilever plate is shown in Figure 2.
The plate is clamped to a rigid hub of radius r.There are a presetting angle θ r in the fixed end and a pretwist angle θ R at the free end.i, j, k is the inertial frame and the origin of it is in the center of the hub.At the edge of the hub, there is the rotating coordinate system e 0 x , e 0 y , e 0 z .We define that the spanwise direction is e 0 x , the chordwise direction is e 0 y , and the thickness direction is e 0 z .Another local triad system êx , êy , êz attached to the free end of the plate is introduced on the plate, which is called the sectional coordinate system.
The twist rate θ ′ of the blade is expressed as and the twist angle θ x of the blade is written as Based on the Kirchhoff hypothesis, the first-order shear deformation theory including the warping effect is considered to establish the displacement field.Displacements of any point along x, y, and z directions can be expressed by the displacement of the neutral plane of the plate as follows [35]: where ε 0 is the membrane strain and ε 1 is the flexural strain.
For the advanced fiber-reinforced composite material blade, the constitutive relation of the composite plate is expressed as follows:  21 , where where k represents the shear coefficient, A ij is called the extensional stiffness, D ij is the bending stiffness, and B ij is the bending-extensional coupling stiffness, which are defined as I 0 , I 1 , and I 2 are mass moments of inertia, which are defined as The rotating plate mounted in the hub is rotating with the variable speed.The displacement of a random point P without deformation on the plate is expressed as r x, y, z, t = R 0 + x e x + y e y + z e z , 7a and the deformation displacement of a random point P on the rotating plate is written as where u, v, and w denote displacement components along x, y, and z directions, respectively.
When the blade rotates with a constant angular velocity ω, the transient angle Θ can be written as The instantaneous direction of the local unit vector e x , e y , e z for any typical point on the blade with respect to the stationary global Cartesian unit vector i, j, k is given by The speed and the acceleration of any random point P on the plate are described as follows: International Journal of Aerospace Engineering The variation of the kinetic energy is given by A large centrifugal force will be generated when the compressor blade of the aircraft engine operates with the high rotating speed.Components of centrifugal forces along the spanwise and chordwise directions are described by F s and F c , which can be expressed as The variation of the potential energy δU consists the variation of the strain energy δU 1 and the variation of the centrifugal potential energy δU 2 , namely, δU = δU 1 + δU 2 , where Aerodynamic forces of the blade derived by the firstorder piston theory are written as where ρ ∞ is the air density, C ∞ is the speed of sound, P is the speed of the airflow, and U t y = P sin θ, U t z = P cos θ.The variation of the virtual work of aerodynamic forces is expressed as follows: Substituting ( 11), (13a), (13b), and ( 15) into ( 16), the nonlinear partial governing equations for nonlinear vibration of the rotating blade expressed in terms of displacement variables can be derived as follows: The boundary conditions are given as follows: The dimensionless variables and parameters are introduced as follows: Substituting ( 23) into ( 17), ( 18), ( 19), (20), and ( 21), we obtain the dimensionless governing equations of motion for the blade.
Based on the practical working condition of the blade, and theoretical and numerical studies given by Dowell et al. [36,37], it is known that vibrations of the first two order modes for the blade play an important role during vibration.The Galerkin approach is applied to obtain ordinary differential equations of motion in the dimensionless form.The Galerkin approach is derived by the Taylor expansion method, which is a mathematically convergent method.The first two order mode functions can produce convergent results.Mode functions of the plate u x, y, t , v x, y, t , w x, y, t , ϕ x x, y, t , and ϕ y x, y, t are given as follows based on a combination method of beam functions where where k i are coefficients related to frequencies, k 4 i = k 4 j = ω 2 ρA/EJ, κ i and λ j are coefficients of beam functions as follows: Substitute (24a), (24b), (24c), (24d), and (24e) into the dimensionless governing equations of motion and express u, v, φ x , and φ y with w.Ordinary differential equations of transverse vibration for the first mode and the second mode of the rotating blade in the dimensionless form are derived through the Galerkin method as follows: 8 International Journal of Aerospace Engineering where expressions of coefficients α i , β i , and f ij are described in the appendix.

Perturbation Analysis
There are both square nonlinear terms and cubic terms in (25a) and (25b), so the asymptotic perturbation method is exploited to derive average equations.The small parameter ε is introduced.The scale transformation is performed as follows: The case of 1 : 3 internal resonance-1/2 sub-harmonic resonance for the rotating plate is considered.Relationship between ω 1 and ω 2 is expressed as where σ 1 and σ 2 are two detuning parameters.Substituting ( 26) and ( 27) into (25a) and (25b), we obtain the following equation The time scale transformation is introduced as τ = ε q t, 29 9 International Journal of Aerospace Engineering where q is the positive rational number, which will be determined in the following deduction.
Functions x t and y t are the solutions of (28a) and (28b), which can be expressed as a power series of small parameter ε.
In (30a) and (30b), if n ≠ 0, then r n = n , otherwise r 0 = r.In the following derivation, we will fix r n .
The real part of x t and y t are written as follows: where the item with the asterisk on the right-hand side is the complex conjugate of the item on the left-hand side.Solutions of (28a) and (28b) can be rewritten more explicitly as follows: where cc denotes terms of the complex conjugate of functions on the right-hand side.
Solutions of (28a) and (28b) can be expressed as harmonic functions, whose coefficients are related to small parameter ε, as shown in (32a) and (32b).Functions ψ n τ, ε and φ n τ, ε are given by We suppose when ε → 0, the limits of ψ n τ, ε and ϕ n τ, ε exist.Meanwhile, we define ψ 0 n = ψ n and ϕ 0 n = ϕ n for n ≠ 1 and ψ 0 1 = ψ and ϕ 0 1 = ϕ for n = 1.When n = 2, we obtain the derivative results as follows: In order to solve the coefficients of ψ n τ, ε and ϕ n τ, ε , substituting (30a) and (30b) into (28a) and (28b), we obtain equations for each harmonic with order n and a fixed order of approximation on the perturbation parameter ε.
When n = 0, we get the following expressions: The solvability of (35a) and (35b) requires r = 2. Thus, we obtain the following expressions: When n = 2, considering derivative forms of (34a) and (34b) yields Then, the following equations are obtained When n = 1, for the balance of the nonlinear and linear terms, we must make q = 2. Equations can be derived as follows: International Journal of Aerospace Engineering Substitute (36a) and (36b) and (38a) and (38b) into (39a) and (39b), differential equations about ψ 1 and ϕ 1 are derived as Make ψ 1 and ϕ 1 as follows: Substitute (41a) and (41b) into (40a) and (40b), averaged equation in the Cartesian form is obtained as where expressions of coefficients σ i , g i , h i , and k i are described in the appendix.

Numerical Simulation of Nonlinear Vibrations
4.1.Comparison about Frequencies.In this section, based on the kinetic energy and the potential energy of the blade deduced in the previous section, the Chebyshev polynomial and Ritz method are used to calculate the frequency of the model, which are compared with the frequency obtained by finite element method.
To validate the model proposed, we compare frequencies theoretically calculated with frequencies obtained by the finite element method.Geometrical parameters and physical properties of the compressor blade are selected as the spanwise length L = 0 0534 m, the chordwise length C = 0 03773 m, the mounting angle θ = 34 49 o , the hub radius r = 0 0845 m, the density ρ = 7800 kg/m 3 , the Poisson's ratio γ = 0 3, the elastic modulus E = 1 96 GPa, the inlet velocity v in = 39 0496 m/s, the inlet mass flow M in = 2 557 kg/s, the real flow velocity V in = 78 472 m/s.Table 1 shows frequencies of different orders obtained by two methods above.When the rotating speed increases from 0 rpm to 15000 rpm, frequencies calculated by those two methods are in good agreement.So, the model is reasonable.

Numerical Simulation.
Based on (42a), (42b), (42c), and (42d), we consider the influences of dimensionless parameters on nonlinear vibration behaviors of the rotating cantilever blade.Numerical simulation is utilized to investigate nonlinear dynamics of the rotating cantilever composite plate subjected to the aerodynamic force and the centrifugal force.Nonlinear oscillation of the system is investigated by choosing the perturbation rotating speed Ω 1 as the controlling parameter.When other parameters and the initial condition do not vary, we only change Ω 1 to detect the influence of the periodic perturbation rotating speed on vibration of the rotating cantilever composite blade.The Runge-Kutta algorithm and the Poincare map theory are used to construct numerical results of the bifurcation diagram, which describes the vibration law of the displacement x 1 , when Ω 1 changes in a certain region.
The chaotic and periodic responses can be identified by several conventional criteria.Thus, based on the bifurcation diagram, the waveform, phase portraits, and the power spectrum are utilized to further verify the existence of the chaotic and periodic motion of the blade.Figures 4 and 5 illustrate that the system occurs multiperiodic motion when the  (b) the waveform on the plane t, x 1 ; (c) the phase portrait on the plane x 3 , x 4 ; (d) the waveform on the plane t, x 3 ; (e) the phase portrait in the three-dimensional space x 1 , x 2 , x 3 ; and (f) the power spectrum on the plane f , x 1 .

12
International Journal of Aerospace Engineering perturbation rotating speed Ω 1 is chosen as Ω 1 = 0 076 and Ω 1 = 0 412, respectively.Figures 4(a) and 4(c) represent the phase portraits on the planes x 1 , x 2 and x 3 , x 4 , respectively.Figures 4(b) and 4(d) give the waveform on the plane t, x 1 and t, x 3 , respectively.Figure 4(e) represents the three-dimensional phase portrait in the space x 1 , x 2 , x 3 .Figure 4(f) describes the power spectrum on the plane f , x 1 .It can be shown from Figure 4 that the amplitude of the first-order mode is larger than that of the secondorder mode.Figure 6 indicates that the chaotic motion of the system appears when Ω 1 = 0 8. Figures 7 and 8 illustrate that multiperiodic motion of the rotating composite plate occur, when the perturbation rotating speed Ω 1 continues to increase and is selected as Ω 1 = 0 93 and Ω 1 = 1 1,  (b) the waveform on the plane t, x 1 ; (c) the phase portrait on the plane x 3 , x 4 ; (d) the waveform on the plane t, x 3 ; (e) the phase portrait in the three-dimensional space x 1 , x 2 , x 3 ; and (f) the power spectrum on the plane f , x 1 .
13 International Journal of Aerospace Engineering respectively.Figure 9 shows that the chaotic motion of the rotating cantilever composite plate exists when Ω 1 = 1 201.Figure 10 illustrates that the system exhibits multiperiodic motion when Ω 1 increases to Ω 1 = 1 33.
Numerical simulation indicates that the system performs complex dynamic behaviors, such as multiperiodic motion and chaotic motion, when Ω 1 changes in a certain region.
The perturbation rotating speed Ω 1 affects the dynamic behaviors of the system significantly.

Conclusions
In this paper, nonlinear vibration behaviors of the high rotating cantilever plate subjected to the centrifugal force    The bifurcation diagram, phase portrait, and power spectrum are conducted to demonstrate that periodic motions and chaotic motions occur in nonlinear vibrations of the rotating blade under certain conditions.In this paper, the perturbation rotating speed is selected as the controlling parameter.Numerical simulation shows that nonlinear dynamic motions of the blade are sensitive to the perturbation rotating speed.The periodic motion and chaotic motion of the blade appear alternately.Periodic motions and chaotic motions exist in the averaged equations.It is well known that  16 International Journal of Aerospace Engineering periodic motions and chaotic motions in the averaged equations can lead to amplitude-modulated periodic and chaotic vibrations in the original system under certain conditions.Therefore, amplitude-modulated periodic and chaotic motions occur in the rotating blade.Occurrence of the chaotic motion means that the rotating cantilever blade   International Journal of Aerospace Engineering may perform large amplitude nonlinear vibration, which leads to the damage of the system.We can control the responses of the system from the chaotic motions to the periodic motions by changing the perturbation rotating speed.Finally, the large-amplitude vibrations of the blade can be controlled.Thus, study on nonlinear behavior of the blade is of great significance.
Analytical results of this paper can be widely used in the practical aircraft engines.Large amplitude nonlinear vibration produced by the forced vibration and selfexcited vibration occasionally leads to the undesirable disaster of the blade.Since chaotic motion is the large amplitude nonlinear vibration, we can adjust the perturbation rotating speed to control the responses of the system in order to avoid the appearance of the chaotic motion.Moreover, the possibility of aircraft engine failure caused by nonlinear vibration is reduced to ensure the safety of the aeroplane.

Appendix Coefficients in Equations
Coefficients represented in (25a) and (25b) are given as follows:

Figure 1 :Figure 2 :
Figure 1: The model of the pretwist rotating cantilever rectangular plate.

S−S
Δp z δwd x d y + S Δp y δvd x d z = Δp z δw + ∂z ∂y Δp y δv d x d y 15 Nonlinear dynamic equations of motion for the rotating blade are established by using Hamilton's principle.t 0 δK − δU + δW d t = 0 16

Figure 4 :
Figure 4: The multiperiodic motion of the blade with varying rotating speeds when Ω 1 = 0 076.(a) The phase portrait on the plane x 1 , x 2 ; (b) the waveform on the plane t, x 1 ; (c) the phase portrait on the plane x 3 , x 4 ; (d) the waveform on the plane t, x 3 ; (e) the phase portrait in the three-dimensional space x 1 , x 2 , x 3 ; and (f) the power spectrum on the plane f , x 1 .

Figure 5 :
Figure 5: The multiperiodic motion of the blade with varying rotating speeds when Ω 1 = 0 412.(a) The phase portrait on the plane x 1 , x 2 ; (b) the waveform on the plane t, x 1 ; (c) the phase portrait on the plane x 3 , x 4 ; (d) the waveform on the plane t, x 3 ; (e) the phase portrait in the three-dimensional space x 1 , x 2 , x 3 ; and (f) the power spectrum on the plane f , x 1 .

Figure 6 :
Figure 6: The chaotic motion of the blade with varying rotating speeds when Ω 1 = 0 8. (a) The phase portrait on the plane x 1 , x 2 ; (b) the waveform on the plane t, x 1 ; (c) the phase portrait on the plane x 3 , x 4 ; (d) the waveform on the plane t, x 3 ; (e) the phase portrait in the three-dimensional space x 1 , x 2 , x 3 ; and (f) the power spectrum on the plane f , x 1 .

Figure 7 :
Figure 7: The multiperiodic motion of the blade with varying rotating speeds when Ω 1 = 0 93.(a) The phase portrait on the plane x 1 , x 2 ; (b) the waveform on the plane t, x 1 ; (c) the phase portrait on the plane x 3 , x 4 ; (d) the waveform on the plane t, x 3 ; (e) the phase portrait in the three-dimensional space x 1 , x 2 , x 3 ; (f) the power spectrum on the plane f , x 1 .

Figure 8 :
Figure 8: The multiperiodic motion of the blade with varying rotating speeds when Ω 1 = 1 1.(a) The phase portrait on the plane x 1 , x 2 ; (b) the waveform on the plane t, x 1 ; (c) the phase portrait on the plane x 3 , x 4 ; (d) the waveform on the plane t, x 3 ; (e) the phase portrait in the three-dimensional space x 1 , x 2 , x 3 ; and (f) the power spectrum on the plane f , x 1 .

Figure 9 :
Figure 9: The chaotic motion of the blade with varying rotating speeds when Ω 1 = 1 201.(a) the phase portrait on the plane x 1 , x 2 ; (b) the waveform on the plane t, x 1 ; (c) the phase portrait on the plane x 3 , x 4 ; (d) the waveform on the plane t, x 3 ; (e) the phase portrait in the three-dimensional space x 1 , x 2 , x 3 ; and (f) the power spectrum on the plane f , x 1 .

Figure 10 :
Figure 10: The multiperiodic motion of the blade with varying rotating speeds when Ω 1 = 1 33.(a) The phase portrait on the plane x 1 , x 2 ; (b) the waveform on the plane t, x 1 ; (c) the phase portrait on the plane x 3 , x 4 ; (d) the waveform on the plane t, x 3 ; (e) the phase portrait in the three-dimensional space x 1 , x 2 , x 3 ; and (f) the power spectrum on the plane f , x 1 .
E 1 is the longitudinal modulus of the fiber, E 2 is the transverse modulus, γ 12 and γ 21 are the major Poisson's ratio, respectively, and G 12 , G 23 , and G 13 are the shear modulus, respectively.In-plane force resultants N xx , N yy , and N xy ; the moment resultants M xx , M yy , and M xy ; and the transverse forces Q x and Q y are written in matrices as follows: