Selected Problems of Random Free Vibrations of Rectangular Thin Plates with Viscoelastic Dampers

The main motivation of this work was to present a semi-analytical extension of the correspondence principle in stochastic dynamics. It is demonstrated for the stochastic structural free vibrations of Kirchhoff–Love elastic, isotropic and rectangular plates supported by viscoelastic generalized Maxwell dampers. The ambient temperature of the plate affects the dampers only and is included in a mathematical model using the frequency–temperature correspondence principle. The free vibration problem of the plate–viscoelastic damper system is solved using the continuation method and also the Finite Element Method (FEM). The stochastic approach begins with an initial deterministic sensitivity analysis to detect the most influential parameters and numerical FEM recovery of the polynomial representation for lower eigenfrequencies versus these parameters. A final symbolic integration leads to the first four basic probabilistic characteristics, all delivered as functions of the input uncertainties.


Introduction
Probabilistic mechanics is a topic that has been extensively studied, e.g., in [1,2], and one of its fundamental numerical methods-the Stochastic Finite Element Method (SFEM)was invented and has been applied in the context of thin rectangular plate bending [3]. Considering multiple geometric scale uncertainties, the bending analysis of thin plates may also be performed using the Wavelet-based Stochastic Finite Element Method [4]. The SFEM was efficiently applied for the stochastic dynamic response analysis of graphite-epoxy composite plates [5]; other studies in this area can be found in [6,7]. An interesting scientific and engineering problem, the modeling of propeller blades, is presented in [8]. The authors did not use the typical random approach; however, they did perform an analysis of the deviation histogram of machining errors that may be random in nature.
The probabilistic structural response for the free damped vibrations of thin elastic and isotropic plates resting on viscoelastic supports is considered in this work. The nonlinear eigenproblem is solved here to determine the natural frequencies of the plate-viscoelastic damper system, and its solution is obtained thanks to the iterative continuation method presented by Lewandowski et al. [9][10][11]. The authors considered different types of viscoelastic dampers based on the generalized Maxwell model of a damper. An experimental study considering a generalized Maxwell model for nonlinear viscoelastic dampers was comprehensively performed by Lu et al. [12]. A comprehensive overview of some other deterministic and stochastic methods used in dynamics can be found in [13]. Several uncorrelated Gaussian random design variables are considered in this study, with an initial sensitivity study leading to the selection of the most influential parameters. The Least-Squares Method enables the determination of the random polynomials of eigenfrequencies, whose further integration with the Gaussian kernel finally returns the probabilistic characteristics. It should be underlined that the FEM experiments were entirely programmed in the Octave environment, whereas sensitivity and probabilistic analyses were all implemented in the computer algebra system MAPLE. The most important novelty of this work is the common application of a probabilistic numerical apparatus for a solution of an eigenproblem obtained using the continuation method.

Eigenvibration Analysis Methodology
The main purpose of this analysis is to determine the first few natural frequencies of thin rectangular elastic and isotropic plates supported by viscoelastic dampers and also their first four probabilistic characteristics. The numerical analysis of this problem is based upon the Finite Element Method (FEM) with a regular discretization including 4-node-plate finite elements with linear approximation functions, as shown in Figures 1 and 2 [14]. The deformation vector w e i of the i-th node within the finite element e can be written as

Eigenvibration Analysis Methodology
The main purpose of this analysis is to determine the first few natural frequencies of thin rectangular elastic and isotropic plates supported by viscoelastic dampers and also their first four probabilistic characteristics. The numerical analysis of this problem is based upon the Finite Element Method (FEM) with a regular discretization including 4-nodeplate finite elements with linear approximation functions, as shown in Figures 1 and 2 [14]. The deformation vector of the -th node within the finite element can be written as The displacement field within the element is expressed as a linear combination of shape functions ( , ): where = [ 1 , 2 , 3 , 4 ] and = [ 1 , 2 , 3 , … 12 ]. The element stiffness matrix and the consistent mass matrix are defined in the traditional way. The stiffness matrix is derived here analytically, and the mass matrix is derived numerically using 16point Gaussian quadrature. Further, it is known that the equation of motion of a structure with viscoelastic dampers can be written in the following form [9,10]: where denotes the global plate-damping matrix.  The application of the Laplace transform with zero initial conditions leads to the following transform of Equation (3): where is the ℒ-transform of , and ̅ can be expressed as The displacement field within the element e is expressed as a linear combination of shape functions N e k (x, y): where w e = w e 1 , w e 2 , w e 3 , w e 4 T and N e = N e 1 , N e 2 , N e 3 , . . . N e 12 . The element stiffness matrix Ke and the consistent mass matrix M e are defined in the traditional way. The stiffness matrix is derived here analytically, and the mass matrix is derived numerically using 16-point Gaussian quadrature. Further, it is known that the equation of motion of a structure with viscoelastic dampers can be written in the following form [9,10]: where C denotes the global plate-damping matrix. The application of the Laplace transform with zero initial conditions leads to the following transform of Equation (3): where q(s) is the L-transform of q(t), and f (s) can be expressed as In this formula, n d is the total number of dampers attached to the plate at selected nodes of a finite element mesh, and L r is a global matrix indicating the location of the r-th damper in the plate. A viscoelastic damper is represented graphically in Figure 3, and it consists of m spring-dashpot elements and an additional spring element. Each of the Maxwell elements contains a viscous part with the constant c j and an elastic part with the constant k j , where j = 1, 2, . . . , m. The application of the Laplace transform with zero initial conditions leads to the following transform of Equation (3): where ̅( ) is the ℒ-transform of ( ), and ̅ ( ) can be expressed as In this formula, is the total number of dampers attached to the plate at selected nodes of a finite element mesh, and is a global matrix indicating the location of theth damper in the plate. A viscoelastic damper is represented graphically in Figure 3, and it consists of spring-dashpot elements and an additional spring element. Each of the Maxwell elements contains a viscous part with the constant and an elastic part with the constant , where = 1, 2, … , . All quantities appearing in Equation (5) can be expressed as follows: = ⁄ is the quotient of the stiffness and damping coefficients of the -th Maxwell element. Obviously, the stiffness and damping parameters, and , of the individual elements constituting the viscoelastic damper attached to the structure additionally depend upon the temperature. Let the damper parameters be known for a certain reference temperature 0 and have the values and . The damper model constants at temperature can be expressed using the frequency-temperature correspondence principle as [9,10] All quantities appearing in Equation (5) can be expressed as follows: where ν j = k j /c j is the quotient of the stiffness and damping coefficients of the j-th Maxwell element. Obviously, the stiffness and damping parameters, k j and c j , of the individual elements constituting the viscoelastic damper attached to the structure additionally depend upon the temperature. Let the damper parameters be known for a certain reference temperature T 0 and have the values k j and c j . The damper model constants at temperature T can be expressed using the frequency-temperature correspondence principle as [9,10] k j = k j ; j = 0, 1, 2, . . . , m, c j = c j α T ; j = 1, 2, . . . , m.
The shift factor α T is a function of temperature T and can be expressed by the William-Landel-Ferry formula as [15,16] where C 1 and C 2 are experimental constants. After substituting Equation (5) into Equation (4), the L-transform of Equation (3) of the motion of a plate with viscoelastic dampers takes the following form: where Equation (10) represents a nonlinear eigenproblem that is solved for the eigenvalue s and the eigenvector q(s) using the continuation method comprehensively described by Lewandowski in, e.g., [9,10]. Below, the foundations of this method are quoted.
In the case of Equation (10), the components containing the variable s in the first power are multiplied by the parameter κ ∈ [0; 1]. The equation can then be re-written as where In order for the elements of the eigenvector q corresponding to the eigenvalue s to be determined unambiguously, an additional normalizing equation of the following form is introduced into the matrix in Equation (13): where a has a given value.
In the first step of the continuation method, in Equation (13), the parameter κ 1 = 0 is assumed, and the generalized eigenproblem is solved.
This problem was solved in the Octave program using the built-in command 'eig', which allows both standard and generalized eigenproblems to be solved. As a result of solving this problem, the first approximations of eigenvalues s 2 , . . . , q (1) 3n are obtained. On their basis, the parameter a where j = 1, 2, . . . , 3n, is determined. In the l-th step (l = 2, 3, 4, . . .), the increment ∆κ l is assumed, and the Newton method is used to solve the system of Equation (13) with the additional Equation (15). For this purpose, the system of incremental equations of the Newton method is solved using . This system of equations takes the following form: where The derivatives in Equations (18b)-(18d) are calculated as follows: The increments δ q and δs are obtained from the system of Equation (17), and the following are calculated: Successive approximations of the j-th eigenvalue and the j-th eigenvector in the l-th step of the algorithm are calculated until the desired accuracies ε 1 and ε 2 , of the final results are achieved, that is, until the following inequalities are satisfied: The final values of s j and a (k) j obtained in the l-th step are taken as starting values for step l + 1 and the new parameter κ l+1 = κ l + ∆κ l+1 .
The procedure described above is carried out up to the value of the parameter κ = 1, when the final eigenvalues and eigenvectors for the nonlinear eigenproblem in Equation (10) are obtained. The obtained eigenvalues of the problem in Equation (10) are complex numbers of the form s j = µ j + iη j . On this basis, the j-th natural frequency ω j of the structure and the non-dimensional damping ratio γ j of the j-th mode of vibration are determined from the formulas: The continuation method makes it possible to calculate the first few natural frequencies and the corresponding non-dimensional damping ratios without having to solve the entire nonlinear eigenproblem. Using proprietary algorithms written in the Octave programming language, the natural frequencies of a plate equipped with viscoelastic vibration dampers are determined. In the program, the user can independently select the positions of the selected number of dampers in the FEM mesh nodes. With the help of proprietary software, the matrices occurring in the nonlinear eigenproblem in Equation (10) are determined. Using the built-in command, the program solves the generalized eigenproblem in Equation (16). Then, the program uses the above-described iterative algorithm of the continuation method, which allows the unknown natural frequencies of a plate equipped with viscoelastic vibration dampers to be determined.

Sensitivity and Uncertainty Analyses
The so-called normalized sensitivity gradients are determined using the following standard definition [1]: where v j denotes the mean value of the given parameter v j . The following design parameters affecting the natural frequencies of the rectangular plate are checked at the initial stage: (i) geometric dimensions of the plate l x × l y × H, (ii) material constants of the plate E, ν p and ρ p , (iii) damper parameters k 0 , k 1 and c 1 , (iv) the ambient temperature of the plate, namely, T, and (v) the reference temperature T 0 . Further probabilistic analysis is carried out for two parameters exhibiting the highest positive and negative sensitivity coefficients. The eigenfrequencies ω i of the plate under consideration are all found via the polynomial basis via the Least-Squares Method fittings made on the basis of several FEM solutions for varying values of the parameter v [1,2]. Statistical optimization of this basis order is employed through the common minimization of the fitting variance and the maximization of the correlation factor. Finally, the basic probabilistic characteristics, i.e., expected values, standard deviations, coefficients of variation, skewness and kurtosis, are computed. The following integral definitions are applied:

Numerical Experiment
A square isotropic plate fixed on one edge was discretized using the 14 × 14 plate rectangular finite element mesh, whose material properties are E = 205 GPa, ν p = 0.3 and ρ p = 7850 kg/m 3 . The plate dimensions are equal to l x × l y × H = (2.0 × 2.0 × 0.01) m.
Three viscoelastic dampers are attached in the middle and at both ends of the free edge of the plate (see Figure 4).  The influence of temperature on the values of the above-mentioned parameters is taken into account by applying the frequency-temperature correspondence principle. In order to calculate the value of the shift function from Equation (9), the values of the constants 1 = 19.5 and 2 = 80.2 were adopted. The initial sensitivity analysis results were computed analytically using polynomial responses and are compared in Table 1 below. Quite expectedly, the two most influential parameters for the given plate are its thickness (the minimum value of the gradient) and the edge length (the maximum gradient). So, these two parameters were further selected for stochastic analysis. They were treated as Gaussian variables having expected values equal to the mean values given above and a coefficient of variation belonging to the interval ( ) ∈ [0.00, 0.20]. Such a wide interval was assumed to check theoretical variations in all characteristics being computed, and it includes the statistical scattering of all possible measurement techniques. It should be noted that the time effort and computer power required for the semi-analytical probabilistic solutions were only a little bit greater than the deterministic origin. Polynomial approximation was used to describe the response function. The degree of the  The influence of temperature on the values of the above-mentioned parameters is taken into account by applying the frequency-temperature correspondence principle. In order to calculate the value of the shift function from Equation (9), the values of the constants 1 = 19.5 and 2 = 80.2 were adopted. The initial sensitivity analysis results were computed analytically using polynomial responses and are compared in Table 1 below. Quite expectedly, the two most influential parameters for the given plate are its thickness (the minimum value of the gradient) and the edge length (the maximum gradient). So, these two parameters were further selected for stochastic analysis. They were treated as Gaussian variables having expected values equal to the mean values given above and a coefficient of variation belonging to the interval ( ) ∈ [0.00, 0.20]. Such a wide interval was assumed to check theoretical variations in all characteristics being computed, and it includes the statistical scattering of all possible measurement techniques. It should be noted that the time effort and computer power required for the semi-analytical probabilistic solutions were only a little bit greater than the deterministic origin. Polynomial approximation was used to describe the response function. The degree of the The influence of temperature on the values of the above-mentioned parameters is taken into account by applying the frequency-temperature correspondence principle. In order to calculate the value of the shift function from Equation (9), the values of the constants C 1 = 19.5 and C 2 = 80.2 were adopted. The initial sensitivity analysis results were computed analytically using polynomial responses and are compared in Table 1 below. Quite expectedly, the two most influential parameters for the given plate are its thickness (the minimum value of the gradient) and the edge length (the maximum gradient). So, these two parameters were further selected for stochastic analysis. They were treated as Gaussian variables having expected values equal to the mean values given It should be noted that the time effort and computer power required for the semi-analytical probabilistic solutions were only a little bit greater than the deterministic origin. Polynomial approximation was used to describe the response function. The degree of the polynomial was assumed to be 4 or 5 depending on the necessary accuracy of matching the response function. Examples of polynomials obtained for a random plate side length and a random plate thickness are given in Chapter 5 by the relations in (28) and (29), respectively. Figures 6 and 7 show the graphs of the dependence of the expected value, variance, skewness and kurtosis on the coefficient of variance when the random variable is the plate side length and its thickness, respectively; each of the graphs in Figures 6a-d and 7a-d shows these characteristics for the first five natural frequencies, ω 1 -ω 5 . These probabilistic coefficients were all found from their integral definitions, and they can be treated as exact in the probabilistic context (no convergence studies are necessary). polynomial was assumed to be 4 or 5 depending on the necessary accuracy of matching the response function. Examples of polynomials obtained for a random plate side length and a random plate thickness are given in Chapter 5 by the relations in (28) and (29), respectively. Figures 6 and 7 show the graphs of the dependence of the expected value, variance, skewness and kurtosis on the coefficient of variance when the random variable is the plate side length and its thickness, respectively; each of the graphs in Figures 6(a-d) and 7(a-d) shows these characteristics for the first five natural frequencies, 1 -5 . These probabilistic coefficients were all found from their integral definitions, and they can be treated as exact in the probabilistic context (no convergence studies are necessary).  It is seen in Figure 6 that Gaussian uncertainty in the plate length causes some nonlinear increases even for the expected values of the fundamental free vibrations. This is in contradiction to the case illustrated in Figure 7, where they are simply constant. Moreover, the plate length randomness greatly amplifies the uncertainty in this problem, because its output-to-input ratio equals almost 3; the plate thickness shows a direct interrelation between the input and output CoVs (Figure 7b). Interestingly, the largest statistical dispersion is always associated with the first eigenvalue. Finally, it is seen that both the skewness and kurtosis monotonously increase together with an additional increase in the input CoV in Figure 6. Hence, the positive non-symmetry and concentration about the expected values remarkably increase together with the input uncertainty level. One notices that the differences between probabilistic characteristics for various eigenfrequencies are rather small. Figure 7c,d show that higher-order probabilistic coefficients relevant to the plate thickness oscillate about (for the 1st and the 2nd) or are almost equal to 0 (for higher eigenfrequencies). It is seen in Figure 6 that Gaussian uncertainty in the plate length causes some nonlinear increases even for the expected values of the fundamental free vibrations. This is in contradiction to the case illustrated in Figure 7, where they are simply constant. Moreover, the plate length randomness greatly amplifies the uncertainty in this problem, because its output-to-input ratio equals almost 3; the plate thickness shows a direct interrelation between the input and output CoVs (Figure 7b). Interestingly, the largest statistical dispersion is always associated with the first eigenvalue. Finally, it is seen that both the skewness and kurtosis monotonously increase together with an additional increase in the input CoV in Figure 6. Hence, the positive non-symmetry and concentration about the expected values remarkably increase together with the input uncertainty level. One notices that the differences between probabilistic characteristics for various eigenfrequencies are rather small. Figure 7c,d show that higher-order probabilistic coefficients relevant to the plate thickness oscillate about (for the 1st and the 2nd) or are almost equal to 0 (for higher eigenfrequencies).

Comparative Analysis-Results Validation
In order to better illustrate the previously obtained results, calculations were performed using three probabilistic approaches: the Semi-Analytical Method (SAM), the Stochastic Perturbation Technique (SPT) and Monte Carlo simulations (MCSs).
These analyses were applied for a plate exhibiting Gaussian uncertainty in the first natural frequency, uniquely defined by its mean value and the specific range of its coefficient of variation, i.e., α(b) ∈ [0.00, 0.025]. The global response function for random plate side length l x = l y = l was obtained in the form of the following fourth-order polynomial: The choice of the degree of the approximating polynomial to the random quantity was dictated by the sufficient accuracy of matching the response function. The number of trials for the Monte Carlo simulations was equal to 150,000. The expected values E(l), coefficients of variation α(l), skewness β(l) and kurtosis κ(l) of the first natural frequency are presented in turn in Figure 8. Considering the large variations in the resulting statistics, the expected values and coefficients of variation for the first natural frequency are presented in Tables 2 and 3. was dictated by the sufficient accuracy of matching the response function. The number of trials for the Monte Carlo simulations was equal to 150,000. The expected values ( ), coefficients of variation ( ), skewness ( ) and kurtosis ( ) of the first natural frequency are presented in turn in Figure 8. Considering the large variations in the resulting statistics, the expected values and coefficients of variation for the first natural frequency are presented in Tables 2 and 3.   A very good convergence of the results of the SAM and SPT approaches can be observed. On the other hand, MCS results are in good agreement only for the coefficient of variations α(l).
Another random parameter for which the validation of the calculation results was performed is the plate thickness. In this case, the global response function for a random plate thickness h was obtained in the form of the following fifth-order polynomial: Similar to the above, the expected values E(h), coefficients of variation α(h), skewness β(h) and kurtosis κ(h) of the first natural frequency are presented in turn in Figure 9. Additionally, the expected values and coefficients of variation for the first natural frequency are presented in Tables 4 and 5. plate thickness h was obtained in the form of the following fifth-order polynomial: Similar to the above, the expected values (ℎ), coefficients of variation (ℎ), skewness (ℎ) and kurtosis (ℎ) of the first natural frequency are presented in turn in Figure  9. Additionally, the expected values and coefficients of variation for the first natural frequency are presented in Tables 4 and 5. (a) (b) (c) (d) Figure 9. Expected values (a), coefficients of variation (b), skewness (c) and kurtosis (d) for a random plate thickness. Similar to the previous case, a very good convergence of the results of the SAM and SPT approaches can be observed. The MCS results are in quite good agreement only for kurtosis and α(h) values from 0.0 to about 0.1.
The Monte Carlo simulation technique is the most time-consuming in terms of numerical calculations when performed on a typical PC computer using the MAPLE v.21 computational package. The computation time ratio can be expressed simply by the quotient t SAM t MCS or t SPT t MCS , which ranges from 1/50 to 1/100, depending mainly on the number of MCS trials.

Conclusions
(1) The theoretical and computational studies presented in this work clearly show that the common application of the semi-analytical stochastic approach with continuation methods allows for the fast and accurate determination of the probabilistic coefficients of free vibrations. It is demonstrated that the output randomness in rectangular elastic thin isotropic plate eigenfrequencies is usually not larger than the input statistical scattering of their design parameters. The only exception is in the plate dimension statistics; however, successful measuring techniques are so accurate now that despite the huge sensitivity to this parameter, the realistic coefficient of variation is less than a few percent. The fact that the largest resulting statistical dispersion is obtained for the first eigenfrequency may be important in the reliability assessment of various dynamical systems. This is due to the fact that a limit function, whose probability serves as the basis of the reliability index, is introduced as the difference between the induced vibrations and the lowest eigenfrequency.
(2) A comparative analysis for the first eigenfrequency was performed using three probabilistic approaches: the Semi-Analytical Method (SAM), the Stochastic Perturbation Technique (SPT) and Monte Carlo simulations (MCSs). For the first two methods, results with very high accuracy were obtained. The Monte Carlo simulation showed convergence only for selected random moments-the coefficient of variation for the random plate side length and kurtosis for the random plate thickness.
(3) A continuation of this research can include stochastic extensions of the continuation method in the numerical analysis of forced vibrations, possibly with the use of non-Gaussian random design parameters, too. In the case of any mathematical difficulties with computer algebra integration, the iterative generalized stochastic perturbation technique implemented as the SFEM is recommended. Further uncertainty analysis may be alternatively completed by the application of probabilistic entropy or its relative version presented recently in the literature [17,18].