Optimal value of amount of water in a water rocket to fly to the highest altitude

This is easily done (as long as you're willing to accept some assumptions) by numerical simulation. I have done it, and can post the Matlab code if that is helpful. Incidentally, for the parameters I chose (rocket volume 2 L, drag coefficient 0.4, initial pressure 100 psi, etc.) the maximum height of 112 m is reached for a fill volume of 0.63 L.

Some interesting/useful facts:

The thrust is twice the internal pressure times the nozzle area (yes, counterintuitively, twice).

As the air expands, its pressure drops more than it would at constant temperature due to adiabatic cooling.

Once all the water's gone, you get an additional boost as the remaining air is expelled (you can verify this empirically by launching the rocket with no water--if you haven't tried it, you'll be surprised by how high it goes). This bit requires (for me) some fudging, as the air jet, which will remain sub-sonic, loses a lot of its energy to friction.

Edit--here's the code as requested in a comment:

function h_max = apex_height(V_w)

% INPUT: V_w, volume of water (L)
% OUTPUT: h_max, maximum height achieved (m)

% ROCKET DESIGN PARAMETERS
m_r = 135;                %mass of empty rocket (g)
V_r = 2;                  %volume of rocket (L)
D = 7/8;                  %diameter of nozzle (in)
D_r = 2.55;               %diameter of rocket (in)
p = 100;                  %initial pressure (psi)

% EMPIRICAL PARAMETERS
C = .4;                   %drag coefficient
e = .3;                   %energy efficiency of air burst

% SIMULATION PARAMETERS
dt = .0001;               %thrust phase time increment (s)
dT = .01;                 %coast phase time increment (s)
N = 5000;                 %max number of time steps

% ENVIRONMENTAL PARAMETERS
atm = 101325;             %atmospheric pressure (Pa)
gamma = 1.4;              %ratio of specific heats for air
rho = 1000;               %density of water (kg/m^3)
rho_a = 1.2;              %density of air (kg/m^3)
g = 9.8;                  %gravity (m/s^2)

% DERIVED QUANTITIES
V_r = V_r/1000;           %volume of rocket (m^3)
V_w = V_w/1000;           %volume of water (m^3)
V_A = V_r - V_w;          %initial air volume (m^3)
V_a = V_A;                %instantaneous air volume (m^3)
D = D*.0254;              %diameter of nozzle (m)
D_r = D_r*.0254;          %diameter of rocket (m)
A = pi*D^2/4;             %area of nozzle (m^2)
A_r = pi*D_r^2/4;         %rocket cross section area (m^2)
P = atm*p/14.7;           %initial relative pressure (Pa)
p = atm + P;              %absolute pressure (Pa)
m_r = m_r/1000;           %mass of empty rocket (kg)
m = m_r + V_w*rho;        %mass of rocket + water (kg)
m_a = V_A*rho_a*p/atm;    %mass of air in rocket (kg)
M = m_a - V_r*rho_a;      %mass of expelled air (kg)
thr = 2*A*P;              %thrust (N)
C = C/2*A_r*rho_a;        %drag coefficient (kg/m)
v_r = 0;                  %rocket velocity (m/s)
v_w = sqrt(2*P/rho);      %exhaust (water) velocity (m/s)
time = 0;                 %time since launch (s)
I = zeros(N,1);           %initialize vectors
p = [p; I];               %rocket pressure vector
v_r = [v_r; I];           %rocket velocity vector
time = [time; I];         %time vector

% THRUST PHASE CALCULATION
i = 1;
while V_w>0         %while there's still water
    i = i+1;
    a = (thr - C*v_r(i-1)^2)/m - g; %acceleration
    time(i) = time(i-1) + dt;       %STEP time
    v_r(i) = v_r(i-1) + dt*a;       % and rocket velocity
    V_a = V_a + dt*A*v_w;           % and air volume
    V_w = V_r - V_a;                %COMPUTE water volume
    p(i) = p(1)*(V_A/V_a)^gamma;    %    and pressure
    thr = 2*A*(p(i)-atm);           %    and thrust
    v_w = sqrt(2*(p(i)-atm)/rho);   %    and water speed
    m = m_r + rho*V_w;              %    and mass
end
p = p(1:i);
h = cumsum(v_r(1:i))*dt;        %compute height
h = [h; zeros(N+1-i,1)];        %initialize height vector

% AIR BURST CALCULATION         (this part is not perfect)
E = (p(i)-atm)*V_r;             %energy stored at burnout
V = sqrt(2*e*E*M/m_r/(M+m_r));  %air burst velocity boost
v_r(i) = v_r(i) + V;

% COAST PHASE CALCULATION
while h(i)>0        %while the rocket is still in the air
    i = i+1;
    s = sign(v_r(i-1));               %up or down
    a = -C*v_r(i-1)^2/m*s - g;        %acceleration
    time(i) = time(i-1) + dT;         %STEP time
    h(i) = h(i-1) + v_r(i-1)*dT;      % and height
    v_r(i) = v_r(i-1) + dT*a;         % and velocity
end
h = h(1:i); time = time(1:i); v_r = v_r(1:i);

% % PLOT THE TRAJECTORY
% plot(time,h*3.281)
% xlabel 'time (s)', ylabel 'height (ft)'

h_max = max(h);

Accurate math/physics modeling of a water rocket is probably very, very difficult but all models have humble beginnings, can be refined in time, supplemented with empirical data and may yet provide some useful insights.

The first port of call would probably be the Rocket Equation.

Its main insight is that:

$$\Delta v=V_{exh}\ln\Big(\frac{m_0}{m_f}\Big)$$

so that for a constant exhaust velocity $V_{exh}$ and mass ratio $\frac{m_0}{m_f}$, the velocity increase $\Delta v$ is given by the above formula. This suggests that higher fill leads to higher final velocity.

But for a real water bottle, the exhaust speed will not be constant because the air pressure actually rapidly decreases. It can be (crudely) modeled as an adiabatic expansion.

Flow of water (presumed here water only) through the nozzle could be approximated with Bernoulli's Equation (ignoring the gravitational term, as well as water viscosity).

These three elements could (note caveat) be combined into a simple model for $v(t)$, as a function of initial mass of water and initial air pressure.

Because $v(t)=\frac{\rm{d}h(t)}{\rm{d}t}$, then:

$$h(\Delta t)=\int_0^{\Delta t}v(t)\rm{d}t$$

where $\Delta t$ is the total 'firing time'. To that height then has to be added the 'free flight' height to calculate the zenith: $h_z=h(\Delta t)+h_{ff}$

So far we have ignored drag, usually modeled as:

$$F_D=\frac12 \rho C_DAv^2$$

The drag coefficient $C_D$ could be strongly reduced by mounting a shallow cone on the rockets upside.

Let me try a simple model. Assume the cross-section of the bottle $A cm^2$, the mass is the mass of water $m(x) = m_0 + A x$, $m_0$ is the mass of the bottle (in CGS unit, the density of water equals to 1); and the pressure $P(x) = P_0 \frac{L-x_0}{L-x}$, where $L$ is the height of bottom, 30 cm, and $x_0$ the intial water level, the parameter of optimization. $P_0$ denotes the initial pressure of the bottle. Assume that the ejection speed a constant $u_0$. All in CGS unit, $\rho_{water} = 1$.

• the mass $m(x) = m_0 + A x$,
• pressure $P(x) = P_0 \frac{L-x_0}{L-x}$,
• ejection speed $u_0$, a constant.
• ejection speed and water level: $u_0 B = -A \frac{dx}{dt}$, where $B$ is the area of ejection nozzle (the bottle's opening).

Rocket motion from conservation of momentum (with a constant ejection speed $u_0$):

$$ \tag{1} dv = - u_0 \frac{dm}{m(x)} - g dt. $$

$dv$ is the change of rocket's speed $v(t+dt) - v(t)$. The relation between time $t$ and water level $x$ can be find from volume flux $u_0 B$:

$$ \tag{2} t = \frac{A}{u_0 B} \left( x_0 - x \right). $$

integral the above equation Eq.(1):

$$ \tag{3} v(x) = u_0 \ln \frac{m_0 + A x_0}{m_0 + A x} - \frac{g A}{u_0 B}\left( x_0 - x \right) $$

The second term is indeed $ g t $ by replacing $t$ in term of $x$.

At $x = 0$ running out of water, the speed of rocket in this moment is:

$$ \tag{4} v(x= 0) = u_0 \ln \frac{m_0 + A x_0}{m_0} - \frac{g A}{u_0 B} x_0. $$

At this moment $t(0) =\frac{A x_0}{u_0 B} $, the height of rocket $h(0)$:

$$ h(x=0) = \int_0^{t(0)} v(t) dt = \int_{x_0}^{0} v(x) \frac{dt}{dx} dx = \int_{0}^{x_0} v(x) \frac{A}{B u_0} dx $$

This integral can be done by using $\int \ln x dx = x \ln x - x$:

$$ \tag{5} h(x=0) = \frac{m_0}{B} \ln\frac{m_0 + A x_0}{m_0} - \frac{A x_0}{B} - \frac{1}{2} g \left(\frac{A x_0}{u_0 B}\right)^2. $$

The last term is $\frac{1}{2} g t^2$ in free-fall motion.

After the water run out, the rocket is still moving upward. There still a final trip to reach the top position where the speed is zero. In this moment, the rocket has speed $v(x=0)$ from Eq.(4), and may add an addition impulse from the high pressure in the bottle burst ( refer to Ben51, one of other answers to this problem,) say $v_i$. The total top distance the rocket travel is then: $$\tag{6} h_{total} = \frac{ \left( v(x=0) + v_i \right)^2}{2g} + \frac{m_0}{B} \ln\frac{m_0 + A x_0}{m_0} - \frac{A x_0}{B} - \frac{1}{2} g \left(\frac{A x_0}{u_0 B}\right)^2. $$

where $v(x=0)$ is given in Eq.(4).

Estimation of $u_0$ and $v_i$

The ejection velocity is a critical parameter. In Eq.(5) the pull-down of gravity is divided by $u_0^2$. If $u_0$ is too small, the gravity will pull down the rocket. The parameter $u_0$ has to contain the information of $x_0$, to counter balance the constant $u_0$. I suggest a reasonable estimate of $u_0$ as:

$$ P_f = P_0 \frac{L - x_0}{L}; $$ and estimate $u_0$ using this final pressure (when $x=0$) from Bernoulli Equation:

$$ u_0 = \lambda \sqrt{2 P_0 \frac{L - x_0}{L-0.5 x_0}} $$ Where $\lambda$ is an adjustable parameter to fit the real situation.

Also $v_i$ is final burst impulse, considering the work done as the pressure lower down to $1 atm$ under free expansion, I suggest:

$$ v_i = \lambda' \sqrt{\left(P_f - P_{atm}\right) * A * L / m_0} $$

Again, the parameter, $\lambda'$ is employed to fit the real situation.

An example of speed as function of time ($x_0 = 15 cm, L = 30 cm$):

And the height as function of $x_o$ using $\lambda = 0.5$ and setting $v_i = 0$: