CFD Online URL
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Programming questions

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

Like Tree26Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   December 15, 2011, 03:57
Default
  #21
Senior Member
 
Claus Meister
Join Date: Aug 2009
Location: Wiesbaden, Germany
Posts: 241
Rep Power: 8
idrama is on a distinguished road
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/OpenFOAM-1.7.x/src/OpenFOAM/lnInclude/.
idrama is offline   Reply With Quote

Old   December 15, 2011, 04:04
Default
  #22
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Rotterdam, The Netherlands
Posts: 1,553
Rep Power: 23
ngj will become famous soon enoughngj will become famous soon enough
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 tmp-status, 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
Bufacchi, Tobi and pfhan like this.
ngj is offline   Reply With Quote

Old   December 15, 2011, 11:02
Default
  #23
Member
 
Rui Vizinho de Oliveira
Join Date: Sep 2011
Posts: 30
Rep Power: 5
RuiVO is on a distinguished road
Thank you so much too both of you

Best regards

Rui
RuiVO is offline   Reply With Quote

Old   February 14, 2012, 16:08
Default
  #24
Member
 
Matthew J. Churchfield
Join Date: Nov 2009
Location: Boulder, Colorado, USA
Posts: 49
Rep Power: 8
mchurchf is on a distinguished road
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 Rhie-Chow interpolation, but it is not. It is some sort of time correction, but also compares the old Rhie-chow 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
mchurchf is offline   Reply With Quote

Old   February 21, 2012, 22:23
Default
  #25
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
santiagomarquezd will become famous soon enough
Matrix manipulation

fvc::ddtPhiCorr(rUA, U, phi)

Momentum predictor


Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   April 7, 2013, 19:37
Default
  #26
New Member
 
Jialin Su
Join Date: Mar 2013
Location: Loughborough
Posts: 29
Rep Power: 3
callumso is on a distinguished road
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 non-linearity 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 right-hand side of the UEqn minus the product of (the off-diagonal 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
callumso is offline   Reply With Quote

Old   April 8, 2013, 02:21
Default
  #27
Senior Member
 
Dr. Alexander Vakhrushev
Join Date: Mar 2009
Posts: 210
Rep Power: 9
makaveli_lcf is on a distinguished road
Send a message via ICQ to makaveli_lcf
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

Franz-Josef-Str. 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 02:23. Reason: Additions
makaveli_lcf is offline   Reply With Quote

Old   April 8, 2013, 07:37
Default
  #28
New Member
 
Jialin Su
Join Date: Mar 2013
Location: Loughborough
Posts: 29
Rep Power: 3
callumso is on a distinguished road
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
callumso is offline   Reply With Quote

Old   April 8, 2013, 07:44
Default
  #29
Senior Member
 
Dr. Alexander Vakhrushev
Join Date: Mar 2009
Posts: 210
Rep Power: 9
makaveli_lcf is on a distinguished road
Send a message via ICQ to makaveli_lcf
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.
sharonyue likes this.
__________________
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; April 8, 2013 at 07:45. Reason: Mods
makaveli_lcf is offline   Reply With Quote

Old   April 8, 2013, 10:45
Default
  #30
New Member
 
Jialin Su
Join Date: Mar 2013
Location: Loughborough
Posts: 29
Rep Power: 3
callumso is on a distinguished road
Danke sehr, Alex. Jetz verstehe ich.

Freundliche Gruesse
Jialin Su
callumso is offline   Reply With Quote

Old   April 13, 2013, 11:05
Default
  #31
Member
 
Avdeev Evgeniy
Join Date: Jan 2011
Location: Togliatty, Russia
Posts: 48
Blog Entries: 1
Rep Power: 11
j-avdeev will become famous soon enough
Send a message via Skype™ to j-avdeev
Hello,
I have one more question around terms .A(), .H(), .T().

As said before
Quote:
Originally Posted by mkraposhin View Post
.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
So what return terms .lower(), .upper(), .diag()?

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();
...
And I checked output: .A() and .diag() are not the same. What the difference?

Thanks.
j-avdeev is offline   Reply With Quote

Old   April 13, 2013, 14:44
Default
  #32
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
santiagomarquezd will become famous soon enough
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 }
Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   May 28, 2013, 05:32
Default
  #33
New Member
 
Emad Tandis
Join Date: Sep 2010
Posts: 26
Rep Power: 6
EmadTandis is on a distinguished road
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?
EmadTandis is offline   Reply With Quote

Old   July 4, 2013, 06:42
Default
  #34
Senior Member
 
Dongyue Li
Join Date: Jun 2012
Location: Torino, Italy
Posts: 549
Rep Power: 6
sharonyue is on a distinguished road
Quote:
Originally Posted by j-avdeev View Post
Hello,
I have one more question around terms .A(), .H(), .T().

As said before


So what return terms .lower(), .upper(), .diag()?

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();
...
And I checked output: .A() and .diag() are not the same. What the difference?

Thanks.
Hi, your problem is totally the same with mine,Ueqn.A() extract the diag. term from the Ueqn??, Do you find out why?

Thanks
sharonyue is offline   Reply With Quote

Old   July 4, 2013, 06:45
Default
  #35
Senior Member
 
Dongyue Li
Join Date: Jun 2012
Location: Torino, Italy
Posts: 549
Rep Power: 6
sharonyue is on a distinguished road
Quote:
Originally Posted by EmadTandis View Post
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?
I think only when entering another runtimeloop, UEqn.A() will change. So in PISO loop, it will not change.

If I was wrong correct me.
sharonyue is offline   Reply With Quote

Old   July 9, 2013, 10:21
Default
  #36
Member
 
Avdeev Evgeniy
Join Date: Jan 2011
Location: Togliatty, Russia
Posts: 48
Blog Entries: 1
Rep Power: 11
j-avdeev will become famous soon enough
Send a message via Skype™ to j-avdeev
Santiago, propably nooby question, anyway

from fvMatrix.C
Code:
00742     tAphi().internalField() = D()/psi_.mesh().V();
This string accord to
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.
j-avdeev 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 09:23
Programming in OpenFOAM vinu OpenFOAM 2 July 11, 2009 11:16
two questions about udf programming Vincent FLUENT 1 October 5, 2008 01: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 08:58.