En esta sección, intentamos trasmitir algunas ideas volcadas en el foro de A.C.E.M.U., sobre simulaciones de vuelo resueltas a través del Método Runge-Kutta de resolución de ecuaciones diferenciales.

Trabajando con las curvas de empuje del motor (ACEMU MX#001), se obtuvieron unos datos interesantes en cuanto a predicción de las posibles velocidades y alturas del vector.

Para comenzar, aclaremos cuales fueron las aproximaciones utilizadas:

1. 1° se tomó la curva de empuje obtenida, pasando los datos a Newton.
2. 2° se asumieron dos etapas de vuelo:
   • Etapa 1: con el motor haciendo empuje sobre el vector.
   • Etapa 2: vuelo libre (sólo actúa la fuerza de la gravedad)

También se asumió que no hay rozamiento con el aire, o sea que la ecuación utilizada fué:

y''(t)=(E(t)-P)/m

donde,

y''(t) es la aceleración (derivada segunda de la altura y(t))
E(t) es la fuerza de empuje del motor al tiempo t
P es el peso del vector
m es la masa del vector

Se asumió que la masa del cohete no varía con el tiempo, la cual se fijó en 1.8 kg, por lo tanto  el peso es P=m*9.8

Como E(t) se encuentra bajo la forma de una tabla de valores, no cabría otra opción que resolver numéricamente el problema, para lo cual se utilizó Runge-Kutta de orden 2.
Todo fue programado empleando el lenguaje AWK Referencia al lenguaje AWK en wikipedia.

También se asumió que lleva paracaídas, pero no se consideró su apertura.

Tiempo de vuelo total: aprox 20 s
Velocidad máxima alcanzada: 90 m/s a los 1.3 segundos de vuelo (en ese momento se encuentra a unos 65 mts. de altura)
Apogeo: 485 m a los 10 s
Velocidad final del vector: -97 m/s (recordemos que no se considera apertura de paracaídas)

1) Usar Runge-Kutta de orden 4
2) Modificar la ecuación para introducir fricción con el aire
3) Incluir la apertura del paracaídas en la simulación
4) Arrancar la simulación cuando el Empuje supera el Peso del cohete, a los efectos de no recortar a mano la tabla de datos de empuje, y evitar consideraciones tales como que al principio el cohete “cae” (obviamente la simulación no tiene en cuenta que el cohete está apoyado)

Código: [Seleccionar]
# Programa AWK para calcular altura y velocidad de un cohete a partir de
# los datos de empuje del motor.
# El método empleado es Runge-Kutta de orden 2
# Inicializamos algunos datos
BEGIN { 
  n=0; m=1.8; P=m*9.8 ; h=1/60; v=0; y=0; t=0 
}
# Primero leemos todos los datos. El archivo de datos tiene 2 columnas: tiempo
# y empuje. Sólo nos interesa el empuje.
{ E[n]=9.8*(($2+0.449)*19.198-0.1311);
  n=n+1;
}
# Ahora realizamos la simulación del vuelo
END {
# Primera etapa: motor encendido, ejerciendo empuje sobre el cohete
  for (i = 0; i < n; i++){
    K1=h*v
    L1=h*(E[i]-P)/m
    K2=h*(v+L1)
    L2=h*(E[i]-P)/m
    y=y+(K1+K2)/2
    v=v+(L1+L2)/2
    t=(i+1)*h
    printf "%6.3f %8.2f %8.2f\n",t,y,v;
  }
# Segunda etapa: motor apagado, vuelo libre
# La única diferencia con el caso anterior es que no hay empuje de motor (E), 
# la única fuerza que actúa sobre el cohete es la fuerza de la gravedad (P).
# Seguimos iterando hasta que la altura pasa a ser negativa (tocó tierra)
#print "cambiamos a vuelo libre";
  while (y>0) {
    K1=h*v
    L1=h*(-P)/m
    K2=h*(v+L1)
    L2=h*(-P)/m
    y=y+(K1+K2)/2
    v=v+(L1+L2)/2
    i=i+1
    t=i*h
    printf "%6.3f %8.2f %8.2f\n",t,y,v;
  }
}

En cuanto a la ecuación la misma se basa en la Ley de Newton: F=ma

En este caso F(t)=E(t)-P(t), es la fuerza resultante que actúa sobre el cohete, que sería la resta del empuje del motor, menos el peso del cohete (suponemos además que el cohete despega y vuela de manera completamente vertical).

Todas son funciones del tiempo, porque el empuje del motor varía a medida que transcurre el tiempo y lo mismo ocurre con la masa del cohete, de ahí que el peso P también varíe con el tiempo.
Una aproximación fue suponer, que P(t)=constante, por eso

F(t)=E(t)-P

Del otro lado de la ecuación tenemos la masa del cohete multiplicada por la aceleración que sufre el mismo.

La aceleración, no es otra cosa que la derivada segunda de la distancia, (en nuestro caso la altura 'y' del cohete), y que suele representarse con dos apóstrofes

(y'')

para indicar que es la derivada segunda.

Por eso es que podemos despejar la aceleración que sufre el cohete en la ecuación de la siguiente manera:

y(t)''=F/m

o sea

y(t)''=(E(t)-P)/m

donde los valores de E(t) son los que obtuvimos en la toma de datos del ensayo del motor.
Teniendo en cuenta que el peso P=mg, en realidad podríamos escribir:

y''=E/m - g

Esa es la ecuación diferencial que estamos resolviendo numéricamente.

Lo único que determina la solución serán las condiciones iniciales (altura y velocidad inicial) que asumimos valían cero.

Una mejora que puede agregarse es tener en consideración la fricción del aire. Generalmente la fricción es una función de la velocidad: a mayor velocidad, mayor fricción.

Aparecería entonces una nueva fuerza:

F(t)=E(t) - P - f(y')

donde f(y') es la fuerza de fricción, que depende de la velocidad y' (derivada primera de la distancia).

Así la ecuación diferencial quedaría:

y''=E(t) - P - f(y')

y se resuelve de manera similar a la anterior.\\ Si la altura que alcanza el cohete es MUY alta, entonces la fricción también dependerá de la altura (ya que varía la densidad del aire y ese es uno de los factores que influye en la fricción), pero para bajas alturas se puede suponer constante y por lo tanto independiente de la altura que alcance el cohete.

Hasta este punto la simulación resulta bastante sencilla, porque los datos o los tenemos, o los va dando la simulación (como la altura y la velocidad en cada instante).
Se complica un poco la estimación del ritmo de variación de masa del motor, para tenerlo en cuenta de manera de mejorar la simulación.
Seguramente se podrá estimar la pérdida de masa en cada instante, mediante el cociente de la fuerza de empuje en el instante t, dividido la integral de la fuerza de empuje, multiplicado por la masa de combustible.
Eso daría la fracción de masa perdida en ese delta t…, suponiendo que la fuerza de empuje es proporcional a la masa eyectada, lo cual no es estrictamente cierto, porque también es proporcional a la velocidad de eyección de esa masa, pero a los efectos de este análisis lo tomaremos como un dato válido.

Igual no creemos que la influencia de la variación de masa sea demasiado importante. Creemos que el término de de fricción debe ser más importante y sería muy sencillo incluir.

Hasta ahí la simulación es 1D… suponemos que el cohete sube derecho y baja derecho (y de cola),habría que, al menos llegar a una simulación 2D, que debe incluir el efecto del viento, o por lo menos un ángulo de despegue.

Con el método de Runge-Kutta, en lugar de resolver la ecuación uno lo que obtiene, son aproximaciones numéricas a la solución de la ecuación para ciertos valores de la variable independiente, en este caso el tiempo.

Por supuesto que las condiciones iniciales (altura y velocidad inicial) determinan la solución.

Lo que no aclaramos, es que al resolver la ecuación diferencial de segundo orden, la solución es la altura y(t).
Como el método de Runge-Kutta se puede aplicar solamente a EDOs (Ecuaciones Diferenciales Ordinarias) de primer orden, lo que tenemos que hacer para aplicarla a una de segundo orden, es aplicar un cambio de variable y transformar la EDO de 2do orden en un sistema de 2 EDOs de 1er orden.

En nuestro caso el cambio de variable lo podríamos escribir así:

y'=v

(la derivada primera de la altura es la velocidad)

Por lo tanto y''=v' (la aceleración es la derivada primera de la velocidad)

y eso da lugar al sistema de 2 EDOs de 1er orden:

y'=v
v'=E(t)-P

La solución de la primera ecuación es la altura y(t) y le corresponde en nuestro caso la condición inicial y(0)=0

La solución de la segunda ecuación es la velocidad v(t) y le corresponde la condición inicial v(0)=0

Podríamos haber puesto otras condiciones iniciales, pero lo normal en nuestro caso sería que el cohete arranque desde el piso con velocidad cero.

Ahora que tenemos eso, le aplicamos la “receta Runge-Kutta” que es muy simple.

Supongamos que tenemos una Ecuación Diferencial Ordinaria genérica,  expresada de la siguiente manera:

y'=f(t,y)

Esto es que el término de la derecha, f(t,y) es una función que depende exclusivamente de la variable independiente (el tiempo), y la solución de la ecuación diferencial, y(t) que es lo que queremos encontrar.

Entonces, si conocemos un valor de la solución a tiempo t, aplicando la receta de Runge-Kutta, podemos encontrar una aproximación a la solución a tiempo t+h, donde h es el incremente, en nuestro caso 1/60 s.

Para eso, hacemos el siguiente cálculo :

K1=h*f(t0,y0)
K2=h*f(t0+h,y0+K1)
y1=y0+(K1+K2)/2
t1=t0+h

t0 e y0 son las condiciones iniciales, con eso encontramos un valor aproximado a la solución para t1=t0+h, y a ese valor le llamas y1.

OK, ¿qué hacemos ahora?

Pues muy simple, ahora conocemos la solución a tiempo t1, podemos volver a aplicar la receta Runge-Kutta para encontrar la solución a tiempo t2=t1+h.

A esto lo denominamos y2.

Ahora con t2 e y2, volvemos a aplicar la receta y calculo y3 para el tiempo t3.

… y así seguimos, todo lo necesario.

Ahora bien, la receta que vimos, es para aplicar a un EDO de 1er orden, pero nosotros tenemos un sistema de 2 EDOs de 1er orden.

No pasa nada, supongamos que el sistema de EDOs es:

y'=f(t,y,v)
v'=g(t,y,v)

con condiciones iniciales t0, y0 y v0 (veamos que a cada ecuación debe corresponderle una condición inicial de la correspondiente solución).

Le aplicamos la receta a CADA ecuación en ORDEN.
Para no confundir, identificaremos con K las Kes de la 1er ecuación y con L las Kes de la segunda ecuación. Así que al aplicar la receta queda así:

K1=h*f(t0,y0,v0)
L1=h*g(t0,y0,v0)
K2=h*f(t0+h,y0+K1)
L2=h*g(t0+h,y0+K1,v0+L1)
y1=y0+(K1+K2)/2
v1=v0+(L1+L2)/2
t1=t0+h

Tenemos un sistema de N EDOs; aplicamos la receta N veces, una a cada EDO, realmente sencillo.
Claro que si tenemos N EDOs lo más probable es que convenga trabajar matricialmente.

A esta altura, si miramos el programa de AWK que vimos más arriba, reconoceremos la “receta de Runge-Kutta” por ahí.
No pretendamos entender al detalle el funcionamiento del programa, porque para quien sabe algo de programación le va a parecer que le faltan cosas.
Por ejemplo la parte que se lee la tabla de datos debe parecer muy oscura (es una sóla linea), pero eso es porque AWK ya tiene implícito la lectura de datos: quienes lo diseñaron sabían bien lo que hacían. AWK es un programa para procesar datos de texto, así que para qué molestarse en que el usuario tenga que programar eso: se da por defecto.

Luego que termina la parte de lectura de datos, el programa está dividido en 2 secciones: mientras tenemos datos de la tabla de empuje, se calculan las ecuaciones considerando el empuje, cuando se termina la tabla de valores de empuje (el motor se apagó), arranca la segunda sección de vuelo libre.
Vemos que en el cálculo desapareció el empuje E, y sólo queda el peso P (en realidad como va dividido por m, alcanzaría con poner g es lugar de P/m).
Por esto consideramos al programa tan sencillo, al igual que la receta.

Volviendo al método numérico, Runge-Kutta de orden 2 no es tan bueno como las recetas de Runge-Kutta de orden superior, pero para comenzar está bien.

Veamos ahora la receta de orden 4:

K1=h*f(t,y)
K2=h*f(t+h/2,y+K1/2)
K3=h*f(t+/2,y+K2/2)
K4=h*f(t+h,y+K3)
t=t+h
y=y+(K1+2*K2+2*K3+K4)/6

La fórmula de fuerza de fricción que estamos usando es:

f(v)=0.5*Cd*da*A*v²

donde:

Cd: es el coeficiente de fricción del cuerpo, que en el caso del cohete consideramos valía 0.7, según una publicación de NAR.
Este coeficiente es muy variable, depende de varios factores, no sólo de la forma del cuerpo, sino también de la velocidad (a través del número de Mach y el número de Reynolds, que mide la turbulencia del fluido que atraviesa), pero para bajas velocidades (subsónicas) se puede considerar independiente de Re y Ma.
Acá es donde toma importancia la terminación del cohete: que tan áspera o suave sea la superficie, que elementos extraños perturben el flujo del aire, etc.

da: es la densidad del aire, tomada como constante y de magnitud 1.29 kg/m³

A: es el área de la sección transversal del cohete, que es Pi*r², medido en m² , por lo tanto r=0.025 m

v: es la velocidad del vector.

Es evidente que esta fuerza varía también en el tiempo, a medida que varía la velocidad del vector.

Así que ahora las ecuaciones quedan así:

y'=v
v'=(E(t)-P-f(v))/m

A continuación el link del fuente AWK para la simulación en cuestión, utilizando Runge-Kutta de orden 4 Simulación - Programa escrito en AWK para Runge-Kutta 4.

Retorno a página anterior