[pi06jl6]

Question:

Why do I get pressure oscillation (PISO) when I use a parallell staggered FSI with implicit coupling?

Background:

Using the names from icoFsiFoam to illustrate the procedure

while (time.run()) {

    while (!DONE) {

         #include "solveFluid.H"
         #include "setTraction.H"
         #include "solveSolid.H"
         #include "updateMesh.H"
         ..... convergence check (DONE) ....

    }

}

This work nicely in serial form, old code using OF1.4, currently i am using OF2.1.x and I have a parallel form that cannot, ironic, be run in serial mode.

If I restrict to one iteration ie putting DONE=true, then I get reasonable result, but I obtain pressure oscillations when DONE=residual criteria [10-20 iterations]. Note however, in both cases the mean pressure is often OK. I can see difference in implementation of movePoints() between OF version, and I have noticed different strategies from other FSI solvers regarding relative flux/meshPhi, but from the development forum, from the tickets of reported bugs i cannot notice anything that points to my problem. [I use pressure/velocity from previous time step as starting field, I don't cycle the either of these in loop, only the "boundary" is the free variable here]

I always assume the error is behind the wheel, i don't complain about the vehicle, I sincerely hope that is some people that recognize what I have missed.

I have noticed few threads pointing out that you cannot "loop" PISO within same time stamp but i cannot honestly see "why not" since in the mesh.update() one grabs the old flux. Further someone has tried change icofsiFoam by setting this into a loop and failed, but no further discussion about "why not" than that (SIMPLE) algorithm cannot be in a loop. To some sense, implicit FSI is pimpleDyMFoam with the exception that mesh.update() runs inside the loop.

This issue has been brought up to the forum few times now (strongly coupled FSI). Further since parallell codes (using OF) exist with published paper, i clearly understand that this issue has been encountered by others, so anyone have any idea what is going on?

Bruno implicated that it could be my solvers (defined in fvSolution), but although I have noticed strong sensitivity but it wont change the high-frequency of the oscillation (most often the amplitude). Further, it is the pressure that goes wrong, not velocity and general while changes in U/meshPhi/V0 is in the third digit - the pressure changes at first (sign shifts)!

fvScheme: Using linearupwind grad(U) for div(phi,U) , celllimited gauss linear 1 for grad(), and gauss linear uncorrected for displacement. Interpolation: linear. temporal:backward
fvsolution: GAMG/Smoothsolver/GaussSeidel

Note: I get significant less pressure oscillations when more diffusion is added to the schemes, which is to some sense logical.

Best Regards

Johan

[wyldckat]

Hi Johan,

I might be wrong here, but AFAIK, the (official) OpenFOAM version is not able to handle FSI simulations, or at least not in the sense of deforming a solid region as a response to forces imposed by the fluid.

You might have to switch to the variant OpenFOAM 1.6-ext or the latest FOAM-Extend 3.0, because from what I know, only these variants by the Extend Project will perform FSI simulations.
Specifically, I'm guessing that the "icoFsiFoam" you're referring to as coming from an old 1.4 version, is actually 1.4-dev, which is from which the Extend variants come from.

Beyond this, regarding the official OpenFOAM version: there were several advances in dynamic meshes, both on 2.2.x and 2.3.x, which were never implemented in 2.1.x. Therefore, you might be experiencing some problems that have already been fixed in these newer versions.

Best regards,
Bruno

[pi06jl6]

Originally Posted by wyldckat

Hi Johan,

I might be wrong here, but AFAIK, the (official) OpenFOAM version is not able to handle FSI simulations ...

Bruno, thank you for your response, FYI, i have written an own solver. First a version that works excellent in serial mode, then a second version, for parallel.

A question, i have modified the pimple and what seems to cause significant amount of oscillation is the following part in pimpleDyMFoam (OF.2.1.x),

// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);


U -= rAU*fvc::grad(p);
U.correctBoundaryConditions();
sources.correct(U);

IMHO, by moving the relative flux to the end seems more correct. Otherwise you modify the velocity and then in the beginning of the pimple in next time step you makes it absolute, which is inconsistent since the U fields are not the same.

Does it sound right or do I owe someone an apology?

Anyway, by moving thay statement (makeRelative) to after source.correct(U) removed the larger proportion of the pressure oscillation.

I still do have oscillation which originates from that fact that the norm is weaker for parellel cp to serial, and further I get proper velocity fields and mean fields but instant pressure fields are non-physical. I can see a clear correlation of the pressure frequency to the flow speed/(cell distance), the amplitude of the oscillation to the variance of the boundary velocity field by the dynamic pressure equation.

I am not finished but it seems that i am currently leaving the BUG-zone [thanks to the change in code above] and now entering the FEATURE-zone, need to take a close look in thresholds, mesh quality and so forth.

My Best Regards

Johan

[pi06jl6]

Pardon for this late update, I have many different projects,

More detailed information on my issue:

A. The pressure oscillation is non-physical, ie same oscillation pattern but different frequency. The amplitude is related to the mesh speed, hence boundary velocity. The larger boundary velocity the larger amplitude.

B. The frequency seems now more related to correction scheme rather than cell size of the mesh [maybe on both], it is a regular pattern based upon consequent steps, not time. (~5 alternative negative gradient of pressure followed by ~5 positive)

C. In the PISO, when such oscillation occurs, a characteristic feature is that for the first correction (nNonOrthogonalCorrectors), it has a small initial residual (1e-3) and then in the second correction its goes up. The mesh is tetrahedral and well optimized.

D: Several different test domain is used, current used: a 2M Cells Mesh with Re 400, Domain size 3x2x6 m, cell size range from 2cm-4cm and the cantilever is of 0.2x0.2x1m with 20 MPa as Young Modulus. The deformation is around 0,2 m.

E. Thresholds For Velocity/pressure: (final): 1e-9 (tried different combination, no change). 1e-6 for motion solver (Linfty norm).

F. In each time step, I perform subiteration in order to find fixed point between fluid/solid. Average convergence ~20 iterations. I always use the previous time step fields as starting in each subiteration, the only free variable is the displacement field.

The study clearly implicates that it is the meshPhi in a loop of same time step causes this.This since explicit technique has no pressure oscillation.

I can see in different implementations and version of OF for pimpleDyMFoam substantial changes in the code layers for mesh movement, but I cannot draw any conclusion how to proceed.

The main reason I ask is that I dont wanna change OF itself, further, i I assume this is feature in the mesh motion solver that easily can be solved but since updates of OF contain different changes around the mesh motion zones with sparse documentation, I wish to hear within this forum if someone else knows the solution.

Therefore I have two questions,

Question One:

Are mesh.update() designed only to be called once during each time step?

To the second, which concerns with my strategy to solve this, namely, the easiest way is to change OF that has verified correct results (), rather than change the few lines which seems odd (namely the pointer switch in movePoints for displacementLaplacian).

Question Two:

For which stable version of OF do there exist a working FSI solver for strongly-coupled motion?

[wyldckat]

Hi Johan,

I'm not an expert on this topic, but here's what I know:

Originally Posted by pi06jl6

Question One:

Are mesh.update() designed only to be called once during each time step?

I believe so. If we run this command:

find $FOAM_APP . -name "*.[CH]" | xargs grep "\.update("

it seems that all solvers only use it once per time iteration. If you look at the file indicated by this command:

echo $FOAM_SOLVERS/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C

you'll see how it handles the corrections to the "phi" and "U" fields!
Oh, and you might want to look at the source code for 2.3.x, since there might have been some improvements made to the source code!

Originally Posted by pi06jl6

Question Two:

For which stable version of OF do there exist a working FSI solver for strongly-coupled motion?

As I've mentioned before:

Originally Posted by wyldckat

You might have to switch to the variant OpenFOAM 1.6-ext or the latest FOAM-Extend 3.0, because from what I know, only these variants by the Extend Project will perform FSI simulations.

I'm not certain if they use strongly-coupled motion, but AFAIK the Extend variants are the current reference on the topic of FSI running on OpenFOAM technology.

Best regards,
Bruno

[pi06jl6]

Thank you again Bruno, I know what you told me earlier, by all respect, I simply hope that someone more experienced with FSI would reply ( ). Your answer fits to my own conclusion . After I have finished my current project i will get my hands dirty again, i havent done anything with the code since my first post in February. I had hoped someone else could point in the right direction .

[pi06jl6]

I have now re-open the case, i hope that i can obtain a solution to this issue.

I do want to make clear on the following: I dont think there is any bug in update(), it is rather a feature.

What i am looking for is a solution to the pressure oscillation I get when I am running mesh update in a loop without update time stamp.

This is what i know:

1. everything runs fine whenever I run the case by calling Mesh.update() only once per time step.
2. oscillation occurs when I place Mesh.update() inside the FSI loop.
3. in Mesh.update(), the movePoints is called, in that routine there is a swap in pointers which i do believe cause my problem.

I would appreciate an answer to this problem. I will now start to code and try to work on a solution that don't require rewrite in the OF package.

[pi06jl6]

Originally Posted by pi06jl6

I have now re-open the case, i hope that i can obtain a solution to this issue.

Case solved.

[wyldckat]

Hi Johan,

Originally Posted by pi06jl6

Case solved.

Any chance you can share the solution?

Best regards,
Bruno

[pi06jl6]

Originally Posted by wyldckat

Hi Johan,


Any chance you can share the solution?

Best regards,
Bruno

Ah sorry my friend, of course, there is one exact solution given by Christophe Kassiotis, he wrote a letter to me after I figured out that in his work he must have encountered the same issue so I contacted him, he suggested that by moving mesh back to previous point field in each subiteration followed by another update with the new point field, the pointers were switched back into the right position.

However, i didnt like the idea of running ALE twice. So what did I do?

I ran statistic on the mesh, figured out that due to CFL condition the time step always lead to that the actual change of the volume where so small that meshPhi hardly changed [checked that on raw data as well], so I updated [mesh.update()] only at the last FSI subcycle and for the other subcycles I moved the boundary points using same mesh topology.

So silly simple!

This need of course to be validated, this was done by calculate added mass effect for cantilever, as result, the calculated shift in time period 0.38 to 0.52, the theoretical value were 0.51, the reduced velocity is 1.94 and thus on lower branch and therefore I can conclude the result to be reasonable.

Now i am calculating the reduced velocity plots to further validate the solver.

And Bruno, thank you for helping me, I will when have time migrate the code to the latest OF version to see if there is any change, since as you pointed out, the routines has been changed. The reason I havent looked into this during last two months is my focus has been on an article I have been working with, a DMD study on twisted cantilever.

Best Regards

Johan

[bieshuxuhe]

Hi ,
I am a Newbie in openfoam and fsi .

Ah sorry my friend, of course, there is one exact solution given by Christophe Kassiotis, he wrote a letter to me after I figured out that in his work he must have encountered the same issue so I contacted him, he suggested that by moving mesh back to previous point field in each subiteration followed by another update with the new point field, the pointers were switched back into the right position.

I am interested in the red sentence ! As we know , the code " mesh.update(); " is located in the main function ( int main(int argc, char *argv[]) ) .
So I guess I should move mesh back to previous point field just before the code " mesh.update(); " in the main function .
But I am confused about that whether the point field of the mesh could be recognized in main function , and we usually operate the point field in src/fvMotionSolver/pointPatchFields for dynamic mesh .
could you give me some advice on moving mesh back to previous point field in each subiteration followed by another update with the new point field ?

Best Regards
xuhe

[pi06jl6]

Originally Posted by bieshuxuhe

... could you give me some advice on moving mesh back to previous point field in each subiteration followed by another update with the new point field ?

There is several code snippets in how to store previous fields and this question has been posted several times on this forum so i hope you excuse me if I only give you the link

http://www.cfd-online.com/Forums/ope...time-step.html
http://www.cfd-online.com/Forums/ope...elaxation.html

Perhaps you could look into icoFsiFoam, icostructFoam, i have peeked into these but i haven't used it: I have used for my FSI solver with respect to OF package the interDymFoam for my Master Thesis and for my Phd work pimpleDymFoam solver, and nothing else. Then i have used the program package DEAL.II for the structure and as AMI (Arbitrary Mesh Interface) between the mesh of the structure and fluid, my own code.

I am about to finalize the validation, a piece of advice, be very careful with the error in AMI, i have noticed high sensitivity with respect to dynamic pressure in the tolerance settings for convergence in displacement field.

[bieshuxuhe]

dear Johan Lorentzon,
thanks for your reply ! I have look into icoFsiFoam(weak coupling) and icoFsiElasticNonLinULSolidFoam(strong coupling) in FOAM-Extend-3.0 .
after I read the source code of icoFsiElasticNonLinULSolidFoam , I found that there is no need to "move mesh back to previous point field in each subiteration followed by another update with the new point field" in FOAM-Extend-3.0 . maybe that is because icoFsiElasticNonLinULSolidFoam is in serial form ?
I am creating my FSI solver with respect to interDyMFoam(which looks like in serial form) , and FEM for solid .
when I just run openfoam in serial form with (official) OpenFOAM version , do I need to "move mesh back to previous point field in each subiteration followed by another update with the new point field" ???

[pi06jl6]

Originally Posted by bieshuxuhe

do I need to "move mesh back to previous point field in each subiteration followed by another update with the new point field" ???

First, this solver you are mentioning is in the extended version of OF, i don't use it and there is no such support for that version on the HPCC I work against. Secondly, the switch in pointer lies in the implementation of the motion solver, the code of update is correct. You specify the motion solver in your case settings and I am using pointDisplacement. I hope this helps you further. If I remember it correctly, they are using another motion solver, tetDecompositionMotionSolver. OF is not easy to read even for programmers like me, you must follow cookie trail careful to find the exact lines that do your bidding. Good luck with your work!

[pi06jl6]

I implemented the "exact" solution in the following fashion,

pointField prev=mesh.points()
subiterate:
solveFluid
solveSolid
mesh.movePoints(prev)
mesh.update()

Will update shortly with more information.

[bieshuxuhe]

F. In each time step, I perform subiteration in order to find fixed point between fluid/solid. Average convergence ~20 iterations. I always use the previous time step fields as starting in each subiteration, the only free variable is the displacement field.

dear Johan,
I have two question to consult you
1. about " use the previous time step fields as starting in each subiteration "
though there may be many subiterations in each time step , I think we will use the previous time step fields as starting in each subiteration by running "runtime++" only once in each time ！just like this

while (time.run()) {

    runtime++;

    while (!DONE) {

         #include "solveFluid.H"
         #include "setTraction.H"
         #include "solveSolid.H"
         #include "updateMesh.H"
         ..... convergence check (DONE) ....

    }

}

am I right ?

2. about "the only free variable is the displacement field"
I think maybe you mean the point displacement of the fluid mesh especially between fluid/solid !
but if we did as you told us that

pointField prev=mesh.points()
subiterate:
solveFluid
solveSolid
mesh.movePoints(prev)
mesh.update()

we will use the previous time step point fields as starting in each subiteration !
am I right ?

thanks

and well done for implementing the "exact" solution ! if everything goes well ，I will try that next month . then I will talk about it .

best wishes

[pi06jl6]

Answering your Q shortly:

1. Correct! The runtime++ is prior the subcycle.

2. Correct! However, this is a matter of stability, in the end it would not affect the solution.


And regarding the "exact" solution, I am not sure whether that is the right procedure, i will return on this matter, currently i spending my time with implementation of finite strain implementation.

[pi06jl6]

Originally Posted by pi06jl6

And regarding the "exact" solution, I am not sure whether that is the right procedure, i will return on this matter, currently i spending my time with implementation of finite strain implementation.

Current status is that moving the points at subcycle without movePoints and apply only at the final subcycle iteration is the correct procedure, you can see similar implementation of FSI solvers in the extended version of OF (mentioned above). I can conclude that the suggestion with moving the pointfield twice is not working and this is now clear from the implementation of the movePoints.

For this reason, the update() routine is not meant to be used in a subcycle iteration such as strong coupling.

The only concern at the moment it works only for Euler (implicit, first order) discretization - at least for 2.1.x distribution. I have one important reservation to my conclusion, I have not changed the package itself, as I would like to, but there is no point creating an own solver when you must distribute a modified package. I would appreciate if someone with experience could provide a clarification on this matter - the correct procedure in the fsi coupling for the strong case. [choice of discretization schemes and so forth]

Edit:

A Collegaue told me to upgrade to OF 2.3 since at 2.1 and 2.2 he couldnt run pimpleDymFoam with backward. For current OF 2.1.x the backward scheme gives correct results only for sufficient small dt and low boundary velocity, for higher velocities I get pressure oscillation.

A minor detail: dont blindly use the snippets in forum. This concern is obvious but I realized that the subject above in this thread [off-topic] "how to restart from previous time step" is worthy a post of its own, since the link above is misleading. The most simple and from my point of view the only way is to duplicate definition of the fields in createFields.H and then assign whenever necessary to save the values since functions such as StorePrevIter() is often triggered in sublayers of the solver, in this case, the relaxation in the PIMPLE algorithm.

[pi06jl6]

The final itch: the small pressure oscillation for backward was a feature [input file: updated the face centers every ten iteration for purpose of interpolation, when setting it to every iteration the spike was gone] in my own implemented of sliding mesh (AMI) where the interpolation error where not visible until I switched from Euler scheme to backward and it didnt affect the total results, it merely introduced a small spike in the pressure plot. I knew it was there, feels very good to have found it, looked for that error in over a month [it was a default setting that were supposed to work together with extended version, forgot the switch, or rather, i underestimated its influence of the result ].

Validation of my solver is finished with good results.

Summation:

1. Be sure to place makeRelative in the right place since different distribution place them differently and if you blindly cuts the loop and place the solid solver in wrong order it plainly goes wrong.
2. update is not designed for subcycle as for strong coupling and this "exact" solution should be implemented, to as least in the way it was proposed, I cannot see how it is done without modification of the OF package. Instead, apply movePoints in subcycle and at the final stage apply update.
3. If sufficient small time step, it doesnt matter if you restart from previous time field in each subcycle or take currently, but in general better stability is achieved, under the same assumption for PISO. Hence, what is stable for PISO is stable fo FSI if the mass ratio is in the stable regime.
4. the restart of meshPhi is not loaded correctly, as reported in the 2.3.0 release.
5. use smoothSolver for stability and check the literature for stability of FSI regarding parameter.

Regarding GLC i have studied the code in more detail and in the sublayer and I can conclude that a recent article where it is stated that it works only for Euler is plainly wrong, just remember to activate dtPhiCorr.

Everything works fine! The engine purrs now like a kitten.

Case closed!

Edit:

Special thanks to wyldckat for the advice to switch to OF.2.3, it was although a difficult upgrade since the packages did not compile under same configuration without some bug fixes of Cmake files (very annoying) and other minor stuff regarding shared libraries.

[wuwenbo]

Johan,thanks for your work for this topic about developing a useful cod to simulate the fsi in strong coupling.
I can understand most of your conclusion,but forth. Can you discrible why you say "the restart of meshPhi is not loaded correctly" detailed.
thank you!