Clear Filters
Clear Filters

Quivers/streamlines not perpendicular to contours

6 views (last 30 days)
Trav
Trav on 23 Sep 2024
Commented: Trav on 23 Sep 2024
I've created a mesh representing a superelevation transition along a highway - basically a rotating plane around X=0 as Y increases.
I'm trying to find the longest raindrop flow line, but the problem I'm having is that currently the contours and streamlines/quivers do not align. THat is, the streamlines are not perpendicular/normal to the contours. I'm not sure if the contours are incorrect or if the streamlines are incorrect, or if I'm interpreting the use of a function incorrectly. Some help would be much appreciated!
The code is below:
% Hardcoded Input Variables
left_width = 13.5; % Left width of the road in meters
right_width = 3; % Right width of the road in meters
grade = 0.02; % Longitudinal grade (e.g., 2%)
normal_crossfall = -0.03; % Normal crossfall (e.g., -3%)
superelevated_crossfall = 0.06; % Superelevated crossfall (e.g., 6%)
transition_length = 70; % Length of transition in meters
contour_interval = 0.20; % Contour interval (e.g., 200 mm)
% Mesh Parameters
num_segments = 50; % Number of segments along the transition length
num_width_points = 20; % Number of points across the road width (ensure this is even)
% Generate Mesh Points
x = linspace(-left_width, right_width, num_width_points); % Evenly spaced points across the total width
y = linspace(0, transition_length, num_segments + 1); % Only transition length
% Initialize Elevation Matrix
elevation = zeros(length(x), length(y)); % Adjust size based on unique x
% Compute Elevations
for i = 1:length(y)
segment_length = y(i);
for j = 1:length(x)
% Longitudinal grade component
longitudinal_elevation = grade * segment_length; % Extend grade throughout
% During transition
start_elevation = longitudinal_elevation + (normal_crossfall * (-x(j)));
end_elevation = longitudinal_elevation + (superelevated_crossfall * (-x(j)));
elevation(j, i) = interp1([0, transition_length], [start_elevation, end_elevation], segment_length);
end
end
% Create 2D Contour Plot
[X, Y] = meshgrid(x, y);
figure;
contourf(X, Y, elevation', 'LineColor', 'none'); % Filled contour plot
hold on;
% Add contour lines using the specified contour interval
contour_levels = min(elevation(:)):contour_interval:max(elevation(:));
contour(X, Y, elevation', contour_levels, 'LineColor', 'k'); % Contours in black
colorbar;
title('2D Contour Plot of Superelevation Transition');
xlabel('Width of Road (m)');
ylabel('Length of Transition (m)');
axis equal; % Set equal scaling for both axes
% Calculate gradients for flow direction
[grad_x, grad_y] = gradient(elevation'); % Calculate gradients of elevation
% Normalize gradients for better visualization
magnitude = sqrt(grad_x.^2 + grad_y.^2);
grad_x = -grad_x ./ magnitude; % Reverse direction for upstream to downstream
grad_y = -grad_y ./ magnitude; % Reverse direction for upstream to downstream
% Create Streamline starting from the highest values of Y along the Y-axis
lowest_x = min(x); % Find the lowest value of X
startY = transition_length:-5:0; % Start from downstream to upstream
streamlines = streamline(X, Y, grad_x, grad_y, lowest_x * ones(size(startY)), startY);
% Add quiver plot for direction visualization
quiver(X, Y, grad_x, grad_y, 'Color', 'r', 'AutoScale', 'off'); % Add quivers
hold off;
  1 Comment
Trav
Trav on 23 Sep 2024
This is a snapshot of the figure I get - the thick blue lines are drawn on by me manually and are what I would expect the streamlines and quivers to follow, according to the contours (in black).

Sign in to comment.

Answers (1)

Bruno Luong
Bruno Luong on 23 Sep 2024
Add the command
axis equal
Different aspect ratio in x and y can skew the perpendicularity visually.
  1 Comment
Trav
Trav on 23 Sep 2024
axis equal is already being used. I think I figured it out - had some incorrect calcs which were skewing the results.
Thanks

Sign in to comment.

Categories

Find more on Contour Plots in Help Center and File Exchange

Products


Release

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!