[idrama]

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/.

[ngj]

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

[RuiVO]

Thank you so much too both of you

Best regards

Rui

[mchurchf]

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

[santiagomarquezd]

http://www.cfd-online.com/Forums/ope...ipulation.html

http://www.cfd-online.com/Forums/ope...rua-u-phi.html

http://www.cfd-online.com/Forums/ope...predictor.html


Regards.

[callumso]

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

[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().

[callumso]

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

[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.

[callumso]

Danke sehr, Alex. Jetz verstehe ich.

Freundliche Gruesse
Jialin Su

[j-avdeev]

Hello,
I have one more question around terms .A(), .H(), .T().

As said before

Originally Posted by mkraposhin

.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

...
    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.

[santiagomarquezd]

Hi, A()=diag()/mesh.V()

See fvMatrix.C

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.

[EmadTandis]

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?

[sharonyue]

Originally Posted by j-avdeev

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

...
    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,http://www.cfd-online.com/Forums/ope...tml#post437736, Do you find out why?

Thanks

[sharonyue]

Originally Posted by EmadTandis

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.

[j-avdeev]

Santiago, propably nooby question, anyway

from fvMatrix.C

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.