# Initial Conditions not Generating on Graph

1 view (last 30 days)
Matthew Charles on 5 Jul 2022
Edited: Matthew Charles on 25 Aug 2022
=
B = 1e-19;
% Vol = nd*(4*pi*(D/2)^3 /3);
theta_0 = nd*pi*D^3 /(6*Vol);
theta_m = nd*pi*(D+1)^3 /(6*Vol);
ftheta_0 = (1 - theta_0)^5.3;
delta_c = 0.567*(B*D^2 /sigma)^(1/4);
beta = 32*delta_c^2 /(9*D);
K3 = 180*beta*theta_d^2 /(1-theta_d)^3;
% Solve for Delta_p (pressure change) and dhw/dt
x0 = [1,1];
% optio
y = fsolve(@(y)[y(1) - 180*y(2)*mu*theta_d^2 *(hd-hw)/(D^2 *(1-theta_d)^3);-lambda*(beta*y(1)*g*theta_d/mu)*(hd-hw)/(1 + K3*(hd-hw)/D^2)], x0);
Delta_p = y(1);
dhw = y(2);
% Delta_p = 994-833;
Vst0 = Delta_p*g*D0^2 /(18*mu);
K1 = abs((2/3)*alpha*Vst0^2 *ftheta_0^2 /(D0 * ((theta_m/theta_0)^(1/3)-1)));
K2 = beta*Delta_p*g*theta_d/mu;
dhs = lambda*(K1*t + Vst0*ftheta_0);
% dhw = lambda*K2*(hd - hw)/(1 + (K3/D^2)*(hd - hw));
dhd = -theta_0/(theta_d - theta_0) *dhs - (1 - theta_d)/(theta_d - theta_0) *dhw;
dD = lambda*alpha/3 *(Vst0*ftheta_0)/((theta_m/theta_0)^(1/3) - 1) *D0/D;
f = [dhs;dhw;dhd;dD];
end
Main function shown above
clc, clear all, close all
alpha = 0.8;
nd = 900;
sigma = 5e-3; % Interfacial tension N/m
rho = 833; % density kg/m^3
mu = 5e-3; % Viscosity mPa/s
theta_d = 0.74;
Vol = 1e-6; % Volume
g = 9.81; % Gravity
D0 = 300e-6;
H = 225e-3;
tspan = 0:4:3600;
hw0 = 0;
hd0 = 0;
hs0 = H;
code for graph to be generated.
As can be seen in the code above, the initial conditons declared did not generate on the graph, only hs0 = H, which I am unsure about. Any assistance will be greatly appreciated. thanks

Torsten on 5 Jul 2022
You plot "hs" for different values of "nd". So the initial condition for hs is always the same, namely 225 mm.
I don't understand your problem with the graph you gry to generate.
alpha = 0.8;
nd = 900;
sigma = 5e-3; % Interfacial tension N/m
rho = 833; % density kg/m^3
mu = 5e-3; % Viscosity mPa/s
theta_d = 0.74;
Vol = 1e-6; % Volume
g = 9.81; % Gravity
D0 = 300e-6;
H = 225e-3;
tspan = 0:4:3600;
hw0 = 0;
hd0 = 0;
hs0 = H;
figure
hold on
nd_vec = [100 150 180 200 220 240 250 260 270 300 600 900];
for i = 1:length(nd_vec)
nd = nd_vec(i);
x0 = [hs0, hw0, hd0, D0];
[t,x] = ode45(@(t,x)df(t, x, nd, alpha, mu, sigma, D0, theta_d, Vol, g), tspan, x0);
% figure
plot(t, x(:,1)*1e3)
end
title('{h_{s}(t)}')
ylabel('mm')
xlabel('Time')
function f = df(t, x, nd, alpha, mu, sigma, D0, theta_d, Vol, g)
hs = x(1);
hw = x(2);
hd = x(3);
D = x(4);
lambda = 1;
xx = (hd-hw);
if xx <= 0
lambda = 0;
end
B = 1e-19;
% Vol = nd*(4*pi*(D/2)^3 /3);
theta_0 = nd*pi*D^3 /(6*Vol);
theta_m = nd*pi*(D+1)^3 /(6*Vol);
ftheta_0 = (1 - theta_0)^5.3;
delta_c = 0.567*(B*D^2 /sigma)^(1/4);
beta = 32*delta_c^2 /(9*D);
K3 = 180*beta*theta_d^2 /(1-theta_d)^3;
% Solve for Delta_p (pressure change) and dhw/dt
x0 = [1,1];
% optio
y = fsolve(@(y)[y(1) - 180*y(2)*mu*theta_d^2 *(hd-hw)/(D^2 *(1-theta_d)^3);-lambda*(beta*y(1)*g*theta_d/mu)*(hd-hw)/(1 + K3*(hd-hw)/D^2)], x0);
Delta_p = y(1);
dhw = y(2);
% Delta_p = 994-833;
Vst0 = Delta_p*g*D0^2 /(18*mu);
K1 = abs((2/3)*alpha*Vst0^2 *ftheta_0^2 /(D0 * ((theta_m/theta_0)^(1/3)-1)));
K2 = beta*Delta_p*g*theta_d/mu;
dhs = lambda*(K1*t + Vst0*ftheta_0);
% dhw = lambda*K2*(hd - hw)/(1 + (K3/D^2)*(hd - hw));
dhd = -theta_0/(theta_d - theta_0) *dhs - (1 - theta_d)/(theta_d - theta_0) *dhw;
dD = lambda*alpha/3 *(Vst0*ftheta_0)/((theta_m/theta_0)^(1/3) - 1) *D0/D;
f = [dhs;dhw;dhd;dD];
end

### Categories

Find more on Networks in Help Center and File Exchange

R2022a

### Community Treasure Hunt

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

Start Hunting!