I'm having issues to get expected results from this script. It is simulating a bubble rising through water using ode45. The expected result is for velocity is to gradually change to a stable point. However, the velocity instantly changes to the stable point which makes no sense. I tried increasing the iterations and debugging but wanted to be certain the problem isn't with the code before moving to the physics side.
The code:
clc
for i = 1:1
velocity = SimData(i,1);
velocity = table2array(velocity);
sparger_dynamics_distance(velocity)
end
function sparger_dynamics_distance(v_int)
tspan = linspace(0,15,500); % Time span for the simulation (from 0 to 10 seconds)
v0 = v_int;
x0 = 0;
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-12);
[t,x] = ode45(@ode_function_x, tspan, [x0;v0], options);
matrix2=[t, x];
% Plot the velocity as a function of time
hold on
plot(t, x)
xlabel('Time (s)');
ylabel('Displacment (m)');
title('Time vs. Displacment');
grid on;
[t, x]
function dx = ode_function_x(t, x)
dx = zeros(2,1);
dx(1) = x(2); % dx(1) is the velocity
density_air = 1.225; % Density of air (kg/m^3)
density_fluid = 1000; % Density of fluid (kg/m^3)
bubble_radius = 0.0005; % Bubble radius (m)
hole_diameter = bubble_radius*2;
hole_area = ((pi/4)*hole_diameter^2); % Area of the hole
V_B = (4/3)*pi*bubble_radius^3; % Volume of the body (m^3)
g = 9.81; % Gravitational acceleration (m/s^2)
m_B = density_air * V_B; % Mass of the body (kg)
visc_W = 1.002*(10^-3); % Dynamic viscosity of the fluid (N*s/m^2)
% Calculate Reynolds number (Re)
Re = (density_fluid * dx(1) * (bubble_radius*2)) / visc_W;
% Calculate drag coefficient (C_D) using the given formula
C_D = (24/Re) + (2.6 * (Re/5) / (1 + (Re/5)^1.52)) + (0.41 * (Re/263000)^-7.94 / (1 + (Re/263000)^-8)) + (Re^0.8 / 461000);
% Calculate forces
F_B = density_fluid * V_B * g;
F_D = (0.5 * density_fluid * hole_area * C_D) * (dx(1))^2;
F_W = m_B * g;
% Calculate acceleration (dv/dt)
dx(2) = (((F_B - F_D - F_W))/m_B);
end
end
If anyone knows where the problem might be, I would greatly appreciate it.
The results I get:
https://docs.google.com/spreadsheets/d/1ejg9SU9IcGXnRQUKNouk9oPKXKuuo0VU9Ur-XbrMQAc/edit?usp=sharing
Tried increasing ode45's tolerance,increasing iterations, and debugging but nothing worked. I'm expecting velocity/dv(1) to gradually decrease and not instantly decrease.