El siguiente es el fuente para la simulación de altura y velocidad de un cohete, con los datos de empuje del motor, utilizando el método de Runge-Kutta de 4to orden.

La base de este programa, es la vista en la sección anterior Simulación por Método Runge Kutta

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 4
# Inicializamos algunos datos
BEGIN {
  n=0; m=1.8; P=m*9.8 ; h=1/60; v=0; y=0; t=0;
  Cd=0.7;              #Coeficiente de Fricción del Cohete (drag Coeff)
  da=1.29;             #Densidad del aire
  A=0.00196            #Area de sección transversal del cohete: radio=2.5cm
  C=0.5*A*Cd*da        #coeficiente para calculo de drag del cohete
}
# 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
# Peso (P) y fricción del aire ofreciendo resistencia al avance
  for (i = 0; i < n; i++){
    sig=1
    if (v < 0) sig=-1
    K1=h*v
    L1=h*(E[i]-P-C*v*v*sig)/m
    K2=h*(v+L1/2)
    v2=(v+L1/2)
    v2=v2*v2
    L2=h*(E[i]-P-C*v2*sig)/m
    K3=h*(v+L2/2)
    v2=(v+L2/2)
    v2=v2*v2
    L3=h*(E[i]-P-C*v2*sig)/m
    K4=h*(v+L3)
    v2=(v+L3)
    v2=v2*v2
    L4=h*(E[i]-P-C*v2*sig)/m
    y=y+(K1+2*K2+2*K3+K4)/6
    v=v+(L1+2*L2+2*L3+L4)/6
    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),
# Las únicas fuerzas que actúan sobre el cohete son el Peso (P) y la fricción del aire
# Seguimos iterando hasta que la altura pasa a ser negativa (tocó tierra)
#print "cambiamos a vuelo libre";
  while (y>0) {
    sig=1
    if (v<0) sig=-1
    K1=h*v
    L1=h*(-P-C*v*v*sig)/m
    K2=h*(v+L1/2)
    v2=(v+L1/2)
    v2=v2*v2
    L2=h*(-P-C*v2*sig)/m
    K3=h*(v+L2/2)
    v2=(v+L2/2)
    v2=v2*v2
    L3=h*(-P-C*v2*sig)/m
    K4=h*(v+L3)
    v2=(v+L3)
    v2=v2*v2
    L4=h*(-P-C*v2*sig)/m
    y=y+(K1+2*K2+2*K3+K4)/6
    v=v+(L1+2*L2+2*L3+L4)/6
    i=i+1
    t=i*h
    printf "%6.3f %8.2f %8.2f\n",t,y,v;
  }
}

En este caso al incluir el término de fricción del aire, se obtienen los siguientes resultados :

Tiempo de vuelo total: aprox. 18.7 seg.
Velocidad máxima alcanzada: 87.97 m/s a los 1.283 segundos de vuelo (en ese momento se encuentra a unos 64.6 mts. de altura)
Apogeo: 401 mts. a los 9.35 seg.
Velocidad final del vector al llegar a tierra: -80.6 m/s (no se considera apertura de paracaídas)

La diferencia con el vuelo sin fricción es de nada menos que 84 mts.

Una cosa en la que hubo que improvisar fue en la manera que actúa la fuerza de fricción, porque siempre se opone al movimiento, eso quiere decir que cambia de signo cuando el cohete empieza a caer.
Lamentablemente AWK es un lenguaje para procesar textos, no para cálculo científico, y si bien cuenta con una muy pequeña biblioteca de funciones matemáticas, no tiene la función valor absoluto.

Para poder solucionar este problema, incluimos una variable “sig” que registra el signo de la velocidad, de manera de cambiar el sentido en que se aplica la fuerza de fricción del aire en la ecuación diferencial.

Evidentemente con esto estamos llegando al límite de lo que puede darnos AWK para resolver este tipo de problemas.

Realizamos algunos ajustes al modelo de cohete con OpenRocket, para lo cual elegimos un motor con un impulso total similar al valor del MX#001, de 185 Ns.
Se ajustó el peso del cohete con una más cercana a la que creemos sea la real, y se efectuó una simulación con un motor similar (o lo más parecido posible) al MX#001.

La altura obtenida con estos datos utilizando OpenRocket, fue de 760 mts.

Utilizando el programa de simulación antes visto, usando los datos de empuje reales del motor, obtuvimos una altura de 790 metros, que coincide bastante bien con la simulación del OpenRocket.

Somos conscientes que no estamos comparando las mismas cosas, porque hay una discrepancia en la masa y en los motores, pero para ser datos similares, reconforta que la simulación propuesta arrojo resultados parecidos.

Respecto a las simulaciones vistas en Simulación por Método Runge Kutta, se corrigió la masa del cohete, reduciendo la misma en 300 grs, y ajustando el coeficiente de fricción a lo que maneja OpenRocket: 0.61.
Ajustando ambos valores las simulaciones se parecen bastante. 

De todos modos creemos que pueden haber ocurrido algunas compensaciones de errores, que nos acercaron a los resultados de OpenRocket.

Las dos diferencias fundamentales entre los cálculos son:

1. no usamos el mismo motor, sino motores con igual impulso total.
2. en el cálculo en sí, usamos un coeficiente de fricción constante y no variable como creemos usa Openrocket, pero como las velocidades no son tan altas, y el período de aceleración es muy corto, las discrepancias no serían tan grandes.

Aunque el programa es absolutamente casero y está a años luz de estos otros programas, al menos parece reproducir razonablemente bien la altura máxima que puede alcanzar el vector.

Claro que eso no es lo mismo que simular el vuelo; lo que estamos haciendo, es asumir que el vector se mueve en una sola dimensión (la vertical), y sabemos que el cohete en realidad puede volar para cualquier lado, pero al menos este programa sirve para estimar la altura máxima que podría alcanzar el cohete.

Retorno a página anterior