Add files via upload

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

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'])
+

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
+
+