CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

Trouble understanding the terms of UEqns.H in twoPhaseEulerFoam

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   February 12, 2019, 19:49
Default Trouble understanding the terms of UEqns.H in twoPhaseEulerFoam
  #1
Member
 
Join Date: Apr 2016
Posts: 30
Rep Power: 6
shanvach is on a distinguished road
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);
    }
}
shanvach is offline   Reply With Quote

Reply

Tags
fvm, momentum, source terms, twophaseeulerfoam, ueqn.h

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Understanding Terms in Compressible pEqn AMK53 OpenFOAM Running, Solving & CFD 5 September 18, 2019 02:50
TwoPhaseEulerFoam: trying to understand the terms in UEqns.H ARUN K RAJ OpenFOAM Programming & Development 1 February 19, 2015 06:43
MUSCL scheme: having trouble understanding what total variation diminishing means JMDag2004 OpenFOAM Running, Solving & CFD 2 January 16, 2015 09:10
Question in definition of terms in solve titio OpenFOAM Running, Solving & CFD 0 March 19, 2009 17:02
K-Epsilon model? Brindaban Ghosh Main CFD Forum 2 June 24, 2000 05:22


All times are GMT -4. The time now is 03:58.