[mechkween]

Hi all,

I've been trying to add a spatial source term to the pimpleFoam application but struggling with the compilation. Hoping someone here can help me out. I need to add a source that ramps from an initial value to a final value. Here's what I did so far:

Defining information I will read from transportProperties:

// Read in ramp data from transport properties
Info<< "Reading ramp properties\n" << endl;

IOdictionary transportProperties
(
IOobject
(
    "transportProperties",
    runTime.constant(),
    mesh,
    IOobject::MUST_READ_IF_MODIFIED,
    IOobject::NO_WRITE
)
);

dimensionedVector ramp_beg
(
transportProperties.lookup("ramp_beg")
);
dimensionedVector ramp_end
(
transportProperties.lookup("ramp_end")
);
dimensionedVector ramp_loc
(
transportProperties.lookup("ramp_loc")
);

// End of ramp reading

Defining the source term and it's addition to UEqn.H:

volVectorField argument = (mesh.C() - ramp_loc);
volVectorField force
(
ramp_beg + 0.5*(1.0+Foam::tanh(argument.value()))*(ramp_end - ramp_beg)
);

// Solve the Momentum equation

tmp<fvVectorMatrix> UEqn
(
    fvm::ddt(U)
  + fvm::div(phi, U)
  + turbulence->divDevReff(U)
  + force // This is the ramped up pressure gradient: +dp/dx
 ==
    fvOptions(U)
);

The source (pressure gradient) is only acting in x direction so I intend to define the following:

• ramp_beg [0 1 -2 0 0 0 0] (0, 0, 0)
• ramp_end [0 1 -2 0 0 0 0] (0.04, 0, 0)
• ramp_loc [0 1 0 0 0 0 0] (10, 0, 0)

At the moment, compiling this says that value (for argument.value()) is not a member of argument.

Any suggestions how to fix it?

Thanks!

[Cyp]

argument.value() only works if argument is a dimensionedVector.

Can you compile if you remove .value() ?

[mechkween]

Originally Posted by Cyp

argument.value() only works if argument is a dimensionedVector.

Can you compile if you remove .value() ?

I tried with dimensionedVector and this is the error message:

UEqn.H:3:40: error: conversion from ‘Foam::tmp<Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> >’ to non-scalar type ‘Foam::dimensionedVector {aka Foam::dimensioned<Foam::Vector<double> >}’ requested
 dimensionedVector argument = (mesh.C() - ramp_loc); 
                                        ^
UEqn.H:6:48: error: no matching function for call to ‘tanh(Foam::Vector<double>&)’
 ramp_beg + 0.5*(1.0+Foam::tanh(argument.value()))*(ramp_end - ramp_beg)
                                                ^

If I stick with volVectorField definition and remove value(), I again get an error message:

UEqn.H:6:40: error: no matching function for call to ‘tanh(Foam::volVectorField&)’
 ramp_beg + 0.5*(1.0+Foam::tanh(argument))*(ramp_end - ramp_beg)
                                        ^

My understanding is that the argument of tanh() must be dimensionless which is why I chose .value() since the actual argument is a vector of dimensions of length i.e. mesh.C() should give me actual coordinates and ramp_loc would be the start of the ramp.

[Cyp]

of course, because "(mesh.C() - ramp_loc) " returns a volVectorField and not a dimensionedVector.

The other problem in your code is the function Foam::tanh( ) that is not defined for a volVectorField (what is an hyperbolic tangent of a vector ?).

Try

volVectorField argument = (mesh.C() - ramp_loc);
volVectorField force
(
	"force",
	ramp_beg + 0.5*(1.0+Foam::tanh(mag(argument)))*(ramp_end - ramp_beg)
);

[mechkween]

Originally Posted by Cyp

of course, because "(mesh.C() - ramp_loc) " returns a volVectorField and not a dimensionedVector.

The other problem in your code is the function Foam::tanh( ) that is not defined for a volVectorField (what is an hyperbolic tangent of a vector ?).

Try

volVectorField argument = (mesh.C() - ramp_loc);
volVectorField force
(
    "force",
    ramp_beg + 0.5*(1.0+Foam::tanh(mag(argument)))*(ramp_end - ramp_beg)
);

Okay, thanks! That did end up compiling but the run crashes with the following:

--> FOAM FATAL ERROR: 
Argument of trancendental function not dimensionless

    From function trans(const dimensionSet&)
    in file dimensionSet/dimensionSet.C at line 430.

which looks like there's still an issue with the tanh() function as that's the only trancendental function in the solver.

[Cyp]

it probably means that your formula is not correct. "argument" should be dimensionless to be apply to tanh().

[mechkween]

Originally Posted by Cyp

it probably means that your formula is not correct. "argument" should be dimensionless to be apply to tanh().

Yes, I managed to non dimensionalise it by dividing by a scalar one with matching dimensions. The ramp is now working.

I was thinking if it is possible to absorb this pressure gradient into the fvc::grad(p) term to calculate for an effective p. The current approach seems like the obtained pressure does not reflect the actual pressure in the system. Would it also affect what the U values would be?

[Cyp]

If you remove fvc::grad(p) from the momentum equation, you can not define an equation for the pressure that satisfies the conservation of the mass (see the description of the PISO algorithm for the pressure-velocity coupling). Basically, in what you did, your velocity is not divergence free.

If you want a ramp of pressure, I will suggest two strategies:

- first, you apply a pressure difference at the boundaries of your domain and you apply a ramp to one end of the domain (this is what people do usually, and also what make sense if you want to compare with some lab experiments)

- Or, you decompose the pressure gradient into a perturbation and a mean value. The mean value is your source term and the pressure gradient correspond to gradient of the perturbation. In this case, you keep fvc::grad(p) (corresponding to the gradient of the perturbation) in the momentum equation, you add your body force (the mean pressure gradient that applies to every cells of the domain) and you solve the pressure equation (you don't have to modify anything). This strategy is usually used for simulations with periodic boundary conditions.

[mechkween]

Originally Posted by Cyp

If you remove fvc::grad(p) from the momentum equation, you can not define an equation for the pressure that satisfies the conservation of the mass (see the description of the PISO algorithm for the pressure-velocity coupling). Basically, in what you did, your velocity is not divergence free.

I haven't removed the fvc::grad(p) term. It's within the solve (Ueqn == -fvc::grad(p)) line, which I didn't show in the first post. I also checked that the divergence of the velocity in the test case is indeed zero everywhere.

Originally Posted by Cyp

If you want a ramp of pressure, I will suggest two strategies:

- first, you apply a pressure difference at the boundaries of your domain and you apply a ramp to one end of the domain (this is what people do usually, and also what make sense if you want to compare with some lab experiments)

- Or, you decompose the pressure gradient into a perturbation and a mean value. The mean value is your source term and the pressure gradient correspond to gradient of the perturbation. In this case, you keep fvc::grad(p) (corresponding to the gradient of the perturbation) in the momentum equation, you add your body force (the mean pressure gradient that applies to every cells of the domain) and you solve the pressure equation (you don't have to modify anything). This strategy is usually used for simulations with periodic boundary conditions.

This part is a little confusing. So you suggest for the first point to define p_inlet and p_outlet as boundary conditions (as uniformFixedValues) such that for the given section length, the gradient is what I need? That would mean no addition of a forcing term in the code?

For the second point, is it the same as I have been doing i.e. adding force to UEqn and solving with fvc::grad(p) term using the PISO algorithm? Or you are suggesting to add the force on the right side of the solve() such that solve (UEqn == -fvc::grad(p) + force) and modify the pEqn to be consistent, like the treatment of gravity in buoyantPimpleFoam?

[Cyp]

Your right for the first point. You consider the ramp only at the inlet boundary.

For the second point, it is important to combine this approach with cyclic boundary condition (so, may be it is not relevant in your case, meaning that you should probably deal with the first option).

Here, again, there are two ways for adding a source term in the momentum equation:

- you add your source term when you build the matrix:

fvVectorMatrix UEqn
(
    fvm::ddt(U) 
+ ...
+ force
);
solve (UEqn == -fvc::grad(p));

in this case, "force" will be consider in the second member of your object "UEqn" and you don't have to change anything in the pressure equation.

- you write

fvVectorMatrix UEqn
(
    fvm::ddt(U) 
+ ...
);
solve (UEqn == -fvc::grad(p)+force);

In this case, you have to modify your pressure-velocity coupling algorithm, i.e. change the pressure equation and the correction stage of the PISO.

The first option is probably the simplest.

Cheers,

[irwin]

I am starting with 2-d pipe flow simulation . I have added source(pressure gradient) with cyclic inlet and outlet in icoFoam solver . I checked it with adding same pressure gradient without cyclic boundary condition at the inlet . Both of them gives completely different results . In fact solver with added source term (pressure gradient ) gives completely different results compared to analytical solution. Can you please suggest something in case I am doing something wrong .




thanks