Sorry for bothering you, I know this has been up for discussion before. I'm using quite an old version of waves2Foam and it doesn't do well with point motions within the relaxation zones (solver just freezes on relaxing.correct) I believe you said this was fixed in newer versions but I'm reluctant to upgrade since I have modified stuff and I don't like to fix things that are working in every other aspect. Do you remember what part of the code is troubled by point/cell motions in the relaxation zones?

Alternatively, I think it would be ideal to actually keep point motion at zero in the relaxation zones to improve accuracy. I have experimented a bit with custom motionDiffusivity classes to achieve this with moderate success. I can greatly reduce the point motion inside the relaxation zones but some motion always "leaks" in for some reason. Is this something you have experimented with as well?

What I fixed was related to topological changes, i.e. a change in the number of computational cells originating from e.g. mesh refinement. I do not believe that I have changed anything related to point motion in the relaxation zones; actually I always try to eliminate that, though since I specify motions at the bed and only have vertical movement, it is fairly easy for me.

I am interested in knowing the reason for the freeze of the relaxation.correct(). Do your cells become severely distorted/non-convex during the mesh motion?

What are your changes? Are they of common interest? Once the SVN is running again we could discuss, if you should contribute in order to bring you (latest version) and all the rest of us (new functionality) a step ahead.

I do not believe you can avoid mesh motion in the relaxation zone through the mesh diffusivity, I would rather suggest that you simply remove those points/cells from the mesh motion matrix.

That is probably why I couldn't find anything when I searched the forum, I must have remembered wrongly.

I have unsuccessfully tried to narrow down the reason for freezing. The cells are not distorting much. The distortion is hardly visible in the relaxation zones since they are far away from the moving object. Nevertheless as soon as some distortion reaches the region, the solver freezes (no FPe or anything, it just stops responding.) If I comment out relaxing.correct, it all works fine and equally, if I turn off the mesh motion. I might try to pick out the points in the relaxation region from the motion matrix as you suggest. It should work but I'm not sure how to do it practically.

As for my changes: not really of common interest I think, they are quite old. Just things that I have stuck with since the first release of waves2Foam. For example, one of them is implementing the "air velocity" which you already released ages ago. I should probably just update and it would out-date most of them. I just don't like messing with things that work

The only ones that are unique for me are relating to linking the waves2Foam libraries with a custom self propulsion library and solver. Whenever that is working properly I will definitely release it but that is not in the near future I fear.

Then it is "quite" easy, even though there will be some tweaking of the source code, both OF and waves2Foam.

What you will need to do is to make a new class, which to begin with is merely a copy of your mesh motion method. Then, in some way or the other find all cells in all relaxation zones from inside the mesh motion method [1].

Somewhere in the mesh motion technique, there will be a point, where you construct the matrix system. After this step and before solving for the mesh motion, you prescribe a fixed displacement of (0 0 0) in all the cells in all relaxation zones. You do this by invoking the command:

Code:

<vectorEqnName>.setValues( cellList, vectorField );

This will then enforce the motion to be zero inside the relaxation zones and the inner limit to the relaxation zones will begin to work as pseudo boundaries.

Kind regards

Niels

[1]: This is not straight forward, since relaxationZone is not registered in the database, but inspired by your question I will see, if I can find the time to make in happen.

That is how I accessed it in my custom motion diffusion model. Of course it would be easier if they were just saved as point/cell/face zones or something. Speaking of that, I found this bit in the displacementLaplacian motion solver

I will have a go and see if actually just giving the points in the relaxation zones as pointZones to this will solve it. Otherwise I will try something a bit more advanced, I'll let you know how it goes.

All you will then need is to put a small method in the relaxation zone, which gathers (unique) cell/point indices from all of the relaxation zones. This should be fairly easy and you could benefit from using the

Code:

labelHashSet

class. Please post if it works and I will put the above changes into the release.

Ok, so here is an easy fix for anyone using dynamic meshing and want the mesh to be frozen in the relaxation regions:

1. Create a pointSet with all the points in the relaxation regions.
2. Convert to a pointZone using setsToZones
3. Add the following line in the dynamicMeshDict

Code:

frozenPointsZone <name of your pointZone>;

This is implemented in the displacementLaplacian and componentLaplacian motion solvers (as of 1.7.1) but It is pretty easy to put in a similar "hands on approach" as the one I showed above in any solver.

Take heed though, this basically just ignores any point displacement added to these points by the motion solver and sets it to zero. If you have a lot of motion close to the relaxation regions this might create very skewed cells on the border to these. If your motion "source" is relatively far away it should be fine. A better fix would probably be to follow Niels' suggestion above to remove these points from the motion matrix which should ensure that the motion gets diffused properly.

//Björn

PS. If you look at the point displacement in Paraview, it will be non zero in the area you specified. This is because the motion gets set to zero manually so the calculated pointDisplacement gets ignored but still printed as a field. The points will be stationary though if you look at them.

I have a question for you. I am having problems with my simulation because of pressure instabilities.

I am trying to model a ship traveling in the ocean through the waves with constant speed.

At first, everything seems to be alright in the simulation. But after some point, the oscillations of the ship in the direction of gravity seems like getting higher and higher (the vehicle sinks more and jumps more on the water as the simulation progress further). You can find a video of my simulation below:

I checked everything (mass of the vehicle, density of the fluid-air-vehicle, wave parameters, etc.) but nothing seems wrong.

I tried to run the simulation without any water movement or waves. But, again the oscillations of the ship gets higher and higher because of the unstable pressure variations under the ship until the analysis gets blown up.

The only thing I know is that the system was solved for total pressure in 1.7, though in the following version this strategy was abandoned, because of issues with obtaining correct forces.

I have never tried to work with floating objects, so I am virtually clueless. However, are you sure that the 6DOF solver, which must be based on some ODE-system, is free of instabilities. It could remind me of some explicit Euler time integration?

I have made some tweaks to the build scripts to allow the top level Allwmake script set the location of the GSL includes and libs, and the location to install the waves2foam libs and applications. I routinely have to change both of these when I build a new version for our cluster.

While I was at it, I updated the solver220 for the head of the OF2.2.x git trunk and included the waveDyMFoam while I was at it. The top level Allwmake also was modified to build waveDyMFoam when building for 220.

I plan on having Kevin create 2d and 3d tutorials for a floating object when I can free him up from his current projects.

I have attached a patch file for these changes. Let me know if you are willing to accept these proposed changes.

Thanks a lot; the already posted and coming contributions are greatly appreciated.

I am currently without my "waves2Foam" computer and also without access to the SVN (still), though I will make sure that these contributions are added in the next release.
The date for the next release, however, is a bit uncertain, as it depends on when the shipping company finds its way to my new home

As I still struggle with the breaking wave problem, I have two questions about wave2Foam:

I assume that other solver options can also be used besides GAMG which is currently selected in wave2Foam. The damBreak tutorial uses PCG solver. In your experience what advantage does GAMG solver has? In the simulations in your thesis, did you try different discretization schemes and/or pimple control parameters?

I noticed that const (zero) wind vector is specified in the source file. If wants to examine the the wind effects on waves, it is relatively straightfoward to implement the desired wind velocity profile. But as you mentioned, by including the air-phase (no wind velocity) in the solver, the air flow driven by wave motion already may make the solution unstable? So what would you expect if non-zero wind velocity is specified?

For a large number of cells GAMG will outperform PCG, and furthermore I have once been told that PCG is not well suited for parallel computing (only statement from my part, as I do not know for certain). For small meshes you will find that PCG is faster than GAMG, however, then the speed hardly matters.

With respect to wind, then I really do not know anything about that. I can only advice you to test it. Only reason for the addition of the wind vector was due to requests from people doing moving ships in a translated coordinate system.

Kind regards

Niels

June 2, 2013, 07:11

Effect of temporal schemes in waveFlume tutorial and steepness limitations

First of all thank you for such a wonderful (and easy to use!) wave generation and absorption toolbox, it has made life so easy for us fellow coastal engineers! I really appreciate your continuous effort in helping everyone out through this forum and for painstakingly answering each and everyone's questions.

I have been using your toolbox for about a couple of months now and my ultimate research goal is to use it for fluid structure interaction, particularly to study wave dissipation and resulting turbulence generation through coastal vegetation.

Since I am interested in simulations based on laboratory scale flumes I started with the waveFlume tutorial and was doing some sensitivity analysis with a bunch of different parameters (resolution, PCG/GAMG pressure solvers, convection and temporal schemes, alpha schemes, turbulence models, Co number, length of relaxation zone, type of relaxation i.e, exponential/free poly/3rd order poly, wave steepness, Ursell number, etc) and happened to stumble upon a phenomenon where there is an un-physical growth of wave height as the waves propagate through the domain in the cases when any of the OF temporal scheme other than the Euler is used (Fig. 1 in attachment). It appears that the extra diffusion due to its first order nature balances whatever is causing the growth. Also noticeable is the standing wave signature when Crank Nicholson/backward is used with reducing 'phi' values for CN yielding better and better results due to progression towards pure Euler. Pure CN (phi = 1) appears to be unstable as the simulation blows up a little after 4.5s. Regarding the sponge layer relaxation technique the free polynomial performs marginally better in reducing the standing wave problem (with increasing exponents till about 6 giving good results). Why this bothers me is that a first order in time will be too diffusive for LES models and since my final goal is to use LES to study fluid structure interactions I am wandering whether waveFoam will be accurate for LES applications. I did not find any second order temporal schemes in any of the available waveFoam tutorials and was wandering if you could provide some insights into this wave height growth or point out to me if I am missing something. I am also not sure if this is a VOF-MULES thing as I did not find any of the interFoam tutorials also using anything other than a Euler. May be OpenFOAM VOF experts working in fields other than waves can also comment on this or provide me with a reference which gives an idea why MULES should be unstable with a second order in time scheme.

I also tested a bunch of different wave cases (all using the streamFunction theory) with the goal of understanding any limitations as to steepness and Ursell number. I found that steepness of > 0.05 causes instability in the free surface and the waves lose their monochromatic nature (Fig 2, all cases have T=0.75s, h=0.3m with at least 6 cells/wave ht) and this also corroborates the findings of Afshar (2010)*. If I understand correctly in the thesis one of the causes was attributed to the un-physical high velocities in the air phase just at the interface. I used the technique as in Patterson (2008)** to take convection out of the air phase but it appears to have no effect on the weird behaviour (only benefit is the time stepping and somewhat reduced computation time). Inspection of the spectra suggest (Figs. 4 & 6) that the interaction among the higher harmonics for >0.05 steepness waves increase rapidly as waves exit the inlet sponge and progress towards the sponge layer. The Ursell number seems to have almost no effect (provided the boundary inlet wave is of sufficiently high streamFunction order) so long the steepness is less than 0.05. I wanted to know whether you have faced similar problems during your runs and if yes if you found a solution to it.

The sample case files are attached (Btw I am using waveFoam with OF 2.1.x). If you need the full dataset and all the files please let me know I will email you the link.

Thank you once again for such a wonderful tool and wishing you all the very best for the continued development of the open-source CFD community.

A warm welcome to the forum. I will try to address your two main questions namely (i) higher order time stepping and (ii) steepness of the waves individually. However, first let me thank you for the analyses that you have performed, as it does reflect a lot of the issues in terms of VOF and waves.

ad. i: If you look into my PhD-thesis ([1]) I considered the Beji and Battjes experiment over a submerged bar; a test case, which has recently been added due to a kind contribution from Jorge Cadelho. Here I considered the effect of using either implicit Euler or Crank-Nicholson time integration and I found a considerable improvement in the predictions by using Crank-Nicholson.
But, but and but! The CN time integration scheme was not the native scheme in OF, since I have never quite understood how the Crank-Nicholson was defined over three (3) time steps, since the ddt0(rho,U) terms are stored and used for the time integration. The solution to be revealed below. I has occurred to me at this instance (incorrectly?) that the second order time integration schemes in OF assumes constant time steps, however, most often variable time steps are used for VOF-simulations. At least I derived the third order explicit Adam-Bashforth time integration for variable time steps once, and the coefficients became exceptionally ugly, and I do not recall seeing such ugliness in the OF-source. If I used the AB3 with variable time steps under the assumption that I could use the constant time step formulation did give mass conservation errors.
Sorry, rambling a bit. The solution was to follow the Ferziger and Peric book and simply define the CN directly based on two time instances. This is done by correcting the implicit terms, e.g. convection, from

Code:

fmv::div(rhoPhi,U)

to

Code:

0.5 * fvm::div(rhoPhi,U)
+ 0.5 * fvc::div(rhoPhi,U)

If I remember correctly this also gave rise to increased stability compared to the native Crank-Nicholson. I have not been using this formulation for wave breaking, simply because it is not robust enough.

ad. ii: Please read our journal article, which accompanied the release of waves2Foam. Here, you will see a nice example on the effect of cell aspect ratio on the stability of the interface. Also, as an interesting analysis see Ed Ransley's et al. poster ([2]), where they found that it is much faster to reach convergence, if you use an aspect ratio of 1 around the free surface.

I hope that this has answered some of your questions and that you will be able to progress from here.

Hi Niels, waves2Foam continues to be an excellent tool for my studies.

I'm currently trying to create a new wave theory to generate focused waves with a second order correction. I'm not familiar with C++ and have done little programming in the past. I've come across the newWaveTheory files in the waveTheories directory, can these be used to generate a new theory easily or should I make a copy of an existing theory, modify it and recompile. Would I have to recompile the whole of waves2foam? Is there any documentation or tutorials on how to generate a new wave theory?

No, the newWaveTheory.C is a runTime selection method, which allows you to specify the wave theories in a dictionary by a type name. Furthermore, at the moment there are no tutorial in creating a new wave theory, however, please fill the gaps on the Wiki, if you are so kind.

For a new wave theory merely choose any of the wave theories and fill out the functions with the correct algebraic functions, e.g. stokesFirst.[C,H], and remember to change to class name throughout. Since you are doing second order focused waves, I suspect that you will find the bichromaticSecond.[C,H] interesting, simply because the transfer functions for the super- and subharmonics are implemented (for two waves). However, if you consider the second order irregular wave description by Fuhrman and Madsen (or is it Madsen and Fuhrman?), then the transfer functions are the same.