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

Programming questions

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

Like Tree27Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   November 24, 2009, 16:28
Default Programming questions
  #1
Member
 
Sylvain Martel
Join Date: Apr 2009
Location: University of Sherbrooke/Quebec/Canada
Posts: 51
Rep Power: 8
smart is on a distinguished road
Hi all,

I would like to know what is the meaning of the following terms .A(), .H(), .T():

volScalarField AU = UEqn().A()
U = UEqn().H()/AU
fvc::div(nueff()*dev(fvc::grad(U)().T()))

In addition, could somebody tell me what is the meaning of the "dev" function?

Finally, what is the meaning of "&" in the following line:
phi = fvc::interpolate(rho*U) & mesh.Sf();

Is this "phi" in the UEqn is the same as the one in the pEqn?

Thank you very much to help me understand more about these think.

Sylvain
smart is offline   Reply With Quote

Old   November 24, 2009, 23:25
Default
  #2
Senior Member
 
Eelco van Vliet
Join Date: Mar 2009
Location: The Netherlands
Posts: 122
Rep Power: 9
eelcovv is on a distinguished road
Hi Sylvain,

I am not an OpenFoam expert, but just started digging into the code and came op with the same questions. For the first question I can refer to

http://physics.ucsd.edu/students/cou.../Tutorial2.pdf

slide 9, where they summarise

- A() returns the central coefficients of an fvVectorMatrix
- H() returns the H operation source of an fvVectorMatrix
- Sf() returns cell face area vector of an fvMesh
- flux() returns the face flux field from an fvScalarMatrix

This terms are related to the PISO solution algorithm, a short explanation is given in
http://www.ifdmavt.ethz.ch/education...ols_slides.pdf
slide 48.

The & is the mathematic cross product vector operator (usually denoted with X). If you want calculate the flux field through a surface (in this case the mass flow), you should multiply the velocity vector with the vector representing the cell face using a cross product, i.e. phi_f = rho u_f X S_f
As far as I understood phi is always the mass flux field through the cell faces. Please, anybody correct me if I am wrong

Hope this helps,

Regards

Eelco
eelcovv is offline   Reply With Quote

Old   November 24, 2009, 23:48
Default correction
  #3
Senior Member
 
Eelco van Vliet
Join Date: Mar 2009
Location: The Netherlands
Posts: 122
Rep Power: 9
eelcovv is on a distinguished road
Sorry, the mass flux through a surface is of course given by the inner product between the velocity and the surface face vector, therefore & should represent the inner product vector operator. That by the way leaves me with a question, because with the inner product you get a scalar, not a vector. Is phi a scalar field of a vector field ? Perhaps somebody could explain.
cheers,

Eelco
eelcovv is offline   Reply With Quote

Old   November 25, 2009, 00:19
Default
  #4
Member
 
Mathieu Olivier
Join Date: Mar 2009
Location: Quebec City, Canada
Posts: 76
Rep Power: 8
mathieu is on a distinguished road
Hi,

A few fast answers:

1. phi is a surfaceScalarField that represents the mass flux through each faces. It is used in convective terms.

2. The .T() function returns the transpose of a tensor (grad(U) in this case).

3. The dev function returns the deviatoric part of a tensor.

4. Generally, phi is not exactly the same in both equations (momentum and pressure). In the momentum equation, phi is the corrected flux from the previous iteration/timestep. In the pressure equation, phi is the mass flux without the "pressure contribution":

U = rUA*UEqn.H();
phi = fvc::interpolate(U) & mesh.Sf();


Hope this helps,

Mathieu
Bufacchi, Matt_B and Thamali like this.
mathieu is offline   Reply With Quote

Old   November 25, 2009, 03:55
Default
  #5
Senior Member
 
Kathrin Kissling
Join Date: Mar 2009
Location: Besigheim, Germany
Posts: 134
Rep Power: 8
kathrin_kissling is on a distinguished road
Mathieu,

you also wanna try to use the Programmers Guide, which can be found under
http://foam.sourceforge.net/doc/Guid...mmersGuide.pdf

This might help a lot.

Best

Kathrin
kathrin_kissling is offline   Reply With Quote

Old   November 25, 2009, 15:49
Default
  #6
Member
 
Sylvain Martel
Join Date: Apr 2009
Location: University of Sherbrooke/Quebec/Canada
Posts: 51
Rep Power: 8
smart is on a distinguished road
Thank you for answers!

What is the meaning of A() and other therms? Now, I know that A() returns the central coefficients of an fvVectorMatrix, but can someone tell me more about each term to have a better understanding of each one?

Thank you all,

Sylvain
smart is offline   Reply With Quote

Old   November 25, 2009, 15:59
Default
  #7
Member
 
Mathieu Olivier
Join Date: Mar 2009
Location: Quebec City, Canada
Posts: 76
Rep Power: 8
mathieu is on a distinguished road
Hi Sylvain,

You may find answers in many finite volume CFD books. Furthermore, let me suggest you prof Jasak's thesis where the nomenclature is very similar to the one used in OpenFOAM:

http://powerlab.fsb.hr/ped/kturbo/Op...jeJasakPhD.pdf

Regards,

Mathieu
mathieu is offline   Reply With Quote

Old   December 3, 2009, 12:28
Default
  #8
Senior Member
 
mkraposhin's Avatar
 
Matvej Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 172
Rep Power: 8
mkraposhin is on a distinguished road
.A() - diagonal elements of matrix
.H() - off-diagonal

ddt(rho U) + div(rho*U*U) - div(nu grad(U)) = - grad(p)

A() U - H() = - grad(p)

.T() - transpose of matrix

& - scalar product
Bufacchi, songwukong and Thamali like this.
mkraposhin is offline   Reply With Quote

Old   January 6, 2010, 09:45
Default
  #9
Senior Member
 
Dr. Alexander Vakhrushev
Join Date: Mar 2009
Posts: 213
Rep Power: 10
makaveli_lcf is on a distinguished road
Send a message via ICQ to makaveli_lcf
Hi all!

Could some one please give more detailed description for flux correction in PISO loop:

phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi);

How do we get part denoted as "fvc::ddtPhiCorr(rUA, U, phi)"?

My question is regarding implicit pressure resistance for porous media, because there for class porousZone we use tmp<volTensorField> trTU representing rUA in explicit case.
How ddtPhiCorr method should be modified to be able to use volTensorField instead of volScalarField?

Ok, I already have found in http://foam.sourceforge.net/doc/Doxy...8H_source.html

Quote:
00013 // ddtPhiCorr not well defined for cases with porosity
00014 phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
So still waiting for your comments concerning this point.
__________________
Best regards,

Dr. Alexander VAKHRUSHEV

Christian Doppler Laboratory for "Advanced Process Simulation of
Solidification and Melting"

Simulation and Modelling of Metallurgical Processes
Department of Metallurgy
University of Leoben

Franz-Josef-Str. 18
A - 8700 Leoben
Österreich / Austria
Tel.: +43 3842 - 402 - 3125
http://smmp.unileoben.ac.at

Last edited by makaveli_lcf; January 6, 2010 at 10:14.
makaveli_lcf is offline   Reply With Quote

Old   January 11, 2010, 16:27
Default
  #10
Senior Member
 
mkraposhin's Avatar
 
Matvej Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 172
Rep Power: 8
mkraposhin is on a distinguished road
Quote:
Originally Posted by makaveli_lcf View Post
Hi all!

Could some one please give more detailed description for flux correction in PISO loop:

phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi);

How do we get part denoted as "fvc::ddtPhiCorr(rUA, U, phi)"?

My question is regarding implicit pressure resistance for porous media, because there for class porousZone we use tmp<volTensorField> trTU representing rUA in explicit case.
How ddtPhiCorr method should be modified to be able to use volTensorField instead of volScalarField?

Ok, I already have found in http://foam.sourceforge.net/doc/Doxy...8H_source.html

So still waiting for your comments concerning this point.
Hi, i'm not sure i can help, but i'll try.
Suggest your flow to be incompressible, take look at this files (numbered by order of calling in program):

1) OpenFOAM-1.6/src/finiteVolume/finiteVolume/fvc/fvcDdt.C, function ddtPhiCorr(rA, U, phi), lines 106-125

2) OpenFOAM-1.6/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C, procedure EulerDdtScheme<Type>::fvcDdtPhiCorr(rA, U, phiAbs). lines 372 - 402

3) OpenFOAM-1.6/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C, procedure ddtScheme<Type>::fvcDdtPhiCoeff(U, phi), lines 138-233

1) Function #1 calls appropriate type of function #2.

2) Function #2 evaluates difference between old value of flux and flux, calculated with old values of velocity. Then, Function #3 is called to calculate value of Kc (see attached file for more explanation)

3) Function #3 evaluates proportion between flux difference (calculated in #2 and given (current) flux values). Proportion should be between some small epsilon and 1 (unity). Evaluated field returned to #2

4) In Function #2, calculated value phi_d multiplied by 1/A, 1/dt and Kc (previous step) and returned.

For more explanations, see attached file
Attached Files
File Type: pdf ddtPhiCorr.pdf (35.6 KB, 737 views)
mkraposhin is offline   Reply With Quote

Old   January 12, 2010, 02:51
Default
  #11
Senior Member
 
Dr. Alexander Vakhrushev
Join Date: Mar 2009
Posts: 213
Rep Power: 10
makaveli_lcf is on a distinguished road
Send a message via ICQ to makaveli_lcf
Thank you Matvej for such detailed explanations! Now this point is clear for me.
All the best!
__________________
Best regards,

Dr. Alexander VAKHRUSHEV

Christian Doppler Laboratory for "Advanced Process Simulation of
Solidification and Melting"

Simulation and Modelling of Metallurgical Processes
Department of Metallurgy
University of Leoben

Franz-Josef-Str. 18
A - 8700 Leoben
Österreich / Austria
Tel.: +43 3842 - 402 - 3125
http://smmp.unileoben.ac.at
makaveli_lcf is offline   Reply With Quote

Old   January 29, 2010, 05:44
Default
  #12
Senior Member
 
Claus Meister
Join Date: Aug 2009
Location: Wiesbaden, Germany
Posts: 241
Rep Power: 8
idrama is on a distinguished road
Hey foamers!

Could anybody tell me please, how to store the matrices A and H into a file? I would like to have a look. Or is there even a better possibility so see them?

Cheers,

Claus
idrama is offline   Reply With Quote

Old   January 29, 2010, 06:34
Default
  #13
Senior Member
 
Kathrin Kissling
Join Date: Mar 2009
Location: Besigheim, Germany
Posts: 134
Rep Power: 8
kathrin_kissling is on a distinguished road
Hello Claus,

have a look into fvMatrix.H its defined there:

volScalarField matrix.A()
tmp< GeometricField< Type,
fvPatchField, volMesh > > matrix.H()


Hope this helps!

Best

Kathrin
kathrin_kissling is offline   Reply With Quote

Old   February 5, 2010, 03:15
Default help
  #14
Member
 
mohsen kh
Join Date: Nov 2009
Posts: 41
Rep Power: 6
mohsenkh599 is an unknown quantity at this point
Hi Foamers
I am an amateur and I want to simulate viscoelastic fluids flow, But I am completely confused with OpenFOAM. And a lot of questions arise for me working with this software. Would you possibly answer some of my questions?
1.For implementing a new model (like LPPT) is it sufficient to write a new application in (OpenFOAM/applications) or it is necessary to modify OpenFOAM/src and/or lib and/or etc.
2. Could you possibly send me one or more new application file which you have made it by yourself (a viscoelastic model for instance). Or a detailed help about implementing new models in OpenFOAM?
I do need your help. Please help. That’s for my master thesis.
Best wishes
Mohsen – m.kh.599@gmail.com
mohsenkh599 is offline   Reply With Quote

Old   February 5, 2010, 08:30
Default
  #15
Senior Member
 
Dr. Alexander Vakhrushev
Join Date: Mar 2009
Posts: 213
Rep Power: 10
makaveli_lcf is on a distinguished road
Send a message via ICQ to makaveli_lcf
Mohsen,

first try to read at least User and Programmer guides, that would help you to avoid such questions.

Good luck!
__________________
Best regards,

Dr. Alexander VAKHRUSHEV

Christian Doppler Laboratory for "Advanced Process Simulation of
Solidification and Melting"

Simulation and Modelling of Metallurgical Processes
Department of Metallurgy
University of Leoben

Franz-Josef-Str. 18
A - 8700 Leoben
Österreich / Austria
Tel.: +43 3842 - 402 - 3125
http://smmp.unileoben.ac.at
makaveli_lcf is offline   Reply With Quote

Old   March 20, 2010, 16:27
Default
  #16
Member
 
YS
Join Date: Jan 2010
Posts: 66
Rep Power: 7
Ya_Squall2010 is on a distinguished road
Quote:
Originally Posted by mkraposhin View Post
Hi, i'm not sure i can help, but i'll try.
Suggest your flow to be incompressible, take look at this files (numbered by order of calling in program):

1) OpenFOAM-1.6/src/finiteVolume/finiteVolume/fvc/fvcDdt.C, function ddtPhiCorr(rA, U, phi), lines 106-125

2) OpenFOAM-1.6/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C, procedure EulerDdtScheme<Type>::fvcDdtPhiCorr(rA, U, phiAbs). lines 372 - 402

3) OpenFOAM-1.6/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C, procedure ddtScheme<Type>::fvcDdtPhiCoeff(U, phi), lines 138-233

1) Function #1 calls appropriate type of function #2.

2) Function #2 evaluates difference between old value of flux and flux, calculated with old values of velocity. Then, Function #3 is called to calculate value of Kc (see attached file for more explanation)

3) Function #3 evaluates proportion between flux difference (calculated in #2 and given (current) flux values). Proportion should be between some small epsilon and 1 (unity). Evaluated field returned to #2

4) In Function #2, calculated value phi_d multiplied by 1/A, 1/dt and Kc (previous step) and returned.

For more explanations, see attached file


Hi, mkraposhin

I was wondering where can we find the mathematical derivation of this treatment or justification? Thanks a lot!
Ya_Squall2010 is offline   Reply With Quote

Old   April 18, 2010, 08:46
Default help
  #17
Member
 
mohsen kh
Join Date: Nov 2009
Posts: 41
Rep Power: 6
mohsenkh599 is an unknown quantity at this point
hi
how can I implement following equation in OpenFoam:
tau_=-2*C*(d(phi)/d(C))
I mean the derivative of phi (a variable) relative to C (a matrix)
the best
mohsenkh599 is offline   Reply With Quote

Old   April 18, 2010, 08:53
Default help
  #18
Member
 
mohsen kh
Join Date: Nov 2009
Posts: 41
Rep Power: 6
mohsenkh599 is an unknown quantity at this point
hello again
when compiling a code I get the following message :

mohsen@mohsen-laptop:~/OpenFOAM/OpenFOAM-1.6/src/transportModels$ ./Allwmake
+ wmake libso incompressible
'/home/mohsen/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libincompressibleTransportModels.so' is up to date.
+ wmake libso interfaceProperties
'/home/mohsen/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libinterfaceProperties.so' is up to date.
+ wmake libso viscoelastic
wmakeLnInclude: linking include files to ./lnInclude
make: *** No rule to make target `viscoelasticLaws/Conf/Conf.dep', needed by `Make/linuxGccDPOpt/dependencies'. Stop.
mohsen@mohsen-laptop:~/OpenFOAM/OpenFOAM-1.6/src/transportModels$


whats the problem?
the code is as follows:

#include "Conf.H"
#include "addToRunTimeSelectionTable.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

defineTypeNameAndDebug(Conf, 0);
addToRunTimeSelectionTable(viscoelasticLaw, Conf, dictionary);

// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

// from components
Conf::Conf
(
const word& name,
const volVectorField& U,
const surfaceScalarField& phi,
const dictionary& dict
)
:
viscoelasticLaw(name, U, phi),
C_
(
IOobject
(
"C" + name,
U.time().timeName(),
U.mesh(),
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
U.mesh()
),
tau_
(
IOobject
(
"tau" + name,
U.time().timeName(),
U.mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
U.mesh(),
dimensionedSymmTensor
(
"zero",
dimensionSet(1, -1, -2, 0, 0, 0, 0),
symmTensor::zero
)
),
I_
(
dimensionedSymmTensor
(
"I",
dimensionSet(0, 0, 0, 0, 0, 0, 0),
symmTensor
(
1, 0, 0,
1, 0,
1
)
)
),
rho_(dict.lookup("rho")),
etaS_(dict.lookup("etaS")),
etaP_(dict.lookup("etaP")),
M_(dict.lookup("M")),
lambda_(dict.lookup("lambda")),
{}


// * * * * * * * * * *t * * * * * Member Functions * * * * * * * * * * * * * //

tmp<fvVectorMatrix> Conf::divTau(volVectorField& U) const
{

dimensionedScalar etaPEff = etaP_;

return
(
fvc::div(tau_/rho_, "div(tau)")
- fvc::laplacian(etaPEff/rho_, U, "laplacian(etaPEff,U)")
+ fvm::laplacian( (etaPEff + etaS_)/rho_, U, "laplacian(etaPEff+etaS,U)")
);

}


void Conf::correct()
{
// Velocity gradient tensor
volTensorField L = fvc::grad( U() );
// Transpose of grad(U)
volTensorField LT = fvc::L.T();
// Twice the rate of deformation tensor
volSymmTensorField twoD = twoSymm( L );
// C . gamma
volTensorField moh1 = twoD & C_;
// gamma . C
volTensorField moh2 = C_ & twoD;
// vorticity
volTensorField vor = L - LT;
//
volTensorField moh3 = vor & C_;
//
volTensorField moh3 = C & vor;


// Evolution of orientation
tmp<fvSymmTensorMatrix> tauEqn
(
fvm::ddt(tau_)
==
1/4*(moh1+moh2)
-1/4*(moh3+moh4)
);

tauEqn().relax();
solve(tauEqn);

// Evolution of the backbone stretch
tmp<fvScalarMatrix> pheeEqn
(
phee
==
(lambda/2)*Foam::sqr(tr((I - inv(C))/2))
+ M*Foam::sqr(tr((I - inv(C))/2))
);

pheeEqn().relax();
solve(pheeEqn);

// Viscoelastic stress
tau_ = -2*C_*fvm::ddt(C);

}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// ************************************************** *********************** //

thank you very much
mohsenkh599 is offline   Reply With Quote

Old   June 16, 2010, 12:44
Default
  #19
Senior Member
 
Join Date: Apr 2010
Posts: 151
Rep Power: 7
flowris is on a distinguished road
A little suggestion: sometimes people forget to run "wclean" before "wmake", provoking the error you had.
flowris is offline   Reply With Quote

Old   December 14, 2011, 22:39
Smile Why grad(U)().T() and not grad(U).T()
  #20
Member
 
Rui Vizinho de Oliveira
Join Date: Sep 2011
Posts: 31
Rep Power: 5
RuiVO is on a distinguished road
Hello Foamers. This thread is a very usefull one indeed. I still have two questions though.
Why is it written :

dev2(fvc::grad(U)().T()) instead of dev2(fvc::grad(U).T())

Why the extra pair of parenthesis ??

The other question is what is the meaning of "twoSymm" in the following expression:

dev(twoSymm(fvc::grad(U_)))

Best Regards

Rui
RuiVO is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
OpenFoam programming prapanj OpenFOAM 10 March 18, 2010 08:23
Programming in OpenFOAM vinu OpenFOAM 2 July 11, 2009 10:16
two questions about udf programming Vincent FLUENT 1 October 5, 2008 00:24
new CFD Programming Forum Thinker Main CFD Forum 14 November 19, 2002 17:03
Programming in C Tony Main CFD Forum 5 March 7, 2002 13:40


All times are GMT -4. The time now is 20:35.