CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

Problem with icoDyMFoam

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By dmoroian

Reply
 
LinkBack Thread Tools Display Modes
Old   January 14, 2007, 16:38
Default Hello, A wonderful Sunday t
  #1
Senior Member
 
Philippose Rajan
Join Date: Mar 2009
Location: Germany
Posts: 530
Rep Power: 16
philippose will become famous soon enough
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 is offline   Reply With Quote

Old   January 18, 2007, 17:02
Default Hello, I have been trying
  #2
Senior Member
 
Philippose Rajan
Join Date: Mar 2009
Location: Germany
Posts: 530
Rep Power: 16
philippose will become famous soon enough
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
philippose is offline   Reply With Quote

Old   January 18, 2007, 18:15
Default Hi, It is wrong to use phi
  #3
New Member
 
Zeljko Tukovic
Join Date: Mar 2009
Posts: 22
Rep Power: 8
ztukovic is on a distinguished road
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 is offline   Reply With Quote

Old   January 19, 2007, 03:50
Default Correction: phi += meshPhi(U)
  #4
New Member
 
Zeljko Tukovic
Join Date: Mar 2009
Posts: 22
Rep Power: 8
ztukovic is on a distinguished road
Correction: phi += meshPhi(U) can not be used to make the flux ABSOLUTE (not relative) at the begining of next time step.

Regards,
Zeljko
ztukovic is offline   Reply With Quote

Old   January 20, 2007, 13:01
Default Hello, Good day! Thank
  #5
Senior Member
 
Philippose Rajan
Join Date: Mar 2009
Location: Germany
Posts: 530
Rep Power: 16
philippose will become famous soon enough
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
philippose is offline   Reply With Quote

Old   March 2, 2007, 12:03
Default I agree on that Philippose. Th
  #6
Senior Member
 
Frank Bos
Join Date: Mar 2009
Location: The Netherlands
Posts: 338
Rep Power: 9
lr103476 is on a distinguished road
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
__________________
Frank Bos
lr103476 is offline   Reply With Quote

Old   March 2, 2007, 15:00
Default Hi Frank, Here is new icoDy
  #7
New Member
 
Zeljko Tukovic
Join Date: Mar 2009
Posts: 22
Rep Power: 8
ztukovic is on a distinguished road
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
ztukovic is offline   Reply With Quote

Old   March 2, 2007, 16:34
Default Thanks Zeljko, I begin to unde
  #8
Senior Member
 
Frank Bos
Join Date: Mar 2009
Location: The Netherlands
Posts: 338
Rep Power: 9
lr103476 is on a distinguished road
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
__________________
Frank Bos
lr103476 is offline   Reply With Quote

Old   March 2, 2007, 16:49
Default It is hard to say because I do
  #9
New Member
 
Zeljko Tukovic
Join Date: Mar 2009
Posts: 22
Rep Power: 8
ztukovic is on a distinguished road
It is hard to say because I don't know which application you are actually using? Can you send me your code?

Zeljko
ztukovic is offline   Reply With Quote

Old   March 2, 2007, 21:40
Default The solver I use is basically
  #10
Senior Member
 
Frank Bos
Join Date: Mar 2009
Location: The Netherlands
Posts: 338
Rep Power: 9
lr103476 is on a distinguished road
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
__________________
Frank Bos
lr103476 is offline   Reply With Quote

Old   March 3, 2007, 06:39
Default Hi Frank, The following par
  #11
New Member
 
Zeljko Tukovic
Join Date: Mar 2009
Posts: 22
Rep Power: 8
ztukovic is on a distinguished road
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
ztukovic is offline   Reply With Quote

Old   March 3, 2007, 19:06
Default Hi Zeliko, I tried your versi
  #12
Senior Member
 
dmoroian's Avatar
 
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11
dmoroian is on a distinguished road
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
mm.abdollahzadeh likes this.
dmoroian is offline   Reply With Quote

Old   March 4, 2007, 16:19
Default Just for the record: this actu
  #13
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,763
Rep Power: 21
hjasak will become famous soon enough
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
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   March 4, 2007, 16:37
Default Thats true, I am using those d
  #14
Senior Member
 
Frank Bos
Join Date: Mar 2009
Location: The Netherlands
Posts: 338
Rep Power: 9
lr103476 is on a distinguished road
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
__________________
Frank Bos
lr103476 is offline   Reply With Quote

Old   March 4, 2007, 16:43
Default Yup, Zeljko has found a very i
  #15
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,763
Rep Power: 21
hjasak will become famous soon enough
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
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   March 4, 2007, 16:55
Default This is very strange, since I'
  #16
Senior Member
 
dmoroian's Avatar
 
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11
dmoroian is on a distinguished road
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
dmoroian is offline   Reply With Quote

Old   March 4, 2007, 17:54
Default Uhh, it appears that there is
  #17
Senior Member
 
Frank Bos
Join Date: Mar 2009
Location: The Netherlands
Posts: 338
Rep Power: 9
lr103476 is on a distinguished road
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
__________________
Frank Bos
lr103476 is offline   Reply With Quote

Old   March 7, 2007, 12:43
Default Anyone? What's wrong with
  #18
Senior Member
 
Frank Bos
Join Date: Mar 2009
Location: The Netherlands
Posts: 338
Rep Power: 9
lr103476 is on a distinguished road
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
__________________
Frank Bos
lr103476 is offline   Reply With Quote

Old   March 7, 2007, 12:58
Default Works for me... Hrv
  #19
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,763
Rep Power: 21
hjasak will become famous soon enough
Works for me...

Hrv
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   March 7, 2007, 13:06
Default :-) of course. I think that al
  #20
Senior Member
 
Frank Bos
Join Date: Mar 2009
Location: The Netherlands
Posts: 338
Rep Power: 9
lr103476 is on a distinguished road
:-) 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
__________________
Frank Bos
lr103476 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Mesh Problem with icoDyMFoam yuhai OpenFOAM Running, Solving & CFD 5 January 14, 2009 15:57
Mesh Problem with icoDyMFoam yuhai OpenFOAM Running, Solving & CFD 0 January 12, 2009 18:53
Problem with icoDyMFoam olivier OpenFOAM Running, Solving & CFD 13 December 19, 2008 10:03
Problem with icoDyMFoam tutorial matlie OpenFOAM Bugs 10 April 26, 2007 04:51
Problem starting with icoDyMFoam kassiotis OpenFOAM Running, Solving & CFD 1 March 12, 2007 12:12


All times are GMT -4. The time now is 06:47.