Mathematical Modeling Numerical simulation of the resonance effect at Re-entry of a rigid body with low inertial and aerodynamic asymmetries into the atmosphere Lyubimov V.V. Samara State Aerospace University Abstract. We consider a motion relative to centre of mass of a rigid body with low inertial and aerodynamic asymmetries at re-entry into the atmosphere. Numerical and mathematical simulation provide a means to study the resonance effect that results in changes of the direction of the rigid body rotation. It is shown that the realization of the effect in question is observed when selecting the initial conditions of the integration and defined values of the asymmetry parameters. Keywords: resonance effect, numerical simulation, asymmetry, rigid body, atmosphere Citation: Lyubimov V.V. Numerical simulation of the resonance effect at Re- entry of a rigid body with low inertial and aerodynamic asymmetries into the atmosphere. Proceedings of Information Technology and Nanotechnology (ITNT-2015), CEUR Workshop Proceedings, 2015; 1490: 198-210. DOI: 10.18287/1613-0073-2015-1490-198-210 1. Introductions It is known [1]-[2] that the presence of the small asymmetry can lead to realization of the resonance that is observed in re-entry of a spacecraft viewed as a rigid body (RB) into the atmosphere. A continuous resonance leads to violation of the technological limitations on the RB’s angle of attack or angular velocity. Besides the resonance itself, some secondary resonance effects can lead to violation of the technological limitations. The mentioned effects were found by Sadov Y.A. in dynamical systems with slow and fast variables [3]. The essence of these effects is: in the area close to the resonance, characteristic evolutions of slow variables of the system are observed, caused by this resonance. From the mathematical point of view, the secondary resonance effects are explained with presence of the resonant frequency mistuning in the denominators of averaging method’s highest approximations that are received in the non-resonant case. Secondary resonance effects were studied particularly in the spherical rotation of the heavy asymmetrical rigid body [4]. Applying to the task of re-entry of an asymmetrical RB in the atmosphere, the mentioned resonance effects were studied in detail in the work [5]. In particular, it has been revealed that the secondary resonance effect itself can lead to the realization of 198 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Lyubimov V.V. Numerical simulation of the… the RB’s strong spin-up that was received earlier in the work [1]. In the following, using of traditional mathematical analysis methods allowed to extend the classification of the secondary resonance effects applying to the task of atmospheric re-entry of RB [6]. Other features of the task of atmospheric re-entry that affect the behaviour of the descent craft, are random character of the initial conditions on the separatrix (that lead or don’t lead to the transfer into the resonance area) [7], and also the possibility of losing part of data on the craft’s movement [8]. In the work [9], it is shown that the negative consequences of the resonance effects in respect to the task of atmospheric re-entry of the spacecraft with low mass and aerodynamic asymmetries can be solved with introducing the control over the amount of asymmetry. In the process of numerical simulation of movement of the RB with low mass and aerodynamic asymmetries, a secondary resonance effect has been considered, leading to changes in the direction of the craft’s rotation [10]. Numerical simulation and analytical study of the similar resonance effect in the task of atmospheric re-entry of RB with low aerodynamic and inertial asymmetries are of practical interest. 2. Mathematical models Initial non-linear equation of the asymmetrical RB’s movement relating to the centre of mass in the atmosphere appear, for example, in the work [11]. However, the significant non-linearity of these equations and instability of frequencies of the system significantly complicate the application of such studies for detection and detailed study of secondary resonance effects. It is known that application of the integral manifold method [12] allow to decrease the order of differential equation system. For instance, in the work [13] it is shown that the application of one of the variants of integral manifold method makes it possible to receive from initial non-linear system of asymmetrical RB’s movement relating to the centre of mass an equation system that is approximated to the non-linear. With allowance only for aerodynamic and inertial asymmetries, this system takes the following form: m sin 2  x   2 cos 2(  3 ) , (1) I x 1,2 2m a 1,2 sin   1,2 2 sin 2        x   cos(2  23 )  Fa  2a  2m Aa 4mzп 2a SL dq  cos(  1 )  3 , (2) Fa IFa dt   x  1,2 cos  . (3) Here, ε is a small parameter that characterizes the value of RB’s small inertial and aerodynamic asymmetries, and the slowness of changing the dynamic pressure q ( dq / dt  O( ) ), x is the RB’s angular velocity relating to axis with least moment 3 of inertia; α is the space angle of attack;   mzп qSLctg / I ; S and L are 199 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Lyubimov V.V. Numerical simulation of the… characteristic area and size of the RB; I x and I= I y  I z are the moments of RB’s inertia relating to axes of XYZ body-fixed coordinate system; I x  I x / I ,   п   / 2; φ is the aerodynamic roll angle; m A and m are generalized parameters of aerodynamic in inertial asymmetries, m A  m   m  , 1 A 2 A 2 2 Ix 1,2  x   , 2 m1A   1  I x  x  31,2 * 2 (mф  C z)tg, sin   m A / m A , y x 1 1 2a mzп m2A   1  I x  x  31,2 * 2 (mф  C y)tg, cos   m A / m A , m A  m A / 2 , z x 1 2 2a mzп m  ( I yz )2  (I )2 , sin 23  I / m , cos 23   I yz / m , mфy , mzф - aerodynamic shape asymmetry coefficients, I  ( I y  I z ) / 2, I  ( I z  I y ) / 2, I  I / I , I yz  I yz / I - inertial (dynamic) asymmetry, mzп - coefficient of 2 restoring moment, a  I x 2x / 4  2 , x  1,2 - resonant frequency mistuning, Fa (x , , q) is the known function of the slow variables. Equation system (1)-(3) is solved together with three differential equations that describe change in characteristics of the RB’s centre of mass movement [1]: airspeed, flight-path inclination angle and flight altitude. Mentioned three variables are considered as slow. Equation system (1)-(3) takes into account the possibility of implementing the main resonance: x  1,2  0 . When passing through the main resonance, significant perturbations of movement parameters (compared to multiple resonances of higher orders) are observed [14]. By solving the equation x  1,2  0 , we find the resonance values of the angular velocity:  rx   . (4) I 1 x I The sign in the expression (4) matches the sign of angular velocity x . The equation system (1)-(3) contains in its right side the dependency on the phase of rapid movement θ, which complicates the analysis of non-resonant evolutions at secondary resonance effects. For further simplification of this system, we will use averaging method in non-resonant case [15]. System (1)-(3) applies to the class of systems with single rapid phase  and several slow variables. It has a standard form for application of averaging method: u  U (u, , ), (5) 200 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Lyubimov V.V. Numerical simulation of the…   (u). (6) Here u  x ,  is the vector of slow variables, (u)  x  1,2 , U (u, , ) is the vector function of right side of equations (1) and (2). Equations for centre of mass movement parameters do not depend on  , so they are not taken into account in process of averaging. System (1)-(3) is averaged on the non-resonant areas of movement under the rapid phase  . Resulted averaging of equations for slow variables have the following form: du o  A1 (u o )  2 A2 (u o )  3 A3 (u o )  ... , (7) dt where the functions Ai , i=1,2,3,... are determined through standard averaging method [15]. After averaging the system (1)-(3) in non-resonant case we obtain: A1  0 , A2  0 A3  0 . Hence, the evolution of slow variables x ,  , caused by the secondary resonance effects, is determined by the members of the third approximation of averaging method. The equations for slow variables x ,  , being averaged on non- resonant areas of movement and taking into account first three approaches of averaging method, are d x  A  m g 2 g3  A  A m g3  A g  3  ( m g )  ( m g3 2 )  dt   3  3   2        A   m cos(21  22 ) , 2 3(m g3 )2  g 2  (   g   )  (8)       2 4 8  Averaged equations (8) and (9) contain generalized parameters of asymmetry m A , m , 1 , 2 in numerator. The parameters m A , m define the value of RB’s  asymmetry, and parameters 1 , 2 characterize the location of asymmetry on the rigid body. From the equations (8)-(9) it appears the main resonance  =0 contributes to the realization of secondary resonance effects, for all the members of mentioned averaged equations contain frequency mistuning x  1,2 in their denominators. In particular, while approaching the resonance in the boundaries of non-resonant area, a decrease of x  1,2 value occurs, which leads to an increase of speed in averaged values of variables x ,  . Such behaviour of slow variables is typical for secondary resonance effects. 201 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Lyubimov V.V. Numerical simulation of the…  A d 3  m g3 g1   A  A 2     ( m g )  m g  3    3 3 dt   2   A 2 m g3  A g g  A   (m g3) 1  12  (m g3 )     3      A m g2 g3   A A    ( m g 3 )   2m g 3  3 x     A A (m )2 g3 g1 g3 (m ) g3   2 2  g     7 g2  4 2    2   2 4 x     (m )2 g32  2  2 g1      A 2   m cos(21  22 )   g4 . 3   2   g1   (9) 2 4   2      8  2a 1,2 sin  1,2 2 sin 2  1,2 2 sin 2  2a 2 Here g1  ( x  ) , g2  , g3  , Fa 2a Ix Fa 4mzп 2 SL dq g4  . IFa dt To analyse the secondary resonance effects, we also have to obtain the equation for the second derivative of the averaged angular velocity: d 2 x  x d x  x d  x d     . (10) dt 2 x dt  dt  dt 3. Numeric simulation with different initial conditions The numerical integration of the equation system (1)-(3) (with allowance of the equation for derivative dq/dt) is implemented via four-staged Runge-Kutta method with adaptive stepsize. The secondary resonance effect is considered, when the angular velocity x (t ) , changing from initial positive value, envelopes the x (t ) r ascending curve and reaches significant values with further keeping its sign in the area of descending curve x (t ) . At certain decrease of initial conditions of r integration or decrease of asymmetry parameters, the situation of angular velocity behaviour changes dramatically. If in the area of ascending curve x (t ) the angular r velocity x (t ) also increases, in the area of descending curve x (t ) , the realization r of continuous resonance occurs that can be followed with transition to the negative range of values x (t ) . Let us assume that the time of integration is 300 s. We consider the influence of various initial condition of integration on the realization of the aforementioned resonance effects. Here, two cases shall be pointed. 202 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Lyubimov V.V. Numerical simulation of the… Case 1. The resonance effect followed by change of rotation direction is observed in various values of initial angular velocity of the body. Let us assume that the condition x (0)  x (0) has been met. In the first case, generalized parameters of r asymmetry take the following values: m A  0.011 , m  0.04 , 1  2  0 . Initial values of the angle of attack: (0)  5 grad, and initial value of aerodynamic roll angle is * (0)  0 . In the process of numerical integration of equations (1)-(3), a narrow interval of 1 positive values for initial angular velocity x (0)  [17, 23]s , has been found. With boundary values of x (0) within this interval, the resonance effect in question is observed. Let us have a more detailed review of these numerical results. 1 1. We assume that the initial value of angular velocity x (0)  23s . In the Fig.1, the upper curve describes the increase of angular velocity x (t ) , and the curve with maximum is associated with the change of resonance values of angular velocity rx (t ) . At that, in the area of ascending branch of resonance curve rx (t ) , a monotonic non-resonance increase of x (t ) occurs. Further, the withdrawal from resonance values is observed, and angular velocity x (t ) is stabilized. From the Fig.1 it appears that the craft acquires the terminal positive angular velocity which exceeds 1 the initial velocity x (0) more than twice. Thus, with x (0)  23s , a slow increase of spacecraft angular velocity positive values is observed. The angle of attack also increases slowly from the initial value of 5 grad, but does not reach the value in 90 angle degree, typical for realization of main resonance. The direction of RB’s rotation remains unchanged and positive. The increase of angular velocity x (t ) in the area of ascending branch of rx (t ) curve is explained by positivity of the derivate (8). The change of signature of the x (t ) curve’s convexity happens due to flip of the second derivate’s sign (10). The increase of the x (t ) fluctuation area after passing the maximum of the resonance curve can be explained with the following: in the right side of the equation (8) there are terms directly proportional to sinα. After 220 seconds of flight, the angle of attack takes values from 40 to 50 grades, which leads to a certain perturbation on the angle of attack in the right side of the equation (8). 1 Fig. 1. – Change of angular velocity at x (0)  23s 203 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Lyubimov V.V. Numerical simulation of the… 1 2. Let the initial angular velocity be x (0)  20s . In this case value of the initial 1 angular velocity is within the x (0)  17 s range. Also, from the Fig.2 it appears that the angular velocity x (t ) in the area of ascending curve x (t ) changes r similarly to the previously considered variant. But in the area of the descending resonance curve, angular velocity x (t ) passes a bit lower than on the Fig.1. Also, after passing the maximum of the resonance curve, the fluctuation amplitude x (t ) increases significantly, which also promotes the angular velocity to reach its resonance values x (t ) . In its turn, an increase of the resonance curve fluctuation r area is observed. The expansion of the angular velocity fluctuation area occurs due to increase of angle of attack that reaches 75 grads on the 200 th second. It should be noticed that the continuous resonance does not occur nevertheless. 1 3. Let the initial angular velocity be x (0)  17 s . From the Fig.3 it appears that in the variant under consideration (similarly to the first one), a growth of non- resonance RB spin-up in the area of ascending resonance curve is realized. 1 Fig. 2. – Change of angular velocity at x (0)  20s Yet, the spin-up is slightly smaller than in variant 2, so the increase of amplitude x (t ) in the area of the curve maximum leads to meeting the resonance values x (t ) r . Further, the main resonance is realized, and the angle of attack quickly reaches 90 grades. After the resonance realization, angular velocity returns to the non-resonance movement in the negative range of its values. Comparing the results of the Fig.1 and Fig.3, it may be noticed that with values of 1 the initial angular velocity that are the boundaries of the x (0)  [17, 23]s interval, in the process of evolution x (t ) , opposite directions of the RB rotating motion are reached. Case 2. Resonance effect followed by change of the rotation direction is also observed in various values of the initial angle of attack (0) . In the second case, generalized parameters of asymmetry take the following values: m A  0.011 , m  0.002 , 1  2  0 . Initial value of the angular velocity: x (0)  25s 1 , and 204 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Lyubimov V.V. Numerical simulation of the… initial value for the aerodynamic roll angle is (0)  0 . In the process of numerical integration of the equations (1)-(3), we found a narrow range of values for initial angle of attack (0)  [3; 5.7] grad, with boundary values of which, the studied resonance effect is observed. The results obtained in the second case are shown on the Fig. 4-6. 4. Let the initial value for the angle of attack be (0)  3 grad. On the Fig.4, in the area of ascending branch of the resonance curve, a monotonic non-resonant increase x (t ) occurs, similar to variant 1 shown on Fig.1. Further, a withdrawal of x (t ) from resonance values x (t ) is also observed, and the angular velocity of the body is r 1 stabilized at x (0)  47s . Yet, different from variant 1, the angular speed changes in monotonic way throughout the whole integration interval. It occurs due to fact that the angle of attack at the slow increase from initial value of 5 grad reaches the value of 38 grad in the end of integration interval. As a result, no perturbations of angular velocity is observed taking place with angles of attack higher than 50 grades and more. The direction of the craft’s rotation also remains positive throughout the whole integration interval. The increase of the angle velocity x (t ) in the area of ascending branch of the curve x (t ) is also explained by positivity of the derivative (8), and the r change of the character of the x (t ) curve’s convexity occurs due to flip of the second derivative’s sign (10). 1 Fig. 3. – Change of angular velocity at x (0)  17s Fig. 4. – Change of angular velocity in the second case at α(0)=3 grad 5. Let us assume that the initial value for the angle of attack is α(0)=4.7 grad. Here, the value for the initial angle of attack is within the section of (0)  [3;5.7] grad. At that, the angular velocity x (t ) in the area of ascending curve x (t ) changes r 205 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Lyubimov V.V. Numerical simulation of the… similarly to variant 4. Yet in the area of descending resonance curve, angular velocity x (t ) passes slightly lower than on the Fig.4. Compared to variant 2, in this variant the angular speed perturbation in the area of descending resonance curve is not significant. This is due to fact that the angle of attack in the considered area takes value of about 40 grades. Comparing to variant 4, a slightly increase takes place in the fluctuation are of the resonance curve. Also, in the area of maximum x (t ) on the r curve x (t ) , another inflexion point arrives. As a result, on the curve x (t ) there are two inflexion points, where the mentioned curve changes the character of the convexity twice. From the Fig.5 it appears that in the end of the integration interval, angular velocity x (t ) doesn’t reach its resonance values x (t ) . r Fig. 5. – Change of angular velocity in the second case at α(0)=4.7 grad 6. Let the initial angle of attack be α(0)=5.7 grad. From the Fig.6 it appears that in the considered variant, similarly to variants 4 and 5, a growths of the non-resonant RB spin-up occurs in the area of the ascending resonance curve. Still, the mentioned spin- up is slightly lower than in the variant 5. As a result, the expansion of range of the resonance values leads angular velocity x (t ) to reaching the resonance values rx (t ) Further, main resonance is realized, and angle of attack reaches 90 grades. After leaving the continuous resonance, RB’s angular velocity x (t ) takes non- resonance values in the negative range. Comparing the results presented on Fig.4 and Fig.6, a conclusion can be made that with the initial values for angle of attack that are the boundaries of section (0)  [3; 5.7] grad, in the process of the evolution the angular velocity x (t ) reaches opposite values of the rotating motion. 4. Numerical simulation at different values of asymmetry parameters We perform a numerical simulation that allow us to study the issue of the influence of different values of asymmetry parameters on the realization of the considered resonance effect. We conduct an individual analysis of influence of inertial and aerodynamic parameter values on the considered resonance effect. Hence, two cases are found. Case 3. Realization of the resonance effect with changing the rotation direction occurs at different values of the generalized parameter of m inertial asymmetry. In the process of the RB rotating motion numerical simulation, a range of small values 206 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Lyubimov V.V. Numerical simulation of the…  was found for m  [0.012, 0..040] parameter, for which three types of body’s angular velocity evolution, similar to the first case, were revealed. Let us give them more detailed overview. At values of the parameters m  0.040 , m A  0.011 , 1 1  2  0 and initial condition of integration x (0)  23s , (0)  5 grad, (0)  0 we obtain the result shown in Fig.1. When integrating the system (1)-(3) with asymmetry parameter m  0.012 (all other parameters and initial conditions are the same as in variant 1), we obtain the result shown in Fig.7. This result confirms that with parameter value m  0.012 angular velocity changes the rotation direction. Fig. 6. – Change of angular velocity in the second case at α(0)=5.7 grad In numerical simulation with parameter m  0.013 we obtain result that conforms the result on Fig.2 in the qualitative manner. Case 4. Also, the realization of the resonance effect with change of the rotation direction occurs with different values of m A aerodynamic asymmetry generalized parameter. In numerical simulation of RB rotating motion, a narrow range of small values for parameter m  [0.008, 0.011] was found. With three values from this A range we obtain three types of spacecraft’s angular speed evolution which are similar to the first case. Let us give them more detailed overview. The initial variant conforms the case shown on the Fig.1. Indeed, the parameters of the craft and initial Fig. 7. – Change of angular velocity at m  0.012 data are identic to the variant 1: m∆ =0.040, mᴬ=0.011, θ1 =θ2 =0, ωx (0)=23c -1 , (0)  5 grad, φ(0)=0. If the aerodynamic asymmetry parameter is considered equal to m A  0.008 (leaving other parameters and initial integration conditions unchanged), we obtain the result of integration presented on the Fig.8. This result is approximate to the results shown in Fig.3 and Fig.7. Here, the angular velocity x (t ) 207 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Lyubimov V.V. Numerical simulation of the… also transits to non-resonant negative range of its values after the realization of the continuous resonance and withdrawal from it. Numerical simulation of angular speed x (t ) at the parameter value m A  0.0082 allow us to obtain a result identical to the result of simulation shown on Fig.2. Fig. 8. – Change of angular velocity at m A  0.008 5. Conclusion In this work, application of numerical simulation together with analytical research allows us to study in detail the new resonance effect observed in the process of movement in relation to centre of mass of the rigid body with low inertial and aerodynamic asymmetries at re-entry into the atmosphere. An unusual feature of the studied resonance effect is the following: while changing in the boundaries of the ascending branch of the resonance curve, angular velocity x (t ) in absolutely most cases “glides” with delay of relatively resonant curve x (t ) without reaching the r resonant values. At the transition to the area of the descending part of the curve x (t ) r , two typical cases were peculiar to angular velocity x (t ) behaviour:  angular velocity kept increasing the positive values (practically by the linear law) with further gradual non-linear transition to some large constant value;  angular velocity reached resonant value, resonance realization occurred with further withdrawal from it and transition to non-resonant negative range of values. In numerical simulation, narrow intervals of initial integration conditions x (0) , (0) and a small interval of asymmetry of m , m A parameters values were found. Taking into account the values on the boundaries of intervals in the process of numerical integration led us to realization of the two described typical cases. While accounting intermediate values from aforementioned intervals in the process of numerical simulations, we found a number of cases where an increase of amplitudes of x (t ) and x (t ) oscillations was observed after passing the maximum by the r resonance curve. At that, change of the x (t ) curve convexity character occurred. The conducted numerical simulation and results of analytical research allow us to draw the following conclusions: 208 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Lyubimov V.V. Numerical simulation of the… 1. The values of initial integration conditions x (0) , (0) and values of asymmetry parameters m , m A , 1  2  0 together have an effect on realization of the considered resonance event. 2. One may talk of a common pattern of influence by the values x (0) and m , m A on the realization of the found resonance effect that takes place in case when mentioned initial conditions and asymmetry parameters are chosen from the intervals found in this study. 3. The comparison of the numerical results obtained in process of integration of the approximated system (1)-(3) with the results of numerical integration of initial non-linear system has shown good qualitative agreement with these results. 4. The further study of the considered dynamical event of x (t ) sign flip is of practical interest, for instance in the process of movement of a RB with low inertial and aerodynamic asymmetries in the atmosphere. References 1. Yaroshevsky VA. The movement of the body in an uncontrolled atmosphere. Moscow: Mechanical engineering, 1978; 168. [in Russian] 2. Platus DH. Roll Resonance Control of Angle of Attack for Re-entry Vehicle Drag Modulation. Journal of Guidance, Control, and Dynamics, 1981; 4(5): 632-636. 3. Sadov YuA. Secondary Resonance Effects in Mechanical Systems. Mechanics of Solids, 1990; 4: 20-24. 4. Zabolotnov YuM, Lyubimov VV. Secondary Resonance Effects in the Rotation of a Rigid Body about a Fixed Point. Mechanics of Solids, 2002; 1: 49-59. 5. Zabolotnov YuM, Lyubimov VV. Secondary Resonance Effect in the Motion of a Spacecraft in the Atmosphere. Cosmic Research, 1998; 36(2): 194-201. 6. Lyubimov VV. Asymptotic Analysis of the Secondary Resonance Effects in the Rotation of a Spacecraft with Small Asymmetry in the Atmosphere. Russian Aeronautics (Iz VUZ), 2014; 57(3): 245-252. 7. Zabolotnov YuM. Statistical Analysis of Movement of Light Capsule around of the Centre of Mass at Re-Entry into the Atmosphere. Cosmic Research, 2013; 51: 1-12. 8. Sedelnikov AV. The Usage of Fractal Quality for Microacceleration Data Recovery and for Measuring Equipment Efficiency Check. Microgravity Science and Technology, 2014; 26(5): 327-334. 9. Lyubimov VV. Dynamics and Control of Angular Acceleration of a Re-Entry Spacecraft with a Small Asymmetry in the Atmosphere in the Presence of the Secondary Resonance Effect. International Sibefian Conference on Control and Communications (SIBCON), 2015; 1-4. 10. Lyubimov VV, Lashin VS. Numerical simulation of resonance effect with the change of the direction of rotation at atmospheric descent asymmetric nanosatellite. Samara: SamNC RAN, 2015; 250-253. [in Russian] 11. Belokonov VM, Belokonov IV, Zabolotnov YuM. Fast calculation of the trajectories in the atmosphere reduce uncontrolled spacecraft with regard to their motion relative to the center of mass. Space activities, 1983; 21(4): 512-521. [in Russian] 12. Strygin VV, Sobolev VA. Separation of Motions by the Integral Manifolds Method. Moscow: Nauka, 1988. [in Russian] 209 Information Technology and Nanotechnology (ITNT-2015) Mathematical Modeling Lyubimov V.V. Numerical simulation of the… 13. Zabolotnov YuM. The method of investigation of resonant vibrational motion of a nonlinear system. Izvestiya RAN: Mechanics of Solids, 1999; 1: 33-45. [in Russian] 14. Zabolotnov YuM. The asymptotic analysis of quasi-linear equations of motion in the atmosphere the spacecraft with a small asymmetry and. Space Exploration, 1994; 32(2): 22-23. [in Russian] 15. Moiseev NN. Asymptotic Methods in Nonlinear Mechanics. Moscow: Nauka, 1986. [in Russia] 210 Information Technology and Nanotechnology (ITNT-2015)