Error using odearguments?

Hi, I tried to run the code below but I was greeted with the below error message. Any assistance will be greatly appreciated

function f = emu(t, nd, alpha, mu, D0, Vol, g, L)

%% Oil/Water System at 30C (Paraffin)

% Concentration for this system is 200ppm

lambda = 1;

theta_0 = nd*pi*D0^3 /(6*Vol);

theta_m = nd*pi*(D0+L)^3 /(6*Vol);

ftheta_0 = (1 - theta_0)^5.3;

Delta_rho = 9.98*10^-9 - 7.97*10^-9;

Vst0 = Delta_rho*g*D0^2 /(18*mu); % Settling Velocity in mm/s

D = sqrt((abs((2/3)*alpha*Vst0*ftheta_0/((theta_m/theta_0)^(1/3)-1)*D0*t +D0^2)));

K1 = abs((2/3)*alpha*Vst0^2 *ftheta_0^2 /(D*((theta_m/theta_0)^(1/3)-1)));

dhs = lambda*(K1*t + Vst0*ftheta_0);

dD = lambda*alpha/3 *(Vst0*ftheta_0)/((theta_m/theta_0)^(1/3) - 1) *D0/D;

f = [dhs;dD];

end

USing the following to generate the figure:

clc, clear all, close all

%% Oil/Water System at 30C (Paraffin)

Vol = 900; % Volume for emulsion injected

g = 9810; % mm/s^2 - gravity

L = 0.5e-6; % distance between the surfaces of neighboring bubbles (guessed value)

nd = 1000; % guessed value for no. of bubbles

mu = 3.8; % Viscosity in mPa.s-1

alpha = 0.08;

D0 = 0.3; % converted 300umVst0 = Delta_rho*g*D0^2 /(18*mu);

tspan = [0 180];

hs0 = 225; % in millimeters - initial height

figure

[t,dhs] = ode45(@(t,dhs)emu(t, nd, alpha, mu, D0, Vol, g, L), tspan, hs0);

hs = x(:,1);

plot (t, dhs(:,1), '-o')

grid on

ylabel ('H(mm)')

xlabel ('Time (s)')

title ('Simulation of Sedimentation-Coalesence Experiment')

Unfortunately, I got the following error message:

Error using odearguments

@(T,DHS)EMU(T,ND,ALPHA,MU,D0,VOL,G,L) returns a vector of length 2, but the length of initial conditions vector is 1. The vector returned

by @(T,DHS)EMU(T,ND,ALPHA,MU,D0,VOL,G,L) and the initial conditions vector must have the same number of elements.

Error in ode45 (line 107)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in fig (line 22)

[t,dhs] = ode45(@(t,dhs)emu(t, nd, alpha, mu, D0, Vol, g, L), tspan, hs0);

Any assistance please?

[Torsten]

The error message seems pretty clear.

You give an initial value hs0 = 225. This signals to ODE45 that you want to solve 1 differential equation.

But in emu, you return f = [dhs;dD] which signal to ODE45 that you want to solve 2 differential equations.

So ODE45 asks itself: Which indicator is the correct one ?

What I find very surprising in your ODE(s) is that they don't seem to depend on the solution variable, but only on time t. Is this really correct ? If yes, use "integral" instead of ODE45. It will be ways faster.

[Matthew Charles]

Hi @Torsten it largely deends on the change in droplet diameter (D) with time (t). Should i still switch to "integral" instead of the ODE45 given this case? and if so, I am not sure how to go about that? Any assistance will be appreciated

Answer 1

[Torsten]

syms t
Vol = 900;              % Volume for emulsion injected
g = 9810;               % mm/s^2 - gravity
L = 0.5e-6;             % distance between the surfaces of neighboring bubbles (guessed value)
nd = 1000;              % guessed value for no. of bubbles

mu = 3.8;                % Viscosity in mPa.s-1
alpha = 0.08;
tspan = [0 18000];
dt = 0.1;
hs0 =  225; % in millimeters - initial height
D0 = 0.3; % converted 300umVst0 = Delta_rho*g*D0^2 /(18*mu);

lambda = 1;

theta_0 = nd*pi*D0^3 /(6*Vol);
theta_m = nd*pi*(D0+L)^3 /(6*Vol);

ftheta_0 = (1 - theta_0)^5.3;
Delta_rho = 9.98*10^-9 - 7.97*10^-9;
Vst0 = Delta_rho*g*D0^2 /(18*mu); % Settling Velocity in mm/s
D = sqrt((abs((2/3)*alpha*Vst0*ftheta_0/((theta_m/theta_0)^(1/3)-1)*D0*t +D0^2)));
K1 = abs((2/3)*alpha*Vst0^2 *ftheta_0^2 /(D*((theta_m/theta_0)^(1/3)-1)));
dhs = lambda*(K1*t + Vst0*ftheta_0);
dD = lambda*alpha/3 *(Vst0*ftheta_0)/((theta_m/theta_0)^(1/3) - 1) *D0/D;

% Integrate dhs
hs = int(dhs,t);
hs = hs0 + hs - subs(hs,t,0);
hs = matlabFunction(hs);
% Use D directly
D = matlabFunction(D);
% Integrate dD to get back D
DD = int(dD,t);
DD = D0 + DD - subs(DD,t,0); 
DD = matlabFunction(DD);
% Plot quantities
t = tspan(1):dt:tspan(2)

t = 1×180001

0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000 2.1000 2.2000 2.3000 2.4000 2.5000 2.6000 2.7000 2.8000 2.9000

figure(1)
plot(t,hs(t))

% Compare D and DD
figure(2)
plot(t,D(t))
hold on
plot(t,DD(t))