Warning: Imaginary parts of complex X and/or Y arguments ignored.

3 次查看（过去 30 天）
Matthew Charles
Matthew Charles 2022-6-29
评论： Sam Chak 2022-6-29
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.
I set the initial value as 240 but even that is not reflected on the graph, Can you assist please?

请先登录，再回答此问题。

采纳的回答

Paul
Paul 2022-6-29
编辑：Paul 2022-6-29
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 个评论
显示 3更早的评论
Dennis
Dennis 2022-6-29
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
Sam Chak 2022-6-29
Ran in:
Hi @Matthew Charles
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 .
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
Pr0 = 1.5708e+07
Prm = ((nd*pi*((D0 + l)^3))/6)/Vol % Maximum Volume Fraction of droplet
Prm = 1.5787e+07
delrho1 = rho_B2 - rho_O2 % difference between the dispersed water and continuous oil phase
delrho1 = 161
Vsto = (delrho1*g*(D0^2))/18*mu2 % Settling Velocity of Hard Spheres (stoke's velocity)
Vsto = 39485250
fPr = (1 - Pr0)^5.3 % Dimensionless
fPr = -8.1032e+37 - 1.1153e+38i
x1 = 60 % (hs - hd) in mm, guessed value
x1 = 60
K1 = ((2/3)*alpha*((Vsto^2)/D0))*((fPr^2)/((Prm/fPr)^(1/3)) - 1)
K1 = 4.4121e+97 + 9.9098e+97i

请先登录，再进行评论。

更多回答（1 个）

Steven Lord
Steven Lord 2022-6-29
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.

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by