[talop]

Hi all,

Long time reader, first time poster

I have been trying to simulate convective heat losses from an isothermal, horizontal cylinder under cross flow in the mixed and natural convection regimes at steady state. This steady state solution will be used to initialize a transient case. Upon successfully simulating this reference case, I plan to move on to rotating bodies and detached-eddy simulations. However, I was unable to get a solution, let alone comparing it with the literature.

The problem geometry consits of a horizontal cylinder with a diameter of 1 meter placed in a 10x10 m square cross section duct. The cylinder is placed 6.5 m downstream of the inlet and the outlet is 12.5 m downstream from the cylinder.

The mesh is generated in ANSYS Fluent using the poly-hexcore method with quality above 0.5. Running checkMesh does not give any errors. Maximum skewness is 3.65 and non-orthogonality is maximum 59.85 with an average of 8.10. For comparison two more meshes were generated in ANSYS Meshing, one using a mapped mesh with hexahedral elements, other using tetrahedrons outside the boundary layer.

So far, I have been able to run some cases with significant forced convection, using the upwind schemes only, with k-epsilon and k-omega SST turbulence models. The air was flowing in at a rate 1 m/s. The inlet temperature was 288.15 K while hot cylinder was at 600 K (Re ~ 24000, Gr ~ 3.97e+09 and Ri ~ 6.9).

1) If a case dominated by natural convection is to be simulated (obtained by reducing the inlet velocity), solution crashes during p_rgh iterations after a few time steps.

2) Similarly, if the divergence scheme for the velociy term "div(phi,U)" is not "bounded Gauss upwind", the simulation crashes regardless of the inlet velocity at the same stage but at a different "time step" (anywhere between 5 to 500).

3) The same mesh performs well in ANSYS Fluent while using both segregated (SIMPLE) and coupled solver with second order upwind schemes except for the k-equation with the same relaxation factors or higher default values. SIMPLE solver takes ~1000 iterations to reach a solution.

I have tried to reduce the relaxation factors, use more correctors, use "type coupled" keyword for solvers and initialize the internal field with the inlet values. Using the coupled solver and reducing the relaxation factors helped slightly but the stability problem is still there.

Using the blockMesh is not possible for me because I will move towards more complicated geometries, which will be difficult , if not impossible, to mesh with blockMesh. I have tried using the snappyHexMesh in a different setting but I was unable to obtain good boundary layer meshes around internal corners. Therefore I will stick with the ANSYS Mechanical or Fluent as a mesher for now.

My guess is that I am using sub-optimal settings at best, leading to divergent behavior for higher order schemes or highly transient caeses. I will be glad if more experienced users can suggest some better settings for running mixed convection simulations using two-equation models and/or DES.

You can find the case files (except mesh because of size limitations) in the attached .zip file.

Best

[talop]

Update: Using a hex mesh gave slightly better results but simulations for the natural and mixed convection cases still blow up even with upwind schemes.


Best,
talop

[Alczem]

Hey!


Natural convection has been a long time nemesis of mine too!


From what I have read and tried, transient or pseudotransient are better at this kind of simulation. If you are using the ESI version, it is also possible to use hydrostaticPressure function object to initialize p and p_rgh with a hydrostatic distribution and a ph_rgh field. I found it can help.


Mesh wise, hexaedral is indeed better from my experience.


How low did you get with the relaxation factors?

[talop]

Hi, Alczem

I have tried initializing the flow field with some velocity to designate the flow direction and it did work partially. At least the simulations do not crash now. But it did not solve the negative p_rgh and back flow problem completely.

I have seen that the fluid "collapses" under gravity and this leads to a back-flow. To solve this, I tried initializing the pressure field (0/p) with overpressure but this did not solve the back flow either. I will try the hydrostaticPressure to see if it helps.

The fvSchemes uses "stadyState" for the ddtSchemes and therefore it runs in the SIMPLE mode even if it reads PIMPLE in the fvSolution file below.

solvers
{
    p_rgh
    {
        // type             coupled; // segregated;
        solver           GAMG; // PBiCICG; // GAMG;
        tolerance        1e-7;
        relTol           0.01;

        smoother         DICGaussSeidel;
        preconditioner   DILU;

    }

    "(U|h|k|epsilon|omega)"
    {
	// type         coupled; // segregated;
        solver          PBiCGStab; // PBiCICG; // PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-8;
        relTol          0.1;
    }
}

PIMPLE
{
    momentumPredictor no;

    nCorrectors               5; // 3; 
    nNonOrthogonalCorrectors  5; // 3;
    nOuterCorrectors          2; // 1; 

    pRefCell        0;
    pRefValue       0;

    residualControl
    {
        p_rgh           1e-4;
        U               1e-4;
        h               1e-4;

        // possibly check turbulence fields
        "(k|epsilon|omega)" 1e-3;
    }
}

relaxationFactors
{
    fields
    {

        rho             0.85; // 1.0;
        p_rgh           0.5; // 0.7;
    }
    equations
    {
        U               0.5; // 0.3;
        h               0.3;
        "(k|epsilon|omega)" 0.4; // 0.7;
    }
}


// ************************************************************************* //

[Alczem]

Hey,


Regarding your fvSchemes, I would use limited 0.33 instead of orthogonal since your mesh is not perfectly orthogonal. You could also try to add a limiter to the gradient calculation.


Why does your p_rgh values are set to 0 instead of something like 101325? I noticed you set pRefValue to 0 in fvSolution and there is a pRef file with 101325 in constant. I am not familiar with OF10, but are you sure the pressure is properly initialized? Again, I don't know exactly how it should be configured.


Cheers!

[talop]

Hi,

I have solved some of the issues and implemented better schemes. Getting there slowly

The p_rgh = p - rho*g*(h - h_ref) is the pressure without hydrostatic effects w.r.t. a reference point and gravity vector. As far as I know, it is defined as gauge pressure referencing the pRef file. In the tutorial cases this definition was used in the tutorial cases.

After tinkering with the case more, I have identified two problems:

1) constant/physicalProperties: equationOfState perfectGas;
The perfectGas uses variable pressure in rho = P/RT. FOr this case, pressure variation is small so incompressiblePerfectGas can be used with a reference pressure of 101325 Pa.

2) system/fvSolution: nOuterCorrectors 3;
By increasing the number of outer correctors to 10-15, I have been able to use higher order schemes for the divergence terms (limitedLinear, vanLeer ... etc.)