[kk415]

Hello All,

I am trying to develop a multiphase flow solver where the concept is to solve the advection equation of VOF on a finer mesh than the Navier Stokes mesh. So as a very basic problem I am solving the RayleighPlateau Instability problem with zero surface tension coefficient.

The finer mesh and coarser mesh are also kept identical currently. So all I am doing is copying the coarse mesh field of Velocity 'U' and Flux 'phi' from coarse mesh to Finer Mesh. 'UF' is identical to 'U' and 'phiF' is identical to 'phi'.

{
Info<< "Interpolating Flux Field from Coarser to Finer Mesh using Gradient method" << endl;
volTensorField gradU(fvc::grad(U));
forAll(LSMesh.cells(), cellI2){ 
UF[cellI2].component(0)=U[tmpList[cellI2]].component(0)+gradU[tmpList[cellI2]].component(0)*r[cellI2].component(0)+gradU[tmpList[cellI2]].component(1)*r[cellI2].component(1)+gradU[tmpList[cellI2]].component(2)*r[cellI2].component(2);
UF[cellI2].component(1)=U[tmpList[cellI2]].component(1)+gradU[tmpList[cellI2]].component(3)*r[cellI2].component(0)+gradU[tmpList[cellI2]].component(4)*r[cellI2].component(1)+gradU[tmpList[cellI2]].component(5)*r[cellI2].component(2);
UF[cellI2].component(2)=U[tmpList[cellI2]].component(2)+gradU[tmpList[cellI2]].component(6)*r[cellI2].component(0)+gradU[tmpList[cellI2]].component(7)*r[cellI2].component(1)+gradU[tmpList[cellI2]].component(8)*r[cellI2].component(2);
}
//phiF = linearInterpolate(UF) & LSMesh.Sf();
//Info << phi << endl;
forAll(phiF(), iFace){     //Method2
//Info << LSMesh.Cf()[iFace] << tab << mesh.Cf()[iFace] << endl;
    //const label& leftcellFine = LSMesh.owner()[iFace];
    //const label& rightcellFine = LSMesh.neighbour()[iFace];
    //if (WVolList[leftcellFine] < 1){
    //    phiF[iFace] = 1e-20;
    //}
    //else {
        phiF[iFace] = phi[iFace];
    //}
}

/*forAll (LSMesh.boundaryMesh() , ipatch ) {
    if( (LSMesh.boundaryMesh()[ipatch].type() == "empty") || (LSMesh.boundaryMesh()[ipatch].type() == "symmetry")){
    }
    else{
        const polyPatch& cPatch = LSMesh.boundaryMesh()[ipatch];
        forAll(cPatch, facei){
        label iFaceGlobal = facei + LSMesh.boundaryMesh()[ipatch].start();
        phiF[iFaceGlobal] = phi[iFaceGlobal];
        //Info << iFaceGlobal << tab << phi[iFaceGlobal] << tab << phiF[iFaceGlobal] << endl;
        }
    }
}*/

}

r is a zero vector currently. The resultant alpha is field is brought back from the fine mesh to coarse mesh as volume average using the following code

if(LSMesh.cells().size()>=mesh.cells().size())
{
//gradPsi == 0*gradPsi;
alpha1 == 0*alpha1;
Info<< "Interpolation from finer to coarser mesh" << endl;
forAll(LSMesh.cells(), cellI){
alpha1[tmpList[cellI]]=alpha1[tmpList[cellI]]+WVolList[cellI]*alpha1F[cellI];}
//C[tmpList[cellI]]=C[tmpList[cellI]]+WVolList[cellI]*C_f[cellI];}
}

WVolList is a unit scalar here and tmpList[cellI]=cellI.

Still the result is different in the current solver from both isoAdvector and MULES. The result picture is attached. It seems there is issue with the boundary treatment. Can anyone suggest what possibly be wrong here.

Boundary Condition for alpha:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      alpha.water;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    /*bottom
    {
        type            cyclic;
    }

    top
    {
        type            cyclic;//
    }

    left
    {
        type            cyclic;//
    }
    right
    {
        type            cyclic;//
    }*/
    symmetry
    {
        type            symmetry;
    }
    atmosphere
    {
        type            zeroGradient;
    }
    front
    {
        type            empty;
    }
    back
    {
        type            empty;
    }
}

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

Needless to say the case setup works fine on single grid solver.