3. Numerical Method

[23]  There are many numerical methods for solution of fluid dynamics equations. For a numerical solution of equation system (2) one have to choose an optimal approximation method suitable for this problem. Problem solutions of computing geophysical fluid dynamics are very complicated by the fact that here writing the algorithm one has to take into account the boundary conditions, abilities of computers (rate and random access memory), calculation methods, etc. Generally the problem depends on many parameters and one has to choose and adapt a numerical method for the only particular problem.

[24]  The processes considered in this problem are nonstationary. Since the disturbances are very quick variable the approximated equations should be integrated with a small time step. The wave propagation time at large distances is tens of minutes and so thousands of iterations are required. Taking into account this fact and the nonlinearity of the equation system, one can conclude that in this case the explicit finite difference method of integration of fluid dynamics equations is the most suitable. It is known that the background atmosphere density r₀ decreases very sharply with height and large gradients can lead to nonphysical oscillations in the numerical solution. This fact should be taken into account choosing a numerical method [Oran and Boris, 1987]. Taking into account this fact, after testing of various difference methods we created a numerical method, using the some properties of the LCPFCT algorithm [Oran and Boris, 1987]. During the recent 10-15 years this method has been successfully used for solution of various types of fluid dynamics problems.

[25]  At first for numerical solution of equation (2), one has to construct a uniform orthogonal difference grid in which the physical boundaries are located along the grid lines. The difference grid is chosen in such a way that its lower boundary coincides with the rigid flat Earth surface along the OX axis (see Figure 2). For example, the value of r_ik^tn is the value of the density perturbation in the cell center with coordinates x_i and z_k at the time moment tⁿ. Indices i and k denote the number of the calculation cell ( i along the horizontal, i = 1,..., I, and k along the vertical, k = 1,..., K ) and n is the number of the time layer ( n = 0,1,...,N ). Choosing parameters of the difference grid, the Courant-Friedrichs-Lewi conditions should be fulfilled [Durran, 1999], that is for the stability of the solution of equations system (2) the following condition D t [c+u²+w²/ min (D x, D z)] < 1 should be fulfilled. Here D t, D x, and D z are the sizes of the grid steps in time and in the horizontal and vertical directions, respectively. The following parameters of the difference grid were chosen: the steps in altitude, horizontal coordinate, and time are 5 km, 10 km, and 0.5 s, respectively. The height and width of the grid are 400 km and 2000 km, respectively.

[26]  We split the system (2) to two systems, where the first equation system contain in the right-hand side only the vertical components and the second one contains only the components along the horizontal direction:

(9)

(10)

This method increases accuracy of the numerical scheme in respect of time. The obtained systems of one-dimensional equations are separately solved using the one-dimensional module developed for equations system solution in the time interval from t to D t. Numerical integration of the equations is performed in the horizontal and vertical directions alternately. Each couple of sequential integrations is an advancement of the solution by one complete step in time. The application of the numerical method for the continuity equation in (9) gives

Here r_ik^tn denotes an intermediate value of the disturbance in the density. A diffusion part, which provides the stability, should be included in this equation:

where r_ik^tn is also an intermediate value of the disturbance in the density, n is the dimensionless diffusion coefficient. At the next stage this large numerical diffusion is minimized using an adding of antidiffusion fluxes. The final value of the density disturbance in the new time moment n+1 is

where f_k+1/2^ad= m (r_ik+1^tn - r_ik^tn) is the antidiffusion flux, m is the dimensionless coefficient of the antidiffusion. The values of the velocities and temperature disturbances are found by the same method.

[27]  Unlike in the fluxes correction method LCPFCT, we did not use the special limitation applied to the antidiffusion fluxes for conservation of a positivity of solution, because the variables in (2) can be positive and negative values. Thus we get rid of such problems as fluxes synchronization [Oran and Boris, 1987] etc. That is why we wrote the equation system (2) in terms of variations in density and temperature, but not of their absolute values.

[28]  For the sake of stability the nonlinear components in (9) are approximated in the following manner [Durran, 1999]:

Equations system (10), where the components related to the horizontal direction are present in the right-hand side, is solved by the same method.

[29]  Solving equations systems (2) one has to assign initial and boundary conditions. The initial conditions for velocity, density disturbances, and temperature were chosen to be zero. Four walls exist in this model. There are two horizontal and two vertical walls, for which the boundary conditions should be assigned. Points at the boundaries may be considered as inner due to introduction of fictitious cells. The values within these sells are assigned by the boundary conditions. The same boundary conditions as in the case of tangential breaks are applied at the lower boundary [Landau and Lifshits, 1988] that is the variables undertake no jump passing through the Earth surface. At the upper and right-hand side boundaries we applied usual boundary conditions for provision of waves passage through these walls without any considerable reflection [Durran, 1999]. At the left-hand side boundary we used the boundary conditions in form (3), which provide penetration of disturbances into the calculations region. After the pulse passage through the boundaries the same boundary conditions are applied there as for the right-hand side boundary.