[JEBland]

I'm attempting to get a square channel toy model working for pressure-driven laminar flow using simpleFoam, and the pressure is quite wild for the first few hundred iterations. The pressure is initialized to (approximately) the average pressure of the inlet & outlet. Is there an efficient way to limit the initial iterations to encourage stability in pressure-driven simulations?



I understand that there might be "shocks" in the flow when setting up pressure based BCs, but I would expect that for an initially still fluid with the interior pressure set to the average pressure of the boundaries, that the peak & min pressures should also be approximately bounded by inlet and outlet pressures - instead the peak min/max values of the pressure at orders of magnitude larger than the final resulting values, even with the relaxation factors set to 0.1 for p and U.



Explicit limitation of p doesn't seem to be respected by simpleFoam, either.



While the pressure-driven BCs are on the list of the -listScalarBCs and -listVectorBCs for potentialFoam, my attempts to initialize with potentialFoam were fruitless. Apparently this isn't new, but to my knowledge, pressure driven flow has not been incorportated into potentialFoam (https://bugs.openfoam.org/view.php?id=1845).



I don't mind that it can take a while to converge in simpleFoam, but it's frustrating that even with the relaxation set to 0.1 that there are still domains moving through and not being smoothed by the solver.



A couple of shallow attempts to bound the solver settings in fvSchemes didn't seem to do much, either.

https://www.openfoam.com/documentati...erical-schemes
https://cfd.direct/openfoam/user-guide/v6-fvSchemes/


Attached are a couple sample images, which represent some of the unstable initial behavior where the inlet has a pressure of 7.3 m^2/S^2 and the outlet has a pressure of 0.52 m^2/s^2. I'm using the totalPressure boundary condition for the pressure and the pressureInletOutlet - from the guide, this combination should have "very good" stability.



p:

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

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

internalField   uniform 3.9 ;

boundaryField
{
    Leg1InOut // inlet
    {
        type            totalPressure;
        p0              uniform 7.3;
        value           7.3; 
    } 
    Leg2InOut // outlet
    {
        type            totalPressure;
        p0              uniform 0.52;
        value           0.52; 
    }
    
    walls
    {
        type            fixedFluxPressure; 
        value           0;    
        gradient        0;
    }
    
}

U:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     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
{

    Leg1InOut
    { 
        type         pressureInletUniformVelocity; 
        value        uniform (0 0 0);
    } 
    Leg2InOut
    { 
        type         pressureInletOutletVelocity; 
        value        uniform (0 0 0);
    }
       
    walls
    {
        type            noSlip; // zeroGradient;
    }
}

blockMeshDict:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}

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


length  1.0;
width   0.1;
height  0.1;

xMin #calc "-0.5 * $length";
xMax #calc " 0.5 * $length";

yMin #calc "-0.5 * $height";
yMax #calc " 0.5 * $height";

zMin #calc "-0.5 * $width";
zMax #calc " 0.5 * $width";

vertices
( 
// 
//                 Leg1
//                  |
//                  |
//                 Leg2   
// 
    //
    ( $xMin $yMin $zMin) // 0  
    ( $xMax $yMin $zMin) // 1 
    ( $xMax $yMax $zMin) // 2
    ( $xMin $yMax $zMin) // 3
    
    ( $xMin $yMin $zMax) // 4  
    ( $xMax $yMin $zMax) // 5 
    ( $xMax $yMax $zMax) // 6
    ( $xMin $yMax $zMax) // 7  
    //   
);

 

blocks
(
//  
    hex (0 1 2 3 4 5 6 7) (100 10 10) simpleGrading (1 1 1) // Leg1 - Leg2
// 
);

 

edges
(
);

boundary
(
    Leg1InOut
    {
        type patch;
        faces
        (
        (0 4 7 3) // nX 
        );
    } 
    
    Leg2InOut
    {
        type patch;
        faces
        (
          (2 6 5 1) // pX 
        );
    } 
      
    walls
    {
        type wall;
        faces
        (
          (2 3 6 7) // pY
          (0 1 5 4) // nY bottom
          (1 0 3 2) // nZ 
          (4 5 6 7) // pZ 
        );
    } 
);


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

controlDict:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     simpleFoam;

startFrom       startTime;

startTime       0;  

stopAt          endTime;

endTime         10000; 

deltaT          1; 

writeControl    timeStep;

writeInterval   10;

purgeWrite      0;

writeFormat     binary; // ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


functions
{
    #includeFunc flowRatePatchpatches
        (
             Leg1InOut
             Leg2InOut
        ); 
}

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

fvSolution:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{    
     p
    {
        solver           GAMG;
        tolerance        1e-7;
        relTol           0.01;
        smoother         DICGaussSeidel;
    } 

    U
    {
        solver          smoothSolver; 
        smoother        GaussSeidel;
        tolerance       1e-05;
        relTol          0;
    }
    "(p|U)Final"
    { 
        relTol          0;
    } 
}
 

SIMPLE
{
    momentumPredictor no;  
    nNonOrthogonalCorrectors 0; 
    consistent       no;

    residualControl
    {
        p               1e-4; // 4
        U               1e-5; // 5
    }
}

relaxationFactors
{
    fields
    {
        p               0.1;  
    }
    equations
    { 
         U               0.1; 
    }
}

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

fvSchemes:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
}

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

laplacianSchemes
{
    default         Gauss linear corrected;
//    default         Gauss linear limited corrected 0.4; // Doesn't do much?
}

interpolationSchemes
{
    default         linear; 
}

snGradSchemes
{
    default         corrected;
} 

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

Any assistance would be appreciated.

[JEBland]

I checked against a sample with pore-based modelling, and the sample had the pressure on the order of 1e-3 (whereas my difference is 1e0 - 1e1). Doing a basic test, I decreased the pressure difference by introducing an order of magnitude so 7.3 becomes 7.3e-3, 3.9 becomes 3.9e-3, and 0.52 becomes 0.52e-3. Then repeated with e-2.. I readily get good convergence for e-3, e-2. But at *e-1, things get substantially slower.



Does the mesh size need to become finer in order to handle the larger pressure differences / faster velocities?

The relaxationFactors are already quite small at 0.05 (reduced from the 0.1 above).


I'm tinkering with just mapping the internal fields a la this thread:
MapFields only on internal fields
Initial test seems to work, but will see if I'm able to "turn on" the pressure bit by bit.


Edit: Looks like it needs a finer mesh, too. Will play with refining the mesh as well as mapping the previous fields in order to get more stability.