
[Sponsors] 
December 15, 2011, 03:57 

#21 
Senior Member
Claus Meister
Join Date: Aug 2009
Location: Wiesbaden, Germany
Posts: 241
Rep Power: 10 
For the 2nd question:
template <class Cmpt> inline SymmTensor<Cmpt> twoSymm(const Tensor<Cmpt>& t) { return SymmTensor<Cmpt> ( 2*t.xx(), (t.xy() + t.yx()), (t.xz() + t.zx()), 2*t.yy(), (t.yz() + t.zy()), 2*t.zz() ); } Take a look into Tensorl.H in OpenFOAM/OpenFOAM1.7.x/src/OpenFOAM/lnInclude/. 

December 15, 2011, 04:04 

#22 
Senior Member
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,702
Rep Power: 27 
Good morning Rui,
With respect to the first question, then fvc::grad(U) will return a tmp<volTensorField>, however, the transpose ( .T() ) is only defined on volTensorField and not tmp<volTensorField>, thus in order to remove the tmpstatus, you write fvc::grad(U)(). Now you have a volTensorField, and you can perform the operations you would like to do with it. I hope that it clarified things, Niels 

December 15, 2011, 11:02 

#23 
Member
Rui Vizinho de Oliveira
Join Date: Sep 2011
Posts: 31
Rep Power: 7 
Thank you so much too both of you
Best regards Rui 

February 14, 2012, 16:08 

#24 
Member
Matthew J. Churchfield
Join Date: Nov 2009
Location: Boulder, Colorado, USA
Posts: 49
Rep Power: 9 
Matvej,
Thanks for your description of ddtPhiCorr. I went into the code myself and saw the same formulation you see. Do you know what the purpose of ddtPhiCorr is? No one explains its purpose. A few people think it is the RhieChow interpolation, but it is not. It is some sort of time correction, but also compares the old Rhiechow interpolated flux to the linearly interpolated flux. I have removed the ddtPhiCorr from my code and found the code runs similarly, but takes more iterations of the pressure solve to converge. Thank you, Matt 

February 21, 2012, 22:23 

#25 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 16 
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

April 7, 2013, 18:37 

#26 
New Member
Jialin Su
Join Date: Mar 2013
Location: Loughborough
Posts: 29
Rep Power: 5 
Hi Foamer,
I find this thread fairly informative. But I would like to seek clarification for some of my understanding for the following code of piso/simple/ico foam: //set up the linear algebra for the momentum equation. The flux of U, phi, is treated explicity //using the last known value of U. fvVectorMatrix UEqn ( fvm::ddt(U) + fvm::div(phi, U)  fvm::laplacian(nu, U) ); // solve using the last known value of p on the RHS. This gives us a velocity field that is // not divergence free, but approximately satisfies momentum. See Eqn. 7.31 of Ferziger & Peric solve(UEqn == fvc::grad(p)); //  PISO loop take nCorr corrector steps for (int corr=0; corr<nCorr; corr++) { // from the last solution of velocity, extract the diag. term from the matrix and store the reciprocal // note that the matrix coefficients are functions of U due to the nonlinearity of convection. volScalarField rUA = 1.0/UEqn.A(); // take a Jacobi pass and update U. See Hrv Jasak's thesis eqn. 3.137 and Henrik Rusche's thesis, eqn. 2.43 // UEqn.H is the righthand side of the UEqn minus the product of (the offdiagonal terms and U). // Note that since the pressure gradient is not included in the UEqn. above, this gives us U without // the pressure gradient. Also note that UEqn.H() is a function of U. U = rUA*UEqn.H(); I would like to know if rUA and UEqn.H() use the values of U from the previous iteration or the values of U obtained after the momentum predictor? i.e. are UEqn.A() and UEqn.H() automatically updated after "solve(UEqn == fvc::grad(p))" is done? I read a post from Santiago saying that they will be updated. But I would like to confirm again. Thank you. Regards, Callum 

April 8, 2013, 01:21 

#27 
Senior Member

Before you have not assigned a new values for velocity with:
1) U = ... or 2) UEqn.solve(....) you have old iteration values. One can store them explicitly as well with U.storePrevIter() and then call again later with U.prevIter().
__________________
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 FranzJosefStr. 18 A  8700 Leoben Österreich / Austria Tel.: +43 3842  402  3125 http://smmp.unileoben.ac.at Last edited by makaveli_lcf; April 8, 2013 at 01:23. Reason: Additions 

April 8, 2013, 06:37 

#28 
New Member
Jialin Su
Join Date: Mar 2013
Location: Loughborough
Posts: 29
Rep Power: 5 
Hi Alex,
Thank you for the reply. I am new to OF (and C++) and the question below might be a bit dumb. But could you or anyone tell me the difference between "solve(...." alone and "UEqn.solve(..."? Are they indeed the same? (I'd thought so until I saw the line "pEqn.solve(..." in the code)...Thanks. Regards, Jialin Su 

April 8, 2013, 06:44 

#29 
Senior Member

The difference is following:
1) UEqn ( ddt(U) + div(phi, U) == laplacian(nu, U)  grad(p) ); here you already include pressure in you matrix UEqn. And you can solve momentum equation with UEqn.solve(). 2) But suppose you would like to keep pressure influence away from UEqn matrix (like it is done in PISO algorithm): UEqn ( ddt(U) + div(phi, U) == laplacian(nu, U) ); and to solve proper momentum equation (where the pressure gradient should exist) you have to call solve(UEqn ==  grad(p)) thereby you include pressure in momentum equation solution, but keep it away from your linear system UEqn.
__________________
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 FranzJosefStr. 18 A  8700 Leoben Österreich / Austria Tel.: +43 3842  402  3125 http://smmp.unileoben.ac.at Last edited by makaveli_lcf; April 8, 2013 at 06:45. Reason: Mods 

April 8, 2013, 09:45 

#30 
New Member
Jialin Su
Join Date: Mar 2013
Location: Loughborough
Posts: 29
Rep Power: 5 
Danke sehr, Alex. Jetz verstehe ich.
Freundliche Gruesse Jialin Su 

April 13, 2013, 10:05 

#31  
Member
Avdeev Evgeniy
Join Date: Jan 2011
Location: Togliatty, Russia
Posts: 58
Blog Entries: 1
Rep Power: 13 
Hello,
I have one more question around terms .A(), .H(), .T(). As said before Quote:
I found them in article http://openfoamwiki.net/index.php/Co...lockedCellFoam Code:
... scalarField& lower_ = UEqn.lower(); scalarField& upper_ = UEqn.upper(); scalarField& diag_ = UEqn.diag(); ... Thanks. 

April 13, 2013, 13:44 

#32 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 16 
Hi, A()=diag()/mesh.V()
See fvMatrix.C Code:
00721 template<class Type> 00722 Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const 00723 { 00724 tmp<volScalarField> tAphi 00725 ( 00726 new volScalarField 00727 ( 00728 IOobject 00729 ( 00730 "A("+psi_.name()+')', 00731 psi_.instance(), 00732 psi_.mesh(), 00733 IOobject::NO_READ, 00734 IOobject::NO_WRITE 00735 ), 00736 psi_.mesh(), 00737 dimensions_/psi_.dimensions()/dimVol, 00738 zeroGradientFvPatchScalarField::typeName 00739 ) 00740 ); 00741 00742 tAphi().internalField() = D()/psi_.mesh().V(); 00743 tAphi().correctBoundaryConditions(); 00744 00745 return tAphi; 00746 }
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

May 28, 2013, 04:32 

#33 
New Member
Emad Tandis
Join Date: Sep 2010
Posts: 27
Rep Power: 8 
Hello
How does OpenFOAM work when This function is called: UEqn.A() I mean if phi or U was changed, UEqn.A() changes? or this function depend on mesh and transportProperty? 

July 4, 2013, 05:42 

#34  
Senior Member
Dongyue Li
Join Date: Jun 2012
Location: Torino, Italy
Posts: 742
Rep Power: 9 
Quote:
Thanks 

July 4, 2013, 05:45 

#35  
Senior Member
Dongyue Li
Join Date: Jun 2012
Location: Torino, Italy
Posts: 742
Rep Power: 9 
Quote:
If I was wrong correct me. 

July 9, 2013, 09:21 

#36 
Member
Avdeev Evgeniy
Join Date: Jan 2011
Location: Togliatty, Russia
Posts: 58
Blog Entries: 1
Rep Power: 13 
Santiago, propably nooby question, anyway
from fvMatrix.C Code:
00742 tAphi().internalField() = D()/psi_.mesh().V(); A()=diag()/mesh.V() or another string? And I tryed on practice and got something more like A()~diag()/mesh.V(), i.e. approximately. Is it right? Thank you. 

Thread Tools  
Display Modes  


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 