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