Numerical simulation of wave phenomena using FreeFEM

In this research free surface motion governed by the shallow water equations is considered. A numerical scheme based on the finite element method, which is incorporated in the open source FreeFEM, was used to simulate several wave phenomena. By carefully setting the corresponding initial condition as well as boundary conditions, several numerical computations were conducted. Numerical simulations presented here are standing wave in a closed basin, progressive wave over a flat bottom, as well as wave shoaling over a decreasing depth and wave refraction. In all cases above, the existing analytical formula are used to validate the numerical results. These computations suggest that explicit-implicit scheme is appropriate for varying water wave simulations.


Introduction
Water wave has been becoming an interesting topic to study. Many academicians try to figure out the behaviors of water wave in many phenomena. For example, the behavior of free surface wave investigated theoretically or by doing some experiments, see for instance [5], [6]. In connection with the development of computing tools extending very fast since few decades ago, some of academicians also try to find out the wave dynamic by using numerical simulation that commonly involves high computation.
FreeFEM is an open source software for solving various partial differential equations (PDEs) arising in many area of physics. The software adopts a numerical scheme based on the finite element method. In this article several 2D wave simulations such as standing wave, wave shoaling and wave refraction using FreeFEM were reported. Simulation results were validated using analytical formula.

Governing equations
The shallow water equations (SWEs), also known as Saint Venant's equations, is named after a French mathematician, Adhémar Jean Claude Barré de Saint-Venant (1797-1886). The equations are commonly used to describe free surface dynamic for incompressible fluid flow, like water in rivers, lakes, channels or coastal areas. In these applications horizontal length scale is much greater than the vertical length scale, therefore shallow water description is adequate. its amplitude is still low, and wavelength is large, as discussed in [1], [4], [7].
Let η denote the free surface elevation measured from the undisturbed water level, and let (u, v) T be the vector velocity of fluid particles in the x and y directions, and d 0 be the depth of fluid, see Figure 1. Next, consider the two-dimensional linear SWEs, which are In the equations above unknown variables are η, u, and v, which all are functions of spatial variables x, y and time variable t, while d 0 and gravitational acceleration g are constants. From physical point of view the first equation of (1) is derived from the principle of mass conservation, while the second is obtained from the momentum balance equation.

Numerical method
The finite element method (FEM) is one of numerical method which is based on grid. Generally, to implement FEM we start with making a grid by triangulating a domain into a number of triangles. The next step is rewrite the equation into its weak formulation. The weak formulation is obtained by multiplying the original equation (strong form) with a test function. In many cases weak form is usually yielded after an integration and implementation of Gauss's divergence theorem, in which we apply boundary conditions.

Weak formulation for two-dimensional linear SWEs
Consider SWEs (1) written in vector formulation, where u = (u, v) T and zero vector 0 ∈ R 2 . First, we apply the Euler time integration scheme to (2), followed with a multiplication with a test function ψ, and taking an integration over the domain Ω, a weak formulation of (2) reads Simplifying the second term on the left hand side of (4), after applying partial integration, one obtains By using the Gauss's divergence theorem applied into the third term on the left hand side of (5), results in where Γ = ∂Ω is a boundary and n is the outward pointing unit normal vector of the boundary.
Meanwhile, in similar fashion the weak form of (3) reads as where each term on the left hand side of (7) contains dot product operation with another test

Standing wave simulation
In this section, we implement FEM to solve SWEs as given in (1) to perform standing wave simulation with initial condition, Computational domain for this simulation is a square, L × L = 8 m × 8 m, with hardwall boundary conditions along with its four sides as shown in Figure 2. This hardwall boundary is formulated as homogeneous Neumann boundary condition as ∂η ∂n = ∇η · n = 0.  Weak formulation for SWEs that we will use are (6) and (7), which we call it explicit-implicit scheme. By expanding the third term on the left hand side of (6) yields Consider that Further, • for Γ 2 , since ∇η · n = η x = 0 and from (3) u t = −gη x , then In the same way, for Γ 3 and Γ 4 , respectively, we get v = 0 and u = 0. Therefore, Γ (dψu n h ) · n ds = 0.

Results on standing wave simulation
We conduct simulation for observational time 10 seconds. Here we set the water depth d 0 = 3 m, the gravitational acceleration g = 9.81 m/s 2 , mesh size h = 0.2667, and time step ∆t = 0.025 s. In FreeFEM we use the structured triangulation of a domain and finite element type in continuous piecewise linear P 1 -element. Then, FreeFEM will generate the data recording the (x, y) particle position along with free surface at its position. We will replot the generated data using gnuplot. Finally, simulation results at subsequent times are depicted in

Error analysis
For standing wave simulation, the exact solutions of equations (1) with initial condition, Error analysis for this computation was conducted by calculating the relative error r , which is an L 2 norm for a two-dimensional domain area as follows To analyze error resulted in standing wave simulation, we also conducted the same simulation, but using other time integration schemes. The results are wave experiences damping for implicitimplicit scheme, blows up for explicit-explicit scheme, and does not undergo any damping for implicit-explicit scheme. After doing calculations using r formula above, we get the relative error in L 2 norm from each scheme as tabulated in Table 1. From Table 1, it is noted that the combination of explicit Euler and implicit Euler scheme is stable from all schemes that we conducted in these standing wave simulations.

Wave shoaling simulation
In this section we will conduct a simulation of wave shoaling. When a wave propagates from deep water towards a shoreline over a shoal of decreasing depth, it experiences shoaling phenomena; because of a decreasing depth, wave speed will be decreasing, amplitude increasing and wavelength decreasing. For shoaling simulation we adopt the SWEs as given in (1) subjecting to zero initial condition, η(x, y, 0) = 0, u(x, y, 0) = 0, v(x, y, 0) = 0.  7 We will conduct wave shoaling simulation by setting a rectangular L × M = 6 m × 2 m as a computational domain with left wave influx boundary condition, η(0, y, t) = 0.25 sin(5πt), (8) right absorbing boundary condition, η(L, y, t) = 1 and homogeneous Neumann boundary condition on the rest of boundaries. Note that this pair of boundary conditions, (8) and (9), are important for having left running wave influx, see [8] for detail explanation. Weak formulations of SWEs model that we use are (6) and (7) and for boundary conditions we use schemes as follows where Γ 2 and Γ 4 are corresponding to boundary on (x = L, y) and (x = 0, y), respectively. In simulation we set the mesh size h = 0.066667, time step ∆t = 0.004038 s, the gravitational acceleration g = 9.81 m/s, and bathymetry is formulated as follows where d 1 = −5 m and d 2 = −1 m. To make sure that implementation of right absorbing boundary works well we conducted simulation till time T = 1.6 s, that is when the wave has reached the boundary x = L. In FreeFEM we use the structured triangulation of a domain and finite element type in continuous piecewise linear P 1 -element. The result on wave shoaling simulation is depicted in Figure 5, which is the snapshot of free surface at time t = 1.59905 s. It clearly shows a shoaling phenomena; wavelength decreasing, amplitude increasing, as the effect of shallower depth. In the next subsection, the amplification of wave amplitude will be examined and be compared with the analytical formula.

Validation of results on wave shoaling simulation
In this section, we conduct wave shoaling simulation. The amplification of wave in shoaling phenomena should follow the Wentzel-Kramer-Brillouin (WKB) formula, This formula is obtained from energy conservation during wave propagation, detail discussion about it can be obtained in many standard textbook about wave, see for instance [2], [3].
We carried out several times of simulation with various of water depths. From simulation result, we recorded wave amplitude at x ∈ [4.5, 6.0] in y = 0 at a certain time. We compare this numerical results of A 2num with exact amplitude A 2exact . Here, we use A 1 = 0.25 as well as the exact amplitude corresponding to water depth d 1 and find the absolute error at y = 0 as formulated in The comparison between numerical and analytical wave amplitude along with the value of other variables is tabulated in Table 2. From that Table we conclude that the error resulted in wave shoaling simulation are all below 1.5%. Therefore, we can conclude that our FEM executed using FreeFEM can simulate wave shoaling accurately. Table 2. Comparison between numerical and analytical wave amplitude in wave shoaling simulation at time t where d 1 , d 2 are water depth with respect to topography, ∆t is time step and A 2 is wave amplitude in y = 0.

Wave refraction simulation
In this section we will perform another wave phenomenon, namely wave refraction. This is a bending wave phenomenon commonly observed near the coast. This phenomenon happens mainly because the phase speed of the wave is proportional to the square root of water depth. Therefore, when a wave approaches a shore, due to refraction mechanism its wave front bends so that it becomes parallel to the coastline.
For conducting wave refraction simulation, we set the initial and boundary condition that is the same with wave shoaling simulation, but with left wave influx boundary condition, η(0, y, t) = 0.15 sin(5πt). The result on simulation is depicted in Figure 6 showing a snapshot of free surface at t = 1.99525 s. To observe refraction phenomena we plot free surface contour resulting on wave refraction simulation at that time as shown in Figure 7. In that figure it is shown that wave bends when it approaches shallower depth.

Conclusions
In this article, we have conducted several wave simulations by utilizing FreeFEM that solves the shallow water equations using the finite element method. Generally, combination of explicit and implicit schemes in time integration would give us simulation results that are free from numerical damping error. Several simulations of benchmark tests have been conducted, such as, standing wave simulation, wave shoaling phenomena as well as wave refraction phenomena. Comparison with analytical formula the numerical results show good agreement. It should be noted that in this study we focus on implementing FreeFEM software for simulating linear wave phenomena in shallow water areas. For simulations that are directly related to the application of real problems, an in-depth and comprehensive study is still needed, specifically the model must take into account non-linear effects.