1 Darcy Fochheimer Explanation

The Darcy Forchheimer model allows us to simply add a porosity zone into our fluid domain without any expense. In order to use the model, you have to put a fvOptions file into the constant folder including the following content (OpenFOAM-v6):

/*--------------------------------*- C++ -*----------------------------------*\ 
| =========                 |                                                 | 
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           | 
|  \\    /   O peration     | Version:  dev                                   | 
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      | 
|    \\/     M anipulation  |                                                 | 
\*---------------------------------------------------------------------------*/ 
FoamFile                                                                        
{                                                                               
    version     2.0;                                                            
    format      ascii;                                                          
    class       dictionary;                                                     
    location    "constant";                                                     
    object      fvOptions;                                                      
}                                                                               
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 
                                                                                
porosity1                                                                       
{                                                                               
    type            explicitPorositySource;                                     
                                                                                
    explicitPorositySourceCoeffs                                                
    {                                                                           
        selectionMode   cellZone;                                               
        cellZone        cat1;                                                   
                                                                                
        type            DarcyForchheimer;                                       
                                                                                
        f   (0.63 1e6 1e6);                                                     
        d   (80.25 1e6 1e6);                                                    
                                                                                
        coordinateSystem                                                        
        {                                                                       
            type    cartesian;                                                  
            origin  (0 0 0);                                                    
            coordinateRotation                                                  
            {                                                                   
                type    axesRotation;                                           
                e1  (1 0 0);                                                    
                e2  (0 1 0);                                                    
            }                                                                   
        }                                                                       
    }                                                                           
}

2 The Darcy-Forchheimer Equation

The Darcy Forchheimer acts in the momentum equation as a sink term $S_m$ . Considering the momentum equation, it follows:

$\frac{\partial \rho \textbf{U}}{\partial t} + \nabla (\rho \textbf{U}\textbf{U}) = \nabla \boldsymbol \sigma + S_m$

Here, the Cauchy stress tensor $\boldsymbol \sigma$ is not split into its deviatoric and hydrostatic part (shear-rate and pressure). The main important term is the source term $S_m$ which is given as:

$S_m = \left(\mu \textbf{D} + \frac{1}{2}\rho\operatorname{tr}(\textbf{U}\bullet\textbf{I})\textbf{F}\right)\textbf{U}$

While the coefficients $\textbf{D}$ and $\textbf{F}$ have to be specified in the fvOptions file (see code above). The souce term $S_m$ represents a sink as the sign is negative.

2.1 Calculate the Coefficients

In order to use the Darcy-Forchheimer equation, the coefficients D and F have to be estimated. For porous media we need the velocity dependent pressure p(u) function of the porous media. E.g., the above equation for $S_m$ can be written as:

$S_m = \mu \textbf{D}\textbf{U} + \frac{1}{2}\rho \operatorname{tr}(\textbf{U}\bullet\textbf{I})\textbf{F}\textbf{U}$

As the source term $S_m$ can be expressed as a pressure gradient it follows $S_m = \nabla p$ . Therefore, we can write:

$\nabla p = \mu \textbf{D}\textbf{U} + \frac{1}{2}\rho \operatorname{tr}(\textbf{U}\bullet\textbf{I})\textbf{F}\textbf{U}$

While switching to the Cartesian coordinate system, we can achieve something similar to:

$\nabla p = \mu \textbf{D}_{i}u_i + \frac{1}{2}\rho \textbf{F}_{i} |u_{kk}| u_i$

If we know the velocity dependent pressure we can calculate a polynomial function such as:

$\nabla p = \frac{\Delta p}{\Delta x} = A u + B u^2$

Comparing both equations, we can say:

$A = \mu \textbf{D}$

$B = \frac{1}{2}\rho\textbf{F}$

Therefore, knowing the velocity dependent pressure values, we can calculate A and B and thus, the coefficients $\textbf{D}$ and $\textbf{F}$ can be recalculated. It should be clear, that the vectors $\textbf{D}$ and $\textbf{F}$ can be direction dependent.

For an easy and fast calculation, one can use the online calculator tool provided by Tobias Holzmann:

https://holzmann-cfd.com/community/blog-and-tools/darcy-forchheimer

2.2 Special treatment in OpenFOAM

If we specify a negative number to a direction of the Darcy and/or Forchheimer coefficient in the fvOptions file, this value gets multiplied by the largest value of the tensor and -1. Hence, the negative number gets a large positive one. Hence, it is easier to distinguish if a direction is valid or not and we still can have the blockage. The code snippet for this treatment is given in the porisityModel.C file

void Foam::porosityModel::adjustNegativeResistance(dimensionedVector& resist)
{
    scalar maxCmpt = cmptMax(resist.value());

    if (maxCmpt < 0)
    {
        FatalErrorInFunction
            << "Cannot have all resistances set to negative, resistance = "
            << resist
            << exit(FatalError);
    }
    else
    {
        vector& val = resist.value();
        for (label cmpt = 0; cmpt < vector::nComponents; cmpt++)
        {
            if (val[cmpt] < 0)
            {
                val[cmpt] *= -maxCmpt;
            }
        }
    }
}

2.3 Porous Media such as Honeycombs

Porous media that only have one flow direction, e.g., honeycombs, has to block two directions of the flow. To achieve that, we can set high values for $\textbf{D}$ and $\textbf{F}$ .

DarcyForchheimer

Contents

1 Darcy Fochheimer Explanation

2 The Darcy-Forchheimer Equation

2.1 Calculate the Coefficients

2.2 Special treatment in OpenFOAM

2.3 Porous Media such as Honeycombs

3 References