Error using odearguments?

2 次查看（过去 30 天）
Matthew Charles
Matthew Charles 2022-7-12
编辑： Torsten 2022-7-12
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?
  2 个评论
Torsten
Torsten 2022-7-12
编辑：Torsten 2022-7-12
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
Matthew Charles 2022-7-12
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

请先登录，再进行评论。

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

回答（1 个）

Torsten
Torsten 2022-7-12
编辑：Torsten 2022-7-12
Ran in:
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))

类别

在 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