Force_Calculation_Capillary_Hydrostatic_Weight.m (4410 bytes)

clear

clc

format compact



%% Force Calculation Information:

% Capillary force:

% Fc = pi*dp*gamma*(1-cos(theta))/2

    % Dp is particle diameter

    % gamma is fluid surface tension

    % theta is critical contact angle

% Hydrostatic force:

% Fhs = 0.25*pi*Dp^2*(1-cos(theta))*(2*gamma/Db - rhow*g*Db/2)

    % Dp, gamma, theta - see above

    % Db is bubble diameter

    % rhow is water density

% Net weight force:

% Fw = 1/6*pi*Dp^3*rhop*g - 1/8*pi*Dp^3*rhow*g*[2/3 + cos(theta/2)-1/3*cos(theta/2)^3]

    % Dp, rhow, theta - see above

    % rhop is particle density

% Force due to local turbulence:

% Ft = 4/3*pi*Rp^3*delrho*bm

% bm = 1.9*epsilon^(2/3)/(r^1/3)

    % Rp is particle radius

    % delrho is density difference between particle and water

    % bm is acceleration due to local turbulence

    % epsilon is fluid energy dissipation per unit volume

    % r is vortex radius, taken as bubble diameter

% Drag force:

% Fd = 3*pi*Dp*eta*u

    % Dp - see above

    % eta is dynamic viscosity of fluid

    % u is particle velocity

% Hydrodynamic force:

% Fhd = 6*pi*eta*Rp*V*f1

% f1 = [1+(Rp/(4*H))^p]^(1/p)

    % eta, Rp - see above

    % V is particle0bubble relative velocity

    % f1 is modification factor

    % H is particle-bubble separation

    % p is empirical factor



%% To get Capillary force:

syms Dp

gamma = 0.073;  %unit: N/m

theta = 91; %contact angle of a flat glass surface, unit: degree



Fc = pi*Dp*gamma*(1-cos(theta))/2;



% figure(1)    %currently useful

% fplot(1000*Dp, Fc, [0 0.002])

% xlabel('Sand Particle Diameter,D_p (mm)')

% ylabel('Capillary Force,F_c (N)')

% set(gca,'box','off');

% set(gca,'TickDir','out');



%% To get Hydrostatic force:

syms Db_1

rhow = 1000;  %unit: kg/m3

g = 9.8;  %unit: m/s2



Db_1 = 0.0005; %unit: m

Db_2 = 0.001;

Db_3 = 0.002;

Db_4 = 0.0025;

Db_5 = 0.005;



Fhs_1 = 0.25*pi*Dp^2*(1-cos(theta))*(2*gamma/Db_1 - rhow*g*Db_1/2);

Fhs_2 = 0.25*pi*Dp^2*(1-cos(theta))*(2*gamma/Db_2 - rhow*g*Db_2/2);

Fhs_3 = 0.25*pi*Dp^2*(1-cos(theta))*(2*gamma/Db_3 - rhow*g*Db_3/2);

Fhs_4 = 0.25*pi*Dp^2*(1-cos(theta))*(2*gamma/Db_4 - rhow*g*Db_4/2);

Fhs_5 = 0.25*pi*Dp^2*(1-cos(theta))*(2*gamma/Db_5 - rhow*g*Db_5/2);



% figure(2)   %currently useful

% fplot(1000*Dp, Fhs_1, [0  0.002])

% hold on

% fplot(1000*Dp, Fhs_2, [0  0.002])

% fplot(1000*Dp, Fhs_3, [0  0.002])

% fplot(1000*Dp, Fhs_4, [0  0.002])

% fplot(1000*Dp, Fhs_5, [0  0.002])

% xlabel('Sand Particle Diameter,D_p (mm)')

% ylabel('Net Hydrostatic Force,F_h_s (N)')

% legend('D_b=0.5 mm','D_b=1 mm','D_b=2 mm','D_b=2.5 mm','D_b=5 mm')

% set(gca,'box','off');

% set(gca,'TickDir','out');

% legend('boxoff');



%% To get Net weight force:

rhop = 2600;  %unit: kg/m3

Fw = 1/6*pi*Dp^3*rhop*g - 1/8*pi*Dp^3*rhow*g*[2/3 + cos(theta/2)-1/3*cos(theta/2)^3];



% figure(3)   %currently useful

% fplot(1000*Dp, Fw, [0 0.002])

% xlabel('Sand Particle Diameter,D_p (mm)')

% ylabel('Net Weight Force,F_w (N)')

% set(gca,'box','off');

% set(gca,'TickDir','out');



%% Combine all static forces in one figure:



figure(4)

fplot(1000*Dp, 1000*Fc, [0 0.002],'k','LineWidth',2)

hold on

fplot(1000*Dp, 1000*Fhs_1, [0  0.002],'--k','LineWidth',2)

%%%%%fplot(1000*Dp, 1000*Fhs_2, [0  0.002],'--k','LineWidth',1)

fplot(1000*Dp, 1000*Fhs_3, [0  0.002],'--k','LineWidth',3)

%%%%%fplot(1000*Dp, 1000*Fhs_4, [0  0.002],'--b')

fplot(1000*Dp, 1000*Fhs_5, [0  0.002],'--k','LineWidth',4)

fplot(1000*Dp, 1000*Fw, [0 0.002],':k','LineWidth',4)

% % % ax=gca;

% % % ax.YAxis.Exponent = 0;



l=legend(' F_c',' F_e_s, with D_b=0.5 mm', ...

    ' F_e_s, with D_b=2 mm', ...

    ' F_e_s, with D_b=5 mm', ...

    ' F_w');

xlabel('Sand Particle Diameter, D_p (mm)')

ylabel('Static Forces (mN)')



ylim([0   1])

set(gca,'box','off');

set(gca,'TickDir','out');

set(gca,'FontName','calibri','FontSize',23)

set(gca,'LineWidth',2)

legend('boxoff');

set(l,'FontName','calibri','FontSize',23);



%% Reynols number:



% % syms n

% % muw = 0.001; % viscosity of water

% % dimp = 0.05; % diameter of impeller, unit m

% % Re = n*dimp^2*rhow/muw;

% % 

% % figure(5)

% % fplot(n,Re, [0 8])

% % yline(10000,'--k')

% % 

% % ylim([0 20000])

% % set(gca,'box','off');

% % set(gca,'TickDir','out');

% % xlabel('Impeller Velocity, v_i_m_p (m/s)','FontSize',14)

% % ylabel('Impeller Reynolds Number, Re_i_m_p (-)','FontSize',14)