Description of flux() method

fumiya · March 20, 2013, 05:19

I try to describe the role of flux() method and how it achieves it's role in OpenFOAM.
This topic is difficult for me to comprehend, so I'm grad if the following explanations could be of some help

Role of flux() method
To construct the conservative face flux
How to achieve this
If we discretize the pressure poisson equation(e.g., pEqn in simpleFoam), we can get

$\sum_{f} \left \{\left ( \frac{1}{A_{P}} \right )_{f}\left | \mathbf{\Delta } \right |\frac{p_{N}-p_{P}}{\left | \mathbf{d} \right |} +\left ( \frac{1}{A_{P}} \right )_{f}\mathbf{k}\cdot \left ( \triangledown p \right )_{f} \right \}=\sum_{f}\left ( \frac{H}{A_{P}} \right )_{f}\cdot \mathbf{S}_{f}$ ---(1)

where $\mathbf{S}_{f} = \mathbf{\Delta } + \mathbf{k}$ (See the attached picture).

The second term of the l.h.s of the eqn. (1) is the explicit non-orthogonal correction term
(source term) and the first term is the implicit orthogonal part(matrix coefficients).
If we arrange the above equation using OpenFOAM notation(v2.2), we get

$\sum_{f} \left \{ \mathrm{phiHbyA} -\left ( \left ( \frac{1}{A_{P}} \right )_{f}\left | \mathbf{\Delta } \right |\frac{p_{N}-p_{P}}{\left | \mathbf{d} \right |} +\left ( \frac{1}{A_{P}} \right )_{f}\mathbf{k}\cdot \left ( \triangledown p \right )_{f} \right ) \right \}=0$ ---(2)

Then we can construct the conservative face flux phi by

$\mathrm{phi} = \mathrm{phiHbyA} - \left ( \left ( \frac{1}{A_{P}} \right )_{f}\left | \mathbf{\Delta } \right |\frac{p_{N}-p_{P}}{\left | \mathbf{d} \right |} +\left ( \frac{1}{A_{P}} \right )_{f}\mathbf{k}\cdot \left ( \triangledown p \right )_{f} \right )$ ---(3)

In the eqn. (3), the pEqn.flux() calculates

$\left ( \frac{1}{A_{P}} \right )_{f}\left | \mathbf{\Delta } \right |\frac{p_{N}-p_{P}}{\left | \mathbf{d} \right |} +\left ( \frac{1}{A_{P}} \right )_{f}\mathbf{k}\cdot \left ( \triangledown p \right )_{f}$ ---(4)

,so we can get the conservative face flux using the following code
Code:
```
if (simple.finalNonOrthogonalIter())
{
    phi = phiHbyA - pEqn.flux();
}
```

Where is flux() function implemented?
As you can see in the eqn. (4), the flux() is the sum of

1st term
the off-diagonal coefficients of the pressure poisson equation multiplied by pressure values at cell centers

in fvMatrix.C

Code:

00895     for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
00896     {
00897         fieldFlux.internalField().replace
00898         (
00899             cmpt,
00900             lduMatrix::faceH(psi_.internalField().component(cmpt))
00901         );
00902     }

and lduMatrix::faceH is implemented in lduMatrixTemplates.C

Code:

00092         for (register label face=0; face<l.size(); face++)
00093         {
00094             faceHpsi[face] =
00095                 Upper[face]*psi[u[face]]
00096               - Lower[face]*psi[l[face]];
00097         }

2nd term
the non orthogonality explicit corrections

in fvMatrix.C

Code:

00937     if (faceFluxCorrectionPtr_)
00938     {
00939         fieldFlux += *faceFluxCorrectionPtr_;
00940     }

Please correct the mistakes, if any.

Best regards,
Fumiya

sharonyue · August 26, 2013, 00:40

Quote:

Originally Posted by fumiya

I try to describe the role of flux() method and how it achieves it's role in OpenFOAM.
This topic is difficult for me to comprehend, so I'm grad if the following explanations could be of some help

..........................

Best regards,
Fumiya

Excellent!!Fumiya now Its a bit clear regarding the p.Eqn.flux().

BTW, do you think

Code:

fvScalarMatrix pEqn
            (
                fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
            );

fvScalarMatrix pEqn
            (
                fvm::div(A*fvc::grad(p)) == fvc::div(phiHbyA)
            );

is the same?

longyun000 · December 21, 2016, 11:39

Quote:

Originally Posted by fumiya

I try to describe the role of flux() method and how it achieves it's role in OpenFOAM.
This topic is difficult for me to comprehend, so I'm grad if the following explanations could be of some help

Role of flux() method
To construct the conservative face flux
How to achieve this
If we discretize the pressure poisson equation(e.g., pEqn in simpleFoam), we can get

$\sum_{f} \left \{\left ( \frac{1}{A_{P}} \right )_{f}\left | \mathbf{\Delta } \right |\frac{p_{N}-p_{P}}{\left | \mathbf{d} \right |} +\left ( \frac{1}{A_{P}} \right )_{f}\mathbf{k}\cdot \left ( \triangledown p \right )_{f} \right \}=\sum_{f}\left ( \frac{H}{A_{P}} \right )_{f}\cdot \mathbf{S}_{f}$ ---(1)

where $\mathbf{S}_{f} = \mathbf{\Delta } + \mathbf{k}$ (See the attached picture).

The second term of the l.h.s of the eqn. (1) is the explicit non-orthogonal correction term
(source term) and the first term is the implicit orthogonal part(matrix coefficients).
If we arrange the above equation using OpenFOAM notation(v2.2), we get

$\sum_{f} \left \{ \mathrm{phiHbyA} -\left ( \left ( \frac{1}{A_{P}} \right )_{f}\left | \mathbf{\Delta } \right |\frac{p_{N}-p_{P}}{\left | \mathbf{d} \right |} +\left ( \frac{1}{A_{P}} \right )_{f}\mathbf{k}\cdot \left ( \triangledown p \right )_{f} \right ) \right \}=0$ ---(2)

Then we can construct the conservative face flux phi by

$\mathrm{phi} = \mathrm{phiHbyA} - \left ( \left ( \frac{1}{A_{P}} \right )_{f}\left | \mathbf{\Delta } \right |\frac{p_{N}-p_{P}}{\left | \mathbf{d} \right |} +\left ( \frac{1}{A_{P}} \right )_{f}\mathbf{k}\cdot \left ( \triangledown p \right )_{f} \right )$ ---(3)

In the eqn. (3), the pEqn.flux() calculates

$\left ( \frac{1}{A_{P}} \right )_{f}\left | \mathbf{\Delta } \right |\frac{p_{N}-p_{P}}{\left | \mathbf{d} \right |} +\left ( \frac{1}{A_{P}} \right )_{f}\mathbf{k}\cdot \left ( \triangledown p \right )_{f}$ ---(4)

,so we can get the conservative face flux using the following code
Code:
```
if (simple.finalNonOrthogonalIter())
{
    phi = phiHbyA - pEqn.flux();
}
```

Where is flux() function implemented?
As you can see in the eqn. (4), the flux() is the sum of

1st term
the off-diagonal coefficients of the pressure poisson equation multiplied by pressure values at cell centers

in fvMatrix.C

Code:

00895     for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
00896     {
00897         fieldFlux.internalField().replace
00898         (
00899             cmpt,
00900             lduMatrix::faceH(psi_.internalField().component(cmpt))
00901         );
00902     }

and lduMatrix::faceH is implemented in lduMatrixTemplates.C

Code:

00092         for (register label face=0; face<l.size(); face++)
00093         {
00094             faceHpsi[face] =
00095                 Upper[face]*psi[u[face]]
00096               - Lower[face]*psi[l[face]];
00097         }

2nd term
the non orthogonality explicit corrections

in fvMatrix.C

Code:

00937     if (faceFluxCorrectionPtr_)
00938     {
00939         fieldFlux += *faceFluxCorrectionPtr_;
00940     }

Please correct the mistakes, if any.

Best regards,
Fumiya

hi,
I am studing icofoam now, and I found the code:
[/CODE]
77 surfaceScalarField phiHbyA
78 (
79 "phiHbyA",
80 fvc::flux(HbyA)
81 + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
82 );
[/CODE]

I can't understand why there is an item "fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)"? According to the equation above, there should be only the first item.

fumiya · December 24, 2016, 13:40

Hi,

The term "ddtCorr" is considered in the transient solvers.
The following references could be of some help:

Best regards,
Fumiya

yuan_neu · April 3, 2017, 21:36

Quote:

Originally Posted by sharonyue

Excellent!!Fumiya now Its a bit clear regarding the p.Eqn.flux().

BTW, do you think

Code:

fvScalarMatrix pEqn
            (
                fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
            );

fvScalarMatrix pEqn
            (
                fvm::div(A*fvc::grad(p)) == fvc::div(phiHbyA)
            );

is the same?

i do not think so. fvc::grad(p) returns a volscalarfield. so the second pEqn is actually based on grad(p).

bbita · January 2, 2018, 12:05

Dear Foamers,

I am trying to add a source term to alpha using mass flux. I am starting with a circle and while using .flux(), I can see my circle changes to a square. I am wondering if .flux() can be used in 2D the same as 1D.

Thanks,

December 24, 2016, 13:40		#4
fumiya Senior Member Fumiya Nozaki Join Date: Jun 2010 Location: Yokohama, Japan Posts: 266 Blog Entries: 1 Rep Power: 19	Hi, The term "ddtCorr" is considered in the transient solvers. The following references could be of some help: DdtPhiCorr http://openfoam.org/release/2-3-0/numerics/ Best regards, Fumiya __________________ [Personal] Blog: https://caefn.com/ Twitter: https://twitter.com/FumiyaNozaki SlideShare: http://www.slideshare.net/fumiyanozaki96

January 2, 2018, 12:05	Flux in 2D Model	#6
bbita Member Join Date: May 2017 Posts: 44 Rep Power: 9	Dear Foamers, I am trying to add a source term to alpha using mass flux. I am starting with a circle and while using .flux(), I can see my circle changes to a square. I am wondering if .flux() can be used in 2D the same as 1D. Thanks,

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Sign of Heat Flux at wall	Kyung	FLUENT	2	February 26, 2016 17:25
Comparison: Finite Volume Method vs. Analytic Method	m-fry	Main CFD Forum	1	April 20, 2010 15:40
About the Flux splitting method	maximus	Main CFD Forum	0	February 24, 2006 01:05
Replace periodic by inlet-outlet pair	lego	CFX	3	November 5, 2002 21:09
total mass flux correction for compressible fluid?	Francesco Di Maio	Main CFD Forum	0	August 21, 2000 05:23