CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Wave Action Equation Solver for OpenFOAM

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   August 29, 2013, 14:55
Default Wave Action Equation Solver for OpenFOAM
  #1
New Member
 
Tom Chyczewski
Join Date: Mar 2009
Location: Bethpage, New York, USA
Posts: 15
Rep Power: 8
chyczewski is on a distinguished road
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:

\frac{\partial N}{\partial t}+\frac{\partial C_xN}{\partial x}+\frac{\partial C_yN}{\partial y}+\frac{\partial C_kN}{\partial k}+\frac{\partial C_\theta N}{\partial \theta} = S

where N is the wave action which is a function of geographical space (x,y), wavenumber magnitude k and wavenumber direction \theta.

(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:

\frac{N^{n+\frac{1}{2}}-N^{n}}{\Delta t}+\frac{\partial C_kN}{\partial k}+\frac{\partial C_\theta N}{\partial \theta}=0

\frac{N^{n+1}-N^{n+\frac{1}{2}}}{\Delta t}+\frac{\partial C_xN}{\partial x}+\frac{\partial C_yN}{\partial y}=S

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 \tilde N(x,y,z) = N(x,y,k_i,\theta_j), where z is a third dimension added to fit the OpenFOAM framework. So, actually, \tilde N is solved in the second set for all (i,j).

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. (i,j)) 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
chyczewski is offline   Reply With Quote

Old   August 30, 2013, 00:06
Default
  #2
lin
Senior Member
 
Hua Zen
Join Date: Mar 2009
Posts: 114
Rep Power: 8
lin is on a distinguished road
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.
lin is offline   Reply With Quote

Old   September 3, 2013, 07:21
Default
  #3
egp
Senior Member
 
egp's Avatar
 
Eric Paterson
Join Date: Mar 2009
Location: Blacksburg, VA
Posts: 197
Blog Entries: 1
Rep Power: 9
egp is on a distinguished road
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
egp is offline   Reply With Quote

Old   September 3, 2013, 14:41
Default
  #4
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,758
Rep Power: 21
hjasak will become famous soon enough
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
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   September 3, 2013, 16:11
Default
  #5
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Rotterdam, The Netherlands
Posts: 1,595
Rep Power: 24
ngj will become famous soon enoughngj will become famous soon enough
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
__________________
Please note that I do not use the Friend-feature, so do not be offended, if I do not accept a request.
ngj is online now   Reply With Quote

Old   September 3, 2013, 16:38
Default
  #6
egp
Senior Member
 
egp's Avatar
 
Eric Paterson
Join Date: Mar 2009
Location: Blacksburg, VA
Posts: 197
Blog Entries: 1
Rep Power: 9
egp is on a distinguished road
Niels, nice summary!

Quote:
Originally Posted by ngj View Post
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
egp is offline   Reply With Quote

Old   September 4, 2013, 03:19
Default
  #7
New Member
 
Alfons Smale
Join Date: Jun 2009
Location: Delft, Netherlands
Posts: 2
Rep Power: 0
asmale is on a distinguished road

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.

asmale is offline   Reply With Quote

Old   September 4, 2013, 11:06
Default
  #8
New Member
 
Tom Chyczewski
Join Date: Mar 2009
Location: Bethpage, New York, USA
Posts: 15
Rep Power: 8
chyczewski is on a distinguished road
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.
chyczewski is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Transport Equation for a scalar - Compressible Solver alessio.nz OpenFOAM Programming & Development 5 April 11, 2013 09:04
Problem with a leapfrog scheme for wave equation Shiranui Main CFD Forum 0 June 22, 2010 09:19
Wave equation, wave velocity components mztcu CFX 1 May 4, 2010 02:14
Wave equation for hfss in fluent jpcfd FLUENT 0 March 16, 2010 11:45
Assistance with one-way wave equation Frank Main CFD Forum 5 March 14, 2007 11:56


All times are GMT -4. The time now is 12:00.