Issue during for loop.

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)

[Rik]

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?

[Shawn Alcorn]

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.

Answer 1

[Cris LaPierre]

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?

[Shawn Alcorn]

dy is not an angle, so it can't be converted to radians. I have already verified that my theta is correct. dy is the derivative of the initial parabolic function which becomes a linear function. By definition dy is equivalent to the rise over run at each point along the parabolic function, then by using the inverse of tangent you get the angle at that specific x location. Looking into the theta_save variable, the theta follows the curve and approaches 0 radians as it should.

[Cris LaPierre]

Good point on atan.

Sounds like your code is working as designed then.