Repositorio Institucional de Documentos

A Finite Element Approach To Inverse Scattering Problems PROYECTO FIN DE CARRERA ESCUELA DE INGENIERÍA Y ARQUITECTURA DE ZARAGOZA INGENIERÍA DE TELECOMUNICACIÓN AUTOR: ENRIQUE MASGRAU RITE DIRECTOR: GIUSEPPE VECCHI PONENTE: ENRIQUE MASGRAU GÓMEZ 22 de noviembre de 2012 Agradecimientos Durante este último año he vivido una de las mejores experiencias vitales de mi vida. Torino es el lugar que ha sido mi casa durante todo este tiempo, y de donde yo me llevaré fantásticos amigos para siempre. En este país, me he sentido como uno más, además de recibir un gran afecto por parte de toda su gente. Realmente, Torino estará siempre dentro de mí cuando esta aventura termine. Torino además es el lugar donde yo decidí terminar mis estudios universitarios, incluido este trabajo. Quiero dar las gracias al Politecnico di Torino por darme esta posibilidad. Estudiar en el extranjero me ha hecho crecer realmente como persona. Ahora es el momento de dar las gracias a todas las personas que me han apoyado a lo largo de este tiempo a realizar mi Proyecto Fin de Carrera. Gracias también al Istituto Superiore Mario Boella, lugar donde he desarrollado mi proyecto final de carrera. En el ISMB yo he podido disfrutar de un ambiente agradable de trabajo. Me gustaría dar mis más sínceras gracias al pofesor Giuseppe Vecchi por darme la oportunidad de realizar el PFC con él y a Elia A.Attardo por todo su apoyo, esfuerzo y paciencia durante todos estos meses. Del mismo modo, me gustaría mostrar mi gratitud a mi familia por su guía y apoyo a lo largo de estos años. Sin ellos, esta experiencia no podría haber sido posible. Finalmente, quisiera dar las gracias a mis amigos de España, Italia y Erasmus en general que me han ayudado a disfrutar de mi trabajo, estudios y de los buenos momentos. De verdad, gracias a todos. 3 Resumen En este proyecto se desarrolla una aplicación basada en "Microwave Imaging" (MWI), implementada con el propósito de obtener las propiedades materiales desconocidas de un determinado objeto. Este método consta de dos partes fundamentales: El Problema Directo y el Algoritmo Inverso. La primera etapa, el Problema Directo, implica la medición del campo eléctrico disperso a lo largo de un dominio material. El Objeto de Interés (OI) representa un cuerpo material desconocido que se encuentra dentro de un tanque, rodeado por un conjunto de antenas que iluminan el escenario y almacenan los datos experimentales de las mediciones relacionadas con el campo eléctrico. Además, se introduce un medio adaptador dentro del recinto cerrado, el cual es denominado como "background". Con el objetivo de simular los datos experimentales, se desarrolla el Método de los Elementos Finitos. FEM representa una técnica matemática e ingenieril muy potente que nos permite resolver un conjunto de ecuaciones lineales que describen el comportamiento electromagnético. De este modo, podremos generar los datos sintéticos referidos a las variables nodales, que definen el escenario de imagen simulado mediante un simulador de FEM, denominado GiD y desarrollado por la UPC. Después de resolver el Problema Directo, se aborda el "Contrast Source Inversion Method" (CSIM) con el propósito de reconstruir los parámetros físicos originales que definen el OI. Haciendo uso de este algoritmo de inversión será viable alcanzar el error mínimo global entre los datos reales y los reconstruidos. Cuando este método iterativo converja, los resultados reconstruidos serán analizados con el objeto de identificar los materiales implicados en el "Imaging Domain". En este trabajo se describen los diferentes experimentos relacionados con el Problema Directo yAlgoritmo Inverso, obteniendo diversas conclusiones sobre el funcionamiento de FEM-CSIM. En concreto, se analizan los conductores eléctricos perfectos, la distribución de las fuentes de corriente, las própiedades dieléctricas del "background" y la influencia de la frecuencia. Del mismo modo, los resultados de reconstrucción serán comparados en diferentes experimentos, obteniendo información sobre la resolución del método y las limitaciones del algoritmo. Finalmente es importante destacar las simulaciones realizadas con medios con y sín pérdidas, y los experimentos de biomedicina que tratan de representar posibles experimentos reales de imagen médica. Observaremos como CSIM proporciona una calidad alta en los resultados cómo para poder detectar la posicion y características de los objetivos. En las mejores situaciones de reconstrucción obtendremos errores en torno al 25%, que aunque puedan parecer discretos, son suficientes en muchas aplicaciones de imagen médica. 5 Índice 1 Introducción 15 1.1 Ámbito de la Aplicación ..................................... 15 1.2 Motivación ............................................ 16 1.3 Definición del problema ..................................... 16 1.4 Solución propuesta ........................................ 17 1.5 Contenido del Proyecto ..................................... 17 2 Conceptos Electromagnéticos 19 2.1 Ecuaciones de Maxwell ...................................... 19 2.2 Condiciones de Contorno .................................... 20 2.3 Ecuación de Helmholtz ...................................... 21 2.4 Ecuaciones Para Ondas Escalares ................................ 21 3 El Problema Directo 23 3.1 Punto de Partida ......................................... 23 3.2 Ondas Electromagnéticas .................................... 24 3.3 Métodos de Simulación ...................................... 25 4 El Método De Los Elementos Finitos 27 4.1 Conceptos Básicos de FEM ................................... 27 4.2 Discretización del Dominio ................................... 28 4.3 Funciones Base .......................................... 29 4.3.1 Funciones Base Nodales ................................. 30 4.3.2 Funciones Base Vectoriales ............................... 31 4.4 Problemas Escalares ....................................... 32 7 4.4.1 Problema de Valores en la Frontera en 2D ....................... 33 4.4.2 Resolviendo BVP Usando el Método de Garlerkin .................. 33 4.5 Problemas Vectoriales ...................................... 36 4.5.1 La Ecuación "Curl-Curl" y Los Elementos Vectoriales ................ 36 4.6 Absorbiendo Condiciones de Contorno ............................. 37 4.6.1 "The Perfectly Matched Layer" en 2D ......................... 38 4.7 Implementación de la Matriz FEM ............................... 41 4.8 Simulaciones ........................................... 41 4.8.1 Recinto Circular ..................................... 41 4.8.2 Recinto Triangular .................................... 46 4.8.3 Recinto Cuadrado .................................... 47 5 Resolución de Problemas Inversos de Dispersión Electromagnética 49 5.1 El Algoritmo de Inversión .................................... 49 5.2 Operadores Matriciales de Inversión .............................. 50 5.3 The Contrast Source Inversion Method ............................. 51 5.3.1 Normas y Productos Internos .............................. 52 5.3.2 Estimación Inicial de FEM-CSIM ............................ 53 5.4 Simulaciones ........................................... 53 6 Resultados y Conclusiones 61 6.1 Limítaciones del CSIM ...................................... 61 6.1.1 Influencía del Contraste ................................. 61 6.1.2 Frecuencia ......................................... 64 6.2 Recintos PEC ........................................... 67 6.3 Elección del "Background" ................................... 69 6.3.1 Medios Con Pérdidas ................................... 69 6.3.2 Medios Sin Pérdidas ................................... 72 6.4 Aplicaciones Biomédicas ..................................... 74 6.5 Conclusiones ........................................... 78 6.6 Trabajos Futuros ......................................... 79 Apéndice 85 A Documento original de la tesis 85 8 Lista de Figuras 1.1 Resultados de una aplicación MWI ............................... 16 1.2 Descripción del proceso ..................................... 17 3.1 Permitividad relativa (línea azul) y conductividad (línea roja) para glycerine-water 80:20 de 300MHz hasta 3GHz ..................................... 24 3.2 Sistema de Imagen Cerrado. ................................... 24 3.3 (a) Protótipo de un sístema MWT (b) Lectura de tomografía ................ 25 3.4 (a) modelo MWI 2D (b) modelo MWI 3D ......................... 26 4.1 Elementos diferentes para discretización del dominio en 1D, 2D y 3D ............ 28 4.2 Error de discretización derivado del uso de triángulos o cuadrados. ............. 29 4.3 Malla adaptativa de una geometría en 2D. ........................... 29 4.4 Numeración de los nodos locales de cada elemento e. ..................... 30 4.5 Elemento triangular basado en elementos vectoriales ..................... 32 4.6 MWI application results ..................................... 37 4.7 PML encerrando un dominio infinito .............................. 40 4.8 Problema Directo resuelto con FEM usando funciones nodales: 0 r= 20, (Xs,Ys)=(0.5,0), (a)(b) |E|y∠Epara σ= 0.0005 S/m (c)(d) |E|y∠Epara σ= 0.5S/m. ......... 42 4.9 Problema Directo resuelto con FEM usando funciones nodales: 0 r= 20 (a)(b) |E|y∠E para σ= 0.0005 [S/m] con fuentes (Xs,Ys)=([0.05 0],[0 -0.05]), (c)(d) |E|y∠Epara σ= 0.0005 [S/m] con fuentes (Xs,Ys)=([-0.05 0.05],[-0.05 0.05]). ............... 42 4.10 Problema Directo resuelto con FEM usando funciones nodales: 0 r= 80 (a)(b) |E|y∠E para σ= 0.0005 [S/m] con (Xs,Ys)=(0.05,0), (c)(d) |E|para σ= 0.0005 [S/m], con (Xs,Ys)=(0.02,0) y (Xs,Ys)=(0.07,0). .............................. 43 9 16 CAPÍTULO 1. INTRODUCCIÓN 1.2 Motivación El Problema Directo y el Método Inverso serán implementados y evaluados separadamente para obtener un método global que nos permite conocer las propiedades materiales de un blanco desconocido, que puede estar localizado dentro de un cuerpo conocido. Conociendo estas propiedades del material es posible entender qué tipo de material se está analizando. Esta situación implica un gran beneficio en la ingenería biomédica, y más exactamente, en el tratamiento de imagen en el ámbito de la medicina. Por ejemplo, la tomografía nos permite establecer elcamino a la detección de tumores de mama, diagnosis de huesos y otro tipo de enfermedades. Estos problemas reflejan las razones de la elección de este proyecto. MWI tiene un importante impacto en la ingenería médica y es innegable lo adecuado de este tipo de estudios para el desarrollo de la medicina. Por ello, la ingenería debe investigar en el campo de la medicina. Es un orgullo participar en el estudio e investigación de aplicaciones de bioingenería que resultan tan importantes para la mejora de la vida. 1.3 Definición del problema Esta sección presenta una descripción del escenario en el que se desarrolla este trabajo. Supóngase que existe un cuerpo con unas cualidades determinadas y desconocidas, y estas propiedades son el objeto del estudio e investigación. Es necesario entender el comportamiento de este material para poder establecer conclusiones que pueden ser de gran importancia para el éxito de la investigación. Por ejemplo para diagnosticar un posible problema de salud (cáncer de mama). En la Fig.1.1 se muestra un escenario médico, una detección de cáncer de mama. En este caso las propiedades materiales del pecho son conocidas, y si se aplica un correcto análisis, será posible detectar la existencia de cuerpos extraños que podrían representar un pequeño tumor. Hay que tener muy en cuenta que el análisis del OI requiere que se aplique un método que nos permita descubrir las propiedades dieléctricas del cuerpo desconocido, sin ser alterado y mucho menos, ser dañado durante el proceso. Figura 1.1: Resultados de una aplicación MWI 1.4 Solución propuesta 17 1.4 Solución propuesta Existen diferentes caminos y métodos para solventar este problema. En este proyecto el método usado para resolver este problema biomédico es conocido como FEM-CSIM (FEM-Contrast Source Imaging Method). FEM será uno de los principales objetos de estudio, ya que permite computar de forma muy eficaz los parámetros electromagnéticos, además de representar un método muy popular para analizar escenarios en infinidad de aplicaciones ingenieriles como diseño de estructuras y análisis de temperatura. Si FEM-CSIM es implementado correctamente, se obtendrán las propiedades dieléctricas del OI con un pequeño error. Más tarde, estos resultados tendrán que ser analizados e interpretados por especialistas para determinar posibles problemas de salud. Figura 1.2: Descripción del proceso Como ya fue mencionado anteriormente, podemos dividir la aplicación en dos partes: Problema Directo, resuelto mediante la implementación de FEM, y el Algoritmo Inverso. Gracias al método iterativo CSIM podremos obtener las propiedades eléctricas del OI. En la Fig.1.2, se describe el esquema general del proceso. 1.5 Contenido del Proyecto Esta memoria esta compuesta por seis capítulos, siendo el objetivo del autor describir correctamente el método, los conceptos implicados en el proceso y, obviamente, los resultados y correspondientes conclusiones de la aplicación. El Capítulo 2 aborda los conceptos electromagnéticos que han sido considerados durante este proyecto, principalmente los principios del electromagnetismo y la formulación de las ecuaciones de Maxwell. Además es necesario explicar otros conceptos imprescindibles en este trabajo como el estudio de las 18 CAPÍTULO 1. INTRODUCCIÓN condiciones de contorno. El Problema Directo es presentado en el Capítulo 3, describiendo el escenario de MWI y cómo este problema es resuelto en experimentos reales, comparándolo con los escenarios simulados mediante FEM, que es el objeto de estudio del capítulo siguiente. En el Capítulo 4 se describe detalladamente el Método de los Elementos Finitos (FEM). La complejidad de este método implica un estudio en profundidad y una importante comprensión del funcionamiento e implementación. También se describen sus aspectos matemáticos y geométricos centrados en los escenarios 2D. Sin embargo, se incluye una explicación breve sobre los problemas en 3D. El problema de dispersión inversa es objeto de estudio en el Capítulo 5, donde se desarrolla e implementa FEM-CSIM y se evalúan los resultados obtenidos. Durante el Capítulo 6 se presentarán más resultados y conclusiones de determinados experimentos con el propósito de comprender más detalladamente el comportamiento de la aplicación. Finalmente, tambíén en el Capítulo 6, se abordan los posibles trabajos futuros que pueden ser llevados a cabo con el objeto de observar la amplitud del campo de estudio que suponen las aplicaciones MWI basadas en problemas inversos. Cap´ ıtulo 2 Conceptos Electromagnéticos El electromagnetismo ha representado una parte indispensable de muchos estudios científicos e ingenieriles desde que J.C.Maxwell completase la teoría del electromagnetismo en 1873 [1]. Como es bien sabido, el problema del análisis electromagnético consiste en la resolución de un conjunto de ecuaciones de Maxwell sujetas a unas condiciones de contorno determinadas. En este capítulo se revisa la formulación matemática necesaria para implementar la aplicación de MWI. De este modo, revisaremos brevemente algunos conceptos básicos de la teoría electromagnética que se usarán en este trabajo. La formulación descrita considera tanto los casos en dos dimensiones (2D) para problemas escalares como los correspondientes a configuraciones vectoriales. 2.1 Ecuaciones de Maxwell Las ecuaciones de Maxwell explican el comportamiento del fenómeno electromagnético y pueden ser expresadas mediante formas diferenciales e integrales. El conocimiento de estas ecuaciones es esencial para una correcta compresión del estudio desarrollado en este proyecto. Las ecuaciones de Maxwell, en su correspondiente forma diferencial, se derivan de su versión integral mediante el uso de los teoremas de Strokes y Gauss. ∇ × ~ E(~r, t) = −∂~ B(~r, t) ∂t (2.1) ∇ × ~ H(~r, t) = ∂~ D(~r, t) ∂t +~ J(~r, t)(2.2) ∇ · ~ D(~r, t) = ρ(~r, t)(2.3) ∇ · ~ B(~r, t)=0 (2.4) ∇ · ~ J(~r, t) = −∂ρ ∂t (2.5) 19 20 CAPÍTULO 2. CONCEPTOS ELECTROMAGNÉTICOS donde los vectores espaciales ~ E,~ H,~ D,~ Band ~ Json respectivamente, el campo eléctrico [V/m], la intensidad del campo magnético [A/m], el flujo de densidad eléctrica [C/m2], el flujo de densidad magnética [Wb/m2]y finalmente la densidad de corriente eléctrica [A/m2], en la ecuacione (2.3) la variable ρse refiere a la carga eléctrica. La densidad de corriente eléctrica ~ Jes la suma de dos contribuciones diferentes: la densidad de corriente de conducción ~ Jc, que implica la capacidad del medio para conducir corriente eléctrica, y la densidad de corriente impuesta ~ Ji, debida a las fuentes de corriente impuestas sobre el medio. Así, ~ J(~r, t) = ~ Jc(~r, t) + ~ Ji(~r, t)(2.6) Además es necesario introducir algunas ecuaciones adicionales para presentar correctamente la teoría desarrollada por Maxwell. Para medios isotrópicos, homogéneos y no dispersivos algunas relaciones de interés a considerar son las siguientes: ~ D(~r, t) = or(~r)~ E(~r, t)(2.7) ~ B(~r, t) = µoµr(~r)~ H(~r, t)(2.8) ~ Jc(~r, t) = σo(~r)~ E(~r, t)(2.9) En estas expresiones aparecen parámetros relacionados con las propiedades materiales como la permitividad en el vacío o, la permitividad relativa compleja r, la permeabilidad en el vacío µo, la permeabilidad relativa µry la conductividad σo. 2.2 Condiciones de Contorno Las ecuaciones diferenciales expuestas en la Sección 2.1 pueden ser resultas si se consideran las correspondientes condiciones de contorno de los medios presentes. En otras palabras, es necesaria una completa descripción de las interfaces entre medios con el objeto de obtener soluciones reales. Algunas condiciones de frontera usadas en la práctica son descritas: ˆn×(~ E1(~r)−~ E2(~r)) = 0 (2.10) ˆn×(~ H1(~r)−~ H2(~r)) = 0 (2.11) ˆn×(~ D1(~r)−~ D2(~r)) = 0 (2.12) 2.3 Ecuación de Helmholtz 21 ˆn×(~ B1(~r)−~ B2(~r)) = 0 (2.13) En las expresiones (2.11) y(2.12), tanto la densidad de corriente eléctrica superficial como la carga eléctrica superficial son consideradas nulas. En el caso de que este supuesto no se cumpla, las expresiones son, ˆn×(~ H1(~r)−~ H2(~r)) = ~ Js(2.14) ˆn×(~ D1(~r)−~ D2(~r)) = ρs(2.15) 2.3 Ecuación de Helmholtz En esta aplicación es necesario resolver ecuaciones diferenciales parciales (PDE) del campo eléctrico ~ E. Para obtener la siguiente expresión es necesario eliminar la componente magnética ~ H, ∇ × ∇ × ~ E(~r)−w2µoo(~r)~ E(~r) = −jwµo~ Ji(~r)(2.16) En la expresión (2.16) se han asumido dos condiciones: (i) no existen cargas eléctricas (ρ=0), (ii) la permeabilidad relativa es nula (problemas no magnéticos) µr= 1. La ecuación (2.16) es conocida como la ecuacion de Helmholtz. Durante la descripción de FEM para problemas vectoriales veremos como se refiere a ella con el nombre de "curl-curl equation". 2.4 Ecuaciones Para Ondas Escalares En el análisis electromagnético, siempre que sea posible, se utiliza una formulación simplificada de los problemas haciendo uso de modelos en 2D como aproximación a los problemas en 3D. Podemos definir el escenario escalar de la siguiente manera, ∇2Ez+k2 orEz=jwµoJz(2.17) también conocida como ecuación para ondas escalares inhomogéneas. Esta expresión será usada durante el Capítulo 4 para describir los problemas escalares en 2D correspondientes al caso de una polarización TM-z. 22 CAPÍTULO 2. CONCEPTOS ELECTROMAGNÉTICOS Cap´ ıtulo 3 El Problema Directo "Microwave Imaging" (MWI) es de interés para diversas aplicaciones, como estudios geofísicos e imágenes médicas. En la aplicación de MWI considerada en este proyecto se trata de reconstruir cuantitativamente las propiedades eléctricas (i.e. permitividad y conductividad), en su mayor parte desconocidas, de un objeto de interés (OI) que está sumergido en un medio "background" de propiedades diléctricas conocidas [9]. Como se comentó en el Capítulo 1, para resolver una aplicación de MWI es necesario implementar el Problema Directo con el propósito de obtener las distribución del campo eléctrico. En este capítulo, describiremos como se realiza la resolución del Problema Directo en aplicaciones experimentales reales en comparación a las apliaciones basadas en problemas simulados. 3.1 Punto de Partida Antes de explicar como se implementa el Problema Directo mediante simulación, es interesante y beneficioso para una mayor compresión del problema explicar cómo se realiza dicho proceso en aplicaciones médicas reales como detección de tumores o diagnósticos de huesos. Se considera un escenario de microondas en donde la región de imagen esta rodeada por una superficie eléctrica conductora. Esta superficie trabaja como una interfaz protectora de las interferencias exteriores y como contenedor de un medio adaptador que debe evitar las reflexiones producidas en las paredes. En la Fig.3.1. (extraida de [6]) se puede observar un liquído común en muchas aplicaciónes MWI basada en una solución de glicerina-agua al 80:20 por ciento. 23 24 CAPÍTULO 3. EL PROBLEMA DIRECTO Figura 3.1: Permitividad relativa (línea azul) y conductividad (línea roja) para glycerine-water 80:20 de 300MHz hasta 3GHz Un conjunto de antenas que rodean el OI, que trabajan tanto como emisores como receptores, iluminan el escenario con el objetivo de captar las mediciones EM. El OI se encuentra localizado en la región D junto con parte del "background" homogéneo (permitividad denotada como b). Figura 3.2: Sistema de Imagen Cerrado. En la Fig.3.2 [10] se puede visualizar el esquema general de la aplicación. En este caso, el recinto cerrado es un tanque cilíndrico. Nótese que durante este proyecto denotaremos las propiedades dieléctricas como =0+j00 =or=o(0 r+j00 r) = o0 r+jσ wo(3.1) Algunas posibles aplicaciones basadas en obtención de imágenes se observan en la Fig.3.3. 3.2 Ondas Electromagnéticas Cuando hablamos sobre el análisis EM en el Problema Directo, necesitamos diferenciar entre los posibles casos que pueden darse. En un medio, la propagación de una onda puede ser onda transversal electromagnética caracterizada por Ez=Hz= 0, lo cual implica que las componentes transversales son 3.3 Métodos de Simulación 25 (a) (b) Figura 3.3: (a) Protótipo de un sístema MWT (b) Lectura de tomografía nulas. Otra posiblidad es el caso de ondas TE. Esta propagación implica que Ez= 0, de forma que la unica componente transversal existente es la magnética, Hz6= 0. Finalmente, existe un caso EM basado en la presencia de una componente eléctrica no nula, Ez6= 0, mientras que la componente tranversal magnética cumple, Hz= 0. En el ambiente de aplicaciones de imagen para microondas, la formulación 2D corresponde a la polarización TM. Este problema será resuelto con solvencia mediante el uso de elementos nodales. 3.3 Métodos de Simulación El objeto de esta sección es la obtención de técnicas de simulación capaces de proporcionar un modelo sintético de posibles experimentos reales basados en MWI. De ellos se podría extraer los parámetros EM del escenario. La técnica elegida en este proyecto es denomidada como método de los elementos finitos, que implica la resolución de una ecuación matricial donde los coeficientes eléctricos que definen el campo eléctrico en nuestro escenario, pueden ser computados mediante la generación de una matriz global, conocida como matriz FEM. Esta matriz relaciona la geometría del escenario con las correspondientes propiedades eléctricas propias de éste. Esta ecuación matricial viene dada por la siguiente expresión, Az =b(3.2) Nuestra simulación directa involucra un método iterativo en el que una determinada antena funciona como emisor mientras el resto de antenas trabajan como receptores. Posteriormente, otra antena será la transmisora, de forma que el Tx de la anterior medición trabaja como Rx junto al resto de antenas. Obviamente, a mayor número de antenas mayor cantidad de información será recolectada por lo que la solución será más proxima a la real. En la Fig.3.4 (imágenes obtenidas de [3]), se definen los escenarios en 2D y 3D. 32 CAPÍTULO 4. EL MÉTODO DE LOS ELEMENTOS FINITOS Tabla 4.1: Numeración de aristas para un elemento triangular Edge No.i Node i1Node i2 1 1 2 2 2 3 3 3 1 Cada uno de esas aristas esta asociada a una función base vectorial: ~ Ne j(~r) = le iϕe i1~r∇ϕe i2~r −ϕe i2~r∇ϕe i1~r(4.7) donde le idenota la longitud de la arista i, ϕe i1~r and ϕe i2~r son las funciones nodales descritas en la expresión (4.4) de cada uno de los nodos que forman la arista e. Figura 4.5: Elemento triangular basado en elementos vectoriales 4.4 Problemas Escalares Es el momento de introducir FEM en la resolución de algunos posibles problemas en dos dimensiones relacionados a problemas escalares. Básicamente, un problema escalar se encuentra definido por una variable que debe ser determinada, por ejemplo, una polariazación TM-z. Esta sección se centra en este caso particular ya que define nuestro marco de estudio en este trabajo. Introduciendo las condiciones de contorno descritas en el Capítulo 2 surge el problema de valores de frontera (BVP) en nuestro enfoque. En definitiva, en esta sección, la formulación FEM asociada al desarrollo de las ecuaciones matriciales será utilizada para la resolución de problemas escalares. 4.4 Problemas Escalares 33 4.4.1 Problema de Valores en la Frontera en 2D En la Sección 2.4 se introdujo la expresión (2.17), que denominamos como ecuación de onda inhomogénea. Ahora tratamos de analizarla para geometrías en 2D. Nótese que las aplicaciones MWI están sujetas a un recinto basado en un PEC. Esto requiere que las condiciones de frontera seán consideradas en la formulación y computación de FEM. De esta forma, nuestro escenario escalar en 2D quedará correctamente definido por la siguiente ecuación −∇·(α∇u) + βu =g(4.8) y las condiciones de contorno u=p on Γ1(4.9) ˆn·(α∇u) + γu =q on Γ2(4.10) donde u es la función desconocida, αyβson los parámetros conocidos asociados a las propiedades físicas del dominio Ω, y g es la fuente o función de excitación. Acorde con las condiciones de frontera, Γ1define la condición de contorno de "Dirichlet" mientras Γ2se refiere a la de Robin, con Γ1+ Γ2= Γ. Γdenota la interfaz que encierra el área Ω.γ, p y q son los parámetros conocidos asociados con las propiedades físicas de la frontera. Cuando γ= 0 la condición de frontera de Robin se convierte en un caso especial denominado como la frontera Neumman. 4.4.2 Resolviendo BVP Usando el Método de Garlerkin La expresión general de BVP en 2D que aparece en (4.8) puede ser reescrita de la siguiente forma equivalente ∂ ∂x αx ∂u ∂x+∂ ∂y αy ∂u ∂y +βu =g(4.11) y las condiciones de contorno u=p on Γ1(4.12) αx ∂u ∂x ·ˆx+αy ∂u ∂y ·ˆy·ˆn+γu =q on Γ2(4.13) donde αx,αyyβson constantes. La formulación necesaria para este problema puede obtenerse construyendo el residuo ponderado de la expresión (4.11) para un simple elemento con dominio Ωe 34 CAPÍTULO 4. EL MÉTODO DE LOS ELEMENTOS FINITOS re=∂ ∂x αx ∂u ∂x+∂ ∂y αy ∂u ∂y +βu −g(4.14) Este elemento residual sería idealmente cero, lo cual implicaría que la solución numérica obtenida de u es igual a la solución exacta. Sin embargo, este no es el caso, de forma que por general, el elemento residual rees no nulo. Para minimizar el residuo, primero debemos multiplicar todos los términos de la expresión por los pesos w, y posteriormente integrar la expresión resultante a lo largo del área del elemento en cuestión, y finalmente fijarla a cero. Z ZΩe w∂ ∂x αx ∂u ∂x+∂ ∂y αy ∂u ∂y +βu −gdxdy = 0 (4.15) Haciendo uso de la Formulación Variacional para elementos finítos y del Teorema de Divergencia podemos desarrollar la expresión (4.15) hasta reducir la formulación débil de la ecuación diferencial a la siguiente expresión (Para más detalles en la obtención de la expresión (4.16) analizar ecuaciones (4.19)-(4.27) del elemento original [A.1]), −Z ZΩeαx ∂w ∂x ∂u ∂x +αy ∂w ∂y ∂u ∂y dxdy +Z ZΩe β w u dxdy =Z ZΩe w g dxdy −IΓe wαx ∂u ∂x ˆx+αy ∂u ∂y ˆydl (4.16) Asumiendo el planteamiento de Garlerkin, la función de pesos wdebe ser relacionada al mismo conjunto de funciones espaciales que interpolan la incognita principal, u. Esta incognita es interpolada usando los polinomios de Lagrange que aparecen en la expresión (4.3). Es decir, u= n X j=1 ue jNj(4.17) donde Njes la correspondiente función espacial basada en elementos nodales y nes el número total de nodos que forman el dominio Ωe(caso general, en nuestro caso n= 3). Imponiendo w=Nicon i= 1,2,3...Ney substituyendo en (4.16) obtenemos, 4.4 Problemas Escalares 35 −Z ZΩe αx∂Ni ∂x   Ne X j=1 ue j ∂Nj ∂x  +αy∂Ni ∂y   Ne X j=1 ue j ∂Nj ∂y   dxdy +Z ZΩe β Ni  Ne X j=1 ue jNj dxdy =Z ZΩe Nig dxdy −IΓe Niαx ∂u ∂x ˆx+αy ∂u ∂y ˆydl, for i = 1,2,3...Ne (4.18) La ecuación (4.18) en forma matricial resulta         Me 11 Me 12 · · · Me 1n Me 21 Me 22 · · · Me 2n . . .. . ..... . . Me n1Me n2· · · Me nn                        ue 1 ue 2 . . . ue n                +         Te 11 Te 12 · · · Te 1n Te 21 Te 22 · · · Te 2n . . .. . ..... . . Te n1Te n2· · · Te nn                        ue 1 ue 2 . . . ue n                =                fe 1 fe 2 . . . fe n                +                pe 1 pe 2 . . . pe n                (4.19) donde Me ij =−Z ZΩeαx∂Ni ∂x ∂Nj ∂x +αy∂Ni ∂y ∂Nj ∂y dxdy (4.20) Te ij =Z ZΩe β NiNjdxdy (4.21) fe i=Z ZΩe Nig dxdy (4.22) pe i=−IΓe Niαx ∂u ∂x ˆx+αy ∂u ∂y ˆydl (4.23) Es posible implentar un sístema matricial más compacto,         Ke 11 Ke 12 · · · Ke 1n Ke 21 Ke 22 · · · Ke 2n . . .. . ..... . . Ke n1Ke n2· · · Ke nn                        ue 1 ue 2 . . . ue n                =                be 1 be 2 . . . be n                (4.24) donde Ke ij =Me ij +Te ij be i=fe i+pe i (4.25) 36 CAPÍTULO 4. EL MÉTODO DE LOS ELEMENTOS FINITOS La matriz Kees la matriz FEM de cada elemento. Para implementar la matriz FEM es necesario sumar la matriz de Stiffness (4.20) y la matriz de Mass (4.21), mientras que Medepende de las condiciones de frontera del dominio, y Tedepende de las propiedades dieléctricas del dominio elemental Ωe. Tras el desarrollo llevado a cabo en el documento original (A.1), ecuaciones (4.36)-(4.41), podemos obtener una nueva versión de la ecuación (4.23) dada por, pe i=−ZΓ2 Ni(q−γ u)dl (4.26) La integral de línea en (4.26) existe solamente para los elementos de contorno, y ésta debe ser evaluado solo para las aristas localizadas a lo largo de Γ2. Para las interiores, la contribución es nula. 4.5 Problemas Vectoriales Aunque los problemas escalares definen el escenario de nuestro proyecto, creemos que es beneficioso introducir la formulación asociada a la teoría vectorial en electromagnetismo. A diferencia de los problemas escalares, los vectoriales se encuentran definidos por dos o más variables. 4.5.1 La Ecuación "Curl-Curl" y Los Elementos Vectoriales El campo electromagnético es presentado como la ecuación de Helmholtz, generando una ecuación particular que es conocida como la ecuación "curl-curl" para el campo ~ E. ∇ × µ−1∇ × ~ E(~r)−w20−jwσ~ E(~r) = −jw ~ Js(~r)on S (4.27) y las condiciones de contorno ˆn×~ E(~r) = ~p(~r)on Γ1(4.28) ˆn×(µ−1∇ × ~ E(~r)) + γ(~r)ˆn×(ˆn×~ E(~r)) = ~q(~r)on Γ2(4.29) El siguiente paso es seguir el planteamiento de Garlerkin pero utlizando elementos vectoriales. La formulación débil de integración viene dada por ZShµ−1∇ × ~ Wi·∇ × ~ E−w20−jwσ~ Wi·~ EidS +ZΓ2 ~ Wi·~q −γˆn׈n×~ Edl =−jw ZS ~ Wi·~ JsdS (4.30) 4.6 Absorbiendo Condiciones de Contorno 37 Expandiendo la solución ~ E(~r)en términos de las nuevas funciones base, ~ E(~r) = Ne X j=1 EjNj(4.31) donde Ejes el campo tangencial a lo largo de la arista j-simo, Ne es el número de elementos vectoriales e imponiendo ~ Wi=Ni, substituimos en (4.30) obteniendo un sístema lineal de ecuaciones Az =bcon Aij =ZSµ−1(∇ × Ni)·(∇ × Nj)−w20−jwσNi·NjdS +ZΓ2 γ(ˆn×Ni)·(ˆn×Nj)dl (4.32) zj=Ej(4.33) bi=−jw ZS Ni·JsdS −ZΓ2 Ni·qdl (4.34) 4.6 Absorbiendo Condiciones de Contorno La mayor parte de los problemas de dispersión descritos durante este trabajo se resuelven en espacios finitos de espacio físico. Sin embargo, existen aplicaciones de electromagnetismo donde el emplazamiento se realiza en espacio abierto. Imagínese una aplicación de medición de los parámetros de una determinada antena o la interacción de una onda incidente en una determinada estrucura. En ambos casos el campo radiado se propaga a través del espacio libre. En otras palabras, las condiciones de frontera deberían encontrarse en el infinito. Figura 4.6: MWI application results En la Sección 4.3 vimos como la principal desventaja del FEM son los costes computacionales necesarios para computar los datos electromagnéticos. De forma que el caso de una aplicación en espacio abierto implicaría unos requerimientos computacionales no aceptables. Por ello, será necesario truncar mediante barreras artificiales la region exterior con el propósito de limitar el tamaño del dominio computacional. Esta frontera deberá aparecer tan transparente como sea posible para el campo dispersivo/radiado. Es 38 CAPÍTULO 4. EL MÉTODO DE LOS ELEMENTOS FINITOS decir, las reflexiones producidas en la frontera artificial deberán ser minimizadas. Este cambio de estudio es conocido como Absorbiendo Condiciones de Contorno (ABCs). ABCs simulan o reemplazan el espacio infinito que rodea un dominio computacional finito. Esta operación nunca es perfecta, de forma que la solución computada dentro del ABC es solo una estimación de la solución real. Algunas propiedades de estas regiones artificiales deben ser: 1. ABC debe estar tan cercano como sea posible a la estructura, de forma que el dominio computacional sea lo más pequeño posible. Esto permitirá computar una cantidad de datos suficientemente aceptables. 2. Los cambios del escenario original deben ser prácticamente inapreciables con el objeto de obtener resultados lo más reales posibles. 3. Si ABC sólo puede absorber ondas planas homogéneas, ésta deberá ser colocada fuera de la región evanescente que rodea la fuente electromagnética (antena, guía de ondas). 4. Si ABC puede absorber campos evanescentes, esta deberá ser colocada lo más cercana posible a la fuente para reducir el dominio computacional. La teoría tradicional sobre el fenómeno de ABC analiza una serie de posibles implementaciones basadas en aproximaciones de Taylor de orden n-simo. Básicamente, se describe la curva del coeficiente de reflexión de nuestra interfaz y se manipula de forma que este coeficiente sea nulo para determinados ángulos de incidencia. Sin embargo, este método presenta una gran desventaja, ya que sólo para una única frecuencia obtendremos el propósito que buscamos. Este problema fue solucionando cuando en 1994 Berenger [5] propusó el concepto del "perfectly matched layer" (PML). 4.6.1 "The Perfectly Matched Layer" en 2D Una "perfectly matched layer" (PML) es una interfaz que no refleja una onda plana para ninguna frecuencia, ni ángulo de incidencia. Al mismo tiempo, permite truncar el dominio infinito a un tamaño de coste computacional aceptable. El planteamiento principal se basa en la introducción de unos parámetros conocidos como factores de estiramiento: sx,syand sz(Coordenadas de estiramiento). Además, es importante presentar un nuevo operador, ∇s, que se usa para la correspondiente modificación de las ecuaciones de onda. ∇s= ˆx1 sx ∂ ∂x + ˆy1 sy ∂ ∂y + ˆz1 sz ∂ ∂z (4.35) 4.6 Absorbiendo Condiciones de Contorno 39 de forma que ks= ˆxkx sx + ˆyky sy + ˆzkz sz (4.36) la nueva formulación de las ecuaciones de Maxwell viene dada por, ∇s×~ E(~r) = −jwµ ~ H(~r)(4.37) ∇s×~ H(~r) = jw ~ E(~r)(4.38) Definamos las coordenadas de estiramiento, (formulación extraida de [5]) sx= 1 + σx jwo s∗ x= 1 + σ∗ x jwµo (4.39) sy= 1 + σy jwo s∗ y= 1 + σ∗ y jwµo (4.40) En conclusion, es posible eliminar las reflexiones entre un medio y el PML para todos los ángulos de incidencia y rango de frecuencias imponiendo igualdad en las coordenadas de estiramiento transversales de los dos medios y forzando las coordendas longitudinales (Ver documento original para más detalles del proceso) a ser iguales a la siguiente expresión [7]: sx(x) = so(x)1−jσx(x) w0with 0=0 ro(4.41) so(x) = 1 + smπ δ2 (4.42) σx(x) = sin2πx 2δ(4.43) donde δes la espesor del amortiguador y smes el coeficiente que depende de la longitud de onda del medio. Otra posibilidad se obtiene de [1] como sx(x) = 1 −jx−δ δ2 δmax (4.44) δmax =σmax wo or δmax ≈w−1(4.45) y una buena elección para la espesor del amortiguador podría ser δ=λ/4. 40 CAPÍTULO 4. EL MÉTODO DE LOS ELEMENTOS FINITOS En la Fig.4.7 puede visualizarse un ejemplo de implementación de interfaz PML para limitar el dominio computacional de una aplicación en región abierta. Como podemos observar, es necesario colocar un región basada en PML en cada una de las posibles fronteras. Figura 4.7: PML encerrando un dominio infinito Finalmente, las distintas regiones descritas en la figura anterior quedarían definidas por los siguientes factores de estiramiento: 1. Region Ω.sx=s∗ x= 1 ∧sy=s∗ y= 1 2. Region I. sx= (4.44) ∧sy= 1 3. Region II. sx= (4.44) ∧sy= 1 4. Region III. sx= 1 ∧sy= (4.44) 5. Region IV. sx= 1 ∧sy= (4.44) 6. Region V. sx=sy= (4.44) 7. Region VI. sx=sy= (4.44) 8. Region VII. sx=sy= (4.44) 9. Region VIII. sx=sy= (4.44) Es necesario mencionar que durante las simulaciones se comprobó que el escenario descrito en la Fig.4.7 no funciona correctamente implementándose con funciones nodales. Es decir, debido a la complejidad del escenario se requiere el uso de elementos vectoriales para obtener soluciones reales. 4.7 Implementación de la Matriz FEM 41 4.7 Implementación de la Matriz FEM En el Capítulo 3 comprobamos cómo el Problema Directo podía ser resuelto evaluando el campo eléctrico dispersado por el OI al ser iluminado por un conjunto de antenas. En el caso de escenarios simulados como el nuestro, es necesario resolver la ecuación matricial descrita en (3.2). Esta puede ser desarrollada acorde a la formulación presentada en la Sección 4.4.2 para problemas escalares o la formulación descrita en la Sección 4.5.1 para problemas vectoriales. En ambos casos, es imprescindible implementar la matriz global FEM. Para la implementación de dicha matriz se debe determinar la contribución de cada uno de los elementos que forman la malla. Una vez calculadas cada una de las matrices locales FEM, debemos relacionar los índices locales de los nodos con sus correspondientes índices globales para construir con exito la matriz global. Detalles del proceso: documento original [A.1] en la Sección 4.7 y la obtención de los elementos que forman la matriz FEM puede analizarse en el apéndice [A.1] y [A.2] del documento original. 4.8 Simulaciones Finalmente, esta sección presenta diversas simulaciones basadas en geometrías en 2D con el objetivo de visualizar los resultados de computación del FEM. Todos los casos que se tratan son problemas escalares, de forma que el Problema Directo se ha ejecutado mediante el uso de funciones nodales. Estudiaremos el método para diferentes escenarios en cuanto a sus propiedades dieléctricas así como la influencia de las cargas de corriente y de la frecuencia de trabajo utilizada. 4.8.1 Recinto Circular En esta sección analizamos un PEC circular cuyo radio del tanque es igual a 0.1m. El propósito de estas simulaciones es analizar las consecuencias que implican diferentes propiedades diléctricas, frecuencias y localización de las fuentes de corriente. Primero, fijamos la frecuencia a 500MHz. Las siguientes simulaciones se caracterizan por un medio con 0 r= 20 y diferentes valores de conductividad. Consideramos una simple carga localizada en las coordenadas (Xs,Ys)=(0.5,0) Los resultados se muestran en la Fig.4.8: dependiendo de la relación entre la conductividad y la permitividad relativa (tangente de pérdidas), el número de reflexiones existentes en el medio es mayor o menor. En el caso donde la conductividad es considerablemente grande, y por tanto la relación es más grande, se observa cómo la intensidad del campo computado se localiza en el punto de la fuente, mientras que en el otro caso, donde la conductividad es suficientemente pequeña y la tangente de pérdidas es demasiado pequeña, que las reflexiones no se atenuan y provocan picos de intensidad en zonas donde no hay cargas eléctricas. 48 CAPÍTULO 4. EL MÉTODO DE LOS ELEMENTOS FINITOS Tabla 4.2: Costes Computacionales Frequency [GHz] Permittivity Edge Size [m] Num.Nodes Num.Elem. Time [s] 0.4 10 0.0158 210 366 0.220609 0.8 20 0.0056 1534 2922 4.004716 1.5 12 0.0038 3295 6376 15.41539 2 30 0.0018 14356 28266 270.684261 2.5 40 0.0013 27370 54122 1000.769933 3.2 40 9.8821e-004 47375 93940 3349.545703 del FEM-CSIM realizadas en los siguientes capítulos, el coste computacional se incrementa rápidamente debido a que el Problema Directo deberá ser resuelto para cada uno de los transmisores además de la correspondiente posterior resolución del algoritmo inverso. Y ello sin olvidar que los escenarios analizados son de una complejidad mayor. Cap´ ıtulo 5 Resolución de Problemas Inversos de Dispersión Electromagnética En este capítulo se presenta el objetivo fundamental de este trabajo, el denominado Algoritmo Inverso. Para la realización de imágenes médicas se requiere la inversión de la ecuación de onda que describe nuestro escenario electromagnético. Es decir, los resultados numéricos obtenidos durante la resolución del Problema Directo mediante FEM deben ser aplicados a un proceso de inversión con el propósito de reconstruir el origen físico que los ha generado. El método de inversión desarrollado durante este proyecto se denomina como "Contrast Source Inversion Method" (CSIM), el cual trabaja iterativamente actualizando dos variables: la fuente de contraste y las variables de constraste. Minimizando iterativamente una función, que describiremos más adelante, llegaremos a obtener las propiedades dieléctricas del OI analizado. 5.1 El Algoritmo de Inversión Una vez obtenidos los datos asociados a la dispersión electromagnética de un determinado escenario en 2D, almacenados para cada uno de los posibles transmisores y receptores que iluminan el OI (Fig. 3.2), estos se utilizan como inicialización del proceso de inversión. Pero antes, definamos la correspondiente variable de contraste para un determinado escenario, χ(r) = r(r)−b(r) b(r)(5.1) donde b(r)es la permitividad relativa compleja del "background", de forma que, χ(r)=0∀r /∈ D. Como se comentó en el Capítulo 3, el Dominio Imagen Des iluminado por uno de los Tx de los Ntx trans- 49 50 CAPÍTULO 5. RESOLUCIÓN DE PROBLEMAS INVERSOS DE DISPERSIÓN ELECTROMAGNÉTICA misores/receptores existentes, generando el campo incidente Einc t, modelado por la siguiente ecuación escalar de Helmholtz, ∇2Einc t(r) + k2 bEinc t(r) = jwµoJt(r)(5.2) donde kb(r) = wpµoob(r)es el vector de ondas del "background". La siguiente medición se basa en introducir el OI en el recinto cerrado, generando el campo eléctrico total Et, para cada transmisor. Satisface la siguiente ecuación escalar, ∇2Et(r) + k2Et(r) = jwµoJt(r)(5.3) En este caso, k(r) = wpµoor(r). Entoces podemos definir el campo disperso como Esct t(r)≡ Et(r)−Einc t(r), referido a ∇2Esct t(r) + k2 bEsct t(r) = −k2 b(r)wt(r)(5.4) donde la fuente de constraste, wt(r) = χ(r)Et(r), describe la variable de contraste en términos del campo eléctrico total generado por las fuentes. La siguiente ecuación matricial relaciona las matrices FEM referidas a cada uno de los nodos que forman la malla con el campo disperso, [S−Tb]Esct t,Ω=Tbwt,Ω(5.5) mientras S∈CN×Nes la matriz de Stiffness que depende de las BCs, Tb∈CN×N, la matriz de Mass que define las propiedades del "background", los vectores Esct t,Ω∈CNywt,Ω∈CNson, respectivamente, los valores de dispersión EM y de las fuentes de contraste para cada Tx en el dominio Ω. Hay que aclarar que, debido al funcionamiento del FEM-CSIM, fue necesario modificar la implementación de la matriz FEM, ya que el método requiere la definición de las propiedades diléctricas nodo a nodo en lugar de elemento a elemento (Documento original [A.1] apéndice [B.1]. 5.2 Operadores Matriciales de Inversión Para la correcta definición del FEM-CSIM es necario introducir algunos operadores matriciales. El primer operador es denotado como Ms∈CR×N, donde R es el número de receptores. Este operador nos permite transformar los N valores nodales del campo disperso relacionados con el dominio Ωcon los 5.3 The Contrast Source Inversion Method 51 R fuentes puntuales. El segundo operador es representado como MD∈CI×N, donde I es el número de nodos definidos en el Dominio Imagen Dy selecciona los valores de campo de los I nodos del dominio imagen. De esta forma podemos obtener el vector fuente de constraste para toda la malla como wt,Ω=MT Dwt(5.6) Usando (5.6) en la ecuación matricial (5.5), obtenemos un nuevo operador, L ∈ CN×I, dado como Esct t,Ω=L[wt]=[S−Tb]−1TbMT D[wt](5.7) 5.3 The Contrast Source Inversion Method CSIM reliza la actualización de las variables de contraste y de las fuentes de contraste, minimizando la siguiente función de coste, F(χ, wt) = FS(wt) + FD(χ, wt) =Ptkft− MsL[wt]k2 S Ptkftk2 S +Pt χEinc t−wt+χ MDL[wt]  2 D Pt χEinc t  2 D (5.8) Aquí, ft∈CRes el campo disperso medido en las R posiciones de las fuentes para cada transmisor. Recuérdese que necesitamos convertir los datos medidos para todo el dominio Ω, a los equivalentes en el Dominio Imagen Dusando Einc t=MDEinc t,Ω. Como primer paso, en el algoritmo CSI actualizamos las fuentes de contraste wtcon el método del gradiente conjugado (CG) con las direcciones de búsqueda Polak Ribière (formulación desarrollada en [9]), asumiendo constantes las variables de contraste χ. En el segundo paso, wtse considera constante, y la función de coste FD(χ, wt)es minimizada. Ambas operaciones se realizan secuencialmente en cada iteración. wt,n =wt,n−1+αt,ndt,n (5.9) donde αt,n es el paso de actualización para el transmisor t en la iteración n, dt,n son las direcciones Polak Ribière. El paso viene definido como [9] αt,n =ηShρt,n−1,MsL[dt,n]iS+ηD nhrt,n−1, dt,n −χn−1 MDL[dt,n]iD ηSkMsL[dt,n]k2 S+ηD nkdt,n −χn−1 MDL[dt,n]k2 D (5.10) donde ηSyηD nson los factores de normalización, 52 CAPÍTULO 5. RESOLUCIÓN DE PROBLEMAS INVERSOS DE DISPERSIÓN ELECTROMAGNÉTICA ηS= X t kftk2 S!−1 ηD n= X t χEinc t  2 D!−1(5.11) y los términos de error ρt,n−1yrt,n−1son ρt,n−1=ft− MsL[wt,n−1] rt,n−1=χEinc t−wt,n−1+χn−1 MDL[wt,n−1] (5.12) Las direcciones de búsqueda de Polak Ribière dt,n son calculadas por la siguiente fórmula [9], dt,n =−gt,n +hgt,n, gt,n −gt,n−1iD kgt,n−1k2 D dt,n−1(5.13) donde dt,0es fijada a cero y gt,n es el gradiente de la función de coste F(χ, wt)respecto a wt,n, y puede ser aproximado por, (formulación extraída de [9]) gt,n =−2ηST−1 DLHMH sρt,n−1−2ηD nT−1 DI− LHMT DXH n−1TDrt,n−1(5.14) Donde I∈RI×Ies una matriz de identidad y Xn−1es la matriz diagonal cuyos elementos de la diagonal son los valores de χn−1. Después de actualizar las fuentes de contraste, χes evaluada minimizando FD(χ, wt)donde wtes considerada constante. χnpuede obtenerse resolviendo la siguiente igualdad matricial, X t EH t,nTDEt,n!χ=X t EH t,nTDwt,n (5.15) donde Et,n ∈CI×Ies la matriz diagonal definida por el vector Et,n =Einc t+MDL[wt,n]. 5.3.1 Normas y Productos Internos En este sección definimos algunos operadores y símbolos que aparecen durante la formulación anterior. La norma-L2 y el producto interior son calculados como, kxk2 D=xHTDxand hx, yiD=yHTDx(5.16) 5.4 Simulaciones 53 siendo x e y vectores de tamaño I, y TD∈RI×Ies la matriz de Mass para los nodos pertenecientes al Dominio Imagen D, definida como TD=ZD λi·λjds (5.17) Para la norma-L2 y el producto interior en Ω, con vectores de tamaño R, se define como kxk2 S=xHxand hx, yiS=yHx(5.18) 5.3.2 Estimación Inicial de FEM-CSIM No podemos incializar el método iterativo fijando la iteración inicial a cero porque la función de coste quedaría indefinida tras la primera iteración. De [3] y[9] extraímos la formulación óptima para la estimación inicial. La expresión (5.19) se obtiene aplicando el método de "steepest descent" a FS(wt). wt,0=Re MsL[GSft], ftS kMsL[GSft]k2 S GSft(5.19) donde el operador GS=−2ηST−1 DLHMH s. 5.4 Simulaciones En esta sección se presentan brevemente algunos resultados de simulaciones basadas en FEM-CSIM. Dado un determinado perfil, se utilizan los resultados electromagnéticos obtenidos tras la ejecucion de FEM para procesar el método de inversión. Se analizan una serie de experimentos donde determinaremos algunas conclusiones sobre la problemática de la aplicación. Para más detalles de cada simulación y experimento leer el documento original [A.1]. (a) (b) Figura 5.1: (a) Parte real permitividad relativa 0 r(b) Conductividad σ 54 CAPÍTULO 5. RESOLUCIÓN DE PROBLEMAS INVERSOS DE DISPERSIÓN ELECTROMAGNÉTICA Comenzaremos presentando un escenario simple basado en un OI circular de r= 0.02m con r= 5.0000 −j1.2839. El OI se encuentra iluminado por 16 Txs a fo= 700MHz equiespaciados con un r= 0.3m, mientras que el PEC circular tiene un radio igual a 0.36m, El "background" elegido está caracterizado por b= 3.0000 −j1.2839. Se muestra el escenario en la Fig.5.1. Para visualizar más fácilmente la calidad de las reconstrucciones obtenidas, calculamos el error entre el perfil reconstruido y el original. El error se define en términos relativos como, Err =kχexact −χreconstkD kχexactkD (5.20) Para comenzar el método iterativo hacemos uso del método de "backpropagation" [9] con el propósito de obtener la éstimación inicial. En la siguiente figura podemos observar χt,0y el correspondiente perfil dieléctrico. (a) (b) Figura 5.2: Estimación Inicial: (a) Prop.Dieléctricas (0 r,σ) (b) Contraste Inicial χt,0 En las siguientes figuras se presentan los resultados de la reconstrucción después de 400 iteraciones. Podemos observar como las reconstrucciones son bastante cercanas al perfil exacto. El OI es localizado e identificado. Los mayores problemas se visualizan a lo largo de las interfaces entre regiones debido a la asignación de parámetros físicos entre nodos pertenecientes a ambas regiones. En la Fig.5.4 se puede observar cómo la función de coste es minimizada hasta aproximadamente 10−4. Considerando que el resultado óptimo podría encontrase en torno al 10−5, el algoritmo ha funcionado correctamente. Se alcanza un error del 45.2%. La función FS(wt), hace referencia al error normalizado de los datos electromagnéticos mientras que FD(χ, wt), indica el error normalizado correspondiente a la 5.4 Simulaciones 55 reconstrucción del Dominio Imagen D. (a) (b) (c) Figura 5.3: (a) Reconst.prop.dieléctricas (b) True Re(χ)and Im(χ)(c) Reconst. Re(χ)and Im(χ) Aunque el error obtenido puede considerarse alto, para este tipo de aplicaciones es suficiente. El propósito de la imagen médica es el de detectar cuerpos extraños o posibles patologías. No se busca obtener resultados con la máxima definición, si no determinar la posición y naturaleza de las distintas regiones involucradas. Figura 5.4: Total, state and domain functional cost; Err: Error vs Iterations En el documento original [A.1] se realiza un experimento con este mismo escenario, tratando de determinar las limitaciones del algoritmo [13]. El siguiente experimento trata de definir la resolución de FEM-CSIM. Para ello, simulamos un escenario donde hay presentes dos OIs circulares de permitividad relativa compleja, r= 10 −j1.7975 a f0= 500MHz, separados por una distancia, d= 0.09m. Se introduce en un recinto circular PEC, de radio igual a 0.31m, con un "background" de bk = 6 −j2.876. Posteriormente, reducimos la distancia hasta d= 0.02m. Comparemos resultados. 56 CAPÍTULO 5. RESOLUCIÓN DE PROBLEMAS INVERSOS DE DISPERSIÓN ELECTROMAGNÉTICA (a) (b) Figura 5.5: d= 0.09m (a) Parte real permitividad relativa 0 r(b) Conductividad σ (a) (b) (c) Figura 5.6: d= 0.09m ((a) Reconst.prop.dieléctricas (b) True Re(χ)and Im(χ)(c) Reconst. Re(χ) and Im(χ) (a) (b) Figura 5.7: d= 0.02m (a) Reconst. diel. properties (b) Reconst. Re(χ)and Im(χ) 5.4 Simulaciones 57 Miestras que para una distancia mayor (Fig.5.9), ambos OIs son identificados con éxito, para la segunda simulación (Fig.5.10), se observa un solapamiento entre las superficies. Este hecho ocurre debido a que la segunda distancia es inferior a la resolución de la aplicación, definida como la mínima distancia entre dos posibles objetivos para poder ser identificados sin solapamiento. La distancia mínima viene dada por dmin =λg/4. Consideremos un escenario basado en un OI de características diléctricas inhomogéneas. Existe un primer círculo de radio igual a 0.25m con permitividad relativa, r= 12 −j0.674. Dentro de este cuerpo circular hay presentes dos más pequeños de permitividad relativa, r= 16 −j1.573 a 800MHz. El PEC circular de radio 0.312m contiene un "background" definido por b= 8 −j2.246. Ejecutamos el CSIM durante 850 iteraciones, usando 16 antenas localizadas a lo largo de un círculo de 0.25m. (a) (b) Figura 5.8: d= 0.09m (a) Parte real permitividad relativa 0 r(b) Conductividad σ (a) (b) (c) Figura 5.9: d= 0.06m (a) Reconst.prop.dieléctricas (b) True Re(χ)and Im(χ)(c) Reconst. Re(χ)and Im(χ) 64 CAPÍTULO 6. RESULTADOS Y CONCLUSIONES (a) (b) (c) Figura 6.5: |χ|>1: (a) Reconst.prop.dieléctricas (b) True Re(χ)and Im(χ)(c) Reconst. Re(χ)and Im(χ). Nótese que la elección del medio es fundamental para una correcta inversión, debiéndose adaptar al objeto de estudio. Para visualizar las distintas funciones de error y coste puede consultarse el documento original [A.1]. 6.1.2 Frecuencia Se deben considerar otros parámetros de simulación que implican consecuencias en los resultados, como la frecuencia de trabajo. Como ya se vio en el Capítulo 4, una frecuencia suficientemente grande puede acarrear una costosa discretización en términos de computación debido a la presencia de un mayor número de variables nodales. Sin embargo, tratamos de determinar las consecuencias que afectan a la calidad de la reconstrucción al usar una frecuencia u otra. Para ello, se ejecutan 3 simulaciones diferentes. Primero, un determinado dominio imagen es introducido en dos recintos circulares de diferente dimensionado, pero igual "background", b= 8 −j0.719. El primero de los PEC presenta un r= 0.21m, mientras que el segundo, r= 0.28m. Las propiedades dieléctricas del OI quedan definidas por r= 11.6−j0.539 afo= 1GHz. Nótese las diferencias del radio para cada uno de los dos casos. Si se denota el radio R en términos de la longitudes de onda, la primera computación es igual a 2λbk mientras que en la segunda es de 2.75λbk. 6.1 Limítaciones del CSIM 65 (a) (b) (c) Figura 6.6: (R/λbk) análisis (a) Prop.Dieléctricas (0 r,σ) (b) Resultados Reconst. r= 0.21m (c) Resultados Reconst. r= 0.28m. En la Fig.6.6(a) el perfil original es comparado con los resultados de la reconstrucción para ambos casos. De la misma manera, en la Fig.6.7 se puede contemplar la variable contraste original junto a las reconstrucciones realizadas para ambas simulaciones. (a) (b) (c) Figura 6.7: (R/λbk) análisis (a) True χ(b) Reconst.χr= 0.21m (c) Reconst.χr= 0.28m Finalmente, se realiza la tercera simulación, consistente en el recinto circular de r= 0.21m donde se introduce el mismo perfil dieléctrico que en la Fig.6.6(a) pero a una frecuencia menor, fo= 500MHz. De este modo, las permitividades relativas presentan un valor de r= 11.6−j1.0785 yb= 8 −j1.438 para el "background". 66 CAPÍTULO 6. RESULTADOS Y CONCLUSIONES (a) (b) (c) Figura 6.8: fo= 500MHz: (a) Reconst. prop. dieléctricas (b) True Re(χ)and Im(χ)(c) Reconst. Re(χ)and Im(χ) En la tabla 6.1 se muestran los términos de error para cada simulación después de 200 iteraciones. Se puede observar como la tercera computación genera los peores resultados mientras que la segunda presenta la mejor calidad de reconstrucción. En el caso de la primera y segunda simulación se está utilizando un medio con las mismas pérdidas. Sin embargo, la relación R/λbk es mayor para el segundo caso, lo que implica una mayor atenuación en las ondas reflectadas en el PEC, y en consecuencia, presenta una mejor calidad. En cuanto al tercer caso, a pesar de que presenta la mayor tangente de pérdidas, R/λbk es demasiado baja, lo que implica una baja atenuación de las reflexiones, provocando fallos durante la resolución del Problema Directo. La conclusión que se puede extraer de las simulaciones descritas en esta sección es que para realizar una correcta medición durante la realización del Problema Directo, es necesario establecer un compromiso entre la tangente de pérdidas que define a nuestro medio adaptador y las dimensiones del recinto para permitir una atenuación lo suficientemente grande, y de esta forma, evitar las reflexiones que representan ruido en nuestra aplicación. Tabla 6.1: Comparación entre las 3 simulaciones fo[GHz] |χ|R_PEC [m] R/λbk 00 b 0 berr(%) 1 0.448 0.21 2 0.089 43.8 1 0.448 0.28 2.75 0.089 42 0.5 0.445 0.21 1 0.179 55 6.2 Recintos PEC 67 6.2 Recintos PEC En esta sección se pretende determinar la mejor opción posible para el recinto conductor. Hasta el momento, hemos utilizado siempre PEC circulares; ahora analizaremos un mismo perfil para tres formas distintas del PEC. Imaginemos que se dispone de tres cilindros de dimensiones iguales pero fabricados con diferentes materiales. Sin embargo, desconocemos la identidad de cada uno de ellos. Para determinar la identidad del material de cada cilindro, se realiza el estudio de la sección transversal de los tres OIs, localizándolos equiespaciadamente. (a) (b) Figura 6.9: Estudio del PEC: (a) Prop.Dieléctricas (b) Contraste Verdadero (a) (b) (c) Figura 6.10: Circular: (a) Escenario (b) reconst.prop.dieléc. (c) reconst.contraste 68 CAPÍTULO 6. RESULTADOS Y CONCLUSIONES En la Fig.6.9 se describen las propiedades dieléctricas y la variable contraste del Dominio Imagen D. Las permititividades relativas de cada unade las regiones son: sc1= 14 −j0.8988;sc2= 12 −j1.0913; sc3= 16 −j1.4123 afo= 1.4GHz. (a) (b) (c) Figura 6.11: Cuadrada: (a) Escenario (b) reconst.prop.dieléc. (c) reconst.contraste (a) (b) (c) Figura 6.12: Triangular: (a) Escenario (b) reconst.prop.dieléc. (c) reconst.contraste Como se comentó anteriormente, el dominio imagen es introducido en tres recintos con formas distintas: un círculo de radio 0.12m Fig.6.10(a), un cuadrado de lado 0.24m Fig.6.11(a) y un triángulo equilatero de lado 0.42m Fig.6.12(a). Para todas las simulaciones, se utilizó un "background" con pérdidas (b= 9 −j2.568 a 1.4GHz). El OI fue iluminado por 24 antenas emplazadas a un radio de 0.1m. El dominio de inversión se basa en un cuadrado de lado 0.15m. En los tres casos el número de nodos era entre 5000 y 6000 nodos. En la Tabla 6.2 se muestran los resultados numéricos tras 550 iteraciones. 6.3 Elección del "Background" 69 Tabla 6.2: Comparación entre PECs Simulación err(%)log10(F)log10(Fs)log10(Fd) Circular 51 -3.61 -5.5 -1.8 Cuadrada 52 -3.375 -5.4 -1.7 Tirangular 52.5 -3.65 -5.7 -1.5 Se puede observar como los resultados de reconstrucción son similares para todas las simulaciones, de forma que podemos entender que no exiten consecuencias importantes en utilizar un tipo de recinto u otro, siempre y cuando las áreas en términos de longitud de onda seán más o menos equivalentes. 6.3 Elección del "Background" Hemos comprobado como la elección del medio donde el OI es sumergido requiere la consideración de las limitaciones del contraste. Igualmente, el "background" implica un mejor o peor funcionamiento de FEM-CSIM en cuanto a las pérdidas que este medio presenta, esencial en la atenuación de las reflexiones. 6.3.1 Medios Con Pérdidas Presentamos un escenario inhomogéneo basado en un OI cuadrado de lado 0.07m dentro de otro más grande de lado, l= 0.12m. El más pequeño tiene una permitividad compleja igual a r= 12 −j0.449, mientras la permitividad del grande, r= 8 −j1.123. Introducimos el dominio imagen dentro de un "background" definido por, b= 6−j1.7975 a 800MHz. A continuación se muestra el perfil original junto a las correspondientes reconstrucciones y estimaciones iniciales. (a) (b) (c) Figura 6.13: Cuadrados: a) True prop.dieléc. (b) Reconst. prop.diléc. (c) Estimación Inicial. 70 CAPÍTULO 6. RESULTADOS Y CONCLUSIONES (a) (b) (c) Figura 6.14: Cuadrados: a) True χ(b) Reconst. χ(c) Estimación Inicial Los resultados obtenidos tras 220 iteraciones son óptimos, con un error del 25%, lo que significa que el error ha decrecido rapidamente desde la iteración inicial. Se ha podido minizar la función hasta 10−4.5, un valor muy cercano a la minimización óptima, 10−5. En la Fig.6.15 se muestran los resultados en términos de error. Figura 6.15: Cuadrados: Total, state and domain functional cost; Err: Error vs Iterations La siguiente simulación se basa en un perfil muy popular en imagen conocido como E-phantom. En este caso suponemos un medio con bajas pérdidas con una permitividad relativa, b= 16 −j1.8 afo= 1.2GHz. El E-phantom esta basado en superficies inhomogéneas donde, la mayor parte de la superficie presenta r= 33−j2, mientras que la pequeña inclusión circular y el rasgo de más a la derecha quedan definidos por r= 33 −j8.33. El radio elegido para el PEC es igual a 0.17m. En las imágenes expuestas a continuación se muestran las propiedades dieléctricas del dominio imagen, y junto a ellas, se presentan los correspondientes resultados de reconstrucción. Del mismo modo, podemos observar tanto los datos como los resultados obtenidos para la variable contraste en la Fig.6.17. 6.3 Elección del "Background" 71 (a) (b) (c) Figura 6.16: E-phantom: a) True prop.dieléc. (b) Reconst. prop.diléc. (c) Estimación Inicial. (a) (b) (c) Figura 6.17: E-phantom: a) True χ(b) Reconst. χ(c) Estimación Inicial Es posible identificar cada una de las superficies que constituyen el OI, lo que indica alta calidad en el proceso de inversión. Sin embargo, debido a la resolución del algoritmo, en la zona donde las líneas se encuentran más próximas, los resultados nos muestran una zona difusa donde las partes no se identifican correctamente. El error obtenido tras 751 iteraciones, estimación inicial incluida, es del 38.1%. Si se analiza las funciones de estado y datos, se puede apreciar que los datos sintéticos presentan poco ruido, lo que ha permitido obtener un buen resultado de la reconstrucción en comparación con los datos reales. 72 CAPÍTULO 6. RESULTADOS Y CONCLUSIONES Figura 6.18: E-phantom: Total, state and domain functional cost; Err: Error vs Iterations 6.3.2 Medios Sin Pérdidas A continuación se presentan los resultados de inversión usando medios sin pérdidas, donde se observa cómo la calidad de la reconstrucción es claramente inferior a los resultados equivalentes con medios con pérdidas. La razón reside en la ausencia de pérdidas durante la propagación de las ondas reflejadas en las paredes, lo cual provoca que aparezcan problemas en la sintetización de los datos simulados tras el Problema Directo debido a las reflexiones que iluminan el dominio imagen. Este hecho genera soluciones no reales; es decir, FEM-CSIM genera reconstrucciones de objetivos que no representan el auténtico OI. Esta situación podría entenderse como otra limitación de la aplicación en escenarios simulados. (a) (b) (c) Figura 6.19: Medios sin pérdidas I: a) True prop.dieléc. (b) Reconst. prop.diléc. (c) Estimación Inicial. Ahora, se simula el mismo perfil dieléctrico que en Fig.6.13(a) dentro de un "background" sin pérdidas, b= 6 a 800MHz. Recordad que el OI inhomogéneo se basaba en dos superficies cuadradas de permitividad, respectivamente, r= 12 −j0.449 yr= 8 −j1.123. 6.3 Elección del "Background" 73 Podemos observar como los resultados de la Fig.6.19, obtenidos presentan valores numéricos cercanos a los originales pero las localizaciones de las regiones y las formas no están bien definidas. (a) (b) (c) Figura 6.20: Medios sin pérdidas I: a) True χ(b) Reconst. χ(c) Estimación Inicial Analicemos ahora los resultados asociados a la variable constraste. Estos resultan ser similares a los anteriores, apareciendo extrañas formas que no permiten identificar correctamente la geometría del dominio imagen. Esos errores se deben a los problemas asociados a las reflexiones y soluciones no reales. Figura 6.21: Medios sin pérdidas I: Total, state and domain functional cost; Err: Error vs Iterations El error en la variable de contraste es del 49%. Observese la escasa minimización realizada por el algoritmo en las funciones de coste en comparación con su equivalente con pérdidas en el medio. Vemos como la inversión obtenida en la función de datos ha sido demasiado pobre. Realicemos de nuevo una comparación respecto a los "backgrounds" utilizados, basada en la misma geometría del E-phantom definida en la Fig.6.16(a). Consideremos un OI homogéneo con una permitividad relativa igual a r= 5.9−j1.7975 , sumergido en un "background" de permitividad b= 4.2a fo= 800MHz. 80 CAPÍTULO 6. RESULTADOS Y CONCLUSIONES Bibliografía [1] Jianming Jin, The Finite Element Method in Electromagnetics. A Wiley-Interscience Publication, John wiley & Sons, Inc. 2nd Edition, 2002. [2] A.Bondenson, T.Rylander, P.Ingelström, Computational Electromagnetics. Springer. June 27,2005. [3] Amer Zakaria, The Finite-Element Contrast Source Inversion Method for Microwave Imaging Applications. A Thesis submitted to the Faculty of Graduate Studies of The University of Manitoba. 2012. [4] Anastasis C.Polycarpou, Introduction to the Finite Element Method in Elecromagnetics. Morgan & Claypool Publishers. First Edition, 2006. [5] Jean-Pierre Berenger, Perfectly Matched Layer (PML) for Computational Electromagnetics. Morgan & Claypool Publishers. First Edition, 2007. [6] Elia A.Attardo, Computational Methods for Microwave Imaging: Biomedical Applications. PhD Thesis, Politecnico di Torino. May 2011. [7] Jiayuan Fang, Zhonghua Wu, Generalized Perfectly Matched Layer-An Extension of Berenger’s Perfectly Matched Layer Boundary Condition. IEEE Microwave And Guided Wave Letters, VOL. 5, NO. 12. December 1995. [8] David M.Pozar, Microwave Engineering. John wiley & Sons, Inc. Second Edition, 1998. [9] A.Zakaria, C.Gilmore and J.LoVetri, Finite-Element Contrast Source Inversion Method for Microwave Imaging. Article IOPscience. September 2010. 81 82 BIBLIOGRAFÍA [10] C.Gilmore and J.LoVetri, Enhancement of Microwave Tomography Through The Use of Electrically Conducting Enclosures. Article IOPscience. April 2008. [11] International Center for Numerical Methods in Engineering, http://gid.cimne.upc.es/support/ manuals CINME,UPC. [12] Microwave Imaging Laboratory Website, http://www.ece.umanitoba.ca/EM_Imaging_Lab/index. html University of Manitoba. [13] M.D’Urso, T.Isernia, and Andrea F.Morabito, On the Solution of 2-D Inverse Scattering Problems via Source-Type Integral Equations. IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL.48, NO.3, March 2010. [14] A.Zakaria and J.LoVetri, Application of Multiplicative Regularization to the Finite-Element Contrast Source Inversion Method. IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 59, NO. 9, September 2011. [15] C.J.Reddy, Manohar D.Deshpande, C.R.Cockrell, and Fred B.Beck Finite Element Method for Eigenvalue Problems in Electromagnetics. NASA Technical Paper 3485, December 1994. [16] Jianming Jin and W.C. Chew Combining PML and ABC for Finite Element Analysis of Scattering Problems. Article of University of Illinois at Urbana-Champaign, [17] J.Fang and Z.Wu Generalized Perfectly Matched Layer for the Absorption of Propagating and Evanescent Waves in Lossless and Lossy Media. IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO.12 December 2012 [18] P.Mojabi and J.LoVetri Eigenfunction Contrast Source Inversion for Circular Metallic Enclosures. IOP PUBLISHING. INVERSE PROBLEMS. 12 January 2010 Apéndice 83 Ap´ endice A Documento original de la tesis 85 A Finite Element Approach To Inverse Scattering Problems MASTER’S DEGREE PROJECT POLITECNICO DI TORINO TELECOMMUNICATIONS ENGINEERING AUTOR: ENRIQUE MASGRAU RITE TUTOR: GIUSEPPE VECCHI November 22, 2012 Acknowledgements During this last year I have lived one of the best vital experiences of my life. Torino has been my home for all this time, and from where I will take amazing friends forever. In this country, I have felt like one more, and I have received a great affection by all its people. I really will have Torino inside me when this adventure will finish. Furthermore, Torino is the place where I decided to finish my university studies, including this work. I want to thank Politecnico di Torino for giving me this possibility. Study abroad has made me grow as a person. Now, I want to thank all the people that support me along this time to perform this research work. Thank to Istituto Superiore Mario Boella, the place where I have been doing my master’s project. In the ISMB I’ve been able to enjoy a friendly atmosphere of work. I would like to sincerely show my gratitude to professor Giuseppe Vecchi for letting me make the thesis with him and for Elia A.Attardo for all his support, effort and patience during these months. In the same way, I would like to show gratitude to my family. For years of guidance and support. Without them, this experience would not have been possible. Finally, I would like to thank my friends from Spain, Italy and Erasmus in general that have helped me to enjoy the work, the studies and good moments. Grazie a tutti, davvero. 3 10 Lista de Figuras 1.1 MWI application results ................................ 19 1.2 Process description ................................... 20 2.1 Interface between two media. ............................. 26 3.1 Relative permittivity (blue line) and conductivity (red line) for glycerine-water 80:20 from 300MHz up to 3GHz ............................ 30 3.2 The enclosed imaging system. ............................. 30 3.3 (a) MWT prototype system (b) Tomography reading ................ 31 3.4 (a) 2D MWI model (b) 3D MWI model ..................... 32 4.1 Different elements to discretize a certain 1D, 2D or 3D geometry ......... 37 4.2 Discretization error obtained using triangular or quadrilateral elements . . . . . . 38 4.3 2D geometry based on two different regions ..................... 38 4.4 The numbering of local nodes for the element e. ................... 40 4.5 Triangular element divided into three sub-triangles. ................. 41 4.6 Triangular edge element. ................................ 43 4.7 (a) Interior edge shared by two neighboring triangles. (b) The outward unit vectors normal to the common edge point in opposite directions .............. 49 4.8 Fictitious Boundary .................................. 52 4.9 Incident, reflected and trasmitted waves ....................... 56 4.10 The PML ABC on a plane boundary ......................... 57 4.11 Concave surface enclosing the domain ........................ 58 11 4.12 2D meshing for 0 r= 20................................. 62 4.13 Forward Problem solved by FEM using nodal functions: 0 r= 20, (Xs,Ys)=(0.5,0), (a)(b) |E|and ∠Efield for σ= 0.0005 S/m (c)(d) |E|and ∠Efield for σ= 0.5 S/m. ........................................... 63 4.14 Forward Problem solved by FEM using nodal functions: 0 r= 20 (a)(b) |E|and ∠Efield for σ= 0.0005 [S/m] and for current sources (Xs,Ys)=([0.05 0],[0 -0.05]), (c)(d) |E|and ∠Efield for σ= 0.0005 [S/m] and for current sources (Xs,Ys)=([- 0.05 0.05],[-0.05 0.05]). ................................. 64 4.15 Forward Problem solved by FEM using nodal functions: 0 r= 80 (a)(b) |E|and ∠Efield for σ= 0.0005 [S/m] and for current sources (Xs,Ys)=(0.05,0), (c)(d) |E| field for σ= 0.0005 [S/m], current sources (Xs,Ys)=(0.02,0) and (Xs,Ys)=(0.07,0). 65 4.16 (a)(b) Inhomogeneous 2D Scenarios (c)(d) |E|fields for current source at (Xs,Ys)=(0,- 0.05). ........................................... 66 4.17 Forward Problem solved by FEM using nodal functions: 0 r= 10,σ= 0.03 [S/m], (a) |E|field for fo= 700MHz, (b) |E|field for fo= 1GHz, (c) |E|field for fo= 1.5GHz, (d) |E|field for fo= 2GHz ...................... 68 4.18 (a)(b) |E|field for fo= 700MHz and fo= 2GHz respectively. .......... 68 4.19 Triangular 2D: 0 r= 20, (a)(b) |E|and ∠Efor σ= 0.02 [S/m] and (Xs,Ys)=(0,0), (c)(d) |E|and ∠Efor σ= 0.5[S/m] and (Xs,Ys)=(0,0) ............... 69 4.20 (a) One current source (b) Three current sources ................ 70 5.1 (a) Real relative permittivity 0 r(b) Conductivity σ................ 79 5.2 Initial Guess: (a) Diel.Properties (0 r,σ) (b) Initial Constrast χt,0......... 80 5.3 (a) Reconst. diel. properties (b) True Re(χ)and Im(χ)(c) Reconst. Re(χ) and Im(χ)........................................ 81 5.4 (a) Total, state and domain functional cost (b) Error vs Iterations ...... 82 5.5 |χ|>1: (a) True Re(χ) and Im(χ) (b) Reconstructed Re(χ) and Im(χ)..... 84 5.6 Re(χ)<0:(a) True Re(χ) and Im(χ) (b) Reconstructed Re(χ) and Im(χ). . . 85 5.7 d= 0.09m (a) Real relative permittivity 0 r(b) Conductivity σ.......... 86 5.8 d= 0.09m (a) Reconst. diel. properties (b) True Re(χ)and Im(χ)(c) Reconst. Re(χ)and Im(χ).................................... 86 12 5.9 d= 0.02m (a) Reconst. diel. properties (b) Reconst. Re(χ)and Im(χ)...... 87 5.10 Inhomogeneous Profile (a) Real relative permittivity 0 r(b) Conductivity σ. . . 87 5.11 Inhomogeneous Profile (a) Reconst. diel. properties (b) True Re(χ)and Im(χ) (c) Reconst. Re(χ)and Im(χ)............................. 88 5.12 Total, state and domain functional cost; Err: Error vs Iterations ......... 88 5.13 Initial Guess: (a) Diel.Properties (0 r,σ) (b) Initial Constrast χt,0......... 89 5.14 U-Shape profile: (a) Real relative permittivity 0 r(b) Conductivity σ...... 90 5.15 U-Shape profile: (a) Reconst. diel. properties (b) True Re(χ)and Im(χ)(c) Reconst. Re(χ)and Im(χ)............................... 90 5.16 U-Shape profile: Total, state and domain functional cost; Err: Error vs Iterations 91 6.1 |χ|<1: (a) Real relative permittivity 0 r(b) Conductivity σ........... 94 6.2 |χ|<1: (a) Reconst. diel. properties (b) True Re(χ)and Im(χ)(c) Reconst. Re(χ)and Im(χ).................................... 95 6.3 re(χ)<0: (a) Reconst. diel. properties (b) True Re(χ)and Im(χ)(c) Reconst. Re(χ)and Im(χ).................................... 95 6.4 |χ|>1: (a) Real relative permittivity 0 r(b) Conductivity σ........... 96 6.5 |χ|>1: (a) Reconst. diel. properties (b) True Re(χ)and Im(χ)(c) Reconst. Re(χ)and Im(χ).................................... 96 6.6 |χ|<1(a) Total, state and domain functional cost (b) Error vs Iterations . . 97 6.7 re(χ)<0(a) Total, state and domain functional cost (b) Error vs Iterations . 98 6.8 (R/λbk) analysis (a) Diel.Properties (0 r,σ) (b) reconst.results r= 0.21m (c) reconst.results r= 0.28m................................. 99 6.9 (R/λbk) analysis (a) True χ(b) Reconst.χr= 0.21m (c) Reconst.χr= 0.28m. . 100 6.10 fo= 500MHz: (a) Reconst. diel. properties (b) True Re(χ)and Im(χ)(c) Reconst. Re(χ)and Im(χ)............................... 101 6.11 PEC study profile: (a) Dielec.Properties (b) True Contrast Source ........ 103 6.12 Circle: (a) Enclosure shape (b) reconst.dielectric prop. (c) reconst.contrast source103 6.13 Square: (a) Enclosure shape (b) reconst.dielectric prop. (c) reconst.contrast source .......................................... 104 13 6.14 Triangle: (a) Enclosure shape (b) reconst.dielectric prop. (c) reconst.contrast source .......................................... 104 6.15 Squares: a) True diel. properties (b) Reconst. diel. properties (c) Initial guess 106 6.16 Squares: a) True χ(b) Reconstructed χ(c) Initial guess ............ 107 6.17 Squares: Total, state and domain functional cost; Err: Error vs Iterations . . . . 107 6.18 E-phantom: a) True diel. properties (b) Reconst. diel. properties (c) Initial guess ........................................... 108 6.19 E-phantom: a) True χ(b) Reconstructed χ(c) Initial guess .......... 108 6.20 E-phantom: Total, state and domain functional cost; Err: Error vs Iterations . . 109 6.21 Lossless I: a) True diel. properties (b) Reconst. diel. properties (c) Initial guess 110 6.22 Lossless I: a) True χ(b) Reconstructed χ(c) Initial guess ........... 110 6.23 Lossless I: Total, state and domain functional cost; Err: Error vs Iterations . . . 111 6.24 Lossless II: (a) Reconst.Dielectric Properties (b) Reconst.Contrast Source . . . . 112 6.25 Lossless II: Total, state and domain functional cost; Err: Error vs Iterations . . . 112 6.26 Forearm: (a) Real relative permittivity 0 r(b) Conductivity σ.......... 113 6.27 Forearm: (a) True Contrast (b) Reconst.Contrast ................. 114 6.28 Forearm: Total, state and domain functional cost; Err: Error vs Iterations . . . . 114 6.29 Brain: (a) Dielectric properties (b) True Contrast ................. 115 6.30 Brain: (a) Reconst.Diel.Properties (b) Reconst.Contrast ............. 116 6.31 Brain: Total, state and domain functional cost; Err: Error vs Iterations . . . . . 116 A.1 (a) Linear triangular element in the xy-plane. (b) Linear triangular element (master element) in the ξη-plane .............................. 124 14 Lista de Tablas 4.1 Edge Numbering for a Triangular Element ...................... 42 4.2 Computational Costs .................................. 71 6.1 Comparision between 3 simulations .......................... 101 6.2 Comparision between PEC shapes .......................... 105 6.3 Forearm properties ................................... 113 6.4 Brain Model ....................................... 115 6.5 CSIM Results ...................................... 118 15 16 Chapter 1 Introduction In this research work a biomedical application is described and designed, to do this, a mathemathical and engineering method is implemented with the purpose of being evaluated to extract conclusions about the results. On the one hand, this method presents important electromagnetic concepts that are used to solve this medical scenario (among many other possibilities), and on the other hand engineering techniques are mixed with those concepts to optimize the method that is researched. 1.1 Scope Nowadays Microwave Imaging (MWI) is one of the most important research fields in biomedial imaging (tomography and breast detection are some biomedial applications based on MWI), it can be defined as discovering the internal structure of an object by illuminating it with electromagnetic methods. Inverse Scattering problem is a MWI application, the aim of this research work is based on describing the method to solve it in the correct way. MWI can be based on several methods according to the medical or engineering application, in this case the method works with electromagnetic concepts mixed with engineering and mathemathical techniques to know the dielectric properties of a certain target that is unknown but can be reconstructed, this target is known as Object of Interest (OI) and it will be described more deeply in Chapter 3 during the description of the Forward Problem. 17 18 CHAPTER 1. INTRODUCTION The Forward Problem is one of the two main steps that form this method, the target is radiated by electromagnetic devices (antennas) and the results are used like an input of the Inverse Scattering method, that is the second part of the process. Along this research work the inverse method and the way to compute the electromagnetic results (Finite Element Method) will be described deeply. 1.2 Motivation Forward Problem and Inverse method will be implemented and evaluated separately to obtain a global method that lets us to known the dielectric properties of a unknown object that can be placed inside of a known body, knowing these material properties is possible to understand what type of material is analyzing, this idea has a great benefit in biomedical engineering, more exactly in medical imaging, for example, the tomography gives us a way to detect breast tumors, problems in bones and other type of diseases. These problems are the reasons why this research work was choosen and correctly developed, MWI has an important impact in medical engineering, it is undeniable that such studies are necessary in the development of medicine, so that, engineering must research about medical applications and methods to improve it. Is proud to be part of the study of these bio-engineering applications so important in the real life. 1.3 Problem Definition This section presents a description of the scenario that is involved in this research work. Imagine that there is a body with certain unknown properties, these unknown properties are the target of a study or investigation, It’s necessary to understand the behaviour of this material to be able to determine some conclusions that are very important from the point of view of the investigation, for example a medical analysis to diagnosticate a possible proplem (breast cancer detection). Then, we have a target that part is known and a part is unknown, as mentioned in section before, this unknown body is called as Object of Interest, in the future we will always refer to this target as OI. In Fig.1. a possible medical scenario is shown, a breast cancer detection, in this case the material properties of the breast are known, so with the correct anaylis is possible to detect some strange 1.4 Proposed Solution 19 bodies in the breast that could be a little tumor. In final results of the diagnostic would be possible to observe a total surface or volume where different dielectric properties are involved and obviously some of them means a dangerous object. Other possibility could be an analysis of an object formed by two parts, firstly the known region called as background medium whose dielectric porperties are known, secondly the Imaging Domain that is formed by the OI and a part of the total background medium (known material) , the total region is the sum of both regions although the Imaging Domain is the aim of this MWI application. For the study of the OI is necessary to implement and design a method that lets us to obatin those properties of the unknown body such that the target must not be altered, much less damaged by the method. Figure 1.1: MWI application results 1.4 Proposed Solution There are different ways and methods to solve this problem, in this project the global method that performs this type of biomedical applications is known as FEM-CSIM (FEM-Contrast Source Imaging Method). FEM means Finite Element method, it will be one of the targets of study more important in this research work, this method lets to compute very efficiently electromagnetic parameters in a certain scenario, it is a popular method to analyze an object or target that is used in a lot of engineering fields such as structure designs, temperature analysis of objects... FEM will be described in Chapter 4 to understand how this method works and how is implemented. If FEM-CSIM is implemented correctly, it will provide the dielectric properties of the OI (object of interest) with a little error like in all engineering applications, later these results could be 26 CHAPTER 2. ELECTROMAGNETIC CONCEPTS Figure 2.1: Interface between two media. In expressions (2.16) and (2.17) both surface electric current density and surface electric charge density have been considered zero, in the case where this assumption is not considered the expresions would be, ˆn×(~ H1(~r)−~ H2(~r)) = ~ Js(2.19) ˆn×(~ D1(~r)−~ D2(~r)) = ρs(2.20) 2.3 Helmholt Equation In the application that is developed during this work is important to solve partial difference equations (PDE) where the ~ Efield is involved, remove ~ His necessary to obtain the next expression, obviously this assumption implies consequences in most of the equations commented until this moment. ∇×∇× ~ E(~r)−w2µo(~r)~ E(~r) = −jwµo~ Ji(~r)(2.21) Two conditions have been assumed, (i) there aren’t electric charges (ρ=0), (ii) the relative permeability is null (non-magnetic problems) µr= 1. Equation (2.21) is known as Helmholtz Equation, in particular, FEM is generally used to solve the frequency domain form of the curl-curl equation, sometimes referred to as the vector Helmholtz equation. 2.4 Scalar Wave Equations 27 2.4 Scalar Wave Equations In electromagnetic analysis, whenever possible problems are simplified by using a 2D model to approximate a 3D problem. Assume that the fields has not variation with respect to one Cartesian coordinate (z-coordinate), we treat a scalar scenario that is defined by ∇2Ez+k2 orEz=jwµoJz(2.22) also called as inhomogeneous scalar wave equation, this expression will be used in Chapter 4 to describe the 2D scalar problem given by TM and z-polarization case. 28 CHAPTER 2. ELECTROMAGNETIC CONCEPTS Chapter 3 The Forward Problem Microwave imaging (MWI) is of interest for various applications such as geophysical surveying and medical imaging. In the form of MWI considered herein, one attempts to quantitatively reconstruct the, mostly unknown, electrical properties (i.e. permittivity and/or conductivity) of an object of interest (OI) which is immersed in a background medium of known electrical properties [9]. As mentioned in Scope Section in Chapter 1, to solve a MWI application as is developed during this research work is necessary firstly implement the Forward Problem with the aim of obtain the electric distribution of the electromagnetic field along the OI, after the electromagnetic field is computed, the next step is based on the Inverse Algorithm that is presented in Chapter 5. This chapter represents a briefly introduction to the Forward Problem in the context of our biomedical application, the chosen method used to compute the electric field is totally defined for 2D geometries in Chapter 4. However, along this chapter we talk about this method, describing how it is executed in real experimental applications in order to understand the differences between experimental and synthetic. 3.1 Starting Point Before explaining how is implemented the Forward Problem based on simulating, we think that could be profitable to understand how MWI is performed in real medical applications such as breast detection or bone tissue analysis. To get this, a description of real experiments based on MWI is developed during this section. We consider microwave scenarios where the imaging region is surrounded by an electrically conducting surface, this surface is denoted as PEC. 29 30 CHAPTER 3. THE FORWARD PROBLEM In real MWI applications the OI is placed inside a tank, the conductive surface serves as both the container for any possible matching fluid and a shield from outside interference as commented before. Obviously, the inclusion of the conducting enclosure considerably changes the distribution of the EM energy as compared to an open region. The matching medium should be as optimal as possible to avoid possible wave reflection from the walls, in Fig.3.1 extracted from [6] we can observe a common used liquid medium based on 80:20 percent glycerine-water solution. Figure 3.1: Relative permittivity (blue line) and conductivity (red line) for glycerine-water 80:20 from 300MHz up to 3GHz Surrounding the OI inside the tank, a set of antennas, that work both Tx and Rx electromagnetic devices, illuminate the OI with the aim of storing the scattering electric field. The scatterer is located entirely within the region Dand is embedded in a homogenous background medium (matching medium) with a background permittivity whose dielectric properties are denoted as b. Figure 3.2: The enclosed imaging system. 3.2 Electromagnetic Waves 31 In Fig.3.2[10] a general scheme is shown. In this case the enclosure is based on a cylinder tank. Notice that during this project the dielectric properties are denoted as =0+j00 =or=o0 r+j00 r=o0 r+σ jwo(3.1) Now, some possible real MWI scenarios are shown in the next figure. (a) (b) Figure 3.3: (a) MWT prototype system (b) Tomography reading 3.2 Electromagnetic Waves When we talk about analyse the electromagnetic wave in the Forward Problem, we need to differentiate between all possible EM cases. In this section a general discussion of the different types of wave propagation that can exist is presented. Maybe in a medium the wave propagation could be based on tranverse electromagnetic waves characterized by Ez=Hz= 0, this implies that tranverse field components are zero. Other possible case is the TE waves, this propagation means that Ez= 0, so the unique transverse component is the magnetiz field Hz6= 0. Finally, there is a wave propagation case based on a non zero electric tranverse component, Ez6= 0, instead the magnetic tranverse field is null, Hz= 0. In microwave terminology, the 2D formulation corresponds to TM polarization. This 2D problem is readily solved using nodal elements and is the object of this performed electromagnetic study along the research work. TM scenarios are described by the Maxwell’s equations defined in Chapter 2 and precisely in Section 2.4. 32 CHAPTER 3. THE FORWARD PROBLEM 3.3 Simulation Methods The purpose of this research work is to simulate a real MWI experiment using a method that let us to compute the scattered electric field, there are different used techniques along the history to get it. One of them is called as the finite element method that will be described deeply in Chapter 4, it is the chosen method to solve the Forward Problem in our MWI application. Basically, this method implies to solve a matrix equation where the electric coefficients that define the electric field can be obtained generating with FEM, a global matrix called as FEM matrix that is related to the geometry and properties of our scenario (OI + background). This matrix equation will be detailed developed in Chapter 4 and is given by Az =b(3.2) Our forward simulation involve an iterative method where a given antenna from the set of antennas, that are surronding the OI, works as a trasnsmitter, i.e., in each iteration one different of the total antennas is the trasmitter while the others work as receivers, later that antenna will work as receiver while other antenna will be the new transmitter. Obviously, if the number of antennas used in the method is higher, the quantity of information will be wider. In Fig.3.4 (images from [3])a 2D and 3D approach is defined, observe how one of the antennas works as Tx while the others are Rx devices. (a) (b) Figure 3.4: (a) 2D MWI model (b) 3D MWI model 3.3 Simulation Methods 33 As we work with synthetic data, it is also necessary to generate a possible 2D or 3D geometry that describe a certain scenario as the defined in last section. There are a lot of simulators that let us to create several surfaces or volumes. In our case, the software is called as GiD, obviously it is based on the FEM (more information on GiD website http://gid.cimne.upc.es/). Finally, it is important to clarify that FEM represents a different way to compute the EM field than the traditional theory based on Green’s function that is called as the Intregral Formulation, in [9] a deeply comparision between both methods is performed. 34 CHAPTER 3. THE FORWARD PROBLEM Chapter 4 The Finite Element Method This chapter presents a finite element method modeling for Microwave Imaging tomography, due to the importance of FEM in this project it represents one of the main targets of study during this research work. We will talk about the basic concepts of this technique, possible applications for using this method, the requirements that the FEM design need to perform correct results, the elements used to build the possible meshing, the assembling matrix and obviously how FEM is used to solve the Forward Problem. Furthermore other problems are solved such as the boundary-value problems (BVP) in 2D geometries for scalar and vector problems (3D), and scenarios where a certain type of structures called as Absorbing Boundary Conditions like Perfect Matched Layers (PMLs) are considered in the design and computational simulation. The finite element method to solve the Forward Problem has been implemented using MatLab and GiD programmes that are respectively a very popular mathemathical simulator based on matrix processing and a powerful structure designer and simulator using the finite element method. 4.1 FEM basic concepts This section presents the basic principles of the finite element method (FEM), the first practical use for FEM began in the 1950s for aircraft design; however it’s first use in electrical engineering was not until 1965 when Winslow used it to solve for the magnetic field on an irregular mesh [3]. The finite element method is a standard tool for solving differential equations in many disciplines, e.g., electromagnetics, solid and structural mechanics,fluid dynamics, acous- 35 42 CHAPTER 4. THE FINITE ELEMENT METHOD Tabla 4.1: Edge Numbering for a Triangular Element Edge No.i Node i1Node i2 1 1 2 2 2 3 3 3 1 4.3.2 Edge Basis Functions For solving vectorial problems using FEM, the use of nodal-based elements exhibits shortcomings as spurious solutions, fortunately, a revolutionary approach was discovered. This approach uses so-called vector basis or vector elements that assign degrees of freedom to the edges rather than to the nodes of the elements. Edge elements eliminate spurious modes that can introduce errors in field calculations for near-field problems. In this section we introduce edge elements in two dimensions. Thanks to edge elements the electric field in each triangular element can be approximated by the next expression: ~ Ee(~r) = 3 X j=1 Ee j~ Ne j(~r)(4.11) where Ee jdenotes the tangential field along the jth edge. The nodes of the triangular element are joined together by three edges. Each edge in the mesh is identified with two labels: a local number to indicate its location in a given triangle and a global number to indicate its location with respect to the entire mesh. Each edge of the triangle is associated with a vector-basis function: ~ Ne j(~r) = le iϕe i1~r∇ϕe i2~r −ϕe i2~r∇ϕe i1~r(4.12) where le idenotes the length od the edge i, ϕe i1~r and ϕe i2~r are the nodal basis functions described in expression (4.4) of each node that forms edge e. 4.4 Scalar Problems 43 Figure 4.6: Triangular edge element. 4.4 Scalar Problems Once the finite element method has been described is necessary to introduce the method to solve possible problems, in Chapter 2 scalar and vector wave equations were commented. Now we develop FEM for two-dimensional scenarios in electromagnetics that are defined by scalar problems, basically a scalar problem is represented by a variable wich must be measured, for example, a TM z-polarization case. This section is focus on TM, z-polarization and two-dimensional case that are the main characteristics that define our stage in this research work. In Chapter 2 the boundary conditions in electromagnetic theory were also presented, here introducing these concepts in scalar equations the Boundary-Value Problem arises in our approach, a particular method called as Garlerkin’s method is deeply described to solve scalar problems and obviously how nodal basis functions are used in these applicatios. In conlusion, during this section FEM formulation based on processing matrix equations is used to solve scalar problems 4.4.1 2D Boundary-Value Problem In Section 2.4 a second-order PDE that defines scalar problems was introduced, expresion (2.22) was named as inhomogeneous scalar wave equation, now we analyse it for 2D geometries. It is important to notice that MWI applications are focused on PEC enclosure, this means that BC must be considered during the performed formulation and computation using the finite element method. Then some boundary conditions are involved in this second-order PDE to define 44 CHAPTER 4. THE FINITE ELEMENT METHOD correctly the scalar problem of our application, in 2D geometry this approach is known as 2D Boundary-Value Problem (BVP). These new conditions must be considered in order to implement with accuracy the formulation that will determine the assembling process, so the boundary-value problem under consideration is defined by the next second order differential equation −∇·(α∇u) + βu =g(4.13) and the boundary conditions u=p on Γ1(4.14) ˆn·(α∇u) + γu =q on Γ2(4.15) where u is the unknown function, αand βare the known parameters associated with the physical propoperties of the domain Ω, and g is the source or excitation function. According to boundary conditions, Γ1defines the Dirichlet boundary while Γ2refers to the Robin boundary, with Γ1+ Γ2= Γ,Γdenotes the contour enclosing the area Ω.γ, p and q are the knwon parameters associated with the physical propperties of the boundary, when γ= 0 the Robin BC is a special case called as Neumman boundary. As we explained in Section 2.4, in scalar problems we assume that the fields has not variation with respect to one Cartesian coordinate, for example z coordinate. Our problem depends only on one variable, this is the case of TM and z-polaraization wave, so the expression (4.13) can be written as given in (2.22). 4.4.2 Solving BVP via Garlerkin’s Method The generic 2D BVP considered in expression (4.13) in this section could be expressed by a second-order partial differential equation given by ∂ ∂x αx ∂u ∂x+∂ ∂y αy ∂u ∂y +βu =g(4.16) and the boundary conditions 4.4 Scalar Problems 45 u=p on Γ1(4.17) αx ∂u ∂x ·ˆx+αy ∂u ∂y ·ˆy·ˆn+γu =q on Γ2(4.18) where αx,αyand βare constants. The weak formulation of this problem can be obtained by first constructing the weighted residual of (4.16) for a single element with domain Ωe re=∂ ∂x αx ∂u ∂x+∂ ∂y αy ∂u ∂y +βu −g(4.19) This element residual is ideally zero, provided that the numerical solution u to be obtained is identical to the exact solution. However, this is not the case, and therefore, the element residual reis, in general, nonzero. Our objective is to minimize this element residual in a weighted sense. To achieve this, we must first multiply rewith a weight function w, then integrate the result over the area of the element, and finally, set the integral to zero. Z ZΩe w∂ ∂x αx ∂u ∂x+∂ ∂y αy ∂u ∂y +βu −gdxdy = 0 (4.20) A new formulation called as The Variational Formulation is introduced to develop the mathematical expression (4.20), this new formulations is given by the next identities ∂ ∂x wαx ∂u ∂x=∂w ∂x αx ∂u ∂x+w∂ ∂x αx ∂u ∂x(4.21) rearranging the identity as w∂ ∂x αx ∂u ∂x=∂ ∂x wαx ∂u ∂x−∂w ∂x αx ∂u ∂x=∂ ∂x wαx ∂u ∂x−αx ∂w ∂x ∂u ∂x (4.22) 46 CHAPTER 4. THE FINITE ELEMENT METHOD obviously with the term relationed to αyin integral (4.20) occurs the same that with αx, substituting expression (4.22) into (4.20) results Z ZΩe∂ ∂x wαx ∂u ∂x+∂ ∂y wαy ∂u ∂y dxdy −Z ZΩeαx ∂w ∂x ∂u ∂x +αy ∂w ∂y ∂u ∂y dxdy +Z ZΩe β w u dxdy =Z ZΩe w g dxdy (4.23) using the divergence theorem Z ZΩe∇· ~ AdA =IΓe ~ A·ˆandl (4.24) that means Z ZΩe∂Ax ∂x +∂Ay ∂y dxdy =IΓe ( ˆaxAx+ ˆayAy)·ˆandl (4.25) and applying the divergence theorem to the first integral of expression (4.23), results to Z ZΩe∂ ∂x wαx ∂u ∂x+∂ ∂y wαy ∂u ∂y dxdy =IΓe wαx ∂u ∂x ˆx+αy ∂u ∂y ˆydl (4.26) Substituting this expression into (4.23), the weak form of the differential equation reduces to −Z ZΩeαx ∂w ∂x ∂u ∂x +αy ∂w ∂y ∂u ∂y dxdy +Z ZΩe β w u dxdy =Z ZΩe w g dxdy −IΓe wαx ∂u ∂x ˆx+αy ∂u ∂y ˆydl (4.27) According to the Garlerkin approach, the weight function wmust be relationed to the same set of shape functions that are used to interpolate the main unknown, u. This unknown is interpolated using a set of Lagrange polynomials as in (4.3). Thus, u= n X j=1 ue jNj(4.28) 4.4 Scalar Problems 47 where Njis the corresponding shape functions based on nodal elements (in expression (4.3) is denoted as ϕe i(x, y)) and nis the total number of nodes that form the element domain Ωe. Imposing w=Niwith i= 1,2,3...Neand substituting into (4.27) −Z ZΩe αx∂Ni ∂x   Ne X j=1 ue j ∂Nj ∂x  +αy∂Ni ∂y   Ne X j=1 ue j ∂Nj ∂y   dxdy +Z ZΩe β Ni  Ne X j=1 ue jNj dxdy =Z ZΩe Nig dxdy −IΓe Niαx ∂u ∂x ˆx+αy ∂u ∂y ˆydl, for i = 1,2,3...Ne (4.29) Notice that the primary unknown uin the contour integral of (4.29) has not been replaced by interpolation functions given by (4.28), this integral will be treated separately. Equation (4.29) can be expressed in a matrix form given by         Me 11 Me 12 ··· Me 1n Me 21 Me 22 ··· Me 2n . . .. . ..... . . Me n1Me n2··· Me nn                        ue 1 ue 2 . . . ue n                +         Te 11 Te 12 ··· Te 1n Te 21 Te 22 ··· Te 2n . . .. . ..... . . Te n1Te n2··· Te nn                        ue 1 ue 2 . . . ue n                =                fe 1 fe 2 . . . fe n                +                pe 1 pe 2 . . . pe n                (4.30) where Me ij =−Z ZΩeαx∂Ni ∂x ∂Nj ∂x +αy∂Ni ∂y ∂Nj ∂y dxdy (4.31) Te ij =Z ZΩe β NiNjdxdy (4.32) fe i=Z ZΩe Nig dxdy (4.33) pe i=−IΓe Niαx ∂u ∂x ˆx+αy ∂u ∂y ˆydl (4.34) 48 CHAPTER 4. THE FINITE ELEMENT METHOD Is possible to implement a more compact matrix system,         Ke 11 Ke 12 ··· Ke 1n Ke 21 Ke 22 ··· Ke 2n . . .. . ..... . . Ke n1Ke n2··· Ke nn                        ue 1 ue 2 . . . ue n                =                be 1 be 2 . . . be n                (4.35) where Ke ij =Me ij +Te ij be i=fe i+pe i (4.36) Matrix Keis the FEM Matrix in each element. FEM Matrix is based on the basis functions that our method is using to represent the corresponding contribution and its relation respect the electromagnetic field. To assemble the FEM Matrix is neccesary to sum Stiffness Matrix (4.31) and Mass Matrix Matrix (4.32), on one hand Medepends on the boundary conditions of the domain, and on the other hand Tedepends on the dielectric properties of the element domain Ωe. The contour integral in (4.34) must be evaluated along the closed boundary of each element in the domain, imagine that the finite element mesh is based on triangular elements, in this case the contour integral should be evaluated along the three edges of each triangle in a counterclockwise direction. However, it is important to realize that a nonboundary edge belongs to two neighboring triangles, as shown in Fig.4.7(a). Evaluating the line integral in (4.34) for element e1in Fig.4.7(b) along the edge from node 1 to node 2 is exactly the same result but opposite sign than evaluating the same line integral for element e2along the edge from node 3 to node 1. 4.4 Scalar Problems 49 The opposite sign stems from the fact that the unit vector normal of the common edge point in opposite directions. ˆan1=−ˆan2(4.37) Figure 4.7: (a) Interior edge shared by two neighboring triangles. (b) The outward unit vectors normal to the common edge point in opposite directions To give an example, we compute the contribution to the global entry of the vector ~p by local node 1 of element e1and local node 1 of element e2, as the integration in (4.34) is evaluated along the common edge to the two triangles: the contribution to the global entry by local node 2 of element e1and local node 3 of element e2is also evaluated. p1−1=−Z1→2 N(e1) 1αx ∂u ∂x ˆx(e1)+αy ∂u ∂y ˆy(e1)dl −Z3→1 N(e2) 1αx ∂u ∂x ˆx(e2)+αy ∂u ∂y ˆy(e2)dl (4.38) p2−3=−Z1→2 N(e1) 2αx ∂u ∂x ˆx(e1)+αy ∂u ∂y ˆy(e1)dl −Z3→1 N(e2) 3αx ∂u ∂x ˆx(e2)+αy ∂u ∂y ˆy(e2)dl (4.39) Substituting (4.37) in the expressions above, and using the fact that N(e1) 1=N(e1) 1∧N(e1) 2=N(e2) 3(4.40) along the path of integration, it is evident that the two integrals cancel each other out. Consequently, the contribution of the line integral in (4.34) to the global right-hand-sive vector 50 CHAPTER 4. THE FINITE ELEMENT METHOD ~p is zero for all interior edges. It is nonzero for edges that belong to the domain boundary Γ. For boundary edges that belong to Γ1, where a Dirichlet boundary condition is to be imposed, the contribution of the line integral in (4.34) will be discarded. Thus, the only contribution of the line integral (4.34) is attributed only to boundary edges that reside on Γ2. Substituting the boundary condition in expression (4.15) into (4.34), the line integral becomes pe i=−ZΓ2 Ni(q−γ u)dl (4.41) The line integral in (4.41) exists only for boundary elements, it must be evaluated only along boundary edges that reside on Γ2. For interior edges, as commented before, the contribution is zero. 4.5 Vector Problems In this section vector problems are briefly described, due to scalar problems define our application they represent most of this research work, nevertheless we consider that it is important to develop vector theory in electromagnetics. While scalar problems are defined by a single magnitude, vector equations describe the behaviour of two or more variables. In Section 2.3 the Helmholtz Equations was introduced for electric fields, now we use that formulation as the curl-curl equation of electromagnetics. 4.5.1 The Curl-Curl Equation And Edge Elements To deal with vector quantities, such as the electric field, a first attempt might be to expand each vector component separately in nodal basis functions. It turns out that such an approach leads to nonphysical solutions, referred to as spurious modes. This can be avoided by using edge elements (Section 4.3.2), which are very well suited for approximating electromagnetic fields. The electromagnetic field is introduced as Helmholtz equation, performing a particular partial equation that is known as the curl-curl equation for ~ Efield. ∇×µ−1∇× ~ E(~r)−w20−jwσ~ E(~r) = −jw ~ Js(~r)on S (4.42) 4.5 Vector Problems 51 and the boundary conditions ˆn×~ E(~r) = ~p(~r)on Γ1(4.43) ˆn×(µ−1∇× ~ E(~r)) + γ(~r)ˆn×(ˆn×~ E(~r)) = ~q(~r)on Γ2(4.44) The next step is to follow the Garlerkin’s method developed during the last section, introducing the test or weigthed function into the residual expression and using the edge basis functions, the integral weak form is given by ZShµ−1∇× ~ Wi·∇× ~ E−w20−jwσ~ Wi·~ EidS +ZΓ2 ~ Wi·~q −γˆn׈n×~ Edl =−jw ZS ~ Wi·~ JsdS (4.45) We expand the solution ~ E(~r)in terms of the basis functions, the edges are labeled by integers 1, 2,..., Ne. ~ E(~r) = Ne X j=1 EjNj(4.46) where Ejis the tangential electric field along the jth edge, choosing ~ Wi=Niand substituting in (4.45) we obtain a linear system of equations Az =bwith Aij =ZSµ−1(∇×Ni)·(∇×Nj)−w20−jwσNi·NjdS +ZΓ2 γ(ˆn×Ni)·(ˆn×Nj)dl (4.47) zj=Ej(4.48) bi=−jw ZS Ni·JsdS −ZΓ2 Ni·qdl (4.49) Notice that during this section we refer to 0as the real part of the permittivity, assuming that =0+j00 =or=o(0 r+j00 r). 58 CHAPTER 4. THE FINITE ELEMENT METHOD various PML media in such a way that the reflection is zero from all the inner interfaces in the domain. This is the case at the vacuum and PML interfaces. Figure 4.11: Concave surface enclosing the domain There is a second way to describe the PML problem, it is known as Anisotropic Formulation which presents (4.58) and (4.59) as ∇× ~ E(~r) = −jwµ ~ H(~r)(4.71) ∇× ~ H(~r) = jw ~ E(~r)(4.72) where =     xx 0 0 0yy 0 0 0 zz      µ=     µxx 0 0 0µyy 0 0 0 µzz      (4.73) and these tensors can be expressed as a function of stretching coordinates Sx,Syand Szlike =orΛµ=µoµrΛ(4.74) 4.6 Absorbing Boundary Conditions 59 where Λ = ˆxˆxsysz sx + ˆyˆysxsz sy + ˆzˆzsxsy sz (4.75) Now we consider a TM and z polarization, so the Helmholz equations is given by ∂ ∂x 1 µxx ∂Ez ∂x +∂ ∂y 1 µyy ∂Ez ∂y −k2 ozzEz=−jwµoJz(4.76) Then, the tensors are defined by =orsxsyµ=µoµr  sy/sx0 0sx/sy (4.77) The next step is how to implement the longitudinal stretching coordinates in each of the possible vacuum and PML interfaces, there are several theories about this object of study. A possibility is given in [7], as the next form, sx(x) = so(x)1−jσx(x) w0with 0=0 ro(4.78) so(x) = 1 + smπ δ2(4.79) σx(x) = sin2πx 2δ(4.80) where δis the thickness of the absorber and smis a coefficient that depends on the wavelength. Other possibily is given by [1] as sx(x)=1−jx−δ δ2 δmax (4.81) δmax =σmax wo or δmax ≈w−1(4.82) a good choice for the absorber thickness could be δ=λ/4. 60 CHAPTER 4. THE FINITE ELEMENT METHOD Finally, we analyze the practical case described in Fig.4.11 which has been implemented during this research work: 1. Region Ω.sx=s∗ x= 1 ∧sy=s∗ y= 1 2. Region I. sx= (4.81) ∧sy= 1 3. Region II. sx= (4.81) ∧sy= 1 4. Region III. sx= 1 ∧sy= (4.81) 5. Region IV. sx= 1 ∧sy= (4.81) 6. Region V. sx=sy= (4.81) 7. Region VI. sx=sy= (4.81) 8. Region VII. sx=sy= (4.81) 9. Region VIII. sx=sy= (4.81) 4.7 Assembling FEM Matrix In Chapter 3 we described how the Forward Problem can be solved measuring the electromagnetic field generated by a set of antennas using the finite element method, to obtain the coefficients that describe the electric field is necessary to solve the matrix equation presented in (3.2), it can be developed according to the formulation described in Section 4.4.2 for scalar problems using nodal basis functions or the formulation described in Section 4.5.1 for vector problems using edge elements. Both of them can be expressed as [1] explain, Ku =b(4.83) where K represents the Global FEM Matrix, u is the vector of electric coefficients and b is a vector called as right-hand side formed by the sum of vector f and p. First we analyse the assemble procedure related to nodal basis functions that represent the implementation of our application. In practice, the FEM matrix is computed by assembling contributions from all 4.7 Assembling FEM Matrix 61 elements, this means that compute the local FEM matrix for each element of the total domain Ωis necessary to implement the global matrix. Matrix expression (4.35) is the general matrix equation for any type of elements used during the domain discretization, in our application the discretization is based on a triangular meshing, so the FEM matrix for a certain element will present the next form Ke=     Ke 11 Ke 12 Ke 13 Ke 21 Ke 22 Ke 23 Ke 31 Ke 32 Ke 33      (4.84) After the assembly of the local matrices, the global numbering scheme is used to build the global FEM matrices. The global matrix K is filled using the following scheme: the first local element Ke i,i in Keis added to the ith row and the ith column of global matrix K, the second local element Ke i,j is added to the ith row and the jth column of global matrix K, and so on. This is repeated for each triangular element. Assembling contributions from all elements, we obtain Kij = Ne X e=1 ZSeα∇ϕe i·∇ϕe j+βϕe iϕe jdS (4.85) In Appendix A.1 a detailed analytical evaluation of FEM matrix elements is given for triangular meshing based on nodal basis functions. A very important aspect about the implementation of FEM matrix is the concept of sparse matrix, as mentioned in Section 4.1 the main disadvantage of the finite element method against other methods is the cost of computation, sparse matrix is based on processing matrices on the way that only are stored the non-zero values of FEM matrix, this method let us to saves us a considerable cost of time and computation. For vector problems where edge elements are involved to improve results the FEM matrix is defined by expression (4.47) with γ= 0. The assembling of Stiffness and Mass matrix is developed in Appendix A.2, formulation obtained in [1]. 62 CHAPTER 4. THE FINITE ELEMENT METHOD 4.8 Simulations Finally, this section presents several simulations based on 2D geometries, the purpose is to visualize the results of using FEM to compute the electric field. Due to scalar problems, the basis functions that define our method are the Nodal Elements. According to several combinations of current sources, a electric distribution is calculated. In addition, to understand the behaviour of the EM field, different material properties defined by the dielectric parameters will be considered. 4.8.1 Circular Enclosure In this section we analyse simulations where the PEC enclosure is a circle, that is an essential factor due to wave reflections in PEC interface. The radius of the circular enclosure is equal to 0.1m. It is necessary to comprise that in real biomedical MWI applications, normally the target to analyse presents a ”small” size assuming a common tumor. In addition, to minimize the computation cost of simulations is better to implement ”small” scenarios. The purpose is to observe different results for different dielectric properties, frecuency and placement for current sources. Firstly, several dielectric cases are computed with fixed frequency at 500MHz, an essential point, commented during Section 4.2, is the edge length that is implemented as a function of the wavelength, which means that for each permittivity case the number of elements in the corresponding mesh changes. For the next simulations 0 r= 20, so the designed meshing is given by Figure 4.12: 2D meshing for 0 r= 20. 4.8 Simulations 63 We consider a single current source allocated in cartesian coordinates (Xs,Ys)=(0.5,0) for different values of conductivity. (a) (b) (c) (d) Figure 4.13: Forward Problem solved by FEM using nodal functions: 0 r= 20, (Xs,Ys)=(0.5,0), (a)(b) |E|and ∠Efield for σ= 0.0005 S/m (c)(d) |E|and ∠Efield for σ= 0.5S/m. We can extract important information from results above, according to the chosen dielectric properties the number of wave reflections that appear in the computed electric field can be higher or lower. On one hand, when the permittivity and the conductivity are more similar, the peaks of electric intensity are focused at the placement of the given source. On the other hand, if the conductivity is enough small, there will be electric duplicities due to the low attenuation that performs the medium which involves reflections in the PEC. 64 CHAPTER 4. THE FINITE ELEMENT METHOD We talk about lossy mediums (attenuate travel waves) when the conductivity can be considered high. A given medium is a well conductor when the relation σ/(wo0)>> 1is satisfied. Now, we present some simulations based on Fig.4.13.(a) with different current source position and different number of sources to understand better the behaviour of the given 2D scenario. (a) (b) (c) (d) Figure 4.14: Forward Problem solved by FEM using nodal functions: 0 r= 20 (a)(b) |E|and ∠E field for σ= 0.0005 [S/m] and for current sources (Xs,Ys)=([0.05 0],[0 -0.05]), (c)(d) |E|and ∠E field for σ= 0.0005 [S/m] and for current sources (Xs,Ys)=([-0.05 0.05],[-0.05 0.05]). In Fig.4.14.(a) is possible to check that appears an overlap between the two charges due to the ¨small¨ distance among both, this problem has a reason that will be evaluated during the next Chapter 5. 4.8 Simulations 65 Let’s consider a new scenario characterized by a higher permittivity, 0 r= 80, obviously the behaviour of the electric field in the medium will be different. It’s not necessary to show the new surface because presents the same structure but with only more number of elements. In Fig.4.15.(a-b), it is possible to observe how the number of performed reflections in PEC has been more numerous. In addition, in Fig.4.15.(c-d) the current source has been placed respectively, farther and nearer to the PEC interface to visualize how the wave reflections decrease or increase because there is more or less space to attenuate the wave propagations. (a) (b) (c) (d) Figure 4.15: Forward Problem solved by FEM using nodal functions: 0 r= 80 (a)(b) |E|and ∠Efield for σ= 0.0005 [S/m] and for current sources (Xs,Ys)=(0.05,0), (c)(d) |E|field for σ= 0.0005 [S/m], current sources (Xs,Ys)=(0.02,0) and (Xs,Ys)=(0.07,0). 66 CHAPTER 4. THE FINITE ELEMENT METHOD We have analysed circular surfaces where the dielectric properties are homogeneous, but in real applications as thomography where the electric field should be computed, dielectric properties variate along the surface, i.e., they are inhomogeneous. In Fig.4.16 inhomogeneous 2D scenarios are presented, they are composed by three different parts, whose permittivities and conductivities are given below the figures. (a) (b) (c) (d) Figure 4.16: (a)(b) Inhomogeneous 2D Scenarios (c)(d) |E|fields for current source at (Xs,Ys)=(0,-0.05). •Green Surface ⇒0 r= 20;σ= 0.0005; •Yellow Surface ⇒0 r= 80;σ= 0.02; •Blue Surface ⇒0 r= 6;σ= 0.1; 4.8 Simulations 67 We can observe how in Fig.4.16.(a)(c) the reflections appear vividly around the first region where the difference between permittivity and conductivity is greater (low losses), however in blue region no electric intensity is involved due to the high level of losses. The analysis in Fig.4.16.(b)(d) is different because the electric charge is placed in the blue region where the behaviour is similar to Fig.4.13.(c), and the reflections waves are more attenuated. Now let’s fix the dielectric properties to compute some EM results for circular enclosure as a function of the working frequency, thereby we may understand the implications of variations in frequency. Consider a circular enclosure with radius equals to 0.3m, where the surface is homogeneous in terms of dielectric parameters (0 r= 10,σ= 0.03), the electric charge is based on a single current source placed in (Xs,Ys)=(0,0.1m). In Fig.4.17 it is possible to visualize how the increase of the working frequency generates a greater number of field lines, the question that we should arise is why reflections due to PEC interface not appear in the shown results. The reason is double, on one hand, the tangent of losses is lower due to high frequencies. On the other hand, the circular enclosures is three times bigger, so the space to attenuate the wave propagations due to the electric charge is higher. Then, if the second hypothesis was correct we would obtain more wave reflections in results using a smaller PEC enclosure. In order to perform this experiment, we use the PEC enclosure defined in the first circular scenario (radius=0.1m). In Fig.4.18 two different cases are depycted. The dielectric properties used in Fig.4.17 (0 r= 10,σ= 0.03) are considered. Then, we can compare the results against Fig.4.17(a)(d). For a working frequency equals to 700MHz, we obtain a higher level of electric field in the peaks where reflections are focused. With the frequency equals to 2GHz happens the same, so we can conclude that the hypotheis approached above is certain. Finally, in the next sections other type of enclosures based on triangles or squares are developed briefly, the results will only change due to the new geometry and the corresponding new behaviour of the wave reflections between the medium and the PEC interface. The aim of these simulations is to show some results with different enclosures and compare some computational cost as a function of material and frequency choice.