Force_Calculation_HydrodynamicAngle.m (7970 bytes)

clear

clc

format compact



%% Case 1: Particle is at the top of bubble:

gamma = 0.073;  %fluiunit: N/m

rhop = 2600; 

rhow = 1000;

g = 9.8;

muw = 0.001;

theta_f = 120; %unit: degree

theta_m = 112;

theta_c = 109;

Db = 0.0025; % A bubble size that can be observed for all types of sand





Dp_f1 = 0.0005;  %unit: m

Dp_m1 = 0.00125;

Dp_c1 = 0.0025;



up_1 = 1;  %unit: m/s

up_2 = 4;

up_3 = 7;

up_4 = 10;



Fc_f1 = 0.5*pi*Dp_f1*gamma*(1-cos(theta_f));

Fc_m1 = 0.5*pi*Dp_m1*gamma*(1-cos(theta_m));

Fc_c1 = 0.5*pi*Dp_c1*gamma*(1-cos(theta_c));





Fw_f1 = 1/6*pi*Dp_f1^3*rhop*g - 1/8*pi*Dp_f1^3*rhow*g*(2/3+cos(theta_f/2)-1/3*(cos(theta_f/2))^3);

Fw_m1 = 1/6*pi*Dp_m1^3*rhop*g - 1/8*pi*Dp_m1^3*rhow*g*(2/3+cos(theta_m/2)-1/3*(cos(theta_m/2))^3);

Fw_c1 = 1/6*pi*Dp_c1^3*rhop*g - 1/8*pi*Dp_c1^3*rhow*g*(2/3+cos(theta_c/2)-1/3*(cos(theta_c/2))^3);





Fhs_f1 = 0.25*pi*Dp_f1^2*(1-cos(theta_f))*(2*gamma/Db-rhow*g*Db/2);

Fhs_m1 = 0.25*pi*Dp_m1^2*(1-cos(theta_m))*(2*gamma/Db-rhow*g*Db/2);

Fhs_c1 = 0.25*pi*Dp_c1^2*(1-cos(theta_c))*(2*gamma/Db-rhow*g*Db/2);





% % % Fd_f11 = 3*pi*Dp_f1*muw*up_1;

% % % Fd_f12 = 3*pi*Dp_f1*muw*up_2;

% % % Fd_f13 = 3*pi*Dp_f1*muw*up_3;

% % % Fd_f14 = 3*pi*Dp_f1*muw*up_4;

% % % 

% % % Fd_m11 = 3*pi*Dp_m1*muw*up_1;

% % % Fd_m12 = 3*pi*Dp_m1*muw*up_2;

% % % Fd_m13 = 3*pi*Dp_m1*muw*up_3;

% % % Fd_m14 = 3*pi*Dp_m1*muw*up_4;

% % % 

% % % Fd_c11 = 3*pi*Dp_c1*muw*up_1;

% % % Fd_c12 = 3*pi*Dp_c1*muw*up_2;

% % % Fd_c13 = 3*pi*Dp_c1*muw*up_3;

% % % Fd_c14 = 3*pi*Dp_c1*muw*up_4;



u_slip = (2/9)*(rhop-rhow)*g*(0.5*Dp_f1)^2/muw;

f1factor = (1+(0.5*Dp_f1/(4*0.1*Dp_f1))^0.719)^(-0.719);

f2pfactor = (1.707+(0.1*Dp_f1)/(0.5*Dp_f1))/(0.836+(0.1*Dp_f1)/(0.5*Dp_f1));

f2ppfactor = (2.656+(0.1*Dp_f1)/(0.5*Dp_f1))/(1.440+(0.1*Dp_f1)/(0.5*Dp_f1));

hfactor = (0.1*Dp_f1+0.5*Dp_f1)/(13*0.5*Dp_f1);

f2factor = (f2pfactor-hfactor*f2ppfactor)/(1-hfactor);

Fhd_f11 =  6*pi*muw*(Dp_f1/2)*(u_slip+up_1)*f1factor + 6*pi*muw*(Dp_f1/2)*up_1*f2factor;

Fhd_f12 =  6*pi*muw*(Dp_f1/2)*(u_slip+up_2)*f1factor + 6*pi*muw*(Dp_f1/2)*up_1*f2factor;

Fhd_f13 =  6*pi*muw*(Dp_f1/2)*(u_slip+up_3)*f1factor + 6*pi*muw*(Dp_f1/2)*up_1*f2factor;

Fhd_f14 =  6*pi*muw*(Dp_f1/2)*(u_slip+up_4)*f1factor + 6*pi*muw*(Dp_f1/2)*up_1*f2factor;



% % Fhd_m11 = 6*pi*muw*(Dp_m1/2)*up_1*(1+(0.5*Dp_m1/(4*Dp_m1))^0.719)^(-0.719);

% % Fhd_m12 = 6*pi*muw*(Dp_m1/2)*up_2*(1+(0.5*Dp_m1/(4*Dp_m1))^0.719)^(-0.719);

% % Fhd_m13 = 6*pi*muw*(Dp_m1/2)*up_3*(1+(0.5*Dp_m1/(4*Dp_m1))^0.719)^(-0.719);

% % Fhd_m14 = 6*pi*muw*(Dp_m1/2)*up_4*(1+(0.5*Dp_m1/(4*Dp_m1))^0.719)^(-0.719);

% % 

% % Fhd_c11 = 6*pi*muw*(Dp_c1/2)*up_1*(1+(0.5*Dp_c1/(4*Dp_c1))^0.719)^(-0.719);

% % Fhd_c12 = 6*pi*muw*(Dp_c1/2)*up_2*(1+(0.5*Dp_c1/(4*Dp_c1))^0.719)^(-0.719);

% % Fhd_c13 = 6*pi*muw*(Dp_c1/2)*up_3*(1+(0.5*Dp_c1/(4*Dp_c1))^0.719)^(-0.719);

% % Fhd_c14 = 6*pi*muw*(Dp_c1/2)*up_4*(1+(0.5*Dp_c1/(4*Dp_c1))^0.719)^(-0.719);





syms alpha

diff11 = Fhd_f11*cos(alpha) / (Fc_f1 + Fw_f1 + Fhs_f1);

diff12 = Fhd_f12*cos(alpha) / (Fc_f1 + Fw_f1 + Fhs_f1);

diff13 = Fhd_f13*cos(alpha) / (Fc_f1 + Fw_f1 + Fhs_f1);

diff14 = Fhd_f14*cos(alpha) / (Fc_f1 + Fw_f1 + Fhs_f1);



% % % diff11 = -(Fhd_f11-Fd_f11)*cos(alpha) / (Fc_f1 + Fw_f1 + Fhs_f1);

% % % diff12 = -(Fhd_f12-Fd_f12)*cos(alpha) / (Fc_f1 + Fw_f1 + Fhs_f1);

% % % diff13 = -(Fhd_f13-Fd_f13)*cos(alpha) / (Fc_f1 + Fw_f1 + Fhs_f1);

% % % diff14 = -(Fhd_f14-Fd_f14)*cos(alpha) / (Fc_f1 + Fw_f1 + Fhs_f1);

% % % 

% % % diff21 = -(Fhd_m11-Fd_m11)*cos(alpha) / (Fc_m1 + Fw_m1 + Fhs_m1);

% % % diff22 = -(Fhd_m12-Fd_m12)*cos(alpha) / (Fc_m1 + Fw_m1 + Fhs_m1);

% % % diff23 = -(Fhd_m13-Fd_m13)*cos(alpha) / (Fc_m1 + Fw_m1 + Fhs_m1);

% % % diff24 = -(Fhd_m14-Fd_m14)*cos(alpha) / (Fc_m1 + Fw_m1 + Fhs_m1);

% % % 

% % % diff31 = -(Fhd_c11-Fd_c11)*cos(alpha) / (Fc_c1 + Fw_c1 + Fhs_c1);

% % % diff32 = -(Fhd_c12-Fd_c12)*cos(alpha) / (Fc_c1 + Fw_c1 + Fhs_c1);

% % % diff33 = -(Fhd_c13-Fd_c13)*cos(alpha) / (Fc_c1 + Fw_c1 + Fhs_c1);

% % % diff34 = -(Fhd_c14-Fd_c14)*cos(alpha) / (Fc_c1 + Fw_c1 + Fhs_c1);



bound = 0;



%% Case 2: Particle is at the bottom of bubble:



diff41 = Fhd_f11*cos(alpha) / (Fw_f1 - Fc_f1 - Fhs_f1);

diff42 = Fhd_f12*cos(alpha) / (Fw_f1 - Fc_f1 - Fhs_f1);

diff43 = Fhd_f13*cos(alpha) / (Fw_f1 - Fc_f1 - Fhs_f1);

diff44 = Fhd_f14*cos(alpha) / (Fw_f1 - Fc_f1 - Fhs_f1);



% % % diff41 = -(Fhd_f11-Fd_f11)*cos(alpha) / (Fw_f1 - Fc_f1 - Fhs_f1);

% % % diff42 = -(Fhd_f12-Fd_f12)*cos(alpha) / (Fw_f1 - Fc_f1 - Fhs_f1);

% % % diff43 = -(Fhd_f13-Fd_f13)*cos(alpha) / (Fw_f1 - Fc_f1 - Fhs_f1);

% % % diff44 = -(Fhd_f14-Fd_f14)*cos(alpha) / (Fw_f1 - Fc_f1 - Fhs_f1);

% % % 

% % % diff51 = -(Fhd_m11-Fd_m11)*cos(alpha) / (Fw_m1 - Fc_m1 - Fhs_m1);

% % % diff52 = -(Fhd_m12-Fd_m12)*cos(alpha) / (Fw_m1 - Fc_m1 - Fhs_m1);

% % % diff53 = -(Fhd_m13-Fd_m13)*cos(alpha) / (Fw_m1 - Fc_m1 - Fhs_m1);

% % % diff54 = -(Fhd_m14-Fd_m14)*cos(alpha) / (Fw_m1 - Fc_m1 - Fhs_m1);

% % % 

% % % diff61 = -(Fhd_c11-Fd_c11)*cos(alpha) / (Fw_c1 - Fc_c1 - Fhs_c1);

% % % diff62 = -(Fhd_c12-Fd_c12)*cos(alpha) / (Fw_c1 - Fc_c1 - Fhs_c1);

% % % diff63 = -(Fhd_c13-Fd_c13)*cos(alpha) / (Fw_c1 - Fc_c1 - Fhs_c1);

% % % diff64 = -(Fhd_c14-Fd_c14)*cos(alpha) / (Fw_c1 - Fc_c1 - Fhs_c1);





%% PLOTS 



figure(1)

fplot(diff11,[0, 2*pi],'k','LineWidth',2)

hold on

fplot(diff12,[0, 2*pi],'--k','LineWidth',2)

fplot(diff13,[0, 2*pi],'--k','LineWidth',3)

fplot(diff14,[0, 2*pi],'--k','LineWidth',4)

fplot(bound,[0 2*pi],':k','LineWidth',2)

ax = gca;

ax.XTick = 0:pi/2:2*pi;

ax.XTickLabel = {'0', '\pi/2','\pi','3\pi/2','2\pi'};

xlabel('Flow Direction with Respect to Gravity, \alpha')

ylabel('F_h_d / ( F_w + F_c + F_e_s )')

l=legend(' D_p=0.15 mm, u_p_b=1 m/s',' D_p=0.15 mm, u_p_b=4 m/s', ...

    ' D_p=0.15 mm, u_p_b=7 m/s',' D_p=0.15 mm, u_p_b=10 m/s')

title('Case 1: Sand Particle Above Air Bubble','FontWeight','normal','FontSize',18)

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

set(gca,'TickDir','out','FontName','calibri','FontSize',25);

set(gca,'LineWidth',2)

legend('boxoff');

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



figure(2)

fplot(diff41,[0, 2*pi],'k','LineWidth',2)

hold on

fplot(diff42,[0, 2*pi],'--k','LineWidth',2)

fplot(diff43,[0, 2*pi],'--k','LineWidth',3)

fplot(diff44,[0, 2*pi],'--k','LineWidth',4)

fplot(bound,[0 2*pi],':k','LineWidth',2)

ax = gca;

ax.XTick = 0:pi/2:2*pi;

ax.XTickLabel = {'0', '\pi/2','\pi','3\pi/2','2\pi'};

xlabel('Flow Direction with Respect to Gravity, \alpha')

ylabel('F_h_d / ( F_c + F_e_s - F_w)')

l=legend(' D_p=0.15 mm, u_p_b=1 m/s',' D_p=0.15 mm, u_p_b=4 m/s', ...

    ' D_p=0.15 mm, u_p_b=7 m/s',' D_p=0.15 mm, u_p_b=10 m/s')

title('Case 2: Sand Particle Below Air Bubble','FontWeight','normal','FontSize',18)

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

set(gca,'TickDir','out','FontName','calibri','FontSize',25);

set(gca,'LineWidth',2)

legend('boxoff');

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









% figure(2)

% fplot(diff11,[0 2*pi],'k')

% hold on

% fplot(diff21,[0 2*pi],'r')

% fplot(diff31,[0 2*pi],'b')

% fplot(bound,[0 2*pi],'--k')

% ax = gca;

% ax.XTick = 0:pi/2:2*pi;

% ax.XTickLabel = {'0', '\pi/2','\pi','3\pi/2','2\pi'};

% xlabel('Flow Direction w.r.t. Gravity, \alpha')

% ylabel('Hydrodynamic / (Gravity + Capillary + Hydrostatic)')

% legend('u_p=1 m/s, D_p=0.15 mm','u_p=1 m/s, D_p=0.27 mm','u_p=1 m/s, D_p=1.3 mm')

% title('Case 1: Sand Particle Above Air Bubble [Fixed Particle Velocity]')



% figure(4)

% fplot(diff41,[0, 2*pi],'k')

% hold on

% fplot(diff51,[0, 2*pi],'r')

% fplot(diff61,[0, 2*pi],'b')

% fplot(bound,[0 2*pi],'--k')

% ax = gca;

% ax.XTick = 0:pi/2:2*pi;

% ax.XTickLabel = {'0', '\pi/2','\pi','3\pi/2','2\pi'}

% xlabel('Flow Direction w.r.t. Gravity, \alpha')

% ylabel('Hydrodynamic / (Gravity - Capillary - Hydrostatic)')

% legend('u_p=1 m/s, D_p=0.15 mm','u_p=1 m/s, D_p=0.27 mm','u_p=1 m/s, D_p=1.3 mm')

% title('Case 2: Sand Particle Below Air Bubble [Fixed Particle Velocity]')