February 12, 2019, 19:49
|
Trouble understanding the terms of UEqns.H in twoPhaseEulerFoam
|
#1
|
Member
Join Date: Apr 2016
Posts: 30
Rep Power: 6
|
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
Code:
//#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);
}
}
|
|
|