2. The Model

2.1. Fluid Dynamics Equations

[14]  We are interested in the wave-like disturbances of the ionosphere electron density observed at large distances from the rocket flight trajectory. These waves (AGW) may be observed from other types of disturbance sources such as explosions, earthquakes, meteors etc. We should simulate generation and propagation of these waves in the lower atmosphere and the ionosphere.

[15]  IGW are the most intense part of the AGW spectrum. Since the 1960s very many papers were dedicated to study of the atmosphere IGW properties. In these papers a solution of the fluid dynamics equations by analytical or numerical methods is one of the widely used methods [Francis, 1975]. New direction in physics of the atmosphere waves was born in the recent years due to the increase of computer processing rate and computational fluid dynamics development. This direction is the study of IGW propagation using numerical solution of the nonlinear geophysical fluid dynamics equations [Zhang and Yi, 2002]. Application of such numerical methods makes it possible to solve numerous problems such as simulation of the intense atmosphere waves excitation at strong impacts of supersonic rocket flights on the environment etc. Thus in order to solve the direct problem, i.e., simulation of an atmosphere disturbance, we have to solve the fluid dynamics equations system with corresponding initial and boundary conditions. Strong inhomogeneous medium exists in the atmosphere around the rocket, being formed by the release of exhausts out of the rocket engine and by the supersonic rocket flight. Therefore it is very complicated and inconvenient to solve the atmosphere disturbances problem during rocket flight in a general form. So the atmosphere region simulated by us is located at some distance from the rocket trajectory, that is, there where the strong SW generated by the rocket flight is attenuated and is transformed into an intense acoustic pulse. Further, there is a nonlinear stretching of this pulse and generation of the AGW wide spectrum with the following disturbance of the ionosphere plasma as result of collision to neutral particles etc. Acoustic pulse is introduced into our model by determination of corresponding boundary conditions. We will come back to the parameters of this pulse in the next paragraph.

[16]  The propagation of AGW in the atmosphere is described by the solution of the fluid dynamics equations system:

(1)

The first, second, third, and fourth equations are continuity equation, momentum conservation equation, energy conservation equation, and equation of the ideal gas state, respectively. The Coriolis force is insignificant for such relatively rapid motions, therefore we neglect by it. Here r, T, p, and v are the density, temperature, pressure, and velocity of environment particles motion, respectively. The values of g, F^d, and Q^d are the gravitational acceleration, viscosity force, and the heat absorbed due to wave dissipation, respectively. The values of c_v, m₀, and R are the specific heat of gas at constant volume, relative molecular mass of air, and the universal gas constant, respectively.

[17]  The viscosity force in this model is introduced as the resistance force in the Rayleigh form F^d = - a v [Sedunov, 1991]. This form of a simple parameterization of the viscous friction forces is often used in complex fluid dynamics calculations. The kinematics friction coefficient c = a /r increasing with height was chosen after the testing of the model at various values. The heat conductivity play the main role at the energy attenuation of such large wave motions [Zhang and Yi, 2002], that is, the dissipative part in the temperature equation is equal to Q^d = k D T, where k is the coefficient of air heat conductivity.

[18]  Each thermodynamic parameter in equations system (1) is split to two parts: the stationary part denoted by index 0 and the disturbed one denoted by a dash:

where U₀ is the horizontal velocity of the background zonal wind in OX direction (the meridional wind is not taken into account). Assuming a hydrostatic equilibrium for the background atmosphere, for the two-dimensional, plane-parallel, compressible atmosphere we obtain after some transformations a system consisting of equations in partial derivatives and the equation of the ideal gas state.

(2)

Figure 2

Here u and w are the horizontal and vertical components of atmosphere particles motion velocity in wave, respectively, and A_1-5 are constants. The frame of reference used here is shown in Figure 2, where the OZ axis is directed vertically upward, and the OX axis is horizontal and lies on the Earth surface.

2.2. Shock Waves

[19]  In order to introduce into our model the boundary conditions describing the disturbance source, one should know the exact parameters of SW generated by the rocket flights. Launching of CR of the space shuttle type is considered for determination of the Mach cone location and orientation in the atmosphere generated by the rocket flight. All launchings of the space shuttle CR have the following stages of flight in the Earth atmosphere. The launching of CR, its ascent to some altitude where the stages are separated and the main engine burns and works during about 7 min [Jacobson and Carlos, 1994]. Only during this stage of the flight (called "the main engine burn" (MEB)), the shuttle and its external tank of liquid fuel continue ascending up to 105-110 km and then accelerate in the horizontal flight. The horizontal flight of MEB extends from 300 km up to 1400 km from the launching point. During the horizontal flight the shuttle velocity increases from 2.5 km h^-1 to 7.5 km h^-1 (up to the end of MEB). Thus the horizontal flight is a supersonic one and whichever waves are generated (acoustic or internal gravity ones), they have to propagate perpendicularly to the flight trajectory. Moreover many other observations show that the generation of AGW occurs mainly during the horizontal rocket flight [Afraimovich et al., 2001].

[20]  Taking the above into account, we introduce the disturbance source into the simulation model in the following manner. Let the vertical coordinate plane XOZ (see Figure 2) be perpendicular to the trajectory of the horizontal rocket flight. One can assumed that a single strong nonlinear acoustic pulse which is further transformed into an AGW packet enters into the computed domain from the left-hand boundary (region AB in Figure 2). To take parameters of this pulse, we consider the general properties of SW from rockets. At the supersonic streamline of the rocket at large distances from it the disturbances caused by SW are weak and therefore may be considered as a cylindrical sound wave divergent from the axis passing through the rocket and parallel to the direction of the streamline [Landau and Lifshits, 1988]. Two shock waves are formed in the cylindrical sound pulse. The velocity in the front break increases by a jump from zero, then there follows the region of a gradual compression decrease changed by a rarefaction, and after that pressure again increases by a jump in the second break. However, the cylindrical sound pulse is specific (as compared with both flat and spherical cases) and can have no rear front: the tendency of the particles motion velocity to zero occurs only asymptotically. From here one can conclude that the cylindrical sound pulse emitted by a rocket flight have a complicated asymmetric form. However, in the first approximation one may take this wave in a sinusoidal form. As it has been noted in section 1, this is also proved by the fact that the initial acoustic wave at large distances is so leveled that its form influences weakly the response form. Running ahead, one may say that the results of our calculations also confirm this fact. As the diameter of the cylindrical pulse is small as compared to the vertical size of the considered region (XOZ), we accepted the front of the wave passing through the region boundary as flat. The acoustic pulse is introduced into the model in such a form that the expression for the horizontal component of the oscillation velocity of air particles on the OZ axis would be

(3)

where u_m, t₀, P, z_m, and D_z are the amplitude, the moment of the pulse arrival, period, altitude of the axis of the sound cylinder from the Earth surface, and the Gauss scale characterizing the cylinder transverse dimension, respectively.

[21]  The values of the density and pressure in the wave can be calculated using formulae describing simple nonlinear acoustic waves [Rudenko and Soluyan, 1975]:

(4)

where g is the adiabatic constant, c is the speed of sound. The value of the temperature perturbation may be determined from the state equation. After that one can easily introduce the boundary conditions into equation (2) for all values at the OZ left-hand side boundary of the simulation region. These boundary conditions provide a transmission of the disturbance from the medium around the rocket trajectory to the atmosphere region considered by us. As for the simulation of the second atmosphere region located to the left from the rocket trajectory, in this case the same boundary conditions will be at the right-hand boundary of the computation domain.

2.3. Ionosphere Disturbances

[22]  Solving equations system (2) one can found the spatial and time distribution of all parameters looked for in the calculation region. As the final goal of our work is a comparison of the simulation results with the data of observations, disturbances in the electron concentration in the ionosphere should be calculated. For this purpose we consider the continuity equation for charged particles:

(5)

Here N_e is the concentration of electrons in undisturbed ionosphere, v_e is the motion velocity of free electrons in the ionosphere, and P_e and L_e are the production and loss rates due to chemical processes, respectively. Andreeva et al. [2001] showed that due to collisions with neutral particles during AGW propagation the ionosphere plasma obtains the velocity (this is most correct for the F layer):

(6)

where b = B/| B| is a unit vector along the Earth magnetic field (see Figure 2). We assume that the magnetic field is homogeneous in the simulation domain. Taking into account equation (6) and neglecting formation and loss of charged particles, after integration of (5) we obtain

(7)

Figure 3

where b_x and b_z are the components of a unit vector of the Earth magnetic field. Here the first integral characterizes the variation in the concentration due to the motion of some volume of the ionosphere and the second integral is a consequence of the processes of compression or rarefaction in the plasma. Using this formula one can calculate the variations in the electron concentration in the ionosphere at AGW propagation at this point in the given moment of time. In our calculations we took the profile (see Figure 3) consisting of two layers as the background ionosphere electron concentration. We used the IGRF model (http://nssdc.gsfc.nasa.gov/space/model/models/igrf.html) for determination of the geomagnetic field components. If one neglects the magnetic field influence, that is everywhere b v, then

(8)

As noted in section 1, the determination of TEC variation ( D TEC) in the ionosphere in various directions, i.e., between the receiver and GPS satellite (see Figure 2), is of a great importance for observations of the ionosphere state:

Integrating, we assume that above the simulation domain D N_e 0. It should be noted that the real ionosphere is horizontally inhomogeneous but we confine our consideration by a parallel-sided ionosphere.