Finite Element Solution of a Problem for Gravity Gyroscopic Equation in the Time Domain

To solve the equation for gravity-gyroscopic waves in a rectangular domain, the distinguished algorithm for the solution of the Cauchy problem for a second-order transient equation is proposed. This algorithm is developed by using the time-varying finite element method. The space derivatives in the gravity-gyroscopic wave equation are approximated with finite differences. The stability and accuracy of the method are estimated. The procedure for the implementation of the method is developed. The calculations were performed for determining the steady-state modes of fluctuations of the solutions of the gravity-gyroscopic wave equation depending on the problem parameters.


Introduction
The strengthening of requirements for solving complex problems in scientific and technological researches results in the improvement of existing numerical computation algorithms.This improvement is of particular importance in studying complex nonstationary processes, for example, the propagation of gravity-gyroscopic waves.
The mathematical statement of problems in the theory of internal waves, including gravity-gyroscopic waves, and methods for solving such problems are discussed in publications [1]- [6].In publications [1]- [3], the equations for describing the propagation of internal waves are derived using the Bussinesk approximation.Various analytical methods of solving problems in the linear theory of gravity-gyroscopic waves are proposed in publications [3]- [6].Specifically, in publications [4] [5], the modes of steady-state oscillations associated with gravity-gyroscopic waves are defined depending on the coefficients of equations describing the propagation of gravity-gyroscopic waves in response to a harmonic disturbance in the unbounded medium.
Since analytical methods not always allow solutions of problems for arbitrary data and solution regions to be obtained, numerical methods are used.In solving nonstationary problems, semidiscrete methods are used, that is, the methods with initial approximation only for spatial variables followed by the solution of the Cauchy problem by using finite-difference algorithms.These methods provide an order of accuracy not higher than the second one [7] [8].
To enhance the calculation accuracy, it is necessary to use numerical algorithms providing a higher order of accuracy [9].This is possible due to the finite element method.Among numerous publications for the finite element method, we refer to publications [10]- [12] which specify the methods and algorithms for solving stationary and nonstationary problems in the mathematical physics.In publication [13], an algorithm for solving second-order nonstationary equations, which provides the fourth order of accuracy, is proposed.According to publications [14]- [18], this algorithm was used and substantiated in solving various problems in the mathematical physics which are described by equations for oscillations.The methods for the implementation of this algorithm were developed, as well as the stability and convergence of the algorithm in the classes of Soboliev functions were proved.

Problem Statement
Consider the following problem: where Equation ( 1) is obtained from the set of hydrodynamic equations in Boussinesq approximation [1]- [5].This equation describes small oscillations of stratified fluid with density ( ) x ρ and α being twice the speed of fluid rotation around 3 ox axis.The quantity ( ) ( ) ( ) denotes the square of the Brent-Väisälä frequency.The fluid velocity is calculated using the formula We will also assume that the Brent-Väisälä frequency 0 ω is constant.Equation ( 1) is not classical and belongs to the Sobolev equations of composite type [1]- [3].
Let us define the generalized solution of problem (1)-(3) as a function ( ) ( ) , that has a derivative ( ) and satisfies the following relationships for [ ] with values in H.The issues of the solution existence and its properties were discussed in the works [1]- [3].

Spatial and Time Approximation
Let us construct a subspace h H H ⊂ that approximates H. Let us also introduce a mesh . In this case ( ) -is a space of mesh functions ( ) with the norm ( ) ( ) ( ) , where constant M is independent of 1 3 , h h [7].
With the use of the corresponding quadrature formulae ( ) , m a u ϑ on the mesh and derive an approximated mesh solution from Equation ( 4): ( The problem defined by Equation ( 5) corresponds to the following Cauchy problem for the function where ( ) ( ) ( ) , , .
Let us seek an approximate solution of problem (6) in the form of the Hermitian spline of the third order si-milarly to [10] [11]: Selecting the various weight functions ( ) t θ , namely: ( ) , the following two parameter vector difference scheme can be obtained from Equation (9) [13] [14]:  The parameters , α β are subject to the fourth order approximation condition for Scheme (10) [14].

Scheme Implementation Algorithm
For implementation of Scheme (10), it is necessary to solve the system of two equations with respect to unknown variables , y y  : Let us assume that condition (11) holds.Excluding ŷ  from Equations ( 12) we obtain an equation to find ŷ : The first algorithm of implementation of Scheme ( 12) is a direct solving Equation ( 13) with further calculation of ŷ  from the first equation in System (12): ( ) The second algorithm of implementation of Scheme ( 12) is reduced to solving the following two equations [13] [14] for each time increment: For instance, if 1 8, 1 24 For inversion of the matrices 1 2 , C C the direct square root method was used once at the initial time.For all other layers the solutions were obtained by multiplying the matrix
then the solution ( ) y t of scheme (10) converges to the solution of Problem (5) and the following estimate holds:

( ) ( ) ( ) ( )
The proof of the statement is based on the separate transformation of two layer scheme (10) to a three layer scheme for both solution y and its derivative y  .Condition (14) for the chosen values of the parameters ( 1 8 α = , 1 24 and the corresponding value 1 8 δ = ) gives the following restriction on the time increment: For evaluation of the accuracy of Scheme (10) the error h z u u = − shall be estimated.Applying the procedure used in the theory of finite difference schemes for such estimates [8] the final result can be formulated.

Computational Experiments
Let us proceed to the description of the computational experiment.The asymptotic amplitude problem for the solution of Equation ( 1) under the uniform initial conditions on the plane was studied in the works [2] [3].As a source of excitation the oscillations in the vicinity of a thin plate of 2 in length inclined at angle θ with 1 ox axis were assumed.The law of plate oscillations as a rigid body im- mersed in a liquid is as follows: ( ) ( ) ( ) > Equation ( 16) is an elliptic one.Under the condition of ( ) ( ) 16) is of the hyperbolic type.It means that there are differences in the behavior of the solution ( ) , , u x x t in its way to the steady state.Furthermore, it has been shown that when t → ∞ the solution decays according the following law Let us consider the right-hand member of Equation ( 1) as an excitation factor of the fluid motion that is acting on a small segment of 2∆ in length which is small compared to the overall length of the domain:

(
) ( ) ( ) x and during the finite period of time: The initial conditions are uniform: The dimensions of the domain are 1 2 10 l l = = , 1 ∆ = and 0 10 t ≤ ≤ .The mesh parameters are as follows: For the first case 0 ω α ω = < and so 2 1 b = .This result is trivial in a certain sense.In this case, Equation (1) can be written as an ordinary differential equation: ( ) Figure 2 shows the plot of the amplitude oscillations at the middle point of the domain Ω , i.e. at 1 1 2 x l = ,  The next plots (Figures 5-7) show the solutions for a hyperbolic case (the problem parameters are described in the caption of Figure 5).
, , y x x t for 50 t = .
, , y x x t for 47.5 t = .
, , y x x t for 49.5 t = .
, , y x x t for 36.5 t = .
, , y x x t for 36.5 t = .
, , y x x t for 43 t = .

Conclusions
The following conclusions can be drawn on the basis of the comparative analysis of the simulation results: • for all parameter combinations the solution achieves a steady state regime (Figures 3-16); • there are no significant differences between the parameter sets of hyperbolic or elliptic types (Figure 6, Figure 7 and , x x .However, along the entire plane the energy has no constraints in the form of the domain boundaries thus it is dissipated along the infinite plane ( ) , x x leading to the following behavior of the solution as a function of t, as was shown in [3]:

Theorem 2 .
If stability condition(15) of Scheme(10) holds the solution ( ) y t converges to a sufficiently smooth solution of Problem (1)-(3) and the following estimate holds:

ωFigure 1 .
the simulations are shown in Figures1-16.For all simulations the frequency of induced osciland α such values were used at which the type of Equation (16) was el- liptic or hyperbolic or singular due to fulfillment of the following equalities: either 0 ω ω = , or α ω = .The plot of the time variation of the amplitude The plot shows that the time of the external influence on the medium is limited by the segment 0 8 t ≤ ≤ for all alternative simulations.The simulation plots for five different combinations of the parameters 0 , , ω ω α that have resulted in the in- teresting solutions with a typical behavior (as follows from the results not included in this work) are shown below.Let us introduce the parameter the solution of the Dirichlet problem for the Poisson elliptic equation:

Figure 3 and
Figure4show the simulation plots for two instants of time with the positive and negative solution values.

Figure 9 ,Figure 15 ,•
Figure 10): in both cases the form of the solution form is determined by the shape of spatial domain; • the singular cases show that the steady state regime corresponds to equation of asymptotic amplitude (14): the solution shown in Figure 12, Figure 13 is "stretched" along ox 1 axis in the middle region and corresponds to the asymptotic equation ( ) Figure 16 the solution is"stretched" along the ox 3 axis and corresponds to the limiting equation, for t → ∞ no decay of the solution amplitude ( ) , u x t is observed.The last observation can be explained as follows.For the solution of problem (1)-(3) the following law of energy conservation holds: can be obtained from identity (4) by simple manipulations assuming u ϑ = .The same conservation law holds for the Cauchy problem solution of Equation (1) in the entire plane ( )