[markc]

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.
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]

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]

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]

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]

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]

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]

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:


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:


And the next pic the erronoues p-field at a different color scale:


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]

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]

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]

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.
cellDisplacement
U
p
pd
gamma
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]

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]

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.
shipFoam.zip
Casefolder.zip

Hope someone helps us a step further again,

Brgds,

Mark

[miliante]

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]

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]

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]

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.

shipFoam.zip
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]

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]

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]

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]

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