clc
clear all
%initial coordinate
x=linspace(0,12.5,125);
y1=0.0875*x+0.169;
plot(x,y1,'LineWidth',3)
hold on
y2=-0.0875*x+0.169;
plot(x,y2,'LineWidth',3)
axis([9.5 12.5 -1.5 1.5])
N=5;%number of points of inlet
gamma=1.4;
x1=9.5; y11=1; x2=9.5; y22=0.5; x3=9.5; y3=0; x4=9.5; y4=-0.5; x5=9.5; y5=-1;
% alpha=6;
% theta1=alpha;
% theta2=alpha/2;
% theta3=alpha*0;
% theta4=-alpha/2;
% theta5=-alpha;
lpha=sqrt ((gamma+1) /(gamma-1)) ;
for q=1:20
if q<N+1
alpha=6; theta(1)=alpha;theta(2)=alpha/2;theta(3)=alpha*0;theta(4)=-alpha/2;theta(5)=-alpha;
Ma=1.3;
beta = sqrt (Ma^2 - 1 ) ;
nu(1)=lpha*atan ( beta / lpha )-atan ( beta ) ;nu(2)=nu(1);nu(3)=nu(1);nu(4)=nu(1);nu(5)=nu(1);
M(q)=1.348;
end
if q>N
theta(q)=1/2*(theta(q-N)+theta(q-(N-1))+nu(q-N)+nu(q-(N-1)));
nu(q)=1/2*(theta(q-N)+theta(q-(N-1))+nu(q-N)+nu(q-(N-1)));
yw(q)=(nu(q)/130.45) ^(2/3) ;
M(q)=(1+1.3604*yw(q)+0.0962*(yw(q)^2) -0.5127*(yw(q)^3) ) /(1 -0.6722*yw(q)-0.3278*(yw(q)^2) ) ;
miu(q)=asin(1/M(q));
endK_minus(q)=theta(q)+nu(q) ;
K_plus(q)=theta(q) -nu(q) ;
slope(q) =0.5*( theta(q)+theta(q+N) -miu(q)-miu(q+N)) ;
slope(q+1) =0.5*( theta(q+1)+theta(q+N)+miu(q+1)+mu(q+N)) ;
x(q)=(y(q-(N-1))-y(q-N)+slope(q,q+N)*x(q-N)-slope(q+1,q+N)*x(q-(N-1)))/(slope(q,q+N)-slope(q+1,q+N));
y(q)=slope(q,q+N)*x(q)+y(q-N)-slope(q,q+N)*x(q-N) ;
end
评论4