[shanvach]

Hi all,

I am trying to add variable source terms in the twoPhaseEulerFoam. However I am unable to understand the various terms of the momentum equation.

I am trying to add two force terms fsid and fl in the momentum equation. However don't know where to add the terms. Could somebody explain what the bold terms physically mean? I also cannot find the laplacian term for viscous stresses.

Thanks and regards,
Shantanu

//#include "Random.H"
//#include "StochasticDispersionRAS.H"
#include "contErrs.H"
Info<< "Constructing momentum equations" << endl;

MRF.correctBoundaryVelocity(U1);
MRF.correctBoundaryVelocity(U2);
MRF.correctBoundaryVelocity(U);

fvVectorMatrix U1Eqn(U1, rho1.dimensions()*U1.dimensions()*dimVol/dimTime);
fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime);
dimensionedScalar dimen 
( 
"dimen", 
dimensionSet(0,0,1,0,0,0,0), 
scalar(1.0) 
); 
dimensionedScalar dimen1 
( 
"dimen1", 
dimensionSet(0,0,-1,0,0,0,0), 
scalar(1.0) 
); 
dimensionedScalar dimen2 
( 
"dimen1", 
dimensionSet(0,0,-1,0,0,0,0), 
scalar(1.0) 
); 
volVectorField gr_alpha = fvc::grad(alpha2);
volTensorField grU_2 = (fvc::grad(U2))*dimen;
volScalarField sc = grU_2 && (grU_2 + grU_2.T());
const scalar csid = 25;

const scalar a  = 4e-6;
//const scalar k = .1;
const scalar eps = 0.206;
//Lift force
//dimensionedScalar up = ("up", [0 0 -1 0 0 0 0], 1.0);
//const scalar up  =1.0;
volScalarField up = .1629+ (.2088/(.1476+(Foam::sqrt(sc))));
//volScalarField up = /*1.0*sc/(sc+ROOTSMALL)*/ sc*dimen1;
const scalar cl = 1.2;

volVectorField fsid = U2*dimen1; // copy of the field;
volVectorField fl = U2*dimen1;
volScalarField dis = wallDist(mesh).y();
volVectorField ncap = wallDist(mesh).n();
//Info<<"data="<<dis<<"\n"<<endl;
//Info<<"ncap="<<ncap<<"\n"<<endl;

Random rnm;
forAll(fsid, cellI)
{
   // Random :: Random(){}
   // Random rnm;
    //Random rnm Random.rndGen();
    //cachedRandom& rnm = this->owner().rndGen();
    
    scalar rn = rnm.GaussNormal<scalar>();
	//Info<<rn<<endl;
    //volScalarField dis = wallPoint.data_(); 
    //wallPointData<scalar>   wpd;
    //scalar dis = wpd.updateCell();
    //wallDist wpd;
   // volScalarField dis = wallDist(mesh).y();
   // Info<<"data="<<dis<<"\n"<<endl;
    
    fsid[cellI][0] = (Foam::sqrt(sc[cellI]))*alpha2[cellI]*a*eps*rn*csid*-((gr_alpha[cellI][0])/(mag((gr_alpha[cellI]))+ROOTSMALL));//*dimen;
    fsid[cellI][1] = (Foam::sqrt(sc[cellI]))*alpha2[cellI]*a*eps*rn*csid*-((gr_alpha[cellI][1])/(mag((gr_alpha[cellI]))+ROOTSMALL));//*dimen;
    fsid[cellI][2] = (Foam::sqrt(sc[cellI]))*alpha2[cellI]*a*eps*rn*csid*-((gr_alpha[cellI][2])/(mag((gr_alpha[cellI]))+ROOTSMALL));//*dimen;
}
forAll(fl,cellI)
{
	fl[cellI][0] = ((cl*(up[cellI])*(sc[cellI])*a*a*a)/((dis[cellI])+ROOTSMALL))*(-ncap[cellI][0]);//*dimen1;
	fl[cellI][1] = ((cl*(up[cellI])*(sc[cellI])*a*a*a)/((dis[cellI])+ROOTSMALL))*(-ncap[cellI][1]);//*dimen1;
	fl[cellI][2] = ((cl*(up[cellI])*(sc[cellI])*a*a*a)/((dis[cellI])+ROOTSMALL))*(-ncap[cellI][2]);//*dimen1;
}
//Info<< "fsid = "<< fsid << " s\n\n" << endl;

volScalarField Kd(fluid.Kd());

{
    volScalarField Vm(fluid.Vm());

    {
        U1Eqn =
        (
            fvm::ddt(alpha1, rho1, U1) + fvm::div(alphaRhoPhi1, U1)
          - fvm::Sp(contErr1, U1)
          + MRF.DDt(alpha1*rho1 + Vm, U1)
          + phase1.turbulence().divDevRhoReff(U1)
         ==
		    //+ fvm::Sp(fluid.Kd()/rho1, U1)

          - Vm
           *(
                fvm::ddt(U1)
              + fvm::div(phi1, U1)
              - fvm::Sp(fvc::div(phi1), U1)
              - DDtU2 
            )
          + fvOptions(alpha1, rho1, U1) 
        );
        U1Eqn.relax();
        U1Eqn += fvm::Sp(Kd, U1);
        fvOptions.constrain(U1Eqn);
        U1.correctBoundaryConditions();
        fvOptions.correct(U1);
    }

    {
        U2Eqn =
        (
            fvm::ddt(alpha2, rho2, U2) + fvm::div(alphaRhoPhi2, U2)
          - fvm::Sp(contErr2, U2)
          + MRF.DDt(alpha2*rho2 + Vm, U2)
          + phase2.turbulence().divDevRhoReff(U2)
         ==
          - Vm
           *(
                fvm::ddt(U2)
              + fvm::div(phi2, U2)
              - fvm::Sp(fvc::div(phi2), U2)
              - DDtU1 //- fsid*dimen*dimen1 - fl*dimen*dimen1
            //-  fvm::Sp(fluid.Kd()/rho1, U1)
            )
          + fvOptions(alpha2, rho2, U2) //- fsid
        );
        U2Eqn.relax();
        U2Eqn += fvm::Sp(Kd, U2);
        fvOptions.constrain(U2Eqn);
        U2.correctBoundaryConditions();
        fvOptions.correct(U2);
    }
}