Tesis Doctoral Modelización Proyectiles Base Burn Fernando Aguirre Estevez

Transcript

ESCUELA POLITÉCNICA SUPERIOR DEL EJÉRCITO TESIS DOCTORAL MODELIZACIÓN DE PROYECTILES BASE BURN Director: Dr. Coronel Francisco Cucharero Pérez Doctorando: Comandante Fernando Aguirre Estévez I.-1 I.-2 ESCUELA POLITÉCNICA SUPERIOR DEL EJÉRCITO TESIS DOCTORAL CAPÍTULO I ÍNDICE MODELIZACIÓN DE PROYECTILES BASE BURN Director: Dr. Coronel Francisco Cucharero Pérez Doctorando: Comandante Fernando Aguirre Estévez I.-3 I.-4 I.-ÍNDICE I. Índice II. Introducción III. Estado del arte IV. Aerodinámica de un proyectil base burn V. Análisis de un proyectil base burn VI. Cálculo numérico de un proyectil base burn VII. Evaluación del centro de masas y del momento de inercia de un proyectil asistido VIII. Optimización IX. Algoritmo de recocido simulado X. Algoritmos genéticos XI. Recocido simulado frente a algoritmo genético XII. Modelización de proyectiles base burn método 1 XIII. Modelización de proyectiles base burn método 2 XIV. Conclusiones XV. Bibliografía XVI. Programas I.-5 I.-6 ESCUELA POLITÉCNICA SUPERIOR DEL EJÉRCITO TESIS DOCTORAL CAPÍTULO II INTRODUCCIÓN MODELIZACIÓN DE PROYECTILES BASE BURN Director: Dr. Coronel Francisco Cucharero Pérez Doctorando: Comandante Fernando Aguirre Estévez II.-1 II.-2 II.-INTRODUCCIÓN 1. Resumen El aumento del alcance sin el aumento del calibre ha supuesto la incorporación de determinados dispositivos a los proyectiles de artillería. En este punto, se han formulado dos tendencias fundamentales; en primer lugar, el incorporar sistemas de propulsión adicional a los proyectiles lo cual ha provocado el nacimiento de los cohetes y de los proyectiles de propulsión adicional (RAP); y en segundo lugar, se han incorporado dispositivos que permiten disminuir la resistencia de los proyectiles inyectando algún fluido por el culote del proyectil lo que se ha denominado como proyectiles base burn o proyectiles base bleed (BB). Los fundamentos teóricos que sostienen la modelización numérica en la que se basa un proyectil base burn se encuentra en forma de STANAG´s, artículos y los códigos de software de NABK (NATO Armaments Ballistic Kernel). Estos documentos contienen la información necesaria para modelizar el vuelo de un proyectil base burn pero para ello es necesario el conocimiento de una serie de coeficientes balísticos. En este trabajo se intentará clarificar la tecnología asociada a un proyectil base burn, su modelización y tratamiento matemático. La primera dificultad que aparece es que la modelización de proyectiles base burn se ha incluido en el STANAG 4355 de dos formas diferentes. Por un lado, se tiene el método 1 (método empleado por los países francófonos) y por el otro, el método 2 (método desarrollado por los EE.UU. y empleado por los países anglófonos). Estos dos métodos difieren totalmente entre si, siendo las bases de datos muy diferentes. La utilización de uno u otro modelo va a depender del país propietario de la munición y cada uno de ellos tiene sus inconvenientes y sus ventajas. El método 1 resulta más intuitivo al utilizar conceptos como velocidad de combustión, densidad del propulsante y superficie de quemado sumando al sistema de siete ecuaciones del Modelo Modificado de Masa Puntual una ecuación de variación de la masa mientras que el método 2 incorpora tres ecuaciones más: la ecuación de variación de la masa, la del pseudotiempo y la del tiempo de apagado del motor. Aún existe otro modelo más del que se tiene noticia, modelo II.-3 desarrollado por Suecia, y que no aparece en el STANAG 4355. Se comenzará este estudio considerando la resistencia aerodinámica de un proyectil así como una descripción física de la unidad base burn. A continuación se analizará detalladamente el Modelo Modificado de Masa Puntual para proyectiles base burn, donde se construirá una aplicación informática en lenguaje C++ que mostrará las trayectorias balísticas de diferentes proyectiles. El cálculo de trayectorias balísticas efectuado numéricamente de un proyectil base burn necesita del conocimiento de los coeficientes aerodinámicos, los coeficientes base burn y los factores de ajuste. Los coeficientes aerodinámicos pueden obtenerse a través de programas informáticos, de ensayos en túneles de viento o bien de bases de datos suministradas por el fabricante del proyectil. No obstante, al aplicar el Modelo Modificado de Masa Puntual para proyectiles base burn existe un conjunto de coeficientes adicionales relativos al proceso de quemado del propulsante así como a la disminución de la resistencia producida por el aumento de la presión en la base. Estos coeficientes pueden obtenerse de ensayos en laboratorio o bien de bases de datos suministradas por el fabricante del proyectil. La obtención del paquete aerodinámico no supone un problema en la actualidad a través de aplicaciones informáticas, si bien la obtención de los coeficientes adicionales base burn constituye un problema superior. La principal aportación de esta tesis ha sido la determinación de estos coeficientes a través de resultados de alcance, deriva y tiempo de extinción del motor base burn obtenidos de experiencias realizadas mediante la resolución de un problema de optimación complejo, desarrollando para ello una serie de herramientas numéricas. Finalmente, y a efectos de ajustes finales en el cálculo de trayectorias se deben encontrar cuatro parámetros denominados factores de ajuste, como se indica en los STANAG 4355 y STANAG 4144, para la correcta modelización numérica según el Modelo Modificado de Masa Puntual para proyectiles base burn. II.-4 2. Palabras clave Ácido desoxirribonucleico (ADN) Deoxyribonucleic acid (DNA) Acuerdo de estandarización Standardization Agreement (STANAG) Alcance Range Aleatorio Random Algoritmo de descenso Descent algorithm Algoritmo genético Genetic algorithm (GA) Altura de quemado Height of burst (HOB) Ángulo de ataque Angle of attack Ángulo de tiro Quadrant elevation (QE) Búsqueda en zonas próximas Neighbourhood search Búsqueda local Local search Búsqueda taboo Taboo search Capa de cortadura Shear layer Capa límite Boundary layer Centro de masas Center of mass Coeficiente balístico Ballistic coefficient Combustión Combustion Condición de parada Stopping criterion Condición de parada Termination criterion Corriente libre Freestream Cruce Crossover Deriva Drift Encendido Igniter Entradas de control de fuego Fire control inputs (FCI) Estela Wake Estocástico Stochastic Finalización de quemado Burnout Factor de ajuste Fitting factor Flujo másico Mass flow Fuerza de resistencia Drag force II.-5 Fuerza lateral Lift force Fuerza Magnus Magnus force Función objetivo Fitness Función objetivo Objective function Generador de gas Gas generator Heurística Heuristic Integración Integration Inyección de propulsante Fuel injection Masa de propulsante Fuel mass Mecánica de fluidos computacional Computational Fluid Dynamics (CFD) Mecánica estadística Statistical mechanics Mínimo global Global minima Mínimo local Local minima Momento de giro Overturning moment Momento de inercia Axial moment of inertia Mutación Mutation Núcleo Balístico de Armamento de la OTAN NATO Armaments Ballistic Kernel (NABK) Número de Mach Mach number Número de Reynolds Reynolds number Obús Howitzer Onda de choque Shock wave Onda de expansión Expansion wave Optimización Optimization Optimización combinatoria Combinatorial optimization Parámetro de inyección Injection parameter Presión de base Base pressure Presión en la cámara Chamber pressure Proyectil con resistencia de culote atenuada Base bleed projectile (BB) Proyectil con resistencia de culote atenuada Base burn projectile (BB) Proyectil de propulsión adicional Rocket assisted projectile (RAP) Punto de remanso Stagnation point Recocido simulado Simulated annealing (SA) II.-6 Reproducción Reproduction Retraso a la ignición Ignition delay Runge Kutta Fehlberg Runge Kutta Fehlberg Solución óptima Optimal solution Soplado Bleeding Temperatura adaptativa Adaptive temperature Temperatura del propulsante Propellant temperature Trayectoria Trajectory Tubo Gun Turbulencia Turbulence Vector de ataque Yaw of repose Velocidad de enfriamiento Cooling schedule Velocidad de rotación axial Axial spin rate Velocidad del sonido Speed of sound Velocidad en boca Muzzle velocity II.-7 II.-8 AGRADECIMIENTOS Ante todo, debo expresar mi más profundo agradecimiento al Coronel del Cuerpo de Ingenieros Politécnicos D. Francisco Cucharero Pérez, gracias al cual se ha podido llevar a efecto esta Tesis. II.-9 II.-10 ESCUELA POLITÉCNICA SUPERIOR DEL EJÉRCITO TESIS DOCTORAL CAPÍTULO III ESTADO DEL ARTE MODELIZACIÓN DE PROYECTILES BASE BURN Director: Dr. Coronel Francisco Cucharero Pérez Doctorando: Comandante Fernando Aguirre Estévez III.-1 III.-2 III.-ESTADO DEL ARTE 1. Estado actual Aunque en la actualidad se conoce la tecnología asociada a un proyectil base burn, su modelización e interpretación desde un punto de vista matemático es una cuestión aparte. La mayoría de la información existente se encuentra en forma de STANAG´s y artículos siendo la principal fuente de información los STANAG 4355 y 4144 si bien el estudio de los códigos de software de NABK (NATO Armaments Ballistic Kernel) constituye una importante referencia a la hora de completar la información de los propios STANAG´s. Estos documentos proporcionan la información necesaria para modelizar el vuelo de un proyectil base burn, siempre que se hayan podido obtener los coeficientes balísticos necesarios para ello. El STANAG 4355, en su cuarta edición, trabaja con el Modelo Modificado de Masa Puntual y con el Modelo de Masa Puntual, los cuales han sido completados para modelar los proyectiles base burn. El modelo de seis grados de libertad resulta demasiado complicado, siendo su principal inconveniente, más que su complejidad, la ausencia de información, y la dificultad de obtenerla, que permita alimentar su extensísima a la vez que espinosa base de datos. Este hecho, unido a que las ventajas de este método en cuanto mejoras son prácticamente nulas, ha provocado que, en estos momentos, se trabaje con el Modelo Modificado de Masa Puntual en un noventa por ciento dejando el restante diez por ciento para el Modelo de Masa Puntual. A esto hay que añadir que la modelización de proyectiles base burn se ha incorporado en el STANAG 4355 de dos maneras. Por un lado, se tiene el método 1 (método empleado por los países francófonos) y por otro, el método 2 (método desarrollado por los EE.UU. y empleado por los países anglófonos). Estos dos métodos son completamente diferentes siendo la comunalidad de las bases de datos prácticamente cero. La utilización de uno u otro modelo depende del desarrollador y cada uno de ellos tiene sus inconvenientes y sus ventajas. El método 1 resulta más intuitivo al utilizar conceptos como velocidad de combustión, densidad del propulsante y superficie de quemado adicionando al sistema de siete ecuaciones del Modelo Modificado de Masa III.-3 Puntual una ecuación de variación de la masa mientras que el método 2 incorpora tres ecuaciones más: la ecuación de variación de la masa, la del pseudotiempo y la del tiempo de apagado del motor. El cálculo de trayectorias balísticas efectuado numéricamente de un proyectil base burn necesita del conocimiento de los coeficientes aerodinámicos, los coeficientes base burn y los factores de ajuste. Los coeficientes aerodinámicos pueden obtenerse a través de programas informáticos existentes, de ensayos en túneles de viento o bien de bases de datos suministradas por el fabricante del proyectil. Existen aplicaciones informáticas como PRODAS o McDrag que proporcionan los coeficientes aerodinámicos a partir de la geometría del proyectil en base a ciertos ajustes fenomenológicos obtenidos a partir de ensayos; o bien a partir de la integración numérica de las ecuaciones de Navier Stokes, proceso conocido abreviadamente como CFD (Computational Fluid Dynamics). La modelización de un modo pormenorizado de la combustión del propulsante así como de la reducción de la resistencia al aumentar la presión en la base se materializa en el STANAG 4355 a través de una serie de funciones que caracterizan el comportamiento del proyectil. Esto implica un conocimiento muy exhaustivo de la reacción química que ocurre al quemarse el propulsante en la cámara de combustión y que provoca una disminución de la resistencia al incrementarse la presión en la base por la inyección de los gases procedentes de la combustión del propulsante. Además, para describir apropiadamente el comportamiento de un proyectil base burn es necesario conocer cuatro parámetros denominados "fitting factors" o factores de ajuste, como se indica en los STANAG 4355 y STANAG 4144, en orden a crear correspondencia entre los resultados obtenidos computacionalmente y las experiencias observadas. Actualmente la tecnología puntera en la modelización de proyectiles base burn está representada para el método 1 por los estudios realizados por Chargelegue y Couloumy de la Direction de Armaments Terrestres y para el método 2 por los trabajos desarrollados por Gunners, Andersson y Hellgren, y por Danberg. III.-4 En el método 1, Chargelegue y Couloumy implementaron un programa informático conectado a un radar doppler de seguimiento de proyectiles que obtenía el coeficiente de reducción de resistencia base burn en función del número de Mach de los proyectiles. El coeficiente de reducción de resistencia base burn CxBB fue obtenido mediante la diferencia de dos coeficientes de reducción de resistencia base burn medidos con el radar en dos situaciones diferentes. Estos trabajos fueron expuestos en el 11º Simposio Internacional de Balística celebrado en Bruselas en 1989. En el método 2, los trabajos desarrollados por Gunners, Andersson y Hellgren, así como por Danberg mostraron los mecanismos de reducción de resistencia de base para proyectiles base burn. A partir de ciertas hipótesis de trabajo se llegó a la conclusión final de que el coeficiente de resistencia de un proyectil base burn puede considerarse igual al coeficiente de resistencia de un proyectil con la unidad base burn desactivada menos la diferencia en la componente de resistencia de base entre un motor inerte y un motor con la unidad base burn activa C D0 = C D0 − ΔC D0 . A partir de esta teoría, se bb bb desarrolló un procedimiento que los investigadores aplicaron a la modelización del proyectil M864 US CBL/DUP BB2 utilizando para ello una serie de complejos ensayos en túnel de viento. Estos estudios fueron mostrados en el 13º Simposio Internacional de Balística celebrado en Estocolmo en 1992; teniendo en cuenta, además la contribución en lo relativo a la integración numérica de las ecuaciones de Navier Stokes, realizada por Nietubicz y Sahu (1988), y expuesta en el 1º Simposio Internacional en Propulsión Química celebrado en Grecia en 1989. En este apartado también cabe destacar los trabajos actualmente en desarrollo realizados por General Dynamics Santa Bárbara Sistemas (GDSBS) y Explosivos Alaveses S.A. (EXPAL) en relación con la caracterización balística del proyectil ER-02/BB SP HEA BB1 disparado por un obús 155/52 APU SBT, empleando el método 1; y utilizando para ello una combinación de ensayos en laboratorio y experiencias. Determinando la velocidad de quemado del propulsante en el laboratorio, Explosivos Alaveses calculó el área de combustión del propulsante base burn SC a partir de la geometría de la unidad base burn. También calculó el exponente n que afecta a la velocidad de combustión de propulsante base burn VC a partir de ensayos en laboratorio. El coeficiente de reducción III.-5 de resistencia base burn CxBB fue obtenido por General Dynamics SBS computando el coeficiente de la fuerza de resistencia total con un radar. El coeficiente CxBB se obtuvo comparando la diferencia existente en el coeficiente de la fuerza de resistencia entre un proyectil con la unidad base burn activa y el mismo con la unidad base burn desactivada. Los factores de ajuste se obtuvieron mediante un procedimiento matemático a partir de los datos de alcance y deriva obtenidos de experiencias. Finalmente, se deben reseñar los trabajos de investigación desarrollados en esta tesis. No es objeto de este trabajo la determinación de los coeficientes aerodinámicos dado que este problema se encuentra resuelto en la actualidad, de una forma u otra, como se ha descrito antes. El Modelo Modificado de Masa Puntual para proyectiles base burn necesita un conjunto de coeficientes adicionales relativos al proceso de quemado del propulsante así como a la disminución de la resistencia producida por el aumento de la presión en la base. Estos coeficientes base burn pueden obtenerse de ensayos en laboratorio o bien de bases de datos suministradas por el fabricante del proyectil. Las aportaciones de esta tesis abarcan la implementación de un método capaz de obtener los coeficientes base burn necesarios para el cálculo a partir de alcance, deriva y tiempo de extinción del motor base burn obtenidos de experiencias realizadas. III.-6 2. Aportaciones Ya se ha comentado anteriormente que para la modelización de una forma detallada del proceso de quemado del propulsante así como de la reducción de la resistencia al aumentar la presión en la base aplicando el Modelo Modificado de Masa Puntual para proyectiles base burn es necesario conocer una serie de funciones que caracterizan el comportamiento del proyectil. Esto conlleva un estudio muy preciso de la reacción química que ocurre al quemarse el propulsante en la cámara de combustión y que ocasiona una disminución de la resistencia al incrementarse la presión en la base por la inyección de los gases procedentes de la combustión del propulsante. Generalmente, el conocimiento de esta información se efectúa mediante ensayos en laboratorios especializados. Estos ensayos conllevan el disponer de un equipamiento muy costoso y únicamente se lo pueden permitir países con gran capacidad económica. Una de las aportaciones de esta tesis es la construcción de un método más asequible económicamente. Este método permitirá el cálculo de estas funciones de caracterización a través de resultados de alcance, deriva y tiempo de extinción del motor base burn obtenidos de experiencias. Posteriormente podrá comprobarse cómo con los datos extraídos de experiencias se construirá una función objetivo cuyo óptimo será la solución buscada. De este modo, el problema de determinación de los coeficientes balísticos necesarios para modelizar numéricamente un proyectil base burn de acuerdo al STANAG 4355 del Modelo Modificado de Masa Puntual se reducirá a un problema de optimización que aunque complicado y laborioso podrá resolverse empleando una serie de herramientas numéricas desarrolladas a los efectos. El incremento del número de parámetros de los que dependen las funciones a optimizar confieren un entorno de resolución muy árido por lo que se han desarrollado una serie de procedimientos de búsqueda de soluciones en un proceso iterativo de cálculo. Como la modelización de proyectiles base burn se ha incorporado en el STANAG 4355 de dos maneras; método 1 (método empleado por los países francófonos) y método 2 (método desarrollado por los EE.UU. y empleado por los países anglófonos); siendo III.-7 ambos métodos completamente diferentes, es necesario desarrollar dos metodologías que, aunque nacen del mismo origen común, calculan coeficientes de caracterización diferentes. Otra de las aportaciones de esta tesis es un procedimiento para calcular los coeficientes de uno cualquiera de los dos métodos a partir de los coeficientes del otro método. 2.1. Método 1 En el caso del método 1; inicialmente, será preciso calcular para diferentes ángulos de tiro y cargas, el parámetro de inyección de eficiencia óptima I0, el área de combustión del propulsante base burn SC, el coeficiente de reducción de resistencia base burn CxBB y el exponente n que afecta a la velocidad de combustión de propulsante base burn VC a partir de alcance, deriva y tiempo de extinción del motor base burn obtenidos de resultados de experiencias. Una vez configurados adecuadamente los coeficientes base burn, se deben determinar a efectos de ajustes finales los factores de ajuste para la correcta modelización numérica según el Modelo Modificado de Masa Puntual para proyectiles base burn método 1; a partir de alcance, deriva y tiempo de extinción del motor base burn obtenidos de experiencias. Para ello, deberá tenerse en cuenta que para cada carga; únicamente el factor base burn f(iBB,MT) y el factor de sustentación fL son funciones del ángulo de tiro siendo el tiempo de retraso a la ignición del motor base burn tDI función de la temperatura de propulsante del motor MT y el factor de quemado con la velocidad de rotación axial del proyectil K(p) una función lineal de la velocidad de rotación axial del proyectil p. 2.2. Método 2 En el caso del método 2; inicialmente, será preciso calcular para diferentes ángulos de tiro y cargas, el tiempo de referencia de apagado del motor t B* , el tiempo de referencia * de retraso a la ignición del motor base burn t DI y la derivada de referencia de la masa de propulsante del motor m& *f . A continuación, para cada carga, se calcula el tiempo de retraso a la ignición del motor base burn tDI y la duración del motor base burn tB - tDI. III.-8 Todos ellos a partir de alcance, deriva y tiempo de extinción del motor base burn obtenidos de resultados de experiencias. Una vez configurados adecuadamente los coeficientes base burn, se deben determinar a efectos de ajustes finales los factores de ajuste para la correcta modelización numérica según el Modelo Modificado de Masa Puntual para proyectiles base burn método 2; a partir de alcance, deriva y tiempo de extinción del motor base burn obtenidos de experiencias. Para ello, deberá tenerse en cuenta que para cada carga; únicamente el factor base burn f(iBB,MT) y el factor de sustentación fL son funciones del ángulo de tiro siendo el factor del tiempo de quemado de la velocidad de rotación axial del proyectil base burn fBTp y el factor del tiempo de quemado de la presión del aire atmosférico local del proyectil base burn fBTP constantes. III.-9 III.-10 ESCUELA POLITÉCNICA SUPERIOR DEL EJÉRCITO TESIS DOCTORAL CAPÍTULO IV AERODINÁMICA DE UN PROYECTIL BASE BURN MODELIZACIÓN DE PROYECTILES BASE BURN Director: Dr. Coronel Francisco Cucharero Pérez Doctorando: Comandante Fernando Aguirre Estévez IV.-1 IV.-2 IV.-AERODINÁMICA DE UN PROYECTIL BASE BURN Existen diferentes tipos de resistencia aerodinámica. Básicamente, hay resistencias de origen potencial y de origen viscoso, no calculables mediante teoría potencial. El movimiento es irrotacional o potencial en una región del campo fluido cuando el rotacional de la velocidad es nulo en todos los puntos de la región. El teorema de Stokes asegura, entonces, la ausencia de circulación para toda línea cerrada, interior al dominio de irrotacionalidad, reducible a un punto por deformación continua de la misma. La circulación de la velocidad a lo largo de cualquier línea que une dos puntos cualesquiera del dominio es independiente de la línea, para líneas reducibles unas a otras por deformación continua. Si el dominio de irrotacionalidad es simplemente conexo, toda línea cerrada interior al dominio es reducible y el potencial de velocidades es unívoco. Si el dominio de irrotacionalidad es múltiplemente conexo (como es el caso de un cuerpo sumergido en un fluido), existirán líneas cerradas no reducibles que pueden tener circulación distinta de cero. El potencial de velocidades es, en ese caso, multiforme. La multiplicidad de conexión del dominio fluido irrotacional está ligada a la existencia de regiones de vorticidad no nula, como capas de torbellinos, que aparecen en la estela en el movimiento alrededor de cuerpos a grandes números de Reynolds. El teorema de Bjerkness-Kelvin dice que la variación en la unidad de tiempo de la circulación alrededor de una línea fluida cerrada, es igual a la circulación de la aceleración. Si la aceleración deriva de un potencial, la circulación de la aceleración a lo largo de la línea cerrada es nula y la circulación de la velocidad se mantiene constante. Si inicialmente, el movimiento es irrotacional, se mantiene irrotacional como es el caso de una corriente uniforme. Las resistencias de origen potencial se dividen en resistencia inducida, calculable mediante el plano de Trefftz, y la debida a las ondas de presión producidas por un móvil que se desplaza en régimen supersónico (resistencia de onda). IV.-3 El plano de Trefftz es un plano perpendicular a la corriente incidente no perturbada que se encuentra en zonas muy alejadas aguas abajo. El movimiento inducido por los torbellinos de la estela, en puntos muy alejados corriente abajo, de una superficie sustentadora es plano y la teoría potencial demuestra que la vorticidad se encuentra concentrada en esta estela. La resistencia inducida está relacionada con la huella que deja la superficie sustentadora en el plano de Treffz y se calcula, conocida la distribución de sustentación a lo largo de la envergadura. La interpretación física se debe a herraduras de torbellinos, cuyas cabezas están situadas en la superficie sustentadora y cuyas colas se extienden hasta el infinito formando la estela. La siguiente expresión refleja, con carácter general, la dependencia de la sustentación con la resistencia inducida, donde Λ es el alargamiento de la superficie sustentadora y τ un factor mayor o igual que cero (es cero para una distribución de fuerza lateral elíptica). C Di = C L2 (1 + τ ) π ⋅Λ De entre las resistencias de origen viscoso, la resistencia de rozamiento se debe a la viscosidad del fluido que desliza sobre la superficie del obstáculo y se adhiere a él, y es calculable usando la aproximación de capa límite mientras que la resistencia de forma o presión requiere un análisis más complejo, ya que en este caso, la capa límite se encuentra parcialmente desprendida, lo que obliga a recurrir a ensayos experimentales o análisis numéricos. El espacio del movimiento puede dividirse en dos: uno, de pequeño espesor, formado por el fluido que rodea al obstáculo, denominado, capa límite donde las fuerzas de inercia y las de viscosidad son del mismo orden, y otro formado por el fluido exterior al propio obstáculo, llamada región exterior, donde la viscosidad no cuenta y el fluido es ideal. En los fluidos perfectos o ideales, se desprecia la viscosidad al ser el número de Reynolds muy grande. Cuando se aplica esta teoría al cálculo de la sustentación o fuerza normal a la dirección del movimiento de un obstáculo, en el seno de una corriente fluida IV.-4 uniforme, los resultados obtenidos coinciden con los experimentales. Sin embargo, si se quiere calcular la resistencia o componente paralela a la dirección del movimiento, ésta teoría es inaplicable, ya que de acuerdo con la paradoja de D´Alembert la resistencia es nula para cualquier obstáculo finito. Es necesario tener en cuenta la capa límite. La presión sufre variaciones muy pequeñas a través de la capa límite, por lo que las presiones sobre el objeto pueden analizarse sin la capa límite. Pero la capa límite se encuentra sometida a gradientes de presión adversos fuertes que pueden provocar su desprendimiento. La distribución de presiones sobre el obstáculo se aproxima bien salvo en la zona de presión mínima, en la parte posterior. Las presiones medidas experimentalmente son inferiores a las teóricas. La resistencia debida a este fenómeno, es la denominada resistencia de forma, ya que depende de la forma del obstáculo. Esto se debe a que la corriente se desprende de la superficie del obstáculo dando lugar a una fuerte estela de torbellinos. La presión en la base del obstáculo influye fuertemente en la frecuencia de desprendimiento de torbellinos. Generalmente, al aumentar la presión de base, disminuye la frecuencia de desprendimiento. La estela tiene su origen en la capa límite desprendida. Analizando la corriente exterior, se presenta como una capa de torbellinos (discontinuidad tangencial de velocidades en una superficie fluida). La capa límite al separarse, se enrolla formando anillos o líneas de torbellinos que se desprenden alternativamente. En el límite, cuando el número de Reynolds crece hasta el infinito, el espesor de estas regiones es nulo. IV.-5 A continuación se va a considerar la resistencia de un proyectil sin resistencia inducida, es decir, el ángulo de ataque es el apropiado para que la sustentación sea nula. IV.-6 En las siguientes figuras puede observarse la zona de baja presión, en la parte posterior de un proyectil, responsable en gran parte de la resistencia de forma. El aumento de la presión en la base proporciona energía a la capa límite, retrasando su desprendimiento y disminuyendo la resistencia de forma o de base. IV.-7 La separación entre la zona turbulenta y potencial es clara, sugiriendo una forma inclinada en la parte posterior del proyectil. Una inclinación suave en la parte posterior del proyectil da como resultado el retraso en el desprendimiento de la capa límite, disminuyendo la zona turbulenta y a su vez la resistencia de base. Sin embargo, una inclinación muy grande en la parte posterior del proyectil puede ocasionar el desprendimiento de la capa límite, al no tener la corriente energía suficiente para acometer esa inclinación, lo cual resulta en un aumento de la región turbulenta y consecuentemente de la resistencia, perdiendo la inclinación por el culote la ventaja inicial. IV.-8 No obstante, si esta forma inclinada, en la parte posterior, se continuará indefinidamente terminando en un cono puntiagudo, se provocaría una reducción en la resistencia pero una fuerte inestabilidad del proyectil, dado que la capa límite terminaría desprendiéndose en algún punto, al estar sometida a un gradiente adverso de presión y los torbellinos se desprenderían alternativamente, produciendo sacudidas en el proyectil. De este modo, existe una configuración óptima de diseño de un proyectil: Base burn o base bleed constituye un procedimiento, mediante el cual, se expulsa algún tipo de gas, sin empuje, por la base del proyectil durante el vuelo. Este gas reduce la resistencia de base, disminuyendo el vacío creado en la parte posterior del proyectil e incrementando el alcance de una forma importante. La ausencia de propulsión adicional tiene la ventaja de que no aumenta la dispersión del proyectil como ocurre en los proyectiles RAP. IV.-9 De acuerdo a la reducción de resistencia, característica de un dispositivo base bleed, este proceso puede dividirse en tres fases: fase inicial transitoria, fase de quemado estable y fase de quemado final. El coeficiente de resistencia correspondiente a cada fase puede verse en la siguiente figura: IV.-10 IV.-11 IV.-12 ESCUELA POLITÉCNICA SUPERIOR DEL EJÉRCITO TESIS DOCTORAL CAPÍTULO V ANÁLISIS DE UN PROYECTIL BASE BURN MODELIZACIÓN DE PROYECTILES BASE BURN Director: Dr. Coronel Francisco Cucharero Pérez Doctorando: Comandante Fernando Aguirre Estévez V.-1 V.-2 V.-ANÁLISIS DE UN PROYECTIL BASE BURN La separación del flujo en la base de un proyectil conduce a la formación de una región de recirculación a baja velocidad, cerca de la base. La presión en esta región es inferior que la presión de la corriente libre. La resistencia de base causada por esta diferencia de presión puede ser superior a dos terceras partes de la resistencia total de un cuerpo de revolución. Técnicas como inducir una ligera inclinación por el culote, el base bleed o el base burn se han usado para reducir la resistencia de base; sin embargo, se han aplicado tradicionalmente de una forma empírica debido a la carencia de datos detallados o del conocimiento suficiente de las interacciones fluidodinámicas que ocurren en la base. La inyección de gases a baja velocidad, cerca de la estela del proyectil, reduce la resistencia de base de un proyectil. Esto provoca un incremento del alcance en los proyectiles de artillería. Esta técnica se denomina base bleed cuando no se considera reacción en la transferencia de masa. El concepto base burn se refiere a la inyección de masa por combustión de un propulsante situado en la base del proyectil. El propulsante se aloja en una cámara y la inyección de masa ocurre a través de un orificio de la cámara. Este agujero no es una tobera de forma que el empuje proporcionado por este dispositivo es muy pequeño. Los modelos numéricos son poco comunes en la determinación de flujos por la base, debido a su complejidad. Las medidas en laboratorio no son muy exactas, por lo que los sistemas que reducen la resistencia de base son difíciles de diseñar. Un esquema de una unidad base burn puede verse en las siguientes figuras, donde se esquematiza el dispositivo de ignición. Activado su encendido, el propulsante arde a través de una membrana y los gases son expulsados por un orificio de gran diámetro, practicado en la parte posterior del proyectil. Los gases arden en la cámara, durante un tiempo del orden de treinta segundos, y el chorro generado sale por el orificio a velocidad subsónica. V.-3 En la figura siguiente se muestra un esquema fluidodinámico de un flujo supersónico sobre un proyectil romo cilíndrico con base bleed. La corriente supersónica libre se expande en la esquina de la base, al producirse la separación de la capa límite, y apareciendo una capa libre de cortadura. Esta capa de cortadura se recomprime V.-4 realineándose y desarrollándose alrededor del eje de simetría. La capa de cortadura arrastra fluido de la zona posterior de la base, que es acelerado, retornando este fluido a la región de la base, formando una zona de recirculación en el proceso. La inyección de fluido a baja velocidad, por la base, desplaza los puntos de remanso aguas abajo, correspondiendo su localización a un equilibrio entre la cantidad de movimiento de gas inyectado y la cantidad de movimiento de fluido recirculado. La cantidad de fluido inyectado se cuantifica a través del parámetro de inyección adimensional I I= 4 ⋅ m& f π ⋅ d b2 ⋅ ρ ⋅ v donde m& f es la derivada de la masa de propulsante del motor, db el diámetro de la base del proyectil, ρ la densidad del aire y v la velocidad aerodinámica. Esta definición del parámetro de inyección no tiene en cuenta el espesor de la capa límite y la cantidad de movimiento del fluido inyectado, los cuales afectan a la presión en la base. La dependencia de la relación de presión en la base respecto a la presión de la corriente libre (Pb/P1), con el parámetro de inyección I muestra, desde un punto de vista experimental, tres regímenes distintos de operación determinados por la cantidad de masa inyectada. La relación Pb/P1 se incrementa casi linealmente para bajos valores de I V.-5 (régimen 1). Un pico en Pb/P1 aparece para un valor intermedio de I (I = 0.01 en el aire), el cual depende de varios factores, como el número de Mach de la corriente libre, el tamaño y geometría del orificio de salida del fluido de sangrado, el parámetro de inyección I, el peso molecular y la temperatura del gas soplado. Es posible obtener incrementos en el valor máximo Pb/P1 de entre 10 y 90 %. Sobrepasado el valor óptimo de Pb/P1, este disminuye (régimen 2) hasta alcanzar un mínimo relativo. Si se continua aumentando el flujo de sangrado se alcanzaría la zona supersónica (régimen 3), resultando en un incremento de Pb/P1. De los resultados experimentales, la efectividad de los dispositivos base-bleed se incrementa con el número de Mach de la corriente libre. A altos números de Mach, el pico en la presión de base ocurre para un parámetro de inyección I inferior, siendo el incremento en la presión de base superior. Para bajos regímenes de sangrado, el incremento en la presión de base con el sangrado es casi independiente del área de salida del gas, pero para regímenes altos de sangrado, la efectividad del dispositivo base bleed aumenta con el área de salida del gas. La inyección con bases porosas es también más efectiva. Análisis efectuados con aire, hidrógeno, helio, argón y nitrógeno han mostrado que el base bleed es más eficaz en la medida que el gas expelido tiene un peso molecular más V.-6 bajo. El pico en la presión de base es más alto y ocurre para valores inferiores del parámetro de inyección I, usando un gas más liviano. También se ha observado un incremento en la presión de base calentando el gas sangrado. Para bajos regímenes de inyección, la elevación de la presión de base es casi proporcional a la entalpía del gas. El pico de la presión de base es mayor y ocurre para valores inferiores del parámetro de inyección I, que para el caso correspondiente al sangrado del gas en frío. La inyección de hidrógeno en combustión ha mostrado presiones en la base, mayores que con sangrado en caliente, y el sangrado de gases con partículas combustibles sólidas es aún más efectivo. Aunque es conocida la efectividad del base bleed como una técnica de reducción de la resistencia, los detalles exactos en cuanto a interacciones fluidodinámicas causadas por el sangrado por el culote no son claramente comprendidos. A continuación pueden observarse fotografías del flujo de inyección, para diferentes parámetros de inyección I. Para sangrado nulo, se produce una fuerte onda de choque de recompresión cerca del punto de remanso. De I = 0 a I = 0.0033 el ángulo de la capa de cortadura se hace más pequeño, la expansión en la esquina de la base se debilita, la estela se hace más ancha y las ondas de recompresión se debilitan. La onda de recompresión se debilita, aún más, cerca de I = 0.0131 , cuando el flujo de sangrado proporciona más fluido que el que requiere la capa de cortadura. Cuando el flujo de sangrado se incrementa más allá de I = 0.0199, la onda de recompresión se desplaza ligeramente aguas arriba. Para I = 0.0279, la garganta se bloquea al alcanzar la velocidad de sangrado las condiciones sónicas. Una onda, que forma un disco de Mach, surge del orificio de sangrado e interactúa con las ondas de choque oblicuas de recomprensión, que vienen de la corriente exterior, formando una compleja interacción de ondas de choque. Este sistema es altamente inestable. Para I = 0.0148 se produce el sangrado óptimo, desde un punto de vista de presión de base, donde la región de recirculación casi desaparece. En este punto, la presión de base se maximiza, la estela se ensancha, el ángulo de la capa de cortadura se hace más pequeño y la velocidad de reversa, a lo largo de la línea central, se hace muy pequeña. Para I = 0.0226, no se detecta velocidad aguas arriba, a lo largo de la línea central, lo que indica la penetración del chorro de sangrado hacia la zona de readherencia. V.-7 V.-8 A continuación puede observarse un resultado computacional, resolviendo las ecuaciones de Navier Stokes, para el proyectil M864, para flujo másico caliente y frío, con valores del parámetro de inyección I superior a 0.04 . V.-9 V.-10 ESCUELA POLITÉCNICA SUPERIOR DEL EJÉRCITO TESIS DOCTORAL CAPÍTULO VI CÁLCULO NUMÉRICO DE UN PROYECTIL BASE BURN MODELIZACIÓN DE PROYECTILES BASE BURN Director: Dr. Coronel Francisco Cucharero Pérez Doctorando: Comandante Fernando Aguirre Estévez VI.-1 VI.-2 VI.-CÁLCULO NUMÉRICO DE UN PROYECTIL BASE BURN 1. Índice 1. Índice 2. Introducción 3. Modelo Modificado de Masa Puntual 4. Términos adicionales para proyectiles base burn método 1 5. Términos adicionales para proyectiles base burn método 2 6. Movimiento del proyectil referido al centro de gravedad 7. Vector de ataque 8. Cálculo numérico 8.1. Datos aerodinámicos 8.2. Datos base burn método 1 8.3. Datos base burn método 2 8.4. Factores de ajuste 8.5. Datos varios 8.6. Algoritmo de integración 8.7. Ordenada máxima 8.8. Condición de parada 8.9. Atmósfera 9. Presentación de resultados 9.1. Proyectil OE155F2RTC FR HEA BB1 9.2. Proyectil OE155BONUS FR CBL/DUP BB1 9.3. Proyectil M864 US CBL/DUP BB2 VI.-3 VI.-4 2. Introducción Para la confección de este trabajo es necesario disponer de un algoritmo de cálculo numérico, que sea capaz de describir matemáticamente la trayectoria balística de un proyectil de artillería base burn. En este punto, se mostrará el desarrollo matemático necesario para el cálculo de las diferentes trayectorias balísticas. Para ello se utilizará el Modelo Modificado de Masa Puntual, según se describe en el STANAG 4355, por el método 1 (método empleado por los países francófonos) y por el método 2 (método desarrollado por los EE.UU. y empleado por los países anglófonos). Con las bases matemáticas consideradas, se procederá a construir una serie de aplicaciones informáticas de cálculo en lenguaje C++, empleando un compilador Microsoft Visual Studio.NET. Las gráficas se construirán utilizando MATLAB R2006a. Estas herramientas informáticas se aplicarán al cálculo de las trayectorias balísticas de una serie de proyectiles de artillería. Proyectil calibre tubo carga MV QE OE155F2RTC FR HEA BB1 155 155AUF1 DLE-CE-155F1/7 C7 825.0 878.0 OE155BONUS FR CBL/DUP BB1 155 155 DLE-CE-B0/7 C7 DM52/8 CH8 808.0 698.0 891.3 OEF3BB NO HEA BB1 155AUF1 M109A3G M864 US CBL/DUP BB2 155 DM662 NO CBL BB2 155 M198 M109A3G M232/5H CH5H DM52/8 CH8 810.0 670.0 En el cuadro se describe la designación del proyectil, indicando la nacionalidad, el tipo de proyectil y el método de cálculo que emplea el país fabricante, el calibre, el tubo empleado en el disparo, la carga, la velocidad en boca del arma MV correspondiente a la carga referenciada y el ángulo de tiro QE que proporciona el máximo alcance. Aunque hay dos proyectiles que son portadores de submunición, el objeto de este apartado es el desarrollo de una capacidad matemática que permita calcular trayectorias balísticas individuales de diferentes proyectiles, por lo que a efectos académicos se considera, en estos proyectiles, una altura de quemado (HOB) de cero metros, considerándose como proyectiles rompedores a los efectos de interés. VI.-5 814.0 907.0 840.0 La experiencia acumulada en estos cálculos permite afirmar que los fenómenos de interés aparecen a medida que las trayectorias son más amplias en tiempo, por lo que el análisis presentado en este apartado se ha efectuado con la velocidad en boca más alta con la que se pueden disparar dichos proyectiles, y con el ángulo de tiro QE con el que se obtiene el máximo alcance. El proyectil noruego OEF3BB NO HEA BB1 se corresponde con el proyectil francés OE155F2RTC FR HEA BB1, pero disparado desde un tubo norteamericano M109A3G y con una carga norteamericana DM52/8 CH8. Análogamente, el proyectil noruego DM662 NO CBL BB2 se corresponde con el proyectil norteamericano M864 US CBL/DUP BB2. Aunque físicamente son los mismos proyectiles con los mismos coeficientes balísticos, algunos coeficientes se han adaptado a las estrictas condiciones climatológicas de Noruega. VI.-6 3. Modelo Modificado de Masa Puntual Las siguientes ecuaciones constituyen un modelo matemático representativo del vuelo de un proyectil estabilizado por rotación, dinámicamente estable (proyectil de artillería convencional) y con simetría axial. La modelización matemática incluye solamente las fuerzas y momentos esenciales despreciando los efectos transitorios. Todos los vectores se referirán a un sistema de ejes a derechas, ortonormal y situado fijo en el suelo. El sistema de coordenadas cartesianas con vectores unitarios se muestra en la siguiente figura. Las leyes de Newton referidas al centro de masas del proyectil son: r r r r r r r F = mu&& = DF + LF + MF + mg + mΛ siendo u la velocidad del proyectil. La aceleración debida a la fuerza de resistencia es: VI.-7 ( )) ⎛π ⋅ ρ ⋅i ⋅ d 2 ⎞ r DF ⎟⎟ C D0 + C D 2 QD ⋅ α e2 ⋅ vv = −⎜⎜ α m 8m ⎠ ⎝ ( siendo i y QD factores de ajuste, ρ la densidad del aire, d el diámetro del proyectil, m la masa del proyectil, αe el ángulo de ataque, v la velocidad aerodinámica, y, CD0 y CDα2 coeficientes aerodinámicos. La aceleración debida a la fuerza lateral es: ( ) r LF ⎛ π ⋅ ρ ⋅ d 2 ⋅ f L ⎞ ⎟⎟ C Lα + C L 3 α e2 ⋅ v 2α e = ⎜⎜ α m ⎝ 8m ⎠ r siendo fL un factor de ajuste denominado factor de sustentación, α e el vector de ataque, y, CLα y CLα3 coeficientes aerodinámicos. La aceleración debida a la fuerza Magnus es: 3 MF d ⋅ QM ⋅ p ⋅ ρ ⋅ Cmag − f r r (α e ∧ v ) = 8m m siendo QM un factor de ajuste denominado factor de la fuerza Magnus, p la velocidad de rotación axial del proyectil y Cmag-f un coeficiente aerodinámico. La aceleración debida a la gravedad es: ⎡ X1 ⎤ ⎢ R ⎥ ⎢ 2X ⎥ r r 2 g = − g 0 R 2 r 3 ⋅ r = − g 0 ⎢1 − ⎥ R ⎥ ⎢ ⎢ X3 ⎥ ⎢⎣ R ⎥⎦ ( ) g 0 = 9.80665[1 − 0.0026 cos(2 ⋅ lat )] VI.-8 r r r r = X −R ⎡ 0 ⎤ r ⎢ ⎥ R = ⎢− R ⎥ ⎢⎣ 0 ⎥⎦ y R = 6.356766 106 metros con X1, X2 y X3 las coordenadas del proyectil y lat la latitud. La aceleración debida a los efectos de la fuerza de Coriolis es: v r Λ = −2(ω ∧ u ) ⎡ Ω cos(lat ) cos( AZ ) ⎤ ⎥ Ω sin (lat ) ω = ⎢⎢ ⎥ ⎢⎣− Ω cos(lat ) sin ( AZ )⎥⎦ r Ω = 7.292115 10-5 rad/s siendo ω la velocidad de rotación de la tierra y AZ el acimut. Para el hemisferio Sur, lat es negativa. La magnitud de la aceleración de rotación axial del proyectil viene dada por: p& = π ⋅ ρ ⋅ d 4 ⋅ p ⋅ v ⋅ Cspin 8I x siendo IX el momento de inercia axial del proyectil y Cspin un coeficiente aerodinámico. El momento de inercia axial puede aproximarse por la siguiente expresión: VI.-9 I X = 0.14 ⋅ m ⋅ d 2 La velocidad de rotación axial del proyectil en el instante inicial viene dada por: p0 = 2 ⋅ π ⋅ u0 tc ⋅ d donde el proyectil da una vuelta en tc calibres. El vector de ataque viene dado por: ( ) r r 8 ⋅ I x ⋅ p v ∧ u& αe = − π ⋅ ρ ⋅ d 3 CM α + CM 3 ⋅ α e2 ⋅ v 4 r ( α ) ⎡0 ⎤ α e0 = ⎢⎢0⎥⎥ ⎢⎣0⎥⎦ r r siendo α e0 el valor inicial del vector de ataque, y, CMα y CMα3 coeficientes aerodinámicos. La velocidad del proyectil en el instante inicial viene dada por: ⎡u0 cos(QE ) cos( AZ )⎤ r ⎢ ⎥ u0 = ⎢ u0 sin (QE ) ⎥ ⎣⎢ u0 cos(QE )sin ( AZ )⎦⎥ La velocidad aerodinámica se obtiene de la expresión r r r v =u−w r siendo w la velocidad del viento. VI.-10 La posición del proyectil en el instante inicial se toma el origen del sistema de coordenadas. El número de Mach M=v/a siendo a la velocidad del sonido que se puede tomar como: a = 20.046 a (m/s) La corrección por esfericidad terrestre se obtiene de: X1 ⎡ ⎤ 2 2 ⎥ r ⎢ (X + X 3 )⎥ E = ⎢X2 + 1 2R ⎢ ⎥ X3 ⎢⎣ ⎥⎦ VI.-11 Los coeficientes aerodinámicos vienen dados como funciones del número de Mach. Esas funciones están dadas en forma de polinomios de cuarto grado o inferior. Se definen en regiones de números de Mach, desde M MAX i−1 hasta M MAX inclusive. Cada coeficiente aerodinámico se describe como polinomios de la forma: Ci = a0 + a1 ⋅ M + a2 ⋅ M 2 + a3 ⋅ M 3 + a4 ⋅ M 4 donde Ci es un coeficiente aerodinámico particular y M el número de Mach. Símbolo CD0 CDα2 CLα CLα3 CMα CMα3 Cmag-f Cspin Descripción Unidad Coeficiente de la fuerza de resistencia para ángulo de ataque cero Coeficiente cuadrático de la fuerza de 1/rad2 resistencia Coeficiente de la fuerza de 1/rad sustentación Coeficiente cúbico de la fuerza de 1/rad3 sustentación Coeficiente del momento de giro 1/rad Dependencia Mach Coeficiente cúbico del momento de giro Coeficiente de la fuerza Magnus 1/rad3 Mach 1/rad2 - Mach Coeficiente del momento amortiguamiento de rotación de Mach Mach Mach Mach Mach En orden a crear correspondencia entre los resultados obtenidos computacionalmente y las experiencias observadas, se definen los factores de ajuste. En cuanto a los factores de ajuste dependientes del ángulo de tiro, se tienen el factor de forma i y el factor de sustentación fL, que vienen dados para cada carga como: VI.-12 f = a0 + a1 ⋅ (QE ) + a 2 ⋅ (QE ) + a3 ⋅ (QE ) 2 3 siendo QE el ángulo de tiro. El factor de resistencia de ataque QD y el factor de la fuerza Magnus QM se toman como constantes. Parámetro Símbolo i Unidad - Factor de sustentación fL - 0.8 E j , se recurre al proceso de recocido. Se genera una perturbación aleatoria llamada p, que es un número real comprendido entre cero y uno, y se calcula la distribución de probabilidad para una determinada energía: ⎛ ΔE ⎞ ⎟⎟ P(ΔE ) = exp⎜⎜ − ⎝ kB ⋅T ⎠ Si P(ΔE ) > p , se acepta la nueva solución; en caso contrario, se desecha. Para la temperatura se ha tomado la ley de enfriamiento propuesta por Lundy y Mees: Ti = Ti −1 1 + βTi −1 La forma de descender la temperatura es la siguiente: ¾ Si se produce innovación, b se incrementa en una unidad y al valor de a se le asigna cero. ¾ Si no se produce innovación, b se mantiene y a se incrementa en una unidad. Cuando ocurre que a < A ó b > B, valores prefijados inicialmente, se produce la innovación en T y al valor de b se le asigna cero. El proceso se repite hasta alcanzar la condición de parada. Fijando la temperatura inicial T0, la temperatura final TM, y el número de iteraciones M, se define ß según la ecuación obtenida anteriormente β≈ T0 − TM TM ⋅ T0 ⋅ M Sin embargo, esta ecuación se ha obtenido descendiendo la temperatura en cada iteración y la temperatura desciende cuando se cumple la condición a < A ó b > B . IX.-18 Teniendo en cuenta los valores asignados a A y B, se tiene que A B = 10 con lo que estadísticamente el descenso de la temperatura ocurrirá con mayor frecuencia cuando se produzca el suceso B que cuando se produzca el A . A esto hay que añadirle que cuando no se produzca innovación, b se mantiene y no cambia de valor. Esta complicación debida a la aleatoriedad del método ha llevado a establecer un argumento aproximado para la estimación de ß. Si se asume que la temperatura siempre desciende cuando ocurre el suceso B y como A >> B se puede corregir la expresión anterior: β ≈ fac T0 − TM TM ⋅ T0 ⋅ M Como la temperatura desciende cuando b > B, el parámetro introducido debe tomar un cierto valor tal que fac > B . Tomando una cota conservativa, el parámetro ß se puede escribir β = (B + 1) T0 − TM TM ⋅ T0 ⋅ M De este modo, se fija la temperatura inicial T0 = 1, la temperatura final TM un cierto porcentaje de la temperatura inicial tal como 0.1-0.3, el número de iteraciones M, y el parámetro ß se obtiene de la ecuación anterior. La temperatura final obtenida tras el proceso de cálculo TF constituye una cota superior de TM, muy cercana a ella, lo cual impide que la temperatura final alcance un valor nulo lo que provocaría la aparición de infinitos. IX.-19 5. Constante de Boltzmann El punto más importante para que el algoritmo funcione de manera eficaz es la determinación de la constante de Boltzmann. Si la constante de Boltzmann es muy grande, P(ΔE ) tiende a la unidad y todos los cambios se aceptan con lo que el algoritmo pasa a ser un mero explorador de soluciones, pero si la constante de Boltzmann es muy pequeña, P(ΔE ) tiende a cero y no se acepta ningún cambio, con lo que se tiene un algoritmo de descenso incapaz de salir de mínimos locales. Un valor apropiado de la constante de Boltzmann es aquella que permite al algoritmo salir de óptimos locales para proceder a la búsqueda de otros óptimos mejores. En el primer momento de la iteración, se tiene la siguiente la ecuación ⎛ ΔE ⎞ ⎟⎟ P0 = exp⎜⎜ − ⋅ k T ⎝ B 0⎠ En concreto, si se retiene únicamente el primer valor de la iteración, se tiene: kB = − E0 ln P0 ⋅ T0 Si se retienen los dos primeros valores kB = − ΔE ln P0 ⋅ T0 Parece lógico que un valor apropiado de la constante de Boltzmann es aquella tal que P(ΔE ) es del orden de 0.5 , a fin de permitir al algoritmo salir de óptimos locales lo cual establece P0 = 0.5 . Generalmente, el problema es que el valor obtenido para la constante de Boltzmann en la primera iteración no resulta útil conforme el proceso avanza. El hecho de considerar una u otra, de las ecuaciones expuestas anteriormente, no representa gran diferencia. La adecuada elección de la constante de Boltzmann permite la exploración del espacio de configuraciones en los momentos iniciales, pues IX.-20 facilita que los saltos de una configuración a otra estén permitidos, pero cuando la temperatura se hace lo suficientemente baja, la probabilidad de salto a una configuración peor será prácticamente nula con lo que el proceso quedará atrapado en el mejor óptimo local. Esto ha provocado que se haya empleado un método más elaborado consistente en fijar un valor para la constante de Boltzmann y ejecutar el proceso de iteración un número fijo de iteraciones, almacenando P(ΔE ) para calcular ) P= ∑ P(ΔE ) número⋅total ⋅de⋅valores⋅ P ( ΔE ) número ⋅ total ⋅ de ⋅ valores ⋅ P (ΔE ) ) ) Como P ≈ 0.5 y dado que P(k B ) , es posible establecer un método simple que permita ) encontrar la constante de Boltzmann tal que P sea del orden de 0.5 , ya sea el método de la bisección, el de la secante, el de Newton, o cualquier otro. Dada la aleatoriedad que conlleva el método de recocido simulado no se necesita mucha precisión en la determinación de la constante de Boltzmann siendo un valor suficiente aquel que ) proporciona P ≈ 0.5 ± 0.1 , aunque si la precisión obtenida es únicamente del 10 %, el óptimo obtenido se puede considerar válido. El número de iteraciones no debe ser inferior a treinta aunque para tener una mayor confianza se recomienda un valor entre cincuenta y cien, si bien cuanto mayor sea este número más seguridad se tendrá en la localización del óptimo. Como valor inicial para empezar la búsqueda de la constante de Boltzmann, se puede emplear la ecuación mostrada anteriormente kB = − E0 ln P0 ⋅ T0 habiendo encontrado buenos resultados introduciendo P0 = 0.8-0.95 aunque no es una regla fija. Finalmente, una vez ejecutado el método y localizado el óptimo, una ) comprobación final es obtener el factor P ≈ 0.5 . En la inmensa mayoría de los casos analizados, el valor obtenido es satisfactorio. Un inconveniente que presenta el algoritmo de recocido simulado es que, además del IX.-21 problema que supone la determinación de la constante de Boltzmann, la variación de cualquier parámetro conlleva tener que volver a recalcularla. En ocasiones, con una ligera modificación es suficiente pero en otras es necesario casi empezar desde el principio. La experiencia acumulada en este campo permite afirmar que la correcta determinación de la constante de Boltzmann permite llegar a obtener valores óptimos con una precisión muy elevada. IX.-22 6. Condición de parada La condición de parada determina cuando se llega al equilibrio térmico y finaliza el algoritmo. En teoría, el algoritmo debería finalizar cuando T = 0. En la práctica, se para cuando T alcanza un cierto valor final, previamente especificado, o después de cumplir alguna de las condiciones que se enumeran a continuación. Una buena opción es parar cuando el número de movimientos aceptados es inferior a un cierto porcentaje. En ese caso, es muy probable que el algoritmo se haya estancado y no vaya a mejorar la solución obtenida. Se van a establecer cuatro condiciones de parada que detienen el algoritmo si se cumple una cualquiera de ellas. 6.1. Condición 1 Si la solución no ha mejorado tras K1 series consecutivas de L pasos, una cierta cantidad definida por ∈1 > 0 . E (K1 ⋅ L + n ) > E (n ) ⋅ (1+ ∈1 ) ⇒ FIN Lundy (1985) se ha basado en este criterio, demostrando que el número de iteraciones para llegar a la temperatura final es polinomial. 6.2. Condición 2 Si el número de movimientos aceptados en K2 series consecutivas de L pasos, es inferior que ∈2 ⋅L ⇒ FIN , con ∈2 > 0 . 6.3. Condición 3 Si la temperatura desciende por debajo de la temperatura final. IX.-23 6.4. Condición 4 Si el número de iteraciones alcanza un valor superior al valor máximo establecido M. IX.-24 7. Cálculo del siguiente paso El siguiente valor se obtendría de una forma aleatoria. Hay investigadores que proponen que se modifiquen todas las componentes del vector de estado, de forma que sean aleatorias puras, aunque también se pueden modificar solamente algunas de ellas. Para estudiar este fenómeno se va a calcular un problema tipo con diferentes configuraciones, y se van a comparar las soluciones obtenidas. Así, se van a localizar los factores de ajuste de un problema balístico para un proyectil base burn, empleando el método 1 del Modelo Modificado de Masa Puntual. Utilizando el algoritmo de recocido simulado se van a determinar los cuatro factores de ajuste tDI, K(p) (función lineal) y f(iBB,MT) para un ángulo de tiro concreto. Se ha escogido para ello el proyectil OE155F2RTC FR HEA BB1, disparado desde un tubo 155AUF1 con una carga C7, por ser un proyectil muy estudiado del que se dispone abundante información calculada mediante el método 1. Dado que en este algoritmo la velocidad de cálculo es muy importante para obtener resultados de forma rápida, se ha recurrido al método 1 por permitir una mayor velocidad en la computación. Inicialmente, se realizaron una serie de cálculos empleando un ángulo de tiro de 300 milésimas pero se ha podido comprobar que la dispersión de resultados aumenta en la medida que se emplean ángulos de tiro superiores. Posteriormente, se podrá ver que para mejorar el proceso de optimización del algoritmo, interesa tener una mayor dispersión en los resultados que permitan ver claramente las diferencias de una configuración a otra. De igual forma, se ha escogido la carga más elevada por presentar una mayor dispersión en los datos obtenidos. Aumentando el ángulo de tiro, se han encontrando muy pocas diferencias entre 700 y 900 milésimas por lo que finalmente los cálculos se han efectuado con un ángulo de tiro de 700 milésimas, al objeto de evitar efectos asociados al segundo sector que podrían enmascarar y confundir los resultados obtenidos, y también para reducir el tiempo de cálculo. Para cada carga, de los cuatro factores de ajuste, únicamente el factor base burn f(iBB,MT) es función del ángulo de tiro siendo el tiempo de retraso a la ignición del motor base burn tDI función de la temperatura de propulsante del motor MT, y el factor de quemado con la velocidad de rotación axial del proyectil K(p) una función lineal de IX.-25 la velocidad de rotación axial del proyectil p. En este momento se va a trabajar con un único ángulo de tiro y no es de aplicación pero se debe tener en cuenta que mientras que f(iBB,MT) es función del ángulo de tiro, tDI y K(p) no lo son. Obtenido el alcance y la deriva de experiencias realizadas, se van a calcular los factores de ajuste tDI, K(p) y f(iBB,MT) que mejor aproximen el alcance y la deriva. Se han escogido estos parámetros por ser los más representativos del tiro de artillería. La adición a los anteriores de la ordenada máxima ha permitido constatar un empeoramiento en el cálculo de los factores de ajuste. Si se obtienen unos factores de ajuste que mejor ajusten alcance y deriva y unos segundos factores de ajuste que mejor ajusten alcance, deriva y ordenada máxima; se ha observado que si se calcula alcance y deriva con los segundos factores de ajuste, se obtienen peores valores que usando los primeros factores de ajuste. Esto es debido a que la ordenada máxima es un parámetro que ante variaciones en los factores de ajuste, experimenta variaciones mucho mayores que las apreciadas en alcance y deriva. En concreto, las variaciones en la ordenada máxima pueden llegar a ser del orden de diez veces superiores a las encontradas para alcance y deriva. También se ha constatado la dificultad por parte del Modelo Modificado de Masa Puntual de proporcionar valores acertados de la ordenada máxima, una vez ajustada la deriva y el alcance, con datos extraídos de experiencias reales. Finalmente, se ha optado por ajustar únicamente alcance y deriva. Si a través de experiencias realizadas se obtiene el alcance y la deriva, que se designaran como alcanceexp y derivaexp, se define la función objetivo como: E= (alcance − alcance ) + (deriva − deriva ) 2 exp 2 exp siendo los términos alcance y deriva los obtenidos como resultado de haber realizado todo el proceso de cálculo en función de los factores de ajuste tDI, K(p) y f(iBB,MT). Una característica muy importante de este problema es que se conoce, a priori, el valor óptimo buscado de la función objetivo. Como condiciones iniciales, se considera que para el tiempo de retraso a la ignición del motor base burn tDI se introducirá el valor correspondiente, en el supuesto de que se conozca. Si no se dispone de información adicional, se considerará tDI = 0 siendo el resto de valores K(p) = 1 y f(iBB,MT) = 1 . IX.-26 Como ya se comentó anteriormente, el algoritmo de integración utilizado ha sido el Runge Kutta 45. Los parámetros eRK, hmin y hmax se han elegido para obtener una precisión inferior a un metro en el cálculo de la deriva y el alcance: hmin = 0.04 s, hmax = 6.0 s y eRK = 0.1 . Se ha considerado apropiada para este problema la siguiente condición de parada; L = 100, K1 = 4, ∈1 = 0.5 , K2 = 4, ∈2 = 0.5 y M = 1000. En el peor de los casos, cada proceso de cálculo independiente finaliza cuando se alcanza el número máximo de iteraciones M, que teniendo en cuenta el tiempo de ejecución de ordenador, se ha considerado suficiente fijarlo en 1000. A fin de comprobar la calidad de los resultados obtenidos se va a establecer un sistema en el cual, tras alcanzar la condición de parada y extraer el óptimo correspondiente, se va a continuar calculando hasta volver a alcanzar una nueva condición de parada con su correspondiente nuevo óptimo, y así sucesivamente. Así, se han ejecutado veinte de estos procesos de cálculo diferentes, comprobando como cuando el algoritmo alcanza el óptimo puede hacerlo con diferentes valores de factores de ajuste; es decir, aunque el valor del óptimo sea muy similar, existen diferentes combinaciones de los factores de ajuste que alcanzan ese valor del óptimo. A fin de investigar este fenómeno, se ha considerado realizar otros diez procesos de cálculo fijando tDI = 0 , no obteniendo grandes diferencias con respecto al caso anterior. Finalmente, se han calculado otros diez procesos de cálculo fijando tDI = 0 y K(p) una función constante, sin considerar la dependencia con la velocidad de rotación axial del proyectil p, mostrando estos últimos resultados, como en los diez procesos de cálculo se ha alcanzado aproximadamente el mismo óptimo con similares factores de ajuste. Resumiendo, se ha establecido un número de procesos cálculo de cuarenta, los cuales se utilizaran posteriormente a efectos de valorar los resultados proporcionados por el algoritmo. A partir del valor inicial es necesario definir el tamaño de paso para cada una de las variables, que permita saltar de una configuración a otra. Con carácter general, este paso se ha definido como Δ = f esc ⋅ Pre ⋅ T ⋅ Rnd ⋅10 IX.-27 siendo fesc el factor de escala, Pre la precisión, T la temperatura y Rnd un número aleatorio. Este tamaño de paso debe ser un orden de magnitud menor que el tamaño característico de óptimo local. Un tamaño mayor no permitirá al algoritmo funcionar como algoritmo de optimización convirtiéndose, en un mero explorador de soluciones. Esta es la idea inherente al introducir el parámetro precisión. En este problema se ha encontrado una buena aproximación para Pre = 10e-4 . Este tamaño puede ser diferente dependiendo de la variable a la que vaya a afectar ya que la escala de cada una de las dimensiones no tiene porque ser la misma. Esto conlleva a que es necesario definir el factor de escala que afecta a cada una de las variables. En cada iteración se va reduciendo la temperatura y por tanto las probabilidades son cada vez más pequeñas conforme avanza el procedimiento y se acerca a la solución óptima, con lo que el tamaño del paso también disminuye. Inicialmente se realiza una diversificación de la búsqueda sin controlar demasiado el coste de las soluciones visitadas. Cuando la temperatura se haga lo suficientemente baja la probabilidad de salto a una configuración peor será prácticamente nula, con lo que el proceso quedará atrapado en el mejor óptimo local conseguido hasta el momento y el tamaño de paso será muy pequeño. Para obtener el siguiente valor es necesario modificar alguna de las componentes del vector de estado. Un método clásico de hacer esto es tomar el vector de estado, que en este caso tiene cuatro componentes; tDI, K(p) y f(iBB,MT), realizar un sorteo entre todas ellas y decrementar la elegida una cantidad dada Δ , con Rnd comprendido entre 0 y 1 .A continuación se realiza un nuevo sorteo entre las componentes restantes, que no han variado, y se incrementa la elegida un nuevo valor Δ , con Rnd comprendido entre 0 y 1 . Esto permite conocer la calidad de soluciones que puede proporcionar el algoritmo y comparar con otro método para obtener el siguiente valor, como puede ser elegir aleatoriamente una componente del vector de estado y adicionarle a ésta una cantidad dada Δ , con Rnd comprendido entre -1 y 1 . IX.-28 En las siguientes gráficas se van a comparar las dos configuraciones descritas anteriormente, mostrando alcance y deriva. Se ha designado como SA2 aquella que modifica dos componentes del vector de estado y SA1 aquella que modifica una sola componente. Definiendo el área existente entre el valor calculado y el valor exacto para SA2 y SA1, y calculando la diferencia entre las dos áreas, puede estimarse que valor calculado se acerca más al valor exacto. Las gráficas de alcance y deriva muestran como es posible obtener valores realmente próximos al valor buscado, encontrándose precisiones del orden de un metro. IX.-29 Otra opción es modificar todas las componentes del vector de estado adicionando a cada una de ellas una cantidad dada Δ , con Rnd comprendido entre -1 y 1, y compararla con la descrita anteriormente que modifica únicamente dos componentes del vector de estado. Se ha designado SA2 aquella que modifica dos componentes del vector de estado y SA4 aquella que modifica todas las componentes. Un problema observado cuando se modifican todas las componentes, es que una vez introducida la constante de Boltzmann, ejecutado el método y localizado el óptimo, se ) ) pasa a comprobar el factor P encontrándose que en algunas ocasiones P resultaba del orden de 0.99 , como es el caso del punto número 11 de la iteración donde puede observarse un error enorme. Este punto se ha incluido a efectos de observar lo que ocurre cuando la constante de Boltzmann no se ha elegido apropiadamente. El problema es que esto ha ocurrido sobre un 15 % de las veces aunque en la gráfica solo se ha incluido en el punto número 11. Este inconveniente se agudiza al aumentar el tamaño del vector de estado. Así, se puede afirmar que el efecto de modificar todas las componentes del vector de estado provoca una mayor dificultad en la determinación de la constante de Boltzmann. IX.-30 IX.-31 Una cuestión que podría plantearse es que ocurriría si solo se retienen aquellos puntos ) tales que P ≈ 0.5 . En este caso, el alcance ha salido ligeramente superior para SA4 mientras que la deriva resulta algo mejor para SA2. La dificultad en la determinación de la constante de Boltzmann supone un grave inconveniente añadido y la mejora obtenida solo se ha manifestado muy levemente en la determinación del alcance, por lo que el efecto de modificar todas las componentes del vector de estado supone más un inconveniente que una ventaja. IX.-32 Como resumen, se puede afirmar que este algoritmo requiere una labor de determinación (y por tanto de prueba y error) muy exhaustiva de los parámetros. IX.-33 IX.-34 ESCUELA POLITÉCNICA SUPERIOR DEL EJÉRCITO TESIS DOCTORAL CAPÍTULO X ALGORITMOS GENÉTICOS MODELIZACIÓN DE PROYECTILES BASE BURN Director: Dr. Coronel Francisco Cucharero Pérez Doctorando: Comandante Fernando Aguirre Estévez X.-1 X.-2 X.-ALGORITMOS GENÉTICOS 1. Índice 1. Índice 2. Introducción 3. Definición 4. Codificación 5. Selección 6. Mutación 7. Cruce 8. Reducción 9. Población 10. Condición de parada 11. Implementación del algoritmo X.-3 X.-4 2. Introducción Esta técnica se basa en los mecanismos de selección que utiliza la naturaleza, de acuerdo a los cuales los individuos más aptos de una población son los que sobreviven, al adaptarse más fácilmente a los cambios que se producen en su entorno. Hoy en día se sabe que estos cambios se efectúan en los genes de un individuo (unidad básica de codificación de cada uno de los atributos de un ser vivo), y que sus atributos más deseables (los que le permiten adaptarse mejor a su entorno) se transmiten a sus descendientes cuando éste se reproduce sexualmente. Una definición bastante completa de un algoritmo genético es la propuesta por Koza: "Es un algoritmo matemático altamente paralelo que transforma un conjunto de objetos matemáticos individuales con respecto al tiempo usando operaciones modeladas de acuerdo al principio Darwiniano de reproducción y supervivencia del más apto, y tras haberse presentado de forma natural una serie de operaciones genéticas de entre las que destaca la recombinación sexual. Cada uno de estos objetos matemáticos suele ser una cadena de caracteres (letras o números) de longitud fija que se ajusta al modelo de las cadenas de cromosomas, y se les asocia con una cierta función matemática que refleja su aptitud. " Los primeros desarrollos en el campo de los algoritmos genéticos aparecieron a finales de los cincuenta y principios de los sesenta, programados en computadoras por biólogos evolutivos que buscaban explícitamente realizar modelos de aspectos de la evolución natural. En 1962, investigadores como Box, Friedman, Bledsoe y Bremermann habían desarrollado independientemente algoritmos inspirados en la evolución para optimización de funciones, pero sus trabajos generaron poca reacción. En 1965, Rechenberg de la Universidad Técnica de Berlín introdujo una técnica que él denominó estrategia evolutiva. En esta técnica no había población ni cruzamiento; un padre mutaba para producir un descendiente, y se conservaba el mejor de los dos, convirtiéndose en el padre de la siguiente ronda de mutación. Versiones posteriores introdujeron la idea de población. X.-5 En 1966, Fogel, Owens y Walsh introdujeron una técnica que llamaron programación evolutiva. En este método, las soluciones candidatas para los problemas se representaban como máquinas de estado finito. Al igual que en la estrategia evolutiva de Rechenberg, su algoritmo funcionaba mutando aleatoriamente una de estas máquinas simuladas y conservando la mejor de las dos. Sin embargo, todavía faltaba introducir la idea de cruce. En 1962, el trabajo de Holland sobre sistemas adaptativos estableció las bases para desarrollos posteriores, y fue el primero en proponer explícitamente el cruce y otros operadores de recombinación. El trabajo fundamental en el campo de los algoritmos genéticos apareció en 1975, con la publicación del libro “Adaptación en Sistemas Naturales y Artificiales”. Basado en investigaciones y artículos anteriores del propio Holland y de colegas de la Universidad de Michigan, este libro fue el primero en presentar la mutación, la selección y el cruce, simulando el proceso de la evolución biológica como estrategia para resolver problemas. Al principio, estas aplicaciones eran principalmente teóricas. Sin embargo, al seguir proliferando la investigación, los algoritmos genéticos se aplicaron en otros sectores en la medida que los ordenadores incrementaban su potencia de cálculo. Actualmente, la computación evolutiva es un campo emergente, y los algoritmos genéticos están resolviendo problemas de interés. No hay más que recordar las palabras de Charles Darwin: “que el azar en la variación, junto con la ley de la selección, es una técnica de resolución de problemas de inmenso poder y de aplicación casi ilimitada”. X.-6 3. Definición Los algoritmos genéticos (GA´s) son métodos adaptativos que pueden usarse para resolver problemas de búsqueda y optimización. Están basados en el proceso genético de los organismos vivos. A lo largo de las generaciones, las poblaciones evolucionan en la naturaleza de acuerdo con los principios de la selección natural y supervivencia de los más aptos, postulados por Darwin, donde solo los mejor adaptados tienen mayores posibilidades de reproducirse. Por imitación de este proceso, los algoritmos genéticos son capaces de ir creando soluciones para problemas del mundo real. Dado un problema específico a resolver, la entrada del GA es un conjunto de soluciones potenciales a ese problema, codificadas de alguna manera, y una métrica llamada función de aptitud que permite evaluar cuantitativamente a cada candidata. Estas candidatas pueden ser soluciones que ya se sabe que funcionan, con el objetivo de que el GA las mejore, pero que suelen generarse aleatoriamente. Inmediatamente, el GA evalúa cada candidata de acuerdo con la función de aptitud. En el conjunto de candidatas, generadas aleatoriamente, la mayoría no funcionarán en absoluto, y serán eliminadas. Sin embargo, algunas pueden resultar interesantes o pueden mostrar actividad, aunque sólo sea actividad débil e imperfecta, hacia la solución del problema. Estas candidatas prometedoras se conservan y se les permite reproducirse. Se realizan múltiples copias de ellas, pero las copias no son perfectas; se introducen cambios aleatorios durante el proceso de copia. Luego, esta descendencia digital prosigue con la siguiente generación, formando un nuevo conjunto de soluciones candidatas, y son sometidas a una ronda de evaluación de aptitud. Las candidatas que han empeorado o no han mejorado con los cambios en su código son eliminadas otra vez; pero, de nuevo, las variaciones aleatorias introducidas en la población pueden haber mejorado a algunos individuos, convirtiéndolos en mejores soluciones del problema, más completas o más eficientes. Nuevamente, se seleccionan y copian estos individuos vencedores hacia la siguiente generación con cambios aleatorios, y el proceso se repite. Las expectativas son que la aptitud media de la población se incrementará en cada ronda y, por tanto, repitiendo este proceso cientos o miles de rondas, pueden descubrirse soluciones muy buenas del problema. X.-7 Los algoritmos genéticos han demostrado ser una estrategia enormemente poderosa para resolver problemas, demostrando de manera espectacular el poder de los principios evolutivos. Aunque no se garantiza que el algoritmo genético encuentre la solución óptima del problema, se encuentran soluciones de un nivel aceptable en un tiempo competitivo comparado con otros algoritmos. Una particularidad de los algoritmos genéticos es que son paralelos. La mayoría de los otros algoritmos son serie y sólo pueden explorar el espacio de soluciones hacia una solución en una dirección al mismo tiempo, y si la solución que descubren resulta subóptima, no se puede hacer otra cosa que abandonar todo el trabajo hecho y empezar de nuevo. Sin embargo, dado que los GA´s tienen descendencia múltiple, pueden explorar el espacio de soluciones en múltiples direcciones a la vez. Si un camino resulta ser poco útil, se puede eliminar fácilmente y continuar el trabajo en otras direcciones. Esta idea puede verse desde otro punto de vista. Considérese el concepto de esquema o patrón que describe un conjunto cromosomas. Para ser preciso un esquema es una cadena compuesta de ceros, unos y asteriscos. Cada esquema representa al conjunto de cromosomas que contienen los ceros y los unos en la misma posición en la que se encuentran en el esquema y cualquier valor en la posición de los asteriscos. Sean todas las cadenas binarias (cadenas de ceros y unos) de n dígitos que forman un espacio de búsqueda, que puede representarse como ******** (donde * significa 0 o 1). La cadena 01101010 es un miembro de este espacio pero también es un elemento del espacio 0*******, del espacio 01******, del espacio 0******0, del espacio 0*1*1*1*, del espacio 01*01**0, etc. Evaluando esta cadena concreta, un algoritmo genético estaría explorando cada uno de los espacios a los que pertenece. Tras muchas evaluaciones, iría obteniendo un valor cada vez más preciso de la aptitud media de cada uno de estos espacios, cada uno de los cuales contiene muchos elementos. Por tanto, un algoritmo GA que evalúe explícitamente un número pequeño de individuos, está evaluando implícitamente un grupo de individuos mucho más grande. Este concepto se conoce como teorema del esquema en el campo de los algoritmos evolutivos. Debido al paralelismo que les permite evaluar implícitamente muchos esquemas a la vez, los algoritmos genéticos funcionan particularmente bien resolviendo problemas cuyo espacio de soluciones potenciales es muy amplio como para hacer una búsqueda X.-8 exhaustiva en un tiempo razonable. La mayoría de los problemas que caen en esta categoría se conocen como no lineales. En un problema lineal, la aptitud de cada componente es independiente, por lo que cualquier mejora en alguna parte dará como resultado una mejora en el sistema completo. La no linealidad es más habitual en la vida real, donde cambiar una componente puede tener efectos en todo el sistema, y donde cambios múltiples que por separado son perjudiciales, en combinación pueden conducir hacia mejoras en la aptitud. Los algoritmos genéticos son especialmente útiles en problemas con un paisaje adaptativo complejo, aquéllos en los que la función de aptitud es discontinua, ruidosa, cambia con el tiempo, o tiene muchos óptimos locales. La mayoría de los problemas reales tienen un espacio de soluciones enorme que es imposible de explorar exhaustivamente. En este caso, se deben tener herramientas que sean capaces de obviar los óptimos locales y no quedar atrapados en ellos, para acometer la búsqueda del óptimo global. El cruce es el factor esencial de los algoritmos genéticos. Sin el cruce, cada solución individual explora el espacio de búsqueda sin información de lo que el resto de individuos puedan haber descubierto. El cruce permite transmitir información entre los mejores individuos para producir una descendencia que tenga las virtudes de sus padres y ninguna de sus debilidades. En ocasiones es necesario modelizar un problema en términos de múltiples objetivos. Los GA´s presentan una estructura que se adapta muy bien a estos problemas donde una solución optimiza un parámetro y otra optimiza otro distinto. Cuando una solución particular a un problema, con múltiples objetivos, optimiza un parámetro pero ese parámetro no puede mejorarse más sin causar una pérdida de calidad en otro parámetro, esa solución se denomina óptimo paretiano o no dominado. Al crear un algoritmo genético se debe realizar una representación del problema. Aunque se pueden idear diferentes formas de hacer esto, lo más común es definir a los individuos como cadenas de números binarios, enteros o reales; donde cada número representa algún aspecto determinado de la solución. Si los individuos son cadenas binarias, un 0 o 1 podría significar la ausencia o presencia de una cierta característica. Si X.-9 son listas de números, estos números podrían representar diferentes cosas. Así, la mutación implica cambiar estos números, cambiar bits o sumar o restar valores aleatorios. La definición de la función de aptitud debe ser analizada cuidadosamente para conseguir una mejora en la resolución del problema en cuestión. Si se elige mal la función de aptitud o se define de manera inexacta, puede que el algoritmo genético sea incapaz de encontrar una solución al problema, o puede acabar resolviendo un problema equivocado. El algoritmo genético únicamente maximiza, pero la minimización puede realizarse fácilmente utilizando el recíproco de la función maximizante. Además de elegir bien la función de aptitud, también deben elegirse cuidadosamente los otros parámetros de un GA; el tamaño de la población, el ritmo de mutación, el cruce y el tipo de selección. Si el tamaño de la población es demasiado pequeño, puede que el algoritmo genético no explore suficientemente el espacio de soluciones para encontrar buenas soluciones consistentemente. Si el ritmo de cambio genético es demasiado alto o el sistema de selección se escoge inadecuadamente, puede alterarse el desarrollo de esquemas beneficiosos y la población cambia demasiado deprisa antes de que la selección llegue a producir convergencia. Si el tamaño de una población es pequeño, los ritmos de mutación son muy altos o la presión selectiva es demasiado fuerte (una situación así podría ser resultado de un cambio ambiental drástico), entonces la especie puede extinguirse. La solución ha sido la evolución de la evolutividad, las adaptaciones que alteran la habilidad de una especie para adaptarse. La mayoría de los seres vivos comprueba y corrige errores, a nivel celular, durante el proceso de replicación del ADN, manteniendo su ritmo de mutación a unos niveles aceptablemente bajos. En otro sentido; en tiempos de fuerte presión ambiental, algunas especies de bacterias entran en un estado de hipermutación en el que el ritmo de errores en la replicación del ADN aumenta bruscamente, aumentando la probabilidad de que se descubra una mutación compensatoria. Un problema que puede surgir con un GA se conoce como convergencia prematura (los biólogos la llaman deriva genética). Si un individuo, que es más apto que la mayoría de sus competidores, emerge muy pronto en el curso de la ejecución, se puede reproducir X.-10 tan abundantemente que merme la diversidad de la población demasiado pronto, provocando que el algoritmo converja hacia el óptimo local que representa ese individuo, en lugar de rastrear el paisaje adaptativo lo bastante a fondo para encontrar el óptimo global. Esto es un problema especialmente común en las poblaciones pequeñas, donde incluso una variación aleatoria en el ritmo de reproducción puede provocar que un genotipo se haga dominante sobre los otros. Los métodos más comunes implementados por los investigadores en GA´s para solucionar este problema implican controlar la fuerza selectiva para no proporcionar tanta ventaja a los individuos excesivamente aptos. La selección escalada, por rango o por torneo, son tres de los métodos principales para conseguir esto. Algunos métodos de selección escalada son el escalado sigma, en el que la reproducción se basa en una comparación estadística de la aptitud media de la población y la selección obtenida por la función de Boltzmann, en la que la fuerza selectiva aumenta durante la ejecución controlando de una manera similar a como lo hace la variable temperatura en el recocido simulado. X.-11 4. Codificación Los individuos (posibles soluciones del problema), pueden representarse como un conjunto de parámetros (que se denominan genes), los cuales agrupados forman un conjunto de valores (a menudo referida como cromosoma). El conjunto completo de material genético se denomina genoma. La evolución es un proceso que opera a nivel de cromosomas y no a nivel de individuos. Cada individuo es codificado como un conjunto de cromosomas por lo que se necesita un método para codificar las soluciones potenciales del problema, de forma que un ordenador pueda procesarlos. Un enfoque común es codificar las soluciones como cadenas binarias: secuencias de 1s y 0s, donde el dígito de cada posición representa el valor de algún aspecto de la solución. Otro método similar consiste en codificar las soluciones como cadenas de enteros o números decimales, donde cada posición, de nuevo, representa algún aspecto particular de la solución. Este método permite una mayor precisión y complejidad que el método comparativamente restringido de utilizar sólo números binarios, y a menudo está intuitivamente más cerca del espacio de problemas; sin embargo, la teoría en la que se fundamentan los algoritmos genéticos se comprende mejor desde un punto de vista binario. Un tercer método consiste en representar a los individuos de un GA como cadenas de letras, donde cada letra, otra vez, representa un aspecto específico de la solución. Esta técnica se ha usado para generar redes neuronales. La virtud de estos tres métodos es que facilitan la definición de los operadores que provocan los cambios aleatorios en los candidatos seleccionados: cambiar un 0 por un 1 o viceversa, sumar o restar al valor de un número una cantidad elegida al azar, o cambiar una letra por otra. Otra estrategia, desarrollada principalmente por Koza y denominada programación genética, representa a los programas como estructuras de datos ramificadas llamadas árboles. En este método, los cambios aleatorios pueden generarse cambiando el operador o alterando el valor de un cierto nodo del árbol, o bien sustituyendo un subárbol por otro. En términos biológicos, el conjunto de parámetros representando un cromosoma particular se denomina fenotipo. El fenotipo contiene la información requerida para construir un organismo, el cual se refiere como genotipo. Los mismos términos se X.-12 utilizan en el campo de los algoritmos genéticos. La adaptación al problema de un individuo depende de la evaluación del genotipo que puede obtenerse a partir del fenotipo. La función de adaptación debe ser diseñada para cada problema de manera específica. Dado un cromosoma particular, la función de adaptación le asigna un número real, que se supone refleja el nivel de adaptación al problema del individuo representado por ese cromosoma. X.-13 5. Selección Durante la fase reproductiva se seleccionan los individuos de la población para cruzarse y producir descendientes, que constituirán, una vez mutados, la siguiente generación de individuos. Una vez seleccionados los padres (generalmente dos), sus cromosomas se combinan utilizando habitualmente los operadores de cruce y mutación. Un algoritmo genético puede utilizar muchas técnicas diferentes para seleccionar a los individuos que deben copiarse hacia la siguiente generación, pero abajo se indican algunos de los más comunes. Algunos de estos métodos son mutuamente exclusivos, pero otros pueden utilizarse en combinación. ¾ Selección elitista: los miembros más aptos de cada generación son seleccionados (La mayoría de los GA´s no utilizan elitismo puro, sino que usan una forma modificada por la que el individuo mejor, o algunos de los mejores, son copiados hacia la siguiente generación, en caso de que no surja nada mejor). Esta técnica tiene el problema de que no permite al algoritmo salir de óptimos locales para buscar óptimos globales. ¾ Selección proporcional a la aptitud: los individuos más aptos tienen más probabilidad de ser seleccionados, pero no la certeza. ¾ Selección por rueda de ruleta: forma de selección proporcional a la aptitud en la que la probabilidad de que un individuo sea seleccionado es proporcional a la diferencia entre su aptitud y la de sus competidores. Conceptualmente, esto puede representarse como un juego de ruleta; cada individuo obtiene una sección de la ruleta, pero los más aptos obtienen secciones mayores que las de los menos aptos. Luego la ruleta se hace girar y cada vez se elige el individuo que posea la sección en la que se para la ruleta. ¾ Selección escalada: al incrementarse la aptitud media de la población, la fuerza de la presión selectiva también aumenta y la función de aptitud se hace más discriminadora. La técnica propuesta por Goldberg traza una recta entre el X.-14 menor valor de aptitud y el mayor valor de aptitud de la población, y luego prolonga o trunca los valores de aptitud de cada individuo al valor de la recta que le corresponda. Este método puede ser útil para seleccionar más tarde, cuando todos los individuos tengan una aptitud relativamente alta y sólo les distingan pequeñas diferencias en la aptitud. ¾ Selección por torneo: se eligen subgrupos de individuos de la población al azar, y los miembros de cada subgrupo compiten entre ellos en un torneo del que resultará ganador el que tenga la mejor aptitud. Sólo se elige a un individuo de cada subgrupo para la reproducción. Generalmente, suele inferirse un torneo binario. ¾ Selección por rango: a cada individuo de la población se le asigna un rango numérico basado en su aptitud y la selección se basa en este ranking, en lugar de las diferencias absolutas en aptitud. La ventaja de este método es que puede evitar que individuos muy aptos ganen dominancia al principio a expensas de los menos aptos, lo que reduciría la diversidad genética de la población y podría obstaculizar la búsqueda de una solución aceptable. ¾ Selección generacional: la descendencia de los individuos seleccionados en cada generación se convierte en toda la siguiente generación. No se conservan individuos entre generaciones. ¾ Selección por estado estacionario: la descendencia de los individuos seleccionados en cada generación vuelven al conjunto genético preexistente, reemplazando a algunos de los miembros menos aptos de la siguiente generación. Se conservan algunos individuos entre generaciones. ¾ Selección jerárquica: los individuos atraviesan múltiples rondas de selección en cada generación. Las evaluaciones de los primeros niveles son más rápidas y menos discriminatorias mientras que los que sobreviven hasta niveles más altos son evaluados más rigurosamente. La ventaja de este método es que reduce el tiempo total de cálculo al utilizar una evaluación más rápida y menos selectiva para eliminar a la mayoría de los individuos que se muestran poco o nada X.-15 capacitados, y sometiendo a una evaluación de aptitud más rigurosa y computacionalmente más costosa sólo a los que sobreviven a esta prueba inicial. Una dificultad en el comportamiento del algoritmo genético puede ser la existencia de gran cantidad de óptimos locales, así como el hecho de que el óptimo global se encuentre muy aislado. Un problema habitual en las ejecuciones de los algoritmos genéticos surge debido a la velocidad con la que el algoritmo converge. En algunos casos la convergencia es muy rápida lo que suele denominarse convergencia prematura, en la cual el algoritmo converge hacia óptimos locales mientras que en otros casos el problema es justo el contrario; es decir, se produce una convergencia lenta del algoritmo. El problema de la convergencia prematura surge a menudo cuando la selección de individuos se realiza de manera proporcional a su función objetivo. En tal caso, pueden existir individuos con una adaptación al problema muy superior al resto, que a medida que el algoritmo avanza dominan a la población y acaban con la diversidad genética del problema. Se debe efectuar algún tipo de artificio para que dichos superindividuos no lleguen a dominar a la población o lo que es lo mismo salir del óptimo local para encontrar otros óptimos que faciliten la localización del óptimo global. La función de selección de padres más utilizada es la denominada función de selección proporcional a la función objetivo, en la cual cada individuo tiene una probabilidad de ser seleccionado como padre proporcional al valor de su función objetivo. Esto coincide generalmente con la suposición de que si en el conjunto generado, la aptitud de cada individuo viene dada por la función objetivo, los individuos con un valor mayor de la función objetivo tienen más probabilidad de ser seleccionados para el cruce. Denominando E a la función objetivo y λ el número de individuos que constituyen la población, la probabilidad de que el individuo j sea seleccionado como padre es: pj = Ej λ ∑E j =1 X.-16 j Esta función de selección es invariante ante un cambio de escala, pero no ante una traslación. Una de las maneras de superar el problema relacionado con la rápida convergencia proveniente de los superindividuos, que surge al aplicar la anterior función de selección, es efectuar la selección proporcional al rango del individuo, con lo cual se produce una repartición más uniforme de la probabilidad de selección. Los individuos de la población han sido ordenados de menor a mayor; es decir, al peor individuo se le asigna el rango uno mientras que el individuo con mejor función objetivo tiene rango el número de individuos de la población, y la probabilidad de que un individuo sea seleccionado como padre cuando la selección se efectúa proporcional al rango del individuo, se tiene como pj = rango(E j ) λ (λ + 1) 2 La suma de los rangos, λ (λ + 1) 2 constituye la constante de normalización. La función de selección basada en el rango es invariante frente a translaciones y a cambios de escala. Cuando las probabilidades de selección varían de generación a generación (por ejemplo, la selección proporcional a la función objetivo), es lo que se denominan métodos de selección dinámicos frente a métodos de selección estáticos, en los cuales dichas probabilidades permanecen constantes (por ejemplo, la selección basada en rangos). Otro posible refinamiento es el modelo de selección del valor esperado, el cual ) ) introduce un contador cj inicializado en E j E (t ) siendo E (t ) la media de la función objetivo en la generación t. Cada vez que el individuo resulta seleccionado para el cruce, dicho contador decrece en una cantidad Δc , con Δc ∈ [0.1,0.5] . El individuo en cuestión dejará de poder ser seleccionado cuando su contador sea negativo. La dificultad de este método es la elección del valor del contador que además puede ser diferente de unos individuos a otros. X.-17 Si se asegura que todos los individuos tienen asignada una probabilidad de selección distinta de cero, el método de selección se denomina preservativo. En caso contrario se denomina extintivo. X.-18 6. Mutación Una vez que la selección ha elegido a los individuos aptos, éstos deben ser alterados aleatoriamente con la esperanza de mejorar su aptitud para la siguiente generación. Existen dos estrategias básicas para llevar esto a cabo. La primera y más sencilla se llama mutación. En la naturaleza, una mutación es un suceso bastante poco común (sucede aproximadamente una de cada mil replicaciones). En la mayoría de los casos las mutaciones son letales, pero en promedio, contribuyen a la diversidad genética de la especie. Al igual que una mutación en los seres vivos, se realiza un cambio a uno de los genes de un cromosoma elegido aleatoriamente. Cuando se usa una representación binaria, un bit se sustituye por su complemento (un cero cambia en uno y viceversa). Este operador permite la introducción de nuevo material cromosómico en la población tal y como sucede en situaciones biológicas. En un algoritmo genético, una mutación causa pequeñas alteraciones en puntos concretos del código de un individuo. El operador mutación se aplica a cada descendiente de manera individual y consiste en la alteración aleatoria (normalmente con probabilidad pequeña) de cada gen componente del cromosoma. La mutación es la causante de que los cromosomas de los hijos puedan ser diferentes a los de los padres. No conviene abusar de la mutación ya que aunque es un mecanismo generador de diversidad, y por tanto proporciona una salida cuando un algoritmo genético está estancado, reduce el algoritmo genético a una búsqueda aleatoria. Siempre es más conveniente usar otros mecanismos de generación de diversidad como aumentar el tamaño de la población o garantizar la aleatoriedad de la población inicial. fenotipo genotipo gen mutado Descendiente ↓ 1 0 1 0 0 1 0 0 164 Descendiente mutado 1 0 1 0 1 1 0 0 172 Como ocurre en el cruce, la mutación se maneja como un porcentaje que indica con qué frecuencia se producirá, aunque a diferencia de la primera sucede de forma más X.-19 esporádica (el porcentaje de cruce es normalmente de más del 60 % mientras que el de mutación nunca suele superar el 5 %). Se puede considerar diferentes técnicas de mutación dependiendo de si la mutación afecta a un bit, un grupo de bits, a un gen o a un grupo de genes: ¾ Mutación de bit: existe una única probabilidad de que se produzca una mutación de algún bit. De producirse, el algoritmo toma aleatoriamente un bit y lo invierte. ¾ Mutación multibit: cada bit tiene una probabilidad de mutarse que es calculada por el operador de mutación multibit. ¾ Mutación de gen: igual que la mutación de bit, solo que en vez de cambiar un bit, cambia un gen completo. Puede sumar un valor aleatorio, un valor constante o introducir un gen aleatorio nuevo. ¾ Mutación multigen: igual que la mutación de multibit, solo que en vez de cambiar un conjunto de bits, cambia un conjunto de genes. Puede sumar un valor aleatorio, un valor constante, o introducir un gen aleatorio nuevo. ¾ Mutación de intercambio: existe una probabilidad de que se produzca una mutación. De producirse, toma dos bits/genes aleatoriamente y los intercambia. ¾ Mutación de barajado: existe una probabilidad de que se produzca una mutación. De producirse, toma dos bits o dos genes aleatoriamente y baraja de forma aleatoria los bits o genes. X.-20 7. Cruce El segundo método de alteración aleatoria se llama cruce e implica elegir a dos individuos para que intercambien segmentos de su código, produciendo una descendencia cuyos individuos son combinación de sus padres. Este proceso simula el proceso de recombinación que se da en los cromosomas durante la reproducción sexual. Las formas comunes de cruzamiento incluyen el cruzamiento de un punto, en el que se establece un punto de intercambio en un lugar aleatorio del genoma de los dos individuos y uno de los individuos contribuye con todo su código anterior a ese punto y el otro individuo contribuye con todo su código a partir de ese punto para producir una descendencia; y el cruzamiento uniforme, en el que el valor de una posición dada en el genoma de descendencia corresponde al valor en esa posición del genoma de uno de los padres o al valor en esa posición del genoma del otro padre, elegidos todos ellos con una cierta probabilidad. Si se emparejan dos descendientes de los mismos padres; ello garantiza la perpetuación de un individuo con buenas características (es una práctica utilizada en la cría de ganado y destinada a potenciar ciertas características frente a otras). Sin embargo, si esto sucede demasiado a menudo puede crear problemas, toda la población puede aparecer dominada por los descendientes de algún gen que además puede tener caracteres no deseados. Esto se corresponde con un estancamiento en un mínimo local. Considerando el operador de cruce basado en un punto, se toman dos padres seleccionados y se cortan sus ristras de cromosomas en una posición escogida al azar para producir dos subristras iniciales y dos subristras finales. Después se intercambian las subristras finales produciendo dos nuevos cromosomas completos. Ambos descendientes heredan genes de cada uno de los padres. fenotipo genotipo punto de cruce↓ Progenitores fenotipo genotipo punto de cruce↓ 1 1 1 1 0 0 0 1 241 0 0 1 0 1 1 1 1 47 Descendientes 1 1 1 1 1 1 1 1 255 0 0 1 0 0 0 0 1 33 X.-21 Considerando el ejemplo anterior visto desde el punto de vista del genotipo, se tiene un descendiente dados dos progenitores: 255 = Rnd1 ⋅ 241 + Rnd 2 ⋅ 47 siendo Rnd1 y Rnd2 números generados aleatoriamente comprendidos entre 0 y 1. Es evidente que es posible encontrar dos valores Rnd1 y Rnd2 que cumplan la relación anterior. Mas restrictivo sería 255 = Rnd ⋅ 241 + (1 − Rnd ) ⋅ 47 Se han investigado otros operadores de cruce, habitualmente teniendo en cuenta más de un punto de cruce. De Jong investigó el comportamiento del operador de cruce basado en múltiples puntos concluyendo que el cruce basado en dos puntos representaba una mejora, mientras que añadir más puntos de cruce no beneficiaba el comportamiento del algoritmo. La ventaja de tener más de un punto de cruce radica en que el espacio de búsqueda puede ser explorado de manera más fácil siendo la principal desventaja el hecho de aumentar la probabilidad de romper buenas estructuras. fenotipo punto de cruce↓ gen. ↓ fenotipo punto de cruce↓ gen. ↓ Progenit. 1 0 1 0 0 0 1 1 1 0 654 0 0 1 1 0 1 0 0 1 0 210 Descend. 1 0 1 0 0 1 0 1 1 0 662 0 0 1 1 0 0 1 0 1 0 202 Tomando el punto de vista del genotipo; se tiene un descendiente dados dos progenitores: 662 = Rnd1 ⋅ 654 + Rnd 2 ⋅ 210 En el operador de cruce uniforme de Syswerda, cada gen, en la descendencia se crea copiando el correspondiente gen de uno de los dos padres, escogido de acuerdo a una máscara de cruce generada aleatoriamente. Cuando existe un 1 en la máscara de cruce, X.-22 el gen es copiado del primer padre, mientras que cuando exista un 0, el gen se copia del segundo padre. Máscara de cruce 1 0 0 1 0 0 1 Progenitor 1 1 1 0 ↓ 1 0 0 ↑ ↑ 0 0 0 Descendiente Progenitor 2 1 1 0 ↓ 1 1 1 ↑ ↑ 1 1 1 1 ↓ 1 0 Existen diversas variantes a la idea de máscara de cruce al introducir el concepto de probabilidad de herencia. Si esta probabilidad es alta se heredan mayores características de los padres. Esta probabilidad puede depender de una gran variedad de parámetros. Aunque podría pensarse que el operador cruce es más importante que el operador mutación ya que proporciona una exploración rápida del espacio de búsqueda, éste último se considera un operador básico que proporciona un pequeño elemento de aleatoriedad en la vecindad de los individuos de la población. Si bien se admite que el operador de cruce es el responsable de efectuar la búsqueda a lo largo del espacio de posibles soluciones, varios investigadores piensan que el operador mutación va ganando en importancia a medida que la población va convergiendo. Utilizando la denominada evolución primitiva, en la cual el proceso evolutivo consta tan sólo de selección y mutación, Schaffer encontró que dicha evolución primitiva supera con creces a una evolución basada exclusivamente en selección y cruce. Es importante introducir la definición de convergencia introducida por De Jong. Si el algoritmo genético ha sido correctamente implementado, la población evolucionará a lo largo de generaciones sucesivas de tal manera que la adaptación media de todos los individuos de la población, así como la adaptación del mejor individuo se irán incrementando hacia el óptimo global. El concepto de convergencia está relacionado con la progresión hacia la uniformidad, un gen ha convergido cuando un elevado porcentaje de los individuos de la población comparten el mismo valor para dicho gen. Se dice que la población converge cuando todos los genes han convergido. X.-23 Normalmente el cruce se emplea dentro de la implementación del algoritmo genético como un porcentaje que indica con qué frecuencia se efectuará. Esto significa que no todas las parejas de cromosomas se cruzarán sino que habrá algunas que pasen intactas a la siguiente generación. De hecho existe una técnica, denominada elitismo, en la que el individuo más apto a lo largo de las distintas generaciones no se cruza con nadie y se mantiene intacto hasta que surge otro individuo mejor que él, que lo desplazará. X.-24 8. Reducción Una vez obtenidos los individuos descendientes de una determinada población, el proceso de reducción al tamaño original consiste en escoger la nueva generación de población entre los progenitores y los descendientes de los mismos. Existen dos planteamientos fundamentales a este respecto, en el primero los individuos generados sustituyen a sus progenitores con lo que un individuo no convive nunca con sus antecesores (reducción simple). En el segundo planteamiento se sustituyen los individuos peor adaptados de toda la población (por eso se dice que este modelo sigue una política elitista), con lo que existirá convivencia entre un individuo y sus ancestros. No obstante, se han planteado otras consideraciones sustitutivas a fin de evitar la aparición de elitismo como sustituir una tasa de la población o incluso otros. El concepto de reducción está ligado con el de tasa de reemplazamiento generacional; es decir, el porcentaje de descendientes generados con respecto del tamaño de la población. Aunque habitualmente suele considerarse un tamaño fijo de población, puede establecerse un aumento de población o una reducción como sucede en biología pudiendo llegar a extinguirse la población en cuestión. X.-25 9. Población Una interesante pregunta es cual es el tamaño idóneo de la población. Las poblaciones pequeñas corren el riesgo de no cubrir adecuadamente el espacio de búsqueda mientras que el trabajar con poblaciones de gran tamaño puede acarrear problemas relacionados con excesivo costo computacional. Goldberg efectuó un estudio teórico obteniendo como conclusión que el tamaño óptimo de la población para ristras de longitud dada, con codificación binaria, crece exponencialmente con el tamaño de la ristra. Este resultado traería como consecuencia un tamaño enorme de la población con lo que la aplicabilidad de los algoritmos genéticos en problemas reales sería muy limitada, ya que resultarían poco competitivos con otros métodos de optimización combinatoria. Alander, basándose en evidencias empíricas sugirió que un tamaño de población no superior a 21 individuos sería suficiente para atacar con éxito los problemas. Habitualmente la población inicial se escoge generando datos al azar, pudiendo contener cada gen uno de los posibles valores del alfabeto, con probabilidad uniforme. En los pocos análisis acerca de que sucedería si los individuos de la población inicial se obtuviesen como resultado de alguna técnica heurística o de optimización local, se ha constatado que esta inicialización no aleatoria de la población inicial puede acelerar la convergencia del algoritmo genético, aunque en algunos casos la desventaja resulta ser la prematura convergencia del algoritmo o lo que es lo mismo la convergencia hacia óptimos locales. Generalmente, salvo poblaciones iniciales realmente degeneradas en los que el operador de mutación va a tener mucho trabajo, la convergencia es poco sensible a la población inicial si esta se escoge de forma aleatoria y es lo suficientemente grande. X.-26 10. Condición de parada Básicamente, la condición de terminación debe ser la convergencia del algoritmo genético hacia una solución estable o alcanzar un número prefijado de generaciones. Como el algoritmo genético opera con una población en cada iteración, se espera que el método converja de manera que al final del proceso la población sea muy similar y en el infinito se reduzca a un sólo individuo. Para cadenas binarias existe una compleja teoría para demostrar la convergencia de los algoritmos genéticos. Básicamente, esta teoría se basa en considerar que una cadena es un representante de una clase de equivalencia o esquema focalizando la búsqueda entre esquemas en vez de entre cadenas. A partir de este resultado el teorema de esquemas demuestra que la población converge a unos esquemas que cada vez son más parecidos y en el límite teórico a una única cadena. En el caso de cadenas no binarias, el problema está menos estudiado y aparecen los conceptos de forma y conjunto de similitud que intentan generalizar al de esquema. Exigiendo que al cruzar dos cadenas de la misma clase se obtenga otra de esta misma clase, y asumiendo la selección de los progenitores bajo ciertas condiciones, es posible demostrar la convergencia del algoritmo. La realidad es que las condiciones impuestas suelen ser muy estrictas y no se suelen cumplir, por lo que en ocasiones los algoritmos genéticos resuelven satisfactoriamente un problema y otras se quedan más lejos del óptimo. La curva de convergencia asociada al algoritmo suele presentar una convergencia rápida al principio, para casi enseguida bloquearse. Esto se debe a que el algoritmo genético es excelente descartando subespacios realmente malos. Cada cierto tiempo, la población vuelve dar el salto evolutivo y se produce un incremento mayor en la velocidad de convergencia. La razón de esto es que algunas veces aparece una mutación altamente beneficiosa, o un individuo excepcional, que propaga algún conjunto de cromosomas excepcional al resto de la población. X.-27 Se van a establecer tres condiciones de parada que detienen el algoritmo si se cumple una cualquiera de ellas. 10.1. Condición 1 Si la solución no ha mejorado, tras K1 series consecutivas de L pasos, una cierta cantidad definida por ∈1 > 0 . E (K1 ⋅ L + n ) > E (n ) ⋅ (1+ ∈1 ) ⇒ FIN 10.2. Condición 2 La convergencia de la población se espera en la medida que al final del proceso la población sea muy similar, y en el infinito se reduzca a un sólo individuo. Con esta condición, se comprobará el grado de semejanza de la población para lo cual se va a emplear un análisis estadístico. Cada individuo se codifica como un conjunto de cromosomas. El conjunto de parámetros representando un cromosoma particular se denomina fenotipo y a partir del fenotipo se obtiene el genotipo. Finalmente, un individuo concreto se representa como un conjunto de genotipos. El análisis que se va a realizar a continuación se reseña para un genotipo concreto. La condición de parada deberá cumplirse para cada uno de los genotipos que representa un individuo. Considérese el genotipo a analizar de una población concreta compuesta por N individuos e imagínese que esta población concreta es una muestra extraída aleatoriamente de una población madre compuesta por una colección de genotipos mucho mayor. Con esta hipótesis de partida, es posible aplicar la Teoría de la Estimación. Considerando el estadístico media y el estadístico varianza: a= x1 + K + x N N X.-28 s 2 2 2 ( x1 − a ) + K + ( x N − a ) = N Se considera que se ha alcanzado la condición de parada cuando un alto porcentaje de la población se encuentra dentro del intervalo de confianza elegido. Se ha tomado un intervalo de confianza del 90 %, aunque el análisis puede efectuarse para cualquier otro intervalo de confianza. De este modo, se considera que la población converge cuando el 90 % de los individuos se encuentran dentro del intervalo de confianza de la distribución estadística. Asumiendo que la población madre presenta una distribución normal, se puede estimar la media y la desviación típica. Para una gran muestra (N > 30), la media se estima como: α = a ± w90 [N (0,1)] s N siendo w90 [N (0,1)] el valor que define el intervalo de confianza de la distribución normal N(0,1), de tal forma que el 90 % de los individuos se encuentran dentro del intervalo de confianza de la campana de Gauss. Para una pequeña muestra (N < 30): α = a ± w90 (t N −1 ) s N −1 siendo w90(tN-1) el valor que define el intervalo de confianza de la distribución t de Student con N-1 grados de libertad, de tal forma que el 90 % de los individuos se encuentran dentro del intervalo de confianza de la distribución. La desviación típica se estima de la siguiente expresión: χ N2 −1 = N ⋅ s2 σ2 X.-29 siendo χ N2 −1 la distribución chi-cuadrado con N-1 grados de libertad. N ⋅ s2 σ2 = ( ) ( ) w0.05 χ N2 −1 N ⋅ s2 w0.95 χ N2 −1 Los dos valores de σ 2 definen el intervalo de confianza de la distribución chi-cuadrado con N-1 grados de libertad, de tal forma que el 90 % de los individuos se encuentran dentro del intervalo de confianza de la distribución. Esta condición se comprobará tras K2 series consecutivas de L pasos. 10.3. Condición 3 Si el número de iteraciones alcanza un valor superior al valor máximo establecido M. X.-30 11. Implementación del algoritmo Los algoritmos genéticos son métodos de optimización que tratan de hallar (xi,...,xn) tales que E(xi,...,xn) sea máximo. En un algoritmo genético, tras parametrizar el problema en una serie de variables, (xi,...,xn), se codifican en forma de cromosomas. Todos los operadores utilizados por un algoritmo genético se aplicarán sobre estos cromosomas o sobre poblaciones de ellos. Un algoritmo genético es independiente del problema a resolver lo cual lo hace un algoritmo robusto, por ser útil para cualquier problema, pero a la vez débil pues no está especializado en ninguno. Se elige un punto inicial del espacio de soluciones que cumpla las condiciones del problema. En el sistema se debe definir el tamaño de paso que permite definir el resto de individuos de la población inicial. Este tamaño debe ser un orden de magnitud menor que el tamaño característico de óptimo local y está íntimamente ligado a la suavidad de la función objetivo. Un grado de discretización excesivamente fino penaliza el rendimiento mientras que un grado de discretización muy grande va ha hacer que los pasos sean inútiles, ya que no se consideran soluciones, limitándose el algoritmo a funcionar como un mero explorador de soluciones. Este hecho condiciona el espacio de soluciones que será capaz de alcanzar el método. En todo problema existe una relación entre el espacio de soluciones, el tamaño de paso, el número de iteraciones a realizar y la velocidad de cálculo del computador. Para un problema concreto, la velocidad de cálculo del computador condiciona el número de iteraciones a realizar, por lo que definido el tamaño del paso se tiene el espacio de soluciones. Evidentemente, para explorar un determinado espacio de soluciones se debe tomar la solución inicial en ese espacio. Todo esto si el espacio de soluciones es lo suficientemente extenso, pero puede ocurrir que el espacio de soluciones se encuentre limitado lo que condiciona los individuos de la población. En el problema que se está analizando, la obtención de la función objetivo requiere un elevado tiempo de cálculo de ordenador ya que implica el cálculo de una trayectoria balística completa. Este inconveniente condiciona el concepto de tasa de reemplazamiento generacional; es decir, la reducción de la población. De este modo, al X.-31 efecto de reducir el tiempo de ordenador sustancialmente, solo se va a considerar un único descendiente de dos progenitores. En cada generación, se localiza el individuo menos apto que desaparecerá siendo sustituido por un descendiente de dos padres. Así, la población estará formada por todos los ancestros y el único descendiente que se ha generado. La gran cantidad de opciones que se pueden elegir al utilizar un algoritmo genético hace necesario establecer algún criterio que permita identificar que características son más apropiadas para el problema concreto a resolver. La forma de elegir la población inicial, el número de individuos de los que consta la población, el tipo de selección, etc. son cuestiones difusas a priori. Para estudiar este fenómeno se va a calcular un problema tipo con diferentes configuraciones y se van a comparar las soluciones obtenidas de forman análoga a como ya se realizó al aplicar el método de recocido simulado Así, se van a localizar los factores de ajuste de un problema balístico para un proyectil base burn, empleando el método 1 del Modelo Modificado de Masa Puntual. Utilizando el algoritmo de recocido simulado se van a determinar los cuatro factores de ajuste tDI, K(p) (función lineal) y f(iBB,MT) para un ángulo de tiro concreto. Se ha escogido para ello el proyectil OE155F2RTC FR HEA BB1, disparado desde un tubo 155AUF1 con una carga C7, por ser un proyectil muy estudiado del que se dispone abundante información calculada mediante el método 1. Dado que en este algoritmo la velocidad de cálculo es muy importante para obtener resultados de forma rápida, se ha recurrido al método 1 por permitir una mayor velocidad en la computación. Inicialmente, se realizaron una serie de cálculos empleando un ángulo de tiro de 300 milésimas pero se ha podido comprobar que la dispersión de resultados aumenta en la medida que se emplean ángulos de tiro superiores. Posteriormente, se podrá ver que para mejorar el proceso de optimización del algoritmo, interesa tener una mayor dispersión en los resultados que permitan ver claramente las diferencias de una configuración a otra. De igual forma, se ha escogido la carga más elevada por presentar una mayor dispersión en los datos obtenidos. Aumentando el ángulo de tiro, se han encontrando muy pocas diferencias entre 700 y 900 milésimas por lo que finalmente los cálculos se han efectuado con un ángulo de tiro de 700 milésimas, al objeto de evitar efectos asociados al segundo sector que podrían enmascarar y confundir los resultados obtenidos, y también para reducir el tiempo de cálculo. X.-32 Para cada carga, de los cuatro factores de ajuste, únicamente el factor base burn f(iBB,MT) es función del ángulo de tiro siendo el tiempo de retraso a la ignición del motor base burn tDI función de la temperatura de propulsante del motor MT, y el factor de quemado con la velocidad de rotación axial del proyectil K(p) una función lineal de la velocidad de rotación axial del proyectil p. En este momento se va a trabajar con un único ángulo de tiro y no es de aplicación pero se debe tener en cuenta que mientras que f(iBB,MT) es función del ángulo de tiro, tDI y K(p) no lo son. Obtenido el alcance y la deriva de experiencias realizadas, se van a calcular los factores de ajuste tDI, K(p) y f(iBB,MT) que mejor aproximen el alcance y la deriva. Se han escogido estos parámetros por ser los más representativos del tiro de artillería. La adición a los anteriores de la ordenada máxima ha permitido constatar un empeoramiento en el cálculo de los factores de ajuste. Si se obtienen unos factores de ajuste que mejor ajusten alcance y deriva y unos segundos factores de ajuste que mejor ajusten alcance, deriva y ordenada máxima; se ha observado que si se calcula alcance y deriva con los segundos factores de ajuste, se obtienen peores valores que usando los primeros factores de ajuste. Esto es debido a que la ordenada máxima es un parámetro que ante variaciones en los factores de ajuste, experimenta variaciones mucho mayores que las apreciadas en alcance y deriva. En concreto, las variaciones en la ordenada máxima pueden llegar a ser del orden de diez veces superiores a las encontradas para alcance y deriva. También se ha constatado la dificultad por parte del Modelo Modificado de Masa Puntual de proporcionar valores acertados de la ordenada máxima, una vez ajustada la deriva y el alcance, con datos extraídos de experiencias reales. Finalmente, se ha optado por ajustar únicamente alcance y deriva. Si a través de experiencias realizadas se obtiene el alcance y la deriva, que se designaran como alcanceexp y derivaexp, se define la función objetivo como: E= (alcance − alcance ) + (deriva − deriva ) 2 exp 2 exp siendo los términos alcance y deriva los obtenidos como resultado de haber realizado todo el proceso de cálculo en función de los factores de ajuste tDI, K(p) y f(iBB,MT). Una X.-33 característica muy importante de este problema es que se conoce, a priori, el valor óptimo buscado de la función objetivo. Como condiciones iniciales, se considera que para el tiempo de retraso a la ignición del motor base burn tDI se introducirá el valor correspondiente, en el supuesto de que se conozca. Si no se dispone de información adicional, se considerará tDI = 0 siendo el resto de valores K(p) = 1 y f(iBB,MT) = 1 . Como ya se comentó anteriormente, el algoritmo de integración utilizado ha sido el Runge Kutta 45. Los parámetros eRK, hmin y hmax se han elegido para obtener una precisión inferior a un metro en el cálculo de la deriva y el alcance: hmin = 0.04 s, hmax = 6.0 s y eRK = 0.1 . Inicialmente, se va a efectuar una selección de padres al azar favoreciendo a los individuos mejor adaptados, asignando a cada individuo una probabilidad de ser seleccionado proporcional a la función objetivo. Este procedimiento está basado en la ruleta sesgada. Según dicho esquema, los individuos bien adaptados se escogerán probablemente varias veces por generación mientras que los pobremente adaptados al problema no se escogerán más que de vez en cuando. Esta función de selección de padres es más utilizada y se denomina función de selección proporcional a la función objetivo. Siendo E la función objetivo y λ el número de individuos que constituyen la población, la probabilidad de que el individuo j sea seleccionado como padre es: pj = Ej ; 0 ≤ pj <1 λ ∑E j =1 j Definiendo para cada individuo el parámetro qj como i qi = ∑ p j j =1 Generando un número aleatorio Rnd, comprendido entre 0 y 1; si qk −1 < Rnd ≤ qk el X.-34 individuo seleccionado es el k siendo q0 = 0 . El segundo progenitor se elige de igual modo aunque se debe tener en cuenta que si un Ej fuese muy grande comparado con el resto de individuos, este individuo sería seleccionado siempre y sería imposible elegir un segundo progenitor, por lo que en este caso el segundo progenitor se elige de forma aleatoria pura. Una vez elegidos los progenitores se va a proceder a combinar el material genético para obtener un descendiente. Considerando la dificultad del problema balístico elegido para aplicar el algoritmo genético, se va a emplear una codificación mediante números decimales. Este enfoque permite una mayor precisión y complejidad que el método de utilizar sólo números binarios, y está intuitivamente más cerca del espacio de problemas. Recordando la expresión de cruce analizada anteriormente descendiente = Rnd1 ⋅ progenitor1 + Rnd 2 ⋅ progenitor 2 siendo Rnd1 y Rnd2 números generados aleatoriamente comprendidos entre 0 y 1. Una vez seleccionada la población inicial, el operador de cruce es el encargado de efectuar la búsqueda a lo largo del espacio de posibles soluciones, pero el operador mutación es el responsable de introducir nuevo material genético en la población tal y como sucede en biología. De ahí, la gran importancia del operador mutación. En los algoritmos genéticos, una mutación es un suceso muy poco habitual aunque su frecuencia es mayor que en la naturaleza. Se recomienda que dicho porcentaje no sea superior al 5 %. En biología, la mayoría de las mutaciones no sobreviven aunque el pequeño porcentaje de las que si lo hacen introducen información genética adicional que permite mejorar la especie. En los algoritmos genéticos, la situación es análoga ya que si la mutación representa una mejora para alcanzar el óptimo, el individuo pasa a formar parte de la población pero si no lo hace, el individuo termina por desaparecer. Este razonamiento permite justificar la utilización del operador mutación en un porcentaje no superior al máximo definido, que es el que se ha empleado finalmente. Se ha considerado apropiada para este problema la siguiente condición de parada; L = X.-35 100, K1 = 4, ∈1 = 0.5 , K2 = 4, ∈2 = 0.5 y M = 1000. En el peor de los casos, cada proceso de cálculo independiente finaliza cuando se alcanza el número máximo de iteraciones M, que teniendo en cuenta el tiempo de ejecución de ordenador, se ha considerado suficiente fijarlo en 1000. A fin de comprobar la calidad de los resultados obtenidos se va a establecer un sistema en el cual, tras alcanzar la condición de parada y extraer el óptimo correspondiente, se va a continuar calculando hasta volver a alcanzar una nueva condición de parada con su correspondiente nuevo óptimo, y así sucesivamente. Así, se han ejecutado veinte de estos procesos de cálculo diferentes, comprobando como cuando el algoritmo alcanza el óptimo puede hacerlo con diferentes valores de factores de ajuste; es decir, aunque el valor del óptimo sea muy similar, existen diferentes combinaciones de los factores de ajuste que alcanzan ese valor del óptimo. A fin de investigar este fenómeno, se ha considerado realizar otros diez procesos de cálculo fijando tDI = 0 , no obteniendo grandes diferencias con respecto al caso anterior. Finalmente, se han calculado otros diez procesos de cálculo fijando tDI = 0 y K(p) una función constante, sin considerar la dependencia con la velocidad de rotación axial del proyectil p, mostrando estos últimos resultados, como en los diez procesos de cálculo se ha alcanzado aproximadamente el mismo óptimo con similares factores de ajuste. Resumiendo, se ha establecido un número de procesos cálculo de cuarenta, los cuales se utilizaran posteriormente a efectos de valorar los resultados proporcionados por el algoritmo. A partir del valor inicial es necesario definir el conjunto de individuos que componen la población inicial. Como cada individuo se encuentra representado por un conjunto de genotipos, es necesario producir variaciones en estos genotipos que permitan obtener la población. A tenor de esta idea, se va a definir un tamaño de paso que con carácter general permitirá extender la población: Δ = f esc ⋅ Pre ⋅ Rnd ⋅10 siendo fesc el factor de escala, Pre la precisión y Rnd un número aleatorio. Este tamaño de paso debe ser un orden de magnitud menor que el tamaño característico de óptimo local. Un tamaño mayor no permitirá al algoritmo funcionar como algoritmo X.-36 de optimización convirtiéndose, en un mero explorador de soluciones. Esta es la idea inherente al introducir el parámetro precisión. En este problema se ha encontrado una buena aproximación para Pre = 10e-4 . Este tamaño puede ser diferente dependiendo de la variable a la que vaya a afectar ya que la escala de cada una de las dimensiones no tiene porque ser la misma. Esto conlleva a que es necesario definir el factor de escala que afecta a cada una de las variables. A partir del valor inicial es necesario especificar el conjunto de individuos que componen la población inicial. Precisado el tamaño de paso, se va a proceder a definir cada uno de los individuos que componen la población inicial para lo cual se va a emplear un método puramente aleatorio consistente en tomar el vector de estado o lo que es lo mismo, el conjunto de genotipos que en este caso son cuatro genotipos; tDI, K(p) y f(iBB,MT), e introducir una pequeña perturbación a partir de las condiciones iniciales sumando una cantidad dada Δ a cada genotipo, con Rnd comprendido entre -1 y 1 . Esto se efectuara para cada uno de los individuos que componen la población inicial. En cuanto al número de individuos que componen la población, haciendo referencia a consideraciones estadísticas que establecen como infinito estadístico el valor de treinta, se ha considerado inicialmente una población constituida por treinta individuos. Una vez definido el tamaño de paso, es posible definir el concepto de mutación de acuerdo a la forma de codificación elegida mediante números decimales. Δ m = f esc ⋅ Pre ⋅ Rnd ⋅10 n siendo fesc el factor de escala, Pre la precisión, Rnd un número aleatorio comprendido entre -1 y 1, y n un número aleatorio comprendido entre 1 y 3. Los elementos fesc y Pre toman los mismos valores que los introducidos en Δ .Como ya se indico anteriormente la frecuencia con la que aparece la mutación no debe superar el valor reseñado. De este modo, cuando surge el individuo candidato a sufrir una mutación, se elige al azar uno de X.-37 los genotipos que compone dicho individuo y se le adiciona una cantidad dada Δ m . La gran cantidad de condicionantes que tiene este algoritmo permite imaginar tantas posibilidades que resulta materialmente impensable mostrarlas todas. No obstante, una primera pregunta que surge y que ya se comentó con anterioridad es que incidencia puede tener la utilización del método de cruce. Si se comparan las expresiones escritas a continuación, se puede observar a priori que la segunda es más restrictiva que la primera; no obstante, para el problema a resolver cabe preguntarse cual proporcionará mejores resultados. descendiente = Rnd1 ⋅ progenitor1 + Rnd 2 ⋅ progenitor 2 descendiente = Rnd ⋅ progenitor1 + (1 − Rnd ) ⋅ progenitor 2 En las siguientes gráficas se van a comparar las dos configuraciones descritas anteriormente mostrando alcance y deriva. Se ha designado como GA2 aquella que introduce dos valores Rnd y GA1 aquella que introduce un único valor Rnd. Definiendo el área existente entre el valor calculado y el valor exacto para GA2 y GA1, y calculando la diferencia entre las dos áreas, puede estimarse que valor calculado se acerca más al valor exacto. Así, se ha obtenido que en alcance GA2 resulta mucho mejor que GA1 ya que un 50 % de los puntos que se obtienen empleando GA1 presentan un error bastante elevado. En cuanto a la deriva se tiene que GA1 presenta una ligera mejoría frente a GA2; sin embargo, porcentualmente hablando, esta mejoría es muy pequeña y no compensa el empeoramiento que se obtiene en alcance. Esto condiciona la elección de GA2, empleando dos valores de Rnd diferentes, como método más eficaz de cálculo debido a que aunque en deriva GA1 presenta ligeras mejoras, al comparar el alcance se puede observar que GA1 muestra errores muy elevados. X.-38 Otra opción sería variar el número de individuos que componen la población inicial para X.-39 conocer si existe alguna mejora por el hecho de emplear una población de un número de individuos determinado. Anteriormente, se estableció una población constituida por treinta individuos. Una población inicial mayor introducirá un mayor contenido genético lo cual constituye una ventaja a la hora de intentar salir de mínimos locales, por lo que se va a comparar el tamaño de población original de treinta individuos con una población formada por cien individuos. Este valor viene condicionado por el número máximo de iteraciones a realizar establecido en las condiciones de parada ya que en este método la población debe evolucionar hacia la convergencia y se debe permitir que las generaciones se crucen y tengan descendencia. De este modo y teniendo en cuenta el número máximo de iteraciones establecido anteriormente, se ha considerado como apropiado poder llegar a tener una población del orden del 10 % el número máximo de iteraciones. Analizando los resultados obtenidos, puede observarse como en general se obtiene una ligera mejoría al utilizar una población mayor. X.-40 Anteriormente se definió la población inicial tomando el vector de estado o lo que es lo mismo, el conjunto de genotipos; tDI, K(p) y f(iBB,MT), e introduciendo una pequeña perturbación a partir de las condiciones iniciales sumando una cantidad dada Δ a cada genotipo, con Rnd comprendido entre -1 y 1 , para cada uno de los individuos que componen la población inicial. No obstante, se plantea definir la población inicial empleando otro método de claras influencias en las técnicas de recocido simulado, consistente en tomar el vector de estado para cada individuo y realizar un sorteo entre todas las componentes para decrementar la elegida una cantidad dada Δ , con Rnd comprendido entre 0 y 1 .A continuación se realizará un nuevo sorteo entre las componentes restantes, que no han variado, y se incrementa la elegida un nuevo valor Δ , con Rnd comprendido entre 0 y 1 . Esto se hará para cada uno de los individuos que componen la población inicial. Esta última forma de definir la población se ha designado como GA pob2 mientras que la primera se ha designado como GA pob1. Analizando los resultados obtenidos, puede observarse como en general se obtiene una X.-41 ligera mejoría con GA pob1. X.-42 Un problema usual que puede surgir es la velocidad con que el algoritmo converge pudiendo aparecer la convergencia prematura, en la cual el algoritmo converge hacia óptimos locales. El problema de la convergencia prematura, nace a menudo cuando la selección de individuos se realiza de manera proporcional a su función objetivo. En tal caso, pueden existir individuos con una adaptación al problema muy superior al resto, que a medida que avanza el algoritmo dominan a la población y acaban con la diversidad genética del problema. Se debe efectuar algún tipo de ardid para que dichos superindividuos no lleguen a dominar a la población; es decir, salir del óptimo local para encontrar otros óptimos que faciliten la localización del óptimo global. Una de las maneras de soslayar el problema relacionado con la rápida convergencia proveniente de los superindividuos, que surge al aplicar la función de selección proporcional a la función objetivo, es efectuar la selección mediante torneo. Se eligen dos subgrupos de individuos de la población al azar, y los miembros de cada subgrupo rivalizan entre ellos en un torneo del que resultará ganador el que tenga mejor aptitud. Sólo se elige un individuo de cada subgrupo para la reproducción. En las siguientes gráficas se va a comparar las dos configuraciones relatadas anteriormente. Se ha designado como GA PFO aquella que incorpora la selección proporcional a la función objetivo y GA torneo aquella que utiliza la selección como su propio nombre indica. Analizando los resultados, puede observarse como GA torneo presenta mejorías respecto a GA PFO en cuanto a un mayor agrupamiento de resultados en torno al valor exacto. X.-43 X.-44 Otra forma de superar el problema relacionado con la rápida convergencia proveniente de los superindividuos, que surge al aplicar la función de selección proporcional a la función objetivo, es efectuar la selección proporcional al rango del individuo con lo que se produce una repartición más uniforme de la probabilidad de selección. Los individuos de la población han sido ordenados de menor a mayor; es decir, al peor individuo se le asigna el rango uno, mientras que el individuo con mejor función objetivo tiene rango el número de individuos de la población, y la probabilidad de que un individuo sea seleccionado como padre cuando la selección se efectúa proporcionalmente al rango del individuo, se tiene como pj = rango(E j ) λ (λ + 1) 2 La suma de los rangos, λ (λ + 1) 2 , constituye la constante de normalización. La función de selección basada en el rango es invariante frente a translaciones y a cambios de escala. En las siguientes gráficas se va a comparar las dos configuraciones descritas anteriormente. Se ha designado como GA rango y GA torneo los que utilizan la selección como su propio nombre indica. Analizando los resultados, puede observarse como la diferencia entre ambos es muy pequeña aunque puede asegurarse sin ningún género de dudas que son superiores a la selección proporcional a la función objetivo. X.-45 X.-46 X.-47 X.-48 ESCUELA POLITÉCNICA SUPERIOR DEL EJÉRCITO TESIS DOCTORAL CAPÍTULO XI RECOCIDO SIMULADO FRENTE A ALGORITMO GENÉTICO MODELIZACIÓN DE PROYECTILES BASE BURN Director: Dr. Coronel Francisco Cucharero Pérez Doctorando: Comandante Fernando Aguirre Estévez XI.-1 XI.-2 XI.-RECOCIDO SIMULADO FRENTE A ALGORITMO GENÉTICO Se van a comparar los resultados obtenidos anteriormente empleando el algoritmo de recocido simulado frente al algoritmo genético. Para ello se va a continuar el estudio iniciado en los capítulos anteriores, analizando el problema balístico para un proyectil base burn empleando el método 1, en el cual se van a determinar los factores de ajuste tDI, K(p) (función lineal) y f(iBB,MT). Se ha utilizado el proyectil OE155F2RTC FR HEA BB1, disparado desde un tubo 155AUF1, con una carga C7 y un ángulo de tiro de 700 milésimas. Obtenido el alcance y la deriva de experiencias realizadas, se han calculado los factores de ajuste tDI, K(p) y f(iBB,MT) que mejor aproximen alcance y deriva. Supóngase que a través de experiencias realizadas se obtiene el alcance y la deriva, que se designaran como alcanceexp y derivaexp, con lo que la función objetivo se define como E= (alcance − alcance ) + (deriva − deriva ) 2 exp 2 exp siendo los términos alcance y deriva los obtenidos como resultado de haber realizado todo el proceso de cálculo en función de los factores de ajuste tDI, K(p) y f(iBB,MT). Se ha considerado apropiada para este problema la siguiente condición de parada; L = 100, K1 = 4, ∈1 = 0.5 , K2 = 4, ∈2 = 0.5 y M = 1000. En el peor de los casos, cada proceso de cálculo independiente finaliza cuando se alcanza el número máximo de iteraciones M, que teniendo en cuenta el tiempo de ejecución de ordenador, se ha considerado suficiente fijarlo en 1000. A fin de comprobar la calidad de los resultados obtenidos se va a establecer un sistema en el cual, tras alcanzar la condición de parada y extraer el óptimo correspondiente, se va a continuar calculando hasta volver a alcanzar una nueva condición de parada con su correspondiente nuevo óptimo, y así sucesivamente. Así, se han ejecutado veinte de estos procesos de cálculo diferentes, comprobando como cuando el algoritmo alcanza el óptimo puede hacerlo con diferentes valores de factores de ajuste; es decir, aunque el valor del óptimo sea muy similar, existen diferentes combinaciones de los factores de ajuste que alcanzan ese valor del óptimo. A fin de investigar este fenómeno, se ha considerado realizar otros diez procesos de cálculo XI.-3 fijando tDI = 0 , no obteniendo grandes diferencias con respecto al caso anterior. Finalmente, se han calculado otros diez procesos de cálculo fijando tDI = 0 y K(p) una función constante, sin considerar la dependencia con la velocidad de rotación axial del proyectil p, mostrando estos últimos resultados, como en los diez procesos de cálculo se ha alcanzado aproximadamente el mismo óptimo con similares factores de ajuste. Resumiendo, se ha establecido un número de procesos cálculo de cuarenta, los cuales se utilizaran posteriormente a efectos de valorar los resultados proporcionados por el algoritmo. En las siguientes gráficas se van a comparar los resultados obtenidos en los capítulos anteriores utilizando el algoritmo de recocido simulado frente al algoritmo genético. En cada uno de ellos, ya se efectuó un extenso proceso de investigación y mejora de manera independiente, para ser capaz de obtener los mejores resultados aprovechando la infraestructura de cálculo disponible. En este punto, ha llegado el momento de comparar uno contra el otro mostrando el alcance y la deriva. En el caso del algoritmo genético, se ha utilizado la configuración de torneo. Analizando los resultados obtenidos, puede observarse como el algoritmo de recocido simulado presenta una mejoría enorme frente al algoritmo genético. XI.-4 Si bien el algoritmo de recocido simulado resulta mejor que el algoritmo genético, el primero presenta el inconveniente de que requiere una labor de determinación (y por tanto de prueba y error) muy exhaustiva de ciertos parámetros, en concreto la constante de Boltzmann debe reajustarse en cuanto se cambia cualquier otro parámetro que afecte al problema. Esto implica un mayor tiempo de cálculo al tener que realizar ajustes constantes para obtener buenos resultados. Este inconveniente no lo presenta el algoritmo genético que es más robusto en ese aspecto, si bien menos preciso. XI.-5 XI.-6 ESCUELA POLITÉCNICA SUPERIOR DEL EJÉRCITO TESIS DOCTORAL CAPÍTULO XII MODELIZACIÓN DE PROYECTILES BASE BURN MÉTODO 1 MODELIZACIÓN DE PROYECTILES BASE BURN Director: Dr. Coronel Francisco Cucharero Pérez Doctorando: Comandante Fernando Aguirre Estévez XII.-1 XII.-2 XII.-MODELIZACIÓN DE PROYECTILES BASE BURN MÉTODO 1 1. Índice 1. Índice 2. Introducción 3. Proyectiles base burn método 1 3.1. Área de combustión del propulsante base burn SC 3.2. Coeficiente de reducción de resistencia base burn CxBB 3.3. Parámetro de inyección de eficiencia óptima I0 3.4. Implementación del método de cálculo XII.-3 XII.-4 2. Introducción El cálculo de trayectorias balísticas efectuado numéricamente de un proyectil base burn necesita del conocimiento de los coeficientes aerodinámicos, los coeficientes base burn y los factores de ajuste. Los coeficientes aerodinámicos pueden obtenerse a través de programas informáticos existentes, de ensayos en túneles de viento o bien de bases de datos suministradas por el fabricante del proyectil. Existen aplicaciones informáticas como PRODAS o McDrag que proporcionan los coeficientes aerodinámicos a partir de la geometría del proyectil en base a ciertos ajustes fenomenológicos obtenidos a partir de ensayos; o bien, a partir de la integración numérica de las ecuaciones de Navier Stokes, proceso conocido abreviadamente como CFD (Computational Fluid Dynamics). La modelización de una forma detallada del proceso de quemado del propulsante así como de la reducción de la resistencia al aumentar la presión en la base, se materializa en el STANAG 4355 a través de una serie de funciones que caracterizan el comportamiento del proyectil. Esto implica un conocimiento muy detallado de la reacción química que ocurre al quemarse el propulsante en la cámara de combustión, y que ocasiona una disminución de la resistencia al incrementarse la presión en la base por la inyección de los gases procedentes de la combustión del propulsante. De este forma, al aplicar el Modelo Modificado de Masa Puntual para proyectiles base burn existe un conjunto de coeficientes adicionales relativos al proceso de quemado del propulsante, así como a la disminución de la resistencia producida por el aumento de la presión en la base. Estos coeficientes base burn pueden obtenerse de ensayos en laboratorio o bien de bases de datos suministradas por el fabricante del proyectil. No es objeto de este trabajo la determinación de los coeficientes aerodinámicos, dado que este problema se encuentra resuelto en la actualidad de una forma u otra, como se ha descrito antes. No obstante, la obtención de los coeficientes base burn representa un problema importante y, constituye el objeto y las aportaciones de esta tesis. No obstante, en pocas ocasiones se tiene conocimiento de esta información que el XII.-5 fabricante se reserva como oro en paño siendo la única alternativa posible el obtener tales funciones de ensayos en laboratorios especializados. Estos ensayos son muy costosos y únicamente se lo pueden permitir países con elevados presupuestos. Una alternativa a esto, mucho más barata y que puede resultar adecuada para Unidades Militares con bajos presupuestos, es la determinación de un método que permita el cálculo de estos factores de ajuste a través de los resultados obtenidos de experiencias. Con los datos obtenidos de experiencias realizadas, se van a calcular unas funciones que mejor aproximen estos mismos datos de experiencias. El alcance y la deriva constituyen los parámetros prioritarios a ajustar en el proceso de optimización, los cuales han sido escogidos por ser los más representativos del tiro de artillería. Como ya se comentó con anterioridad, la adición a los anteriores de la ordenada máxima ha permitido constatar un empeoramiento en el cálculo de los coeficientes. Si se obtienen unos coeficientes que mejor ajusten alcance y deriva, y unos segundos coeficientes que mejor ajusten alcance, deriva y ordenada máxima; se ha observado; que si se calcula alcance y deriva con los segundos coeficientes, se obtienen peores valores que usando los primeros coeficientes. Esto es debido a que la ordenada máxima es un parámetro que ante variaciones en los coeficientes, experimenta variaciones mucho mayores que las apreciadas para alcance y deriva. También se ha constatado la dificultad por parte del Modelo Modificado de Masa Puntual de proporcionar valores acertados de la ordenada máxima, una vez ajustada la deriva y el alcance, con datos extraídos de experiencias reales. El alcance y la deriva son los parámetros más representativos del tiro de artillería y por tanto el principal objetivo en este trabajo es su ajuste; no obstante, con posterioridad se podrá comprobar como va a ser necesario otro parámetro obtenido de experiencias, el tiempo de extinción del motor base burn. Sin este dato no se podrá determinar apropiadamente la velocidad de combustión del propulsante base burn y habría que recurrir a algún otro tipo de artificio. Así, el Modelo Modificado de Masa Puntual para proyectiles base burn involucra el conocimiento de una serie de funciones que modelizan el comportamiento del proyectil, y que implican un conocimiento preciso del proceso de quemado del propulsante así como de la reducción de la resistencia al incrementar la presión en la base. La determinación de un método que permita el cálculo de estas funciones a través de los resultados obtenidos de experiencias realizadas constituye un problema de optimación XII.-6 complejo para lo cual se han desarrollado una serie de herramientas numéricas. El incremento del número de parámetros de los que dependen las funciones a optimizar hace de ardua aplicación los algoritmos clásicos, lo que ha dado lugar a la utilización de algoritmos heurísticos. De este modo, se han implementado diversos procedimientos de búsqueda de soluciones que tenían un objetivo común; resolver un determinado problema para el que era muy costoso obtener su solución exacta mediante los métodos de optimización puramente analíticos. La dificultad de encontrar el mejor óptimo mediante métodos clásicos en situaciones donde se tienen múltiples óptimos resulta evidente. Si a esto se le añade el hecho de que la función a optimizar tiene un número elevado de variables de las que depende, la complicación del problema crece substancialmente. Existen muchos problemas de optimización que son muy difíciles de resolver de manera exacta. En estos casos, se plantean algoritmos que permitan obtener soluciones cercanas al valor óptimo de la función objetivo, invirtiendo un tiempo razonable. Este tipo de algoritmos se denominan heurísticos, siendo útiles en contextos tales que no se conoce un procedimiento exacto de resolución del problema, o bien éste requiere un gran esfuerzo computacional. En ocasiones no es necesario obtener la solución óptima global del problema, siendo suficiente con conocer una solución cercana a dicho óptimo. Finalmente, y a efectos de ajustes finales en el cálculo de trayectorias se deben encontrar cuatro parámetros denominados factores de ajuste, como se indica en los STANAG 4355 y STANAG 4144, para la correcta modelización numérica según el Modelo Modificado de Masa Puntual para proyectiles base burn. Este ejercicio se resolverá aplicando los mismos algoritmos comentados anteriormente pero particularizados a esta situación aunque, como puede preverse, constituye un problema de rango inferior. XII.-7 3. Proyectiles base burn método 1 Se enumeran a continuación los coeficientes base burn para proyectiles base burn método 1. Estos coeficientes se han separado en dos grupos atendiendo a su dificultad, los primeros de ellos son los más simples y se reflejan a continuación. Los segundos merecen un estudio más exhaustivo. Símbolo IXB mCB0 mf ρp XCG0 XCGB Descripción Momento de inercia axial al finalizar el quemado del propulsante Masa de propulsante del motor quemada en el tubo Masa de propulsante del motor Densidad del propulsante del motor base burn Distancia inicial del centro de masas a la ojiva Distancia del centro de masas a la ojiva, al finalizar el quemado del propulsante Texto IXB Unidad Dependencia kgm2 MCB0 kg MFUEL kg RHOP kg/m3 XCG0 m XCGB m carga Si no se conocen los momentos de inercia, tanto el momento de inercia axial IX como el momento de inercia axial al finalizar el quemado del propulsante IXB, pueden estimarse a partir de proyectiles asistidos semejantes, como ya se estudió en análisis previos. La masa de propulsante del motor quemada en el tubo mCB0 puede estimarse de proyectiles asistidos semejantes; no obstante, su valor es muy pequeño y puede considerarse nulo en primera aproximación. El ajuste del resto de coeficientes absorberá la diferencia mínima que se produzca. Usualmente, la masa de propulsante del motor mf es dato pero si no lo fuera, podría estimarse dentro del proceso de optimización con el resto de coeficientes; no obstante, la masa de propulsante del motor depende del tiempo de extinción del motor base burn. Este dato resulta fácilmente medible mediante diversos dispositivos empleados en las experiencias. XII.-8 La densidad del propulsante del motor base burn ρp es dato una vez conocido el tipo de propulsante. Si no se conociera, se puede estimar de proyectiles base burn semejantes o bien de tablas: XII.-9 La distancia inicial del centro de masas a la ojiva XCG0 y la distancia del centro de masas a la ojiva al finalizar el quemado del propulsante XCGB, pueden estimarse a partir de proyectiles asistidos semejantes, como ya se estudió en análisis previos Los coeficientes base burn que se relatan a continuación son los más complejos de determinar y ya en si mismos constituyen un estudio de relevancia mayor, por lo que se examinaran cada uno por separado a fin de efectuar un estudio exhaustivo. Símbolo CxBB I0 Descripción Texto Coeficiente de reducción de CXBB resistencia BB Parámetro de inyección de I0 eficiencia óptima del motor base burn XII.-10 Unidad - Dependencia Mach Mach Símbolo SC Descripción Texto Área de combustión del SC propulsante base burn Velocidad de combustión VC del propulsante base burn VC Texto CXBB Unidad 2 m m/s Dependencia propulsante quemado MT, p, presión Nota Polinomio función del número de Mach a0 + a1M + a2M2 + a3M3 + a4M4 para diferentes intervalos del número de Mach. Los intervalos deben introducirse en orden ascendente con el formato I0 CXBB Mmin Mmax a0 a1 a2 a3 Polinomio lineal función del número de Mach a4 a0 + a1M para diferentes intervalos del número de Mach. Los intervalos deben introducirse en orden ascendente con el formato I0 SC Mmin Mmax a0 a1 Polinomio lineal función de la masa de propulsante quemada en el motor (mCB) a0 + a1 mCB para diferentes intervalos de mCB. Los intervalos deben introducirse en orden ascendente con el formato VC SC mCB,min a0 a1 mCB,max La velocidad de combustión viene dada por VC0eβ(MT-21) kPnK(p) donde MT es la temperatura de propulsante del motor (°C), P es la presión ambiente (Pascal) y p es velocidad de rotación axial del proyectil (rad/s). K(p) tiene una entrada separada. El rango de temperaturas validas puede introducirse con el formato VC VC0 β k n MTmin MTmax Para establecer correspondencia entre los resultados obtenidos computacionalmente y las experiencias observadas, se aplican los factores de ajuste. XII.-11 Símbolo f(iBB, MT) Descripción Factor base burn K(p) Factor de quemado con la KP velocidad de rotación axial Tiempo de retraso a la TDI ignición del motor base burn tDI Texto FIBB Texto Unidad - FIBB - Dependencia MT, QE, carga carga, p s MT, carga Nota Polinomio función de la temperatura de propulsante del motor y del ángulo de tiro, para cada carga a0 + a1QE + a2QE2 + a3QE3 + b1(MT-21) + b2(MT-21)2 + b3(MT-21)3 para diferentes intervalos de MT (temperatura de propulsante del motor en °C) y QE (ángulo de tiro en milésimas), para una carga dada. Los intervalos de temperatura deben introducirse en orden ascendente con el formato KP FIBB MTmin MTmax a0 a1 a2 a3 b1 b2 b3 Polinomio lineal función de la velocidad de rotación axial del proyectil, para cada carga a0 + a1p Donde p es la velocidad de rotación axial (rad/s). El formato es TDI KP a0 a1 Polinomio función de la temperatura de propulsante del motor, para cada carga a0 + a1 (MT - 21) + a2 (MT - 21)2 + a3 (MT - 21)3 para diferentes intervalos de MT (temperatura de propulsante del motor en °C), para una carga dada. Los intervalos de temperatura deben introducirse en orden ascendente con el formato TDI MTmin MTmax a0 a1 a2 a3 Parámetro Símbolo Tiempo de retraso a la tDI ignición del motor base burn Factor de quemado con la K(p) velocidad de rotación axial Factor base burn f(iBB,MT) XII.-12 Unidad Límites típicos s - - - - - 3.1. Área de combustión del propulsante base burn SC El área de combustión del propulsante base burn consiste en el conjunto de todas las superficies del propulsante que estén expuestas a la combustión (y que no están inhibidas de alguna manera). El área de combustión del propulsante base burn depende de la geometría del grano y del uso de inhibidores. El propulsante utilizado puede ser una composición simple con dos constituyentes principales; combustible y oxidante, aunque pueden encontrarse propulsantes compuestos con una composición más compleja y contener oxidantes de varios tipos. Pueden existir otros aditivos en pequeños porcentajes que pueden controlar la velocidad de quemado, acelerándola o ralentizándola. Sin importar su composición, todos los propulsantes son procesados en una forma geométrica similar llamada grano propulsante. Como regla, los granos propulsantes son de forma cilíndrica para encajar perfectamente dentro de la unidad base burn con el fin de maximizar la eficiencia volumétrica. El núcleo puede tener una amplia variedad de formas como circular, estrella, cruz, hueso, etc. aunque la forma más común es la circular. La forma del núcleo tiene una profunda influencia en la forma del perfil de la derivada de la masa de propulsante del motor en función del tiempo. La derivada de la masa de propulsante del motor m& f es proporcional al área de combustión del propulsante base burn en cualquier instante de tiempo. El área de combustión del propulsante en cualquier punto sigue la dirección normal a la superficie en ese punto y ocurre de manera uniforme a lo largo de toda el área superficial expuesta a la combustión, siendo el resultado una relación entre el área de combustión del propulsante y la distancia de avance quemada, la cual depende casi exclusivamente de la forma inicial del grano y de los inhibidores. Así, conocida la geometría de la unidad base burn y la velocidad de quemado del propulsante, es posible calcular el área de combustión del propulsante base burn en XII.-13 función de la masa de propulsante quemado SC(mCB): SC = ai + bi ⋅ mCB mB = m0 − m f mCB = m0 − m mCBi < mCB ≤ mCBi+1 donde ai y bi se definen sobre regiones de mCB desde mCBi=0 hasta mCBi=n inclusive. Si no se conoce la geometría de la unidad base burn, puede estimarse el área de combustión del propulsante base burn SC(mCB) de proyectiles base burn semejantes, partiendo de la hipótesis de que las unidades base burn son similares. Dado que el objetivo prioritario es encontrar los parámetros que mejor ajusten alcance y deriva dentro de un proceso de optimización, cabe preguntarse hasta que punto el hecho de sustituir virtualmente una unidad base burn por otra de un proyectil base burn semejante afectará al alcance y a la deriva. Dado que la masa de propulsante del proyectil mf varía de unos proyectiles a otros y que el área de combustión del propulsante vendrá afectado por la masa de propulsante del proyectil, es necesario aplicar la teoría de análisis dimensional a la ecuación que expresa el área de combustión del propulsante base burn como función de la masa de propulsante quemado. Aplicando el teorema π de Vaschy y Buckingham, lo más lógico es adimensionalizar la masa de propulsante quemada en el motor mCB con la masa de propulsante del proyectil mf. En cuanto al área de combustión del propulsante base burn, se va a adimensionalizar con una superficie de referencia que se analizará posteriormente. ⎛ b ⎞ m SC a = i + ⎜ i m f ⎟ ⋅ CB ⎟ m S ref S ref ⎜⎝ S ref f ⎠ XII.-14 Reescribiendo la ecuación anterior en términos de magnitudes adimensionales: S C = S C* ⋅ S ref , ai* = b ai m * , bi* = i m f , mCB = CB S ref S ref mf ( ) * * S C* mCB = ai* + bi* ⋅ mCB Si la unidad base burn se concentra en el culote del proyectil en forma de disco, la superficie de referencia Sref debería ser el diámetro de la base del proyectil al cuadrado. S ref = d b2 Aplicando la anterior ecuación a los proyectiles base burn semejantes, se tiene: La separación mostrada por el proyectil OE155BONUS FR CBL/DUP hace pensar que la unidad base burn no se concentra en un disco en el culote. Dado que no se conoce XII.-15 como es la geometría de las unidades base burn y que se desea obtener una expresión que, con carácter general, sirva para todos los proyectiles base burn, se va a repetir el proceso anterior suponiendo que la unidad base burn no tiene una dimensión predominante presentando un aspecto de cuasiesfera. Así, considerando la superficie de referencia Sref como S ref ⎛ mf =⎜ ⎜ρ ⎝ p ⎞ ⎟ ⎟ ⎠ 2 3 y aplicando esta superficie de referencia a los proyectiles base burn semejantes, se tiene: Donde se puede comprobar que el resultado obtenido presenta una agrupación más uniforme. El siguiente paso es investigar como afectaría el sustituir virtualmente una unidad base burn por otra de un proyectil base burn semejante en el cálculo del alcance y la deriva, dado que el objetivo prioritario es determinar los parámetros que mejor ajusten alcance y deriva al aplicar un algoritmo de optimación. Considerando todas las unidades base burn disponibles (se ha descartado la del proyectil noruego OEF3BB NO XII.-16 HEA BB1 ya que es la misma que la del proyectil francés OE155F2RTC FR HEA BB1) y disparando con todas las cargas y a cinco ángulos de tiro distintos (300, 500, 700, 900 y 1200 milésimas), se ha calculado la diferencia porcentual en alcance y deriva al intercambiar las unidades base burn: (alcance – alcanceBB)/alcance y (deriva – derivaBB)/deriva. Estos resultados se han resumido en la tabla que se anexa a continuación reflejando únicamente el mayor valor encontrado para todas las cargas y ángulos de tiro analizados: Alcance Proyectil OE155F2RTC Unidad base burn OE155BONUS OE155F2RTC FR HEA BB1 0.00 % 0.68 % OE155BONUS FR CBL/DUP BB1 1.12 % 0.00 % ER-02/BB SP HEA BB1 - 0.79 % 0.34 % Deriva Proyectil OE155F2RTC Unidad base burn OE155BONUS OE155F2RTC FR HEA BB1 0.00 % 1.29 % OE155BONUS FR CBL/DUP BB1 -1.25 % 0.00 % ER-02/BB SP HEA BB1 - 1.69 % 0.48 % A la vista de las diferencias obtenidas y del aspecto que presentan las diferentes áreas de combustión del propulsante base burn SC(mCB), se plantea reducir aún más estas diferencias efectuando una homotecia de dirección el eje vertical sobre la curva del área de combustión del propulsante base burn adimensional de un proyectil base burn semejante introduciendo un nuevo coeficiente que deberá ser ajustado. ( ) ( ) * * S C* mCB = S C0 ⋅ S C* mCB proyectil⋅semejante Aunque pudiera parecer que es necesario disponer de información proporcionada por proyectiles base burn semejantes; sin embargo, la robustez del método permite trabajar XII.-17 sobre otro tipo de curva que se plantee como área de combustión de propulsante base burn. Avanzando sobre este concepto y recordando la definición de la derivada de la masa de propulsante del motor m& f m& = − m& f m& = −VC ⋅ ρ p ⋅ SC (mCB ) VC = VC0 ⋅ f (MT ) ⋅ g (P ) ⋅ K ( p ) f (MT ) = e β ( MT − 21) g (P ) = k ⋅ P n Rearreglando, se tiene m& = −VC0 ⋅ e β ( MT − 21) ⋅ k ⋅ P n ⋅ K ( p ) ⋅ ρ p ⋅ S C (mCB ) donde no se conocen los factores VC0, β, k, y n. K(p) es un factor de ajuste y se tratará posteriormente y ρp ya se trató anteriormente. El factor β tiene en cuenta la desviación producida cuando la temperatura de propulsante del motor MT se separa de los 21 º C y su análisis requiere del conocimiento de experiencias efectuadas a diferentes temperaturas de propulsante del motor, por lo que no ha sido objeto de este trabajo, pudiendo estimarse de proyectiles base burn semejantes aunque si se dispara a temperaturas cercanas a 21 º C el efecto es nulo. Si se considera, análogamente a como se hizo antes, introducir dos nuevos coeficientes basados en proyectiles base burn semejantes XII.-18 ( ) VC0 = VC00 ⋅ VC0 proyectil ⋅semejante k = k0 ⋅ k proyectil ⋅ semejante Obsérvese que esto únicamente constituye un cambio en la escala, dado que VC0 y k son escalares y se introducen únicamente para identificar más apropiadamente el factor SC0. Introduciendo las expresiones anteriores, se tiene ( ) m& = −VC 00 VC 0 proyectil ⋅ semejante ( ) * ⋅ e β ( MT − 21) ⋅ k0 ⋅ k proyectil ⋅ semejante ⋅ P n ⋅ K ( p ) ⋅ ρ p ⋅ SC0 ⋅ S ref ⋅ SC* mCB proyectil ⋅ semejante Agrupando los términos a ser ajustados en uno único S C0 = VC00 ⋅ k 0 ⋅ S C0 {( ) m& = − SC0 VC 0 p⋅s [ ] } [ ( ) ] * ⋅ e β ( MT − 21) ⋅ k p ⋅ s ⋅ P n ⋅ K ( p ) ⋅ ρ p ⋅ S ref ⋅ SC* mCB p⋅s Así, finalmente se tiene que el área de combustión del propulsante base burn depende directamente de los términos SC0 y n. Estos dos elementos constituyen los parámetros a ser ajustados en el proceso de optimización posterior que han sido aportados por el área de combustión del propulsante base burn. Los análisis efectuados con la información existente muestran que el término SC0 puede sufrir variaciones del orden del 10 % mientras que n toma valores del orden de 0.76 , si bien estos valores son valores típicos únicamente ya que al aplicar el algoritmo de optimización correspondiente se debe permitir un mayor movimiento para aumentar la probabilidad del éxito. La velocidad de combustión de propulsante base burn VC es muy sensible al efecto del término n, condicionando enormemente el tiempo de extinción del motor base burn. Este efecto será necesario tenerlo en cuenta a la hora de seleccionar apropiadamente la función objetivo en el proceso de optimación posterior. XII.-19 3.2. Coeficiente de reducción de resistencia base burn CxBB El cambio en aceleración, debida a la reducción de resistencia de base de un motor base burn BB& , durante la fase de quemado ( t DI ≤ t y m ≥ mB , siendo tDI el tiempo de retraso a la ignición del motor base burn y mB la masa del proyectil al finalizar el quemado del propulsante), se adiciona a la ecuación del movimiento del centro de masas del proyectil r r F = mu&& : ⎤ ⎡⎛ π ⎞ 2 2 ⎢ ⎜ 8 ⎟ ⋅ ρ ⋅ d ⋅ v ⋅ C xBB ⋅ f (I ) ⋅ f (iBB , MT ) ⎥⎛ vr ⋅ cos α r ⎞ e ⎥⎜ BB& = ⎢ ⎝ ⎠ +αe ⎟ m v ⎥⎝ ⎢ ⎠ ⎥ ⎢ ⎦ ⎣ El coeficiente CxBB es el coeficiente de reducción de resistencia durante la fase de quemado. Se conceptúa como cualquier otro coeficiente aerodinámico, los valores de estos coeficientes vienen dados por funciones polinómicas de orden cuatro (o inferior) del número de Mach. Resultan de gran interés los estudios realizados por Chargelegue y Couloumy de la Direction de Armaments Terrestres a finales de los ochenta. Empleando un radar doppler de seguimiento se dispararon setenta y cuatro proyectiles base burn con diferentes cargas, ángulos de tiro y temperaturas de propulsante del motor. Mediante un programa informático se obtenía el coeficiente de reducción de resistencia base burn en función del número de Mach. Estos trabajos fueron expuestos en el 11º Simposio Internacional de Balística celebrado en Bruselas en 1989. En este aspecto, debe tenerse en cuenta que estos estudios se realizaron en 1986 y 1987, siendo precursores del STANAG 4355, y empleando una expresión ligeramente simplificada del cambio en aceleración debida a la reducción de resistencia de base del motor base burn durante la fase de quemado BB& : BB& = r 1 ρ ⋅ S ⋅ v ⋅ v ⋅ C xBB (M , I , P, MT ,K) 2m XII.-20 Designando como CxBB(M,I,P,MT,…) el coeficiente de reducción de resistencia base burn manejado por Chargelegue y Couloumy y CxBB el coeficiente de reducción de resistencia base burn que aparece en el STANAG 4355, se pueden comparar ambos términos. r r v ⎛ v ⋅ cos α e r ⎞ C xBB (M , I , P, MT ,K) = C xBB ⋅ f (I ) ⋅ f (iBB , MT )⎜ + αe ⎟ v v ⎝ ⎠ siendo M el número de Mach, I el parámetro de inyección de propulsante en el motor base burn, P la presión del aire atmosférico local y MT la temperatura de propulsante del motor en °C. A efectos simplificatorios se ha supuesto que la superficie de referencia S se corresponde con π d 2 4 . Si se considera el ángulo de ataque αe suficientemente pequeño, se tiene en primera aproximación C xBB (M , I , P, MT ,K) ≈ C xBB ⋅ f (I ) ⋅ f (iBB , MT ) Puede observarse como el coeficiente de reducción de resistencia base burn empleado por Chargelegue y Couloumy CxBB(M,I,P,MT,…) constituye un coeficiente “global” de reducción de resistencia base burn. Estos investigaron ya indicaron entonces la existencia de un flujo de masa óptimo, por encima del cual, este coeficiente “global” de reducción de resistencia base burn CxBB(M,I,P,MT,…) solo depende del número de Mach. Chargelegue y Couloumy obtuvieron la siguiente figura que presenta el coeficiente global de reducción de resistencia base burn CxBB(M,I,P,MT,…), para un proyectil DTC, con temperatura de propulsante del motor de 21 º C, mostrando la dependencia que presenta CxBB(M,I,P,MT,…) con el número de Mach. XII.-21 Chargelegue y Couloumy pensaron que como la velocidad de combustión de propulsante base burn VC está afectada por la temperatura de propulsante del motor, podría ocurrir que el coeficiente global de reducción de resistencia base burn CxBB(M,I,P,MT,…) también estuviese afectado por la temperatura de propulsante del motor. Por ello, representaron el coeficiente global de reducción de resistencia base burn CxBB(M,I,P,MT,…) a diferentes temperaturas de propulsante del motor, mostrando que el efecto de la temperatura modifica muy débilmente la dispersión del coeficiente global de reducción de resistencia base burn. Separando este efecto, Chargelegue y Couloumy sugirieron que existe un flujo de masa óptimo por encima del cual, el coeficiente global de reducción de resistencia base burn CxBB(M,I,P,MT,…) solo depende del número de Mach. XII.-22 Finalmente, Chargelegue y Couloumy obtuvieron el coeficiente global de reducción de resistencia base burn CxBB(M,I,P,MT,…) a través de un procedimiento mediante la diferencia de dos términos obtenidos de experimentos; Cx(M,I,P,MT,…) con bloque propulsante a 21 º C (CxABB) y Cx(M,I,P,MT,…) sin bloque propulsante (CxSBB) C x BB (M , I , P, MT ,K) = C x SBB − C x ABB Las tres gráficas CxBB(M,I,P,MT,…), CxABB y CxSBB, como función del número de Mach, pueden verse en la siguiente figura. CxBB(M,I,P,MT,…) se expresa como un conjunto de funciones polinómicas del número de Mach con grado inferior a cuatro. A continuación se puede examinar el coeficiente de reducción de resistencia base burn CxBB de los proyectiles de artillería OE155F2RTC FR HEA BB1, OE155BONUS FR CBL/DUP BB1 y OEF3BB NO HEA BB1. XII.-23 XII.-24 En las anteriores gráficas puede observarse el incorrecto acoplamiento de las curvas en la zona intermedia próxima al número de Mach unidad con el consiguiente problema de valores anómalos. Con el método desarrollado en esta tesis, este problema no se producirá. Si no se conoce el coeficiente de reducción de resistencia base burn CxBB, puede efectuarse una estimación de esta función a partir de alguna función similar conocida. El método de Rayleigh Ritz plantea sustituir la función buscada por una combinación lineal de funciones linealmente independientes bien condicionadas y conocidas, que están ponderadas por unos coeficientes a determinar mediante métodos variacionales ψ = c1 ⋅ f1 + K + cn ⋅ f n Para la determinación de los n coeficientes hay que imponer n condiciones. Primeramente debe decidirse acerca de cual es el valor más adecuado para el número n. Téngase en cuenta que si bien incrementar este valor mejorará la definición de la curva, XII.-25 también impondrá unas penalizaciones excesivas al necesitar más ecuaciones para determinar los coeficientes. El incremento del número de parámetros de los que depende la función a optimizar, provoca de una manera casi indefectible la aplicación de algoritmos heurísticos, dada la enorme dificultad que supone el obtener la solución exacta mediante métodos de optimización puramente analíticos. No obstante, debe tenerse en cuenta que aunque estas herramientas tienen una gran potencia, el entorno de búsqueda es tremendamente hostil con múltiples óptimos cercanos. Esto permite asegurar que el procedimiento de resolución será mucho más favorable en la medida que el número de parámetros de los que depende la función a optimizar sea lo más reducido posible. Así, se debe establecer una solución de compromiso de forma que la curva se encuentre suficientemente bien definida y representada con un número no muy elevado de parámetros. Finalmente se ha optado por una representación definida por cuatro puntos; (M0, CxBB0), (M1, CxBB1), (M2, CxBB2), y (M3, CxBB3) como se indica en la figura. La elección de los números de Mach M0, M1, M2 y M3 dependerá del problema concreto a resolver; no obstante, M0 se toma un número de Mach ligeramente inferior a la unidad y M1 un número de Mach ligeramente superior a la unidad. La curva se ajusta entre M1 y M3 con una parábola que pasa por los puntos 1, 2 y 3. Entre M0 y M1 se emplea otra parábola que pasa por los puntos 0 y 1 y tenga la misma pendiente en el punto 1 que la primera parábola. Para números de Mach inferiores a M0 se toma un valor constante y en el punto 0 se utiliza un radio de acuerdo que suavice el cambio de tramo. El objetivo prioritario es encontrar los parámetros que mejor ajusten alcance y deriva, dentro de un proceso de optimización. No obstante, el número de parámetros a ajustar en el proceso de optimización posterior puede resultar ligeramente alto lo que va a provocar que la función de optimización presente un entorno complicado con múltiples óptimos cercanos. XII.-26 XII.-27 El incremento del número de parámetros de los que depende la función de optimización configura un entorno de búsqueda complejo con múltiples óptimos cercanos. Esto direcciona a la necesidad de utilizar algoritmos heurísticos, ante la gran dificultad de aplicar métodos clásicos analíticos. En este sentido, aunque las herramientas heurísticas tienen una gran capacidad para resolver problemas, es necesario establecer medidas de apoyo a la resolución. Se debe limitar cuidadosamente el campo de actuación de cada uno de los parámetros de los que depende la función objetivo, así como las condiciones iniciales de iteración, para lo cual es necesario disponer de una extensa experiencia, disponiendo de un grupo de proyectiles base burn semejantes convenientemente bien caracterizados. En este sentido, se han instalado una serie de filtros a fin de definir más apropiadamente la naturaleza del problema; CxBB3 < CxBB2, CxBB2 < CxBB1, CxBB0 < CxBB1, la pendiente del segmento CxBB1CxBB2 debe ser inferior a la pendiente del segmento CxBB2CxBB3 (parábola 123 cóncava) y la pendiente entre el número de Mach cero y el M0 es nula. Elegidos los números de Mach M0, M1, M2 y M3 de acuerdo al problema concreto a resolver, que vendrá influenciado por la propia experiencia que se tenga en la XII.-28 configuración proyectil arma carga, se deben encontrar los parámetros CxBB0, CxBB1, CxBB2, y CxBB3 que mejor se adapten a los datos extraídos de experiencias en un proceso de optimización. El siguiente punto es investigar hasta que punto el hecho de sustituir la curva original correspondiente al coeficiente de reducción de resistencia base burn CxBB por una curva parabólica por tramos similar, afectará al alcance y a la deriva. Considerando los proyectiles OE155F2RTC FR HEA BB1 y OE155BONUS FR CBL/DUP BB1 (no se ha considerado el proyectil noruego OEF3BB NO HEA BB1 ya que es el mismo que el proyectil francés OE155F2RTC FR HEA BB1) y disparando con todas las cargas y a cinco ángulos de tiro distintos (300, 500, 700, 900 y 1200 milésimas), se ha calculado la diferencia porcentual entre utilizar el coeficiente de reducción de resistencia base burn CxBB por la curva parabólica por tramos indicada anteriormente; (alcance – alcancepol)/alcance y (deriva – derivapol)/deriva. Estos resultados se han resumido en la tabla que se anexa a continuación reflejando únicamente el mayor valor encontrado para todas las cargas y ángulos de tiro analizados: Proyectil Alcance Deriva OE155F2RTC FR HEA BB1 0.08 % 0.22 % OE155BONUS FR CBL/DUP BB1 - 0.49 % - 0.40 % Con respecto a las diferencias obtenidas en la tabla anterior, se debe indicar que son meramente ilustrativas, dado que la curva parabólica por tramos con la que se han efectuado los cálculos es una curva obtenida geométricamente a partir de las curvas originales. Realmente, la determinación de los parámetros CxBB0, CxBB1, CxBB2, y CxBB3 que mejor ajusten alcance y deriva, empleando el algoritmo de optimización desarrollado, mostrará posteriormente que estas diferencias son prácticamente nulas. 3.3. Parámetro de inyección de eficiencia óptima I0 La función de variación del flujo del motor base burn viene dada por: XII.-29 f (I ) = I I 0 si I ≤ I 0 f (I ) = 1 si I ≥ I 0 donde I es el parámetro de inyección de propulsante en el motor base burn I= 4 ⋅ m& f π ⋅ d b2 ⋅ ρ ⋅ v donde m& f es la derivada de la masa de propulsante del motor y db el diámetro de la base del proyectil. I0 es el parámetro de inyección de eficiencia óptima del motor base burn, que es función del número de Mach y que se puede tomar como un polinomio de un grado o inferior. Se define sobre regiones de números de Mach, desde M MAX i−1 hasta M MAX i inclusive. Como ya se comentó anteriormente, los estudios realizados por Chargelegue y Couloumy apuntaron que existe un flujo de masa óptimo por encima del cual, el coeficiente global de reducción de resistencia base burn CxBB(M,I,P,MT,…) solo depende del número de Mach. Estos investigadores utilizan el primer grupo de disparos para calcular el parámetro de inyección de eficiencia óptima del motor base burn I0. XII.-30 Cuando el parámetro de inyección de propulsante en el motor base burn I es mayor que el parámetro de inyección de eficiencia óptima I0, el coeficiente “global” de reducción de resistencia base burn solo depende del número de Mach. Cuando empieza a agotarse el propulsante, el parámetro de inyección de propulsante en el motor base burn I se vuelve más pequeño que el parámetro de inyección de eficiencia óptima I0, y este coeficiente “global” de reducción de resistencia base burn cae rápidamente. La función de variación del flujo del motor base burn f(I) se usa para modelizar este comportamiento de caída brusca del coeficiente “global” de reducción de resistencia base burn: C xBB (M , I , P, MT ,K) ≈ C xBB ⋅ f (I ) ⋅ f (iBB , MT ) A continuación se puede examinar el parámetro de inyección de eficiencia óptima I0 de los proyectiles de artillería OE155F2RTC FR HEA BB1 y OE155BONUS FR CBL/DUP BB1 (no se ha considerado el proyectil noruego OEF3BB NO HEA BB1 ya que es el mismo que el proyectil francés OE155F2RTC FR HEA BB1). XII.-31 Donde puede apreciarse como el parámetro de inyección de eficiencia óptima I0 toma la forma de una poligonal con pendientes negativas. Al principio de la trayectoria, inmediatamente después del momento del disparo, el parámetro de inyección de propulsante en el motor base burn I es mayor que el parámetro de inyección de eficiencia óptima I0, de forma que la función de variación del flujo del motor base burn es la unidad siendo el coeficiente global de reducción de resistencia base burn una función dependiente únicamente del número de Mach. Cuando el propulsante empieza a agotarse, el parámetro de inyección de propulsante en el motor base burn I se vuelve más pequeño que el parámetro de inyección de eficiencia óptima I0, y el coeficiente global de reducción de resistencia base burn cae rápidamente, siendo la función de variación del flujo del motor base burn f(I) la que refleja este comportamiento de caída brusca del coeficiente global de reducción de resistencia base burn. Esta forma de actuar viene regulada por el parámetro de inyección de eficiencia óptima I0 que debe proceder de la manera que se ha relatado anteriormente para todas las cargas y todos los ángulos de tiro. Una mala elección del parámetro de inyección de XII.-32 eficiencia óptima I0 puede provocar que al comienzo de la trayectoria el parámetro de inyección de propulsante en el motor base burn I se vuelva más pequeño que el parámetro de inyección de eficiencia óptima I0 para alguna carga o algún ángulo de tiro, produciendo una caída abrupta al comienzo de la trayectoria totalmente contraria a lo que se pretendía. Al separar el coeficiente “global” de reducción de resistencia base burn CxBB(M,I,P,MT,…) en tres términos; el coeficiente de reducción de resistencia base burn CxBB, la función de variación del flujo del motor base burn f(I) y el factor base burn f(iΒΒ,MT) (factor de ajuste que será calculado en un proceso posterior), se debe calcular primeramente la función de variación del flujo del motor base burn f(I) cumpliendo las condiciones que se han descrito en el párrafo anterior, para posteriormente determinar el coeficiente de reducción de resistencia base burn CxBB. C xBB (M , I , P, MT ,K) ≈ C xBB ⋅ f (I ) ⋅ f (iBB , MT ) De este modo, la función de variación del flujo del motor base burn f(I) sirve para expresar la caída brusca que sufre el coeficiente global de reducción de resistencia base burn cuando el propulsante del motor empieza a agotarse. Así, el coeficiente de reducción de resistencia base burn CxBB es el encargado de expresar un ajuste exacto a lo largo de la trayectoria del proyectil a partir de la función de variación del flujo del motor base burn f(I). Esto no significa que la determinación de la función de variación del flujo del motor base burn f(I) sea un trabajo simple ya que la caída abrupta que soporta el coeficiente global de reducción de resistencia base burn debe ocurrir para todas las cargas y todos los ángulos de tiro, y nunca al comienzo de la trayectoria. Una vez analizado el comportamiento del parámetro de inyección de eficiencia óptima I0, la acepción óptima no parece la más apropiada para el parámetro I0 dado su funcionamiento de regulación que ejerce sobre el coeficiente global de reducción de resistencia base burn CxBB(M,I,P,MT,…). Este calificativo puede provenir de una traducción poco apropiada dado que Chargelegue y Couloumy lo designan como parámetro de inyección característico I0, mucho más apropiado para la función realizada. XII.-33 A continuación se puede observar la comparativa entre el parámetro de inyección de propulsante en el motor base burn I y el parámetro de inyección de eficiencia óptima I0 a lo largo de la trayectoria. En estas figuras puede apreciarse el fenómeno comentado anteriormente. Al principio de la trayectoria, el parámetro de inyección de propulsante en el motor base burn I es superior que el parámetro de inyección de eficiencia óptima I0, de forma que la función de variación del flujo del motor base burn es la unidad. Cuando el propulsante empieza a terminarse, el parámetro de inyección de propulsante en el motor base burn I se vuelve más pequeño que el parámetro de inyección de eficiencia óptima I0, y el coeficiente global de reducción de resistencia base burn cae velozmente. Estas gráficas se han impreso para los proyectiles de artillería OE155F2RTC FR HEA BB1 y OE155BONUS FR CBL/DUP BB1 para la carga más alta, la más baja y una intermedia, así como para cinco ángulos de tiro distintos (300, 500, 700, 900 y 1200 milésimas). XII.-34 XII.-35 XII.-36 XII.-37 XII.-38 XII.-39 XII.-40 XII.-41 XII.-42 XII.-43 XII.-44 XII.-45 XII.-46 XII.-47 XII.-48 3.4. Implementación del método de cálculo A continuación se va ha habilitar un método original de cálculo que constituye una de las aportaciones de esta tesis, con el objeto de determinar los diferentes parámetros necesarios para caracterizar apropiadamente el vuelo de un proyectil base burn empleando el método 1. Aunque el algoritmo de recocido simulado (SA) resulta mejor que al algoritmo genético (GA) en cuanto a la precisión de los resultados obtenidos, presenta el inconveniente de que requiere una labor de determinación (y por tanto de prueba y error) muy exhaustiva de los parámetros. Esto implica un mayor tiempo de cálculo al tener que realizar ajustes constantes para obtener buenos resultados. Este inconveniente no lo presenta el algoritmo genético que aunque es menos preciso, es más robusto en este sentido, lo cual condiciona la aplicación de ambos algoritmos. El algoritmo genético es más apropiado para una fase incipiente del problema, cuando se tiene un mayor desconocimiento de los resultados a obtener. Entonces, las soluciones obtenidas utilizando el algoritmo genético se utilizan como valor de entrada para la aplicación del algoritmo de recocido simulado. XII.-49 La sinergia de ambos algoritmos produce los efectos más potentes en la determinación y resolución del problema. El proceso de determinación de los diferentes parámetros que modelan el proceso de quemado del propulsante así como de la reducción de la resistencia al aumentar la presión en la base, constituye un procedimiento multipaso e iterativo donde determinados unos parámetros (A), se debe acometer el cálculo de los siguientes (B), con el consiguiente problema de que cuando se comienza el cálculo a la búsqueda de los parámetros A es necesario presuponer los parámetros B, los cuales obviamente no se disponen. He aquí donde es necesario acudir en ayuda de los proyectiles base burn semejantes para estimar estos parámetros B previos. Cuanto más semejante sean estos proyectiles, menor será el número de iteraciones a realizar. Este proceso debe repetirse hasta alcanzar la precisión deseada. Si no fuera dato la masa de propulsante del motor mf, podría estimarse dentro del proceso de optimización con el resto de coeficientes; no obstante, la masa de propulsante del motor depende del tiempo de extinción del motor base burn; es decir, de la duración de la fase asistida, hecho que está íntimamente relacionado con la velocidad a la que se consume el propulsante en el motor. La velocidad de combustión de propulsante base burn VC es muy sensible al efecto del término n, condicionando enormemente el tiempo de extinción del motor base burn. Un valor grande de n provoca el apagado prematuro del motor base burn mientras que un valor pequeño produce un apagado tardío. Este hecho habrá que tenerlo en cuenta a la hora de seleccionar apropiadamente la función objetivo en el proceso de optimización. En definitiva, la masa de propulsante del motor mf, el tiempo de extinción del motor base burn y el término n que aparece en la expresión de la velocidad de combustión de propulsante base burn VC se encuentran acoplados. Afortunadamente, el tiempo de extinción del motor base burn es un dato fácilmente medible mediante diversos dispositivos empleados en experiencias. Se va a esquematizar todo el proceso de resolución en los siguientes pasos, suponiendo que la masa de propulsante del motor mf es dato: 1. Inicialmente, se introduce el parámetro de inyección de eficiencia óptima XII.-50 I0, teniendo en cuenta que al comienzo de la trayectoria el parámetro de inyección de propulsante en el motor base burn I es mayor que el parámetro de inyección de eficiencia óptima I0, lo cual ocurre hasta que el propulsante del motor empieza a agotarse. En ese momento, el parámetro de inyección de propulsante en el motor base burn I se vuelve más pequeño que el parámetro de inyección de eficiencia óptima I0 y el coeficiente global de reducción de resistencia base burn cae rápidamente. Este comportamiento debe ocurrir para todas las cargas y todos los ángulos de tiro. El resto de parámetros se pueden estimar de un proyectil base burn semejante en la primera etapa del cálculo, si no se disponen de otros datos. 2. Aplicar el algoritmo de optimación para determinar los términos SC0 (para ello se ha estimado previamente un área de combustión del propulsante base burn S C* adimensional), n y 4CxBB (CxBB0, CxBB1, CxBB2, y CxBB3). Obtenido el alcance, la deriva y el tiempo de extinción del motor base burn de experiencias realizadas, se calculan los parámetros SC0, n y 4CxBB que mejor aproximen el alcance, la deriva y el tiempo de extinción del motor base burn. A través de experiencias realizadas se obtiene el alcance, la deriva y el tiempo de extinción del motor base burn, que se designaran como alcancei-exp, derivai-exp y tBi-exp para diferentes ángulos de tiro y cargas, definiendo la función objetivo como i=m ( E = ∑ walc (alcancei − alcancei −exp ) + wder (derivai − derivai −exp ) + wt B t Bi − t Bi−exp 2 2 i =1 alcancei −exp = alcancei −exp (QEi , MVi ) derivai −exp = derivai −exp (QEi , MVi ) t Bi−exp = t Bi−exp (QEi , MVi ) ( alcancei = alcancei QEi , MVi , S C0 , n,4C xBB ( derivai = derivai QEi , MVi , S C0 , n,4C xBB XII.-51 ) ) ) 2 ( t Bi = t Bi QEi , MVi , SC0 , n,4C xBB ) siendo los términos alcancei, derivai y tBi los obtenidos como resultado de haber realizado todo el proceso de cálculo en función de los parámetros SC0, n y 4CxBB. Los parámetros walc, wder y wtB son factores de peso que pueden ser diferentes dependiendo del paréntesis al que vayan a afectar y que permiten favorecer un requisito en perjuicio de los otros. Esto conlleva a que es necesario seleccionar el peso que afecta a cada uno de los paréntesis. Una característica muy importante de este problema es que se conoce a priori el valor óptimo buscado de la función objetivo. 3. Volver al paso 1. y finalizar el proceso cuando se haya alcanzado la precisión requerida. La función objetivo se calcula como un sumatorio que depende de las triada de valores alcancei-exp, derivai-exp y tBi-exp, función a su vez del ángulo de tiro y la carga, siendo m el número total de triadas. Aunque puede preverse a simple vista que cuantos más valores se dispongan de las experiencias realizadas, mejor será el procedimiento de ajuste, también se debe tener en cuenta que mayor será el tiempo de cálculo. No obstante se debe tener en cuenta que el coeficiente de reducción de resistencia base burn CxBB es una función del número de Mach y que alcancei-exp, derivai-exp y tBi-exp dependen del ángulo de tiro y la velocidad en boca del arma. Si se dispara con velocidades iniciales elevadas y el ángulo de tiro no es muy grande, el proyectil no alcanzará números de Mach bajos con lo que se habrá ajustado la parte de la curva CxBB correspondientes a número de Mach altos pero no la parte correspondientes a números de Mach bajos; o dicho de otra forma, en el proceso de optimización los números altos de Mach tendrán un peso de ponderación mucho mayor que los números de Mach bajos pudiendo llegar el caso que el peso de los números bajos de Mach sea muy pequeño. En esta situación, los coeficientes CxBB0 y CxBB1 serán los más afectados ya que el valor obtenido no tendrá mucha precisión en los resultados. Así, se debe intentar que la batería de triada de valores alcancei-exp, derivai-exp y tBi-exp correspondan a ángulos de tiro y cargas tales que XII.-52 los números de Mach altos tengan un peso de ponderación similar a los números de Mach bajos. Esto también conlleva un inconveniente que es preciso evitar; si la trayectoria es corta porque la carga y el ángulo de tiro son bajos, la unidad base burn puede no llegar a su extinción ya que el impacto se produce antes de que se haya consumido todo el propulsante con lo que el tiempo de extinción del motor base burn será erróneo y el procedimiento de optimización funcionará de manera inexacta. Ya se comentó anteriormente que la inclusión de la ordenada máxima en la función objetivo introduce un empeoramiento en los coeficientes ajustados de esta forma. Si con estos segundos coeficientes se calcula alcance y deriva, estos resultados son inferiores a los obtenidos empleando coeficientes ajustados prescindiendo de la ordenada máxima en la función objetivo. La razón es que la ordenada máxima es un parámetro que ante variaciones en los coeficientes, experimenta variaciones mucho mayores que las advertidas para alcance y deriva. Otro factor es que el Modelo Modificado de Masa Puntual presenta dificultades para determinar la ordenada máxima, una vez ajustada la deriva y el alcance, con datos obtenidos de experiencias reales. Si no fuera dato la masa de propulsante del motor mf, el problema se abordaría incorporando a los términos SC0, n y 4CxBB (CxBB0, CxBB1, CxBB2, y CxBB3) el propio mf y realizando todo el proceso de la misma forma. La experiencia adquirida en la realización de este trabajo advierte que el número de parámetros a ajustar de los que depende la función objetivo resulta ligeramente alto (67). En la medida que se incrementan el número de parámetros de los que depende la función a optimizar, el entorno de búsqueda se vuelve más complejo con múltiples óptimos cercanos. Esto establece la necesidad de utilizar nuevos algoritmos ante la casi imposibilidad de aplicar métodos clásicos tipo analíticos. En este sentido, aunque las herramientas heurísticas tienen una gran potencia es necesario establecer normas de apoyo en la resolución. Se debe limitar cuidadosamente el campo de actuación de cada uno de los coeficientes de los que depende la función objetivo así como las condiciones iniciales de iteración, para lo cual es necesario disponer de una amplia experiencia, por lo que es tan importante disponer de una panoplia de proyectiles base burn semejantes convenientemente bien caracterizados. Una vez abordado el problema y en un estadio incipiente de resolución, es más conveniente atacar su resolución empleando el XII.-53 algoritmo genético que aunque no es tan preciso como el algoritmo de recocido simulado requiere una menor carga de trabajo. Es entonces cuando la solución obtenida empleando el algoritmo genético se utiliza como valor de entrada en el algoritmo de recocido simulado, con el objetivo de obtener una mayor precisión en la resolución del problema, ya que aunque este algoritmo es más preciso requiere una mayor labor de determinación (y por tanto de prueba y error) de determinados parámetros, lo que implica mayor tiempo de cálculo al tener que realizar ajustes constantes para obtener buenos resultados. Esta fase es previa a la determinación de los factores de ajuste pero es necesario introducir valores para proceder al cálculo, por lo que para el tiempo de retraso a la ignición del motor base burn tDI se introducirá el valor correspondiente en el supuesto de que se conozca. Sino se conociera, se considerará tDI = 0 siendo el resto de valores K(p) = 1 , f(iBB,MT) = 1 y fL = 1 . Téngase en cuenta que el tiempo de retraso a la ignición del motor base burn tDI es un valor fácilmente medible en experiencias. Como ya se comentó anteriormente, el algoritmo de integración utilizado ha sido el Runge Kutta 45, los parámetros eRK, hmin y hmax se han elegido para obtener una precisión inferior a un metro en el cálculo de la deriva y el alcance: hmin = 0.04 s, hmax = 6.0 s y eRK = 0.1 . Se ha considerado apropiado para este cálculo la siguiente condición de parada; L = 200, K1 = 5, ∈1 = 0.5 , K2 = 5, ∈2 = 0.5 y M = 5000. Una vez configurados adecuadamente los coeficientes base burn, se deben determinar a efectos de ajustes finales los factores de ajuste, para la correcta modelización numérica según el Modelo Modificado de Masa Puntual para proyectiles base burn método 1. Aplicando el algoritmo de optimación, se deben determinar los factores de ajuste tDI, K(p), f(iBB,MT) y fL utilizando para ello datos de experiencias. En este trabajo no se ha tenido en cuenta la dependencia con la temperatura de propulsante del motor MT, utilizando como valor de referencia MT = 21 . La obtención de valores de temperatura de propulsante del motor extraídos de experiencias es una labor ardua, dado que estos valores están sujetos a las condiciones climatológicas del XII.-54 momento, y es necesario disponer de un amplio abanico de valores, por lo que no ha sido objeto de este trabajo, pudiendo estimarse de proyectiles base burn semejantes, aunque si se dispara a temperaturas de propulsante del motor cercanas a 21 º C, el efecto es nulo. Para cada carga; únicamente el factor base burn f(iBB,MT) y el factor de sustentación fL son funciones del ángulo de tiro, siendo el tiempo de retraso a la ignición del motor base burn tDI función de la temperatura de propulsante del motor MT y el factor de quemado con la velocidad de rotación axial del proyectil K(p) una función lineal de la velocidad de rotación axial del proyectil p. El alcance disminuye al aumentar tDI, al aumentar K(p) o/y al disminuir f(iBB,MT), aunque la contribución más importante es esta última por lo que puede decirse que el alcance viene afectado fundamentalmente por el factor base burn. La deriva es muy poco sensible a tDI (aunque disminuye al aumentar tDI) disminuyendo al aumentar K(p) o/y al disminuir f(iBB,MT). De cualquier forma, numéricamente, todos los parámetros se encuentran acoplados. Dado que K(p) depende de cada carga, es necesario disponer de un factor de ajuste que sea capaz de ajustar la deriva para cada ángulo de tiro. El STANAG 4355 indica que i debe considerarse la unidad mientras que fL se utilizará para ajustar la deriva, siendo QM y QD parámetros constantes. En cuanto al tiempo tras la finalización del quemado del propulsante tB; si K(p) aumenta, disminuye el tiempo extinción del motor siendo ésta la contribución más importante. Si f(iBB,MT) aumenta, disminuye tB, pero de manera imperceptible. Si tDI aumenta, tB crece prácticamente como una traslación y si fL aumenta, tB crece pero sin percepción aparente. Así, para cada carga, se aplica el algoritmo de optimación para determinar los factores de ajuste tDI y K(p) (si se conoce el tiempo de retraso a la ignición del motor base burn tDI, se fijará este valor). En este punto, los análisis efectuados han mostrado como los mejores resultados se obtienen definiendo la función objetivo únicamente como función del tiempo de extinción del motor base burn. XII.-55 A continuación, para cada carga y ángulo de tiro, se vuelve a aplicar el algoritmo de optimación para determinar los factores de ajuste f(iBB,MT) y fL definiendo aquí la función objetivo función del alcance y la deriva únicamente, ya que se ha determinado que así se obtienen los mejores resultados. Se va a esquematizar todo el proceso de resolución en los siguientes pasos: 1. Para cada carga, se aplica el algoritmo de optimación para determinar los factores de ajuste tDI y K(p) (si se conoce el tiempo de retraso a la ignición del motor base burn tDI, se fijará este valor). Obtenido el tiempo de extinción del motor base burn de experiencias realizadas, se calculan los factores de ajuste tDI y K(p) que mejor aproximen el tiempo de extinción del motor base burn. A través de experiencias realizadas se obtiene el tiempo de extinción del motor base burn, que se designará como tBi-exp para diferentes ángulos de tiro, definiendo la función objetivo como i=m E = ∑ t Bi − t Bi −exp i =1 t Bi −exp = t Bi −exp (MVi ) t Bi −exp = t Bi −exp (MVi , t DI , K ( p )) siendo los términos tBi los obtenidos como resultado de haber realizado todo el proceso de cálculo en función de los factores de ajuste tDI y K(p). En este paso se tomará f(iBB,MT) = 1 y fL = 1 . 2. A continuación, para cada carga y ángulo de tiro, se aplica el algoritmo de optimación para determinar los factores de ajuste f(iBB,MT) y fL. Obtenido el alcance y la deriva de experiencias realizadas, se calculan los XII.-56 factores de ajuste f(iBB,MT) y fL que mejor aproxime el alcance y la deriva. A través de experiencias realizadas se obtiene el alcance y la deriva, que se designaran como alcanceexp y derivaexp, para cada ángulo de tiro, definiendo la función objetivo como E = walc (alcancei − alcancei −exp ) + wder (derivai − derivai −exp ) 2 2 alcancei −exp = alcancei −exp (QEi , MVi ) derivai −exp = derivai −exp (QEi , MVi ) alcancei = alcancei (QEi , MVi , f (iBB , MT ), f L ) derivai = derivai (QEi , MVi , f (iBB , MT ), f L ) siendo los términos alcance y deriva los obtenidos como resultado de haber realizado todo el proceso de cálculo en función de los factores de ajuste f(iBB,MT) y fL. Los parámetros walc y wder son factores de peso que pueden ser diferentes dependiendo del paréntesis al que vayan a afectar y que permiten favorecer un requisito en perjuicio de los otros. Esto conlleva a que es necesario seleccionar el peso que afecta a cada uno de los paréntesis. Los factores de ajuste tDI y K(p) se han calculado en el apartado anterior para cada carga. Obtenido los factores de ajuste f(iBB,MT) y fL para cada ángulo de tiro, se pueden determinar los siguientes polinomios mediante un ajuste por mínimos cuadrados f (iBB , MT ) = iBB( MT =21) + b1 ⋅ (MT − 21) + b2 ⋅ (MT − 21) + b3 ⋅ (MT − 21) 2 iBB( MT = 21) = a0 + a1 ⋅ (QE ) + a2 ⋅ (QE ) + a3 ⋅ (QE ) 2 XII.-57 3 3 f L = a0 + a1 ⋅ (QE ) + a2 ⋅ (QE ) + a3 ⋅ (QE ) 2 3 Se ha considerado apropiado para este cálculo la siguiente condición de parada: L = 100, K1 = 4, ∈1 = 0.5 , K2 = 4, ∈2 = 0.5 y M = 1000. Se ha aplicado el anterior procedimiento para localizar los términos SC0, n y 4CxBB de los proyectiles de artillería OE155F2RTC FR HEA BB1 y OE155BONUS FR CBL/DUP BB1, y compararlos con los correspondientes coeficientes de reducción de resistencia base burn CxBB de dichos proyectiles. El proyectil OE155F2RTC FR HEA BB1 se ha disparado desde un tubo 155AUF1 y el proyectil OE155BONUS FR CBL/DUP BB1 desde un tubo 155AUF1, con tres cargas (la más alta, la más baja y una intermedia) y a cinco ángulos de tiro diferentes (300, 500, 700, 900 y 1200 milésimas). XII.-58 El siguiente punto es investigar hasta que punto el hecho de sustituir la curva original correspondiente al coeficiente de reducción de resistencia base burn CxBB por una curva parabólica por tramos similar, afectará al alcance y a la deriva. Disparando con todas las cargas y a cinco ángulos de tiro diferentes (300, 500, 700, 900 y 1200 milésimas), se ha obtenido una diferencia máxima de unos quince metros al intercambiar el coeficiente de reducción de resistencia base burn CxBB por la curva parabólica por tramos calculada. Estos valores suponen un error porcentual de alrededor de un 0.05 %. No obstante, no debe olvidarse que aún no se ha hecho uso de los factores de ajuste. Una vez ajustados los factores de ajuste para cada carga y ángulo de tiro se obtienen precisiones del orden de un metro, aplicando este método. Se ha aplicado el anterior procedimiento al proyectil de artillería M864 US CBL/DUP BB2, el cual había sido calculado originalmente por el método 2, y ahora se resolverá aplicando el método 1. El proyectil M864 US CBL/DUP BB2 se ha disparado desde un tubo M198 con tres cargas (la más alta, la más baja y una intermedia) y a cinco ángulos de tiro diferentes (300, 500, 700, 900 y 1200 milésimas). XII.-59 Inicialmente se define un parámetro de inyección de eficiencia óptima I0, teniendo en cuenta que al comienzo de la trayectoria el parámetro de inyección de propulsante en el motor base burn I debe ser mayor que el parámetro de inyección de eficiencia óptima I0, lo cual ocurre hasta que el propulsante del motor empieza a agotarse. En ese momento, el parámetro de inyección de propulsante en el motor base burn I se vuelve más pequeño que el parámetro de inyección de eficiencia óptima I0, y el coeficiente global de reducción de resistencia base burn cae rápidamente. Este comportamiento debe ocurrir para todas las cargas y todos los ángulos de tiro, por lo que el resto de parámetros necesarios para el cálculo deben estimarse de un proyectil base burn semejante, en esta primera etapa del cálculo. Como proyectil base burn semejante se ha considerado en esta fase el proyectil OE155F2RTC FR HEA BB1, considerando el área de combustión del propulsante base burn SC adimensional, n y el coeficiente de reducción de resistencia base burn CxBB. XII.-60 A continuación, se aplica el algoritmo de optimación para determinar los términos SC0, n y 4CxBB (CxBB0, CxBB1, CxBB2, y CxBB3). XII.-61 semejante seleccionado, podría ocurrir que se produjese este efecto, aunque como se verá a continuación los resultados han sido plenamente satisfactorios. Si por desgracia hubiese ocurrido esta circunstancia resultando en una caída abrupta al comienzo de la trayectoria, sería necesario volver a definir un nuevo parámetro de inyección de eficiencia óptima I0 y repetir todo el cálculo de nuevo. La habilidad para elegir el parámetro de inyección de eficiencia óptima I0 es crucial, dado que una mala elección conlleva un cálculo computacional muy importante. XII.-65 XII.-66 XII.-67 XII.-68 XII.-69 XII.-70 XII.-71 XII.-72 Se ha aplicado el procedimiento definido anteriormente al proyectil de artillería ER (Extended Range)-50/BB SP HEA BB1, empleando el método 1. Para aumentar las prestaciones en el campo de la artillería de 105 mm, EXPAL (Explosivos Alaveses S.A.), ha desarrollado este nuevo proyectil disparado con una carga L36, obteniendo alcances del orden de 20.1 kilómetros. La unidad de corrimiento baja es un conjunto atornillado en la base del proyectil actuando como generador del gas sin empuje. La ignición es activada por la carga propulsiva en el momento del disparo. El gas generado llena el volumen normal de la base donde hay turbulencia, aumentando el alcance hasta un 17%. Las bandas de forzamiento han sido adaptadas para poder disparar el proyectil desde un obús L118 Light Gun. Las experiencias fueron realizadas en el Centro de Ensayos Torregorda los días 21 y 22 de mayo de 2002. Se disparó con una carga L36, desde un obús L118 y a cinco ángulos de tiro diferentes (400, 525, 650, 750 y 850 milésimas). XII.-73 XII.-75 XII.-76 XII.-77 Inicialmente se define un parámetro de inyección de eficiencia óptima I0, teniendo en cuenta que al comienzo de la trayectoria el parámetro de inyección de propulsante en el motor base burn I debe ser mayor que el parámetro de inyección de eficiencia óptima I0, lo cual ocurre hasta que el propulsante del motor empieza a agotarse. En ese momento, el parámetro de inyección de propulsante en el motor base burn I se vuelve más pequeño que el parámetro de inyección de eficiencia óptima I0, y el coeficiente global de reducción de resistencia base burn cae rápidamente. Este comportamiento debe ocurrir para todas las cargas y todos los ángulos de tiro, por lo que el resto de parámetros necesarios para el cálculo deben estimarse de un proyectil base burn semejante, en esta primera etapa del cálculo. Como proyectil base burn semejante se ha considerado en esta fase el proyectil OE155F2RTC FR HEA BB1, considerando el área de combustión del propulsante base burn SC adimensional, n y el coeficiente de reducción de resistencia base burn CxBB. XII.-78 A continuación, se aplica el algoritmo de optimación para determinar los términos SC0, n y 4CxBB (CxBB0, CxBB1, CxBB2, y CxBB3). XII.-79 QM 1.0 Finalmente se debe comprobar que al comienzo de la trayectoria, el parámetro de inyección de propulsante en el motor base burn I no se vuelva más pequeño que el parámetro de inyección de eficiencia óptima I0 para alguna carga o algún ángulo de tiro, como se explicó en detalle anteriormente. Dado que inicialmente se había definido un parámetro de inyección de eficiencia óptima I0 basándose en el proyectil base burn semejante seleccionado, podría ocurrir que apareciese este efecto, aunque como se verá a continuación los resultados han sido enteramente satisfactorios. Si por desventura hubiese ocurrido esta circunstancia resultando en una caída abrupta al comienzo de la trayectoria, sería necesario volver a definir un nuevo parámetro de inyección de eficiencia óptima I0 y repetir todo el cálculo de nuevo. XII.-81 XII.-82 XII.-83 XII.-84 ESCUELA POLITÉCNICA SUPERIOR DEL EJÉRCITO TESIS DOCTORAL CAPÍTULO XIII MODELIZACIÓN DE PROYECTILES BASE BURN MÉTODO 2 MODELIZACIÓN DE PROYECTILES BASE BURN Director: Dr. Coronel Francisco Cucharero Pérez Doctorando: Comandante Fernando Aguirre Estévez XIII.-1 XIII.-2 XIII.-MODELIZACIÓN DE PROYECTILES BASE BURN MÉTODO 2 1. Índice 1. Índice 2. Proyectiles base burn método 2 2.1. Cambio en la presión de base adimensional para cambios en el parámetro de inyección de propulsante δBP δI 2.2. Derivada de referencia de la masa de propulsante del motor m& *f y tiempo de referencia de apagado del motor t B* 2.3. Duración del motor base burn tB - tDI, tiempo de referencia de * , presión de retraso a la ignición del motor base burn t DI referencia del aire atmosférico local Pr y velocidad de rotación axial de referencia pr 2.4. Implementación del método de cálculo XIII.-3 XIII.-4 2. Proyectiles base burn método 2 Se enumeran a continuación los coeficientes base burn para proyectiles base burn método 2. Estos coeficientes se han separado en dos grupos atendiendo a su dificultad, los primeros de ellos son los más simples y se reflejan a continuación. Los segundos merecen un estudio más exhaustivo. Símbolo IXB mDI mf tDI XCG0 XCGB Texto TDI Descripción Momento de inercia axial al finalizar el quemado del propulsante Masa de retraso a la ignición. Masa de propulsante del motor Tiempo de retraso a la ignición del motor base burn Distancia inicial del centro de masas a la ojiva Distancia del centro de masas a la ojiva, al finalizar el quemado del propulsante Texto Unidad IXB kgm2 MDI kg MFUEL kg TDI s XCG0 m XCGB m Dependencia MT, carga Nota Polinomio función de la temperatura de propulsante del motor, para cada carga a0 + a1 (MT - 21) + a2 (MT - 21)2 + a3 (MT - 21)3 para diferentes intervalos de MT (temperatura de propulsante del motor en °C), para una carga dada. Los intervalos de temperatura deben introducirse en orden ascendente con el formato TDI MTmin MTmax a0 a1 a2 a3 Si no se conocen los momentos de inercia, tanto el momento de inercia axial IX como el momento de inercia axial al finalizar el quemado del propulsante IXB, pueden estimarse a partir de proyectiles asistidos semejantes, como ya se estudió en análisis previos. El tiempo de retraso a la ignición del motor base burn tDI viene dado como una función de la temperatura de propulsante del motor para cada carga: XIII.-5 t DI = a0 + a1 ⋅ (MT − 21) + a2 ⋅ (MT − 21) + a3 ⋅ (MT − 21) 2 3 Es un valor pequeño por lo que si no se conociera, lo más simple es considerarlo nulo. El ajuste del resto de coeficientes absorberá la diferencia mínima que se produzca. La masa de retraso a la ignición mDI es la masa de las componentes del proyectil que se han desprendido hasta el tiempo de retraso a la ignición. Su valor es muy pequeño y puede considerarse nula en primera aproximación. El resto de coeficientes absorberá la pequeña diferencia que se produzca. Usualmente, la masa de propulsante del motor mf es dato pero si no lo fuera, podría estimarse dentro del proceso de optimización con el resto de coeficientes; no obstante, la masa de propulsante del motor depende del tiempo de extinción del motor base burn. Este dato resulta fácilmente medible mediante diversos dispositivos empleados en las experiencias. La distancia inicial del centro de masas a la ojiva XCG0 y la distancia del centro de masas a la ojiva al finalizar el quemado del propulsante XCGB, pueden estimarse a partir de proyectiles asistidos semejantes, como ya se estudió en análisis previos Los coeficientes base burn que se relatan a continuación son los más complejos de determinar y ya en si mismos constituyen un estudio de relevancia mayor, por lo que se examinaran cada uno por separado a fin de efectuar un estudio exhaustivo. Símbolo ∂BP ∂I m& *f Pr Descripción Cambio en la presión de base adimensional para cambios en el parámetro de inyección de propulsante Derivada de referencia de la masa de propulsante del motor Presión de referencia del aire atmosférico local XIII.-6 Texto Unidad - Dependencia Mach, parámetro de inyección MFDOTST kg/s pseudotiempo PRESREF Pa DBPDI Símbolo pr t B* , tBST tB - tDI * t DI , tDIST Texto DBPDI Descripción Velocidad de rotación axial de referencia Tiempo de referencia de apagado del motor Duración del motor base burn Tiempo de referencia de retraso a la ignición del motor base burn Texto SPINREF Unidad rad/s TBST s TBURN s TDIST s Dependencia MT, carga Nota Polinomio función del número de Mach a0 + a1M + a2M2 + a3M3 + a4M4 para diferentes intervalos del número de Mach y para un valor dado (b) del parámetro de inyección de propulsante en el motor base burn (I). Los números de Mach deben introducirse en orden ascendente con el formato MFDOTST DBPDI Mmin I b Mmax a0 a1 a2 a3 a4 Esta secuencia puede repetirse para diferentes valores del parámetro de inyección de propulsante que debe introducirse en orden ascendente. Polinomio función del pseudotiempo t* a0 + a1t* + a2t*2 para diferentes intervalos de t*. Los intervalos deben introducirse en orden ascendente con el formato TBURN MFDOTST t*min t*max a0 a1 a2 Polinomio función de la temperatura de propulsante del motor, para cada carga a0 + a1 (MT - 21) + a2 (MT - 21)2 + a3 (MT - 21)3 para diferentes intervalos de MT (temperatura de propulsante del motor en °C), para una carga dada. Los intervalos de temperatura deben introducirse en orden ascendente con el formato TBURN MTmin MTmax a0 a1 a2 a3 Para establecer correspondencia entre los resultados obtenidos computacionalmente y las experiencias observadas, se aplican los factores de ajuste. XIII.-7 Símbolo fBTP fBTp f(iΒΒ,MT) Texto Descripción Factor del tiempo de quemado de la presión del aire atmosférico Factor del tiempo de quemado de la velocidad de rotación axial Factor base burn Texto FBTPRES Dependencia carga FBTSPIN - carga FIBB - MT, QE, carga FBTPRES Nota Es una constante para cada carga, con el formato FBTSPIN FBTPRES a Es una constante para cada carga, con el formato FIBB Unidad - FBTSPIN a Polinomio función de la temperatura de propulsante del motor y del ángulo de tiro, para cada carga a0 + a1QE + a2QE2 + a3QE3 + b1(MT-21) + b2(MT-21)2 + b3(MT-21)3 para diferentes intervalos de MT (temperatura de propulsante del motor en °C) y QE (ángulo de tiro en milésimas), para una carga dada. Los intervalos de temperatura deben introducirse en orden ascendente con el formato FIBB MTmin MTmax a0 a1 a2 a3 b1 Parámetro Símbolo Factor del tiempo de quemado de la velocidad fBTp de rotación axial del proyectil base burn Factor del tiempo de quemado de la presión del fBTP aire atmosférico local del proyectil base burn Factor base burn f(iBB,MT) b2 b3 Unidad Límites típicos - − 0.5 < f BTp < −0.1 - − 0.9 < f BTP < −0.5 - 0.9