Submarine Explosion
-----------------------------
A large mass of incompressible, inviscid fluid contains a spherical bubble obeying Boyle's Law:
p V = constant
At great distances from the bubble, the pressure is zero.
Neglecting body forces, show that the radius R(t) of the bubble at time t satisfies the equation:
(d/dt) (R^2 (d/dt)(R) ) - (1/2) R ((d/dt)(R))^2 = k/R^2
for a constant k
Solution 1.
------------
Here spherical symmetry applies and so:
phi(r, t) = (F(t)/r) + G(t)
Then we consider our boundary condition.
A unit normal to the boundary between the bubble and the fluid is e_r, and so:
grad(phi) dot n = (d/d r)(phi)
The expanding bubble represents a moving boundary, given by R = R(t), with velocity:
U = (d/dt)( R(t)) e_r
Applying the boundary condition at r=R, we have:
u(R,t) = (d/dt)(phi(R,t)) = - F(t)/R(t)^2 = (d/dt)(R(t))
so that F(t) = - R^2 (d/dt)(R) and
phi(r, t) = - (R(t)^2)/r (d/dt)(R(t)) + G(t), v(r, t) = R(t)^2/(r^2) (d/dt)(R(t)) e_r
Bernoulli's equation for unsteady irrotational flow under zero body forces takes the form:
p(r, t)/rho + (1/2)( R(t)^4 / r^4) (d/dt)(R(t)) - (1/r)(d/dt)( R(t)^2) (d/dt) (R(t)) + G'(t) = H(t)
Writing this as r goes to infinity, where the fluid is at zero pressure, leads to:
G'(t) = H(t).
Thus Bernoulli's equation reduces to:
p(r, t) / rho + (1/2) (R(t)^4 / r^4) ((d/dt)(R(t)))^2 - (1/r) (d/dt)( R(t)^2 (d/dt)(R(t)) ) = 0,
to hold everwhere.
In particular when r == R, i.e. at the boundary, we have:
p(R, t) / rho + (1/2)*((d/dt)(R(t)))^2 - (1/R(t)) (d/dt)(R(t)^2 (d/dt)(R(t)) ) = 0
Now at r == R, Boyle's law holds: pV = C, a constant:
p(R, t) = C/( (4/3) pi R^3 ),
so that
p(R, t)/ rho = k/ R^3,
for some constant k, and:
k / R(t)^3 = (1/R(t)) (d/dt) (R(t)^2 (d/dt)(R(t))) - (1/2) ( (d/dt)(R(t)) )^2.
Multiplying across by R, we obtain the required nonlinear equation, which we may solve numerically.
equation1 := (3/2) R(t) ( (d/dt)(R(t)) )^2 + R(t)^2 ((d/dt)^2(R(t)) ) - 10/R(t)^2 = 0
-----------------------
Maple code:
-----------------------
> restart: with(plots) : with(Detools):
k := 10
eq1 := diff( R(t)^2* diff(R(t), t) ,t) - 1/2*R(t)* (diff(R(t),t))^2 - k/R(t)^2 = 0;
p := dsolve( {eq1, R(0)=1, D(R)(0) = 0, R(t), type=numeric, range = 0..10} ):
odeplot(p);