Force_Calculation_Shear_Rate_vs_Bubble_Diameter.m (11958 bytes)

clear

clc

format compact



%% Input Parameters:



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

gamma = 0.073;  %fluid, unit: N/m

theta = 91; %unit: degree

rhow = 1000;  %unit: kg/m3

rhop = 2600;  %unit: kg/m3

rhodiff = rhop-rhow; 

g = 9.8;  %unit: m/s2



Dp_f = 0.00015;  % FINE, unit: m

Dp_m = 0.00027;  % MEDIUM, unit:m

Dp_c = 0.001;  % COARSE, unit: m



ep1 = 10;  %unit: kW/m3

ep2 = 50;

ep3 = 100;



%% CASE 1： 

% Particle at the bottom of the bubble:

% drag and turbulence are opposite to each other:

% ------------ Case not analyzable ----------------



% % % syms W

% % % Rp_f = Dp_f/2;

% % % H_f = 0.1*Dp_f;

% % % 

% % % vs = (2/9)*(Rp_f)^2*g*rhodiff/muw;

% % % f1 = (1+((Rp_f)/(4*H_f))^0.719)^(1/0.719);

% % % f2_prime = (1.707+H_f/Rp_f)/(0.836+H_f/Rp_f);

% % % f2_pprime = (2.656+H_f/Rp_f)/(1.44+H_f/Rp_f);

% % % h = (H_f+Rp_f)/(10*Rp_f);

% % % f2 = (f2_prime-h*f2_pprime)/(1-h);



% % % RHS1t = (1/3*pi*(Dp_f)^2*rhodiff*1.9*ep1^(2/3))/(-3*pi*muw*Dp_f*(W+vs)*f1+3*pi*muw*Dp_f*W*f2);

% % % RHS2t = (1/3*pi*(Dp_f)^2*rhodiff*1.9*ep2^(2/3))/(-3*pi*muw*Dp_f*(W+vs)*f1+3*pi*muw*Dp_f*W*f2);

% % % RHS3t = (1/3*pi*(Dp_f)^2*rhodiff*1.9*ep3^(2/3))/(-3*pi*muw*Dp_f*(W+vs)*f1+3*pi*muw*Dp_f*W*f2);



% % % figure(1)

% % % fplot(-W, 2*(RHS1t)^(1/3), [-7.8   -4],'-k')

% % % hold on

% % % fplot(-W, 2*(RHS2t)^(1/3), [-7.8   -4],'--k')

% % % fplot(-W, 2*(RHS3t)^(1/3), [-7.8   -4],'-.k')

% % % xlim([0   7.8])

% % % ylim([0   20])





%% Case 2-A and 2-B:

% Particle at the left:

% A drag and turbulence are in opposite direction:

% B laminar, no turbulence

syms Db



X_max = 5.5;



% --------- fine sand case

Rp_f = Dp_f/2;

H_f = 0.1*Dp_f;



vs = (2/9)*(Rp_f)^2*g*rhodiff/muw;

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

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

f2_pprime = (2.656+H_f/Rp_f)/(1.44+H_f/Rp_f);

h = (H_f+Rp_f)/(Db);

f2 = (f2_prime-h*f2_pprime)/(1-h);



STAR1t = 0.5*pi*Dp_f*gamma*(1-cos(theta))+0.25*pi*(Dp_f)^2*(1-cos(theta))*(2*gamma/Db-rhow*g*Db/2)+(4/3)*pi*(Dp_f/2)^3*rhodiff*(1.9*ep1^(2/3))/(Db^(1/3));

Vf1t = (vs*f1+STAR1t/(3*pi*muw*Dp_f))/(f1+f2);



STAR2t = 0.5*pi*Dp_f*gamma*(1-cos(theta))+0.25*pi*(Dp_f)^2*(1-cos(theta))*(2*gamma/Db-rhow*g*Db/2)+(4/3)*pi*(Dp_f/2)^3*rhodiff*(1.9*ep2^(2/3))/(Db^(1/3));

Vf2t = (vs*f1+STAR2t/(3*pi*muw*Dp_f))/(f1+f2);



STAR3t = 0.5*pi*Dp_f*gamma*(1-cos(theta))+0.25*pi*(Dp_f)^2*(1-cos(theta))*(2*gamma/Db-rhow*g*Db/2)+(4/3)*pi*(Dp_f/2)^3*rhodiff*(1.9*ep3^(2/3))/(Db^(1/3));

Vf3t = (vs*f1+STAR3t/(3*pi*muw*Dp_f))/(f1+f2);



STAR3l = 0.5*pi*Dp_f*gamma*(1-cos(theta))+0.25*pi*(Dp_f)^2*(1-cos(theta))*(2*gamma/Db-rhow*g*Db/2);

Vf3l = (vs*f1+STAR3l/(3*pi*muw*Dp_f))/(f1+f2)-2;



figure(2)

fplot(1000*Db, Vf3t/0.01,[0.00015 0.0008],'--k','LineWidth',1.5)

hold on

fplot(1000*Db, Vf2t/0.01,[0.00015 0.0008],'-.k','LineWidth',1.5)

fplot(1000*Db, Vf1t/0.01,[0.00015 0.0008],':k','LineWidth',1.5)

fplot(1000*Db, Vf3l/0.01,[0.0008 0.0055],'k','LineWidth',1.5)

xline(5.5,'--k')

yline(400, '--k')

xlabel('Bubble Diameter for a Stable Liquid Marble, D_b (mm)')

ylabel({'Particle Shear Rate for Stable Liquid'; 'Marble Shape, d\gamma/dt (s^-^1)'})

l=legend(' D_p=0.15 mm, tubulent \epsilon=100 kW/m^3', ...

    ' D_p=0.15 mm, turbulent \epsilon=50 kW/m^3', ...

    ' D_p=0.15 mm, turbulent \epsilon=1 kW/m^3', ...

    ' D_p=0.15 mm, laminar')

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

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

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

set(gca,'LineWidth',2)

legend('boxoff');

xlim([0  8])

ylim([0  1000])

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



% --------- medium sand case

Rp_m = Dp_m/2;

H_m = 0.1*Dp_m;



vsm = (2/9)*(Rp_m)^2*g*rhodiff/muw;

f1m = (1+((Rp_m)/(4*H_m))^0.719)^(1/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.44+H_m/Rp_m);

hm = (H_m+Rp_m)/(Db);

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



STAR4t = 0.5*pi*Dp_m*gamma*(1-cos(theta))+0.25*pi*(Dp_m)^2*(1-cos(theta))*(2*gamma/Db-rhow*g*Db/2)+(4/3)*pi*(Dp_m/2)^3*rhodiff*(1.9*ep1^(2/3))/(Db^(1/3));

Vf4t = (vsm*f1m+STAR4t/(3*pi*muw*Dp_m))/(f1m+f2m);



STAR5t = 0.5*pi*Dp_m*gamma*(1-cos(theta))+0.25*pi*(Dp_m)^2*(1-cos(theta))*(2*gamma/Db-rhow*g*Db/2)+(4/3)*pi*(Dp_m/2)^3*rhodiff*(1.9*ep2^(2/3))/(Db^(1/3));

Vf5t = (vsm*f1m+STAR5t/(3*pi*muw*Dp_m))/(f1m+f2m);



STAR6t = 0.5*pi*Dp_m*gamma*(1-cos(theta))+0.25*pi*(Dp_m)^2*(1-cos(theta))*(2*gamma/Db-rhow*g*Db/2)+(4/3)*pi*(Dp_m/2)^3*rhodiff*(1.9*ep3^(2/3))/(Db^(1/3));

Vf6t = (vsm*f1m+STAR6t/(3*pi*muw*Dp_m))/(f1m+f2m);



STAR6l = 0.5*pi*Dp_m*gamma*(1-cos(theta))+0.25*pi*(Dp_m)^2*(1-cos(theta))*(2*gamma/Db-rhow*g*Db/2);

Vf6l = (vsm*f1m+STAR6l/(3*pi*muw*Dp_m))/(f1m+f2m)-2;



figure(3)

fplot(1000*Db, Vf6t/0.01,[0.00025 0.0014],'--k','LineWidth',1.5)

hold on

fplot(1000*Db, Vf5t/0.01,[0.00025 0.0014],'-.k','LineWidth',1.5)

fplot(1000*Db, Vf4t/0.01,[0.00025 0.0014],':k','LineWidth',1.5)

fplot(1000*Db, Vf6l/0.01,[0.0014 0.0055],'k','LineWidth',1.5)

xline(5.5,'--k')

yline(400, '--k')

xlabel('Bubble Diameter for a Stable Liquid Marble, D_b (mm)')

ylabel({'Particle Shear Rate for Stable Liquid'; 'Marble Shape, d\gamma/dt (s^-^1)'})

l=legend(' D_p=0.27 mm, tubulent \epsilon=100 kW/m^3', ...

    ' D_p=0.27 mm, turbulent \epsilon=50 kW/m^3', ...

    ' D_p=0.27 mm, turbulent \epsilon=1 kW/m^3', ...

    ' D_p=0.27 mm, laminar')

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

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

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

set(gca,'LineWidth',2)

legend('boxoff');

xlim([0  8])

ylim([0  1000])

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



% --------- coarse sand case

Rp_c = Dp_c/2;

H_c = 0.1*Dp_c;



vsc = (2/9)*(Rp_c)^2*g*rhodiff/muw;

f1c = (1+((Rp_c)/(4*H_c))^0.719)^(1/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.44+H_c/Rp_c);

hc = (H_c+Rp_c)/(Db);

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



STAR7t = 0.5*pi*Dp_c*gamma*(1-cos(theta))+0.25*pi*(Dp_c)^2*(1-cos(theta))*(2*gamma/Db-rhow*g*Db/2)+(4/3)*pi*(Dp_c/2)^3*rhodiff*(1.9*ep1^(2/3))/(Db^(1/3));

Vf7t = (vsc*f1c+STAR7t/(3*pi*muw*Dp_c))/(f1c+f2c);



STAR8t = 0.5*pi*Dp_c*gamma*(1-cos(theta))+0.25*pi*(Dp_c)^2*(1-cos(theta))*(2*gamma/Db-rhow*g*Db/2)+(4/3)*pi*(Dp_c/2)^3*rhodiff*(1.9*ep2^(2/3))/(Db^(1/3));

Vf8t = (vsc*f1c+STAR8t/(3*pi*muw*Dp_c))/(f1c+f2c);



STAR9t = 0.5*pi*Dp_c*gamma*(1-cos(theta))+0.25*pi*(Dp_c)^2*(1-cos(theta))*(2*gamma/Db-rhow*g*Db/2)+(4/3)*pi*(Dp_c/2)^3*rhodiff*(1.9*ep3^(2/3))/(Db^(1/3));

Vf9t = (vsc*f1c+STAR9t/(3*pi*muw*Dp_c))/(f1c+f2c);



STAR9l = 0.5*pi*Dp_c*gamma*(1-cos(theta))+0.25*pi*(Dp_c)^2*(1-cos(theta))*(2*gamma/Db-rhow*g*Db/2);

Vf9l = (vsc*f1c+STAR9l/(3*pi*muw*Dp_c))/(f1c+f2c);





figure(4)

fplot(1000*Db, Vf9t/0.011,[0.0012 0.0043],'--k','LineWidth',1.5)

hold on

fplot(1000*Db, Vf8t/0.011,[0.0012 0.0043],'-.k','LineWidth',1.5)

fplot(1000*Db, Vf7t/0.011,[0.0012 0.0043],':k','LineWidth',1.5)

fplot(1000*Db, Vf9l/0.015,[0.0043 0.0055],'k','LineWidth',1.5)

xline(5.5,'--k')

yline(400, '--k')

xlabel('Bubble Diameter for a Stable Liquid Marble, D_b (mm)')

ylabel({'Particle Shear Rate for Stable Liquid'; 'Marble Shape, d\gamma/dt (s^-^1)'})

l=legend(' D_p=1.33 mm, tubulent \epsilon=100 kW/m^3', ...

    ' D_p=1.33 mm, turbulent \epsilon=50 kW/m^3', ...

    ' D_p=1.33 mm, turbulent \epsilon=1 kW/m^3', ...

    ' D_p=1.33 mm, laminar')

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

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

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

set(gca,'LineWidth',2)

legend('boxoff');

xlim([0  8])

ylim([0  1000])

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





%% Case 3-A and 3-B: 

% Particle at the left-45deg position:

% A drag and turbulence are in opposite direction:

% B laminar, no turbulence

% not the worst scenario

% syms Db

% % ---------- fine sand case

% Rp_f = Dp_f/2;

% H_f = 0.1*Dp_f;

% 

% vsf = (2/9)*(Rp_f)^2*g*rhodiff/muw;

% f1f = (1+((Rp_f)/(4*H_f))^0.719)^(1/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.44+H_f/Rp_f);

% hf = (H_f+Rp_f)/(Db);

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

% 

% A_f = 0.5*pi*Dp_f*gamma*(1-cos(theta));

% B_f = (1/6)*pi*(Dp_f)^3*rhop*g-(1/8)*pi*(Dp_f)^3*rhow*g*(2/3+cos(theta/2)-(1/3)*(cos(theta/2))^3);

% C_f = 0.25*pi*(Dp_f)^2*(1-cos(theta))*(2*gamma/Db-rhow*g*Db/2);

% 

% Box_f1t = (A_f-sin(45)*B_f+sin(45)*(4/3)*pi*(Dp_f)^3*rhodiff*1.9*ep1^(2/3)/((Db)^(1/3))+C_f)/sin(45);

% W1t = (vsf*f1f+Box_f1t/(3*pi*muw*Dp_f))/(f1f+f2f)-1;

% 

% Box_f2t = (A_f-sin(45)*B_f+sin(45)*(4/3)*pi*(Dp_f)^3*rhodiff*1.9*ep2^(2/3)/((Db)^(1/3))+C_f)/sin(45);

% W2t = (vsf*f1f+Box_f2t/(3*pi*muw*Dp_f))/(f1f+f2f)-1;

% 

% Box_f3t = (A_f-sin(45)*B_f+sin(45)*(4/3)*pi*(Dp_f)^3*rhodiff*1.9*ep3^(2/3)/((Db)^(1/3))+C_f)/sin(45);

% W3t = (vsf*f1f+Box_f3t/(3*pi*muw*Dp_f))/(f1f+f2f)-1;

% 

% Box_f1l = (A_f-sin(45)*B_f+C_f)/sin(45);

% W1l = (vsf*f1f+Box_f1l/(3*pi*muw*Dp_f))/(f1f+f2f)-3;

% 

% figure(5)

% fplot(1000*Db, W1l/0.01, [0.0008  0.005],'-k')

% hold on

% fplot(1000*Db, W1t/0.01, [0.0001  0.0008],'-.b')

% fplot(1000*Db, W2t/0.01, [0.0001  0.0008],'-.r')

% fplot(1000*Db, W3t/0.01, [0.0001  0.0008],'-.k')

% 

% xlim([0   5])

% ylim([0 1000])

% 

% % ---------- medium sand case

% Rp_m = Dp_m/2;

% H_m = 0.1*Dp_m;

% 

% vsm = (2/9)*(Rp_m)^2*g*rhodiff/muw;

% f1m = (1+((Rp_m)/(4*H_m))^0.719)^(1/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.44+H_m/Rp_m);

% hm = (H_m+Rp_m)/(Db);

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

% 

% A_m = 0.5*pi*Dp_m*gamma*(1-cos(theta));

% B_m = (1/6)*pi*(Dp_m)^3*rhop*g-(1/8)*pi*(Dp_m)^3*rhow*g*(2/3+cos(theta/2)-(1/3)*(cos(theta/2))^3);

% C_m = 0.25*pi*(Dp_m)^2*(1-cos(theta))*(2*gamma/Db-rhow*g*Db/2);

% 

% Box_m1t = (A_m-sin(45)*B_m+sin(45)*(4/3)*pi*(Dp_m)^3*rhodiff*1.9*ep1^(2/3)/((Db)^(1/3))+C_m)/sin(45);

% W4t = (vsm*f1m+Box_m1t/(3*pi*muw*Dp_m))/(f1m+f2m)-1;

% 

% Box_m2t = (A_m-sin(45)*B_m+sin(45)*(4/3)*pi*(Dp_m)^3*rhodiff*1.9*ep2^(2/3)/((Db)^(1/3))+C_m)/sin(45);

% W5t = (vsm*f1m+Box_m2t/(3*pi*muw*Dp_m))/(f1m+f2m)-1;

% 

% Box_m3t = (A_m-sin(45)*B_m+sin(45)*(4/3)*pi*(Dp_m)^3*rhodiff*1.9*ep3^(2/3)/((Db)^(1/3))+C_m)/sin(45);

% W6t = (vsm*f1m+Box_m3t/(3*pi*muw*Dp_m))/(f1m+f2m)-1.5;

% 

% Box_m1l = (A_m-sin(45)*B_m+C_m)/sin(45);

% W4l = (vsm*f1m+Box_m1l/(3*pi*muw*Dp_m))/(f1m+f2m)-3;

% 

% figure(6)

% fplot(1000*Db, W4l/0.01, [0.0014  0.005],'-k')

% hold on

% fplot(1000*Db, W4t/0.01, [0.0002  0.0014],'-.b')

% fplot(1000*Db, W5t/0.01, [0.0002  0.0014],'-.r')

% fplot(1000*Db, W6t/0.01, [0.0002  0.0014],'-.k')

% 

% xlim([0   5])

% ylim([0 1000])

% 

% % ---------- coarse sand case

% Rp_c = Dp_c/2;

% H_c = 0.1*Dp_c;

% 

% vsc = (2/9)*(Rp_c)^2*g*rhodiff/muw;

% f1c = (1+((Rp_c)/(4*H_c))^0.719)^(1/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.44+H_c/Rp_c);

% hc = (H_c+Rp_c)/(Db);

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

% 

% A_c = 0.5*pi*Dp_c*gamma*(1-cos(theta));

% B_c = (1/6)*pi*(Dp_c)^3*rhop*g-(1/8)*pi*(Dp_c)^3*rhow*g*(2/3+cos(theta/2)-(1/3)*(cos(theta/2))^3);

% C_c = 0.25*pi*(Dp_c)^2*(1-cos(theta))*(2*gamma/Db-rhow*g*Db/2);

% 

% Box_c1t = (A_c-sin(45)*B_c+sin(45)*(4/3)*pi*(Dp_c)^3*rhodiff*1.9*ep1^(2/3)/((Db)^(1/3))+C_c)/sin(45);

% W7t = (vsc*f1c+Box_c1t/(3*pi*muw*Dp_c))/(f1c+f2c);

% 

% Box_c2t = (A_c-sin(45)*B_c+sin(45)*(4/3)*pi*(Dp_c)^3*rhodiff*1.9*ep2^(2/3)/((Db)^(1/3))+C_c)/sin(45);

% W8t = (vsc*f1c+Box_c2t/(3*pi*muw*Dp_c))/(f1c+f2c);

% 

% Box_c3t = (A_c-sin(45)*B_c+sin(45)*(4/3)*pi*(Dp_c)^3*rhodiff*1.9*ep3^(2/3)/((Db)^(1/3))+C_c)/sin(45);

% W9t = (vsc*f1c+Box_c3t/(3*pi*muw*Dp_c))/(f1c+f2c);

% 

% Box_c1l = (A_c-sin(45)*B_c+C_c)/sin(45);

% W7l = (vsc*f1c+Box_c1l/(3*pi*muw*Dp_c))/(f1c+f2c)-3.5;

% 

% figure(7)

% fplot(1000*Db, W7l/0.01, [0.00385  0.005],'-k')

% hold on

% fplot(1000*Db, W7t/0.01, [0.0009  0.00385],'-.b')

% fplot(1000*Db, W8t/0.01, [0.0009  0.00385],'-.r')

% fplot(1000*Db, W9t/0.01, [0.0009  0.00385],'-.k')

% 

% xlim([0   5])

% ylim([0 1000])