Reparto óptimo de combustible entre las dos fases

La masa inicial m₀ es la suma de la carga útil, más el combustible y más la masa de los recipientes que contienen el combustible. Para calcular esta última cantidad, se ha supuesto que los recipientes metálicos tiene una masa que es el factor r multiplicado por la masa de combustible. Donde r es del orden del 5% ó 0.05.

masa inicial, m₀=carga útil+(1+r) · combustible total.

m₀=M+(1+r)c

La cantidad de combustible en la primera fase c₀ es igual al producto del combustible total, por el tanto por ciento p y dividido por cien.

combustible en la primera fase c₀ =combustible total · tanto por ciento/100;

c₀=p·c

Una vez que ha transcurrido un tiempo t₀ igual al cociente entre el combustible en la primera fase c₀ y la cantidad D que se quema por segundo, t₀=c₀/D se alcanza una velocidad máxima v₀

${\displaystyle v_0=u\ln \left(\frac{m_0}{m_0-c_0}\right)=u\ln \left(\frac{M+(1+r)c}{M+(1+r-p)c}\right)}$

El cohete se desprende de la primera fase disminuyendo la masa inicial del cohete m₀ en una cantidad igual a la suma de la masa del combustible quemado c₀ y la masa del recipiente que lo contiene

masa inicial al encenderse la segunda fase, m₁=carga útil+(1+r) · c₁

m₁=M+(1+r)c₁

Siendo c₁ la masa de combustible de la segunda fase, que es igual a la masa del combustible total menos la masa de combustible de la primera fase c₀ ya quemado.

combustible en la segunda fase c₁ =combustible total - combustible en la primera fase c₀

c₁=(1-p)c

En el instante t₁ se agota el combustible de la segunda fase, es igual al cociente entre la masa de combustible total y la cantidad D que se quema por segundo

t₁=combustible total/D.

Cuando se agota el combustible, el cohete alcanza la velocidad máxima v₁, continuando con la misma velocidad ya que no actúan fuerzas sobre el mismo.

${\displaystyle v_1=v_0+u\ln \left(\frac{m_1}{m_1-c_1}\right)=u\ln \left(\frac{m_0}{m_0-c_0}\right)+u\ln \left(\frac{m_1}{m_1-c_1}\right)}$

Teniendo en cuenta que

${\displaystyle \begin{array}{l}{m}_0=M+(1+r)c\\ m_1=M+(1+r)c_1=M+(1+r)(1-p)c\\ m_0-c_0=M+(1+r-p)c\\ m_1-c_1=M+r(1-p)c\end{array}}$

La velocidad final v₁ es una función de p

${\displaystyle \begin{array}{l}{v}_1=u\ln \left(\frac{m_0}{m_0-c_0}\frac{m_1}{m_1-c_1}\right)=u\ln \left(\frac{M+(1+r)c}{M+(1+r-p)c}\frac{M+(1+r)(1-p)c}{M+r(1-p)c}\right)\\ v_1=u\ln \left(\frac{m+(1+r)}{m+(1+r-p)}\frac{m+(1+r)(1-p)}{m+r(1-p)}\right),m=\frac{M}{c}\end{array}}$

Para calcular el valor máximo de la velocidad v₁, derivamos con respecto de p. De este modo determinamos la proporción óptima de combustible en las dos fases p·c en la primera y (1-p)c en la segunda

Ya que m₀ no depende de p calculamos el extremo de la función

${\displaystyle f(p)=\frac{m+(1+r)(1-p)}{\left(m+1+r-p\right)\left(m+r(1-p)\right)}}$

Obtenemos una ecuación de segundo grado en p

${\displaystyle \begin{array}{l}\frac{df(p)}{dp}=\frac{-(1+r)\left(m+1+r-p\right)\left(m+r(1-p)\right)-\left\{-\left(m+r(1-p)\right)-r\left(m+1+r-p\right)\right\}\left(m+(1+r)(1-p)\right)}{{\left(m+1+r-p\right)}^2{\left(m+r(1-p)\right)}^2}=0\\ -(1+r)\left(m+1+r-p\right)\left(m+r(1-p)\right)+\left(m(1+r)+r(2+r-2p)\right)\left(m+(1+r)(1-p)\right)=0\\ -(1+r)\left\{m^2+m\left(1+2r-p-rp\right)+r(1+r-p)(1-p)\right\}+m^2(1+r)+m\left\{{\left(1+r\right)}^2(1-p)+r(2+r-2p)\right\}+r(2+r-2p)(1+r)(1-p)=0\\ m\left\{\left(1+4r+2r^2-p-4pr-r^{2}p\right)-\left(1+2r-p-rp\right)(1+r)\right\}+r(1+r)(1-p)\left\{2+r-2p-\left(1+r-p\right)\right\}=0\\ mr(1-2p)+(1+r)(1-2p+p^2)r=0\\ p^2-2p\left(1+\frac{m}{1+r}\right)+\left(1+\frac{m}{1+r}\right)=0\end{array}}$

La solución es la raíz p<1

${\displaystyle p=k-\sqrt{k^2-k},k=1+\frac{m}{1+r}}$

u=2000; %velocidad de escape de los gases (respecto al cohete)
c=9000; %combustible
M=800;  %carga útil
r=0.05;
m=M/c;
%velocidad final del cohete
vf=@(p) u*log((m+(1+r))./(m+(1+r-p)))+u*log((m+(1+r)*(1-p))./(m+r*(1-p)));
fplot(vf,[0,1])
k=1+m/(1+r);
p_m=k-sqrt(k^2-k); %máximo
v_m=vf(p_m); %máxima velocidad
line([p_m,p_m],[4200,v_m],'lineStyle','--')
line([0,p_m],[v_m,v_m],'lineStyle','--')
disp(p_m) 
disp(vf(p_m)) %velocidad máxima
grid on
ylim([4200,4650])
xlabel('p');
ylabel('v_f(m/s)')
title('Proporción de combustible en la primera fase')

La proporción óptima es p=0.78 (78%) de combustible en la primera fase y 0.22 (22%) en la segunda. La velocidad máxima del cohete es 4 637 m/s

v_m =   4.6374e+03
>> p_m
p_m =    0.7816

Posición y velocidad del cohete en función del tiempo

Velocidad y posición del cohete en la primera etapa

Masa inicial del cohete, m₀=M+(1+r)·c.

Donde M es la carga útil y c es el combustible total

El combustible en la primera fase c₀=c·p

La velocidad y posición del cohete en el instante t

${\displaystyle \begin{array}{l}v=u\ln \frac{m_0}{m_0-Dt}\\ x=ut\ln m_0+\frac{u}{D}\left(\left(m_0-Dt\right)\ln \left(m_0-Dt\right)+Dt-m_0\ln m_0\right)\end{array}}$

El combustible de la primera fase se agota en el instante t₀=c₀/D, en ese momento la posición x₀ y velocidad v₀ del cohete en este instante son, respectivamente

${\displaystyle \begin{array}{l}{v}_0=u\ln \frac{m_0}{m_0-c_0}\\ x_0=u\frac{c_0}{D}\ln m_0+\frac{u}{D}\left(\left(m_0-c_0\right)\ln \left(m_0-c_0\right)+c_0-m_0\ln m_0\right)\end{array}}$

Velocidad y posición del cohete en la segunda etapa

Agotado el combustible de la primera fase queda c₁=c-c₀. La masa inicial del cohete es ahora m₁=M+(1+r)·c₁

La velocidad y posición del cohete en el instante t>t₀ es

$$\begin{gathered} {\displaystyle \begin{array}{l}v=v_0+u\ln \frac{m_1}{m_1-D(t-t_0)}\\ x=x_0+v_0(t-t_0)+u(t-t_0)\ln m_1+ \\ \frac{u}{D}\left(\left(m_1-D(t-t_0)\right)\ln \left(m_1-D(t-t_0)\right)+D(t-t_0)-m_1\ln m_1\right)\end{array}} \end{gathered}$$

El combustible se agota en el instante t=t₀+c₁/D=c/D

La velocidad final y la posición del cohete son, respectivamente

${\displaystyle \begin{array}{l}{v}_1=v_0+u\ln \frac{m_1}{m_1-c_1}\\ x_1=x_0+v_0\frac{c_1}{D}+u\frac{c_1}{D}\ln m_1+\frac{u}{D}\left(\left(m_1-c_1\right)\ln \left(m_1-c_1\right)+c_1-m_1\ln m_1\right)\end{array}}$

• Combustible total en el cohete c=9000 kg
• Carga útil que transporta M=800 kg
• Combustible quemado por segundo D=1000 kg/s
• Tanto por ciento de combustible en la primera fase, 70%, p=0.7
• La velocidad constante de salida de los gases es u=2000 m/s respecto del cohete
• La masa del recipiente que contiene el combustible es el 5%, r=0.05 de la masa del combustible

• Primera fase

Masa total y combustible

m₀=800+9000+0.05·9000=10250 kg,
c₀=0.7·9000=5600 kg

Velocidad final y la posición del cohete al agotarse el combustible en la primera fase son, respectivamente.

$$\begin{gathered} {\displaystyle \begin{array}{l}{v}_0=2000\ln \frac{10250}{10250-5600}=1907\text{m/s}\\ x_0=2000\frac{5600}{1000}\ln 10250+ \\ \frac{2000}{1000}\left(\left(10250-5600\right)\ln \left(10250-5600\right)+5600-10250\ln 10250\right)=5067\text{m}\end{array}} \end{gathered}$$

• Segunda fase

Estas son la velocidad inicial y la posición inicial para la segunda fase.

El combustible en la segunda fase es c₁=c-c₀=9000-5600=3400 kg
La masa inicial del cohete en esta fase es m₁=800+3400+0.05·3400=4370 kg

El tiempo que tarda en agotarse el combustible es c₁/D=3.4 s. La velocidad y posición final del cohete en el instante t=5.6+3.4=9 s es

$$\begin{gathered} {\displaystyle \begin{array}{l}v=v_0+2000\ln \frac{4370}{4370-3400}=4623\text{m/s}\\ x=x_0+v_0·3.4+2000·3.4\ln 4370+ \\ \frac{2000}{1000}\left(\left(4370-3400\right)\ln \left(4370-3400\right)+3400-4370\ln 4370\right)=13077\text{m}\end{array}} \end{gathered}$$

u=2000; %velocidad de escape de los gases (respecto al cohete)
combustible=9000; %combustible
c0=0.7*combustible;  %combustible en la primera fase
c1=combustible-c0; %combustible en la segunda fase
carga=800;  %carga útil
m0=carga+1.05*combustible; %masa inicial de la primara fase
m1=carga+1.05*c1; %masa inicial cuando empieza la segunda fase
D=1000; %kg de combustible quemado por segundo
t0=c0/D; %tiempo que tarda en agotarse el combustible de la primera fase
t1=c1/D; %tiempo que tarda en agotarse el combustible de la segunda fase
v0=u*log(m0/(m0-D*t0)); %velocidad final en la primera fase
%posición final primera fase
x0=u*t0*log(m0)+u*((m0-D*t0).*log(m0-D*t0)+D*t0-m0*log(m0))/D; 
v1=v0+u*log(m1/(m1-D*t1)); %velocidad final
%posición final
x1=x0+v0*t1+u*t1*log(m1)+u*((m1-D*t1).*log(m1-D*t1)+D*t1-m1*log(m1))/D; 
fprintf('Velocidad final %4.1f, posición final %5.1f en la 
primera fase\n',v0,x0);
fprintf('Velocidad final %4.1f, posición final %5.1f en la 
segunda fase\n',v1,x1);
 
figure
hold on
t=0:0.05:t0;
v=u*log(m0./(m0-D*t));
plot(t,v,'b');
t=0:0.05:t1;
v=v0+u*log(m1./(m1-D*t));
plot(t0,v(1),'ro','markersize',4,'markerfacecolor','r')
plot(t+t0,v,'r')
hold off
grid on
xlabel('t(s)')
ylabel('v(m/s)')
title('Velocidad')
 
figure
hold on
t=0:0.05:t0;
x=u*t*log(m0)+u*((m0-D*t).*log(m0-D*t)+D*t-m0*log(m0))/D;
plot(t,x,'b')
t=0:0.05:t1;
x=x0+v0*t+u*t*log(m1)+u*((m1-D*t).*log(m1-D*t)+D*t-m1*log(m1))/D;
plot(t0,x(1),'ro','markersize',4,'markerfacecolor','r')
plot(t+t0,x,'r')
hold off
grid on
xlabel('t(s)')
ylabel('x(m)')
title('Posición')

Velocidad final 1907.1, posición final 5066.9 en la primera fase
Velocidad final 4622.8, posición final 13077.0 en la segunda fase

Actividades

Se introduce

• El combustible total en ambas fases, en el control titulado Combustible total en el cohete
• El tanto por ciento del combustible total en la primera fase, en el control titulado Tanto por ciento de combustible en la primera fase
• La carga útil que transporta el cohete, en el control titulado Carga útil que transporta
• La cantidad D de combustible que se quema por segundo, en el control titulado Combustible quemado por seg.

Se pulsa el botón titulado Nuevo y a continuación, Empieza

Ahora se tratará de comprobar, que un cohete de dos fases que transporta la misma cantidad de combustible y la misma carga útil, es más ventajoso que el mismo cohete de una sola fase.

En segundo lugar, se tratará de investigar la dependencia de la velocidad final del cohete con el reparto de combustible total entre las dos fases. Manteniendo fijas la cantidad total de combustible y la carga útil, se tratará de modificar el tanto por ciento de combustible en la primera fase y anotar la velocidad final una vez agotado todo el combustible de la segunda fase. Comprobar la distribución óptima de combustible, aquella que da lugar a una mayor velocidad final, completando una tabla como la siguiente.

Cohete de tres etapas

Sea un cohete de tres etapas. La cantidad de combustible de la primera etapa es, c₀=p₀c, de la segunda, c₁=p₁c y de la tercera, c₂=p₂c. Donde c es la cantidad total de combustible y p_i es la proporción de combustible en cada etapa, naturalmente, p₀+p₁+p₂=1

Vamos a determinar el reparto óptimo de combustible en cada una de las tres etapas

La velocidad final es

${\displaystyle v_f=u\ln \left(\frac{m_0}{m_0-c_0}\right)+u\ln \left(\frac{m_1}{m_1-c_1}\right)+u\ln \left(\frac{m_2}{m_2-c_2}\right)}$

Teniendo en cuenta que

${\displaystyle \begin{array}{l}{m}_0=M+(1+r)c\\ m_1=M+(1+r)(c_1+c_2)=M+(1+r)(p_1+p_2)c\\ m_2=M+(1+r)c_2=M+(1+r)p_{2}c\\ m_0-c_0=M+(1+r-p_0)c=M+(r+p_1+p_2)c\\ m_1-c_1=M+(1+r)p_{2}c+r{p}_{1}c\\ m_2-c_2=M+r{p}_{2}c\end{array}}$

La velocidad final, v_f, cuando se ha agotado el combustible, es

${\displaystyle \begin{array}{l}{v}_f=u\ln \left(\frac{m_0}{m_0-c_0}\frac{m_1}{m_1-c_1}\frac{m_2}{m_2-c_2}\right)=u\ln \left(\frac{M+(1+r)c}{M+(r+p_1+p_2)c}\frac{M+(1+r)(p_1+p_2)c}{M+(1+r)p_{2}c+r{p}_{1}c}\frac{M+(1+r)p_{2}c}{M+r{p}_{2}c}\right)=\\ u\ln \left(\frac{m+1+r}{m+r+p_1+p_2}\frac{m+(1+r)(p_1+p_2)}{m+(1+r)p_2+r{p}_1}\frac{m+(1+r)p_2}{m+r{p}_2}\right),m=\frac{M}{c}\end{array}}$

Ya que m₀ no depende de p_i calculamos el extremo de la función

${\displaystyle f(p_1,p_2)=\frac{\left(m+(1+r)(p_1+p_2)\right)\left(m+(1+r)p_2\right)}{\left(m+r+p_1+p_2\right)\left(m+(1+r)p_2+r{p}_1\right)\left(m+r{p}_2\right)}}$

Se trata de una función de dos variables p₁ que denominamos x y p₂ que denominamos y

${\displaystyle f(x,y)=\frac{\left(m+(1+r)(x+y)\right)\left(m+(1+r)y\right)}{\left(m+r+x+y\right)\left(m+(1+r)y+rx\right)\left(m+ry\right)}}$

El extremo dela función f(x,y) se calcula

${\displaystyle \frac{\partial f(x,y)}{\partial x}=\mathrm{0,}\frac{\partial f(x,y)}{\partial y}=0}$

Math Symbolic de MATLAB nos ayuda a calcular estas derivadas parciales y a resolver mediante el comando solve el sistema no lineal de dos ecuaciones con dos incógnitas x e y

syms m r x y;
z=(m+(1+r)*(x+y))*(m+(1+r)*y)/((m+r+x+y)*(m+(1+r)*y+r*x)*(m+r*y));
z_x=diff(z,x);
z_x=simplify(z_x);
z_y=diff(z,y);
z_y=simplify(z_y);

c=9000; %combustible
M=800;  %carga útil
z_x=subs(z_x,{m,r},{M/c,0.05});
z_y=subs(z_y,{m,r},{M/c,0.05});
eqns=[z_x==0, z_y==0];
s=solve(eqns,[x,y]);
disp(double(s.x))
disp(double(s.y))

   0.2654 + 0.0000i
  -0.1327 + 0.5730i
  -0.1327 - 0.5730i

   0.1134 + 0.0000i
  -0.1837 - 0.1716i
  -0.1837 + 0.1716i

• La primera raíz real es el valor de x_m=p₁=0.2654
• La segunda raíz real es el valor de y_m=p₂=0.1134
• Naturalmente, p₀=1-p₁-p₂=0.6212

El reparto óptimo de combustible para este cohete de tres etapas es, 62.1% en la primera etapa, 26.5% en la segunda y 11.3% en la tercera

Vamos a comprobarlo representando la superficie f(x,y) y el valor de su extremo, z_m=f(x_m,y_m)

c=9000; %combustible
M=800;  %carga útil
r=0.05;
m=M/c;
%velocidad final del cohete
f=@(x,y) (m+(1+r)*(x+y)).*(m+(1+r)*y)./((m+(r+x+y)).*(m+(1+r)*y+r*x).*(m+r*y));
x=0:0.01:1;
y=0:0.01:1;
[X,Y]=meshgrid(x,y);
hold on
mesh(X,Y,f(X,Y));
x_m=0.2654;
y_m=0.1134;
z_m=f(x_m,y_m);
plot3(x_m, y_m, z_m,'ko','markersize',4,'markeredgecolor','k',
'markerfacecolor','k')
line([x_m,x_m], [y_m,y_m], [7,z_m],'lineStyle','--' )
line([0,x_m], [y_m,y_m], [7,7] )
line([x_m,x_m], [0,y_m], [7,7])
hold off
grid on
xlabel('x')
ylabel('y')
zlabel('z')
title('Cohete de tres etapas')
view(27,34)

Calculamos la velocidad (máxima posible) del cohete de tres etapas, para las proporciones óptimas

u=2000; %velocidad de escape de los gases (respecto al cohete)
c=9000; %combustible
M=800;  %carga útil
r=0.05;
m=M/c;
%velocidad final del cohete
p1=0.2654; %proporciones óptimas
p2=0.1134;
vm=u*log((m+1+r)*(m+(1+r)*(p1+p2))*(m+(1+r)*p2)/((m+r+p1+p2)*(m+(1+r)*p2+r*p1)
*(m+r*p2)));
disp(vm)

   4.7297e+03

Obtenemos con tres etapas una velocidad de 4 730 m/s un poco mayor que con dos etapas 4 637 m/s