Force_Calculation_Turbulence.m (8298 bytes)
clear
clc
format compact
%% Force Calculation Information:
% 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
%% Fundamental Input Parameters:
rhow = 1000; %unit: kg/m3
rhop = 2600; %unit: kg/m3
rhodiff = rhop-rhow;
%% To get Local Turbulence force:
syms Db
%------------------1) x-axis: bubble diameter
%---------------------fixed particle diameter, varying dissipation
Dp = 0.00015;
ep1 = 1;
ep2 = 50;
ep3 = 100;
bm1 = 1.9*(ep1)^(2/3)/(Db^(1/3));
Ft_bm1 = (4*3)*pi*(Dp/2)^3*(rhodiff)*bm1;
bm2 = 1.9*(ep2)^(2/3)/(Db^(1/3));
Ft_bm2 = (4*3)*pi*(Dp/2)^3*(rhodiff)*bm2;
bm3 = 1.9*(ep3)^(2/3)/(Db^(1/3));
Ft_bm3 = (4*3)*pi*(Dp/2)^3*(rhodiff)*bm3;
% % % figure(1)
% % % fplot(1000*Db, Ft_bm1, [0 0.005],'k')
% % % hold on
% % % fplot(1000*Db, Ft_bm2, [0 0.005],'b')
% % % fplot(1000*Db, Ft_bm3, [0 0.005],'r')
% % % legend('D_p=0.15 mm, \epsilon = 1 kW/m^3', ...
% % % 'D_p=0.15 mm, \epsilon = 50 kW/m^3', ...
% % % 'D_p=0.15 mm, \epsilon = 100 kW/m^3')
% % % xlabel('Center Bubble Diameter,D_b (mm)')
% % % ylabel('Local Turbulence Force,F_t (N)')
% % % set(gca,'box','off');
% % % set(gca,'TickDir','out');
% % % set(gca,'YScale','log');
% % % legend('boxoff');
%------------------2) x-axis: bubble diameter
%---------------------fixed dissipation, varying particle diameter
Dp1 = 0.00015;
Dp2 = 0.00027;
Dp3 = 0.00133;
bm4 = 1.9*(ep1)^(2/3)/(Db^(1/3));
Ft_bm4 = (4*3)*pi*(Dp1/2)^3*(rhodiff)*bm4;
bm5 = 1.9*(ep1)^(2/3)/(Db^(1/3));
Ft_bm5 = (4*3)*pi*(Dp2/2)^3*(rhodiff)*bm5;
bm6 = 1.9*(ep1)^(2/3)/(Db^(1/3));
Ft_bm6 = (4*3)*pi*(Dp3/2)^3*(rhodiff)*bm6;
figure(2)
fplot(1000*Db, 1000*Ft_bm4, [0 0.005],'k','LineWidth',2)
hold on
fplot(1000*Db, 1000*Ft_bm5, [0 0.005],'--k','LineWidth',2)
fplot(1000*Db, 1000*Ft_bm6, [0 0.005],':k','LineWidth',2)
legend(' D_p=0.15 mm, \epsilon = 1 kW/m^3', ...
' D_p=0.27 mm, \epsilon = 1 kW/m^3', ...
' D_p=1.33 mm, \epsilon = 1 kW/m^3')
xlabel('Center Bubble Diameter,D_b (mm)')
ylabel('Local Turbulence Force,F_t (mN)')
set(gca,'box','off');
set(gca,'TickDir','out');
set(gca,'YScale','log');
set(gca,'FontName','calibri','FontSize',23)
set(gca,'LineWidth',2)
legend('boxoff');
%------------------3) x-axis: particle diameter
%---------------------fixed dissipation, varying bubble diameter
syms Dp
Db_1 = 0.0005;
Db_2 = 0.0025;
Db_3 = 0.005;
bm7 = 1.9*(ep1)^(2/3)/(Db_1^(1/3));
Ft_bm7 = (4*3)*pi*(Dp/2)^3*(rhodiff)*bm7;
bm8 = 1.9*(ep1)^(2/3)/(Db_2^(1/3));
Ft_bm8 = (4*3)*pi*(Dp/2)^3*(rhodiff)*bm8;
bm9 = 1.9*(ep1)^(2/3)/(Db_3^(1/3));
Ft_bm9 = (4*3)*pi*(Dp/2)^3*(rhodiff)*bm9;
figure(3)
fplot(1000*Dp, 1000*Ft_bm7, [0 0.002],'k','LineWidth',2)
hold on
fplot(1000*Dp, 1000*Ft_bm8, [0 0.002],'--k','LineWidth',2)
fplot(1000*Dp, 1000*Ft_bm9, [0 0.002],':k','LineWidth',2)
legend(' D_b=0.5 mm, \epsilon = 1 kW/m^3', ...
' D_b=2.5 mm, \epsilon = 1 kW/m^3', ...
' D_b=5 mm, \epsilon = 1 kW/m^3')
xlabel('Sand Particle Diameter,D_p (mm)')
ylabel('Local Turbulence Force,F_t (mN)')
set(gca,'box','off');
set(gca,'TickDir','out');
set(gca,'YScale','log');
set(gca,'FontName','calibri','FontSize',23)
set(gca,'LineWidth',2)
legend('boxoff');
%------------------4) x-axis: particle diameter
%---------------------fixed bubble diameter, varying dissipation
bm10 = 1.9*(ep1)^(2/3)/(Db_2^(1/3));
Ft_bm10 = (4*3)*pi*(Dp/2)^3*(rhodiff)*bm10;
bm11 = 1.9*(ep2)^(2/3)/(Db_2^(1/3));
Ft_bm11 = (4*3)*pi*(Dp/2)^3*(rhodiff)*bm11;
bm12 = 1.9*(ep3)^(2/3)/(Db_2^(1/3));
Ft_bm12 = (4*3)*pi*(Dp/2)^3*(rhodiff)*bm12;
% % % figure(4)
% % % fplot(1000*Dp, Ft_bm10, [0 0.002],'k')
% % % hold on
% % % fplot(1000*Dp, Ft_bm11, [0 0.002],'b')
% % % fplot(1000*Dp, Ft_bm12, [0 0.002],'r')
% % % legend('D_b=2 mm, \epsilon = 1 kW/m^3', ...
% % % 'D_b=2 mm, \epsilon = 50 kW/m^3', ...
% % % 'D_b=2 mm, \epsilon = 100 kW/m^3')
% % % xlabel('Sand Particle Diameter,D_p (mm)')
% % % ylabel('Local Turbulence Force,F_t (N)')
% % % set(gca,'box','off');
% % % set(gca,'TickDir','out');
% % % set(gca,'YScale','log');
% % % legend('boxoff');
%------------------5) x-axis: dissipation
%---------------------fixed bubble diameter, varying particle diameter
syms ep
bm13 = 1.9*(ep)^(2/3)/(Db_2^(1/3));
Ft_bm13 = (4*3)*pi*(Dp1/2)^3*(rhodiff)*bm13;
bm14 = 1.9*(ep)^(2/3)/(Db_2^(1/3));
Ft_bm14 = (4*3)*pi*(Dp2/2)^3*(rhodiff)*bm14;
bm15 = 1.9*(ep)^(2/3)/(Db_2^(1/3));
Ft_bm15 = (4*3)*pi*(Dp3/2)^3*(rhodiff)*bm15;
figure(5)
fplot(ep, 1000*Ft_bm13, [1 100],'k','LineWidth',2)
hold on
fplot(ep, 1000*Ft_bm14, [1 100],'--k','LineWidth',2)
fplot(ep, 1000*Ft_bm15, [1 100],':k','LineWidth',2)
legend(' D_b=2.5 mm, D_p=0.15 mm', ...
' D_b=2.5 mm, D_p=0.27 mm', ...
' D_b=2.5 mm, D_p=1.33 mm')
xlabel('Fluid Energy Dissipation,\epsilon (kW/m^3)')
ylabel('Local Turbulence Force,F_t (mN)')
set(gca,'box','off');
set(gca,'TickDir','out');
set(gca,'YScale','log');
set(gca,'FontName','calibri','FontSize',23)
set(gca,'LineWidth',2)
legend('boxoff');
%------------------6) x-axis: dissipation
%---------------------fixed particle diameter, varying bubble diameter
% % % bm16 = 1.9*(ep)^(2/3)/(Db_1^(1/3));
% % % Ft_bm16 = (4*3)*pi*(Dp1/2)^3*(rhodiff)*bm16;
% % % bm17 = 1.9*(ep)^(2/3)/(Db_2^(1/3));
% % % Ft_bm17 = (4*3)*pi*(Dp1/2)^3*(rhodiff)*bm17;
% % % bm18 = 1.9*(ep)^(2/3)/(Db_3^(1/3));
% % % Ft_bm18 = (4*3)*pi*(Dp1/2)^3*(rhodiff)*bm18;
% % %
% % % figure(6)
% % % fplot(ep, Ft_bm16, [1 100],'k')
% % % hold on
% % % fplot(ep, Ft_bm17, [1 100],'b')
% % % fplot(ep, Ft_bm18, [1 100],'r')
% % % legend('D_b=0.5 mm, D_p=0.15 mm', ...
% % % 'D_b=2.5 mm, D_p=0.15 mm', ...
% % % 'D_b=5 mm, D_p=0.15 mm')
% % % xlabel('Fluid Energy Dissipation,\epsilon (kW/m^3)')
% % % ylabel('Local Turbulence Force,F_t (N)')
% % % set(gca,'box','off');
% % % set(gca,'TickDir','out');
% % % set(gca,'YScale','log');
% % % legend('boxoff');
%% To PLOT ALL THREE GRAPHS IN ONE FIGURE:
% % % % % figure(7)
% % % % % subplot(3,1,1)
% % % % % fplot(1000*Db, Ft_bm4, [0 0.005],'k')
% % % % % hold on
% % % % % fplot(1000*Db, Ft_bm5, [0 0.005],'b')
% % % % % fplot(1000*Db, Ft_bm6, [0 0.005],'r')
% % % % % legend( 'D_p=0.15 mm, \epsilon = 1 kW/m^3', ...
% % % % % 'D_p=0.27 mm, \epsilon = 1 kW/m^3', ...
% % % % % 'D_p=1.33 mm, \epsilon = 1 kW/m^3')
% % % % % xlabel('Bubble Diameter,D_b (mm)')
% % % % % %%% ylabel('Local Turbulence Force,F_t (N)')
% % % % % ylim([0.0000001 1])
% % % % % set(gca,'box','off');
% % % % % set(gca,'TickDir','out','FontSize',12);
% % % % % set(gca,'YScale','log');
% % % % % legend('boxoff');
% % % % %
% % % % % subplot(3,1,2)
% % % % % fplot(1000*Dp, Ft_bm7, [0 0.002],'k')
% % % % % hold on
% % % % % fplot(1000*Dp, Ft_bm8, [0 0.002],'b')
% % % % % fplot(1000*Dp, Ft_bm9, [0 0.002],'r')
% % % % % legend('D_b=0.5 mm, \epsilon = 1 kW/m^3', ...
% % % % % 'D_b=2.5 mm, \epsilon = 1 kW/m^3', ...
% % % % % 'D_b=5.0 mm, \epsilon = 1 kW/m^3')
% % % % % xlabel('Sand Particle Diameter,D_p (mm)')
% % % % % ylabel('Local Turbulence Force,F_t (N)')
% % % % % ylim([0.00000000001 0.001])
% % % % % set(gca,'box','off');
% % % % % set(gca,'TickDir','out','FontSize',12);
% % % % % set(gca,'YScale','log');
% % % % % legend('boxoff');
% % % % %
% % % % % subplot(3,1,3)
% % % % % fplot(ep, Ft_bm13, [1 100],'k')
% % % % % hold on
% % % % % fplot(ep, Ft_bm14, [1 100],'b')
% % % % % fplot(ep, Ft_bm15, [1 100],'r')
% % % % % legend('D_b=2.5 mm, D_p=0.15 mm', ...
% % % % % 'D_b=2.5 mm, D_p=0.27 mm', ...
% % % % % 'D_b=2.5 mm, D_p=1.33 mm')
% % % % % xlabel('Fluid Energy Dissipation,\epsilon (kW/m^3)')
% % % % % %%% ylabel('Local Turbulence Force,F_t (N)')
% % % % % ylim([0.0000001 1])
% % % % % set(gca,'box','off');
% % % % % set(gca,'TickDir','out','FontSize',12);
% % % % % set(gca,'YScale','log');
% % % % % legend('boxoff');
- Log in to post comments