Can you modify my code?

clear all;
close all;
clc;

%% 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 = load('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 * ho_1) ./ height1;
C2 = (Co2 * ho_2) ./ height2;
C3 = (Co3 * ho_3) ./ height3;

%% Voidage Calculation
e1 = 1 - C1;
e2 = 1 - C2;
e3 = 1 - C3;

%% Gradient Calculation
dy1 = diff(height1);
dy2 = diff(height2);
dy3 = diff(height3);

dt = mean(diff(time));

grad1 = dy1/dt;
grad2 = dy2/dt;
grad3 = dy3/dt;

%% Reynolds Number Calculation
d = 0.0002; % Particle diameter in m
Re1 = abs(grad1) .* d ./ nu;
Re2 = abs(grad2) .* d ./ nu;
Re3 = abs(grad3) .* d ./ nu;

%% 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);
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
Up1 = abs(grad1) .* d;
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
Up2 = abs(grad2) .* d;
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
Up3 = abs(grad3) .* d;
Utz3 = Up3 ./ (e3.^n3);

%% Average Terminal Velocity Calculation
Ut_avg1 = mean([Ut1, Utz1]);
Ut_avg2 = mean([Ut2, Utz2]);
Ut_avg3 = mean([Ut3, Utz3]);

I created the above code, but still need some help. I want to know how to find Up1, Up2 and Up3, which are different gradients from different concentrations. After that, I want to find Ut_averages using Utz and Ut_R at different times. Can you modify this? If you still need more clarification, let me know.

Answer 1

[Mathieu NOE]

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');