Introducción
La tráquea, bronquios y bronquiolos son, según [12], las principales partes que componen la zona inferior de las vías respiratorias humanas. La carina es la zona donde la tráquea se bifurca, y da comienzo a los bronquios principales derecho e izquierdo. Este punto morfológico ha sido de especial atención para estudios de fenómenos asociados a dispersión y acumulación de partículas, como los realizados por [6] y [8], y el análisis de condiciones patológicas, como lesiones traqueobronquiales [7].
El estudio numérico de la carina implica su respectiva modelación geométrica. Dos parámetros morfológicos claves en este proceso son el ángulo subcarinal y la forma de la carina. Investigaciones sobre el ángulo, como las realizadas por [9], [10] y [11], muestran cómo este valor oscila entre 30 ° y 130 °. Respecto a la forma, el estudio llevado a cabo por [1] especifica cuatro configuraciones: romo (blunt), parabólica (parabolic), punto de silla (saddle) y asimétrica (asymmetric). En ese estudio se presentan claramente los efectos sobre la aceleración del flujo, siendo el punto de silla la configuración que presenta regiones con magnitud de velocidad más elevada.
Algunos modelos de representación de vías respiratorias no consideran directamente la forma de la carina. Uno de los modelos más usados, el modelo de Weibel [13], describe su sistema de ramificaciones como tuberías, de lo que se obtiene una unión recta en las bifurcaciones. Sin embargo, algunos trabajos de investigación, como el realizado por [14], muestran que existe una relación entre el radio de redondeo de la carina y el diámetro de la rama principal, relación que no supera 0,2. Modelos de bifurcación simple, como el implementado por [2], se basan en tales consideraciones.
El principal objetivo de esta investigación es estudiar y analizar el efecto que tiene la variación del radio de redondeo de la carina en el desarrollo del flujo de aire a través de las vías respiratorias. La metodología y las ecuaciones son presentadas en la primera parte; los resultados y su respectivo análisis, en la segunda parte, y las conclusiones y referencias, en la parte final.
Métodos
Modelo computacional
El modelo sintético usado en las simulaciones computacionales está basado en el modelo de bifurcación simétrico propuesto por [2], el cual se observa en la figura 1. Las letras L y D hacen referencia a la longitud y al diámetro, respectivamente. A su vez, los subíndices p y d están asociadas a la rama principal (parent) y a las ramas hijas (daughter). Los parámetros R y a describen la curvatura de la bifurcación. El parámetro morfológico en que se enfoca esta investigación es el radio de redondeo de la carina (rc). Algunas investigaciones, como la realizada por [3], sugieren una versión adimensional de este radio, de la forma:
|
|
(1) |
Este valor de r no es mayor, generalmente, a 0,2. En la tabla 1 se presentan los valores numéricos usados en este estudio, para cada uno de los parámetros que definen la geometría. Todas las unidades son en milimetros, excepto el ángulo a (que está en grados) y el radio adimensional r. Las longitudes de las ramas se modificaron, respecto al modelo original, a cinco veces el diámetro de la ramificación, con el fin de evitar problemas numéricos asociados al desarrollo del flujo cerca de las salidas.
Figura 1. Modelo de bifurcación simple propuesto por [2]. Fuente: elaboración propia.
Tabla 1. Parámetros morfológicos del modelo sintético de vías respiratorias.
|
Dp |
Lp |
Dd |
Ld |
Rd |
a |
r |
|
16 |
80 |
14 |
70 |
64 |
35 |
0, 0.07 y 0.14 |
Fuente: elaboración propia.
La solución numérica de las ecuaciones gobernantes se obtuvo mediante el método de los volúmenes finitos, específicamente usando el software libre y de código abierto OpenFOAM. Los dominios computacionales se discretizaron mediante el enmallador no estructurado nativo de este software, snappyHexMesh. La cantidad promedio de volúmenes por cada geometría es de 2,5 millones. Como se puede observar en la figura 2, se realizó un refinamiento enfocado en el punto de unión, el cual nos permite analizar en mayor detalle el efecto del redondeo de carina.
Figura 1. Radios de redondeo. Fuente: elaboración propia.
Las ecuaciones que gobiernan el comportamiento para este tipo de flujos son las ecuaciones de Navier-Stokes. El fluido es aire, considerado como incompresible e isotérmico, teniendo así como ecuaciones gobernantes (masa y momento):
|
|
(2) |
|
|
(3) |
Donde u y v son los componentes de la velocidad, p es la presión y ν = μ/ρ es la viscosidad cinemática. Se usaron valores de propiedades del aire a temperatura ambiente (15°C): ρ = 1,23 [kg/m3] y μ = 1,8204x10-05 [Pa.s].
El solucionador nativo de OpenFoam para flujo incompresible y en estado transitorio (pimpleFoam) fue empleado para resolver las ecuaciones gobernantes. El acoplamiento entre presión y velocidad es realizado a través del algorítmo pimple, nativo de OpenFOAM. La discretización temporal fue llevada a cabo usando el esquema Crank-Nicolson.
Condiciones de frontera
La condición respiratoria analizada en este estudio es el proceso de inhalación en estado tranquilo, el cual, según [5], oscila entre 200 y 2000, números de Reynolds (Re) pertenecientes al régimen laminar. El número de Reynolds está definido como:
|
|
(4) |
Tomando como longitud característica el diámetro de la rama principal. Una condición de velocidad uniforme en la entrada es adoptada. Dos números de Reynolds dentro de la zona laminar son seleccionados (Re=500 y Re=2000). Asimismo, en las salidas se asume una condición de presión constante, así como condiciones de Neumann nulo de velocidad.
En las paredes de las vías se considera la condición de no deslizamiento.
Resultados
El monitoreo de los resultados se realiza en las líneas y secciones que se observan en la figura 3, y sobre los planos morfológicos coronal y sagital. Aquí, LLong hace referencia a la longitud total de la rama hija, la sección a se encuentra a 0,1 veces esta longitud, y la sección b a 0,5 veces, medidas a partir del inicio de la bifurcación.
Figura 3. Secciones de monitoreo. Fuente: elaboración propia.
Perfiles de velocidad
La figura 4 muestra los perfiles de velocidad promedio (Umean) contra el diámetro adimensional (d/D), medidos sobre el plano coronal para Re=500. Trabajos experimentales y numéricos, como los desarrollados por [15] y [16], evidencian cómo el perfil se deforma al pasar por la bifurcación. Se observa cómo el redondeo del radio de carina tiene una influencia directa sobre la aceleración del fluido, de lo que se obtienen mayores picos de velocidad a medida que el radio aumenta, además de influir en el estrechamiento del perfil hacia la pared interna de la rama.
Este comportamiento se mantiene para diferentes velocidades de entrada, ya que, como se observa tanto en la figura 4 como en la figura 5, para ambos números de Reynolds el fenómeno de estrechamiento se rige por el mismo patrón. A medida que se avanza a través de la rama, este efecto se va atenuando hasta casi desaparecer.
A diferencia del plano coronal, sobre el plano sagital se conserva la simetría de los perfiles de velocidad respecto
Figura 4. Perfiles de velocidad en la sección a. Re=500. Plano coronal. Fuente: elaboración propia
Figura 5. Perfiles de velocidad en la sección a. Re=2000. Plano coronal. Fuente: elaboración propia
al eje axial de las ramificaciones. La figura 6 y la figura 7 muestran cómo los perfiles toman una forma de M. Este fenómeno de deformación fue descrito por [4], y es debido al efecto de aceleración del flujo cerca a las paredes y a una desaceleración localizada en eje axial de la rama. Como se puede observar, estas variaciones de la aceleración son más pronunciadas a medida que el radio de redondeo disminuye, de lo que se obtiene una mayor depresión de la curva de velocidad cuando el radio es 0.
Vorticidad y patrones de flujo secundario
Para cuantificar el efecto del radio de redondeo sobre la vorticidad, se toma la magnitud de esta sobre una línea diametral perteneciente al plano coronal. La figura 8 muestra los perfiles de vorticidad (ω) adimensional respecto al pico (ωmax) contra el diámetro adimensional (d/D). En la figura 8-(a) y 8-(e) se observa que, cerca al centro de la rama, la magnitud de la vorticidad se aumenta a medida que el radio es mayor. Este comportamiento se va atenuando a lo largo de la ramificación, y llega a ser imperceptible en la mitad de la longitud, como se ve en la figura 8-(c) y 8-(g).
Figura 6. Perfiles de velocidad en la sección b. Re=500. Plano sagital. Fuente: elaboración propia
Figura 7. Perfiles de velocidad en la sección b. Re=2000. Plano sagital. Fuente: elaboración propia
Investigaciones relacionadas a la dispersión de partículas, como las realizadas por [8] y [17], describen cómo estas estructuras vorticales observadas a través de los flujos secundarios tienen una gran influencia en las tasas de mezcla y en el comportamiento de la deposición de las partículas. Los patrones de flujo secundario no presentan un cambio notable respecto al radio de carina. Sin embargo, se puede observar cómo el desarrollo y evolución de las estructuras vorticales dependen tanto de la longitud de las ramas como de la velocidad de entrada. Los patrones de flujo observados en la figura 8-(b) y 8-(d), pertenecientes a las secciones a y b, respectivamente, muestran cómo a medida que se avanza a través de la ramificación, los vórtices tienden a acercarse a la línea del plano coronal (guardando la simetría respecto a este plano). A su vez, al analizar la figura 8-(f) y 8-(h), se evidencia cómo los vórtices se acercan a la pared externa de la vía, es decir, se alejan del punto de velocidad máxima.
(a) Vorticidad en sección a. Re=500 (b) Patrones de flujo en sección a. Re=500
(c) Vorticidad en sección b. Re=500 (d) Patrones de flujo en sección b. Re=500
(e) Vorticidad en sección a. Re=2000 (f) Patrones de flujo en sección a. Re=2000
(g) Vorticidad en sección a. Re=2000 (h) Patrones de flujo en sección a. Re=2000
Figura 8. Comportamiento de la vorticidad y patrones de flujo secundario. Fuente: elaboración propia.
Esfuerzos cortantes
Algunos estudios, como los llevados a cabo por [18], [19] y [20], describen cómo los esfuerzos cortantes tienen relación con la deformación de las células epiteliales de las paredes de las vías respiratorias, además de influir sobre la barrera epitelial, la cual protege las vías de infecciones. En este estudio, los esfuerzos cortantes son tomados sobre la pared interna de la ramificación. La figura 9 muestra los esfuerzos cortantes sobre la pared (WSS) contra la longitud adimensional (x/L) para Re=500. Se puede observar cómo el efecto del redondeo de carina es desplazar la ubicación del esfuerzo máximo, lo que lo aleja del punto de bifurcación. Este fenómeno se conserva para distintas velocidades de entrada, como se evidencia en la figura 10, donde, para Re=2000, se mantiene ese desplazamiento. Los esfuerzos son graficados solo hasta la mitad de la rama (0,5 veces la longitud total), ya que después de este punto los esfuerzos convergen al mismo valor, sin importar el radio.
Figura 9. Esfuerzos cortantes de pared. Re=500. Fuente: elaboración propia.
Figura 10. Esfuerzos cortantes de pared. Re=2000. Fuente: elaboración propia.
Caídas de presión
En las vías respiratorias, la caída de presión está sujeta a efectos de disipación viscosa y de energía cinética [21]. Esta caída a través de la longitud de las ramas se ve ligeramente alterada por el radio de carina. En la figura [11] se observa la caída de presión (P) contra la longitud adimensional respecto a la longitud total de la rama (x/L) para Re=500. Se puede ver que en el punto exacto de unión el valor de la presión es mayor para radios más grandes de redondeo, y, a su vez, la caída de presión es más pronunciada al avanzar por la rama. Este comportamiento se mantiene para diferentes números de Reynolds, como se observa en la figura 12 (con Re=2000). Al recorrer una distancia de aproximadamente el 20 % de la longitud de la rama, las curvas de presión convergen a un mismo valor (razón por la cual se grafica solo hasta 0,2 la longitud adimensional).
Figura 11. Caída de presión. Re=500. Fuente: elaboración propia.
Figura 12. Caída de presión. Re=2000. Fuente: elaboración propia.
Conclusiones
El estudio numérico del efecto del radio de redondeo de carina sobre el flujo a través de un modelo sintético de bifurcación simple fue desarrollado con la finalidad de analizar la influencia de este parámetro morfológico en el proceso respiratorio. Tres radios de redondeo y dos velocidades de entrada fueron comparados y estudiados.
Los resultados muestran cómo el aumento del radio de redondeo de la carina tiene relación con el aumento del pico de velocidad (tanto en el plano coronal como en el sagital), así como con la magnitud de la vorticidad. A su vez, ese aumento influye en el desplazamiento del valor máximo de esfuerzo cortante de pared, y se aleja del punto de bifurcación, y en las caídas de presión, las cuales también crecen con el aumento del radio.
Es importante resaltar que, como se observó en la mayoría de los resultados, el efecto del radio de redondeo es visible en zonas cercanas a la carina, y que, al avanzar a través de la rama, la atenuación de esos efectos es notoria.
Dentro de los efectos fisiológicos que se pueden derivar a partir del efecto del radio, cabe resaltar:
- El desplazamiento de la posición del máximo esfuerzo cortante de pared indica que la deformación de las células epiteliales se aleja del punto de bifurcación a medida que aumenta el radio de carina.
- Para cualquier radio, las estructuras vorticales tienden hacia la línea coronal a medida que el flujo avanza a través de la ramificación, con lo que fenómenos de mezcla y dispersión de partículas se concentran hacia el centro de la rama.
La versatilidad de las técnicas de dinámica de fluidos computacional (CFD) para el análisis de biofluidos es notoria, desde el estudio de enfermedades y condiciones patológicas hasta la variación de parámetros morfológicos (como el aquí estudiado). La robustez de los códigos y la reducción del costo computacional que se ha venido mostrando en la actualidad ubican a la simulación numérica CFD como una de los principales métodos de investigación biomédica, para fenómenos asociados al transporte de fluidos.
Referencias
[1] Martonen, T. B., Yang, Y. and Xue, Z. Q. “Effects of carinal ridge shapes on lung airstreams”, Aerosol science and technology, vol. 21, no. 2, pp. 119-136, 1994.
[2] Lee, D., Park, S. S., Ban-Weiss, G. A., Fanucchi, M. V., Plopper, C. G. and Wexler, A. S. “Bifurcation model for characterization of pulmonary architecture”, The Anatomical Record, vol. 291, no. 4, pp. 379-389, 2008.
[3] Kang, M. Y., Hwang, J. and Lee, J. W. “Effect of geometric variations on pressure loss for a model bifurcation of the human lung airway”, Journal of biomechanics, vol. 44, no. 6, pp. 1196-1199, 2011.
[4] Schroter, R. C. and Sudlow, M. F. “Flow patterns in models of the human bronchial airways”, Respiration physiology, vol. 7, no. 3, pp. 341-355, 1969.
[5] Liu, Y., So, R. M. C. and Zhang, C. H. “Modeling the bifurcating flow in a human lung airway”, Journal of biomechanics, vol. 35, no. 4, pp. 465-473, 2002.
[6] Comer, J. K., Kleinstreuer, C. and Zhang, Z. “Flow structures and particle deposition patterns in double-bifurcation airway models. Part 1. Air flow fields”. Journal of Fluid Mechanics, vol. 435, pp. 25-54, 2001.
[7] Chu, C. P. W. and Chen, P. P. “Tracheobronchial injury secondary to blunt chest trauma: diagnosis and management”. Anaesthesia and intensive care, vol. 30, no. 2, pp. 145, 2002.
[8] Van Ertbruggen, C., Hirsch, C. and Paiva, M. “Anatomically based three-dimensional model of airways to simulate flow and particle transport using computational fluid dynamics”, Journal of applied physiology, vol. 98, no. 3, pp. 970-980, 2005.
[9] Haskin, P. H. and Goodman, L. R. “Normal tracheal bifurcation angle: a reassessment”. American Journal of Roentgenology, vol. 139, no. 5, pp. 879-882, 1982.
[10] Sahni, D., Batra, Y. K. and Rajeev, S. “Anatomical dimensions of trachea, main bronchi, subcarinal and bronchial angles in fetuses measured ex vivo”. Pediatric Anesthesia, vol. 18, no. 11, pp. 1029-1034, 2008.
[11] D'Souza, A., Ankolekar, V. H., Das, A., Padmashali, S., D’Souza, A. S. and Mamatha, H. “Determination of inter-bronchial and subcarinal angles in fetuses of different gestational age and their clinical implication”. Muller Journal of Medical Sciences and Research, vol. 6, no. 2, pp. 129-132, 2015.
[12] West, J. B. “Respiratory physiology: the essentials”. Lippincott Williams & Wilkins, MLA. 2012
[13] Weibel, E. R. “Geometric and dimensional airway models of conductive, transitory and respiratory zones of the human lung”. Morphometry of the human lung. pp. 136-142, 1963.
[14] Horsfield, K., Dart, G., Olson, D. E., Filley, G. F.and Cumming, G. “Models of the human bronchial tree”. Journal of applied physiology, vol. 31, no. 2, pp. 207-217, 1971.
[15] Zhao, Y. and Lieber, B. B. “Steady inspiratory flow in a model symmetric bifurcation”, Transactions of the ASME-K-Journal of Biomechanical Engineering, vol. 116, no. 4, pp. 488-496, 1994.
[16] Calay, R. K., Kurujareon, J. and Holdo, A. E. “Numerical simulation of respiratory flow patterns within human lung”, Respiratory physiology & neurobiology, vol. 130, no. 2, pp. 201-221, 2002.
[17] Yu, G., Zhang, Z. and Lessmann, R. “Computer simulation of the flow field and particle deposition by diffusion in a 3-D human airway bifurcation”, Aerosol Science and Technology, vol. 25, no. 3, p. 338-352, 1996.
[18] Xia, G., Tawhai, M. H., Hoffman, E. A. and Lin, C. L. “Airway wall stiffening increases peak wall shear stress: a fluid–structure interaction study in rigid and compliant airways”, Annals of biomedical engineering, vol. 38, no. 5, pp. 1836-1853, 2010.
[19] Dailey, H. L., Yalcin, H. C. and Ghadiali, S. N. “Fluid-structure modeling of flow-induced alveolar epithelial cell deformation”, Computers & structures, vol. 85, no. 11, pp. 1066-1071, 2007.
[20] Sidhaye, V. K., Schweitzer, K. S., Caterina, M. J., Shimoda, L. and King, L. S. “Shear stress regulates aquaporin-5 and airway epithelial barrier function”, Proceedings of the National Academy of Sciences, vol. 105, no. 9, pp. 3345-3350, 2008.
[21] Pedley, T. J., Schroter, R. C. and Sudlow, M. F. “The prediction of pressure drop and variation of resistance within the human bronchial airways”, Respiration physiology, vol. 9, no. 3, p. 387-405, 1970.