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

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)]);

Answer 1

[Shishir Reddy]

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.

Answer 2

[Moreno, M.]

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