[DaveR]

Hi all!

I'm currently running a simulation to model the water flow in a swimming pool. There are two inlets with velocity at ~2.5 m/s, a "skim" outlet (top of the domain) and a "sump" outlet (bottom). Figure 1 illustrates the domain. The aim has been to use this model to learn a range of OpenFOAM (v.4) techniques such as steady-state, transient, turbulent and eventually free surface modelling.

Figure 1 - Domain

I have calculated the flow to be turbulent and I have attempted to simulate with both a k-e and k-Omega-SST model. The model runs fine in simpleFoam, albeit with a residuals output higher than i'd like (Figure 2).

residuals_simpleFoam&k-w-SST.png Figure 2 - residuals for SimpleFOAM (steady-state) & k-w-SST model (Click to enlarge)

From what I saw here, due to the magnitude and relative unsteadiness of the residuals I'm assuming transient modelling is more appropriate for this simulation. I have tried to simulate turbulence in transient using pisoFoam, however this quickly spits out floating point errors. I calculated my turbulence parameters from here, here and here. I initially used a k-Omega-SST model based on this article , however since I was having issues I changed to a k-e model inline with the cavity tutorial in the OpenFOAM user guide.

Considering the list of articles here, from watching the logs I get the impression my epsilon value is destabilising (rapidly diverges to a high negative value) and hence is the reason for the crashes. This seems to be backed up by the fact that my simulation runs fine in transient if I turn the turbulence off in the transportProperties dictionary. Doing this gives me residual values that, to my knowledge, are much more acceptable (Figure 3). As such, I'm assuming my problem lies within the definition of turbulence somewhere and have switched back to a k-w-SST model as this article suggests it to be better for nozzle flows, which I'm assuming is what my inlets represent.

residuals_pisoFoam&laminar.png Figure 3 - residuals for pisoFoam (transient) but with turbulence off (laminar flow) (Click to enlarge)

I've attached my case if anyone is interested in taking a look at my simulation setup. https://drive.google.com/open?id=0B3...WxTMU1aNDktR0E

I think the problem arises from my state of knowledge, however having scoured the internet for a few weeks now I can't seem to find a explanation that sticks, particularly;

• In practical engineering applications, can anyone suggest an appropriate thought flow on determining whether a simulation is transient or steady state?

• My turbulence calculations rely heavily on guesses such as the turbulent intensity and length scale. While I assume I'll develop confidence in my guesses as I become more experienced, can anyone suggest a good list of guidelines while I'm still developing said experience? Although, from the information I've looked at, perhaps these initial guesses are not particularly critical to the simulation?

• My turbulence calculations seem dissimilar in their order of magnitude compared with other tutorial cases (k=0.00384, e=0.154) Is there a glaring issue here?

• The fvSchemes and fvSolutions dictionaries are still a little alien to me. I've looked through the OpenFOAM user guide and my interpretation is for real engineering applications these do not seem to deviate far from the tutorial cases. Is my understanding here correct, or do my issues with floating point errors perhaps lie in these files?

Can anyone shed any light on these points/my simulation & understanding in general? Your help would be greatly appreciated!


Apologies if my post is a little lengthy - trying to provide all the details!

Best regards,

Dave

[tomf]

Hi Dave,

I have had a quick look at your case and the most important thing I can see is the erroneous definition of your patches (see attached pictures). Overall I think the quality of the mesh is very good, but the mesh may be too much refined.

For your settings I have seen some errors:

First of all, your timestep still gets you slightly above a Courant number of 1 (from the log file), so I guess it would be best to lower it a bit more.

Second of all, I would suggest the following modifications to the fvSchemes file:

gradSchemes
{
    default          Gauss linear;
    grad(k)          cellLimited Gauss linear 1.0;
    grad(omega)  cellLimited Gauss linear 1.0;
    grad(epsilon) cellLimited Gauss linear 1.0;
}

If this does not solve the issue, you may also want to change the divSchemes for the turbulence as:

divSchemes
{
    default         none;
    div(phi,U)      Gauss limitedLinearV 1;
    div(phi,k)      Gauss upwind; 
    div(phi,epsilon) Gauss upwind;
    div(phi,omega)  Gauss upwind;
    div(phi,R)      Gauss limitedLinear 1;
    div(R)          Gauss linear;
    div(phi,nuTilda) Gauss limitedLinear 1;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

For the fvSolution file, I would suggest to change the number of nonOrthogonalcorrectors to 1, the log file shows most of them do not do anything (No Iterarions 0) anyway and maybe also the nCorrectors to 2. This is not making the solution more stable, but just cleaning up your settings a bit.

    nCorrectors     2;
    nNonOrthogonalCorrectors 1;

For the solvers section in fvSolution, maybe it helps to add a minimal number of iterations for the momentum equations:

    "(U|k|epsilon|omega|R|nuTilda)"
    {
        solver          smoothSolver;
        smoother    GaussSeidel;
        tolerance    1e-05;
        relTol          0;
        minIter       1;
    }

I am not sure how to fix your mesh, but it is probably due to snappyHexMesh, which I do not use regularly.

Regards,
Tom

[DaveR]

Hi Tom,

Thanks for your comprehensive reply! Apologies in the lateness of mine I haven't been very well.

Yesterday evening I made the adjustments to my fvSchemes and Solution as suggested and the simulation ran! Perhaps you could explain a little about how you knew to do this? I still find fvSchemes and Solutions are my poorest areas of knowledge with regards to CFD/OpenFOAM and i'm struggling to adjust them too far away from the tutorial cases for my specific runs.

With regards to the mesh, my approach was to model the geometry in solidworks with inlets & outlet openings, then "cap" them with individual thin plates to represent my patches. I could then apply boundary conditions to these patches. I think that is why the images are showing a strange overlap?

I'm not convinced this is the right way to approach modelling custom geometries as the reason the mesh count is so high is that the features (inlet/outlet patches) are lost without very high resolution in these areas. Could you suggest a better way of approaching this?

I'm also a little confused as to how to approach calculating the minimum cell length for use with the Courant number. Currently i've been assuming an aspect ratio of ~1 and taking the cube root of the minimum cell volume as shown in the checkMesh results. i'm not sure this is correct though, as my time step is calculated for Co=0.5 yet my simulation runs at Co = 1.2 like you say. Is there a more efficient approach to this?

Currently i'm getting a time step so miniscule that on my machine it would take approximately 18 years to complete! On the transient simulation my Courant number settles to ~0.2 so I guess increasing the time step again could aid this.

Many thanks for your help!

Best regards,
Dave

[tomf]

Hi Dave,

No problem. Part of the changes come from experience, but some of it is also from knowledge about the finite volume discretisation method and some particular issues with the numerical solution of the Navier-Stokes equations. A lot of information can be found in "Computational Methods for Fluid Dynamics" by Peric and Ferziger.

For the meshing I would assume it may be better to separate the inlets/outlets from the main geometry as diferent STLs and have openings in the main geometry, however I do not use snappy that often, so I am not sure.

The high resolution for the patches is indeed necessary, but it is not in the bulk of your domain. So you could use a coarse base mesh and refine up to higher levels near these patches, and perhaps use some refinement boxes near these jets.

For the Courant number your approach is reasonable, but there may be some trial and error involved to get the optimal timestep.

Good luck,
Tom

[DaveR]

Hi Tom,

Thanks for the reference - I've downloaded the book and i'll be sure to work my way through it.

Apologies for the confusion - but the method of meshing you describe is what I did in this case. I'm not sure if it's just a visualisation error rather than anything erroneous about the mesh? I'm going to look at reducing the mesh density and hopefully maintain an acceptable order of magnitude for the residuals, however this seems difficult to achieve at present

Thanks for clarifying that about the Courant number, I've never been sure if this method was a reasonable approach.

I've modified a couple of other transient simulations using your modified fvSchemes and Solutions and they also run - I can't thank you enough for your help!

Best regards,

Dave

[tomf]

Hi Dave,

The errors also pop up in the result of checkMesh, where there are a couple of patches reported with a very low number of faces:

Checking patch topology for multiply connected surfaces...
                   Patch    Faces   Points                  Surface topology
                    body  1808760  1864011  ok (non-closed singly connected)
                  inlet1       54      213  ok (non-closed singly connected)
                  inlet2       35      138  ok (non-closed singly connected)
                    skim        4       16  ok (non-closed singly connected)
                    sump       53      219  ok (non-closed singly connected)

So I think there is somehow something incorrect with the meshing.

Good luck,
Tom