# Tutorial 20 – Non-Newtonian Flow Simulation

## Code

```% Clear MATLAB Command Window and Workspace
clc;
clear;

[p,e,t] = importMeshGmsh('uBend.msh');
displayMesh2D(p,t);
``` ```% Convert mesh to second order
[p,e,t,nVnodes,nPnodes,indices] = convertMeshToSecondOrder(p,e,t);

% Define fluid
nu0 = 10;
nuinf = 0;
n = 0.8;
lambda = 2;

% Init solution and convergence criteria
[u, convergence] = initSolution(p,t,[1 0],0);
maxres = 1e-3;
maxiter = 25;

% Define inlet velocity profile
vel = [1 0];

% Iterate nonlinear terms
for iter = 1:maxiter
% Compute shear rate and viscosity for non-Newtonian fluids if needed
gamma = shearRate(p,t,u);
nu = nuinf + (nu0-nuinf)*(1+(lambda*gamma).^2).^((n-1)/2);

% Assemble matrix and right hand side
[NS, F] = assembleNavierStokesMatrix2D(p,e,t,nu,u(indices.indu),u(indices.indv),'nosupg');
% [NS, F] = assembleNavierStokesMatrix2D(p,e,t,nu,u(indices.indu),u(indices.indv),'supgDoublyAsymptotic');

% Apply boundary conditions
[NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,18,'inlet', [1 0]);
[NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,19,'wall');

% Compute and plot residuals
[stop, convergence] = computeResiduals(NS,F,u,size(p),convergence,maxres);
plotResiduals(convergence,2);

% Break if solution converged
if(stop)
break;
end

% Solve equations
u = NS\F;
end
% Display solution
% We limit range of displayed shear rate in order to reduce impact of
% artefacts that are caused by homogenous boundary condition on inlet. To
% reduce it you can prescribe parabolic velocity profile on inlet.
figure(3);
displaySolution2D(p,t,gamma,'Shear rate [1/s]', [0 1200]);
``` 