Finding axial, radial velocity and core radius of a 3D vortex from PTV data

I have obtained PTV data and would like to calculate various components of the vortex, is there a better code than the one i have used below.Not certain about the accuracy of the code.

% select a CSV file to import
[file, path] = uigetfile('*.csv', 'Select the CSV file to import');
if file == 0
    return; % user canceled file selection
end
%Load the data from the CSV file
data = readmatrix(fullfile(path, file), 'HeaderLines', 1);
% x,y coordinate vectors (mm) and Vx and Vy fields (m/s)
x = data(:, 1);
y = data(:, 2);
z = data(:, 3);
u = data(:, 4);
v = data(:, 5);
w = data(:, 6);
vorticity =data(:,20);
% Example data (replace these with your actual data)
% Example data (replace these with your actual data)
y = linspace(-5, 5, 100) / 1000;    % y-coordinates (converted to meters)
z = linspace(-5, 5, 100) / 1000;    % z-coordinates (converted to meters)
x = linspace(-5, 5, 100) / 1000;    % x-coordinates (converted to meters)
[Y, Z, X] = meshgrid(y, z, x);      % Create meshgrid for y, z, x


% Step 1: Calculate the vorticity magnitude
vorticity_magnitude=max(abs(vorticity(:)));

% Step 2: Define the threshold for vortex core (e.g., 10% of max vorticity)
max_vorticity = max(vorticity_magnitude(:));
threshold = 0.1 * max_vorticity;  % Threshold for vortex core (10% of max vo

% Step 3: Find the vortex core radius (where vorticity exceeds threshold)
% Compute the radial distance from the center (0, 0, 0)
radius = sqrt(X.^2 + Y.^2 + Z.^2);  % Radial distance in 3D

% Identify the region where vorticity exceeds threshold
core_radius = max(radius(vorticity_magnitude >= threshold));  % Vortex core radius

% Step 4: Find the location of maximum vorticity (core of the vortex)
[~, core_idx] = max(vorticity_magnitude(:));  % Find the index of max vorticity
[core_x, core_y, core_z] = ind2sub(size(vorticity_magnitude), core_idx);  % Get the coordinates of the core

% Step 5: Radial and Axial Velocities at the vortex core
% The radial velocity is the component of velocity in the direction of the radial vector
% The axial velocity is the component of velocity along the axis of the vortex (e.g., the x-axis)


% Calculate the radial velocity at the core
r_vec = [X(core_x, core_y, core_z), Y(core_x, core_y, core_z), Z(core_x, core_y, core_z)];  % Position vector at core
r_mag = norm(r_vec);  % Magnitude of position vector

% Calculate the radial velocity at core
radial_velocity = (u(core_x, core_y, core_z) * r_vec(1) + ...
                   v(core_x, core_y, core_z) * r_vec(2) + ...
                   w(core_x, core_y, core_z) * r_vec(3)) / r_mag;  % Dot product of velocity and radial direction

% The axial velocity is the velocity component along the x-axis (assuming the vortex is aligned with the x-axis)
axial_velocity = u(core_x, core_y, core_z);  % Axial velocity along x-axis at the core

% Step 6: Find the velocity at the outer boundary (where radius is maximum)
% Find the index at the outermost boundary (max radial distance)
[~, boundary_idx] = min(abs(radius(:) - max(radius(:))));  % Find the outer boundary index
[boundary_x, boundary_y, boundary_z] = ind2sub(size(radius), boundary_idx);  % Get boundary coordinates

boundary_velocity = sqrt(u(boundary_x, boundary_y, boundary_z)^2 + v(boundary_x, boundary_y, boundary_z)^2 + w(boundary_x, boundary_y, boundary_z)^2);

% Display Results
fprintf('Vortex core radius: %.4f meters\n', core_radius);
fprintf('Radial velocity at the core: %.4f m/s\n', radial_velocity);
fprintf('Axial velocity at the core: %.4f m/s\n', axial_velocity);


%

Answer 1

[Himanshu]

Hi Merina,

I see that you are trying to calculate the axial, radial velocity, and core radius of a 3D vortex from PTV data. You can try the following changes in your MATLAB code:

1. Ensure X, Y, Z are created using the unique values from the imported data to match the actual data dimensions.
2. Use the velocity components to compute the vorticity magnitude accurately as part of the 3D vector field analysis.
3. Ensure the velocity at the outer boundary is calculated correctly by considering the maximum radial distance.

I hope this helps.

[Merina Amon]

Hi Himanshu,

Thanks, may you elaborate thius on the code itself, i have refined the code so many times to no avail.