Chaotic Behavior of the Biharmonic Dynamics System

Motion of a biharmonic system under action of small periodic force and small damped force is studied. The biharmonic oscillator is a physical system acting under a biharmonic force like a sin θ b sin 2θ. The article contains biharmonic oscillator analysis, phase space research, and analytic solutions for separatrixes. The biharmonic oscillator performs chaotic motion near separatrixes under small perturbations. Melnikov method gives analytical criterion for heteroclinic chaos in terms of system parameters. A transition from chaotic to regular motion of the biharmonic oscillator was found as the heteroclinic chaos can be removed by increasing the coefficient of a damping force. The analytical results obtained using Melnikov method have been confirmed by a good match with numeric research.


Reduction of the Biharmonic Dynamics System to the Duffing Equation
As is well known, nonlinear systems can perform a chaotic motion under the action of periodic forces 1-4 .Frequently the Duffing equation is used to illustrate chaos 1-4 , and the chaotic behavior of various forms of the Duffing equation 5 , some of which exhibit twofrequency excitation 6 as well as the chaotic motion of Duffing system with bounded noise 7 have been investigated.However the Duffing equation in the expanded generalized form has many mechanical applications and it can be interesting to researchers.
In this paper, we will study periodically driven biharmonic dynamic system with a damping force: θ a sin θ b sin 2θ ε a 1 sin θ b 1 sin 2θ c cos ωt − δ θ, 1.1 where ε and δ are assumed to be small positive parameters; ω > 0 is the frequency of the International Journal of Mathematics and Mathematical Sciences external force; and a, a 1 , b, b 1 , and c are coefficients.The terms in ε and in δ in 1.1 can be considered as small perturbations.
If ε 0 and δ 0, then the periodic and the damping force are absent, and we have the conservative system describing the motion of the undisturbed biharmonic oscillator as θ a sin θ b sin 2θ. 1.2 We will show that 1.1 is equivalent to the Duffing equation for small values of θ.It is known that and 1.1 can be written as If the variable θ is a small of order ε, then we have the Duffing equation where λ − a 2b , μ a 8b /3!.

Mechanical Applications
The disturbed system 1.1 has applications in space flight mechanics when studying a problem of a spacecraft motion about its center of mass in the atmosphere.Atmospheric reentry is a critical phase for space vehicles.Dynamic stability issues play a crucial role for the success of their mission.For effective breaking while descending in the rarefied atmosphere of Mars or Titan we can compensate blunt-shaped spacecrafts by a small increase in length 8-11 .We will consider blunt-shaped spacecraft.An aerodynamic restoring moment strongly influences the motion of the spacecraft relative to the center of mass.The aerodynamic restoring moment coefficient for an axisymmetric rigid body can be written as 12, 13 where θ is the spatial angle of attack θ is defined as the angle between symmetry axis and the velocity of the spacecraft V , Figure 1 , m O is the pitching moment coefficient concerning the leading edge body, c n is the normal force coefficient, x c x c /L is the nondimensional coordinate of the center of mass, and L is aerodynamic reference length shield diameter .Such spacecrafts can have three positions of equilibrium according to the angle of attack θ: stable position at the points θ * 0 and θ * π; and unstable in the third intermediate point θ * ∈ 0, π 12, 13 .The restoring moment can be approximated as the biharmonic dependence Figure 1 m θ a sin θ b sin 2θ. 1.7 The presence of the second harmonic in 1.7 causes the possibility of appearance of an additional equilibrium position-saddle point on a phase portrait.For the considered spacecrafts position, θ 0 is stable; therefore, a derivative of the function m θ θ with respect to the angle of attack θ at this point is negative: and if there exists an intermediate position of equilibrium inside the interval of 0, π , then It is obvious that 1.9 and 1.10 are valid simultaneously when b < 0. Note that the dependence of m θ θ given in Figure 1 satisfies conditions 1.9 and 1.10 .The stable position occurs not only in the point of θ 0 but also in the point of θ π when 1.9 for the spacecraft is fulfilled.The motion of the spacecraft in a neighborhood of International Journal of Mathematics and Mathematical Sciences θ π cannot be allowed, because in this case the back part of the spacecraft will move toward to an approach flow.
The aerodynamic restoring moment can be written by where S is reference surface and q is the dynamic pressure.Disturbance will simulate periodic change of position of the center of mass where x c0 is an initial position of the center of mass.Using expression 1.6 , we represent the aerodynamic restoring moment as Thus, the planar motion of the spacecraft about the center of mass can be described by where I is the transverse moment of inertia, M d −δ 1 θ is small damping moment.Obviously, the equation of spacecraft motion 1.15 corresponds to the biharmonic dynamics system 1.1 at c 0 up to constants.Let us observe, that the undisturbed equation 1.2 describes the motion of a known mechanical system-a heavy material point on a circle, rotating about a vertical axis 14 where g is the gravitational acceleration, l is the radius of the circle and Ω is the angular velocity of the circle.

Aim and Structure of the Article
Let conditions 1.9 and 1.11 be satisfied, and let there be three positions of equilibrium.
In this case, we find three regions in the phase portrait separated by separatrixes.Under the effect of disturbances the phase trajectory θ θ θ can repeatedly intersect the separatrixes, thus moving from one an area to another, is accompanied by a jump change of the variable θ.We can observe chaos.It is an accepted fact in the theory of nonlinear dynamic systems that knowledge of the stable and unstable manifolds of hyperbolic equilibrium or hyperbolic periodic orbits may play a crucial role in understanding many issues of dynamics 1, 2 .For many dynamic systems, the only general way of studying such stable and unstable manifolds is computing them numerically.However, in some cases we can obtain analytic solutions.
The aim of this paper is to analyze the motion of the disturbed system 1.1 near the undisturbed separatrixes and to define the boundaries of chaos.We will carry out the theoretical studies by means of the Melnikov method 15 .The Melnikov method is an analytical tool used to determine to first order, the existence of homo/heteroclinic intersections and so chaotic behavior.The Melnikov method allows us to obtain a necessary condition for the existence of chaos, therefore, numerical simulation is needed to confirm the predicted behavior and to give a deeper understanding of the global dynamics of the system.We will show that our theoretical results are in good agreement with the results of numerical calculations.
The present paper is structured in the following way.Section 2 gives the analysis of the unperturbed motion of the biharmonic dynamics system 1.1 and the phase portrait.On the phase portrait characteristic regions of the possible motions are found and for these areas the analytical solutions of the equation of the unperturbed motion are obtained.The main features of the phase space of the unperturbed system are defined.In Section 3 the Melnikov criterion for the perturbed system is analytically calculated for various areas of the phase portrait with help of the theory of residues.In Section 4 the Melnikov criterion is numerically calculated for various areas of a phase portrait for the case of disturbed motion of the spacecraft.By means of computer numerical simulations of the disturbed motion, we use several numerical techniques to check the validity of the analytical criterion for chaos obtained using Melnikov method.

The Unperturbed Solutions
The ε and δ terms in 1.1 are considered as small perturbations.If ε 0 and δ 0, then periodic and dissipation forces are absent, and we have the conservative system 1.2 .It is obvious that if b 0 or a 0 replacement of variables ϕ 2θ we have the equation of a mathematical pendulum.However, if 1.9 and 1.11 are satisfied b < 0 , then biharmonic dynamic system has a more complicated phase portrait in comparison with the mathematical pendulum or with the heavy material point on a circle 1.16 .The equilibrium positions of the system 1.2 are defined from 1.10 .If 1.9 and 1.11 are satisfied, then the undisturbed system 1.2 has four equilibrium positions at θ ∈ −π, π : two stable-center type θ 0, π 2.1 and two unstable-saddle type where b < 0. The center θ * −π coincides with the center θ * π.At θ * → −π and at θ * → π the speeds θ coincide, therefore we can say, that phase trajectories are closed on a cylindrical phase space.We will consider the evolution of the cylindrical space in the range θ ∈ −π, π .We will separate two regions A 0 and A 1 , divided by the two saddles s 1 and s −1 Figure 2 .It is necessary to note, that the region A 1 of the development of the cylinder undergoes a break at θ π, −π.From 2.2 it follows, that if the coefficient a is equal 0, the saddle s 1 is in the position: θ * π/2.At positive values of the coefficient a > 0 the saddle s 1 belongs to the interval: θ * ∈ 0, π/2 , and at negative values a < 0 the saddle s 1 belongs to the interval: The following energy integral corresponds to 1.2 :

International Journal of Mathematics and Mathematical Sciences
where E is total energy.The biharmonic oscillator as well as the mathematical pendulum can perform oscillations and rotation.The shape of the phase portrait depends on the potential energy: W θ a cos θ b cos 2 θ.

2.4
The centers 2.1 correspond to the minimum of the potential energy 2.4 , and the saddles 2.2 -to the maximum of the potential energy 2.4 .If E > W * , where W * W θ * , then the motion is possible in the outer regions Figure 2 .In the opposite case E < W * the motion can occur in any of the inner regions, depending on initial conditions.The equality E W * corresponds to the motion along separatrixes.In this case, the two saddles s 1 and s −1 are connected by four heteroclinic trajectories.First of all, we will consider the separatrixes, limiting the region A 0 .Separating the variables in the energy integral 2.3 , and taking into account 2.2 and 2.4 , the equation of the motion on the separatrixes can be written as in the integrated form where Substituting the variables: we can rewrite 2.5 as where The integral 2.8 can be simplified 16 where x * tan θ * /2 and P 2 W θ * a − b − a − 2b 2 /2b > 0 and conditions 1.9 and 1.11 are satisfied.
Finally, using the change of variables 2.7 and 2.9 , the solution of 1.2 for the heteroclinic orbits, for the region A 0 Figure 2 , can be written as where Now we will consider the region A 1 , including the center c 1 .Let us make a substitution to the new variable in the equation of undisturbed motion 1.2 and obtain the following equation β −a sin β b sin 2β.

2.13
Computing as in 2.3 -2.9 , we obtain the solution of this equation in the form where β * π − θ * .Then, coming back to the variable θ with help of 2.11 , we will receive the equation of the heteroclinic orbits, bounding the region A 1 , including the center c 1 Figure 2 : 2.15

General Positions
Now we set the stage for our study of the disturbed system 1.1 .The stable and unstable manifolds do not necessarily coincide and it is possible that they can cross transversally leading to an infinite number of new heteroclinic points.Then, a heteroclinic tangle is generated.In this case, because of the perturbation, the motion of the system 1.1 , near the unperturbed separatrices, becomes chaotic.Inside this chaotic layer small isolated regions of regular motion with periodic orbits can also appear.The existence of heteroclinic intersections may be proved by means of the Melnikov method 15 .We present a more convenient form the for application of Melnikov method to the disturbed nonautonomous equation of the second order 1.1 as three differential autonomous equations of the first order 2 : φ ω,

3.2
The Melnikov function 2 for system 3.1 is given by where q 0 ± t θ ± t , σ ± t are the solutions of the undisturbed heteroclinic orbits 2.10 or 2.15 for the areas A 0 or A 1 .

Case 1 a 1 0, b 1 0, c 1
The disturbed system 3.1 in this case takes the form where Substituting 3.5 into 3.3 gives where M ε and M δ are the functions corresponding to both perturbations: the external periodic force ε cos ωt and the damping force −δ θ , respectively.The Melnikov function describes the splitting of the stable and unstable manifolds of the disturbed hyperbolic fixed points defined on the cross-section.Thus, there are transverse intersections between the stable and unstable trajectories, if M ± t 0 0. Firstly we consider the functions M

3.7
So using the tabulated integrals 16, 17 , for the integrals 3.7 we obtain the following expressions

International Journal of Mathematics and Mathematical Sciences
Similar expressions can be obtained for M 1 δ and M 1 ε for the region A 1 , including the center c 1 Figure 2 , using the solutions 2.15 where M 1 ε max and M 2 ε max are measures of the maximum splitting of the stable and unstable manifolds, when the disturbed system 3.4 is only under the action of the one perturbation the external periodic force ε cos ωt for the regions A 0 and A 1 , respectively.
Obviously, at a 0 the undisturbed biharmonic oscillator 1.2 is transformed to the simpler system: θ sin 2θ.The regions A 0 and A 1 are equal.From 2.2 , 2.10 , and 2.15 , we obtain

3.11
Following the expressions 3.8 -3.10 , the Melnikov function becomes identical for the regions A 0 and A 1 : From 3.8 -3.10 , it is easy to see that the conditions for the manifolds to intersect in terms of the parameters δ, ε is given by

3.13
Let us define a new parameter of the damping force, divided into amplitude of external force: The criteria Δ j and the variable θ * of parameter a for b −1 and ω 1.
The criteria Δ j as functions of the frequency ω. then conditions 3.13 are given by

3.15
Let us note, that θ * and λ, according to 2.2 and 2.11 , depend on coefficients a and b, therefore criteria 3.15 are functions of the parameters a, b, and ω Δ j Δ j a, b, ω , j 0, 1.

3.16
The criteria 3.16 define chaotic behavior of the perturbed system 3.4 in the regions A 0 and A 1 .In Figure 3, we graph these criteria and the variable θ * as functions of parameter a for the fixed parameter b −1 and ω 1. Figure 4 shows the criteria 3.16 as functions of the frequency ω.

Case 2 a
1 / 0, b 1 / 0, c 0 The disturbed system 3.1 in this case takes the form of φ ω,

3.17
where or where

3.20
For the two regions A 0 and A 1 the functions 3.20 can be represented as where

3.23
It is obvious that using 3.8 integrals 3.23 can be rewritten as The improper integral 3.22 in view of 2.10 and 2.15 is calculated numerically.For parameter 3.14 the conditions for the manifolds to intersect are given by Δ 1 for the area A 1 .

3.25
Criteria 3.25 define behavior of the perturbed system 3.17 in a vicinity of separatrixes.

Numerical Analysis
We have analyzed the evolution of the dynamical behavior of the disturbed system 1.1 as the parameters vary, studying the time histories of the variable θ and its derivative θ.Numerical techniques are based on the numerical integration of the equation of the disturbed motion 1.1 implementing a fixed-step fourth-order Runge-Kutta algorithm.For all numerical calculations the following biharmonic force parameters were used: a 1, b −1 and the frequency of the perturbed force was ω 1.For the numerical analysis of the disturbed system 3.1 we use the P Poincaré cross-section method, examining manifolds with plane sections, perpendicular to the phase axis φ in the two-dimensional space θ, θ , divided with an interval of 2π.It allows us to study the disturbed system 3.1 using a discrete phase instead of examining the continuous dynamics of the system.At ε 0, δ 0 the regular structure of phase space is observed, trajectories have no intersections, and Poincaré sections coincide with undisturbed phase portrait Figure 5 .
Disturbances ε / 0 result in the complication of phase space and the occurrence of a chaotic layer near the undisturbed separatrixes Figures 6-9 .Figures 6 and 7 contain Poincaré sections for the case considered in Section 3.2, and in Figures 8 and 9-in Section 3.3.The growth of disturbances there leads to an increase in the width of the chaotic layer, and the new oscillatory modes determined by closed curves, uncharacteristic for the undisturbed case are observed.
In the presence of damping phase trajectories eventually tended to reach steady positions of equilibrium of the undisturbed system Figures 10, 11 .
In order to check in a quantitative way the validity of the analytic criteria 3.15 we focus on the evolution of the stable and unstable manifolds associated to the saddle fixed points.Parameters a 1, b −1 and the frequency ω 1 in expressions 3.15 give a critical value for the regions A 0 and A 1 : or critical values of the coefficients of a damping force for ε 0.02:            Figure 12 shows numerical simulations of the phase space with initial conditions close to the undisturbed separatrix θ 0 −1.0572, θ0 0.01, φ 0 π/10 for the region A 0 .Now, we reset the value of δ from δ 0 0.01823 to greater ones see Figure 12 .It can be observed clearly that, for δ < δ 0 δ 0.018 , the stable and unstable manifolds transversally intersect each other Figure 12 a .However, when δ > δ 0 δ 0.020 , the invariant manifolds do not intersect Figure 12 b .

Conclusion
This work attempts to describe the transient cases occurring during a spacecraft descent in a planet atmosphere using methods of chaotic mechanics, in particular, the Melnikov method.We have suggested to introduce the concept of a biharmonic system 1.1 which reflects the behavior of the spacecrafts, and also, probably, more general mechanical systems.We have established the existence of transient heteroclinic chaos by means of the Melnikov method.Moreover, this method has provided an analytic criterion for the existence of chaotic behavior in terms of the system parameters.We have found a transition from chaotic to regular regime in the motion of the biharmonic oscillator, as the heteroclinic chaos can be removed by increasing the coefficient of a damping force.The analytic results given by the Melnikov method have been confirmed by a good match with the numeric research.

Figure 2 :
Figure 2: The potential energy W θ a cos θ b cos 2 θ and the phase space for a 1, b −1.

0 δ and M 0 ε
for the area A 0 , including the center c 0 Figure 2 .Substituting 2.10 into 3.6 gives

Figure 13 b
Figure 13 b δ > δ 1 δ 0.0113 .Thus the description, based on numerical simulations for some certain parameter values, makes a good match with the analytic criteria 3.15 provided by Melnikov method.