[zhangxc0223]

Hi, all

I am struggle in a two-phase pipe model with cyclic boundary conditions. I would like to solve it using DPMFoam. The problem is that it doesn't work if I add a momentumSource in 'fvOptions' file like the following code.

momentumSource
{
    type            meanVelocityForce;
    active          true;

    meanVelocityForceCoeffs
    {
        selectionMode   all;

        fields          (U.air);
        Ubar            (0 -24 0);
    }
}

It gives the warning message as below when running, and the bulk velocity are continuously decreasing from the initial velocity to 0.

--> FOAM Warning : 
    From function virtual void Foam::fv::option::checkApplied() const
    in file cfdTools/general/fvOptions/fvOption.C at line 118
    Source momentumSource defined for field U.air but never used

So, I have two questions about that.

1. is there any way to fix the problem of using 'fvOptions'?
2. is that possible to add a certain value in the 'UcEqn.H' file, to drive the flow with a constant bulk velocity? If so, please give me a hit how to modify it.

Very appreciate for any help.
Regards.

[RobertHB]

Idea 1: Change active true; to active yes; and try again. At least that is my fvOptions setup and its working.

Originally Posted by zhangxc0223

It gives the warning message as below when running, and the bulk velocity are continuously decreasing from the initial velocity to 0.

--> FOAM Warning : 
   [...]
     Source momentumSource defined for field U.air but never used

Idea 2: Are you sure your velocity file is called U.air? Your error clearly states that the field you define your fvOptions for does not exist.

[zhangxc0223]

Originally Posted by RobertHB

Idea 1: Change active true; to active yes; and try again. At least that is my fvOptions setup and its working.

Idea 2: Are you sure your velocity file is called U.air? Your error clearly states that the field you define your fvOptions for does not exist.

Thanks Robert for your reply. I have checked to type 'true', 'yes', and 'on', all of them don't work for it. May I ask what version of OpenFOAM are you using for your fvOptions?

In my DPMFoam solver, the velocity field is U.air for my air phase. That really puzzles me for not exist. I have also tried to use 'U','Uc', still doesn't work.

[RobertHB]

Originally Posted by zhangxc0223

Thanks Robert for your reply. I have checked to type 'true', 'yes', and 'on', all of them don't work for it. May I ask what version of OpenFOAM are you using for your fvOptions?

I'm using OF 5.0 and 6.0. But i think i have found your "problem". Take a look at the UEqn.H (Link) used by e.g. simpleFoam, with which fvSolution works:

tmp<fvVectorMatrix> tUEqn
(
    fvm::div(phi, U)
    + MRF.DDt(U)
    + turbulence->divDevReff(U)
    ==
    fvOptions(U)
 );

and the UcEqn.H (Link) used by DMPFoam

fvVectorMatrix UcEqn
(
    fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
    - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
    + continuousPhaseTurbulence->divDevRhoReff(Uc)
    ==
    (1.0/rhoc)*cloudSU
 );

There is no connection to the fvOptions file. Maybe you can add this connection by adding "+ fvOptions(U)". So that it becomes

fvVectorMatrix UcEqn
(
    fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
    - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
    + continuousPhaseTurbulence->divDevRhoReff(Uc)
    ==
    (1.0/rhoc)*cloudSU
    + fvOptions(U)
);

But i'm not sure about this.

[zhangxc0223]

Originally Posted by RobertHB

So that it becomes

fvVectorMatrix UcEqn
(
    fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
    - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
    + continuousPhaseTurbulence->divDevRhoReff(Uc)
    ==
    (1.0/rhoc)*cloudSU
    + fvOptions(U)
);

But i'm not sure about this.

Thanks Robert, I have try to do following your suggestion, but it gives the error that:

UcEqn.H: In function 'int main(int, char**)':
UcEqn.H:7:38: error: 'fvOptions' was not declared in this scope
     (1.0/rhoc)*cloudSU + fvOptions(Uc)
                                      ^

Here I use the 'Uc' instead of 'U' to avoid that 'U' was not declared in this scope. So What else can I do to fix this? Sorry that I am pretty new on the C++ in OpenFOAM.

[RobertHB]

Try adding #include "fvOptions.H" to your DPMFoam.C

[farzadmech]

Dear Matthew
Did you find your answer? I mean how to use fvOption in DPMfoam Solver?


Thanks,
Farzad

Originally Posted by zhangxc0223

Thanks Robert, I have try to do following your suggestion, but it gives the error that:

UcEqn.H: In function 'int main(int, char**)':
UcEqn.H:7:38: error: 'fvOptions' was not declared in this scope
     (1.0/rhoc)*cloudSU + fvOptions(Uc)
                                      ^

Here I use the 'Uc' instead of 'U' to avoid that 'U' was not declared in this scope. So What else can I do to fix this? Sorry that I am pretty new on the C++ in OpenFOAM.

[dscian]

fvOptions is already added in DPMFoam/MPPICFoam by the following commit:
https://github.com/OpenFOAM/OpenFOAM...714fa03dc67e10

Though it is not in the cell momentum equation, and I don't know why. I guess it should have been added but I don't really know how fvOptions works code-wise.

The options I have encountered in other solvers are; fvOptions(U), fvOptions(rho, U), fvOptions(alpha, rho, U). But I didn't see fvOptions(alpha, U). You might try that one or alphac*fvOptions(U).

Please share the solution if you find out.

Utkan

[farzadmech]

Dear Utkan Caliskan
Thanks for your response. Which OpenFOAM version you are using? I am using version 6.0 and in version 6.0, in DPMFOAM.c we have;

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "PhaseIncompressibleTurbulenceModel.H"
#include "pimpleControl.H"

Also, I am thinking in the Uceqn.C there must be fvOption in the momentum equation, but I cannot find such a thing.



Thanks,
Farzad

Originally Posted by dscian

fvOptions is already added in DPMFoam/MPPICFoam by the following commit:
https://github.com/OpenFOAM/OpenFOAM...714fa03dc67e10

Though it is not in the cell momentum equation, and I don't know why. I guess it should have been added but I don't really know how fvOptions works code-wise.

The options I have encountered in other solvers are; fvOptions(U), fvOptions(rho, U), fvOptions(alpha, rho, U). But I didn't see fvOptions(alpha, U). You might try that one or alphac*fvOptions(U).

Please share the solution if you find out.

Utkan

[farzadmech]

Dear Utkan Caliskan
I did everything to make fvOption work in my simulation. I exactly copied what was in OF8(I mean I have add extra lines to my code). my code compiles fine, but it does not have any effect on the momentum of the specific region which I need.
I think, I will stop using fvOption and I will code instead.



Thanks,
Farzad

Originally Posted by dscian

fvOptions is already added in DPMFoam/MPPICFoam by the following commit:
https://github.com/OpenFOAM/OpenFOAM...714fa03dc67e10

Though it is not in the cell momentum equation, and I don't know why. I guess it should have been added but I don't really know how fvOptions works code-wise.

The options I have encountered in other solvers are; fvOptions(U), fvOptions(rho, U), fvOptions(alpha, rho, U). But I didn't see fvOptions(alpha, U). You might try that one or alphac*fvOptions(U).

Please share the solution if you find out.

Utkan

[dscian]

Can you tell me which modification you made in the solver?

I suggest you add fvOptions in the cell momentum equation (UcEqn).

fvVectorMatrix UcEqn
(
    fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
  - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
  + continuousPhaseTurbulence->divDevRhoReff(Uc)
 ==
    (1.0/rhoc)*cloudSU
  + alphac*fvOptions(U)
);

I've studied its code a bit and think that the way the source added (addSup) depends on the fvOptions options. Moreover, there is no fvOptions(alpha, U). Depending on what you are trying to add, the term above may solve your problem.

Utkan

[farzadmech]

Dear Utkan Caliskan
I did like this;

fvVectorMatrix UcEqn
(
    fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
  - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
  + continuousPhaseTurbulence->divDevRhoReff(Uc)
 ==
    (1.0/rhoc)*cloudSU 

);

UcEqn.relax();

fvOptions.constrain(UcEqn);   //new for momentum purpose

volScalarField rAUc(1.0/UcEqn.A());
surfaceScalarField rAUcf("Dp", fvc::interpolate(rAUc));

surfaceScalarField phicForces
(
   fvc::flux(rAUc*cloudVolSUSu/rhoc) + rAUcf*(g & mesh.Sf())
);

if (pimple.momentumPredictor())
{
    solve
    (
        UcEqn
     ==
        fvc::reconstruct
        (
            phicForces/rAUcf - fvc::snGrad(p)*mesh.magSf()
        )
    );

fvOptions.correct(Uc);    //new for momentum purpose

}

By adding fvOption directly to the momentum equation, it gives me error, so I used what I have found for OF8 for DPMFoam as I wrote by RED in the code, but it did not worked. Also, I have included fvOptions.H to the code.


Thanks,
Farzad

Originally Posted by dscian

Can you tell me which modification you made in the solver?

I suggest you add fvOptions in the cell momentum equation (UcEqn).

fvVectorMatrix UcEqn
(
    fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
  - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
  + continuousPhaseTurbulence->divDevRhoReff(Uc)
 ==
    (1.0/rhoc)*cloudSU
  + alphac*fvOptions(U)
);

I've studied its code a bit and think that the way the source added (addSup) depends on the fvOptions options. Moreover, there is no fvOptions(alpha, U). Depending on what you are trying to add, the term above may solve your problem.

Utkan

[farzadmech]

Dear Utkan Caliskan
Finally I decided to code instead of using fvOption, and I did the coding and I have add externalMomentum to the momentum equation in that specific region like this;

forAll(mesh.C(),i)
{
	scalar x = mesh.C()[i].component(0);
	scalar y = mesh.C()[i].component(1);
	scalar z = mesh.C()[i].component(2);

	
	if(x>0.096 && x<0.404 && y>0.10 && y<0.20 && z>0.0357 && z<0.0643)
	{
	externalMomentum[i] = (1.0 ,0.0 ,0.0);
    }
	
}

But it gives me below error;

ExMom.H:10:22: error: no match for ‘operator=’ (operand types are ‘Foam::Vector<double>’ and ‘double’)
  externalMomentum[i] = (1.0 ,0.0 ,0.0);

by the way, externalMomentum is a volVectorField. Do you know why this problem happens?


Thanks,
Farzad

[dscian]

Sorry for the late reply. The system didn't allow me to send a reply a few days ago.


I was able to recompile it with that addition.

For that one try the following:

forAll(mesh.C(),i)
{
    scalar x = mesh.C()[i].component(0);
    scalar y = mesh.C()[i].component(1);
    scalar z = mesh.C()[i].component(2);

    
    if(x>0.096 && x<0.404 && y>0.10 && y<0.20 && z>0.0357 && z<0.0643)
    {
    externalMomentum[i] = vector(1.0 ,0.0 ,0.0);
    }
    
}

[ktron]

Hi, the approach dscian suggested works for me, with the minimal modification of usig Uc instead of U.

fvVectorMatrix UcEqn
(
    fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
  - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
  + continuousPhaseTurbulence->divDevRhoReff(Uc)
 ==
    (1.0/rhoc)*cloudSU
  + alphac*fvOptions(Uc)
);

[ktron]

I tried to do the same in OpenFOAM 9, but was somewhat confused.


In OF9, the solver DPMFoam is replaced by the solver denseParticleFoam.
Also, fvOptions is replaced by fvModels and fvConstraints.


I am using the momentum source correction "meanVelocityForce", which used to be in fvOptions in OF8. Now, In OF9, it seems to be in fvConstraints, according to:


https://cpp.openfoam.org/v9/meanVelo...8H_source.html


I wanted to try and modify denseParticleFoam in the way that the source correction is done in simpleFoam, which ist

     tmp<fvVectorMatrix> tUEqn
     (
         fvm::div(phi, U)
       + MRF.DDt(U)
       + turbulence->divDevSigma(U)
      ==
         fvModels.source(U)
     );

But of course, this does not include the corrections from fvConstraints. Trying to do something like



+ fvConstraints.source(Uc)



in denseParticleFoam gives the error, that 'class Foam::fvConstraints' has no member 'source'. Simply


+ fvConstraints(Uc)



does not work either, and putting 'meanVelocityForce' in fvModels does not work either. I have to admit, I was just blindly trying some things here, without much knowledge about the background.


However, when running the case with simpleFoam the source term corrections ARE taken into account. I suspect, that there is something to these fat printed parts in denseParticleFoam, which has been mentioned before by farzadmech:

 fvVectorMatrix UcEqn
 (
     fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
   - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
   + continuousPhaseTurbulence->divDevTau(Uc)
  ==
     (1.0/rhoc)*cloudSU
 );
 
 UcEqn.relax();
 
 fvConstraints.constrain(UcEqn);
 
 volScalarField rAUc(1.0/UcEqn.A());
 volScalarField rASpUc(1.0/(UcEqn.A() - cloudSUp/rhoc));
 surfaceScalarField rASpUcf("Dp", fvc::interpolate(rASpUc));
 
 surfaceScalarField phicSUSu
 (
     fvc::flux(rASpUc*cloudSUu/rhoc)
   + rASpUcf*(g & mesh.Sf())
 );
 surfaceScalarField phicSUSp
 (
     fvc::interpolate(rASpUc*cloudSUp/rhoc)
 );
 
 if (pimple.momentumPredictor())
 {
     solve
     (
         UcEqn
      ==
         fvc::reconstruct
         (
             (phicSUSu + phicSUSp*phic)/rASpUcf
           - fvc::snGrad(p)*mesh.magSf()
         )
       + (1.0/rhoc)*(fvm::Sp(cloudSUp, Uc) - cloudSUp*Uc)
     );
 
     fvConstraints.constrain(Uc);
 }

But, since it does not work, I suspect something in 'constain' is somehow not implemented yet for denseParticleFoam. Please excuse the long post. Can anyone provide some clarity?

[ktron]

It also works to do this:

fvVectorMatrix UcEqn
(
    fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
  - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
  + continuousPhaseTurbulence->divDevRhoReff(Uc)
 ==
    (1.0/rhoc)*cloudSU
  + fvOptions(alphac, Uc)
 );

Originally Posted by ktron

Hi, the approach dscian suggested works for me, with the minimal modification of usig Uc instead of U.

fvVectorMatrix UcEqn
(
    fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
  - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
  + continuousPhaseTurbulence->divDevRhoReff(Uc)
 ==
    (1.0/rhoc)*cloudSU
  + alphac*fvOptions(Uc)
);