[hugol]

Dear community,

I started learning OpenFOAM a few months ago to solve the internal flow of an aerosol virtual impactor that I am designing (passive particle size discriminator). This technique consists in accelerating the air flow through an acceleration nozzle at a specified velocity, and splitting it into a minor flow (1 L/min here) facing the acceleration nozzle and a major flow (7.5 L/min) at 90 degrees so that small particles tend to follow the streamlines and go into the major flow and big particles into the minor flow.

The current set up aims to have a cutoff diameter of 1 um, so the max mach number should be about 0.7. In consequence, rhoSimpleFoam should be the most appropriate solver for steady-state, using kOmegaSST for RANS model (no particles to be modeled for the moment).

Throughout my many trials to establish mesh convergence, I was not able to reach convergence with 2nd order schemes as the solver would crash, especially for finer meshes. I noticed that the pressure would be constantly bounded.
With 1st order schemes (bounded Gauss upwind), convergence could be reached but velocity and pressure monitoring showed huge discrepencies between the meshes so the problem was not solved.

Following recommandations of other threads stating that rhoSimpleFoam could be highly unstable sometimes, I tried rhoPimpleFoam but the results were roughly similar. The solver would crash when schemes were switched to 2nd order (after 1st order "convergence"), with a higher tolerance for coarser meshes. This time, omega and k were bounded, not p.

I have tried to play with the URF, ramp flow or to fix temporary outlets pressure but no success.

I am using OpenFOAM v2412, the simulation is 3D and the mesh is generated with cfMesh. checkMesh returns OK for all meshes even if the finest mesh (1.6 million cell) has a few non ortho over 70 (max 77).

0/U:

/*--------------------------------*- C++ -*----------------------------------* \
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2412                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];
internalField   uniform (0 0 0);          

boundaryField
{
    // ---------- INLET ----------
    inlet
    {
        type        pressureInletOutletVelocity;
        value       uniform (0 0 0);      
    }

    // ---------- OUTLETS ----------
    outlet_min
    {
        type            flowRateOutletVelocity;
        massFlowRate    table  
        (
            (0 0.53e-5)
        );   // [kg/s] final 1/4*2.12e-5
        value           uniform (0 0 0);    // dummy, overridden
    }

    outlet_maj
    {
        type            flowRateOutletVelocity;
        massFlowRate    table  
        (
            (0 0.38e-4)
        );    // [kg/s] final 1/4*1.70e-4
        value           uniform (0 0 0);    // dummy, overridden
    }

    // ---------- WALLS ----------
    walls
    {
        type            noSlip;
    }
    front
    {
        type            noSlip;
    }
    // ---------- SYMMETRY ----------
    sym
    {
        type            symmetry;
    }
    back
    {
        type            symmetry;
    }
}


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

0/p:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2412                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [1 -1 -2 0 0 0 0];
internalField   uniform 95000;       // initial guess (Pa)

boundaryField
{
    // ---------- INLET ----------
    inlet
    {
        type        totalPressure; //totalPressure
        p0          uniform 95000;       // ambient static pressure
        value       uniform 95000;
    }

    // ---------- OUTLETS ----------
    outlet_min
    {
        type    fixedFluxPressure;
        rho     rho;
        value   uniform 0;
    }
    outlet_maj
    {
        type    fixedFluxPressure;
        rho     rho;
        value   uniform 0;
    }

    // ---------- WALLS ----------
    walls
    {
        type            zeroGradient;
    }
    front
    {
        type            zeroGradient;
    }
    // ---------- SYMMETRY ----------
    sym
    {
        type            symmetry;
    }
    back
    {
        type            symmetry;
    }
}


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

0/T:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2412                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 1 0 0 0];
internalField   uniform 300;          // K

boundaryField
{
    inlet      
	{ 
		type 	fixedValue;
		value 	uniform 300;
	}
    outlet_min 

	{ 
		type 	zeroGradient; 
	}
    outlet_maj 
	{ 
		type 	zeroGradient; 
	}
    walls
    {
        type      zeroGradient; 
	}   
    sym        
	{ 
		type 	symmetry; 
	}
    front      
	{ 
		type 	zeroGradient; 
	}
    back       
	{ 
		type 	symmetry; 
	}
}


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

I would greatly appreciate your help.
Thank you for your time.

[Tobermory]

When you say that you have played with URFs, what have you tried?

Also, does it run better if you switch one of the outlets off (turn to a wall)? It would be interesting to know whether its a numerics issue (most likely, esp if your mesh is fine - you essentially have less numerical dampening) or a boundary condition issue.

Finally, what is your mesh quality (orthoginality, skewness) like? These can sometimes impart instability unless controlled in the fvSchemes.

[hugol]

Dear Tobermory,

Thank you for your inetrest to my post.

Regarding the relaxation factors, I've tried to decrease them a lot (up to 0.1 for most of them, 0.05 for p, and even 0.01 for rho), sometimes alltogether, sometimes only a few at once, but in the end it is always crashing.

I am currently trying to run with a closed outlet, the calculation has not converged yet in first order but it shows roughly the same behaviour than before, i.e. bounding p and k/omega with rhoSimpleFoam and only bounding k/omega with rhoPimpleFoam.
Currently I am revising my boundary layer specification as omega blows up near the walls, but I am not very confident.

Please find below the checkMesh of my intermediate mesh that converges in 1st order but not in 2nd order (linearUpwind limited):

Mesh stats
    points:           530889
    faces:            1507343
    internal faces:   1423133
    cells:            488156
    faces per cell:   6.00315
    boundary patches: 7
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     472476
    prisms:        7216
    wedges:        0
    pyramids:      3336
    tet wedges:    0
    tetrahedra:    1376
    polyhedra:     3752
    Breakdown of polyhedra by number of faces:
        faces   number of cells
            9   2880
           12   548
           15   324

Checking topology...
    Boundary definition OK.
    Cell to face addressing OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
                   Patch    Faces   Points  Surface topology
                    back    27000    27468  ok (non-closed singly connected)
                   front    27000    27468  ok (non-closed singly connected)
                   inlet      855      920  ok (non-closed singly connected)
              outlet_maj      874      940  ok (non-closed singly connected)
              outlet_min      741      800  ok (non-closed singly connected)
                     sym     4237     4480  ok (non-closed singly connected)
                   walls    23503    24800  ok (non-closed singly connected)
                    ".*"    84210    83556      ok (closed singly connected)


Checking faceZone topology for multiply connected surfaces...
    No faceZones found.

Checking basic cellZone addressing...
    No cellZones found.

Checking basic pointZone addressing...
    No pointZones found.

Checking geometry...
    Overall domain bounding box (0.00468622 0.04 0.0075) (0.016 0.0615149 0.008)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (1.4486e-15 -9.01761e-16 1.49343e-14) OK.
    Max cell openness = 1.24363e-15 OK.
    Max aspect ratio = 79.3439 OK.
    Minimum face area = 1.71327e-12. Maximum face area = 1.06241e-08.  Face area magnitudes OK.
    Min volume = 6.35694e-17. Max volume = 1.19149e-12.  Total volume = 1.58973e-08.  Cell volumes OK.
    Mesh non-orthogonality Max: 58.8981 average: 5.20859
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 2.73945 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

End

My finer mesh (1.6 M cells) also converges fine in first order but blows up instantly when switching schemes, non-orthogonality is a bit higher with ~80 faces over 70 (but average at 6). Other parameters are similar to the intermediate.