En este artículo, intentamos trasmitir en forma de artículo, un posteado en el Foro de A.C.E.M.U., realizado por el maestro Kenneth Irving, sobre las simulaciones antes vistas, pero resueltas a través del Método Runge-Kutta de resolución de ecuaciones diferenciales.
Para no dejar ideas fuera de contexto, dejamos el posteado tal cual fué realizado, por eso el diálogo en primera persona.

Estimados: aprovechando el feriado, me puse a jugar con la curva de empuje del motor (la “normal”), y llegué a unos número interesantes en cuanto a predicción de velocidad y altura del vector.

Para empezar, aclaremos cuales fueron las aproximaciones que usé:

1. usé la curva de empuje obtenida, pasando los datos a Newton.
2. asumí 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 asumí que no hay rozamiento con el aire, o sea que la ecuación que usé es:

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

Asumí que la masa del cohete no varía con el tiempo y la fijé en 1.8 kg, por lo tanto  el peso es P=m*9.8

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

Asumí que el cohete lleva paracaídas… pero que no abre.

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 m de altura)
Apogeo: 485 m a los 10 s
Velocidad final del vector: -97 m/s (dejó un lindo agujerito en la tierra)

1) Usar RK 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, así no hay que andar recortando a mano la tabla de datos de empuje, y evitamos tonterías 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 no hay misterio alguno, es física de liceo, la vieja 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, eso me olvidé de aclararlo).

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 que hice fue suponer que P(t)=constante, por eso

F(t)=E(t)-P

Del otro lado de la ecuación tenés la masa del cohete multiplicada por la aceleración que sufre el cohete.

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 podés 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 estoy resolviendo numéricamente.
OK, cómo hago eso no lo voy a explicar acá, porque el fundamento es más complicado, pero es una simple receta, una vez que entendiste de donde sale la receta, la aplicás a ojos cerrados y no tiene ningún misterio.

Lo único que determina la solución serán las condiciones iniciales (altura y velocidad inicial) que asumí que ambas 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 aqui es realmente sencillo de simular, porque los datos o los tenés, o te los va dando la simulación (como la altura y la velocidad en cada instante).
Se me 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, pero ya le voy a encontrar una vuelta.
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 me daría la fracción de masa perdido 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 más o menos sirve).

Igual no creo que la influencia de la variación de masa sea demasiado importante. Calculo que el término de fricción debe ser más importante y es realmente una pavada de incluir.

Hasta ahí la simulación es 1D… suponemos que el cohete sube derechito y baja derechito (y de culete)… 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 aclaré 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 tenés 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

(lo cual es obvio, pues uno lo que lee es: 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 tenés eso, le aplicás la “receta Runge-Kutta” que es MUY simple.

Supongamos que tenés una Ecuacion Ddiferencial 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 conocés un valor de la solución a tiempo t, aplicando la receta de RK podés 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, hacés el siguiente cálculo (esto tomalo como viene, otro día te explico de donde sale):

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

¡Listo! t0 e y0 son las condiciones iniciales, con eso encontraste un valor aproximado a la solución para t1=t0+h, y a ese valor le llamás y1.

OK, ¿qué hago ahora?

Pues muy simple, ahora conozco la solución a tiempo t1, puedo volver a aplicar la receta RK para encontrar la solución a tiempo t2=t1+h.

Llamale a eso y2.

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

… y así sigo, todo lo que quiera. Es realmente una papa  8)

Ahora bien, la receta que te pasé es para aplicar a 1 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 condciones iniciales t0, y0 y v0 (fijate que a cada ecuación debe corresponderle una condición inicial de la correspondiente solución).

Le aplicás 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

¡y listo! Si tenés un sistema de N EDOs… aplicás la receta N veces, una a cada EDO… una verdarera chotada. Claro que si tenés N EDOs lo más probable es que te convenga trabajar matricialmente, porque te podés llegar a hacer un entrevero lindo  

A esta altura, si mirás el programa de AWK que escribí, vas a empezar a reconocer la “receta de RK” por ahi. No pretendas entender al detalle el funcionamiento del programa, porque si sabés algo de programación te va a aprecer 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: los tipos que 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 tengo datos de la tabla de empuje, calculo las ecuaciones considerando el empuje, cuando se me terminó la tabla de valores de empuje (el motor se apagó) arranca la segunda sección, de vuelo libre. Viste 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).  Pero viste que el programa es rechoto, igual que la receta.

OK, volviendo al método numérico, RK de orden 2 no es tan bueno como las recetas de RK de orden superior, pero para empezar está bien.

Te paso 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

al final soy un nabo, demoro más en escribir toda esta explicación que en reescribir el programa para hacerlo con RK4 e incluir la fricción del aire

Bueno, acabo de hacer la simulación incluyendo la fricción del aire.

La fórmula de fuerza de fricción que estoy 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 consideré que 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, que tomé 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 tomé 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

Retorno a página anterior