[AMR96]

Hello,

I'm trying to simulate flow through a curved wind tunnel using the k-omega SST turbulence model. The plan is to simulate a car inside the tunnel but for now I want to just simulate the flow through the empty tunnel to make sure the boundary conditions and other parameters are working properly.

Unfortunatly, I have not been able to get a simulation with sesnsible residuals. The residuals tend to stay in the order of 10^-1. Also my y+ value explodes and increases up to 1e12 which is completely wrong. I have not been able to pin-point where the problem is. I tried the same simulation with Realizable k-epsilon and got better residuals and y+ values.

I have generated the mesh in Fluent and calculated the first layer height of the inflation layer using y+ calculators online. I also used the CFD forum calculator to get the values of k and omega. The turbulence intensity and viscosity ratio were both set to 1% and 1 respectively.

So far, I have tried adjusting the schemes to see if anything changes. I tried using first order scheme for U and second order for omega and k. I also setting everything to first order scheme but nothing changed.

The boundary conditions for U and p come from the experimental tunnel.

I would like to know where the problem is and how I can solve it. Any guidance is greatly appreciated.

Here are the Boundary condition codes:

U

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0.0 0.0 0.0);

boundaryField
{

    inlet
        {
            type            groovyBC;
			variables
			(
			"xp2=pos().x;"
			"zp2=pos().z;"
			"z0=-3;"
			"r2=sqrt(pow(xp2,2)+pow(zp2-z0,2));"

			"w=0.88;"
	
			"U_tx=w*(zp2-z0);"
			"U_tz=w*xp2;"
			"U_ax=0;"
			);
			valueExpression "vector(U_tx, U_ax, U_tz)";
			value uniform (0 0 0); 
        }

    ground
    {
        
        type		rotatingWallVelocity;
		origin		(0 0 -3);
		axis		(0 1 0);
		omega		0.88;
		value		uniform (0 0 0);
    }

    outlet
    {
        type            pressureInletOutletVelocity;
		value           $internalField;

        
    }
	
	tunnel-walls
    {
        type            noSlip;
        
    }

}

p

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0.0;

boundaryField
{

    inlet
    {
        type        zeroGradient;
    }

    ground
    {
        type        zeroGradient;
    }

    outlet
    {
        type			zeroGradient;
		type            fixedValue;
        value           uniform 0;
    }

    tunnel-walls
    {
        type            zeroGradient;
        
    }
	

}

k

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0.001;

boundaryField
{

    inlet
    {
        type            fixedValue;
        value           $internalField;
    }

    ground
    {
        type        kqRWallFunction;
        value       $internalField;
    }

    outlet
    {
        type            inletOutlet;
        inletValue      $internalField;
        value           $internalField;
    }

    tunnel-walls
    {
        type        kqRWallFunction;
        value       $internalField;
        
    }

}

Omega

dimensions      [0 0 -1 0 0 0 0];

internalField   uniform 25 ;

boundaryField
{
    
	inlet
    {
         type            fixedValue;
         value           $internalField;
    }

    outlet
    {
        type            inletOutlet;
        inletValue      $internalField;
        value           $internalField;
    }

    tunnel-walls
    {
        type            omegaWallFunction;
        value           $internalField;
    }

    ground
    {
        type            omegaWallFunction;
        value           $internalField;
    }

   
}

nut

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{

    inlet
    {
        type        calculated;
        value       $internalField;
    }

    ground
    {
        type            nutkWallFunction;
        value           uniform 0;
    }

    outlet
    {
        type        calculated;
        value       $internalField;
    }

    tunnel-walls
    {
        type            nutkWallFunction;
        value           uniform 0;
        
    }

}

fvscheme

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         leastSquares;
   
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss upwind;
    div(phi,k)      Gauss linearUpwind grad(k);
    div(phi,omega)  Gauss linearUpwind grad(omega);
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method meshWave;
}

[Alczem]

Hey


You have defined two types of condition for your outlet pressure condition, which one is it? From the rest of the setup, I guess it should be fixedValue.


For fvSchemes, if you want the most robust settings, I think you can try Gauss upwind for U, k, and omega (as you said) but also cellLimited Gauss linear 1 for the gradSchemes. You could also try Gauss linear uncorrected for the laplacianSchemes and uncorrected for the snGradSchemes (it removes the non-orthogonal correction but also improves stability if I remember correctly)



Other than that, maybe try without boundary layers first, with a coarse mesh. Are your mesh metrics ok?


Good luck!

[AMR96]

Originally Posted by Alczem

Hey


You have defined two types of condition for your outlet pressure condition, which one is it? From the rest of the setup, I guess it should be fixedValue.


For fvSchemes, if you want the most robust settings, I think you can try Gauss upwind for U, k, and omega (as you said) but also cellLimited Gauss linear 1 for the gradSchemes. You could also try Gauss linear uncorrected for the laplacianSchemes and uncorrected for the snGradSchemes (it removes the non-orthogonal correction but also improves stability if I remember correctly)



Other than that, maybe try without boundary layers first, with a coarse mesh. Are your mesh metrics ok?


Good luck!

Thank you for your reply Alczem.

You are correct, I have fixed my outlet condition and only used fixed value.

Regarding the fvScheme suggestions, would they also improve the robustness of the Realizable k-epsilon method?

The mesh quality of the empty tunnel is good but if I add the car, I will get high skewness of around 7 and non-orthogonality of around 80. That is why I prefer to have use the corrected method for gradient and laplacien schemes.

Best,

Amirmohammad

[Alczem]

Regarding the fvScheme suggestions, would they also improve the robustness of the Realizable k-epsilon method?

Yep it should, although I would stick with the kOmegaSST even with a coarser mesh. The wall laws should take care of the larger cells for a first simulation.

The mesh quality of the empty tunnel is good but if I add the car, I will get high skewness of around 7 and non-orthogonality of around 80. That is why I prefer to have use the corrected method for gradient and laplacien schemes.

A skewness of 7 is not ideal but is still ok in my experience depending on where the skewed cells are located. A non-orthogonality of 80 is more concerning in my opinion. Is it due to the boundary layers? That is why I suggest to try without boundary layers first, and make sure you have nNonOrthogonalCorrectors > 1 in your fvSolution, it might help too. But I would try to improve the mesh first if possible.


Out of curiosity, can you include a picture of a crinkle slice of your mesh near the car so we can get a proper idea of the geometry?

[AMR96]

Originally Posted by Alczem

A skewness of 7 is not ideal but is still ok in my experience depending on where the skewed cells are located. A non-orthogonality of 80 is more concerning in my opinion. Is it due to the boundary layers? That is why I suggest to try without boundary layers first, and make sure you have nNonOrthogonalCorrectors > 1 in your fvSolution, it might help too. But I would try to improve the mesh first if possible.


Out of curiosity, can you include a picture of a crinkle slice of your mesh near the car so we can get a proper idea of the geometry?

I create the mesh using Ansys Fluent Mesh. I am not exactly sure what causes the high non-orthogonality. The smallest non-orthogonality that have managed to get is 78 which is still high.

I haven't tried anything with no boundary layers yet, but I will test it to see how the orthogonal quality changes.

As for the geometry, it is a F1 car. I am trying to simulate it in a curved tunnel. I have attached some pictures. As you can see the front wing profiles are very close to each other. Unfortunately, this causes some complications with the boundary layers.

Do you have any idea how I can improve the mesh quality, specially the orthogonality aspect of it?

Best,

Amirmohammad

[Alczem]

Hey!


Aah such a satisfying mesh how did you manage to import it though? I remember trying to import a polyhex-core mesh in OpenFOAM but it would not read some cells. What does checkMesh say about your mesh?


Concerning the non-orthogonality, try to locate the bad cells before importing the mesh. I am guessing they should be near the layers, but they might also be at some rough edge or some other complex part of the geometry.


Have you tried snappyHexMesh or cfMesh by the way? sHM is a bit difficult to learn at first, but it usually handle complex geometries quite well although adding layers can be hard. cfMesh is ridiculously simple and is one of my favorite meshing tools for OpenFOAM when the simulation involves only one region of fluid and no baffles. It can also generate layers on almost anything and it is really fast. Both sHM and cfMesh require stl files of your geometry. You should check them out

[AMR96]

Hi,

It looks good from the outside, it's not so good when you run the simulation though )

I followed what was mentioned in the last message of this thread to export it in ASCII and import it in OF.

Error in polyhedral-mesh (Fluent Meshing) conversion to OpenFOAM.

Normally, everything is fine when I do the checkMesh. But sometimes the ''point usage'' check fails and it creates a set called unused points. Is this something that I should be concerned about?

Thanks for the tip, I'll look into it before importing the mesh.

I haven't worked with sHM and cfMesh. I would like to try them but my geometry has so many surfaces that I'm not sure how I can export the stl file with proper named selections.

[Alczem]

I followed what was mentioned in the last message of this thread to export it in ASCII and import it in OF.

Error in polyhedral-mesh (Fluent Meshing) conversion to OpenFOAM.

I remember reading this, I will have to give it a try then.

Normally, everything is fine when I do the checkMesh. But sometimes the ''point usage'' check fails and it creates a set called unused points. Is this something that I should be concerned about?

Not sure what this message is about sorry. Maybe there are some duplicate points that are automatically put aside by OpenFoam.

I haven't worked with sHM and cfMesh. I would like to try them but my geometry has so many surfaces that I'm not sure how I can export the stl file with proper named selections.

That is often the problem, you need good stl files. Get creative to generate your STLs! For instance, if you managed to import your mesh in OpenFOAM, it means you can open it in Paraview, which means that you can extract your patches and then save them invidually as STL files.


Quick guide if you go with cfMesh for a first try, you need one stl file that includes all your groups. So once you have exported every patch as a STL file via Paraview, you can concatenate them into one master STL file. This master file will be used by cfMesh and the mesh will include all the patches you defined. For the patch names to be properly created, the individual STL files saved through Paraview should be edited: when you open an ASCII STL file, the first line is "solid NAME_OF_THE_STL" and then the facets. You should replace NAME_OF_THE_STL by the appropriate name you need in OpenFOAM. Then you can concatenate all your STL files into one. cfMesh uses a file called meshDict and you can find templates on this forum or check out the user guide

[AMR96]

Hi Alczem,

Thank you very much for the detailed explanation, I'll look into using cfmesh.

I also wanted to know your opinion on the simulation that I am running now. I have placed the car in the tunnel and have modified the fvScheme according to your suggestions. However, I am still using realizable k-epsilon for now.

Right now I keep getting a divergent behaviour. For about 200 iterations, the residuals seem to be converging but after a while the epsilon residual starts to increase again and as a result the velocity residuals keep increasing with it until they reach a magnitude of 1e-1. The values of the lift and drag coefficient start to explode to really high magnitudes. Meanwhile the pressure residual decreases to a magnitude of about 1e-8.

Do you know what could be causing such behaviour and how I can resolve it?

Best,

Amirmohammad

[Alczem]

Hey,


You can try to start from scratch with the most simple case and add complexity step by step. For instance, try the most simple boundary conditions (no groovy or rotatingWall ) with a very slow air speed, turn off the turbulence to see if it helps for the first iterations. You could also try to replace your F1 by a simple geometry like a cylinder or a cube just to see if it is the mesh (btw did you try without layers ? ), the geometry or the numerical settings that cause the divergence. You can also try different values of relaxation but simpleFoam is a really basic and robust solver, so I would still say that the issue lies within the mesh if you are using the same fvSchemes. Last option, you could try to run a transient simulation since the flow is probably unsteady. This might be why the residuals don't go down as much, so I would rely on a stable Cl and Cd to monitor convergence.



These are the main things off the top of my head, there could be others. Good luck!

[AMR96]

Thank you very much for the tips Alcezm.

I have already tried some of your suggestions. For example, I tried my groovy BC and rotating wheels BC without any car in the tunnel and the simulations were fine.

I actually looked into the mesh and fixed some of the regions. Now I have a max non-orthogonality of 62 and the simulation is at least stable. However, the pressure residual is still not converging to lower values.

I just have one question regarding evaluating residuals in OpenFoam. Currently, I am using 2 nNonOrthogonalCorrectors in my fvSolution. As a result I run through the inner loop 3 times for the pressure, my question is that in this case, which residual should I consider as the final value. For example, here is the residuals for the latest iteration:

Time = 2123

smoothSolver: Solving for Ux, Initial residual = 0.00032114322, Final residual = 2.8402765e-05, No Iterations 3
smoothSolver: Solving for Uy, Initial residual = 0.003573726, Final residual = 0.00030546567, No Iterations 3
smoothSolver: Solving for Uz, Initial residual = 0.00097934705, Final residual = 9.1617658e-05, No Iterations 3
GAMG: Solving for p, Initial residual = 0.013409324, Final residual = 8.4451575e-05, No Iterations 4
GAMG: Solving for p, Initial residual = 0.0015543383, Final residual = 7.6809216e-06, No Iterations 4
GAMG: Solving for p, Initial residual = 0.00011836172, Final residual = 1.1426608e-06, No Iterations 4
time step continuity errors : sum local = 1.8951391e-06, global = 1.1958294e-08, cumulative = -0.00023340107
smoothSolver: Solving for epsilon, Initial residual = 3.0484591e-05, Final residual = 1.3867492e-06, No Iterations 3
bounding epsilon, min: -612.72797 max: 549252.92 average: 647.93591
smoothSolver: Solving for k, Initial residual = 1.2121451e-13, Final residual = 8.7842305e-15, No Iterations 3
bounding k, min: -0.061413721 max: 5.5413289 average: 0.13781884
ExecutionTime = 231.45 s ClockTime = 233 s

In this case is my pressure residual 0.013 (from the first run) or 0.0001 (from the last run)?

Best,

Amirmohammad

[Alczem]

Hey!


Good to hear that there have been at least some improvement!


If I am not mistaken, the residual you should consider is the first initial one, so 0.013. But seeing the other inner iterations lower this value is a good thing If it remains stable at around 1e-2 and you get stable coefficients for the drag and lift, I would say it is getting close to a good first run, but the bounding values on k and epsilon are not that ideal. Maybe try to switch back to kOmegaSST, and use omegaWallFunction with the blending option to binomial for an y+ independant wall function (and nutUBlendingWallFunction for nut).

[AMR96]

Originally Posted by Alczem

blending option to binomial for an y+ independant wall function (and nutUBlendingWallFunction for nut).

Hi Alcezm,

Thanks for your help and advice. I just have one question regarding your suggestion. What does setting the blending option to binomial do?

Also, should I use boundary layers with specific first layer heights in order to keep the y+ value below 1?

For the previous simulation, I was using the aspect ratio boundary layer (AR=5) and number of layers set to 5. However, the y+ values on the car were not correct since they were below 30. On the tunnel walls and the ground I have better values.

[Alczem]

The blending option allows you to choose different ways to manage the transition zone between y+>5 and y+>30. The binomial setting makes it smooth and continuous. By default, it is using stepwise which is disconstinuous. Try different settings, there are a few. There is also nutUSpaldingWallFunction which is able to handle every y+. Check out their description on the OF documentation, it is explained better than I do



Using these boundary conditions with omegaWallFunction and kqRWallFunction "should" make the wall treatment independant from your y+ so you don't have to worry as much about layers and first cell thickness (and that's why I was asking you to try with a coarse mesh with no layers). Of course, if your results are far too different from the expected values, you can start working a better and more accurate wall modelling, but make sure the differences do not come frome a wrong setup of the boundary conditions.

[AMR96]

Thanks for your explanation.

I tried the binomial option with the nutUBlending wall function but I got floating point exception after 30 iterations. If I use first order schemes for omega and k, the simulation works however the pressure residual is still around 0.05 after 700 iterations.

I will try to see which configurations work best.

[AMR96]

Hi Alcezm,

Sorry for my numerous messages!

I have been trying to make k-omega SST work during the past week. I finally have something stable. What I did was that I decided to not use any wall functions since I have very small y+ values. I also set up faceLimited gradient limiters to make sure omega does not explode, as a result the bounding message for omega has disappeared.

After letting the simulation run for a while, I decided to play around with the relaxation factors to see if I can get better residuals. Since the pressure residual was changing really slowly, I decided to increase it up to 0.8 and set the velocity residual to 0.2. This helped with lowering the pressure residual to around 0.009. However the pressure residual is taking quite some time to work go down. I wanted to see if you have any advice on how I can accelerate the pressure residual convergence.

Here is the latest iteration of the simulation:

Time = 3969

smoothSolver: Solving for Ux, Initial residual = 0.0001159485, Final residual = 2.1835152e-06, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.001512918, Final residual = 2.9030925e-05, No Iterations 2
smoothSolver: Solving for Uz, Initial residual = 0.00035526131, Final residual = 7.0361856e-06, No Iterations 2
GAMG: Solving for p, Initial residual = 0.0095298769, Final residual = 3.7831935e-05, No Iterations 4
GAMG: Solving for p, Initial residual = 0.00069856619, Final residual = 4.0308205e-06, No Iterations 4
GAMG: Solving for p, Initial residual = 6.9457521e-05, Final residual = 8.809016e-07, No Iterations 4
time step continuity errors : sum local = 4.3554646e-07, global = -3.1564626e-08, cumulative = -0.00011683361
smoothSolver: Solving for omega, Initial residual = 2.2712009e-05, Final residual = 2.2212845e-06, No Iterations 1
smoothSolver: Solving for k, Initial residual = 0.0023921728, Final residual = 0.00021069256, No Iterations 2
bounding k, min: -0.39160315 max: 6.3462664 average: 0.058743534
ExecutionTime = 76720.02 s ClockTime = 77256 s

forceCoeffs force_coefficient write:
Coefficient Total Pressure Viscous Internal
Cd: 0.93310885 0.85012659 0.082982262 0
Cd(f): 0.49064226 0.44965085 0.040991402 0
Cd(r): 0.44246659 0.40047573 0.04199086 0
Cl: -0.75699953 -0.75746573 0.00046619877 0
Cl(f): -0.45293633 -0.45375812 0.00082179428 0
Cl(r): -0.3040632 -0.30370761 -0.00035559551 0
CmPitch: -0.074436563 -0.075025258 0.0005886949 0
CmRoll: 0.024087831 0.02458756 -0.0004997291 0
CmYaw: 0.087711516 0.083018467 0.0046930485 0
Cs: 0.59764956 0.5623019 0.035347657 0
Cs(f): 0.3865363 0.36416942 0.022366877 0
Cs(r): 0.21111326 0.19813248 0.01298078 0

Sorry for the long message!

Best,

Amirmohammad

[Alczem]

Hey!


Glad you managed to make it stable! I have very little experience using no wall functions so I can't tell you why it works although it should have worked even with the wall functions since they can tell which Y+ approximation they should use depending on the cell size.



The faceLimited limiter is the most aggressive one, make sure it suits your needs concerning accuracy. But I guess it is okay if it is used only for turbulence.


And again, concerning convergence, do not rely only on residuals. Check if your coefficients are stable. If they are oscillating, it might indicate an unsteady flow which could be why the residuals do not go down much farther.


Keep us posted on the results and the expected values you should be getting, this should be interesting

[AMR96]

Hi,

Thanks for your reply.

I tried the cellLimited but it didn't work so I had to use faceLimited.

I actually looked at the drag and lift coefficient curves and drag has a really good convergence however lift has some oscillations which I believe could be due to somewhat higher residual of velocity in the lift direction. I also looked at the velocity flow field and it seems like it behaves as expected.

I believe the flow is unsteady since my domain is quite small and the car is placed at a short distance from the inlet and outlet (compared to what is simulated for a car on the road). I think the latter causes some undesired effects as the flow doesn't have enough time to become stable once it enters the domain.

I will keep you posted on what I find.

For now I am trying to see if I can impose a boundary condition at the outlet where the pressure will not be a fixed value, this can be closer to reality of what is happening. I am using a fixedGradient boundary at the outlet. Any insight regarding this would be really helpful.

Best,

Amirmohammad

[AMR96]

Hi Alcezm,

I hope you are doing well!

I took your advice and started looking into the lift and drag coefficients. The drag coefficient seems to be more or less stable and only oscillates between small values of around 0.91 and 0.92. I studied 4 different meshes and averaging the drag coefficient over the last 500 iterations gave me a value of around 0.92 for all the cases except the coarsest mesh which was 0.919.

The main problem lies in the lift coefficient, there is considerable oscillations between 0.68 up to 0.75. I noticed that as I refined the mesh the oscillation got worse. I tired averaging the values from the coarsest mesh to the most refined and this is what I get:
0.70 0.725 0.728 0.688

I have read on the forum that if the nature of the simulation is unsteady, refining the mesh could result in more oscillations, could this be the case here or is there another problem?

I would appreciate your insight on this.

I have attached graphs of the lift coefficient to the post. (mesh2 is the coarsest and mesh5 is the most refined)

Best,

Amirmohammad

[Alczem]

Hey!


I don't do a lot of Cl/Cd studies but it does not seem strange. I would try to refine the mesh one more time to see if you can get an average lift coefficient that remains stable (so far it keeps decreasing). If you are not satisfied with the average value, try to run a transient simulation, it will be more accurate, and you will be able to see if your flow detaches from your geometry which could be the source of the oscillations.


How to intepret this kind of results is slightly out of my "expertise" so maybe try to find similar threads about lift and drag (if you have not already) and get your hands on some book about CFD and aerodynamics