hello
I can make following suggestions :
- you can use readmatrix instead of load to read excel files
- also gradient exist now , so no need to use the older diff (also the resulting array will have the same size as the input array which is not the case with diff)
- reynold's number: the denominator should be mu (and not nu)
- solving the non linear equation can be done by hand , no need for fsolve here (and that is fine for me as I don't have the Optimisation toolbox
% solving (4.8 - n)/(n - 2.4) - 0.43*Ar^(0.57)*(1 - 2.4*(dp/D)^2.4) = 1
% is equivalent to :
% (4.8 - n)/(n - 2.4) = A , with A = 1 + 0.43*Ar^(0.57)*(1 - 2.4*(dp/D)^2.4) ; (A = 1.0303e+05)
% we have then n_R = (4.8 - 2.4*A)/(1+A);
- Up1,Up2,Up3 should be computed before the lines :
Ut_R1 = Up1 ./ (e1.^n_R);
Ut_R2 = Up2 ./ (e2.^n_R);
Ut_R2 = Up3 ./ (e3.^n_R);
otherwise the code fails
- minor bugs corrected like ho1/ho2/ho3 in the init section becomes ho_1/ho_2/ho_3 in the main code , corrected
I could plot the Up1,Up2,Up3 results vs time, but I don't know where Ut1/Ut2/Ut3 are computed - they show up in those lines
Ut_avg1 = mean([Ut1, Utz1]);
Ut_avg2 = mean([Ut2, Utz2]);
Ut_avg3 = mean([Ut3, Utz3]);
%% Constants
g = 9.81; % Acceleration due to gravity in m/s^2
mu = 8.9E-4; % Kinematic viscosity of air in m^2/s
rho_p = 2930; % Density of particles in kg/m^3
rho_f = 1000; % Density of fluid (water) in kg/m^3
dp = 125E-6 % particle size in m
D = 5.05 %
m1 = 40 % mass of CaCO3 in g
m2 = 60 % mass of CaCO3 in g
m3 = 80 % mass of CaCO3 in g
Vf = 1.5 % volume of the fluid in L
%% Import data
data = readmatrix('data.xlsx');
time = data(:,1); % time (s)
height1 = data(:,2); % height for concentration 40g (cm)
height2 = data(:,3); % height for concentration 60g (cm)
height3 = data(:,4); % height for concentration 80g (cm)
%% Concentration Calculation
Co = 3; % Initial concentration in kg/m^3
ho1 = max(height1); % initial heigh in com
ho2 = max(height2);
ho3 = max(height3);
Vs1 = m1/(rho_p*1000);
Vs2 = m2/(rho_p*1000);
Vs3 = m3/(rho_p*1000);
Co1 = Vs1/(Vs1+Vf);
Co2 = Vs2/(Vs2+Vf);
Co3 = Vs3/(Vs3+Vf);
C1 = (Co1 * ho1) ./ height1;
C2 = (Co2 * ho2) ./ height2;
C3 = (Co3 * ho3) ./ height3;
%% Voidage Calculation
e1 = 1 - C1;
e2 = 1 - C2;
e3 = 1 - C3;
%% Gradient Calculation
dt = mean(diff(time));
grad1 = gradient(height1,dt);
grad2 = gradient(height2,dt);
grad3 = gradient(height3,dt);
%% Reynolds Number Calculation
d = 0.0002; % Particle diameter in m
Re1 = abs(grad1) .* d ./ mu;
Re2 = abs(grad2) .* d ./ mu;
Re3 = abs(grad3) .* d ./ mu;
%% Terminal Velocity Calculation - Richardson's Correlation
Ar = (rho_p - rho_f) * g * D^3 ./ (mu);
% n_R = fsolve(@(n) (4.8 - n)/(n - 2.4) - 0.43*Ar^(0.57)*(1 - 2.4*(dp/D)^2.4), 1);
% solving (4.8 - n)/(n - 2.4) - 0.43*Ar^(0.57)*(1 - 2.4*(dp/D)^2.4) = 1
% is equivalent to :
% (4.8 - n)/(n - 2.4) = A , with A = 1 + 0.43*Ar^(0.57)*(1 - 2.4*(dp/D)^2.4);
% we have then n_R = (4.8 - 2.4*A)/(1+A);
A = 1 + 0.43*Ar^(0.57)*(1 - 2.4*(dp/D)^2.4);
n_R = (4.8 - 2.4*A)/(1+A);
Up1 = abs(grad1) .* d;
Up2 = abs(grad2) .* d;
Up3 = abs(grad3) .* d;
Ut_R1 = Up1 ./ (e1.^n_R);
Ut_R2 = Up2 ./ (e2.^n_R);
Ut_R2 = Up3 ./ (e3.^n_R);
%% Terminal Velocity Calculation - Richardson-Zaki Equation
if Re1 < 0.3
n1 = 4.65;
elseif Re1 > 500
n1 = 0.4;
else
n1 = 2.78 - 2.72*log10(Re1);
end
Utz1 = Up1 ./ (e1.^n1);
if Re2 < 0.3
n2 = 4.65;
elseif Re2 > 500
n2 = 0.4;
else
n2 = 2.78 - 2.72*log10(Re2);
end
Utz2 = Up2 ./ (e2.^n2);
if Re3 < 0.3
n3 = 4.65;
elseif Re3 > 500
n3 = 0.4;
else
n3 = 2.78 - 2.72*log10(Re3);
end
Utz3 = Up3 ./ (e3.^n3);
% %% Average Terminal Velocity Calculation
% Ut_avg1 = mean([Ut1, Utz1]);
% Ut_avg2 = mean([Ut2, Utz2]);
% Ut_avg3 = mean([Ut3, Utz3]);
figure
plot(time,Up1,time,Up2,time,Up3);
legend('Up1','Up2','Up3');
