Numerical comparison models

Numerical schemes

Our goal is to compare three models, the symmetric Boussinesq one, the usual KdV approximation (M) and its topographically modified version (Mb). The comparison is made for a solitary wave propagating to the right above two topographies : a step and a slowly varying sinusoidal bottom. We use for the Boussinesq system (Σ) and the KdV equations (ΣKdV ) a Crank-Nicholson scheme combined with a relaxation method coming from BesseBruneau in [10] and justified by Besse in [8]. This type of scheme is of order two in space and time, which is appropriate for our purpose. 6.1.1 Numerical scheme for the KdV equations Due to the identical structure of the two KdV equations of (ΣKdV ), we only present the numerical scheme for the first equation. Defining u(t, x) = U0(T, x−t), we can reformulate this equation as follows ∂tu + ∂xu + ε · 3 4 u∂xu + 1 6 ∂ 3 xu ¸ = 0 . (6.1.1) 75 Numerical comparison models 76 6. NUMERICAL COMPARISON OF THE MODELS We use a Crank-Nicholson scheme and the relaxation method introduced by Besse-Bruneau in [10] and justified by Besse in [8] which replace the costly numerical treatment of the nonlinear term by a predictive step. This provides us with the following semi-discretized in time equation : u n+1 − u n dt + ∂x µ u n+1 + u n 2 ¶ + ε · 3 4 µ αun+ 1 2 ∂x µ u n+1 + u n 2 ¶ + (1 − α) u n+1 + u n 2 ∂xu n+ 1 2 ¶ + 1 6 ∂ 3 x u n+1 + u n 2 ¸ = 0 , where the predictive term u n+ 1 2 is defined as follows u n = u n+ 1 2 + u n−1/2 2 . The discretization of the nonlinear term u∂xu here takes advantage of the two possible discretizations u n+ 1 2 ∂x µ u n+1 + u n 2 ¶ and u n+1 + u n 2 ∂xu n+ 1 2 by introducing a parameter α ∈ [0, 1] and taking a convex combination of these possibilities. Keeping in mind that we want to preserve the semi-discrete L 2 norm, an easy integration by parts gives us the appropriate value α = 2/3. We then choose the spatial discretization so that the discrete L 2 norm is preserved by the complete scheme, which gives the final discretization of (6.1.1) : u n+1 i − u n i δt + µ D1 u n+1 + u n 2 ¶ i + ε   1 4  u n+ 1 2 i + u n+ 1 2 i+1 + u n+ 1 2 i−1 2   µ D1 u n+1 + u n 2 ¶ i + 1 4 u n+1 i + u n i 2 ³ D1u n+ 1 2 ´ i + 1 6 µ D3 u n+1 + u n 2 ¶ i ¸ = 0 , (6.1.2) where the matrix D1 and D3 are to the classical centered discretizations of the derivatives ∂x and ∂ 3 x . 6.1.2 Numerical scheme for the Boussinesq system As far as the discretization of the Boussinesq system (Σ) is concerned, we consider the same ideas. Using a Crank-Nicholson scheme and the same relaxation method, we aim here at preserving the specific norm |(v, η)| 2 H1 ε = |v| 2 L2 +|η| 2 L2 +εa2|∂xv| 2 L2 +εa4|∂η| 2 L2 . This quantity is indeed conserved by (Σ) (see [14] for more details). To this end, the nonlinear terms v∂xv, η∂xη, η∂xv and v∂xη are discretized in order to preserve both this specific discrete norm and their symmetric structure. Remarking that the equalities (v∂xv, v)L2 = 0 ; (η∂xη, v)L2 + (η∂xv, η)L2 + (v∂xη, η)L2 = 0 , 6.1. Numerical schemes 77 hold for (Σ) and using the same kind of method as for the KdV equation leads to the following semi-discretization of the nonlinear terms : . The matrix D1 and D3 are as defined in the KdV scheme, and the matrix D2 is the classical centered discretization of the derivative ∂ 2 x .

Initial data

Let us now talk about the initialization of the two schemes. First, all the prevision terms are initialized with a simple explicit integration of the equations on a half-step in time. Then, the initial conditions are chosen such that the simulated wave is unidirectional and propagating to the right. To this end, we first take the initial data of the second KdV equation to be zero. Then the system (ΣKdV ) reduces to the equation (6.1.1) for which we know the existence of solitary waves expressed as follows : u(t, x) = α cosh2 ¡ k(x − ct + l) ¢ , (6.1.4) with c = 1 + εα 4 , k = r 3α 8 and α, l being arbitrary. It is hence natural to specify the initial condition for the KdV equation (6.1.1) as follows : u(t = 0, x) = u0(x) = α cosh2 ¡ k(x + d) ¢ . (6.1.5) Finally, and because of the way the KdV approximation was constructed from the Boussinesq model, we specify the initial conditions for this latter as follows : v(t = 0, x) = η(t = 0, x) = 1 2 u0(x) . 6.2. Numerical results and comments 79 6.1.4 Validation of the numerical method With the initial data (6.1.5), the KdV scheme is expected to propagate the corresponding solitary wave to the right, without any deformation for all time. In order to validate this scheme, the numerical results obtained with the initial data (6.1.5) have been compared with the analytical solution (6.1.4). The following relative errors on the free surface have been computed in the L ∞ norm for several values of epsilon and for computation times T = 1/ε : ε T L δx δt relative error 0.05 20 80 0.03 0.03 1.5546.10−3 0.1 10 80 0.04 0.04 1.3717.10−3 0.2 5 80 0.05 0.05 1.0534.10−3 where L is the length of the computational domain and δx,δt are respectively the spatial and time discretization steps. These results allow to validate the scheme proposed for the KdV equations.

Numerical results and comments 

Numerical results

As specified in Chapter 2, the choice of the parameters a1, a2, a4 is very interesting in a numerical point of view. Indeed, the parameter a1 controls the presence of the dispersive terms ∂ 3 x v and ∂ 3 xη whereas the parameters a2 and a4 correspond to the terms ∂ 2 x∂tv and ∂ 2 x∂tη. These last terms have the main advantage of being regularizing terms analytically and numerically speaking, they smooth in some way the solution because they provide a control of the quantities ∂xv and ∂xη in the L 2 norm. We decided to use here the system (Σ) corresponding to a1 = a2 = a4 = 1/12 because it is likely to provide the better results. All the forthcoming results are expressed in non-dimensionalized variables. We recall that both the free surface and the bottom are of size ε : z = εη for the free surface and z = −1 + εb for the bottom. However, in order to get clear and readable results, we have plotted a rescaled free surface z = η and a rescale bottom z = −1 + b. A quick word on the duration T of the simulations : the previous chapters provided us with a justification of the models on large time scales of order O(1/ε). We have decided – only in the first example of the step – to overtake this large time scale and simulate the models on the very large time T = 1/ε3/2 , in order to see if the model remains stable on such time scales. The three models have been tested on two different examples of bottom. The first one 80 6. NUMERICAL COMPARISON OF THE MODELS correspond to a step at the bottom, defined similarly to [25] by b(x) =    0 , ∀x ∈ · 0, L 2 − 3 2 ¸ , β0 2 µ 1 + sin ³π 3 (x − L 2 ) ´¶ , ∀x ∈ · L 2 − 3 2 , L 2 + 3 2 ¸ , β0 , ∀x ∈ · L 2 + 3 2 , L¸ , (6.2.1) where β0 is an arbitrary constant of order O(1) and L is the length of the computation domain. The second example corresponds to a slowly varying sinusoidal bottom, defined as follows : b(x) = b0 sin µ π 2 + 2π l x ¶ , ∀x ∈ R , (6.2.2) where l is defined by l = 1 + εα/4 ε and α is the amplitude of the initial data defined in (6.1.5). The following results show the snapshots of the simulations at different times – so that the time evolution is relatively visible – and the evolution of the relative L ∞ error between the free surfaces obtained with the Boussinesq model and respectively the KdV approximation and the topographically modified approximation. The three models have been systematically plotted together in the same pictures in order to compare efficiently their respective behaviours. The numerical simulations have been performed for different values of ε in the case (6.2.1) of a step : ε = 0.05, ε = 0.1 and ε = 0.2, which are typical values of the upper part of the range of validity of the long waves approximation. As far as the case (6.2.2) is concerned, we simulated the models for the values ε = 0.05 and ε = 0.1. For all the simulations, the amplitude α of the initial free surface and the constant β0 linked to the bottom have been taken equal to 0.5. Here is a global tabular precising all the values of interest used in the simulations.