Modelling and simulation of a wave energy converter

In this work we present the mathematical model and simulations of a particular wave energy converter, the so-called oscillating water column. In this device, waves governed by the one-dimensional nonlinear shallow water equations arrive from offshore, encounter a step in the bottom and then arrive into a chamber to change the volume of the air to activate the turbine. The system is reformulated as two transmission problems: one is related to the wave motion over the stepped topography and the other one is related to the wave-structure interaction at the entrance of the chamber. We finally use the characteristic equations of Riemann invariants to obtain the discretized transmission conditions and we implement the Lax-Friedrichs scheme to get numerical solutions.


General setting
This work is devoted to model and simulate an on-shore oscillating water column (OWC), which is a particular type of wave energy converter (WEC) that transforms the energy of waves reaching the shore into electric energy. The structure is installed at the shore in such a way that the water partially fulfills a chamber, which is connected with the outside through a hole where a turbine is placed (see Figure 1). Incoming waves collide with the exterior part of the immersed wall and, after the collision, one part of the wave is reflected while the other part passes below the fixed partially immersed wall and enters the chamber. This increases the water volume inside the chamber and consequently, it creates an airflow that actives the turbine by passing through it and the same occurs when the volume of water reduces inside the chamber. The perpetuation of the incoming waves makes the water inside the chamber oscillate and act as a liquid piston, whose oscillations create electric energy. In this work the wave energy converter is deployed with stepped bottom, which means that incoming waves encounter a step in the bottom topography just before reaching the structure. The influence of such step in the OWC device will be discussed later in Section 4. The present research is essentially motivated by a series of works by Rezanejad and collaborators on the experimental and numerical study of nearshore OWCs, in particular, we refer to Rezanejad and Soares [14], where the authors used a linear potential theory to do simulations and showed the improvement of the efficiency when a step is added. Our goal is to numerically study this type of WEC considering as the governing equations for this wave-structure interaction the nonlinear shallow water equations derived by Lannes in [8], whose local well-posedness was obtained by Iguchi  [7] in the one-dimensional case and by Bocchi in [1] in the two-dimensional axisymmetric case. In the Boussinesq regime and for a fixed partially immersed solid similar equations were studied by Bresch, Lannes and Métivier in [2] and in the shallow water viscous case by Maity, San Martín, Takahashi and Tucsnak in [11] and by Vergara-Hermosilla, Matignon, and Tucsnak in [15]. We consider an incompressible, irrotational, inviscid and homogeneous fluid in a shallow water regime, which occurs in the region where the OWC is installed. Following [8], the motion of the fluid is governed by the 1D nonlinear shallow water equations where ζ(t, x) is free surface elevation, h(t, x) is the fluid height, ρ is the fluid density, P is the surface pressure of the fluid and q(t, x) is the horizontal discharge defined by where u(t, x, z) is the horizontal component of the fluid velocity vector field. Let us first give the boundary conditions related to (1). The relevance of these boundary conditions will be explained in Section 2.2. The boundary conditions on the horizontal discharge are q is continuous at x = 0, x = l 0 ± r, and the boundary conditions on the surface elevation are where f is a prescribed function depending only on time. The surface pressure is given by the constant atmospheric pressure where the fluid is directly in contact with the air, i.e.
and no surface tension is considered here. On the other hand, under the partially immersed structure, the fluid surface elevation is constrained to be equal to the parametrization of the bottom of the solid ζ w , i.e. ζ = ζ w in (l 0 − r, l 0 + r).
To complete the system, we consider an initial configuration where the fluid is at rest,

Organization of the paper
In Section 2, we derive the model used in the numerical simulations following [1,2,7]. In particular, we show that the equations (1) can be reformulated as two transmission problems, one related to the step in the bottom topography and one related to the wave-structure interaction at the entrance of the chamber. Furthermore, the equations in the exterior domain are written as two transport equations on Riemann invariants. In Section 3, we discretize the equations in conservative form using the Lax-Friedrichs scheme and use the Riemann invariants to derive the discretization of the entry condition and boundary conditions. In Section 4, we give several computations showing the numerical solutions of the model and compare the OWC device with and without stepped bottom. At the end of this section, we show the accuracy of the numerical scheme to validate our computations and we discuss the absorbed power and the efficiency of the OWC.

Notations
We divide the domain of the problem (−l, l 1 ) into two parts. The interval I = (l 0 − r, l 0 + r) is called interior domain, which is the projection onto the line of the wetted part of the structure, and its complement E = (−l, l 1 ) \ I, called exterior domain, which is the union of three intervals where l 1 is the position of the end of the chamber and l 0 and r are respectively the position of the center and the half length of the partially immersed structure. From the nature of the problem, l 1 > l 0 > r. Moreover, the boundary of I is formed by the contact points {l 0 ± r}, which are the projections on the real line of the triple contact points between fluid, solid and air. For any function f defined in the real line, its restrictions on the interior domain and the exterior domain are respectively denoted by

Governing equations
In this section, we present the mathematical model that describes the oscillating water column process considered in this work. The model can be essentially divided in three parts: the wave motion over a discontinuous topography represented by the step, the wave-structure interaction at the entrance of the chamber and the wave motion in the chamber. In the exterior domain E, where the fluid is in contact with the air, the surface pressure P e is constrained and is assumed to be equal to the constant atmospheric pressure P atm , while the surface elevation ζ e is not known. Contrarily, in the interior domain I, that is the region under the partially immersed structure, the surface elevation ζ i is constrained to coincide with the parametrization of the wetted surface, which is assumed to be the graph of some function ζ w . The surface pressure P i is unknown and it turns out to be a Lagrange multiplier associated with the constraint on ζ i . For more details on this approach for the study of wave-structure interaction, we refer to [8]. In this work we consider a partially immersed fixed structure with vertical side walls, the parametrization ζ w is a constant both in time and space. Summing up, we have an opposite behaviour for the surface elevation and the surface pressure under the structure and elsewhere, that is ζ i = ζ w , P i is unknown and ζ e is unknown, P e = P atm .
For the exterior domain, we distinguish the region before the step, denoted by E 0 and the region after the step, denoted by E 1 ∪ E 2 . The fluid heights are defined respectively by where h s and h 0 are the fluid heights at rest before the step and after the step respectively. Denoting by s the height of the step, we have h s = h 0 + s. Therefore the nonlinear shallow water equations (1) can be written as the following three systems: and h e = h s + ζ e ,

Derivation of the transmission conditions
The following section is devoted to showing that the motion over the stepped bottom and the wavestructure interaction can be reduced to two transmission problems for the nonlinear shallow water equations. To do that, we derive the transmission conditions relating the different parts of the model, respectively at the step in the bottom topography and at the side walls of the partially immersed structure.

At the topography step
We consider the problem before the entrance of the chamber not as one shallow water system with a discontinuous topography but rather as a transmission problem between two shallow water systems with flat bottoms where the fluid heights are respectively h s + ζ e and h 0 + ζ e . The first transmission condition is given by the continuity of the surface elevation at the step, namely where the traces at x = 0 − and at x = 0 + are the traces at x = 0 of the unknowns before the step and after the step respectively. The second transmission condition is given by the continuity of the horizontal discharge at the step, namely

At the structure side-walls
The transmission conditions at the side-walls of the partially immersed structure are derived from the continuity of the horizontal discharge at the side-walls and the assumption that the total fluid-structure energy is equal to the integral in time of the energy flux at the entry of the domain. The continuity of the horizontal discharge at x = l 0 ± r together with the fact that ∂ x q i = 0 gives the first transmission condition between the E 1 and E 2 , which reads Let us now derive the second transmission condition at x = l 0 ± r. To do that, we show the local conservation of the fluid energy in the exterior domain and in the interior domain as in [2].
Exterior domain. Considering the nonlinear shallow water equations in E, multiplying the first equation in (7)- (8) by ρgζ e and the second equation by ρ q e h e , and considering the fact that Adding both equations in (13), we obtain We compute that and, denoting by e ext and by f ext respectively the local fluid energy and the local flux we obtain the local conservation of the fluid energy in the exterior domain Interior domain. Let us remark that from the first equation in (9) one gets that q i ≡ q i (t) in the interior domain. Multiplying the second equation in (9) by and, denoting by e int and by f int respectively the local fluid energy and the local flux we obtain the local conservation of the fluid energy in the interior domain, Now we assume that the total fluid-structure energy at time t is equal to the integral between 0 and t of the sum between the energy flux at the entry of the domain and the difference of the energy fluxes at the step, i.e.
with the fluid energy defined by This assumption is an adaptation to a bounded domain case of the conservation of total fluid-structure energy assumed in [2]. We remark that the difference of the energy fluxes at the step f ext| x=0 + − f ext| x=0 − does not vanish due to the discontinuity of the fluid height at x = 0 in the presence of the step. The fact that the structure is fixed ( d dt E solid = 0) yields From (14) and (15) we have Using the boundary conditions (2) and (3) we get where the brackets · are defined as in (12). By definition of the fluxes it follows + gζ e q e and from (2) and (12) we obtain Integrating on (l 0 + r, l 0 − r), the second equation in (9) yields Combining the last two equalities, we get the following transmission condition

Reformulation as two transmission problems
Coupling the governing equations (7)-(9) with the conditions derived in the previous section, we have therefore reduced the problem of the OWC essentially to two transmission problems. The first one in E 0 ∪ E 1 reads: with transmission conditions at x = 0 The second transmission problem in E 1 ∪ E 2 reads: with transmission conditions at x = l 0 ± r where α = 2r h w and h w = h 0 + ζ w .

Riemann invariants
Let us now rewrite the nonlinear shallow water equations (7) and (8) in the exterior domain E in a compact form by introducing the couple U = (ζ e , q e ) T : where Notice that λ + > 0 and λ − > 0. Taking the scalar product of (21) and eigenvectors, we obtain ∂ t 2 gh e ± q e h e ± gh e ± q e h e ∂ x 2 gh e ± q e h e = 0.
Let us introduce the right and the left Riemann invariant R and L associated to the nonlinear shallow water equations, respectively Hence we can write the 1D nonlinear shallow water equations in the exterior domain as the two following transport equations on R and L: We will see that these two transport equations of Riemann invariants are helpful when we solve our model by numerical method. More details about Riemann invariants of the nonlinear shallow water equations can be found in [9].

Discretization of the model
We have reformulated in the previous section the mathematical model of the oscillating water column as two transmission problems. This section is devoted to discretize the nonlinear shallow water equations (7)-(9) at the level of the numerical scheme. More precisely, we will use the Lax-Friedrichs scheme to solve our main equations and use Riemann invariants to address the entry conditions and all boundary conditions.

Discretization of the equation
The finite difference method is a standard discretization approach for partial differential equations, especially those that arise from conservation laws. We first rewrite equation (21) as the following conservative form : with By means of a finite difference approach, equation (24) can be discretized as where the flux F is discretized with cell centres indexed as i and cell edge fluxes indexed as i ± 1/2. The choice of F m i±1/2 depends on the numerical scheme. We consider here the well-known Lax-Friedrichs scheme proposed by Lax [10] to get the discrete flux where i ≥ 1 and F m i = F (U m i ).

Discretization of the entry condition
At the entrance of our system, the surface elevation is given by a prescribed function f depending only on time, In order to express the entry condition for the horizontal discharge, let us first recall that from (22) one has where R and L are respectively the right and the left Riemann invariant associated to the nonlinear shallow water equations. We get Hence, the value of q e at x = −l is given by On the right-hand side of the relation above, L| x=−l is unknown. First we have to determine L| x=−l in order to determine q e | x=−l . This can be achieved by the transport equation for L in (23). After discretizing it as in [12], we get where L m 0 is the value of L at x = −l at time t m and λ − is computed as a linear interpolation between λ −,0 and λ −,1 following [9], namely with 0 ≤ β ≤ 1 such that λ − δ t = βδ x . Moreover, we can compute λ − as .

Thus, we have
which gives L m 0 in terms of its values at the previous time step and in terms of interior points.

Discretization of the boundary conditions
Since our system is composed by four parts, it remains three boundary conditions should be taken into consideration besides the entry condition at x = −l. When wave arrives from the offshore, it will encounter a step in the bottom and then arrive into a chamber, and finally arrive to the wall (see the configuration 1). More precisely, the first boundary condition is at the discontinuity of the topography located at x = 0 and the second is at the partially immersed structure side-walls located at x = l 0 ± r. The last boundary condition is at the end of the chamber, located at x = l 1 .

At the topography step
Let us first consider the shallow water wave equations with discontinuous topography, namely, it is a system with depth h s on R − = {x < 0} and depth h 0 on R + = {x > 0}. Our equation turns out to be From transmission conditions (10) and (11), we have the continuity of the surface elevation ζ e and of the horizontal discharge q e at x = 0: Let us denote the right Riemann invariant in the domain R − by R l and the left Riemann invariant in the domain R + by L r . We then find two expressions of q e describing q l e | x=0 and q r e | x=0 , respectively, According to the relations (28), we observe that (29) is a system of two nonlinear equations on the two unknowns ζ l e | x=0 (respectively ζ r e | x=0 ) and q l e | x=0 (respectively q r e | x=0 ). We write it in the compact form where x 1 = ζ l e | x=0 , x 2 = q l e | x=0 and the vector F = (F 1 , F 2 ) is given by In the case h s = h 0 (without step) we can derive from (29) a third degree equation on h 0 + ζ l e | x=0 and take the unique solution that gives ζ l e | x=0 = 0 when R l | x=0 , L r | x=0 = 0 (we refer to [8] for this case). Here, since h s = h 0 , we use MATLAB nonlinear system solver fsolve with initial point (0, 0) to solve (30). Before doing that, we have to determine the values of the two Riemann invariants R l | x=0 and L r | x=0 . The transport equations for R l and L r are the following: where the corresponding eigenvalue λ l + in the domain R − is given by and the corresponding eigenvalue −λ r − in the domain R + is given by Let us emphasize that we use here the same interpolation for λ + and λ − as in [12]. After discretization of equations (31), we get where λ l + , λ r − are as in (32)-(33) and we recall that (R l ) m n1,x is the value of R l at x n1,x and t m (see Notations 3.0.1). Hence, we have Gathering the relations (28), (29) and (34), we can solve ζ l e | x=0 (respectively ζ r e | x=0 ) and q l e | x=0 (respectively q r e | x=0 ), which give us the boundary conditions at the step.

At the structure side-walls
Compared with the derivation of the boundary conditions near the step, the idea to derive the boundary condition near the fixed partially immersed structure is almost the same. There are two differences between them. The first one is that, since the depth is always h 0 , the eq. (29) becomes where we denote the horizontal discharge in the exterior domain on the left-hand side of the object by q l e and on the right-hand side of the object by q r e . Let us recall that q i is the horizontal discharge in the interior domain I. From the first transmission condition in (20), we know that The second difference is that, unlike in the previous subsection, we do not have the continuity condition of ζ e at the structure side-walls. Nevertheless, we consider the discretization of the second transmission condition in (20), hence we get where for the sake of clarity (q e ) m l0−r = (q e ) m n1,x+n2,x and (q e ) m l0+r = (q e ) m n1,x+n2,x+n3,x (analogously for (ζ e ) m l0−r and (ζ e ) m l0+r ). Then, q e at x = l 0 − r is expressed as which gives (q e ) m l0−r in terms of its values at the previous time step and in terms of interior points. Now we can solve (q e ) m l0−r immediately. Once the value of (q e ) m l0−r is obtained, we can find the values of ζ l e | x=l0−r and ζ r e | x=l0+r by using equations (35) and the transport equations for the Riemann invariants as the strategy in Section 3.3.1.

At the end of the chamber
The corresponding boundary condition at the end of the chamber, located at x = l 1 , is given by Hence, recalling the definition of the right-going Riemann invariant R, we recover the surface elevation ζ e at x = l 1 , namely

Numerical validations
In this section, we use the scheme introduced in Section 3 to simulate our model. For the fluid, we always consider the density of water ρ = 1000 kg/m 3 and the gravitational acceleration g = 9.81 m/s 2 . The entry of the domain is set at x = −l = −30 m and the prescribed function f is given by where T = 1.5 s is the period. Using the notations as before, we consider l 0 = 11 m, r = 1 m and l 1 = 17 m and the fluid height at rest before the step h s = 15 m. We compute the solution by using the Lax-Friedrichs scheme in the exterior domain [−30, 10] ∪ [12,17], with a refined mesh with N x = 2300 and a time step δ t = 0.7 √ ghs δ x with space step δ x = 0.02 m. Here, the CFL number is 0.7, which is commonly used to prescribe the terms of the finite-difference approximation of a PDE (see for instance [13]). In the interior domain, the solution can be computed using the transmission conditions (20) with h w = h 0 + ζ w and ζ w = −7.5 m.

Numerical solutions
In real applications, an OWC device can be deployed on a stepped sea bottom in order to improve its performance. It is important then to have a good understanding of the impact of a step in the topography. We find that, before the waves encounter the step, there is no significant difference between the OWC model without stepped bottom and with stepped bottom (see (a) and (b)). But when the waves encounter the step in the bottom and arrive into the chamber, we can see that, the waves in the OWC model without stepped bottom move significantly faster than the waves in the OWC model with stepped bottom. In particular, at t = 3.3 s the waves in the OWC model without stepped bottom has already arrived to the chamber and will begin to change the water level in the chamber, while the waves in the OWC model with stepped bottom have not reached yet and the water will rise inside the chamber later (see (c) and (d)). As the step at bottom is a sort of obstacle for the incoming wave, this phenomenon is reasonable.
As one may expect, the incoming wave split into two parts when it touches the left wall of the partially immmersed structure. One part enters the chamber and changes the volume of the air that makes the turbine rotate. The other part is reflected and becomes an outgoing wave, as we can see in Figure 2. At t = 5 s, the reflected wave in the OWC model without stepped bottom already reaches x = −10 m, while the reflected wave in the OWC model with stepped bottom has not reached x = −10 m (see (e) and (f)). This shows that the reflected waves in the OWC model with stepped bottom move slower than the waves in the OWC model without stepped bottom.
This difference can be explained by the fact that more incident wave energy is converted when a step is added. In other words, the OWC with stepped bottom would be more efficient than the one without stepped bottom, which is in agreement with the result by Rezanejad and Soares in [14].

Accuracy analysis
In numerical validations, accuracy analysis is of importance. As we can see in Figure 1, the configuration of OWC device is essentially constituted from three parts: the domain before the step in the sea bottom, the domain after the step and the chamber. We implement our algorithms by gathering together the three parts. It is worth mentioning that one compact algorithm is also actionable.
In order to make it possible to verify our algorithm, we do the following accuracy analysis. Under the same initial wave and physical parameters, we compare the free surface elevation ζ e of the classical nonlinear shallow water wave model with our model without discontinuous topography. Figure 3 shows that there is no significant difference between the two cases. Moreover, we also find that the error is of order 10 −3 (see Figure 4), which is acceptable since the Lax-Friedrichs method is first-order accurate in space.

Absorbed power and efficiency
Designing a WEC of high efficiency is nowadays a hot topic in all regions and countries over the world. In this regard, we present in this section the method to calculate the absorbed power as well as the efficiency of the OWC considered in this work. The primary efficiency η Reg of the device is defined by the ratio of the absorbed power from the waves to the incident wave power. From the seminal work of Evans in [6], we know that in the linear timeharmonic theory the volume flux Q(t) = Re{qe −iωt } is assumed linearly proportional to the pressure in the chamber P (t) = Re{pe −iωt }. Using this assumption, the average power absorbed from regular waves over one wave period, denoted by P Reg , is given by where p is the time independent and λ is a positive constant associated with linear air turbine characteristics. On the other hand, following [14] in experiments the average power absorbed from regular waves can be determined by: where T is the duration of the test. The incident wave power P inc is defined as the product of total energy per wave period E inc and the group velocity c g (see [3]): 1 + 2kh s sinh(2kh s ) , and the dispersion relation given by where ω is the frequency, k is the wave number, h s is the fluid height at rest before the step, ρ is the density of the fluid, g is the gravitational acceleration and L is the projected width of the WEC perpendicular to the incident wave direction, A is the amplitude of the wave. In the shallow water regime kh s 1 and the group velocity reduces to c g = √ gh s . Thus, the primary efficiency of the device in regular wave is given by We notice that in both (37) and (38) the absorbed power (hence the primary efficiency) strongly depends on the air pressure in the chamber. In our model, it is considered to be a constant, namely the atmospheric pressure P atm . However, when the waves arrive into the chamber and change the volume of the air, the air pressure in the chamber will certainly change as well. In this case, the pressure will no more be a constant, but depends on time. Hence, to study more rigorously the absorbed power and the primary efficiency of the OWC, this fact must be taken into account in the model. This will be addressed in our future work. Analogously, the improvement of the efficiency of an OWC device deployed on a stepped sea bottom can be also investigated with a better knowledge of the air pressure in the chamber. From the results in Section 4.1, we can expect that significant improvements in the efficiency can be achieved by adding a step at the bottom of the sea.