[mrishi]

Hi all,


I am trying to setup an algebraic body force distribution with pimpleFoam. Previously, this code had worked quite well when I tried with ANSYS Fluent UDFs. However, when I use fvOptions to add momentum source, I find that velocity generated is in opposite direction to the source Term. A simple case to illustrate my point:

1. mesh.png shows the cylindrical bath I am simulating




2. The region marked for adding source (shown in the cross section)





3. The sourceTerm vector





4. and finally, the resultant velocity field that is generated:



There is a clear inversion of direction going from direction of source to direction of velocity, whose reason I can't quite pinpoint.

Things which seem important:

a. the direction of g is (0 -9.81 0), so -ve y direction. Could this be interfering? The solver is pimpleFoam, where createFields creates a sense of height opposite to gravity through through ghRef.H and gh.H.

b. this body force is a toy case for the actual force distribution, which has gaussian decay. I have used simpler expression to verify this "misdirection". Here I have applied a basic constant source term as vectorCodedSource as follows:

codeAddSup
#{
  vectorField& momSource = eqn.source(); //momentum source pointer

    const Time& runTime = mesh_.time();
        const vectorField& C = mesh_.C();
        const scalarField& V = mesh_.V();

//termSource creation with r/w
       volVectorField  termSource
       (
           IOobject
           (
               "termSource",
               runTime.timeName(),
               mesh_,
               IOobject::NO_READ,
               IOobject::AUTO_WRITE
           ),
           mesh_,
           dimensionedVector
           (
               "source",
               dimVelocity,
               vector::zero
           )
        );


         volScalarField plumeMarker
        (
             IOobject
             (
                "plumeMarker",
                runTime.timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            mesh_,
            dimensionedScalar
            (
                 "marker",
                 dimless,
                 scalar(0)
            )
        );

//cylindrical polar coordinates
        const scalarField cellRadius = sqrt(pow(C.component(0),2)+pow(C.component(2),2)); 
        const scalarField cellTan = C.component(1)/cellRadius;
        const scalarField cellHeight = C.component(1);

scalarField Rp = (cellHeight)*0.176; //threshold to mark region of force action

//loop to set values
  forAll(C, i)
  {

      if(cellRadius[i]<Rp[i])
      {
          sourceLocal[i].component(0) = 0.0;
          sourceLocal[i].component(1) = 998.1*9.81; //use of density is redundant, gets divided out later, but keeping it in spirit of the actual algebraic Body force term.
          sourceLocal[i].component(2) = 0.0;


          plumeMarker[i] = 0; //IOobject for postProcess
          termSource[i] = sourceLocal[i]/998.1;
      }
      else
      {
          sourceLocal[i] = vector(0,0,0);
          plumeMarker[i] = 1;
      }
   momSource += sourceLocal/998.1;
   if(runTime.outputTime())
   {
       plumeMarker.write();
       termSource.write();
   }


  #};

}

So effectively, in the central region there is zero net gravity while everywhere else, it is downward. And yet the flow seems to go the other way round. To me it seems an issue of wrong direction specification when the fields are declared, but if i think in terms of global coordinate system, the terms appear to be oriented alright.

I really appreciate help on figuring out what went wrong in the process.

[mrishi]

I must add that with reference to above, even with changing the direction of g vector, the reversal from source to actual velocity persists.
The way this source term is added to UEqn is:

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

so sign reversal of source should not happen considering that it is on RHS.

[mrishi]

In the documentation of codedSource,



https://cpp.openfoam.org/v6/classFoa...e.html#details



the example shows a sample energySource code, where the sourceTerm is SUBTRACTED instead of adding. Is that really the key to understanding it? It seems very counter-intuitive.

Perhaps I should try and play with some scalar and see how this affects things, does positive source reduce and vice versa?

 codeAddSup
         #{
             const Time& time = mesh().time();
             const scalarField& V = mesh_.V();
             scalarField& heSource = eqn.source();
             heSource -= 0.1*sqr(time.value())*V;
         #};