[Kellis]

Good afternoon everyone,

I am currently working on a project simulating a type of turbine, using Foam Extend v3.2, and am hoping to use an FSI solver to model flexible blades. As the "first step" of the solution procedure, we have created an initialization solver which accounts for the centrifugal force acting on the flexible portion of the blade, and moves the solid and fluid meshes accordingly. Below are some pictures to help give you an idea of what's going on.

blade_overview.png

Above: a picture of the blade, attached to the hub (blue). The front half (grey) is rigid, while the trailing half (red) is flexible.

blade_pre.png

Above: this picture is a slice of the domain through the blade, before running the solver. We are solving 1/8th of an annulus with cyclic BC's on the edges. You can see there is a small tip gap (5% of chord length) between the edge of the blade and the stator.

blade_post.png

Above: this is what the blade looks like after running the initialization solver. As expected, the trailing edge deforms outward slightly, decreasing the tip gap.

However, this deformation has caused a big headache, as it collapses the cells in the tip gap area, forming negative volumes and other abnormalities. I have tried both the FVM (velocityLaplacian) and FEM (laplace) solvers, each with an array of diffusivity methods, none of which seem to cure the issue. Additionally, I have tried increasing the mesh resolution in the tip gap area, which only seems to make the problem worse.

Another thread has been posted recently, with a similar issue, but after trying the suggestions there I am still having errors:

https://www.cfd-online.com/Forums/op...imulation.html

Below are some screenshots of the problem area. I've included the ones which I think are the most important, but I can include more if needed. The screenshots are of the fluid mesh in the tip gap area (top right of blade in the two slice pictures above).

VeloLaplacianQuadInvFSIFrozDiffOff.png

The dynamicMeshDict for that one is as follows:

dynamicFvMesh   dynamicMotionSolverFvMesh;
twoDMotion      no;

solver          velocityLaplacian;
diffusivity     quadratic inverseDistance 2 (FSI);

frozenDiffusion no;

distancePatches
(
    FSI // the trailing edge patch
);

You can see some of the cells are pushed out of the domain entirely, even though the stator (outer patch) isn't supposed to be moving.

Next I tried using a directional diffusivity, still with the velocityLaplacian solver:

VeloLaplacianDirecV1.png

While the mesh looks *much* better here, unfortunately a quick checkMesh still returned a number of negative volume cells. The diffusivity vector here was (1 10 1), with the y-direction pointing radially outward in this case.

Additionally I tried a number of things with the laplace motion solver, but I will have to post a reply with the pictures from those, as I've reached the max number of attachments on this post.

Any help is much appreciated!

Thank you,
Kellis

[Kellis]

As I mentioned above, I also tried using the laplace solver, but reached the max number of attachments to I had to make another post.

First, I tried the laplace with the standard quadratic inverseDistance 2 diffusivity method:

LaplaceQuadInvFSIFrozDiffOff.png

There was some improvement over the FVM method, as the cells didn't move out past the edges, but they still bunched up severely.

Next, per a recommendation I read on the thread linked in the first post, I increased the number of cells in the tip gap, using the same diffusivity:

LaplaceQuadInvFSIHiRes.png

I'm not really sure how to describe that one, but it looks pretty ugly

Additionally, I tried using the deformation/distortionEnergy methods. Below is an example using distortionEnergy with the diffusivity exponent set to 3.

LaplaceDistortion3.png

This looks pretty good except the area right next to the blade.

Other than what I mentioned above, I also tried:

• Turning frozenDiffusivity on/off
• Making and using a quartic diffusivity method
• patchEnhanced method
• Many other things...

I am starting to think there just isn't a good way to make this work. If anyone has dealt with something similar before, please let me know!

Thank you,
Kellis

[wyldckat]

Greetings Kellis,

I've come looking for posts of yours after the PM you sent me earlier today and found this thread of yours.
My apologies if my post seems like I'm an old geezer complaining, but I'm probably too tired right now to properly understand the geometry/meshes, hence the seemly complaining post below

----

Unfortunately I'm unable to make heads or tails of what should move to where in your mesh, so I'm guessing that other forum members will have a hard time understanding it, unless they've dealt with very similar geometries...

Because I can see the wing-like shape on the first image, but then in the two images following that it is not so clear as to where that section cut really is and then the meshes seem completely different from the first 3 images...

Then there is the mesh itself, which seems a bit odd: why does it look like they were hexahedral cells that were cut into tetrahedral cells?
Is it really a tetrahedral mesh, or is it a rendering issue... see here for what I mean, namely the options to render in polyhedral mode, instead of decomposing to tetrahedrals: https://openfoamwiki.net/index.php/F...is_in_ParaView


I used the forum's image viewer to try and compare how each mesh moves, but it seems like the mesh is shown in an inverted orientation, when compared with the 2nd and 3rd images... assuming that we are seeing the blue region in meshed form... what I mean is that there is no clear reason as to why the mesh would move outside of the domain like that... it looks like it implodes/explodes, as if both top and bottom surfaces were held stiff in place?



I do have an idea on what you can do to better study how each dynamic mesh algorithm works and so that you can run it faster than using the FSI toolkit.... er, well, my idea went down the drain, because the "SnakeRiverCanyon" tutorial case that is in OpenFOAM, is not in foam-extend
If you are willing to try it out, use in almost any version of OpenFOAM the tutorial case "mesh/moveDynamicMesh/SnakeRiverCanyon" and you might see what I mean:

1. There is only one moving patch, that works in function of an STL surface that is moving against the mesh (or vice versa).
2. It's relatively fast to run and gives a step-by-step output of the mesh, to see it in ParaView.
3. It allows you full control of moving speed and the dynamic meshing settings, to seewhat does what exactly and in function of what.

However, there are two interesting tutorial cases in the folder "mesh/moveDynamicMesh" in foam-extend 3.2:

• circCylinder3d
• circCylinder3dHex

Both use the Mesquite motion solver, which allows for remeshing the mesh, when it gets too distorted, while it's running. I have no idea if it will work with the FSI toolkit...


What I mean by all of this is that in the few cases I helped out in using the FSI toolkit, using these settings (example case "run/fsiFoam/HronTurekFsi3/fluid" from the FluidStructureInteraction toolkit):

dynamicFvMesh dynamicMotionSolverFvMesh;

solver velocityLaplacian;

diffusivity quadratic inverseDistance 2(plate cylinder);

mostly worked without problems... although care was needed when meshing, given that inconsistent mesh resolution would lead to sudden distortions to the mesh... or whenever the wrong time step or solid properties were given (i.e whenever the plate would bend way too much and too fast).
Oh, and only hexahedral cells were used in the cases I assisted... this to say that tetrahedral cells can potentially be more constraining in some situations, at least in the mesh you're showing, given that it's bounding the diagonals of the hexahedral cells...


Wait, wait wait... hold on... this is wrong:

diffusivity     quadratic inverseDistance 2 (FSI);

it should be this:

diffusivity     quadratic inverseDistance 1(FSI);

no space before the parenthesis... and the number refers to the number of elements inside the list, i.e. the number of entries inside the parenthesis.
In the example case I mentioned:

2(plate cylinder);

is because there are two elements...

So, perhaps, somehow, the "2" in your case was actually applying a power of 2 to all distortions, hence the weird jumps outside of bounds...?

Best regards,
Bruno

[Kellis]

Bruno,

Thank you for taking the time to type out that response! I'm sorry that the pictures weren't especially illustrative of the issue at hand. I've taken a few more to hopefully help give you a better idea of what's going on.

First, here's an overview of the domain. You can see it's 1/8th of an annulus, with the blade showing up in the clear section in the middle:

extendeddomainpic.png

Next, here's a view of the blade (red/grey) and the slice from the pictures in my previous post (light blue). I've redone the pictures using the "Extract Cells by Region" filter, to avoid the rendering issue where the hex mesh appears to be a tet mesh.

BladeWithSlice.png

Looking at this same blade and slice from above, now with the cells shown. You were correct to assume it was a hexahedral mesh, just with a rendering issue from paraView.

BladeWithSliceMeshPre.jpg

Next, the same area but after the solver has been run. You can see the trailing edge of the blade has moved outward, causing the cells to deform. The outer edge of the domain shouldn't move at all, but some of the cells appear to be pushed past it.

BladeWithSliceMeshPost.jpg

And, finally, a close-up of the same situation as described above.

BladeWithSliceMeshPostZoom.jpg

I did fix the syntax issue with the "1(FSI)" entry in the dynamicMeshDict, but altering it didn't seem to have an effect.

Additionally, I tried using the Mesquite motion solver, but I don't think it's compatible with the FSI cases, as it spits out the following error when it tries to solve the mesh motion:

--> FOAM FATAL ERROR: 
Problem with mesh motion solver selection

    From function initFoam
    in file moveFluidMesh.H at line 54.

FOAM aborting

Looking at the moveFluidMesh.H file, it appears to be checking for either a fvMotionSolver or feMotionSolver, and apparently the Mesquite is neither.

Assuming the problem can be solved with an existing FSI mesh solver, perhaps there is just too much mesh motion happening at once? I have tried ramping up the angular velocity of the solver slowly from 0, but this didn't help. Maybe starting with a high Young's Modulus and slowly lowering it to the desired level would help?

Again, thanks for your reply! Looking forward to any more ideas you or anyone else may have.

Thanks,
Kellis

[wyldckat]

Quick'ish answer: This is a lot clearer! And very interesting!

Originally Posted by Kellis

Next, the same area but after the solver has been run. You can see the trailing edge of the blade has moved outward, causing the cells to deform. The outer edge of the domain shouldn't move at all, but some of the cells appear to be pushed past it.

Attachment 68162

Holy... why... why did it not compress the mesh? It's as if it is abiding to a minimum cell thickness, which resulted in not compensating for the mesh distortion properly and resulting in the cells staying at the minimum size and pushing everything else outward...

OK, I have a few... several questions:

1. Does the outer wall really have to be so close to the blade?
2. Is the cell thickness between the blade and the outer wall have to be so small?
3. Why is the blade deforming so much?
   1. Is it made of some kind of plastic?
   2. But even worse, what forces resulted in such a distortion?
4. What is the time step that you have in your simulation?
   1. And if it's fixed time step or self-adjusting?
   2. If it's fixed step, what is the Courant number at as the simulation progresses?
5. What are the definition that you have in the "fvSchemes" files? I'm mostly concerned about the "ddt" schemes, but the others are also suspect.
6. While we are at it, the settings you have in "fvSolution" is also important.

Suggestions:

1. Try running the fluid flow without solid deformation and calculate the forces being imposed on the blade. Or at least check the pressure values... perhaps the forces are off the scale and that would explain the massive deformations.
2. Triple check if you have the correct dimensions (in meters) on the resulting mesh.
3. Triple check if you have the correct dimensions for the other solid properties, because... something feels off the scale... as if the fluid or the blades are made of a gummy substance...
4. Write more time steps to disk, to see when it starts messing up the mesh distortion, because it is possible that there could have occurred a non-physical burst of pressure somewhere that ruptured the simulation and the subsequent damage.
   • And look at the min-max pressure values and the velocity values... and where they are located... something doesn't look right to me...

[Kellis]

Originally Posted by wyldckat

Holy... why... why did it not compress the mesh? It's as if it is abiding to a minimum cell thickness, which resulted in not compensating for the mesh distortion properly and resulting in the cells staying at the minimum size and pushing everything else outward...

That's what it looks like to me as well. I've been trying to poke around the quadratic motionDiff source code to see if it keeps the first layer size constant, but I'm not good enough with the code to figure much out.

Does the outer wall really have to be so close to the blade

Unfortunately, yes. The turbine we are trying to simulate (Wells turbine) has a tip gap between the blade and stator which the performance is heavily dependent on. The current case is 5% of the blade's chord length - I plan to look at 10% as well, which will hopefully be a bit more forgiving.

Is the cell thickness between the blade and the outer wall have to be so small?

Yes - we need a few layers in the tip gap area to get good results from the fluid mesh. At least 5, but more is preferable according to the grid size tests we've run.

Why is the blade deforming so much? Is it made of some kind of plastic? But even worse, what forces resulted in such a distortion?

Confession time: the current case is not a "true" FSI simulation. We've commented out the fluid solver section of the code, and are simply trying to solve the case for the turbine rotating "in vacuum" - i.e. the blade is just experiencing a centrifugal force, no traction from pressure or viscous forces. However, the FSI structure is maintained and the fluid and solid meshes are moved after the solid is solved.

We're hoping this initialization step will shorten the time it takes the case to reach a quasi-steady state once we add in the fluid solver as well.

For what it's worth, the Young's Modulus in the test cases is 1E+07. This was chosen somewhat arbitrarily to give a decent amount of deflection, but not so much as to contact the wall. We want a good amount of "flex" in the blade for what we're doing.

What is the time step that you have in your simulation?

Currently the ddtScheme is set to steadyState. We have retained the original updated lagrangian framework of the stock FSI solver, but our centrifugal body force term is total lagrangian (I think):

volVectorField FCentSolid = -(OMEGA ^ (OMEGA ^ stressMesh.C()))

...

        fvVectorMatrix DUEqn
            (
                fvm::d2dt2(rho,DU)
                ==
                fvm::laplacian(2*muf + lambdaf, DU, "laplacian(DDU,DU)")
                + divDSigmaExp
                + divDSigmaLargeStrainExp
                + rho*FCentSolid
              );

So I think we are restricted to solving the solid with a single timestep or else the centrifugal force will be re-added with each new time. This could be solved by using the UL form of the body force but I'm not sure what that looks like. The other option is recompiling the initialization solver to use a TL framework, but that is a pain since all of the FSI stuff relies on the DU variable.

What are the definition that you have in the "fvSchemes" files? I'm mostly concerned about the "ddt" schemes, but the others are also suspect. While we are at it, the settings you have in "fvSolution" is also important.

fvSchemes and fvSolution are attached.

I will try your suggestions #2 and 3. For suggestions #1 and 4 - It's my fault for not giving you the full story, but the fluid velocity and pressure are both 0 everywhere since we're just solving the solid side of things to get started.

Again, thank you for the help! This case is pretty convoluted, so if I need to explain anything more clearly, just ask!

Cheers,
Kellis

[wyldckat]

Quick answers:

Originally Posted by Kellis

Confession time: the current case is not a "true" FSI simulation. We've commented out the fluid solver section of the code, and are simply trying to solve the case for the turbine rotating "in vacuum" - i.e. the blade is just experiencing a centrifugal force, no traction from pressure or viscous forces. However, the FSI structure is maintained and the fluid and solid meshes are moved after the solid is solved.

We're hoping this initialization step will shorten the time it takes the case to reach a quasi-steady state once we add in the fluid solver as well.

Mmm... that's a fairly dangerous assumption. Did you test this with a unit test case with which you checked that things are working as intended?

If not, then that might explain why things look weird...

Originally Posted by Kellis

For what it's worth, the Young's Modulus in the test cases is 1E+07. This was chosen somewhat arbitrarily to give a decent amount of deflection, but not so much as to contact the wall. We want a good amount of "flex" in the blade for what we're doing.

Uh... so it's made of rubber...?
0.01 GPa is on the low end of rubbers, according to the Engineering Toolbox website...

Originally Posted by Kellis

Currently the ddtScheme is set to steadyState. We have retained the original updated lagrangian framework of the stock FSI solver, but our centrifugal body force term is total lagrangian (I think):

Oh, oh no no no... OK, I don't have proof, but I have very strong doubts you can use "steadyState" solvers for this formulation of FSI, due to how the equations are formulated. Specially if it's on the d2t2 scale...

Well, if it does happen to work, then the question is how many time steps (should have deltaT = 1.0) does it take to reach this mesh distortion? Does it do it in a single time step?

Because if the distortion happens in a single time step, then that is because most mesh distortion algorithms were not designed to do the complete distortion in a single time step. Even snappyHexMesh uses several iterations when smoothing and distorting the mesh during the snapping phase.

Furthermore, I don't remember ever seeing a solver with 'DyM' in the name that worked in steady-state...

Originally Posted by Kellis

Again, thank you for the help! This case is pretty convoluted, so if I need to explain anything more clearly, just ask!

Uhm... actually, one of the problems is me having the time to ask But I'm certainly curious here, hence the somewhat frequent feedback...

[Kellis]

Originally Posted by wyldckat

Mmm... that's a fairly dangerous assumption. Did you test this with a unit test case with which you checked that things are working as intended?

If not, then that might explain why things look weird...

I think the solid motion itself seems pretty expected - the blade flexes outward, with highest DU at the tip. One of the things on my to-do list is to create a "SRFElasticSolidFoam" solver to verify the displacements the FSI solver is giving.

Uh... so it's made of rubber...?

More or less - we test with a lot of silicone materials in the lab, so that's more or less what we're trying to simulate as well.

Oh, oh no no no... OK, I don't have proof, but I have very strong doubts you can use "steadyState" solvers for this formulation of FSI, due to how the equations are formulated. Specially if it's on the d2t2 scale...

Well, if it does happen to work, then the question is how many time steps (should have deltaT = 1.0) does it take to reach this mesh distortion? Does it do it in a single time step?

Yes, the pictures above are from a single timestep of displacement... from 0 to 1s, but set to steadyState. I have suspected for a while that it might be better to do in increments, so after hearing this from you, I've started to work on a transient version of the solver.

To keep things in the "UL" frame, by my reckoning, I need to divide the total desired angular velocity (say, 100 rad/s), by the number of timesteps. So, for 100 rad/s in 10 timesteps (0 to 1s, with a deltaT of 0.1s), adding the centrifugal force associated with 10 rad/s should work, correct?

Some example code:

In the createFields, create the incremental omega:

    // Total desired angular velocity
    scalar OM = readScalar(transportProperties.lookup("omega"));
    
    // Subdivide the omega into however many timesteps you have
    scalar tEnd = readScalar(controlDictionary.lookup("endTime"));
    scalar deltaT = readScalar(controlDictionary.lookup("deltaT"));
    
    scalar nSteps = tEnd/deltaT;

    OM = OM/nSteps;
    
    // Axis of rotation (vector).  Must already be normalized.
    dimensionedVector OMEGA (transportProperties.lookup("omegaaxis"));

    // "Incremental" angular velocity vector
     OMEGA = OM*OMEGA;

Then, in the solveSolid file, the same as the post above:

volVectorField FCentSolid = -(OMEGA ^ (OMEGA ^ stressMesh.C()));

        fvVectorMatrix DUEqn
            (
                fvm::d2dt2(rho,DU)
                ==
                fvm::laplacian(2*muf + lambdaf, DU, "laplacian(DDU,DU)")
                + divDSigmaExp
                + divDSigmaLargeStrainExp
                + rho*FCentSolid
                );

I think these will add up to give the right "total" angular velocity, no?

Uhm... actually, one of the problems is me having the time to ask But I'm certainly curious here, hence the somewhat frequent feedback...

Maybe I should prepend every question with a misguided PM to attract your attention then? Kidding - thank you for your help!

[wyldckat]

Originally Posted by Kellis

I think these will add up to give the right "total" angular velocity, no?

Mmm... the pseudo code doesn't add up properly... wait, OK, all that is missing is using the current time value to act as a multiplier factor into "FCentSolid", or have an intermediate variable, e.g.:

OmegaNow = OMEGA*runTime.value();

This is assuming that the simulation stops when it reaches the end time... although it could miss the max speed by 1 time step...

Originally Posted by Kellis

Maybe I should prepend every question with a misguided PM to attract your attention then?

"Please leave a message after the beep" came to mind... this was because I do have a template message for PMs that I received and don't have time to look into them in the near future

[Kellis]

Bruno (and anyone else who may be watching),

After some trial and error I've managed to resolve the issue with the mesh in the tip gap area (thankfully) by breaking the FCent up into 10 equal increments. Pictures of the mesh are below.

Overview of slice through the blade before running the solver, same view as previous post. Fluid mesh is blue, solid is brown.
overview.png

Tip gap mesh before (pink circled area above)
tipgappre.png

Tip gap mesh after:
tipgappost.png

It has worked wonderfully! I've also graded the cells to be smaller towards the blade, which helped with the apparent "minimum cell size" issue, and slightly increased the Young's Modulus.

These results were with the "quadratic inverseDistance 1(FSI)" diffusivity.

However, there is now an issue with the mesh on the blade where the FSI and non-FSI patches meet. Here is a zoomed in picture of the red circled area above, after running the solver:
interfacepost.png

It appears that the point shared by the FSI and non-FSI patches is experiencing a large displacement, similar in magnitude to the point to its right. However, the solid mesh doesn't move at all at this point. Is it possible to enforce the motionU to be zero here? Perhaps that point is "owned" by the FSI patch, rather than the blade (non-FSI) patch?

There are also some negative volume / wrong oriented cells on the surface of the blade at the FSI/nonFSI patch interface, as you can see below. The FSI patch is red, non-FSI is grey, and the negative/zero volume cells are bright green. The slice from above has been included for reference:
bigoverview.png

I am hoping fixing the point motion issue will resolve this as well.

As always, I appreciate any insight!

Thanks,
Kellis