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

## Loading and preparing the mesh

First, we load and display the mesh.

**foam.mat 2.63 MB**

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

**Tutorial 29 _ results.mat 5.54 MB**

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