[i.sabahi]

TLDR: If you use pimpleDyMFoam with sixDoFRigidBodyMotion, will it be considered as loosely coupled FSI or closely coupled FSI?

Hi there

I've nearly done the simulations for my thesis project using OpenFOAM and I'm at the final stages, actually writing down the thesis and etc. Thesis is about the flutter analysis of a 2DoF airfoil in low Re numbers, nothing fancy. I used pimpleDyMFoam as my solver and dynamicMotionSolverFvMesh as the dynamicFvMesh class, and sixDoFRigidBodyMotion as my solver.

Now I'm getting mixed up regarding what the actual algorithm does. Mainly, I understand that the fluid part is done by solving simplified Navier-Stokes equations using the predictor-corrector loop(pimple), and the solid part solves the familiar equations of motion for a rigid airfoil, and treats it as a point of mass: at each time step the acceleration is derived and displacements after delta(t) are calculated.

But 1.Does it remember the airfoil's speed at the last iteration(v_0)?
2. Is my following understanding of the algorithm correct?

a.Start
b.Based on BCs and ICs, the PIMPLE algorithm calculates the pressure distribution
c.The forces on the airfoil surface corresponding to the pressure are calculated
d.Once the forces are obtained, the equations of motion are solved and the displacements are calculated
e.Domain is re-meshed according to the new position of the airfoil
f.Next step

Or the algorithm has some inner loops to make sure we have a consistent displacement? In other words, will the Added mass effect happen on this algorithm?

Please help me out if you can, i have to defend my thesis in 40 days

[i.sabahi]

Would anyone be kind enough to enlighten me regarding my question?

I understand it may sound dumb, but you tend to question all the basic stuff as you approach your defending time. And I think the basic stuff actually matter the most. Help me out please, it's very appreciated even if you point me to the direction of your unsure idea, or guess

Regards

[i.sabahi]

bump......

[i.sabahi]

Dear forum members

It would mean a lot to me if you could help me out. I understand that you have a limited time and answering (probably stupid) questions might not be the best usage of it, but it would surely be appreciated by me. I'm sure you've all experienced the confusion and stress that befalls each student before they defend their thesis.

Thanks in advance.

[fredicenci]

Hi.



Idk about the coding.. but I believe the added mass is intrinsic to the modification of the pressure field on the body as the mass on the dynamicMeshDict is the structural mass (weight in air).

I am trying free decay tests on a cylinder with no damping (structural). When in air, there is a small amplification of the response (dont know why yet). When in water, there is a more evident damiping, using f = 1/2pi sqrt(k/(m+ma)) the result is closer to the simulation, that means added mass is being calculated.

[i.sabahi]

Originally Posted by fredicenci

Hi.



Idk about the coding.. but I believe the added mass is intrinsic to the modification of the pressure field on the body as the mass on the dynamicMeshDict is the structural mass (weight in air).

I am trying free decay tests on a cylinder with no damping (structural). When in air, there is a small amplification of the response (dont know why yet). When in water, there is a more evident damiping, using f = 1/2pi sqrt(k/(m+ma)) the result is closer to the simulation, that means added mass is being calculated.

Hey there Fredi, thanks for your reply!

Do you use pimpleDyMFoam with rigidBodyMotion like me? Or are you using some other solver?

Allow me to express my understanding. The amplitude of the response is higher in the air, and in the water you observe more damping and lower response amplitude. Am I right? And I didn't understand the last part of your sentence, are you suggesting that the added mass effect is being calculated for, or considered in the solution procedure(an inner loop or something that calculates the "displaced" mass by the movement of airfoil, adds it to the airfoil and runs the structural equations with the old forces and the new mass, and sees if there's any change)?

I guess the difference in viscosity and density might have something to do with the extra damping in water, what do you think?

thanks again for sharing your experience
Other forum members are welcome to do the same, I'd be more than glad.

[fredicenci]

Originally Posted by i.sabahi

Hey there Fredi, thanks for your reply!

Do you use pimpleDyMFoam with rigidBodyMotion like me? Or are you using some other solver?

Allow me to express my understanding. The amplitude of the response is higher in the air, and in the water you observe more damping and lower response amplitude. Am I right? And I didn't understand the last part of your sentence, are you suggesting that the added mass effect is being calculated for, or considered in the solution procedure(an inner loop or something that calculates the "displaced" mass by the movement of airfoil, adds it to the airfoil and runs the structural equations with the old forces and the new mass, and sees if there's any change)?

I guess the difference in viscosity and density might have something to do with the extra damping in water, what do you think?

thanks again for sharing your experience
Other forum members are welcome to do the same, I'd be more than glad.

Yes, I am using pimpleDyMFoam with rigidBodyMotion. Actually, I found amplification in the decay test on air, which is not physical. The results on water seems to be better.

The term Added Mass is the modification on the pressure field and phase with velocity and acceleration. Frequently we translate this phenomena putting mass on the system. This isnt needed on pympleDyMFoam.. that is what I meant.

[aakie]

Hi,

I am curious if you managed to find all the information you wanted to have on the pimpleDyMSolver scheme (weakly vs strongly coupled, is it implicit or explicit etc.)?

Soon I will also have to defend my thesis and I am also using pimpleDyMFoam, but so far I was not able to find some reliable information on the solver.

Thanks.

BR, A

[Bloerb]

This depends. If you are using nOuterCorrectors >1, which means PIMPLE instead of PISO you are solving the system more than once per time step. Without specifing moveMeshOuterCorrectors though you are only solving the fluid side more than once. This would hence be weak coupling. Strong coupling is only achieved with nOuterCorrectors >1 and moveMeshOuterCorrectors yes; inside your fvSolution file.

[aakie]

Hi Bloerb,

Thank you for your quick response. Is the moveMeshOuterCorrectors option also available in pimpleDyMFoam within OpenFOAM-plus (OpenFOAM.com, v1712) ? Based on this link, I would think so: https://develop.openfoam.com/Develop...eDynamicMesh.C.

And would there be a way to monitor this strongly-coupled solver during run time, i.e. see something back in the log file? Because I just switched it on but I dont see an update from the structural solver in between the outer corrector (only after PIMPLE converged). Or is there another way to really see if this option works?

And if the option works, then I presume that the pimple tolerances (pressure/momentum and turb. if desired) determine the conservation pressure-momentum coupling AND the conservation at the fluid-structure interface. Or is there a way to impose a tolerance on the displacement computed by the structural solver?

Thanks a lot! Cheers

[Bloerb]

The option is available in v1712 but I just tested it and it does not work in OF1712, I have tested it with OF1806 though as well and it does what it is supposed to. So you might need to switch versions.

For monitoring you can look into the log file and you can use the sixDoF function object to log them separately to file.

To answer your second question. No, you can impose residual checks on the fluid side for the nOuterCorrectors, but there is no residual control for the solid side. It should however be fairly easy to implement this on your own.

Some example settings:

PIMPLE
{
    momentumPredictor        yes;
    correctPhi               yes;
    nOuterCorrectors         1000;
    nCorrectors              1;
    nNonOrthogonalCorrectors 1;
    pRefCell                 0;
    pRefValue                0;
    moveMeshOuterCorrectors  yes;
    consistent               yes;
    
    residualControl
    {
        "(U|p)"
        {
            relTol          0;
            tolerance       1e-6;
        }
    }  
}

relaxationFactors
{
    equations
    {
        U           0.95;
    }
}

[aakie]

Originally Posted by Bloerb

The option is available in v1712 but I just tested it and it does not work in OF1712, I have tested it with OF1806 though as well and it does what it is supposed to. So you might need to switch versions.

For monitoring you can look into the log file and you can use the sixDoF function object to log them separately to file.

To answer your second question. No, you can impose residual checks on the fluid side for the nOuterCorrectors, but there is no residual control for the solid side. It should however be fairly easy to implement this on your own.

Some example settings:

PIMPLE
{
    momentumPredictor        yes;
    correctPhi               yes;
    nOuterCorrectors         1000;
    nCorrectors              1;
    nNonOrthogonalCorrectors 1;
    pRefCell                 0;
    pRefValue                0;
    moveMeshOuterCorrectors  yes;
    consistent               yes;
    
    residualControl
    {
        "(U|p)"
        {
            relTol          0;
            tolerance       1e-6;
        }
    }  
}

relaxationFactors
{
    equations
    {
        U           0.95;
    }
}

Thanks a lot for your help, this is really helpful. Indeed, in OF1806 the moveMeshOuterCorrector seems to work (in the log file, I now see the motion status update in between the outer Correctors while this was not the case for OF1712). In addition, the sixDofState functionObject did not work for OF1712 and it does for OF1806

Would you say that a only a fluid tolerance (instead of having fluid AND structural tolerance) criteria suffices energy conservation at the fluid-structure interface? And btw, what does the setting 'consistent' mean in your Pimple Settings?

Thanks again! BR, A

[Bloerb]

consistent implies SIMPLEC algorithm for relaxed inner loops. Which allows for less relaxation in steady state simulations. If this is beneficial may depend on your problem and if you are using relaxation at all. Here i relax only U, but not in the final loop of the outerCorrectors (this would be UFinal).

The standard mesh motion solver is an explicit motion of the points using slerp interpolation. Hence no matrix system that needs a solution. This shouldn't have any dependence on time step. The solution of the rigid body on the other hand depends on your damping factor. A highly damped rigid body can be solved with basically any time step. In theory the time step for the fluid solution should be lower if Co<1. I am however uncertain if this is really the case for high frequency oscillations without damping. The best way to test this is to model a simple pendulum or something similar and test the maximal time step you can choose. Simply deactivate the fluid solution by setting nCorrectors 0; In that case you are only solving the rigid body. Once you know which time step gives you a stable solution choose that as a minimum.

[aakie]

Originally Posted by Bloerb

consistent implies SIMPLEC algorithm for relaxed inner loops. Which allows for less relaxation in steady state simulations. If this is beneficial may depend on your problem and if you are using relaxation at all. Here i relax only U, but not in the final loop of the outerCorrectors (this would be UFinal).

The standard mesh motion solver is an explicit motion of the points using slerp interpolation. Hence no matrix system that needs a solution. This shouldn't have any dependence on time step. The solution of the rigid body on the other hand depends on your damping factor. A highly damped rigid body can be solved with basically any time step. In theory the time step for the fluid solution should be lower if Co<1. I am however uncertain if this is really the case for high frequency oscillations without damping. The best way to test this is to model a simple pendulum or something similar and test the maximal time step you can choose. Simply deactivate the fluid solution by setting nCorrectors 0; In that case you are only solving the rigid body. Once you know which time step gives you a stable solution choose that as a minimum.

Thanks for the tip. I will see if the model is stable now. I am using quite a low time-step since I am running for Re > 3*10^6, which requires time-step size of dT = 0.001. I hope this time-step is sufficient (since I am only solving a 1DOF structural system with spring and damper).

And I read in one of your other posts (Flow induced 3D motion OpenFoam) that it is also possible to specify nIter in the dynamicMeshDict for the implicit structural solver. Whith nIter > 1 it should be possible to have a stable structural solver with higher a time-step size right? And do you know if nIter is a possible option within OFv1806 by any chance?

Thanks again! Super helpful

BR, Adriaan

[Bloerb]

No nIter is only available inside the rigidBodyDynamics library. Since you are however only solving one body this is already sufficiently achieved by the outer correctors. What you can do is add the following to your dynamicMeshDict:

accelerationRelaxation 1.0;
accelerationDamping 1.0;

damping the acceleration should only be done to reach a steady state. Relaxation on the other hand can be beneficial if your solid solver crashes. Make sure though that you have enough outerCorrectors so that you are not changing the solution.

[aakie]

Originally Posted by Bloerb

No nIter is only available inside the rigidBodyDynamics library. Since you are however only solving one body this is already sufficiently achieved by the outer correctors. What you can do is add the following to your dynamicMeshDict:

accelerationRelaxation 1.0;
accelerationDamping 1.0;

damping the acceleration should only be done to reach a steady state. Relaxation on the other hand can be beneficial if your solid solver crashes. Make sure though that you have enough outerCorrectors so that you are not changing the solution.

Oke that makes sense. accRelax is being used every OuterCorrection for the updated acceleration where aPrevIter is compared against the newly computed aNew. If aNew is close to aPrevIter, the relaxation should not have impact the aNew and hence the solution.

But what does aDamp do? It seems to me that aDamp basically cuts a piece off the acceleration every outer correction if I am correct (numerical damping). So this could affect the solution? (Btw, I am running transient).

And in addition, what is your opinion on symplectic (explicit) vs. Newmark (implicit)? So far Ive read that symplectic is better for energy conservation between time-steps, because its based on leap frog which is conservative. On the other hand, Newmark with multiple predictor-corrector loops should also provide energy conservation. Since there is no option in dynamicMeshDict to control the implicit iterations of Newmark/CrankNicolson, it is possible that the symplectic could be better in that case?

Sorry to bother you with all these questions, but you seem to be on point, thanks again!!

BR, Adriaan

[Magistrane]

Hello Aakie and Bloerb (again),

I join this subject because it is exactly were I am now (after the topic : Flow induced 3D motion OpenFoam). As I told, my case was running but... I am suspicious about the convergence : I have pressure bubbles that appears around my rotating mesh. So I suppose that it is because my coupling was poor :

PIMPLE
{
    correctPhi          no;
    nOuterCorrectors    2;
    nCorrectors         1;
    nNonOrthogonalCorrectors 0;
}

In DyM :

solver              sixDoFRigidBodyMotion;
accelerationRelaxation 0.3;

So I changed fvSolution with Bloerb's code (but with a relaxation factor at 0.8 because of divergence) aaaaaand I don't get how this is supposed to work : flow is now iterating 1000 times before changing time AND Pimple still writing :

PIMPLE: iteration 1000
smoothSolver:  Solving for Ux, Initial residual = 1.71972e-15, Final residual = 1.71972e-15, No Iterations 0
smoothSolver:  Solving for Uy, Initial residual = 6.10606e-17, Final residual = 6.10606e-17, No Iterations 0
smoothSolver:  Solving for Uz, Initial residual = 6.69439e-17, Final residual = 6.69439e-17, No Iterations 0
GAMG:  Solving for p, Initial residual = 3.68232e-06, Final residual = 3.68232e-06, No Iterations 0
GAMG:  Solving for p, Initial residual = 3.68232e-06, Final residual = 1.71712e-07, No Iterations 1
time step continuity errors : sum local = 5.51833e-12, global = 6.57687e-15, cumulative = -2.75267e-08
smoothSolver:  Solving for omega, Initial residual = 2.80963e-05, Final residual = 5.05917e-08, No Iterations 1
bounding omega, min: -0.000540369 max: 684.726 average: 24.866
smoothSolver:  Solving for k, Initial residual = 0.0467692, Final residual = 5.842e-08, No Iterations 3
PIMPLE: not converged within 1000 iterations
ExecutionTime = 2879.78 s  ClockTime = 2906 s

How is it not converged ? Residuals are all below 10^-6. So I don't get why it is not converge and why it didn't go to the next iteration...
Btw, I continu to follow the post that match with what I need to compute

Good evening !

[aakie]

Originally Posted by Magistrane

Hello Aakie and Bloerb (again),

I join this subject because it is exactly were I am now (after the topic : Flow induced 3D motion OpenFoam). As I told, my case was running but... I am suspicious about the convergence : I have pressure bubbles that appears around my rotating mesh. So I suppose that it is because my coupling was poor :

PIMPLE
{
    correctPhi          no;
    nOuterCorrectors    2;
    nCorrectors         1;
    nNonOrthogonalCorrectors 0;
}

In DyM :

solver              sixDoFRigidBodyMotion;
accelerationRelaxation 0.3;

So I changed fvSolution with Bloerb's code aaaaaand I don't get how this is supposed to work : flow is now iterating 1000 times before changing time AND Pimple still writing :

PIMPLE: iteration 1000
smoothSolver:  Solving for Ux, Initial residual = 1.71972e-15, Final residual = 1.71972e-15, No Iterations 0
smoothSolver:  Solving for Uy, Initial residual = 6.10606e-17, Final residual = 6.10606e-17, No Iterations 0
smoothSolver:  Solving for Uz, Initial residual = 6.69439e-17, Final residual = 6.69439e-17, No Iterations 0
GAMG:  Solving for p, Initial residual = 3.68232e-06, Final residual = 3.68232e-06, No Iterations 0
GAMG:  Solving for p, Initial residual = 3.68232e-06, Final residual = 1.71712e-07, No Iterations 1
time step continuity errors : sum local = 5.51833e-12, global = 6.57687e-15, cumulative = -2.75267e-08
smoothSolver:  Solving for omega, Initial residual = 2.80963e-05, Final residual = 5.05917e-08, No Iterations 1
bounding omega, min: -0.000540369 max: 684.726 average: 24.866
smoothSolver:  Solving for k, Initial residual = 0.0467692, Final residual = 5.842e-08, No Iterations 3
PIMPLE: not converged within 1000 iterations
ExecutionTime = 2879.78 s  ClockTime = 2906 s

How is it not converged ? Residuals are all below 10^-6. So I don't get why it is not converge and why it didn't go to the next iteration...
Btw, I continu to follow the post that match with what I need to compute

Good evening !

Hi Magistrane,

It seems im not the only one struggling with this About your 1000 iterations; it can mean a couple of things in my opinion:
- what are the tolerances of your linear solver? make sure they are tight enough so that the desired pimple tolerance of Bloerb is achieved (1e-06). Put your linear tolerances to maybe 1e-10 or higher. Your linear tolerance is concerned with the inner Corrector while the Pimple tolerance is concerned with your outerCorrector. Once the linear tolerance is achieved of your velocity for example, it moves on to the pressure and then turb. quantities etc. At one point your PIMPLE solver compares the final values of your present outer loop to the previous one and if these are close enough to each other (i.e. if your pressure-momentum coupling is correctly solved), then it moves on to the next time -step.
- other possible solutions:
- use initial flow fields to give your solve a good starting point
- check your time-step/Courant number
- check your discretization schemes
- check your mesh quality

I am sure you've checked already a couple of these suggestions Btw, these suggestions are more linked to your fluid solver, which is the problem I guess since your still stuck in your first time step ~ right? Gr. A

[Bloerb]

I have attached a simple test case for your convenience. It is a simple undamped single mass oscillator. The solid equation is hence




the analytical solution to this is





I have set nCorrectors to zero. Hence this only solves the solid. You should be able to test the performance of the different time integrators and time steps with this, as well as the influence of the outer correctors. A simple plot script is added as well so you can visualize the solution. From my experience symplectic is the most stable version. the implicit solvers can be more accurate, and allow for higher time steps for damped systems but require enough outer correctors.



Magistrane you have bounding in your fluid solution for the turbulence. This is not good and probably means that there is something wrong with it. You should add relaxation for your turbulence quantities. You are currently doing the full 1000 loops because my settings only exit the pimple loop if the residual for p and U falls below 1e-6. This is quite strict. your p residual seems to be at 3e-6 so above that. Hence you are not leaving the loop.


And yes the accelerationDamping will induce numerical damping and hence falsify your transient run. The same applies to relaxation if you do not choose enough iterations.

[aakie]

Oke thanks a lot Bloerb, I will start to play around with the NouterCorrector/deltaT! it seems indeed the symplectic is more robuust (as long as the time-step is sufficiently small).

Besides the 2-way coupled FSI model of the 6DoFSolver, I am currently running a 1-way coupled model (dynamicMotionSolverFvMesh) where the mesh diffusion was solved by Laplace's method (with inverse diffusion distance). I noticed that when I put the moveMeshOuterCorrectors parameter on, the cellDisplacement was solved for every PIMPLE outer loop, while when moveMeshOuterCorrectors was switched off, the cellDipsplacement was only solved one time per time-step (not every Pimple OUTER loop). Would you say when the laplace equation is not being solved every PIMPLE Outer loop together with the fluid, the solution might be affected? In this model, the structure is not being solved but only the dynamic mesh for every time-step (1-way coupled).

Also, I noticed in your example script, you are using correctPhi with pcorr. What does this function do and is this setting especially necessary for a moving mesh? I read somewhere that this setting might correct the flux while the mesh moves.

Thanks in advance! BR, A