## Copyright (C) 2010 Simone Zuccher ## ## This program is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 2, or (at your option) ## any later version. ## ## This program is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ## General Public License for more details. ## ## ## function [U, V] = blasius(y,Uinf,nu,x) ## ## INPUT ## y: vector of wall-normal distance on a generally uneven grid ## Uinf: external streamwise velocity ## nu: kinematic viscosity ## x: distance from the leading edge, ## OUTPUT ## U(y) stream-wise velocity ## V(y) wall-normal velocity ## ## blasius.m solves Blasius' equation ## f f'' + 2f''' = 0 ## with BC ## f(0) = f'(0) = 0 and f'(ymax) = 1 ## by introducing u = f', so as to recast the third-order ODE into system ## _ ## | ## | f u' + 2 u'' = 0 ## < ## | f' - u = 0 ## |_ ## ## with BC ## f(0) = u(0) = 0 and u(ymax) = 1. ## The last equation for f at ymax is ## f'(ymax) - u(ymax) = 0. ## If you experience convergence problems, make sure to use grid nodes ## clustered close to y=0. ## Uinf=1, nu=1 and x=1 return U and V as functions of eta (similarity ## variable. ## Author: Simone Zuccher ## Created: 18 June 2010 ## Last modified: function [U, V] = blasius(y,Uinf,nu,x) %--------------------------------- Preamble -----------------------------------% N=length(y); y=y(:); eta=y*sqrt(Uinf/(nu*x)); % Sanity-check on eta if (min(eta)<0|max(eta)<8|x<=0) disp(['error: blasius.m: eta negative or max(eta) < 8 or x<=0. '... 'Exiting blasius.m...']); U=[]; V=[];return end %-------------------------------- Constants -----------------------------------% % These constants are use in the FD code nv=2; ne=2; uv=0; fv=1; e1=0; e2=1; % A temporary rectangular matrix which will become a sparse diagonal matrix smin=-1; smax=1; sizeJac=[N*ne,nv*(smax-smin+2)]; % md: column of rectangular matrix corresponding to the main diagonal md = -(smin-1)*nv; %--------------------------------- Newton IC ----------------------------------% % Initial guess for Newton: use U'(0) = 0.332 and solve a Cauchy problem F=zeros(size(eta)); U=zeros(size(eta)); j=(2:N-1); d10(j) = 1./(eta(j+1)-eta(j-1)); d10=d10(:); d2p(j) = 2*d10(j)./(eta(j+1)-eta(j)); d2p=d2p(:); d2m(j) = 2*d10(j)./(eta(j)-eta(j-1)); d2m=d2m(:); d20(j) = -d2p(j)-d2m(j); d20=d20(:); U(2)=0.332*eta(2); % key choice for guess for j=2:N-1 F(j)=F(j-1)+(U(j-1)+U(j))/2.*(eta(j)-eta(j-1)); U(j+1)=-(2*d20(j).*U(j)+(2*d2m(j)-d10(j).*F(j)).*U(j-1))./ ... (2*d2p(j)+d10(j).*F(j)); end F(N)=F(N-1)+(U(N-1)+U(N))*(eta(N)-eta(N-1))/2; f=zeros(2*N,1); j=(1:N); eq=(j-1)*ne+1; f(eq+uv)=U; f(eq+fv)=F; j=(1:N); eq=(j-1)*ne+1; f(eq+uv)=U; f(eq+fv)=F; %-------------------------------- Newton loop ---------------------------------% err=1; while(err>1e-10) % Set everything back to zero Jac=zeros(sizeJac); rhs=zeros(sizeJac(1),1); % % Boundary conditions at y=0 j=1; eq = (j-1)*ne+1; rhs(eq+e1)=f(eq+uv)-0; rhs(eq+e2)=f(eq+fv)-0; Jac(eq+uv,md-e1+uv)=1; Jac(eq+fv,md-e2+fv)=1; % % Internal points from 2 to N-1 j=(2:N-1); eq = (j-1)*ne+1; rhs(eq+e1)=f(eq+fv).*d10(j).*(f(eq+uv+nv)-f(eq+uv-nv)) + ... 2*(f(eq+uv+nv).*d2p(j) + f(eq+uv).*d20(j) + f(eq+uv-nv).*d2m(j)); rhs(eq+e2)=d10(j).*(f(eq+fv+nv)-f(eq+fv-nv)) - f(eq+uv); Jac(eq+uv-nv,md-e1+uv-nv)= -f(eq+fv).*d10(j) + 2.*d2m(j); Jac(eq+uv,md-e1+uv)= 2*d20(j); Jac(eq+uv+nv,md-e1+uv+nv)= f(eq+fv).*d10(j) + 2.*d2p(j); Jac(eq+fv,md-e1+fv)= d10(j).*(f(eq+uv+nv)-f(eq+uv-nv)); Jac(eq+fv-nv,md-e2+fv-nv)= -d10(j); Jac(eq+uv,md-e2+uv)= -1; Jac(eq+fv+nv,md-e2+fv+nv)= d10(j); % % Boundary conditions at y=ymax j=N; eq = (j-1)*ne+1; rhs(eq+e1)=f(eq+uv)-1; rhs(eq+e2)=(f(eq+fv)-f(eq+fv-nv))/(eta(N) - eta(N-1))-f(eq+uv); Jac(eq+uv,md-e1+uv)=1; Jac(eq+fv,md-e2+fv)=1/(eta(N) - eta(N-1)); Jac(eq+uv,md-e2+uv)=-1; Jac(eq+fv-nv,md-e2+fv-nv)=-1/(eta(N) - eta(N-1)); % % Make the Jacobian a sparse diagonal matrix from a rectangular one and solve Jac=spdiags(Jac,(1:sizeJac(2))-md,sizeJac(1),sizeJac(1)); ftmp = -Jac\rhs; % Update solution f = ftmp + f; % Check the absolute error err=norm(ftmp); endwhile %---------------------------------- Output ------------------------------------% j=(1:N); eq=(j-1)*ne+1; U=Uinf*f(eq+uv); V=Uinf/2*((y/x).*U-sqrt(nu/(Uinf*x))*f(eq+fv)); endfunction