# Tutorial 29 – Simulating Flow Through a 3D Porous Medium

This tutorial is intended for the full version of the toolbox.

## Introduction

This tutorial shows how to perform a simulation of fluid flow through a 3D porous medium reconstructed from CT imaging. To learn how to create a geometric model from raster images, see tutorial 28

First, we load and display the mesh.

```load('foam.mat', 'p', 'e', 't');
viewCase

``` Next, we convert the mesh to the second order.

```[p, e, t, no_vnodes, no_pnodes, indices] = convertMeshToSecondOrder(p, e, t);
```

For convenience, we define the boundary IDs as variables.

```inlet_ID = 2;
wall_ID = 1;
outlet_ID = 3;

```

Now we extract the IDs of the nodes lying on the inlet.

```inlet_nodes = extractNodeIdsOnFaces(e, inlet_ID);

```

## Generating the inlet velocity profile

We compute the distance between the nodes and the walls.

```node_distance = computeWallDistance(p, e, t, wall_ID);

```

Generating wall distance field…
gmres(75) converged at outer iteration 1 (inner iteration 56) to a solution with relative residual 9.5e-11.
Done

We want to specify the inlet velocity so that it is proportional to the square of the distance from the wall (so that it is in a sense paraboloid).

```inlet_velocity = node_distance.^2;

```

We will now scale it to achieve the desired flow rate.

```computed_flow_rate = boundaryIntegral(p, e, inlet_velocity, inlet_ID);
desired_flow_rate = 5e-7;
inlet_velocity = desired_flow_rate/computed_flow_rate*inlet_velocity;

```

Let us display the computed inlet velocity profile

```clf
displaySolution(p, t, inlet_velocity, 'Inlet velocity')

``` We set the remaining inlet velocity components to zero.

```inlet_velocity = [inlet_velocity, zeros(length(inlet_velocity), 2)];

```

## Setting up the solver

First, we initiate the convergence settings and the velocity field.

```[u, convergence] = initSolution(p, t, [0, 0, 0], 0);

```

Next, we define the fluid properties.

```rho = 1000;
nu = 1e-6;

```

We define the solver settings.

```max_residual = 1e-10;
max_no_iterations = 20;
underrelaxation_parameter = 0.8;
solverOptions.p = p;
solverOptions.e = e;
solverOptions.t = t;
solverOptions.diffusivity = nu;
solverOptions.maxres = 1e-12;
solverOptions.method = 'ns';
solver = linearSolver(solverOptions);
u_old = u; % copy of the velocity field for under-relaxation purposes

```

## Computations

We are ready to solve the problem. We iterate the nonlinear terms in a loop. Note that due to the size of the problem the simulation takes a long time, if you do not wish to wait, the results have been provided in

```for iter = 1:max_no_iterations
% Under-relax velocity field
u = (1 - underrelaxation_parameter)*u_old + underrelaxation_parameter*u;
% Assemble the problem
[NS, F] = assembleNavierStokesMatrix(p, e, t, nu, u(indices.indu),...
u(indices.indv), u(indices.indw), 'supgDoublyAsymptotic');
% Impose boundary coditions
[NS, F] = imposeCfdBoundaryCondition(p, e, t, NS, F, inlet_ID, 'inlet', inlet_velocity);
[NS, F] = imposeCfdBoundaryCondition(p, e, t, NS, F, wall_ID, 'wall');
% Compute and plot residuals
[stop, convergence] = computeResiduals(NS, F, u, size(p), convergence, max_residual);
plotResiduals(convergence, 2);
% Break if solution is converged
if(stop)
break;
end

% Solve eqations
u_old = u;
u = solver.Solve(NS, F, u);
end

```

0 4.2908 0 0 0
gmres(75) converged at outer iteration 2 (inner iteration 65) to a solution with relative residual 9.3e-13.
1.0000 0.8582 0.0000 0.0000 0.0000
gmres(75) stopped at outer iteration 2 (inner iteration 75) without converging to the desired tolerance 1e-12
because the maximum number of iterations was reached.
The iterate returned (number 2(75)) has relative residual 4.4e-12.
2.0000 0.1716 0.0000 0.0000 0.0000
gmres(75) stopped at outer iteration 2 (inner iteration 75) without converging to the desired tolerance 1e-12
because the maximum number of iterations was reached.
The iterate returned (number 2(75)) has relative residual 7.9e-12.
3.0000 0.0343 0.0000 0.0000 0.0000
gmres(75) stopped at outer iteration 2 (inner iteration 75) without converging to the desired tolerance 1e-12
because the maximum number of iterations was reached.
The iterate returned (number 2(75)) has relative residual 3.8e-11.
4.0000 0.0069 0.0000 0.0000 0.0000

because the maximum number of iterations was reached.
The iterate returned (number 2(75)) has relative residual 2.9e-11.
5.0000 0.0014 0.0000 0.0000 0.0000
gmres(75) stopped at outer iteration 2 (inner iteration 75) without converging to the desired tolerance 1e-12
because the maximum number of iterations was reached.
The iterate returned (number 2(75)) has relative residual 1.5e-11.
6.0000 0.0003 0.0000 0.0000 0.0000
gmres(75) stopped at outer iteration 2 (inner iteration 75) without converging to the desired tolerance 1e-12
because the maximum number of iterations was reached.
The iterate returned (number 2(75)) has relative residual 6.9e-12.
7.0000 0.0001 0.0000 0.0000 0.0000
gmres(75) stopped at outer iteration 2 (inner iteration 75) without converging to the desired tolerance 1e-12
because the maximum number of iterations was reached.
The iterate returned (number 2(75)) has relative residual 1.7e-12.
8.0000 0.0000 0.0000 0.0000 0.0000
gmres(75) stopped at outer iteration 2 (inner iteration 75) without converging to the desired tolerance 1e-12
because the maximum number of iterations was reached.
The iterate returned (number 2(75)) has relative residual 1e-12.
9.0000 0.0000 0.0000 0.0000 0.0000
gmres(75) converged at outer iteration 2 (inner iteration 69) to a solution with relative residual 1e-12.
10.0000 0.0000 0.0000 0.0000 0.0000
gmres(75) converged at outer iteration 2 (inner iteration 62) to a solution with relative residual 9.1e-13.
11.0000 0.0000 0.0000 0.0000 0.0000
gmres(75) converged at outer iteration 2 (inner iteration 60) to a solution with relative residual 9.9e-13.
12.0000 0.0000 0.0000 0.0000 0.0000
gmres(75) converged at outer iteration 2 (inner iteration 57) to a solution with relative residual 8.8e-13.
13.0000 0.0000 0.0000 0.0000 0.0000
gmres(75) converged at outer iteration 2 (inner iteration 53) to a solution with relative residual 8.9e-13.
14.0000 0.0000 0.0000 0.0000 0.0000
gmres(75) converged at outer iteration 2 (inner iteration 48) to a solution with relative residual 8.7e-13.
15.0000 0.0000 0.0000 0.0000 0.0000
gmres(75) converged at outer iteration 2 (inner iteration 43) to a solution with relative residual 8.6e-13.
16.0000 0.0000 0.0000 0.0000 0.0000 ## Postprocessing

We can now display the results. First, the x velocity field.

```displaySolution(p, t, u(indices.indu), 'x velocity')

``` Now we obtain the physical pressure by multiplying the normalized pressure from the solver by the density.

```pressure = rho*u(indices.indp);

```

Finally, we display the pressure field.

```displaySolution(p, t, pressure, 'Pressure')

``` Author: Jakub Gałecki