Wave Action Equation Solver for OpenFOAM
 

I need to develop a solver for the wave action equation using the OpenFOAM library and am very interested in any insights veteran users might have on the best way to get this done.

The wave action equation is of the following form:



where is the wave action which is a function of geographical space , wavenumber magnitude and wavenumber direction .

(The wave action is related to the ocean wave height spectra and this equation governs the evolution of this spectra in the presence of nonuniform currents and winds).

My first question is, does anyone know of any implementations of a solver for an equation like this in OpenFOAM?

My suspicion is that the answer to my first question is 'no' (though I hope I'm wrong). So, I've been thinking of how I will do this and this is what I've come up with:

This equation can be solved using a fractional step method:





The first step is local in geographic space, including only derivatives in wavespace, which is discretized in a simple structured manner. So the first step can be completed sequentially for each point in geographical space.

Then the second step can be solved sequentially in wavenumber space using the features of the OpenFOAM library to solve the equation that contains only geographical space terms. This would require, for each point in wavenumber space, copying , where is a third dimension added to fit the OpenFOAM framework. So, actually, is solved in the second set for all .

Assuming this approach makes sense, can someone outline how I would create a volume scalar field array with two added dimensions to represent the wavenumber space so that I can simply specify a point in wavenumber space (i.e. ) to be copied to a volume scalar field that can be processed in OpenFOAM? Once that is accomplished I think the only other thing to contend with is how to apply the boundary conditions to the second step equation, which are a function of the position in wavenumber space. Though it doesn't seem to daunting, I haven't completely thought it through yet.

Any insights or suggestions would be greatly appreciated.

Tom

From my limited knowledge of openfoam and ww3, in my opinion, unless you are very familiar with openfoam, it may be better to use the programming language that you are familiar with, such as C, Fortran, Maltab etc.

Hi Tom,

This sounds like an interesting problem! Ignoring the wave space issue for a moment, I think you need to first cast the problem as a 2D surface using either the finiteArea machinery of 1.6-ext, or the thin film model in 2.2.x (see reactingParcelFilmFoam).

The wave space issue is a bit more complex and I don't have a good recommendation for an elegant solution. Brute force for a fixed number of wave-space bins would be doable, but ugly. What is the typical wave-space discretization in other codes such as WaveWatch III http://polar.ncep.noaa.gov/waves/wavewatch/

Good luck,

Eric

Hi Tom,

Sounds interesting - I have done waves, but only in 3-D. Some clarifications if you don't mind:
- you are solving on a 2-D surface, right? Is it curved or flat?

For the wavenumber space, are we talking 1-D, where each wavenumber carries its own Ck, Ctheta pair?

- now mane Ck and Ctheta fields will you have? Will it just be N of each, or NxN as you seem to imply? I will call it NxM for the code below...

If N, you need

PtrList<volScalarField> ks(N);

PtrList<volScalarField> thetas(N);

You will access them as ks[i] and thetas[i].

If NxM, make a PtrList<PtrList<volScalarField> > bigOne(N);

forAll (bigOne, i)
{
bigOne.set(i, new PtrList<volScalarField>(M));
}

and then fill the lot. You will access them as bigOne[i][j];

I have done quite a bit of this sort of thing (OpenFOAM) :) and would be happy to help. I guess you will find my details on the web; if not, please send me a message and we can catch up.

Yours,

Hrv

Hello all,

I have been following this thread out of pure curiosity, so allow me to follow up on the equation.

The main equation above it indeed in a four dimensional space, where two of the directions are spatial coordinates and the two remaining "directions" is the wave number (or frequency) and the wave propagation direction.

Using the above suggested decomposition it is easy to understand the physical system. In each spatial point (x_i, y_i) you have a wave spectrum, where the energy is distributed on a range of frequencies/wave numbers and a range of directions (optimally, the domain will have to be cyclic along the angular axis). The spectrum is represented by the wave action N. This is also why this type of model is called "spectral wave model".

Therefore, the solution to the equation is not a deterministic description of the behaviour of the free surface as e.g. obtained with interFoam/waveFoam, but rather a statistical representation in a given point. E.g. in a simple 1D spatial domain, you would have the JONSWAP spectrum (http://www.google.nl/url?sa=i&rct=j&...78322546123846), since the wave can only propagate in one direction (or rather positive and negative directions).

The propagation speeds (C) for x and y are simple propagation speeds typically computed from the linear dispersion relation including the effect of wave-current interaction. The C_theta is a propagation speed related to refraction, i.e. advection of wave action from one propagation direction to another due to changes in the bathymetry. I am not quite sure with respect to the last velocity.

The source terms on the right hand side include effects from bottom friction, wind drag, dissipation due to wave breaking and non-linear transfer functions from one frequency to another. The latter involves either three or four frequencies/wave numbers (triad and quad interactions of the wave components).

Finally, to answer the question of the lay-out of the computational domain. The well know commercial and open-source tools available (Wave Watch III, SWAN and MIKE SW) all solve these equations on flat meshes, i.e. local conditions around a harbour, or on spherical meshes for the global prediction of wave conditions.

This type of models are then typically coupled to large scale hydrodynamic models/oceanographics models, such that predictions of storm surge, wave induced circulations, etc can be simulated. All of these parameters have importance for a wide variety of topics: swimmers safely, coastal protection, disaster management (should we or should we not evacuate this area), navigation, etc.

Kind regards

Niels

Niels, nice summary!

Hi all,
In addition to the summary of Niels, I would like to add the following:

Presently I’m working on solver that solves the above mentioned equation, neglecting the variation in frequency space from the above shown wave action equation (so that results in the dimensions time, x, y and direction), similar to Holthuijsen et al (1989). The implementation of that wave action equation was possible using readily available components from OpenFOAM and assigning the directional space to the z-axis. Results are unfortunately not yet ready for publication, that will take some more time, but it seems to be working.

The next step is to extend the above mentioned solver with the frequency space. For this, a “trick” is needed to allow for the extra dimension within OpenFOAM. One way is, as above mentioned by Mr. Jasak, to create additional fields for the different frequencies to be considered (or even extra fields for each frequency and direction component if you use the FaMesh approach from Mr. Patterson). This would subsequently also require a connectivety table that can be used for solving energy transport from one frequency to another (either by means of the propagation or by means of source terms). There are also other “work arounds” to obtain the same result: either by splitting up the action equation (as stated in the first post), using multiple meshes, or using a VectorN approach.

I will be working on the next step (including the frequency space) in a few months time, using approaches from the spectral wave models WaveWatch and SWAN. In this next step I will use specific schemes for discretisation of the wave action equation that we are presently investigating as part of our long term research and development of SWAN.

I’ll try and keep you posted on the progress/results.

Regards,
Alfons Smale


Holthuijsen, L.H., Booij, N., Herbers, T.H.C., 1989. A prediction model for stationary short-crested waves in shallow water with ambient currents. Coast. Eng. 13, 23–54.

Thank you all for your interest in this problem and the suggestions.

Eric, Thanks for the suggestions on finiteArea and the thin film model. I'll look into them. As far as wave-space discretization schemes, here's a summary of what I've seen: WaveWatch III has both a simple first order upwind scheme as well as a third order / limited scheme available (which they call ULTIMATE QUICKEST). Explicit methods are used to solve the equation. SWAN uses a hybrid central difference - upwind scheme. In this case the equation is solved with an implicit method. (In both WaveWatch III and SWAN the same methods are applied to both wavenumber magnitude and direction space.) There are also a couple of unstructured implementations of SWAN; one using a finite element method (Hsu, Ou and Liau, 'Hindcasting nearshore wind waves using a FEM code for SWAN', Coastal Engineering 52 (2005) 177-195) and another using a finite volume method (Qi, Chen, Beardsley, Perrie, Cowles and Lai, 'An unstructured-grid-finite-volume surface wave model (FVCOM-SWAVE): Implementation, validations and applications', Ocean Modelling 28 (2009) 153-166). Both use the Flux Corrected Transport scheme in the wavenumber magnitude space and a central differencing / Crank Nicholson scheme in wavenumber direction space. In all cases, a logarithmic distribution discretization is used for the wavenumber magnitude and a uniform distribution in wavenumber direction. My initial impression is in line with yours; a brute force method should be doable. My plan is to start with that and work towards something more elegant.

Hrvoje, yes, it is a 2-D surface. Both WaveWatch III and SWAN solve a spherical coordinate form of the equations allowing for curved surfaces. My plan is to solve the Cartesian form and just consider flat surfaces.

"For the wavenumber space, are we talking 1-D, where each wavenumber carries its own Ck, Ctheta pair?"
Yes, each point in wave-space (i.e., each wavenumber magnitude and direction pair) has a distinct Ck, Ctheta pair. However, Ck and Ctheta are also functions of your position in geographical space since they are a function of the surface current speeds and ocean depth. This is also true for C_x and C_y. You can store k and theta with 1-D arrays since the discretization does not vary with geographical space (as you suggest in the code contained in your response). Also note that all of the C's are only functions of the surface currents, ocean depth and wave-space position. So they need only be calculated up front and stored for use through the solution procedure. I believe I need to store these things with 3-D arrays; one for geographical space and two for wave-space.

For ks and thetas your suggestion works:

PtrList<volScalarField> ks(N);

PtrList<volScalarField> thetas(M);

But, for the C's and the wave action I need C(I,N,M) where I identifies the position in geographical space. Is there a way to generate an array that can be accessed bigOne[i][n][m], where i is the mesh index, n is the wavenumber magnitude index and m is the wavenumber direction index? If so, for each position in wave-space [n][m], I can copy bigone[i][n][m] into smallOne[i] to solve the second step equations (as outlined in my first post).

Your offer to help is greatly appreciated. I will in all likelihood take you up on it as I move forward.

Niels and Alfons, thanks for your posts. I am interested in following your progress Alfons. I didn't know there was a version of the equations that included just the wave direction and not the frequency.