[miliante]

Good day for you OF Users!!!

I am merging several solvers in order to simulate a 2d roll damping in ships. I have successfully integrated a 3DOF movement for my ship with moving mesh by calculatin and integrating the forces and momentum. Gravity is playing its part in the calculation of these forces, but not in the fluid. Water is just flying away when it hits the ship.
I noticed that when I merged the solvers I had forgotten the gravity (- ghf*fvc::snGrad(rho)). But now there is no Hydrostatic Pressure to counter act gravity and the fluid just keeps practically static due the action of gravity pressing down.
Can someone help me on how to add hydrostatic pressure (inside the domain and as a boundary conditon). I read on the posts that OF uses a "modified pressure" which should account for everything, but then how should I use it in my case?

Thanks in advance.

T.

[egp]

Hi Thiago,

We have recently validated a rasInterFoam derivative for role damping (both free roll and roll in beam seas). The total pressure (p) is recovered by added the hydrostatic pressure back into the specific pressure (pd). For surface tracking methods, this can be done in a brute force way, p = pd + rho*g*z. However, for a surface capturing VOF method with a complex free surface (e.g., one that includes wave breaking), a more elegant approach is required. We are reconstructing p by solving an additional Laplacian at the end of each time step.

{
rho = gamma*rho1 + (scalar(1) - gamma)*rho2;

// Add gravity term, re-using the momentum equation
UEqn += interface.sigmaK()*fvc::grad(gamma) - rho*g;

volScalarField rUAPres = 1.0/UEqn.A();

// Pick up velocity boundary conditions
volVectorField Upres = U;
Upres = rUAPres*UEqn.H();

surfaceScalarField phiPres = fvc::interpolate(Upres) & mesh.Sf();

for(int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rUAPres, p) == fvc::div(phiPres)
);

pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
}
}

I will be presenting some of this work at the OpenFOAM workshop in Milan. Also, I am preparing a tutorial to distribute to the participants of the workshop.

Regards, Eric

[miliante]

Hello Eric,

All be damned!!! It worked. But if it worked as it should, there is something wrong with my calculation of forces, because it is too dissipative now. It should die within 20 secs and the ship stops rolling in less then 5 secs.
When I calculate the pressure forces I sum:
p.boundaryField()[something] + (g & mesh.Cf()) * rho.boundaryField()[something]

Perhaps it is something else. I'll keep looking into it and post back my results.

P.S.: It is a pity Milan is too far to go, wish I could your presentation.

Sincerely, T.

[nicasch]

HiEric and Thiago!

A couple of questions, to clear all misunderstandings.

Once p is reconstructed as described above (by solving an additional Laplacian at the end of each time step), shouldn't then values of p be exactly equal to pd + rho*g*z?

If not, why not?
If yes, why bothering with solving an additional Laplacian, why not simply calculate p as: p = pd + rho*g*z?

Regards.

[egp]

nicasch,

Yes, you could reconstruct p = pd + rho*g*z. This would be fine for small-amplitude clean waves and free-surface disturbances. However, it is not so clear when you have breaking waves, splashes, and sprays. Therefore, solving for p is a more general approach.

Eric

[nicasch]

Hi Eric, thank You for Your answer!

It seems that solving for p is generally a safer approach.
I have only one more question regarding boundary conditions for p, it would be nice if You shared Your opinion. Looking at interFoamPressure application, p is defined as:

volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
pd
);

This means that boundary conditions for p are same as for pd. Now, if You solve for p during Your simulation, before the end of the time step (after pd and U are solved), and You define p as in interFoamPressure, this would mean that pd and p would have the same values at boundaries. Is this correct?
What I am worried about is that the solution for p might not yield values which are equal to p=pd+rho*g*z. Is this correct?
Does this mean that, after U is solved, You are only seeking for pressure field p which exactly correspond to calculated U field? If this is so, than one may specify an open boundary for pd (and p) with same boundary values?
For example, in damBreak-case there is an open boundary with prescribed pressure pd. If You have same boundary value for p at this open boundary, then p will not be equal to p=pd+rho*g*z in the domain.
Is this phisically justified? Is it obligatory to obtain values for p equal to p=pd+rho*g*z, or not?

Thanks in advance!

[miliante]

Hello Eric and nicasch,

nicasch questions were very interesting and it allowed me to better understand what is going on. I thought everything was working properly, but I finally decided to check the p field. I isn't ok and I think it is because of boundary conditions. I am using zeroGradient for the sides and bottom and totalPresure for the top. What kind of BC are you using?

Thanks in advance.

T.

[kev4573]

Zero gradient should be okay for the sides however try using symmetryPlane for the bottom.

[nicasch]

Kevin,

do You mean symmetryPlane should be set only for p equation or for all fields? If symmetryPlane is set for pd and, even worse, for U and gamma at bottom wall, this would not be right!

If symmetryPlane should be set only for p equation, how can that be done? If You set base type patch for bottom wall in blochMeshDict, then it is not allowed to specify symmetryPlane at that patch. symmetryPlane has to be defined already in blochMeshDict, and therefore it holds for all fields.
How can this be separated, i.e. how can You set the following B.C. at one and the same boundary:
U - fixedValue
pd and gamma - zeroGradient
p - symmetryPlane?

If someone knows the answer, please share it.

Regards.

[kev4573]

You are right that symmetryPlane must be set for all variables if you use it on a patch. I'm unclear myself on why this boundary condition works so I'd be interested in any insight as well.

[miliante]

Well people,

I am also looking for a "correct" solution for this problem. But right now i cheated...now I have a field just for p with a fixedValue for the bottom with the calculated (as in my pocket calculator) value for the pressure. I is now giving me what seems to be the correct pressure p (p = pd + rho*g*h) and the simulation is much better, but I am still fighting with gravity since the fluid still "flyies" away.

Thanks people.

T.

[miliante]

People...

I'll give up for a while. I never worked with free surface and moving mesh before. I'll spend some time reading more. I am leaving my files on my website, in case someone want to check, or compare, or...well, have fun!!!

http://reis.miliante.googlepages.com/software2

T.

[miliante]

hello people!!!

Has someone found out how to solve the p field correctly? Even "cheating", I am loosing pressure thru the simulation leading to wrong results.

Sincerely,

T.

[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

[miliante]

Thank you Hrv,

but I have already seen and downloaded the "interFoamPressure" file. My problem is that, probably duo to boundary conditions, the field does not converge correctly and the p field is useless during simulation. The use of simmetryPlane on the bottom does solve correctly the p field, but the overall simulation diverges. And Using the fixed values for the pressure on the bottom and top for the p field give me much better results, but there is some kind of "pressure loss", leading to error further in the simulation, which I have not looked into very deeply.

Thanks again for the reply.

T.

[miliante]

Never mind it. It is working, there is no "pressure loss". I was so used on working without the freesurface that I forgot that energy is dissipated just by making these waves. I just have to find out why in the beginning it is very close to the experimental result and then everything screws up.

Thanks people.

[hjasak]

In those famous words from Sreenath: "It Works For me!"

In any case, you have a possibility of defining your own boundary conditions on p if the ones from pd are not OK. Also, I believe Eric did some tests and the method seemed to work.

Please keep me posted if you can confirm the pressure field is indeed the problem.

Hrv

[miliante]

Hello!!!

The problem is in the part proposed by Eric, not that it isn't correct, it should work. The problem is probably on my coding. Perhaps boundary condition or the fact that the mesh is moving.
What happens is: I start the simulation with a 10 degrees displacement and let it roll. But it only goes 'till -7 degrees and back to 10 degrees. I turned off dynamic pressure to prevent losses. It should oscilate between 10 and -10.

If anyone got it working on a moving mesh and could give me some help, I would appreciate it.

Thanks in advance.
T.

[egp]

I'll post both a tutorial (free roll of a box barge) and a rasInterDyMFOAM6DOF Solver as soon as I finish my slides for the Milan workshop. Since they are due tomorrow, I'll post something in a day or two, and check it in to the SVN repository.

Eric

[miliante]

Thanks Eric, that would be of great help.
Best on your presentation.