From ee5f31fb57eb81f6ab0e9127104b2cfb44f73918 Mon Sep 17 00:00:00 2001 From: "Eshraghi, Javad" Date: Tue, 16 Feb 2021 13:54:13 -0500 Subject: [PATCH] Add files via upload --- Step1.m | 138 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Step2.m | 106 +++++++++++++++++++++++++++++++++++++++++++ Step3.m | 107 +++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 351 insertions(+) create mode 100644 Step1.m create mode 100644 Step2.m create mode 100644 Step3.m diff --git a/Step1.m b/Step1.m new file mode 100644 index 0000000..d83920e --- /dev/null +++ b/Step1.m @@ -0,0 +1,138 @@ +close all; +clc; +clear; + +%% Directory +Sphere_Speed = '20'; % Impact velocity +Diameter = '0750'; % Sphere Diameter +Material = '4'; % density ratio +Trial = 1; % trial # +Sphere_Diameter_Cavity_Im = 85; % Sphere diameter in Pixel +fps = 5000; % the rate that the camera was filming at +Thresh = 50000; + +cFrame = 10; % the frame where contact with the surface occurs +sFrame = 15; % the frame where the sphere completely submerged +pFrame = 334; % the frame that pinch off occurs at +derivative_order = 4; % The order of numerical derivative to find air flow rate + +imageName_contact = sprintf('S%s_D%s_M%s_%02d%05d.tif',Sphere_Speed,Diameter,Material,Trial,cFrame); +imageName_pinchoff = sprintf('S%s_D%s_M%s_%02d%05d.tif',Sphere_Speed,Diameter,Material,Trial,pFrame); + +folderName = ['S' num2str(Sphere_Speed) '_D' num2str(Diameter) '_M' num2str(Material) '_0' num2str(Trial)]; + +% cd('F:\..\Sphere_experiment\Data_set_images\Cavity_camera') +% cd(num2str(folderName)) + + +%% Body +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% User selects the following points on the submerged (one point) and +% pinch-off (three points): +% i. free surface +% ii. pinch-off location +% iii. bottom of sphere +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +figure(1) +Im1 = imread(imageName_contact); +imshow(Im1) +[x1, y1] = ginput(1); + +figure(2) +Im2 = imread(imageName_pinchoff); +imshow(Im2) +[x2, y2] = ginput(2); + +freesurface_Y = y1; +sphere_center_Y = (y2(2)-Sphere_Diameter_Cavity_Im/2); +pinchoff_Y = y2(1); + + +%% Initializations +calib = Sphere_Diameter_Cavity_Im/(0.75*2.54); %pix/cm +R = Sphere_Diameter_Cavity_Im/2; % radius, pi +rho_air = 1.205; +Cavity_volume = zeros; + +%% Main Time Loop, Frame by Frame increments +for ii = sFrame:pFrame + Im_Name = sprintf('S%s_D%s_M%s_%02d%05d.tif',Sphere_Speed,Diameter,Material,Trial,ii); + Im = imread(Im_Name); % reads the image + Im_Size_orig = size(Im); + %% Image Processing + BWcrop = imcrop(Im,[1 y1 Im_Size_orig(2) Im_Size_orig(1)]); % crops the image + %%Volume calculation + left = zeros; + right = zeros; + Cavity_radius = zeros; + Im_Size = size(BWcrop); + + for jj = 1:Im_Size(1) + for kk = 1:Im_Size(2) + if BWcrop(jj,kk)>Thresh + left (jj) = kk; + break + end + end + for kk = 1:Im_Size(2) + if BWcrop(jj,Im_Size(2)-kk+1)>Thresh + right (jj) = Im_Size(2)-kk; + break + end + end + end + % Find velocity of the spheres + y_max (ii) = max(length(left),length(right)); + + imshow(BWcrop) + hold on + plot(left,'LineWidth',2) + plot(right,'LineWidth',2) + hold off + pause(0.005) + Cavity_radius = (right-left)/2; + Cavity_volume (ii) = sum(pi.*Cavity_radius.^2)-(4/3*pi*R^3); +end +Vol_Smoothed = smooth(0:1:(pFrame-1),Cavity_volume(1:pFrame),0.35,'rloess'); % Smoothing +Vol_Smoothed_P = Vol_Smoothed; +% if Vol_Smoothed(1)<0 + Vol_Smoothed = Vol_Smoothed - min(Vol_Smoothed); +% end +Volume_Dimensionless = Vol_Smoothed/(4/3*pi*R^3); + +% 2nd order derivative +if (derivative_order==2) + for jj = 2:(pFrame-1) + Vol_Flow_Rate (jj) = 1/2 * (Vol_Smoothed(jj+1) - Vol_Smoothed(jj-1)); % cubic px per frame + end +end + +% 4th order derivative +if (derivative_order==4) + for jj = 2:(pFrame-3) + Vol_Flow_Rate (jj-1) = 1/12 * (-3 * Vol_Smoothed(jj-1) - 10 * Vol_Smoothed(jj) + 18 * Vol_Smoothed(jj+1) - 6 * Vol_Smoothed(jj+2) + Vol_Smoothed(jj+3)); + end +end + +Vol_Flow_Rate_size=size(Vol_Flow_Rate); +Vol_flow_rate_Smoothed = smooth(1:1:(Vol_Flow_Rate_size(2)),Vol_Flow_Rate(1:Vol_Flow_Rate_size(2)),0.35,'rloess'); % Smoothing +Volume_flow_rate = Vol_flow_rate_Smoothed; % cubic px per frame +Air_velocity = Volume_flow_rate/(pi*R^2); % px per frame +Air_velocity_Dimensionless = Air_velocity / ((str2double(Sphere_Speed)*calib/fps*10)); +Pressure_diff = rho_air * 0.5 * Air_velocity.^2; % kg/(px.frame squared) +Pressure_diff_Dimensionless = Pressure_diff / (rho_air * 0.5 * (str2double(Sphere_Speed)*calib/fps*10).^2); %Dimensionless +Sphere_Vel = diff (y_max); % px per frame +Sphere_Pinchoff_Vel = mean(Sphere_Vel(end-5:end)); + +figure,plot(Volume_Dimensionless) +figure,plot(Cavity_volume) +hold on +plot(Vol_Smoothed) +hold off +figure,plot(Air_velocity_Dimensionless) +figure,plot(Pressure_diff_Dimensionless) +% Save Functions +cd('F:\Processing\Sphere_experiment\Data_set_processing\Cavity_volume_info') +% save([folderName '_CavityVolumeinfo' '.mat'],'Sphere_Pinchoff_Vel','Sphere_Vel','Pressure_diff_Dimensionless','Cavity_volume','Volume_flow_rate','Air_velocity') +save([folderName '_Cavityinfo' '.mat']) diff --git a/Step2.m b/Step2.m new file mode 100644 index 0000000..f082754 --- /dev/null +++ b/Step2.m @@ -0,0 +1,106 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% either bofore or after of this step we need to determine the sphere radius +% and x position of the center of the sphere to be able to dimensionalize +% and find the mapping function to compare the results +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +close all; +clc; +clear; + +%% Directory +Sphere_Speed = '35'; % Impact velocity +Diameter = '0750'; % Sphere Diameter +Material = '2'; % the height it was dropped at +Trial = 3; % the number of the trial to be examined +Sphere_Diameter_Splash_Im = 180; % Sphere diameter in Pixel + +ssFrame = 22; % the frame where the sphere completely appeared +ccFrame = 88; % the frame that the sphere contacted the surface +closureFrame = 265; % the frame that the splash is closed + +imageName_submerged = sprintf('S%s_D%s_M%s_%02d%05d.tif',Sphere_Speed,Diameter,Material,Trial,ssFrame); +imageName_contact = sprintf('S%s_D%s_M%s_%02d%05d.tif',Sphere_Speed,Diameter,Material,Trial,ccFrame); +folderName = ['S' num2str(Sphere_Speed) '_D' num2str(Diameter) '_M' num2str(Material) '_0' num2str(Trial)]; +cd(num2str(folderName)) + +figure(1) +imshow(imageName_submerged) +[x_right_side, y_right_side] = ginput(1); +Sphere_right_side = x_right_side; + +figure(2) +imshow(imageName_contact) +[x_level, y_level] = ginput(1); + + + +Sphere_Radius = Sphere_Diameter_Splash_Im/2; +Sphere_Center_X = Sphere_right_side - Sphere_Radius; +Im_Size = size(imageName_submerged); +Image_height = Im_Size(2); +Free_Surface_Y = y_level; +set_height = 1; + +% cd('F:\..\Sphere_experiment\Data_set_images\Splash_camera') + + +%% +calib = Sphere_Diameter_Splash_Im/(str2double(Diameter)*2.54/1000); %pix/cm +rho_air = 1.205; + +% %% Initializations +% Cavity_volume = zeros; + +%% Select the surface height +% figure(3) +% Im1 = imread(imageName_contact); +% imshow(Im1) +% [x1, y1] = ginput(1); + +% %% Main Time Loop, Frame by Frame increments +% for ii = ssFrame:ccFrame +% +% Im_Name = sprintf('S%s_D%s_M%s_%02d%06d.tif',Sphere_Speed,Diameter,Material,Trial,ii); +% Im = imread(Im_Name); % reads the image +% Im_Size_orig = size(Im); +% +% %% Image Processing +% BWcrop = imcrop(Im,[1 1 Im_Size_orig(2) y1]); % crops the image +% BW = im2bw(BWcrop,.7); % converts image to black and white +% % only pixels that are connected to another # pixels are kept +% BWred = bwareaopen(~BW, 100); +% BWfill=imfill(BWred,'holes'); % filles in holes of the image +% B = bwboundaries(BWfill); % finds the boundaries of the image +% sizeB = size(B); % size of the boundary matrix +% +% %%%% Sometimes, B will have empty cell array (sizeB = [0 1]) +% if sizeB ~= 0 +% bound = B{1}; +% else +% bound = [0 0]; +% end +% +% % Stores the boundary in variables +% boundY = bound(:,1)'; +% boundX = bound(:,2)'; +% +% % Find velocity of the spheres +% S_y_max (ii) = max(boundY); +% end +% +% Sphere_Impact_Vel = diff (S_y_max); % px per frame +% Sphere_Impact_Velocity = Sphere_Impact_Vel(end); + + +%% Find Trajectory and initial velocity + +% rr_right_dimensionless = (rim_x_right - Sphere_Center_X) / Sphere_Radius; +rr_left_dimensionless = -(rim_x_left - Sphere_Center_X) / Sphere_Radius; +% zz_right_dimensionless = ((Free_Surface_Y) - rim_y_right) / Sphere_Radius ; +zz_left_dimensionless = ((Free_Surface_Y) - rim_y_left) / Sphere_Radius ; + + +cd ('F:\Processing\Sphere_experiment\Data_set_processing\Experimental_trajectory_info') +save([folderName '_ExpTrajectory' '.mat']) + diff --git a/Step3.m b/Step3.m new file mode 100644 index 0000000..20d7b94 --- /dev/null +++ b/Step3.m @@ -0,0 +1,107 @@ +clear all; clc; close all; +cd('F:\Processing\Sphere_experiment\Data_set_processing\Experimental_trajectory_info\375') +%% Directory +Sphere_Speed = '20'; % Impact velocity +Diameter = '0375'; % Sphere Diameter +Material = '1'; % the height it was dropped at +Trial = 2; % the number of the trial to be examined + +folderName = ['S' num2str(Sphere_Speed) '_D' num2str(Diameter) '_M' num2str(Material) '_0' num2str(Trial)]; +load([folderName '_ExpTrajectory' '.mat']) +%% Initialization +% Splash Parameters +R = 0.00238125; +a = 0.14 * R; +Cd = 0.001; +deltaP = 150; +left = 1; +right = 0; +save_model_trajectory_right = 0; +save_model_trajectory_left = 0; + +save_modified_trajectory = 0; +% Initial conditions +u_0 = 1.4; +v_0 = 1.7; + +% Show comparing plots of two models (with and without 2nd surface tension) +show_comp_plot = 0; + +% fluid properties and Constants +g = 9.81; +sigma = 0.0728; +rho = 998.21; +nu_air = 15.11e-6; +rhoair = 1.205; + +%% Considering Pressure with constant Cd and delta P +f1 = @(t,x) [x(3);x(4);-(g*(x(3)*x(4))/(2*(x(3)^2+x(4)^2))+2*sigma*x(3)/(rho*pi*a^2*sqrt(x(3)^2+x(4)^2))+Cd*x(3)*sqrt(x(3)^2+x(4)^2)/(pi*a)+(deltaP*x(4))/(rho*pi*a*sqrt(x(3)^2+x(4)^2)));-(g*(x(3)^2+2*x(4)^2)/(2*(x(3)^2+x(4)^2))+2*sigma*x(4)/(rho*pi*a^2*sqrt(x(3)^2+x(4)^2))+Cd*x(4)*sqrt(x(3)^2+x(4)^2)/(pi*a)-(deltaP*x(3))/(rho*pi*a*sqrt(x(3)^2+x(4)^2)))]; +[t1,xa1] = ode45(f1,[0:10e-5:0.0086],[R 0 u_0 v_0]); +r1 = xa1(:,1)/R; +z1 = xa1(:,2)/R; + +f2 = @(t,x) [x(3);x(4);-(g*(x(3)*x(4))/(2*(x(3)^2+x(4)^2))+sigma*x(4)*(sqrt(x(3)^2+x(4)^2))*(2*x(3)^2+x(4)^2)/(2*a*x(1)*rho*pi*(x(3)^2+x(4)^2)^2)+2*sigma*x(3)/(rho*pi*a^2*sqrt(x(3)^2+x(4)^2))+Cd*sqrt(x(3)^2+2*x(4)^2)*x(3)/(pi*a)+(deltaP*x(4))/(rho*pi*a*sqrt(x(3)^2+x(4)^2)));-(g*(x(3)^2+2*x(4)^2)/(2*(x(3)^2+x(4)^2))+sigma*x(4)*(sqrt(x(3)^2+x(4)^2))*(x(3)*x(4))/(2*a*x(1)*rho*pi*(x(3)^2+x(4)^2)^2)+2*sigma*x(4)/(rho*pi*a^2*sqrt(x(3)^2+x(4)^2))+Cd*sqrt(x(3)^2+2*x(4)^2)*x(4)/(pi*a)-(deltaP*x(3))/(rho*pi*a*sqrt(x(3)^2+x(4)^2)))]; +[t2,xa2] = ode45(f2,[0:10e-5:0.0086],[R 0 u_0 v_0]); +r2 = xa2(:,1)/R; +z2 = xa2(:,2)/R; + +%% Plotting +if show_comp_plot == 1 + figure('units','normalized','outerposition',[0 0 0.75 0.75]) + set(gca,'fontsize',14,'FontName','Garamond','FontWeight','bold','Color','w'); + set(gcf,'color','white'); + hold on + plot(r1, z1, 'co-') + hold on + plot(r2, z2, 'ro-') + xlabel('r/R','fontsize',16,'FontName','Garamond','FontWeight','bold') + ylabel('z/R','fontsize',16,'FontName','Garamond','FontWeight','bold') + legend('Neglecting 2nd surface tension','Considering 2nd surface tension','location','southeast') + legend boxoff + hold off +end + +if left == 1 + figure('units','normalized','outerposition',[0 0 0.75 0.75]) + set(gca,'fontsize',14,'FontName','Garamond','FontWeight','bold','Color','w'); + set(gcf,'color','white'); + hold on + plot(r2, z2, 'b-') + hold on + plot(rr_left_dimensionless,zz_left_dimensionless, 'rX') + xlabel('r/R','fontsize',16,'FontName','Garamond','FontWeight','bold') + ylabel('z/R','fontsize',16,'FontName','Garamond','FontWeight','bold') + legend('Model Prediction','Experiment','location','southeast') + legend boxoff + hold off +end +if right == 1 + figure('units','normalized','outerposition',[0 0 0.75 0.75]) + set(gca,'fontsize',14,'FontName','Garamond','FontWeight','bold','Color','w'); + set(gcf,'color','white'); + hold on + plot(r2, z2, 'b-') + hold on + plot(rr_right_dimensionless,zz_right_dimensionless, 'rX') + xlabel('r/R','fontsize',16,'FontName','Garamond','FontWeight','bold') + ylabel('z/R','fontsize',16,'FontName','Garamond','FontWeight','bold') + legend('Model Prediction','Experiment','location','southeast') + legend boxoff + hold off +end + +% Saving + +if save_model_trajectory_right == 1 + % Save Info + cd('F:\Processing\Sphere_experiment\Data_set_processing\Model_trajectory_info') + save([folderName '_RightModelTrajectoryinfo' '.mat']) +end + +if save_model_trajectory_left == 1 + % Save Info + cd('F:\Processing\Sphere_experiment\Data_set_processing\Model_trajectory_info') + save([folderName '_LeftModelTrajectoryinfo' '.mat']) +end + +