Troubleshooting my code ( ~15 lines)

I am using the following code:

syms B1 B2 SWR Cv 

i = 1;
Cv(1,2) = 0.39;
B1(1,2) = 55;
B2(1,2) = 29;

while i<61
    i = i + 1;

      Cv(i,2)= Cv(i-1) + 0.01;

      R = 0.5 == Cv(i,2)./2*(tan(B1*pi/180) + tan(B2*pi/180));
      Cp = 0.5 == 1- (1+tan(B2*pi/180).^2)/(1+tan(B1*pi/180).^2); 

      [B1(i,2),B2(i,2)] = vpasolve([R, Cp], [B1,B2]); 
      SWR = Cv(i,2)*(tan(B1(i,2)*pi/180) - tan(B2(i,2)*pi/180)); 

      X(i,2) = Cv(i,2);
      Y(i,2) = SWR;

end

And it is giving me the following errors:

Warning: The system is inconsistent. Solution does not
exist. 
> In symengine
  In sym/privBinaryOp (line 908)
  In  /  (line 309)
  In gasturbine_hw7_try2 (line 16) 
Error using sym.getEqnsVars>checkVariables (line 92)
The second argument must be a vector of symbolic
variables.

Error in sym.getEqnsVars (line 54)
    checkVariables(vars);

Error in sym/vpasolve (line 132)
    [eqns,vars] = sym.getEqnsVars(varargin{1:N});

Error in gasturbine_hw7_try2 (line 18)
    [B1(i,2),B2(i,2)] = vpasolve([R, Cp], [B1,B2]);

I am having trouble understanding why the errors are occurring. What I am trying to do is solve equation R and Cp for increasing values of Cv. So I am simply making an array for the variables B1, B2 and Cv and using the values in each index to solve for a specific value of SWR. I know you can do this using anonymous functions, but I am having trouble doing it that way and so am doing it this way. I seek help in this method, not the other.

[Walter Roberson]

Why do you sometimes access Cv with one index and sometimes access it with two indices?

 Cv(i,2) = Cv(i-1) + 0.01;

Answer 1

[Star Strider]

What was the problem with my fsolve approach? Using fsolve is much more efficient than using the Symbolic Math Toolbox for iterative problems.

Note that:

R = 0.5 == Cv(i,2)./2*(tan(B1*pi/180) + tan(B2*pi/180));

is equivalent to:

R = 0.5 == Cv(i,2)*(tan(B1*pi/180) + tan(B2*pi/180))./2;

You need to be sure your parentheses are nested correctly.

[Vidhan Malik]

For some reason it gave me incorrect values even when I initialized it properly and I still don't completely understand it all to debug it, so that's why I went back to using this method.

[Star Strider]

I would prefer that we stay with fsolve and discover the problem with your assignments and anonymous functions. If you want to find the zeros of your system (that seems to be your intent), fsolve is much more efficient.

I know almost nothing about jet engines, so I can’t help you with that. (I’m an Instrument-Rated Private Pilot, but with one exception, have flown only recips. My areas of expertise are in clinical medicine and biomedical engineering.) I can help you use the Symbolic Math Toolbox to derive your equations correctly if you lead me through what you need to do, and then use fsolve to calculate the result.

[Vidhan Malik]

Thank you for putting so much effort in! So currently, the equations are in fact correct, I've double checked with my professor about it. Hence, I can only see that there is something wrong with the code. The graph should be like the right side of a upside-down parabola (non-linear negative function) but instead we are getting a positive linear line. And I haven't been able to determine where the problem lies when using your method.

So what I have been wanting matlab to do is solve equations R and Cp simultaneously for B1 and B2 without having to initialize it (I have used the same equations without initializing and manually varied the value of Cv and I do get the correct values of B1 and B2, but just can't seem to be able to figure out how to do it in a loop).

[Star Strider]

My pleasure.

You will likely need to experiment with various initial estimates for your parameters to get the set that will give you the result you want. Since your function is periodic, there are an infinite number of possible solutions, both positive and negative. Using the Symbolic Math Toolbox to estimate them may provide a usable set, but I cannot guarantee that approach will be any more successful.

[Vidhan Malik]

I am trying to change the value of the following:

 Cvv = Cvv = [0.3:0.1:101];

So that the graph shows from 0.4 to 1 but I keep getting the error that

Error using plot
Vectors must be the same length.

Error in s2 (line 12)
plot(Cvv, SWR)

How do you change that?

[Vidhan Malik]

Cvv = [0.4:0.01:1];
R = @(B1,B2,Cv) -0.5 + Cv/2*(tand(B1) + tand(B2)); %%Eq 1
Cp = @(B1,B2,Cv) -0.5 + (1+tand(B2).^2)/(1+tand(B1).^2); %%Eq2
% MAPPING: A1 = p(1),    B2 = p(2)
Eqn = @(p,Cv) [R(p(1),p(2),Cv), Cp(p(1),p(2),Cv)];
for k1 = 1:length(Cvv)
    Cv = Cvv(k1);
    [P(:,k1),Fval] = fsolve(@(p)Eqn(p,Cv), [58.667766321328347724797896587002; 42.659761905600685286031836767658]);
    SWR(k1) = Cv*(tand(P(1,k1)) - tand(P(2,k1)));
end
figure(1)
plot(Cvv, SWR)
grid
xlabel('C_v')
ylabel('SWR')

What I am currently using that gives that error.

[Star Strider]

That code runs for me without error.

I can’t explain the reason it’s giving you a positive rather than a negative parabolic curve, unless the angles are reversed in the ‘SWR’ call. Perhaps it should be:

SWR(k1) = Cv*(tand(P(2,k1)) - tand(P(1,k1)));

instead?

Guessing.

[Vidhan Malik]

What changes would I have to make so that the X-axis goes from 0.4 to 1?

[Star Strider]

It does with that code. At least it did when I ran it.

[Vidhan Malik]

That's really weird that we would get two different graphs and values with the same code

[Star Strider]

In R2016a, this code (note the argument switch in ‘SWR’):

Cvv = [0.4:0.01:1];
R = @(B1,B2,Cv) -0.5 + Cv/2*(tand(B1) + tand(B2)); %%Eq 1
Cp = @(B1,B2,Cv) -0.5 + (1+tand(B2).^2)/(1+tand(B1).^2); %%Eq2
% MAPPING: A1 = p(1),    B2 = p(2)
Eqn = @(p,Cv) [R(p(1),p(2),Cv), Cp(p(1),p(2),Cv)];
for k1 = 1:length(Cvv)
    Cv = Cvv(k1);
    [P(:,k1),Fval] = fsolve(@(p)Eqn(p,Cv), [58.667766321328347724797896587002; 42.659761905600685286031836767658]);
    SWR(k1) = Cv*(tand(P(2,k1)) - tand(P(1,k1)));
end
figure(1)
plot(Cvv, SWR)
grid
xlabel('C_v')
ylabel('SWR')

gives me this plot:

[Vidhan Malik]

I'm really perplexed how that worked! Initially I was getting the linear positive function, but by switching the part in the SWR function, it gives the correct graph! I'm definitely going to have to spend some time understanding this code.

THANK YOU SO MUCH FOR YOUR HELP! I honestly really appreciate it. And I have learned more in the last days than what I knew all together before, thank you!

[Star Strider]

As always, my pleasure!

I learn as well, although I wish I knew a bit more about jet engines just now! I’ve just been guessing.

[Vidhan Malik]

Ask any questions about them, I should be able to answer most of them!

In this case, B1 represents the angle of attack of the air as it comes into the rotor blade and B2 is the angle that the air leaves the rotor blade (angles are relative to the blade itself). This allows you to create velocity triangles and those shape the design of the blade itself!

[Vidhan Malik]

This is what a singular compressor stage looks like where the blade between 1 and 2 is the rotor blade which revolves. The blade between 2 and 3 is the stator blade which remains stationary and corrects the flow coming out of the stator blade and gets the flow of air it ready to go into the next stage.

This is a crude drawing of the velocity triangle for the air entering the compressor stage (like wise you will have an additional triangle for the air leaving the blade). In the triangle, W represents the speed of air relative to the blade, C is the absolute velocity of the air, and U is the speed at which the rotor blades revolve!

Finally, this is the velocity triangle that you get for just the rotor blade (The triangle on the bottom, triangle1 is the velocity triangle of the air entering the rotor blade. The triangle on the top, triangle2 is the velocity triangle of the air leaving the rotor blade. Like wise, Beta1 is the angle at which the air enters the rotor stage, and Beta2 is the angle at which the air leaves the rotor stage, the angles being relative to the rotor blade)! Cz simply represents the horizontal component of the velocity of air (axial velocity). The flow coefficient with applied assumptions and simplification turns out to be Cv=Cz/u. As you guess, this is the parameter I was changing in my equations! As the value of Cv is changed, you can see that the geometry of the triangles would change too, giving you new Beta1 and Beta2 values so that is exactly what I am finding out using this code! Fun stuff right? :)

[Star Strider]

Definitely fun stuff! I understand the basics because of the aerodynamics I learned on my own in order to understand aircraft stability and control. I now know more than I ever though I would about the compressor stage of jet engines. Before, I just marveled at them and wondered what went into their design. Thank you for the explanation!

[Vidhan Malik]

I usually spend any free time I have on this channel haha

Youtube Channel devoted to Jet Engines!

It's always amazing to know that there is still so much to learn!

[Vidhan Malik]

I think the best part about this is that there is still something wrong with the graph :D. Having a negative stage work ratio kind of doesn't make sense haha, I guess its back to figuring out what's wrong!

[Star Strider]

I’ll take a look at the YouTube channel. I just never took time to explore jet engines in more depth than in my aviation books. (My interest is in control theory, primarily with respect to analysing physiologic homeostasis, but I got interested in aircraft control when I started flying, and engines are part of that.)

I can’t determine what a stage work ratio should be, but the sign could be simply a matter of the reference frame. For example if something is coming toward you, its sign would be positive, if going away from you, negative. The maths would be the same, except for the sign.

[Vidhan Malik]

Really interesting! I would imagine the nerve endings throughout out body are the probes which pick up on things like temperatures and pressures (bad analogy but something of the sort?)

[Star Strider]

That’s part of it, particularly with muscles, but even those turn out to be extremely difficult to model. Most of my control research involved hormones, energy metabolism, and such. (I trained as an endocrinologist before I went back to get my M.S. in biomedical engineering.) It turns out to be a virtually impossible task because of the complexity and adaptability of physiological control systems. They are inherently both uncontrollable and unobservable, making them extremely difficult to model with any accuracy. I take some comfort in that others have the same problems I had (and have), so it’s the nature of the systems themselves, and not any particular inability on my part, that I also have problems with them.

[Vidhan Malik]

Weird idea but you could flood the particular system with a fluid that allows you to map its surroundings in real life, which would allow you to see and make a model of it potentially?

[Star Strider]

Unfortunately, no. (Besides, that’s already part of all living systems.) With respect to muscles, the problem is separating the transducer dynamics from the muscle dynamics. This works for isometric models (such as the Hill model), but muscles don’t work isometrically, they work isotonically. That’s the problem.

With respect to endocrine control, the problem is that there are too many small but important systems (such as intracellular regulatory systems) to measure reliably, and those that are measurable tend to require very financially expensive and physically invasive instrumentation.

I think about this a lot in my spare time, and have as yet to come up with a solution. (I was an undergraduate Chemistry major, so I have that additional background, but as yet no solutions.)

[Vidhan Malik]

Lost me in the second sentence haha, but no worries. Sounds like you have quite the challenge on your hands!

[Star Strider]

It’s relatively esoteric physiology and biochemistry, likely on par with the detailed fluid dynamics of jet engines. The point is that there are a lot of components to physiological regulatory systems, not all of which are even known, and most of which are important but are both uncontrollable and unobservable, at least in real time.

[Vidhan Malik]

Only until you figure it out! :)