Issue during for loop.
1 次查看(过去 30 天)
显示 更早的评论
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 个评论
Rik
2022-9-9
You have extremely long expressions. I never trust myself to get all the braces correct. How did you verify that you did?
And did you use the debugger to see step by step what happens to your variables? Did you intend to overwrite M(1) every iteration of j?
回答(1 个)
Cris LaPierre
2022-9-9
编辑:Cris LaPierre
2022-9-9
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
- deltaTheta(j) = abs(theta_save(j) - theta_save(j+1));
- theta_save(i) = theta;
- theta = atan(dy);
- dy = 0.609425 - .185186*x;
- 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 个评论
Cris LaPierre
2022-9-9
编辑:Cris LaPierre
2022-9-9
Good point on atan.
Sounds like your code is working as designed then.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!