Computer simulation of oscillatory processes of viscoelastic elements of thin-walled structures in a gas flow

Results of numerical investigation of dynamic behavior of deformed wing aircraft in a gas flow are presented in the paper. Vibrations with respect to deflections are described by a system of integro-differential equations in partial derivatives. Using the Bubnov-Galerkin method, the problem is reduced to a system of ordinary integro-differential equations, where time is an independent variable. The solutions of integro-differential equations are determined by a numerical method based on the use of quadrature formulas. Computational algorithms and a package of applied programs have been created to solve problems on nonlinear flutter of viscoelastic elements of an aircraft. The reliability of the solution of the problem is confirmed by comparison with known numerical and analytical results. The effect of different boundary conditions on critical flutter velocity is studied. Critical velocity and critical flutter time of viscoelastic plates are determined. It is shown that the singularity parameter  affects not only the vibrations of viscoelastic systems, but also critical time and critical flutter velocity. It is stated that consideration of viscoelastic properties of plate material leads to 40 60% decrease in critical flutter velocity.


Introduction
Enormous scale of the development in aviation industry and shipbuilding necessitates the further development of the theory and practice of mathematical simulation. The study of structure material with viscoelastic and nonlinear properties, the consideration of which has great theoretical and practical importance, approximates the theory of calculation to the actual conditions of structure operation. Therefore, the problems of the theory of hereditary elasticity attract serious attention of researchers.
Of particular interest are the nonlinear problems of the theory of hereditary elasticity, which, apart from their practical importance, are of considerable scientific interest in the spheres of mechanics, mathematical physics and computational mathematics. This is due to the fact that nonlinear problems of the theory of hereditary elasticity are reduced to boundary and initial-boundary value problems for nonlinear weakly singular integral-differential equations withpartial derivatives, the complete investigation of which is connected with the development of new methods for solving weakly singular integraldifferential equations.
The above-mentioned scientific problem gives grounds to assert that the development of adequate mathematical models, numerical methods and algorithms for solving nonlinear integral-differential equations of dynamic problems of the hereditary theory of viscoelasticity is actual.
In connection with this, the development of mathematical models of individual elements of aircraft made of composite material is becoming very important.
One of the main difficulties for a complete understanding of the supersonic flutter phenomenon is the fact that critical velocity of the flutter depends on a large number of parameters. At present, the difficulty in isolating many of these factors in experimental study does not allow us to obtain a satisfactory agreement between experimental and theoretical results.In literature there are numerous reviews of the problem under investigation. An extensive bibliography is given by Marco Amabili [1,2], Farbod Alijani and Marco Amabili [3], Y. Wang and Z. M. Wang [4]. The development of problems on the flutter-plate strip, plates and panels with account of viscoelastic properties of structure material is reflected in publications by Mojtaba [34,35] and others. In [5][6][7][8][9], the Kelvin-Voigt model is used to describe the strain processes occurringin viscoelastic materials.
As is known, exact solutions to the problem of flutter, even in the simplest cases, are non-existent to this day. Therefore, there are different opinionsregarding the effect of viscoelastic properties of structure material on critical velocity of flutter [12][13][14][15]. In [14,15], theoretically (by average method), and in [12,13] by computational experiment it is shown that the effect of viscoelastic properties of structure material on critical flutter velocity in both steady and transient processes leads to a significant decrease in value of critical velocity of a flutter.In [16,17] solving the problem of the flutter of viscoelastic strip (in the case of the exponential kernel of heredity), the conclusion is drawn that critical flutter velocity, in both ideally elastic and viscoelastic cases, does not differ much, and the "viscous" properties of material affect the character of strip motion in subcritical regiononly. Similar conclusions have been established earlier in [12,18,19], and this phenomenon is explained by the fact that the exponential kernels do not correctly describe the hereditarily deformable properties of structure material.
This drawback lies in the fact that the relaxation kernel, proportional to the strain rate, has a finite value at initial time, however the experiment shows an arbitrarily high strain rate, i.e. at t 0, R(t), that contradicts the experiment [20][21][22], and as a result, when solving any dynamic problem (not only the flutter problem), errors accumulate over time and the results will be distorted in comparison with reality processes.In fact, in [23] it is exactly stated that the vibrations of any viscoelastic system (beams, plates and shells) under constant load occur along the creep curve and attenuate over time along this curve. This law is fulfilled if any weakly singular kernel of heredity of Abel-type is used as the kernel of heredity; the use of an exponential kernel over time due to errors accumulation leads to a severe violation of this law of motion of viscoelastic systems.
In [12,18,19] it is shown that if to solve viscoelastic flutter problems with weakly singularkernels of heredity, then there immediately appears a significant effect of viscosity parameter and singularity on the values of critical velocity and critical time of both linear and nonlinear flutter problems. Therefore, the development of a new method for solving and demonstrating the reliability of results of the study of dynamic stability problem, both in ideally elastic and in hereditary deformable systems in a gas flow, is a very urgent problem. The present paper is devoted to the above problem. The accuracy and convergence of the method is tested on known results [24][25][26] related to the flutter of elastic plates and viscoelastic strips [16,17].

Formulation of the problem 2.1 Nonlinear flutter of viscoelastic plates
Consider a rectangular viscoelastic plate with sides a and b, which is flowed over from one side by a supersonic gas flow of velocity V. Aerodynamic pressure is taken into account by the piston theory of A.A. Ilyushin [27].
For the case of finite deflections of a plate commensurate with its thickness h, strains of viscoelastic plate are described by equations: In accordance with the boundary conditions: G 2 ) -hinged support on two edges and fixing on the other two: (1) is taken in the form: Where w nm =w nm (t) andФ nm =Ф nm (t)are the sought for time functions;  nm (x,y) -known functions, depending on boundary conditions: ).
Integrating equation (1) by the Bubnov-Galerkin method with relationships (3) -(5), the following system of integro-differential equations (IDE) areobtainedwith respect to w nm (t) andФ nm (t): ,  (7) are written in general form and are called the basic resolving IDE oftwo-dimensional problems of dynamics of viscoelastic systems. Many problems of vibrations and stability of viscoelastic plates arereduced to equations of the form (7).
Comparison of different cases of plate fixation shows that with an increase in a number of fixed sides of the plate,flutter critical velocity increases.
It can be seen from the obtained results that if the exponential kernel ( = 1) is used, flutter velocity decreases by approximately 5%, and when the Koltunov-Rzhanitsyn kernel is used this velocity decreases by 45% relative to critical velocity of the flutter of ideally elastic plates. Therefore, when using exponential kernels, flutter velocity of viscoelastic plate practically coincides with critical flutter velocity for ideally elastic plates. These conclusions and results fully agreewith the conclusions and results given in [16,17], where critical flutter velocities are determined by a numerical-analytical method.  According to the results obtained by many researchers [18,19,28,29] using the integral stress-strain law with weakly singular kernelof heredity, it is obvious that the viscosity parameter leads to a decrease in critical velocity and an increase in critical time (Figs. 1, 2). With vanishingly small internal friction, the velocity of the panel flutter is approximately 2 times less than the one calculated under the assumption that internal friction is completely absent [12,13]. The results obtained by the authors fully agreewith the conclusions and results given in [13,16,17]. Table 2 shows the effect of the shape of plate deflection on critical velocity of the flutter. AtN=5, critical velocity is 1540m/s, and for N=6 this velocity is 1595m/s. The difference is 3.5%.
6 Table 2. Effect of the shape of plate deflection on critical velocity of the flutter   4 shows the graphs of vibration amplitude as a function of time for studying the role of geometric and aerodynamic nonlinearity. In a linear statement, the amplitude of vibrations increases rapidly (curve 1), flow velocity is greater than critical one, and with account of nonlinearity, the amplitude of vibrations decreases (curve 2).

The problem of vibrations and stability of viscoelastic strip
Consider the problem of vibrations and stability of elastic and viscoelastic strip in a gas flow [16,17] and present a comparative analysis of the results of solution with the ones obtained by the proposed method.
In a rectangular coordinate system, the strip occupies region 0yl, x0. On one side it is flown overby a gas flow with velocity vector V=Vn 0 , n 0 = {cos, sin}.
Strip vibrations are described by equation [16,17]   Assuming that given in [16,17], and retaining the previous notations, equation (8) is written in the form: Here the notations are introduced,and other notations correspond to the ones assumed in [16,17].
Solution of equation (9) is chosen in the form: (10) Substituting (10) into (9) at=0 we get 0 ) 1 ( Equation (11) with initial conditions (12) is solved by a numerical method. This method is based on the use of various analytical transformations that allow initial systems to be reduced to the systems of integral equations with regular kernels and stable numerical integration ensuring the solution of problems with a high degree of accuracy.
According to numerical method, we regularize the IDE system (11) with singular kernels. Substituting the variables the integral with the Koltunov-Rzhanitsyn kernel with a singularity of the following form has the form Note that after substitution of variables, the integrand with respect to z becomes regular. Then, assuming thatt=t i , t i =it, i=1,2,…(t=const-is the interpolation step) and replacing the integrals by some quadrature formulas (in particular, of trapezoids), we get where the coefficients are On the basis of this method, an algorithm for numerical solution of system (11) is described. Integrating system (11) twice with respect to t, it can be written in integral form; with rational transformation we eliminate the singular features of integral operator Thus, according to numerical method with respect to the unknowns, a system of algebraic equations is obtained. The Gauss method is used to solve the system. Based on the developed algorithm, a package of applied computer programs is created. Results of calculations are given in Tables 3.  Table 3 shows the results of specific calculations for the parameters values [16,17]: The third column shows the results [16,17] when critical flutter velocities are determined by a numerical-analytical method.
As seen from theresults obtained, for ideally elastic and viscoelastic strips (in the case of exponential kernel of heredity) the critical velocities of the flutter exactly coincide with the results given in [16,17]. For viscoelastic strip with a weakly singular heredity kernel, this velocity decreases.

Results reliability of solving the dynamic stability problem of hereditary deformable systems
Conclusions on results reliability of solving the dynamic stabilityproblem of hereditary deformable systems in a gas flow, obtained by eliminating weakly singular features of integral and integro-differential equations.
Thus, the computational experiments carried out according to the algorithm of the proposed method for solving the flutter problem of viscoelastic systems completely refute some intuitive conclusions and natural dissatisfaction with the effect of viscoelastic properties on critical velocity and critical flutter time.Therefore, when solving dynamic problems of mechanics of a deformable rigid body, it is necessary to use the integral stressstrain law with weakly singular heredity kernels of Abel-type. Numerical experiment of the dynamics of corresponding structures has shown the presence of a significant effect of this feature on the nature of their vibrations, for example, an account of weakly singular feature of heredity kernel leads to a significant decrease in the value of critical velocity and an increase in critical time (Figs. 5, 6).This fundamentally new mechanical effect may be of interest to the specialists in the field of designing such structures [33].

Conclusions
It should be noted that the algorithm of the proposed method makes it possible to investigate in detail the effect of rheological parameters on the character of vibrational stability of hereditary deformed systems, in particular, in the study of flutter problem of ideally elastic systems.
As seen from Tables 1 and 3 the reliability of study results is proved by testing with known results related to the flutter of elastic plates [24,26] and viscoelastic strip [16,17]. In both cases, a satisfactory agreement of the solutions obtained is shown;that shows the reliability and high accuracy of the proposed calculation procedure.