# Contrib icoBlockedCellFoam

## Contents

## 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

that for the "blocked" cell (indexed as
*k*) the diagonal element equals 1, and all off-diagonal
coefficients are zeroed along with the source term .
The consequence of the described operations with the system (1) is that

In accordance to (2) the solution for the cell value 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)