[Martin Ramirez]

Hi all,

I thank you for entering this publication.

The context of my concern is as follows, I want to simulate the discharge of a pressure pipe on a hydraulic structure and observe the interaction that occurs with a pipe at the outlet for open flow conditions, initially create a mesh of the structure and pipe in conjunction with an elbow as accessory for the discharge at the entrance, in the SALOME 9.3 software, creating the mesh of the faces of 'wall-elbow', 'Atmosphere', 'wall-solid', 'Inlet' and 'Outlet' In geometry, two objects were merged, which were an inspection box with their respective pipes and an elbow. The mesh was imported into OpenFOAM through the ideasUnvToFoam intermediary which allows to bring the unv file to a representation of 'polyMesh'. For which the 'boundary' file is built, which contains 6 faces, taking the output and input of the elbow as 'defaultFaces', which are modified with an 'empty' value so that they do not have an interaction within the simulation.

FoamFile
{
        version            2.0;
        format             ascii;
        class                polyBoundaryMesh;
        location            "constant/polyMesh";
        object               boundary;
}

6
(
        wall-elbow
        {                      
              type           wall;
              nFaces        290;
              startFace     581153;
         }
         Atmosphere
        {                      
              type           patch;
              nFaces        290;
              startFace     581153;
         }
          wall-solid
        {                      
              type           wall;
              nFaces        290;
              startFace     581153;
         }
          Inlet
        {                      
              type           patch;
              nFaces        290;
              startFace     581153;
         }
          Outlet
        {                      
              type           patch;
              nFaces        290;
              startFace     581153;
         }
           defaultFaces
        {                      
              type           empty;
              nFaces        290;
              startFace     581153;
         }
)

Once the analysis mesh was imported, we proceeded to configure 'system', for which 'setFieldsDict' was edited, which expressed the regions with an alpha value of 1 for a section of cylindrical geometry, taking the radius and position of the entry.

FoamFile
{
        version            2.0;
        format             ascii;
        class                dictionary;
        location            "system";
        object              setFieldsDict;
}

defaultFieldValues
(
       volScalarFieldValue alpha.water 0
);

regions
(
        cylinderToCell
        {
             p1 (0.40 -1 0.40);         // center of first circle
             p2 (0.40 0 0.40);          // center of second circle
             radius 0.04021;            // radius of cylinder
             fieldValues
             (
                  volScalarFieldValue alpha.water 1
             );
         }
);

The solver used to analyze the boundary conditions, properties of the fluid and the physics of the flow, is carried out by means of 'interFoam', taking the inlet pipe with a full tube analysis (that is, it contains only water), the discharge of the fluid on the bottom of the structure and its drainage to the outlet pipe it is taken as a relation of multiphase fluids where air and water interact. The physics of the flow is observed with the appreciations of turbulent flow, for this reason the K-Omega model is used, it solves the RAS equation by means of 'PIMPLE'.

FoamFile
{
        version            2.0;
        format             ascii;
        class                dictionary;
        location            "system";
        object              setFieldsDict;
}

solvers
{
       "alpha.water.*"
       {
             nAlphaCorr             1;
             nAlphaSubCycles     1;
             cAlpha                    1;

             MULESCorr              yes;
             nLimiterIter             3;

             solver                      smoothSolver;
             smoother                 symGaussSeidel;
             tolerance                 1e-8;
             relTol                       0;
        }
        
         pcorrFinal
       {
             solver             PCG;
             preconditioner
             {
                    preconditioner         GAMG;
                    tolerance                  1e-5;
                    relTol                        0;
                    smoother                GaussSeidel;
                    nPreSweeps              0;
                    nPostSweeps             2;
                    nFinestSweeps           2;
                    cacheAgglomeration    true;
                    nCellsInCoarsestLevel       10;
                    agglomerator              faceAreaPair;
                    mergeLevels               1;
              }
               tolerance      1e-5;
               relTol            0;
               maxIter        50;
         }
         p_rgh
         {
                    solver                    GAMG;
                    tolerance                 5e-9;
                    relTol                       0.01;

                    smoother                GaussSeidel;
                    nPreSweeps            0;
                    nPostSweeps           2;

                    cacheAgglomeration      true;
                    
                    nCellsInCoarsestLevel      10;
                    agglomerator           faceAreaPair;
                    mergeLevels            1;

                    maxIter                   50;
            };

            p_rghFinal
            {
                    $p_rgh;
                    tolerance                  5e-9;
                    relTol                        0;
             }

             "(U|k|omega).*"
             {
                    solver                       smoothSolver;
                    smoother                  symGaussSeidel;
                    nSeewps                   1;
                    tolerance                   1e-6;
                    relTol                         0.1;
             };
}

PIMPLE
{
             momentumPredictor           no;
             nCorrectors                          2;
             nNonOrthogonalCorrectors    0;
}

relaxationFactors
{
             equations
             {
                    ".*"  1;
              }
}

The determination of the VOF started from the methodology that establishes alpha.water 1 for the value of water and 0 for the value of air, for this purpose the entrance was taken as a uniform value of 1, the walls were established with the type 'zeroGradient'.

FoamFile
{
        version            2.0;
        format             ascii;
        class                volScalarField;
        location            "0";
        object               alpha.water;
}

dimensions                [0 0 0 0 0 0 0]

internalField               uniform 0;
boundaryField
{
        Inlet
        {
              type             fixedValue;
              value            uniform 1;
        }

        wall-solid
        {
              type              zeroGradient;
         }
        wall-elbow
        {
              type               zeroGradient;
        }

        Outlet
        {
              type               zeroGradient;
              value              uniform 0;
        }
        Atmosphere
        {
               type               inletOutlet;
               inletValue       uniform  0;
               value             uniform   0;
         }

}

The velocity value corresponds to an inlet flow to pressure pipe of 0.20 m3 / s.

FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 2.29 0);

boundaryField
{
    Inlet
    {
	type            flowRateInletVelocity;
        volumetricFlowRate constant 0.20;
    }

    wall-solid
    {
	type            fixedValue;
        value           uniform (0 0 0);
    }

    Atmosphere
    {
        type            pressureInletOutletVelocity;
        value           uniform (0 0 0);
    }

    Outlet
    {
	type            inletOutlet;
        inletValue      uniform (0 0 0);
        value           $internalField;
    }
    
    wall-elbow
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
}

For the pressure two moments are established, the first is the discharge pressure of the inlet pipe and the second is the atmospheric pressure.

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p_rgh;
}

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

internalField   uniform 78451.04;

boundaryField
{
    Atmosphere
    {
        type            totalPressure;
        p0              uniform 101325;
        U               U;
        phi             phi;
        rho             rho;
        psi             none;
        gamma           1;
        value           uniform 101325;
    }
        
    Outlet
    {
         type            zeroGradient;
    }

    Inlet
    {
        type            fixedFluxPressure;
        value           uniform 78451.04;
    }

    ".*"
    {
        type            fixedFluxPressure;
        value           uniform 101325;
    }

}

The turbulence energy value is taking into account the turbulence intensity for the hydraulic diameter.

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      k;
}

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

internalField   uniform 0.00969;

boundaryField
{
    Inlet
    {
        type            fixedValue;
        intensity       0.0351;
        value           $internalField;
    }

    wall-solid
    {
        type            kqRWallFunction;
        value           $internalField;
    }
    
    wall-elbow
    {
        type            kqRWallFunction;
        value           $internalField;
    }


    ".*"
    {
        type            inletOutlet;
        inletValue      $internalField;
        value           $internalField ;
    }
}

The dissipation of energy is taken using the omega model, for this the length of the turbulence is used.

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      omega;
}

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

internalField   uniform 32.22;

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

    wall-solid
    {
        type            omegaWallFunction;
        value           $internalField;
    }

    wall-elbow
    {
        type            omegaWallFunction;
        value           $internalField;
    }


    ".*"
    {
        type            inletOutlet;
        inletValue      $internalField;
        value           $internalField;
    }
}

The nut value is set for the input and output conditions.

{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      nut;
}

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

internalField   uniform 0;

boundaryField
{
    wall-solid
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
    
    wall-elbow
    {	
	type		nutkWallFunction;
	value		uniform 0;
    }

    ".*"
    {
        type            calculated;
        value           uniform 0;
    }
}

The turbulence properties are analyzed by means of the RAS equations solved with the k-omega model.

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      turbulenceProperties;
}

simulationType  RAS;

RAS
{
    RASModel        kOmegaSST;

    turbulence      on;

    printCoeffs     on;
}

The properties of the fluids follow the density, viscosity, and surface tension conditions for air and water.

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}

phases (water air);

water
{
    transportModel  Newtonian;
    nu              [0 2 -1 0 0 0 0]  1e-06;
    rho             [1 -3 0 0 0 0 0] 1000;
}

air
{
    transportModel  Newtonian;
    nu              [0 2 -1 0 0 0 0]  1.48e-05;
    rho             [1 -3 0 0 0 0 0] 1;
}

sigma           [1 0 -2 0 0 0 0] 0.07;

The gravity value is specified as 9.81 m / s2

FoamFile
{
    version     2.0;
    format      ascii;
    class       uniformDimensionedVectorField;
    location    "constant";
    object      g;
}

dimensions      [0 1 -2 0 0 0 0];
value           (0 0 -9.81);

To establish the Courant number, a deltaT value of 0.0000001 was used.

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}

application     interFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         10;

deltaT          0.0000001;

writeControl    adjustableRunTime;

writeInterval   0.05;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression compressed;

timeFormat      general;

timePrecision   6;

runTimeModifiable yes;

adjustTimeStep  yes;

maxCo           6;
maxAlphaCo      6;
maxDeltaT       1;

functions
{
    inletFlux
    {
        type            faceSource;
        functionObjectLibs ("libfieldFunctionObjects.so");
        outputControl   timeStep;
        log             true;
        // Output field values as well
        valueOutput     false;
        source          patch;
        sourceName      Inlet;
        operation       sum;

        fields
        (
            rhoPhi
        );
    }

    OutletFlux
    {
        $inletFlux;
        sourceName      Outlet;
    }
    
    AtmosphereFlux
    {
        $inletFlux;
        sourceName      Atmosphere;
    }
}

Once the 'interFoam' solver is executed, an error is generated in the model's journey based on the loaded conditions and characteristics, the stability of the system is damaged.

Courant Number mean: 0.00202408 max: 0.0216337
Interface Courant Number mean: 0 max: 0
deltaT = 7.19942e-06
Time = 7.19942e-06

smoothSolver:  Solving for alpha.water, Initial residual = 1.79803e-06, Final residual = 4.27735e-12, No Iterations 1
Phase-1 volume fraction = 0.00838449  Min(alpha.water) = 0  Max(alpha.water) = 1
MULES: Correcting alpha.water
Phase-1 volume fraction = 0.00838449  Min(alpha.water) = -1.61933e-13  Max(alpha.water) = 1
GAMG:  Solving for p_rgh, Initial residual = 1, Final residual = 1.24433e+37, No Iterations 50
time step continuity errors : sum local = 1.71056e+34, global = -5.28206e+32, cumulative = -5.28206e+32
GAMG:  Solving for p_rgh, Initial residual = 1, Final residual = 4.3059e-09, No Iterations 19
time step continuity errors : sum local = 1.42133e+59, global = 7.45096e+58, cumulative = 7.45096e+58
smoothSolver:  Solving for omega, Initial residual = 0.811838, Final residual = 0.0685701, No Iterations 7
smoothSolver:  Solving for k, Initial residual = 0.999997, Final residual = 0.0973059, No Iterations 40
bounding k, min: 8.78028e-62 max: 0.0785186 average: 0.00967384
ExecutionTime = 12.45 s  ClockTime = 14 s

Courant Number mean: 1.10307e+67 max: 2.31367e+69
Interface Courant Number mean: 0 max: 0
deltaT = 1.86701e-74
Time = 7.19942e-06

smoothSolver:  Solving for alpha.water, Initial residual = 1.2844e-08, Final residual = 2.55005e-13, No Iterations 1
Phase-1 volume fraction = 0.00838449  Min(alpha.water) = -7.91099e-14  Max(alpha.water) = 1
MULES: Correcting alpha.water
Phase-1 volume fraction = 0.00838449  Min(alpha.water) = -2.03375e-12  Max(alpha.water) = 1
GAMG:  Solving for p_rgh, Initial residual = 1, Final residual = 0.00547086, No Iterations 5
time step continuity errors : sum local = 2.88229e-05, global = -1.51418e-06, cumulative = 7.45096e+58
GAMG:  Solving for p_rgh, Initial residual = 0.310287, Final residual = 4.13946e-09, No Iterations 37
time step continuity errors : sum local = 4.48349e-11, global = 4.75626e-12, cumulative = 7.45096e+58
smoothSolver:  Solving for omega, Initial residual = 0.00265864, Final residual = 3.78146e-05, No Iterations 1
smoothSolver:  Solving for k, Initial residual = 2.73075e-06, Final residual = 9.66491e-11, No Iterations 1
ExecutionTime = 16.4 s  ClockTime = 18 s

Courant Number mean: 0.0109246 max: 3.60706
Interface Courant Number mean: 0 max: 0
deltaT = 2.17757e-74
Time = 7.19942e-06

smoothSolver:  Solving for alpha.water, Initial residual = 5.45264e-09, Final residual = 5.45264e-09, No Iterations 0
Phase-1 volume fraction = 0.00838449  Min(alpha.water) = -2.03375e-12  Max(alpha.water) = 1
MULES: Correcting alpha.water
Phase-1 volume fraction = 0.00838449  Min(alpha.water) = -2.03375e-12  Max(alpha.water) = 1
GAMG:  Solving for p_rgh, Initial residual = 0.399741, Final residual = 0.00178865, No Iterations 5
time step continuity errors : sum local = 2.75265e-05, global = -1.45885e-06, cumulative = 7.45096e+58
GAMG:  Solving for p_rgh, Initial residual = 0.296846, Final residual = 0.00128416, No Iterations 50
time step continuity errors : sum local = 1.17308e-05, global = -3.1497e-07, cumulative = 7.45096e+58
smoothSolver:  Solving for omega, Initial residual = 0.00213625, Final residual = 2.79442e-05, No Iterations 1
smoothSolver:  Solving for k, Initial residual = 0.000438583, Final residual = 9.42752e-08, No Iterations 1
bounding k, min: 2.22045e-16 max: 0.0785185 average: 0.00967385
ExecutionTime = 21.27 s  ClockTime = 23 s

Bearing in mind the context of the simulation and its errors, I would like to know your opinions regarding this situation and the possible bibliographic reviews that I must carry out to change, improve or configure the solution.

Thank you.

Martin Ramirez