[markc]

Hello All,

The BC is great but: not perfect yet. I found a bug and some other strange things.

1.
A clear bug:
interfaceWeight.H contains the following line:
>>>
const label patchIDInlet = mesh.boundaryMesh().findPatchID("inlet");
<<<
Having a case which does not have a boundary called "inlet", the BC still continuous without even any suspective behaviour. the label patchIDInlet gets the value -1 and everything runs smooth. Strange? Some lines further in this file the label is used:
>>>
label patchID = mesh.boundaryMesh().whichPatch(cProp[cJ]);
if (patchID == patchIDInlet)
<<<
I have some simple 2 D cases where it does not matter what the name of the inlet boundary is.

2.
Strange behaviour:
As told I have some 2D cases which run fine. Meshes are generated with blockMesh and bodies are added using snappyHexMesh. However I also have a larger case (appr 80k cells), which is 3D, generated at the same way (blockMesh, sHM), everything looks the same, but this case fails. The last thing which is writen in the log file is:
>>>
This routine has not yet been implemented
<<<
which is writen by evaluateFace.H. A few lines further it crashes due to division by zero (sigFpe).

I played around with wave settings in environmentalProperties (e.g. sealevel, waterDepth etc), without succes.
Well, it is clear that the specific part is not implemented yet, but it seems strange to me that other cases which are almost the same do run fine. I am not sure if this has anything to do with the "bug" as explained under (1).
Furthermore I played around with the X position of the inlet boundary. This does not have any influence, it can also be placed at a negative position.

Any help is appreciated,

Brgds,

Mark

[ngj]

Hi Mark

Yes, that is true, I did hardcode that some time ago. I am actually working on another branch of this BC now, so it is some time ago I looked into the code.

I believe you should be able to replace "inlet" by this->patch().name().

Regarding 2: I wouldn't call it a bug, as it is deliberate from my side, as I didn't needed it. The reason for successful run in 2D is due to the fact that dy/dz=0 for the "horizontal" edged in 2D, whereas they might be modified by snappyHex and changed to a non-zero gradient.
This "bug" would take some time to fix, and as I am working with another branch, it does not have sufficient priority at the moment. But if you sit down and go through all of the geometrical considerations leading to that, you should be able to modify the code.

Could you elaborate of the X-position of the inlet boundary? I do not understand what you mean? If you mean, that it does not change the way the waves are generated, then yes, as X is not part of the formulation for gamma, U, pd (assumed X=0).

I do know it is not perfect, but I wanted to contribute with something, as I didn't have success with those versions floating around at the forum prior to my posting.

I hope you will be able to succeed with your work.

Niels

[markc]

Hello Niels,

First of all don't get me wrong. I am very enthousiastic about your work.
Next: somewhere earlier in this thread I read something about X position of boundary, that's why I investigated that matter. And sure: this does not have any influence, as you explained already.

Regarding 1: I will try the "this->patch().name()
Regarding 2: Ok, that makes things clear. So it is indeed true that the BC is yet onlt for 2D cases and that the "not yet implemented" part is meant for faces adjoining non-empty boundaries. That explains a lot.
Finally: implementing code to calculate velocities at the interface gives a lot of extra work. Isn't it sufficient to simply distinghuish whether a cell centre is above or below the so-called eta value? Apparantly not, otherwise you would not take all that effort.

Brgds,

Mark

[ngj]

Hi Mark

My apologies, after submitting I regretted the tone in my post. I had been having problems for days with the wave BC with a non-consistent error. I found the bug and it is definitely also in your version:

add "evaluate()" in the end of the constructors taking "dict".

That solves it. Otherwise you would end up working with non-initialized patches, and if you are lucky, you get reasonable values, otherwise everything blows up.

@2: Well not quite, I have been running 3D cases successfully, however my mesh was purely blockMesh, where the inlet is a rectangle aligned with the y-z coordinatesystem and no non-uniform edge stretching, thus still dy/dz == 0 for the horizontal edges on the patch.
I suspect that the snappyHex modifies the points on your patch to adapt the mesh to your STL-surfaces, and even though you might have generated a blockMesh as described above, the dy/dz criterion is invalidated.

@3: Using the approach and having a wet face only soforth y-coord on the face is smaller than eta and a dry face otherwise results in a binary approach. My experience is that it generated secondary surface features such as small wave length wiggles. This is unwanted.

I hope it explains your doubts,

Brgds

Niels

[markc]

Hello Niels,

Thanks for your clarifications again. However, my next question again (sorry for that):
Adding your "bug fix" as you explained results in a runtime error on the cases which previously ran fine. The error given is:
>>>
request for surfaceScalarField phi from objectRegistry region0 failed
available objects of type surfaceScalarField are

0
(
)
<<<
Running the solver on a case which already contains a phi file from previous runs results in a similar error.

Am I doing something wrong? What I did was:
add "evaluate();" at the end of constructors which are called with dict (clear enough) in the 3 .C files (surfaceWave...C). Is this correct? Should it only be in the ...velocity...C?

Brgds,

Mark

[ngj]

Hi Mark

I went into the version you are playing around with. After a little grepping I found that the line:

const fvsPatchField<scalar> & phip = patch().lookupPatchField<surfacescalarfield,scalar >("phi");

in interfaceWeights.H is unnecessary. If you remove it, everything should be working as before.

Did it work changing "inlet" to "this->patch().name()"?

Best regards,

Niels

[markc]

Hello Niels,

1. Thanks, I will try asap to remove that line
2. Forget to mention but YES it worked!

Brgds,

Mark

[markc]

Hello Niels,

Removed the phip - part and also had to remove the next line, where U is getting read. This worked out!
But it still did not work on my other 3D case, though the specific boundary is made by blockMesh and not touched by sHM. Therefore I started with a simplification. For every cell face on the boundary I only look at the points with the largest and smallest Z value. Determining whether the cell is wet or dry is simple. When eta is in between both than we have an interface cell. Next I simply do (eta-minZ) / (maxZ-minZ) to determine some kind of interfaceRatio (comparable to your Awet/(Awet+Wdry). This approach seems to work. Actually on my 2D cases I have exactly the same answers as with the more sophisticated version but the simple version also works on my 3D case.
I am not certain about what kind of accuracy I am throughing away with this approach. I suspect that it is perfectly valid for rectangular and orthogonal faces. Do you have ideas about this?

Other question: Playing around with the BC people may have noticed reflection of waves. This can occur on the outlet but setting zeroGradient on U, pd and gamma effectively cures this for my cases so far.



However reflection (or similar kinds of unphysical secondary waves) also seem to occur in the middle of the domain, especially when the domain contains objects. Others must have seen this as well. Any experience about this issue is welcome. Things which come to my mind are:

• mesh size / mesh grading
• mesh is exactly aligned with wave travel direction
• PISO settings (?)
• discretisation schemes (?)

Any ideas here?

Brgds,

Mark

[ngj]

Hi Mark

I am glad it worked out for you! Though it is pretty weird that you 3D-cases doesn't work, as I have been running those several times without any problems.

With respect to the discussion of reflection, then of course structures in the interior will cause reflection, and one would like to avoid this reflection to come in contact with the inlet-boundary.
Further I have experienced reflection when large gradients in the computational cells occur.
However it is a relevant and interesting discussion.

Currently I am working on a completely new structure / approach to the wave-making including a library which should treat the spongelayers/relaxation zones. In that context I have a couple of questions for the use of the BC:
1. Has it been a limitation that the vertical component is the y-coordinate? I am planing to hold on to that.
2. Are you exclusively using plane patches? In terms of the relaxation zones it would make the implementation significantly easier and more cost effective, if your domain boundaries to a large extent are plane/retangular.

Best regards,

Niels

P.S. Have you succeeded in setting up the new Forum-profile, so you can receive an email for all post made in the OpenFoam context?

[markc]

Hello Niels,

1. Actually in my modified versions of your BC I swapped Y with Z as this is more intuitive for the cases we use over here. I intend to make this variable however, by looking at the gravity vector. Also the positive Z direction could be taken into account. A little extra overhead but not a big deal. Once I succeed in this I will post it here.

2. So far we only use square and orthogonal inlet boundaries.

Brgds,

Mark

[ngj]

Thanks for that, it might not be to much of an effort to make it generic. I will give it some thoughts.

Best regards,

Niels

[ngj]

Hi Mark

I reread your post, and it is clear that an implementation, which takes the direction of the gravity vector (i.e. X, Y or Z direction) into account, is the most favorable.

However, I wondered about your note on the direction of the Z-vector!?! Is it such that you sometimes work with g(0,0,+9.81), i.e. having the vertical axis pointer toward the center of the Earth?
If so, it is a good point, which I had not considered at all.

If your are having any ideas for the new implementation I am currently working on, they should be most welcome.

Best regards,

Niels

[markc]

Hello Niels,

So far I am not in need of the g vector having a positive value but I am of the opinion that it is better to make it as general as possible. After all it is a relatively simple task, compared to what you have achieved already.
I will try to attach my version here. It is generalised in terms of the wave elevation direction, looking at the gravity vector and also in the wave travel direction, which is set by an extra input in the environmentalPropertiesDict. This input is yet not really used but is intended for a future option when I want to state waves coming in from a certain direction.
Also it contains an input "refU" which is a velocity vector which adds to the wave vectors. I need this to simulate objects traveling at a certain speed against the waves.
The BC is much more basic than Niels's. It has a very simple approach for determining the interface, as explained earlier in this thread. I suspect that this approach works well on orthogonal faces but will be less accurate on (e.g. tet meshes. Any comments are welcome.
Furthermore I removed the 2nd order stokes, just to make it somewhat easier to implement and test the other new features.
My version is not meant to become a sophisticated BC, it is only meant to play around a little and test some features. So one can use it e.g. to see how the direction of gravity is implemented. Only partly tested, so be careful.

Brgds,

Mark

[ngj]

Hi Mark

Thanks for sharing your implementations.

Yes, I agree wih what you have done in your code. Just to summarize for everyone else (using Stokes first order theory as an example, developed for origo at still water level):







where is wave number vector, , the depth and a vertical coordinate related to the direction of the gravity vector as



This leads to a velocity field defined as:



The only requirement for this to work is obviously that



With the inherent inner-product operations in OF, the above seems to be quite straight forward to implement.

Have a nice day,

Niels

P.S. The above actually ensures that you are not limited to be having the gravity vector pointing in the direction of any of main coordinate axes. It is completely free to choose the direction, but it complicates the expression of the wave number vector

[ngj]

Hi all

I have been in the process of reimplementing the above mentioned wave boundary conditions into a completely new structure, which for instance make use the runTime selectable tables.

In the process of this reimplementation I realized that some of the approaches I did take was erroneous, and thus the results from the current available wave BC must be considered with some care and tested carefully.

Unfortunately I am not able to deliver the necessary patches to the current BC and the new implementation will not be made available for some time.

This message should be considered as a warning, and users are encouraged to make thorough testing before actual use.

My hope is to release the new version at some point.

Best regards,

Niels

[icingfish]

Originally Posted by ngj

Hi Charlotte

My initials have been changed, so please replace "s001581" with "ngja" and you should be able to access the files.

Best regards,

Niels

i've replaced "s001581" with "ngja", but i still can't follow the link,
would you please Niels tell me how i can get and learn from your bcs, thanks!

[ngj]

Hi

The source was deleted some time ago, due to the above mentioned errors. No new version is yet publicly available.

Best regards

Niels

[icingfish]

Originally Posted by ngj

Hi

The source was deleted some time ago, due to the above mentioned errors. No new version is yet publicly available.

Best regards

Niels

OK！thank you, Niels
i'm new and haven't looked over all the posts. now i'll continue do it.

[tfuwa]

Originally Posted by ngj

Hi all

I have been in the process of reimplementing the above mentioned wave boundary conditions into a completely new structure, which for instance make use the runTime selectable tables.

In the process of this reimplementation I realized that some of the approaches I did take was erroneous, and thus the results from the current available wave BC must be considered with some care and tested carefully.

Unfortunately I am not able to deliver the necessary patches to the current BC and the new implementation will not be made available for some time.

This message should be considered as a warning, and users are encouraged to make thorough testing before actual use.

My hope is to release the new version at some point.

Best regards,

Niels

Hi Niels and All,

Is new wave BC available now? I am interested at internal wavemaker in OF. If would be great that any of you could suggest any related reading materials. Your help would be greatly appreciated.

Kind regards,
Albert