CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Hydrostatic Pressure and Gravity (https://www.cfd-online.com/Forums/openfoam-solving/57774-hydrostatic-pressure-gravity.html)

miliante May 29, 2008 10:10

Good day for you OF Users!!!
 
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 May 29, 2008 10:38

Hi Thiago, We have recently
 
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 May 29, 2008 14:23

Hello Eric, All be damne
 
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 June 3, 2008 12:17

HiEric and Thiago! A couple
 
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 June 4, 2008 06:38

nicasch, Yes, you could rec
 
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 June 5, 2008 06:05

Hi Eric, thank You for Your an
 
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 June 10, 2008 15:49

Hello Eric and nicasch,
 
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 June 10, 2008 19:32

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

nicasch June 11, 2008 02:46

Kevin, do You mean symmetry
 
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 June 11, 2008 16:31

You are right that symmetryPla
 
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 June 11, 2008 18:29

Well people, I am also l
 
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 June 12, 2008 14:08

People... I'll give up f
 
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 June 19, 2008 12:03

hello people!!! Has some
 
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 June 19, 2008 12:20

I wrote a utility for this aaa
 
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 June 19, 2008 14:17

Thank you Hrv, but I ha
 
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 June 19, 2008 16:13

Never mind it. It is working,
 
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 June 20, 2008 08:59

In those famous words from Sre
 
In those famous words from Sreenath: "It Works For me!" http://www.cfd-online.com/OpenFOAM_D...part/happy.gif

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 June 30, 2008 17:21

Hello!!! The problem is
 
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 June 30, 2008 18:59

I'll post both a tutorial (fre
 
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 July 1, 2008 07:37

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

markc July 13, 2008 07:13

Thiago, Two other questions
 
Thiago,

Two other questions about your case:
1.
Looking at your model setup in paraview it shows that the initial "corpo" angle of 10 degrees is simulated by having a liquid level which is higher on the left side than on the right side. Only very close to the corpo the watersurface actually "slopes" but for the remainder of the model the level is just different, but constant, left from right.
http://www.cfd-online.com/OpenFOAM_D...s/mime_pdf.gif snapshot.pdf
Is this correct?
2.
Absolute pressure, p:
If I am right, the hydrostatic pressure is set in the main solver, just before the time loop starts. Isn't it better to set the initial hydrostatic pressure in the case/0/p file before starting the solver?

Brgds,

Mark

markc July 14, 2008 16:33

Hello All, I made some prog
 
Hello All,

I made some progress, solved the issue with gamma values. What went wrong:
1.
Mixed up meaning of variables. Now:
>>>
cellDisplacement.boundaryField()[patchIndex] == vector(0,0,dZ);
pointDisplacement.boundaryField()[patchIndex] == vector(0,0,dZ);
<<<
As stated some days ago I was not aware that some unphysical values were added.
2.
Had to put movingWallVelocity BC on U.

But: we are not there yet. Next issue:
If I run the solver without actually moving the mesh (by changing dZ in the above two lines into "0"), the solution runs fine. If I print values for calculated dZ (calculated change of body), reasonable and stable values are given every time step.
However switch on meshmovement again, the solution for pressure (and consequently dZ, change in Z direction) immediately starts oscilating from +much to minus even more to plus yet even more until the solver stops.
I suspect the body (Hull) is moved but pressurefield is not correctly updated, causing either a "vacuum" or a very high pressure, which is being released the next time step, causing an even larger p value towards the other side. It does not matter how small dZ values I allow (by simply limiting them). It explodes sooner or later.
So it looks like a continuityError. We keep on investigating and any help is welcome!
Btw: I am very curious about EricÅ› tutorial.

Thanks in advance for any help,

Mark

nicasch July 15, 2008 01:57

Hello Eric and others. I am
 
Hello Eric and others.

I am very interested in your tutorial and your solver too, it would of great help to me.

Thanks in advance!

miliante July 15, 2008 10:32

Hey Mark, Good to know
 
Hey Mark,

Good to know people working on similar problems...

Here we go:

Post 1

1. Yes, in the end it is just: acceleration = Force / mass [m/s²] or acceleration = Moment / inertia [rad/s²], C = 0 and K = 0.
2. It is a 2D solver, so it doesn't move in the Z direction, only X (horizontal) and Y (vertical).
3. Yes, But there is only one patch, I am just too lazy to correct that.

I also noriced that. Perhaps the dZ is not small enough? Since it always skips the correct solution due to the moving mesh, it tends to accumulate an unrealistic acceleration. I haven't solved that problem yet, but I tried limiting the maximum allowed acceleration and it works better.

Post 2

1. That is odd. That happens because the boat should be rotated 10 degrees.
2. No, after much consideration, hydrostatic pressure should be disconsidered except when calculating the forces. So "pEqnGravity.H" is responsible for solving the hydrostatic pressure and assigning it to the variable p which is only used in "solveTotalBalance.H". The "pd" variable should only contain the dynamic part. But there is where my problem lies. It is not correctly solving the pressure "p".

I know I haven't helped much. I have been strugling with these problems the last few weeks without much success...and much time. But I am working on it. I'll post an updated version of the solver later today, preferably with everything in english.

Sincerely,

T.

markc July 15, 2008 15:00

Thiago and others, Also som
 
Thiago and others,

Also some progress here. I solved the issue with oscilating p by simply updating the meshmotion once in (e.g.) 10 time steps while the p is solved every timestep. Doing so the p-solver gets the opportunity to obtain stable values before it is used as input in the forceSolver and consequently as input for meshmovement. Not the nicest way but for the moment it works. Hopefully someone gets some better ideas or I will see the light some later.
Btw, in my case the movement will be relatively small so in fact updating the meshmotion less frequently makes the solver somewhat faster.

We continue working, as time allows.

Awaiting others comments / results,

Mark

markc July 24, 2008 02:32

Hello All, Some advice abou
 
Hello All,

Some advice about movingMesh needed.
I want my patch "Hull" to translate in Z and rotate around Y.
Translate in Z I managed, it works now. However next step is the rotation.
It should do something like:
dZ as a result of Y-rotation=Xi*sin(theta)
dX as a result of Y-rotation=Zi*cos(theta)
with Xi and Zi position of cellface on hull patch.
I would like to use a transformation matrix and the "transform" function.
For the movingMesh I am only using celldisplacement and pointdisplacement, no velocities. So I have vectorField celldisplacement and pointField pointdisplacement.
Next is to add dZ (which worked) and dZ and dX as a result of rotation around Y axis. Can anyone explain how to perform this task, prefereably in a generalised way (transformation matrix)?
Thiago, you are doing something like this in your forceSolver, can you elaborate about what is exactly happening?
Btw, can anyone explain the difference between and the exact purpose of the fields pointDisplacement and cellDisplacement?

Brgds,

Mark

markc July 30, 2008 03:29

Hello All, The struggle goe
 
Hello All,

The struggle goes on. At present my solver is able to calculate the required 2 DOF (trim and sinkage) for the ship problem. The mesh is updated correctly, everything seems to run fine, with sensible solutions. However sooner or later the solver crashes on floating point exception.
I update my free body solver every 10 timesteps and than the mesh is being updated.
One thing which might cause the problem (or give a clue) is the fact that maximum gamma becomes larger than 1 (e.g. 1.005 or so) after the first meshmotion and it does not go down again. Often the logfile shows max gamma of ...e299 the last time step before crash. I wonder if this is a root cause or a just a secondary result...
So far I tried to play with tolerances, relaxationFactors and divergence schemes. Tried the following schemes:
div(rho*phi,U) Gauss upwind
div(phi,gamma) Gauss vanLeer
div(phirb,gamma) Gauss interfaceCompression

and
div(rho*phi,U) Gauss limitedLinearV 1.0;
div(phi,gamma) Gauss Gamma01 1.0;
div(phirb,gamma) Gauss Gamma01 1.0;

With both I get max gamma which slowly creeps higher and higher.

Next:
I need hydrostatic pressure to calculate free body forces. I use the exact code as given by Eric in this thread. It works:
http://www.cfd-online.com/OpenFOAM_D...ges/1/8557.jpg

But every now and then some strange p-fields become visible in paraview. The picture below shows a p-field ploted in the same color scale as the previous picture:
http://www.cfd-online.com/OpenFOAM_D...ges/1/8558.jpg

And the next pic the erronoues p-field at a different color scale:
http://www.cfd-online.com/OpenFOAM_D...ges/1/8559.jpg

The pictures are with a solution time difference of 0.1 second. The next step everything is (almost) ok again.
Something similar also happens with pd field.
Does this have something to do with continuity errors? Or boundary conditions? As far as I can find, the latter must be correct, which is being confirmed by the fact that usually the results seem to make sense.

Hope that someone comments,

Brgds,

Mark

hjasak July 30, 2008 03:44

You probably have a boundary c
 
You probably have a boundary condition problem - it seems the solver completely blew up. I have seen this in cases where I have fixed boundary conditions, moving mesh and incompressible flow, when the mesh volume changes (this is physically impossible).

Definitely check your b.c.-s.

Hrv

markc July 30, 2008 04:30

Right! My boat consists out of
 
Right! My boat consists out of 2 patches (Hull and Deck). I only move the hull while the deck remains in a fixed position (deck is well above water, I am not interested in air resistance). So if the hull moves down, the total boat volume increases and the domain volume decreases.
BC's are same as I used for the case without movingMesh, except that the moving patch has movingWallvelocity. (velocity inlet, zerogradient outlet with fixed pd, atmosphere on top etc.
cell- and pointDisplacement all have fixedValue(0 0 0)
The solver allows the hull to move (and thus the domain volume to change) upto a certain extend: quite some mesh-changes have been performed before crash. But apparently there is a restriction?
I will work on moving both patches and report what happens,

Brgds,

Mark

markc July 30, 2008 08:00

Ok, no improvement. Both Hull
 
Ok, no improvement. Both Hull and Deck were moving correctly with tiny little steps. Simulation went a little further but in principle the same problems:
-max gamma above 1 after first mesh change and slightly increasing (from 1.0002 upto 1.05)
-some timesteps erroneous p and pd fields

Attached my boundary conditions.
http://www.cfd-online.com/OpenFOAM_D...s/mime_txt.gif cellDisplacement
http://www.cfd-online.com/OpenFOAM_D...s/mime_txt.gif U
http://www.cfd-online.com/OpenFOAM_D...s/mime_txt.gif p
http://www.cfd-online.com/OpenFOAM_D...s/mime_txt.gif pd
http://www.cfd-online.com/OpenFOAM_D...s/mime_txt.gif gamma
http://www.cfd-online.com/OpenFOAM_D...s/mime_txt.gif pointDisplacement

Note that my mesh is tet, on which checkMesh gives nil errors/warnings.
I read about the discussion on BC's in this thread but I do not see what is wrong on mine. Hopefully anyone can give advice?

Thanks in advance,

Mark

miliante July 30, 2008 09:24

Mark, the problem that yo
 
Mark,
the problem that you are having is the same I have. I was hopping Erik's tutorial would help me with that.
As a test, would you mind putting a "slip" boundary condition for U, p and pd for the bottom and run again. Let me know what happens.

Sincerely,

T.

markc July 30, 2008 13:25

Regret to inform you... Apply
 
Regret to inform you...
Applying slip BC on bottom did not work. Something interesting tough: the absolute pressure value on the bottom is now better, with zerogradient BC I missed about 1 m static head (out of 16) on the hydrostatic pressure on the bottom.
During the simulation however the p-solution goes all the way. Some timesteps it seems ok, then suddenly you get those fields as plotted above.
Strange: during some time steps p and pd were (almost) similar to each other.
Also: some timesteps p-fields more looks like gamma, though not exact.
Question: What is the difference between continuityErrors and movingMeshContinuityErrors.H?

For those who are interested I attach my solverfiles and case files, unfortunately without mesh data, as this is too large.
http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif shipFoam.zip
http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif Casefolder.zip

Hope someone helps us a step further again,

Brgds,

Mark

miliante July 30, 2008 17:22

Mark, That is really strang
 
Mark,

That is really strange. When I use slip, The solver runs much better, without the strange p-field. Nevertheless, p and pd should not be similar, if pd is only the dynamic part. But p like gamma should not be a problem, since p = rho*g*h and rho is proportional to gamma.

Don't know the answer to the question, but I haven't checked the files yet.

I'll check your solver later.

Continuously trying to solve the problem, sincerely,

T.

markc August 1, 2008 09:26

Hello All, Did some trouble
 
Hello All,

Did some trouble shooting (without finding the trouble).
Recall: modified interFoam solver. Included free body force solver. Forces are calculated using hydrostatic pressure (p). Forces are used to (slowly) move the free body (hull). After meshupdate solution becomes unpredictable.


1.
Ran solver without mesh update by limiting maximum mesh change to 0 (dictionary input).
BC's for bottom on U, pd and p: slip
divSchemes: rho*phi,U: Gauss limitedLinearV 1.0;
phi,gamma: Gauss Gamma01 1.0;
phirb,gamma: Gauss Gamma01 1.0;

Max Courant# 0.25
Results: floating point exception after 0.2 runtime seconds. upto then everything looked perfect, no strange p-fields, gamma bounded between 0 and 1 etc. But suddenly crash.

2.
As 1 but BC;s on bottom:
fixedValue on U
zerogradient on p and pd
Result: similar as 1.

3.
as 2 but with max Courant 0.5
Result: as 2

4.
as 3 but with divSchemes:
rho*phi,U: Gauss upwind;
phi,gamma: Gauss vanLeer;
phirb,gamma: Gauss interfaceCompression;

Result:
Solution runs perfect, exactly as with standard interFoam solver.

5.
as 4 but with mesh update every 10 timesteps.
See previous posts: funny p fields, unbounded gamma.

Conclusions:
Not very many I'm afraid but let's try some reasoning:
-BC's?
-p field calculation seems to work well without mesh update.
-divSchemes are important but this has nothing to do with meshmotion, it's interFoam in general.
-after the mesh motion, gamma is solved earlier than p. Gamma becomes larger than 1 immediately after the mesh motion. p depends on gamma and not the other way around. Now: do we have to look after p-solution? In my opinion it makes more sense to look after mesh motion and/or how gamma becomes unbounded.

Hopefully someone can give us some advice again,

Brgds,

Mark

markc August 1, 2008 16:14

One more test: Played with BC
 
One more test:
Played with BC on movingPatches. Found out that slip BC worked well. At least: the solver runs fine, p and pd fields are stable etc. Though with slip gamma still becomes slightly larger than 1, it does not continuously grow.
Not satisfying me to use a BC which does not make sense, but it may give a clue to experts...
Brgds,

Mark

markc August 10, 2008 12:09

Hello All, I can still use
 
Hello All,

I can still use some help....
The above described ship solver; updated to OpenFOAM-1.5; created a new solver based on the interDyMFoam and my old shipFoam solver.
So now everything running in OF1.5.
Unfortunately no changes on my results. In my opinion the main cause is the fact that gamma becomes larger than 1 immediately after the 1st mesh change. The solver still continues for quite a while but slowly gamma increases until the solution explodes.
BC's are suspected causes. However I can not find anything wrong. The moving patch (Hull) has movingWallVelocity in U, uniform(0 0 0). phi stays 0 at this patch during the solution, e.g.:
>>>

0.08/phi:
...
boundaryField
{
Hull
{
type calculated;
value uniform 0;
}
...
>>>

Fot those interested the solver files and case files (without mesh, which is too large), which shows BC's and solver settings etc.

http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif shipFoam.zip
http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif SFtest_2.zip

Tried to change the BC for Hull to wall, which actually did not change the situation very much, surprisingly: also here phi kept 0 for the moving patch (???)

For reference a part of the logfile just before and after the mesh displacement is applied the first time:

>>>
DICPCG: Solving for cellDisplacementx, Initial residual = 0, Final residual = 0, No Iterations 0
DICPCG: Solving for cellDisplacementy, Initial residual = 0, Final residual = 0, No Iterations 0
DICPCG: Solving for cellDisplacementz, Initial residual = 0, Final residual = 0, No Iterations 0
Execution time for mesh.update() = 0.65 s
time step continuity errors : sum local = 3.63498e-11, global = -1.35965e-11, cumulative = -2.57582e-09
GAMG: Solving for pcorr, Initial residual = 1, Final residual = 7.39126e-07, No Iterations 17
time step continuity errors : sum local = 2.7088e-17, global = -9.96824e-18, cumulative = -2.57582e-09
MULES: Solving for gamma
Liquid phase volume fraction = 0.890809 Min(gamma) = -4.86844e-20 Max(gamma) = 1
MULES: Solving for gamma
Liquid phase volume fraction = 0.890809 Min(gamma) = -2.44672e-19 Max(gamma) = 1
MULES: Solving for gamma
Liquid phase volume fraction = 0.890809 Min(gamma) = -4.72367e-20 Max(gamma) = 1
MULES: Solving for gamma
Liquid phase volume fraction = 0.890809 Min(gamma) = -6.21236e-20 Max(gamma) = 1
smoothSolver: Solving for Ux, Initial residual = 0.000501712, Final residual = 2.71806e-08, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.00171506, Final residual = 9.19563e-08, No Iterations 2
smoothSolver: Solving for Uz, Initial residual = 0.00186756, Final residual = 9.82137e-08, No Iterations 2
GAMG: Solving for pd, Initial residual = 0.0279634, Final residual = 6.05306e-07, No Iterations 8
time step continuity errors : sum local = 5.13561e-11, global = 3.64385e-12, cumulative = -2.57218e-09
GAMG: Solving for pd, Initial residual = 0.0044137, Final residual = 4.9019e-07, No Iterations 7
time step continuity errors : sum local = 3.97271e-11, global = 1.03753e-11, cumulative = -2.5618e-09
GAMG: Solving for pd, Initial residual = 0.00106885, Final residual = 7.32548e-07, No Iterations 5
time step continuity errors : sum local = 5.85887e-11, global = -6.32875e-12, cumulative = -2.56813e-09
DICPCG: Solving for p, Initial residual = 0.00022892, Final residual = 9.88458e-07, No Iterations 17
ExecutionTime = 77.03 s ClockTime = 77 s

Courant Number mean: 0.0125713 max: 0.143702
deltaT = 0.002
forceBalance

Fp = 961351
Mp = -2.93254e+06
XCB = -3.05043
forceBalance
dZ = 0.01
cumulative sinkage = 0.01
trim [m]= -0.946567
dPitch [radians] = -0.000626566
cumulative trimangle [radians]= -0.000626566
cumulative trimangle [degrees]= -0.0359007
Time = 0.022

DICPCG: Solving for cellDisplacementx, Initial residual = 1, Final residual = 8.43641e-07, No Iterations 61
DICPCG: Solving for cellDisplacementy, Initial residual = 0, Final residual = 0, No Iterations 0
DICPCG: Solving for cellDisplacementz, Initial residual = 1, Final residual = 8.55758e-07, No Iterations 60
Execution time for mesh.update() = 2.3 s
time step continuity errors : sum local = 5.8587e-11, global = -6.32857e-12, cumulative = -2.57446e-09
GAMG: Solving for pcorr, Initial residual = 1, Final residual = 5.84293e-07, No Iterations 17
time step continuity errors : sum local = 3.44426e-17, global = -1.21738e-17, cumulative = -2.57446e-09
MULES: Solving for gamma
Liquid phase volume fraction = 0.89079 Min(gamma) = -1.80021e-10 Max(gamma) = 1.00149
MULES: Solving for gamma
Liquid phase volume fraction = 0.890694 Min(gamma) = -1.25075e-09 Max(gamma) = 1.00297
MULES: Solving for gamma
Liquid phase volume fraction = 0.890599 Min(gamma) = -2.38389e-12 Max(gamma) = 1.00444
MULES: Solving for gamma
Liquid phase volume fraction = 0.890503 Min(gamma) = -2.34219e-12 Max(gamma) = 1.00591
smoothSolver: Solving for Ux, Initial residual = 0.0174879, Final residual = 2.69513e-07, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.0025134, Final residual = 9.39992e-08, No Iterations 2
smoothSolver: Solving for Uz, Initial residual = 0.00235204, Final residual = 9.58132e-08, No Iterations 2
GAMG: Solving for pd, Initial residual = 0.070277, Final residual = 6.11576e-07, No Iterations 11
time step continuity errors : sum local = 5.03159e-11, global = 5.16696e-12, cumulative = -2.5693e-09
GAMG: Solving for pd, Initial residual = 0.00953424, Final residual = 8.32796e-07, No Iterations 7
time step continuity errors : sum local = 1.14229e-10, global = 1.06366e-11, cumulative = -2.55866e-09
GAMG: Solving for pd, Initial residual = 0.00257332, Final residual = 7.13764e-07, No Iterations 6
time step continuity errors : sum local = 9.75252e-11, global = 2.52487e-12, cumulative = -2.55613e-09
DICPCG: Solving for p, Initial residual = 0.000279455, Final residual = 8.03508e-07, No Iterations 44
ExecutionTime = 85.96 s ClockTime = 86 s

Courant Number mean: 0.0130403 max: 0.477204
deltaT = 0.001
Time = 0.023
>>>

Hrv, is your case / solver which is shown in the video (dropping cylinder) publicly available somewhere? This example does everyhting I need, I probably overlook something small.

Btw, I was positively surprised by the new possiblities of OF1.5.

Hope someone can help me out here...

Brgds,

Mark

markc August 12, 2008 16:13

Hi All, For everyone who is
 
Hi All,

For everyone who is not on holiday and is interested here some small steps forward.
1.
Found that gamma stays bounded if number of gamma subcycles is set to 1 (and no more than 1).
Herewith gamma stays bounded but erroneous pressurefields continue to occur.
2.
Confirmed that the erroneous p-fields always occured immediately after mesh update.
3.
Changed equation for solving hydrostatic pressure from codesnippet mentioned in this thread (solving a laplacian) into default method incorporated in interDyMFoam, which is simply pd + rho*gh.
Problem with strange p-fields persisted.
4.
Removed pd from the above equation. Now p is only hydrostatic pressure, not total pressure.
p field now remains correct. pd however persisted erroneous after mesh update.
5.
Checked U fields. Found that U field shows sharp increase 1st timestep after mesh update, even with strong relaxationfactors (0.1).
6.
Focussed on U: tried different div schemes:
-Vanleer, upwind, linear limited01. Without significant improvements

At this moment I think I tracked down the problem to U solution.
Note that all the time BC on moving patch was movingWallVelocity.

Any help and/or working examples are greatly appreciated.

Brgds,

Mark

markc August 24, 2008 08:27

Hello All, We are desperate
 
Hello All,

We are desperately in need of some expert advice / help about the above problems with moving mesh and interFoam:
-The timestep following the mesh motion, pressure fields tend to explode, possibly due to velocity explosion.
-This also occurs with very small mesh motion (taking into account motion Courant number).
-We have tried all kinds of BC's, fvSchemes, solver settings, diffusivity settings, without success.
-Solver without motion runs correctly, motion itself also runs correctly but I can still not exclude programming errors in my solver.
-Tried the new interDyMFoam from OF 1.5. The tutorials run fine but are different from my problem, where some body (ship) is being moved IN the domain, instead of the complete domain being moved.

I would be grateful if anyone could give me access to some working example. E.g. the falling cylinder from mr. Hrv seems to do everything I am trying to do. Also Eric presented some tutorial about this subject in Milan, where I was not able to attend.
Is any of this material available somewhere?

Thanks in advance,

Mark

hjasak August 27, 2008 05:29

Hi Mark, Why don't you drop
 
Hi Mark,

Why don't you drop me an E-mail and we can discuss this further. It seems the solver I am working on already does all you need.

Hrv

egp August 28, 2008 06:22

Hi Mark, Unfortunately, I h
 
Hi Mark,

Unfortunately, I have been really swamped since the Milan meeting, and haven't had much chance to read the Discussion Board, or to push out the tutorial. However, I see that Hrv has offered to help you, so you should be good to go.

Eric


All times are GMT -4. The time now is 08:12.