How to create a shape in MATLAB
 

Hi guys. I'm a student and I've a homework. If you're gonna help me about that I'll be thankful.

I'm studying mechanical engineering. In the lecture of CFD, everybody got a shape of pipe (Actually it's not a pipe. A structure which lets fluid flows inside) We'd to make a 2-d simulation at Star CCM first. Than we must make the same thing in MATLAB and than compare the results. I've already made the Star CCM simulation and got results:

Scalar Scene 1_Velocity Magnitude
Scalar Scene 2_Pressure Coefficient
Vector Scene 1_Cell Relative Velocity

http://www.urimage.net/images/2013/01/12/E538.jpg
http://www.urimage.net/images/2013/01/12/RQln.jpg
http://www.urimage.net/images/2013/01/12/H8N9.jpg
http://www.urimage.net/images/2013/01/12/ShA1e.jpg



Now I've to do the same thing at the MATLAB but I've really no idea how to start. Can you guys please help me?

My pipe is:

http://www.urimage.net/images/2013/01/12/xPuEo.jpg


The exercise is: The flow is frictionless, incompressible, stable and irrotational. Fluid is water and entrance velocity is 8 m/s, and d= 1000 Kg/m3

1) The flow function (ψ) for the water passage must be calculated point by point with Laplace-Equation

2) The velocity compenents U and V must be calculated for all the corners. U = dψ/dy, V = -d/dx.

3) The pressure distribution must be calculated with Bernoulli-equation. Especially the pressure coefficient (Cp) must be calculated for all the points.

I've some old MATLAB codes from kind of related project

%Streamfunction etc.
function [psi]=msteinmetz_a()
dy=1;
dx=1;
ymax=40;
xmax=2000;
yb=25;
xb=100;
eps=0.00001;
dif=eps+1;
Uin=0.8;
Uout=(ymax-1)*dy/((ymax-yb)*dy)*Uin;
psi=zeros(ymax,xmax);

%initial distribution
for i=1:1:xb-1
for j=2:1:ymax
psi(j,i)=psi(j-1,i)+dy*Uin;
end
end
for i=xb:1:xmax
for j=yb+1:1:ymax
psi(j,i)=psi(j-1,i)+dy*Uout;
end
end


% Finding Stable distribution
while (dif > eps)
psibefore=psi;
for j=2:1:yb
for i=2:1:xb-1
psi(j,i)=0.25*(psi(j,i-1)+psi(j,i+1)+psi(j-1,i)+psi(j+1,i));
end
end

for j=yb+1:1:ymax-1
for i=2:1:xmax-1
psi(j,i)=0.25*(psi(j,i-1)+psi(j,i+1)+psi(j-1,i)+psi(j+1,i));
end
end
dif=abs(max(max(psibefore-psi)));

end

%U
U=zeros(ymax,xmax);

for j=2:1:yb
for i=2:1:xb-1
U(j,i)=(psi(j+1,i)-psi(j-1,i))/2*dy;
end
end

for j=yb+1:1:ymax-1
for i=2:1:xmax-1
U(j,i)=(psi(j+1,i)-psi(j-1,i))/2*dy;
end
end

%V
V=zeros(ymax,xmax);

for j=2:1:yb
for i=2:1:xb-1
V(j,i)=(psi(j,i-1)-psi(j,i+1))/2*dx;
end
end

for j=yb+1:1:ymax-1
for i=2:1:xmax-1
V(j,i)=(psi(j,i-1)-psi(j,i+1))/2*dx;
end
end

%Velocity Magnitude
Umag=zeros(ymax,xmax);

for j=2:1:yb
for i=2:1:xb-1
Umag(j,i)=sqrt((U(j,i)^2+V(j,i)^2));
end
end

for j=yb+1:1:ymax-1
for i=2:1:xmax-1
Umag(j,i)=sqrt((U(j,i)^2+V(j,i)^2));
end
end

%cp Field
cp=zeros(ymax,xmax);

for j=2:1:yb
for i=2:1:xb-1
cp(j,i)=1-(Umag(j,i)^2)/2*Uin^2;
end
end

for j=yb+1:1:ymax-1
for i=2:1:xmax-1
cp(j,i)=1-(Umag(j,i)^2)/2*Uin^2;
end
end


%plotting

subplot(2,1,1)
hold on
title({'Velocitiy Vectors'})
quiver(U,V);
hold off

subplot(4,2,5)
hold on
title({'Streamlines'})
contour(psi,20);
hold off

subplot(4,2,7)
hold on
title({'Pressure Coefficient'})
surf(cp);
colormap jet;
hold off

subplot(2,2,4)
hold on
title({'Streamfunction'})
surf(psi);
hold off

Please hellp