Area Mach Number Relation

I need to plot Mach number (M) as a function of Area ratios (A/A*) for subsonic and supersonic cases. I am trying to use newton-raphson method to iterate and find a converging solution for Mach number (M) at specific area ratios (A/A*), however i would like to be able to call in an array of A/A* inputs. The values for this array would from 0.1 to 10 (or 0.1:0.1:10) .

The function in question is

@(M) (1/M^2)*(((2+gm1*M^2)/gp1)^(gp1/gm1))-ARatio^2;

the variables for this function are;

g = 1.4;

gm1 = g-1;

gp1 = g+1;

Any help would be greatly appreciated!!!

Answer 1

[darova]

Here is what i reached using polyxpoly

g   = 1.4; 
gm1 = g-1; 
gp1 = g+1;
F = @(M) 1./M.^2.*((2+gm1*M.^2)/gp1).^(gp1/gm1);%-ARatio^2;
M = linspace(0.1,3.5);                    % Mach number
A = sqrt( F(M) );                         % A ratio

plot(M,A)                                 % draw function
hold on
    % find Mach number of each A 
for a = linspace(0.1,5,10)
    mm = [0 4];                           % just horizontal line
    aa = [a a];
    [xm,ya] = polyxpoly(mm,aa,M,A);
    plot(xm,ya,'.-r')
end
hold off

[Steven Castrillon]

i get an error for the line 13

[xm,ya] = polyxpoly(mm,aa,M,A);

Also, although this gives me a correct plot i do not want to plot the equation as is, i want to reach the converged Mach number values and plot them against the A/A* ratios. Here is an example

a1.Create a main routine and another external function.

a2. Create a for loop inside the main function. This for loop will iterate the values of A/A*

a3. Call the external function with arguments A/A* and an initial Mach number (M0). Remember that to generate the subsonic side the initial M0 has to be subsonic and for the supersonic solution, the initial guess has to be supersonic.

a4.The external function should return the converged Mach number (M_conv) value for that particular A/A*

a5.

Store the values of A/A*’s and the returned M_conv’s in two arrays.

a6.

Plot M_conv vs. A/A*

[darova]

I used polyxpoly to find Mach number at specific points Aratio

There is also intersections

[Steven Castrillon]

The following code shows me how to get solutions for both sub and supersonic, however, how can i input a range of ARatio instead of just one value?





clear;
clc;

%% INPUTS

% Define some paramters
g   = 1.4;
gm1 = g-1;
gp1 = g+1;


% Define anonymous function with two inputs (M and ARatio)
% - Will be used in the methods below
% - Pass M and ARatio as arguments to AM_EQN to get function value
% -     funVal = AM_EQN(M,ARatio)
AM_EQN = @(M,ARatio) sqrt((1/M^2)*(((2+gm1*M^2)/gp1)^...
                            (gp1/gm1)))-ARatio;

% Solve for Msub and Msup using this area ratio (A/A*)
ARatio = 1.5;

% Error tolerance
errTol = 1e-4;

% Flags for printing iterations to screen
verboseBisection   = 0;
verboseIncremental = 0;


%% SUBSONIC INCREMENTAL SEARCH

% Initial values
dM       = 0.1;                                                             % Initial M step size
M        = 1e-6;                                                            % Initial M value
iConvSub = 0;                                                               % Initial converge index

if (verboseIncremental == 1)
    fprintf('Incremental Search Method: Subsonic\n');
    fprintf('-----------------------------------\n');
end

% Iterate to solve for root
iterMax = 100;                                                              % Maximum iterations
stepMax = 100;                                                              % Maximum step iterations
for i = 1:1:iterMax
    for j = 1:1:stepMax
        
        % Evaluate function at j and j+1
        fj   = AM_EQN(M,ARatio);
        fjp1 = AM_EQN(M+dM,ARatio);
        
        % Print iterations to command window
        if (verboseIncremental == 1)
            fprintf('fj | fjp1: %3.4f\t%3.4f\n',fj,fjp1);
        end
        
        % Update M depending on sign change or not
        % - If no sign change, then we are not bounding root yet
        % - If sign change, then we are bounding the root, update dM
        if (fj*fjp1 > 0)
            M = M + dM;                                                     % Update M
        elseif (fj*fjp1 < 0)
            dM = dM*0.1;                                                    % Refine the M increment
            break;                                                          % Break out of j loop
        end
        
    end % END: j Loop
    
    % Check for convergence
    if (abs(fj-fjp1) <= errTol)                                             % If converged
        iConvSub = i;                                                       % Set converged index
        break;                                                              % Exit loop
    end
    
end % END: i Loop

% Set subsonic Mach number to final M from iterations
Msub = M;

%% SUPERSONIC INCREMENTAL SEARCH

% Initial values
dM       = 1;                                                               % Initial M step size
M        = 1+1e-6;                                                          % Initial M value
iConvSup = 0;                                                               % Initial converge index

if (verboseIncremental == 1)
    fprintf('\nIncremental Search Method: Supersonic\n');
    fprintf('-------------------------------------\n');
end

% Iterate to solve for root
iterMax = 100;                                                              % Maximum iterations
stepMax = 100;                                                              % Maximum step iterations
for i = 1:1:iterMax
    for j = 1:1:stepMax
        
        % Evaluate function at j and j+1
        fj   = AM_EQN(M,ARatio);
        fjp1 = AM_EQN(M+dM,ARatio);
        
        % Print iterations to command window
        if (verboseIncremental == 1)
            fprintf('fj | fjp1: %3.4f\t%3.4f\n',fj,fjp1);
        end
        
        % Update M depending on sign change or not
        % - If no sign change, then we are not bounding root yet
        % - If sign change, then we are bounding the root, update dM
        if (fj*fjp1 > 0)
            M = M + dM;                                                     % Update M
        elseif (fj*fjp1 < 0)
            dM = dM*0.1;                                                    % Refine the M increment
            break;                                                          % Break out of j loop
        end
        
    end % END: j Loop
    
    % Check for convergence
    if (abs(fj-fjp1) <= errTol)                                             % If converged
        iConvSup = i;                                                       % Set converged index
        break;                                                              % Exit loop
    end
    
end % END: i Loop

% Set supersonic Mach number to final M from iterations
Msup = M;

% Print solutions to command window
fprintf('==== INCREMENTAL SEARCH SOLVER ====\n');
fprintf('Msub: %3.4f after %i iterations\n',Msub,iConvSub);
fprintf('Msup: %3.4f after %i iterations\n',Msup,iConvSup);
fprintf('===================================\n\n');

[darova]

I like it

[Steven Castrillon]

yes but can you please help me to introduce an array of values for ARatio?

in the code i provided, ARatio is set as : ARatio = 1.5

when i set it as: ARatio = [0.1:0.1:10] i get an error at

Error in AREAMACH2 (line 62)

if (fj*fjp1 > 0)

Please help