[philippose]

Hello,
A good Sunday to everyone!

Over the past week or so (more actually) I have been trying to write a solver for turbulent incompressible flow with mesh motion, and then extend it to include rigid body dynamics.

I based the solver on a combination of icoDyMFoam and turbFoam. While doing so, a couple of questions have risen, and it would be nice if someone could shed some light.

1. Some of the solvers (interFoam, icoFsiFoam, etc) do the conversion from absolute to relative flux only within the PISO loop, whereas in icoDyMFoam, the flux is converted to absolute before mesh motion, back to relative after the mesh motion, and again within the PISO loop.... how does one decide which is the correct approach?

2. In the line: phi = (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, U, phi); in the PISO loop

In icoDyMFoam, the latter half "+ fvc::ddtPhiCorr(rUA, U, phi)" is commented out. Is this always required for solvers with mesh motion?

3. icoDyMFoam uses an additional step "correctPhi" outside the PISO loop... is this required? Though, I havent seen it in any of the other solvers which handle moving meshes.

Basically, I currently have a turbulent solver with mesh motion which has all the code used in icoDyMFoam, but with the additional code required to handle turbulence, and the system seems to run fine. Since I am not too deep into the theory behind CFD, I dont know if there are any conceptual errors (anyone would like to take a look at the code??)

When I extend this to include rigid body dynamics, the pressure forces are calculated from the surfaces correctly, and I get a resulting velocity for the moving part of the mesh (by solving the motion equation: ma + cv + kx = F).... the only problem being, that in the first few iterations the pressure swings wildly, and that causes the calculated velocity to switch signs, and shoot up. This in turn messes up the fluid solver, and within about 8 to 10 iterations, the whole thing blows up.

I tried with very small time steps of around 1e-08 too... but there seems to be no difference. Is there anything special I need to take into account when modifying the "motionU" matrix based on the flow solution in each iteration?

Has anyone else tried coupling openFOAM solvers with rigid body dynamics? I know that Fluid structure interaction works fine with openFOAM, which is why I am surprised I havent got this relatively simple case of rigid body dynamics to work yet.... I am sure I am missing something simple... but what?

Any form of help... even wild ideas and shots in the dark are welcome :-)!

Have a nice day!

Philippose

[hjasak]

That would be because your initial field does not satisfy your equations. Thus, in the first few iterations, the pressure will be there to enforce continuity and will have little to do with the actual pressure.

You can eithe start from a potential flow solution (potentialFoam) or solve the flow on a static mesh + use the static solution as the initial guess.

You can safely throw out ddtPhiCorr - this is to do with a minor accuracy correction (div(U) or div(flux)).

So much for now...

Hrv

[philippose]

Hi Hrv,
Thanks for clearing up the ddtPhiCorr bit. So thats one thing done :-)!

I was thinking about using a converged solution as the starting point for the rigid body dynamics case. Now that you have confirmed that, I shall get back to you on the result of that soon.

The only thing that kept me back from trying that out so far is that from what I have seen, you dont seem to do that for your Fluid Structure interaction simulations. In theory, this problem would also be present in that case, or?

Shall be back with more and for more :-)!

Philippose

[msrinath80]

Hi Philippose,

I am very interested in trying out your rigid body dynamics code. Would you mind posting it here so I could fool around with it a bit.

Thanks!

[lr103476]

Hi Philippose,

Interesting stuff! I am also eager to try your solver, it would be nice if you could post it here.

Frank

[philippose]

Hi,
Its really nice to hear that people are interested in the concept I am trying out... on the other hand, now it increases pressure on me to try and get it to work :-)!

Anyway... as of now, I have a working framework for the system. Right from its conception, I have tried to make it as general as possible.

It is already dictionary based, and uses a modular approach (like the normal openFOAM solvers) with the code split into logically contiguous ".H" files which are included into the top level solver which I call "turbForceFoam".

The solver itself, is an amalgamation of icoDyMFoam and turbFoam to provide an incompressible turbulent solver capable of mesh motion.

With the rigid dynamics part, currently I solve the very basic equation:

m*x'' + c*x' + k*x = sigma(F) + m*g

The motion right now is constrained to only one dimension which can be specified by a vector in the dictionary, along with a host of other parameters:

1. Direction vector for gravity
2. Magnitude of the acceleration due to gravity
3. Mass of moving part
4. Spring stiffness
5. Spring pretension (as a distance)
6. Damping constant
7. Patches on which the forces are calculated

And in addition, the capability to switch on or off the force balance bit, and since I plan to include the effect of viscous forces (wall shear stress) too, I can individually switch on and switch off what factors are used for calculating the resultant force (currently):

1. Pressure forces
2. Wall shear forces (not yet implemented)

All these are dictionary based, and most of them are run-time modifiable.

All that said, though the code seems to be fine, I must say.... I have NOT got it working yet. When I enable the pressure forces bit, within a few iterations the simulation diverges beyond repair. It seems to be something to do with the fact that the pressure takes a couple of iterations to converge, and calculating the forces before it converges causes extremely large variations in the mesh velocity, which causes the whole simulation to spiral out of control.

Hrvoje suggested that I try to run the system on a case which has already converged for the static case.... I hope to try that tomorrow.

Are you people still interested in the code? If so, I can zip up the directory of the solver and post it here.... though I must warn you, that the code may look untidy, and amateurish :-)! Its still very much a pre-alpha version, and I am still learning true C++ programming :-)!

Let me know if there is still interest :-)!

Philippose

[msrinath80]

This is a great start. I am even more interested to try it out now. I was planning to work on this area and I think I can start by looking into your version and adding enhancements in addition to looking for bugs. I have a good number of experimental LDA/PIV results on settling spheres and I think it is an ideal way to test and enhance capabilities of your version.

Much appreciated! Please do post it

[lr103476]

Yeah, me too. Although, at this moment I am quite busy with my own problems :-), I really like to take a look at your code. Maybe I'll learn something from it and other people can help you to get it working.

Pui, what kind of experiments did you do? Plunging spheres? What are the Reynolds numbers? I am searching some validation cases for low Reynolds 3D moving objects (in general).

Regards, Frank

[msrinath80]

Hi Frank,

Currently we are studying two types of experiments. First is a simple settling sphere in a 3D box (closed); Second is the more tricky situation where we adjust the density of the fluid to match the sphere (through temperature control) so that it is neutrally buoyant at the center of the box. Then we release a heavier sphere from above. In both cases, the fluid velocity at a specific point in the box is probed using LDA (for e.g.). So we have the time-series to validate. Re varies from 1 to about 120. Of course we are at the very early stage in the experiments which are notoriously difficult to control. For the present we plan to use Prof. Derkson's PIV data for preliminary validation.

[philippose]

Hello again,
Great :-)! In that case, it would be a pleasure to post the code....it would be quite interesting to get this working :-)!

I hit upon an idea which I tried out, and it seems to help a great deal. I basically decided to put in a "bandwidth limiter" in the form of a limit on the maximum acceleration possible on the body.

My reasoning is, that under normal circumstances, you dont expect to see large accelerations on a body once the system has attained realistic values, so I have added an adjustable acceleration limit (run time modifiable via the dictionary).

So far I ran a test simulation for 0.00025 seconds before my hard disk ran out of space :-) It looked pretty much under control... I need to reduce the write interval and try again.

Anyway... here is the code... Any form of feedback (even however bad it is), is welcome :-)! And if anyone hits upon some really cool stuff, or new ideas, way of improving, bugs, etc., I would love to know, and be involved

Enjoy and thanks for the help / interest :-)!

Philippose

turbForceFoam.tar.gz

[msrinath80]

Thank you Philippose I'll keep you updated.

[philippose]

Hello,
I have a quick question.... in the solver that I am currently writing (turbForceFoam), I calculate the velocity and distance of the moving body by numerical integration at each time step within the solver.

If I stop the simulation in the middle, and restart, this history is lost, and the velocity and distance start from zero.

Is there anyway I can save two vectors (one for velocity, and the other for the distance moved in each direction from start of simulation) into a data file, which is written into the time corresponding time folder only when the other data (such as U, phi, etc...) are saved (i.e., only at every "writeInterval") ?

Saving it for every time step can be done with the OFStream class, but how do I connect this with the write interval used to write the IOObjects?

In theory, I could get the distance moved by taking a point on the relevant patch, and finding the difference between the co-ordinates at time t=0, and at time t=T (since the points file is saved at each writeInterval, and I can read the motionU file to get the velocity that was last saved, but I would like to know if the other option is possible.

Have a nice weekend!

Philippose

[gschaider]

Hi Philippose!

I think the easies way to accomplish this (saving non-field data for restarting) is this: create a IOdictionary with READ_IF_PRESENT, AUTO_WRITE (just like the IOobject for field data) and update the data (in your case the two vectors) there. With the AUTO_WRITE they will get written at the same time as all the other fields. If at a solver start the data is not found in that dictionary, then you've got a "fresh" run, otherwise (data is found) it's a restart.

[philippose]

Hello Bernhard,
Thanks for the response..... Thats very interesting :-)! I didnt consider the possibility of creating a dictionary... rather, I didnt know that was possible :-)!

I was always looking at the IOobject, but found that the IOobject seems to work only with field data as you have mentioned.

Shall look into the IOdictionary option, and let you know what happens....

Thanks once again, and have a nice day!!

Philippose

(P.S.: For the people interested in the solver itself, I completed a simulation over the weekend which showed that it is working as it is supposed to.

I have also added wall shear stresses to the force equation, but am yet to test it out... shall do check it once I have this restart problem solved :-)!)

[msrinath80]

Hi Philippose,

Could you elaborate the physical meaning for each of the following terms:

m*x'' + c*x' + k*x = sigma(F) + m*g

Thanks

PS: It would be nice if you could post the updated source along with a sample case.

[philippose]

Hello,
A Good day to you!

Its always difficult to decide how deep into the basics one needs to go when asked to explain a concept :-)!

Since you are involved in numerical solutions of engineering problems, I shall assume that you are technically quite advanced.

In the equation:

m*x'' + c*x' + k*x = sigma(F) + m*g

m (kg): Mass of the moving body

c (Ns/m): Damping constant of a virtual damping system you might want to add

k (N/m): Spring constant of a virtual spring you may want to add

g (m/s^2): acceleration due to gravity (note.. if you use an acceleration due to gravity, you also need to choose the right direction vector in which it is supposed to act with respect to your system)

sigma(F) (N): Sum of all the external forces acting on the body...ofcourse..excluding gravity since it is accounted for separately (resolved into the direction of motion of the body as specified using the "aVector" dictionary term - currently only 1-D motion is allowed). Currently the external forces are basically the pressure forces and the wall shear forces.

x'' (m/s^2): acceleration of the body = d^2(x)/dt^2

x' (m/s): Velocity of the body = d(x)/dt. This term is calculated by: x'[n] = x'[n-1] + x''[n]*(dt[n])

x (m): Distance moved by the body (with position at start of simulation taken as zero).. this is calculated by: x[n] = x[n-1] + x'[n]*(dt[n])

In the above three equations, [n] is the current time step, and [n-1] the previous time step.... basically a simple numerical integration.

As for an updated version of the solver.... An important addition I am yet to implement is the ability to stop and restart the simulation. As of now, if I stop the simulation and then restart it later, the velocity and position variables start from zero because I dont store these values on disk.

I would like to finish that before putting up the source again.... this might take up to Saturday or Sunday.

Have a nice day, and feel free to get back in case of any doubts or questions :-)!

Philippose

[msrinath80]

Thanks Philippose for that really verbose explanation. It was very useful

One last question though: When you say that currently only 1-D motion is allowed, you mean that movement in the other two directions is explicitly restricted, right? Is there any special reason for this (e.g. numerical issues etc.)?

[philippose]

Hi again,
Yet another lovely day to you...

Usually, in any such mechanical system, the more degrees of freedom you allow in your system, the more complex it becomes, and the more vulnerable it becomes to instabilities... oscillations, and divergence.

This combined with the fact that I am more or less new to computational fluid dynamics (and a new born baby when it comes to coupled solvers :-)!), I decided that I would initially begin with a simple system with only 1 DOF (with the other 5 explicitly restricted).... so that I can get more acquainted with the way the fluid solver reacts to the coupling, the various problems which arise due to the coupling, etc.

On the other hand, since I intend to use this solver initially for the analysis and simulation of hydraulic valves, I currently dont need any more than 1 DOF. A normal hydraulic spool is almost always a one DOF system.

Having said this much, I must also say, that its very simple to change the current solver to support 3 DOFs.... Just remove the "aVector" which constrains the motion, and solve the motion equation for x, y and z directions, and pass on the resulting velocities to the motionU file.

One problem I have is that I currently work with meshes generated by Netgen... pure Tet meshes... and I very often land up with an invalid mesh after a certain amount of mesh motion, due to high skew, and non-orthogonality. (Any ideas how to solve that?? An Open Source Hex Mesh generator exists??)

As usual..... questions / doubts / bugs / ideas are totally welcome :-)!

Philippose

[msrinath80]

"...Just remove the "aVector" which constrains the motion, and solve the motion equation for x, y and z directions, and pass on the resulting velocities to the motionU file..."

I will give this a try later this week. Thanks for your help

As for your skewness problem, you might want to check this out:

http://www.aero.lr.tudelft.nl/~frank/index.php?id=research/cfd/OpenFOAM/movingme sh/2D_movingmesh_wing

[msrinath80]

Hi Philippose,

Did you have any luck getting the save data after n time steps to work. If so could you post your latest version (preferably one which also calculates the fluid shear force on the rigid body)? Many thanks!