[pi__sec]

Dear Foamers,

As I have read countless times going over the forum, I am new to OpenFoam. A friend of mine and I have started learning OpenFoam ourselves a while ago for a project we have started. We want to simulate a small hydro turbine half submerged in a flowing river.

The case we have set up is based on the many tutorials/threads/example cases we have read and studied, and managed to get quite far in the process, starting from the very basics slowly working our way up. However, because recent developments have been slow and we could really use some help.

The case is set up to use interDyMFoam with LES turbulence modelling. Because the case is a turbine in a river it seemed obvious to us to use interDyMFoam, and as far as we have read LES would provide us with far better results in this type of scenario than RAS. The case included in this threat is a section of our setup to reduce the size of the meshing, the bounding box is also reduced. The case contains a ramp, the blades of a turbine that rotate and a rotating AMI around them.

We have long spent on the creation of the mesh. The ramp and turbine blades (again, only a section is provided) were constructed using CAESES, and then loaded in OpenFoam. The meshing itself was done using snappyHexMesh. In the case the meshing it set up in parallel for 8 cores.

Two warnings are still present in the sHM log:
- Several times: Displacement at mesh point
- At end: [...] doSnap [...] Did not succesfully snap mesh. Continuing to snap to resolve easy surfaces but the resulting mesh will not satisfy your quality constraints
I mention this because I am not fully sure this is not the cause of the problem scenario described below. A checkMesh log is also provided in the case. The rotational movement in the mesh was seperatelty tested without simulation to make sure this works as well.

In the simulation phase itself we encountered several problems such as floating point errors in different scenarios. In the case I have for now also replaced the atmosphere by a wall. In first instance we try to get a fully working simulation. In the uploaded case parallelism is not activated for the simulation run.

Current scenario: the courant number seems to steadily increase in each timestep. The adjustTimeStep option is activated. The below observations seem to be the most obvious.

• Running the case with a larger deltaT of 1e-4 or 1e-6 causes an increasing courant number and quickly ends up in causing a floating point exception.
• Running the case with a deltaT of 1e-8 causes an increasing courant number and when the max courant limit is reached ends up in deltaT reducing in size every step (to e-64 and beyond) with only minimal changes in courant number

I am assuming that if I would decrease the maxCo parameter in the larger timesteps scenarios, then deltaT would just decrease in size. When using the kEqn model k also seems to explode.

I have also been experimenting with changes to fvSchemes, fvSolutions and BCs, but so far without success. I am not sure which solver would be most suitable for this particular case.

The full case including STLs can be found at the following link as the upload limit was too low: https://drive.google.com/open?id=0B7...kR4LVdidGowOFU

Any suggestions would be greatly appreciated! Even a suggestion of where the problem would probably originate from or what best to change in the case would be a great help. Thank you in advance! If any more information is needed please let me know.

tl;dr: interDyMFoam LES simulation of turbine in river flow (no atmosphere currently, just wall), the courant number increases and either causes a crash or causes the time step to become ever smaller

[pi__sec]

Already a small update with some more information. I ran the case as a laminar simulation (for just 25ms). This time there were no stability problems or explosion of Courant number. Note that the simulation was only ran for a very short time, maybe additional problems will arrive down the line.

The simulation would stabilise at a deltaT of 0.0002 and a maximum Courant number of around 0.35 - 0.36.

The (shortened) interDyMFoam log has been attached.

[hokhay]

Hi,
I am a new foamer and I am curious of your project, as I recently start dealing with LES. The time step also bother me a lot as well. Have you find a way to get around it?

The only two ways I know to reduce courant number is to increase mesh size or reduce flow velocity, which can be explain by the courant number formula. I wonder if there is a third way.

The error message suggest that the snapping process cannot be completed because of poor mesh quality. The poor quality mesh could make the max courant number lower, so probably you sure get a good mesh before starting simulation.

Sent from my LG-H818 using CFD Online Forum mobile app

[pi__sec]

Hey hokhay,

Sorry for the late reply; I have been quite busy lately. Our project has advanced a bit since my last post on the forum, but because no replies came to my post here I did not bother to update it until now.

In this post I won't give an elaborate description of the progress we have made, but in summary: we have managed to get stable LES simulations running using different turbulence models (and also a solver for one phase simulations). The time step and crash issues have all been resolved. The next steps are determining which models provide the most accurate results for our case. So far we have been comparing the lift and drag coefficients that OpenFoam calculates with a naca model, but also here the results vary over time and depend on a range of factors. We are in the process now in running simulations on a benchmark case to compare the different configuration settings. If anyone has any experience with determining the accuracy of the lift and drag coefficients of simulations by comparing it to the naca data then feel free to share your insights. If there are any pitfalls we should avoid or certain aspects we should take into consideration we are keen to know.

To come back to your own post and remarks (in reverse order):
The snapping warning that we got in the snappyhexmesh log is still present, but does not really seem to impact the quality of the mesh in any major way. Your remark about the quality of the mesh affecting the Courant number does not seem correct to me. A poor quality mesh would increase the (maximum) Courant number rather than lowering it.

The Courant number:
Below the one-dimentional formula taken quickly from wikipedia as an easy reference. The logic extends to the multidimensional case by summation of the dimensions.

With
C: Courant number
: velocity magnitude
: time step
: length interval

Answer to your remark: in addition to increasing mesh size (increasing
) and reducing flow velocity (), a reduction of the time step () can be used to reduce the Courant number. This is usually the way the Courant number is controlled, as is the case in our simulation. By setting the maxCo to a maximum value and adjustTimeStep to yes in the controlDict, OpenFoam will adjust the time step for you in order to keep the maximum Courant number under the given limit. This is what was also happening in the case given in our original post: the time step was being reduced so much in order to preserve the value to which the maximum Courant number was set. However, if your case is set up properly (which ours was not), I think this is a great way to automatically have the time step be adjusted during the simulation. This makes sure that the time step is increased if the limit is not yet reached, thus saving computation time, and that the time step is decreased if the limit would be surpassed, thus preserving accuracy. This should also provide an answer to your first question then. I would recommend to use the automatic adjustment of the time step; this way OpenFoam does the hard work for you. Please note that the original value of the time step should not be put on an unrealistic value, otherwise it could still cause the simulation to blow up right away.

I hope that you find my reply useful and that it can help you to make progress in your own simulation. I am not an expert in CFD simulations myself though, and I am sure that in these forums there are a lot of persons with much more knowledge about the subject, so in advance my apologies if any of my insights are not completely correct.

Kind regards,
A

[sheaker]

Hello dear Aäron Thijs.

I wonder if snappyhexmesh creates arcs on the AMI or just straight edges / flat planes? If there are no arcs to perfectly fit the AMI, rotation could cause some troubles.

Have a nice day.
sheaker

[pi__sec]

Hey Sheaker,

Thank you for your reply.

Snappyhexmesh creates straight edges / flat planes on the AMI, but that is always the case I think. This could introduce a small error into the results, but depending on the settings I believe it will not impact the results too much. However, we do currently still have a water phase inaccuracy that carries over into the air phase at the location where the rotation "comes out of the water and goes into the air". This is an issue that still needs to be solved as well.

A nice day to you as well!

[JNSN]

Hi Aäron,
I had a short look at your case. First, I think you should not use LES but RAS for this problem. With LES you can of course get more accurate results but only with MUCH MORE effort. If using LES in a wrong way (mesh size) the results will be rubbish. For your application RAS should provide very good results. kOmegaSST or similar should be fine. Second, as along as you do not expect a significant influence due to the free surface, e.g. Ventilation, you should go with a single phase solver. This makes everything much simpler. Third, the mesh: for testing it is fine, but even using a RAS model it is way too coarse. You should use prism layer cells around the blades and maybe also one layer around the AMI patches. Not sure if this works well with SHM, I have not much experience with this tool. Also try to extend your domain downstream to ensure the outlet BC is far enough from the turbine. You may use stretching for this.

Best regards,
Jan

Sent from my D5803 using CFD Online Forum mobile app

[pi__sec]

Hey Jan,

First of all, thanks a lot for the feedback!

I'll reply to your remarks in order:

1. The reason we opted to use LES instead of RAS is, besides the improved accuracy (despite the higher computational cost), that we wanted to have a realistic view on the large eddy creation and transition throughout the domain. If a large eddy would be caused by a turbine wing "in the front", then would this eddy have any effect on any wing further downstream. Maybe the effect of this eddy would be minimal and have no impact at all, but we do not know. Perhaps we could also set up a RAS case to compare, although so far we did try hybrid models (SpalartAllmarasIDDES and variations; although for computational cost I guess this is a not a good alternative anyway ).
2. I fully agree that using a single phase solver makes everything a lot easier. This is also why for now we have switched to pimple for further testing. We still want to do some final two phase simulations as well, as to see the effect of this on the lift, drag and momentum coefficients. There should at least be an impact.
3. The mesh has also changed quite a bit since I first posted here. If you say the mesh is too coarse based on my original case I understand your point as it has quite huge cells in it.
   1. Do you think at least the region around the wings is fine enough? We did do some research regarding mesh size and are now also testing with meshes of varying course/fineness. Do you have a rule of thumb we could follow with regard to both RAS and LES simulations? More specifically for our case as I can imagine it can differ quite a bit depending on the situation. Currently we are trying to do determine the best settings based on a naca wing in a situation with Reynolds number equal to 50000.
   2. We could try to add some extra layers surrounding the blades and AMI patches. We did try to do more fancy meshing with SHM before, but that has been a while ago. We'll give it a try.
   3. The outlet was originally also further from the turbine, as to not have it influence the results, but for testing purposes the domain was made smaller. If you say stretching I am assuming you mean extending the domain in the horizontal x direction and then stretching the added cells progressively in the same x direction? Does it matter much if this is done progressively and the rate in which the stretching is applied?

Again thank you for the advice; it is much appreciated!

Kind regards,
Aäron

[JNSN]

Originally Posted by pi__sec

Hey Jan,

First of all, thanks a lot for the feedback!

I'll reply to your remarks in order:

1. The reason we opted to use LES instead of RAS is, besides the improved accuracy (despite the higher computational cost), that we wanted to have a realistic view on the large eddy creation and transition throughout the domain. If a large eddy would be caused by a turbine wing "in the front", then would this eddy have any effect on any wing further downstream. Maybe the effect of this eddy would be minimal and have no impact at all, but we do not know. Perhaps we could also set up a RAS case to compare, although so far we did try hybrid models (SpalartAllmarasIDDES and variations; although for computational cost I guess this is a not a good alternative anyway ).
2. I fully agree that using a single phase solver makes everything a lot easier. This is also why for now we have switched to pimple for further testing. We still want to do some final two phase simulations as well, as to see the effect of this on the lift, drag and momentum coefficients. There should at least be an impact.
3. The mesh has also changed quite a bit since I first posted here. If you say the mesh is too coarse based on my original case I understand your point as it has quite huge cells in it.
   1. Do you think at least the region around the wings is fine enough? We did do some research regarding mesh size and are now also testing with meshes of varying course/fineness. Do you have a rule of thumb we could follow with regard to both RAS and LES simulations? More specifically for our case as I can imagine it can differ quite a bit depending on the situation. Currently we are trying to do determine the best settings based on a naca wing in a situation with Reynolds number equal to 50000.
   2. We could try to add some extra layers surrounding the blades and AMI patches. We did try to do more fancy meshing with SHM before, but that has been a while ago. We'll give it a try.
   3. The outlet was originally also further from the turbine, as to not have it influence the results, but for testing purposes the domain was made smaller. If you say stretching I am assuming you mean extending the domain in the horizontal x direction and then stretching the added cells progressively in the same x direction? Does it matter much if this is done progressively and the rate in which the stretching is applied?

Again thank you for the advice; it is much appreciated!

Kind regards,
Aäron

Hi Aäron,

To check for the influence between the blades also RAS will give you valuable information. Main difference of course is that you do not resolve the anisotropic turbulent structures. But if you are interested in these, some hybrid models (DES) should do the trick. In my experience no big difference between the different models. With DES you can use still a RANS mesh around the blades, which saves a lot cells. But even with DES you should first start with RANS to determine appropriate mesh size. I can forward a nice paper from Menter where this approach is nicely described.
Regarding mesh setup such as layer etc I get back to you later. My train just arrived

Sent from my D5803 using CFD Online Forum mobile app

[pi__sec]

Hey Jan,

Alright, then maybe I should first do some RAS simulations. Please forward me the paper if you have it handy!

Thanks again! Have a nice evening!

[JNSN]

Hi,

as promised, here's a link to the paper:

http://enu.kz/repository/2011/AIAA-2011-3474.pdf

Best regards,
Jan

P.S.: As soon as I find the time, I'll try to give you some (hopefully) helpfull tips for the meshing.

[pi__sec]

Thank you very much!

I have come to the conclusion that the friend doing the project with me already added it to the list of papers to read... Sorry for the extra work. I will definitely go through it!

I look forward to your advice!

Kind regards,
Aäron

[JNSN]

Hi Aäron,

some thouhgts/suggestions for your mesh setup:

1. In order to capture the flow around the blades, you should add prism layer cells around. Cell height of cells next to the blade depends weather you use wall functions or not. With wall function the non-dimensional distance y+ can be something > 35, without y+ should be < 1. For your Re-Number it should be possible to go without wall function. To get an idea about the required cell height, you can use the tool: https://www.cfd-online.com/Tools/yplus.php
   
2. Prism layer height and number of prism layers: simply spoken, use as much as possible for both :-) layer generation is unfortunately a known limitation of snappyHexMesh, so you need to try and see, how far you can get.
3. Cell count along the chord lenght: maybe something around 30 should be a good starting point.
4. When running the simulation, you should check the y+ distribution at your blades

   postProcess -func yPlus

   if I remember correct. If required, adjust your mesh to get correct y+ values.
5. Prism layer at the AMI patch: This is not required but nice to have. In order to ensure the interpolation at the AMI works fine, try to make cells on both sides with similar size. Easiest would be just to put one or max two layers on both sides, but as said, I am not sure if having different layer settings at different patches works well in sHM.
6. For the downstream region near the outlet you can e.g. use the extrudeMesh tool to increase the cell size. Just avoid large cell volume jumps between neighbour cells.

Regarding time step size: I avoid using the dynamic time step adjustment. Often you do not get smooth forces with. With reasonable limiters you can get stable and accurate results with fixed time step size, even with max Co-Numbers far above 1. Important is reasonalbe mean Co distribution in the flow regions of interesst. A time step senitivity study should help to see the required time step size. I would think, a time step size corresponding to something between 0.25 - 1 degree rotation per time step should give reliable results.

Just out of curiosity: are you doing this for your own (hobby) or with some professional background? And do you have some experience with CFD in general?

Have a nice weekend,
Jan

[pi__sec]

Hey Jan,

Sorry for my late response; I have had a very busy past few days.

1. We did try to play around with the y+ value and the calculators before. We'll definitely take it into account when trying to optimise the mesh.
2. When trying to generate boundary layers with sHM in the past it did not always give the best results, so indeed we have to see how far we can get.
3. Alright
4. You remembered correct, although you should give it as an argument to the application otherwise it will complain it does not find the turbulence model in the database.
5. Yeah, this will require some experimenting with sHM. I remember trying to make the transitional cells finer, but generating layers on both sides of a region was not possible I believe. It is either inside, outside or distance. Maybe working with a negative distance could do the trick. Or another way could be to create a slightly larger/smaller cylinder and to work with distance from there. I'll update you if we work it out.
6. Thanks for the tip! I had in mind to use blockMesh to slowly increase the cell size towards the outlet. I'll look into both solutions. Or is using blockMesh a bad idea for some reason I overlook?

Regarding time step size: I was not aware that the dynamic time step size had a negative effect on the force calculation. A time step sensitivity study is definitely desirable then.


The next steps to be taken on my end will be:

• Do some more reading of papers
• Set up the RAS case: single phase to start with; I was thinking pimpleFoam/pimpleDyMFoam
• Optimise the mesh with all the pointers you have given me
• Afterwards I will update the thread again with how far we have come

I will satisfy your curiosity in a pm, as to not make this post even longer.

As always, your feedback is greatly appreciated!

Have a nice week,
Aäron