Tutorial 20 – Non-Newtonian Flow Simulation

Code

% Clear MATLAB Command Window and Workspace
clc;
clear;


% Read and display mesh
[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]);