[velan]

Hi Mark,

I compiled those sources in 1.5 version. But i am facing problem in 1.5-dev version.

Thanks for your sources.

- Velan

[norman1981]

Hi Velan,

you are right. I succeeded in compiling this nice solver some weeks ago in this way (if I remember right):

1- compile ODE v1.5 library in v1.5-dev (give it another target name to avoid a replacement of the original one ;-) )
2- modify the Make/options file in order to make the solver link to the 1.5 ODE library
3- wclean and then wmake

Maybe there is a cleverer way to solve this problem but in this way it works for me.

Kind Regards

Norman

[albcem]

Hi Mark,

I managed to modify your code so that multi-patch objects can be passed as moving bodies. Now I assign a motionPatches list under each object to move in the dictionary.

As I worked on the code, I realized some corrections had to be made to account for cases when a simulation is stopped & resumed from the last saved time step. Essentially the center of gravity CoG variable has to be initialized from the last saved time step as:

if (CoGtrack->found("totalDisplacement"))
{
totalDisplacement = CoGtrack->lookup("totalDisplacement");
CoG += totalDisplacement;
}

else the moments computed are relative to the CoG at time-step 0 and wrong. I remember observing a jerking type motion on the cube example when I resumed say at time=3s which is advanced enough for the CoG position to be significantly different than at time=0s. As a result the block is subjected to some extreme moment causing the jerky rotation.

To preserve exact continuity for resumed runs, you may also want to preserve and restore Fstore[0], Mstore[0] & Fstore[1], Mstore[1] variables in the CoGtrack dictionary.

Also the limiter on the mesh motion will need to be made ineffective with something like:

if (counter > body0.startUpdate() || mesh.time().value() > mesh.time().deltaT().value() * body0.startUpdate() )

Lastly, while they are small, I do not see the viscous forces and moments being lumped into Fstore[0] and Mstore[0] to affect the object motion. Is this intentional?

Cem

[markc]

Hello Cem,

Thanks for the constructive remarks. With respect to viscous forces: I deliberately excluded them as I do not have a lot of experience with these forces and turbulence models. The best way is probably to make a switch in the shipDict, saying if you want to include those forces or not.
Further more: you are fully right with respect to the issues regarding a restart. I tried to modify the solver but ran into a few problems. First for reference how it is right now:
A IOdictionary is initialized:
CoGtrack = autoPtr<IOdictionary>
(
new IOdictionary
(
IOobject
(
fileName_,
mesh_.time().timeName(),
mesh_,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
)
)
);
This file is written in every time directory and will obtain some relevant simulation parameters. The writing is done as follows (for a few random variables):
>>>
CoGtrack->add("velocity",U_cog_old);
CoGtrack->add("rotationSpeed",Omega_cog_old*180/pi_);
CoGtrack->add("Fstore",Fstore);
<<<
Where U_cog_old and Omega_... are vectors and Fstore is a pointField, defined as:
>>>
Fstore = pointField(3,vector::zero);
<<<
Now the problems:
1.
Writing of Fstore goes fine. However reading gives problems. I try to read as follows:
>>>
if (CoGtrack->found("Fstore"))
{
Fstore = CoGtrack->lookup("Fstore");
}
<<<
The compile error:
>>>
error: ambiguous overload for ‘operator=’ in ‘((Foam::bodyMotion*)this)->Foam::bodyMotion::Fstore = ((Foam::bodyMotion*)this)->Foam::bodyMotion::CoGtrack.Foam::autoPtr<T>:per ator-> [with T = Foam::IOdictionary]()->Foam::IOdictionary::<anonymous>.Foam::dictionary: :lookup(((const Foam::word&)(& Foam::word(((const char*)"Fstore")))), 0)’
<<<
Apparrently reading a pointField from a dictionary needs some other approach than reading a vector?

2.
The CoGtrack file is written in every time directory at write-time. This works but when running in parallel it is only written in one of the decomposed time directories. And reconstructing the decomposed case, the track file is not put in the reconstructed time directories. So writting the CoGtrack file in the time directories is not the good approach. One would like to write it in a dedicated directory, e.g. "CoGtrack" which will obtain subdirectories at write-time with the name as timeName(). Unfortunately I can't get this to work. I tried to initialze the IOdictionary like:
>>>
CoGtrack = autoPtr<IOdictionary>
(
new IOdictionary
(
IOobject
(
fileName_,
"CoGtrack"/mesh_.time().timeName(),
mesh_,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
)
)
);
<<<
This seems not to be a very difficult task but yet it fails during runTime. Anyone knows how to perform this task?

Brds,

Mark

[albcem]

Hi Mark,

The work-around I found for your first issue, ie writing/reading of the Fstore and Mstore variables (tensors), was to write and read the vector components:

Fstore[0] is written to fStore0
Fstore[1] is written to fStore1
Mstore[0] is written to mStore0
Mstore[1] is written to mStore1

I don't think we need the third component as it is directly over-written.

For the second issue, I actually did not resume a run from a reconstructed parallel run, so I did not realize there was an issue there. I will see if I can produce something useful about it.

Now that you raised the turbulence issue, all I did to activate it was to change the UEqn.H file by adding the following to the top

surfaceScalarField muEff
(
"muEff",
twoPhaseProperties.muf()
+ fvc::interpolate(rho*turbulence->nut())
);

and replacing mu with muEff. I think this is the only difference between the rasInterFoam and InterFoam solvers...

You are right about putting some special consideration to the compuation of viscous forces. With the turbulent flow field solution based on muEff and utilizing the original mu to compute the resulting viscous stresses, we have some solution for viscous stresses, but I can't judge its completeness. I have seen an old liftDrag.C code where laminar and turbulent viscous stresses are computed separately. However, I don't think the naming convention is right there as the velocity field is a turbulent solution in both parts.

Cem

[velan]

Hello Mark,

Can you please tell me, which file reads the contents of environmental properties. Today i was changing the value for wave number, depth....etc in environmentalproperties file. But i havent find any difference. I would like to know, is your solver reading those values ?

Thanks
- Velan

[ngj]

Hi Mark and Velan

@Velan: I can easily answer your question, as I can recognize the input, and I believe Mark has been using my wave generation tool. Unfortunately I have removed the tool from my website as I "found" a set of errors while implementing it again in a rather different approach (see disclaimer elsewhere on the forum). Furthermore the tool is not included in shipFoam, and hence the parameters in environmentalProperties are not read.
As soon as the new wave boundary conditions are mature for release they will be so through the svn under openfoam-extend. The time horizon however, is not known at present.

@Mark: I am trying your tutorials with shipfoam and it runs smoothly. Thank you very much for sharing, it is amazing On the Forum you did ask at some point about the gamma-quantity exceeding 1. I can see from your tutorials that this does not occur anymore. At present I am trying to update the bed level (bathymetrical change) due to the sediment transport induced by wave action. However I experience the exact same problems with diverging gamma.
Would it be possible for you to tell me what you did to achieve that?

Thank you very much for any help.

Best regards,

Niels

[markc]

Hello Niels,

It's a pity you was not in Montreal on the conference. We formed a SIG on the field of 2-phase flows. We defined our common interests. You can find some information via the wiki, where we have a dedicated area. You might be interested in joining.

With respect to boundedness of gamma:

• I found that the gamma compression term is important for stability. The lower, the more stable. But at the cost of a more "soft" interface. I think however that in the demo's I still used cGamma=3. Nowadays I use maximum 1.
• using upwind provides more stability, though accuracy may be impaired.
• I do not know from when you found my question but in shipFoam I added the possiblity to use relaxation on U and pd. I found this is not implemented in the OF1.5 interFoam solvers. I am still wondering whether I am right here, or if I overlook something or that there may be a fundamental reason to not use relaxation in VOF solvers. However it greatly enhances stability, especially with the moving meshes.

Hope to helped you a bit further,

Mark

[ngj]

Hi Mark

Thanks for your quick reply. I actually found the SIG earlier, and I have now joined the mailing-list.
Just to help narrowing down the problems, I have used cGamma = 1 up to now, thus that seems not to be of significant importance for the stability (?), hence the solution to the problem might be due to the relaxation. I will try it as soon as possible and when I have progressed I will report back to you.

Thanks again,

Niels

P.S. I hope you all had a good workshop

[matteoL]

Dear all,
I am trying to run some simple steady S60 boat cases with OF-1.5-dev and I have strange pressure issues:

Please check if my following statements are correct:
1) if the box domain is big enough, theorically the pd (pressure minus hydrostatic part should be zero on the boundary (but the symmetry plane ofc))

2) p is useful only for the evaulation of forces, otherwise it is never really used in the solver. In fact in the dev version of rasinterfoam the line p=pd+rho*gz has been cancelled (p is never updated).
This is a problem because when i use the force lib it uses p, and thus says that the pressure forces are always constant (I have thus readded the line..)

-Now, if i put as bc for the pd: zero gradient at the inlet and outlet, I obtain that the pressure at teh outlet and inlet are quite big, mainly at the outlet,like if there was some sort of blocking..
After a while this make the freesurface rise all over the domain!!!

-I have thus tried to impose pd= fixed value=0 at inlet-outlet -side and at the beginning it seems to work much better, but then after a while it creates high values to the cells just before (after) the outlet (inlet)...
and the explodes..
I can understand this if i think that my bc are too strongly imposed, and I am completely blocking any sort of perturbation (numerical or physical due to the wave arriving from the boat) and thus i should impose something weaker..

I have thsu tried the BC:fixedmeanvalue=0 strategy but this is even worse and teh code explods very quickly.

Any ideas?

Finally, Mark I cannot download your shipfoam code posted before (it appears as a .unk file rather than .tar.gz ), Am i doing somethign wrong or is there some problem?Could you post it again pls?

Thank you very much,
best regards,
Matteo

[norman1981]

Dear Matteo,

regarding your last question, rename the .unk extension to .tar.gz.

Kind regards

Norman

[matteoL]

Thank you very much Norman. i should have suspected that it was something easy...

Another point: Has anyone managed to compile it under OF-1.5-dev? Also with the changes proposed in this Thread, I still have problems running it.. It seems they have changed something in the mehm motion (in particular in the laplace problem)
and the code crashes at the first iteration...
Any idea?

thanks,
Matteo

[matteoL]

To answer myself :
-I has strange issues with the od variable because I was using the groovyBC BC and thus had one single inlet domain. thus I had no reason to place a cell layer at exactly z=0 (i.e. freesurface level). For some reason which are not completely clear to me, this was creating problem with the hydrostatic pressure subtraction...

Anyway, I have now set a layer of cell at exactly z=0 and that problem now seems to be solved...

ciao,
matteo

[matteoL]

Hello,
Just to let you know that I think I have found a bug int he evaluation of the forces when usingthe libforces library with rasInterfoam.

In fact the viscous terms are evaluated multiplying the dynamic viscosity with rhoRhef which is correct only if you have one single phase. In the case of rasinterfoam, I think you have to multiply it with the real rho.

I have thus added some line the forces.C and forces.H so that it "load" also the rho field and multiply it correctly (mainly the importat change is line 321 that becomes vectorField vf = (Sfb[patchi] & devRhoReffb[patchi])*rho.boundaryField()[patchi]; )


There are some more lines added changed to load files, modify constructor etc..
but i think i have put enough comments in the files to make it easily readable.

Here is the changed lib:

To use it, just go in it and type "wmake libso". it will thus create a lib called libforcesrasInterFoam.so in your $foam_usr_lib directory.

Then as usual in the system/controlDict you have to load it instead of the usual libforces.so and set rhoref=1.

p.s: I have just modifed the forces class, not the forcesCoeffs (anyway it is straightforward to modify that one as well)

p.p.s: this libs works only for rasInterfoam, it may not work for other solvers..(in that case use the original one).

I hope this helps and please let me knwo if you think my modification is not correct,

ciao
Matteo

[ngj]

Hi Mark

I have made some progress, and I have found a method to update the mesh without the need of any relaxation and better behaviour in gamma. Up to now I have been working with periodic mesh updating, i.e. updating every elapsed N seconds. It seems that not updating every time step is the cause to all of the problems, as mesh.V0() is only updated when the mesh is moved, and hence gamma correction is carried based on an erroneous dV/dt in the time steps, where no mesh movement is computed.

I have tried to run your test with the dropping of a sphere, and as you can see the effect of the relaxation is significant both with respect to displacement amplitude and with respect to the oscillating frequency.

However there are still some effects on gamma becoming larger than 1, but I have not analysed whether or not this could have something to do with the shear magnitude of the displacement or another reason. The positive thing is that the solution is still stable, i.e. gamma overshoots and converges toward 1 afterward.

Best regards,

Niels

[markc]

Hello Niels,

Nice work again. However, I do not fully understand what you're meaning. Are you now updating the mesh every time step or only periodically?
I played with updating only every N steps as well. At the moment of updating the mesh, peakpressures are introduced which get the time to settle down in the next steps. I was not in favor of this method, as you still introduce instabilities. Updating every time step seems more natural to me. We still intend to do some investigation about the amount of relaxation which is needed. I think a decent study object would be a falling cylinder or wedge, from which experimental data is available. You can than study the influence of relaxation both on motions and on (peak) pressures. Note that the amount of relaxation which is used in the "tutorials" are just based on the experience that with those values you get stable results.

Brgds,

Mark

[ngj]

Hi Mark

I see that I forgot to include some informations.

The curve "Relaxation" is based on your tutorial with updates every 5 time steps and the standard relaxation parameters, whereas "No relaxation" updates every time step and uses no relaxation.

The reason for not pursuing the relaxation approach was that I experienced that the relaxation changed the propagation speed of surface water waves, which actually does make sence, however that changes the entire response, and hence it was not a feasible way.

I have achieved to bridge the two methods, i.e. only periodic updating and updating every other time step with small time consumption. This is done by setting the moving patch velocities to 0 and multiply the internal field, e.g. motionU, by 0 in these time steps. It leads to a method where the only time used is to set up the motion matrix and verify that the old initial solution, i.e. the constructed 0-field, fulfills the requirements. Actually I think it would be sufficient to do this trick twice after the actual update, as mesh.V0() and mesh.V00() would be filled with the correct values. It does, however, still lead to pressure peaks, but for some applications this might be accepted.

Best regards,

Niels

[kjetil]

Originally Posted by hjasak

I wrote a utility for this aaaages ago and it is documented on this Forum. Have a look at:

http://openfoam-extend.svn.sourceforge.net/viewvc/openfoam-extend/trunk/Core/Ope nFOAM-1.4.1-dev/applications/utilities/postProcessing/stressField/interFoamPress ure/

As for the pd + rho*g*h, how do you plean to include surface tension pressure jump or deal with the situation where you have a multiple layers of air and water for example. So, it is not rho*g*h but integral of rho*(g & x), which is MUCH more difficult.

Enjoy,

Hrv

Has readEnvironmentalProperties.H been substituted with another file in the 1.6 release? I try to rebuild this utility, but this file is not present in the file structure - but is so in the previous one. Any hints how to replace/fix this?

[ngj]

Hi Kjetil

As an easy fix, there are three possible solutions to this:

1. If readEnvironmentalProperties.H is still present in 1.6, then find this new location and source it.
2. Copy readEnvironmentalProperties.H to the location of the utility, you are building and try again.
3. Look in the source to interFoam and see if readEnvironmentalProperties.H is read. If not, find out where the gravity vector is read from and substitute #include "readEnvironmentalProperties.H" in your utility with the current file.

Good luck,

Niels

P.S. However, as I understand interFoam v. 1.6, there are being solved for the total pressure, so what is it exactly you are after?

[kjetil]

1. File is not found. However, it now seems like it has been replaced by readGravitationalAcceleration.H.

About what I'm after ... please see attached screenshot. I would like my flow to start from zero and accelerate up to a steady state. However, pressure turns out to be zero in entire field, though I have a varying velocity field. So I read that there is a utility to find hydrostatic pressure when using interFoam (though I'm using the new 1.6 version called compressibleInterFoam). But this utility was not included in 1.6 ...only in the "extended" version which is based on 1.5.