[Ya_Squall2010]

Hello,

I wrote a custom b.c. named "regularWaveMakerFixedValue" which is used to generate regular waves at the inlet. It can be easily extended to generate other kind of waves. Basicall it applies solid wall (U=0) at the air part of the patch and wave velocities computed from wave theory at the water part of the patch, just the same way, at least I believe so, as what GroovyBC does at the boundary. But as can be seen from the attached pictures, due to the strong vortex generated in the air phase near the free surface within the domain, the wave shape has been seriously corrupted, and soon after quarter of wave period, negative value of VOF appears at this location (that is the big visible unphysical ripple right next to the boundary) which leads to the crash of the solver. Could anybody kindly advise what's going wrong here or what should I do to work it around? many thanks! I've attached the source code of the b.c. as well as my initialization & setup files.

b.t.w, I also tried to apply inflow&outflow boundary condition at the air part of the patch but got quite similar results.

[ngj]

Hi

Are you also applying the boundary condition on the VOF-ratio, alpha? If not, doing so might solve your problem.

Best regards,

Niels

[Ya_Squall2010]

yes, I applied the bc to both alpha1 and U.

[Ya_Squall2010]

so far, I've changed different kind of waves and all got similar error. Could any one who's done wave simulations give me some advice on the right chioce of schemes and parameter setup? Many thanks!

[Ya_Squall2010]

now I've tried to force the air phase velocity equal to zero at the end of each time step, and also got same problem. A 'warm' start can not solve it as well.

[Ya_Squall2010]

attached is the wave variables time histories. 2nd Stokes wave was used. Nothing strange in the profiles.
export.png

[ngj]

Hi

I do remember similar problems, however I do not remember what solved it in the end. Have you tried ramping the solution over one or several periods by multiplying with a smooth function?

The problem seems to occur where the velocity gradient is the largest. Are you also applying the vertical velocity component from the wave theory on the boundary?

Best regards,

Niels

[Ya_Squall2010]

Thanks for attention. If you are talking about the "warm" start, that is multiply the wave variables by a damping function in the initial couple of waves, yes, I did try, but no use. As can be seen from the snapshot in my latest post, the solver crashes just when the surface elevation is about to change its moving direction, but the velocity components don't. I wonder if this is due to the special shape of the wave shape being chosen. I shall try another waves next. I did apply both horizontal & vertical velocity components on the boundary.

I considered the possibility that the mesh size is not enough to resolve the wave profiles between two wave crests. Please see attached for a fine grid result, which seems even worse and the solver crashes a little bit earlier.
export.png

[Ya_Squall2010]

Originally Posted by ngj

Hi

I do remember similar problems, however I do not remember what solved it in the end. Have you tried ramping the solution over one or several periods by multiplying with a smooth function?

The problem seems to occur where the velocity gradient is the largest. Are you also applying the vertical velocity component from the wave theory on the boundary?

Best regards,

Niels

it would also help if you could let me know you fvscheme and fvsolution setup for wave simulations. Thanks!

[Ya_Squall2010]

update:

Today I tried following wave profile, the small hump at the thought was removed. But the code still cannot proceed after the instant shown. Big nagetive alpha1 appears at the same place as before, and as can be seen from the snapshot, abnormal big velocity was resulted. I was wondering if this is due to the inaccurate solving of the velocity field, therefore, I tried to change the default PISO scheme in interFoam to PIMPLE and hoped accurate velocity field can be predicted by introducing more outer iterations. But the experiments told me that it didn't help. I got pretty much the same error and the only 'improvement' was that the crash happened a little bit late than before. That's all.

results of pimple+vof
1.jpg
2.jpg

[ngj]

Here are my fvSchemes and fvSolution. I have made a run using your parameters and it works in my end using my implementation.

One thought: Are you sure, that you compute the wave number correctly? I get 1.358925 m^-1. I have noticed that if omega and k does not fulfill the dispersion relation, the solution typically diverges.

Best regards,

Niels

[Ya_Squall2010]

Thanks，I shall have a look at my wave number and will post mine soon.

[Ya_Squall2010]

Hi, I found my wave number is 1.358925, which is the same as yours. I've attached my implementation of wave maker B.C., is it possible for you to compile it and have a go on your computer to see what you get? Alternatively, you could post your files and I will have a try. Many thanks!

just uncompress the files to .../OpenFOAM-1.6/src/finiteVolume/fields/fvPatchFields/derived/, and change the make file accordingly.
regularWaveMakerFixedValue.tar.gz

[ngj]

The only reason I can come up with for your problems is that you for some reason or the other have a force imbalance, which is tried to be corrected through the presence of large velocities, i.e. large shear stresses.

Hence your gamma distribution or momentum flux at the boundary seems not to model what you are looking at. This could be either due to missing pressure gradient at the boundary or the way you determine wet/dry cells. In my implementation I consider "continuous" variation in alpha instead of binary variation. Adding a pressure gradient seems to have limited impact.

Good luck,

Niels

[wavemaster]

Hi,
Hope I am not too late here to get some help...
I have implemented a waveBC similar to the one shown above, and seem to face the same problem.
What puzzles me most is that the same case and BC implemented in groovyBC runs fine. In the picture You can see that the same velocities and alpha conditions are applied, but inside the domain the surface elevation develops different.
Ramping etc. has no effect either.
Is there something special in groovyBC which is not obvious to understand when setting time dependant values on the Patch?
Are waveBC so sensitive that a slightly different value for solving the same equation in groovyBC or in Foam might give wrong results in one and work perfectly in another?
Any other idea?
Thanks in advance for any hint...
Pal

[ngj]

Hi Pal

You need to give more details. The waveBC implementation linked to above is for infinite deep water, so that might explain your problems, if you are comparing with a groovyBC version on finite water depth.

Best regards,

Niels

[wavemaster]

Hi,
OK, I will try to be more precise.
I implemented and used a simple waveBC in groovyBC, it works fine.
alpha1

        type            groovyBC;
        valueExpression "(pos().z<=Depth+A*sin(w*time())) ? 1 : 0";
        variables       "Depth=1.2;l=1;A=0.05;Endtime=2;rampfactor= ((time()/Endtime)< 1 ) ? (time()/Endtime) : 1;g=vector(0,0,-9.81);k=2*pi/l;w=sqrt(k*mag(g));";

U

     type            groovyBC;
     valueExpression "(pos().z<=Depth+A*sin(w*time())) ? vector( A*w*exp(k*(pos().z-Depth))*sin(w*time()), 0,A*w*exp(k*(pos().z-Depth))*cos(w*time())) : vector(0,0,0)";      
          
     variables       "Depth=1.2;l=1;A=0.05;Endtime=2;rampfactor=((time()/Endtime)< 1 ) ? (time()/Endtime) : 1;g=vector(0,0,-9.81);k=2*pi/l;w=sqrt(k*mag(g));";

Then I wanted to start coding and learn more and implemented it in a waveBC
alpha

        forAll(result,i)
        {
            if (ZCo[i] <= waterLevel_+*amplitude_*sin(omegaTime))
            {
                result[i] = scalar(1);
            }
            else
            {
                result[i] = scalar(0);
            }

and U

        forAll(result,i)
        {
            if ( ZCo[i] <= waterLevel_+*amplitude_*sin(omegaTime) )
            {
                Depthexpress = exp( wavenumber*(ZCo[i]-waterLevel_) );

                result[i] = vector(
                        amplitude_*omegawave*Depthexpress*sin(omegaTime),
                        0,
                        amplitude_*omegawave*Depthexpress*cos(omegaTime)
                                  );
            }
            else
            {
                result[i] = vector(0,0,0);
            }

When I compare in paraFoam the boundary patches (alpha1 and U ) over time look exactly the same (Image), but after some timesteps the surface behaves differently, the formula in groovy looks nice, the case with my waveBC creates strange velocities.
If I visualize only the patch, the values shown in paraview are the ones created by the BC, or?
So, on the one hand the BC seem to be equal, but on the other hand the results are not. Thus the question if there is some trick in groovyBC...
Thanks for Your input,
Pal

[ngj]

The velocity component in the x-direction for the 2 implementations differ by a "-".

Best regards,

Niels

[wavemaster]

Sorry, this was a copy error. It did not affect my simulations.
I had just tried waves running in different directions, with groovyBC all is fine and again my brilliant piece of code fails.
If I set the amplitude to zero, there is no velocity or alpha movement visible on the patch , but I still get crazy velocities with my stuff...
Thanks for looking at it...
Pal