Anales AFA - XVI Meeting on Recent Advances of Physics of Fluids and its Applications 41-45
SIMULACIÓN NUMÉRICA DEL LLENADO DE UNA CÁMARA ELÁSTICA: UNA
PRIMERA APROXIMACIÓN AL FLUJO DIASTÓLICO EN EL VENTRÍCULO
IZQUIERDO
NUMERICAL SIMULATION OF THE FILLING OF AN ELASTIC CHAMBER: A FIRST
APPROACH TO THE DIASTOLIC FLOW IN THE LEFT VENTRICLE
J. P. G. Marí*1, M. N. Pizarro Ricciotti1y C. A. Perazzo1,2
1Facultad de Ingeniería y Ciencias Exactas y Naturales, Universidad Favaloro,
Sarmiento 1853, (1044) CABA, Argentina.
2IMeTTyB, Universidad Favaloro-CONICET,
Solís 453, (1078) CABA, Argentina.
Recibido: 31/10/21; Aceptado: 25/02/22
Gracias a la tecnología de imágenes médicas es posible conocer el campo de velocidades de la sangre en el ventrícu-
lo izquierdo con una resolución espacial y temporal cada vez mayor. El objetivo final de dicha tecnología es obtener
información de interés médico a partir del campo de velocidades. Entre otras cosas, se busca inferir las propiedades
elásticas diastólicas del músculo cardíaco, las cuales se ven alteradas en distintas patologías. Motivado por este proble-
ma y, como una primera aproximación, ignorando los detalles anatómicos y estructurales del ventrículo izquierdo, en
este trabajo se estudia numéricamente el flujo de llenado de una cámara elástica circular. Se resuelven las ecuaciones
de Navier–Stokes para el fluido dentro de la cavidad y las ecuaciones de la elasticidad para la pared de la cámara,
acopladas por las condiciones de borde en el contorno interior de la cámara, por lo que el problema a resolver es del
tipo interacción fluido–estructura. Se muestran distintas resoluciones numéricas variando sólo el módulo de Young de
la pared elástica y se analiza cómo se modifica el flujo dentro de la cavidad.
Palabras Clave: flujo diastólico, rigidez miocárdica, simulación numérica.
Owing to medical imaging technology it is now possible to ascertain the blood velocity field inside the left ventricle with
a growing spatial and temporal resolution. The long–term goal of this technology is to obtain information of medical
interest from the velocity field. Among other aims, it seeks to infer the elastic diastolic properties of the cardiac muscle,
which are altered in a variety of pathologies. Driven by this issue and, as a first approach, ignoring the anatomical
and structural details of the left ventricle, this work numerically studies the filling flow of a circular elastic chamber.
The Navier–Stokes equations for the fluid within the cavity are solved and coupled with the solutions of the elasticity
equations for the chamber wall by the boundary conditions in the internal border of the chamber, resulting in a fluid-
structure interaction problem. Several numerical resolutions are shown, adjusting merely the Young modulus of the
elastic wall and analysing how the flow inside the cavity alters.
Keywords: diastolic flow, miocardical rigidity, numerical simulation.
https://doi.org/10.31527/analesafa.2022.fluidos.41 ISSN 1850-1168 (online)
I. INTRODUCCIÓN
El aumento de la rigidez del miocardio, y la consecuen-
te pérdida de elasticidad, generan disfunción diastólica, es
decir, patologías que generan un menor llenado ventricular,
pudiendo o no estar afectada la fracción de eyección. Como
tal, la posibilidad de detectar precozmente este aumento, an-
tes de que la enfermedad esté completamente desarrollada,
podría contribuir a una mayor sobrevida de los pacientes.
Nuevas técnicas de diagnóstico por imágenes, como la
resonancia magnética 4D–flow, permiten cuantificar de ma-
nera no invasiva el campo de velocidades del flujo en todo
el ciclo cardíaco. Esto es un gran avance con respecto a las
mediciones con ecografía Doppler, ya que esta sólo brinda
información de la velocidad sobre una línea, y no sobre todo
el espacio.
Este avance en las técnicas de visualización del flujo ven-
* juannpmari@gmail.com
tricular impulsó su estudio numérico. A partir del trabajo de
Peskin et al. [1] existe una gran cantidad de publicaciones
en las que se estudia numéricamente el flujo dentro del ven-
trículo izquierdo con un grado creciente de detalle (ver por
ejemplo los trabajos de Meschini et al. [2] y Xu [3] y la
bibliografía allí citada). En todos ellos se concluye que la
característica central del flujo es un vórtice anular generado
en la válvula mitral. En la actualidad existe cierto consenso
acerca de que la dinámica de dicho anillo juega un rol fisio-
lógico y puede modificarse en ciertas patologías. De hecho,
existen estudios que analizan distintos descriptores del vór-
tice anular [4] en pacientes sanos y con cardiomiopatías. Sin
embargo, la mayoría de la literatura previa centra su aten-
ción en la dinámica de la válvula mitral, dado que es donde
se genera el vórtice, ya sea en condiciones normales, pato-
lógicas e incluso cuando es artificial. Por otra parte, a pesar
de la cantidad de simulaciones numéricas existentes del flu-
©2022 Anales AFA 41
FIG. 1: Geometría utilizada. Gris oscuro representa al miocardio
y gris claro a la sangre. Los dos segmentos rectos representan las
valvas.
jo ventricular [5], no hay estudios que investiguen el modo
en que cambios en la rigidez miocárdica altera al flujo. Da-
da esta situación, el objetivo de este trabajo es precisamen-
te estudiar la influencia de la rigidez del miocardio sobre
el flujo diastólico. Para ello, se generó un modelo simple
del ventrículo izquierdo durante la diástole. Como primera
aproximación, el modelo generado es en dos dimensiones y
con una geometría simplificada, pero basado en parámetros
estructurales y dinámicos reales. Se realizaron varias simu-
laciones que difieren sólo en el módulo de Young del sólido
y se analizan las modificaciones en el flujo. Los resultados
aquí presentados muestran que para el modelo utilizado es
posible inferir la rigidez de la pared sólida de la cámara a
partir del campo de velocidades durante su llenado. Esto su-
giere por primera vez que quizás también pueda ser posible
hacerlo en un ventrículo a partir de imágenes de resonan-
cia, y así proveer de alguna herramienta de utilidad clínica
para detectar precozmente el aumento de rigidez asociado a
patologías como dilatación o hipertrofia del miocardio.
II. PROBLEMA A RESOLVER
Geometría y materiales
La geometría utilizada se presenta en la Fig. 1. Consis-
te de dos circunferencias concéntricas, cuya diferencia de
radio representa el espesor del miocardio. Hay una entrada
rectangular, que emula el orificio de la válvula mitral por
donde ingresa la sangre, y dos segmentos rectos que repre-
sentan las valvas (asumidas rígidas). Esta forma y dimen-
siones corresponden al ventrículo al final de la relajación
isovolumétrica, antes de comenzar el llenado.
Como aproximación, se asumió al miocardio como un
material linealmente elástico e isótropo, con un módulo de
Young Econstante a lo largo de toda la diástole. Se asu-
me además a la sangre como Newtoniana e incompresible.
Los valores de los parámetros utilizados pueden verse en la
Tabla 1.
Ecuaciones que gobiernan el problema
El flujo está gobernado por la ecuación de Navier–Stokes
incompresible
ρs
∂u
∂t+ρs(u·∇)u=−∇p+µs∇2u,(1)
∇·u=0,(2)
TABLA 1: Parámetros geométricos y estructurales.
Parámetro Valor
Radio interior[6]ri0.0325 m
Espesor[7]e0.01091 m
Diámetro válvula mitral[2]Dv0.024 m
Longitud de las valvas[8]Lv0.023 m
Densidad de la sangre[9]ρs1.06 ×103kg/m3
Densidad del miocardio[9]ρm1.37 ×103kg/m3
Viscosidad de la sangre[9]µs4.71 ×10−3Pa s
Módulo de Poisson del mio-
cardio νm
0.4
FIG. 2: La velocidad a la entrada de la cavidad es U f (t). En esta
figura se muestra f (t).
donde uyprepresentan el campo de velocidades y la pre-
sión del fluido, respectivamente.
A la entrada de la cavidad se asumió un perfil de velo-
cidades uniforme cuyo módulo está dado por U f (t), con
U=0,8 m/s, donde f(t)es un pulso Gaussiano [10] cu-
yo máximo está en tp=0,09 s y tiene desvío estándar
sd =0.027 s para 0 ≤t≤2tp, y f(t) = 0 si t≥2tp(ver
Fig. 2). Este flujo de entrada simula la onda E (early filling)
del llenado diastólico.
Por otra parte, la mecánica del sólido elástico está gober-
nada por la ecuación
ρm
∂2v
∂t2=∇(FS),(3)
dondeves el desplazamiento del sólido elástico, F=I+∇v
es el gradiente de deformación, Ses el segundo tensor de
Piola–Kirchoff y ρmes la densidad inicial del sólido. Notar
que la divergencia ∇(FS)está calculada con respecto a las
coordenadas materiales.
Dado que se asume que el sólido es lineal, homogéneo e
isótropo, la relación entre Sy el tensor de deformaciones de
Green–Lagrange, dado por
ε=1
2∇v+ (∇v)T+ (∇v)T∇v,(4)
está determinada por la ley de Hooke
S11
S22
S33
S23
S13
S12
=E
(1+νm)(1−2νm)D
ε11
ε22
ε33
2ε23
2ε13
2ε12
,(5)
Marí et al. / Anales AFA - XVI Meeting on Recent Advances of Physics of Fluids and its Applications 41-45 42
FIG. 3: Módulo de la velocidad para E =10 kPa y E =70 kPa
para tiempos tp,2tpy3tp. Para facilitar la comparación entre los
flujos a diferente E, la escala de colores es la misma para ambos
E a un tiempo, pero es diferente en cada tiempo. Las escalas es-
paciales son las mismas para todas las figuras.
D=
1+νmνmνm000
νm1−νmνm000
νmνm1−νm000
0001−2νm
20 0
0 0 0 0 1−2νm
20
0 0 0 0 0 1−2νm
2
.(6)
Como condiciones iniciales se asumióu=0 y p=0 para
el fluido y ∂v
∂t=0 y ∂2v
∂2t=0 para el sólido, y como condición
de contorno se utilizó la condición de no deslizamiento, es
decir, la velocidad del fluido en contacto con la pared inte-
rior de la cavidad es igual en todo momento a la velocidad
de dicha pared en todos sus puntos.
Las Ecs. (1) - (6) y las condiciones iniciales y de contorno
constituyen el modelo a resolver. Este problema es del tipo
interacción fluido–estructura, en el que el fluido afecta a la
estructura y viceversa. La interacción entre el fluido y la
estructura ocurre en la pared interior de la cavidad.
Todas las longitudes fueron adimensionalizadas con ri,
las velocidades con U, el tiempo t,tpysd con ri/Uy
los esfuerzos y Econ ρsU2. Las cantidades adimensiona-
les se representan con el mismo símbolo que las respectivas
con unidades, dado que de ahora en más todas las canti-
dades informadas son adimensionales, excepto Edel cual
se retendrán sus unidades. El número de Reynolds resultó
Re =ρsUri
µs=5851.
Método numérico
Es posible plantear las ecuaciones del problema basadas
en un sistema de coordenadas espacial, que se mantiene fijo
en el espacio, o en un sistema material, que se mueve junto
con el material al deformarse. Otra opción es usar una com-
binación de ambos sistemas de referencia, lo que se conoce
como método Arbitrary Lagrangian–Eulerian. Este combi-
na las ventajas de ambos métodos por separado, y es el que
se utilizó en el presente trabajo, ya que es ideal para apli-
caciones de interacción fluido–estructura y es ampliamente
utilizado en esta clase de problemas [10].
Se recurrió a un enfoque fully coupled, que implica que
se resuelvan todas las ecuaciones simultáneamente para to-
das las variables, incluyendo todos los acoplamientos entre
estas, en una sola iteración. Debido a su mayor robustez se
eligió está opción frente al enfoque segregado, que consiste
en dividir el problema en distintos pasos a resolver sucesi-
vamente.
Para la discretización espacial se utilizó el método finite
element method, con una red de 19661 elementos triangula-
res y 10340 nodos, y se resolvieron las ecuaciones usando el
método de Newton–Raphson. Debido a la no linealidad del
sistema de ecuaciones, se requieren múltiples pasos de ite-
ración hasta llegar a una solución aceptable, tomando como
criterio de terminación un error relativo menor a 0.005.
En cada una de estas iteraciones se resuelve un sistema
linealizado, para lo cual hay dos algoritmos que se pueden
usar, uno directo y uno iterativo [10]. El directo requiere
más memoria pero es más robusto en geometrías complejas,
por lo que fue el elegido.
Los cálculos fueron realizados en un procesador Intel i7
de dos núcleos y 8 Gb de RAM, tomando alrededor 30 mi-
nutos cada simulación.
Los datos de ambas componentes de la velocidad y de
presión de cada solución numérica para cada paso temporal
fueron posteriormente recalculados por interpolación sobre
una grilla rectangular con un paso espacial adimensional
de 0,005 en ambas direcciones (la de entrada del flujo y la
transversal). Para el cálculo de la vorticidad, las derivadas
fueran calculadas por diferencia finitas a segundo orden.
III. RESULTADOS
Las ecuaciones de la sección anterior se resolvieron nu-
méricamente para un tiempo entre 0 y 3tp, y se registró la
solución en 75 instantes de tiempo equidistantes. Como re-
sultado, se obtuvieron los valores de p,uy la forma del
ventrículo.
Dado que el objetivo del trabajo fue analizar los cambios
en los patrones del flujo al modificar la rigidez miocárdica,
se repitió esta resolución para nueve valores de Eentre 10
kPa y 70 kPa (alrededor de 20 kPa se considera normal [2]),
manteniendo todos los demás parámetros constantes. En las
Figs. 3y4se muestran el módulo de la velocidad y la vor-
ticidad respectivamente para E=10 kPa y E=70 kPa para
tiempos (adimensionales) tp, 2tpy 3tp.
IV. DISCUSIÓN Y CONCLUSIONES
Las simulaciones muestran que lo más destacado del flu-
jo es un par de vórtices de diferente signo generado en los
extremos de la valva y que se desplaza en la cavidad (esto
se corresponde con un vórtice anular en casos con simetría
axial [9] o tridimensionales [2,3], pero con diferente di-
námica). El análisis no muestra cambios apreciables en su
evolución al variar E. Por otra parte, la evolución del má-
ximo de |u|muestra grandes diferencias a diferentes Esi
t≳3 (Fig. 5).
Marí et al. / Anales AFA - XVI Meeting on Recent Advances of Physics of Fluids and its Applications 41-45 43
FIG. 4: Vorticidad para E =10 kPa y E =70 kPa para tiempos tp,
2tpy3tp. Para facilitar la comparación entre los flujos a diferente
E, la escala de colores es la misma para ambos E a un tiempo,
pero es diferente en cada tiempo. Las escalas espaciales son las
mismas para todas las figuras.
FIG. 5: Máximo de la magnitud de la velocidad como función del t
(adimensionalizado con ri/U) para diferentes valores del módulo
de Young. La línea discontinua es la velocidad a la entrada de la
cavidad.
Las diferencias más notables al variar Ese observan en la
longitud de la cavidad en la dirección en que ingresa el flu-
jo. En la Fig. 6se aprecia que todas las curvas inicialmente
crecen hasta alcanzar un máximo, pero este valor máximo y
el tiempo en el que se llega a él decrecen con E. A partir del
máximo, la longitud para los valores más chicos de Edis-
minuye monótonamente, pero para valores más grandes de
Ela longitud disminuye hasta alcanzar un mínimo a partir
del cual vuelve a aumentar.
Las Figs. 5y6muestran que en el modelo aquí adoptado
es posible inferir el valor de Ea partir de información del
flujo de llenado de la cámara. Por ejemplo, en la Fig. 7se
muestran como función de Eel valor sde la pendiente en
t=3.2 de cada curva de la Fig. 5, el valor Lmdel primer
máximo relativo de cada curva de la Fig. 6y el tiempo tmen
FIG. 6: Longitud en la dirección en la que entra el flujo para dife-
rentes valores del módulo de Young.
FIG. 7: s, Lmy tm(definidos en el texto) vs. E.
el cual se alcanza Lm. Dado que estas curvas son monótonas
(al menos para el rango de Eaquí explorados), si se tuviera
una simulación numérica del modelo del cual se desconoce
E, a partir de s,Lmytmse puede obtener una buena estima-
ción de su valor. Lo destacable de esto es que sugiere que
es posible detectar cambios en la rigidez de un ventrículo a
partir del flujo de llenado. Sin embargo, para confirmar esto
es necesario un modelo de ventrículo mucho más detallado.
El modelo aquí usado es el más sencillo posible de acuer-
do al objetivo de este trabajo, por lo que, aunque esté moti-
vado por el flujo diastólico ventricular, no es aceptable co-
mo modelo de un ventrículo real, por lo cual es inviable
Marí et al. / Anales AFA - XVI Meeting on Recent Advances of Physics of Fluids and its Applications 41-45 44
comparar estas simulaciones con imágenes médicas o con
simulaciones sumamente realistas [2,3]. Por otra parte, la
sencillez extrema de modelo permite mantener un bajo cos-
to computacional y reconocer fácilmente la modificaciones
en el flujo al cambiar la rigidez de la cámara.
Como continuación de este trabajo y con el objetivo de
trasladar sus resultados a un ventrículo, se debe modificar
al modelo para hacerlo gradualmente más realista y anali-
zar de qué modo estas modificaciones alteran sus resulta-
dos. Algunas de estas modificaciones serían, por ejemplo,
valvas de diferente longitud, que el flujo ingrese a la ca-
vidad con una pequeña inclinación, adoptar una geometría
más realista, pasar a un modelo tridimensional, mejorar el
modelo mecánico del miocardio, etc. El objetivo final y a
largo plazo, del cual el presente trabajo es sólo la etapa ini-
cial, es identificar patrones en el flujo durante el llenado
diastólico que se vean modificados al variar la rigidez del
miocardio, para luego tratar de hallar la misma correlación
en imágenes 4D–flow de pacientes. Esto sería de gran utili-
dad para permitir al médico detectar precozmente posibles
patologías, y así plantear un tratamiento más eficaz desde el
comienzo.
AGRADECIMIENTOS
JPGM y MNPR agradecen a la Universidad Favaloro sus
Becas de Investigación para Estudiantes. CAP agradece los
subsidios PICT 2019-2674 y 2019-1563 de la ANPCyT.
REFERENCIAS
[1] C. S. Peskin y D. M. McQueen. Modeling prosthetic heart
valves for numerical analysis of blood flow in the heart. J
Comput Phys 37, 113-132 (1980).
[2] V. Meschini, M. De Tullio, G. Querzoli y R. Verzicco.
Flow structure in healthy and pathological left ventricles
with natural and prosthetic mitral valves. J Fluid Mech 834,
271-307 (2018).
[3] F. Xu y S. Kenjereš. Numerical simulations of flow patterns
in the human left ventricle model with a novel dynamic
mesh morphing approach based on radial basis function.
Comput Biol Med 130, 104184 (2021).
[4] J. Mangual, E. Kraigher-Krainer, A. De Luca, L. Tonce-
lli, A. Shah, S. Solomon, G. Galanti, F. Domenichini y G.
Pedrizzetti. Comparative numerical study on left ventricu-
lar fluid dynamics after dilated cardiomyopathy. J Biomech
46, 1611-1617 (2013).
[5] S. N. Doost, D. Ghista, B. Su, L. Zhong e Y. S. Morsi.
Heart blood flow simulation: a perspective review. Biome-
dical engineering online 15, 1-28 (2016).
[6] M. Di Donato, P. Dabic, S. Castelvecchio, C. Santambro-
gio, J. Brankovic, L. Collarini, T. Joussef, A. Frigiola, G.
Buckberg y L. Menicanti. Left ventricular geometry in nor-
mal and post-anterior myocardial infarction patients: sphe-
ricity index and ‘new’ conicity index comparisons. Eur Jour
Cardio-Thorac 29, S225-S230 (2006).
[7] M. Lorente, C. Escalona, M. Zabalza-Cerdeiriña y J.
Álvarez-Moro. Left ventricle morphometry in healthy hu-
mans. Long axis, contrast enhanced CT study. Scientific
Medical Data (2017).
[8] K. Nomura, Y. Ajiro, S. Nakano, M. Matsushima, Y. Ya-
maguchi, N. Hatakeyama, M. Ohata, M. Sakuma, T. No-
naka, M. Harii, M. Utsumi, K. Sakamoto, K. Iwade y N.
Kuninaka. Characteristics of mitral valve leaflet length in
patients with pectus excavatum: A single center cross-
sectional study. PLOS ONE 14, e0212165 (2019).
[9] H. Watanabe, S. Sugiura, H. Kafuku y T. Hisada. Multiphy-
sics simulation of left ventricular filling dynamics using
fluid-structure interaction finite element method. Biophys
J87, 2074-2085 (2004).
[10] Y. Cheng, H. Oertel y T. Schenkel. Fluid–structure coupled
CFD simulation of the left ventricular flow during filling
phase. Ann Biomed Eng 33, 567-576 (2005).
Marí et al. / Anales AFA - XVI Meeting on Recent Advances of Physics of Fluids and its Applications 41-45 45
