|
[Sponsors] |
February 1, 2019, 18:59 |
Negative Cell Volumes with FSI Solver
|
#1 |
Member
Kellis
Join Date: Mar 2017
Posts: 39
Rep Power: 9 |
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: 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 |
|
February 1, 2019, 19:08 |
|
#2 |
Member
Kellis
Join Date: Mar 2017
Posts: 39
Rep Power: 9 |
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:
Thank you, Kellis |
|
February 1, 2019, 20:33 |
|
#3 |
Retired Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,981
Blog Entries: 45
Rep Power: 128 |
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:
However, there are two interesting tutorial cases in the folder "mesh/moveDynamicMesh" in foam-extend 3.2:
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): Code:
dynamicFvMesh dynamicMotionSolverFvMesh; solver velocityLaplacian; diffusivity quadratic inverseDistance 2(plate cylinder); 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: Code:
diffusivity quadratic inverseDistance 2 (FSI); Code:
diffusivity quadratic inverseDistance 1(FSI); In the example case I mentioned: Code:
2(plate cylinder); 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
__________________
|
|
February 4, 2019, 16:27 |
|
#4 |
Member
Kellis
Join Date: Mar 2017
Posts: 39
Rep Power: 9 |
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: Code:
--> FOAM FATAL ERROR: Problem with mesh motion solver selection From function initFoam in file moveFluidMesh.H at line 54. FOAM aborting 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 |
|
February 4, 2019, 18:06 |
|
#5 | |
Retired Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,981
Blog Entries: 45
Rep Power: 128 |
Quick'ish answer: This is a lot clearer! And very interesting!
Quote:
OK, I have a few... several questions:
|
||
February 4, 2019, 18:57 |
|
#6 | ||||||
Member
Kellis
Join Date: Mar 2017
Posts: 39
Rep Power: 9 |
Quote:
Quote:
Quote:
Quote:
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. Quote:
Code:
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 ); Quote:
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 |
|||||||
February 5, 2019, 17:51 |
|
#7 | |||
Retired Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,981
Blog Entries: 45
Rep Power: 128 |
Quick answers:
Quote:
If not, then that might explain why things look weird... Quote:
0.01 GPa is on the low end of rubbers, according to the Engineering Toolbox website... Quote:
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... Uhm... actually, one of the problems is me having the time to ask But I'm certainly curious here, hence the somewhat frequent feedback... |
||||
February 5, 2019, 19:50 |
|
#8 | ||||
Member
Kellis
Join Date: Mar 2017
Posts: 39
Rep Power: 9 |
Quote:
Quote:
Quote:
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: Code:
// 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; Code:
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 ); Quote:
|
|||||
February 6, 2019, 19:54 |
|
#9 | |
Retired Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,981
Blog Entries: 45
Rep Power: 128 |
Quote:
Code:
OmegaNow = OMEGA*runTime.value(); "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 |
||
February 19, 2019, 19:32 |
|
#10 |
Member
Kellis
Join Date: Mar 2017
Posts: 39
Rep Power: 9 |
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 |
|
Tags |
diffusivity, dynamic, dynamicmesh, fsi, mesh |
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
[blockMesh] blockMesh error - Negative Volume Block | adoledin | OpenFOAM Meshing & Mesh Conversion | 2 | June 22, 2016 11:44 |
negative volumes | farhadkhari | Fidelity CFD | 0 | May 15, 2016 13:38 |
thobois class engineTopoChangerMesh error | Peter_600 | OpenFOAM | 4 | August 2, 2014 10:52 |
having negative volumes in FSI simulation | sasankheirandish | ANSYS | 0 | August 2, 2014 09:17 |
Cells with t below lower limit | Purushothama | Siemens | 2 | May 31, 2010 22:58 |