[gaza]

Hi Foamers,
According to Prof. Hrvoje Jasak's PhD thesis:

"The H(U) term consists of two parts: the “transport part”, including the matrix coefficients for all neighbours multiplied by corresponding velocities and the “source part” including the source part of the transient term and all other source terms apart from the pressure gradient".

My question is on implementation of .H() operator:
If I change the basic momentum equation and add my source term
and then I use the .H() function in the code, is it contain my source term also or I have to implement this by hand?

[gaza]

more clearly:

If I add source term to the momentum equation, is it already included in the .H() function?

[Zeppo]

Yes, when you add your (user-difined) source to the equation, you are actually modifying the source coefficients (source_) of fvMatrix object to be solved. All the modifications are incorporated into the matrix. The function fvMatrix<Type>::H() const does the following:

// "src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C"
...
788     Hphi.primitiveFieldRef() += lduMatrix::H(psi_.primitiveField()) + source_;
...

[gaza]

Hi Zeppo,
Thank you very much for your help

[mszeto715]

Originally Posted by Zeppo

Yes, when you add your (user-difined) source to the equation, you are actually modifying the source coefficients (source_) of fvMatrix object to be solved. All the modifications are incorporated into the matrix. The function fvMatrix<Type>::H() const does the following:

// "src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C"
...
788     Hphi.primitiveFieldRef() += lduMatrix::H(psi_.primitiveField()) + source_;
...

Hi,

Is lduMatrix::H(psi_.primitiveField()) negative?

There is a command that "negates" a couple lines above that. I just want to confirm that that makes the code consistent with other explanations, e.g., https://openfoamwiki.net/index.php/O...ide/H_operator

Thanks.
-Mimi

[mAlletto]

Originally Posted by gaza

Hi Foamers,
According to Prof. Hrvoje Jasak's PhD thesis:

"The H(U) term consists of two parts: the “transport part”, including the matrix coefficients for all neighbours multiplied by corresponding velocities and the “source part” including the source part of the transient term and all other source terms apart from the pressure gradient".

My question is on implementation of .H() operator:
If I change the basic momentum equation and add my source term
and then I use the .H() function in the code, is it contain my source term also or I have to implement this by hand?

I just had a look at the source code in fvMatrix.C:

1523 template<class Type>
1524 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1525 (
1526  const fvMatrix<Type>& A,
1527  const DimensionedField<Type, volMesh>& su
1528 )
1529 {
1530  checkMethod(A, su, "==");
1531  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1532  tC.ref().source() += su.mesh().V()*su.field();
1533  return tC;
1534 }
1535 
1536 template<class Type>
1537 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1538 (
1539  const fvMatrix<Type>& A,
1540  const tmp<DimensionedField<Type, volMesh>>& tsu
1541 )
1542 {
1543  checkMethod(A, tsu(), "==");
1544  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1545  tC.ref().source() += tsu().mesh().V()*tsu().field();
1546  tsu.clear();
1547  return tC;
1548 }
1549 
1550 template<class Type>
1551 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1552 (
1553  const fvMatrix<Type>& A,
1554  const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu
1555 )
1556 {
1557  checkMethod(A, tsu(), "==");
1558  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1559  tC.ref().source() += tsu().mesh().V()*tsu().primitiveField();
1560  tsu.clear();
1561  return tC;
1562 }
1563 
1564 template<class Type>
1565 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1566 (
1567  const tmp<fvMatrix<Type>>& tA,
1568  const DimensionedField<Type, volMesh>& su
1569 )
1570 {
1571  checkMethod(tA(), su, "==");
1572  tmp<fvMatrix<Type>> tC(tA.ptr());
1573  tC.ref().source() += su.mesh().V()*su.field();
1574  return tC;
1575 }
1576 
1577 template<class Type>
1578 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1579 (
1580  const tmp<fvMatrix<Type>>& tA,
1581  const tmp<DimensionedField<Type, volMesh>>& tsu
1582 )
1583 {
1584  checkMethod(tA(), tsu(), "==");
1585  tmp<fvMatrix<Type>> tC(tA.ptr());
1586  tC.ref().source() += tsu().mesh().V()*tsu().field();
1587  tsu.clear();
1588  return tC;
1589 }
1590 
1591 template<class Type>
1592 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1593 (
1594  const tmp<fvMatrix<Type>>& tA,
1595  const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu
1596 )
1597 {
1598  checkMethod(tA(), tsu(), "==");
1599  tmp<fvMatrix<Type>> tC(tA.ptr());
1600  tC.ref().source() += tsu().mesh().V()*tsu().primitiveField();
1601  tsu.clear();
1602  return tC;
1603 }
1604 
1605 template<class Type>
1606 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1607 (
1608  const fvMatrix<Type>& A,
1609  const dimensioned<Type>& su
1610 )
1611 {
1612  checkMethod(A, su, "==");
1613  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1614  tC.ref().source() += A.psi().mesh().V()*su.value();
1615  return tC;
1616 }
1617 
1618 template<class Type>
1619 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1620 (
1621  const tmp<fvMatrix<Type>>& tA,
1622  const dimensioned<Type>& su
1623 )
1624 {
1625  checkMethod(tA(), su, "==");
1626  tmp<fvMatrix<Type>> tC(tA.ptr());
1627  tC.ref().source() += tC().psi().mesh().V()*su.value();
1628  return tC;
1629 }

From the source code I would judge that everything on the right hand side of the operator == of the equation e.g. the Ueqn will be put to the source_ attribute of the fvMatrix. Or better said the operator takes the left side (a fvMatrix object) adds the right hand side of operator (a fvField) to the source_ attribute and returns the reference to the fvMatrix:

1  // Construct the Momentum equation
2 
3  MRF.correctBoundaryVelocity(U);
4 
5  tmp<fvVectorMatrix> tUEqn
6  (
7  fvm::div(phi, U)
8  + MRF.DDt(U)
9  + turbulence->divDevReff(U)
10  ==
11  fvOptions(U)
12  );
13  fvVectorMatrix& UEqn = tUEqn.ref();
14 
15  UEqn.relax();
16 
17  // Include the porous media resistance and solve the momentum equation
18  // either implicit in the tensorial resistance or transport using by
19  // including the spherical part of the resistance in the momentum diagonal
20 
21  tmp<volScalarField> trAU;
22  tmp<volTensorField> trTU;
23 
24  if (pressureImplicitPorosity)
25  {
26  tmp<volTensorField> tTU = tensor(I)*UEqn.A();
27  pZones.addResistance(UEqn, tTU.ref());
28  trTU = inv(tTU());
29  trTU.ref().rename("rAU");
30 
31  fvOptions.constrain(UEqn);
32 
33  volVectorField gradp(fvc::grad(p));
34 
35  for (int UCorr=0; UCorr<nUCorr; UCorr++)
36  {
37  U = trTU() & (UEqn.H() - gradp);
38  }
39  U.correctBoundaryConditions();
40 
41  fvOptions.correct(U);
42  }
43  else
44  {
45  pZones.addResistance(UEqn);
46 
47  fvOptions.constrain(UEqn);
48 
49  solve(UEqn == -fvc::grad(p));
50 
51  fvOptions.correct(U);
52 
53  trAU = 1.0/UEqn.A();
54  trAU.ref().rename("rAU");
55  }

So I would guess that the pressure gradient is also included in the H() function. Do you agree ore am I completely wrong.