Nonlinear analysis of an actuated seafloor-mounted carpet for a high-performance wave energy extraction

Mohammad-Reza Alam


It is known that muddy seafloors can extract significant energy from overpassing surface waves via engaging them in strong interaction processes. If a synthetic seabed can respond to the action of surface gravity waves similar to the mud response, then it too can take out a lot of energy from surface waves. Analysis of the performance of a mud-resembling seabed carpet in harvesting ocean wave energy is the subject of this article. Specifically, and on the basis of the field measurements and observations of properties/responses of seafloor mud, we focus our attention on an artificial viscoelastic seabed carpet composed of (vertically acting) linear springs and generators. We show that the system of sea/synthetic-carpet admits two propagating wave solutions: the surface mode and the bottom mode. The damping of a surface-mode wave is proportional to its wavelength and hence is classic. However, the damping of a bottom-mode wave is larger for shorter waves, and is in general stronger than that of the surface-mode wave. To address the effect of (high-order) nonlinear interactions as well as to investigate the performance of our proposed carpet of wave energy conversion (CWEC) against a spectrum of waves, we formulate a direct simulation scheme based on a high-order spectral method. We show, by taking high-order nonlinear interactions into account, that the CWEC efficiency can be significantly higher for steeper waves. We further show that the bandwidth of high performance of the CWEC is broad, it yields minimal wave reflections and its theoretical efficiency asymptotically approaches unity within a finite and (relatively) short extent of deployment.

1. Introduction

Gade (1958) reports a place in the Gulf of Mexico, known to locals as mud hole, where owing to the accretion of mud banks has turned into, for the local fishermen, a safe haven against strong waves during storms. Within the mud hole, the interaction of surface waves with the mud is very strong such that waves are completely damped out within a few wavelengths. Observations of strong wave damping owing to the coupling with the bottom mud are not limited to the Gulf, but almost anywhere with a very muddy seafloor (Silvester 1974; Macpherson & Kurup 1981; Sheremet & Stone 2003). If mud can take a substantial energy out of incident surface gravity waves, an artificial carpet deployed on the seabed that responds to the action of the overpassing waves in the same way (as the response of a mud layer) must be able to extract the same amount of energy. Analysis of performance of this synthetic seabed-carpet wave energy conversion technique (CWEC) is the subject of this article.

The complicated nature of the seafloor mud and the wide range of its mechanical/material properties that may be location and even time dependent, has aroused a great deal of research on this subject. For understanding wave–mud interactions, several models have been incorporated (e.g. Dalrymple et al. 1978; Macpherson 1980; Liu et al. 1993; Lin & Jeng 2003) with the correct model in general yet a matter of dispute (Winterwerp & Kesteren 2004). However, under the periodic forcing, a viscoelastic model has been shown to be a very good approximation and is now widely used (Jiang & Mehta 1995; Mei et al. 2010). While the general idea presented here can incorporate any mud model, for specificity, we focus our attention on a linear viscoelastic mud model.

In this paper, we consider a seabed carpet composed of sets of vertically acting linear springs and generators, with the generator's action modelled to be linearly proportional to the vertical speed. We show that the coupled governing equation of the waves/carpet system admits two modes of propagating wave solution (§2): a surface-mode wave and a bottom-mode wave. For a surface mode, the higher the wavelength, the higher the energy extraction by the bottom. This is in agreement with most observations of long waves damped by the bottom mud (Liu et al. 1993). For a bottom-mode wave, however, the energy extraction increases as the wavelength decreases (§§2 and 3). Besides the attractiveness of this feature of a bottom-mode wave for the wave energy community (Farley 2012; Falnes & Hals 2012), it may provide an explanation for the recent (and yet unexplained) observation of strong short-wave damping by muddy seafloors (Sheremet & Stone 2003; Alam et al. 2011).

To address the effect of nonlinear interactions on the efficiency of the proposed wave energy extraction technique, and also to show the performance of the scheme against a powerful spectrum of waves, we formulate a computational high-order spectral technique (§4). It is shown that the efficiency of the viscoelastic carpet in extracting the wave energy is higher for steeper waves (more than 30% in some cases). We also show that our energy-harnessing scheme acts (almost) uniformly across the frequency band of a broadband spectrum, extracting most of the energy while yielding (almost) no reflected waves.

The presented wave energy conversion device is completely under the water surface, hence imposes minimal danger to boats and to the sea life (i.e. no mammal entanglement). The carpet is survivable against high momentum of storm surges and, in fact, can perform even better under very energetic (e.g. stormy) sea conditions, when most existing wave energy devices are needed to shelter themselves by going into an idle mode. The proposed idea and its variations may also be used to create localized safe havens for fishermen and sailors in open seas, or if implemented in large scales to protect shores and harbours against strong storm waves.

2. Governing equations

We consider the irrotational motion of a homogeneous inviscid incompressible fluid with a free surface over a viscoelastic bottom with the mean depth z=−h (figure 1). Governing equations in terms of velocity potential ϕ and surface/bottom elevations ηs,ηb, ignoring surface tension, read Embedded Image 2.1a Embedded Image 2.1b Embedded Image 2.1c Embedded Image 2.1d Embedded Image 2.1e Embedded Image 2.1f where Pb is the pressure on the seabed, b* and k* are, respectively, the viscous damping and stiffness coefficient of the viscoelastic bottom per unit area, ρ is the density and g is the gravity acceleration. In governing equations (2.1), (2.1a) is the continuity equation in the fluid domain; (2.1b), (2.1c) are, respectively, kinematic and dynamic boundary conditions on the free surface; (2.1d), (2.1e) are, respectively, kinematic and dynamic boundary conditions on the viscoelastic bottom and (2.1f) represents the dynamics of the viscoelastic bottom.

Figure 1.

Schematic of a synthetic viscoelastic carpet on the seafloor for extracting energy from surface gravity waves. The carpet is composed of linear springs (stiffness coefficient k*) that provide the restoring force, and generators (modelled by dash-pot type dampers with damping coefficient b*) that extract energy. The distance between each module of the spring damper is assumed much smaller than the typical wavelength of overpassing waves such that the assumption of continuously distributed spring dampers is valid. (Online version in colour.)

The linearized form of governing equations ((2.1) admits a propagating wave solution in the form (cf. Alam et al. 2009a) Embedded Image 2.2a Embedded Image 2.2b Embedded Image 2.2c with any arbitrary phase (it is assumed zero here), and with Embedded Image 2.3a Embedded Image 2.3b where as,ab are, respectively, surface and bottom amplitudes. Variables k,ω are, respectively, the wavenumber and frequency of the propagating wave and must satisfy the dispersion relation, which, in dimensionless form, reads Embedded Image 2.4 with Embedded Image 2.5 in which the variable Ω is the dimensionless frequency, ζ is the dimensionless damping ratio, γ is the dimensionless restoring force (the ratio of hydrostatic restoring coefficient to the bottom elasticity) and μ is shallowness. Note that the limiting case of Embedded Image (i.e. Embedded Image) is equivalent to a rigid bottom and (2.4) readily reduces to the rigid and flat bottom dispersion relation.

In the absence of damping (i.e. ζ=0), the necessary condition for having a stable (i.e. non-growing) solution to (2.4) is γ<1. In physical space, this requirement is equivalent to k*>ρg, which means that, if the bottom is perturbed, the added force on the bottom owing to the weight of the water column perturbation must be compensated by the restoring force of bottom elasticity. Under this necessary condition and for any given μ, (2.4) with ζ=0 attains two real solutions for Ω. The two undamped solutions of (2.4) in the limit of deep water (i.e. μ≫1) are further simplified to Embedded Image 2.6 and Embedded Image 2.7 Equation (2.6) is clearly the dispersion relation of deep-water waves and is not affected by the bottom elasticity. Note that equation (2.6) (surface mode in the limit of deep water) results in abas (carpet amplitude is much smaller than the surface) and hence is of less physical significance.

Undamped solutions of (2.4) along with the ratio of bottom amplitude to the surface amplitude (2.3b) are plotted in figure 2 (dashed-dotted-dotted and solid lines) for three values of γ=0.1,0.5,0.9. Also plotted in this figure are asymptotic curves (2.6) and (2.7) and their respective amplitude ratios represented by, respectively, dashed and dashed-dotted lines. For ease of reference, we refer to the branch of (2.4) that at the limit of μ≫1 tends to (2.6) as the surface-mode solution, and the branch that, at this limit, tends to (2.7), as the bottom-mode solution. It can be shown that for a surface-mode solution as>ab, and for a bottom-mode solution as<ab, except for γ=γcr=0.5, for which as μ increases, the surface and bottom solutions asymptotically tend to each other and towards the deep-water dispersion relation (2.6), as is shown in figure 2b and can be confirmed by equations (2.6) and (2.7).

Figure 2.

Plot of branches of the undamped dispersion relation (equation (2.4) with ζ=0) and ratio of bottom elevation to the surface elevation (2.3b). (ac) Surface-mode (dashed-dotted-dotted-line) and bottom-mode (solid line) solutions to the dispersion relation (2.4) for, respectively, γ=0.1, 0.5, 0.9. Asymptotic relations (2.6) (dashed line) and (2.7) (dashed-dotted line) are also plotted for comparison. Note that for γ=γcr=0.5, all curves asymptotically tend to the deep-water dispersion relation (2.6). (df) Ratios of bottom elevation amplitude ab to the surface elevation amplitude as (2.3b) for, similarly, γ=0.1, 0.5, 0.9. Line styles of each of panels (df), respectively, follow their top figure counterparts. (Online version in colour.)

For a sub-critical 0<γsub<γcr, any given wavenumber is associated with one surface-mode wave and one bottom-mode wave, with the frequency of the bottom-mode wave being much higher than that of the surface mode (figure 2a). However, the frequency of bottom-mode waves cannot be below the threshold of Embedded Image. The inverse scenario holds for the supercritical Embedded Image where now the frequency of a surface-mode wave is limited by the Embedded Image. Note that the substitution of values of k,ω from figure 2 into the assumed form of the solution (2.2) results in a right-going wave. For left-going waves, the values of k,ω are to be picked from the second half of plots in figure 2 (not shown), which are left or bottom mirrors of presented plots (Phillips 1960; Ball 1964).

If damping is non-zero (ζ>0) and in the limit of deep water (μ≫1), the dispersion relation (2.4) is simplified to Embedded Image 2.8 where, similar to the undamped case, the first/second parenthesis is referred to as the surface/bottom-mode solution. Clearly, bottom damping has no effect on the surface-mode waves but affects the bottom mode. In the presence of damping, the frequency may gain an imaginary part that corresponds to exponential decay of the amplitude and accordingly the loss of energy. Real and imaginary parts of the solution Ω of (2.4) for a fixed γsup=0.9 and ζ=0.35,0.60 are shown in figure 3, along with the bottom-to-surface-amplitude ratios and deep-water asymptotes.

Figure 3.

(a,b) Real and (c,d) imaginary branches of the solution Ω(μ) to the dispersion relation (2.4) for γ=0.9, and (a,c) ζ=0.35 and (b,d) ζ=0.60. (e,f) Also plotted is the ratio of the bottom amplitude to the surface amplitude. Plotted curves are branches associated with the surface mode (dashed-dotted-dotted line) and the bottom mode (solid line). Note the bifurcation at μcr=3.62 and μcr=1.23 (cf.(2.9)). Branches of deep-water asymptotes (2.8) are also plotted for comparison (dashed and dashed-dotted lines). (Online version in colour.)

The real part of the frequency of surface-mode waves, ℜ(Ωs), experiences minimal change compared with the undamped case (figure 3a,b, dashed-dotted-dotted line). The imaginary part of a surface-mode wave, ℑ(Ωs), starts from a finite negative value (corresponds to finite decay rate) and asymptotically decreases to zero as waves get shorter (figure 3c,d, dashed-dotted-dotted line). This is expected and in agreement with existing theories and observations of stronger long-wave (compared with short-wave) damping.

The behaviour of the bottom-mode waves is however more complex. Roots of the bottom-mode bracket of (2.8) show a bifurcation at a critical (dimensionless) depth1 (cf. figure 3cf, solid line), Embedded Image 2.9 Specifically, if μ<μcr, then Ωb=ix3±x4, where Embedded Image. This is similar to the surface-mode case in which waves propagate with the frequency of ℜ(Ωb) and at the same time, exponentially get damped with the rate Embedded Image. However, |ℑ(Ωb)|increases as the wavenumber increases. This means that for a weak enough damping coefficient, short waves are damped faster than long waves, an opposite of what happens to surface-mode waves. The region of μ<μcr=3.62 in figure 3a,c,e and μ<μcr=1.23 in figure 3b,d,f correspond to this scenario.

For μ>μcr, Ωb=ix1,ix2 with Embedded Image, that is, the waves become overdamped and cannot propagate further. In this case, depending on the value/mode of Ωb (determined by the initial condition), shorter waves may be damped much faster (lower branch) or asymptotically approach a wavelength-independent decay rate. In either case, again very short waves are still damped, and are stronger than the rate of decay of the longest surface waves. We finally note that at the critical shallowness, μ=μcr, Ωb=ix1=ix2 and there is only one damping rate.

In the limit of shallow water (μ≪1), the dispersion relation (2.4), correct to the order O(μ2), reduces to Embedded Image 2.10 If we further assume that ζ<1, then the solution to (2.10) is Embedded Image 2.11 This can be seen in the Embedded Image limit of figure 3ad. Note that Embedded Image is the natural frequency of the water column on a compliant bottom of stiffness k*.

We would like to note that fourth-order dispersion relations, similar to the form presented here but different in details and properties, are also obtained in two-fluid systems such as inviscid two-layer systems (Alam et al. 2009a) or for an inviscid fluid over a viscoelastic fluid (Macpherson 1980).

3. Energy extraction rate

Generators attached to the flexible seabed-carpet extract energy from the surface gravity waves as the waves march forward over, and interact with, the CWEC. Assume a monochromatic progressive surface wave in the system of water over a CWEC carpet. The total energy of the undamped system is composed of three parts: kinetic energy of water particles, potential energy of water and potential energy of the carpet owing to its elasticity. These (kinetic and potential energies per unit length) are, respectively, given by Embedded Image 3.1a and Embedded Image 3.1b and therefore the total energy Embedded Image, where from (3.1), D is clearly a (dimensionless) constant and a function of water depth, wavenumber and the carpet properties. Note that in the limit of rigid bottom when ab=0, the dispersion relation tends to that of the finite depth gravity waves, the kinetic and potential energy of regular gravity waves are retained.

In the presence of damping, the solution to the dispersion relation (2.4) gains an imaginary part ω=ωr+i. The linear solution (2.2a) can therefore be written in the form Embedded Image 3.2 where Embedded Image. If the damping is weak, the rate of the total energy decay is obtained from Embedded Image 3.3 If we define the normalized energy Embedded Image where E0=1/2ρga2s0 and normalize the time with Embedded Image (in accordance with (2.5)), the initial rate of energy decay in dimensionless form turns into Embedded Image 3.4 A plot of Embedded Image as a function of damping coefficient ζ and shallowness μ is shown in figure 4 for the surface and bottom modes (dashed line), along with the numerical results to be discussed in §5. The initial rate of energy damping is finite for the surface mode and attains a maximum at a critical damping coefficient. This is expected because for a small damping coefficient, energy absorption is small, and for a very high damping coefficient, the bottom is basically too stiff to move and hence little energy will be harnessed. For the bottom mode, however, a higher damping coefficient results in a higher energy absorption. This happens because for a given surface amplitude, the bottom has to have a higher amplitude in the bottom mode. Note that for the bottom mode and beyond the bifurcation, no propagating wave will survive and all surface perturbations will simply damp to zero in place.

Figure 4.

The initial rate of energy damping Embedded Image in the (a) surface and (b) bottom mode of a monochromatic wave travelling over the CWEC carpet. The surface-mode damping is finite and (relatively) small, but gains a maximum for an optimum damping coefficient. The bottom-mode damping increases as the damping coefficient increases. Plots are for γ=0.9 and ϵ=0.001. Direct simulation parameters are M=1, N=128 and T/δt=104, where δt is the time step. (a,b) Solid black line, direct simulation (approximate); dashed-dotted line, direct simulation (exact); dashed line, theory. (Online version in colour.)

If an incident wave arrives at a patch of flexible carpet, the quantity of interest is the rate of change of the energy per spatial coordinate (i.e. dE/dx), which is related to (3.3) by Embedded Image 3.5 where Cg is the group velocity of the wave. Another quantity that is widely used for measuring the energy of a wave train is the energy flux FECg (cf. §5b).

4. Direct simulation

If the amplitude of the incident wave is small, linear theory predictions of energy absorption are satisfactory. But quite frequently, real ocean waves are steep, and the spectrum may not be narrow banded. For steep waves, the effects of nonlinearity in the governing equations become significant. Therefore, true evaluation of the total energy absorption requires a precise consideration of all nonlinear interactions as well as careful account of newly generated (free or locked) waves as they interact with the bottom. While it may be possible to follow the leading order nonlinear interaction of a few isolated waves, in general, tracking a system of many waves undergoing a high-order nonlinear interaction is algebraically tedious, if not impossible. To address the problem of the effect of high-order nonlinearities on the performance of the CWEC, as well as the efficiency of the device in extracting energy from a spectrum of waves here, we formulate a direct simulation scheme based on a high-order spectral method. The technique was originally developed independently by Dommermuth & Yue (1987) and West et al. (1987) for the problem of nonlinear wave–wave interactions in deep water. It was then extended to a number of problems including finite depth water (Liu & Yue 1998; Alam et al. 2010) and two-layer density stratified fluid (Alam et al. 2009b,c). Here, we formulate the high-order spectral scheme for the problem of waves propagating over a viscoelastic seabed and discuss the procedure of convergence test and validation.

Consider fully nonlinear governing equations (2.1). Let us define surface and bottom potentials Φs,Φb as the following: Embedded Image 4.1 In terms of (4.1), boundary conditions (2.1b)–(2.1f) conform to Embedded Image 4.2a Embedded Image 4.2b Embedded Image 4.2c Embedded Image 4.2d

In the spectral method presented here, at each step, we start from (assumed to be) known values of variables ηs,Φs,ηb,Φb. We will show that if values of these four variables are known, we can solve the boundary-value problem for ϕ to obtain the vertical velocity on the surface and the bottom, i.e. ϕz(x,z=ηs,t) and ϕz(x,z=−h+ηb,t). If this is achieved, then the right-hand side of equations (4.2a)–(4.2d) are in hand. Numerical integration is then used to integrate (4.2a)–(4.2d) and to provide initial conditions for the next time step.

Let us assume that the potential ϕ(x,z,t) can be expressed in terms of a perturbation series Embedded Image 4.3 where quantity X(m)O(ϵm) with ϵ being a measure of the wave steepness on the surface and bottom. The surface potential Φs(x,t) can then be Taylor expanded about the mean free surface. From (4.3), correct to the order O(ϵM), we obtain Embedded Image 4.4 If Φs(x,t) is known at the initial time t=t0, then equation (4.4) can be solved for f(m)ϕ(m)(x,z=0,t) sequentially, Embedded Image 4.5a and Embedded Image 4.5b Similar expressions are obtained for g(m)ϕ(m)(x,z=−h,t), Embedded Image 4.6a and Embedded Image 4.6b Equations (4.5) and (4.6) together with Laplace's equation at each order (∇2ϕ(m)=0) form a complete linearized boundary-value problem at each order m. These boundary-value problems can be solved sequentially up to an arbitrary order of nonlinearity M.

To construct the spectral method solution, we write the Fourier expansion of ϕ(m) in the form Embedded Image 4.7 where kn=2πn/Lx, Lx is the computational domain, and Embedded Image are the complex modal amplitudes that are determined from Dirichlet boundary conditions on z=0,h. Specifically, from (4.7),(4.5) and (4.6), we obtain Embedded Image 4.8a and Embedded Image 4.8b for n=0,±1,…,±N, where Embedded Image are, respectively, the nth Fourier amplitudes of f(m),g(m). Having Embedded Image, Embedded Image from (4.8), we can calculate the vertical velocity for each order m on the surface and the bottom, Embedded Image 4.9 and Embedded Image 4.10 and the overall vertical velocity at the surface and the interface is Embedded Image 4.11 and Embedded Image 4.12

To provide the right-hand side of the time-evolution equations (4.2), a pseudo-spectral method is used. Specifically, we calculate all spatial derivatives in the spectral space, and all nonlinear terms in the physical space at the discrete points of xi, i=1,2,…,2N. A fast Fourier transform is used to transform between physical and wavenumber spaces. Once the right-hand side of (4.2) is constructed, the fourth-order Runge–Kutta scheme is employed for the integration in time. The operation count for a perturbation of order M with N Fourier modes is Embedded Image.

To validate the accuracy and efficiency of this high-order spectral scheme, we use a fully nonlinear solution to the governing equation (2.1) obtained by Newton's iterative method (see Alam et al. (2009b) for details of the procedure for a similar problem). For the convergence test, this benchmark solution is accurate to 14 decimal places and is considered exact.

5. Numerical results and discussions

The direct simulation technique developed in §4 enables us to consider the general initial-value problem of waves interacting with the bottom. For a monochromatic wave train, for instance, depending on the initial condition and strength of the damping, energy may be distributed (or modulated) between either of the two linear modes of motion (i.e. surface or bottom). The numerical scheme further provides the capability to include (any order of) nonlinearity as well as interactions between a spectrum of waves and our CWEC carpet. Specifically, for relatively steep waves when effects of nonlinearities are important, part of the energy of original waves may flow into higher/lower harmonics, which may be revealed in the form of free propagating (in case of resonance, e.g. Alam (2012)) or locked waves. These newly generated waves themselves interact with the bottom and may pave the path for further depletion of energy from the initial spectrum (e.g. Alam et al. 2011). We first consider the problem of energy extraction from a monochromatic wave train and then study a case of a spectrum of waves losing its energy to our CWEC.

(a) Monochromatic wave train

Consider a surface-mode monochromatic wave train travelling over a viscoelastic seabed. For the initial wave profile, we consider an undamped linear solution of surface-mode waves (2.2). Clearly, if ζ=0 and nonlinear effects are small, the time evolution will result in a steady progressive profile without any change in the shape. However, when the damping effect of the bottom is included, the amplitude of the surface wave will, of course, decay as it loses energy to the seabed, and also because energy may flow into the bottom mode (which will also be damped, but at a higher rate). The latter phenomenon is the counterpart of mode coupling observed in most two-degree-of-freedom vibrating systems with damping (e.g. Rao 2004). For the problem at hand, unless in the limit of deep water (2.8), an incident surface wave exchanges energy with the bottom mode through damping. As will be shown later in this section, the effect of mode coupling on the overall damping increases with increase in damping coefficient.

For direct simulation purposes, consider a surface-mode wave train with the initial surface steepness ϵs=kas=0.001 travelling over a viscoelastic bottom of γ=0.9. The initial rate of energy decay Embedded Image as a function of damping coefficient ζ is plotted in figure 4a for four values of μ. The energy decay rate predicted by the linear theory of (3.4) (dashed line) is plotted against direct simulation results (solid and dashed-dotted lines). Specifically, to highlight the effect of mode coupling, we calculate the rate of energy decay from the numerical simulation in two different ways: (i) by neglecting the effect of mode coupling, we plug the amplitude of the surface-mode wave (obtained from numerical simulation2) into (3.1) (Embedded Image, solid line) and (ii) to obtain the true (or exact) energy decay, we sum the potential and kinetic energy over the entire fluid domain (Embedded Image, dashed-dotted line).

The linear results of (3.4) are expected to be very close to the numerical value of Embedded Image, as is confirmed in figure 4a where the two curves are graphically indistinguishable. These two also well approximate the exact rate of energy decay for (relatively) small values of ζ: all three curves overlap in the left side of figure 4a. However, as ζ increases, the approximate solutions underpredict the true damping. High damping clearly transfers energy from the surface mode to the bottom mode of the same frequency, where energy is damped with a higher rate. The discrepancy is highlighted more for shallower waters (smaller values of μ) as interaction with the bottom is stronger.

The inverse effect happens if we start with a bottom-mode wave (figure 4b). In this case again, Embedded Image and linear solution (3.4) match very well for all values. However, for (relatively) higher damping coefficients, the approximate solution overpredicts the true rate of energy decay. Here, converse to the previous case, part of the energy goes from the initial bottom-mode wave to the same frequency surface-mode wave, where the rate of energy decay is slower. Also note that Embedded Image is finite for a surface-mode wave and gains a maximum for an optimum ζ, while for a bottom mode, it increases as ζ increases and is also higher in value.

Nonlinearity is also affecting the initial rate of energy decay, and the direct simulation scheme developed here enables us to precisely comment on this aspect. To see the effect of order of nonlinearity on the initial rate of energy decay, we consider a (relatively) steep wave of ϵs=0.3. Note that we consider only the energy decay at the very early initial stage of the motion. Figure 5a compares Embedded Image as a function of ζ for different order of nonlinearity M. Linear theory and approximate linear simulation obtain the same normalized rate of energy decay, independent of the steepness (solid and dashed lines). The exact decay rate from linear simulation M=1 predicts a value higher than that of linear prediction, as is seen from figure 5a and was discussed earlier. As the considered order of nonlinearity M increases, the value of Embedded Image changes significantly until M=4, where the solution is converged in this example. The nonlinear damping rate may be as high as 30 per cent more than the prediction of linear theory. The trend is more highlighted in figure 5b, where Embedded Image is compared for different steepnesses of waves. Clearly, as the steepness of the waves increases, the initial rate of damping Embedded Image increases further from the linear theory predictions: the CWEC carpet harnesses energy more efficiently from steeper waves.

Figure 5.

Effect of nonlinear interactions on the initial rate of energy absorption by the CWEC carpet. (a) Linear and nonlinear results for an initial wave of steepness ϵs=0.3. The numerical solution is converged for M=4 and predicts approximately 30 per cent higher damping rate than linear theory predictions. (b) Initial rate of energy damping for a monochromatic wave train with initial steepnesses of ϵs=0.001,0.1,0.2,0.3,0.4. Other parameters are γ =0.9 and μ =1. For all simulations M=5, N=128 and T/δt=512. (a,b) Solid line, direct simulation (approximate); dashed-dotted line, direct simulation (exact); dashed line, theory. (Online version in colour.)

Direct simulation enables us furthermore to investigate the long time evolution of the surface and the bottom. Figure 6 compares the time history of energy of a monochromatic incident surface wave as it interacts with bottoms of ζ=0.1,0.3,0.8 and 1. The simulation starts from a surface wave of ϵs=0.3 and for μ=1. On the basis of initial damping rate predictions (figure 5a), energy decay rate has to (initially) increase as ζ increases until ζcr≈0.8, where the trend reverses. Figure 6 shows that this trend is not limited to the initial stage and continues over time, despite mode couplings and nonlinear interactions. The decay is not uniform and smooth as a result of energy modulating back and forth between higher and lower harmonics. This is more highlighted for the case of ζ=0.1, where the wave steepness stays high for a longer time. As expected and already noted (cf. figures 4 and 5), for smaller values of ζ, our theoretical predictions (3.4) as well as numerically computed Embedded Image and Embedded Image predict (on average) the same damping rate. For higher values of ζ (e.g. ζ=1.5 in figure 6), neglecting the mode coupling underpredicts the overall energy decay. Note that for a moderate damping coefficient of ζ=0.3, more than 80 per cent of the energy of the incident wave is absorbed in a time of about two periods or equivalently in a distance of less than two wavelengths (cf. (3.5)).

Figure 6.

Time evolution of energy of a monochromatic surface wave for different damping coefficient of the bottom ζ. Fixed variables are ϵs=0.3, γ=0.9 and μ=1, and numerical simulation parameters are M=5, N=128 and T/δt=512. Solid line, direct simulation (approximate); dashed-dotted line, direct simulation (exact); dashed line, theory. (Online version in colour.)

(b) Wave spectrum

Real ocean surface is often composed of a multitude of superimposed waves: the so-called ocean spectrum. The spectrum is under constant evolution via (nonlinear) self-interaction and through interacting with its surrounding environment. These continuous interactions and changes in the form can be well captured by the direct simulation scheme of the high-order spectral method. Here, we study the energy absorption from a (relatively) steep spectrum in a (relatively) shallow water as it proceeds and interacts with our CWEC carpet.

For the configuration and to stick to a realistic near-shore spectrum, we adopt one of the recent measurements of Elgar & Raubenheimer (2008), 2 km off the coast of Louisiana. We feed this spectrum as the incident spectrum (dashed-dotted line in figure 7) to a 2 km CWEC carpeted domain and measure the spectrum at each 0.5 km from the incident point. We assume that the water depth is constant and is equal to 3 m everywhere.

Figure 7.

Energy absorption from a spectrum of waves. The initial spectrum (dashed-dotted line) is adopted from measurements of Elgar & Raubenheimer (2008) (their fig. 3b) on the Louisiana continental shelf. (a,b) Linear simulation and (c,d) nonlinear simulation (M=4, for which numerical simulation has converged) for, respectively, γ=0.10 and γ=0.03. The bottom is assumed flat, h=3 m and the spectrum is plotted at distances x=0.5,1,1.5 and 2 km from the initial spectrum (x=0, dashed-dotted line). Nonlinearity flattens the spectrum and significantly contributes to the overall energy absorption. Direct simulation parameters are N=1024 and T/δt=64. (Online version in colour.)

If nonlinearities are not taken into account (figure 7a,b), the energy of each wave is independently damped, i.e. no energy transfer across the spectrum occurs. This is compared for a fixed ζ=0.5 and γ=0.10 (figure 7a) and γ=0.03, i.e. stiffer bottom (figure 7b). The more compliant bottom extracts more energy.

The initial spectrum of figure 7a,b upon evolution over a rigid bottom (i.e. no damping), taking into account all order nonlinear interactions, will evolve into the highest energy curves of figure 7c,d. Results of the interaction of this evolved spectrum over a CWEC carpet are plotted in figure 7c,d. Upon engagement with the CWEC carpet, the spectrum attenuates (almost) uniformly in all frequencies. Figure 7d clearly shows that owing to nonlinearities, higher energy is extracted from the spectrum (cf. figure 7b). Figure 8 shows a snapshot of surface and bottom elevations after a steady state is reached (corresponding to figure 7c). The nonlinear spectrum is fed from the left-hand side, and it interacts with the compliant bottom over the range 0<x<2 km. The interaction is stronger at the beginning of the patch (note the amplitude of bottom variations between x=0−500). Clearly, the energy of the incident spectrum decreases over the patch and asymptotically approaches zero as the interaction distance increases. All the energy is absorbed by the CWEC carpet. For the specific spectrum studied here, the energy production is of the order of approximately 20 kW m−1 of the coast line, which is achieved, as an example in figure 7c, over a distance of 2 km. Note that for the discussed example, the values of γ are not chosen for the optimum damping through which the spectrum will die much faster. For instance, if we choose γ=0.9, then the 20 kW m−1 will be extracted in less than 100 m (compared with 2 km for γ=0.1) along the direction of wave propagation.

Figure 8.

A snapshot of water surface and the bottom elevation for the spectrum studied in §5b. The snapshot corresponds to the nonlinear case of figure 7c. (Online version in colour.)

While engineering aspects of CWEC is beyond the scope of this paper, we would like to briefly comment that γ=0.9 corresponds to a springiness of k*≈11 kN m−3 (note that k* is defined as per unit area, cf. (2.1)). That is, for every square metre of the carpet, a total springiness of 11 kN m−1 is needed. This is within the range of available industrial coil springs. For instance, a typical automobile coil spring has a stiffness coefficient of about 25 kN m−1. A damping ratio of ζ=0.8 corresponds to a damping coefficient of b≈5.5 kNs m−3 in 5 m water deep. A typical automobile shock absorber coefficient is about 2.5 kNs m−1 (Dixon 2007). A typical wave energy on high-energy parts of the west coast of Europe is approximately 90 kW m−1. Noting that a typical swell wavelength may be of the order of 100 m (μ≈0.3 if h=5 m), then it can be shown that the distance required to absorb 50 per cent of the incident wave energy is approximately 7 m, which results in an absorption rate of approximately 6.5 kW m−2 of the carpeted sea (efficiencies are not included). This value is approximately 2.5 kW m−2 for the coasts of the United States.

6. Conclusion

Inspired by the strong attenuation of ocean surface waves by muddy seafloors, we have proposed that a mud-resembling synthetic compliant seabed carpet, composed of linear springs and generators, can be used as an efficient wave energy absorption device. Specifically, we consider and analyse a linear and vertically acting viscoelastic carpet (CWEC) placed on the ocean floor and under the action of (nonlinear) overpassing ocean surface waves. We showed that the coupled system of gravity waves and the CWEC admits two modes of propagating waves: the surface mode and the bottom mode. The major difference between the two modes is that the rate of energy decay of a surface-mode wave is higher for longer waves, whereas for a bottom mode, shorter waves are damped faster. To address the problem of energy transfer between these two modes, which is necessary for an accurate estimation of power absorption, and to investigate the effect of surface-wave nonlinearities on the performance of the CWEC, a direct simulation technique of a high-order spectral method was formulated and cross validated against analytical results and exact numerical solutions. We showed that for strong dampings, a considerable portion of energy flows from the surface mode to the bottom mode, resulting in a different overall rate of energy decay than predicted if the mode coupling is ignored. We also showed that for steep waves, nonlinear interactions result in a significantly higher rate of energy absorption and therefore a higher performance of the CWEC in harvesting ocean wave energy.

The presented CWEC can absorb the entire energy content of incident waves (narrow or broadband), hence it has a theoretical efficiency of unity and a broad bandwidth of high performance. The idea presented here can essentially incorporate any mud model, and its performance under different models/assumptions may be worth further investigation (e.g. frequency-dependent damping and stiffness coefficients, e.g. Jiang & Mehta (1995)). Extension of the numerical scheme to the general three-dimensional case is straightforward (e.g. Alam et al. 2010). In fact, the performance of a patch of such a viscoelastic carpet in extracting energy from an incident wave train, the patch's optimal shape and its diffraction pattern are interesting and worth studying. Finally, we would like to note that the characteristics of the CWEC discussed here (such as strong short-wave damping and bifurcation in the response) may have implications in the inverse problem of wave damping by the bottom mud or by the ice-covered boundaries.


I would like to thank Y. Liu, P. Hassanzadeh, L.-A. Couston, A. Immas and Y. Liang for careful reading of the manuscript and valuable comments. The support from the American Bureau of Shipping is greatly acknowledged.


  • 1 Equivalently, if μ is given and fixed, the expression (2.9) can be re-written for the critical damping ζcr, Embedded Image

  • 2 Note that the amplitude associated with a specific frequency output from a spectral code has energy of both modes (that have the same frequency) embedded in it.

  • Received March 27, 2012.
  • Accepted May 14, 2012.


View Abstract