[unilord]

Hello everyone,

I am an Aeronautical Engineering student working in his Master Thesis about Numerical Analysis of a Centrifugal Pump.

The initial analysis is a steady state. I am using the MRF approach, with 3 domains separate, as you can see from the attached picture. The inlet is in the y direction the dark blue surface from the purple domain and the outlet is radial, represented by the light blue surface from the pink domain.

Thus, and since I am using OF v11, I need to make use of the NCC tool they provide. If you go to the documentationpage, you will see that, for surfaces that fully overlap, you should use movingWallSlipVelocity (or slip if the mesh is stationary) for velocity U; zeroGradient or fixedFluxPressure for pressure p; zeroGradient for other fields.

The mesh has around 720000 elements, the pump has around 30 cm diameter, with 600 rpm. The working fluid is water.

I am performing a laminar analysis, and using the following boundary conditions:

Velocity:

inlet
    {
        type            pressureInletOutletVelocity;
        value           $internalField;	     
    }
    outlet
    {
        type			flowRateInletVelocity;
        volumetricFlowRate	0.01356;
        profile			laminarBL;
    }
    
    "(wall_1|wall_2|bladeWalls|wall_3)" 
    {
        type            noSlip;
    }
    
    "(NCC1_2|NCC2_1)"
    {
        type            slip;
    }
    
    "(NCC2_3|NCC3_2)"
    {
    	type		noSlip;
    }

Pressure

inlet
    {
        type			fixedValue;
        value                   uniform 0;
    }
    outlet
    {
        type			zeroGradient;
    }
    
    "(wall_1|wall_2|bladeWalls|wall_3)" 
    {
        type            zeroGradient;
    }
    
    "(NCC1_2|NCC2_1)"
    {
        type            zeroGradient;
    }
    
    "(NCC2_3|NCC3_2)"
    {
    	type		zeroGradient;
    }

Where NCC1_2 / NCC2_1 are the patches corresponding to the adjacent surfaces of the blue and purple domains (fully overlap), and NCC2_3 / NCC3_2 are the patches corresponding to the adjacent surfaces of the purple and pink domain (partially overlap). The results are presented in the attached pictures.

I can't understand why there are "aureolas" where the non conformal coupling interfaces are. The tutorials are only showing moving components, but I want to simulate a static component. Does anyone have any idea of what I could be doing wrong?

Thank you in advance.

[prud1k]

Can you share your case? I'm student and really don't understand how NCC is working and some waves to move case from ACMI method to NCC

[unilord]

Hey,

I am now using another centrifugal pump that I found in a paper. What do you mean when you say you don't understand how NCC is working an some waves?

[prud1k]

Originally Posted by unilord

Hey,

I am now using another centrifugal pump that I found in a paper. What do you mean when you say you don't understand how NCC is working an some waves?

Dear Unilord,
I'm used OF v2306 to research aerodynamic characteristics of rotating projectile with feather stabilizers with using ACMI method. Now I need to do the same thing using the NCC method. There was a difficulty in translating the case from the AMCI method from v2306 to the NCC method from v11. Unfortunately, so far, as I understand it, the Internet is not yet rich in examples of using a rotating mesh using the NCC method. An additional difficulty lies in the syntax discrepancy between the branches of OF, which makes it difficult to switch from one version to another due to the imposition of errors due to the lack of ACMI support and the lack of perception of directories due to the syntax.
Sorry for the obvious questions, I'm a novice user OF

[unilord]

Originally Posted by prud1k

Dear Unilord,
I'm used OF v2306 to research aerodynamic characteristics of rotating projectile with feather stabilizers with using ACMI method. Now I need to do the same thing using the NCC method. There was a difficulty in translating the case from the AMCI method from v2306 to the NCC method from v11. Unfortunately, so far, as I understand it, the Internet is not yet rich in examples of using a rotating mesh using the NCC method. An additional difficulty lies in the syntax discrepancy between the branches of OF, which makes it difficult to switch from one version to another due to the imposition of errors due to the lack of ACMI support and the lack of perception of directories due to the syntax.
Sorry for the obvious questions, I'm a novice user OF

Dear prud1k,

I am also facing the same issue. There is really not enough data regarding the new feature "Non-conformal coupling" present in OF11 (and OF10 I believe). I am attaching here my case directory. However, I would not use it as a faithful example, since I am facing some issues myself.

For instance. I am using four different regions, thus having six non conformal patches (NCC1_2, NCC2_1, NCC2_3, NCC3_2, NCC4_3, NCC3_4). Four of those six regions are in contact with the MRF zone. One issue I am having is, should I declare them as nonRotatingPatches() in the /constant/MRFProperties file or not?

Anyways, if anyone could help us regarding this issue it would be awesome. I am leaving here the files you need and the methodology I am using. I will keep this thread updated with the results I am getting. I am testing with kEpsilon at the moment, although I was using SST for some weeks. Please comment if you have any doubts


Methodology:

1 - having a separate directory for each region you are using. Each directory should have a constant, system and 0 folder (I am not sure if the 0 folder is necessary).


2 - Meshing each component in the respective directory. I am currently using salome for meshing, since I believe tetraheadron are better for centrifugal pump than cartesian hexa.


3 - After meshing all components, you should use "mergeMesh" command to merge all meshes into one directory. Usage is "> mergeMesh /Domain_2 /Domain_1 - overwrite". Keep doing it untill all domain is in one folder.


4 - Perform a "checkMesh" command inside the folder with all the regions. It will state that there are "x" separate regions, and will write all of them in a "set" folder.


5 - If you are doing steady state, now you must select the MRFzone. I use "topoSet" for it. If you did step nr 4, you will have the MRF region already present in a /set folder. If you compare the nr of elements of each region (given in checkMesh) you will know which one is the impeller/fan.


6 - "createNonConformalCouples" command. Since I have more than one NCC, I have a file defining them, called "createNonConformalCouplesDict", so I just have to run createNonConformalCouples in the terminal.


7 - run.
----------------------------------


NonConformalCouplesDict

// Whether or not to add boundary conditions for the added patch to the fields
fields  no;

// The list of non-conformal couples to be created. Each entry in this section
// creates a single non-conformal coupling.
nonConformalCouples
{
    nonConformalCouple_12
    {
        patches         (NCC1_2 NCC2_1);
        transform       none;
    }
    
    nonConformalCouple_23
    {
        patches         (NCC2_3 NCC3_2);
        transform       none;
    }
    
    nonConformalCouple_34
    {
        patches         (NCC3_4 NCC4_3);
        transform       none;
    }
}

topoSetDict ( first part is commented because it's a geometric way of selecting the cells. If you want to use it instead of steps nr 4 and 5, you can also do so)

actions
(
/*
    {
        name        MRFSet;
        type        cellSet;
        action        new;
        source        cylinderAnnulusToCell;
        p1        (0 -0.01 0);
        p2        (0 0.02301 0);
        outerRadius    0.14301;
        innerRadius    0.05701;
    }
   */
    {
        name        MRFZone;
        type        cellZoneSet;
        action        new;
        source        setToCellZone;
        sourceInfo
        {
            set region2;
        }
    }
    
);

Pressure, p

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"
    
    "NCC.*"
    {
        type            zeroGradient;
    }
    outlet
    {
        type            pressure;
        p        $internalField;
        value        $internalField;
    }
    "wall.*"
    {
        type            zeroGradient;
    }
    "blade.*"
    {
        type            zeroGradient;
    }
    inlet
    {
        type            zeroGradient;
    }
}

Velocity, U

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"
    
    "NCC.*"
    {
        type            slip;
    }
    outlet
    {
        type            pressureInletOutletVelocity;
        value           $internalField;
    }
    "wall.*"
    {
        type            noSlip;
    }
    "blade.*"
    {
        type            MRFnoSlip;
        value           uniform (0 0 0);
    }
    inlet
    {
        type            flowRateInletVelocity;
        volumetricFlowRate 
        {
            type            constant;
            value           0.0063;
        }
        value           uniform (0 0 0);
    }
}

turbulent viscosity, nut

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"
    
    "NCC.*"
    {
        type            zeroGradient;
    }
    outlet
    {
        type            calculated;
        value           uniform 0;
    }
    "wall.*"
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
    "blade.*"
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
    
    inlet
    {
        type            calculated;
        value           uniform 0;
    }
}

Turbulent Kinetic Energy, k

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"
    
    "NCC.*"
    {
        type            zeroGradient;
    }
    outlet
    {
        type            inletOutlet;
        inletValue      $internalField;
        value           $internalField;
    }
    "wall.*"
    {
        type            kqRWallFunction;
        value           uniform 1e-12;
    }
    "blade.*"
    {
        type            kqRWallFunction;
        value           uniform 1e-12;
    }
    inlet
    {
        type            turbulentIntensityKineticEnergyInlet;
        intensity       0.05;
        value           uniform 1;
    }
}

Turbulent Kinetic Energy Dissipation Rate, epsilon

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"
    
    "NCC.*"
    {
        type            zeroGradient;
    }
    outlet
    {
        type            inletOutlet;
        inletValue      $internalField;
        value           $internalField;
    }
    "wall.*"
    {
        type            epsilonWallFunction;
        value           $internalField;
    }
    "blade.*"
    {
        type            epsilonWallFunction;
        value           $internalField;
    }
    inlet
    {
        type            fixedValue;
        value        $internalField;
    }
}

You need a file in /constant called "MRFProperties", with the following format:

MRFImpeller
{
    cellZone    MRFZone;
    active    yes;
    
    nonRotatingPatches (NCC1_2|NCC2_1|NCC2_3|NCC3_2);
    
    origin    (0 0 0);
    axis      (0 0 1);
    omega     constant 64.9262481;
}

[prud1k]

Dear unilord,
Thanks a lot. In addition to my question: Have you come across a GPU connection for calculations? On this forum, I read that the connection takes place using third-party libraries with PETSc and AMGXWrapper, do you know what other methods exist?

[unilord]

Originally Posted by prud1k

Dear unilord,
Thanks a lot. In addition to my question: Have you come across a GPU connection for calculations? On this forum, I read that the connection takes place using third-party libraries with PETSc and AMGXWrapper, do you know what other methods exist?

Hey prud1k,

As far as I know, there is no opensource/free tool that allows the usage of GPU in openfoam. There is one thing that I am aware about, which is called dicehub, but they do not provide GPU for openfoam (not 100% sure). The guys there gave me some hours in exchange of using my simulation to put in their website. However they only gave me in CPU hours. There is a freelancer engineer called Giles Richardson in LinkedIn that is always posting simulations made in UFOCFD (it's another code) ran in a GPU

[unilord]

I am also going to leave here the results of my analysis, in case that someone with more experience than me wants to provide feedback.


I am using kEpsilon at the moment. yPlus between 30 and 200. I am using upwind for divSchemes of k and epsilon (diffusive terms), and linearUpwind for convective terms. However, I want to use linearUpwind for diffusive terms as well. However, when I try to do so, my simulation gives me "bounding omega" value of around E+16, which is obviously wrong. The mesh is tetrahedral, with hexahedral as boundary layer and some polyhedral.


fvSolution

solvers
{
    p
    {
        solver                    GAMG;
        tolerance                 1e-08;
        relTol                    0.0;
        smoother                  GaussSeidel;
        agglomerator    faceAreaPair;
        nCellsInCoarsestLevel     1000;
    }

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-08;
        relTol          0.0;
    }

    k
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-08;
        relTol          0.0;
    }

    epsilon
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-08;
        relTol          0.0;
    }
    
    Phi
    {
        solver          GAMG;
        smoother        DIC;
        tolerance       1e-06;
        relTol          0.01;
    }
    
    yPsi
    {
        solver          GAMG;
        smoother        GaussSeidel;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 100;
        agglomerator    faceAreaPair;
        mergeLevels     1;
        tolerance       1e-6;
        relTol          0;
    }

}

potentialFlow
{
    nNonOrthogonalCorrectors 3;
}

SIMPLE
{
    momentumPredictor         yes;
    nNonOrthogonalCorrectors  2;
    nCorrector                0;
    
    residualControl
    {
        p        1e-6;
        U        1e-7;
        k        1e-7;
        epsilon        1e-7;
    }
}

relaxationFactors
{
    fields
    {
        p               0.3;
    }
    equations
    {
        U               0.7;
        k               0.7;
        epsilon         0.7;
    }
}

fvSchemes

ddtSchemes
{
    default        steadyState;
}

gradSchemes
{
    default        Gauss linear;
    grad(U)             cellMDLimited Gauss linear 0.333;
}

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

laplacianSchemes
{
    default        Gauss linear limited corrected 0.5;
}

interpolationSchemes
{
    default        linear;
}

snGradSchemes
{
    default        limited corrected 0.5;
}

wallDist
{
    method Poisson;
}

checkMesh

Checking geometry...
    Overall domain bounding box (-0.160034 -0.194078 -0.0123) (0.235363 0.409899 0.2409)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (1.79101e-16 -1.00849e-16 -8.91094e-16) OK.
    Max cell openness = 5.34516e-16 OK.
    Max aspect ratio = 12.9508 OK.
    Minimum face area = 3.30405e-08. Maximum face area = 1.60733e-05.  Face area magnitudes OK.
    Min volume = 5.95973e-12. Max volume = 3.00018e-08.  Total volume = 0.00436354.  Cell volumes OK.
    Mesh non-orthogonality Max: 69.9681 average: 16.6896
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 2.23813 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

End