sandbox/easystab/boundary_layer_pyl.m

    Model for the boundary layer in curved geometry

    We seek here to build a model based on the viscous/acceleration growth/decay of a boundary layer to compare to venturi.m.

    The model used is the Interactive Boundary Layer Theory with integral method for the boundary layer. It consists to solve the conservation of flux in the ideal fluid layer: \displaystyle u_e(R-y_w-\delta_1)^2=1 It is as if the radius of the pipe, the initial one minus the shape of the stenosis: (R-y_w) was reduced because of the boundary layer thickness \delta_1. And the second equation is the growth of the boundary layer thickness because of the effect of the pipe geometry and the variation of the ideal fluid velocity u_e, this is the Von K’arm’an boundary integral layer equation: \displaystyle \frac{\partial }{\partial x} (\frac{\delta_1}{H})+\frac{\delta_1}{u_e}(1+\frac{2}{H})\frac{\partial u_e }{\partial x} =\frac{f_2H}{\delta_1 u_e}

    Note that in classical Boundary layer theory u_e is given by u_e(R-y_w)^2=1, here we have a strong coupling.

    H and f_2 are two fixed parameters (for now).

    % if octave
%setenv("GNUTERM","X11")

clear all; clf;

% parameters
L=.5;             % domain length
N=200;            % number of grid points
H=2.59;           % H shape parameter
f2=0.22;          % f2 friction parameter
d0=0.01;          % thickness at inflow
R=1;              % rayon du tuyau
xb=.2;            %position of the bump
wb=0.05 ;         % width of the bump
alpha=0.3;        % degree of closure 

% Differentiation and integration
[D,DD,wx,x]=dif1D('cheb',0,L,N,3);
Z=zeros(N,N); I=eye(N); 
    
% differentiation and integration
scale=-2/L;
[x,DM] = chebdif(N,2); 
D=DM(:,:,1)*scale;    
DD=DM(:,:,2)*scale^2;    
x=(x-1)/scale; 
Z=zeros(N,N); I=eye(N); 

% shape of the pipe is R-yw
yw=alpha*exp(-((x-xb)/wb).^2);

%initial guess
sol=[1+0*x;d0+0*1.7*x.^(.5)];
    
% Newton iterations
disp('Newton loop')
quit=0;count=0;
while ~quit     

    % the present solution and its derivatives
    u=sol(1:N); d=sol(N+1:2*N);
    ux=D*u; dx=D*d;

    The nonlinear function

    The first equation is the conservation of flux in the ideal fluid layer u_e: \displaystyle u_e(R-y_w-\delta_1)^2-1=0 This is not a differential equation, but we keep it explicitely so that the expression of the Jacobian will be easier (we could have removed u_e from the variable by saying \displaystyle u_e=1/(R-y_w-\delta_1)^2 and replaced that in the second equation)

    And the second equation is the growth of the boundary layer thickness because of the effect of the pipe geometry and the variation of the centerline velocity u_e (and as well viscosity): \displaystyle \frac{\partial }{\partial x} (\frac{\delta_1}{H})+\frac{\delta_1}{u_e}(1+\frac{2}{H})\frac{\partial u_e }{\partial x} -\frac{f_2H}{\delta_1 u_e}=0 Note that the shear stress is \tau = u_e\frac{f_2H}{\delta_1 } it is an important quantity

        % shear stress
    tau = (f2*H)*u./d ;
    
    % flux, nonlinear VK function
    f=[u.*(R-yw-d).^2-1; ...
       dx/H+(1+2/H)*d.*ux./u-tau./u./u];

    The analytical Jacobian

    We get it by replacing in the equations u_e by u_e+\tilde{u_e} and the same for the thickness \delta_1+\tilde{\delta_1}. The Jacobian is the linear operator acting on the order one terms of \tilde{u_e} and \tilde{\delta_1}.

        % analytical jacobian
    A=[diag((R-d-yw).^2),   diag(-2*u.*(R-d-yw)); ...
       (1+2/H)*(diag(1./u)*D-diag(ux./u.^2))+f2*H*diag(1./(d.*u.^2)), ...
       D/H+(1+2/H)*diag(ux./u)+f2*H*diag(1./(d.^2.*u))];

    Boundary conditions

    We have two unknowns but just one derivative, so there is one boundary condition to impose. We set at inflow (x=0) the value of the boundary layer thickness to \delta_1(x=0)=\delta_{10}.

        % Boundary conditions
    f(1)=d(1)-d0;
    A(1,:)=[Z(1,:),I(1,:)];

    % convergence test  
    plot(x,u,'b-',x,d,'r-',x,R-yw,'b-',x,tau/10,'k--'); axis([0,L,0,4*R]);grid on; drawnow
    
    res=norm(f);
    disp([num2str(count) '  ' num2str(res)]);
    if count>666|res>1e5; disp('no convergence');conv=0;break; end
    if res<15; quit=1; disp('converged'); conv=1;continue; end
    
    % Newton step
    sol=sol-A\f;   
    count=count+1;
end

    Plot of the results

    we plot the ideal fluid velocity, the boundary layer thickness and the wall shear stress

    xlabel('x'); ylabel('y');title('Boundary layer');
legend('Ideal Fluid velocity','BL thickness','shear stress/10','geometry');
set(gcf,'paperpositionmode','auto');
print('-dpng','-r75','boundary_layer_pyl.png')
    Evolution of the ideal fluid velocity and the boundary layer thickness

    Evolution of the ideal fluid velocity and the boundary layer thickness

    Exercices/Contributions

    • Please use a closure of H(\delta_1^2 du_e/dx) and f_2(H) and modify the code accordingly
    • Please compare to venturi.m.

    Links

    Bibliography