Force_Calculation_Hydrodynamic_Drag.m (3513 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



%% Fixed Parameters: 



gamma = 0.073;  %unit: N/m

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

rhow = 1000;  %unit: kg/m3

g = 9.8;  %unit: m/s2

muw = 0.001;  %unit: kg/m-s

rhop = 2600;  %unit: kg/m3



Dp_f = 0.00015;  %unit: m

Dp_m = 0.00027;  %unit: m

Dp_c = 0.00133;  %unit: m





%% To get Drag force:

syms u_f



%----------- For Fine Sand Case:

Rp_f = Dp_f/2;

H_f = 0.1*Dp_f;

Rb_f = 13*Rp_f;



vs_f = (2/9)*(rhop-rhow)*g*(Rp_f)^2/muw;



f1f = (1+(Rp_f/(4*H_f))^0.719)^-0.719;

f2f_prime = (1.707+(H_f/Rp_f))/(0.836+(H_f/Rp_f));

f2f_pprime = (2.656+(H_f/Rp_f))/(1.440+(H_f/Rp_f));

hf = (H_f+Rp_f)/Rb_f;

f2f = (f2f_prime - hf*f2f_pprime)/(1-hf);



F_dragf = 6*pi*muw*Rp_f*(vs_f+u_f)*f1f+6*pi*muw*Rp_f*u_f*f2f;



%---------- For Medium Sand Case: 

Rp_m = Dp_m/2;

H_m = 0.1*Dp_m;

Rb_m = 7*Rp_m;



vs_m = (2/9)*(rhop-rhow)*g*(Rp_m)^2/muw;



f1m = (1+(Rp_m/(4*H_m))^0.719)^-0.719;

f2m_prime = (1.707+(H_m/Rp_m))/(0.836+(H_m/Rp_m));

f2m_pprime = (2.656+(H_m/Rp_m))/(1.440+(H_m/Rp_m));

hm = (H_m+Rp_m)/Rb_m;

f2m = (f2m_prime - hm*f2m_pprime)/(1-hm);



F_dragm = 6*pi*muw*Rp_m*(vs_m+u_f)*f1m+6*pi*muw*Rp_m*u_f*f2m;



%---------- For Coarse Sand Case: 

Rp_c = Dp_c/2;

H_c = 0.1*Dp_c;

Rb_c = 2*Rp_c;



vs_c = (2/9)*(rhop-rhow)*g*(Rp_c)^2/muw;



f1c = (1+(Rp_c/(4*H_c))^0.719)^-0.719;

f2c_prime = (1.707+(H_c/Rp_c))/(0.836+(H_c/Rp_c));

f2c_pprime = (2.656+(H_c/Rp_c))/(1.440+(H_c/Rp_c));

hc = (H_c+Rp_c)/Rb_c;

f2c = (f2c_prime - hc*f2c_pprime)/(1-hc);



F_dragc = 6*pi*muw*Rp_c*(vs_c+u_f)*f1c+6*pi*muw*Rp_c*u_f*f2c-1.1*10^-5;



%% FIGURE

figure(1)

fplot(u_f, 1000*0.75*F_dragf, [0   8],'k','LineWidth',1.8)

hold on

fplot(u_f, 1000*0.75*F_dragm, [0   8],'--k','LineWidth',2)

fplot(u_f, 1000*0.75*F_dragc, [0   8],':k','LineWidth',3)

ylim([0   0.2])

xlabel('Fluid Velocity, v_f (m/s)')

ylabel('Hydrodynamic Drag Force, F_d (mN)')

l=legend(' D_p=0.15 mm, H=0.1D_p', ...

       ' D_p=0.27 mm, H=0.1D_p', ...

       ' D_p=1.33 mm, H=0.1D_p')

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

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

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

set(gca,'LineWidth',2)

legend('boxoff');

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