Tutorial 23 – Inlet Velocity Profiles in Channel Flows

Posted on July 17, 2017 - Laminar Flow, Tutorials

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

Introduction

In this tutorial we will prescribe nonhomogeneous velocity profile to inlet boundary. It may become very handy as we want to run advanced and more accurate simulations. It can be also done in order to improve convergance of one’s simulation.

We will focus here on two inhomogeneous velocity profiles. First profile will be prescribed in dependence of radius of a pipe, while second one will be prescribed in dependence of two axis.

Both cases are ilustrated below

Cylindrical pipe case:

For the circular pipe, inlet velocity is of course a function of the pipe radius,

    $$U\left(r\right)=U_{0}\left(1-\left(\frac{r}{R}\right)^{2}\right).$$

Here, $U_{0}$ is the maximal (i.e. centerline) velocity and $R$ is the pipe radius. Rectangular pipe case:

The velocity profile in the rectangular pipe will be prescribed in dependence of width and height of a pipe,

    $$U(x,y)=\frac{16U_{0}}{H_{x}^{2}\cdot H_{y}^{2}}\cdot(\frac{H_{x}}{2}-x)\cdot(\frac{H_{x}}{2}+x)\cdot(\frac{H_{y}}{2}-y)\cdot(\frac{H_{y}}{2}+y).$$

Here, U_{0} is the maximal (i.e. centerline) velocity, H_{x} is pipe height and H_{y} is pipe width. See Figure [fig:rec-pipe] for how it looks. Also note that this profile is not an accurate solution of the Navier-Stokes equation in a rectangular channel. You will see that in simulation results. In contrast, the flow in a pipe is an exact solution.

Code

Clear MATLAB Command Window and Workspace

clc;
clear;


% Read and display mesh
% Cylidrical pipe

% [p,e,t] = importMeshGmsh('cylinder_pipe.msh');
% displayMesh(p,t);


% Rectangular pipe

[p,e,t] = importMeshGmsh('rectangular_pipe.msh');
displayMesh(p,t);


% Convert mesh to second order

[p,e,t,nVnodes,nPnodes,indices] = convertMeshToSecondOrder(p,e,t);



% Define fluid

nu = 1.5e-5;


% Init solution and convergence criteria

[u, convergence] = initSolution(p,t,[0 0 0],0);
maxres = 1e-10;
maxiter = 20;



% Define inlet velocity profile
% Paraboloidal velocity profile
% First of all declare radius of cylindrical pipe and maximum velocity.
% We will also defien width and height of rectangular pipe.

R=0.002;
vmax=0.4;
h_x=0.003;
h_y=0.002;

Secondly, we have to specify velocity vectors in all mesh elements. We will only use vectors that are attached to mesh boundary elemnts, but in order to do so we must prescribe velocity field in the whole domain. To prescribe paraboloidal velocity field we have to build vectors that will consist of two zero values in X and Y direction and value specified by location of mesh element i Z direction

% Cylindric pipe inlet velocity profile
% vx= zeros(1,nVnodes);
% vy= zeros(1,nVnodes);
% vz= vmax - (sqrt(p(1,:).^2+p(2,:).^2)).^2 *vmax/((R)^2);
% vel_cylindrical=[vx; vy; vz]';

% Square pipe welocity profile
vx= zeros(1,nVnodes);
vy= zeros(1,nVnodes);
vz= vmax*16/(h_x^2*h_y^2)*(-p(1,:)+h_x/2).*(p(1,:)+h_x/2).*(-p(2,:)+h_y/2).*...
    (p(2,:)+h_y/2);
vel_rectangular=[vx; vy; vz]';

% Iterate nonlinear terms

for iter = 1:maxiter

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


    % Apply boundary conditions

    % Now we assign velocity vectors to inlet

    %[NS, F] = imposeCfdBoundaryCondition(p,e,t,NS,F,28,'inlet',vel_cylindrical);
    [NS, F] = imposeCfdBoundaryCondition(p,e,t,NS,F,28,'inlet',vel_rectangular);
    [NS, F] = imposeCfdBoundaryCondition(p,e,t,NS,F,29,'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


% Plot solution

% Rectangular pipe

[us_c,ps_c,es_c,ts_c]=extractDataOnSurface(p,e,t,u(indices.indw),[0 0 0],[1 0 0]);
figure(3)
displaySolution(ps_c,ts_c,us_c,'Velocity on surface x = 0')


% % Cylindrical pipe
% [us_r,ps_r,es_r,ts_r]=extractDataOnSurface(p,e,t,u(indices.indw),[0 0 0],[1 0 0]);
% figure(3)
% displaySolution(ps_r,ts_r,us_r,'Velocity on surface x = 0')

figure(4)
displaySolution(p,t,u(indices.indw),'Velocity in Z direction')


Tags