A New Approach for Studying Nonlinear Dynamic Response of a Thin Plate with Internal Resonance in a Fractional Viscoelastic Medium

In the previous analysis, the dynamic behaviour of a nonlinear plate embedded into a fractional derivative viscoelastic medium has been studied by the method of multiple time scales under the conditions of the internal resonances two-to-one and one-to-one, as well as the internal combinational resonances for the case when the linear parts of nonlinear equations of motion occur to be coupled. A new approach proposed in this paper allows one to uncouple the linear parts of equations of motion of the plate, while the same method, the method of multiple time scales, has been utilized for solving nonlinear equations. The influence of viscosity on the energy exchange mechanism between interacting nonlinear modes has been analyzed. It has been shown that for some internal resonances there exist such particular cases when it is possible to obtain two first integrals, namely, the energy integral and the stream function, which allows one to reduce the problem to the calculation of elliptic integrals. The new approach enables one to solve the problems of vibrations of thin bodies more efficiently.


Introduction
It is well known that the nonlinear vibrations of plates are an important area of applied mechanics, since plates are used as structural elements in many fields of industry and technology.Different methods such as analytical [1][2][3][4], numerical [5,6], and experimental methods [5,7] could be employed to investigate the nonlinear vibrations of plates.Extensive review of recent research developments in the field could be found in Amabili [8,9] and Sathyamoorthy [10].
The study of free undamped [6,9], as well as damped [2,3,[11][12][13], nonlinear systems is important because the dynamics characteristics of the system-defined by the amplitudefrequency relations and modes of vibrations-are determined [6].Moreover, nonlinear vibrations could be accompanied by such a phenomenon as the internal resonance, resulting in multimode response with a strong interaction of the modes involved [14] accompanied by the energy exchange phenomenon.The internal resonance could be found within some combination of natural frequencies of one and the same type of vibrations.As an example, a 1:3 internal resonance was discovered in [6], when the fourth natural frequency of the out-of-plane vibrations was three times larger than the fundamental frequency of the out-of-plane vibrations of the plate.A 1:1 internal resonance was investigated in [11,15] applying a two-mode Galerkin procedure for the plate deflection in their analysis.
Another type of the internal resonance was investigated in [2,3], when one frequency of the in-plane vibrations was equal to (a 1:1 internal resonance) or twice as large as (a 1:2 internal resonance) some frequency of the out-of-plane vibrations.The combinational resonances of the additive and difference type were also investigated in [3].These types of the internal resonance result in the energy exchange between two or three subsystems.Investigations in this field were initiated by Witt and Gorelik [16], who were among the first to show theoretically and experimentally the twoto-one internal resonance with the energy exchange from one subsystem to another using the simplest two-degree-offreedom mechanical system as an example, and the interest of researchers in the problems of the internal resonance in mechanical systems does not relax.
In the given paper, nonlinear free vibrations of a thin plate described by three coupled nonlinear differential equations are considered when the plate is being under the conditions of the internal resonance resulting in the interaction of modes corresponding to the mutually orthogonal displacements.The displacement functions are determined in terms of eigenfunctions of linear vibrations.The procedure resulting in decoupling linear parts of equations is proposed with further utilization of the method of multiple scales for solving nonlinear governing equations of motion; in so doing the amplitude functions are expanded into power series in terms of the small parameter and depend on different time scales.It is shown that the phenomenon of internal resonance can be very critical, since in a thin plate internal resonance of the types two-to-one and one-to-one as well as additive and difference combinational resonances are always present.
To describe the nonlinear damped vibrations of the thin plate under consideration, the fractional derivative Kelvin-Voigt model is used, since this model has an advantage over the conventional Kelvin-Voigt model, because its prediction is in a good compliance with experimental data.Thus, for example, the experimental data obtained by Abdel-Ghaffar and Housner [17] and Abdel-Ghaffar and Scanlan [18] during ambient vibration studies of the Vincent-Thomas Suspension Bridge and the Golden Gate Bridge, respectively, show that different vibrational modes feature different amplitude damping factors, and the order of smallness of these coefficients highlights the low damping capacity of suspension combined systems, resulting in prolonged energy transfer from one partial subsystem to another.Besides, as natural frequencies of vibrations increase, the corresponding damping ratios decrease.To lead the theoretical investigations in line with the experiment, fractional derivatives were introduced by Rossikhin and Shitikova [19] for describing the processes of internal friction proceeding in suspension combined systems at nonlinear free vibrations, which allowed the authors to obtain the damping coefficients dependent on the natural frequency of vibrations.Moreover, the good agreement between the theoretical results and the experimental data has been found through the appropriate choice of the fractional parameter (the order of the fractional derivative) and the viscosity coefficient [20].

Problem Formulation and the Method of Solution
Let us consider the dynamic behavior of a free supported nonlinear thin rectangular plate (Figure 1), vibrations of which in a viscoelastic medium are described in the Cartesian system of coordinates by the following three differential equations written in the dimensionless form [2,3,21]: subjected to the initial u | =0 =  0 (, ) , V | =0 =  0 (, ) , where  = (, , ), V = V(, , ), and  = (, , ) are the displacements of points located in the plate's middle surface in the -, -, and -directions, respectively; ] is Poisson's ratio;  1 = / and  2 = ℎ/ are the parameters defining the dimensions of the plate;  and  are the plate's dimensions along the and -axes, respectively; ℎ is the thickness;  is the time; an overdot denotes the time-derivative; lower indices label the derivatives with respect to the corresponding coordinates;  0 (, ),  0 (, ), and  0 (, ) are functions describing the distribution of the initial velocities of the points lying in the middle surface of the plate;  is a small value.
In (1)- (6), the dimensionless values are introduced similarly as it has been done in [2]: where dimensional values are marked by asterisks, and  is the Young modulus of the plate material.
In ( 1)-( 3), the terms ae 1   , ae 2   V, and ae 3    represent the viscous resistance forces initiated when the plate vibrates in the surrounding viscous medium, where ae  ( = 1, 2, 3) are damping coefficients and   is the Riemann-Liouville fractional derivative of the -order [22]: But distinct to the traditional modeling the viscous resistance forces via first-order time-derivatives [23], in the present research we adopt the fractional order timederivatives   , which, as it will be shown below, will allow us to obtain the damping coefficients dependent on the natural frequency of vibrations.It has been demonstrated in [19,20] that such an approach for modeling the damped nonlinear vibrations of thin bodies provides the good agreement between the theoretical results and the experimental data through the appropriate choice of the fractional parameter (the order of the fractional derivative) and the viscosity coefficient.
It has been noted in [2,19] that a fractional derivative is the immediate extension of an ordinary derivative.In fact, when  → 1,    tends to ẋ ; that is, at  → 1 the fractional derivative goes over into the ordinary derivative, and the mathematical model of the viscoelastic plate under consideration transforms into the Kelvin-Voigt model, wherein the elastic element behaves nonlinearly, but the viscous element behaves linearly.When  → 0, the fractional derivative    tends to ().To put it otherwise, the introduction of the new fractional parameter along with the parameters ae  allows one to change not only the magnitude of viscosity at the cost of an increase or decrease in the parameters ae  , but also the character of viscosity at the sacrifice of variations in the fractional parameter.
Reference to (5) shows that free vibrations are excited by weak disturbance from the equilibrium position.
We seek the solution of (1)-(3) in the form where  1 (),  2 (), and  3 () are the generalized displacements corresponding to the displacements in the plane of the plate and to its deflection, respectively, but the natural functions satisfying the boundary conditions ( 6) have the form and  and  are integers.Linear undamped natural modes of flexural in-plane and out-of-plane vibrations of the plate are the solution of eigenvalue problem The set of equations ( 11) and ( 12) has the characteristic equation the roots of which are the dimensionless natural frequencies of the flexural in-plane vibrations of the plate where From ( 13), the natural frequency of the flexural out-ofplane vibrations of the plate could be obtained: Shock and Vibration Substituting ( 9) into ( 1)-(3), multiplying (1), (2), and (3) by  1 ,  2 , and  3 , respectively, integrating over  and , and using the orthogonality conditions for linear modes within the domains of 0 ≤ ,  ≤ 1, we are led to an infinite number of sets each involving a system of three coupled nonlinear ordinary differential equations of the second order in   ( = 1, 2, 3): where the summation is carried out over two repeated indices and the elements of the matrix    are defined as (16).The linear parts of nonlinear equations ( 18) and (19) have been studied in time and frequency domains in detail in [24] using the Laplace integral transform method as a method of solution.
The nonlinear parts   of ( 18) and ( 19) have the form where and the coefficients ( = 1, 2, . . ., 8) depending on the combinations of sine and cosine functions entering into the eigenfunctions (10) are presented in the appendix.
Since the second-rank tensor    is symmetric, then it has two real eigenvalues   defined in (15) which are in correspondence with two mutually orthogonal eigenvectors that is, where here and below summation is carried out over the repeated index , while there is no summation over the repeated indices  and .Thus, the matrix    and the generalized displacements   entering in ( 18) and ( 19) could be expanded in terms of the vectors (22) [25] as Substituting ( 24) and ( 25) in ( 18) and ( 19) and then multiplying (18) successively by  I  and  II  with due account for (23), we obtain the following three equations: where  3 =  3 .It should be emphasized that left-hand side parts of (26) are linear and independent of each other, while (26) are coupled only by nonlinear terms in their right-hand sides.
In order to show the influence of the initial conditions ( 4) and ( 5) on the solution to be constructed, let us expand the desired functions   ( = 1, 2, 3) in a series in terms of the small parameter : Substituting (27) in the set of equations ( 26) and restricting ourselves by the terms of the order of , we are led to a linear homogeneous set of differential equations The solution of (28) has the form where   () and   () are complex conjugate functions to be found.Representing functions   () in the polar form the solution ( 27) is reduced to where   () and   () are the amplitudes and phases of nonlinear vibrations, respectively.Differentiating (31) with respect to time  and ignoring the terms of the order of  2 , we obtain Now substituting (25) in (9) with due account for (31), we have Using the initial conditions (4) and ( 5) with due account for relationships (32), from (33) at  = 0 we obtain an infinite set of algebraic equations for determining   (0) and   (0) ( = 1, 2, 3): where All subsequent approximations are determined from an inhomogeneous set of differential equations with known right parts.Since the general solution of such a system is the sum of two solutions, a particular solution of the inhomogeneous system and a general solution of the corresponding homogeneous system, then arbitrary constants could be chosen in such a way that the initial conditions of all subsequent approximations would be zero ones.
It is known [14,[26][27][28] that during nonstationary excitation of thin bodies not all possible modes of vibration would be excited.Moreover, the modes which are strongly coupled by any of the so-called internal resonance conditions are initiated and dominate in the process of vibration, resulting in the energy transfer from one subsystem to another between the coupled modes; in so doing the types of modes to be excited are dependent on the character of the external excitation.As a result of damping, the modes that are not excited directly by an external excitation or through an internal resonance will die out, and only three modes will dominate the motion of the plate in three directions.
Assume hereafter that the vibration process occurs in such a way that only three natural modes corresponding to the generalized displacements  1 1  2 ,  2 1  2 , and  3 1  2 are excited and dominate over other natural modes.In this case, the right parts of (26) are significantly simplified, and equations of free vibrations (26) take the form where and the coefficients ( = 1, 2, . . ., 8) involving in (41)-(44) are presented in the appendix.

Method of Solution
An approximate solution of (45) for small but finite amplitudes weakly varying with time can be represented by a thirdorder uniform expansion in terms of different time scales in the following form [29]: where  = 1, 2, 3,  is a small dimensionless parameter of the same order of magnitude as the amplitudes, and   =    are new independent variables; among them,  0 =  is a fast scale characterizing motions with the natural frequencies and  1 =  and  2 =  2  are slow scales characterizing the modulation of the amplitudes and phases of the modes with nonlinearity.
Recall that the first and the second time-derivatives are defined, respectively, as follows: Further as a fractional derivative we will utilize the Gr ü nvald-Letnikov definition which is defined on a semiaxis [30] where is the th-order finite difference of the function () with the step ℎ and with the point  as a center and (   ) are the binomial coefficients which are determined by the formula Substituting (47) in (49), we obtain following Rossikhin and Shitikova [19] that Shock and Vibration 7 where   = /  and   0 ,  −1 0 ,  −2 0 , . . .are the Gr ü nvald-Letnikov fractional derivatives, which for "reasonably good" functions, such as exponential functions  ± , coincide with the Riemann-Liouville fractional derivatives in time  defined on a semiaxis [30].
Considering that viscosity of the medium surrounding the plate under consideration is small, that is, where   is the relaxation time of the th generalized displacement,   is a finite value, and the choice of  depends on the order of smallness of the viscosity coefficients ae  , substituting (46)-( 52) in (45), after equating the coefficients at like powers of  to zero, we are led to a set of recurrence equations to various orders: to order  3 In order to construct the uniformly valid solution, it is necessary on each step to use the solution from the preceding step and to eliminate secular terms during integration [29].
To solve the sets of ( 55) and ( 56)-( 58), it is necessary to specify the action of the fractional derivative   0 (8) on the functions  1 , that is, to calculate   0     .It has been shown in [31] that The first term in ( 60) is equivalent to the action of the fractional derivative with the low limit tending to −∞: the application of which for the exponential function is reduced to while the second term of (60), as it has been proved in [32], does not influence the solution constructed via the method of multiple time scales restricted to the zeroth-and first-order approximations.
In other words, even the utilization of exact formula (60) in the problem under consideration produces completely equivalent results given by the approximate formula (62) if the solution is constructed via the method of multiple time scales within the considered orders of approximation.Thus, in further analysis we will utilize formula (62).
For further analysis it is also a need to specify the order of weak damping.

Viscosity of the Order of 𝜀.
Let us first consider the case of viscosity of the order of .Then at  = 1 (55) are reduced to

Shock and Vibration
Substituting (59) in the right-hand side of (63) with due account for (62) yields where cc is the complex conjugate part to the preceding terms.
Reference to (64) shows that the following types of the internal resonance could occur on this step: (1) the two-to-one internal resonance, when one natural frequency is twice the other natural frequency, (2) the one-to-one-to-two internal resonance, that is, 3.2.Viscosity of the Order of  2 .Now let us consider vibrations of a nonlinear plate putting  = 2 in its equations of motion (55): Substituting (59) in the right-hand side of (68) yields Assume that the natural frequencies   ( = 1, 2, 3) possess such magnitudes that any of the combinations from (65) to (67) could not occur.Then the functions exp(±   0 ) entering into the right-hand side of (69) produce secular terms.
To eliminate circular terms in (69), it is necessary to vanish to zero the coefficients standing at exp(±   0 ), that is, whence it follows that   ( = 1, 2, 3) are  1 -independent.Considering (72), solution of (69)-( 71) at  = 1, 2, 3 has the form where   ( 2 ) ( = 1, 2, 3) are new functions to be determined and coefficients   ( = 1, 2, . . ., 8) have the form Substituting ( 59) and ( 73) with due account for (72) in ( 56)-( 58) at  = 2, we have Reference to (75) shows that the following types of the internal resonance could occur on this step: (1) the one-to-one internal resonance, when one natural frequency is nearly equal to the other natural frequency, (2) the one-to-one-to-one internal resonance (3) the combinational resonance of the additivedifference type, that is,

Governing Nonlinear Differential Equations Describing Amplitude-Phase Modulation for Different Types of the Internal Resonance of the Order of 𝜀
To deduce the nonlinear differential equations describing the modulation of amplitudes and phases of the nonlinear plate under consideration, we should consider separately each type of the internal resonance which could occur with due account for weak damping of the order .Then eliminating secular terms in (64), we obtain the following solvability equations: Reference to the set of equations ( 83)-(85) shows that its second equation is independent of other two, while the first and third ones represent a set of two nonlinear equations.The similar situation was noted in [2] for a fractionally damped nonlinear plate in the case of the two-to-one internal resonance, when equations describing plate's in-plane motion are coupled.Let us multiply (83)-(85), respectively, by  1 ,  2 , and  3 and find their complex conjugates.Adding every pair of the mutually adjoint equations with each other and subtracting one from another and considering that as a result we have where a dot denotes differentiation with respect to  1 ,  = 2 3 −  1 , and The set of nonlinear equations ( 87)-(88) was solved in [2].It has been found that during free vibrations of the plate being under the conditions of the two-to-one internal resonance three regimes can be observed: stationary (absence of damping at  = 0 and  = 0), quasi-stationary (damping is defined by an ordinary derivative at  = 1), and transient (damping is defined by a fractional derivative at 0 <  < 1).
Representing functions   entering in (89) in the polar form (86) and applying the same procedure as it has been carried out above for the 2:1 internal resonance, we have The nonlinear set of equations ( 90)-( 95) with the corresponding initial conditions completely describes the vibrational process of the mechanical system being investigated under the condition of the internal resonance 1:1:2 and could be solved numerically.This case has been discussed in [33].

Governing Nonlinear Differential Equations
Describing Amplitude-Phase Modulation for Different Types of the Internal Resonance of the Order of  2 To deduce the nonlinear differential equations describing the modulation of amplitudes and phases of the nonlinear plate under consideration, we should consider separately each type of the internal resonance which could occur with due account for weak damping of the order  2 .

Internal Resonance 1:
1.The one-to-one internal resonance could be of two types: when the natural frequencies of two modes of in-plane vibrations are close to each other, that is, when  1 =  2 = , or when the natural frequency of the out-of-plane mode is close to the natural frequency of one of the in-plane modes, that is,  3 =  1 =  or  3 =  2 = .

Internal
Resonance  1 =  2 = .Now let us first consider case (76) of the one-to-one internal resonance, when  1 =  2 = , but  3 ̸ = .Then eliminating secular terms in (75), we obtain the following solvability equations: Let us multiply (98)-(100), respectively, by  1 ,  2 , and  3 and find their complex conjugates.Adding every pair of the mutually adjoint equations with each other and subtracting one from another and considering that the functions   ( 2 ) could be represented in the polar form (86), where   and   ( = 1, 2, 3) are  2 -dependent, as a result we have where the phase difference has the form  =  2 −  1 and an overdot denotes the differentiation with respect to  2 .
It could be shown that the expression in square brackets in (103) and ( 106) is vanished to zero when  1 =  2 = .Then (103) turns out to be uncoupled with respect to all other equations from (101) to (106): the solution of which has the form where  3 is a constant of integration to be determined from the initial conditions, while (106) is reduced to whence it follows that the phase  3 is dependent only on the squared values of all amplitudes of the three interacting modes of vibrations.

Shock and Vibration
Considering (110) and ( 113), (112) could be reduced to where  =  −1 ( 2  6 −  1  5 ) 3 .Equations ( 101) and (112) could be rewritten in another form if we substitute (110) and (108) in all its terms: where The nonlinear set of equations ( 111), (115), and (116), or (114), with the initial conditions completely describes the vibrational process of the mechanical system being investigated under the condition of the oneto-one internal resonance  1 =  2 and could be solved numerically.
It should be noted that the equation defining the stream function could be obtained from (114), if we neglect the first term in (114) which decays rather rapidly as compared with its last term.
Particular Case.In the particular case at Σ = 0 and  1 =  2 =  3 = , (111) has the form whence it follows that where  0 is the initial energy.Considering (120), ( 115) and (116), wherein the term rapidly decaying with time has been neglected, are reduced to the following: while (114), wherein its first term rapidly decaying with time has been neglected, takes the form and it could be integrated, resulting in its first integral The second first integral (124) defines the stream function (, ) such that which describes steady-state vibrations of an elastic plate decaying with time.
Eliminating the variable  from ( 121) and ( 124) and integrating over  2 , we have The integral in the left-hand side of (126) can be transformed to an incomplete elliptic integral of the first kind and can be easily calculated using special tables [34].
The solution of (126) allows one to find the value ( 2 ) and, thus, to solve the problem under consideration.
If we represent δ considering (121) as and retain all terms in (116) with due account for (120), then (116) could be rewritten as follows: Integrating (128) yields Note that relationships (125) and (126) are valid in this case as well.

Internal Resonance
= , could be treated similarly.Then eliminating secular terms in (75), we obtain the following solvability equations:

Shock and Vibration 13
Representing functions   entering in (130) in the polar form and applying the same procedure as it has been carried out above for the other case of the 1:1 internal resonance, as a result we have where the phase difference has the form  = 2( 3 −  1 ).Since (132) is independent of others from (131) to (136), its solution has the form where  2 is a constant of integration to be determined from the initial conditions.Introducing new functions  1 ( 2 ) and  3 ( 2 ) such that and adding (131) and (133) with due account for (138) yield while subtracting (134) from ( 136) we obtain where From ( 131) and (133) we could find that Considering ( 138) and ( 141), ( 131) and (140) could be reduced to where The nonlinear set of equations ( 139) and ( 142) with the initial conditions (117) completely describes the vibrational process of the mechanical system being investigated under the condition of the one-to-one internal resonance  1 =  3 and could be solved numerically.
Particular Case.In the particular case at Σ = 0 and  1 =  2 =  3 = , (139) has the form whence it follows that Considering ( 146), (142) are reduced to while (148), wherein its first term rapidly decaying with time has been neglected, takes the form and it could be integrated, resulting in its first integral This first integral (150) defines the stream function (, ) such that which describes steady-state vibrations of an elastic plate decaying with time.
Eliminating the variable  from (147) and (150) and integrating over  2 , we have The solution of (152) allows one to find the value ( 2 ) and, thus, to solve the problem under consideration.
If we represent δ considering (147) as and retain all terms in (148) with due account for (146), then (148) could be rewritten as follows: Integrating (154) yields Note that relationships (151) and (152) are valid in this case as well.
Thus, during free vibrations of the plate being under the conditions of the one-to-one internal resonance three regimes can be observed: stationary (absence of damping at  = 0), quasi-stationary (damping is defined by an ordinary derivative at  = 1), and transient (damping is defined by a fractional derivative at 0 <  < 1).
Then eliminating secular terms in (75), we obtain the following solvability equations: Representing functions   entering in (151) in the polar form and applying the same procedure as it has been carried out above for the cases of the 1:1 internal resonance, as a result we have The nonlinear set of equations ( 157) with the initial conditions (117) completely describes the vibrational process of the mechanical system being investigated under the condition of the internal resonance 1:1:1 and could be solved numerically.

Combinational Resonances of the Additive-Difference
Type.Reference to (75) shows that three cases of the combinational internal resonance of the additive-difference type (80)-(82) could occur for the plate under consideration.

Combinational Resonance: 2𝜔
Then eliminating secular terms in (75), we obtain the following solvability equations: Let us multiply (158)-(160), respectively, by  1 ,  2 , and  3 and find their complex conjugates.Adding every pair of the mutually adjoint equations with each other and subtracting one from another and considering that the functions   ( 2 ) could be represented in the polar form (86), where   and   ( = 1, 2, 3) are  2 -dependent, as a result we have where the phase difference has the form  = 2 3 −  2 −  1 .

Shock and Vibration
Then (161), (163), and (169) could be rewritten in the following form: where ) . ( The nonlinear set of equations ( 168), ( 171), (172), and (173) with the initial conditions (117) completely describes the vibrational process of the mechanical system being investigated under the condition of the combinational internal resonance 2 3 =  1 +  2 and could be solved numerically.This case has been discussed in [35].
Particular Case.In the particular case at Σ = 0 and  1 =  2 =  3 = , (168) has the form whence it follows that where  1 ,  2 , and  3 are constants of integration such that Considering ( 177), ( 171) and ( 173) are reduced to while ( 180), wherein its first term rapidly decaying with time has been neglected, takes the form and it could be integrated, resulting in its first integral This first integral (182) defines the stream function (, ) such that which describes steady-state vibrations of an elastic plate decaying with time.
Eliminating the variable  from ( 179) and ( 182) and integrating over  2 , we have The solution of (184) allows one to find the value ( 2 ) and, thus, to solve the problem under consideration.
If we represent δ considering (179) as and retain all terms in (180) with due account for (177), then (180) could be rewritten as follows: Integrating (186) yields Note that relationships (183) and ( 184) are valid in this case as well.

Combinational Resonance: 2𝜔
Then eliminating secular terms in (75), we obtain the following solvability equations: Applying to (188) the same procedure as it has been done above for (158)-( 160), as a result we have where the phase difference has the form Introducing new functions  1 ( 2 ),  2 ( 2 ), and  3 ( 2 ) such that and adding ( 189)-( 191) with due account for (195) yield while subtracting (192) from the sum of (193) and the doubled (194) we obtain where From ( 189)-(191) we could find that

Shock and Vibration
Then (189), (191), and (197) could be rewritten in the following form: where ) . ( The nonlinear set of equations ( 196), (199), (200), and (201) with the initial conditions (117) completely describes the vibrational process of the mechanical system being investigated under the condition of the combinational internal resonance 2 3 =  1 −  2 and could be solved numerically.
Particular Case.In the particular case at Σ = 0 and  1 =  2 =  3 = , (196) has the form whence it follows that where  1 ,  2 , and  3 are constants of integration such that Considering ( 205), (199) and ( 201) are reduced to while (208), wherein its first term rapidly decaying with time has been neglected, takes the form and it could be integrated, resulting in its first integral This first integral (210) defines the stream function (, ) such that which describes steady-state vibrations of an elastic plate decaying with time.
Eliminating the variable  from (207) and (210) and integrating over  2 , we have The solution of (212) allows one to find the value ( 2 ) and, thus, to solve the problem under consideration.
If we represent δ considering (207) as and retain all terms in (208) with due account for (205), then (208) could be rewritten as follows: Integrating (214) yields Note that relationships (211) and ( 212) are valid in this case as well.
Then eliminating secular terms in (75), we obtain the following solvability equations: reducing to six differential equations in amplitudes and phases of the three interacting modes: where the phase difference has the form Introducing new functions  1 ( 2 ),  2 ( 2 ), and  3 ( 2 ) such that and adding ( 217)-( 219) with due account for (223) yield while subtracting (221) from the sum of (220) and the doubled (222) we obtain where From (217)-(219) we could find that Then ( 217), (219), and (225) could be rewritten in the following form: where The nonlinear set of equations ( 224), ( 227), (228), and (229) with the initial conditions (117) completely describes the vibrational process of the mechanical system being investigated under the condition of the combinational internal resonance 2 3 =  2 −  1 and could be solved numerically.
Particular Case.In the particular case at Σ = 0 and  1 =  2 =  3 = , (224) has the form whence it follows that where  1 ,  2 , and  3 are constants of integration such that Considering ( 233), ( 227) and ( 229) are reduced to while ( 236), wherein its first term rapidly decaying with time has been neglected, takes the form and it could be integrated, resulting in its first integral This first integral (238) defines the stream function (, ) such that which describes steady-state vibrations of an elastic plate decaying with time.
Eliminating the variable  from (235) and (238) and integrating over  2 , we have The solution of (240) allows one to find the value ( 2 ) and, thus, to solve the problem under consideration.
If we represent δ considering (235) as and retain all terms in (236) with due account for (233), then (236) could be rewritten as follows: Integrating (242) yields Note that relationships (239) and ( 240) are valid in this case as well.

Numerical Investigations
As examples, let us carry out the qualitative analysis of the cases of one-to-one internal resonance (76), when two different modes of in-plane vibrations are coupled, and (77), when one mode of in-plane vibrations is coupled with a certain mode of out-of-plane vibrations.

One-to-One Internal
Resonance  1 =  2 (76).For this case, the stream function (, ) is defined by relationship (129), and the phase portrait to be constructed according to (129) depends essentially on the magnitudes of the coefficient .
6. 1.1. Case (76) at  = 0.For this case, the stream function (, ) is constructed according to (124) and the streamlines of the phase fluid in the phase plane - are presented in Figure 2. Magnitudes of  are indicated by digits near the curves which correspond to the streamlines; the flow directions of the phase fluid elements are shown by arrows on the streamlines.
Reference to Figure 2 shows that the phase fluid flows within the circulation zones, which tend to be located around the perimeter of the rectangles bounded by the lines  = 0,  = 1, and  = ±(/2) ± 2 ( = 0, 1, 2, . ..).As this takes place, the flow in each such rectangle becomes isolated.On all four rectangle sides  = 0 and inside it the value  preserves its sign.The function  attains its extreme magnitudes ±0.5 at the points with the coordinates  = 1/2,  = ± ( = 0, 1, 2, . ..).
Along the lines  = ±(/2) ± 2 ( = 0, 1, 2, . ..) pure amplitude modulated aperiodic motions are realized, since with an increase in time  the value  increases from  0 to 1 (along the line  = −/2) or decreases from  0 to 0 (along the line  = /2), and solution (126) could be written in the form of a soliton-like solution: Streamlines give a pictorial estimate of the connection of  with all types of the energy-exchange mechanism.Thus, in the case of undamped vibrations, that is, when the damping coefficient is equal to zero and  = 0, the points with the coordinates  0 = 1/2,  0 = ± ( = 0, 1, 2, . ..) correspond to the stationary regime, since ξ = 0 and δ = 0 according to (125).The stationary points  0 = 1/2,  0 = ± are centers; as with a small deviation from a center, a phase element begins to move around the stationary point along a closed trajectory.Closed streamlines correspond to the periodic change of both amplitudes and phases.
The  2 -dependence of the values  and 1 −  which are proportional to the square of the amplitudes  1 and  2 , respectively, at  0 = /4 and  0 = 0.5 is presented in Figure 2, wherein the solid and dashed lines correspond, respectively, to the cases of damped and undamped vibrations.Reference to Figure 3 shows the damping of the energy exchange between two subsystems.Moreover, the points of touching of solid curves, that is, the points of tangency of the envelopes of vibrations (  = 1 −   ,  2  ), allow one to determine the value where   are defined by (110) with account for (120).
As in the previous case, the phase fluid flows in an infinitely long channel, the boundaries of which are the straight lines  = 0 and  = 1, corresponding to the phase modulated motions.In one part the streamlines are nonclosed, which corresponds to the periodic change of amplitudes and the aperiodic change of phases; in another part they are closed, which corresponds to the periodic change of both amplitudes and phases.The aperiodic regime lines are the boundaries of the closed and unclosed streamline areas.From the phase portrait in Figure 4 it is seen that the circulation zones are located in a staggered arrangement by the right and left channel sides (this configuration resembles that of von Kármán staggered vortex tracks).
Each zone by the side  = 1 is surrounded by a line with the value  = −0.1.However, only those parts of the line  = −0.1 which bound these zones from below and come closer to the side  = 1 at the points  = 1,  = /2 ±  belong to the domain of the fluid flow.The transition of fluid elements from the points  = 1,  = (3/2) ±  to the points  = 0,  = (/2) ±  proceeds instantly.The soliton-like solution along the separatrix is written as wherein  varies from  min = 0.0099 to  max = 1.
The circulation zones by the side  = 0 are surrounded by the line with the value  = 0.However, only those parts of the line  = 0 which bound these zones from above and come closer to the side  = 0 at the points  = 0,  = /2 ±  belong to the domain of the fluid flow.The transition of fluid elements from the points  = 0,  = (/2) ±  to the points  = 0,  = (−/2) ±  proceeds instantly.The soliton-like solution along the separatrix line  = 0 has the form wherein  varies from  min = 0 to  max = 0.99.Inside the two circulation zones there are points with the extreme values of the stream function: maximal  max = 0.4525 and minimal  min = −0.5525,respectively.These points are the centers corresponding to the stable stationary regimes  =  0 = 0.45,  =  0 = ±2 and  =  0 = 0.55,  =  0 =  ± 2, respectively.
Between the lines corresponding to  = 0 and  = −0.1,unclosed streamlines are located which are in accordance with the periodic change of the amplitudes and the aperiodic change of the phase difference.76) at (1/2) −1 = ±1.In this case, the stream function (, ) (129) is defined as

Case (
and Figures 5 and 6 show the streamlines of the phase fluid in the phase plane, respectively.From the comparison of Figures 5 and 6 it is evident that with the increase of the parameter  the zone of the unclosed streamlines becomes wider while the circulation zones locating in a staggered arrangement by the right and left channel sides  = 0 and  = 1 turn out to be more narrow.The soliton-like solution along the separatrix line  = 0 near the bank  = 0 is found to be as follows: wherein  varies from  min = 0 to  max = 0.5, while along the separatrix  = −1 attached to the side  = 1 it has the form wherein  varies from  min = 0.5 to  max = 1.
Inside the circulation zones there are the stationary points, that is, the centers corresponding to the stable stationary regimes, where the stream function takes the extreme values:  max = 0.207 at the points with the coordinates  =  0 = 0.1465,  =  0 = ±2 and  min = −1.207at the points with the coordinates  =  0 = 0.8535,  =  0 =  ± 2.
At the negative value of parameter  −1 , the phase fluid flow portrait shown in Figure 5 is shifted along the -axis on the value , the flow is opposite in direction (see Figure 6), and the maxima and minima change places.In so doing the solutions along the curve separatrices with  = 0 and  = 1 are written according to (251) and (252), respectively.

One-to-One Internal
Resonance  1 =  3 (77).For this case, the stream function (, ) is defined by relationship (155), and the phase portrait to be constructed according to (155) depends essentially on the magnitudes of the coefficients  2 and  3 .Reference to Figure 7 shows that the phase fluid flows within the circulation zones, which tend to be located around the perimeter of the rectangles bounded by the lines  = 0,  = 1, and  = ±(/2) ± 2 ( = 0, 1, 2, . ..).As this takes place, the flow in each such rectangle becomes isolated.On all four rectangle sides  = 0 and inside it the value  preserves its sign.The function  attains its extreme magnitudes ±0.25 at the points with the coordinates  = 1/2,  = ± ( = 0, 1, 2, . ..).

Case (77) at
Along the boundaries  = 0 and  = 1 the phase modulated regimes are realized because the particles of the phase fluid flow with finite velocities (as distinct to the similar case shown in Figure 2 for  1 =  2 , where velocities of fluid particles tend to infinity as  → 0 and  → 1), while along the lines  = ±(/2) ± 2 ( = 0, 1, 2, . ..) there exist pure amplitude modulated aperiodic motions, since with an increase in time  the value  increases from  0 to 1 (along the line  = −/2) or decreases from  0 to 0 (along the line  = /2), and the solution (152) could be written in the form of a soliton-like solution: 6.2.2.Case (77) at  2 ( 0 ) −1 = −1 and  3  −1 = 1/2.In this particular case, the stream function (, ) (150) is defined as and Figure 8 shows the streamlines of the phase fluid in the phase plane.
As in the previous case, the phase fluid flows in an infinitely long channel, the boundaries of which are the straight lines  = 0 and  = 1, corresponding to the phase modulated motions.In one part the streamlines are nonclosed, which corresponds to the periodic change of amplitudes and the aperiodic change of phases; in another part they are closed, which corresponds to the periodic change of both amplitudes and phases.The aperiodic regime lines are the boundaries of the closed and unclosed streamline areas.From the phase portrait in Figure 8 it is seen that the circulation zones are located in a staggered arrangement by the right and left channel sides (this configuration resembles that of von Kármán staggered vortex tracks).
Each zone by the side  = 1 is surrounded by a line with the value  = 0.5.This line consists of two parts connected with each other at the branch points with the coordinates  = 1,  = /2 ±  ( = 0, 1, 2, . ..).One branch of this line corresponds to the phase modulated regime  = 1, and the other corresponds to the aperiodic regime, wherein  varies from  min = 0.2 to  max = 1.At the branch point itself, the phase fluid flow velocity is equal to zero.The soliton-like solution along the separatrix with  = 0.5 is written as The circulation zones by the side  = 0 are surrounded by the line with the value  = 0.25.This line consists of two parts connected with each other at the branch points with the coordinates  = 0,  = 2/3± ( = 0, 1, 2, . ..).One branch of this line corresponds to the phase modulated regime  = 0, and the other corresponds to the aperiodic regime, wherein  varies from  min = 0 to  max = 2/3.At the branch point itself, the phase fluid flow velocity is equal to zero.The soliton-like solution along the separatrix line  = 0.25 has the form Inside the two circulation zones there are points with the extreme values of the stream function: maximal  max = 0.7 and minimal  min = 1/6, respectively.These points are the centers corresponding to the stable stationary regimes  =  0 = 0.6,  =  0 = ±2 and  =  0 = 1/3,  =  0 =  ± 2, respectively.
Between the lines corresponding to  = 0.25 and  = 0.5, unclosed streamlines are located which are in accordance with the periodic change of the amplitudes and the aperiodic change of the phase difference.and Figure 9 shows the streamlines of the phase fluid in the phase plane.
From Figure 9 it is seen that, unlike the previous case presented in Figure 8, the circulation zones by the side  = 0 and the aperiodic regime disappear.If  → 0, then the streamlines level off and tend to the line  = 0 with  =  min = −0.5,where the phase modulated regime is realized.
Each circulation zone by the side  = 1 is surrounded by a line with the value  = 0.5.This line consists of two parts connected with each other at the branch points with the coordinates  = 1,  = /2 ±  ( = 0, 1, 2, . ..).One branch Inside the circulation zones there are points with the maximal value  max = 0.625.These points are the centers corresponding to the stable stationary regimes  =  0 = 0.75,  =  0 = ±2.The lines corresponding to −0.5 <  < 0.5 are unclosed streamlines which are in accordance with the periodic change of the amplitudes and the aperiodic change of the phase difference.
and Figure 10 shows the streamlines of the phase fluid in the phase plane.Figure 10 illustrates the phase portrait with only unclosed phase fluid streamlines along which the fluid flows in the direction of a decrease in .With  → 0 and  → 1, the streamlines level off and tend, respectively, to the lines  = 0 with  =  max = 0.5 and  = 1 with  =  min = −0.5,where the boundary phase modulated regimes are realized.Reference to Figure 11 shows that the phase fluid flows within the circulation zones, which tend to be located around the perimeter of the rectangles bounded by the lines  = 0,  = 1, and  = ±()±2 ( = 0, 1, 2, . ..).As this takes place, the flow in each such rectangle becomes isolated.On all four rectangle sides  = −0.5, and inside it the value −0.5 <  < 0 attaining its maximal magnitude  max = 0 at the points with the coordinates  = 1/2,  = ±2 ( = 0, 1, 2, . ..).

Conclusion
Free damped vibrations of a nonlinear plate in a fractional derivative viscoelastic medium have been investigated.Nonlinear vibrations are described by coupled three equations with respect to the three orthogonal displacements.
The proposed analytical approach for investigating the damped vibrations of the nonlinear plate subjected to the conditions of the internal resonance has been possible owing to the new procedure suggested in this paper, resulting in decoupling linear parts of equations with further utilization of the method of multiple scales for solving nonlinear governing equations of motion.
It has been shown that the phenomenon of internal resonance could be very critical, since in the thin plate under consideration the internal resonance is always present.Moreover, its type depends on the order of smallness of the viscosity involved into consideration.Thus, at the order, damped vibrations occur within the two-to-one and Shock and Vibration one-to-one-to-two internal resonance.Other types of the internal resonance, such as one-to-one, one-to-one-to-one, and combinational resonances of the additive and difference types could be found at  2 -order; that is, the type of the resonance depends on the order of smallness of the fractional derivative entering in the equations of motion of the plate.
For each type of the resonance, the nonlinear sets of resolving equations in terms of amplitudes and phase differences have been obtained.It has been shown that for the one-to-one and combinational internal resonances there exist such particular cases when it is possible to obtain two first integrals, namely, the energy integral and the stream function, which allows one to reduce the problem to the calculation of elliptic integrals.
The influence of viscosity on the energy exchange mechanism has been analyzed.It has been shown that each mode is characterized by its damping coefficient connected with the natural frequency by the exponential relationship with a negative fractional exponent.Thus, during free vibrations of the plate with internal resonances three regimes can be observed: stationary (absence of damping at  = 0), quasistationary (damping is defined by an ordinary derivative at  = 1), and transient (damping is defined by a fractional derivative at 0 <  < 1).
The phenomenological analysis carried out for two types of the one-to-one internal resonance with the help of the phase portraits constructed for different magnitudes of the plate parameters reveals the great variety of vibrational motions: stationary vibrations, two-sided energy exchange between two subsystems under consideration, and onesided energy interchange resulting in the complete one-sided energy transfer.Under the presence of small viscosity all regimes attenuate with time.
The analysis of the phase portraits for various oscillatory regimes cited above demonstrates that they contain closed as well as nonclosed streamlines, which are separated by the curves separatrixes.Along the separatrixes one succeeds in finding analytic solutions which determine the irreversible energy exchange from one subsystem into another and which are inherently soliton-like solutions in the theory of vibrations.
Comparison of the results obtained in this paper via the new suggested procedure with those for a nonlinear plate, the motion of which is described either by one equation in terms of plate's deflection [6] or also by three coupled nonlinear equations in terms of three displacements but without decoupling the linear parts of governing equations [2,3], reveals the fact that the new approach allows one to discover and analyze much more diversified variety of internal resonances, including the one-to-one-to-one and oneto-one-to-two internal resonances and additional one-toone, two-to-one, and combinational resonances (ten possible cases in all), than all other existing previously in literature approaches.
It should be noted that the proposed analytical approach could be used for other cases if the matrix of the natural modes coupled by the internal resonance is Hermitian matrix; that is, it should be symmetric and all its elements should be real.

Figure 1 :
Figure 1: Scheme of a freely supported rectangular plate.

2 Figure 3 :Figure 4 :
Figure 3: The  * 2 -dependence of the dimensionless amplitudes of vibrations for the case of 1:1 internal resonance.