# Problem with icoDyMFoam

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

 January 14, 2007, 15:38 Hello, A wonderful Sunday t #1 Senior Member   Philippose Rajan Join Date: Mar 2009 Location: Germany Posts: 552 Rep Power: 25 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

 January 18, 2007, 16:02 Hello, I have been trying #2 Senior Member   Philippose Rajan Join Date: Mar 2009 Location: Germany Posts: 552 Rep Power: 25 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

 January 18, 2007, 17:15 Hi, It is wrong to use phi #3 New Member   Zeljko Tukovic Join Date: Mar 2009 Posts: 22 Rep Power: 17 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

 January 19, 2007, 02:50 Correction: phi += meshPhi(U) #4 New Member   Zeljko Tukovic Join Date: Mar 2009 Posts: 22 Rep Power: 17 Correction: phi += meshPhi(U) can not be used to make the flux ABSOLUTE (not relative) at the begining of next time step. Regards, Zeljko

 January 20, 2007, 12:01 Hello, Good day! Thank #5 Senior Member   Philippose Rajan Join Date: Mar 2009 Location: Germany Posts: 552 Rep Power: 25 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

 March 2, 2007, 11:03 I agree on that Philippose. Th #6 Senior Member   Frank Bos Join Date: Mar 2009 Location: The Netherlands Posts: 340 Rep Power: 18 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

 March 2, 2007, 14:00 Hi Frank, Here is new icoDy #7 New Member   Zeljko Tukovic Join Date: Mar 2009 Posts: 22 Rep Power: 17 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

 March 2, 2007, 15:34 Thanks Zeljko, I begin to unde #8 Senior Member   Frank Bos Join Date: Mar 2009 Location: The Netherlands Posts: 340 Rep Power: 18 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

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

 March 2, 2007, 20:40 The solver I use is basically #10 Senior Member   Frank Bos Join Date: Mar 2009 Location: The Netherlands Posts: 340 Rep Power: 18 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

 March 3, 2007, 05:39 Hi Frank, The following par #11 New Member   Zeljko Tukovic Join Date: Mar 2009 Posts: 22 Rep Power: 17 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

 March 3, 2007, 18:06 Hi Zeliko, I tried your versi #12 Senior Member     Dragos Join Date: Mar 2009 Posts: 648 Rep Power: 20 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::operator() ( Field& 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.

 March 4, 2007, 15:19 Just for the record: this actu #13 Senior Member   Hrvoje Jasak Join Date: Mar 2009 Location: London, England Posts: 1,905 Rep Power: 33 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

 March 4, 2007, 15:37 Thats true, I am using those d #14 Senior Member   Frank Bos Join Date: Mar 2009 Location: The Netherlands Posts: 340 Rep Power: 18 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

 March 4, 2007, 15:43 Yup, Zeljko has found a very i #15 Senior Member   Hrvoje Jasak Join Date: Mar 2009 Location: London, England Posts: 1,905 Rep Power: 33 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

 March 4, 2007, 15:55 This is very strange, since I' #16 Senior Member     Dragos Join Date: Mar 2009 Posts: 648 Rep Power: 20 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

 March 4, 2007, 16:54 Uhh, it appears that there is #17 Senior Member   Frank Bos Join Date: Mar 2009 Location: The Netherlands Posts: 340 Rep Power: 18 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

 March 7, 2007, 11:43 Anyone? What's wrong with #18 Senior Member   Frank Bos Join Date: Mar 2009 Location: The Netherlands Posts: 340 Rep Power: 18 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

 March 7, 2007, 11:58 Works for me... Hrv #19 Senior Member   Hrvoje Jasak Join Date: Mar 2009 Location: London, England Posts: 1,905 Rep Power: 33 Works for me... Hrv __________________ Hrvoje Jasak Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk

 March 7, 2007, 12:06 :-) of course. I think that al #20 Senior Member   Frank Bos Join Date: Mar 2009 Location: The Netherlands Posts: 340 Rep Power: 18 :-) 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