Modelling and simulation of a wave energy converter
EESAIM: PROCEEDINGS AND SURVEYS,
Vol. ?, 2021, 1-16
Editors: Will be set by the publisher
MODELLING AND SIMULATION OF A WAVE ENERGY CONVERTER
Edoardo Bocchi , Jiao He and Gast´on Vergara-Hermosilla Abstract.
In this work we present the mathematical model and simulations of a particularwave energy converter, the so-called oscillating water column. In this device, waves governed bythe one-dimensional nonlinear shallow water equations arrive from offshore, encounter a step inthe 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 motionover the stepped topography and the other one is related to the wave-structure interaction atthe entrance of the chamber. We finally use the characteristic equations of Riemann invariantsto obtain the discretized transmission conditions and we implement the Lax-Friedrichs schemeto get numerical solutions. Introduction
General setting
This work is devoted to model and simulate an on-shore oscillating water column (OWC), which isa particular type of wave energy converter (WEC) that transforms the energy of waves reaching theshore into electric energy. The structure is installed at the shore in such a way that the water partiallyfulfills a chamber, which is connected with the outside through a hole where a turbine is placed (seeFigure 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 andenters the chamber. This increases the water volume inside the chamber and consequently, it creates anairflow that actives the turbine by passing through it and the same occurs when the volume of waterreduces inside the chamber. The perpetuation of the incoming waves makes the water inside the chamberoscillate and act as a liquid piston, whose oscillations create electric energy. In this work the wave energyconverter is deployed with stepped bottom, which means that incoming waves encounter a step in thebottom topography just before reaching the structure. The influence of such step in the OWC devicewill be discussed later in Section 4. The present research is essentially motivated by a series of works byRezanejad 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 simulationsand showed the improvement of the efficiency when a step is added. Our goal is to numerically study thistype of WEC considering as the governing equations for this wave-structure interaction the nonlinearshallow water equations derived by Lannes in [8], whose local well-posedness was obtained by Iguchi Departamento de An´alisis Matem´atico & Instituto de Matem´aticas de la Universidad de Sevilla (IMUS), Universidadde Sevilla, Avenida Reina Mercedes, 41012 Sevilla, Espa˜na (e-mail: [email protected]) Laboratoire de Math´ematiques et Mod´elisation d’Evry (LaMME), Universit´e d’Evry Val d’Essonne, 23 Boulevard deFrance, 91037, Evry Cedex, France (e-mail: [email protected]) Institut de Math´ematiques de Bordeaux (IMB), Universit´e de Bordeaux, 351 Cours de la Lib´eration, 33405 TalenceCedex, France (e-mail: [email protected]) © EDP Sciences, SMAI 2021 a r X i v : . [ m a t h . A P ] F e b ESAIM: PROCEEDINGS AND SURVEYS h s h ζ ( x, t ) z h w Incoming Waves Air s − l l − r l l + r l x Figure 1.
Configuration of the OWCand Lannes in [7] in the one-dimensional case and by Bocchi in [1] in the two-dimensional axisymmetriccase. In the Boussinesq regime and for a fixed partially immersed solid similar equations were studied byBresch, Lannes and M´etivier in [2] and in the shallow water viscous case by Maity, San Mart´ın, Takahashiand 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 governedby the 1D nonlinear shallow water equations ∂ t ζ + ∂ x q = 0 ∂ t q + ∂ x (cid:18) q h (cid:19) + gh∂ x ζ = − hρ ∂ x P for x ∈ ( − l, l ) , (1)where ζ ( t, x ) is free surface elevation, h ( t, x ) is the fluid height, ρ is the fluid density, P is the surfacepressure of the fluid and q ( t, x ) is the horizontal discharge defined by q ( t, x ) := (cid:90) ζ ( t,x ) − h u ( t, x, z ) dz, 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 willbe explained in Section 2.2. The boundary conditions on the horizontal discharge are q is continuous at x = 0 , x = l ± r,q = 0 at x = l , (2)and the boundary conditions on the surface elevation are ζ = f at x = − l,ζ is continuous at x = 0 , (3)where f is a prescribed function depending only on time. The surface pressure is given by the constantatmospheric pressure where the fluid is directly in contact with the air, i.e. P = P atm in ( − l, l − r ) ∪ ( l + r, l ) (4)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 SAIM: PROCEEDINGS AND SURVEYS ζ w , i.e. ζ = ζ w in ( l − r, l + r ) . (5)To complete the system, we consider an initial configuration where the fluid is at rest, ζ (0 , x ) = (cid:40) − l, l − r ) ∪ ( l + r, l ) ζ w in ( l − r, l + r ) and q (0 , x ) = 0 . (6) 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 thestep in the bottom topography and one related to the wave-structure interaction at the entrance of thechamber. Furthermore, the equations in the exterior domain are written as two transport equationson 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 andboundary conditions. In Section 4, we give several computations showing the numerical solutions of themodel and compare the OWC device with and without stepped bottom. At the end of this section, weshow the accuracy of the numerical scheme to validate our computations and we discuss the absorbedpower and the efficiency of the OWC.1.2.1.
Notations
We divide the domain of the problem ( − l, l ) into two parts. The interval I = ( l − r, l + r ) iscalled interior domain , which is the projection onto the line of the wetted part of the structure, and itscomplement E = ( − l, l ) \ I , called exterior domain , which is the union of three intervals E ∪ E ∪ E with E = ( − l, , E = (0 , l − r ) and E = ( l + r, l ) , where l is the position of the end of the chamber and l and r are respectively the position of the centerand the half length of the partially immersed structure. From the nature of the problem, l > l > r .Moreover, the boundary of I is formed by the contact points { l ± r } , which are the projections on thereal line of the triple contact points between fluid, solid and air. For any function f defined in the realline, its restrictions on the interior domain and the exterior domain are respectively denoted by f i := f | I and f e := f | E . Presentation of the model
Governing equations
In this section, we present the mathematical model that describes the oscillating water column processconsidered in this work. The model can be essentially divided in three parts: the wave motion over adiscontinuous topography represented by the step, the wave-structure interaction at the entrance of thechamber and the wave motion in the chamber. In the exterior domain E , where the fluid is in contact withthe air, the surface pressure P e is constrained and is assumed to be equal to the constant atmosphericpressure P atm , while the surface elevation ζ e is not known. Contrarily, in the interior domain I , that isthe region under the partially immersed structure, the surface elevation ζ i is constrained to coincide withthe parametrization of the wetted surface, which is assumed to be the graph of some function ζ w . Thesurface pressure P i is unknown and it turns out to be a Lagrange multiplier associated with the constrainton ζ i . For more details on this approach for the study of wave-structure interaction, we refer to [8]. Inthis 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 ESAIM: PROCEEDINGS AND SURVEYS 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 and the region afterthe step, denoted by E ∪ E . The fluid heights are defined respectively by h e = h s + ζ e in E , h e = h + ζ e in E ∪ E , where h s and h 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 + s .Therefore the nonlinear shallow water equations (1) can be written as the following three systems:(1) for x ∈ E , ∂ t ζ e + ∂ x q e = 0 ,∂ t q e + ∂ x (cid:18) q e h e (cid:19) + gh e ∂ x ζ e = 0 and h e = h s + ζ e , (7)(2) for x ∈ E ∪ E ∂ t ζ e + ∂ x q e = 0 ,∂ t q e + ∂ x (cid:18) q e h e (cid:19) + gh e ∂ x ζ e = 0 and h e = h + ζ e , (8)(3) for x ∈ I , ∂ x q i = 0 ,∂ t q i = − h w ρ ∂ x P i and h w = h + ζ w . (9) Derivation of the transmission conditions
The following section is devoted to showing that the motion over the stepped bottom and the wave-structure interaction can be reduced to two transmission problems for the nonlinear shallow water equa-tions. To do that, we derive the transmission conditions relating the different parts of the model, respec-tively at the step in the bottom topography and at the side walls of the partially immersed structure.2.2.1.
At the topography step
We consider the problem before the entrance of the chamber not as one shallow water system with adiscontinuous topography but rather as a transmission problem between two shallow water systems withflat bottoms where the fluid heights are respectively h s + ζ e and h + ζ e .The first transmission condition is given by the continuity of the surface elevation at the step, namely ζ e | x =0 − = ζ e | x =0+ , (10)where the traces at x = 0 − and at x = 0 + are the traces at x = 0 of the unknowns before the step andafter the step respectively.The second transmission condition is given by the continuity of the horizontal discharge at the step,namely q e | x =0 − = q e | x =0+ . (11) SAIM: PROCEEDINGS AND SURVEYS
At the structure side-walls
The transmission conditions at the side-walls of the partially immersed structure are derived from thecontinuity of the horizontal discharge at the side-walls and the assumption that the total fluid-structureenergy is equal to the integral in time of the energy flux at the entry of the domain. The continuity ofthe horizontal discharge at x = l ± r together with the fact that ∂ x q i = 0 gives the first transmissioncondition between the E and E , which reads (cid:74) q e (cid:75) := q e | x = l r − q e | x = l − r = 0 . (12)Let us now derive the second transmission condition at x = l ± r . To do that, we show the local conser-vation 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 equa-tion in (7)-(8) by ρgζ e and the second equation by ρ q e h e , and considering the fact that ∂ t h e = − ∂ x q e , weobtain ∂ t (cid:18) ρ q e h e (cid:19) + ρgζ e ∂ x q e = 0 ,∂ t (cid:18) ρ q e h e (cid:19) − ρ q e h e ∂ x q e + ρ q e h e ∂ x (cid:18) q e h e (cid:19) + ρgq e ∂ x ζ e = 0 . (13)Adding both equations in (13), we obtain ∂ t (cid:18) ρg ζ e ρ q e h e (cid:19) + ρgζ e ∂ x q e + ρgq e ∂ x ζ e − ρ q e h e ∂ x q e + ρ q e h e ∂ x (cid:18) q e h e (cid:19) = 0 . We compute that gζ∂ x q + gq∂ x ζ = ∂ x ( gζq ) and − q h ∂ x q + qh ∂ x (cid:18) q h (cid:19) = ∂ x (cid:18) q h (cid:19) , and, denoting by e ext and by f ext respectively the local fluid energy and the local flux e ext = ρ q e h e + gρ ζ e f ext = ρ q e h e + gρζ e q e , we obtain the local conservation of the fluid energy in the exterior domain ∂ t e ext + ∂ x f ext = 0 . (14) Interior domain . Let us remark that from the first equation in (9) one gets that q i ≡ q i ( t ) in the interiordomain. Multiplying the second equation in (9) by q i h w , we obtain ∂ t (cid:18) ρ q i h w (cid:19) + ∂ x ( q i P i ) = 0 , and, denoting by e int and by f int respectively the local fluid energy and the local flux e int = ρ q i h w + ρg ζ w f int = q i P i , we obtain the local conservation of the fluid energy in the interior domain, ∂ t e int + ∂ x f int = 0 . (15) ESAIM: PROCEEDINGS AND SURVEYS
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 atthe step, i.e. E fluid + E solid = (cid:90) t (cid:16) f ext | x = − l + f ext | x =0+ − f ext | x =0 − (cid:17) , with the fluid energy defined by E fluid = (cid:90) I e int + (cid:90) E e ext . This assumption is an adaptation to a bounded domain case of the conservation of total fluid-structureenergy 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 factthat the structure is fixed ( ddt E solid = 0) yields ddt E fluid = (cid:90) I ∂ t e int + (cid:90) E ∂ t e ext = f ext | x = − l + f ext | x =0+ − f ext | x =0 − . From (14) and (15) we have − (cid:90) I ∂ x f int − (cid:90) E ∂ x f ext = f ext | x = − l + f ext | x =0+ − f ext | x =0 − . Using the boundary conditions (2) and (3) we get (cid:74) f int (cid:75) = (cid:74) f ext (cid:75) , where the brackets (cid:74) · (cid:75) are defined as in (12). By definition of the fluxes it follows (cid:74) q i P i (cid:75) = ρ (cid:115) q e h e + gζ e q e (cid:123) and from (2) and (12) we obtain (cid:74) P i (cid:75) = ρ (cid:115) q e h e + gζ e (cid:123) . Integrating on ( l + r, l − r ), the second equation in (9) yields − ρ rh w ddt q i = (cid:74) P i . (cid:75) Combining the last two equalities, we get the following transmission condition − rh w ddt q i = (cid:115) q e h e + gζ e (cid:123) . (16) Reformulation as two transmission problems
Coupling the governing equations (7)-(9) with the conditions derived in the previous section, we havetherefore reduced the problem of the OWC essentially to two transmission problems. The first one in E ∪ E reads: ∂ t ζ e + ∂ x q e = 0 ,∂ t q e + ∂ x (cid:18) q e h e (cid:19) + gh e ∂ x ζ e = 0 , h e = h s + ζ e in E , h e = h + ζ e in E , (17) SAIM: PROCEEDINGS AND SURVEYS x = 0 ζ e | x =0 − = ζ e | x =0+ , q e | x =0 − = q e | x =0+ . (18)The second transmission problem in E ∪ E reads: ∂ t ζ e + ∂ x q e = 0 ,∂ t q e + ∂ x (cid:18) q e h e (cid:19) + gh e ∂ x ζ e = 0 , h e = h + ζ e , (19)with transmission conditions at x = l ± r (cid:74) q (cid:75) = 0 , − α ddt q i = (cid:115) q e h e + gζ e (cid:123) , (20)where α = 2 rh w and h w = h + ζ w . Riemann invariants
Let us now rewrite the nonlinear shallow water equations (7) and (8) in the exterior domain E in acompact form by introducing the couple U = ( ζ e , q e ) T : ∂ t U + A ( U ) ∂ x U = 0 , (21)where A ( U ) = (cid:32) gh e − q e h e q e h e (cid:33) . The eigenvalues λ + ( U ) and − λ − ( U ) of the matrix A ( U ) and the associated eigenvectors e + ( U ) and e − ( U )are given by λ + ( U ) = q e h e + (cid:112) gh e , − λ − ( U ) = q e h e − (cid:112) gh e ,e + ( U ) = (cid:18)(cid:112) gh e − q e h e , (cid:19) T , e − ( U ) = (cid:18) − (cid:112) gh e − q e h e , (cid:19) T . Notice that λ + > λ − >
0. Taking the scalar product of (21) and eigenvectors, we obtain ∂ t (cid:18) (cid:112) gh e ± q e h e (cid:19) ± (cid:18)(cid:112) gh e ± q e h e (cid:19) ∂ x (cid:18) (cid:112) gh e ± q e h e (cid:19) = 0 . Let us introduce the right and the left Riemann invariant R and L associated to the nonlinear shallowwater equations, respectively R ( U ) := 2 (cid:16)(cid:112) gh e − (cid:112) gh (cid:17) + q e h e , L ( U ) := 2 (cid:16)(cid:112) gh e − (cid:112) gh (cid:17) − q e h e . (22)Hence we can write the 1D nonlinear shallow water equations in the exterior domain as the two followingtransport equations on R and L : ∂ t R ( U ) + λ + ( U ) ∂ x R ( U ) = 0 , ∂ t L ( U ) − λ − ( U ) ∂ x L ( U ) = 0 . (23)We will see that these two transport equations of Riemann invariants are helpful when we solve ourmodel by numerical method. More details about Riemann invariants of the nonlinear shallow waterequations can be found in [9]. ESAIM: PROCEEDINGS AND SURVEYS Discretization of the model
We have reformulated in the previous section the mathematical model of the oscillating water columnas 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 tosolve our main equations and use Riemann invariants to address the entry conditions and all boundaryconditions.3.0.1.
Numerical notations
We use the following notations throughout this section: • in our system, the whole numerical domain [ − l, l ] is composed of four parts: [ − l, , l − r ],[ l − r, l + r ] and [ l + r, l ]. Each interval is divided into cells ( A i ) ≤ i ≤ n x with A i = [ x i − , x i ] ≤ i ≤ n x of size δ x . More precisely, we have x = − l, ..., x i = − l + iδ x , ..., x n ,x = 0; x n ,x +1 = δ x , ..., x n ,x + i = iδ x , ..., x n ,x + n ,x = l − r ; x n ,x + n ,x +1 = l − r + δ x , ..., x n ,x + n ,x + i = l − r + iδ x , ..., x n ,x + n ,x + n ,x = l + r ; x n ,x + n ,x + n ,x +1 = l + r + δ x , ..., x n ,x + n ,x + n ,x + i = l + r + iδ x , ..., x n ,x + n ,x + n ,x + n ,x = l , with l = n ,x δ x , l − r = n ,x δ x , r = n ,x δ x and l − ( l + r ) = n ,x δ x ; • we denote by δ t the time step. According to CFL condition, time step δ t can be specified by δ x ; • for any quantity U , we denote by U mi its value at the position x i at time t m = mδ t . For instance,the variables ζ mi denotes the value of the free surface elevation ζ at the position x i at time t m = mδ t . 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 conser-vative form : ∂ t U + ∂ x ( F ( U )) = 0 , (24)with F ( U ) = (cid:18) q e , g (cid:0) h e − h (cid:1) + q e h e (cid:19) T . By means of a finite difference approach, equation (24) can be discretized as U m +1 i − U mi δ t + F mi +1 / − F mi − / δ x = 0 , where the flux F is discretized with cell centres indexed as i and cell edge fluxes indexed as i ± / F mi ± / depends on the numerical scheme. We consider here the well-known Lax–Friedrichsscheme proposed by Lax [10] to get the discrete flux F mi − / = 12 (cid:0) F mi + F mi − (cid:1) − δ x δ t (cid:0) U mi − U mi − (cid:1) , (25)where i ≥ F mi = F ( U mi ). SAIM: PROCEEDINGS AND SURVEYS Discretization of the entry condition
At the entrance of our system, the surface elevation is given by a prescribed function f depending onlyon time, ζ m | x = − l = f m := f ( t m ) . In order to express the entry condition for the horizontal discharge, let us first recall that from (22) onehas q e = h e R − L ) , R + L = 4 (cid:16)(cid:112) gh e − (cid:112) gh (cid:17) , where R and L are respectively the right and the left Riemann invariant associated to the nonlinearshallow water equations. We get q e = h e (cid:16) (cid:16)(cid:112) gh e − (cid:112) gh (cid:17) − L (cid:17) . Hence, the value of q e at x = − l is given by q e | x = − l = ( h + f ( t )) (cid:16) (cid:16)(cid:112) g ( h + f ( t )) − (cid:112) gh (cid:17) − L | x = − l (cid:17) . 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). Afterdiscretizing it as in [12], we get L m − L m − δ t − λ − L m − − L m − δ x = 0 , (26)where L m is the value of L at x = − l at time t m and λ − is computed as a linear interpolation between λ − , and λ − , following [9], namely λ − = βλ − , + (1 − β ) λ − , with 0 ≤ β ≤ λ − δ t = βδ x . Moreover, we can compute λ − as λ − = λ − , δ t δ x λ − , − δ t δ x λ − , . Thus, we have L m = (cid:18) − λ − δ t δ x (cid:19) L m − + λ − δ t δ x L m − , (27)which gives L m 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 takeninto consideration besides the entry condition at x = − l . When wave arrives from the offshore, it willencounter a step in the bottom and then arrive into a chamber, and finally arrive to the wall (see theconfiguration 1). More precisely, the first boundary condition is at the discontinuity of the topographylocated at x = 0 and the second is at the partially immersed structure side-walls located at x = l ± r .The last boundary condition is at the end of the chamber, located at x = l .0 ESAIM: PROCEEDINGS AND SURVEYS
At the topography step
Let us first consider the shallow water wave equations with discontinuous topography, namely, it is asystem with depth h s on R − = { x < } and depth h on R + = { x > } . Our equation turns out to be ∂ t U + ∂ x ( F ( U )) = 0 , with F ( U ) = (cid:18) q e , g (cid:0) ( h s + ζ e ) − h s (cid:1) + q e h s + ζ e (cid:19) T , in (0 , T ) × R − , (cid:18) q e , g (cid:0) ( h + ζ e ) − h (cid:1) + q e h + ζ e (cid:19) T , in (0 , T ) × R + . From transmission conditions (10) and (11), we have the continuity of the surface elevation ζ e and ofthe horizontal discharge q e at x = 0: ζ le | x =0 = ζ re | x =0 , q le | x =0 = q re | x =0 . (28)Let us denote the right Riemann invariant in the domain R − by R l and the left Riemann invariant inthe domain R + by L r . We then find two expressions of q e describing q le | x =0 and q re | x =0 , respectively, q le | x =0 = (cid:0) h s + ζ le | x =0 (cid:1) (cid:18) R l | x =0 − (cid:18)(cid:113) g ( h s + ζ le | x =0 ) − (cid:112) gh s (cid:19)(cid:19) ,q re | x =0 = ( h + ζ re | x =0 ) (cid:16) (cid:16)(cid:112) g ( h + ζ re | x =0 ) − (cid:112) gh (cid:17) − L r | x =0 (cid:17) . (29)According to the relations (28), we observe that (29) is a system of two nonlinear equations on the twounknowns ζ le | x =0 (respectively ζ re | x =0 ) and q le | x =0 (respectively q re | x =0 ). We write it in the compact form F ( x , x ) = 0 , (30)where x = ζ le | x =0 , x = q le | x =0 and the vector F = ( F , F ) is given by F = ( h s + x ) (cid:16) R l | x =0 − (cid:16)(cid:112) g ( h s + x ) − (cid:112) gh s (cid:17)(cid:17) − x ,F = ( h + x ) (cid:16) (cid:16)(cid:112) g ( h + x ) − (cid:112) gh (cid:17) − L r | x =0 (cid:17) − x . In the case h s = h (without step) we can derive from (29) a third degree equation on (cid:112) h + ζ le | x =0 andtake the unique solution that gives ζ le | x =0 = 0 when R l | x =0 , L r | x =0 = 0 (we refer to [8] for this case).Here, since h s (cid:54) = h , 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: ∂ t R l + λ l + ( U ) ∂ x R l = 0 , ∂ t L r − λ r − ( U ) ∂ x L r = 0 , (31)where the corresponding eigenvalue λ l + in the domain R − is given by λ l + ( U ) = q e h s + ζ e + (cid:112) g ( h s + ζ e ) , (32)and the corresponding eigenvalue − λ r − in the domain R + is given by − λ r − ( U ) = q e h + ζ e − (cid:112) g ( h + ζ e ) . (33) SAIM: PROCEEDINGS AND SURVEYS λ + and λ − as in [12]. After discretizationof equations (31), we get( R l ) mn ,x − ( R l ) m − n ,x δ t + λ l + ( R l ) m − n ,x − ( R l ) m − n ,x − δ x = 0 , ( L r ) mn ,x − ( L r ) m − n ,x δ t − λ r − ( L r ) m − n ,x +1 − ( L r ) m − n ,x δ x = 0 , where λ l + , λ r − are as in (32)-(33) and we recall that ( R l ) mn ,x is the value of R l at x n ,x and t m (seeNotations 3.0.1). Hence, we have( R l ) mn ,x = (cid:18) − λ l + δ t δ x (cid:19) ( R l ) m − n ,x + λ l + δ t δ x ( R l ) m − n ,x − , ( L r ) mn ,x = (cid:18) − λ r − δ t δ x (cid:19) ( L r ) m − n ,x + λ r − δ t δ x ( L r ) m − n ,x +1 , (34)which give ( R l ) mn ,x and ( L r ) mn ,x in terms of their values at the previous time step and in terms of interiorpoints.Gathering the relations (28), (29) and (34), we can solve ζ le | x =0 (respectively ζ re | x =0 ) and q le | x =0 (respectively q re | x =0 ), which give us the boundary conditions at the step.3.3.2. At the structure side-walls
Compared with the derivation of the boundary conditions near the step, the idea to derive the boundarycondition near the fixed partially immersed structure is almost the same. There are two differencesbetween them. The first one is that, since the depth is always h , the eq. (29) becomes q le | x = l − r = ( h + ζ le | x = l − r ) (cid:18) R l | x = l − r − (cid:18)(cid:113) g ( h + ζ le | x = l − r ) − (cid:112) gh (cid:19)(cid:19) ,q re | x = l + r = ( h + ζ re | x = l + r ) (cid:16) (cid:16)(cid:112) g ( h + ζ re | x = l + r ) − (cid:112) gh (cid:17) − L r | x = l + r (cid:17) , (35)where we denote the horizontal discharge in the exterior domain on the left-hand side of the object by q le and on the right-hand side of the object by q re . Let us recall that q i is the horizontal discharge in theinterior domain I . From the first transmission condition in (20), we know that q le | x = l − r = q i = q re | x = l + r . The second difference is that, unlike in the previous subsection, we do not have the continuity conditionof ζ e at the structure side-walls. Nevertheless, we consider the discretization of the second transmissioncondition in (20), hence we get − α ( q e ) ml − r − ( q e ) m − l − r δt = (cid:0) ( q le ) m − l + r (cid:1) (cid:0) h + ( ζ le ) m − l + r (cid:1) + g ( ζ le ) m − l + r − (cid:0) ( q re ) m − l − r (cid:1) (cid:0) h + ( ζ re ) m − l − r (cid:1) − g ( ζ re ) m − l − r . where for the sake of clarity ( q e ) ml − r = ( q e ) mn ,x + n ,x and ( q e ) ml + r = ( q e ) mn ,x + n ,x + n ,x (analogously for( ζ e ) ml − r and ( ζ e ) ml + r ). Then, q e at x = l − r is expressed as( q e ) ml − r = ( q e ) m − l − r − δtα (cid:32) (cid:0) ( q le ) m − l + r (cid:1) (cid:0) h + ( ζ le ) m − l + r (cid:1) − (cid:0) ( q re ) m − l − r (cid:1) (cid:0) h + ( ζ re ) m − l − r (cid:1) (cid:33) − δtα g (cid:0) ( ζ le ) m − l + r − ( ζ re ) m − l − r (cid:1) , (36)which gives ( q e ) ml − r in terms of its values at the previous time step and in terms of interior points. Nowwe can solve ( q e ) ml − r immediately. Once the value of ( q e ) ml − r is obtained, we can find the values of ζ le | x = l − r and ζ re | x = l + r by using equations (35) and the transport equations for the Riemann invariantsas the strategy in Section 3.3.1.2 ESAIM: PROCEEDINGS AND SURVEYS
At the end of the chamber
The corresponding boundary condition at the end of the chamber, located at x = l , is given by q e | x = l = 0 . Hence, recalling the definition of the right-going Riemann invariant R , we recover the surface elevation ζ e at x = l , namely ζ e | x = l = 1 g (cid:18) R | x = l (cid:112) gh (cid:19) − h . Numerical validations
In this section, we use the scheme introduced in Section 3 to simulate our model. For the fluid, wealways consider the density of water ρ = 1000 kg / m and the gravitational acceleration g = 9 .
81 m / s .The entry of the domain is set at x = − l = −
30 m and the prescribed function f is given by f ( t ) = sin (cid:18) πT t (cid:19) , where T = 1 . l = 11 m, r = 1 m and l = 17 m and the fluid height at rest before the step h s = 15 m. We compute the solution by using theLax-Friedrichs scheme in the exterior domain [ − , ∪ [12 , N x = 2300 anda time step δ t = . √ gh s δ x with space step δ x = 0 .
02 m. Here, the CFL number is 0 .
7, which is commonlyused to prescribe the terms of the finite-difference approximation of a PDE (see for instance [13]). In theinterior domain, the solution can be computed using the transmission conditions (20) with h w = h + ζ w and ζ w = − . Numerical solutions
In real applications, an OWC device can be deployed on a stepped sea bottom in order to improve itsperformance. It is important then to have a good understanding of the impact of a step in the topography.Here, we test and compare the case without step s = 0 m . ( h = 15 m) to the case with a step of height s = 5 m ( h = 10 m ) considering the previous physical parameters. The numerical solutions are plottedin Figure 2 at times t = 1 . t = 3 . t = 5 s. The plots (a), (c), (e) show the solutions withoutstepped bottom, while the plots (b), (d), (f) show the solutions with stepped bottom.We find that, before the waves encounter the step, there is no significant difference between the OWCmodel without stepped bottom and with stepped bottom (see (a) and (b)). But when the waves encounterthe step in the bottom and arrive into the chamber, we can see that, the waves in the OWC model withoutstepped bottom move significantly faster than the waves in the OWC model with stepped bottom. Inparticular, at t = 3 . t = 5 s, the reflected wave in the OWC model without stepped bottom already reaches x = −
10 m, whilethe 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 wavesin the OWC model without stepped bottom.This difference can be explained by the fact that more incident wave energy is converted when a stepis added. In other words, the OWC with stepped bottom would be more efficient than the one withoutstepped bottom, which is in agreement with the result by Rezanejad and Soares in [14].
SAIM: PROCEEDINGS AND SURVEYS -30 -20 -10 0 10 X -15-10-50510 Z t = 1.7 q / 5 -30 -20 -10 0 10 X -15-10-50510 Z t = 1.7 q / 5 (a) (b) -30 -20 -10 0 10 X -15-10-50510 Z t = 3.3 q / 5 -30 -20 -10 0 10 X -15-10-50510 Z t = 3.3 q / 5 (c) (d) -30 -20 -10 0 10 X -15-10-50510 Z t = 5 q / 5 -30 -20 -10 0 10 X -15-10-50510 Z t = 5 q / 5 (e) (f) Figure 2.
Comparisons between the numerical results without step (left) and with step(right) at times t = 1 . t = 3 . t = 5 s.4 ESAIM: PROCEEDINGS AND SURVEYS X -1.5-1-0.500.511.5 Z t = 1.7 Classic NSW ModelNSW without step X -1.5-1-0.500.511.5 Z t = 3.3 Classic NSW ModelNSW without step (a) (b) X -1.5-1-0.500.511.5 Z t = 5 Classic NSW ModelNSW without step (c)
Figure 3.
Comparison between classical NSW model and our model without step in differenttimes considering δ x = L/ L = 30 m. Accuracy analysis
In numerical validations, accuracy analysis is of importance. As we can see in Figure 1, the configu-ration of OWC device is essentially constituted from three parts: the domain before the step in the seabottom, the domain after the step and the chamber. We implement our algorithms by gathering togetherthe 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. Underthe same initial wave and physical parameters, we compare the free surface elevation ζ e of the classicalnonlinear shallow water wave model with our model without discontinuous topography. Figure 3 showsthat there is no significant difference between the two cases. Moreover, we also find that the error is oforder 10 − (see Figure 4), which is acceptable since the Lax-Friedrichs method is first-order accurate inspace. 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 theefficiency of the OWC considered in this work.
SAIM: PROCEEDINGS AND SURVEYS Time -1-0.500.51 D i ff e r en c e -3 Figure 4.
Difference between classical NSW model and our NSW model without step consid-ering δ x = L/ L = 30 m. The primary efficiency η Reg of the device is defined by the ratio of the absorbed power from the wavesto the incident wave power. From the seminal work of Evans in [6], we know that in the linear time-harmonic theory the volume flux Q ( t ) = Re { qe − iωt } is assumed linearly proportional to the pressure inthe chamber P ( t ) = Re { pe − iωt } . Using this assumption, the average power absorbed from regular wavesover one wave period, denoted by P Reg , is given by P Reg = 12 λ | p | , (37)where p is the time independent and λ is a positive constant associated with linear air turbine character-istics. On the other hand, following [14] in experiments the average power absorbed from regular wavescan be determined by: P Reg = 1 T (cid:90) T P Qdt, (38)where T is the duration of the test. The incident wave power P inc is defined as the product of totalenergy per wave period E inc and the group velocity c g (see [3]): P inc = E inc c g , with E inc = 12 ρgLA , c g = ω k (cid:18) kh s sinh(2 kh s ) (cid:19) , and the dispersion relation given by ω = gk tanh( kh s ) , where ω is the frequency, k is the wave number, h s is the fluid height at rest before the step, ρ is the densityof the fluid, g is the gravitational acceleration and L is the projected width of the WEC perpendicularto the incident wave direction, A is the amplitude of the wave. In the shallow water regime kh s (cid:28) c g = √ gh s . Thus, the primary efficiency of the device in regular wave isgiven by η Reg = P Reg P inc . We notice that in both (37) and (38) the absorbed power (hence the primary efficiency) stronglydepends on the air pressure in the chamber. In our model, it is considered to be a constant, namely theatmospheric pressure P atm . However, when the waves arrive into the chamber and change the volume of6 ESAIM: PROCEEDINGS AND SURVEYS the air, the air pressure in the chamber will certainly change as well. In this case, the pressure will nomore be a constant, but depends on time. Hence, to study more rigorously the absorbed power and theprimary efficiency of the OWC, this fact must be taken into account in the model. This will be addressedin our future work. Analogously, the improvement of the efficiency of an OWC device deployed on astepped 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 beachieved by adding a step at the bottom of the sea.
Acknowledgements
The authors warmly thank David Lannes for his helpful comments and advises and also the organisersof the summer school CEMRACS 2019 during which this work was done. E.B. is supported by the StartingGrant project “Analysis of moving incompressible fluid interfaces” (H2020-EU.1.1.-639227) operated bythe European Research Council. J. H. is supported by the PostDoc program Sophie Germain of theFondation Math´ematique Jacques Hadamard. G. V-H. is supported by the European Union’s Horizon2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 765579.
References [1]
E. Bocchi , Floating structures in shallow water: local well-posedness in the axisymmetric case , SIAM J. Math. Anal.,52 (1) (2020), pp. 306–339.[2]
D. Bresch, D. Lannes, and G. M´etivier , Waves interacting with a partially immersed obstacle in the boussinesqregime , to appear in Analysis & PDE.[3]
R. G. Dean and R. A. Dalrymple , Water wave mechanics for engineers and scientists , vol. 2, World ScientificPublishing Company, 1991.[4]
A. Elhanafi, A. Fleming, G. Macfarlane, and Z. Leong , Numerical energy balance analysis for an onshoreoscillating water column–wave energy converter , Energy, 116 (2016), pp. 539–557.[5] ,
Numerical hydrodynamic analysis of an offshore stationary–floating oscillating water column–wave energyconverter using cfd , International Journal of Naval Architecture and Ocean Engineering, 9 (2017), pp. 77–99.[6]
D. Evans , Wave-power absorption by systems of oscillating surface pressure distributions , Journal of Fluid Mechanics,114 (1982), pp. 481–499.[7]
T. Iguchi and D. Lannes , Hyperbolic free boundary problems and applications to wave-structure interactions , toappear in Indiana University Mathematics Journal.[8]
D. Lannes , On the dynamics of floating structures , Annals of PDE, 3 (2017), pp. 11–81.[9]
D. Lannes and L. Weynans , Generating boundary conditions for a boussinesq system , Nonlinearity, 33 (2020),pp. 6868—-6889.[10]
P. D. Lax , Weak solutions of nonlinear hyperbolic equations and their numerical computation , Communications onpure and applied mathematics, 7 (1954), pp. 159–193.[11]
D. Maity, J. San Mart´ın, T. Takahashi, and M. Tucsnak , Analysis of a simplified model of rigid structure floatingin a viscous fluid , Journal of Nonlinear Science, 29 (2019), pp. 1975–2020.[12]
F. Marche , Theoretical and numerical study of shallow water models: applications to nearshore hydrodynamics , PhDthesis, Bordeaux 1, 2005.[13]
N. ¨Ozis¸ik, H. R. Orlande, and M. Colac¸o , Finite difference methods in heat transfer , CRC press, 2017.[14]
K. Rezanejad and C. G. Soares , Enhancing the primary efficiency of an oscillating water column wave energyconverter based on a dual-mass system analogy , Renewable Energy, 123 (2018), pp. 730 – 747.[15]