# Programming questions

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

 November 24, 2009, 15:28 Programming questions #1 Member   Sylvain Martel Join Date: Apr 2009 Location: University of Sherbrooke/Quebec/Canada Posts: 51 Rep Power: 16 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 Zhiheng Wang likes this.

 November 24, 2009, 22:25 #2 Senior Member   Eelco van Vliet Join Date: Mar 2009 Location: The Netherlands Posts: 124 Rep Power: 18 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 mm.abdollahzadeh, songwukong, smraniaki and 7 others like this.

 November 24, 2009, 22:48 correction #3 Senior Member   Eelco van Vliet Join Date: Mar 2009 Location: The Netherlands Posts: 124 Rep Power: 18 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 Zhiheng Wang likes this.

 November 24, 2009, 23:19 #4 Member   Mathieu Olivier Join Date: Mar 2009 Location: Quebec City, Canada Posts: 77 Rep Power: 16 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, Thamali and 1 others like this.

 November 25, 2009, 02:55 #5 Senior Member   Kathrin Kissling Join Date: Mar 2009 Location: Besigheim, Germany Posts: 134 Rep Power: 16 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

 November 25, 2009, 14:49 #6 Member   Sylvain Martel Join Date: Apr 2009 Location: University of Sherbrooke/Quebec/Canada Posts: 51 Rep Power: 16 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

 November 25, 2009, 14:59 #7 Member   Mathieu Olivier Join Date: Mar 2009 Location: Quebec City, Canada Posts: 77 Rep Power: 16 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

 December 3, 2009, 11:28 #8 Senior Member     Matvey Kraposhin Join Date: Mar 2009 Location: Moscow, Russian Federation Posts: 355 Rep Power: 20 .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, mgg and 5 others like this.

January 6, 2010, 08:45
#9
Senior Member

Dr. Alexander Vakhrushev
Join Date: Mar 2009
Posts: 250
Blog Entries: 1
Rep Power: 18
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());
__________________
Best regards,

Dr. Alexander VAKHRUSHEV

Christian Doppler Laboratory for "Metallurgical Applications of Magnetohydrodynamics"

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

http://smmp.unileoben.ac.at

Last edited by makaveli_lcf; January 6, 2010 at 09:14.

January 11, 2010, 15:27
#10
Senior Member

Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 355
Rep Power: 20
Quote:
 Originally Posted by 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 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
 ddtPhiCorr.pdf (35.6 KB, 1125 views)

 January 12, 2010, 01:51 #11 Senior Member   Dr. Alexander Vakhrushev Join Date: Mar 2009 Posts: 250 Blog Entries: 1 Rep Power: 18 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 "Metallurgical Applications of Magnetohydrodynamics" Simulation and Modelling of Metallurgical Processes Department of Metallurgy University of Leoben http://smmp.unileoben.ac.at

 January 29, 2010, 04:44 #12 Senior Member   Claus Meister Join Date: Aug 2009 Location: Wiesbaden, Germany Posts: 241 Rep Power: 16 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

 January 29, 2010, 05:34 #13 Senior Member   Kathrin Kissling Join Date: Mar 2009 Location: Besigheim, Germany Posts: 134 Rep Power: 16 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

 February 5, 2010, 02:15 help #14 Member   mohsen kh Join Date: Nov 2009 Posts: 41 Rep Power: 14 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

 February 5, 2010, 07:30 #15 Senior Member   Dr. Alexander Vakhrushev Join Date: Mar 2009 Posts: 250 Blog Entries: 1 Rep Power: 18 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 "Metallurgical Applications of Magnetohydrodynamics" Simulation and Modelling of Metallurgical Processes Department of Metallurgy University of Leoben http://smmp.unileoben.ac.at

March 20, 2010, 15:27
#16
Member

YS
Join Date: Jan 2010
Posts: 91
Rep Power: 15
Quote:
 Originally Posted by mkraposhin 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::fvcDdtPhiCorr(rA, U, phiAbs). lines 372 - 402 3) OpenFOAM-1.6/src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C, procedure ddtScheme::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!

 April 18, 2010, 08:46 help #17 Member   mohsen kh Join Date: Nov 2009 Posts: 41 Rep Power: 14 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

 April 18, 2010, 08:53 help #18 Member   mohsen kh Join Date: Nov 2009 Posts: 41 Rep Power: 14 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 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 tauEqn ( fvm::ddt(tau_) == 1/4*(moh1+moh2) -1/4*(moh3+moh4) ); tauEqn().relax(); solve(tauEqn); // Evolution of the backbone stretch tmp 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

 June 16, 2010, 12:44 #19 Senior Member   Join Date: Apr 2010 Posts: 151 Rep Power: 15 A little suggestion: sometimes people forget to run "wclean" before "wmake", provoking the error you had.

 December 14, 2011, 21:39 Why grad(U)().T() and not grad(U).T() #20 Member   Rui Vizinho de Oliveira Join Date: Sep 2011 Posts: 31 Rep Power: 13 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