Contrib icoBlockedCellFoam

Valid versions:

1 Introduction

The current example is desired to show how the matrix coefficients of the linear system of equations are modified from user-side. The specific aim of the presented test case is to exclude single cell from the flow and pressure field calculations. The application of such approach can be found e.g. when modeling solidification processes where fully solidified regions should not produce neither flow nor pressure sources. Also it can be used for simulating obstacles immersed in the flow etc.

Disclaimer: this example is not claimed to be perfect and also physics behind can be discussed. So please consider it as a guide for the future steps.

To exclude one cell of the domain from calculations a following technique will be used: we set diagonal coefficient to the identity value and remove the influence of the neighboring cells. Latest is done by canceling the off-diagonal coefficients. Thereby matrix coefficients are adjusted in such manner for the system

$\mathbf{A}_{ii}\boldsymbol{\psi}_{i}+\sum_{j}\mathbf{A}_{ij}\boldsymbol{\psi}_{j}=\mathbf{S}_{i}\qquad(1)$

that for the "blocked" cell (indexed as k) the diagonal element $\mathbf{A}_{kk}$ equals 1, and all off-diagonal coefficients $\mathbf{A}_{kj}$ are zeroed along with the source term $\mathbf{S}_{k}$. The consequence of the described operations with the system (1) is that

$\left\{ \begin{array}{c} \mathbf{A}_{kk}=1,\\ \mathbf{A}_{kj}=0,\\ \mathbf{S}_{k}=0.\end{array}\quad\Longrightarrow\quad\mathbf{1}\cdot\boldsymbol{\psi}_{k}=0.\right. \qquad(2)$

In accordance to (2) the solution for the cell value $\boldsymbol{\psi}_{k}$ is trivial one. Hereby if one applies the condition (2) for the momentum and pressure matrices in the PISO algorithm, that results in zero velocity and pressure values automatically for the general solution of the full system (1). To get the solution in the cell different from the zero one, the source term can be adjusted to the corresponding value.

The source code example supplemented with the test case is a modification of the incompressible viscous flow solver icoFoam Further the required additions are described.

2 Mods to createFields.C file

The following modification should be applied to the initial createFields.C file. First for the allocation of the blocked cell it is required to define its position:

 
//  Read the position of the excluded cell
vector position_
(
    transportProperties.lookup("exclude")
);

Then exclude parameter should be defined in the transportProperties file in the case folder.

nu          nu [ 0 2 -1 0 0 0 0 ] 0.01;
exclude     (0.05 0.05 0.005);

In the above example the position vector will correspond to the center of the driven cavity.

3 Solver icoBlockedCell.C file

The main icoFoam.C file of the icoFoam solver is supplemented then with the matrices manipulation commands. First of all the cell containing the position defined by the exclude parameter in transportProperties dictionary is searched for through the whole mesh object:

 
//  Find the cell at the defined location to exclude from calculation
label cell_ = mesh.findCell(position_);
Info<< "Location " << position_ << " is at the cell #"
    <<  cell_ << "." << nl << endl;

Then the diagonal element of the momentum equation is set to the identity value and all neighboring cells coefficients are sequentially set to zero. As the final step here the source term is also eliminated and modified system for the momentum predictor is solved.

 
//  Modify matrix coefficients for momentum equation
{
    //  Access lower, upper triangles and diagonal to modify
    scalarField& lower_ = UEqn.lower();
    scalarField& upper_ = UEqn.upper();
    scalarField& diag_ = UEqn.diag();
 
    //  Get matrix LDU addressing.
    //  Connects face indexing for the off-diagonal elements
    //  with cell indexing numbering of the diagonal
    const unallocLabelList& l_ = UEqn.lduAddr().lowerAddr();
    const unallocLabelList& u_ = UEqn.lduAddr().upperAddr();
 
    //  Set diagonal coefficient to identity value
    diag_[cell_] = scalar(1);
 
    //  Loop over all off-diagonal elements.
    //  Set them to zero value to exclude the influence
    //  of the neighbouring cells
    for (register label face_=0; face_<l_.size(); face_++)
    {
    	if (l_[face_] == cell_)
    	{
        	lower_[face_] = scalar(0);
    	}
 
    	if (u_[face_] == cell_)
    	{
        	upper_[face_] = scalar(0);
    	}
    }
 
    //  Exclude source term, which comes from dU/dt discretization
    UEqn.source()[cell_] = vector(0, 0, 0);
//    Info<< "UEqn sorce term:" << nl
//    	<< UEqn.source() << nl << endl;
 
    //  Copy matrix to manipulate the pressure gradient source term.
    //  Initial UEqn is required for pressure correction.
    fvVectorMatrix UpGradEqn = UEqn;
 
    //  Include pressure for all cells
    UpGradEqn += fvc::grad(p);
 
    //  Exclude any sources for the cell in desired location
    UpGradEqn.source()[cell_] = vector(0, 0, 0);
 
    //  Solve corrected momentum predictor equation
    UpGradEqn.solve();
}

Corresponding step is repeated for the pressure corrector step:

 
//  Modify pressure equation matrix
{
    //  Take care about matrix symmetry.
    //  Due to read-write access to the off-diagonal elements
    //  it can be converted to assymetric one!
    //  Thus either lower or upper elements are used.
    scalarField& off_diag_ =
        (pEqn.hasUpper()) 
            ? pEqn.upper() 
            : pEqn.lower();
 
    //  Get diagonal coefficients
    scalarField& diag_ = pEqn.diag();
 
    //  LDU off-diagonal addressing for symmetric matrix
    const unallocLabelList& o_ =
        (pEqn.hasUpper())
            ? pEqn.lduAddr().upperAddr()
            : pEqn.lduAddr().lowerAddr();
 
    //  Set diagonal to identity
    diag_[cell_] = scalar(1);
 
    //  Remove the influence of the neighbours 
    for (register label face_=0; face_<o_.size(); face_++)
    {
        if (o_[face_] == cell_)
        {
            off_diag_[face_] = scalar(0);
        }
    }
 
    //  Remove source for desired cell
    pEqn.source()[cell_] = scalar(0);
}

4 Results

Here what we get if the central cell is "blocked'" for the driven cavity case (Fig. 1)

Figure 1: Cavity case results: central cell is "blocked"

5 Outlook

Presented algorithm can be improved to take into account arbitrary number of cells.

6 Files

• IcoBlockedCellFoam.tar.gz‎: solver sources (use Allwmake script to build)
• Cavity.tar.gz: example case (supplemented with Run.sh script)

7 History

• 19-04-2011: Initial upload

--Makaveli lcf 07:38, 20 April 2011 (CEST)