[Kishore S]

%From JD Anderson textbook for Computational Fluid Dynamics Chap 8,using
%Maccormack technique of space marching PM wave expansion is solved
%numerically
x=0;
spacestep=0;
y=linspace(0,1,41);%This is the computational y domain

Co=0.5;%Courant no
prompt=input('Enter the coeffcient for artificial viscosity:');
Cy=prompt;
digits(4);
for i=1:41 %Boundary Conditions for x=0 is defined here
M(i)=2;
p(i)=101325;
T(i)=286;
D(i)=1.225;
a=(1.4*287*T(i))^0.5;
u(i)=a*2;
v(i)=0;
F1(i)=D(i)*u(i);
F2(i)=D(i)*(u(i)^2)+p(i);
F3(i)=F1(i)*v(i);
F4(i)=3.5*p(i)*u(i)+(F1(i)*(u(i)^2)/2);
G1(i)=0;
G2(i)=0;
G3(i)=p(i);
G4(i)=0;
end
net=[transpose(D) transpose(u) transpose(v) transpose(p) transpose(T) transpose(M)];
disp(net);
while(spacestep<50)
if x>10
h=40+(x-10)*tan(deg2rad(5.352));%Height of the domain which is a function of x
def=5.352;%Deflection Angle
else
h=40;
def=0;
end
%Loop for calculating the marching step size of the domain across x
%direction
for i=1:41
MA=rad2deg(asin(1/M(i)));
pos=abs(tan(deg2rad(5.352+MA)));
neg=abs(tan(deg2rad(5.352-MA)));
a1=max(pos,neg);
delX(i)=Co*h*0.025/a1;
end

delX;
step=min(delX);%Finding the minimum value of step size of the 41 grid pts along y axis
[dF1_dx,dF2_dx,dF3_dx,dF4_dx,F1_p,F2_p,F3_p,F4_p,p_ p,G1_p,G2_p,G3_p,G4_p]=predictor(F1,F2,F3,F4,G1,G2,G3,G4,def,h,Cy,step,p ,y);
[dF1c_dx,dF2c_dx,dF3c_dx,dF4c_dx,SF1c,SF2c,SF3c,SF4 c]=corrector(F1_p,F2_p,F3_p,F4_p,p_p,G1_p,G2_p,G3_p, G4_p,h,Cy,y,def);
dF1avg_dx=(dF1_dx+dF1c_dx)/2;
dF2avg_dx=(dF2_dx+dF2c_dx)/2;
dF3avg_dx=(dF3_dx+dF3c_dx)/2;
dF4avg_dx=(dF4_dx+dF4c_dx)/2;
for i=1:41
F1(i)=F1(i)+SF1c(i)+dF1avg_dx(i)*step;%F1,F2,F3,F4 are the flux terms required to be found
F2(i)=F2(i)+SF2c(i)+dF2avg_dx(i)*step;
F3(i)=F3(i)+SF3c(i)+dF3avg_dx(i)*step;
F4(i)=F4(i)+SF4c(i)+dF4avg_dx(i)*step;
A=((F3(i)^2)/(2*F1(i)))-F4(i);
B=3.5*F1(i)*F2(i);
C=-3*(F1(i)^3);
D(i)=(-B+(B^2-4*A*C)^0.5)/(2*A);%Here the decoding of primitive variables are done
u(i)=F1(i)/D(i);
v(i)=F3(i)/F1(i);

p(i)=F2(i)-F1(i)*u(i);
T(i)=p(i)/(D(i)*287);
magV=(u(i)^2+v(i)^2)^0.5;
M(i)=magV/((1.4*287*T(i))^0.5);
G1(i)=D(i)*v(i);%G1,G2,G3,G4 are the source terms
G2(i)=F3(i);
G3(i)=D(i)*(v(i)^2)+p(i);
G4(i)=3.5*p(i)*v(i)+F1(i)*(v(i)^2)/2;
end
v

%Here the Abbot Boundary Condition has been employed for the wall grid point
if x>10
ang=rad2deg(atan(abs(v(1))/u(1)));
Rot_ang=def-ang;
else
ang=rad2deg(atan(v(1)/u(1)));
Rot_ang=ang;
end
M
a=M(1);
pmf=prandtlmeyer(a);
pmactual=pmf+Rot_ang;
act=Mach(1.4,pmactual);
act_p=p(1)*((1+0.2*a^2)/(1+0.2*act^2))^3.5;
act_T=T(1)*(1+0.2*a^2)/(1+0.2*act^2);
D(1)=act_p/(287*act_T);
p(1)=act_p;
T(1)=act_T;
v(1)=-u(1)*tan(deg2rad(def));
M(1)=act;
%M(1)=((u(1)^2+v(1)^2)^0.5)/((1.4*287*T(1))^0.5);
F1(1)=D(1)*u(1);
F2(1)=D(1)*(u(1)^2)+p(1);
F3(1)=F1(1)*v(1);
F4(1)=3.5*p(1)*u(1)+F1(1)*(u(1)^2)/2;
G1(1)=D(1)*v(1);
G2(1)=F3(1);
G3(1)=D(1)*(v(1)^2)+p(1);
G4(1)=3.5*p(1)*v(1)+F1(1)*(v(1)^2)/2;

M






%net=[transpose(D) transpose(u) transpose(v) transpose(p) transpose(T) transpose(M)];
%disp(net);
x=x+step;
spacestep=spacestep+1;

end

function [dF1_dx,dF2_dx,dF3_dx,dF4_dx,F1_p,F2_p,F3_p,F4_p,p_ p,G1_p,G2_p,G3_p,G4_p] = predictor(F1,F2,F3,F4,G1,G2,G3,G4,def,h,Cy,step,p, y)


for i=1:41

deta_dx=(1-y(i))*tan(deg2rad(def))/h;
if i~=41
dF1_dy=(F1(i)-F1(i+1))/0.025;
dG1_dy=(G1(i)-G1(i+1))/0.025;
dF1_dx(i)=deta_dx*dF1_dy+(dG1_dy/h);
dF2_dy=(F2(i)-F2(i+1))/0.025;
dG2_dy=(G2(i)-G2(i+1))/0.025;
dF2_dx(i)=deta_dx*dF2_dy+(dG2_dy/h);
dF3_dy=(F3(i)-F3(i+1))/0.025;
dG3_dy=(G3(i)-G3(i+1))/0.025;
dF3_dx(i)=deta_dx*dF3_dy+(dG3_dy/h);
dF4_dy=(F4(i)-F4(i+1))/0.025;
dG4_dy=(G4(i)-G4(i+1))/0.025;
dF4_dx(i)=deta_dx*dF4_dy+(dG4_dy/h);

else
dF1_dy=(F1(i-1)-F1(i))/0.025;
dG1_dy=(G1(i-1)-G1(i))/0.025;
dF1_dx(i)=deta_dx*dF1_dy+(dG1_dy/h);
dF2_dy=(F2(i-1)-F2(i))/0.025;
dG2_dy=(G2(i-1)-G2(i))/0.025;
dF2_dx(i)=deta_dx*dF2_dy+(dG2_dy/h);
dF3_dy=(F3(i-1)-F3(i))/0.025;
dG3_dy=(G3(i-1)-G3(i))/0.025;
dF3_dx(i)=deta_dx*dF3_dy+(dG3_dy/h);
dF4_dy=(F4(i-1)-F4(i))/0.025;
dG4_dy=(G4(i-1)-G4(i))/0.025;
dF4_dx(i)=deta_dx*dF4_dy+(dG4_dy/h);

end

if i==1 || i==41
SF1=0;
SF2=0;
SF3=0;
SF4=0;
else %These are the artificial viscosity terms for predictor step
SF1=(abs(p(i+1)-2*p(i)+p(i-1))*Cy*(F1(i+1)-2*F1(i)+F1(i-1)))/(p(i+1)+2*p(i)+p(i-1));
SF2=(abs(p(i+1)-2*p(i)+p(i-1))*Cy*(F2(i+1)-2*F2(i)+F2(i-1)))/(p(i+1)+2*p(i)+p(i-1));
SF3=(abs(p(i+1)-2*p(i)+p(i-1))*Cy*(F3(i+1)-2*F3(i)+F3(i-1)))/(p(i+1)+2*p(i)+p(i-1));
SF4=(abs(p(i+1)-2*p(i)+p(i-1))*Cy*(F4(i+1)-2*F4(i)+F4(i-1)))/(p(i+1)+2*p(i)+p(i-1));
end
F1_p(i)=F1(i)+SF1+dF1_dx(i)*step;
F2_p(i)=F2(i)+SF2+dF2_dx(i)*step;
F3_p(i)=F3(i)+SF3+dF3_dx(i)*step;
F4_p(i)=F4(i)+SF4+dF4_dx(i)*step;
A=((F3_p(i)^2)/(2*F1_p(i)))-F4_p(i);
B=3.5*F1_p(i)*F2_p(i);
C=-3*(F1_p(i)^3);
D_p(i)=(-B+(B^2-4*A*C)^0.5)/(2*A);
u_p(i)=F1_p(i)/D_p(i);
% v_p(i)=F3_p(i)/F1_p(i);
p_p(i)=F2_p(i)-F1_p(i)*u_p(i);
%T_p(i)=p_p(i)/(D_p(i)*287);
G1_p(i)=D_p(i)*F3_p(i)/F1_p(i);
G2_p(i)=F3_p(i);
G3_p(i)=D_p(i)*(F3_p(i)/F1_p(i))^2+F2_p(i)-(F1_p(i)^2/D_p(i));
te=F3_p(i)/F1_p(i);
G4_p(i)=3.5*(F2_p(i)-(F1_p(i)^2/D_p(i)))*te+0.5*D_p(i)*te*((F1_p(i)/D_p(i))^2+te^2);

end

end
function [dF1c_dx,dF2c_dx,dF3c_dx,dF4c_dx,SF1c,SF2c,SF3c,SF4 c] = corrector(F1_p,F2_p,F3_p,F4_p,p_p,G1_p,G2_p,G3_p,G 4_p,h,Cy,y,def)

for i=1:41
deta_dx=(1-y(i))*tan(deg2rad(def))/h;
if i>1
dF1_dy=(F1_p(i-1)-F1_p(i))/0.025;
dG1_dy=(G1_p(i-1)-G1_p(i))/0.025;
dF1c_dx(i)=deta_dx*dF1_dy+(dG1_dy/h);
dF2_dy=(F2_p(i-1)-F2_p(i))/0.025;
dG2_dy=(G2_p(i-1)-G1_p(i))/0.025;
dF2c_dx(i)=deta_dx*dF2_dy+(dG2_dy/h);
dF3_dy=(F3_p(i-1)-F3_p(i))/0.025;
dG3_dy=(G3_p(i-1)-G3_p(i))/0.025;
dF3c_dx(i)=deta_dx*dF3_dy+(dG3_dy/h);
dF4_dy=(F4_p(i-1)-F4_p(i))/0.025;
dG4_dy=(G4_p(i-1)-G4_p(i))/0.025;
dF4c_dx(i)=deta_dx*dF4_dy+(dG4_dy/h);
else
dF1_dy=(F1_p(i)-F1_p(i+1))/0.025;
dG1_dy=(G1_p(i)-G1_p(i+1))/0.025;
dF1c_dx(i)=deta_dx*dF1_dy+(dG1_dy/h);
dF2_dy=(F2_p(i)-F2_p(i+1))/0.025;
dG2_dy=(G2_p(i)-G1_p(i+1))/0.025;
dF2c_dx(i)=deta_dx*dF2_dy+(dG2_dy/h);
dF3_dy=(F3_p(i)-F3_p(i+1))/0.025;
dG3_dy=(G3_p(i)-G3_p(i+1))/0.025;
dF3c_dx(i)=deta_dx*dF3_dy+(dG3_dy/h);
dF4_dy=(F4_p(i)-F4_p(i+1))/0.025;
dG4_dy=(G4_p(i)-G4_p(i+1))/0.025;
dF4c_dx(i)=deta_dx*dF4_dy+(dG4_dy/h);


end
SF1c(1)=0;
SF2c(1)=0;
SF3c(1)=0;
SF4c(1)=0;
for i=2:40
SF1c(i)=(abs(p_p(i+1)-2*p_p(i)+p_p(i-1))*Cy*(F1_p(i+1)-2*F1_p(i)+F1_p(i-1)))/(p_p(i+1)+2*p_p(i)+p_p(i-1));
SF2c(i)=(abs(p_p(i+1)-2*p_p(i)+p_p(i-1))*Cy*(F2_p(i+1)-2*F2_p(i)+F2_p(i-1)))/(p_p(i+1)+2*p_p(i)+p_p(i-1));
SF3c(i)=(abs(p_p(i+1)-2*p_p(i)+p_p(i-1))*Cy*(F3_p(i+1)-2*F3_p(i)+F3_p(i-1)))/(p_p(i+1)+2*p_p(i)+p_p(i-1));
SF4c(i)=(abs(p_p(i+1)-2*p_p(i)+p_p(i-1))*Cy*(F4_p(i+1)-2*F4_p(i)+F4_p(i-1)))/(p_p(i+1)+2*p_p(i)+p_p(i-1));
end
SF1c(41)=0;
SF2c(41)=0;
SF3c(41)=0;
SF4c(41)=0;
%dF1c_dx(41)=0;
%dF2c_dx(41)=0;
%dF3c_dx(41)=0;
%dF4c_dx(41)=0;

end
end
function [pm]=prandtlmeyer(m)
c=(1.4+1)/(1.4-1);

p2=m^2 -1;
p3=rad2deg(atan(sqrt(p2/c)));
p4=rad2deg(atan(sqrt(p2)));
pm=sqrt(c)*p3-p4;
end
function [t]=Mach(g,p)

b=2.9;
a=1.1;
precision=0.0000001;
c=(g+1)/(g-1);
%To find mach,bisection method has been employed
p2=a^2 -1;
r2=((a+b)/2)^2-1;
p3=rad2deg(atan(sqrt(p2/c)));
r3=rad2deg(atan(sqrt(r2/c)));
r4=rad2deg(atan(sqrt(r2)));
p4=rad2deg(atan(sqrt(p2)));
pm=sqrt(c)*p3-p4;
rm=sqrt(c)*r3-r4;
zero1=pm-p;
zero2=rm-p;
while((b-a)/2>precision)

if zero1*zero2<=0
b=(a+b)/2;
else
a=(a+b)/2;
end
p2=a^2 -1;
r2=((a+b)/2)^2-1;
p3=rad2deg(atan(sqrt(p2/c)));
r3=rad2deg(atan(sqrt(r2/c)));
r4=rad2deg(atan(sqrt(r2)));
p4=rad2deg(atan(sqrt(p2)));
pm=sqrt(c)*p3-p4;
rm=sqrt(c)*r3-r4;
zero1=pm-p;
zero2=rm-p;

end
t=(a+b)/2;
end