Problem $2$:

Run the entrance flow program of Appendix F for Re$=0\mathrm{.}1,3\mathrm{.}0,$ and $35.$ Refer to presented

results in the lecture, generate:

a$.$ Plots for the stream function.

b$.$ Plots for velocity profiles at x$=+-0\mathrm{.}2$

c$.$ Plots for velocity versus x$//$h for the centerline and stagnation line.

d$.$ Plots the vorticity lines from the wall to the centerline.

Here is the Program in MATLAB:

$\%$ MATLAB$$ PROGRAM $7\mathrm{.}11\mathrm{.}0($R$2010$b$)$

$\%$ program to compute entrance flow

$\%$ vort$($i$0+1,$ jj$)$ given by Stokes Flow

$\%$

$\%$ set input data

$\%$

n$=5$;

i$0=$n$*10$;

ii$=$n$*20$;

jj$=$n$*10$;

al$=4\mathrm{.}0$;

re$=1\mathrm{.}0$;

epsi$=1\mathrm{.}0$e$-4$;

evort$=1\mathrm{.}0$e$-4$;

$\%$

$\%$ allocate memory for vortn$1,$ycord

$\%$

vortn$1=$ zeros$($ii$+1,$jj$+1)$;

ycord $=$ zeros$($jj$+1)$;

$\%$

$\%$ compute dx$,$ beta and f

$\%$

dx$=$al$/$ii;

beta$=$al$*$jj$/$ii;

beta$2=$beta$^2$;

beta$2$plus$1=$beta$2+1\mathrm{.}0$;

dt$=0\mathrm{.}5/((1\mathrm{.}0+$beta$)/$dx$+$beta$2$plus$1*4\mathrm{.}0/$re$/$dx$^2)$;

e$=($cos$($pi$/$double$(2*$ii$+1))+$beta$2*$cos$($pi$/$double$($jj$)))/$beta$2$plus$1$;

eta$=$e$^2$;

f$=2\mathrm{.}0*(1\mathrm{.}0-$sqrt$(1\mathrm{.}0-$eta$))/$eta;

$\%$

$\%$ set initial conditions and boundary values

$\%$

$\%$ initial conditions

u $=$ ones$($ii$+1,$jj$+1)$;

v $=$ zeros$($ii$+1,$jj$+1)$;

vort $=$ zeros$($ii$+1,$jj$+1)$;

psi $=($repmat$((1$:jj$+1)-1,$ii$+1,1))./$jj;

$\%$ boundary conditions on plate

u$($i$0+1$:ii$+1,$jj$+1)=0$;

v$($i$0+1$:ii$+1,$jj$+1)=0$;

psi$($i$0+1$:ii$+1,$jj$+1)=1\mathrm{.}0$;

vort$($i$0+1$:ii$+1,$jj$+1)=3\mathrm{.}0$;

$\%$ boundary conditions at outlet

y $=(0$:$1/$jj:$1)$;

v$($ii$+1,$:$)=0\mathrm{.}0$;

u$($ii$+1,$:$)=1\mathrm{.}5*(1\mathrm{.}0-($y$.^2))$;

psi$($ii$+1,$:$)=1\mathrm{.}5*$y$-0\mathrm{.}5*($y$.^3)$;

vort$($ii$+1,$:$)=3\mathrm{.}0*$y;

$\%$

$\%$ solve for vorticity at interior points

$\%$

isw$2=$true;

while isw$2$

isw$2=$false;

for i$=2$:ii

for j$=2$:jj

delsq$=$vort$($i$+1,$j$)+$vort$($i$-1,$j$)+$beta$2*($vort$($i$,$j$+1)...$

$+$vort$($i$,$j$-1))-2\mathrm{.}0*$beta$2$plus$1*$vort$($i$,$j$)$;

if $($u$($i$,$j$)>0\mathrm{.}0)$

conu$=$u$($i$,$j$)*$vort$($i$,$j$)-$u$($i$-1,$j$)*$vort$($i$-1,$j$)$;

else

conu$=$u$($i$+1,$j$)*$vort$($i$+1,$j$)-$u$($i$,$j$)*$vort$($i$,$j$)$;

end

if $($v$($i$,$j$)>0\mathrm{.}0)$

conv$=$v$($i$,$j$)*$vort$($i$,$j$)-$v$($i$,$j$-1)*$vort$($i$,$j$-1)$;

else

conv$=$v$($i$,$j$+1)*$vort$($i$,$j$+1)-$v$($i$,$j$)*$vort$($i$,$j$)$;

end

dvort$=$dt$/$dx$*(-($conu$+$beta$*$conv$)+2\mathrm{.}0/$re$/$dx$*$delsq$)$;

vortn$1($i$,$j$)=$vort$($i$,$j$)+$dvort;

if dvort $>$ evort

isw$2=$true;

end

end

end

$\%$

$\%$ update vorticity matrix

$\%$

vort$(2$:end$-1,2$:end$-1)=$vortn$1(2$:end$-1,2$:end$-1)$;

$\%$

$\%$ solve for stream function

$\%$

isw$1=$true;

while isw$1$

isw$1=$false;

for i$=2$:ii

for j$=2$:jj

dstr$=$psi$($i$+1,$j$)+$psi$($i$-1,$j$)+$beta$2*$psi$($i$,$j$+1)...$

$+$beta$2*$psi$($i$,$j$-1)-2\mathrm{.}0*$beta$2$plus$1*$psi$($i$,$j$)+$vort$($i$,$j$)*$dx$^2$;

psi$($i$,$j$)=$psi$($i$,$j$)+$f$/2\mathrm{.}0/$beta$2$plus$1*$dstr;

if dstr $>$ epsi

isw$1=$true;

end

end

end

end

$\%$

$\%$ calculate u and v velocities at interior points

$\%$

u$(2$:ii$,2$:jj$)=($psi$(2$:ii$,3$:jj$+1)-$psi$(2$:ii$,1$:jj$-1))*$beta$/2\mathrm{.}0/$dx;

v$(2$:ii$,2$:jj$)=($psi$(1$:ii$-1,2$:jj$)-$psi$(3$:ii$+1,2$:jj$))/2\mathrm{.}0/$dx;

$\%$

$\%$ calculate centerline and stagnation streamline

$\%$ values of u

$\%$

u$(2$:ii$,1)=$psi$(2$:ii$,2)*$beta$/$dx;

u$(2$:i$0,$jj$+1)=(1\mathrm{.}0-$psi$(2$:i$0,$jj$))*$beta$/$dx;

$\%$

$\%$ calculate vorticities on the walls

$\%$

for i$=$i$0+2$:ii

vortemp$=$vort$($i$,$jj$+1)$;

vort$($i$,$jj$+1)=(1\mathrm{.}0-$psi$($i$,$jj$))*2*$beta$2/$dx$^2$;

dvort$=$vort$($i$,$jj$+1)-$vortemp;

if$($dvort $>$ evort$)$

isw$2=$true;

end

end

vort$($i$0+1,$jj$+1)=$vort$($i$0+1,$jj$)+(1\mathrm{.}0-$psi$($i$0+1,$jj$))/$dx$^2*5*$beta$2/4\mathrm{.}0$;

end

omega$=$vort

stm$=$psi

xmin$=-$i$0*$al$/$ii;

xmax$=$al$*($ii$-$i$0)/$ii;

x $=($xmin:al$/$ii:xmax