Dynamic Stability of Orthotropic Viscoelastic Rectangular Plate of an Arbitrarily Varying Thickness

Featured Application: The results obtained apply to dynamic analyses of thin-walled plates and shells of variable thickness. Abstract: The research object of this work is an orthotropic viscoelastic plate with an arbitrarily varying thickness. The plate was subjected to dynamic periodic load. Within the Kirchhoff–Love hypothesis framework, a mathematical model was built in a geometrically nonlinear formulation, taking into account the tangential forces of inertia. The Bubnov–Galerkin method, based on a polynomial approximation of the deﬂection and displacement, was used. The problem was reduced to solving systems of nonlinear integrodifferential equations. The solution of the system was obtained for an arbitrarily varying thickness of the plate. With a weakly singular Koltunov–Rzhanitsyn kernel with variable coefﬁcients, the resulting system was solved by a numerical method based on quadrature formulas. The computational algorithm was developed and implemented in the Delphi algorithmic language. The plate’s dynamic stability was investigated depending on the plate’s geometric parameters and viscoelastic and inhomogeneous material properties. It was found that the results of the viscoelastic problem obtained using the exponential relaxation kernel almost coincide with the results of the elastic problem. Using the Koltunov–Rzhanitsyn kernel, the differences between elastic and viscoelastic problems are signiﬁcant and amount to more than 40%. The proposed method can be used for various viscoelastic thin-walled structures such as plates, panels, and shells of variable thickness.


Introduction
Thin-walled plates and shells of variable thickness are widely used in various technology fields due to the requirements for strength, durability, and thin-walled structural elements' exterior appearance. The study of plates and shells of variable thickness made from homogeneous and inhomogeneous materials is a challenging task, and sometimes it encounters insurmountable difficulties. The difficulties are due to the solution of rather cumbersome equations, which are obtained by the desire to reflect, in mathematical modeling, the real mechanical behavior of the thin-walled structure. Variable thickness plate problems involve computational complexity, i.e., the lack of suitable universal numerical methods for solving the obtained equations, and as a consequence, unified computational algorithms.
Studies of parametric vibrations of thin-walled structures have become a separate research area in the mechanics of deformable solids. Thin-walled structures have versatile applications for various mechanical systems, in particular for plates and shells. The solution of dynamic stability problems includes the derivation of the equation of motion, discretization, and determination of areas of dynamic instability of structures.
Several publications [1][2][3] were devoted to studying the behavior of plates and shells of constant thickness under dynamic loads in an elastic formulation, where detailed reviews of these studies can be found. The research of the dynamic stability of elastic systems was started in 1964 [1]. Bolotin's first-order approximation followed by an eigenvalue analysis gave the dynamic behavior including natural frequencies and instability load-frequency margins. The stability of the plane form of the equilibrium position of plates has been studied thoroughly in [2] for the case of compression, typically, in a uniaxial stress field, where it was assumed that the shape of buckling of the plate does not depend on buckling position's coordinate. The article [3] is devoted to reviewing studies carried out from 1987 to 2005 in dynamic stability and parametric vibrations of plate and shell structures.
A new approach [4] based on the Amabili-Reddy higher-order shear deformation theory was implemented to modelling the geometrically nonlinear forced vibrations of laminated circular cylindrical shells. The numerical solution of discretized motion equations showed the one-to-one internal resonance for a complete circular cylindrical shell, giving rise to pitchfork bifurcations of the nonlinear response with quasiperiodic vibrations.
Theoretical and experimental aspects of nonlinear vibration and stability of shells and plates when the vibration amplitude is larger than the thickness are explained in [5]. A very good agreement between experiments and theoretical predictions was obtained for rectangular plates and circular cylindrical panels in the cases considered.
The article [6] provides an overview of publications devoted to studying nonlinear free and forced vibrations of shells made of traditional and modern materials published from 2003 to 2013. In particular, the results of studies of parametric oscillations were considered.
Nonlinear vibrations of viscoelastic thin rectangular plates were investigated in [7]. The plates were subjected to harmonic excitation. The plates' material was modeled as a Kelvin-Voigt viscoelastic solid by retaining all the nonlinear terms. It was shown that the nonlinear response of viscoelastic plates is different from the response of the same plate with viscous damping.
Experiments presented in [8] have shown that the most significant nonlinear effect of geometrically nonlinear vibrations of rectangular plates is the increase in damping with the vibration's amplitude. At the same time, soft materials present an increase in their stiffness with the vibration frequency. These two phenomena appearing together were explained in the framework of viscoelasticity. The nonlinear damping model for a continuous system was presented in [8]. The plate's material was assumed to be linearly viscoelastic and modeled by the classical standard linear solid. Then, the equations of motion describing the rectangular plates' nonlinear vibrations were derived by Lagrange equations. Rectangular plates of constant thickness with nonlinear fractional damping were investigated in [9]. It was underlined that damping increases with the vibration amplitude during nonlinear vibrations and, at the same time, soft materials such as silicone rubber increase their stiffness with the vibration frequency. The fractional linear solid model applies to describe the viscoelastic material behavior. The comparison of numerical and experimental results showed a good match.
The experimental part of the research in [10] was similar to [9]. A square silicone plate of constant thickness made of silicone with fixed edges was experimentally tested. The plate was excited with harmonic force with a frequency lower than the fundamental natural frequency. Linear and nonlinear plate responses were measured by automatic scanning laser Doppler vibrometer. The increase in damping was around 60% when the vibration amplitude was 1.6 times the plate thickness.
Out-of-plane thin membrane actuators have a vital role in micro-electromechanical systems. In [11], the nonlinear vibration control and optimal locations of actuators for a pre-tensioned membrane structure were investigated based on the nonlinear theory of plates. The theoretical and approximate solutions of nonlinear frequencies of the membrane were obtained for a complex parametric study of the behavior of membranes, and the discrepancies between the two solutions under different vibration amplitudes were discussed.
The dynamic instability of thin laminated composite rectangular plates with various length-width and length-thickness aspect ratios subjected to harmonic in-plane loading has been theoretically studied [12]. The equations of motion of the plate were developed using the von Kármán type of plate equation considering geometric nonlinearity. The equations were solved by using Galerkin's technique. Dynamically unstable regions and both stable-and unstable-solution amplitudes of the steady-state vibrations were obtained by applying Bolotin's method.
Dynamic instability and nonlinear parametric vibrations of multilayer composite plates subjected to the action of periodic compressive loads in the plane applied along two opposite edges were considered in [13]. A system of ordinary differential nonlinear equations of the second-order Mathieu-Hill with periodic coefficients was obtained for determining the regions of dynamic instability and nonlinear parametric oscillations based on the Bolotin method. The influence of various parameters on the dynamic stability and nonlinear parametric oscillations was investigated: static load, dynamic load, elongation, boundary conditions, orthotropy, and damping coefficient.
Functionally graded materials, the material properties of which are varied along the thickness direction, are widely used in aerospace structures, nuclear reactors, and other structures with high functional performance [14]. The nonlinear dynamic instability of a functionally graded rectangular plate of constant thickness was investigated [15]. The plate was made by mixing two different materials, metal on one surface and ceramic on the other surface. The plate was subjected to nonuniform loading. The equations of motion used higher-order shear deformation theory considering von Kármán geometric nonlinearity. Galerkin's method was used to reduce the nonlinear governing partial differential equations into ordinary differential equations describing the plate's nonlinear dynamic instability. Bolotin's method was used to obtain the boundaries of dynamic instability regions. Effects of power index, span-to-thickness ratio, aspect ratio, static load factor, boundary conditions, and various types of linearly varying loadings and parabolically distributed loadings were studied.
Analytical studies on the dynamic instability analysis of a functionally graded constant thickness skew plate (in the form of the parallelogram) made by mixing ceramic and metal materials were presented in [16]. The plate was subjected to uniform and linearly varying inplane periodic loadings with four different types of boundary conditions. The total energy functional of the plate was formulated based on Reddy's third-order shear deformation theory. The energy functional was converted into a set of ordinary differential equations (Mathieu-Hill equations) using the Rayleigh-Ritz method in conjunction with BCOPs. The solution of the Mathieu-Hill equations described the dynamic instability behavior of the skew plate. The instability regions were traced using the Bolotin method.
The parametric instability of a functionally graded polymer-based graphene-reinforced nanocomposite plate of constant thickness was investigated in [17]. The plate comprised multiple graphene platelet reinforced composite layers in which graphene platelets are uniformly distributed in each layer with graphene platelet concentration varying layerwise across the plate thickness. The plate was subjected to periodic uniaxial in-plane force and a uniform temperature rise. The basic equations were derived from the theory of first-order shear deformation. The problems were solved using the differential quadrature approach in combination with the Bolotin method.
The instability of rectangular plates due to parametric resonance was considered in [18]. The partial differential equation was derived using the extended Hamilton principle. Using the Galerkin method, the resulting equation was reduced to a system of ordinary differential equations. Semianalytical results were obtained for various plate boundary conditions, confirmed by numerical simulations and other published results.
The boundaries of instability of axially accelerating plates with internal resonance were studied in [19]. A connection was introduced between acceleration and tensions changing in the longitudinal direction. The main equation and the corresponding boundary conditions were derived from the generalized Hamilton's principle. The effects of internal resonances and inhomogeneous boundary conditions at the boundaries of instability were distinguished. Multiple scales established modified conditions of solvability in the main parametric and internal resonances. The Routh-Hurwitz criterion was introduced to determine the boundaries of instability. The influence of the coefficient of viscoelasticity and the coefficient of viscous damping on instability boundaries was investigated. Abnormal instability boundaries were detected upon the application of the internal resonance. The phenomenon of local zigzag and V-shaped boundaries was explained in terms of modal interactions. Numerical calculations of differential quadrature circuits for the first four complex frequencies, the first four complex modes, and stability boundaries were used to confirm the analytical method's results.
The R-functions theory was used for solving dynamic problems of the theory of plates and shells, connected with the construction of the basic function in the case of an arbitrary domain [20][21][22][23][24][25][26]. In [20], dynamic instability and nonlinear parametric vibrations of multilayer plates of complex shapes were investigated. A numerical-analytical method based on a combination of the theory of R-functions and the variational method was proposed. Instability regions and response curves were plotted for various types of materials and load parameters.
In [21], using the theory of R-functions, a method was proposed for studying dynamical instability and nonlinear parametric vibrations of laminated plates. The plates had complex shapes and different cutouts. The structure of the plate was symmetrical concerning its middle surface, consisting of layers of constant thickness. The loads consisted of a periodic in-plane load and static load. The mathematical setting was carried out within the framework of the classical theory based on the Kirchhoff-Love hypothesis. The proposed method was based on the first-order shear deformation theory and the classical plate theory. A plane stress analysis was carried out using the variational Ritz method and the R-functions theory. The obtained results were applied to investigate buckling and parametric vibrations of laminated plates.
In [22], the problem of nonlinear vibrations and analysis of the stability of multilayer symmetric plates of complex shape, loaded with a static or periodic load in the plane, was considered. An approach based on the application of the theory of R-functions and variational methods was proposed. The proposed approach was demonstrated in testing tasks and applied to laminated plates with cutouts. The influence of geometrical parameters, load, boundary conditions on stability regions, and nonlinear oscillations was investigated.
The paper [23] proposes a method for studying parametric vibrations of orthotropic plates of complex shape for various types of boundary conditions, based on variational methods in combination with the theory of R-functions. Specific tasks were solved. Regions of dynamic instability and deflection depending on the frequency and excitation time were plotted for a plate of complex shape.
R-functions theory was used to study the free vibrations and dynamic instability of a symmetrically layered plate [24]. The plate was subjected to combined static and periodic axial forces. The problem was considered within the framework of the classical theory of plates. The solution of problems was based on the variational Ritz method's joint application, the Galerkin procedure, the theory of R-functions, and Bolotin's theory. The influence of the plate's geometric parameters and the mechanical characteristics of the material on the parametric resonance was studied.
The combination of the R-function theory and variational methods was proposed in [25][26][27]. The parametric vibrations of T-shape, rectangular, and rectangular cutout orthotropic plates of constant thickness were investigated for various boundary conditions [25,27]. Using the Kirchhoff-Love hypothesis, the motion equations were built. The proposed hybrid approach is based on the combination of the R-functions method and the variation Ritz method. The numerical analysis showed that one mode of approximation cannot be used in all nonlinear dynamics of the orthotropic plates, particularly during their chaotic dynamics.
A numerical-analytical method was proposed in [26] for studying parametric vibrations of plates under the action of periodic loads applied in the middle plane. Equations of motion were obtained within the framework of classical theory. The developed approach based on applying the theory of R-functions and variational methods was used to solve the equations, making it possible to study plates of arbitrary geometric shapes with various boundary conditions. The authors used the Bolotin method to construct the region of dynamic instability. The results obtained using the developed approach were compared with the results of other authors.
In [28], the problem of dynamic stability of laminated composite plates under a combination of biaxial periodic loads and bending was considered based on the theory of plates, taking into account shear deformation. A system of second-order ordinary differential equations with periodic coefficients of the Mathieu-Hill type was obtained. The Bolotin method was applied to determine the areas of dynamic instability. The influence of modules' ratio, the number of layers, and static and dynamic parameters of the load on the dynamic instability was investigated.
The article [29] considered the dynamic stability of a viscoelastic plate under a periodic load. A differential equation was obtained that describes the dynamic stability of the plate. The solution was sought by the method of separation of variables. Linear differential equations of the third order were obtained for the time function. The prime instability domain position was affected by relaxation time and the correlation between long-term and instantaneous moduli of elasticity.
In [30], using the Bubnov-Galerkin method, the equations of motion of a shallow shell in a mixed form under the action of a longitudinal harmonic load were reduced to a discrete system with finite degrees of freedom. The chaotic behavior of systems with different degrees of freedom was analyzed. It was shown that by increasing the parameter of periodical longitudinal exciting, the investigated system starts to lose the regular behavior, forgetting its initial state, and chaotic dynamics occur.
The parametric instability of rectangular plates made of functionally graded material, subjected to biaxial periodic loading in the plane, was investigated [31]. The equations of motion were constructed, taking into account shear deformation. Hamilton's Principle was used to transform the basic equations into a linear system of Mathieu-Hill equations. With the help of the Floquet theory, areas of dynamic instability were determined. The influence of various parameters on areas of dynamic instability was investigated.
Viscoelastic rectangular plates with constant thickness, shells, and cylindrical panels were studied in a geometrically nonlinear formulation [32][33][34][35]. The models were based on the Kirchhoff-Love hypothesis. The problem was reduced to a system of nonlinear ordinary integrodifferential equations by using the Bubnov-Galerkin method. The resulting system with a weakly singular Koltunov-Rzhanitsyn kernel was solved using a numerical method based on quadrature formulas.
The article [36] was devoted to studying the stability of parametric vibrations of an isolated symmetric cross-ply laminated plate system. The plate system's mechanical model was given as an isolated symmetric cross-ply laminated constant thickness rectangular plate subjected to uniformly distributed stochastic in-plane time-dependent loading. The influence of transverse shear on the laminated plate was taken into account. The equations of motion were derived based on the shear deformation theory and the Kirchhoff-Love theory [37].
As models of viscoelastic material, the Maxwell and Kelvin-Voigt constitutive equations and the standard solid model are the simplest and most used by many researchers. These models are reduced to solving differential equations with initial conditions. The use of differential dependencies between stress and strain in dynamic problems of viscoelasticity theory leads to a specific inaccuracy. The inaccuracy is especially noticeable at the initial Appl. Sci. 2021, 11, 6029 6 of 27 moment. The reason for this inaccuracy is that the strain rate and stress rate according to these models will always be finite.
In contrast, the experiment shows an arbitrarily high rate of deformation. It should be noted that the inaccuracy made at the initial moment significantly affects the final result of the study, i.e., there is an accumulation of errors. There will be no accumulation of errors if the integral relationship between stresses and strains uses a weakly singular Koltunov-Rzhanitsyn kernel with weakly singular features of the Abel type [38][39][40][41]. The use of integral dependence significantly expands the class of problems to be solved in the theory of viscoelasticity. Simultaneously, mathematical models with integral dependence are reduced to solving complex systems of nonlinear integrodifferential equations by numerical methods [42].
The available literature analysis showed a small number of publications on nonlinear vibrations and dynamic stability of thin viscoelastic plates and shells of variable thickness. Reviewed studies show that the calculation of plates and shells of variable thickness faces significant computational difficulties. The calculation of plates and shells of this kind requires the solution of rather cumbersome equations [43][44][45]. There is a lack of suitable universal numerical methods and unified computational algorithms for solving the obtained equations.
It is necessary to carry out research in a geometrically nonlinear formulation and consider the material's viscoelastic properties to obtain a performance of the stress-strain state of thin-walled structural elements of variable thickness.
The research object is an orthotropic viscoelastic plate with an arbitrarily varying thickness subjected to a periodic compressive force in the plate's plane. The research subjects were the plate's deformations, the vibration modes of the plate, the vibration amplitudes, and the change in vibration over time depending on the thickness (geometry) of the plate, the ratio of the frequencies of the periodic force, the natural frequencies of the plate, and the viscoelastic properties of the plate material.
This research aimed to develop the theory of the orthotropic plate's dynamic stability, considering the plate's materials' viscoelastic properties.
The research objectives are as follows: -Development of a mathematical model of the dynamic behavior of plates of variable thickness; -Development of a method for numerical solutions of nondecaying systems of nonlinear integrodifferential equations with weakly singular kernels of an orthotropic viscoelastic plate of variable thickness under the influence of external periodic loads; -Clarification of the orthotropic viscoelastic plate's dynamic stability with varied mechanical and geometric parameters of the plate.

Methods
An analytical-numerical method was used to study the nonlinear vibrations of an orthotropic rectangular plate of variable thickness under the action of a uniformly distributed vibration load leading (at certain combinations of natural vibration frequencies and disturbing forces) to parametric resonance.

Analytical and Numerical Method
Consider an anisotropic viscoelastic rectangular plate of an arbitrarily varying thickness h(x, y) with sides a and b under the action of axial dynamic loads along the side a ( Figure 1). In the general case, the law of change in thickness h can be specified analytically as the function of the x and y coordinates. This function h(x, y) must have continuous second derivatives in each of the x and y coordinates. Further, a certain form of dependence (106) is used in numerical examples, but it is not the only possible one.  The mathematical model of the problem in a geometricall follows the classical Kirchhoff-Love theory. The system of equati framework of the chosen theory has the following form [46,47]: Here,  The mathematical model of the problem in a geometrically nonlinear formulation follows the classical Kirchhoff-Love theory. The system of equations of motion within the framework of the chosen theory has the following form [46,47]: Here, N x and N y are the normal forces per unit length; N xy is the shear force; M x and M y are the bending moments, H is the torque moment; p x and p y are the in-plane loads per unit area; q is the normal distributed load per unit area; P x (t) is the compressive force per unit length applied to the plate edge; ρ is the density of the material; and u, v, and w are the components of the displacement vector {U} in the directions of the axes Ox, Oy, and Oz, respectively. The compressive force P x (t) is applied to the side of dimension b and is a function of time according to the law P x = P 0 + P t cos(θ t), where P 0 is the constant part of the load, P t is the amplitude of the variable part, and θ is the frequency of the load pulsations.
The system of Equations (1)-(3) is supplemented by the corresponding boundary conditions [46], which will be used to solve the problems: 1.
All edges are simply supported: at All edges are clamped: at Two opposite edges are simply supported, the other two edges are clamped: at The initial conditions at t = 0 [46] are as follows: . v 0 (x, y), and . w 0 (x, y) are given functions.

Computational Methods
The solution of the system (1)-(3) is obtained by using the Bubnov-Galerkin method based on a polynomial approximation of the deflection and displacements; discretization is performed in spatial variables at each moment and is reduced to solving systems of nonlinear integrodifferential equations with weakly singular kernels with variable coefficients. The improved numerical method [48] based on quadrature formulas was applied to solve the obtained nondecaying system with a weakly singular Koltunov-Rzhanitsyn kernel [40,41,49].
A relaxation kernel Γ(t) = Ae −βt t α−1 , (0 < α < 1) with three rheological parameters was first proposed in [40]. Publications [40,41] explain a method for determining the parameters A, α, and β. In [41], these parameters are tabulated. The parameters of specific materials are determined by combining the experimental curve and the theoretical curve plotted in logarithmic coordinates.
The number of relaxation kernels for isotropic and orthotropic plates is different: (i) The mechanical parameters of the viscoelastic isotropic are the same in all directions. Therefore, the viscoelastic properties of the plate material are described with one core with three rheological parameters. After applying the Bubnov-Galerkin method, most of the dynamic problems of viscoelastic thin-walled structures are reduced to solving nondecaying systems of integrodifferential equations of the following form: w 0nm , n = 1, 2, . . . , N; m = 1, 2, . . . , M, where u nm = u nm (t), v nm = v nm (t), and w nm = w nm (t) are unknown functions of time; X nm , Y nm , Z nm , φ 1nm , φ 1nm , ψ 1nm are continuous functions in the field of changing arguments; and a klnm , b klnm , c klnm , λ 2 klnm , χ 2 klnm , and ω 2 klnm are given constant numbers. Integration of system (4) twice on t leads to its integral form. Then, t = t i , t i = i∆t, i = 1, 2, . . . (∆t = const is the interpolation step) is assumed. Replacing integrals with quadrature formulas for calculations u nm = u nm (t i ), v nm = v nm (t i ), and w nm = w nm (t i ) results in the following system: The next stage of the numerical method is the regularization of nonlinear integrodifferential Equation (5) with weakly singular kernels. By changing variables the integral at the Koltunov-Rzhanitsyn kernel with a singularity of the following form: takes the form After the change of variables, the integrand concerning z becomes regular. Direct replacement of the integrals in the system (5) by a sum over some quadrature formula yields a numerical solution to the numerical solution of system (5).
The trapezium formula used is where the coefficients are Thus, due to the double integration of the original system (4) over the time t and the use of the quadrature formula, system (5) is obtained to find deflections w nm = w nm (t i ) and displacements u nm = u nm (t i ),v nm = v nm (t i ). The solutions of the system of Equation (5) are found by the Gauss method.

Analytical Solution for an Orthotropic Viscoelastic Rectangular Plate of Variable Thickness
This section presents the obtention of the analytical solution for the dynamic stability of an orthotropic viscoelastic rectangular plate of variable thickness. The solution is based on the integral relationship between stresses and strains using the weakly singular Koltunov-Rzhanitsyn kernel with weakly singular features of the Abel type.
The components of the vector of forces {N} = N x , N y , N xy and moments {M} = M x , M y , M xy for a plate of symmetric structure in matrix form can be written as follows: where ε y = ∂v ∂y The stiffness matrixes [C] and [D] are given as follows [47]: where the coefficients of the stiffness matrixes C ij and D ij (ij = 11, 22, 12, 16, 26, 66), depending on the mechanical characteristics of the viscoelastic material, and the mass per unit area coefficient m are determined as follows: where B ij is the coefficient of the stiffness [47] and Γ * , Γ * ij are integral operators with relaxation kernels Γ(t) and Γ ij (t), respectively: In operator form, the system of equations of motion (1)-(3) is as follows: Then, If the considered plate has orthotropic properties, then for the coefficients C 16 = C 26 = 0 and D 16 = D 26 = 0. Then, following the relations (11) and (12), the stiffness matrices take the following form: In this case of plate's orthotropic properties, the coefficients C ij and D ij are expressed in terms of the elastic constants E 1 , E 2 , G 12 , µ 1 , and µ 2 as follows: Here, E 1 and E 2 are elastic modulus in the direction of the axes x and y, G 12 is shear modulus, and µ 1 and µ 2 are Poisson's ratios.
If the plate has isotropic properties (E 1 = E 2 , µ 1 = µ 2 ), then the elements of the stiffness matrix take on a more straightforward form with two elastic constants E (shear modulus) and µ (Poisson's ratio): The Koltunov-Rzhanitsyn singular kernel [49] is used in the calculations as relaxation kernel Γ(t), Γ ij (t), i, j = 1, 2: If solving problems of the dynamics of viscoelastic plates based on the Kirchhoff-Love hypothesis in an isotropic formulation, then only one relaxation kernel Γ(t) with three different rheological viscosity parameters A, α, and β is involved in the system of integrodifferential equations. Then, in the orthotropic formulation, already 5 different kernels Γ(t), Γ ij (t), i, j = 1, 2 are involved with 15 rheological parameters A ij ,α ij and β ij .
The system of Equations (28)- (30) is the most general. The motion equations of isotropic and orthotropic viscoelastic rectangular plates of variable thickness can be obtained from this system. The corresponding initial and boundary conditions are assumed to be given.
For example, consider a rectangular orthotropic viscoelastic plate of variable thickness h(x,y) with sides a and b that is dynamically loaded in-plane along the side a by a periodic load. The equation of motion, based on the Kirchhoff-Love hypothesis, in a geometrically nonlinear setting has the form (28)-(30), taking into account (38)- (47).
If the following dimensionless variables are introduced into Equations (28)- (30): and the previous notation is retained, then a dimensionless system of nonlinear integrodifferential equations for the problems of parametric vibrations of an orthotropic viscoelastic rectangular plate of variable thickness is obtained. In this case, the operators L ij from the expressions (31)-(37) take the following form: The resulting system of Equations (56)-(63) with the corresponding boundary and initial conditions describes the motion of an orthotropic viscoelastic rectangular plate of variable thickness.
By changing variables t − τ = z 1 α , 0 ≤ z ≤ t α , (0 < α < 1), the integrals over the Koltunov-Rzhanitsyn kernels with a singularity of the following form Note that after the change of variables, the integrands of the system (108)-(110) with respect to z become regular.
Then, assuming t = t i , t i = i∆t, i = 1, 2, . . . (where ∆t is the integration step) and replacing the integrals with quadrature trapezoidal formulas to calculate the unknowns w inm = w inm (t i ), u inm = u inm (t i ), and v inm = v inm (t i ), a system of recurrent formulas is obtained.

Numerical Results
An efficient computational algorithm was developed, and the program was compiled in the algorithmic language Delphi and implemented on a computer. The results of calculations for various physical and geometric parameters of an orthotropic viscoelastic rectangular plate are displayed in graphs.
In the general case, the law of change in thickness can be given analytically. In the numerical examples, the plate thickness changes according to the following law: i.e., the plate thickness decreases linearly. Here, α * is the thickness variability parameter; h 0 is the plate's thickness at α * = 0.

Behavior of the Plate versus Viscoelastic and Inhomogeneous Properties of the Material
The influence of the viscoelastic properties of the material on the behavior of the plate was studied. Figure 2 shows the time history of the deflection w for the following parameters: λ = 1, δ = 25, α* = 0.3, w 0 = 0.01, δ 0 = 0.3, δ 1 = 0.5, θ = 1.1, ∆ = 1.2, µ 1 = 0.32, µ 2 = 0.27, and g = 0.3833. The influence of the viscoelastic properties of the material on the behavior of the plate was studied. Figure 2 shows the time history of the deflection w for the following parameters: λ = 1, δ = 25, α* = 0.3, w0 = 0.01, δ0 = 0.3, δ1 = 0.5, θ = 1.1, Δ = 1.2, μ1 = 0.32, μ2 = 0.27, and g = 0.3833. The results show that taking into account the viscosity of the plate's material leads to a decrease in the vibration amplitude obtained as a result of the calculation. At the initial moment, the amplitudes of the oscillations coincide. Then, over time, the results obtained in the elastic and viscoelastic formulations differ from each other by 40% or more. Figure 3 shows the vibration modes of an elastic and viscoelastic rectangular plate. The results show that taking into account the viscosity of the plate's material leads to a decrease in the vibration amplitude obtained as a result of the calculation. At the initial moment, the amplitudes of the oscillations coincide. Then, over time, the results obtained in the elastic and viscoelastic formulations differ from each other by 40% or more. Figure 3 shows the vibration modes of an elastic and viscoelastic rectangular plate. The results show that taking into account the viscosity of the plate's material leads to a decrease in the vibration amplitude obtained as a result of the calculation. At the initial moment, the amplitudes of the oscillations coincide. Then, over time, the results obtained in the elastic and viscoelastic formulations differ from each other by 40% or more. Figure 3 shows the vibration modes of an elastic and viscoelastic rectangular plate.      Figure 4 shows the effect of inhomogeneous material properties of a viscoelastic plate on deflection at λ = 1, δ = 25, α* = 0.3, w0 = 0.01, δ0 = 0.3, δ1 = 0.5, θ = 1.1, μ1 = 0.32, μ2 = 0.27, and g = 0.3833. Here, curve 1 is in the setting of the viscoelastic properties only along the shear directions, and curve 2 is in the setting of the viscoelastic properties in all directions. The results show that taking into account the plate material's inhomogeneous properties leads to decreased vibration amplitude.     Figure 6 shows the change in the deflection of a viscoelastic plate as a function of time for different values of the thickness variability parameter α * at λ = 1, δ = 25, A ij = 0.05, w 0 = 0.01, δ 0 = 0.3, δ 1 = 0.5, θ = 1.1, ∆ = 1.2, µ 1 = 0.32, µ 2 = 0.27, and g = 0.3833. It is seen that an increase in this parameter, i.e., a decrease in the plate's stiffness, leads to a decrease in the vibration amplitude.  (2) all sides are clamped.

Behavior of the Plate at Various Thicknesses of the Plate
The results obtained show that, for given initial conditions and damping, with an increase in the plate's number of clamped sides, the vibration amplitude decreases and the vibration frequency increases. In general, the results reflect the physical behavior of the considered plate vibrations. They are in qualitative agreement with the results obtained by the authors earlier [7][8][9][10]45,48].
It is of interest to compare the above results with the results obtained by the finite element method. It was noted above that such calculations are associated with significant computational and methodological difficulties. Of the publications found [51][52][53][54][55], only one publication provides a satisfactory possibility for the comparison. In [55], a finite element model based on the four-variable refined plate theory was used for forced vibration analysis of rectangular viscoelastic laminated composite plates integrated with a piezoelectric layer.  The results obtained show that, for given initial conditions and damping, with an increase in the plate's number of clamped sides, the vibration amplitude decreases and the vibration frequency increases. In general, the results reflect the physical behavior of the considered plate vibrations. They are in qualitative agreement with the results ob-