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

Spatial source term

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

Like Tree1Likes
  • 1 Post By Cyp

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   April 9, 2018, 04:24
Question Spatial source term
  #1
Member
 
Join Date: Dec 2012
Posts: 33
Rep Power: 13
mechkween is on a distinguished road
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:

Code:
// 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:

Code:
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!
mechkween is offline   Reply With Quote

Old   April 9, 2018, 06:13
Default
  #2
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 299
Rep Power: 18
Cyp is on a distinguished road
argument.value() only works if argument is a dimensionedVector.

Can you compile if you remove .value() ?
Cyp is offline   Reply With Quote

Old   April 9, 2018, 07:15
Default
  #3
Member
 
Join Date: Dec 2012
Posts: 33
Rep Power: 13
mechkween is on a distinguished road
Quote:
Originally Posted by Cyp View Post
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:
Code:
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:

Code:
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.
mechkween is offline   Reply With Quote

Old   April 9, 2018, 08:46
Default
  #4
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 299
Rep Power: 18
Cyp is on a distinguished road
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
Code:
volVectorField argument = (mesh.C() - ramp_loc);
volVectorField force
(
	"force",
	ramp_beg + 0.5*(1.0+Foam::tanh(mag(argument)))*(ramp_end - ramp_beg)
);
Cyp is offline   Reply With Quote

Old   April 9, 2018, 09:58
Default
  #5
Member
 
Join Date: Dec 2012
Posts: 33
Rep Power: 13
mechkween is on a distinguished road
Quote:
Originally Posted by Cyp View Post
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
Code:
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:

Code:
--> 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.
mechkween is offline   Reply With Quote

Old   April 9, 2018, 10:02
Default
  #6
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 299
Rep Power: 18
Cyp is on a distinguished road
it probably means that your formula is not correct. "argument" should be dimensionless to be apply to tanh().
Cyp is offline   Reply With Quote

Old   April 10, 2018, 21:27
Default
  #7
Member
 
Join Date: Dec 2012
Posts: 33
Rep Power: 13
mechkween is on a distinguished road
Quote:
Originally Posted by Cyp View Post
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?
mechkween is offline   Reply With Quote

Old   April 11, 2018, 05:42
Default
  #8
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 299
Rep Power: 18
Cyp is on a distinguished road
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.
Cyp is offline   Reply With Quote

Old   April 12, 2018, 00:47
Default
  #9
Member
 
Join Date: Dec 2012
Posts: 33
Rep Power: 13
mechkween is on a distinguished road
Quote:
Originally Posted by Cyp View Post
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.

Quote:
Originally Posted by Cyp View Post
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?
mechkween is offline   Reply With Quote

Old   April 12, 2018, 04:43
Default
  #10
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 299
Rep Power: 18
Cyp is on a distinguished road
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:
Code:
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
Code:
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,
Zane likes this.
Cyp is offline   Reply With Quote

Old   April 23, 2021, 15:06
Default
  #11
New Member
 
irwin
Join Date: Apr 2014
Location: New Delhi, India
Posts: 19
Rep Power: 12
irwin is on a distinguished road
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
irwin is offline   Reply With Quote

Reply

Tags
pimplefoam, source term, spatial

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
Source Term due to evaporation in energy transport equation styleworker OpenFOAM Programming & Development 3 September 7, 2022 04:09
[OpenFOAM.org] Patches to compile OpenFOAM 2.2 on Mac OS X gschaider OpenFOAM Installation 136 October 10, 2017 18:25
polynomial BC srv537 OpenFOAM Pre-Processing 4 December 3, 2016 10:07
"parabolicVelocity" in OpenFoam 2.1.0 ? sawyer86 OpenFOAM Running, Solving & CFD 21 February 7, 2012 12:44
UDFs for Scalar Eqn - Fluid/Solid HT Greg Perkins FLUENT 0 October 11, 2000 04:43


All times are GMT -4. The time now is 02:28.