[blais.bruno]

Hello everyone,

I am trying to work a specific subset of the volume averaged navier stokes equation using OpenFOAM. I couple openfoam to a DEM code and use the results from the DEM code (position and velocities of the particles) and use those in the OpenFOAM solver.

However, first, I am trying to verify that i am solving the right set of equation using the method of manufactured solutions.

The equations I want to solve are :

div (U * voidfraction) = 0

dt (voidfraction * U ) + div (voidfraction U U) = - grad P + div tau + G

To use the method of manufactured solution, I inject a velocity field U into mathematica and calculate the necessary source term G to compensate the equation so that the solution works.

My problem is this. These equation can have non-divergence free velocity field, in fact only div (voidfraction * U) must be divergence free. I have manufactured solution for the Navier-Stokes equation and for this equation but with a divergence free velocity and in both case I have managed to recoved second-order convergence for BOTH pressure and velocity.

However, as soon as I input a non-divergence free velocity, I get a wrong pressure field.

The code for UEqn is :

            // Momentum predictor
            fvVectorMatrix UEqn
            (
		fvm::ddt(voidfraction,U) 
                +fvm::div(phi, U)
		+particleCloud.divVoidfractionTau(U, voidfraction)
            );

		if(momentumPredictor)
		    solve(UEqn == - fvc::grad(p)  + G);

The code for the PISO loop is :

            for (int corr=0; corr<nCorrSoph; corr++)
            {
                volScalarField rUA = 1.0/UEqn.A();

                surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA*voidfraction));
                volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);

		U = rUA*UEqn.H();

                phi = (fvc::interpolate(U*voidfraction) & mesh.Sf() )
                     + fvc::ddtPhiCorr(rUA, U, phi);

                // Non-orthogonal pressure corrector loop
                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
                {
                    // Pressure corrector
                    fvScalarMatrix pEqn
                    (
                        fvm::laplacian(rUAvoidfraction, p) ==  fvc::div(phiGes) 
			+fvc::div(rUAvoidfraction * G)
                    );
                    pEqn.setReference(pRefCell, pRefValue);
                    if
                    (
                        corr == nCorr-1
                     && nonOrth == nNonOrthCorr 
                    )
                    {
                        pEqn.solve(mesh.solver("pFinal"));
                    }
                    else
                    {
                        pEqn.solve();
                    }

                    if (nonOrth == nNonOrthCorr)
                    {
                        phiGes -= pEqn.flux();
                     }

                } // end non-orthogonal corrector loop

                #include "continuityErrorPhiPU.H"

                    U -= rUA*fvc::grad(p) - G * rUA ;
        }

Obviously, something is wrong with my PISO loop, but I cannot figure out what is wrong exactly. I must say I am a bit confused where i should go with rUA and rUAvoidfraction. I have read Jasak thesis, but this has not enlighted me clearly on this issue.

Does anyone have a clue?

Thank you very much for both your time and your help.
Best regards,