Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow
Abstract
Departamento de Ingeniería Energética y Fluidomecánica
Full text
UNIVERSIDAD DE VALLADOLID ESCUELA DE INGENIERÍAS INDUSTRIALES GRADO EN INGENIERÍA EN TECNOLOGÍAS INDUSTRIALES DESARROLLO DE UN MODELO DE DINÁMICA DE FLUIDOS COMPUTACIONAL PARA LA CORRECCIÓN DE FLUJOS MEDIDOS POR RESONANCIA MAGNÉTICA - 4D FLOW Autor: Serrano Holgado, Ricardo Tutor: Sierra Pallares, José Ingeniería Energética y Fluidomecánica Valladolid, junio 2023
Resumen Este trabajo trata sobre la complementación de las mediciones de la velocidad de un fluido llevadas a cabo mediante 4D Flow utilizando Dinámica de Fluidos Computacional. Para ello, se ha recopilado información sobre los fundamenentos de 4D Flow, de la mecánica de fluidos y de la Dinámica de Fluidos Computacional. Después, se han seleccionado los modelos numéricos adecuados y se han realizado simulaciones del flujo sanguíneo en un phantom de arteria aorta, así como de trayectorias de partículas en el phantom, utilizando el software Ansys Fluent. Por último, se han analizado los resultados de varias propiedades del fluido y de las trayectorias de las partículas, y se han indicado las conclusiones del trabajo. Palabras clave 4D Flow, Dinámica de Fluidos Computacional (CFD), Ansys Fluent, phantom, arteria aorta. Abstract This paper is about the complementation of the velocity measurements of a fluid carried out by 4D Flow using Computational Fluid Dynamics. To achieve this, information about the fundamentals of 4D Flow, fluid mechanics and Computational Fluid Dynamics was gathered. Then, suitable numerical models were selected, and simulations of blood flow inside a phantom of the aorta artery, as well as of particle trajectories inside the phantom, were performed using the software Ansys Fluent. Lastly, the results of various fluid properties and of the particle trajectories were analysed, and the paper’s conclusions were stated. Keywords 4D Flow, Computational Fluid Dynamics (CFD), Ansys Fluent, phantom, aorta artery.
Agradecimientos A todas las personas que me han ayudado a llegar hasta aquí. En especial, a mi tutor, José, sin cuya ayuda no habría podido realizar este trabajo. Y sobre todo, a mis padres, por la confianza que depositan en mí y el apoyo que me ofrecen en todas mis decisiones.
Índice de Figuras I UNIVERSIDAD DE VALLADOLID Índice Índice de Figuras .................................................................................................... II Índice de Tablas ................................................................................................... VII Índice de Ecuaciones .......................................................................................... VIII Nomenclatura ....................................................................................................... XII Acrónimos ......................................................................................................... XII Símbolos .......................................................................................................... XIII 1 Introducción y objetivos ................................................................................. 1 1.1 Estado del arte......................................................................................... 1 1.1.1 Imagen por Resonancia Magnética ................................................ 1 1.2 Objetivos ................................................................................................ 18 1.3 Estructura del TFG ................................................................................. 18 2 Desarrollo del TFG ........................................................................................ 19 2.1 Modelo fluidodinámico .......................................................................... 19 2.1.1 Ecuaciones de conservación ......................................................... 19 2.1.2 Propiedades físicas ........................................................................ 25 2.2 Modelo numérico ................................................................................... 27 2.2.1 Ventajas de la CFD ......................................................................... 27 2.2.2 Métodos numéricos frente a otros métodos ................................ 28 2.2.3 Pasos para resolver un problema por CFD ................................... 29 2.2.4 Aspectos numéricos de la CFD ...................................................... 30 2.2.5 Condiciones de contorno e iniciales ............................................. 41 2.2.6 Criterios de parada ........................................................................ 42 2.2.7 Detección de vórtices .................................................................... 43 2.3 Caso de estudio ..................................................................................... 45 2.3.1 Caso bidimensional ........................................................................ 46 2.3.2 Caso tridimensional ....................................................................... 52 2.3.3 Resumen del modelo aplicado al caso de estudio ...................... 56 2.4 Resultados ............................................................................................. 57 2.4.1 Simulación bidimensional ............................................................. 58 2.4.2 Simulación tridimensional ............................................................. 82 3 Conclusiones ............................................................................................. 107 4 Bibliografía ................................................................................................. 108
II Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Índice de Figuras Figura 1. Representación de un escáner de MRI [4]. .......................................... 2 Figura 2. Imagen obtenida por MRI [4]. ................................................................ 2 Figura 3. Precesión del momento magnético de los núcleos atómicos ante un campo B0 [6]. .......................................................................................................... 5 Figura 4. Esfera de Bloch ante la aplicación del campo B0 [6]. .......................... 7 Figura 5. Esfera de Bloch ante la aplicación de los campos B0 y B1 [6]. ............ 8 Figura 6. Esfera Bloch ante la aplicación de los campos B0 y B1 en un sistema de referencia móvil [6]. ....................................................................................... 10 Figura 7. Desfases típicos de la densidad de magnetización [6]. ................... 10 Figura 8. Tiempos de relajación de la densidad de magnetización [8]. .......... 12 Figura 9. Campo magnético oscilante linealmente polarizado [6]. ................. 14 Figura 10. Celda estructurada [11]. ................................................................... 31 Figura 11. Geometría de estudio tridimensional. ............................................. 45 Figura 12. Geometría de estudio bidimensional. .............................................. 46 Figura 13. División en áreas para el mallado de la geometría bidimensional. 47 Figura 14. Nomenclatura de las secciones del contorno de la geometría bidimensional. ..................................................................................................... 47 Figura 15. Velocidad en la entrada. ................................................................... 50 Figura 16. División en áreas para el mallado de la geometría tridimensional. .............................................................................................................................. 53 Figura 17. Punto en el que se miden las velocidades. ..................................... 55 Figura 18. Nomenclatura de las secciones del contorno de la geometría tridimensional. ..................................................................................................... 56 Figura 19. Geometría bidimensional. Campo de velocidades. Vista general. t=19 s. .................................................................................................................. 58 Figura 20. Geometría bidimensional. Campo de velocidades. Vista general. t=19,1 s. ............................................................................................................... 59 Figura 21. Geometría bidimensional. Campo de velocidades. Vista general. t=19,2 s. ............................................................................................................... 59 Figura 22. Geometría bidimensional. Campo de velocidades. Vista general. t=19,3 s. ............................................................................................................... 60 Figura 23. Geometría bidimensional. Campo de velocidades. Vista general. t=19,4 s. ............................................................................................................... 60
Índice de Figuras III UNIVERSIDAD DE VALLADOLID Figura 24. Geometría bidimensional. Campo de velocidades. Vista general. t=19,5 s. ............................................................................................................... 61 Figura 25. Geometría bidimensional. Campo de velocidades. Vista general. t=19,6 s. ............................................................................................................... 61 Figura 26. Geometría bidimensional. Campo de velocidades. Vista general. t=19,7 s. ............................................................................................................... 62 Figura 27. Geometría bidimensional. Campo de velocidades. Vista general. t=19,8 s. ............................................................................................................... 62 Figura 28. Geometría bidimensional. Campo de velocidades. Vista general. t=19,9 s. ............................................................................................................... 63 Figura 29. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19 s. ............................................................................................................................ 63 Figura 30. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,1 s. ............................................................................................................... 64 Figura 31. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,2 s. ............................................................................................................... 64 Figura 32. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,3 s. ............................................................................................................... 65 Figura 33. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,4 s. ............................................................................................................... 65 Figura 34. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,5 s. ............................................................................................................... 66 Figura 35. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,6 s. ............................................................................................................... 66 Figura 36. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,7 s. ............................................................................................................... 67 Figura 37. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,8 s. ............................................................................................................... 67 Figura 38. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,9 s. ............................................................................................................... 68 Figura 39. Geometría bidimensional. Contorno de velocidades. t=19 s. ........ 69 Figura 40. Geometría bidimensional. Contorno de velocidades. t=19,1 s. ..... 70 Figura 41. Geometría bidimensional. Contorno de velocidades. t=19,2 s. ..... 70 Figura 42. Geometría bidimensional. Contorno de velocidades. t=19,3 s. ..... 71 Figura 43. Geometría bidimensional. Contorno de velocidades. t=19,4 s. ..... 71 Figura 44. Geometría bidimensional. Contorno de velocidades. t=19,5 s. ..... 72
X Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Ecuación 42. Aplicación del teorema de Gauss-Ostrogradski al término de convección [11]. .................................................................................................. 24 Ecuación 43. Aplicación del teorema de Gauss-Ostrogradski al término de difusión [11]. ........................................................................................................ 24 Ecuación 44. Ecuación integral de conservación de una propiedad genérica del fluido. Forma 2 [11]............................................................................................. 24 Ecuación 45. Ecuación integral de conservación de la masa [12]. ................. 25 Ecuación 46. Ecuaciones integrales de conservación de la cantidad de movimiento [12]. ................................................................................................. 25 Ecuación 47. Ley de los Gases Ideales [11]. ..................................................... 25 Ecuación 48. Relación entre la densidad y la temperatura en líquidos [11]. . 26 Ecuación 49. Esfuerzos tangenciales en fluidos newtonianos [11]. ............... 26 Ecuación 50. Ley constitutiva de los fluidos de Carreau [17]. ......................... 26 Ecuación 51. Ley constitutiva de los fluidos de Casson [19]. .......................... 27 Ecuación 52. Integración del término de acumulación en una celda estructurada [11]................................................................................................. 31 Ecuación 53. Integración del término de convección en una celda estructurada [11]. ...................................................................................................................... 32 Ecuación 54. Integración del término de difusión en una celda estructurada [11]. ...................................................................................................................... 32 Ecuación 55. Integración del término de generación en una celda estructurada [11]. ...................................................................................................................... 32 Ecuación 56. Esquema de diferencias centradas [11]. ................................... 33 Ecuación 57. Estimación del gradiente [11]. .................................................... 33 Ecuación 58. Esquema upwind de primer orden [11]. ..................................... 34 Ecuación 59. Esquema upwind de segundo orden [11]. .................................. 34 Ecuación 60. Esquema van Leer [11]. ............................................................... 35 Ecuación 61. Ecuación para la presión [11]. .................................................... 35 Ecuación 62. Sub-relajación [11]. ...................................................................... 36 Ecuación 63. Integración espacial y temporal de la ecuación de conservación de una propiedad genérica del fluido [11]. ....................................................... 37
Índice de Ecuaciones XI UNIVERSIDAD DE VALLADOLID Ecuación 64. Factor e peso en la integración temporal [11]. .......................... 37 Ecuación 65. Condición de delimitación del método explícito [11]. ................ 38 Ecuación 66. Condición de delimitación del método de Crank-Nicolson [11]. 38 Ecuación 67. Condición de Courant-Friedrichs-Lewy [20]. ............................... 38 Ecuación 68. Número de Courant [20]. ............................................................. 39 Ecuación 69. Número de Courant generalizado para n dimensiones [20]. .... 39 Ecuación 70. Índice de convergencia en caso de utilizar la malla más densa [21]. ....................................................................................................................... 40 Ecuación 71. Ratio de refinado de malla [21]. .................................................. 40 Ecuación 72. Índice de convergencia en caso de utilizar la malla menos densa [21]. ....................................................................................................................... 41 Ecuación 73. GCI asintótico [21]. ....................................................................... 41 Ecuación 74. Tensor gradiente de velocidad. .................................................... 43 Ecuación 75. Tensor de velocidades de deformación [23]. ............................. 43 Ecuación 76. Tensor de velocidades de rotación [23]. ..................................... 44 Ecuación 77. Invariante Q [24]. .......................................................................... 44 Ecuación 78. Tensor S2 + Ω2 [24]. ...................................................................... 44 Ecuación 79. Número de Reynolds [25]. ........................................................... 48 Ecuación 80. Condición de régimen laminar para flujo interno en tuberías de sección circular [25]. ........................................................................................... 48 Ecuación 81. Cálculo de la viscosidad dinámica. .............................................. 49 Ecuación 82. Velocidad en la entrada. ............................................................... 49 Ecuación 83. Cálculo del número de Courant para el área general de la malla. ............................................................................................................................... 52 Ecuación 84. Cálculo del número de Courant para el área afectada por las esferas de influencia de la malla. ....................................................................... 52
XII Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Nomenclatura Acrónimos Acrónimo Descripción MRI NMR MRA PC-MRA PC-MRI CFD CAD GCI Magnetic Resonance Imaging Imagen por Resonancia Magnética Nuclear Magnetic Resonance Resonancia Magnética Nuclear Magnetic Resonance Angiography Angiografía por Resonancia Magnética Phase Contrast Magnetic Resonance Angiography Angiografía por Resonancia Magnética por Contraste de Fase Phase Contrast Magnetic Resonance Imaging Imagen por Resonancia Magnética por Contraste de Fase Computational Fluid Dynamics Dinámica de Fluidos Computacional Computer-Aided Design Diseño Asistido por Ordenador Grid Convergence Index Índice de Convergencia de Malla
Nomenclatura XIII UNIVERSIDAD DE VALLADOLID Símbolos Símbolo Descripción t 𝜌 Tiempo Densidad U T P F A V E R 𝜇 𝜈 Velocidad Temperatura Presión Fuerza Área Volumen Energía Constante de los gases ideales Viscosidad dinámica Viscosidad cinemática
XIV Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID
1 Introducción y objetivos 1 UNIVERSIDAD DE VALLADOLID 1 Introducción y objetivos La MRI (Magnetic Resonance Imaging) es una de las técnicas más avanzadas y efectivas para la obtención de imágenes de los tejidos del cuerpo humano y la realización de diagnósticos de forma no invasiva [1]. La MRI comprende una amplia variedad de técnicas. Una de ellas, denominada PC-MRA (Phase Contrast Magnetic Resonance Angiography) se basa en el hecho de que algunos tejidos, como la sangre, se encuentran en movimiento [2]. Esta técnica es adecuada para la detección de problemas cardiovasculares. La PC-MRA presenta algunos problemas como excesivo ruido, baja resolución o incapacidad para determinar parámetros fluidodinámicos de la sangre aparte de su velocidad que pueden ser de importancia para el diagnóstico [3]. En caso de que el paciente posea algún dispositivo interno metálico como un stent o un marcapasos, estos problemas pueden verse agravados, ya que el dispositivo puede interferir con los campos magnéticos del escáner MRI, provocando una falta de fiabilidad en la información obtenida en la zona donde se sitúa el dispositivo. Una forma de solucionar estos problemas es la complementación de la PC-MRA con la CFD (Computational Fluid Dynamics). La CFD presenta fortalezas en los puntos más débiles de la PC-MRA, lo que hace que la combinaicón de ambas técnicas, si se hace correctamente, proporcione una información precisa, fiable y completa [3]. 1.1 Estado del arte 1.1.1 Imagen por Resonancia Magnética La Imagen por Resonancia Magnética (MRI por sus siglas en inglés, Magnetic Resonance Imaging) es un procedimiento médico no invasivo que permite obtener imágenes de alta resolución de tejidos internos del cuerpo [1]. Para llevarla a cabo, el sujeto debe primero deshacerse de todos los objetos metálicos que porte, ya que pueden producir daños tanto humanos como materiales y alterar el resultado de la prueba. En caso de que el médico lo considere necesario, se le inyecta al sujeto una sustancia conocida como medio de contraste, basada en gadolinio, que sirve para mejorar la resolución de la imagen final. El sujeto debe acostarse en una mesa motorizada que lo introduce dentro de un gran escáner cilíndrico. El escáner produce una serie de campos magnéticos
2 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID ante los que los distintos tejidos del cuerpo del sujeto reaccionan de formas diferentes, dando lugar a distintas señales que el escáner detecta. Se pueden configurar los parámetros del escáner para que en la imagen final aparezcan resaltados uno o varios tejidos. Con esta información, se genera por ordenador la imagen final [4]. Figura 1. Representación de un escáner de MRI [4]. Figura 2. Imagen obtenida por MRI [4]. Si el sujeto tiene algún dispositivo metálico interno, como un stent o un marcapasos, se debe atender a ciertas consideraciones especiales. Se debe consultar cuáles son las propiedades magnéticas del dispositivo y cuál es la
1 Introducción y objetivos 3 UNIVERSIDAD DE VALLADOLID recomendación del fabricante respecto a la realización de la MRI. Actualmente la mayoría de los dispositivos que se implantan tienen un carácter no ferromagnético o débilmente ferromagnético y es posible realizar la MRI sin riesgo [5], aunque si la zona estudiada es cercana al dispositivo es posible que los resultados de la prueba se vean alterados. Resonancia Magnética Nuclear La MRI se fundamenta en el fenómeno físico de Resonancia Magnética Nuclear (NMR por sus siglas en inglés, Nuclear Magnetic Resonance). Este fenómeno es de naturaleza cuántica. Sin embargo, al estudiar objetos macroscópicos, el elevado número de núcleos atómicos involucrados hace que las particularidades cuánticas pierdan relevancia. Por ello, es habitual describir la NMR mediante una analogía con un fenómeno de naturaleza clásica: un objeto en rotación dentro de un campo magnético [6]. Aun así, es importante no olvidar que la descripción clásica no es más que una analogía, y no describe correctamente los fenómenos que tienen lugar a escala microscópica. A continuación se desarrollan ambas perspectivas. NMR según la mecánica cuántica Las partículas como los electrones, protones o neutrones poseen una propiedad denominada espín (𝑆), que está asociada con un momento angular intrínseco. Este momento angular no tiene nada que ver con una rotación. El espín es una cantidad vectorial que presenta algunas particularidades: su módulo es fijo, y está limitada a un conjunto de posibles orientaciones. Para una partícula de espín ½ como los electrones, protones o neutrones, el espín solo puede tomar dos valores: ±ℏ 2, donde ℏ es la constante de Plank dividida entre 2π [6]. Las partículas con espín poseen un momento magnético de valor 𝑚 =𝛾ℏ𝑆 Ecuación 1. Momento magnético cuántico de una partícula [6].
4 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID donde 𝛾 es el ratio giromagnético, que es característico de cada partícula. Los espines de los nucleones que forman un núcleo atómico se combinan cuánticamente para dar lugar al espín total del núcleo 𝐼. Los átomos que poseen un número impar de protones o neutrones tienen espín no nulo y por tanto tienen momento magnético no nulo [6]. Cuando un núcleo con momento magnético no nulo es sometido a un campo magnético externo 𝐵0 se produce una interacción que hace que la energía que posee el núcleo varíe. El espín del núcleo puede adquirir una componente de +ℏ 2 o −ℏ 2 en la dirección del campo magnético externo. A estos dos estados se les llama estados estacionarios y entre ellos existe una diferencia de energía de Δ𝐸=−𝛾ℏ𝐵0 Ecuación 2. Diferencia de energía entre los estados estacionarios [6]. Si ahora se somete al núcleo a un segundo campo magnético externo 𝐵1 perpendicular a 𝐵0 y que rota a su alrededor con una frecuencia angular de 𝜔0=|𝛾𝐵0| Ecuación 3. Frecuencia angular del campo B1 [6]. el espín del núcleo comienza a oscilar a través de una serie de estados de superposición cuántica entre los dos estados estacionarios. Si el campo magnético externo 𝐵1 se retira cuando el espín se encuentra en un estado intermedio, sus componentes transversales (las perpendiculares a la dirección de 𝐵0 ) continúan oscilando entre sus posibles valores a una frecuencia 𝜔0 [6].
1 Introducción y objetivos 5 UNIVERSIDAD DE VALLADOLID NMR según la mecánica clásica La explicación clásica del fenómeno de NMR parte de considerar el espín de las partículas y núcleos atómicos como una rotación real. Esta rotación tiene asociada un momento angular (𝐽) que causa un momento magnético (𝑚 ), relacionados según la siguiente expresión 𝑚 =𝛾𝐽 Ecuación 4. Momento magnético clásico de una partícula [6]. En ausencia de manipulaciones externas, el vector momento magnético presenta una orientación aleatoria. Sin embargo, al someter un núcleo atómico con momento magnético no nulo a un campo magnético externo 𝐵0 , su momento magnético describe un movimiento de precesión alrededor de la dirección del campo 𝐵0 , bien en sentido paralelo al mismo o bien en sentido antiparalelo [6]. Figura 3. Precesión del momento magnético de los núcleos atómicos ante un campo B0 [6]. Sumando los momentos magnéticos de todos los núcleos atómicos presentes en una determinada región del espacio se obtiene nueva magnitud de carácter macroscópico denominada densidad de magnetización (𝑀 ), que representa el momento magnético por unidad de volumen [6]. Cada grupo de núcleos
12 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Figura 8. Tiempos de relajación de la densidad de magnetización [8]. Durante la relajación, el módulo de 𝑀 no permanece constante. Agrupando los términos que describen la precesión de 𝑀 y los términos que describen su relajación se obtienen las ecuaciones completas del movimiento de 𝑀 , denominadas ecuaciones de Bloch [6]. 𝑑𝑀 𝑑𝑡=𝛾𝑀 ×𝐵 −𝑀⊥ 𝑇2−𝑀𝑧−𝑀0 𝑇1 Ecuación 18. Ecuaciones de Bloch [6]. Generación de los campos magnéticos Para que el fenómeno de NMR se produzca es necesaria la generación de varios campos magnéticos controlados. Los escáneres de MRI disponen de conjuntos de bobinas que, en combinación con materiales magnéticos y magnetizables, producen estos campos magnéticos al hacer circular corrientes eléctricas por ellas, según describe la ley de Ampère [6].
1 Introducción y objetivos 13 UNIVERSIDAD DE VALLADOLID Campo estático Para generar el campo 𝐵0 se utiliza un gran solenoide a través del cual se hace circular una corriente eléctrica. La expresión que define el campo generado es 𝐵0=𝜇0𝑛𝐼 √1+(𝐷 𝐿)2 Ecuación 19. Módulo del campo B0 [6]. donde 𝜇0 es la permeabilidad magnética del vacío, 𝜇0=4𝜋∙10−7𝑇∙𝑚 𝐴, n es el numero de espiras de la bobina, I es la intensidad de corriente que circula por la bobina, D es el diámetro de la bobina y L es la longitud de la bobina. En la práctica, el campo magnético producido por el solenoide no tiene la homogeneidad que se necesita para conseguir una imagen de alta resolución. Por ello, los escáneres de MRI disponen también de una serie de bobinas de corrección que ayudan a conseguir la uniformidad requerida [6]. Campos gradientes Para distinguir de qué punto del espacio provienen las señales de detección de NMR, es necesario que señales provenientes de regiones distintas tengan intensidades distintas. Para ello, cada punto de la región de estudio del escáner MRI debe estar sometido a una intensidad de campo magnético distinta. Para conseguir esto, se superponen al campo estático los llamados campos gradientes: campos magnéticos que varían en función de la posición. Los escáneres de MRI disponen de conjuntos de bobinas destinados a este propósito [6]. Campo de radiofrecuencia El campo rotatorio 𝐵1 , también llamado campo de radiofrecuencia o de RF debido a que la frecuencia a la que rota pertenece a este rango del espectro electromagnético, se obtiene a partir de un campo magnético oscilante linealmente polarizado. Este campo se puede descomponer en dos componentes de módulos iguales y constantes y que rotan con velocidades
14 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID angulares iguales en módulo pero de sentido opuesto, como se observa en la Figura 9 [6]. Figura 9. Campo magnético oscilante linealmente polarizado [6]. Una de estas componentes será el campo 𝐵1 . La otra tiene una influencia despreciable sobre la NMR, ya que no gira a la frecuencia de Larmor, y por tanto puede ser ignorada. Para producir el campo magnético oscilante linealmente polarizado se hace circular por una bobina dedicada a este propósito una corriente eléctrica que varía con el tiempo según la expresión 𝑖(𝑡)=𝑖0cos (𝜔0𝑡) Ecuación 20. Intensidad de corriente aplicada a la bobina [6]. Detección La densidad de magnetización está asociada a un pequeño campo magnético. Las variaciones de este campo magnético ante las interacciones con los campos 𝐵0 y 𝐵1 previamente explicadas provocan variaciones en el flujo magnético que atraviesa los conjuntos de bobinas de detección que los escáneres de MRI poseen. Esto, según las leyes de Faraday y Lenz de inducción electromagnética, produce en las bobinas una fuerza electromotriz que se opone a la variación de flujo magnético que las atraviesa.
1 Introducción y objetivos 15 UNIVERSIDAD DE VALLADOLID Normalmente, esta fuerza electromotriz es muy pequeña, y debe ser amplificada tan pronto como sea posible para evitar que la señal se vea alterada debido a interferencias o ruidos externos [6]. Angiografía por Resonancia Magnética Se denomina Angiografía por Resonancia Magnética (MRA por sus siglas en inglés, Magnetic Resonance Angiography) al conjunto de técnicas de MRI aplicadas a la obtención de imágenes de los vasos sanguíneos, especialmente de las arterias. La MRA es útil para detectar aneurismas, estenosis, trombos y otras afecciones del sistema circulatorio, y también para obtener mapas de velocidad del flujo sanguíneo [9]. La MRA comprende técnicas que se pueden dividir en dos grandes familias: MRA dependiente del flujo, que se basa en el hecho de que la sangre está en movimiento para distinguir los vasos sanguíneos de otros tejidos; y MRA independiente del flujo, que no tiene en consideración este hecho. De las técnicas de MRA dependiente del flujo, la más extendida y relevante para este trabajo se denomina Angiografía por Contraste de Fase (PC-MRA por sus siglas en inglés, Phase Contrast Magnetic Resonance Angiography). La PC-MRA, también denominada PC-MRI, consiste en la medición indirecta de la velocidad de una partícula a través de la medición de la variación que se produce en la fase de su espín al someterla a ciertos campos magnéticos. Estos campos, denominados campos gradientes bipolares, consisten en la aplicación de un campo gradiente (variable con la posición) durante cierto periodo de tiempo y posteriormente la aplicación del campo gradiente opuesto durante el mismo periodo de tiempo, en la dirección sobre la que se quiere medir la velocidad [2]. El cambio de fase provocado por los campos gradientes en el espín de una partícula es 𝜃=𝛾∫𝐵0 𝑡 0+𝐺(𝜏)𝑟(𝜏)𝑑𝜏 Ecuación 21. Desfase provocado por un campo gradiente [10]. donde 𝑟(𝜏) es la posición de la partícula en función del tiempo.
16 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID 𝑟(𝜏)=𝑟0+𝑣𝜏𝜏+1 2𝑎𝜏𝜏2+⋯ Ecuación 22. Posición de una partícula en función del tiempo [10]. donde 𝑟0 es la posición inicial de la partícula, 𝑣𝜏 su velocidad, 𝑎𝜏 su aceleración… Tomando en consideración solo las partículas con velocidad en la dirección sobre la que se quiere medir, y considerando despreciables los términos de orden superior a 2 de la Ecuación 22, la Ecuación 21 queda reducida a 𝜃=𝛾𝑥0∫𝐺𝑥(𝜏)𝑑𝜏+𝛾𝑣𝑥∫𝐺𝑥(𝜏)𝜏𝑑𝜏 𝑡 0 𝑡 0 Ecuación 23. Desfase provocado por un campo gradiente. Simplificación [10]. donde la primera integral se conoce como momento 0 del campo gradiente (M0, no confundir con la densidad de magnetización en equilibrio térmico 𝑀0 ) y la segunda integral como primer momento del campo gradiente (M1). 𝜃=𝛾(𝑥0𝑀0+𝑣𝑥𝑀1) Ecuación 24. Desfase provocado por un campo gradiente. Segunda simplificación [10]. Para una partícula estacionaria, la aplicación del primer campo gradiente modificará la fase de su espín en cierto grado, y la aplicación del segundo campo gradiente la modificará en la misma medida, pero en sentido opuesto, resultando una variación en la fase del espín neta nula. Para una partícula en movimiento, la aplicación del primer campo gradiente modificará la fase de su espín en cierto grado. Sin embargo, las posiciones en las que se encontrará la partícula durante la aplicación del segundo campo gradiente estarán sometidas a intensidades de campo magnético distintas, y por tanto la fase de su espín se verá modificada en distinta medida que durante la aplicación del primer campo gradiente. Esto ocasiona una variación en la fase de su espín neta proporcional a su velocidad y al primer momento del campo gradiente aplicado.
1 Introducción y objetivos 17 UNIVERSIDAD DE VALLADOLID Δ𝜃=𝛾𝑣𝑥Δ𝑀1 Ecuación 25. Desfase provocado por un campo gradiente bipolar [10]. Detectando esta variación en la fase del espín es posible conocer la velocidad de la partícula. La PC-MRA permite determinar la velocidad de una partícula solo en una dirección del espacio cada vez que se aplica. Mediante la aplicación repetida de esta técnica en varias direcciones es posible obtener un mapa de velocidades de 4 dimensiones (3 dimensiones espaciales y el tiempo) [2]. Esta técnica se conoce como 4D Flow. El 4D Flow, a pesar de ser una técnica muy potente para conocer la velocidad de un flujo, no está exento de limitaciones. Algunas de ellas son el ruido, la baja resolución espaciotemporal o la incapacidad para determinar otros parámetros hemodinámicos, como el esfuerzo cortante en las paredes de los vasos sanguíneos. Para compensar todas estas carencias, es posible combinar el 4D Flow con la Dinámica de Fluidos Computacional (CFD por sus siglas en inglés, Computational Fluid Dynamics). La CFD es capaz de conseguir una resolución espaciotemporal muy elevada, no tiene ruido y permite calcular todas las características del fluido. A su vez, para conseguir un resultado adecuado mediante CFD, es necesario que la geometría y las condiciones de contorno estén descritas de forma precisa, lo cual se puede conseguir con la ayuda del 4D Flow [3]. Así, cada técnica cubre las carencias de la otra, de modo que la combinación de ambas permite conseguir resultados precisos y de alta resolución. Además, es posible simular la adquisición de resultados por 4D Flow a partir de una simulación de trayectorias de partículas de CFD. Las partículas simuladas modelan los isocromatos del 4D Flow, y permiten resolver las ecuaciones de Bloch (Ecuación 18) para cada isocromato [7].
18 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID 1.2 Objetivos El objetivo global de este TFG es complementar las técnicas de medida de flujos por 4D Flow mediante técnicas de CFD. Para ello, se plantean los siguientes subobjetivos: 1. Determinar un modelo fluidodinámico para la simulación de flujos en arterias grandes. 2. Aplicar métodos de detección de vórtices a la localización de zonas de recirculación y estasis en simulaciones de flujo cardiovascular. 3. Contribuir mediante la simulación de trayectorias de partículas fluidas a la simulación de la obtención de flujos por 4D Flow. 1.3 Estructura del TFG Este trabajo se divide en los siguientes capítulos. Capítulo 1. Introducción y objetivos. En este capítulo se introduce el tema sobre el que va a tratar el TFG, se describen los fundamentos de las técnicas actuales para la obtención de Imágenes por Resonancia Magnética y su relación con la Dinámica de Fluidos Computacional y se plantean los objetivos del trabajo. Capítulo 2. Desarrollo del TFG. En este capítulo se describen las ecuaciones y propiedades de la mecánica de fluidos relacionadas con el TFG, se exponen las principales características, metodologías y modelos utilizados en la Dinámica de Fluidos Computacional y de ellos se escogen los modelos más adecuados para el caso a resolver. Después, se exponen y comentan los resultados obtenidos. Capítulo 3. Conclusiones. En este capítulo se evalúa el cumplimiento de los objetivos planteados, se comentan las repercusiones del trabajo y se sugieren varias líneas de trabajo futuras relacionadas con el TFG. Capítulo 4. Bibliografía. En este capítulo se enuncian las fuentes de las que se ha obtenido la información que se expone en el trabajo, formateadas según el estándar IEEE.
2 Desarrollo del TFG 19 UNIVERSIDAD DE VALLADOLID 2 Desarrollo del TFG 2.1 Modelo fluidodinámico Para realizar un uso correcto de la CFD es necesario tener conocimientos de la mecánica de fluidos subyacente para así poder escoger modelos adecuados, definir condiciones de contorno e iniciales correctamente y analizar los resultados críticamente [11]. Por ello, en este apartado se proporcionan unas nociones básicas de los aspectos de la mecánica de fluidos más relevantes para este trabajo. 2.1.1 Ecuaciones de conservación Ecuaciones diferenciales de conservación Ecuación diferencial de conservación de una propiedad genérica del fluido Antes de empezar a exponer las ecuaciones de conservación de las diferentes propiedades del fluido es conveniente saber identificar los distintos términos que aparecerán en ellas. Para ello, se presenta la ecuación de conservación de una propiedad genérica del fluido, representada como 𝜙 [11]. 𝜌𝜕𝜙 𝜕𝑡+𝜌𝜕(𝑈𝑗𝜙) 𝜕𝑥𝑗=𝜕 𝜕𝑥𝑗(Γ𝜕𝜙 𝜕𝑥𝑗)+𝑆𝜙 Ecuación 26. Ecuación diferencial de conservación de una propiedad genérica del fluido [11]. El primer término corresponde a la acumulación, y es identificable por la derivada temporal. El segundo término corresponde al transporte por convección, y es identificable por la velocidad y la derivada de primer orden. El tercer término corresponde al transporte por difusión, y es identificable por la derivada segunda y el coeficiente de transporte (difusividad, conductividad o viscosidad). El cuarto término corresponde a la generación, y es función únicamente de variables locales [11].
20 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Ecuación diferencial de conservación de la masa o de continuidad Realizando un balance de materia en un volumen infinitesimal de fluido se obtiene la ecuación de continuidad. 𝜕𝜌 𝜕𝑡+𝜕(𝜌𝑈𝑗) 𝜕𝑥𝑗=0 Ecuación 27. Ecuación diferencial de conservación de la masa. Forma 1 [11]. A veces es conveniente describir la conservación de la masa de un elemento fluido que se encuentra en movimiento con el flujo. Para ello, se desarrolla la derivada del producto en la Ecuación 27, obteniendo 𝜕𝜌 𝜕𝑡+𝑈1𝜕𝜌 𝜕𝑥1+𝑈2𝜕𝜌 𝜕𝑥2+𝑈3𝜕𝜌 𝜕𝑥3=−𝜌(𝜕𝑈1 𝜕𝑥1+𝜕𝑈2 𝜕𝑥2+𝜕𝑈3 𝜕𝑥3) Ecuación 28. Ecuación diferencial de conservación de la masa. Forma 2 [11]. Esta ecuación también se puede expresar como 𝐷𝜌 𝐷𝑡+𝜌∙𝑑𝑖𝑣(𝑈 )=0 Ecuación 29. Ecuación diferencial de conservación de la masa. Forma 3 [11]. donde el primer término es la derivada sustancial de la densidad. La derivada sustancial representa la derivada temporal de una propiedad de un elemento fluido que se mueve con el flujo 𝐷 𝐷𝑡≡𝜕 𝜕𝑡+𝑈1𝜕 𝜕𝑥1+𝑈2𝜕 𝜕𝑥2+𝑈3𝜕 𝜕𝑥3 Ecuación 30. Derivada sustancial [11].
2 Desarrollo del TFG 21 UNIVERSIDAD DE VALLADOLID Se puede asumir que el flujo es incompresible cuando se cumple 𝑀= 𝑈 √𝛾𝑅𝑇<0.1 Ecuación 31. Condición de flujo incompresible [11]. En los casos en los que el flujo se puede considerar incompresible, la densidad es constante y uniforme, y por tanto su derivada sustancial es nula. La Ecuación 29 queda reducida a 𝑑𝑖𝑣(𝑈 )=0 Ecuación 32. Ecuación diferencial de conservación de la masa para flujos incompresibles [11]. Ecuaciones diferenciales de conservación de la cantidad de movimiento o de Navier-Stokes Realizando un balance de cantidad de movimiento en un volumen infinitesimal de fluido se obtiene 𝜕(𝜌𝑈𝑖) 𝜕𝑡 +𝜕(𝜌𝑈𝑖𝑈𝑗) 𝜕𝑥𝑗=𝜕𝜏𝑖𝑗 𝜕𝑥𝑗+𝜌𝑓𝑚𝑖 Ecuación 33. Ecuaciones de conservación de la cantidad de movimiento. Forma 1 [11]. Desarrollando las derivadas de los productos 𝜌𝐷𝑈𝑖 𝐷𝑡+𝑈𝑖[𝐷𝜌 𝐷𝑡+𝜌∙𝑑𝑖𝑣(𝑈 )] =𝜕𝜏𝑖𝑗 𝜕𝑥𝑗+𝜌𝑓𝑚𝑖 Ecuación 34. Ecuaciones de conservación de la cantidad de movimiento. Forma 2 [11]. donde el segundo término es nulo por cumplir la ecuación de continuidad, resultando
28 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID El coste de la CFD se limita al del equipo en el que el software se va a ejecutar, el del propio software y el salario del empleado que realiza la simulación. Este coste es reducido en comparación con el de otros métodos de resolución, como los métodos empíricos, que involucran elevados costes materiales. Existe una amplia variedad de softwares comerciales, que pueden ser usados incluso por usuarios inexpertos para obtener soluciones a problemas sencillos, aunque para problemas más avanzados el usuario debe tener conocimientos suficientes en mecánica de fluidos para determinar qué modelo es el más adecuado para el caso de estudio y para detectar errores o falta de precisión en el resultado. La CFD, a pesar de no proporcionar las soluciones exactas de los problemas, es capaz de proporcionar soluciones suficientemente precisas en una amplia gama de situaciones. Incluso en los casos en los que la CFD no es capaz de proporcionar soluciones con suficiente precisión, como en sistemas multifásicos o en sistemas con velocidades de reacción muy elevadas, la CFD puede resultar de utilidad, ya que proporciona cierta información que puede ayudar a definir las pruebas experimentales que se deben realizar, así como sus condiciones, para obtener soluciones más precisas [11]. 2.2.2 Métodos numéricos frente a otros métodos La CFD se basa en la resolución de las ecuaciones de la mecánica de fluidos mediante métodos numéricos. Los métodos numéricos presentan ciertas ventajas sobre otros métodos de resolución. Los métodos analíticos proporcionan una solución exacta del problema en los casos en los que se puede asumir una simplificación (esto es, cuando la viscosidad es despreciable, el movimiento del fluido es unidireccional, etc.). Sin embargo, su mayor inconveniente es que aún no se ha encontrado una solución analítica general para las ecuaciones que gobiernan el movimiento de los fluidos, y por lo tanto su uso queda limitado a unos pocos casos particulares. Los métodos empíricos o semi-empíricos funcionan muy bien para sistemas ampliamente estudiados. Sin embargo, al modificar las condiciones de los procesos estudiados o al desarrollar equipos y procesos nuevos, generalmente no se dispone de una base de datos lo suficientemente amplia como para que un modelo empírico o semi-empírico resulte fiable.
2 Desarrollo del TFG 29 UNIVERSIDAD DE VALLADOLID Los métodos numéricos, además de no contar con ninguno de los inconvenientes que presentan los métodos anteriores, cuentan con la gran ventaja de que su implementación en ordenador se puede llevar a cabo con facilidad [11]. 2.2.3 Pasos para resolver un problema por CFD Los pasos a seguir para la resolución de un problema de mecánica de fluidos mediante CFD son: 1. Modelado de la geometría. Definir la geometría del sistema que se va a estudiar. Para ello, los softwares comerciales cuentan con programas de Diseño Asistido por Ordenador (CAD por sus siglas en inglés, Computer-Aided Design) integrados, aunque también se puede definir la geometría en un programa CAD externo y posteriormente importarla. 2. Mallado. El mallado consiste en la división del volumen de estudio en elementos más pequeños denominados celdas. Cada celda funciona como un elemento discreto y toma un único valor de las propiedades del fluido. 3. Definición de modelos. Mientras que para casos sencillos se puede obtener una solución numérica de las ecuaciones de Navier-Stokes, para flujos más complejos es necesario utilizar otros modelos, empíricos o semi-empíricos, que pueden estar definidos en el programa CFD o ser definidos por el usuario mediante una subrutina. 4. Definición de las propiedades del fluido. Se definen las propiedades físicas del fluido tales como viscosidad, densidad o composición. Los softwares comerciales disponen de polinomios con coeficientes modificables por el usuario para los casos en los que estas propiedades no son constantes. También se pueden definir mediante una subrutina escrita por el usuario. 5. Definición de condiciones iniciales y de contorno. Se especifican las propiedades del fluido en las entradas, salidas, paredes y otras fronteras, así como en el instante inicial.
30 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID 6. Resolución. Se definen las propiedades del solver y se ejecuta la simulación. El simulador resuelve las ecuaciones de la mecánica de fluidos numéricamente en cada celda y para cada instante de tiempo. 7. Análisis de la solución. Se analiza si la solución es razonable y precisa, su dependencia del tamaño de malla y de las condiciones iniciales y de contorno, su relación con los resultados experimentales, etc. [11]. 2.2.4 Aspectos numéricos de la CFD El método de volúmenes finitos Las primeras simulaciones por CFD empleaban el método de diferencias finitas para obtener una aproximación algebraica de las ecuaciones de la mecánica de fluidos. Este método consiste en aproximar las derivadas parciales por su desarrollo en serie de Taylor. El objetivo del método de diferencias finitas es conseguir el mínimo error global entre la solución exacta y la solución aproximada. Sin embargo, en mecánica de fluidos es importante garantizar la conservación local, y el método de diferencias finitas no la garantiza necesariamente, y por tanto no es la mejor elección. El método de volúmenes finitos, por el contrario, tiene por uno de sus principios garantizar la conservación local. Este método consiste en dividir el volumen de estudio en pequeños sub-volúmenes llamados celdas. Después, se aplica el teorema de Gauss-Ostrogradski en cada celda para transformar las ecuaciones diferenciales en ecuaciones integrales, se establecen las superficies limitantes de la celda como límites de integración y se integra, obteniendo un conjunto de ecuaciones algebraicas que se resuelven mediante métodos numéricos. Debido a la discretización del volumen de estudio existe un error numérico. Este error se aproxima a cero al tender a cero el tamaño de la celda. Por tanto, para asegurar un resultado preciso es suficiente con escoger un tamaño de celda lo suficientemente pequeño. Sin embargo, escoger un tamaño de celda demasiado pequeño conlleva aumentar el número de celdas y por lo tanto el tiempo de computación, por lo que se debe encontrar un equilibrio. La forma de las celdas es un factor a tener en cuenta para la aplicación del método de volúmenes finitos. Lo ideal es disponer de una malla que contenga únicamente celdas con ángulos internos de 90 grados. Esto se denomina
2 Desarrollo del TFG 31 UNIVERSIDAD DE VALLADOLID malla estructurada, y es el caso que se va a desarrollar en los siguientes apartados [11]. Figura 10. Celda estructurada [11]. Sin embargo, en la mayoría de los casos prácticos, las geometrías tienen partes que no se pueden dividir mediante este tipo de celdas, y por tanto esta no es una estructura muy común en la práctica [11]. Balances de celda Cada uno de los términos de la Ecuación 44 se debe transformar en una expresión algebraica para poder resolver la ecuación mediante un método numérico. - Por el momento solo se considerarán problemas estacionarios, y por lo tanto el término de acumulación es nulo. ∫ 𝜌𝜕𝜙 𝜕𝑡𝑑𝑉 𝑉𝐶 =0 Ecuación 52. Integración del término de acumulación en una celda estructurada [11]. - El término de convección se integra con las caras de una celda estructurada como límites de integración.
32 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID ∫ 𝜌𝑈𝑗𝑛𝜙𝑑𝐴 𝑆𝐶 =−𝜌[(𝐴𝑈𝜙)𝑤−(𝐴𝑈𝜙)𝑒+ +(𝐴𝑈𝜙)𝑠−(𝐴𝑈𝜙)𝑛+(𝐴𝑈𝜙)𝑡−(𝐴𝑈𝜙)𝑏] Ecuación 53. Integración del término de convección en una celda estructurada [11]. - El término de convección se integra con las caras de una celda estructurada como límites de integración. ∫ Γ𝜕𝜙 𝜕𝑥𝑗𝑛𝑑𝐴 𝑆𝐶 = −[(𝐴Γ𝜕𝜙 𝜕𝑥)𝑤−(𝐴Γ𝜕𝜙 𝜕𝑥)𝑒+ +(𝐴Γ𝜕𝜙 𝜕𝑥)𝑠−(𝐴Γ𝜕𝜙 𝜕𝑥)𝑛+(𝐴Γ𝜕𝜙 𝜕𝑥)𝑡−(𝐴Γ𝜕𝜙 𝜕𝑥)𝑏] Ecuación 54. Integración del término de difusión en una celda estructurada [11]. - El término de generación se promedia en cada celda. ∫ 𝑆𝜙𝑑𝑉 𝑉𝐶 ≈𝑆𝜙𝑉 Ecuación 55. Integración del término de generación en una celda estructurada [11]. Esquemas de discretización espacial Para resolver las ecuaciones algebraicas que resultan de las trasformaciones anteriores es necesario hacer una estimación de los valores de 𝜙 y el gradiente de 𝜙 en las caras de las celdas, Para hacerlo existen varios métodos, llamados esquemas de discretización [11]. Delimitación y transportividad Existen dos propiedades que resultan deseables en un esquema de discretización: delimitación y transportividad.
2 Desarrollo del TFG 33 UNIVERSIDAD DE VALLADOLID Un esquema de discretización está delimitado si las propiedades estimadas están comprendidas entre los valores que se han utilizado para estimarlas. Un esquema de discretización posee transportividad si la dirección de transmisión de la información (la dirección del flujo) resulta obvia con tan solo conocer el esquema [11]. A continuación se describen algunos de los esquemas de discretización más relevantes. Se asume que el flujo es unidimensional y su sentido es de oeste a este. Los subíndices W, P, E, w, e representan la celda al oeste, la celda actual, la celda al este, la cara oeste de la celda actual y la cara este de la celda actual respectivamente. Esquema de diferencias centradas El valor de 𝜙 en una cara se calcula mediante una interpolación lineal entre los valores de 𝜙 en las celdas que comparten la cara. 𝜙𝑤=𝜙𝑊+𝜙𝑃 2 𝜙𝑒=𝜙𝐸+𝜙𝑃 2 Ecuación 56. Esquema de diferencias centradas [11]. El valor del gradiente se calcula mediante la aproximación de primer orden de su desarrollo en serie de Taylor. Esto es así para todos los esquemas de discretización que se describen a continuación, y por tanto no se repetirá en cada uno de ellos. (𝑑𝜙 𝑑𝑥)𝑤=𝜙𝑃−𝜙𝑊 𝑥𝑃−𝑥𝑊 (𝑑𝜙 𝑑𝑥)𝑒=𝜙𝐸−𝜙𝑃 𝑥𝐸−𝑥𝑃 Ecuación 57. Estimación del gradiente [11]. Este esquema está condicionalmente delimitado y no tiene transportividad. El esquema falla en problemas en los que los fenómenos convectivos son de importancia. El error que se comete es del orden de la segunda potencia de la distancia entre celdas [11].
34 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Esquema upwind de primer orden El valor de 𝜙 en una cara es igual al valor de 𝜙 en la celda inmediatamente aguas arriba de la cara. 𝜙𝑤=𝜙𝑊 𝜙𝑒=𝜙𝑃 Ecuación 58. Esquema upwind de primer orden [11]. Este esquema está delimitado y tiene transportividad. El error que se comete es del orden de la distancia entre celdas [11]. Esquema upwind de segundo orden El valor de 𝜙 en una cara se calcula asumiendo que el gradiente entre la cara y la celda inmediatamente aguas arriba de la cara, y el gradiente entre la celda inmediatamente aguas arriba de la cara y dos celdas aguas arriba de la cara, son iguales. 𝜙𝑒−𝜙𝑃 𝑥𝑒−𝑥𝑃=𝜙𝑃−𝜙𝑊 𝑥𝑃−𝑥𝑊→𝜙𝑒=𝜙𝑃+(𝜙𝑃−𝜙𝑊)(𝑥𝑒−𝑥𝑃) 𝑥𝑃−𝑥𝑊 Ecuación 59. Esquema upwind de segundo orden [11]. Este esquema no está delimitado y tiene transportividad. Existen otros esquemas de segundo orden que solucionan el problema de delimitación, como el esquema van Leer [11]. Esquema van Leer Si el gradiente es pequeño, se utiliza una expresión similar a la del esquema upwind de segundo orden. Si el gradiente es grande, se utiliza la expresión del esquema upwind de primer orden. Así, se consigue resolver el problema de delimitación del esquema upwind de segundo orden.
2 Desarrollo del TFG 35 UNIVERSIDAD DE VALLADOLID 𝜙𝑒={𝜙𝑃+(𝜙𝑃−𝜙𝑊)(𝜙𝐸−𝜙𝑃) 𝜙𝐸−𝜙𝑊 𝜙𝑃 𝑠𝑖 |𝜙𝐸−2𝜙𝑃+𝜙𝑊|≤|𝜙𝐸−𝜙𝑊| 𝑠𝑖 |𝜙𝐸−2𝜙𝑃+𝜙𝑊|>|𝜙𝐸−𝜙𝑊| Ecuación 60. Esquema van Leer [11]. Resolución del campo de velocidades Para resolver las ecuaciones de conservación como se ha explicado en el apartado anterior es necesario conocer previamente el campo de velocidades. Para ello se debe resolver un sistema de ecuaciones que comprende las tres ecuaciones de Navier-Stokes y la ecuación de continuidad, siendo las incógnitas las tres componentes de la velocidad y la presión. Las tres ecuaciones de Navier-Stokes se utilizan para hallar las tres componentes de la velocidad. Sin embargo, no es posible hallar la presión directamente a partir de la ecuación de continuidad, por lo que se deben hacer ciertas modificaciones. Se toma la divergencia de las ecuaciones de Navier-Stokes y, tras reescribir algunos términos, asumiendo densidad y viscosidad constantes, se combina con la ecuación de continuidad [11]. 𝜕2𝑃 𝜕𝑥𝑖𝜕𝑥𝑖=−𝜕2(𝜌𝑈𝑖𝑈𝑗) 𝜕𝑥𝑖𝜕𝑥𝑗 Ecuación 61. Ecuación para la presión [11]. Esta ecuación sí puede ser utilizada para hallar la presión directamente. Sin embargo, su resolución mediante métodos numéricos presenta algunas dificultades. Por ello, para resolverla se utilizan algoritmos específicos. Algunos de los más usados son SIMPLE, SIMPLEC, SIMPLER o PISO [11]. Sub-relajación Las ecuaciones de Navier-Stokes contienen términos no lineales que pueden provocar que el método numérico diverja. Para evitarlo se utiliza un método denominado sub-relajación.
36 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID En cada iteración, se tomará como nuevo valor de 𝜙, en vez del valor obtenido directamente por el método numérico, una media ponderada entre el valor obtenido por el método numérico y el valor anterior de 𝜙 [11]. 𝜙𝑛𝑢𝑒𝑣𝑜=𝛼 𝜙𝑚é𝑡𝑜𝑑𝑜+(1−𝛼) 𝜙𝑎𝑛𝑡𝑒𝑟𝑖𝑜𝑟 Ecuación 62. Sub-relajación [11]. 𝛼 se denomina factor de sub-relajación y está comprendido entre 0 y 1. Su valor puede variar entre iteraciones y cada ecuación puede tener su propio factor de sub-relajación. Un valor de 𝛼 demasiado bajo provoca tiempos de computación elevados y un valor demasiado alto puede provocar que el método diverja, por lo que es necesario encontrar un equilibrio [11]. Flujos no estacionarios Antes de comenzar con la explicación sobre la resolución de problemas con flujo no estacionario, es conveniente hacer un comentario sobre su definición. Si se define un flujo estacionario como un flujo que varía en el tiempo, casi todos los flujos serían no estacionarios, ya que la mayoría de los flujos de importancia industrial son turbulentos y la turbulencia implica variación en el tiempo de las propiedades del fluido. Por tanto, es conveniente ajustar esta definición. En el estudio de flujos turbulentos se toman en consideración las propiedades del fluido promediadas en el tiempo en vez de las cantidades instantáneas. Por tanto, en estos casos es conveniente definir la no estacionariedad de un flujo como la variación en el tiempo de dichas propiedades promediadas. El método de resolución de problemas con flujo no estacionario es similar al de problemas con flujo estacionario. Primero, se integra la Ecuación 26, en este caso no solo en el espacio, sino también en el tiempo [11]. ∫ (∫ 𝜌𝜕𝜙 𝜕𝑡𝑑𝑡 𝑡+Δ𝑡 𝑡)𝑑𝑉+∫ [∫ 𝜌𝜕(𝑈𝑗𝜙) 𝜕𝑥𝑗 𝑉𝐶 𝑑𝑉]𝑑𝑡= 𝑡+Δ𝑡 𝑡𝑉𝐶
2 Desarrollo del TFG 37 UNIVERSIDAD DE VALLADOLID =∫ [∫ 𝜕 𝜕𝑥𝑗(Γ𝜕𝜙 𝜕𝑥𝑗) 𝑉𝐶 𝑑𝑉]𝑑𝑡+∫ (∫ 𝑆𝜙𝑑𝑉 𝑉𝐶 )𝑑𝑡 𝑡+Δ𝑡 𝑡 𝑡+Δ𝑡 𝑡 Ecuación 63. Integración espacial y temporal de la ecuación de conservación de una propiedad genérica del fluido [11]. Después, se aplican los esquemas de discretización. Para los términos espaciales se utilizan los esquemas de discretización espacial expuestos anteriormente. Para los términos temporales se utilizan los esquemas de discretización temporal que se expondrán a continuación. Una vez se ha discretizado la ecuación, se realiza un procedimiento iterativo tanto en el tiempo como en el espacio: para cada iteración en el tiempo se realizan sub-iteraciones en el espacio hasta alcanzar la convergencia en ese instante, y se continúa iterando en el tiempo hasta alcanzar el instante de tiempo deseado. Para cada instante se utiliza como aproximación inicial la solución alcanzada en el instante anterior [11]. Esquemas de discretización temporal Para integrar los términos dependientes de 𝜙 en cada paso temporal, es necesario estimar un valor de 𝜙 para el paso completo, ya que su valor solo se conoce en los instantes inicial y final [11]. Para ello, existen varias opciones: tomar el valor de 𝜙 al principio del paso, al final del paso o una mezcla de los dos. Esto se puede expresar utilizando un factor de peso 𝜃 de forma que ∫ 𝜙𝑑𝑡 𝑡+Δ𝑡 𝑡=[θ𝜙(𝑡+Δ𝑡)+(1−θ)𝜙(𝑡)]Δ𝑡 Ecuación 64. Factor e peso en la integración temporal [11]. Método explícito El valor de 𝜃 es 0, es decir, se considera que el valor de 𝜙 durante todo el intervalo de integración es igual a su valor al inicio del intervalo [11]. El error que se comete es del orden del paso de tiempo. Para que este esquema esté delimitado se debe cumplir la condición
44 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Ω=1 2(𝜕𝑢𝑖 𝜕𝑥𝑗−𝜕𝑢𝑗 𝜕𝑥𝑖)=1 2 [ 0𝜕𝑢1 𝜕𝑥2−𝜕𝑢2 𝜕𝑥1𝜕𝑢1 𝜕𝑥3−𝜕𝑢3 𝜕𝑥1 𝜕𝑢2 𝜕𝑥1−𝜕𝑢1 𝜕𝑥20𝜕𝑢2 𝜕𝑥3−𝜕𝑢3 𝜕𝑥2 𝜕𝑢3 𝜕𝑥1−𝜕𝑢1 𝜕𝑥3𝜕𝑢3 𝜕𝑥2−𝜕𝑢2 𝜕𝑥30 ] Ecuación 76. Tensor de velocidades de rotación [23]. Algunos de los métodos más importantes se explican a continuación. Método Q Se denota como Q el segundo invariante del tensor gradiente de velocidad. 𝑄=1 2(‖Ω‖2−‖𝑆‖2) Ecuación 77. Invariante Q [24]. Se define un vórtice como una región conectada para la cual Q es positivo en todos sus puntos, o equivalentemente, una región conectada en la cual la velocidad local de rotación del fluido es mayor que su velocidad local de deformación [24]. Método Lambda2 Si se desprecian los términos no estacionarios y los debidos a la viscosidad de las ecuaciones de Navier-Stokes se cumple que 𝑆2+Ω2=−1 𝜌∇2𝑝 Ecuación 78. Tensor S2 + Ω2 [24].
2 Desarrollo del TFG 45 UNIVERSIDAD DE VALLADOLID 𝑆2+Ω2 es una matriz real y simétrica y por tanto tiene tres valores propios reales. Se denotan dichos valores propios, en orden decreciente, como 𝜆1≥ 𝜆2≥𝜆3. Se define un vórtice como una región conectada para la cual 𝜆2 es negativo en todos sus puntos [24]. 2.3 Caso de estudio La geometría de estudio es un phantom de la arteria aorta diseñado y modelado mediante CAD por el Departamento de Ingeniería Energética y Fluidomecánica. Figura 11. Geometría de estudio tridimensional. El phantom consta, aparte de la geometría que emula la de la arteria aorta y sus intersecciones con otras arterias, de dos depósitos y dos tubos (uno a la entrada y otro a la salida) para la alimentación del fluido. Previamente al estudio de la geometría presentada, se realizará el estudio de una geometría bidimensional obtenida a partir del plano medio de la geometría tridimensional.
46 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Figura 12. Geometría de estudio bidimensional. Así, en caso de cometer algún error al establecer las condiciones de contorno, los parámetros del solver o algún otro aspecto de la simulación, el coste computacional del error, y consecuentemente su coste temporal, serán mucho más pequeños que si se hubiese cometido el mismo error directamente en la simulación tridimensional. El estudio de ambas geometrías se realizará utilizando el software Ansys, concretamente sus módulos SpaceClaim, DesignModeler, Fluent y CFD-Post. 2.3.1 Caso bidimensional Mallado Se comienza realizando el mallado de la geometría bidimensional. Para ello, se establecen cinco áreas: - Área A. Es el área general. Se malla con elementos de 1 mm de lado. - Áreas B, C y D. Son las áreas donde la representación de la aorta interseca con las representaciones de otras arterias del cuerpo humano (el tronco braquiocefálico, la arteria carótida común izquierda y la arteria subclavia izquierda, respectivamente). En estas
2 Desarrollo del TFG 47 UNIVERSIDAD DE VALLADOLID áreas se colocan esferas de influencia de 20 mm de radio en las que el mallado se realiza con elementos de 0,5 mm de lado. - Área E. Es el contorno de la geometría. Se malla mediante dos capas de elementos de 0,1 mm de espesor. Figura 13. División en áreas para el mallado de la geometría bidimensional. Además, se da nombre a la entrada, las salidas y el contorno de la geometría. Figura 14. Nomenclatura de las secciones del contorno de la geometría bidimensional.
48 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Definición de modelos Para escoger el modelo de viscosidad es necesario conocer el régimen de flujo. Para ello, se calcula el número de Reynolds. Antes de ejecutar la simulación, únicamente se conoce la velocidad en la entrada. Se calcula el número de Reynolds tomando su valor máximo. 𝑅𝑒𝑝𝑖𝑐𝑜=𝑢𝑝𝑖𝑐𝑜∙𝐷 𝜈=0,5∙0,02 9,264∙10−6=1079,45 Ecuación 79. Número de Reynolds [25]. La condición de régimen laminar para flujo interno en tuberías de sección circular es: 𝑅𝑒≤2300 Ecuación 80. Condición de régimen laminar para flujo interno en tuberías de sección circular [25]. Por tanto, el flujo es aparentemente laminar. Si tras la ejecución de la simulación existiesen velocidades sustancialmente superiores a la velocidad máxima en la entrada, esto se reconsiderará. Por el momento, se escoge el modelo de viscosidad laminar. Para la simulación de trayectorias de partículas, se escoge el modelo de fase discreta sin masa. Definición de las propiedades del fluido Las propiedades del fluido de estudio (sangre) son: - Densidad. 𝜌=1153𝑘𝑔 𝑚3 - Viscosidad cinemática. 𝜈=9,246∙10−6𝑚2 𝑠 Ansys Fluent requiere la densidad y la viscosidad dinámica para la definición de un fluido. Por tanto, se calcula la viscosidad dinámica como
2 Desarrollo del TFG 49 UNIVERSIDAD DE VALLADOLID 𝜇=𝜌∙𝜈=1,066∙10−2 𝑘𝑔 𝑚∙𝑠 Ecuación 81. Cálculo de la viscosidad dinámica. Definición de las condiciones de contorno Entrada Como condición de contorno en la entrada se establece un perfil de velocidad uniforme, variable en el tiempo y periódico. Cada periodo dura un segundo (emulando el latido de un corazón a 60 pulsos por minuto) y consta de dos partes: - En los primeros 0,436 segundos la velocidad varía con forma de onda sinusoidal, partiendo desde 0,1 m/s, alcanzando su pico en 0,5 m/s y regresando a 0,1 m/s. Esta parte emula la fase sistólica del latido del corazón. - En el resto del periodo la velocidad permanece constante e igual a 0,1 m/s. Esta parte emula la fase diastólica del latido del corazón. 𝑢𝑒𝑛𝑡𝑟𝑎𝑑𝑎(𝑡)={0,5∙𝑠𝑒𝑛(2𝜋(𝑡+0,0320472)) 0,2 𝑠𝑖 𝑡≤0,436 𝑠𝑖 𝑡>0,436 Ecuación 82. Velocidad en la entrada. donde t es el tiempo local dentro de cada periodo.
50 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Figura 15. Velocidad en la entrada. Para poder introducir la condición de contorno en la entrada en Fluent, se implementa la función en lenguaje C, tomando como base una función similar [26]. #include "udf.h"//archivo que contiene las definiciones de las funciones DEFINE y las operaciones de Fluent #define PI 3.141592654 DEFINE_PROFILE(inlet_velocity_60_BPM,th,i) { face_t f; begin_f_loop(f,th) double t = CURRENT_TIME-floor(CURRENT_TIME); //t es el tiempo local en cada periodo { if(t <= 0.436) F_PROFILE(f,th,i) = 0.5*sin(2*PI*(t+0.0320472)); else
2 Desarrollo del TFG 51 UNIVERSIDAD DE VALLADOLID F_PROFILE(f,th,i) = 0.1; } end_f_loop(f,th); } Salidas Como condición de contorno en las salidas se establece una presión manométrica de 13332 Pa. Pared arterial Como condición de contorno en la pared arterial se establece la condición de no deslizamiento (la velocidad de las partículas fluidas próximas a la pared es igual a la velocidad de la pared, en este caso nula). Definición de las propiedades del solver Se configura el solver de la siguiente forma: - Para la resolución del campo de velocidades se utiliza el esquema SIMPLE. - Para la discretización espacial: o Para la discretización del gradiente se utiliza el método de mínimos cuadrados. o Para la discretización de la presión se utiliza el método de segundo orden. o Para la discretización del momento se utiliza el esquema upwind de segundo orden. - Para la discretización temporal se utiliza el método implícito de primer orden. - Los factores de sub-relajación son: o 0,3 para la presión. o 1 para la densidad.
52 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID o 1 para las fuerzas másicas. o 0,7 para el momento. - Los parámetros de avance del tiempo son: o Se realizan 400 pasos temporales. o Cada paso es de 0,05 segundos. o Se realiza un máximo de 50 iteraciones por paso temporal. Para comprobar si el paso temporal es adecuado, se calcula el número de Courant. Para el área general de la malla 𝐶=u∙Δ𝑡 Δ𝑥 =0,5∙0,05 0,0032 =7,81 Ecuación 83. Cálculo del número de Courant para el área general de la malla. Para el área afectada por las esferas de influencia de la malla 𝐶=u∙Δ𝑡 Δ𝑥 =0,5∙0,05 0,002 =12,5 Ecuación 84. Cálculo del número de Courant para el área afectada por las esferas de influencia de la malla. Dado que la discretización temporal se realiza mediante el método implícito, estos valores son aceptables y el paso temporal es adecuado. 2.3.2 Caso tridimensional Mallado Para el mallado de la geometría tridimensional se establecen las mismas cinco zonas (en este caso se denominarán volúmenes desde A hasta E, en vez de áreas) que para el de la geometría bidimensional.
2 Desarrollo del TFG 53 UNIVERSIDAD DE VALLADOLID Figura 16. División en áreas para el mallado de la geometría tridimensional. Sin embargo, los tamaños de los elementos en cada zona difieren de los tamaños descritos en el caso bidimensional, ya que si fuesen iguales resultaría un número de celdas mucho mayor, lo cual conllevaría dos problemas: primero, el coste computacional de la simulación sería inasequiblemente alto; segundo, el programa (Ansys Fluent, con licencia para estudiantes) impone un número máximo de celdas (512.000) que se superaría. Para ajustarse a esta restricción, se tantea con los tamaños de elemento hasta llegar a una malla que posee un número de celdas cercano al máximo. Los tamaños de elemento a los que se llega son: - Volumen A. 3,2 mm de lado. - Volúmenes B, C y D. 2 mm de lado. - Volumen E. 5 capas de 0,1 mm cada una. Se obtiene así una malla de 501.719 elementos (malla 1). Dado que la malla se ha obtenido atendiendo más a las limitaciones impuestas por el programa que a que el tamaño de sus elementes sea el adecuado, es necesario realizar un estudio de convergencia de malla para comprobar su viabilidad. Para ello, es necesario disponer de tres mallas. La malla 1 será la malla de alta densidad. Se define la malla de densidad intermedia con los siguientes tamaños de elemento:
60 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Figura 22. Geometría bidimensional. Campo de velocidades. Vista general. t=19,3 s. Figura 23. Geometría bidimensional. Campo de velocidades. Vista general. t=19,4 s.
2 Desarrollo del TFG 61 UNIVERSIDAD DE VALLADOLID Figura 24. Geometría bidimensional. Campo de velocidades. Vista general. t=19,5 s. Figura 25. Geometría bidimensional. Campo de velocidades. Vista general. t=19,6 s.
62 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Figura 26. Geometría bidimensional. Campo de velocidades. Vista general. t=19,7 s. Figura 27. Geometría bidimensional. Campo de velocidades. Vista general. t=19,8 s.
2 Desarrollo del TFG 63 UNIVERSIDAD DE VALLADOLID Figura 28. Geometría bidimensional. Campo de velocidades. Vista general. t=19,9 s. Figura 29. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19 s.
64 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Figura 30. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,1 s. Figura 31. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,2 s.
2 Desarrollo del TFG 65 UNIVERSIDAD DE VALLADOLID Figura 32. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,3 s. Figura 33. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,4 s.
66 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Figura 34. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,5 s. Figura 35. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,6 s.
2 Desarrollo del TFG 67 UNIVERSIDAD DE VALLADOLID Figura 36. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,7 s. Figura 37. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,8 s.
68 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Figura 38. Geometría bidimensional. Campo de velocidades. Vista detalle. t=19,9 s. Al comienzo del ciclo, el flujo es lento. En las zonas más periféricas de los depósitos de entrada y salida se observa variedad en las direcciones de los vectores de velocidad, lo que indica que el fluido se encuentra en rotación (Figura 19). También es observable un ligero desorden en algunas zonas de las intersecciones, aunque en una medida mucho menor (Figura 29). En el contorno interior de la representación de la aorta descendente y en ambos contornos de la representación de la aorta ascendente puede apreciarse flujo reverso. Al empezar a aumentar la velocidad de la entrada, todo el flujo se acelera. De las zonas que antes presentaban rotación solo lo siguen haciendo las más periféricas, mientras que en las más próximas al centro el fluido adquiere una dirección predominante de avance (Figuras 20 a 22). El flujo reverso en la representación de la aorta ascendente desaparece, aunque se mantiene en el contorno interior de la representación de la aorta descendente. Se aprecia más claramente que en instantes anteriores cómo el fluido que circula por la representación de la aorta ascendente se divide en dos partes: las partículas fluidas que se encuentran más cerca del contorno exterior se dirigen hacia las intersecciones, mientras que las partículas fluidas que se encuentran más cerca del contorno interior continúan por la representación de la aorta descendente (Figuras 30 a 32). Al comenzar a descender la velocidad de la entrada, vuelven a aparecer vectores de velocidad orientados en diversas direcciones en las zonas periféricas de los depósitos y en las intersecciones, indicativo de que el fluido
2 Desarrollo del TFG 69 UNIVERSIDAD DE VALLADOLID en esas zonas comienza a rotar. También vuelve a aparecer flujo reverso en la representación de la aorta ascendente. El fluido en la zona central de la geometría decelera, pero se mantiene ordenado (Figuras 23 a 28 y 33 a 38). La velocidad máxima en la geometría de estudio durante el ciclo es de 0,58 m/s. Aunque esta velocidad es mayor que la velocidad pico en la entrada (0,5 m/s), que es la que se ha utilizado para calcular el número de Reynolds de forma provisional (1079,45), se aprecia sin necesidad de realizar cálculos que el nuevo Reynolds no supera la condición de régimen laminar para flujo interno en tuberías de sección circular (2300). Por tanto, se confirma que la elección del modelo de viscosidad laminar es adecuada para el caso bidimensional. Para complementar la información que se ha obtenido a partir del campo de velocidades, se presentan a continuación imágenes del contorno de velocidades, tomadas durante el último ciclo simulado, cada 0,1 segundos. Estas, aunque proporcionan menos información que las del campo de velocidades, ya que indican únicamente el módulo de la velocidad y no su dirección y sentido, cuentan con la ventaja de que su sencillez permite una interpretación visual mucho más rápida e intuitiva. Figura 39. Geometría bidimensional. Contorno de velocidades. t=19 s.
76 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Figura 51. Geometría bidimensional. Vórtices. t=19,2 s. Figura 52. Geometría bidimensional. Vórtices. t=19,3 s.
2 Desarrollo del TFG 77 UNIVERSIDAD DE VALLADOLID Figura 53. Geometría bidimensional. Vórtices. t=19,4 s. Figura 54. Geometría bidimensional. Vórtices. t=19,5 s.
78 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Figura 55. Geometría bidimensional. Vórtices. t=19,6 s. Figura 56. Geometría bidimensional. Vórtices. t=19,7 s.
2 Desarrollo del TFG 79 UNIVERSIDAD DE VALLADOLID Figura 57. Geometría bidimensional. Vórtices. t=19,8 s. Figura 58. Geometría bidimensional. Vórtices. t=19,9 s. Hay que tener en cuenta que los criterios actuales de detección de vórtices, incluido el criterio Q, no se caracterizan por tener una tasa de acierto especialmente elevada, y tampoco distinguen la intensidad de los vórtices, sino únicamente su presencia o ausencia. Por ello, hay que corroborar la información que proporcionan las imágenes de vórtices con la que proporcionan las imágenes del campo de velocidades. Al comienzo del ciclo se aprecian grandes vórtices en el conducto de entrada, el depósito de entrada, la representación del arco aórtico, parte del conducto
80 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID de salida y el depósito de salida. También se aprecian pequeños vórtices en las intersecciones. Observando las imágenes del campo de velocidades se aprecia que, de estos, los únicos vórtices relevantes son los del depósito de entrada, el depósito de salida, las zonas más próximas al contorno de la representación del arco aórtico y las intersecciones (Figuras 49 y 19). Al aumentar la velocidad en la entrada, los vórtices del depósito de entrada desaparecen (el fluido deja de rotar y adquiere una velocidad predominante de avance). Aparecen dos grandes vórtices en la representación de la aorta ascendente que, observando las imágenes del campo de velocidades, se desestiman. Los vórtices del contorno externo de la representación de la aorta ascendente desaparecen también, y el tamaño los pequeños vórtices de las intersecciones disminuye hasta prácticamente desaparecer. Los únicos vórtices que quedan presentes son los del depósito de salida, que también ven reducido su tamaño (Figuras 40 a 52 y 20 a 22). Cuando la velocidad en la entrada comienza a descender, poco a poco se vuelven a formar vórtices en las zonas en las que habían desaparecido, y los vórtices del depósito de salida recuperan su tamaño inicial (Figuras 53 a 58 y 23 a 28). Aproximación lagrangiana La aproximación lagrangiana a un campo fluido es la descripción del campo mediante la especificación de las propiedades de partículas del fluido que se desplazan a través del campo [28]. En este apartado se describirán las trayectorias de dichas partículas. Se realiza una inyección de partículas sin masa en la entrada de la geometría en el instante t = 10 s. A continuación se presenta una imagen de las trayectorias de las partículas, 10 segundos después de su inyección, coloreadas según su tiempo de residencia en el volumen de estudio.
2 Desarrollo del TFG 81 UNIVERSIDAD DE VALLADOLID Figura 59. Geometría bidimensional. Trayectorias de las partículas. Las partículas parten de la entrada de forma ordenada hasta llegar al depósito de entrada. Allí, incluso aquellas que se posicionan más al centro de la geometría describen una trayectoria rizada hasta abandonar el depósito, en concordancia con el análisis de vórtices realizado anteriormente. Las partículas siguen su camino por la representación de la aorta ascendente de nuevo de forma ordenada. Hasta aquí apenas se aprecian diferencias en la velocidad de las partículas en función de su posición radial: todas las partículas, salvo aquellas situadas muy en los extremos, presentan aproximadamente los mismos tiempos de residencia al alcanzar los mismos planos transversales. Al llegar a las intersecciones, entre 4 y 5 segundos tras su inserción, la mayoría de las partículas abandonan la geometría de estudio. Se puede observar que las partículas que abandonan la geometría de estudio a través de la representación de la arteria carótida común izquierda y de la arteria subclavia izquierda (la segunda y la tercera intersección) describen una trayectoria rizada antes de hacerlo. Las partículas restantes continúan por la representación de la aorta descendente. Ahora sí se aprecian diferencias en la velocidad de las partículas en función de su posición radial: las partículas posicionadas en el centro tienen una velocidad mayor, y abandonan la geometría de estudio unos 7 segundos después de su inserción; mientras que las partículas más periféricas tienen una velocidad menor, y siguen en la geometría de estudio cuando la simulación termina, 10 segundos después de su inserción. Para simular la adquisición de resultados por 4D Flow mediante la simulación de trayectorias de partículas es necesario que las partículas recorran la geometría de estudio en su totalidad, lo cual no se ha conseguido mediante una única inyección en la entrada. Por tanto, es necesario realizar más
82 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID inyecciones de partículas en otros puntos. Para ello, es conveniente tener en cuenta, entre otros factores, el análisis de vórtices que se ha realizado anteriormente: las posiciones en las que hay vórtices permanentemente, como las zonas más alejadas del centro de los depósitos, requieren de una inyección de partículas propia, ya que no se pueden alcanzar inyectando partículas en otros puntos. Ansys Fluent permite realizar inyecciones multitudinarias de partículas de forma sencilla solo en las superficies o contornos delimitantes de la geometría de estudio (entrada, salidas y pared arterial). Como ya se ha visto, la inyección de partículas en la entrada no es suficiente y se requieren inyecciones en puntos intermedios. Por tanto, quedan planteadas como líneas de trabajo futuro la búsqueda de un software que permita realizar las inyecciones necesarias para cubrir la geometría de estudio completa de forma sencilla (o en caso de que no exista, su desarrollo), así como la realización de dichas inyecciones. 2.4.2 Simulación tridimensional Aproximación euleriana En este apartado se describirán el campo de velocidades y la vorticidad del fluido de la geometría tridimensional. Se comienza presentando las imágenes del campo de velocidades del fluido, así como de una vista detalle del mismo en la zona próxima a las representaciones de las intersecciones de la aorta con otras arterias, tomadas durante el último ciclo simulado, cada 0,1 segundos y con un factor de reducción de 75.
2 Desarrollo del TFG 83 UNIVERSIDAD DE VALLADOLID Figura 60. Geometría tridimensional. Campo de velocidades. Vista general. t=19 s. Figura 61. Geometría tridimensional. Campo de velocidades. Vista general. t=19,1 s.
84 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Figura 62. Geometría tridimensional. Campo de velocidades. Vista general. t=19,2 s. Figura 63. Geometría tridimensional. Campo de velocidades. Vista general. t=19,3 s.
2 Desarrollo del TFG 85 UNIVERSIDAD DE VALLADOLID Figura 64. Geometría tridimensional. Campo de velocidades. Vista general. t=19,4 s. Figura 65. Geometría tridimensional. Campo de velocidades. Vista general. t=19,5 s.
92 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Figura 78. Geometría tridimensional. Campo de velocidades. Vista detalle. t=19,8 s. Figura 79. Geometría tridimensional. Campo de velocidades. Vista detalle. t=19,9 s. Al comienzo del ciclo todo el flujo es lento, especialmente en las zonas periféricas del volumen. Se puede apreciar flujo reverso muy lento en la periferia de la representación de la aorta ascendente. En los depósitos se observa una variedad de direcciones en los vectores de velocidad que sugiere que el fluido está rotando (Figuras 60 y 70).
2 Desarrollo del TFG 93 UNIVERSIDAD DE VALLADOLID Al aumentar la velocidad en la entrada, la velocidad del fluido próximo al centro de la geometría aumenta, especialmente entre la entrada y las intersecciones. Sin embargo, apenas se aprecia diferencia en las zonas periféricas: el flujo reverso que existía anteriormente en la representación de la aorta ascendente desaparece, pero todo el fluido sigue circulando con la misma lentitud en estas zonas. La variedad en las direcciones de los vectores de velocidad en el depósito de entrada disminuye, es decir, en algunas de las zonas en las que previamente el fluido se encontraba en rotación ahora adquiere una dirección predominante de avance. Sin embargo, en el depósito de salida la situación no cambia respecto del instante anterior (Figuras 61 a 63 y 71 a 73). Cuando la velocidad en la entrada vuelve a disminuir, el flujo en la parte central decelera, el flujo reverso en la representación de la aorta ascendente vuelve a aparecer y el fluido en el depósito de entrada vuelve a adquirir una velocidad de rotación predominante sobre la de avance (Figuras 64 a 69 y 74 a 79). En resumen, en la simulación tridimensional se observa un comportamiento similar al que se observaba en la bidimensional, pero con una diferencia fundamental: los cambios de velocidad del flujo se producen únicamente en la parte central del volumen de estudio, que está rodeada de una capa de fluido en las zonas periféricas que permanece casi inmóvil incluso en los momentos en los que la velocidad en la entrada está cerca de su valor pico. La velocidad máxima en la geometría de estudio durante el ciclo es de 0,55 m/s. Al igual que en el caso bidimensional, esta velocidad es mayor que la velocidad pico en la entrada (0,5 m/s), que es la que se ha utilizado para calcular el número de Reynolds de forma provisional (1079,45). De nuevo, se aprecia sin necesidad de realizar cálculos que el nuevo Reynolds no supera la condición de régimen laminar para flujo interno en tuberías de sección circular (2300). Por tanto, se confirma que la elección del modelo de viscosidad laminar es adecuada también para el caso tridimensional. Al igual que en el caso anterior, se complementa la información proporcionada por las imágenes del campo de velocidades con imágenes del contorno de velocidades. Sin embargo, su representación en el caso tridimensional resulta difícil, ya que el color de la superficie de la geometría opaca el interior e impide visualizar correctamente todo el volumen de estudio. Lo mismo ocurre con la representación de los vórtices. Por ello, para permitir la visualización de estas magnitudes, se mostrará su representación sobre el plano de simetría de la geometría.
94 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Figura 80. Geometría tridimensional. Plano de simetría. Contorno de velocidades. t=19 s. Figura 81. Geometría tridimensional. Plano de simetría. Contorno de velocidades. t=19,1 s.
2 Desarrollo del TFG 95 UNIVERSIDAD DE VALLADOLID Figura 82. Geometría tridimensional. Plano de simetría. Contorno de velocidades. t=19,2 s. Figura 83. Geometría tridimensional. Plano de simetría. Contorno de velocidades. t=19,3 s.
96 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Figura 84. Geometría tridimensional. Plano de simetría. Contorno de velocidades. t=19,4 s. Figura 85. Geometría tridimensional. Plano de simetría. Contorno de velocidades. t=19,5 s.
2 Desarrollo del TFG 97 UNIVERSIDAD DE VALLADOLID Figura 86. Geometría tridimensional. Plano de simetría. Contorno de velocidades. t=19,6 s. Figura 87. Geometría tridimensional. Plano de simetría. Contorno de velocidades. t=19,7 s.
98 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Figura 88. Geometría tridimensional. Plano de simetría. Contorno de velocidades. t=19,8 s. Figura 89. Geometría tridimensional. Plano de simetría. Contorno de velocidades. t=19,9 s. Al principio del ciclo el flujo es lento. En la zona central del conducto de entrada, el depósito de entrada y la representación de la aorta ascendente el flujo es más rápido que en otras zonas, aunque sigue sin contar con una velocidad especialmente elevada (Figura 80).
2 Desarrollo del TFG 99 UNIVERSIDAD DE VALLADOLID Al aumentar la velocidad en la entrada, aumenta la velocidad del fluido desde la entrada hasta las intersecciones. Desde las intersecciones hasta la salida no se aprecia variación en la velocidad (Figuras 81 a 83). Al disminuir la velocidad en la entrada, disminuye la velocidad del fluido en todas las zonas en las que previamente había aumentado, exceptuando una zona de forma alargada a la salida del depósito de entrada, que mantiene una celeridad relativamente alta durante unos instantes. Después, va perdiendo poco a poco esta celeridad hasta que regresa a la situación inicial (Figuras 84 a 89). En comparación con los resultados obtenidos en la simulación bidimensional, la zona del depósito de entrada en la que el fluido mantiene cierta velocidad tras la caída de la velocidad en la entrada es mucho más excéntrica. Esto indica que en el caso tridimensional la velocidad de avance del fluido en esa zona es más significativa respecto de la velocidad de rotación que en el caso bidimensional. Además, el incremento de velocidad en la entrada en el caso tridimensional tiene un menor efecto sobre el resto del fluido que en el caso bidimensional: la velocidad en la representación de la aorta ascendente aumenta menos, y en la representación de la aorta descendente ni siquiera se aprecia una variación. Para finalizar con el análisis euleriano de la geometría tridimensional, se presentan imágenes, tomadas durante el último ciclo simulado, cada 0,1 segundos, de los vórtices en el plano de simetría del volumen de estudio, analizados según el criterio Q. El color rojo representa la presencia de vórtices, mientras que el color azul representa su ausencia. Figura 90. Geometría tridimensional. Plano de simetría. Vórtices. t=19 s.
100 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID Figura 91. Geometría tridimensional. Plano de simetría. Vórtices. t=19,1 s. Figura 92. Geometría tridimensional. Plano de simetría. Vórtices. t=19,2 s.
2 Desarrollo del TFG 101 UNIVERSIDAD DE VALLADOLID Figura 93. Geometría tridimensional. Plano de simetría. Vórtices. t=19,3 s. Figura 94. Geometría tridimensional. Plano de simetría. Vórtices. t=19,4 s.
108 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID 4 Bibliografía [1] «Wikipedia - Imagen por resonancia magnética,» 12 2022. [En línea]. Available: https://es.wikipedia.org/wiki/Imagen_por_resonancia_magn%C3%A9tica. [Último acceso: 02 2023]. [2] L. Barberini, «Introduction to the Basics of MRI to Introduce the Macroscopic Magnetization M,» de Image Principles, Neck, and the Brain, CRC Press, 2016. [3] E. Ferdian, A. Suinesiaputra, D. J. Dubowitz, D. Zhao, A. Wang, B. Cowan y A. A. Young, «4DFlowNet: Super-Resolution 4D Flow MRI Using Deep Learning and Computational Fluid Dynamics,» Frontiers in Physics, 2020. [4] «Manual MSD - Resonancia magnética nuclear,» 04 2021. [En línea]. Available: https://www.msdmanuals.com/es-es/hogar/temas-especiales/pruebas-dediagn%C3%B3stico-por-la-imagen-habituales/resonancia-magn%C3%A9tica-nuclear-rmn. [Último acceso: 02 2023]. [5] G. N. Levine, A. S. Gomes, A. E. Arai, D. A. Bluemke, S. D. Flamm, E. Kanal, W. J. Manning, E. T. Martin, J. M. Smith, N. Wilke y F. S. Shellock, «Safety of Magnetic Resonance Imaging in Patients With Cardiovascular Devices,» AHA Journals, vol. 116, nº 24, 2007. [6] M. E. Hayden y P.-J. Nacher, «History and Physical Principles of Magnetic Resonance Imaging,» de Image Principles, Neck, and the Brain, CRC Press, 2016. [7] T. Puiseux, A. Sewonu, R. Moreno, S. Méndez y F. Nicoud, «Numerical simulation of timeresolved 3D phase-contrast magnetic resonance imaging,» PLoS ONE, vol. 16, nº 3, 2021. [8] I. Arias Fernández, Desarrollo en GPU de algoritmos paralelos de procesado de imagen aplicados a la simulación del fenómeno de resonancia magnética, 2022. [9] «Wikipedia - Magnetic Resonance Angiography,» 05 2023. [En línea]. Available: https://en.wikipedia.org/wiki/Magnetic_resonance_angiography. [Último acceso: 05 2023]. [10] «Wikipedia - Phase contrast magnetic resonance imaging,» 03 2023. [En línea]. Available: https://en.wikipedia.org/wiki/Phase_contrast_magnetic_resonance_imaging. [Último acceso: 04 2023]. [11] B. Andersson, R. Andersson, L. Hakansson, M. Mortenssen, R. Sudiyo y B. van Wachem, Computational Fluid Dynamics for Engineers, Cambridge University Press, 2012.
4 Bibliografía 109 UNIVERSIDAD DE VALLADOLID [12] A. Crespo Martínez, Mecánica de Fluidos, Ediciones Paraninfo, 2006. [13] «Wikipedia - Teorema de la divergencia,» 03 2023. [En línea]. Available: https://es.wikipedia.org/wiki/Teorema_de_la_divergencia. [Último acceso: 03 2023]. [14] «Wikipedia - Fluido newtoniano,» 02 2023. [En línea]. Available: https://es.wikipedia.org/wiki/Fluido_newtoniano. [Último acceso: 02 2023]. [15] «Wikipedia - Fluido no newtoniano,» 02 2023. [En línea]. Available: https://es.wikipedia.org/wiki/Fluido_no_newtoniano. [Último acceso: 02 2023]. [16] «Wikipedia - Generalized Newtonian fluid,» 02 2023. [En línea]. Available: https://en.wikipedia.org/wiki/Generalized_Newtonian_fluid. [Último acceso: 02 2023]. [17] «Wikipedia - Carreau fluid,» 09 2022. [En línea]. Available: https://en.wikipedia.org/wiki/Carreau_fluid. [Último acceso: 02 2023]. [18] «Wikipedia - Shear rate,» 02 2023. [En línea]. Available: https://en.wikipedia.org/wiki/Shear_rate. [Último acceso: 02 2023]. [19] «Ansys Innovation Courses - 3D Bifurcating Artery - Simulating the Non-Newtonian Fluid Using Casson Fluid Model,» [En línea]. Available: https://courses.ansys.com/index.php/courses/fluent-3d-bifurcatingartery/lessons/exercise-lesson-10/topic/simulating-the-non-newtonian-fluid-usingcasson-fluid-model/. [Último acceso: 02 2023]. [20] «Hmong - Condición de Courant-Friedrichs-Lewy,» [En línea]. Available: https://hmong.es/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition. [Último acceso: 04 2023]. [21] J. W. Slater, «NASA - Examining Spatial (Grid) Convergence,» 02 2021. [En línea]. Available: https://www.grc.nasa.gov/www/wind/valid/tutorial/spatconv.html. [Último acceso: 05 2023]. [22] «Wikipedia - Vortex,» 04 2023. [En línea]. Available: https://en.wikipedia.org/wiki/Vortex. [Último acceso: 04 2023]. [23] Y. Dong, Y. Yan y C. Liu, «New visualization method for vortex structure in turbulence by,» Applied Mathematical Modelling, vol. 40, nº 1, 2016. [24] J.-m. Zhan, Z.-y. Chen, C.-w. Li, W.-q. Hu y Y.-t. Li, «Vortex identification and evolution of a jet in cross,» Engineering Applications of Computational Fluid, vol. 14, nº 1, 2020. [25] «Wikipedia - Número de Reynolds,» 09 2022. [En línea]. Available: https://es.wikipedia.org/wiki/N%C3%BAmero_de_Reynolds. [Último acceso: 06 2023].
110 Ricardo Serrano Holgado –Desarrollo de un modelo de dinámica de fluidos computacional para la corrección de flujos medidos por resonancia magnética - 4DFlow UNIVERSIDAD DE VALLADOLID [26] «Ansys Innovation Courses - 3D Bifurcating Artery - Boundry Conditions,» [En línea]. Available: https://courses.ansys.com/index.php/courses/fluent-3d-bifurcatingartery/lessons/physics-setup-lesson-5/topic/boundary-conditions/. [Último acceso: 03 2023]. [27] T.-R. Teschner, «GitHub - pyGCS (Grid Convergence Study),» 07 2022. [En línea]. Available: https://github.com/tomrobin-teschner/pyGCS. [Último acceso: 06 2023]. [28] «Wikipedia - Lagrangian and Eulerian specification of the flow field,» 10 2022. [En línea]. Available: https://en.wikipedia.org/wiki/Lagrangian_and_Eulerian_specification_of_the_flow_field. [Último acceso: 06 2023].