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

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?

Answer 1

[Paul]

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

[Sam Chak]

Hi @Paul

OP @Matthew Charles started the question originally here:

https://www.mathworks.com/matlabcentral/answers/1750120-too-many-input-arguments-i-am-new-to-matlab-and-not-sure-how-to-solve-this

I also had the same suspicion on K1 because it returns a complex-valued number. Since ,

dhdt = x1*(K1*t - Vsto*fPr)

and K1 has an unbelievably large positive real part, then

will cause the height of the column h to explode.

[Matthew Charles]

HI @Sam Chak & @Paul, I mistakingly forgot the brackets for K1 as its raised to the 1/3 power indeed. However, wrt 1-Pro, that value is squared which, I presume, will give an even greater value. However, even with these changes the resulting graph still explodes. Will review once more

[Paul]

1 - Pr0 is not squared. fPr is squared and fPr is complex number. Is fPr supposed to be complex? If not, there's a problem upstream. If so, maybe some other operation was intended, like abs(fPr)^2 perhaps. But then fPr is used elsewhere in so there will stll be a problem with fPr being complex-valued.

[Dennis]

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]

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

Answer 2

[Steven Lord]

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.