Propagación de ondas trasversales generadas con fuerza de radiación

UNIVERSIDAD POLITECNICA DE VALENCIA ESCUELA POLITECNICA SUPERIOR DE GANDIA I. T. Telecomunicación (Sonido e Imagen) “PROPAGACIÓN DE ONDAS TRANSVERSALES GENERADAS CON FUERZA DE RADIACIÓN” POYECTO FIN DE CARRERA Autor: Luis J. Salmerón Contreras Director/es: D. Rubén Picó Vila D. Noé Jiménez González GANDIA, 2011 2 Propagación de ondas transversales mediante fuerza de radiación Resumen El estudio de la propagación de ondas transversales en el interior de un tejido aporta información acerca de la elasticidad de éste que puede emplearse para diseñar nuevas técnicas de diagnóstico médico. Sin embargo, la alta atenuación de estas ondas dificulta extraer información de capas profundas. En este trabajo se emplea un haz focalizado de ultrasonidos para inducir Fuerza de Radiación en el interior del tejido, generando así un foco de ondas transversales secundarias. El objetivo del presente trabajo es estudiar la propagación de dichas ondas transversales mediante predicción numérica basada en diferencias finitas en dominio temporal (FDTD), introduciendo excitaciones acústicas móviles bajo un sistema de referencia euleriano. Los resultados muestran cómo mediante la adecuada excitación en profundidad del tejido es posible generar ondas transversales supersónicas que aportan nuevas posibilidades a la hora de caracterizar el tejido mejorando así las capacidades de diagnóstico. Autor: Luis J. Salmerón Contreras, email: [email protected] Fecha de entrega: 23-12-2011 PROPAGACIÓN DE ONDAS TRANSVERSALES GENERADAS CON FUERZA DE RADIACIÓN Autor: Luis J. Salmerón Contreras Tutor: Rubén Picó Vila Cotutor: Noé Jiménez González Propagación de ondas transversales generadas con fuerza de radiación 3 Índice I. OBJETIVOS ............................................................................................................................................ 4 II. INTRODUCCIÓN TEÓRICA ............................................................................................................... 5 II. 1. Elastografía ....................................................................................................................................... 5 II. 2. Modelo físico .................................................................................................................................... 7 II. 2. a. Mecánica de Sólidos. Fórmula de Cauchy .............................................................................. 8 II. 2. b. Elasticidad ............................................................................................................................ 12 a) Deformación Axial (Ley de Hooke): ......................................................................................... 13 b) Deformación lateral: .................................................................................................................. 15 c) Deformación angular. ................................................................................................................ 15 II. 2. c. Leyes de Hooke Generalizadas ............................................................................................. 16 II. 2. d. IIª Ley de Newton, ecuación de la Dinámica ........................................................................ 18 II. 2. e. Ecuaciones constitutivas del modelo ..................................................................................... 19 II. 3. Modelo Numérico ........................................................................................................................... 20 II. 3. a. FDTD (Finite Difference Time Domain) ............................................................................... 20 II. 3. b. Series de Taylor ..................................................................................................................... 20 II. 3. c. Simplificación del modelo, simetría axial ............................................................................. 21 II. 3. d. Absorción y PML ................................................................................................................... 22 II. 3. e. Discretización y tresbolillado del modelo ............................................................................. 24 II. 3. f. Consideraciones del modelo ................................................................................................. 25 III. CONDICIONES DE SIMULACIÓN ................................................................................................... 26 IV. RESULTADOS ...................................................................................................................................... 28 IV. 1. Foco estático, transmisión OL ........................................................................................................ 28 IV. 2. Foco estático, transmisión OS ........................................................................................................ 29 IV. 3. Foco estático, transmisión OL+OS ................................................................................................. 31 IV. 4. Foco dinámico, transmisión OL+OS .............................................................................................. 32 V. REFERENCIAS .................................................................................................................................... 35 Anexo I. Ecuaciones del modelo numérico .................................................................................................. 36 4 Propagación de ondas transversales mediante fuerza de radiación I. OBJETIVOS Este proyecto tiene como objetivo principal el estudio de la propagación de éstas ondas transversales generadas mediante fuerza de radiación. La metodología empleada principalmente es la predicción numérica basada en diferencias finitas en dominio temporal. A partir de las deformaciones producidas mediante fuerza de radiación generada por un haz de ultrasonidos focalizados en el interior de un tejido biológico blando se induce en dicho tejido un foco secundario de ondas transversales. El presente trabajo aporta un método de análisis para la propagación de las ondas transversales mediante la introducción de una fuente acústica acorde con las deformaciones producidas por el haz principal de ultrasonidos focalizados. En un segundo paso, variando la posición del foco del haz de ultrasonidos en función del tiempo a una velocidad superior a la velocidad de propagación de las ondas transversales se logra una fuente acústica móvil capaz de producir un frente de ondas transversales supersónicas. A través de las características de propagación de este frente de ondas se extrae información muy precisa acerca de la elasticidad del tejido por el que se propagan las ondas transversales, permitiendo obtener una imagen de diagnóstico médico de la elasticidad del tejido biológico. Propagación de ondas transversales generadas con fuerza de radiación 5 II. INTRODUCCIÓN TEÓRICA II. 1. Elastografía La elastografía es un método de imagen médica no invasivo en el cual se caracteriza la rigidez o la deformación de los tejidos blandos utilizando un sistema basado en ultrasonidos. En los últimos años esta técnica de imagen se está investigando y utilizando para la detección y clasificación de tumores. Un tumor o un cáncer sospechoso normalmente es de 5 a 28 veces más rígido que el fondo de tejido blando normal. Cuando se aplica una compresión mecánica o vibración, el tumor se deforma menos que el tejido circundante. Por lo tanto, una imagen de la deformación puede, bajo determinadas hipótesis simplificadoras, interpretarse como representación de la distribución subyacente del módulo de Young (o módulo de elasticidad). Se ha demostrado que los elastogramas (imágenes de la deformación del tejido) se ven afectados por el grado de adherencia del tumor a su entorno, lo que indica una gran oportunidad para extender la elastografía a la caracterización de la movilidad del tumor para mejorar la precisión diagnóstica y orientación quirúrgica. TECNICAS DE ELASTOGRAFIA ARFI (Acoustic Radiation force impulse) En este método un haz de ultrasonidos focalizados (FUS, focused ultrasound) se utiliza para aplicar una fuerza de radiación localizada en pequeños volúmenes durante poco tiempo y se miden los desplazamientos de tejido resultantes. Los desplazamientos de los tejidos son inversamente proporcionales a la rigidez de los tejidos y, por tanto, una zona de tejido más rígido presenta desplazamientos más pequeños que una zona menos rígida. Debido a la corta duración de la aplicación de la fuerza, este método proporciona información sobre la respuesta mecánica al impulso del tejido, lo que refleja las variaciones en las características viscoelásticas del tejido. HMI (Harmonic Motion Imaging) Esta técnica basada también en la fuerza de radiación, induce desplazamientos oscilatorios (harmónicos) en la zona focal del transductor FUS (focalizado) para la detección de cambios en la rigidez de zonas localizadas. SWEI (Shear wave elasticity imaging) Esta técnica es un nuevo enfoque de la imagen y la caracterización de las estructuras de tejido basado en la utilización de ondas acústicas transversales inducidas a distancia por la fuerza de radiación de un haz de ultrasonido enfocado. SWEI proporciona al médico un "dedo" 6 Propagación de ondas transversales mediante fuerza de radiación virtual para probar la elasticidad de las regiones internas del cuerpo. En SWEI, en comparación con otros métodos de elastografía, la tensión inducida en el tejido puede ser muy localizada, ya que las ondas transversales inducidas de forma remota se atenúan por completo dentro de un área muy limitada de tejido en las proximidades del punto focal del haz focalizado. SWEI puede añadir una nueva calidad a los métodos convencionales de ultrasonido o resonancia magnética. La adición de datos de corte elasticidad ("información palpación") mediante la superposición de código de color de datos sobre la elasticidad de las imágenes de resonancia magnética o ultrasonido puede permitir una mejor diferenciación de los tejidos y aumentar aún más el diagnóstico. Se demuestra que la exposición típica de ultrasonidos de SWEI está significativamente por debajo del umbral de los efectos dañinos de ultrasonidos. SSI (Supersonic Shear Imaging). Es una técnica de imagen de elasticidad de tejidos basada en la estimación de la velocidad de una onda transversal generada con fuerza de radiación. Se genera una fuerza de radiación a distancia (o “pushing beam”, haz de empuje) mediante un transductor focalizado utilizando una sonda ultrasónica convencional. Esta fuerza de radiación da lugar a en un desplazamiento de unos pocos micrómetros de tejido que hace que se propague una onda transversal en el medio. Si sucesivamente focalizamos varios haces en profundidades crecientes obtenemos un frente de ondas transversales cuasi plano. Después de la generación de la onda transversal, el dispositivo de medida cambia al modo de adquisición de datos a una frecuencia de imagen muy elevada (más de 5000 imágenes/seg.) utilizando la misma sonda ultrasónica. Cada ecografía ultrarrápida se consigue transmitiendo una única onda plana ultrasónica formada durante el proceso de recepción en el postprocesado. Luego utilizando un algoritmo de correlación cruzada en una dimensión, las consecutivas imágenes ultrasónicas se comparan para determinar el desplazamiento axial inducido por la propagación la onda transversal dentro del tejido. Propagación de ondas transversales generadas con fuerza de radiación 7 II. 2. Modelo físico Con el objetivo de caracterizar como interacciona una onda a lo largo de un medio, en este apartado se estudiará las propiedades elásticas del medio (Ley de Hooke) y como éste interactúa cuando está sometido a la acción del movimiento (Ley de Newton), ésta parte de la Física se conoce como Mecánica Racional. Éste estudio de los estados de equilibrio y movimiento de los cuerpos, se subdivide a su vez en las siguientes ramas: Mecánica del Punto Material: en esta rama se asume que toda la masa de un cuerpo está concentrada en un punto llamado centro de masas o de gravedad. Con esta teoría se puede explicar, entre otros, el movimiento de los planetas Mecánica del Sólido Rígido: ésta teoría se basa en qué las dimensiones del material del cuerpo no se ven alteradas ante la acción de una fuerza, es decir, su forma no varía y el cuerpo es indeformable. Con ella se explica satisfactoriamente los movimientos que puede sufrir un móvil y caracterizarlos matemáticamente. Mecánica de los Medios Continuos (MMC): en este caso el cuerpo sí puede deformarse ante la acción de una fuerza (incluso llegando al límite de romperse). Dependiendo del estado de la materia que lo compone, esta teoría se subdivide en Mecánica de Fluidos y Mecánica de Sólidos. Esta última, a su vez, se subdivide dependiendo del comportamiento que se asuma para el cuerpo (Elasticidad, Plasticidad, Viscoelasticidad y Viscoplasticidad). Para este estudio es interesante por tanto, definir que propiedades presentará el sólido elástico: − Deformable elástico, puede variar su forma ante la acción de una fuerza y volver a su estado inicial cuando ésta deja de aplicarse. − Continuo, no existen saltos (vacío) entre una porción de material y el inmediatamente siguiente. Esta implicación nos permite tomar elementos infinitesimales del material y que estos están completamente unidos a los circundantes. − Homogéneo, el medio presenta la misma composición y propiedades en cualquier parte del mismo. − Isótropo, las propiedades del material son las mismas en cualquier dirección del espacio, es decir una rotación del cuerpo no implica que los fenómenos que vallan a ocurrir sean diferentes. 8 Propagación de ondas transversales mediante fuerza de radiación II. 2. a. Mecánica de Sólidos. Fórmula de Cauchy Supongamos un sólido elástico que se encuentre en equilibrio estático (o en movimiento de traslación y rotación uniformes, es decir la fuerza y momento torsor resultante es nulo: ∑𝐹𝚤 � �  𝑖= 0 �  y ∑𝑀𝑝 � � � � �  𝑖= 0 �  ), pero que sí esté sujeto a fuerzas (que opuestas vectorialmente se anulan, pero no en módulo, ∑�𝐹𝚤 � �  � 𝑖> 0), éstas tienden a deformar el objeto. La deformación que presenta el cuerpo no sólo depende de la fuerza, también depende inversamente de la superficie donde se aplique la fuerza, es decir la misma fuerza actuando sobre una superficie más grande obtendrá una menor deformación. Por lo tanto, si tomamos un área muy pequeña (pero finita) ∆𝑆 se define vector tensión mecánica o esfuerzo (𝝉), como el vector resultante de calcular la relación entre la fuerza aplicada y la superficie dónde se aplica cuando ésta tienda a cero: 𝜎= lim ∆𝑆→0∆𝐹 ∆𝑆 �𝑃𝑎=𝑁 𝑚2� A efectos prácticos, suponiendo que la distribución de fuerzas es constante para toda la superficie 𝑆𝑛: 𝜎=�𝜎𝑥 𝜎𝑦 𝜎𝑧�=𝐹 𝑆𝑛 El vector esfuerzo dependerá de la superficie 𝑆𝑛 caracterizada por el vector normal 𝑛�  = (𝑛𝑥,𝑛𝑦,𝑛𝑧), perpendicular a ésta y de módulo la unidad. Por lo tanto podemos definir las componentes normal y tangencial del vector esfuerzo: − Esfuerzo normal: 𝜎=𝜏·𝑛�  =𝜏𝑥𝑛𝑥+𝜏𝑦𝑛𝑦+𝜏𝑧𝑛𝑧 − Esfuerzo tangencial: 𝜏=√𝜏·𝜏−𝜎2 A pesar de la sencillez de trabajar con un vector, el cálculo se vuelve complejo si suponemos que por un punto pueden pasar infinitos planos y por consecuente infinitos vectores de tensión referidos a cada superficie. Por tanto otra forma conveniente de trabajar el estado tensional de un cuerpo es tomar un paralepípedo, de aristas diferenciales y caras paralelas a los planos coordenados, correspondiente a la zona infinitesimal en torno al punto de aplicación. Al ser tan pequeño el cubo, se puede suponer que no existen variaciones de tensión entre caras paralelas, por tanto con 3 vectores tensión aplicados en cada una de las 3 caras independientes, podemos caracterizar su estado tensional. Propagación de ondas transversales generadas con fuerza de radiación 9 Reescribiendo cada vector por sus componentes y agrupándolos es una matriz obtenemos lo que se conoce como Matriz ó Tensor de Tensiones ó de Esfuerzos: [𝜏]=�𝜎𝑥𝑥 𝜏𝑥𝑦 𝜏𝑥𝑧 𝜏𝑦𝑥 𝜎𝑦𝑦 𝜏𝑦𝑧 𝜏𝑧𝑥 𝜏𝑧𝑦 𝜎𝑧𝑧� Llegados a este punto es conveniente remarcar que esta matriz posee simetría. Para demostrarlo se aplicará la condición de equilibrio del momento torsor ∑𝑀𝑝 � � � � �  𝑖=𝑟×𝐹= 0 �  y se obtendrá que las componentes tangenciales no son independientes si no que son iguales en parejas 𝜏𝑖𝑗=𝜏𝑗𝑖. Fijándonos en una cara del paralepípedo, se pueden expresar las componentes de esfuerzo en torno al punto de aplicación como, la magnitud de dicho esfuerzo más un pequeño incremento producido por no estar aplicado directamente sobre el punto (serie de Taylor). Éste incremento dependerá de cómo varía esa función con respecto a la dimensión escogida por la distancia recorrida: 𝜎�𝑥+𝑑𝑥 2�≃𝜎(𝑥)+𝜕𝜎 𝜕𝑥𝑑𝑥 2 𝜏𝑧𝑦 𝜎𝑧𝑧 𝜎𝑥𝑥 𝜎𝑦𝑦 𝜏𝑥𝑧 𝜏𝑥𝑦 𝜏𝑦𝑥 𝜏𝑦𝑧 𝜏𝑧𝑥 𝑥 𝑧 𝑦 𝑇(𝑧) 𝑇(𝑦) 𝑇(𝑥) 16 Propagación de ondas transversales mediante fuerza de radiación El módulo de rigidez o segundo coeficiente de Lamé se relaciona con los módulos de Young y de Poisson mediante la expresión: 𝜇=𝐸 2(1 + 𝜈) Por último, para tener una relación con la deformación producida supongamos un sólido elástico al que se le aplica un esfuerzo tangencial en una de sus caras, y para que no obtenga momento torsor, otro de igual magnitud pero de signo contrario. Esto provocará una deformación angular en otra de sus caras que se relaciona con la deformación. 𝜀𝑗𝑖=𝛿ℎ=ℎ2tan 𝛼 ℎ≈𝛼2 𝜏𝑗𝑖=𝜇·𝛼=𝐸 (1 + 𝜈)𝜀𝑗𝑖 II. 2. c. Leyes de Hooke Generalizadas Según el principio de superposición se puede descomponer los esfuerzos que actúan sobre un sólido para analizar sus efectos y luego sumar la aportación de cada análisis al conjunto total. Para ello, se tiene un sólido sometido a esfuerzos en cada una de sus caras, siendo los esfuerzos de caras paralelas iguales en magnitud y opuestos en sentido: = + + 𝜎𝑥 𝜎𝑧 𝜎𝑦 𝜎𝑥 𝜎𝑦 𝜎𝑧 𝜎𝑧 𝜎𝑧 𝜎𝑥 𝜎𝑥 𝜎𝑦 𝜎𝑦 𝜀𝑥=𝜎𝑥 𝐸 𝜀𝑦=−𝜈𝜀𝑥=−𝜈𝜎𝑥 𝐸 𝜀𝑧=−𝜈𝜀𝑥=−𝜈𝜎 𝑥 𝐸 𝜀𝑥=−𝜈𝜎𝑦 𝐸 𝜀𝑦=𝜎𝑦 𝐸 𝜀𝑧=−𝜈𝜎 𝑦 𝐸 𝜀𝑥=−𝜈𝜎𝑧 𝐸 𝜀𝑦=−𝜈𝜎𝑧 𝐸 𝜀𝑧=𝜎 𝑧 𝐸 𝜏𝑗𝑖 𝑗 𝑖 𝛿 ℎ 𝛼 Propagación de ondas transversales generadas con fuerza de radiación 17 Sumando la aportación de cada contribución y añadiendo las posibles contribuciones de esfuerzos tangenciales se obtienen las leyes de Hooke generalizadas: 𝜀𝑥=𝜎𝑥 𝐸−𝜈𝐸�𝜎𝑦+𝜎𝑧� 𝜀𝑦=𝜎𝑦 𝐸−𝜈𝐸(𝜎𝑥+𝜎𝑧) 𝜀𝑧=𝜎𝑧 𝐸−𝜈𝐸�𝜎𝑥+𝜎𝑦� 𝜀𝑥𝑦=(1 + 𝜈) 𝐸𝜏𝑥𝑦 𝜀𝑥𝑧=(1 + 𝜈) 𝐸𝜏𝑥𝑧 𝜀𝑦𝑧=(1 + 𝜈) 𝐸𝜏𝑦𝑧 Estas ecuaciones se pueden resumir de forma matricial, siendo 𝑆[𝑚2/𝑁] el tensor compliancias del material, que en el caso de ser isótropo sólo tiene dos constantes independientes inherentes del material: 𝜀=𝑆:𝜏 (𝜀𝑖𝑗=𝑆𝑖𝑗𝑘𝑙·𝜏𝑘𝑙)                               + + + −− −− −− = E E E E EE E E E EE E S ν ν ν νν νν νν 1 00000 0 1 0000 00 1 000 000 000 000 En ocasiones es necesario conocer la relación inversa, es decir, a partir de la deformación generada conocer la tensión que la ha producido. Las ecuaciones resultantes se conocen como Ecuaciones de Lamé: 𝜏=𝐶:𝜀 (𝜏𝑖𝑗=𝐶𝑖𝑗𝑘𝑙·𝜀𝑘𝑙) Siendo 𝐶 [𝑁𝑚2 ⁄]=𝑆−1, tensor de elasticidad:                     + + + = µ µ µ µλλλ λµλλ λλµλ 200000 020000 002000 0002 0002 0002 C 18 Propagación de ondas transversales mediante fuerza de radiación Donde: Primer coeficiente de Lamé o constante de dilatación, 𝜆=𝜈𝐸 (1+𝜈)(1−2𝜈) Segundo coeficiente de Lamé o módulo de rigidez, 𝜇=𝐸 2(1+𝜈) Las ecuaciones desarrolladas serían: �𝜀𝑇=𝜀𝑥+𝜀𝑦+𝜀𝑧� 𝜎𝑥=𝜆𝜀𝑇+ 2𝜇𝜀𝑥 𝜎𝑦=𝜆𝜀𝑇+ 2𝜇𝜀𝑦 𝜎𝑧=𝜆𝜀𝑇+ 2𝜇𝜀𝑧 𝜏𝑥𝑦= 2𝜇𝜀𝑥𝑦 𝜏𝑥𝑧= 2𝜇𝜀𝑥𝑧 𝜏𝑦𝑧= 2𝜇𝜀𝑦𝑧 II. 2. d. IIª Ley de Newton, ecuación de la Dinámica La segunda ley de Newton explica la relación existente entre la fuerza que se le aplica a un cuerpo y el cambio en la cantidad de movimiento (𝑝) que éste sufre en el tiempo. 𝐹=𝑑𝑝 𝑑𝑡 En la mecánica clásica, la cantidad de movimiento (momento lineal o momentum) se define como el producto de la masa (𝑚) por su velocidad en un instante determinado (𝑣). Si queremos conocer el momento lineal para todo un cuerpo continuo, es suficiente integrar para todo el volumen (𝑉) el producto de la velocidad para cada porción infinitesimal de masa (𝑑𝑚). 𝑝=�𝑣·𝑑𝑚=�𝑣𝜌·𝑑𝑉 𝑉 Si suponemos que la masa es constante, el diferencial de masa es el producto entre la densidad (𝜌 [𝐾𝑔/𝑚3]) y un diferencial de volumen. Por lo tanto, podemos reescribir la ley de Newton como: 𝐹=�𝜌𝜕𝑣 𝜕𝑡𝑑𝑉 𝑉 En la otra parte de la igualdad podemos reescribir la fuerza como el conjunto de todos los esfuerzos que actúan por diferencial de superficie. 𝐹=�𝜎·𝑑𝑆 𝑆 Propagación de ondas transversales generadas con fuerza de radiación 19 Mediante la fórmula de Cauchy y a posterior aplicando el teorema de la divergencia obtenemos: �𝜎·𝑑𝑆 𝑆=�𝜏·𝑛�  𝑆𝑑𝑆=�∇·𝜏 𝑑𝑉 𝑉 Igualando obtenemos �𝜌𝜕𝑣 𝜕𝑡𝑑𝑉 𝑉=�𝛻·𝜏 𝑑𝑉 𝑉 𝜌𝜕𝑣 𝜕𝑡=𝛻·𝜏 II. 2. e. Ecuaciones constitutivas del modelo El modelo propuesto para la propagación de ondas elásticas en sólidos viene dado por las ecuaciones de Lamé extraídas de la Ley de Hooke en el apartado III. 3 y de la ecuación de la dinámica (segunda ley de Newton) en el apartado III. 4. Para poder formar el modelo es necesario relacionar las dos ecuaciones, la estática y la dinámica entre ellas. Para ello, se persigue que ambas expresiones estén en función de las mismas variables, en este caso serán la velocidad de partícula (𝑣) y la matriz de esfuerzos (𝜏). �𝜏=𝐶:𝜀 𝜌𝜕𝑣 𝜕𝑡=𝛻·𝜏 La relación entre la deformación unitaria y la velocidad es la derivada temporal de la primera variable: 𝜕𝜀 𝜕𝑡=1 2(∇𝑣+𝑣∇) Derivando con respecto del tiempo y sustituyendo en la primera ecuación constitutiva se obtiene: 𝜕𝜏 𝜕𝑡=𝐶2(∇𝑣+𝑣∇) Siendo 𝑣∇ la matriz traspuesta de ∇𝑣. Si ahora se reordenan los términos el modelo resultante para la propagación de ondas en un medio elástico resulta: �𝜕𝜏 𝜕𝑡=𝜆𝐼(𝛻𝑣)+𝜇(𝛻𝑣+𝑣𝛻) 𝜌𝜕𝑣 𝜕𝑡=𝛻·𝜏 𝐼, matriz identidad. 20 Propagación de ondas transversales mediante fuerza de radiación II. 3. Modelo Numérico II. 3. a. FDTD (Finite Difference Time Domain) El modelo que se ha obtenido resulta ser un sistema de ecuaciones diferenciales acopladas. Un estudio analítico de estas ecuaciones resultaría complejo y no en todos los casos sería posible una solución, sólo se podrían obtener una colección de soluciones analíticas para casos muy particulares. Para resolver el problema es recomendable, por lo tanto, realizar un estudio numérico del mismo. El método escogido es el de diferencias finitas en el dominio del tiempo (FDTD, Finite Difference Time Domain), este algoritmo se basa en discretizar el sistema de ecuaciones tanto en espacio como en tiempo e ir resolviéndolas explícitamente en dichos pasos temporales en un dominio finito y discreto. Las ventajas de éste método numérico frente a otros es que el hecho de utilizar una diferencia finita como aproximación a una derivada lo hace sencillo y el hecho de trabajar en el dominio del tiempo lo hace muy práctico frente a otros métodos que trabajan en el dominio de la frecuencia y por tanto sería imposible medir cambios en frecuencias que no sean la de estudio, como se desea realizar en este caso. II. 3. b. Series de Taylor Según el Teorema de Taylor una función en un punto se puede aproximar mediante la evaluación de las derivadas de dicha función como: 𝑓(𝑥𝑖)=𝑓|𝑥𝑖+(𝑥−𝑥𝑖)𝜕𝑓 𝜕𝑥�𝑥𝑖+(𝑥−𝑥𝑖)2 2! 𝜕2𝑓 𝜕𝑥2�𝑥𝑖+(𝑥−𝑥𝑖)3 3! 𝜕3𝑓 𝜕𝑥3�𝑥𝑖+ … Para obtener una diferencia finita centrada, nos conviene aproximar la función en un par de puntos 𝑥0±∆𝑥 2: 𝑓�𝑥0+∆𝑥 2�=𝑓|𝑥0+∆𝑥 2𝜕𝑓 𝜕𝑥�𝑥0+1 2! �∆𝑥 2�2𝜕2𝑓 𝜕𝑥2�𝑥𝑖+1 3! �∆𝑥 2�3𝜕3𝑓 𝜕𝑥3�𝑥𝑖+ … 𝑓�𝑥0−∆𝑥 2�=𝑓|𝑥0−∆𝑥 2𝜕𝑓 𝜕𝑥�𝑥0+1 2! �∆𝑥 2�2𝜕2𝑓 𝜕𝑥2�𝑥𝑖−1 3! �∆𝑥 2�3𝜕3𝑓 𝜕𝑥3�𝑥𝑖+ … Ahora obtenemos la diferencia de ambas aproximaciones, y como puede observarse los términos impares se anulan: 𝑓�𝑥0+∆𝑥 2�−𝑓�𝑥0−∆𝑥 2�=∆𝑥𝜕𝑓 𝜕𝑥�𝑥0+∆𝑥3 6𝜕3𝑓 𝜕𝑥3�𝑥𝑖+⋯ El siguiente paso sería dividir la expresión por ∆𝑥: 𝑓�𝑥0+∆𝑥 2�−𝑓�𝑥0−∆𝑥 2� ∆𝑥 =𝜕𝑓 𝜕𝑥�𝑥0+∆𝑥2 6𝜕3𝑓 𝜕𝑥3�𝑥𝑖+⋯ Propagación de ondas transversales generadas con fuerza de radiación 21 Despreciando los términos de orden superior podemos aproximar la primera derivada de una función en un punto dado, cometiendo un error de orden (∆𝑥)2, mediante: 𝜕𝑓 𝜕𝑥�𝑥0≈𝑓�𝑥0+∆𝑥 2�−𝑓�𝑥0−∆𝑥 2� ∆𝑥 En esta aproximación se basará todo el modelo numérico expuesto a continuación. II. 3. c. Simplificación del modelo, simetría axial Para este estudio se ha desarrollado un modelo en coordenadas cilíndricas (𝑟,𝑧,𝜃) y simplificado suponiendo simetría de revolución, por lo que se obtiene un modelo en coordenadas axisimétricas (𝑟,𝑧). Esta particularidad nos permitirá estudiar un dominio tridimensional sin necesidad de calcular 3 dimensiones exactamente, dado que al suponer simetría de revolución el cálculo es mucho menos costoso y más real que utilizando un modelo basado en coordenadas cartesianas y dos dimensiones. Particularizar el modelo para un sistema de coordenadas dado es peculiarmente sencillo, sólo hace falta utilizar los operadores diferenciales de nuestro sistema de coordenadas, en nuestro caso: �𝜕𝑥=𝜕 𝜕𝑥� ∇ϕ=𝜕𝑟𝜙 𝑟+𝜕𝜃𝜙 𝑟𝜃�+𝜕𝑧𝜙 𝑧 div F �  =∇· F �  =𝜕𝑟𝜕𝐹𝑟 𝜕𝑟+1 𝑟(𝐹𝑟+𝜕𝜃𝐹𝜃)+𝜕𝑧𝐹𝑧 El sistema de ecuaciones resultante sería: Esfuerzos normales: 𝜕𝑡𝜏𝑟𝑟=(𝜆+ 2𝜇)𝜕𝑟𝑣𝑟+𝜆� 𝜕𝑧𝑣𝑧+𝜕𝜃𝑣𝜃 𝑟+𝑣𝑟 𝑟� 𝜕𝑡𝜏𝑧𝑧=(𝜆+ 2𝜇)𝜕𝑧𝑣𝑧+𝜆�𝜕𝑟𝑣𝑟+𝑣𝑟 𝑟+𝜕𝜃𝑣𝜃 𝑟� 𝜕𝑡𝜏𝜃𝜃=(𝜆+ 2𝜇) �𝑣𝑟 𝑟+𝜕𝜃𝑣𝜃 𝑟�+𝜆(𝜕𝑟𝑣𝑟+𝜕𝑧𝑣𝑧) Esfuerzos de cizalla: 𝜕𝑡𝜏𝑟𝑧=𝜇 (𝜕𝑟𝑣𝑧+𝜕𝑧𝑣𝑟) 𝜕𝑡𝜏𝑟𝜃=𝜇 �𝜕𝑟𝑣𝜃+𝜕𝜃𝑣𝑟 𝑟−𝑣𝜃 𝑟� 𝜕𝑡𝜏𝜃𝑧=𝜇 �𝜕𝜃𝑣𝑧 𝑟+𝜕𝑧𝑣𝜃� Velocidad de partícula: 𝜕𝑡𝑣𝑟=1 𝜌�𝜕𝑟𝜏𝑟𝑟+𝜕𝜃𝜏𝑟𝜃 𝑟+𝜕𝑧𝜏𝑟𝑧+𝜏𝑟𝑟−𝜏𝜃𝜃 𝑟� 22 Propagación de ondas transversales mediante fuerza de radiación 𝜕𝑡𝑣𝑧=1 𝜌�𝜕𝑟𝜏𝑟𝑧+𝜕𝜃𝜏𝑧𝜃 𝑟+𝜕𝑧𝜏𝑧𝑧+𝜏𝑟𝑧 𝑟� 𝜕𝑡𝑣𝜃=1 𝜌�𝜕𝑟𝜏𝑟𝜃+𝜕𝜃𝜏𝜃𝜃 𝑟+𝜕𝑧𝜏𝜃𝑧+2𝜏𝑟𝜃 𝑟� Para introducir simetría de revolución se impondrá: 𝜕𝜃= 0 𝑣𝜃= 0 𝜏𝑟𝜃= 0 𝜏𝑧𝜃= 0 Esfuerzos normales: 𝜕𝑡𝜏𝑟𝑟=(𝜆+ 2𝜇)𝜕𝑟𝑣𝑟+𝜆� 𝜕𝑧𝑣𝑧+𝑣𝑟 𝑟� 𝜕𝑡𝜏𝑧𝑧=(𝜆+ 2𝜇)𝜕𝑧𝑣𝑧+𝜆�𝜕𝑟𝑣𝑟+𝑣𝑟 𝑟� 𝜕𝑡𝜏𝜃𝜃=(𝜆+ 2𝜇) 𝑣𝑟 𝑟+𝜆(𝜕𝑟𝑣𝑟+𝜕𝑧𝑣𝑧) Esfuerzo de cizalla: 𝜕𝑡𝜏𝑟𝑧=𝜇 (𝜕𝑟𝑣𝑧+𝜕𝑧𝑣𝑟) Velocidad de partícula: 𝜕𝑡𝑣𝑟=1 𝜌�𝜕𝑟𝜏𝑟𝑟+𝜕𝑧𝜏𝑟𝑧+𝜏𝑟𝑟−𝜏𝜃𝜃 𝑟� 𝜕𝑡𝑣𝑧=1 𝜌�𝜕𝑟𝜏𝑟𝑧+𝜕𝑧𝜏𝑧𝑧+𝜏𝑟𝑧 𝑟� II. 3. d. Absorción y PML Debido a que el dominio es finito y por economía de cálculo no puede ser suficientemente grande para no obtener reflexiones dentro de la zona de estudio y poder considerarlo como infinito, se implementa lo que se conoce como PML (Perfect Mach Layer). Esta técnica consiste en crear alrededor del dominio una zona donde se introduce absorción artificial anisótropa y el valor de ésta aumenta cuadráticamente. De esta manera se provoca que la onda se atenúe exponencialmente y por tanto no se produzcan reflexiones en los límites del dominio. Como ya se ha observado anteriormente cada nueva variable que se calcula depende a su vez de otras variables calculadas anteriormente con sus direcciones de propagación. Para añadir absorción en cada dirección será necesario atenuar cada término en función de si éste presenta gradiente en esa dirección. Por lo tanto será necesario descomponer las diferentes componentes de velocidad y esfuerzo en subcomponentes en función de la dirección de gradiente que presenten. El modelo resultante es: Propagación de ondas transversales generadas con fuerza de radiación 23 𝜌(𝛼𝑟𝜕𝑡+𝜔𝑟)𝑣𝑟(𝑟)=𝜕𝑟𝜏𝑟𝑟 𝜌(𝛼𝑧𝜕𝑡+𝜔𝑧)𝑣𝑟(𝑧)=𝜕𝑧𝜏𝑟𝑧 𝜌(Α𝑟𝜕𝑡+Ω𝑟)𝑣𝑟(𝜃)=𝜏𝑟𝑟−𝜏𝜃𝜃�𝑣𝑟=𝑣𝑟(𝑟)+𝑣𝑟(𝑧)+𝑣𝑟(𝜃) 𝜌(𝛼𝑟𝜕𝑡+𝜔𝑟)𝑣𝑧(𝑟)=𝜕𝑟𝜏𝑟𝑧 𝜌(𝛼𝑧𝜕𝑡+𝜔𝑧)𝑣𝑧(𝑧)=𝜕𝑧𝜏𝑧𝑧 𝜌(Α𝑟𝜕𝑡+Ω𝑟)𝑣𝑧(𝜃)=𝜏𝑟𝑧 �𝑣𝑧=𝑣𝑧(𝑟)+𝑣𝑧(𝑧)+𝑣𝑧(𝜃) (𝛼𝑟𝜕𝑡+𝜔𝑟)𝜏𝑟𝑟 (𝑟)=(𝜆+ 2𝜇)𝜕𝑟𝑣𝑟 (𝛼𝑧𝜕𝑡+𝜔𝑧)𝜏𝑟𝑟 (𝑧)=𝜆 𝜕𝑧𝑣𝑧 (Α𝑟𝜕𝑡+Ω𝑟)𝜏𝑟𝑟 (𝜃)=𝜆 𝑣𝑟�𝜏𝑟𝑟=𝜏𝑟𝑟 (𝑟)+𝜏𝑟𝑟 (𝑧)+𝜏𝑟𝑟 (𝜃) (𝛼𝑟𝜕𝑡+𝜔𝑟)𝜏𝑧𝑧 (𝑟)=𝜆 𝜕𝑟𝑣𝑟 (𝛼𝑧𝜕𝑡+𝜔𝑧)𝜏𝑧𝑧 (𝑧)=(𝜆+ 2𝜇)𝜕𝑧𝑣𝑧 (Α𝑟𝜕𝑡+Ω𝑟)𝜏𝑧𝑧 (𝜃)=𝜆 𝑣𝑟�𝜏𝑧𝑧=𝜏𝑧𝑧 (𝑟)+𝜏𝑧𝑧 (𝑧)+𝜏𝑧𝑧 (𝜃) (𝛼𝑟𝜕𝑡+𝜔𝑟)𝜏𝜃𝜃 (𝑟)=𝜆 𝜕𝑟𝑣𝑟 (𝛼𝑧𝜕𝑡+𝜔𝑧)𝜏𝜃𝜃 (𝑧)=𝜆 𝜕𝑧𝑣𝑧 (Α𝑟𝜕𝑡+Ω𝑟)𝜏𝜃𝜃 (𝜃)=(𝜆+ 2𝜇) 𝑣𝑟�𝜏𝜃𝜃=𝜏𝜃𝜃 (𝑟)+𝜏𝜃𝜃 (𝑧)+𝜏𝜃𝜃 (𝜃) (𝛼𝑟𝜕𝑡+𝜔𝑟)𝜏𝑟𝑧 (𝑟)=𝜇 𝜕𝑟𝑣𝑧 (𝛼𝑧𝜕𝑡+𝜔𝑧)𝜏𝑟𝑧 (𝑧)=𝜇 𝜕𝑧𝑣𝑟�𝜏𝑟𝑧=𝜏𝑟𝑧 (𝑟)+𝜏𝑟𝑧 (𝑧) De esta manera cuando estemos en la zona del dominio donde no queramos atenuar 𝛼= 1 y 𝜔= 0, mientras que en las zonas donde queramos conseguir atenuación iremos aumentando esos valores de forma exponencial como se muestra en la siguiente figura: 𝛼 𝜔 020 40 60 80 100 120 140 160 180 200 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 020 40 60 80 100 120 140 160 180 200 0 0.5 1 1.5 2 2.5 3 3.5 4x 10 6 24 Propagación de ondas transversales mediante fuerza de radiación II. 3. e. Discretización y tresbolillado del modelo El último paso es dicretizar nuestro modelo, para ello escogeremos en qué puntos tanto espaciales como temporales deseamos evaluar nuestro modelo. Según el desarrollo de Taylor para evaluar la derivada parcial, con respecto una variable, de una función en un punto de la variable se necesita conocer el valor de dicha función para esa variable en los puntos contiguos. Esto temporalmente implica que para calcular el instante actual se necesita el instante anterior y el instante siguiente, el instante futuro servirá para calcular el nuevo valor que toma la función y el instante pasado será conocido (ya sea una condición inicial, por ejemplo reposo, o el valor de la función calculada en el paso anterior). Esto implica que necesitaremos conocer el instante actual del resto de funciones que necesitamos para el cálculo. Debido a esta peculiaridad se necesita desplazar temporalmente unas funciones de otras, esto se conoce como tresbolillado. Espacialmente sucede el mismo caso, quedando así el mallado de la simulación: En esta representación espacial del mallado cabe imaginar que las variables de esfuerzo están desplazadas también temporalmente. Con estas consideraciones, así quedarían evaluadas las derivadas del modelo: 𝜕𝜏𝑟𝑟 𝜕𝑡�(𝑖,𝑗) 𝑛+12=(𝜆+ 2𝜇)𝜕𝑣𝑟 𝜕𝑟�(𝑖,𝑗) 𝑛+12+𝜆�𝜕𝑣𝑧 𝜕𝑧�(𝑖,𝑗) 𝑛+12+𝑣𝑟|(𝑖,𝑗) 𝑛+12 𝑟� 𝜕𝜏𝑧𝑧 𝜕𝑡�(𝑖,𝑗) 𝑛+12=(𝜆+ 2𝜇)𝜕𝑣𝑧 𝜕𝑧�(𝑖,𝑗) 𝑛+12+𝜆�𝜕𝑣𝑟 𝜕𝑟�(𝑖,𝑗) 𝑛+12+𝑣𝑟|(𝑖,𝑗) 𝑛+12 𝑟� 𝜕𝜏𝜃𝜃 𝜕𝑡�(𝑖,𝑗) 𝑛+12=(𝜆+ 2𝜇)𝑣𝑟|(𝑖,𝑗) 𝑛+12 𝑟+𝜆�𝜕𝑣𝑟 𝜕𝑟�(𝑖,𝑗) 𝑛+12+𝜕𝑣𝑧 𝜕𝑧�(𝑖,𝑗) 𝑛+12� 𝜕𝜏𝑟𝑧 𝜕𝑡�(𝑖+12,𝑗+12) 𝑛+12=𝜇 �𝜕𝑣𝑧 𝜕𝑟�(𝑖+12,𝑗+12) 𝑛+12+𝜕𝑣𝑟 𝜕𝑧�(𝑖+12,𝑗+12) 𝑛+12� 𝜕𝑣𝑟 𝜕𝑡�(𝑖+12,𝑗) 𝑛=1 𝜌�𝜕𝜏𝑟𝑟 𝜕𝑟�(𝑖+12,𝑗) 𝑛+𝜕𝜏𝑟𝑧 𝜕𝑧�(𝑖+12,𝑗) 𝑛+𝜏𝑟𝑟|(𝑖+12,𝑗) 𝑛−𝜏𝜃𝜃|(𝑖+12,𝑗) 𝑛 𝑟� 𝜕𝑣𝑧 𝜕𝑡�(𝑖,𝑗+12) 𝑛=1 𝜌�𝜕𝜏𝑧𝑧 𝜕𝑧�(𝑖,𝑗+12) 𝑛+𝜕𝜏𝑟𝑧 𝜕𝑟�(𝑖,𝑗+12) 𝑛+𝜏𝑟𝑧|(𝑖,𝑗+12) 𝑛𝑟� Propagación de ondas transversales generadas con fuerza de radiación 25 II. 3. f. Consideraciones del modelo El modelo de FDTD es condicionalmente estable, por lo que se deben cumplir ciertas condiciones para la solución no diverja: 𝑆≡ 𝑐𝑚𝑎𝑥·∆𝑡 �(∆𝑟)2+(∆𝑧)2≤1 Donde 𝑐𝑚𝑎𝑥 es la velocidad máxima de propagación dentro del dominio. Como en él se van a propagar tanto ondas longitudinales como trasversales, al ser las primeras las que más velocidad de propagación tienen: 𝑐𝐿=�𝜆+ 2𝜇 𝜌 Por otro lado, para reducir la dispersión numérica es necesario que en el espacio estén bien representadas las ondas que en él se propagan, para ello necesitaremos un buen muestreo espacial (𝑁𝜆), este dato vendrá dado por las veces que podemos muestrear la longitud de onda más pequeña: 𝑁𝜆=𝜆𝑚𝑖𝑛 ∆𝑚𝑎𝑥=𝑐𝑚𝑖𝑛 𝑓𝑚𝑎𝑥�(∆𝑟)2+(∆𝑧)2≥10 Este valor se recomienda que sea de 10 muestras por longitud de onda. Para este caso la velocidad mínima la marcarán las ondas transversales que son las que se propagan con mayor lentitud: 𝑐𝑠=�𝜇𝜌 32 Propagación de ondas transversales mediante fuerza de radiación Realizando el mismo procedimiento, pero ahora teniendo en cuenta que la ventana temporal tiene que estar correctamente escogida para escoger en cada caso el frente de ondas, obtenemos una estimación de la velocidad de propagación: IV. 4. FOCO DINÁMICO, TRANSMISIÓN OL+OS En los primeros instantes de propagación la fuente empieza a desplazarse en sentido del eje y comienza a perturbar el medio: 46810 12 14 x 10 -6 0 0.005 0.01 0.015 0.02 0.025 c S =513.3 m/s (r 2 =0.9872) c L =1972 (r 2 =0.9857) Propagación de ondas transversales generadas con fuerza de radiación 33 A medida que se propaga se empieza a discernir el cono de ondas transversales supersónicas causado por el movimiento de la fuente. En este caso el análisis de los receptores muestra la magnitud del frente de ondas supersónico, este frente de ondas aparece distorsionado debido a que la fuente está viajando mucho más rápido que la velocidad con la que el frente avanza. 00.5 11.5 2 x 10 -5 0 2 4 6 x 10 -3 0 0.01 0.02 0.03 0.04 r (m) t (s) v (m/s) 34 Propagación de ondas transversales mediante fuerza de radiación En este caso aún distorsionado el frente de ondas respecto a la señal original, todavía se puede asegurar cierta correlación entra las señales que captan los receptores. Pudiendo así realizar el ajuste lineal y estimar la velocidad de propagación del frente supersónico. Con estos datos se puede afirmar que el frente de ondas generado es de ondas transversales, dado que se propaga a la velocidad dichas ondas. 0.6 0.7 0.8 0.9 11.1 1.2 1.3 1.4 1.5 x 10 -5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 -3 t (s) r (m) c s =477.5 m/s r 2 =0.996 Propagación de ondas transversales generadas con fuerza de radiación 35 V. REFERENCIAS Auld B.A. Acoustic Fields and Elastic Waves in Solids. New York. Wiley, 1973. Vol. vols. I and II. Bastard Cécile, Remeniéras Jean-Pierre y Callé Samuel Simulation of shear wave propagation in a soft medium using a pseudospectral time domain method. J. Acoust. Soc. Am. 2009. 126:4 . Bayliss A., Gunzburger M. y Turkel E. Boundary conditions for the numerical solution of elliptic equations in exterior regions. J. Appl. Math.. 1982. 42:430-451. Berenger J.-P. A perfectly matched layer for the absorption of electromagnetics waves. J. Comput. Phys. 1994. 114:185–200. Jiménez González, N. Simulación de tejidos vegetales mediante diferencias finitas. Tesis de Máster, Universidad Politècnica de València. Maleke C., Konofagou E E. Harmonic motion imaging for focused ultrasound (HMIFU): a fully integrated technique for sonication and monitoring of thermal ablation in tissues.. Phys. Med. Biol. 53 (2008) 1773– 1793 Salido González, N. Análisis Numérico de parámetros de calidad en haces focalizados. Trabajo final de Carrera, Universidad Politècnica de València. 36 Propagación de ondas transversales mediante fuerza de radiación ANEXO I. ECUACIONES DEL MODELO NUMÉRICO Considerando todas las aportaciones discutidas en el modelo teórico, las ecuaciones resultantes listas para implementar son las siguientes: 𝑣𝑟_𝑟 𝑛+1�𝑖+1 2,𝑗�=��1−𝜔𝑟𝑑𝑡 2𝛼𝑟�𝑣𝑟_𝑟 𝑛�𝑖+1 2,𝑗�+1 𝜌�𝑑𝑡 𝛼𝑟�𝜕𝑟𝜏𝑟𝑟|𝑖+12,𝑗� �1+𝜔𝑟𝑑𝑡 2𝛼𝑟�� 𝑣𝑟_𝑧 𝑛+1�𝑖+1 2,𝑗�=��1−𝜔𝑧𝑑𝑡 2𝛼𝑧�𝑣𝑟_𝑧 𝑛�𝑖+1 2,𝑗�+1 𝜌�𝑑𝑡 𝛼𝑧�𝜕𝑧𝜏𝑟𝑧|𝑖+12,𝑗� �1 + 𝜔𝑧𝑑𝑡 2𝛼𝑧�� 𝑣𝑟_𝜃 𝑛+1�𝑖+1 2,𝑗�=��1−Ω𝑟𝑑𝑡 2Α𝑟�𝑣𝑟_𝜃 𝑛�𝑖+1 2,𝑗�+1 𝜌�𝑑𝑡 Α𝑟�(𝜏𝑟𝑟−𝜏𝜃𝜃)|𝑖+12,𝑗� �1 + Ω𝑟𝑑𝑡 2Α𝑟�� 𝑣𝑧_𝑟 𝑛+1�𝑖,𝑗+1 2�=��1−𝜔𝑟𝑑𝑡 2𝛼𝑟�𝑣𝑧_𝑟 𝑛�𝑖,𝑗+1 2�+1 𝜌�𝑑𝑡 𝛼𝑟�𝜕𝑟𝜏𝑟𝑧|𝑖,𝑗+12� �1 + 𝜔𝑟𝑑𝑡 2𝛼𝑟�� 𝑣𝑧_𝑧 𝑛+1�𝑖,𝑗+1 2�=��1−𝜔𝑧𝑑𝑡 2𝛼𝑧�𝑣𝑧_𝑧 𝑛�𝑖,𝑗+1 2�+1 𝜌�𝑑𝑡 𝛼𝑧�𝜕𝑧𝜏𝑧𝑧|𝑖,𝑗+12� �1 + 𝜔𝑧𝑑𝑡 2𝛼𝑧�� 𝑣𝑧_𝜃 𝑛+1�𝑖,𝑗+1 2�=��1−Ω𝑟𝑑𝑡 2Α𝑟�𝑣𝑧_𝜃 𝑛�𝑖,𝑗+1 2�+1 𝜌�𝑑𝑡 Α𝑟�(𝜏𝑟𝑧)|𝑖,𝑗+12� �1+Ω𝑟𝑑𝑡 2Α𝑟�� 𝜏𝑟𝑟_𝑟 𝑛+1(𝑖,𝑗)=��1−𝜔𝑟𝑑𝑡 2𝛼𝑟�𝜏𝑟𝑟_𝑟 𝑛(𝑖,𝑗)+(𝜆+ 2𝜇)�𝑑𝑡 𝛼𝑟�𝜕𝑟𝑣𝑟|𝑖,𝑗� �1+𝜔𝑟𝑑𝑡 2𝛼𝑟�� 𝜏𝑟𝑟_𝑧 𝑛+1(𝑖,𝑗)=��1−𝜔𝑧𝑑𝑡 2𝛼𝑧�𝜏𝑟𝑟_𝑧 𝑛(𝑖,𝑗)+𝜆�𝑑𝑡 𝛼𝑧�𝜕𝑧𝑣𝑧|𝑖,𝑗� �1 + 𝜔𝑧𝑑𝑡 2𝛼𝑧�� 𝜏𝑟𝑟_𝜃 𝑛+1(𝑖,𝑗)=��1−Ω𝑟𝑑𝑡 2Α𝑟�𝜏𝑟𝑟_𝜃 𝑛(𝑖,𝑗)+𝜆�𝑑𝑡 Α𝑟�(𝑣𝑟)|𝑖,𝑗� �1 + Ω𝑟𝑑𝑡 2Α𝑟�� 𝜏𝑧𝑧_𝑟 𝑛+1(𝑖,𝑗)=��1−𝜔𝑟𝑑𝑡 2𝛼𝑟�𝜏𝑧𝑧_𝑟 𝑛(𝑖,𝑗)+𝜆�𝑑𝑡 𝛼𝑟�𝜕𝑟𝑣𝑟|𝑖,𝑗� �1 + 𝜔𝑟𝑑𝑡 2𝛼𝑟�� 𝜏𝑧𝑧_𝑧 𝑛+1(𝑖,𝑗)=��1−𝜔𝑧𝑑𝑡 2𝛼𝑧�𝜏𝑧𝑧𝑟_𝑧 𝑛(𝑖,𝑗)+(𝜆+ 2𝜇)�𝑑𝑡 𝛼𝑧�𝜕𝑧𝑣𝑧|𝑖,𝑗� �1 + 𝜔𝑧𝑑𝑡 2𝛼𝑧�� 𝜏𝑧𝑧_𝜃 𝑛+1(𝑖,𝑗)=��1−Ω𝑟𝑑𝑡 2Α𝑟�𝜏𝑧𝑧_𝜃 𝑛(𝑖,𝑗)+𝜆�𝑑𝑡 Α𝑟�(𝑣𝑟)|𝑖,𝑗� �1+Ω𝑟𝑑𝑡 2Α𝑟�� 𝜏𝜃𝜃_𝑟 𝑛+1(𝑖,𝑗)=��1−𝜔𝑟𝑑𝑡 2𝛼𝑟�𝜏𝜃𝜃_𝑟 𝑛(𝑖,𝑗)+𝜆�𝑑𝑡 𝛼𝑟�𝜕𝑟𝑣𝑟|𝑖,𝑗� �1+𝜔𝑟𝑑𝑡 2𝛼𝑟�� 𝜏𝜃𝜃_𝑧 𝑛+1(𝑖,𝑗)=��1−𝜔𝑧𝑑𝑡 2𝛼𝑧�𝜏𝜃𝜃_𝑧 𝑛(𝑖,𝑗)+𝜆�𝑑𝑡 𝛼𝑧�𝜕𝑧𝑣𝑧|𝑖,𝑗� �1 + 𝜔𝑧𝑑𝑡 2𝛼𝑧�� 𝜏𝜃𝜃_𝜃 𝑛+1(𝑖,𝑗)=��1−Ω𝑟𝑑𝑡 2Α𝑟�𝜏𝜃𝜃_𝜃 𝑛(𝑖,𝑗)+(𝜆+ 2𝜇)�𝑑𝑡 Α𝑟�(𝑣𝑟)|𝑖,𝑗� �1 + Ω𝑟𝑑𝑡 2Α𝑟�� 𝜏𝑟𝑧_𝑟 𝑛+1�𝑖+1 2,𝑗+1 2�=��1−𝜔𝑟𝑑𝑡 2𝛼𝑟�𝜏𝑟𝑧_𝑟 𝑛�𝑖+1 2,𝑗+1 2�+𝜇�𝑑𝑡 𝛼𝑟�𝜕𝑟𝑣𝑧|𝑖+12,𝑗+12� �1+𝜔𝑟𝑑𝑡 2𝛼𝑟�� 𝜏𝑟𝑧_𝑧 𝑛+1�𝑖+1 2,𝑗+1 2�=��1−𝜔𝑧𝑑𝑡 2𝛼𝑧�𝜏𝑟𝑧_𝑧 𝑛�𝑖+1 2,𝑗+1 2�+𝜇�𝑑𝑡 𝛼𝑧�𝜕𝑧𝑣𝑟|𝑖+12,𝑗+12� �1 + 𝜔𝑧𝑑𝑡 2𝛼𝑧�� Propagación de ondas transversales generadas con fuerza de radiación 37