ODE45 solver not using the previous calculated value in the new time step

Hello,

I am trying to solve the equation of motion for a single air bubble rising in water. When using the ode45 solver on MATLAB, the initial conditions (zero velocity) produces v=0 at t=0. However, at the next time step (t=0.001), the forces are balanced and the bubble has already reached terminal velocity (which is unrealistic). It is shooting up immediatly to the correct terminal velocity, but I am expected a slope as the drag force should gradually increase (As velocity is increasing) until it is balanced. F_bouyancy = F_gravity + F_drag.

I am unsure of where I went wrong in my code and any help will be greatly appreciated.

function [t, x, y, vx, vy] = ode_projectile_with_air_drag_model()
    %% Solves the equation of motions for a an air bubble rising in water
    % r''= F = Fb + Fg+ Fd  
    % where
    % r is the radius vector, Fb is the bouyancy force, Fg is the gravity pull force, and Fd is the air drag force.
    % The above equation can be decomposed to x and y projections
    % x'' = Fd_x (which is going to be zero)
    % y'' = Fb  -Fg + Fd_y
    % Fd =  1/2 * rho * v^2 * Cd * A is the drag force magnitude
    % where v is speed.
    % The drag force directed against the velocity vector
    % Fd_x= - Fd * v_x/v  ; % vx/v takes care of the proper sign of the drag projection
    % Fd_y= - Fd * v_y/v  ; % vy/v takes care of the proper sign of the drag projection
    % where vx and vy are the velocity projections

    % at the first look it does not look like ODE but since x and y depends only on t 
    % it is actually a system of ODEs

    % transform system to the canonical form
    % x  -> y1
    % vx -> y2
    % y  -> y3
    % vy -> y4
    % t  -> x
    %
    % f1 -> y2
    % f2 -> Fd_x
    % f3 -> y4
    % f4 -> -g + Fd_y

    % some constants 
    rho_a=1;  % the density of air kg/m^3
    rho_w=998; % the density of water kg/m^3
    dp = 3.5 * (10^-3);   % diameter of the air bubble [mm] 
    Cd=.6;    % an arbitrary choice of the drag coefficient
    m=rho_a*(4/3)*pi*(dp/2)^3;   % the mass of the air bubble  in kg
    g=9.81;    % the acceleration due to gravity
    A=(pi/4)*(dp^2); % the area of the air bubble in m^2, (Pi*r^2)
    mu_w = 1e-3;         % viscosity of water [Pa*s]
    vol_g = (4/3)*pi*(dp/2)^3; % volume of air bubble
    

    function fvec = projectile_forces(x,y)
        % it is crucial to move from the ODE notation to the human notation
        vx=y(2);
        vy=y(4);
        v=sqrt(vx^2+vy^2); % the speed value

        Fg = m*g;
        Fd= 1/2 * rho_w * (v^2) * Cd * A;
        Fb = (vol_g)*(rho_w)*g; %(m/rho_a)*rho_w*g;
     
        if v == 0
          fvec(1) = y(2);
          fvec(2) = 0
          fvec(3) = y(4);
          fvec(4) = ((Fb - Fg -Fd)/m);
      else
        fvec(1) = y(2);
        fvec(2) = -(Fd*vx/v/m);
        fvec(3) = y(4);
        fvec(4) = Fb/m - Fg/m - Fd*vy/v/m ;
      
      end
        % To make matlab happy we need to return a column vector.
        % So, we transpose (note the dot in .')
        fvec=fvec.'; 
    end

    %% Problem parameters setup:
    %  We will set initial conditions to be zero, bubble released from rest
    tspan=[0, 1];       % time interval of interest
    y0(1)=0;             % the initial x position
    y0(2)=0;             %v0*cos(theta); % the initial vx velocity projection
    y0(3)=0;             % the initial y position
    y0(4)=0;             %:v0*sin(theta); % the initial vy velocity projection

    % We are using matlab solver
    [t,ysol] = ode45(@projectile_forces, tspan, y0);

    % Assigning the human readable variable names
    x =  ysol(:,1);
    vx = ysol(:,2);
    y =  ysol(:,3);
    vy = ysol(:,4);
    v=sqrt(vx.^2+vy.^2); % speed

   
    ax(1)=subplot(2,1,1);
    plot(x,y, 'r-');
    set(gca,'fontsize',14);
    xlabel('Position x component, m');
    ylabel('Position y component, m');
    title ('Trajectory');
    

    ax(2)=subplot(2,1,2);
    plot(t,v*1000, 'r-');
    set(gca,'fontsize',14);
    xlabel('Time, s');
    ylabel('Velocity, mm/s');
    title ('Speed vs. time');
  

    %linkaxes(ax,'x'); % very handy for related subplots
end

Answer 1

[Saarthak Gupta]

Hi Omar,

Looks like you are getting rapid attainment of terminal velocity for an air bubble rising in water.

As per my understanding, the system is modelled correctly, and the equations of motion are consistent with literature. It is unlikely that the ODE solver is behaving incorrectly, and the solution should be expected.

The time it takes for a bubble to achieve terminal velocity depends on fluid properties, such as density and viscosity, as well as bubble characteristics like size and shape. Initial conditions, including the bubble's starting velocity, external forces, and the drag coefficient, also play crucial roles in determining the rate at which the bubble accelerates towards terminal velocity.

For your system, the bubble achieves a terminal velocity of ~ 0.276 m/s around the t = 5.7e-05 second mark.

Refer to the following code. It solves for the velocity of the air bubble, which is a simpler version of your problem since it involves a first-order differential equation:

bubble_motion(3.5e-3);

function bubble_motion(bubbleDiameter)
r = bubbleDiameter/2; % Radius of the bubble in meters
rho_water = 1000; % Density of water in kg/m^3
rho_air = 1; % Density of air in kg/m^3 (at sea level)
g = 9.81; % Acceleration due to gravity in m/s^2
Cd = 0.6; % Drag coefficient for a sphere

m = (4/3) * pi * r^3 * rho_air; % Mass of the air bubble
v0 = 0; % Initial velocity of the bubble
tspan = [0 10]; % 10 seconds simulation

% Solve the ODE using ode45
[t, v] = ode45(@(t, v) dvdt(t, v, r, rho_water, g, Cd, m), tspan, v0);

% Plot the results
plot(t, v);
xlabel('Time (s)');
ylabel('Velocity (m/s)');
xlim([0 1e-4]);
title(sprintf('Velocity of Rising Air Bubble with radius %.5f m\n', r));
end

function dv = dvdt(t, v, r, rho_water, g, Cd, m)
% Drag force using drag coefficient
A = pi * r^2; % Cross-sectional area of the bubble
F_drag = (1/2) * Cd * rho_water * A * v^2;

% Buoyant force
F_buoy = (4/3) * pi * r^3 * rho_water * g;

% Gravitational force on the bubble
F_gravity = m * g;

% Net force
F_net = F_buoy - F_gravity - F_drag;

% Acceleration (dv/dt)
dv = F_net / m;
end

Setting the x-axis limits to [0 1e-4] reveals the downwards concave nature of the velocity profile.

Refer to the following MATLAB documentation for further reference:

• ode45: https://www.mathworks.com/help/matlab/ref/xlim.html
• xlim: https://www.mathworks.com/help/matlab/ref/ode45.html