Question about the Porosity model using DarcyForchheimer Equation

Chaoran · September 12, 2020, 01:21

I'm confused about the implementation of the DarcyForchheimer equation in Openfoam. As I found in openfoam wiki https://openfoamwiki.net/index.php/DarcyForchheimer
the source term was like
$S_m = - (\mu \bf{D}+\frac{1}{2}\rho tr(\bf{U}\cdot\bf{I})\bf{F})\bf{U}$
or
$\Delta p = -(\mu \bf{D}_i \bf{u}_i - \frac{1}{2} \bf{F}_i |u_{kk}|u_i)$
But when I dig into the source code in DarcyForchheimer.H. I found something like

Code:

 forAll(cells, i)
         {
             const label celli = cells[i];
             const label j = this->fieldIndex(i);
             const tensor Cd =
                 mu[celli]*dZones[j] + (rho[celli]*mag(U[celli]))*fZones[j];
 
             const scalar isoCd = tr(Cd);
 
             Udiag[celli] += V[celli]*isoCd;
             Usource[celli] -= V[celli]*((Cd - I*isoCd) & U[celli]);
         }

where user defined values of D and f from the fvOptions dictionary are stored in the diagonal elements of the tensor $\emph{dZones[j]}$ and 0.5*f in $\emph{fZones[j]}$ .
I just can't understand the relations between the source code and the above equation. Why use the $\emph{isoCd}$ trace of tensor $\emph{Cd}$ and the term $\emph{(Cd - I*isoCd)}$ .

And I want to implement a source term like
$S_m = -(A + B*|\bf{u}|^{0.9})\bf{u}$
Can anybody help? Thanks in advance.

Tobi · September 29, 2020, 07:33

Hi Chaoran, well this is not trivial out of the box but is related to numerical stabilization. What is the Darcy-Forchheimer actually? It is a source term in the momentum equation, right? So putting the source term directly to the source (RHS) does make the system more explicit and as you can see here, we extract some part from the source term to add it implicitly to the matrix (it is actually the trace part).

So what we see is that we do create the source Cd. Then we extract the implicit part and put that to the matrix. However, the explicit part does now have to be modified by removing the already added implicitly part which lead to the (Cd - I*isoCd) term.

However, sure, you can re-implement the following:

Code:

             const label celli = cells[i];
             const label j = this->fieldIndex(i);
             const tensor Cd =
                 mu[celli]*dZones[j] + (rho[celli]*mag(U[celli]))*fZones[j];
 
             Usource[celli] -= V[celli]*(Cd  & U[celli]);

but this is 100 % explicit and requires more outer corrections because the source is dependent on the velocity.

That´s it. Even without making the mathematicial proof. The explanation abvoe should be fine.

Chaoran · September 29, 2020, 21:17

Quote:

Originally Posted by Tobi

Hi Chaoran, well this is not trivial out of the box but is related to numerical stabilization. What is the Darcy-Forchheimer actually? It is a source term in the momentum equation, right? So putting the source term directly to the source (RHS) does make the system more explicit and as you can see here, we extract some part from the source term to add it implicitly to the matrix (it is actually the trace part).

Hi, Tobias. Your explanations help me out. Really appreciate it.

September 12, 2020, 01:21	Question about the Porosity model using DarcyForchheimer Equation	#1
Chaoran New Member Chaoran Guan Join Date: Aug 2019 Posts: 5 Rep Power: 7	I'm confused about the implementation of the DarcyForchheimer equation in Openfoam. As I found in openfoam wiki https://openfoamwiki.net/index.php/DarcyForchheimer the source term was like $S_m = - (\mu \bf{D}+\frac{1}{2}\rho tr(\bf{U}\cdot\bf{I})\bf{F})\bf{U}$ or $\Delta p = -(\mu \bf{D}_i \bf{u}_i - \frac{1}{2} \bf{F}_i \|u_{kk}\|u_i)$ But when I dig into the source code in DarcyForchheimer.H. I found something like Code: forAll(cells, i) { const label celli = cells[i]; const label j = this->fieldIndex(i); const tensor Cd = mu[celli]dZones[j] + (rho[celli]mag(U[celli]))fZones[j]; const scalar isoCd = tr(Cd); Udiag[celli] += V[celli]isoCd; Usource[celli] -= V[celli]((Cd - IisoCd) & U[celli]); } where user defined values of D and f from the fvOptions dictionary are stored in the diagonal elements of the tensor $\emph{dZones[j]}$ and 0.5f in $\emph{fZones[j]}$ . I just can't understand the relations between the source code and the above equation. Why use the $\emph{isoCd}$ trace of tensor $\emph{Cd}$ and the term $\emph{(Cd - IisoCd)}$ . And I want to implement a source term like $S_m = -(A + B\|\bf{u}\|^{0.9})\bf{u}$ Can anybody help? Thanks in advance. vitto likes this. Last edited by Chaoran; September 12, 2020 at 01:59. Reason: typing mistake*

September 29, 2020, 07:33		#2
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Bad Wörishofen Posts: 2,711 Blog Entries: 6 Rep Power: 52	Hi Chaoran, well this is not trivial out of the box but is related to numerical stabilization. What is the Darcy-Forchheimer actually? It is a source term in the momentum equation, right? So putting the source term directly to the source (RHS) does make the system more explicit and as you can see here, we extract some part from the source term to add it implicitly to the matrix (it is actually the trace part). So what we see is that we do create the source Cd. Then we extract the implicit part and put that to the matrix. However, the explicit part does now have to be modified by removing the already added implicitly part which lead to the (Cd - IisoCd) term. However, sure, you can re-implement the following: Code: const label celli = cells[i]; const label j = this->fieldIndex(i); const tensor Cd = mu[celli]dZones[j] + (rho[celli]mag(U[celli]))fZones[j]; Usource[celli] -= V[celli]*(Cd & U[celli]); but this is 100 % explicit and requires more outer corrections because the source is dependent on the velocity. That´s it. Even without making the mathematicial proof. The explanation abvoe should be fine. Origami and wudong like this. __________________ Keep foaming, Tobias Holzmann

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Table bounds warnings at: END OF TIME STEP	CFXer	CFX	4	July 17, 2020 00:44
Porosity and the energy equation	shock77	OpenFOAM Running, Solving & CFD	11	May 9, 2020 16:52
[swak4Foam] swakExpression not writing to log	alexfells	OpenFOAM Community Contributions	3	March 16, 2020 19:19
Identification of the porous region in powerLawLopesdaCosta porosity model	Tarak	OpenFOAM	1	February 11, 2020 17:32
Two-Phase Buoyant Flow Issue	Miguel Baritto	CFX	4	August 31, 2006 13:02