How to use contourf with 1D data

Hello,

I'm using the method of characteristics to solve 2D supersonic flows inside a nozzle. The discrete solution is thus computed on n different points inside the nozzle.

At the end of the process, I get :

x % a 1D vector (1 x n) of the x coordinates of the points where a solution is computed
y % a 1D vector (1 x n) of the y coordinates of the points where a solution is computed
M % a 1D vector (1 x n) of the Mach number of the flow at the corresponding [x,y] points

I'd like to depict the solution I have on a 2D plot using contourf to obtain something like that :

My problem is that I don't know how to use meshgrid or reshape to do so.

I also would like to use quiver to show the flow direction, but I logically face the same problem...

Can you help me ?

Thanks a lot :-)

Answer 1

[darova]

Use linspace and meshgrid to create regular mesh

xx = linspace(min(xdata),max(xdata),20);
yy = linspace(min(ydata),max(ydata),20);
[X,Y] = meshgrid(xx,yy);

Use griddata to calculate value at each grid point

Z = griddata(xdata,ydata,zdata,X,Y);

And finally use contour to create contour

[Antoine Laterre]

Seems great now! Thanks a lot!

[Lalit Duseja]

Antoine Laterre can you share the code or guide me how you did this. I have a smilar problem. The diameter varies along x and I have to plot the temperature contours.

[Antoine Laterre]

Dear @Lalit Duseja,

As I worked on this project some times ago, I don't remember exactly everything I did. Thus I'll share with you my full code. I hope this will help you!

Best regards,

Antoine

_____________________________________________________________________________________

• First, I declare the initial conditions of my problem and I define the geometry of the nozzle. I then display the characteristics net.

%% Defigning paramaters
n_init = 5;
display = 1;
gamma = 1.4;
Mi = 1.2;
p_design = 0.0836535;
p_atm_p_0 = 1*p_design;
lambda = (gamma+1)/(gamma-1);

%% Nozzle geometry
n = 2000; % number of points used for the horizontal discretization
h_0 = 0.1; % [m]
delta_h = 0.05; % [m]
l = 0.5; % [m]
x_geo = linspace(0,l,n);
h_up = h_0+delta_h*(1-cos(pi.*x_geo./l))/2;
h_low = -h_up;
theta_up = tan(delta_h/2*sin(pi.*x_geo./l).*pi/l);
theta_low = -theta_up;
h_sym = zeros(1,n);

figure; hold on
plot(x_geo./h_up(1),h_up./h_up(1),'r',x_geo./h_up(1),h_low./h_up(1),'r') 
plot(x_geo./h_up(1),h_sym./h_up(1),'--k','LineWidth',1)
title('Characteristics net');
xlabel('Non dimensional length [x/y_0]');
ylabel('Non dimensional height [y/y_0]');
grid on;  

• Next I initialize my numerical model and compute the solution (the subfunctions I use are at the bottom of the comment).

%% Initialization
data = struct();
% mesh
data.index = 1:n_init; % index of point in "mesh"
data.column = zeros(1,n_init)+1; % index of column in "mesh"
data.line = 1:n_init; % index of line in column
% local flow parameters
data.x = zeros(1,n_init)+0; % x coordinates
data.y = zeros(1,n_init)+linspace(h_up(1),h_low(1),n_init); % y coordinates
data.theta = [theta_up(1), zeros(1,n_init-2), theta_low(1)]; % theta (flow angle)
data.M = zeros(1,n_init)+Mi; % M (Mach number)
data.mu = find_Mach_angle(data.M); % mu (Mach angle)
data.nu = find_PM_angle(data.M,lambda); % nu (Prandt-Meyer angle)
data.p = find_pressure(gamma,data.M); % p/p_0
% characteriscs parameters
data.c_eta = data.nu-data.theta; % c_eta (eta charac.)
data.c_xsi = data.nu+data.theta; % c_xsi (xsi charac.)
data.eta = data.theta+data.mu; % eta (angle)
data.xsi = data.theta-data.mu; % xsi (angle)
% equations parameters for characteristics
data.a_eta = tan(data.eta); % a_eta
data.a_xsi = tan(data.xsi); % a_xsi
data.b_eta = data.y-data.a_eta.*data.x; % b_eta
data.b_xsi = data.y-data.a_xsi.*data.x; % b_xsi 
end_column_loop = 0;
i_column = 1;

%% Solution computation
while end_column_loop ~= 1 % column by column construction of the mesh
   find_column = ismember(data.column,i_column);
   n_elem = sum(find_column);
   index = find(find_column);
   i_line = 1;
   is_part_wall = abs(data.y(index(1))-wall_height(data.x(index(1))))<0.001;      
   if is_part_wall % points on walls
       fprintf('Column n°%d \n',i_column)
       fprintf('From index n°%d to index n°%d\n',index(1),index(end))
       for i = index(1:length(index)-1)  
          % coordinates and labels of new point
           [x, y] = find_inter_charac(data.a_xsi(i),data.a_eta(i+1),data.b_xsi(i),data.b_eta(i+1));
           data.x = [data.x, x]; data.y = [data.y, y];
           if display == 1
           plot([data.x(i),x]./h_up(1),[data.y(i),y]./h_up(1),'b',[data.x(i+1),x]./h_up(1),[data.y(i+1),y]./h_up(1),'b','LineWidth',1)
           end
           data.index = [data.index, data.index(end)+1];
           data.line = [data.line, i_line];
           data.column = [data.column, i_column+1];
           % local flow properties
           c_eta = data.c_eta(i+1); data.c_eta = [data.c_eta, c_eta];
           c_xsi = data.c_xsi(i); data.c_xsi = [data.c_xsi, c_xsi];
           nu = (c_xsi+c_eta)/2; data.nu = [data.nu, nu];
           theta = (c_xsi-c_eta)/2; data.theta = [data.theta, theta];
           M = find_Mach(Mi,nu,lambda); data.M = [data.M, M];
           mu =  find_Mach_angle(M); data.mu = [data.mu, mu];
           p = find_pressure(gamma,M); data.p = [data.p, p];
           % characteriscs parameters
           eta = theta+mu; data.eta = [data.eta, eta];
           xsi = theta-mu; data.xsi = [data.xsi, xsi];    
           % charac. equations
           a_eta = tan(eta);  data.a_eta = [data.a_eta, a_eta];
           a_xsi = tan(xsi);  data.a_xsi = [data.a_xsi, a_xsi];
           b_eta = y-a_eta.*x;  data.b_eta = [data.b_eta, b_eta];
           b_xsi = y-a_xsi.*x;  data.b_xsi = [data.b_xsi, b_xsi];
           % next point
           i_line = i_line+1;
      end
   else % no point on walls
       fprintf('Column n°%d \n',i_column)
       fprintf('From index n°%d to index n°%d\n',index(1),index(end))     
       for i = index(1:length(index))               
           if i == index(1) % intersection with the upper wall
               [x, y] = find_inter_upper_wall(h_0,delta_h,l,data.a_eta(i),data.b_eta(i),data.x(i));
               data.x = [data.x, x]; data.y = [data.y, y];
               if display == 1
               plot([data.x(i),x]./h_up(1),[data.y(i),y]./h_up(1),'b','LineWidth',1)
               end
               data.index = [data.index, data.index(end)+1];
               data.line = [data.line, i_line];
               data.column = [data.column, i_column+1];
               % local flow properties
               theta = wall_angle(x); data.theta = [data.theta, theta];
               nu = data.c_eta(i)+theta; data.nu = [data.nu, nu];
               M = find_Mach(Mi,nu,lambda); data.M = [data.M, M];
               mu =  find_Mach_angle(M); data.mu = [data.mu, mu];
               p = find_pressure(gamma,M); data.p = [data.p, p];
               % characteriscs parameters
               c_eta = nu-theta; data.c_eta = [data.c_eta, c_eta];
               c_xsi = nu+theta; data.c_xsi = [data.c_xsi, c_xsi];
               eta = theta+mu; data.eta = [data.eta, eta];
               xsi = theta-mu; data.xsi = [data.xsi, xsi];
               % charac. equations
               a_eta = tan(eta); data.a_eta = [data.a_eta, a_eta];
               a_xsi = tan(xsi); data.a_xsi = [data.a_xsi, a_xsi];
               b_eta = y-a_eta.*x; data.b_eta = [data.b_eta, b_eta];
               b_xsi = y-a_xsi.*x; data.b_xsi = [data.b_xsi, b_xsi];
               % next point
               i_line = i_line+1;
           elseif i == index(end) % intersection with the lower wall               
               [x, y] = find_inter_lower_wall(h_0,delta_h,l,data.a_xsi(i),data.b_xsi(i),data.x(i));
               data.x = [data.x, x]; data.y = [data.y, y];
               if display == 1
               plot([data.x(i),x]./h_up(1),[data.y(i),y]./h_up(1),'b','LineWidth',1)
               end
               data.index = [data.index, data.index(end)+1];
               data.line = [data.line, i_line];
               data.column = [data.column, i_column+1];
               % local flow properties
               theta = -wall_angle(x); data.theta = [data.theta, theta];
               nu = data.c_xsi(i)-theta; data.nu = [data.nu, nu];
               M = find_Mach(Mi,nu,lambda); data.M = [data.M, M];
               mu =  find_Mach_angle(M); data.mu = [data.mu, mu];
               p = find_pressure(gamma,M); data.p = [data.p, p];
               % characteriscs parameters
               c_eta = nu-theta; data.c_eta = [data.c_eta, c_eta];
               c_xsi = nu+theta; data.c_xsi = [data.c_xsi, c_xsi];
               eta = theta+mu; data.eta = [data.eta, eta];
               xsi = theta-mu; data.xsi = [data.xsi, xsi];
               % charac. equations
               a_eta = tan(eta); data.a_eta = [data.a_eta, a_eta];
               a_xsi = tan(xsi); data.a_xsi = [data.a_xsi, a_xsi];
               b_eta = y-a_eta.*x; data.b_eta = [data.b_eta, b_eta];
               b_xsi = y-a_xsi.*x; data.b_xsi = [data.b_xsi, b_xsi];
               % next point
               i_line = i_line+1;
           end 
           if i < index(end)
               % coordinates and labels of new point
               [x, y] = find_inter_charac(data.a_xsi(i),data.a_eta(i+1),data.b_xsi(i),data.b_eta(i+1));
               data.x = [data.x, x]; data.y = [data.y, y];
               if display == 1
               plot([data.x(i),x]./h_up(1),[data.y(i),y]./h_up(1),'b',[data.x(i+1),x]./h_up(1),[data.y(i+1),y]./h_up(1),'b','LineWidth',1)
               end
               data.index = [data.index, data.index(end)+1];
               data.line = [data.line, i_line];
               data.column = [data.column, i_column+1];
               % local flow properties
               c_eta = data.c_eta(i+1); data.c_eta = [data.c_eta, c_eta];
               c_xsi = data.c_xsi(i); data.c_xsi = [data.c_xsi, c_xsi];
               nu = (c_xsi+c_eta)/2; data.nu = [data.nu, nu];
               theta = (c_xsi-c_eta)/2; data.theta = [data.theta, theta];
               M = find_Mach(Mi,nu,lambda); data.M = [data.M, M];
               mu =  find_Mach_angle(M); data.mu = [data.mu, mu];
               p = find_pressure(gamma,M); data.p = [data.p, p];
               % characteriscs parameters
               eta = theta+mu; data.eta = [data.eta, eta];
               xsi = theta-mu; data.xsi = [data.xsi, xsi];
               % charac. equations
               a_eta = tan(eta); data.a_eta = [data.a_eta, a_eta];
               a_xsi = tan(xsi); data.a_xsi = [data.a_xsi, a_xsi];
               b_eta = y-a_eta.*x; data.b_eta = [data.b_eta, b_eta];
               b_xsi = y-a_xsi.*x; data.b_xsi = [data.b_xsi, b_xsi];
               % next point
               i_line = i_line+1;
           end
       end
   end
   % stopping process
   if n_elem < 2
      end_column_loop = 1; 
   end
   if data.x(end) > x_geo(end)
       end_column_loop = 1;
   end  
   if i_column >= 500
       end_column_loop = 1; 
   end
   % next column
   i_column = i_column+1; 
end

• After I post-process my data to prepare the plot.

%% Post-procesing
line = n_init;
x_mesh = zeros(line,max(data.column))+NaN; % x-coordinates of the nozzle
y_mesh = zeros(line,max(data.column))+NaN; % y-coordinates of the nozzle
M_mesh = zeros(line,max(data.column))+NaN; % Mach number in the nozzle
p_mesh = zeros(line,max(data.column))+NaN; % pressure in the nozzle
theta_mesh = zeros(line,max(data.column))+NaN; % flow angle in the nozzle
for i = 1:max(data.column)
    find_column = ismember(data.column,i);
    n_elem = sum(find_column);
    index = find(find_column);
    for j = 1:n_elem
        x_mesh(j,i) = data.x(index(j));
        y_mesh(j,i) = data.y(index(j));
        M_mesh(j,i) = data.M(index(j));
        p_mesh(j,i) = data.p(index(j));
        theta_mesh(j,i) = data.theta(index(j));
    end
end

• Finally I can make my plots.

%% Plots for first question
scatter(data.x./h_up(1),data.y./h_up(1),10,'b','filled')     
xlim([0,5.5]);
ylim([-1.5,1.5]);

% Points
figure;
subplot(2,1,1)
hold on
plot(x_geo,h_up,'r',x_geo,h_low,'r',x_geo,h_sym,'--k','LineWidth',1)
scatter(data.x,data.y,10,'b','filled')
title('Nozzle profile');
xlabel('Length [m]');
ylabel('Height [m]');
grid on;
subplot(2,1,2)
hold on
plot(x_geo,theta_up,'r',x_geo,theta_low,'r',x_geo,h_sym,'--k','LineWidth',1)
scatter(data.x,data.theta,10,'b','filled')
xlabel('Length [m]');
ylabel('Theta [rad]');
grid on;

% Mach and pressure : contourf
figure; hold on;
subplot(2,1,1); hold on;
contourf(x_mesh./h_up(1),y_mesh./h_up(1),M_mesh,25,'LineWidth',0.1); hold on;
c = colorbar;
c.Label.String = 'Mach number';
plot(x_geo./h_up(1),h_up./h_up(1),'-k',x_geo./h_up(1),h_low./h_up(1),'-k','LineWidth',3)
title('Flow inside the nozzle : Mach number');
xlabel('Non dimensional length [x/y_0]');
ylabel('Non dimensional height [y/y_0]');
xlim([0,5.5]);
ylim([-1.5,1.5]);
subplot(2,1,2); hold on;
colormap('jet'); 
contourf(x_mesh./h_up(1),y_mesh./h_up(1),p_mesh./find_pressure(gamma,Mi),25,'LineWidth',0.1); hold on;
c = colorbar;
c.Label.String = 'pressure (p/p_{in})';
plot(x_geo./h_up(1),h_up./h_up(1),'-k',x_geo./h_up(1),h_low./h_up(1),'-k','LineWidth',3)
title('Flow inside the nozzle : pressure (p/p_{in})');
xlabel('Non dimensional length [x/y_0]');
ylabel('Non dimensional height [y/y_0]');
xlim([0,5.5]);
ylim([-1.5,1.5]);

% Mach and pressure : surf
figure; hold on;
subplot(2,1,1); hold on;
s1 = surf(x_mesh./h_up(1),y_mesh./h_up(1),M_mesh);
s1.EdgeColor = 'none';
c1 = colorbar;
c1.Label.String = 'Mach number';
plot(x_geo./h_up(1),h_up./h_up(1),'-k',x_geo./h_up(1),h_low./h_up(1),'-k','LineWidth',3)
title('Flow inside the nozzle : Mach number');
xlabel('Non dimensional length [x/y_0]');
ylabel('Non dimensional height [y/y_0]');
xlim([0,5.5]);
ylim([-1.5,1.5]);
subplot(2,1,2); hold on;
s2 = surf(x_mesh./h_up(1),y_mesh./h_up(1),p_mesh./find_pressure(gamma,Mi));
s2.EdgeColor = 'none'; hold on;
c2 = colorbar;
c2.Label.String = 'pressure (p/p_{in})';
plot(x_geo./h_up(1),h_up./h_up(1),'-k',x_geo./h_up(1),h_low./h_up(1),'-k','LineWidth',3)
title('Flow inside the nozzle : pressure (p/p_{in})');
xlabel('Non dimensional length [x/y_0]');
ylabel('Non dimensional height [y/y_0]');
xlim([0,5.5]);
ylim([-1.5,1.5]);

% Velocity - Streamlines
figure; hold on;
Mgrid_x = cos(theta_mesh).*M_mesh;
Mgrid_y = sin(theta_mesh).*M_mesh;
quiver(x_mesh./h_up(1),y_mesh./h_up(1),Mgrid_x,Mgrid_y,'b'); hold on;
plot(x_geo./h_up(1),h_up./h_up(1),'-k',x_geo./h_up(1),h_low./h_up(1),'-k','LineWidth',3)
title('Flow inside the nozzle : direction');
xlabel('Non dimensional length [x/y_0]');
ylabel('Non dimensional height [y/y_0]');
xlim([0,5.5]);
ylim([-1.5,1.5]); 

% Comparison with compressible flow 1D
index = (n_init+1)/2:2*n_init-1:data.index(end);
x_2D = data.x(index);
Mx_2D = data.M(index);
x_2Dmean = x_2D(1);
Ai = pi*h_up(1)^2;
Astar = (((gamma+1)/2)/(1+((gamma-1)/2)*Mi^2))^((gamma+1)/(2*(gamma-1)))*Ai;
Mx_1D = zeros(1,length(x_geo));
Mx_2D_mean = zeros(1,length(index));
for i = 1:length(index)
    index_x_i = find(abs(data.x-x_2D(i)) < 0.002);
    Mx_2D_mean(i) = mean(data.M(index_x_i));
end
for i = 1:length(x_geo)
    A = pi*h_up(i)^2;
    Mx_1D(i) = fsolve(@(M) ((Astar/A)-(((gamma+1)/2)/(1+((gamma-1)/2)*M^2))^((gamma+1)/(2*(gamma-1)))),...
        Mi,optimset('Display','off'));
end
figure; hold on;
scatter(x_2D./h_up(1),Mx_2D,'b','filled');
plot(x_2D./h_up(1),Mx_2D_mean,'--k','LineWidth',2);
plot(x_geo./h_up(1),Mx_1D,'-r','LineWidth',2);
title('Flow inside the nozzle : Mach number');
xlabel('Non dimensional length [x/y_0]');
ylabel('Mach number [-]');
legend('2D, symmetry axis','2D, mean','1D','location','best');
xlim([0,5]);
grid on;

figure; hold on;
index_y = find(abs(data.x-x_geo(end)) < 0.002);
y_2D = data.y(index_y);
y_1D = -1.5:0.01:1.5;
My_2D = data.M(index_y);
My_2D_mean = zeros(1,length(y_2D))+mean(My_2D);
My_1D = zeros(1,length(y_1D))+Mx_1D(end);
scatter(y_2D./h_up(1),My_2D,'b','filled');
plot(y_2D./h_up(1),My_2D_mean,'--k','LineWidth',2);
plot(y_1D./h_up(1),My_1D,'-r','LineWidth',2);
title('Flow inside the nozzle : Mach number');
xlabel('Non dimensional height [y/y_0]');
ylabel('Mach number [-]');
legend('2D','2D, mean','1D','location','best');
xlim([-1.5,1.5]);
grid on;

• All the sub-functions I used are hereunder.

%% Subfunctions
% flow paramters
function [nu] = find_PM_angle(M,lambda)
    nu = sqrt(lambda)*atan(sqrt((M.^2-1)/lambda))-atan(sqrt(M.^2-1));
end
function [M] = find_Mach(Mtest,nu,lambda)
   M = fsolve(@(M) nu-(sqrt(lambda)*atan(sqrt((M.^2-1)/lambda))-atan(sqrt(M.^2-1))),...
       Mtest,optimset('Display','off'));
   M = abs(M);
end
function [mu] = find_Mach_angle(M)
    mu = asin(1./M);
end
function [x, y] = find_inter_charac(a_xsi_1,a_eta_2,b_xsi_1,b_eta_2)
    x = -(b_eta_2-b_xsi_1)/(a_eta_2-a_xsi_1);
    y = a_xsi_1*x+b_xsi_1;
end
function [p] = find_pressure(gamma,M)
    p = (1+((gamma-1)/2)*M.^2).^(-gamma/(gamma-1));
end
function [M] = find_Mach_pressure(p_p0,gamma,Mtest)
   M = fsolve(@(M) p_p0-(1+((gamma-1)/2)*M.^2).^(-gamma/(gamma-1)),...
       Mtest,optimset('Display','off'));
   M = abs(M);
end
% geometry, intersection with : walls, contact surface, symmetry surface
function [h] = wall_height(x)
   h_0 = 0.1; % [m]
   delta_h = 0.05; % [m]
   l = 0.5; % [m]
   h = h_0+delta_h*(1-cos(pi*x/l))/2;
end
function [theta] = wall_angle(x)
   delta_h = 0.05; % [m]
   l = 0.5; % [m]
   theta = atan(delta_h/2*sin(pi*x/l)*pi/l); % [rad]
end
function [x, y] = find_inter_upper_wall(h_0,delta_h,l,a_eta,b_eta,xtest)
   x = fsolve(@(x) (h_0+delta_h*(1-cos(pi.*x./l))/2)-(a_eta*x+b_eta),...
       xtest,optimset('Display','off'));
   y = a_eta*x+b_eta;
end
function [x, y] = find_inter_lower_wall(h_0,delta_h,l,a_xsi,b_xsi,xtest)
   x = fsolve(@(x) -(h_0+delta_h*(1-cos(pi.*x./l))/2)-(a_xsi*x+b_xsi),...
       xtest,optimset('Display','off'));
   y = a_xsi*x+b_xsi;
end
function [x, y, theta] = find_inter_wall_spike(x_wall,y_wall,theta_wall,a_xsi,b_xsi)
    diff = abs(y_wall-(a_xsi.*x_wall+b_xsi));
    mins = mink(diff,2);
    i1 = find(diff==mins(1));
    i2 = find(diff==mins(2));
    if i1 == i2
       x = mean([x_wall(i1(1)),x_wall(i2(2))]);
       y = mean([y_wall(i1(1)),y_wall(i2(2))]);
       theta = mean([theta_wall(i1(1)),theta_wall(i2(2))]);
    else
       x = mean([x_wall(i1),x_wall(i2)]);
       y = mean([y_wall(i1),y_wall(i2)]);
       theta = mean([theta_wall(i1),theta_wall(i2)]);
    end
end
function [x, y] = find_inter_free_surf_spike(a_eta,b_eta,x_free,y_free,theta_free)
    x = (y_free-tan(theta_free)*x_free-b_eta)/(a_eta-tan(theta_free));
    y = a_eta*x+b_eta;
end
function [h] = wall_height_spike(x,x_wall,y_wall)
    diff = abs(x_wall-x);
    mins = mink(diff,2);
    i1 = find(diff==mins(1));
    i2 = find(diff==mins(2));
    if i1 == i2
       h = mean([y_wall(i1(1)),y_wall(i2(2))]);
    else
       h = mean([y_wall(i1),y_wall(i2)]); 
    end
end