OpenFOAM solvers

From cpl-wiki
Jump to: navigation, search

The OpenFOAM Solver

OpenFOAM is a collection of numerical libraries, which can be used to write a top level solver to do whatever you want. The coupling of OpenFOAM through CPL library maintains this philosophy, by simply including a CPL::send and CPL::recv command in a custom solver which is designed for the type of coupling you want. The only code needed to coupled in the simplest case is,

    CPL.pack(U, p, nu, mesh, CPL.STRESS);
    CPL.send();
    CPL.recvVelocity();
    CPL.unpackVelocity(U, mesh);

The pack command accepts a range of arguments as bitflags, for example to pack up stress and velocity, we could use "CPL.stress | CPL.vel". All possible values which can be obtained from velocity U and pressure p can be packed (see the source code in CPLSocketFOAM). The user is encouraged to simply start from one of the example solvers in this section (or the OpenFOAM tutorials), manipulate as needed to solve whatever you want and send/receive. If the default pack and unpack command to the socket code are not sufficient, it should be straight forward to add new ones as needed. Note the main restriction that we assume a uniform grid for all exchanged information.


The Coupled IcoFOAM Solver

The simplest solver takes icoFOAM and adds in a CPL_send and CPL_recv operation. Some work is still required in the CPLSocket to pack stress of velocity at the correct location and to unpack the required values into the boundary cells.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "pisoControl.H"
#include "CPLSocketFOAM.H"

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

int main(int argc, char *argv[])
{

    CPLSocketFOAM CPL;
    CPL.initComms(argc, argv);

    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    pisoControl piso(mesh);

    #include "createFields.H"
    #include "initContinuityErrs.H"

    // MPI_Init is called somewhere in the PStream library
    CPL.initCFD(runTime, mesh);

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

	// Initial communication to initialize domains
    CPL.pack(U, p, nu, mesh, CPL.STRESS);
    CPL.send();
    CPL.recvVelocity();
    CPL.unpackVelocity(U, mesh);

    Info<< "\nStarting time loop\n" << endl;
    while (runTime.loop())
    {

        //Pack and send to MD
        if (CPL.Stresscoupling) 
        {
            CPL.pack(U, p, nu, mesh, CPL.VEL);
        } elseif (CPL.velocitycoupling) 
        {
            CPL.pack(U, p, nu, mesh, CPL.STRESS);
        }
        CPL.send();

        //Receive and unpack from MD
        CPL.recvVelocity();
        CPL.unpackVelocity(U, mesh);

        #include "CourantNo.H"

        // Momentum predictor
        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );

        if (piso.momentumPredictor())
            solve(UEqn == -fvc::grad(p));

        // --- PISO loop
        while (piso.correct())
        {
            volScalarField rAU(1.0/UEqn.A());

            volVectorField HbyA("HbyA", U);
            HbyA = rAU*UEqn.H();
            surfaceScalarField phiHbyA
            (
                "phiHbyA",
                (fvc::interpolate(HbyA) & mesh.Sf())
              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
            );

            adjustPhi(phiHbyA, U, p);

            // Non-orthogonal pressure corrector loop
            while (piso.correctNonOrthogonal())
            {
                // Pressure corrector
                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                );

                pEqn.setReference(pRefCell, pRefValue);
                pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));

                if (piso.finalNonOrthogonalIter())
                    phi = phiHbyA - pEqn.flux();

            }

            #include "continuityErrs.H"

            U = HbyA - rAU*fvc::grad(p);
            U.correctBoundaryConditions();
        }

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }
    Info<< "End\n" << endl;
	CPL::finalize();

    return 0;
}

Otherwise the solver proceeds much as if coupling is not applied.

The SediFOAM Solver

Based on lammpsFoam from the SediFOAM project [1]. The system solved is the incompressible fluid phase and a second phase which is the particles in the liquid modeled with DEM (LAMMPS).

The actual coupling all occures through a CPL_Send and CPL_recv interface.

The SediFOAM code is based on bubbleFoam/twoPhaseEulerFoam, described in a forum exchange [1] using the phase intensive formulation of the momentum equation [2]. The details of the implementation in OpenFOAM(r) are described in an internal report of OpenCFD [3] but are apparently summed up in Henrik Rusche PhD thesis [2].

The idea behind the solution algorithm is the following

- momentum predictor without pressure gradient, solve one Jacobi iteration to get a guess of velocity.
- move gravity and part of the drag term to the pressure solver (known as semi-implicit coupling).
- solve the poisson equation to get the velocity correction which gives the required pressure but do not correct the velocity field directly, instead correct the combined flux for both phases of the system.
- obtain the velocity correction from a reconstruction of the flux correction.

The last two steps are the key point, if you want to have a stable solution when sharp interfaces with large density gradients are present.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "CPLSocketFOAM.H"

int main(int argc, char *argv[])
{

    //This Command turns off solver output 
    solverPerformance::debug=0;

    // Create a CPL object (not used if uncoupled)
    // and intialise MPI
    CPLSocketFOAM CPL;
    MPI_Init(&argc, &argv);

    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "readEnvironmentalProperties.H"
    #include "createFields.H"
    #include "readPISO.H"
    #include "initContinuityErrs.H"

    if (coupled)
        CPL.initComms(argc, argv);

    // Also update/create MD-related fields.
    beta = scalar(1) - alpha;
    volScalarField dragCoef = alpha*dimensionedScalar("dum", dimensionSet(1, -3, -1, 0, 0), 0.0);

	// MPI_Init is called somewhere in the PStream library
    if (coupled)
        CPL.initCFD(runTime, mesh);

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.run())
    {
        runTime++;
        
        if (runTime.outputTime())
            Info<< "Time = " << runTime.timeName() << endl;

        if (coupled){
            //Packup and send velocity, gradident of pressure and divergence of stress
            CPL.pack(Ub, p, nub, mesh, CPL.VEL | CPL.GRADPRESSURE | CPL.DIVSTRESS);
            CPL.send();

            //Recieve and unpack particle velocity, force, 
            //sum of force weightings and porosity
            CPL.recv();
            CPL.unpackPorousVelForceCoeff(Ua, F, dragCoef, beta, maxPossibleAlpha, mesh);
            alpha = scalar(1) - beta;
        }

        fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);
        betaf = fvc::interpolate(beta);
        betaPhib = betaf*phib;

        // See H. Xiao and J. Sun / Commun. Comput. Phys., 9 (2011), pp. 297-323 
        // For explaination of various terms omega and A from Cloud
        // \sum DragCoef*U_b where U_b is fluid velocity 
        UbEqn =
        (
            fvm::ddt(beta, Ub)
          + fvm::div(betaPhib, Ub, "div(phib,Ub)")
          - fvm::Sp(fvc::ddt(beta) + fvc::div(betaPhib), Ub)
          // This term is different in Anderson & Jackson or Kafui et al: e*du VS. d(e*u)
          - fvm::laplacian(nub*beta, Ub)
         ==
          - beta*fvm::Sp(dragCoef/rhob, Ub)   // Implicit drag transfered to p-equation
        );


        // E.S. A full solve seems to be needed here in place of relax to give correct answer
        // for plain Couette flow solver
        solve(UbEqn == -fvc::grad(p));

        // Only a Jacobi iteration is done to find a guess of the velocity 
        // before solving for the pressure equation based on the mixture
        //UbEqn.relax();

        // --- PISO loop
        volScalarField rUbA = 1.0/UbEqn.A()*beta;

        // Iterate over number of nCorr specified by PISO input
        for (int corr = 0; corr < nCorr; corr++)
        {
            surfaceScalarField alphaf = fvc::interpolate(alpha);
            surfaceScalarField betaf = scalar(1) - alphaf;
            surfaceScalarField rUbAf = fvc::interpolate(rUbA);
            Ub = rUbA*UbEqn.H()/beta;

            // The gravity and explicit part of drag are moved to the
            // pressure equation (this is known as semi-implicit coupling)
            surfaceScalarField phiDragb = fvc::interpolate(rUbA/rhob) 
                                         *(fvc::interpolate(F) & mesh.Sf())
                                         + rUbAf*(g & mesh.Sf());  
            forAll(p.boundaryField(), patchi)
            {
                if (isA<zeroGradientFvPatchScalarField>(p.boundaryField()[patchi]))
                    phiDragb.boundaryField()[patchi] = 0.0;
            }
            Ua.correctBoundaryConditions();

            // Solve the pressure equation to enforce mass conservation for the mixture
            phia = (fvc::interpolate(Ua) & mesh.Sf());
            phib = (fvc::interpolate(Ub) & mesh.Sf())
                  + rUbAf*fvc::ddtCorr(Ub, phib)
                  + phiDragb;
            phi = alphaf*phia + betaf*phib;
            surfaceScalarField Dp("(rhob*(1|A(U)))", betaf*rUbAf/rhob);
            fvScalarMatrix pEqn(fvm::laplacian(Dp, p) == fvc::div(phi));
            pEqn.setReference(pRefCell, pRefValue);
            pEqn.solve();

            // Do not correct the velocity field directly, instead correct the flux
            surfaceScalarField SfGradp = pEqn.flux()/Dp;
            phib -= rUbAf*SfGradp/rhob;
            phi = alphaf*phia + betaf*phib;
            p.relax();
            SfGradp = pEqn.flux()/Dp;
            Ub += (fvc::reconstruct(phiDragb - rUbAf*SfGradp/rhob));
            Ub.correctBoundaryConditions();
            U = alpha*Ua + beta*Ub;
    
        }
        Ub.correctBoundaryConditions();

        #include "write.H"

    }

    Info<< "End\n" << endl;

    if (! Pstream::parRun())  MPI_Finalize();
    return(0);
}


References

Template:Reflist
  1. https://www.cfd-online.com/Forums/openfoam-solving/58178-twophaseeulerfoam-documentation-2.html
  2. P. J. Oliveira, R. I. Issa, Numerical aspects of an algorithm for the Eulerian simulation, of two-phase flows, International Journal of Numerical Methods in Fluids, 2003; 43:1177–1198 (DOI: 10.1002/fld.508)
  3. NOT AVAILABLE ANYWHERE