I am having a hard time converting fluid dynamics equation into MATLAB code

Question

I converted these two equations into MATLAB code before but I lost the file and now I am trying to write them again but there is a problem and I can't figure it out.

This is the main equation:

main equation

This is the equilibrium radius equation:

equilibrium radius equation

MATLAB figure

This is the figure that I got from the MATLAB code that I lost.

This is the MATLAB code that I am writing now:

clc

%Converting to SI units

Contact_angle =     142*pi/180; 
Volume =            5* 1e-9; 
ST_L =              0.072 ; 
Density =           997 ; 
Gravity =           9.807; 
Shape_factor =      37.1; 
T_zero =            0; Dynamic_viscosity = 8.9e-4 ;

R_e = ((4*Volume)/(pi*Contact_angle))^(1/3);


T = 0:0.1:5.2; 

R = R_e * (1 - exp(((-(2*ST_L)/R_e^12) + ((Density*Gravity)/(R_e^10))) * ((24*Shape_factor*Volume^4 * (T+T_zero))/((pi^2)*Dynamic_viscosity)))).^(1/6);


%Convert R to mm R = R*1000;


plot(T,R); 
grid on  xlim([0 6]); 
ylim([0 6]); 
legend('Water','Location','northwest'); 
ylabel('Radius (mm)'); 
xlabel("Time (s)");

I tried playing around with the shape factor (because it's an arbitrary constant) and playing with the units to get similar results to my previous code. However, when I tried to use a different fluid to test my mathematical model, I realized that changing units or constants is not a proper solution.

please remember to format your code so that folks can easily read it. — Mike 'Pomax' Kamermans, Apr 11 '23 at 15:35
@Mike'Pomax'Kamermans That was probably my fault. I didn't notice that the newlines weren't retained when adding the code fences. I'll fix that. — beaker, Apr 12 '23 at 15:44

score 1 · Answer 1 · answered Apr 13 '23 at 17:46

I get to the same steady result like duffymo :

1.- with the r_e approximation

theta_e =     142*pi/180; 
V =            5* 1e-9; 
gamma_L =              0.072 ; 
rho =           997 ; 
g =           9.807; 
lambda =      37.1; 
t_zero =    0; 
etha = 8.9e-4 ;

% r_e=(6*V*(cos(theta_e/2))^3/(pi*(2+cos(theta_e))*sin(theta_e/2)))^(1/3)
r_e = (4*V/(pi*theta_e))^(1/3);

t = [0:0.1:5]; 

R =r_e*(1 - exp(-(2*gamma_L/r_e^12 + rho*g/(9*r_e^10)) *...
    24*lambda*V^4*(t+t_zero)/(pi^2*etha))).^(1/6)

R =
  Columns 1 through 4
         0    0.0014    0.0014    0.0014
  Columns 5 through 8
    0.0014    0.0014    0.0014    0.0014

...

  Columns 45 through 48
    0.0014    0.0014    0.0014    0.0014
  Columns 49 through 51
    0.0014    0.0014    0.0014

2.- With the exact r_e expression

theta_e =     142*pi/180; 
V =            5* 1e-9; 
gamma_L =              0.072 ; 
rho =           997 ; 
g =           9.807; 
lambda =      37.1; 
t_zero =    0; 
etha = 8.9e-4 ;

r_e=(6*V*(cos(theta_e/2))^3/(pi*(2+cos(theta_e))*sin(theta_e/2)))^(1/3)
% r_e = (4*V/(pi*theta_e))^(1/3);

t = [0:0.1:5]; 

R =r_e*(1 - exp(-(2*gamma_L/r_e^12 + rho*g/(9*r_e^10)) *...
    24*lambda*V^4*(t+t_zero)/(pi^2*etha))).^(1/6)

R =
   1.0e-03 *
  Columns 1 through 4
         0    0.6600    0.6600    0.6600
  Columns 5 through 8
    0.6600    0.6600    0.6600    0.6600

...

  Columns 45 through 48
    0.6600    0.6600    0.6600    0.6600
  Columns 49 through 51
    0.6600    0.6600    0.6600

3.- Changes that bring R to linear change

R =1e6* r_e*(1 - exp(-(2*gamma_L/r_e^12 + rho*g/(9*r_e^10)) *...
    24*lambda*V^4*(1e-12)*(t+t_zero)/(pi^2*etha))) % .^(1/6)

%Convert R to mm R = R*1000;

plot(t,R); 
grid on
ylim([0 6]);xlim([0 6]);  
legend('Water','Location','northwest'); 
ylabel('Radius (mm)'); 
xlabel('Time (s)')

But this R is linear, not exponential, and when doubling t range, R keeps increasing, yet it seems to be the case that it should eventually reach a 'charged capacitor' state.

In my opinion there may be one or more values in the input parameters that may need to be checked in case should be different.

1e6 next to front r_e implies r_e value may need further assessment.

And the 1e-12 in the exponential means that the overall expression in the question pushes it all far away from where the 'capacitor charging' slope takes place.

Hope it helps.

score -1 · Answer 2 · answered Apr 11 '23 at 20:07

I tried reproducing your code using Python and displaying it with matplotlib, without success.

The values in the y array are all O(1.e-30). The resulting values in the r array jump from zero to re = 1.3695297485934952 and stay there, without changing. I'm not sure what I did wrong.

# see https://stackoverflow.com/questions/75985316/i-am-having-a-hard-time-converting-fluid-dynamics-equation-into-matlab-code

import math
import matplotlib.pyplot as plt
import numpy as np

if __name__ == '__main__':
    theta = 142*math.pi/180
    v = 5e-9
    gamma_L = 0.072
    rho = 997
    g = 9.807
    lamb = 37.1
    t0 = 0
    mu = 8.9e-4
    t = np.arange(0.0, 2.0, 0.01)
    re = (4*v/math.pi/theta)**(1/3)
    x = 2*gamma_L/(re**12) + rho*g/9.0/(re**10)
    y = [(24*lamb*(v**4)*(ti+t0)/(math.pi**2)/mu) for ti in t]
    r = [(1000.0 * re * (1 - math.exp(-x*yi))**(1/6))for yi in y]

    fig, ax = plt.subplots()
    ax.plot(t, r, linewidth=2.0)
    ax.set(xlim=(0, 5),
           ylim=(0, max(r)))
    plt.show()

thank you all for the help. But I figured it out – shammas mohamed Apr 16 '23 at 20:08 — shammas mohamed, Apr 16 '23 at 20:08
Tell us what you "figured out". I'd like to know. – duffymo Apr 16 '23 at 23:48 — duffymo, Apr 16 '23 at 23:48

I am having a hard time converting fluid dynamics equation into MATLAB code

2 Answers2