Need to extract boundary layer thickness from 2D planar PIV data.

17 次查看（过去 30 天）
Miguel Olvera
Miguel Olvera 2024-7-5
回答： Moreno, M. 2024-8-12
This is 2D Planar PIV data at a freestream velocity of ~ 13 m/s.
My code currently selects a velocity profile where max(profile_u) = 12.41 m/s. I want to use a velocity profile where max(profile_u) is closer to 13 m/s.
Scatter(u,y) will show all the velocity profiles I am working with.
clc;
close all;
clear all;
%% BASELINE
% Load the data
opts = detectImportOptions('FlatPlate_Baseline_avg0001(in).csv');
opts.VariableNamingRule = 'preserve';
data = readtable('FlatPlate_Baseline_avg0001(in).csv', opts);
% Extract variables
x = data.x / 1000; % Convert to meters
y = data.y / 1000; % Convert to meters
u = data.u;
v = data.v;
%%
% Filter data for x values between 0.05m to 0.08m
filter_indices = (x >= 0.05) & (x <= 0.08);
x_filtered = x(filter_indices);
y_filtered = y(filter_indices);
u_filtered = u(filter_indices);
v_filtered = v(filter_indices);
% Find unique x and y values within the filtered range
unique_x = unique(x_filtered);
unique_y = unique(y_filtered);
nu = 1.5e-5;
kappa = 0.41;
B = 5.0;
profile_spacing = 2 / 1000;
profile_x_positions = unique_x(unique_x >= 0.05);
wall_shear_stress = zeros(size(profile_x_positions));
u_tau_values = zeros(size(profile_x_positions));
delta_99_values = zeros(size(profile_x_positions));
delta_star_values = zeros(size(profile_x_positions)); % Displacement thickness
theta_values = zeros(size(profile_x_positions)); % Momentum thickness
for i = 1:length(profile_x_positions)
profile_indices = find(abs(x - profile_x_positions(i)) < 1e-6);
profile_y = y(profile_indices);
profile_u = u(profile_indices);
valid_indices = profile_y > 0 & profile_y * max(profile_u) / nu > 1e-6;
profile_y = profile_y(valid_indices);
profile_u = profile_u(valid_indices);
if isempty(profile_y) || isempty(profile_u)
fprintf('No valid data points at x = %f meters. Skipping...\n', profile_x_positions(i));
continue;
end
% Estimate u_inf (freestream velocity) as the max value of the profile
u_inf = max(profile_u);
% f(eta) = u / u_inf
f_eta = profile_u / u_inf;
[f_eta_unique, unique_indices] = unique(f_eta);
profile_y_unique = profile_y(unique_indices);
% Find y_99 where f(eta) = 0.99
y_99 = interp1(f_eta_unique, profile_y_unique, 0.99, 'linear', 'extrap');
eta_99 = 0.99;
% delta_99
delta_99 = y_99 / eta_99;
delta_99_values(i) = delta_99;
% displacement thickness (delta_star)
delta_star = abs(trapz(profile_y, (1 - f_eta)));
delta_star_values(i) = delta_star;
% momentum thickness (theta)
theta = abs(trapz(profile_y, f_eta .* (1 - f_eta)));
theta_values(i) = theta;
Re_theta = (u_inf * theta)/nu;
% Calculate wall shear stress and friction velocity
fun_u_tau = @(u_tau) mean(profile_u - (u_tau / kappa) * log(profile_y * u_tau / nu + 1) - B * u_tau);
try
options = optimset('Display', 'off');
u_tau = fzero(fun_u_tau, [0.01, 10], options);
u_tau_values(i) = u_tau;
catch err
fprintf('Cannot converge at x = %f meters: %s\n', profile_x_positions(i), err.message);
continue;
end
rho = 1.225;
tau_wall = rho * u_tau^2;
wall_shear_stress(i) = tau_wall;
% % Plot y/delta vs u/u_inf
% figure;
% y_over_delta = profile_y / delta_99;
% u_over_u_inf = profile_u / u_inf;
% plot(u_over_u_inf, y_over_delta, '-o');
% title(sprintf('Velocity Profile at x = %.3f meters', profile_x_positions(i)));
% xlabel('u / u_{\infty}');
% ylabel('y / \delta');
% grid on;
end
% Display the results
results_table = table(profile_x_positions, delta_99_values, delta_star_values, theta_values);
disp(results_table);
disp(['Momentum Thickness Reynolds Number (Re_theta): ', num2str(Re_theta)]);

请先登录，再回答此问题。

回答（2 个）

Shishir Reddy
Shishir Reddy 2024-7-17
Hi Miguel,
As per my understanding you would like to use a velocity profile where max(profile_u) is closer to 13 m/s for calculating boundary layer thickness.
Upon investigation, I see that the maximum velocity i.e., max(profile_u) does not match the expected freestream velocity of 13m/s because of noise or outliers in the provided data which could be the affecting factor for the maximum value.
This can be solved by normalizing the velocity values, such that the scaled value of max(profile_u) should be 13 m/s.
You can add the following line of code in your script to achieve scaling of the profile_u data.
profile_u = profile_u*(13/max(profile_u));
I hope this helps.
Moreno, M.
Moreno, M. 2024-8-12
Dear Miguel,
What boundary layer thickness do you wish to calculate (displacement, momentum, energy, or boundary layer edge)? If the boundary layer is 'laminar', a typical 99% freestream velocity approach should work to find the boundary layer edge. If you have upstream distortions, you may need to do a little extra work to define the non-integral properties with a degree of robustness.
If you were interested, please find my script for calculating boundary layer thicknesses systematically here: https://uk.mathworks.com/matlabcentral/fileexchange/171234-boundary-layer-thickness-calculator
It works well with many types of flows, even with SBLIs and distortions.
If your profile was very laminar, you may get away with this:
% Find index for DELTA99 based on 'U'-velocity
idx = find(u-0.99*u(end)<0,1,'last')+1;
I hope this helps,
Miguel

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by