3
$\begingroup$

I have to solve $$\left\{\begin{matrix} x\ddot{x}+\frac{3}{2}\dot{x}^2=\frac{1}{\rho}\left (P_v-a_1t^2-a_2t-P_0 \right )\\ x(0)=5\cdot10^{-6}\\ \dot{x}(0)=0 \end{matrix}\right.$$ I do not know if it may be helpful, but note that $a_1=2.488125\cdot 10^{14}$ is very big while the timespan is very small ($t_f=70\cdot 10^{-6}$). I get the result without MATLAB returning any error or warning, but this result does not convince me (I expected the peak to be around $t=4\cdot 10^{-5}$). I have already tried switching to ode15s, but the output is the same. What else can I do to see if the output changes or not and then be sure that this result is correct?

function vdot = no_thermal_effect_98(t,v)
rho = 959.78;
P0 = 101325;
Pvap = 94285.313;
a1 = (P0 - 1800) / (20e-6)^2;    %nondimensional
a2 = 2 * (1800 - P0) / (20e-6);  %nondimensional
%                                       
vdot = zeros(2,1);
vdot(1) = v(2);
vdot(2) = -1.5 * v(2) * v(2) / v(1) + 1 / (v(1) * rho) *...
    (Pvap - P0 - a1 * t^2 - a2 * t);

run

x0 = 5e-6; %meters
tf = 70e-6; %seconds
[t,v] = ode45(@no_thermal_effect_98,[0,tf],[x0,0]);
[t,v(:,1)];
plot(t,v(:,1))

enter image description here

Edit

Let us first consider the case where the pressure outside the bubble is constant: $x\ddot{x}+\frac{3}{2}\dot{x}^2=\frac{1}{\rho}\left (P_v-P_0\right )$.

If $P_0>P_v$ the bubble collapses (i.e. the radius decreases in a monotonic fashion); conversely, if $P_0<P_v$ the bubble grows.

Now let us consider a pressure variation outside the bubble described by the equation $P(t)=a_1t^2+a_2t+P_0$.

enter image description here

At $40 \mu s$ the pressure outside the bubble is greater than the pressure inside the bubble (here in the picture, indicated by the orange line and equal to $Pv=94285.313 Pa$). That is why I expected the peak, that is, the moment where the bubble started collpasing, would be $40 \mu s$.

Improve this question
$\endgroup$
9
  • $\begingroup$ You've used the default (not very accurate) tolerance in ode45(). What happens when you tighten the error tolerance? $\endgroup$
    – Brian Borchers
    Commented Nov 28, 2016 at 16:16
  • $\begingroup$ @BrianBorchers I tried options = odeset('RelTol',1e-13,'AbsTol',1e-13); I got the same output. $\endgroup$
    – Marco
    Commented Nov 28, 2016 at 16:20
  • 3
    $\begingroup$ Why do you expect the output to be different from what you've got? $\endgroup$
    – Brian Borchers
    Commented Nov 28, 2016 at 17:03
  • 1
    $\begingroup$ Sundials also gives the same plot. I don't think this has anything to do with ode45. $\endgroup$
    – Kirill
    Commented Nov 28, 2016 at 20:09
  • 1
    $\begingroup$ I think the comparison with the constant-pressure case is wrong: if $\dot x$ is positive and $\ddot x$ given by the r.h.s. starts decreasing ($\ddot x<0$ for $t>24\mathrm{\mu s}$ for this ODE), it will take some time before $\dot x = 0$, there's no reason it would happen immediately. I think a better comparison might be with a piecewise-constant example, $P_v(t) = \mathrm{const} > P_0$ for $t\in[0,t_1)$ and $P_v(t) = \mathrm{const} < P_0$ for $t\in [t_1,t_f]$. $\endgroup$
    – Kirill
    Commented Nov 28, 2016 at 21:06

1 Answer 1

Reset to default
8
$\begingroup$

I ran your same code in Python

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt


def func(v, t):
rho = 959.78
P0 = 101325
Pvap = 94285.313
a1 = (P0 - 1800) / (20e-6)**2
a2 = 2 * (1800 - P0) / (20e-6) 
v1 = v[1];
v2 = -1.5 * v[1] * v[1] / v[0] + 1 / (v[0] * rho) *\
(Pvap - P0 - a1 * t**2 - a2 * t)
return [v1, v2]


y0 = 5e-6
tf = 70e-6
t = np.linspace(0, tf, 1000)
y = odeint(func, [y0, 0], t)
plt.plot(1e5*t, 1e5*y[:, 0])

and I obtained the same result

enter image description here

If you are worried about the numerical values of your constants, maybe you can change the units or non-dimensionalize your equation. For example you can change the time for $\tau = t/t_0$, with $t_0 = 20\times 10^{-6}$. And use the chain rule to make the changes to the derivatives in the left—hand side to obtain something like

$$y y''+\frac{3}{2}y'^2 = \frac{1}{\rho}\left [(P_v - P_0) - (P_0 - P_1)\tau^2 - 2(P_1 - P_0)\tau\right] $$

with $y = x/t_0$, $y' = dy/d\tau$, and $P_1 = 1800$. I tried this version and gave me the same answer.

Notes

  • In both cases, I scaled the axes by $10^5$.
  • It might be better to non-dimensionalize your equation using some knowledge of the physics in your problem. You can check reference 1 to get some ideas.

enter image description here

References

  1. Kudryashov, Nikolay A., and Dmitry I. Sinelshchikov. "Analytical solutions of the Rayleigh equation for empty and gas-filled bubble." Journal of Physics A: Mathematical and Theoretical 47.40 (2014): 405202.
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.