Estudio analítico y numérico del scattering de una onda electromagnética plana por un cilindro conductor ideal

Estudio analítico y numérico del scattering de una onda electromagnética plana por un cilindro conductor ideal Héctor Bravo Muñoz Supervisado por Dr. Rafael Rodriguez Boix Grado en Física, Universidad de Sevilla 5 de septiembre de 2022 Índice 1. Introducción 1 1.1. Justificación y antecedentes . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2. Breve descripción del trabajo . . . . . . . . . . . . . . . . . . . . . . . . . 2 2. Ecuación integral para la densidad de corriente sobre la superficie del cilindro (caso TMz). Resolución de la ecuación integral mediante el Método de los Momentos 4 2.1. Expresiones del campo eléctrico en términos de funciones de Green y de la densidad de corriente en los casos tridimensional y bidimensional . . . . . . 4 2.2. Ecuación integral para la densidad de corriente en el problema de scattering bidimensional .................................. 11 2.3. Resolución de la ecuación integral mediante el Método de los Momentos con funciones base rectangulares . . . . . . . . . . . . . . . . . . . . . . . . 14 3. Resultados numéricos 22 3.1. Validación para el caso particular del cilindro de sección circular (sección radar bidimensional y densidad de corriente) . . . . . . . . . . . . . . . . . 22 3.2. Análisis de convergencia del MdM . . . . . . . . . . . . . . . . . . . . . . . 26 3.3. Resultados para cilindros con otras secciones transversales (elíptica, ojival yrectangular).................................. 28 4. Conclusiones 43 Agradecimietos Me gustaría agradecer en primer lugar el trabajo de mi tutor Rafael Rodriguez Boix, ya que a pesar de las dificultades de dirigir un TFG a distancia, decidió contar igualmente conmigo para ello. Agradecerle el tiempo dedicado en la redacción de sus apuntes que explicaban toda la teoría tras este trabajo, y el tiempo dedicado a corregirlo de manera exigente. Sin duda, más allá de haberme enseñado electromagnetismo, me ha enseñado cómo ser un buen profesor. Como él diría: “Me gusta hacer las cosas bien”. También quiero dar las gracias a mis compañeros de clase Alex y Pablo, por haberme acompañado durante estos 4 años y haberme ayudado tanto. Por supuesto agradecer también a Fernando por su compañía durante los 3 primeros años. Agradecer también a mis amigos de Montequinto: Abraham, Cristina, Andrea, Álvaro, Mayo y Marta, ya que aunque las cosas se pusieran difíciles, ellos han estado siempre apoyándome. Dar las gracias por supuesto a mis compañeros de ERASMUS, mi segunda familia: Maribel, Juan, Albert, Boni, Emma, María, Laura, Miriam, Ana, Emmanuelle, Martos y Alex, por acompañarme este último año en Alemania. Gracias a ellos, he aprendido mucho más allá de la física. Y por último y más importante, dar las gracias a mis padres Paloma y Alberto, y a mi hermano Alberto. A pesar de que estos últimos años han sido especialmente complicados por muchas razones, ellos lo han dado todo por mí. A ellos se lo debo todo. 1. Introducción 1.1. Justificación y antecedentes Este trabajo se centra en el estudio de los campos electromagnéticos de scattering generados al hacer incidir una onda electromagnética plana sobre un cilindro conductor perfecto -esto es, sobre un cilindro conductor de conductividad infinitade longitud infinita y sección transversal arbitraria. En el caso concreto estudiado en el trabajo, la onda incidente tendrá una polarización del tipo TMz, lo cual significa que el campo eléctrico de la onda incidente es paralelo al eje del cilindro [1]. Adicionalmente, también tiene interés el estudio del scattering en el caso en que la onda tiene polarización TEz (esto es, cuando el campo eléctrico es perpendicular al eje del cilindro), pero el estudio completo de los dos casos es demasiado amplio para ser incluido en un Trabajo Fin de Grado [1]. El considerar infinito el cilindro objeto de estudio permite reducir el problema del scattering a 2 dimensiones. El estudio del scattering por objetos tridimensionales tiene desde luego mayor interés práctico, pero el caso tridimensional es muchísimo más complejo que el bidimensional, y la resolución de ese problema debería abordarse en el marco de un doctorado, y por tanto, queda fuera del ámbito de los Trabajos Fin de Grado. A lo largo de este Trabajo Fin de Grado, llevaremos a cabo cálculos de la sección radar bidimensional [1] de cilindros de sección arbitraria (en el caso tridimensional, la sección radar es la figura de mérito que caracteriza a los objetos cuando sobre ellos incide una onda electromagnética, y es el parámetro relevante del objeto que interviene en la ecuación de radar [2]) y de las corrientes generadas en el conductor. Con vistas al cálculo de la sección radar y de las corrientes en los cilindros, resolveremos numéricamente una ecuación integral de campo eléctrico mediante el Método de los Momentos. Para verificar que el método numérico se ha implementado correctamente, se llevarán a cabo estudios de convergencia de los distintos parámetros involucrados. En un principio, se llevará a cabo el estudio exhaustivo del scattering por un cilindro de sección circular ya que para este caso existe solución analítica. A continuación, se estudiarán cilindros de secciones elípticas, ojivales y rectangulares. Dado que existen casos límites en los que una elipse y una ojiva pasan a ser circunferencias, se estudiarán esos casos límite para validar los resultados obtenidos para cilindros de sección elíptica y ojival. El cálculo de los campos de scattering que se producen cuando una onda electromag- 1 nética incide sobre un objeto es un problema clásico de la Física. Este problema tiene interés práctico en conexión con el funcionamiento de los radares. Estos sistemas se utilizan para medir la distancia a un objeto, para medir la velocidad del objeto si éste está en movimiento (radar de efecto Doppler), y finalmente, para estimar el tamaño del objeto a través de su sección radar (lo que en Física se suele llamar “sección eficaz de dispersión”). Los radares tienen aplicaciones civiles (tráfico en aeropuertos, meteorología, alarmas para ladrones, radares de velocidad en carreteras, etc.), militares (navegación aérea y marina o detección de misiles), y por supuesto, aplicaciones científicas (astronomía, radiodetección de recursos naturales, etc.). El conocimiento de este campo ha sido determinante en algunos momentos de la historia ya que aporta una gran ventaja militar (así ocurrió, por ejemplo, a finales de la Segunda Guerra Mundial). Un ejemplo muy conocido es el del caza F-117, que fue el primer avión militar de la historia con tecnología stealth [3](tecnología de invisibilidad en castellano). En la tabla 12.3 del libro de David M. Pozar Microwave Engineering [2] podemos encontrar valores típicos de secciones radar. Los valores usuales de la sección radar para un caza están entre 3 y 8 m2. Sin embargo, debido a su geometría y a los materiales con los que está construido, la sección radar del F-117 es del orden de 0.003 m2, que es la sección radar de un pájaro del tamaño de un colibrí. Pues bien, este hecho aportó una considerable ventaja al bando estadounidense durante la Guerra del Golfo, y fue crucial para que los estadounidenses acabaran ganando la guerra. 1.2. Breve descripción del trabajo Con vistas a obtener los campos de scattering, vamos a resolver la ecuación integral que se obtiene al imponer que se anule el campo eléctrico tangencial total (suma del incidente y del de scattering) en la superficie de los cilindros. A esta ecuación integral se la suele llamar ecuación integral de campo eléctrico (electric field integral equation o EFIE en inglés). Para deducirla, partiremos de las ecuaciones de Maxwell, expresaremos los campos en función de los potenciales, y haciendo uso de funciones de Green para los potenciales y de la ecuación de continuidad, llegaremos a una ecuación integral que tiene como incógnita la densidad de corriente superficial inducida a los cilindros. Dado que esta ecuación integral se obtiene haciendo uso del potencial vector magnético y del potencial escalar eléctrico, también se la llama a veces ecuación integral de potenciales mixtos. Dicha ecuación integral la vamos a resolver mediante el Método de los Momentos. 2 El Método de los Momentos se basa en aproximar la densidad de corriente inducida en el cilindro como una combinación lineal de funciones base conocidas (en nuestro caso las funciones base serán pulsos rectangulares). Los coeficientes indeterminados de esta combinación lineal se obtienen proyectando la expresión resultante de sustituir la aproximación de la densidad de corriente en la ecuación integral sobre las mismas funciones base (a esta implementación del Método de los Momentos se la conoce como Método de Galerkin), lo cual lleva a que los coeficientes indeterminados sean la solución de un sistema de ecuaciones lineales. Los elementos de la matriz de coeficientes de este sistema de ecuaciones son un conjunto de integrales dobles que es preciso calcular numéricamente. Dado que los integrandos de estas integrales poseen singuaridades logarítmicas, hay que utilizar un método específico de integración numérica capaz de tratar con integrandos con singularidades logarítmicas. En este trabajo se ha utilizado el método de integración numérica basado en las cuadraturas de Ma-Rokhlin-Wandzura [4]. El cálculo numérico de las integrales, la resolución del sistema de ecuaciones para obtener la densidad de corriente de los cilindros y el cálculo adicional de la sección radar bidimensional se llevan a cabo mediante códigos en lenguaje MATLAB. Las magnitudes calculadas son validadas mediante comparación con resultados previamente publicados o mediante comparación con resultados analíticos (concretamente, con los resultados que existen para un cilindro de sección circular). Asimismo, se llevan a cabo estudios de convergencia para chequear el buen funcionamiento del Método de los Momentos. 3 2. Ecuación integral para la densidad de corriente sobre la superficie del cilindro (caso TMz). Resolución de la ecuación integral mediante el Método de los Momentos 2.1. Expresiones del campo eléctrico en términos de funciones de Green y de la densidad de corriente en los casos tridimensional y bidimensional De ahora en adelante, supondremos que todas las magnitudes electromagnéticas escalares y vectoriales de este trabajo dependen del tiempo sinusoidalmente con una frecuencia angular ωde forma que en las ecuaciones que vamos a manejar dichas magnitudes puedan ser sustituidas por sus fasores complejos. En el caso de que la dependencia con el tiempo fuera arbitraria, dichos fasores representarían las transformadas de Fourier de las magnitudes, y para recuperar la dependencia con el tiempo de las correspondientes magnitudes, habría que calcular la transformada inversa de Fourier. Consideremos ahora un volumen τ′en el vacío, en el que existen corrientes y cargas (fuentes de campo electromagnético) variables sinusoidalmente en el tiempo (régimen armónico) que llevan asociadas un fasor densidad de corriente  J( r′)y un fasor densidad de carga ρ( r′)(véase la Figura 1). Debido a la ecuación de continuidad, se tiene que cumplir que:  ∇′· J( r′) + jωρ( r′)=0 ⇒ρ( r′) = j ω ∇′· J( r′).(1) donde jes la unidad imaginaria. Las ecuaciones de Maxwell para los fasores del campo eléctrico  E(r)y del campo magnético  B(r)creados por las fuentes son:  ∇·  E(r) = ρ(r) ε0 (2)  ∇×  E(r) = −jω  B(r)(3)  ∇·  B(r) = 0 (4)  ∇×  B(r) = µ0 J(r) + jωε0µ0 E(r).(5) 4 Figura 1: Volumen arbitrario con densidad de carga y densidad de corriente. Las ecuaciones que relacionan los fasores de los campos con los fasores del potencial vector magnético  A(r)y del potencial escalar eléctrico ϕ(r)son:  B(r) =  ∇×  A(r)(6)  E(r) = −jω  A(r)− ∇ϕ(r).(7) Si hacemos uso del gauge de Lorentz, se debe cumplir que:  ∇·  A(r) + jωε0µ0ϕ(r)=0 ⇒ϕ(r) = j ωε0µ0  ∇·  A(r).(8) Sustituyendo (6) y (7) en (5), se llega a que:  ∇×( ∇×  A) = µ0 J+jωε0µ0·(−jω  A− ∇ϕ)⇒ ⇒ ∇( ∇·  A)− ∇2 A=µ0 J+k2 0 A−jωε0µ0 ∇ϕ, (9) donde k0=ω√µ0ε0=ω ces el número de ondas del vacío. Indroduciendo el gauge de Lorentz (8) en (9) se obtiene:  ∇2 A+k2 0 A=−µ0 J. (10) Esta es la ecuación de Helmholtz para el fasor del potencial vector  A. Como es bien sabido, la solución de la ecuación (10) para  Ase puede escribir mediante la expresión para 5 el potencial retardado:  A(r) = µ0 4π˚τ′  J( r′)e−jk0|r− r′| r − r′ dτ′=˚τ′ GA(r − r′)· J( r′)dτ′,(11) donde GA(r − r′)es la función de Green diádica para el potencial vector, dada por: GA(r − r′) =        µ0e−jk0|r− r′| 4π|r− r′|0 0 0µ0e−jk0|r− r′| 4π|r− r′|0 0 0 µ0e−jk0|r− r′| 4π|r− r′|        .(12) Gxx A(r − r′)ˆx=µ0e−jk0|r− r′| 4π|r− r′|ˆx(véase la ecuación (12)) sería el potencial vector creado por un elemento infinitesimal de densidad de corriente  J(r) = δ(x−x′)δ(y−y′)δ(z−z′)ˆx. A este elemento de corriente le corresponde una función de Green para el potencial escalar Gxx ϕque, de acuerdo con la ecuación (8), vendría dada por: Gxx ϕ(r − r′) = j ωε0µ0  ∇·hGxx A(r − r′)ˆxij ωε0µ0  ∇·hµ0e−jk0|r− r′| 4πr − r′ ˆxi= =−j ωε0µ0 µ0 4π"1 + jk0r − r′e−jk0|r− r′| r − r′ 3#(x−x′).(13) Esta función de Green es el potencial escalar electrodinámico creado por un dipolo eléctrico infinitesimal (dos cargas iguales de signo contrario que existen en los extremos del elemento de corriente infinitesimal debido a la ecuación de continuidad). Imaginemos que queremos calcular el potencial de una sola de esas cargas infinitesimales. De acuerdo con Michalsky [5], ese potencial es una funcion de Green Gxx ϕq (r − r′)que se relaciona con la función de Green del dipolo Gxx ϕ(r − r′)mediante la ecuación: Gxx ϕ(r − r′) = 1 jω ∂Gxx ϕq (r − r′) ∂x′.(14) En general, si el elemento infinitesimal de corriente esta orientado en la dirección dada por el vector unitario ˆu, la relación será: Gϕ(r − r′) = 1 jω  ∇′Gϕq(r − r′)·ˆu. (15) 6 scattering) tiene que ser nulo en la superficie del cilindro conductor perfecto, se va a cumplir que: Ei z(x, y) + ESC z(x, y) = 0 (x, y)∈C. (51) O lo que es lo mismo, que: Ei z(x, y) = ωµ0 4ˆC′ H(2) 0k0p(x−x′)2+ (y−y′)2K(x′, y′)dl′(x, y)∈C. (52) La ecuación (52) es una ecuación integral de campo eléctrico (EFIE en inglés) para la densidad superficial de corriente inducida en el cilindro. En el próximo subapartado, se explicará cómo se puede resolver numéricamente esta ecuación integral mediante el Método de los Momentos (MdM). No obstante, y para facilitar la aplicación del MdM, vamos a parametrizar primero la curva Cque limita la sección transversal del cilindro en función del parámetro angular t∈[0,2π)que se define en la Figura 4. Figura 4: Parametrización usada para la curva Cque limita la sección transversal del cilindro. De acuerdo con la Figura 4, se cumple que en los puntos de la curva Cse cumple que: x(t) = ρ(t) cos(t)(53) y(t) = ρ(t) sin(t).(54) En ese caso, el diferencial de longitud que aparece en la integral de línea de la ecuación (52) viene dado por: 13 dl′=p(dx′)2+ (dy′)2=v u u t dx′ dt′!2 + dy′ dt′!2 dt′=w(t′)dt′.(55) Expresando xeyen función de tyx′ey′en función de t′y utilizando la ecuación (55), la ecuación integral (52) se puede reescribir en términos de las dos variables tyt′ como: Ei z(x(t), y(t)) = ωµ0 4ˆ2π 0 H(2) 0k0p(x(t)−x′(t′))2+ (y(t)−y′(t′))2 ·K(x′(t′), y′(t′))w(t′)dt′0≤t≤2π. (56) De hecho, la ecuación integral (56) se puede reescribir de forma más simplificada como: F(t) = ωµ0 4ˆ2π 0 M(t, t′)I(t′)dt′0≤t≤2π, (57) donde las funciones F(t),M(t, t′)eI(t′)se definen como: F(t) = Ei z(x(t), y(t)) = E0e−jk0(cos φix(t)+sin φiy(t)) (58) M(t, t′) = H(2) 0k0p(x(t)−x′(t′))2+ (y(t)−y′(t′))2w(t′)(59) I(t′) = K(x′(t′), y′(t′)).(60) La ecuación integral (57), que es una expresión más simplificada de la ecuación (56), es la que se va a resolver mediante el MdM. 2.3. Resolución de la ecuación integral mediante el Método de los Momentos con funciones base rectangulares El objetivo que se busca es despejar la función incógnita I(t)de la ecuación integral (57), que básicamente representa la densidad de corriente inducida en el cilindro conductor de la Figura 3 en función del parámetro angular t. Como ya hemos dicho anteriormente, el cálculo de I(t)lo haremos mediante el Método de los Momentos (MdM), para lo cual, aproximaremos I(t)como una combinación lineal de funciones base conocidas. Siguiendo 14 las recomendaciones del libro Computational Methods for Electromagnetics de Peterson, Ray y Mittra [9], utilizaremos pulsos rectangulares como funciones base. Esto es así porque, al no depender la componente zde la densidad de corriente de la coordenada z(en el caso de una onda incidente con polarización TMz), dicha corriente no lleva aparejada densidad de carga de acuerdo con la ecuación de continuidad (véase la ecuación (49)), y por tanto, puede ser aproximada mediante una función discontinua. En el caso del scattering por un cilindro cuando incide una onda con polarización TEz, la densidad de corriente inducida lleva apareada densidad de carga y la densidad de corriente debe ser aproximada mediante una función continua para evitar que la densidad de carga presente singularidades que no tienen sentido físico [9] (téngase en cuenta que la densidad de carga es proporcional a la divergencia de la densidad de corriente de acuerdo con la ecuación de continuidad). Vamos a dividir el intervalo 0≤t≤2πen Nsubintervalos ti−1≤t≤ti(i= 1, ..., N), cumpliéndose que t0= 0 < t1< t2< ... < tN−1< tN= 2π. Definimos el pulso rectangular i-ésimo (i= 1, . . . , N) como (véase la Figura 5): Ri(t) =    1ti−1< t < ti 0en otro caso. (61) Figura 5: Esquema de los pulsos rectangulares. Ahora aproximamos la función I(t′)como una combinación lineal de los pulsos rectangulares Ri(t), esto es: I(t′) = N X j=1 ajRj(t′),(62) donde aj(j= 1, ..., N) son coeficientes indeterminados que hay que determinar. Sustituyendo (62) en (57) se obtiene que: 15 F(t) = N X j=1 ωµ0 4ˆtj tj−1 M(t, t′)dt′aj.(63) A continuación, siguiendo el esquema del MdM, multiplicamos la ecuación (63) por Ri(t)e integramos en el intervalo ti−1<t<ti. Esto equivale a decir que en el marco del MdM, estamos utilizando funciones peso del mismo tipo que las funciones base, y por tanto, que estamos utilizando la versión “Galerkin” del MdM. Al proyectar la ecuación (63) sobre las funciones peso Ri(t)(i= 1, . . . , N), se obtienen las siguientes Necuaciones: ˆti ti−1 F(t)dt = N X j=1 ωµ0 4ˆti ti−1ˆtj tj−1 M(t, t′)dt′dtaj.(64) La expresión (64) se puede reescribir como: bi= N X j=1 Qijaj,(65) que básicamente es un sistema de Necuaciones lineales en los coeficientes indeterminados aj. En este sistema de ecuaciones, los elementos de la matriz de coeficientes Qij y los términos independientes bivienen dados por las siguientes integrales: Qij =ωµ0 4ˆti ti−1ˆtj tj−1 M(t, t′)dt′dt (66) bi=ˆti ti−1 F(t)dt. (67) Las integrales de (66) y de (67) no se pueden hacer analíticamente, y hay que hacerlas mediante integración numérica. El problema que surge al llevar a cabo la integración numérica es que la función de dos variables M(t, t′)tiene una singularidad logarítmica cuando t=t′, y esta singularidad aparece en el dominio de integración de Qij cuando i=j(véase la Figura 6), lo cual puede introducir errores importantes en el proceso de integración numérica de las citadas integrales. Se puede demostrar que en el dominio de integración de Qii (ti−1≤t≤tiyti−1≤ t′≤ti), cuando t→t′la función M(t, t′)se comporta como se indica a continuación: M(t, t′)t→t′≃ {A+Bln |t−t′|} (68) 16 Figura 6: Dominio de integración en tyt′para Qii. Las singularidades se encuentran en la recta t′=t, marcada en rojo. donde AyBson constantes, y donde se ha tenido en cuenta el desarrollo en serie en torno al origen de la función de Hankel H(2) 0(·)cuando su argumento tiende a cero [8]. La ecuación (68) muestra claramente la singularidad logarítmica de M(t, t′)cuando t→t′. Para solventar el problema de las singularidades logarítmicas, al llevar a cabo la integración numérica de las integrales de (66) y de (67), se han utilizado fórmulas de cuadratura de Ma-Rokhlin-Wandzura (MRW) [4], que están especialmente pensadas para integrar funciones con singularidades logarítmicas. Concretamente, si F(x) = h(x) + f(x) ln xes una función de una variable que es una combinación lineal de dos funciones, una h(x)que tiene un comportamiento suave y otra f(x) ln xque tiene una singularidad logarítmica en x= 0 -se supone que f(x)también tiene un comportamiento suave-, las fórmulas de cuadratura MRW permiten integrar F(x)en el intervalo 0≤x≤1con gran precisión mediante la siguiente ecuación: ˆ1 0 F(x)dx ≃ Nc X m=1 wmF(xm),(69) donde los pesos wmy los nodos xmestán tabulados en [4], y donde Nces el orden de la cuadratura. La ecuación (69) es tanto más precisa cuanto mayor es el valor de Nc. En [4] se presentan tablas de wmyxmpara Nc= 5,10,15,20,25,30,35 y40. Debido a la expresión de F(x), las fórmulas de cuadratura MRW permiten integrar, tanto funciones suaves (en ese caso f(x)=0) como funciones con singularidades logarítmicas (en ese caso 17 h(x) = 0)), con lo cual, han sido utilizadas indistintamente para hacer todas las integrales de (66) y (67). Por ejemplo, con vistas a calcular las integrales de (66) cuando i=j, definimos primero: Nj(t) = ωµ0 4ˆtj tj−1 M(t, t′)dt′.(70) Teniendo en cuenta la definición de los coeficientes Qij en (66), se cumple que: Qiji=j =ˆti ti−1 Nj(t)dt = (ti−ti−1)ˆ1 0 Nj((ti−ti−1)u+ti−1)du. (71) donde se ha introducido el cambio de variables t= (ti−ti−1)u+ti−1en (71) para que la integral vaya de 0 a 1 y así poder aplicar la fórmula de cuadratura MRW que aparece en (69). Al utilizar esta fórmula se obtiene: Qiji=j≃(ti−ti−1) Nc X m=1 wmNj(t= (ti−ti−1)xm+ti−1),(72) pero: Nj(t= (ti−ti−1)xm+ti−1) = ωµ0 4ˆtj tj−1 M(t= (ti−ti−1)xm+ti−1), t′)dt′ =ωµ0 4ˆ1 0 M(t= (ti−ti−1)xm+ti−1), t′= (tj−tj−1)v+tj−1)dv′ ≃ωµ0 4(tj−tj−1) Nc X n=1 wnM(t= (ti−ti−1)xm+ti−1), t′= (tj−tj−1)v+tj−1).(73) Sustituyendo ahora (73) en (72,) se obtiene una expresión para los Qiji=j , que está dada por: Qiji=j =ωµ0 4(tj−tj−1)(ti−ti−1) Nc X m=1 Nc X n=1 wmwnM(t= (ti−ti−1)xm+ti−1), t′= (tj−tj−1)xn+tj−1). (74) Haciendo unos cálculos similares, la integral de (67) se puede obtener mediante la ecuación: 18 bi≃(ti−ti−1) Nc X m=1 wmF(t= (ti−ti−1)xm+ti−1).(75) Como se indicó anteriormente, el cálculo de Qii es más problemático debido a la singularidad logarítmica que aparece en la ecuación (68). Por tanto, en este caso, lo que se hace es dividir la integral de (70) en dos sumandos de forma que el punto del intervalo del intervalo de intergración donde aparece la singularidad logarítmica aparezca como límite de integración en cada uno de los sumandos como ocurre en la fórmula de cuadratura de (69), y a continuación, se llevan a cabo en cada sumando los mismos pasos de las ecuaciones (72) a (74). El resultado final es una expresión para Qii basada en las fórmulas de cuadratura MRW, que está dada por: Qii ≃(ti−ti−1)2ωµ0 4·Nc X m=1 Nc X n=1 wmwnxm·Mt= (ti−ti−1)xm+ti−1, t′=−(ti−ti−1)xmxn+ (ti−ti−1)xm+ti−1+ Nc X m=1 Nc X n=1 wmwn(1 −xm)· ·Mt= (ti−ti−1)xm+ti−1), t′= (ti−ti−1)(1 −xm)xn+ (ti−ti−1)xm+ti−1.(76) Sustituyendo (74), (75) y (76) en (65), podemos resolver el sistema de ecuaciones para los coeficientes indeterminados aj(j= 1, . . . , N) de (62), y así llegar a una expresión para I(t). Una vez conocida I(t)en el intervalo 0≤t≤2π, podemos llegar a calcular la sección radar bidimensional (sección eficaz de dispersión por unidad de longitud) del cilindro conductor de la figura 3 para ondas incidentes con polarización TMz. Esta sección radar bidimensional tiene dimensiones de longitud, y de alguna manera, nos da una medida de la anchura que presenta el cilindro a la onda incidente. De acuerdo con la ecuación (1.125) de [9], la sección radar bidimensional vendría dada por: σT M (φ, φi) = l´ım ρ→∞ 2πρ|ESC z(ρ, φ)|2 |E0|2,(77) donde ESC z(ρ, φ)es el campo eléctrico de scatering de la ecuación (50) escrito en coordenadas cilíndricas y E0es la amplitud compleja del campo eléctrico de la onda incidente definido en la ecuación (46). Para poder obtener σT M (φ, φi)mediante la ecuación (77), necesitamos calcular ESC z(ρ, φ)cuando ρ≫. Para ello, reescribimos la ecuación (50) en 19 coordenadas cilíndricas como se indica a continuación: ESC z(ρ, φ) = −ωµ0 4ˆC′ H(2) 0k0p(ρcos φ−x′)2+ (ρsin φ−y′)2K(x′, y′)dl′.(78) Ahora aproximamos la raíz cuadrada que aparece en el argumento de la función de Hankel para ρ≫: p(ρcos φ−x′)2+ (ρsin φ−y′)2ρ≫≃ρ−(x′cos φ+y′sin φ).(79) Teniendo en cuenta la aproximación de (79) en (78), se obtiene: ESC z(ρ, φ)ρ≫≃ −ωµ0 4ˆC′ H(2) 0k0ρ−(x′cos φ+y′sin φ)K(x′, y′)dl′ρ≫ .(80) De acuerdo con la ec (9.2.4) del libro de Abramowitz y Stegun [8], se cumple que: H(2) 0(x)x≫≃r2 πxe−j(x−π 4)=r2j πxe−jx.(81) Haciendo uso de (81) en (80), se obtiene: ESC z(ρ, φ)ρ≫≃ −ωµ0 4r2j πˆC′ e−jk0ρejk0(x′cos φ+y′sin φ) pk0ρ−k0(x′cos φ+y′sin φ)K(x′, y′)dl′ρ≫ ≃ −ωµ0√je−jk0ρ 2√k0√2πρ ˆC′ ejk0(x′cos φ+y′sin φ)K(x′, y′)dl′.(82) Introduciendo (82) en (77) y teniendo en cuenta que ω2µ2 0 k0=k0Z2 0, se llega a que: σT M (φ, φi) = k0Z2 0 4E2 0ˆC′ ejk0(x′cos φ+y′sin φ)K(x′, y′)dl′ 2 .(83) Por analogía con las ecuaciones (52) y (56), se cumple que: ˆC′ ejk0(x′cos φ+y′sin φ)K(x′, y′)dl′=ˆ2π 0 ejk0(x′(t′) cos φ+y′(t′) sin φ)K(x′(t′), y′(t′))w(t′)dt′ =ˆ2π 0 S(t′)I(t′)dt′,(84) 20 donde: S(t′) = ejk0(x′(t′) cos φ+y′(t′) sin φ)w(t′)(85) eI(t′)ha sido definida en la ecuación (60). Sustituyendo ahora (62) en (84) y teniendo en cuenta la definición de los pulsos rectangulares dada en (61), llegamos a que: ˆC′ ejk0(x′cos φ+y′sin φ)K(x′, y′)dl′= N X j=1 ajˆ2π 0 S(t′)Rj(t′)dt′= N X j=1 ajˆtj tj−1 S(t′)dt′.(86) Sustituyendo (86) en (83), podemos expresar σT M (φ, φi)como: σT M (φ, φi) = k0Z2 0 4E2 0 N X j=1 ajˆtj tj−1 S(t′)dt′ 2 .(87) Las integrales de (87) se pueden calcular utilizando cuadraturas MRW como ocurre con las integrales de (66) y (67). siguiendo un procedimiento análogo al seguido en las ecuaciones (72) a (74), se deduce que: ˆtj tj−1 S(t′)dt′≃(ti−ti−1) Nc X m=1 wmS(t= (ti−ti−1)xm+ti−1).(88) Por lo tanto y a modo de resumen, para calcular σT M (φ, φi)basta con utilizar las ecuaciones (74), (75) y (76) para calcular los coeficientes Qij ybi(i, j = 1, ..., N), con resolver el sistema de ecuaciones (65) para obtener ai(i= 1, ..., N), y con utilizar (88) para obtener las integrales de (87). 21 3. Resultados numéricos Teniendo en cuenta los resultados obtenidos anteriormente, podemos calcular distintas magnitudes como la sección radar bidimensional σT M o la densidad de corriente K para cilindros de distintas secciones. Con el fin de validar nuestros códigos en MATLAB, empezaremos por mostrar resultados para el cilindro de sección circular ya que para este caso existen resultados analíticos. 3.1. Validación para el caso particular del cilindro de sección circular (sección radar bidimensional y densidad de corriente) Sea un cilindro de sección circular de radio R. Sin pérdida de generalidad, en todos los casos tratados en este subapartado supondremos que la onda incidente sobre el cilindro viene de izquierda a derecha en la Figura 3 con un ángulo de incidencia φi= 0. La parametrización de curva Cque limita la sección transversal del cilindro (véase la Figura 7) vendrá dada por las ecuaciones: x(t) = Rcos(t)(89) y(t) = Rsin(t).(90) Figura 7: Sección transversal circular. El diferencial de longitud dl′de la ecuación (52) se calcula mediante la ecuación (55): dl′=p((−Rsin(t′))2+ ((Rcos(t′))2)dt′=Rdt′=w(t′)dt′⇒w(t′) = R. (91) 22 cilindros con secciones elíptica, ojival y rectangular. Asimismo, en los tres casos llevaremos a cabo un análisis de convergencia en función del número de funciones base utilizadas en el MdM para un Ncfijo (esto es, haremos análisis de convergencia del tipo mostrado en la Figura 12). En todos los casos asumiremos que la onda incidente viaja de izquierda a derecha con un ángulo de incidencia φi= 0 y, mientras no se indique lo contrario, los cálculos se harán con un orden de cuadratura Nc= 15 y un número de funciones base N= 100. En primer lugar, estudiaremos el caso del cilindro de sección elíptica que se muestra en la Figura 13. La parametrización de la curva Cque limita la sección transversal del cilindro vendrá dada por las ecuaciones: x(t) = acos(t)(97) y(t) = bsin(t),(98) donde aybson los semiejes horizontal y vertical, respectivamente. Figura 13: Sección transversal elíptica. El diferencial de longitud dl′la ecuación (55) viene dado por: dl′=p((−asin(t′))2+ ((bcos(t′))2)dt′=w(t′)dt′.(99) En la Figura 14 se representa la sección radar normalizada a la longitud de onda para cilindros elípticos en los que b/λ = 0.4ya/λ toma valores 0.2, 0.4 y 0.6. En el caso en el que a=b= 0.4λlos valores obtenidos coinciden con los que se representan en la Figura 8 para un cilindro circular en el que R= 0.4λ, lo cual sirve para validar el código implementado para cilindros elípticos. Al igual que en el caso de la Figura 8, en la Figura 29 14 se observa que la sección radar como función de φtiene un máximo en la dirección de propagación de la onda (φ= 0). Figura 14: Representación de la sección radar bidimensional en función del ángulo φpara cilindros de sección elíptica. En la Figura 15 se representa la densidad de corriente inducida en los cilindros elípticos estudiados en la Figura 14. Se observa que el máximo de la densidad de corriente en la dirección φ=πse va haciendo más agudo conforme se va achatando el cilindro elíptico en dirección vertical, esto es, conforme a/b va disminuyendo. Esto es lógico si pensamos que en una tira conductora infinita, que sería un cilindro elíptico en el que a/b = 0, la densidad de corriente en los bordes de la tira es singular [6], con lo cual, lo que ocurre es que el comportamiento del cilindro elíptico va tendiendo al de la tira conductora conforme ava disminuyendo. En la Figura 15 también se observa que la forma de la densidad de corriente en el caso a=b= 0.4λcoincide con la que se aprecia en la Figura 9 para un cilindro circular en el que R= 0.4λ, lo cual sirve una vez más para validar los resultados obtenidos para el cilindro elíptico. En la Figura 16 se representa la sección radar monoestática de cilindros elípticos frente aa/λ, suponiendo que b a= 2. Como era de esperar, al aumentar el tamaño del cilindro, aumentar su sección radar monoestática, y lo hace de forma “cuasi-lineal” como ocurría 30 Figura 15: Módulo de la densidad superficial de corriente normalizada frente al ángulo φ en cilindros elípticos. en el caso de los cilindros de sección circular (véase la Figura 10). En la Figura 17 se lleva a cabo un estudio de convergencia con el número de funciones base Npara la sección radar monoestática de un cilindro elíptico, que es similar al que se llevó a cabo para un cilindro circular en la Figura 12. Se observa que el ritmo de convergencia en la Figura 17 es muy similar al que se observa en la Figura 12. Si tomamos como valores virtualmente exactos los que se obtienen cuando N= 340, cuando N= 10 se está obteniendo una cifra significativa de precisión en los valores de la sección radar monoestática, cuando N= 100, se obtienen cuatro cifras significativas, y cuando N= 260, se obtienen seis cifras significativas. Estos resultados justifican haber utilizado N= 100 en nuestros cálculos para cilindros elípticos. A continuación, vamos a obtener resultados para cilindros de sección ojival con una sección del tipo mostrado en la Figura 18. La ojiva está formada por dos porciones iguales de circunferencia, que son imagen especular una de la otra a través del plano y= 0. La parametrización de la curva Cque limita la sección ojival vendrá dada por las ecuaciones (53) y (54) con la siguiente función ρ(t): 31 Figura 16: Sección radar monoestática para cilindros elípticos en función del tamaño del semieje horizontal apara b a= 2. Figura 17: Análisis de convergencia con Npara cilindros de sección elíptica. Nc= 20, a/λ = 0,4yb= 2a. 32 Figura 18: Sección transversal ojival. ρ(t) =    −(R−c) sin(t) + p(R−c)2sin2(t) + (2Rc −c2) 0 ≤t≤π (R−c) sin(t) + p(R−c)2sin2(t) + (2Rc −c2)π≤t≤2π, (100) donde: R=c2+d2 2c.(101) Aplicando la ecuación (55), se puede demostrar que el diferencial de longitud dl′en este caso viene dado por: dl′=sdρ dt′2 + (ρ(t′))2dt′=w(t′)dt′.(102) donde la derivada de la coordenada radial ρ(t)con respecto al parámetro tviene dada por: dρ dt =     −(R−c) cos(t) + (R−c)2sin(t) cos(t) √(R−c)2sin2(t′)+(2Rc−c2)0≤t≤π (R−c) cos(t) + (R−c)2sin(t) cos(t) √(R−c)2sin2(t)+(2Rc−c2)π≤t≤2π. (103) A la hora de aplicar el MdM al análisis del scattering por una sección ojival, esta sección presenta ciertas dificultades que no nos encontrábamos en los casos de las secciones circular y elíptica. La ojiva tiene esquinas en φ= 0 yφ=π, y en esas esquinas la corriente inducida en la superficie de la ojiva va a ser singular [6]. Si alguno de los pulsos rectangulares 33 elegidos como funciones base muestrea la corriente en una esquina, ello puede a dar a lugar a errores numéricos en la aplicación del MdM. Por tanto, es conveniente escoger los pulsos rectangulares de manera que las esquinas estén en los extremos de los intervalos en los que están definidos. En el caso de la sección ojival, solventar este problema es tan sencillo como elegir un número de funciones base Nque sea par. En la Figura 19 se presentan resultados para la sección radar de cilindros de sección ojival en función del ángulo φ. Se observa que conforme se va achatando la ojiva en dirección vertical (esto es, conforme c/d va disminuyendo), el mínimo de la sección radar para φ=πva dismninuyendo, y esos mínimos son más profundos que los que se observan para cilindros circulares y elípticos de dimensiones similares. También se observa que los valores obtenidos en la Figura 19 para la sección radar de una sección ojival en la que c=d= 0.6λcoinciden con los obtenidos en la Figura 8 para un cilindro circular en el que R= 0.6λ, lo cual sirve para validar el código implementado para cilindros de sección ojival. Figura 19: Representación de la sección radar bidimensional en función del ángulo φpara cilindros de sección ojival. En la Figura 20 se presentan resultados para la densidad de corriente en cilindros de sección ojival. Se observa de nuevo que los resultados obtenidos para el caso c=d= 0.6λ 34 coinciden con los resultados obtenidos en el caso R= 0.6λen la Figura 9. La Figura 20 pone de manifiesto la singularidad de la corriente en la esquina φ=πde las ojivas cuando c= 0.5λy cuando c= 0.4λ. La singularidad es tanto más abrupta cuanto menor es la razón c/d. La densidad de corriente no se puede ir a infinito en esas esquinas porque no hemos utilizado funciones base que den cuenta de la singularidad de borde. Si las hubiéramos utilizado, estaríamos dando cuenta de forma exacta de la singularidad de la corriente en las esquinas [6]. Los resultados obtenidos en la Figura 20 se pueden comparar con los que aparecen en la Figura 2.3 de [9] para una ojiva con una esquina más pronunciada (c= 1.5λyd= 5λ) observándose que la forma cualitativa de la densidad de corriente es esencialmente la misma. Figura 20: Módulo de la densidad superficial de corriente normalizada frente al ángulo φ en cilindros de sección ojival. La Figura 21 muestra resultados para la sección radar monoestática de cilindros de sección ojival frente a c/λ cuando c/d=2/3. De nuevo, la sección radar monoestática crece al crecer las dimensiones del cilindro pero, a diferencia de lo que ocurre en las Figuras 10 y 16, ya no crece de forma “cuasi-lineal”. Es más, se puede comprobar que los valores de la sección radar monoestática de los cilindros de sección ojival son sustancialmente inferiores a los de cilindros de sección circular y elíptica, lo cual significa que la energía 35 Figura 21: Sección radar monoestática para cilindros de sección ojival en función del tamaño del semieje vertical cpara c/d = 2/3. Figura 22: Análisis de convergencia con Npara para cilindros de sección ojival. Nc= 20, c= 0.4λyd= 0.6λ. 36 “rebotada” por el objeto al incidir la onda es muy inferior. De alguna manera, la arista de la ojiva, dispersa la onda en todas las direcciones en lugar de rebotarla hacia atrás. Esto podría explicar por qué los aviones con tecnología stealth tienen muchas aristas en su fuselaje para conseguir reducir la sección radar. La Figura 22 muestra la convergencia del MdM para la sección radar monoestática de cilindros de sección ojival con el número de funciones base N. Es evidente que la convergencia es peor que en el caso de cilindros de sección circular y elíptica, lo cual se atribuye a que las funciones base no son capaces de reproducir con precisión las singularidades de borde de la densidad de corriente en las esquinas. Si tomamos como valores virtualmente exactos los que se obtienen cuando N= 340, cuando N= 20 se está obteniendo una cifra significativa de precisión en los valores de la sección radar monoestática, cuando N= 100, se obtienen aproximadamente tres cifras significativas, y cuando N= 260, se obtienen cuatro cifras significativas. Por último, vamos a estudiar el caso de un cilindro de sección rectangular. La sección va a ser un rectángulo de semilados ayb, tal y como se indica en la Figura 23. Figura 23: Sección transversal rectangular. En el caso del cilindro de sección rectangular, la parametrización de la curva Cque 37 limita la sección rectangular vendrá dada por las ecuaciones: x(t) =                      a0< t < φo bcotg(t)φo< t < π −φo −a π −φo< t < π +φo −bcotg(t)π+φo< t < 2π−φo a2π−φo< t < 2π (104) y(t) =                      atg(t) 0 < t < φo b φo< t < π −φo −atg(t)π−φo< t < π +φo −b π +φo<t<2π−φo atg(t) 2π−φo<t<2π, (105) donde φ0= arctanb a. La función w(t′)que aparece en la ecuación (55) puede ser calculada fácilmente, obteniéndose: w(t′) =                      asec2(t′) 0 < t′< φo bcsc2(t′)φo< t′< π −φo asec2(t′)π−φo< t′< π +φo bcsc2(t′)π+φo< t′<2π−φo asec2(t′) 2π−φo< t′<2π. (106) De manera similar a lo que ocurre en la sección ojival, la aplicación del MdM al cilindro de sección rectangular plantea problemas porque el cilindro presenta cuatro esquinas en los que la densidad de corriente es singular [6]. Para evitar problemas numéricos, los intervalos en los que están definidos los pulsos rectangulares utilizados como funciones base no deben incluir a estas cuatro esquinas. Los cálculos que se presentan en este subapartado se han realizado para una sección cuadrada (φ0=π 4) y para 2 secciones rectangulares distintas en las que φ0=π 6yφ0=π 3. En todos los casos analizados se ha utilizado un número de funciones base N= 96 (ya que este número es próximo a 100, que es el número de funciones base utilizados para los cálculos de los cilindros de secciones elíptica y ojival). El número 96 es divisible por 8 y por 12, lo cual es necesario para que los intervalos de los pulsos rectangulares no incluyan a las esquinas en los casos particulares de secciones 38 Referencias [1] C. A. Balanis, Advanced Engineering Electromagnetics, John Wiley &Sons, Inc., pp. 570-658, New York, 1989. [2] David M. Pozar, Microwave Engineering, 2nd edition, John Wiley &Sons, Inc., pp. 673-679, New York, 1998. [3] Wikipedia. Stealth Technololgy. URL: https://en.wikipedia.org/wiki/Stealth_ technology (Accedido: 29/08/2022). [4] J. Ma, V. Rokhlin and S. Wandzura, “Generalized Gaussian Cuadrature Rules for Systems of Arbitrary Functions", SIAM Journal Num. Anal., Vol. 33, No. 3, June 1996, pp. 971-996. [5] K. A. Michalsky, “On the Scalar Potential of a Point Charge Associated with a Time-Harmonic Dipole in a Layered Medium", IEEE Transactions on Antennas and Propagations, Vol. AP-35, No. 11, November 1987, pp. 1299-1301. [6] J. Van Bladel, Electromagnetic Fields, 2nd edition, IEEE Press/Wiley-Interscience, pp. 172-175, pp. 749-751 and p. 1027, Hooboken, New Jersey, 2007. [7] I. S. Gradshteyn and I. M. Ryzhik, Table of integrals, series and products, 7th edition, Elsevier & Academic Press, p. 482, Amsterdam, 2007. [8] M. Abramowitz, Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables, 9th printing, Dover Publications, pp. 358-389, New York, 1970. [9] A. F. Peterson, S. L. Ray and R. Mittra, Computational Methods for Electromagnetics, IEEE Press and Oxford University Press, pp. 1-93, New York, 1998. 45