Chaotic Behavior of the Biharmonic Dynamics System

Abstract

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 𝑎sin𝜃+𝑏sin2𝜃. 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.

1. Introduction

1.1. 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 two-frequency 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: 𝜃̈=𝑎sin𝜃+𝑏sin2𝜃+𝜀(𝑎₁sin𝜃+𝑏₁sin2𝜃+𝑐)cos𝜔𝑡−𝛿𝜃̇, (1.1) where 𝜀 and 𝛿 are assumed to be small positive parameters; 𝜔>0 is the frequency of the external force; and 𝑎, 𝑎₁, 𝑏, 𝑏₁, and 𝑐 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 𝜃̈=𝑎sin𝜃+𝑏sin2𝜃. (1.2) We will show that (1.1) is equivalent to the Duffing equation for small values of 𝜃. It is known that 𝑥sin𝑥=𝑥−³3!+⋯(𝑥=𝜃,2𝜃), (1.3) and (1.1) can be written as 𝜃̈−𝑎(𝜃−𝜃³3!)−𝑏[2𝜃−(2𝜃)³𝜃3!]=𝜀{𝑎(𝜃−³3!)+𝑏[2𝜃−(2𝜃)³̇3!]+𝑐}cos𝜔𝑡−𝛿𝜃+⋯. (1.4) If the variable 𝜃 is a small of order 𝜀, then we have the Duffing equation 𝜃̈+𝜆𝜃+𝜇𝜃³=𝜀(𝑐−𝜆𝜃)cos𝜔𝑡−𝛿𝜃̇+𝑂(𝜀⁴), (1.5) where 𝜆=−(𝑎+2𝑏), 𝜇=(𝑎+8𝑏)/3!.

1.2. 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 re-entry 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] 𝑚_𝜃(𝜃)=𝑚_𝑂(𝜃)+𝑐_𝑛(𝜃)𝑥_𝑐, (1.6) where 𝜃 is the spatial angle of attack (𝜃 is defined as the angle between symmetry axis and the velocity of the spacecraft 𝑉, Figure 1), 𝑚_𝑂 is the pitching moment coefficient concerning the leading edge body, 𝑐_𝑛 is the normal force coefficient, 𝑥_𝑐=𝑥_𝑐/𝐿 is the nondimensional coordinate of the center of mass, and 𝐿 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) 𝑚_𝜃=𝑎sin𝜃+𝑏sin2𝜃. (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 𝑚_𝜃(𝜃) with respect to the angle of attack 𝜃 at this point is negative: 𝑚_𝜃′|_𝜃=0=(𝑎cos𝜃+2𝑏cos2𝜃)|_𝜃=0<0 (1.8) or 2𝑏<−𝑎, (1.9) and if there exists an intermediate position of equilibrium inside the interval of (0,𝜋), then 𝑚_𝜃(𝜃)=sin𝜃(𝑎+2𝑏cos𝜃)=0, (1.10) which holds true, if |2𝑏|>|𝑎|. (1.11) It is obvious that (1.9) and (1.10) are valid simultaneously when 𝑏<0. Note that the dependence of 𝑚_𝜃(𝜃) given in Figure 1 satisfies conditions (1.9) and (1.10).

Figure 1 

A blunt-shaped spacecraft.

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 𝜃=𝜋 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 𝑀_𝜃=𝑚_𝜃(𝜃)𝑞𝑆𝐿, (1.12) where 𝑆 is reference surface and 𝑞 is the dynamic pressure. Disturbance will simulate periodic change of position of the center of mass 𝑥_𝑐=𝑥_𝑐0+𝜀cos𝜔𝑡, (1.13) where 𝑥_𝑐0 is an initial position of the center of mass. Using expression (1.6), we represent the aerodynamic restoring moment as 𝑀_𝜃=𝐴sin𝜃+𝐵sin2𝜃+𝜀(𝐴sin𝜃+𝐵sin2𝜃)cos𝜔𝑡. (1.14) Thus, the planar motion of the spacecraft about the center of mass can be described by 𝐼𝜃̈=𝐴sin𝜃+𝐵sin2𝜃+𝜀(𝐴sin𝜃+𝐵sin2𝜃)cos𝜔𝑡+𝑀_𝑑, (1.15) where 𝐼 is the transverse moment of inertia, 𝑀_𝑑=−𝛿₁𝜃̇ is small damping moment. Obviously, the equation of spacecraft motion (1.15) corresponds to the biharmonic dynamics system (1.1) at 𝑐=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] 𝜃̈=𝑎sin𝜃+𝑏sin2𝜃(𝑎=−𝑔𝑙,𝑏=Ω²>0), (1.16) where 𝑔 is the gravitational acceleration, 𝑙 is the radius of the circle and Ω is the angular velocity of the circle.

1.3. 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.

2. 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 𝑏=0 or 𝑎=0 (replacement of variables 𝜑=2𝜃) we have the equation of a mathematical pendulum. However, if (1.9) and (1.11) are satisfied (𝑏<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 𝜃_∗𝑎=±arccos(−2𝑏), (2.2) where 𝑏<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 𝐴₀ and 𝐴₁, divided by the two saddles 𝑠₁ and 𝑠₋₁ (Figure 2). It is necessary to note, that the region 𝐴₁ of the development of the cylinder undergoes a break at 𝜃=𝜋,−𝜋. From (2.2) it follows, that if the coefficient 𝑎 is equal 0, the saddle 𝑠₁ is in the position: 𝜃_∗=𝜋/2. At positive values of the coefficient 𝑎>0 the saddle 𝑠₁ belongs to the interval: 𝜃_∗∈(0,𝜋/2), and at negative values 𝑎<0 the saddle 𝑠₁ belongs to the interval: 𝜃_∗∈(𝜋/2,𝜋) (Figure 2).

(a)

(b)

Figure 2 

The potential energy 𝑊(𝜃)=𝑎cos𝜃+𝑏cos2𝜃 and the phase space for 𝑎=1, 𝑏=−1.

The following energy integral corresponds to (1.2): 12𝜃̇²+𝑎cos𝜃+𝑏cos²𝜃=𝐸, (2.3) where 𝐸 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: 𝑊(𝜃)=𝑎cos𝜃+𝑏cos²𝜃. (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 𝐸>𝑊_∗, where 𝑊_∗=𝑊(𝜃_∗), then the motion is possible in the outer regions (Figure 2). In the opposite case (𝐸<𝑊_∗) the motion can occur in any of the inner regions, depending on initial conditions. The equality 𝐸=𝑊_∗ corresponds to the motion along separatrixes. In this case, the two saddles 𝑠₁ and 𝑠₋₁ are connected by four heteroclinic trajectories.

First of all, we will consider the separatrixes, limiting the region 𝐴₀. 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 𝑡−𝑡₀=^𝜃_𝜃0𝑑𝜃2[𝑊(𝜃_∗)−𝑎cos𝜃−𝑏cos²𝜃], (2.5) where 𝑊(𝜃_∗)=𝑎cos𝜃_∗+𝑏cos²𝜃_∗𝑎=−²4𝑏. (2.6) Substituting the variables: 𝜃𝑥=tan2, (2.7) we can rewrite (2.5) as 𝑡−𝑡₀=2^𝑥_𝑥0𝑑𝑥𝜓(𝑥),𝜓(𝑥)=𝑃𝑥⁴+𝐶𝑥²+𝐴, (2.8) where 𝑃=2[𝑊(𝜃_∗)+𝑎−𝑏], 𝐶=4[𝑊(𝜃_∗)+𝑏], 𝐴=2[𝑊(𝜃_∗)−𝑎−𝑏].

The integral (2.8) can be simplified [16] 𝑡−𝑡₀=2√𝑃^𝑥𝑥₀𝑑𝑥𝑥²_∗−𝑥²=2√𝑃𝑥ln|^∗+𝑥𝑥_∗−𝑥||^𝑥_𝑥0, (2.9) where 𝑥_∗=tan(𝜃_∗/2) and 𝑃=2[𝑊(𝜃_∗)+𝑎−𝑏]=−(𝑎−2𝑏)²/2𝑏>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 𝐴₀  (Figure 2), can be written as 𝜃₊𝜃(𝑡)=2arctan[tan^∗2tanh(𝜆𝑡2)],𝜎₊(𝑡)=(𝜃̇)=𝜆sin𝜃_∗cosh(𝜆𝑡)+cos𝜃_∗,[𝜃₋(𝑡),𝜎₋(𝑡)]=[−𝜃₊(𝑡),−𝜎₊(𝑡)], (2.10) where 𝜃𝜆=tan^∗2√𝑃=𝑎²−4𝑏²2𝑏. (2.11) Now we will consider the region 𝐴₁, including the center 𝑐₁. Let us make a substitution to the new variable 𝛽=𝜋−𝜃 (2.12) in the equation of undisturbed motion (1.2) and obtain the following equation 𝛽̈=−𝑎sin𝛽+𝑏sin2𝛽. (2.13) Computing as in (2.3)–(2.9), we obtain the solution of this equation in the form 𝛽𝛽=2arctan{tan^∗2tanh[𝜆(𝑡−𝑡₀)2]}, (2.14) where 𝛽_∗=𝜋−𝜃_∗. Then, coming back to the variable 𝜃 with help of (2.11), we will receive the equation of the heteroclinic orbits, bounding the region 𝐴₁, including the center 𝑐₁ (Figure 2): 𝜃₊𝜃(𝑡)=𝜋−2arctan[cot^∗2tanh(𝜆𝑡2)],𝜎₊(𝑡)=(𝜃̇)=𝜆sin𝜃_∗cosh(𝜆𝑡)−cos𝜃_∗,[𝜃₋(𝑡),𝜎₋(𝑡)]=[2𝜋−𝜃₊(𝑡),−𝜎₊(𝑡)]. (2.15)

3. Chaotic Motion. The Melnikov Criterion

3.1. 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]: 𝜃̇=𝜎=𝑓₁+𝑔₁,𝜎̇=𝑎sin𝜃+𝑏sin2𝜃+𝜀(𝑎₁sin𝜃+𝑏₁sin2𝜃+𝑐)cos𝜙−𝛿𝜎=𝑓₂+𝑔₂,𝜙̇=𝜔, (3.1) where 𝑓₁=𝜎,𝑔₁=0,𝑓₂𝑔=𝑎sin𝜃+𝑏sin2𝜃,₂=𝜀(𝑎₁sin𝜃+𝑏₁sin2𝜃+𝑐)cos𝜙−𝛿𝜎. (3.2) The Melnikov function [2] for system (3.1) is given by 𝑀^±(𝑡₀,𝜙₀)=^∞_−∞{𝑓₁[𝑞⁰_±(𝑡)]𝑔₂[𝑞⁰_±(𝑡),𝜔𝑡+𝜔𝑡₀+𝜙₀]−𝑓₂[𝑞⁰_±(𝑡)]𝑔₁[𝑞⁰_±(𝑡),𝜔𝑡+𝜔𝑡₀+𝜙₀=]}𝑑𝑡_∞−∞{𝑓₁[𝑞⁰_±(𝑡)]𝑔₂[𝑞⁰_±(𝑡),𝜔𝑡+𝜔𝑡₀+𝜙₀]}𝑑𝑡, (3.3) where 𝑞⁰_±(𝑡)=[𝜃_±(𝑡),𝜎_±(𝑡)] are the solutions of the undisturbed heteroclinic orbits (2.10) or (2.15) for the areas 𝐴₀ or 𝐴₁.

3.2. Case 1 (𝑎₁=0,𝑏₁=0,𝑐=1)

The disturbed system (3.1) in this case takes the form 𝜃̇=𝜎=𝑓₁+𝑔₁,𝜎̇=𝑎sin𝜃+𝑏sin2𝜃+𝜀cos𝜙−𝛿𝜎=𝑓₂+𝑔₂,𝜙̇=𝜔, (3.4) where 𝑓₁=𝜎,𝑔₁=0,𝑓₂=𝑎sin𝜃+𝑏sin2𝜃,𝑔₂=𝜀cos𝜙−𝛿𝜎. (3.5) Substituting (3.5) into (3.3) gives 𝑀^±(𝑡₀,𝜙₀)=^∞_−∞𝜎_±{𝜀cos(𝜔𝑡+𝜔𝑡₀+𝜙₀)−𝛿𝜎_±}𝑑𝑡=𝜀^∞_−∞𝜎_±cos(𝜔𝑡+𝜔𝑡₀+𝜙₀)𝑑𝑡−𝛿^∞_−∞(𝜎_±)²𝑑𝑡=𝑀_𝜀+𝑀_𝛿, (3.6) where 𝑀_𝜀 and 𝑀_𝛿 are the functions corresponding to both perturbations: the external periodic force (𝜀cos𝜔𝑡) 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 𝑀^±(𝑡₀)=0.

Firstly we consider the functions 𝑀_𝛿⁽⁰⁾ and 𝑀_𝜀⁽⁰⁾ for the area 𝐴₀, including the center 𝑐₀ (Figure 2). Substituting (2.10) into (3.6) gives 𝑀_𝛿⁽⁰⁾=−𝛿^∞_−∞(𝜎_±)²𝑑𝑡=−𝛿sin²𝜃_𝑆^∞_−∞𝑑𝑡[cosh(𝜆𝑡)+cos𝜃_∗]²,𝑀_𝜀⁽⁰⁾(𝑡₀,𝜙₀)=𝜀^∞_−∞𝜎_±cos(𝜔𝑡+𝜔𝑡₀+𝜙₀)𝑑𝑡=𝜀𝜆sin𝜃_𝑆^∞_−∞cos(𝜔𝑡+𝜔𝑡₀+𝜙₀)cosh(𝜆𝑡)+cos𝜃_∗𝑑𝑡. (3.7) So using the tabulated integrals [16, 17], for the integrals (3.7) we obtain the following expressions 𝑀_𝛿⁽⁰⁾=−2𝛿𝜆(1−𝜃_∗cot𝜃_∗𝑀),_𝜀(0)(𝑡₀,𝜙₀)=2𝜀𝜋sinh(𝜃_∗(𝜔/𝜆))𝜆sin(𝜃_∗)sinh(𝜋(𝜔/𝜆))cos(𝜔𝑡₀+𝜙₀)=𝑀⁽⁰⁾_𝜀maxcos(𝜔𝑡₀+𝜙₀). (3.8) Similar expressions can be obtained for 𝑀_𝛿⁽¹⁾ and 𝑀_𝜀⁽¹⁾ for the region 𝐴₁, including the center 𝑐₁  (Figure 2), using the solutions (2.15) 𝑀_𝛿⁽¹⁾=−𝛿^∞_−∞(𝜎_±)²𝑑𝑡=−𝛿sin²𝜃_∗^∞_−∞𝑑𝑡[cosh(𝜆𝑡)−cos𝜃_∗]²,𝑀_𝜀⁽¹⁾(𝑡₀,𝜙₀)=𝜀^∞_−∞𝜎_±cos(𝜔𝑡+𝜔𝑡₀+𝜙₀)𝑑𝑡=𝜀𝜆sin𝜃_𝑆^∞_−∞cos(𝜔𝑡+𝜔𝑡₀+𝜙₀)cosh(𝜆𝑡)−cos𝜃_∗𝑑𝑡, (3.9) or 𝑀_𝛿⁽¹⁾=−2𝛿𝜆[1+(𝜋−𝜃_∗)cot𝜃_∗𝑀],_𝜀(1)(𝑡₀,𝜙₀)=2𝜀𝜋sinh[(𝜋−𝜃_∗)(𝜔/𝜆)]𝜆sin(𝜃_∗)sinh(𝜋(𝜔/𝜆))cos(𝜔𝑡₀+𝜙₀)=𝑀⁽¹⁾_𝜀maxcos(𝜔𝑡₀+𝜙₀), (3.10) where 𝑀⁽¹⁾_𝜀max and 𝑀⁽²⁾_𝜀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𝜔𝑡) for the regions 𝐴₀ and 𝐴₁, respectively.

Obviously, at 𝑎=0 the undisturbed biharmonic oscillator (1.2) is transformed to the simpler system: 𝜃̈=sin2𝜃. The regions 𝐴₀ and 𝐴₁ are equal. From (2.2), (2.10), and (2.15), we obtain 𝜃_∗=𝜋2,𝜆=√−2𝑏. (3.11) Following the expressions (3.8)–(3.10), the Melnikov function becomes identical for the regions 𝐴₀ and 𝐴₁: 𝑀^±(𝑡₀,𝜙₀)=𝑀_𝛿+𝑀_𝜀(𝑡₀,𝜙₀𝜋)=−2𝛿𝜆+𝜀𝜆𝜋secℎ(2𝜔𝜆)cos(𝜔𝑡₀+𝜙₀). (3.12)

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 𝛿<[𝜋sinh(𝜃_∗(𝜔/𝜆))𝜆²(1−𝜃_𝑆cot𝜃_∗)sin(𝜃_∗)sinh(𝜋(𝜔/𝜆))]𝜀(forthearea𝐴₀),𝛿<{𝜋sinh[(𝜋−𝜃_∗)(𝜔/𝜆)]𝜆²[1+(𝜋−𝜃_∗)cot𝜃_∗]sin(𝜃_∗)sinh(𝜋(𝜔/𝜆))}𝜀(forthearea𝐴₁). (3.13) Let us define a new parameter of the damping force, divided into amplitude of external force: 𝛿Δ=𝜀, (3.14) then conditions (3.13) are given by Δ<𝜋sinh(𝜃_∗(𝜔/𝜆))𝜆²(1−𝜃_∗cot𝜃_∗)sin(𝜃_∗)sinh(𝜋(𝜔/𝜆))=Δ₀(forthearea𝐴₀),Δ<𝜋sinh[(𝜋−𝜃_∗)(𝜔/𝜆)]𝜆²[1+(𝜋−𝜃_∗)cot𝜃_∗]sin(𝜃_∗)sinh(𝜋(𝜔/𝜆))=Δ₁(forthearea𝐴₁). (3.15) Let us note, that 𝜃_∗ and 𝜆, according to (2.2) and (2.11), depend on coefficients 𝑎 and 𝑏, therefore criteria (3.15) are functions of the parameters 𝑎, 𝑏, and 𝜔Δ_𝑗=Δ_𝑗(𝑎,𝑏,𝜔),𝑗=0,1. (3.16) The criteria (3.16) define chaotic behavior of the perturbed system (3.4) in the regions 𝐴₀ and 𝐴₁. In Figure 3, we graph these criteria and the variable 𝜃_∗ as functions of parameter 𝑎 for the fixed parameter 𝑏=−1 and 𝜔=1. Figure 4 shows the criteria (3.16) as functions of the frequency 𝜔.

Figure 3 

The criteria Δ𝑗 and the variable 𝜃∗ of parameter 𝑎 for 𝑏=−1 and 𝜔=1.

Figure 4 

The criteria Δ𝑗 as functions of the frequency 𝜔.

3.3. Case 2 (𝑎₁≠0,𝑏₁≠0,𝑐=0)

The disturbed system (3.1) in this case takes the form of 𝜃̇=𝜎=𝑓₁+𝑔₁,𝜎̇=𝑎sin𝜃+𝑏sin2𝜃+𝜀(𝑎₁sin𝜃+𝑏₁sin2𝜃)cos𝜙−𝛿𝜎=𝑓₂+𝑔₂,𝜙̇=𝜔, (3.17) where 𝑓₁=𝜎, 𝑔₁=0, 𝑓₂=𝑎sin𝜃+𝑏sin2𝜃, 𝑔₂=𝜀(𝑎₁sin𝜃+𝑏₁sin2𝜃)cos𝜙−𝛿𝜎, 𝑀^±(𝑡₀,𝜙₀)=^∞_−∞𝜎_±[𝜀(𝑎₁sin𝜃_±+𝑏₁sin2𝜃_±)cos(𝜔𝑡+𝜔𝑡₀+𝜙₀)−𝛿𝜎_±]𝑑𝑡 (3.18) or 𝑀^±(𝑡₀,𝜙₀)=𝑀_𝜀+𝑀_𝛿, (3.19) where 𝑀_𝜀=𝜀^∞_−∞𝜎_±(𝑎₁sin𝜃_±+𝑏₁sin2𝜃_±)cos(𝜔𝑡+𝜔𝑡₀+𝜙₀𝑀)𝑑𝑡,_𝛿=−𝛿^∞_−∞(𝜎_±)²𝑑𝑡. (3.20) For the two regions 𝐴₀ and 𝐴₁ the functions (3.20) can be represented as 𝑀_𝜀^(𝑘)=−𝜀𝐼_±^(𝑘)sin(𝜔𝑡₀+𝜙₀𝑀),_𝛿(𝑘)=−𝛿𝐽_±^(𝑘),𝑘=0,1, (3.21) where 𝐼_±^(𝑘)=^∞_−∞𝜎_±^(𝑘)(𝑎₁sin𝜃_±^(𝑘)+𝑏₁sin2𝜃_±^(𝑘)𝐽)sin(𝜔𝑡)𝑑𝑡,_±(𝑘)=^∞_−∞(𝜎_±^(𝑘))²𝑑𝑡. (3.22)(3.23) It is obvious that using (3.8) integrals (3.23) can be rewritten as 𝐽_±⁽⁰⁾=2𝜆(1−𝜃_∗cot𝜃_∗),𝐽_±⁽¹⁾=2𝜆[1+(𝜋−𝜃_∗)cot𝜃_∗], (3.24) 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 𝐼Δ<^±(0)𝐽_±⁽⁰⁾=Δ₀(forthearea𝐴₀𝐼),Δ<^±(1)𝐽_±⁽¹⁾=Δ₁(forthearea𝐴₁). (3.25) Criteria (3.25) define behavior of the perturbed system (3.17) in a vicinity of separatrixes.

4. 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: 𝑎=1, 𝑏=−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).

Figure 5 

Poincaré sections for 𝜀=0, 𝛿=0.

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.

Figure 6 

Poincaré sections in the case of 𝑎1=0,𝑏1=0,𝑐=1 for 𝜀=0.01, 𝛿=0.

Figure 7 

Poincaré sections in the case of 𝑎1=0,𝑏1=0,𝑐=1 for 𝜀=0.02, 𝛿=0.

Figure 8 

Poincaré sections in the case of 𝑎1=1,𝑏1=−1,𝑐=0 for 𝜀=0.01, 𝛿=0.

Figure 9 

Poincaré sections in the case of 𝑎1=1,𝑏1=−1,𝑐=0 for 𝜀=0.02, 𝛿=0.

In the presence of damping phase trajectories eventually tended to reach steady positions of equilibrium of the undisturbed system (Figures 10, 11).

Figure 10 

Poincaré sections in the case of 𝑎1=0,𝑏1=0,𝑐=1 for 𝜀=0.02, 𝛿=0.0001.

Figure 11 

Poincaré sections in the case of 𝑎1=1,𝑏1=−1,𝑐=0 for 𝜀=0.02, 𝛿=0.0001.

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 𝑎=1, 𝑏=−1 and the frequency 𝜔=1 in expressions (3.15) give a critical value for the regions 𝐴₀ and 𝐴₁: Δ₀=0.9115,Δ₁=0.4530, (4.1) or critical values of the coefficients of a damping force for 𝜀=0.02: 𝛿₀=𝜀Δ₀=0.01823,𝛿₁=𝜀Δ₁=0.00906. (4.2) Figure 12 shows numerical simulations of the phase space with initial conditions close to the undisturbed separatrix (𝜃₀=−1.0572, 𝜃̇₀=0.01, 𝜙₀=𝜋/10) for the region 𝐴₀. Now, we reset the value of 𝛿 from 𝛿₀=0.01823 to greater ones (see Figure 12). It can be observed clearly that, for 𝛿<𝛿₀  (𝛿=0.018), the stable and unstable manifolds transversally intersect each other (Figure 12(a)). However, when 𝛿>𝛿₀  (𝛿=0.020), the invariant manifolds do not intersect (Figure 12(b)).

(a)

(b)

Figure 12 

Evolution of stable and unstable manifolds as a function of the damping force coefficients 𝛿 for 𝑎=1, 𝑏=−1, 𝜀=0.02, 𝜔=1 and two different values of 𝛿 close to the critical value 𝛿0=0.01823. (a) 𝛿=0.018 (b) 𝛿=0.020, with following initial conditions: 𝜃0=−1.0572, 𝜃̇0=0.01, 𝜙0=𝜋/10 for the area 𝐴0.

Similar results for the region 𝐴₁ are shown in Figure 13 (𝛿₁=0.00906) for the following initial conditions: 𝜃₀=0.9472, 𝜃̇₀=0.2, 𝜙₀=𝜋. Figure 13(a)  𝛿<𝛿₁  (𝛿=0.0090). Figure 13(b)  𝛿>𝛿₁  (𝛿=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.

(a)

(b)

Figure 13 

Evolution of stable and unstable manifolds as a function of the damping force coefficients 𝛿 for 𝑎=1, 𝑏=−1, 𝜀=0.02, 𝜔=1 and two different values of 𝛿 close to the critical value 𝛿1=0.00906. (a) 𝛿=0.0090 (b) 𝛿=0.0113, with following initial conditions: 𝜃0=0.9472, 𝜃̇0=0.2, 𝜙0=𝜋 for the area 𝐴1.

5. 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.

Acknowledgment

The author would like to thank the Russian Foundation for Basic Research (Project no. 06-01-00355a), which has supported this work.