Warning: Imaginary parts of complex X and/or Y arguments ignored.
4 views (last 30 days)
Show older comments
I am using the below code in the below function to generate a graph, however the graph I am trying to replicate is not what matlab is giving.
function [dhdt] = thesis(t,h)
% Oil/Brine Systems Main Characteristics
%% Oil/Water System at 30C (Paraffin)
mu1 = 3.8; % Viscosity in mPa.s-1
rho_O1 = 797; % Oil Density in kg/m3
rho_B1 = 998; % Brine Density in kg/m3
IFT1 = 21; % Interfacial Tension mN/m
g = 9.81; % m2/s
% Concentration for this system is 80ppm
%% Oil/Water System at 40C (Crude Oil A)
mu2 = 5; % Viscosity in mPa.s-1
rho_O2 = 833; % Oil Density in kg/m3
rho_B2 = 994; % Brine Density in kg/m3
IFT2 = 5; % Interfacial Tension mN/m
% Concentration for this system is 200ppm
%% Oil/Water System at 50C (Crude Oil B)
mu3 = 6; % Viscosity in mPa.s-1
rho_O3 = 834; % Oil Density in kg/m3
rho_B3 = 989; % Brine Density in kg/m3
IFT3 = 3; % Interfacial Tension mN/m
% Concentration for this system is 200ppm
%% Oil/Water System at 60C (Crude Oil C)
mu4 =4.9; % Viscosity in mPa.s-1
rho_O4 = 826; % Oil Density in kg/m3
rho_B4 = 985; % Brine Density in kg/m3
IFT4 = 1; % Interfacial Tension mN/m
%Concentration for this system is 200ppm
nd = 1000; % No. of Droplets
Vol = 900; % Liquid volume of emulsion (ml)
l = 0.5; % Mean Distance between droplets
alpha = 0.08; % Empirical Collision Effiency Parameter
D0 = 300; % Initial Droplet Diameter (microns)
%% Sedimentation interface, hs: dhs/dt
Pr0 = ((nd*pi*D0^3)/6)/Vol; % Initial Volume Fraction of droplet
Prm = ((nd*pi*((D0+l)^3))/6)/Vol; % Maximum Volume Fraction of droplet
delrho1 = rho_B2 - rho_O2; % difference between the dispersed water and continuous oil phase
Vsto = (delrho1*g*(D0^2))/18*mu2; % Settling Velocity of Hard Spheres (stoke's velocity)
fPr = (1-Pr0)^5.3; % Dimensionless
x1 = 60; % (hs - hd) in mm, guessed value
K1 = ((2/3)*alpha*((Vsto^2)/D0))*((fPr^2)/((Prm/fPr)^1/3)-1);
dhdt = x1*(K1*t-(Vsto*fPr));
end
Command window code to generate graph:
tspan = [0 180];
[t, dhdt] = ode45(@thesis, tspan, [240]);
% where H = 240mm (Height of column)
% plotting graph
plot (t, dhdt(:,1), '-o')
grid on
ylabel ('H(mm)')
xlabel ('Time (s)')
title ('Simulation of Sedimentation-Coalesence Experiment')
When I run this code in the command window for the above function, I am getting this response :
Warning: Imaginary parts of complex X and/or Y arguments ignored.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1049410/image.png)
I set the initial value as 240 but even that is not reflected on the graph, Can you assist please?
0 Comments
Accepted Answer
Paul
on 29 Jun 2022
Edited: Paul
on 29 Jun 2022
Hi Matthew,
Why is dhdt the output argument of this line?
[t, dhdt] = ode45(@thesis, tspan, [240]);
Shouldn't this be
[t, h] = ode45(@thesis, tspan, [240]);
I guess it doesn't matter what the output variable is called, but it's confusing.
Anyway, the scale on the plot goes up to 1e124. The initial value of dhdt is 240, as it should be.
The solution becomes immediately complex because 1-Pr0 is a large negative number, and then rasing that to the power 5.3 is complex.
Also, is the equation for K1 correct? Should this term ((Prm/fPr)^1/3) really be ((Prm/fPr)^(1/3)) ? The former raise Prm/fPr to the 1 power and then divides the result by 3, whereas the latter raises Prm/fPr to the 1/3 power
5 Comments
Dennis
on 29 Jun 2022
Is there a problems with units? If D0 is in microns and Vol is in mL the calculation of Pr0 looks very odd to me.
Sam Chak
on 29 Jun 2022
Here you can probably investigate what went wrong. Pr0 and Prm are very large values due to the Initial Droplet Diameter, D0 that is measure in microns. Please check the unit conversion.
fPr gives you a complex-valued number. Please check if it should be (Pr0 - 1)^5.3 or unit conversion issue where (1 - Pr0) > 0. If you settle this, you will settle K1, though it is still
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1050330/image.png)
Should x1 be defined as hd - hs instead?
mu1 = 3.8; % Viscosity in mPa.s-1
rho_O1 = 797; % Oil Density in kg/m3
rho_B1 = 998; % Brine Density in kg/m3
IFT1 = 21; % Interfacial Tension mN/m
g = 9.81; % m2/s
mu2 = 5; % Viscosity in mPa.s-1
rho_O2 = 833; % Oil Density in kg/m3
rho_B2 = 994; % Brine Density in kg/m3
IFT2 = 5; % Interfacial Tension mN/m
mu3 = 6; % Viscosity in mPa.s-1
rho_O3 = 834; % Oil Density in kg/m3
rho_B3 = 989; % Brine Density in kg/m3
IFT3 = 3; % Interfacial Tension mN/m
mu4 = 4.9; % Viscosity in mPa.s-1
rho_O4 = 826; % Oil Density in kg/m3
rho_B4 = 985; % Brine Density in kg/m3
IFT4 = 1; % Interfacial Tension mN/m
nd = 1000; % No. of Droplets
Vol = 900; % Liquid volume of emulsion (ml)
l = 0.5; % Mean Distance between droplets
alpha = 0.08; % Empirical Collision Effiency Parameter
D0 = 300; % Initial Droplet Diameter (microns)
%% Sedimentation interface, hs: dhs/dt
Pr0 = ((nd*pi*D0^3)/6)/Vol % Initial Volume Fraction of droplet
Prm = ((nd*pi*((D0 + l)^3))/6)/Vol % Maximum Volume Fraction of droplet
delrho1 = rho_B2 - rho_O2 % difference between the dispersed water and continuous oil phase
Vsto = (delrho1*g*(D0^2))/18*mu2 % Settling Velocity of Hard Spheres (stoke's velocity)
fPr = (1 - Pr0)^5.3 % Dimensionless
x1 = 60 % (hs - hd) in mm, guessed value
K1 = ((2/3)*alpha*((Vsto^2)/D0))*((fPr^2)/((Prm/fPr)^(1/3)) - 1)
More Answers (1)
Steven Lord
on 29 Jun 2022
You probably want to recheck your units, do some dimensional analysis. Note that at the top of the Y axis that 1 unit on the Y axis represents 1e124 mm, and your X axis spans 180 seconds.
Google says that means whatever this represents is traveling around 2e119 kilometers per hour -- 1e110 c. This thing takes off not like a rocket, not like the Starship Enterprise at maximum warp, but more like the TARDIS from Doctor Who.
At a quick glance, should your initial condition be in millimeters or meters? Most of the physical parameters are in units of kg/m^3.
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!