Issue during for loop.

3 views (last 30 days)
Shawn Alcorn
Shawn Alcorn on 8 Sep 2022
Edited: Cris LaPierre on 9 Sep 2022
The issue occurs during the for loop using j. When I calculate deltaTheta it should be pretty constant at the start but eventually decrease in value, problem is my deltaTheta continues to grow by a constant value. The next issue is my M array. The initial value is correct but the while loop is suppose to compare the values until the error is met but my M is decreasing at a constant rate when it should be increasing until a final value of about 2.4. I have verified that all values up to this point are accurate.
%% AEM 624 Hypersonic Flow HW 2
clear all
close all
clc
M1 = 3; %Upstream Mach Number
gamma = 1.4;
P1 = 101325; %Pa
T1 = 298.15;%K
i = 1;
V = M1*343;
q = (gamma/2)*P1*(M1^2);
CpMax = (2/(gamma*(M1^2)))*((((((gamma+1)^2)*(M1^2))/((4*gamma*(M1^2))-(2*(gamma-1))))^(gamma/(gamma-1)))*((1-gamma+(2*gamma*(M1^2)))/(gamma+1))-1);%Modified Newtonian
for x = .01370207:.001:3.291
x_save(i) = x;
y = -0.008333 + (0.609425*x) - (0.092593*(x^2));
dy = 0.609425 - .185186*x;
y_save(i) = y;
theta = atan(dy);
theta_save(i) = theta;
thetaTan = tan(theta);
thetaTan_save(i) = thetaTan;
betaIterDeg = 1;
betaIterRad = betaIterDeg*(pi/180);
thetaTanIter = 2*cot(betaIterRad)*(((M1^2)*(sin(betaIterRad)^2)-1)/(((M1^2)*(gamma+cos(2*betaIterRad)))+2));
while abs(thetaTan - thetaTanIter) > .001
betaIterDeg = betaIterDeg + .001;
betaIterRad = betaIterDeg * (pi/180);
thetaTanIter = 2*cot(betaIterRad)*(((M1^2)*(sin(betaIterRad)^2)-1)/(((M1^2)*(gamma+cos(2*betaIterRad)))+2));
end
beta(i) = betaIterRad;
CpTanWedge = (4/(gamma+1))*((sin(betaIterRad)^2)-(1/(M1^2)));% Tangent Wedge
Cp = 2*((sin(theta))^2); %Newtonian
CpModNewt = CpMax*(sin(theta)^2);% Modified Newtonian
p = (q*Cp)+P1;% Newtonian
pModNewt = (q*CpModNewt)+P1;% Modified Newtonian
pTanWedge = (q*CpTanWedge)+P1;% Tangent Wedge
pnorm(i) = p./P1;% Newtonian
pModNewtNorm(i) = pModNewt/P1;% Modified Newtonian
pTanWedgeNorm(i) = pTanWedge/P1;% Tangent Wedge
i = i + 1;
end
for j = 1:3277
M(1) = (1/sin(beta(1)-theta_save(1)))*sqrt((1+((gamma-1)/2)*(M1^2)*(sin(beta(1))^2))/((gamma*(M1^2)*(sin(beta(1))^2))-((gamma-1)/2)));
PrandtlM(j) = sqrt((gamma+1)/(gamma-1))*(atan(sqrt(((gamma-1)/(gamma+1))*((M(j)^2)-1))))-atan(sqrt((M(j)^2)-1));
deltaTheta(j) = abs(theta_save(j) - theta_save(j+1));
deltaTheta_save(j) =deltaTheta(j);
PrandtlM(j+1) = PrandtlM(j) + deltaTheta(j);
MIter = 1;
PrandtlIter = (sqrt((gamma+1)/(gamma-1))*(atan(sqrt(((gamma-1)/(gamma+1))*((MIter^2)-1)))))-atan(sqrt((MIter^2)-1));
while abs(PrandtlM(j+1)-PrandtlIter) > .001
MIter = MIter + .001;
PrandtlIter = sqrt((gamma+1)/(gamma-1))*(atan(sqrt(((gamma-1)/(gamma+1))*((MIter^2)-1))))- atan(sqrt((MIter^2)-1));
end
j = j + 1;
M(j) = MIter;
end
for k = 1:3277
pipn = ((2*gamma)/(gamma+1))*((M(k)^2)*(sin(beta(k))^2)-1);
CpExpansion = (2/(gamma*(M(1)^2)))*(pipn-1);
pExpansion = (CpExpansion*q)+P1;
pExpansionNorm(k) = pExpansion./P1;
k = k + 1;
end
plot(x_save,y_save)
hold on
plot(x_save,(-1)*y_save)
figure
plot(x_save,pnorm)
hold on
plot(x_save,pModNewtNorm)
hold on
plot(x_save,pTanWedgeNorm)
hold on
plot(x_save(2:3278),pExpansionNorm)
  2 Comments
Shawn Alcorn
Shawn Alcorn on 9 Sep 2022
I verified by using the expressions from a previous working code that I used last week. The main problem is my deltaTheta keeps increasing, the values should be almost constant for a while and then begin to approach 0.

Sign in to comment.

Answers (1)

Cris LaPierre
Cris LaPierre on 9 Sep 2022
Edited: Cris LaPierre on 9 Sep 2022
We are not familiar with the problem you are trying to solve, so can't comment on what the code 'should' be doing. We can only comment on what it is doing.
Follow your calculation steps to understand why deltaTheta is behaving the way it is. Working backwards
  1. deltaTheta(j) = abs(theta_save(j) - theta_save(j+1));
  2. theta_save(i) = theta;
  3. theta = atan(dy);
  4. dy = 0.609425 - .185186*x;
  5. x = .01370207:.001:3.291
So you create an evenly spaced vector, x, and then the equation of a line to calculate dy. You then take atan to get theta.
x = .01370207:.001:3.291;
dy = 0.609425 - .185186*x;
theta = atan(dy);
plot(theta)
I notice you treat your angles as degrees elsewhere in your code, converting them to radians before using trigonometric functions on them (the ones you use expect radians). Should dy also be converted to radians before calculating atan?
  2 Comments
Cris LaPierre
Cris LaPierre on 9 Sep 2022
Edited: Cris LaPierre on 9 Sep 2022
Good point on atan.
Sounds like your code is working as designed then.

Sign in to comment.

Tags

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by