Resumen: El objetivo del presente trabajo es obtener y válidar un modelo y su aplicación para el cálculo de la posición de los controles y las caracteristicas aerodinámicas de cualquier helicóptero convencional en la condición de compensación o vuelo equilibrado, utilizando para ello métodos del análisis numérico. Obtener la posición de los controles, para el vuelo equilibrado en diferentes condiciones es parte de una etapa esencial durante el diseño general del helicóptero, sobretodo si están a disposición los datos de variación de caractristicas aerodinámicas y mecánicas en función del tiempo para las condiciones dadas. Será utilizado para el análisis aerodinámico el método cuasi estacionario de cantidad de movimiento – elemento de pala en su forma diferencial, utilizando el método de Runge-Kutta para resolver las ecuaciones diferenciales de batimiento y retraso.
Palabras clave: Diferencias finitas, Método de Newton, Análisis numérico, Aerodinamica de helicópteros, Runge-Kutta.
Abstract: The objective of this paper is to obtain a method for the calculation of the controls’ position and aerodynamics characteri-tics of conventional helicopters on the trim condition or steady flight, using numerical analysis methods. Obtain the controls’ position for steady flight in different conditions is essential during the general design of an helicopter, overall if there are avaible the aerodynamics and mechanics characteristics data as function of time in the gived condition. It will be use for the aerodynamics analysis the quasi steady blade element momentum theory in their differential form, using the Runge-Kutta method to resolve the flap and lag motion differential equations.
Keywords: Finite differences, Newton’s Method, Numerical analysis, Helicopters’ aerodynamics, Runge-Kutta.
1 Introducción
Este trabajo está orientado a la obtención de la posición de los controles y las caracteristicas aerodinamicas para un vuelo equilibrado en helicópteros convencionales. Un helicóptero de tipo convencional es aquella aeronave de alas giratorias que posee un único rotor principal encargado de generar sustentación y propulsión, cuyo tipo de control antitorque se emplea mediante un rotor que se encuentra situado en la parte trasera de la aeronave (Bramwell y col., 2001).
Los helicópteros son actualmente empleados para distintas tareas de transporte, y durante sus primeras etapas del de diseño preliminar se estima que deben responder adecuadamente a las condiciones establecidas para un vuelo estacionario, crucero, ascendente y descendente. Al momento de realizar una alteración mayor a un helicóptero convencional y modificar su objetivo principal de diseño, es posible que los limites máximos y minimos de los controles sean inadecuados para la nueva tarea a realizar.
El control de un helicóptero convencional es obtenido mediante el accionamiento de los mecanismos que cambien la posición del ángulo de paso colectivo del rotor de cola, el ángulo de paso colectivo del rotor principal, el ángulo de paso cíclico lateral y longitudinal, posicionar el helicóptero al ángulo de ataque y de inclinación lateral adecuado. En las diferentes etapas de vuelo, las condiciones de velocidades relativas en el rotor principal de la aeronave se ven modificadas, desde una condición simétrica, para un vuelo estacionario o axial, hasta una condición de asimetría de velocidades durante un vuelo con componente horizontal.
La idea principal es que cualquier helicóptero convencional, tanto experimental como comercial, pueda ser estudiado a través del método que se propone desarrollar en el trabajo. Esto conlleva a asumir que pueden generarse geometrías de palas poco convencionales, en donde las aproximaciones realizadas por diferentes autores en el campo pueden tener un margen de error considerable.
El método que comúnmente se emplea en la industria de los helicópteros durante su análisis preliminar es el de la teoría de cantidad de movimiento – elemento de pala (Bramwell y col., 2001). Cada una de ellas puede ser utilizada de manera independiente para realizar el análisis, como Cruz, Oropeza y Vigueras, quienes emplearon la teoría del elemento de pala para realizar el análisis de un helicóptero experimental. La combinación de estos métodos para el análisis de un helicóptero UH-60 fue llevado a cabo por Sevinç Çalişkan en su trabajo de maestría, a través de la implementación en un código para realizar cálculos aerodinámicos, que serían utilizados para estimar la condición de trimado en vuelo horizontal de la aeronave mediante el desarrollo analítico de la teoría.
La solución analítica de las ecuaciones es muy compleja de obtener, por lo que será empleado el análisis numérico. En general, el modelo se trata de un sistema de ecuaciones diferenciales acopladas, y la idea es obtener un resultado convergente de las variables de interés en función del cambio de posición de las palas en el acimut del rotor. Será utilizado un método de integración numérica, el método de las diferencias finitas y el método de Runge – Kutta (4RK) para hallar la solución al sistema.
Por otra parte, para la obtención de la posición de los controles, será necesario resolver un sistema de seis ecuaciones con seis incógnitas, en donde se desconoce la forma analítica de las variables en cuestión, por lo tanto, será aplicado el método de Newton para encontrar los valores correspondientes de estas variables que satisfagan el sistema de ecuaciones.
2 Marco Teórico
2.1 Extrapolación de Richardson
Está fundamentada en la expansión de la serie de Taylor. La ecuación de la extrapolación de Richardson está definida (Yang Cao y col., 2005):
(1)Donde n define el orden de truncamiento. Para 𝐧=𝟏 el orden de truncamiento es 4, y la ecuación correspondiente para la primera derivada es
(2)Esta derivada numérica será utilizada para el sistema en derivadas parciales del método de Newton.
2.2 Integración por la regla de Simpson
Con una expansión de tercer grado del polinomio de Taylor se obtiene la expresión expandida de la regla de Simpson, para un grado de truncamiento de 4. La regla de Simpson implica que (Burden & Faires, 2011):
(3)Para un valor de 𝛏 en el intervalo (𝐱𝟎, x𝟐). Este método será utilizado para resolver las integrales del modelo.
2.3 Método de Newton
Si se tiene un sistema de ecuaciones de 6 funciones y 6 incógnitas, este sistema puede escribirse de la siguiente forma
Donde [𝐉] es la matriz Jacobiana, por lo que la so-lución al vector [𝐗] sería
Donde
Entonces, la solución final para las incógnitas sería a través de
(4)La solución a las ecuaciones de equilibrio estatico del helicoptero será obtenida mediante este método.
2.4 Método Runge-Kutta de cuarto orden
A través de la expansión de la serie de Taylor para dos variables, se obtiene un método para la resolución de ecuaciones diferenciales con un grado de truncamiento de 4. Para una ecuación diferencial ordinaria de primer grado, se tiene que la solución puede ser obtenida mediante las siguientes ecuaciones (Burden y col., 2011):
(5)
(6)
(7)
(8)
(9)
(10)Donde 𝐰𝐢+𝟏 es la solución de la variable, 𝐡 es el paso de tiempo, 𝐭𝐢 es el tiempo discreto y 𝛂 el valor inicial. Este método será aplicado a las ecuaciones de batimiento y retraso.
3 El Modelo
Por el modelo de teoría de cantidad de movimiento axial propuesto por Rankine & Froude (Durand 1943), en el cual se asume un volumen de control donde un fluido está pasando por un disco poroso, como se muestra en la Fig. 1. Entonces la tracción generada se define como.
Fig. 1.
Volumen de control
Donde 𝐦̇ es el flujo másico y Δ𝐕 el cambio de velocidad, por lo tanto
(11)Siendo 𝛒 la densidad del fluido y 𝐀 el área del disco. Aplicando el teorema de Bernoulli al volumen de control establecido se obtiene
Obedeciendo la ley de conservación de la masa y cantidad de movimiento, podemos multiplicar por el área del disco para obtener
(12)Asumiendo que la velocidad axial en el disco 𝐕𝐝=𝐕∞𝐬𝐢𝐧𝛂+, e igualando (11) y (12) para despejar 𝐕𝐬 se obtiene
Sustituyendo esta última expresión en (11), quedaría
(13)Donde 𝐯 es la velocidad inducida en el disco. El área del disco es igual al área dentro de su circunferencia 𝐀=𝛑𝐑𝟐. Si se sustituye el área del disco en la de tracción y se deriva en función del radio
(14)Sin embargo, al considerar un movimiento de traslación del disco del rotor, se puede observar que la velocidad verdadera en dicho disco es igual a la resultante de las componentes, es decir
Sustituyendo en la Ec. (14) y despejando 𝐯 se tiene
(15)Haciendo un cambio de variable para ser usado como auxiliar
Entonces
(16)Con la aplicación del método Newton – Raphson
(17)Donde
La convergencia de 𝐯 es obtenida tras cuatro o cinco iteraciones, dependiendo del orden de precisión deseado. La Ec. (13) también fue deducida por Glauert al resolver el problema del arrastre inducido, considerando un ala circular, evaluada por la línea de sustentación de Prandtl asumiendo una forma de distribución de sustentación elíptica.
Mangler y Squire, al igual que Coleman, Feingold y Stempin (Johnson 1980), tras refinar el modelo de cantidad de movimiento utilizando la teoría de vórtices, concluyen que la expresión adecuada para la distribución de velocidad inducida sobre el rotor es
(18)Donde 𝐯𝟎 es la velocidad inducida obtenida por la Ec. 17, 𝐊𝐱 es un factor de corrección, 𝛚 es la velocidad angular del rotor y 𝐭𝐢=𝐭𝐢−𝟏+𝐝𝐭, siendo 𝐭 el tiempo. A su vez
Y para 𝐕∞𝐜𝐨𝐬𝛂≈, entonces 𝐊𝐱=𝟎
Para obtener las fuerzas y momentos generados por el rotor, se utilizará el modelo del elemento de pala, cuyo fundamento está basado en la Figura 2. Además se utilizará la corrección del factor de pérdida en punta de pala de Prandtl.
Fig. 2.
Elemento de Pala
De la Fig. 2 se puede deducir que
Donde 𝐍 es el número de palas y 𝐅 es el factor de corrección de Prandtl. Basados en la Fig 2, el diferencial de tracción en función de un diferencial de 𝐫 sería
(19)Donde 𝐂𝐋 es el coeficiente de sustentación. A su vez, para la ecuación del torque en función del radio se tiene
(20)Podemos decir que la velocidad axial total relativa al elemento de pala es la suma algebraica de todas las velocidades axiales en dicho elemento, como se ilustra en la Fig. 3.
(21)
Fig. 3.
Velocidades relativas
Donde 𝛃 es el ángulo de batimiento y 𝛃̇ representa la velocidad angular de batimiento. Por otra parte, la velocidad tangencial al elemento de pala será la suma algebraica de todas las velocidades aerodinámicas en esta dirección.
(22)
Fig. 4.
Plano del rotor
Donde 𝛇 es el angulo de retraso y 𝛇̇ es la velocidad angular de retraso. Entonces la velocidad total 𝐖 en el elemento de pala sería
(23)El ángulo de paso efectivo 𝛗, y el ángulo de ataque del elemento de pala 𝛂, estarían definidos
(24)
(25)A partir de este modelo se pueden definir las fuerzas y los momentos de manera cuasi estacionaria.
(26)
(27)
(28)
(29)Y en su forma integral, las fuerzas y momentos en función del tiempo para cada pala
(30)
(31)
(32)
(33)
(34)Donde 𝐫𝐡 es el radio del cubo del rotor. Para evaluar el desempeño estático de la aeronave bajo las condiciones deseadas se necesitan los datos del modelo cuasi estacionario hasta un mínimo de 𝚿𝐢=𝐭𝐢𝛚=𝟐𝛑 (Tomando como referencia lo mostrado en Fig. 4), y de esa manera obtener
(35)
(36)
(37)
(38)Por otro lado, los momentos son obtenidos a través del sistema mostrado en la Fig. 5
Fig. 5.
Sistema para momentos en el rotor
(39)
(40)
(41)
(42)Y sus coeficientes son (Leishman, 2006)
(43)
(44)
(45)
(46)
(47)
(48)
(49)Se definirán ahora las ecuaciones que rigen el movimiento de retraso, batimiento y el control de paso. El ángulo de paso es 𝛉, y se define (Johnson 1980)
(50)Donde 𝛉𝟎 es el paso colectivo, 𝛉𝟏𝐬 el paso cíclico longitudinal, 𝛉𝟏𝐜 es paso cíclico lateral y Δ𝛉 el factor del acopla-miento paso-batimiento, cuya ecuación es (Johnson W, 1980)
(51)
Fig. 6.
Sistema de batimiento (Vista lateral del rotor)
Donde 𝛅𝟑 es el ángulo formado por el eje de la articulación de batimiento y el link de paso. A su vez, el batimiento está definido por el sistema mostrado en la Fig. 6 y la ecuación de movimiento correspondiente está expresada en una ecuación diferencial, como sigue
(52)La primera integral del lado derecho de la igualdad corresponde a la contribución aerodinámica, la segunda a la contribución por fuerza centrífuga, el tercer factor es el momento restaurador del muelle, y el último es la contribución de la fuerza de Coriolis. A su vez, la ecuación diferencial que rige el movimiento de retraso es
(53)
Las contribuciones en este caso son las mismas que para el movimiento de batimiento, agregando la contribución del amortiguador. Llegado a este punto, es importante definir ciertas funciones para el análisis de los parámetros 𝛃 y 𝛇 en función del acimut del rotor desde una perspectiva estacionaria. Por las series de armónicas de Fourier se puede definir el ángulo de batimiento y retraso en función del acimut como (Leishman 2006):
(54)
(55)Las armónicas 𝛃𝟎, 𝛇𝟎, 𝛃𝐧𝐜, 𝛇𝐧𝐜, 𝛃𝐧𝐬, 𝛇𝐧𝐬 son obtenidas
(56)
(57)
(58)Donde 𝛃 y 𝛇 se obtienen en función del tiempo al resolver las ecuaciones diferenciales. Como se observa, se tiene un sistema de ecuaciones diferenciales acopladas. El objetivo de la simulación es obtener la respuesta dinámica del sistema para unos valores iniciales establecidos.
Fig. 7.
Sistema de retraso (Vista superior del rotor)
Se considerará un sistema de fuerzas en equilibrio, cuyas incógnitas serán la posición de los controles, el ángulo de ataque del helicóptero y el ángulo de inclinación lateral del helicóptero. De los sistemas mostrados en las Fig. 8 y 9 se deduce que:
Fig. 8.
Fuerzas y momentos vistos desde la vista lateral
Fig. 9.
Fuerzas y momentos vistos desde la vista posterior
La ecuación de equilibrio de fuerzas en el eje vertical
(59)La ecuación de equilibrio de fuerzas en el eje longitudinal
(60)La ecuación de equilibrio de fuerzas en el eje lateral
(61)La ecuación de equilibrio de momento de cabeceo
(62)La ecuación de equilibrio de momento de alabeo
(63)La ecuación de equilibrio de momento de guiñada
(64)Las incógnitas son 𝛉𝟎, 𝛉𝟏𝐬, 𝛉𝟏𝐜, 𝛉𝟎𝐓𝐑, 𝛂𝐇 y 𝛟𝐅. Será utilizado el método de Newton para resolver el sistema de ecuaciones, con derivadas parciales halladas por diferencias finitas.
3.1 Cálculo de la posición de los controles para un helicóptero convencional
El Bo105 S123 es un helicóptero bimotor de 2.5 toneladas, preparado para diversas misiones en transporte, operaciones costa fuera, policía y operaciones de campo de batalla. Es un helicoptero de cuatro palas fabricadas en materiales compuestos de refuerzo de fibra, cuyo rotor no posee articulaciones.
Los datos de campo de trimado, junto con los calculados por Helisim, están disponibles en la literatura (Padfield, 2007) y serán utilizados como referencia y comparación para el método desarrollado. Los valores iniciales
La geometría y otros datos del helicóptero son
Fig. 10.
Bo105 S123
Caracteristicas del Bo105
Caracteristicas del Bo105
Caracteristicas del Bo105
Donde Clα es la pendiente de sustentación del perfil en la pala, rh es el radio del cubo del rotor, δ3 es el ángulo de acoplamiento de paso-batimiento, αE es el ángulo de inclinación del eje del rotor, h es la distancia vertical al cubo del rotor, x es la distancia longitudinal al eje del rotor, A es el área de referencia del componente, αEH es el ángulo de incidencia del estabilizador horizontal, βEV es el ángulo de incidencia del estabilizador vertical.
Los sufijos TR, EV, EH y CG hacen referencia a que las variables pertenecen al rotor de cola, estabilizador vertical, estabilizador horizontal y centro de gravedad correspondientemente. Para el cálculo fue utilizada la herramienta Scilab. Los datos del fuselaje pueden ser encontra-Los resultados se muestran de la Fig. 11 a la 16.
Fig. 11.
Torque contra Velocidad
Fig. 12.
Paso Colectivo de rotor principal contra Velocidad
Fig. 13.
Ángulo de ataque contra Velocidad
Fig. 14.
Paso Colectivo de rotor de cola contra Velocidad
Fig. 15.
Paso Ciclico longitudinal contra Velocidad
Fig. 16.
Paso Ciclico lateral contra Velocidad
En las graficas obtenidas se muestra la comparación de datos de campo del Bo105, datos obtenidos por Padfield mediante Helisim y los datos obtenidos mediante el método y modelo del presente artículo. El comportamiento de las todas las graficas tienen bastante similitud. El torque y el paso colectivo del rotor principal muestran una curva convexa con el minimo entre los 25 y 35 metros por segundo.
El paso colectivo del rotor de cola fue sobreestimado por el presente modelo, sin embargo presenta un comportamiento semejante a los demás datos, siendo una curva convexa. El paso cíclico lateral muestra una curva con un punto de inflexión luego de cierta velocidad, y un máxmo pronunciado para el caso de los datos de la prueba de vuelo. El ángulo de ataque del helicóptero se ve representada por una curva decreciente. Y finalmente el paso ciclico longitudinal muestra una curva con umportamiento igualmente decreciente.
4 Cálculo de las características aerodinámicas del helicóptero convencional
Fueron calculadas las características aerodinámicas del Bo105 bajo su condición de trimado para un vuelo estacionario, como se muestran en las Fig. 17 a 19 y un vuelo a 45 metros por segundo, ilustrado en las Fig. 20 a 22. Se destaca que en las figuras mostradas, la dirección del vuelo del helicóptero es en dirección al eje negativo de las abscisas.
La Fig. 17 muestra una distribución aparentemente simétrica de velocidad inducida en todo el rotor, a su vez, dicha velocidad aumenta radialmene hasta llegar a la punta de pala. En la Fig. 18 el coeficiente de sustentación (que está asociado directamente con los momentos en el cubo del rotor) tiene una concentración de este a estribor del rotor, momentos antes de la posiciòn donde el análisis realizado muestra el máximo de la función de batimiento.
Fig. 17.
Velocidad inducida en el rotor. V = 0 m/s
Fig. 18.
Coeficiente de sustentación en el rotor. V = 0 m/s
Si se eleva la velocidad de traslación de la aeronave, se percibe a través de la Fig. 19 un cambio dramático en la velocidad inducida del rotor. Lo que podría denominarse como el borde de ataque del disco, es donde se presenta el minimo de la velocidad inducida, mientras que en el borde de fuga se encuentra el máximo.
La Fig. 20 por su parte, deja en clara evidencia una zona cercana al cubo, con forma aproximadamente cincunferencial, en donde el coeficiente de sustentación asume muchos valores, tanto superiores como inferiores a la escala mostrada. Alrededor de esta zona, muestra una distribución parecida a la de la Fig. 18, solo que focalizada a bavor y hacia el borde de ataque del rotor. De forma contraria al caso en vuelo estacionario, el máximo en el ángulo de batimiento se presenta en dirección al morro de la aeronave, mientras que el minimo está en dirección a la cola.
Fig. 19.
Velocidad inducida en el rotor. V = 45 m/s
Fig. 20.
Coeficiente de sustentación en el rotor. V = 45 m/s
5 Conclusiones
El comportamiento de las Fig. 11 y 12 claramente representan las diferentes contribuciones del arrastre inducido y el arrastre parasito en la aeronave. El arrastre inducido va decreciendo conforme aumenta la velocidad, caso contrario el arrastre parasito, el cual crece exponencialmente en función de la velocidad, lo que hace que se requiera un paso colectivo superior, y por lo tanto se genere un torque mayor. Estos controles fueron estimados por el modelo con una precisión aceptable, siendo este un primer nivel de modela-do.
La Fig. 13 revela que fue subestimado el ángulo de ataque del helicóptero para altas velocidades, sin embargo, el comportamiento presentado está conforme a la realidad. El helicóptero para el caso estacionario requerirá un angulo de ataque muy cercano al de inclinación del rotor, a medida que aumenta la velocidad de traslación necesita que este angulo sea creciente puesto que el rotor debe proveer la propulsión para estar en equilibrio con la fuerza de arrastre, además del peso de la aeronave.
En la Fig. 14 se muestra una clara sobreestimación del paso de colectivo del rotor de cola. Las causas son las simplificaciones en modelo de rotor de cola que fue utilizado para estimar los datos que posteriormente se emplearían en las ecuaciones de equilibrio de la aeronave. Bastaría con hacer unas consideraciones adicionales en el modelo de rotor de cola para corregir esto. Con ello, muestra la importancia de considerar un modelo igual de complejo para el rotor de cola, donde las velocidades relativas y el número de elementos discretos sean tomados en cuenta correctamente.
El paso cíclico longitudinal y lateral mostrado en las Fig 15 y 16 respectivamente, presentan una apreciación aceptable con un comportamiento adecuado a las condiciones que presenta la aeroanve en la realidad. El paso cíclico longitudinal está directamente relacionado con el momento de cabeceo, y el fuselaje del helicóptero al producir dicho momento que tiende a bajar el morro de la aeronave, este necesita compensar con un momento opuesto.
Por otra parte, el paso cíclico lateral se relaciona al momento de alabeo, y se acciona para compensar la contribución del paso cíclico longitudinal al momento de alabeo y el momento producido por el rotor de cola, esto mientras no existan vientos laterales.
La velocidad inducida en el rotor, mostrado en la Fig. 17 presenta ese comportamiento puesto a que no existe velocidad de traslación, por lo que la distribución es aparentemente simétrica, lo cual simboliza una buena aproximación al caso real. La Fig. 18 revela una focalización del coeficiente de sustentación en la zona mencionada ya que se necesita crear un momento de alabeo y cabeceo para contrarrestar los momentos generados por la excentricidad del centro de gravedad respecto al cubo del rotor, además del momento de alabeo del rotor de cola.
Para una velocidad superior, es decir, 45 metros por segundo, la velocidad inducida cambia por completo su distribución, estando conforme a las ecuaciones desarrolladas por Feingold y Stempin, en donde la velocidad decae en el borde de ataque del disco del rotor. Por tratarse de una superficie en movimiento relativo al aire que se encuentra delante, se crea una zona parecida un punto de estancamiento, tal y como sucede en cualquier solido que enfrenta un flujo. Ademàs de esto, el angulo formado entre el disco del rotor y la dirección del flujo, provocan este efecto en la velocidad inducida.
Finalmente, en la Fig. 20 la zona de flujo inverso está claramente detallada por el exceso de contornos cercanos al cubo del rotor, en donde los valores de coeficiente de sus-tentación alcanzados tienen un comportamiento caótico.
A su vez, la distribución de coeficiente de sustentación en el resto del rotor se encuentra conforme a las fuerzas y momentos que necesita compensar a la aeronave para que se mantenga en equilibrio a altas velocidades (donde el fuselaje posee una contribución considerable, tanto en momento de cabeceo, como en arrastre), considerando el efecto de asimetría de velocidades (la cual contribuye de manera superior en la ecuación de la sustentación).
Referencias
Bramwell ARS, Done George, Balmford David, 2001, Bramwell’s Helicopter Dynamics Second edition. Butter-worth-Heinemann, Inglaterra, Oxford.
Burden Richard, Faires JD, 2011, Numerical Analysis, Brooks/Cole, EUA, Boston.
Cruz F, Oropeza R, Vigueras R, 2011, Análisis aerodinámico de una pala del rotor principal de un helicóptero experimental. Trabajo de Grado, Instituto Politécnico Nacional de Mexico, Mexico, DF.
Durand William, 1943, Aerodynamic Theory, California Institute of Technology, EUA, California.
Yang Won, Cao Wenwu, Chung Tae-Sang, Morris John, 2005, Applied Numerical Methods Using MATLAB®, John Wiley & Sons Inc, EUA, Nueva Jersey.
Johnson Wayne, 1980, Helicopter Theory, Princeton Uni-versity Press, EUA, Nueva Jersey.
Leishman JG, 2006, Principles of Helicopter Aerodynamics, Pess Syndicate of the University of Cambridge, Reino Unido, Cambridge.
Padfield GD, 2007, Helicopter Flight Dynamics, The Theory And Application of Flying Qualities and Simulation Modelling, Second Edition, Blackwell Publishing, Reino Unido, Oxford.
Sevinç Çalişkan, 2009, Development of Forward Flight Trim and Longitudinal Dynamic Stability Codes and Their Application to a UH-60 Helicopter, Tesis de Maestría, Universidad Técnica de Medio Oriente, Turkia, Ankara.
Notas de autor
