[philippose]

Hello,
A wonderful Sunday to everyone!

I was trying a simple moving mesh case which involves a cylindrical tube which opens out to a larger cylindrical tube. At the point where the diameter changes, I have a conical object which can move to increase or decrease the effective area of the orifice.

The primary fluid flow direction is in the +ve y-direction, whereas the conical object moves in the -ve y-direction with a velocity of 0.075 m/s.

In the motionU file, the velocity for the patches which define the conical object is rendered with a velocity: (0 -0.075 0)

I used the "movingWallVelocity" boundary condition for the respective patches in the "U" file.

In addition, I also ran a control case with icoFoam at the initial position of the conical object.

On running the simulation, I found that the massflow (and hence the velocity) was very different between the two solvers, for a given simulation time. Here are the values I got after 0.001 seconds:

icoFoam: 0.002188 m^3/s
icoDyMFoam: 5.897e-05 m^3/s

Additional to this, I had also modified the turbFoam solver to handle mesh motion with turbulent flow. I ran the same system through both turbFoam and turbDyMFoam (The modified solver). The results after 0.001 seconds are:

turbFoam: 0.00238 m^3/s
turbDyMFoam: 3.548e-05 m^3/s

The values I expect are around the values given by icoFoam and turbFoam. The mass flow calculated by the moving mesh variants are way too small.

This difference is not due to the reduction in cross-sectional area because of the moving cone, since in 0.001 seconds, the motion is only 0.075mm, which does not give a reduction of 2 orders of magnitude.

The mesh motion itself seems to be correct as can be seen in Paraview.

Any ideas what might be wrong? I can send in the complete case if required.

Have a nice day!

Philippose

[philippose]

Hello,
I have been trying to get the turbFoam working with moving mesh, but so far, I keep getting flow which is much smaller than expected.

I was wondering.... for cases with mesh motion, the phi is made relative with:

phi -= meshPhi(U)

At the point when the function "runtime.write()" is called, the "phi" is in the relative state. Is this right, or should I convert the "phi" to absolute just before calling "runtime.write()" and convert back to relative?

Currently, before solving the U equation right at the beginning of the runtime loop, I have:

phi += meshPhi(U);
mesh.update();
phi -= meshPhi(U);

and then, within the PISO loop, after the pressure equation, I convert the phi to relative again.

Any suggestions?

Have a nice day!

Philippose

[ztukovic]

Hi,

It is wrong to use phi += meshPhi(U) in order to make flux
relative, because function meshPhi(U) do not return mesh motion flux
from the previous time step which is only appropriate for making flux
relative. I propose following alternative for icoDyMFoam:

Info<< "\nStarting time loop\n" << endl;

while (runTime.run())
{
# include "readPISOControls.H"
# include "readTimeControls.H"
# include "CourantNo.H"

# include "setDeltaT.H"

runTime++;

Info<< "Time = " << runTime.timeName() << nl << endl;

// Make the fluxes absolute
phi += meshFlux;

mesh.update();

// Update mesh motion flux
meshFlux = fvc::meshPhi(U);

// Make the fluxes relative
phi -= meshFlux;

fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
);

if (momentumPredictor)
{
solve(UEqn == -fvc::grad(p));
}

// --- PISO loop

for (int corr=0; corr<nCorr; corr++)
{
volScalarField rUA = 1.0/UEqn.A();

U = rUA*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf());

adjustPhi(phi, U, p);

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

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

if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}

# include "continuityErrs.H"

// Make the fluxes relative
phi -= meshFlux;

U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();
}

runTime.write();

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}


It is also possible to leave the flux absolute at the end of PISO loop an then it is not necessary to make the flux relative at the beginning of next time step.

Regards,
Zeljko Tukovic

[ztukovic]

Correction: phi += meshPhi(U) can not be used to make the flux ABSOLUTE (not relative) at the begining of next time step.

Regards,
Zeljko

[philippose]

Hello,
Good day!

Thank you for the pointers Zeljko... I tried to set up a turbulent solver with the guidelines you provided, but ended up with the same results I had before.

After scouring through the forum, I dont think there is anything wrong with the icoDyMFoam code, because it is used extensively by Prof.Jasak and Frank Bos for testing moving meshes. If there was a mistake it would have definitely been found and corrected.

Looking deeper into the turbFoam, icoFoam and icoDyMFoam code, I found that in icoDyMFoam, the equation in the PISO loop of icoFoam:

phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi);

has been changed to:

phi = (fvc::interpolate(U) & mesh.Sf());

in icoDyMFoam, and there seems to be an additional "correctPhi" step outside the PISO loop.

Could someone give me a short explanation as to why the "+ fvc::ddtPhiCorr" bit was left out in the moving mesh solver? And also, what does the "correctPhi" step do?

Today I downloaded the precompiled development release of OpenFOAM 1.3_11_09 from Frank Bos' website, and modified the turbFoam code to include moving meshes, but this time, I left the "+ fvc::ddtPhiCorr" bit as it is in the original turbFoam code, and did not add the "correctPhi" step from icoDyMFoam.

An additional change I made was to use "movingMeshContinuityErrs" instead of "continuityErrs", though I dont think that affects the solution as such.

It seems to be working now, though I will confirm with further tests... on the other hand... this was more of a "trial and error" solution to the problem, and it would be nice if someone could tell me what those changes in "icoDyMFoam" really imply, so that I can understand what really happens.

Have a nice weekend!

Philippose

[lr103476]

I agree on that Philippose. There are some 'strange' differences between icoFoam and icoDyMFoam. Indeed I use icoDyMFoam a lot, but up to now only to see if OpenFoam is capable of solving 3D flapping wing stuff. For the record, I did not validate the icoDyMFoam solver.

Since Zeljko Tukovic came with a different implementation of the icoDyMFoam solver, I think it is time to ask prof. Jasak for a stepwise description of the icoDyMFoam solver. Especially the relative flux thing and the movingMeshContinuityErrs piece is not clear enough for me.

So is there anybody who can explain the icoDyMFoam solver in a few steps?

Thanks, Frank

[ztukovic]

Hi Frank,

Here is new icoDyMFoam application from development version of OpenFOAM

icoDyMFoam.tgz

You will see that continuityErrs.H header is used instead
movingMeshContinuityErrs.H. This is because the flux phi is absolute at that
place in the code. After that phi is made relative using mesh motion flux
(meshFlux) which corresponds to the mesh displacement performed in the current
time step. The same mesh motion flux is used at the beginning of the next time
step to make the flux absolute. Before solution of the momentum equation, the flux is
made relative again but now with the meshFlux which corresponds to the mesh
displacement performed in the new (current) time step.

I hope this help.

Regards,
Zeljko Tukovic

[lr103476]

Thanks Zeljko, I begin to understand the principle. In fact I should with my cfd background:-S

Do think that this new approach will have a large effect on the solutions that I obtained using the old approach.

Regards,
Frank

[ztukovic]

It is hard to say because I don't know which application you are actually using? Can you send me your code?

Zeljko

[lr103476]

The solver I use is basically the original icoDyMFoam solver which came with release 1.3.

The time loop is as follows:

==================================================

while (runTime.run())
{
# include "readPISOControls.H"
# include "readTimeControls.H"
# include "CourantNo.H"

if (mesh.moving())
{
// Make the fluxes absolute
phi += fvc::meshPhi(U);
}

# include "setDeltaT.H"

runTime++;

Info<< "Time = " << runTime.timeName() << nl << endl;

bool meshChanged = mesh.update();

if (mesh.moving() || meshChanged)
{
# include "correctPhi.H"
}

if (mesh.moving())
{
// Make the fluxes relative
phi -= fvc::meshPhi(U);
}

# include "UEqn.H"

// --- PISO loop

for (int corr=0; corr<nCorr; corr++)
{
rUA = 1.0/UEqn.A();

U = rUA*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf());
//+ fvc::ddtPhiCorr(rUA, U, phi);

adjustPhi(phi, U, p);

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

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

if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}

# include "continuityErrs.H"

if (mesh.moving())
{
// Make the fluxes relative
phi -= fvc::meshPhi(U);
}

U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();
}

runTime.write();

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}

Info<< "End\n" << endl;

return(0);


==================================================

Regards,
Frank

[ztukovic]

Hi Frank,

The following part of the code was problematic for me:

if (mesh.moving())
{
// Make the fluxes absolute
phi += fvc::meshPhi(U);
}

but now I realize that time step shift is actually done after that place at the line runTime++. So fvc::meshPhi(U) function still returns correct mesh motion flux in the above part of code. Sorry, my mistake.

Please note that correctPhi must be done only in the case of topological changes of the mesh:

if (meshChanged)
{
# include "correctPhi.H"
}

Also be sure that dynamicFvMesh::update() function which you are using for the mesh motion returns false.


Additionally, I would like to say something about meshPhi(U) dependence on temporal discretisation scheme. I can guaranty that function ddtScheme::meshPhi() is implemented correctly for ``Euler'' and ``backward'' scheme, but testing on some mesh motion interface tracking cases shows that something is wrong with the CrankNicholson::meshPhi(). Hence, if you want to use second order accurate temporal discretisation scheme with the mesh motion, I would suggest you to use backward instead CrankNicholson scheme.

Regards,
Zeljko

[dmoroian]

Hi Zeliko,
I tried your version of the icoDyMFoam solver and I cannot make it work. It gives the following error:

--> FOAM FATAL ERROR : Incompatible size before mapping. Field size: 0 map size: 1076

From function
void MapInternalField<type,>::operator()
(
Field<type>& field,
const MeshMapper& mapper
) const
in file lnInclude/MapFvSurfaceField.H at line 77.

FOAM aborting

Aborted

The error appears regardless of what mesh I use, even for the mixer2D tutorial.
Do you have any ideea of what may cause the problem?

Dragos

[hjasak]

Just for the record: this actually works and I use it all the time and so does a number of other people.

It is likely that you are trying to use the ancient 1.3 release, which si riddled with bugs - you will need any version of my development sources later than approx April/2006.

Hrv

[lr103476]

Thats true, I am using those development sources since november now, without any problems. In fact it works very well.

The 'alternative' code from Zeljko, described in this topic, is exactly the same as the development sources, but coded slightly different. At least, as far I can tell.....

Regards, Frank

[hjasak]

Yup, Zeljko has found a very interesting problem with the consistent use of mesh fluxes and the version of icoDyMFoam you see is the correct one - you can trust Zeljko with such things :-)

Hrv

[dmoroian]

This is very strange, since I'm using a development version downloaded from the site you indicated in some of the posts (OpenFOAM), and the icoDyMFoam solver that comes with it works without complains.

Dragos

[lr103476]

Uhh, it appears that there is a problem with the 'new' code used by Zeljko. I run in the same problem as Dragos when using topology changes (mixer2D, mixer3D, movingConeTopo).

The error reads:

--> FOAM FATAL ERROR : Incompatible size before mapping

Maybe there is a problem with the 'MapFvSurfaceField.H' that we use, probably it is outdated and you are using an improved version.

For cases using only automated mesh motion, this code from Zeljko works fine for me. I did not discoverd before since I am not (yet) using topological changes.

Regards, Frank

[lr103476]

Anyone?

What's wrong with the 'new' icoDyMFoam from Zeljko when topological changes are used. At least the code uploaded in this topic does not work for mixer2D, mixer3D, movingConeTopo.

Frank

[hjasak]

Works for me...

Hrv

[lr103476]

:-) of course. I think that along with the icoDyMFoam there were some other changes to the code, maybe in MapFvSurfaceField.H as the error message says.

Could you possibly upload you current sources, or do they contain some brand new stuff which need to be published first :-)

Regards,
Frank