Unexpected results from matlab script (using ode45)

39 views Asked by At

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.

0

There are 0 answers