%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%
%% Copyright (C) 2005-2007 Anurag Singh %%
%% %%
%% This program/code snippet/file (hence forth refered as "library") %%
%% is free software; you can redistribute it and/or %%
%% modify it under the terms of the GNU Lesser General Public %%
%% License as published by the Free Software Foundation; either %%
%% version 2.1 of the License, or (at your option) any later version. %%
%% %%
%% This library is distributed in the hope that it will be useful, %%
%% but WITHOUT ANY WARRANTY; without even the implied warranty of %%
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU %%
%% Lesser General Public License for more details. %%
%% %%
%% You should have received a copy of the GNU Lesser General Public %%
%% License along with this library; if not, write to the Free Software %%
%% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA %%
%% %%
%% FILENAME: %%
%% piv_zplus50.m %%
%% %%
%% AUTHOR: %%
%% %%
%% Anurag Singh, %%
%% MS, 2007 %%
%% Aerospace Engineering & Mechanics %%
%% University of Minnesota - Twin Cities. %%
%% Minneapolis, MN 55455 (USA) %%
%% %%
%% (currently working at: Exa Corporation, Burlington, MA 01803) %%
%% %%
%% CONTACT/EMAIL: %%
%% %%
%% anurag@aem.umn.edu %%
%% anurag9@gmail.com %%
%% %%
%% SOURCE CONTROL INFORMATION: %%
%% None (since I was planning on putting it under source control since it has %%
%% reached the critical file system size. Would be a good thing to put it under %%
%% source control while making changes. %%
%% %%
%% DESCRIPTION: %%
%% %%
%% This matlab script loads data for turbulence statistics computation. %%
%% it uses different high PIV dataset and picks only certain images (usually 1/5th of the %%
%% total available data). %%
%% %%
%% Enter the data for the experiment in the following variables %%
%% function piv %%
%% %%
%% %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
% load global variables (the file "piv_variables.m" is a global variable list scipt which declares all the variables
piv_variables;
%------------------------------------------ INPUT Data Initializations -----------------------
load_stat = 1;
save_stat = 0;
calc_swirl = 0;
showplot = 0;
% Define Conastants
R_gas = 287.05; % Gas Constants
mmHg_Pa = 133.322368421; % Conversion Factor, mm Hg to Pa
g = 9.8; % Acceleration due to gravity, m/s^2
nu = 1.51e-5; % Kinematic Viscosity, m^2/s (computed by interpolating for a given temperature, from the web
rho_hg = 13579.04; % Density of Mercury, kg/m3
%-----------------------
% Enter the values here
%-----------------------
% Data file information
NX = 64;
NY = 64;
num_files = 150; % number/range of files to be read, actual number would depend on skip_files value
num_datset = 1; % total number of data sets
skip_files = 0; % number of files to be skipped while loading them in
%
% Note:
% total files read would be = integer(num_files/(skip_files+1))*num_datset
% when no files are skipped it would be
% Path & File details/ Common part in the file names/ Folder names:
folder_name= 'C:\Documents and Settings\Jazz\My Documents\My Research\12272006_175us_3mm';
%folder_name= 'E:\Anurag\Laptop\Research\12272006_175us_3mm';
file_name = 'B'; % Common part in the file names
data_folder= 'C:\Documents and Settings\Jazz\My Documents\My Research\PIV BL Data';
%data_folder= 'E:\Anurag\Laptop\Research\Matlab workspace Data';
blp_flname = 'pivblp.mat';
pivstat_flname = 'pivstat.mat';
% Tunnel & Lab Atmospheric Data
t_inf = 24.8; % Lab Temperature, Measured, Deg C
p_inf = 746.5; % Lab Pressure, measured, mm of Hg
delta_p = 0.159; % Tunnel free stream temperature
% Preston tube data for u_tau computations, usually taken from the HW data
dP_pres = 0.0353; % Preston tube reading for utau computatissons
diapres = 0.001; % Preston Tube diameter, m
% Field of view and PIV data
dt = 175/1.0e6; % Time delay between the two pulses, seconds
fov = 86; % Field of view, mm
piv_z_pos = 0.003; % Laser sheet location from the Wall, m
ccd_siz = 1024; % Camera CCD size, pixels
int_w = 32; % Interrogation Window dimensions, pixels
capfreq = 500; % Image capture frequency, Hz
% Computations
t_inf = t_inf + 273.15; % Lab Pressure, deg K
p_inf = p_inf*mmHg_Pa; % From pitot tube measurement
rho_air = p_inf/(R_gas*t_inf); % Density of air @ , kg/m^3
piv_u_inf = sqrt(2*delta_p*rho_hg*g*1.0e-3/rho_air); % Free stream velocity, m/s
u_exp = piv_u_inf; % piv_u_inf in experiment
vecimfreq = capfreq/2; % frequency of the vector images (it is half the capture frequency)
%Conversion factors- variable name convention - <from units>_<to units>
pix_m = fov/ccd_siz/1000; % pixels to meters conversion factor
% Load Data
% Load hot-wire boundary layer profile data
load([data_folder,'\', blp_flname]);
% Create the variables beforehand to speed up the processing
% Change the initial filenames here
% Assumption is filenames are serially marked and if it breaks no more files will be read
% Code must be changed for the general case, for my case it works fine.
% Create variables
if(skip_files == 0)
NZ = num_files*num_datset;
else
if(mod(num_files*num_datset,(skip_files +1)) ~= 0)
NZ = floor(num_files*num_datset/(skip_files +1)) + 1;
else
NZ = floor(num_files*num_datset/(skip_files +1));
end
end
u_pix = zeros(NZ,NX,NY);
v_pix = zeros(NZ,NX,NY);
sw2D = zeros(NZ,NX,NY); % 2D Swirl
vort2D = zeros(NZ,NX,NY); % Z component of vorticity
coreid = zeros(NZ,NX,NY); % Vortex core signature array
circ = zeros(NX,NY);
swfrm = zeros(NX,NY);
swfrm1 = zeros(NX,NY);
idfrm = zeros(NX,NY);
pidfrm = zeros(NX,NY);
nidfrm = zeros(NX,NY);
markfrm= zeros(NX,NY);
pointcnt= zeros(1,NX*NY);
cores = zeros(NX,NY);
swcores = zeros(NX,NY);
vortcores = zeros(NX,NY);
core = struct('length',{},'width',{},'area',{},'circ',{},'avgvort',{},'maxsw',{},'maxswx',{},'maxswy',{});
regions = struct('length',{},'width',{},'area',{},'minu',{},'minux',{},'minuy',{});
thresh = struct('value',{},'regions',regions);
dudx = zeros(NX);
dudy = zeros(NX);
dvdx = zeros(NX);
dvdy = zeros(NX);
dfdx = zeros(NX);
dfdy = zeros(NX);
ufrm = zeros(NX,NY);
ufrm1 = zeros(NX,NY);
nuidfrm = zeros(NX,NY);
newregion = zeros(NX,NY);
slowregion = zeros(NX,NY);
%display_status; % Print on the screen some statistics about this data reading
disp('===================================================================');
disp(' OPERATIONAL STATISTICS FOR VARIOUS TASKS');
disp('===================================================================');
disp ([' Input data folder = ',folder_name]);
disp ([' Load data folder = ',data_folder]);
disp ([' Num_files = ',num2str(num_files)]);
disp ([' Num_datset = ',num2str(num_datset)]);
disp ([' skip_files = ',num2str(skip_files)]);
disp ([' Total files to be read = ',num2str(num_datset * floor( (num_files + skip_files)/(skip_files+1)))]);
disp(' ---------------------------------------------------------------');
disp ([' BLP file name = ',blp_flname]);
disp ([' PIV Stat file name = ',pivstat_flname]);
disp (' ');
disp ([' load_stat = ',num2str(load_stat),' (1- load blp data, 0 - no blp load)']);
disp ([' save_stat = ',num2str(save_stat),' (1- save blp data, 0 - no blp save)']);
disp ([' calc_swirl = ',num2str(calc_swirl),' (1- calculate swirl data, 0 - no swirl calculation)']);
disp ([' showplot = ',num2str(showplot),' (1- show plots, 0 - disable show plot)']);
disp('===================================================================');
disp('Reading data files, please wait (will prompt after every 25th file) ...');
disp('');
flincr = skip_files + 1;
k = 1;
for j = 1:num_datset
for i = 1:flincr:num_files
fl_name = [folder_name,'\',file_name,num2str(i,'%.5d'),'.txt'];
fid = fopen(fl_name,'r');
if(fid <0 )
disp(['ERROR Opening file ''',fl_name,''' .... skipping']);
continue;
end
firstline = fgetl(fid);
[DNZ NX NY] = strread(firstline, '%*s%*s%*s %f %f %f %*c%*s%*c%*s%*s%*s', 'delimiter', ' ');
temp = fscanf(fid,'%g %g %g %g',[4,inf]);
if i == 1
x_pix = reshape(temp(1,:), [NX NY]);
y_pix = reshape(temp(2,:), [NX NY]);
end
u_pix(k,:,:) = reshape(temp(3,:), [NX NY]);
v_pix(k,:,:) = reshape(temp(4,:), [NX NY]);
fclose(fid);
if rem(k,25) == 0 % Display a message on screen for every 25th file read
disp(['Done past ',file_name,num2str(i,'%.5d'),'.txt',', Data Set=',num2str(j), ' ; Total files read so far = ',num2str(k),]);
end
k = k+1;
end
end
disp('--------------- ');
disp(['Total files read = ',num2str(k-1),' from ',num2str(num_datset),' datasets while skipping ',num2str(skip_files), ' files in each set.']);
disp('--------------- ');
%--------------------------------------------------------------------------------------
%---------------- Process Data
%--------------------------------------------------------------------------------------
u = u_pix.*pix_m./dt;
v = v_pix.*pix_m./dt;
x = x_pix.*pix_m;
y = y_pix.*pix_m;
%u(:,end-1:end,:) = [];
%u(:,:,end-1:end) = [];
%u(:,1:4,:) = [];
%u(:,:,1:4) = [];
%
%v(:,end-1:end,:) = [];
%v(:,:,end-1:end) = [];
%v(:,1:4,:) = [];
%v(:,:,1:4) = [];
%
%
%% shift the scales
%x = x - x(4,1);
%y = y - y(4,1);
%
%x(end-1:end,:) = [];
%x(:,end-1:end) = [];
%x(1:4,:) = [];
%x(:,1:4) = [];
%
%y(end-1:end,:) = [];
%y(:,end-1:end) = [];
%y(1:4,:) = [];
%y(:,1:4) = [];
%
% calculate the mean and the rms
if(load_stat == 0)
piv_ubar = mean(mean(mean(u, 1)));
piv_ustd = mean(mean(std(u, 1)));
piv_vbar = mean(mean(mean(v, 1)));
piv_vstd = mean(mean(std(v, 1)));
piv_uvstd = mean(mean(std(u.*v,1)));
if(save_stat == 1)
save([data_folder,'\',pivstat_flname], 'piv_ubar', 'piv_vbar', 'piv_ustd', 'piv_vstd', ...
'piv_ubarplus', 'piv_vbarplus', 'piv_ustdplus', 'piv_vstdplus');
end
else
load([data_folder,'\',pivstat_flname]);
end
u_norm = u - piv_ubar;
% Compute u_tau only if its not available already from Boundary layer profile measurement
if(size(who('u_tau'))==0)
%% Preston tube and other calculations
dp = rho_hg*g*0.001*dP_pres; %
% Patel(JFM) Method
x_star = log10(dp*diapres^2/(4*rho_air*nu*nu));
%y_star = x_star/2 + 0.037;
y_star = 0.8287 - 0.1381*x_star + 0.1437*x_star^2 - 0.0060*x_star^3;
%y_star = solve('y+2*log(1.95*y + 4.10) - 9.3286');
tau_w = 4*rho_air*nu^2*10^(y_star)/diapres^2;
u_tau = sqrt(tau_w/rho_air);
end
% Parameters in wall units
piv_zplus = piv_z_pos*u_tau/nu;
piv_ubarplus = piv_ubar/u_tau;
piv_ustdplus = piv_ustd/u_tau;
piv_vbarplus = piv_vbar/u_tau;
piv_vstdplus = piv_vstd/u_tau;
tau_w = u_tau*u_tau*rho_air;
u_normplus = u_norm/u_tau;
v_plus = v/u_tau;
lenconvfac = (x(2,1) - x(1,1))*u_tau/nu;
cellarea = (x(2,1) - x(1,1))^2;
areaconvfac = (x(2,1) - x(1,1))^2 * (u_tau/nu)^2;
%save pivdata.mat
if(calc_swirl == 1)
calcsw(1,k-1);
end
%function display_status
%global foldername data_folder num_files num_datset skipfiles blp_flname piv_statname load_stat calc_swirl save_stat showplot
%return;
return;
