Le malattie cardiovascolari sono causa di numerosi decessi. Tra queste risultano molto pericolosi gli aneurismi dell’aorta addominale per i quali ancora non sono state sviluppate tecniche adeguate per la diagnosi e la prevenzione. L’implementazione di analisi numeriche può aiutare ad analizzare l’influenza dei parametri geometrici e materiali sul rischio di rottura dell’aneurisma.
In tale report si considera un modello semplificato dal punto vista geometrico, emodinamico e strutturale e vengono condotte alcune simulazioni di fluidodinamica numerica e analisi strutturale su una geometria rappresentante un vaso aneurismatico.
I risultati dell’analisi strutturale confermano che la sollecitazione massima aumenta con il diametro e con la riduzione dello spessore del vaso.
Il seguente articolo è un estratto di un progetto realizzato per il corso di Modellazione e Simulazione di Sistemi Fisiologici. Tenuto per Ingegneria Medica presso l’Università degli Studi di Roma Tor Vergata.
Introduzione
Negli ultimi decenni gran parte della ricerca sul sistema cardiovascolare si è concentrata sul rapporto tra l’emodinamica locale e l’insorgenza di fenomeni fisiopatologici. Le malattie cardiovascolari sono causa di circa il 48% dei decessi annuali in Europa portando a spese oltre i 100 miliardi di euro [1].
Oggi è ritenuto che l’insorgenza di alcune patologie sia legata a particolari pattern fluidodinamici alterati. Ad esempio, è ritenuto che dove vi sia ristagno di sangue legato a basse velocità questo incide maggiormente sull’insorgenza di processi infiammatori.
L’emodinamica del sistema cardiovascolare è molto complessa legata alla natura turbolenta del flusso dove anche un piccolo cambiamento del flusso sanguigno influenza la biomeccanica di cuore e vasi.
Nel caso degli aneurismi questa possibilità è fortemente sostenuta dal fatto che non sono presenti metodi quantitativi certi utilizzati per validarne concretamente il rischio di rottura e questo rimane ancora un problema aperto. Nonostante i forti sviluppi nelle tecniche chirurgiche e di imaging diagnostico i tassi di mortalità e morbilità rimangono molto elevati.
Attualmente in letteratura si è concordi nel cercare un criterio di valutazione affidabile associato alla potenziale rottura e non alla dimensione massima. Inoltre, essendo un fenomeno che sottopone la parete a forte stress meccanico è fortemente probabile che un criterio basato sulla quantificazione degli stress sia ottimale [2]. Tuttavia, attualmente non c’è un metodo in grado di ottenere misurazioni in vivo valide.
In questo report viene considerato un modello fortemente semplificato dal punto vista geometrico, emodinamico e strutturale e vengono condotte alcune simulazioni numeriche. Viene dunque condotta sia una simulazione fludiodinamica all’interno di una geometria rappresentante un vaso aneurismatico sia un’analisi di interazione fluido struttura disaccoppiata.
![FIG. 1: Ricostruzione 3D dell’aorta (a) e dell’aorta aneurismatica (b)[3]. Approssimazione geometrica dell’anueurisma utilizzata per l’analisi strutturale, il dominio interno è stato invece utilizzato per l’analisi CFD (c)](https://alessandromastrofini.it/wp-content/uploads/2022/06/image-1024x572.png)
Background
Un aneurisma è una dilatazione patologica di un vaso sanguigno, generalmente di un’arteria. Ne derivano diversi problemi e rischi per la salute, tipicamente legati all’indebolimento della parete dilatata che rischia di rompersi portando ad emorragia massiva e conseguente morte.
L’aorta è il vaso arterioso più grande e importante del corpo umano che, originando dal cuore, diffonde il sangue verso ogni organo e distretto. Esistono diversi siti di localizzazione di un aneurisma ma tipicamente si parla di aneurisma addominale, toracico o intracraniale.
Aneurisma addominale
L’aneurisma dell’aorta addominale (AAA) è una dilatazione localizzata dell’aorta infrarenale. L’AAA deriva da cambiamenti della struttura della parete aortica, compresa la perdita di cellule muscolari lisce e la degradazione della matrice extracellulare. La parte rigonfia è fragile e rischia di rompersi con relativa facilità.
In condizioni normali il diametro dell’aorta addominale è circa di 2 cm e si parla di aneurisma quando il rigonfiamento è superiore a 3 cm. L’AAA viene ritenuto operabile quando raggiunge un diametro di 5,5 cm [4].
La causa di origine è tuttora sconosciuta ma la letteratura attuale concorda nel riconoscere il ruolo fondamentale di diversi fattori quali l’invecchiamento, l’aterosclerosi, l’ipertensione, il fumo di sigaretta, la vasculite e la predisposizione genetica.
Nell’ultima decade sono stati proposti numerosi modelli numerici per valutare e predire la rottura dell’AAA mostrando come vanno portate in conto diverse informazioni quali lo spessore di parete, la presenza di trombi intraluminali, l’omogeneità della parete e l’uso di un modello costitutivo appropriato [5].
Analisi
In questo report si affrontano alcuni aspetti della biomeccanica dell’AAA. Seguono diverse semplificazioni, sia dal punto di vista emodinamico che strutturale. Viene quindi considerata un’approssimazione geometrica di un AAA fusiforme rappresentata in fig. 1c e descritta più avanti sulla quale vengono effettuate diverse simulazioni volte a considerarne il comportamento biomeccanico del flusso sanguigno e della parete vascolare (inizialmente considerata rigida).
Tali analisi numeriche vengono implementatenell’ambente simulativo Comsol Multiphysics.
Moto alla Poiseuille
Affrontando la descrizione del moto del sangue possiamo parlare del modello di Navier Stokes per fluido viscoso alla Newton incomprimibile dove, indicato con 𝐮, il campo di velocità, vale:
\begin{equation}
\left\{\begin{array}{l}
\rho \left(\frac{\partial \mathbf{u}}{\partial t}+ \nabla \mathbf{u} \cdot \mathbf{u}\right)=-\nabla p+\mu \Delta_{2} \mathbf{u}+\mathbf{f} \\
\nabla \cdot \mathbf{u}=0
\end{array}\right.
\end{equation}Tali equazioni rappresentano proprio il secondo principio della dinamica e la conservazione della massa.
Cilindro cavo
Nel moto alla Poiseuille si considera un moto unidirezionale lungo l’asse del cilindro del tipo:
\begin{equation}
\mathbf{u}=(u, 0,0)
\end{equation}Tale moto, insieme all’equazione di continuità, esprimono come la derivata della velocità in direzione longitudinale sia nulla:
\begin{equation}
\frac{\partial u}{\partial x}=0
\end{equation}Inoltre, si trascurano le forze di volume 𝐟 e si considera un caso stazionario. Queste ipotesi portano a semplificare notevolmente le equazioni di Navier-Stokes. In particolare, viene meno il contributo convettivo, fonte di non linearità. Si semplificano notevolmente le dipendenze spaziali e, considerando coordinate cilindriche con geometria e carichi simmetrici rispetto l’anomalia, le variabili rimangono:
\begin{equation}
\begin{aligned}
p(x, y, z,) & \rightarrow p(x) \\
\mathbf{u}(x, y, z) & \rightarrow u(r)
\end{aligned}
\end{equation}Ovvero:
\begin{equation}
\mu \Delta_{2} u=\frac{\partial p}{\partial x}
\end{equation}Profilo di velocità
Tale equazione, integrata due volte e considerando la condizione di velocità nulla al bordo, si traduce nella soluzione analitica del moto alla Poiseuille:
\begin{equation}
u(r)=\frac{1}{4 \mu} \frac{P_{L}-P_{0}}{L}\left(r^{2}-R^{2}\right)
\end{equation}Ovvero è presente un andamento parabolico dove la velocità è massima sull’asse del canale:
\begin{equation}
u(r)=\frac{1}{4 \mu} \frac{\Delta p}{L} R^{2}
\end{equation}
Tale problema si affronta sia analiticamente, graficandone l’andamento, sia numericamente. Il confronto è presente in fig. 2a.
Nella modellazione CFD viene considerato un canale largo 2 cm con una lunghezza di 20 cm [2] e una pressione pari a circa 90 mmHg [6], ovvero 12000 Pa, con una caduta di pressione di 2 Pa [7]. Tali parametri sono rappresentativi dell’aorta intrarenale.
Vaso aneurismatico
Si varia il dominio fluido andando ad approssimare un AAA mediante una dominio assialsimmetrico. Tale dominio si definisce a partire dal bordo del precedente dominio fluido aggiungendo, all’interno dell’ambiente Comsol, una sfera di raggio tale da restituire complessivamente il raggio aneurismatico. Ovvero, sia 𝚍 il diametro del vaso e 𝑎𝑛𝑒𝑢𝑟𝑦𝑠𝑚_𝐷 il diametro dell’aneurisma, la sfera si disegna tramite una semicirconferenza di raggio (𝚊𝚗𝚎𝚞𝚛𝚢𝚜𝚖_𝙳 − 𝚍)∕𝟸 centrata sulla parte esterna del vaso, a metà altezza.
Pressione interna
Successivamente le due geometrie sono unite tramite un’operazione booleana e i punti di raccordo vengo smussati mediante un raccordo con raggio pari al 40% del raggio dell’aneurisma. Questo permette sia di avere una geometria più realistica sia di evitare singolarità numeriche legate alla discontinuità.
Valgono condizioni analoghe al caso precedente per l’inlet, l’outlet e la parete con condizione no slip. Tuttavia, l’inlet si considera nella sezione superiore e l’outlet nella sezione inferiore, facendo riferimento alla situazione anatomica.
Le linee di flusso sono presenti in fig. 3 e un render tridimensionale è presente in appendice.
Si estraggono anche ulteriori dati come la pressione lungo il raggio del vaso a differenti altezze, in fig. 4 e la pressione sulla parete, fig. 5. Tali informazioni possono essere estratte sia tramite interfaccia grafica sia all’interno di Matlab (vedi appendice).
Analisi strutturale
Noto il carico pressorio viene condotta un’analisi strutturale riportando tale carico come condizione al bordo. Quindi, le due simulazioni rimangono disaccoppiate ma vengono riportati i carichi pressori.
Per implementare l’analisi è necessario definire un dominio strutturale. Dunque si implementa una procedura, all’interno di Comsol, tale da mantenere il bordo della geometria interna (precedentemente descritta per il dominio fluido) ma restituire un guscio di spessore costante. Viene quindi aperto il dominio lasciando solo la parete e poi quest’ultima viene inspessita verso l’esterno con uno spessore di 1.5 mm [8].
Il carico lo si definisce tramite una funzione per punti a partire dai dati numerici estratti dalla precedente simulazione. Sono caricati i dati tabulati e poi interpolati tramite Comsol fornendo una funzione 𝑝𝑤𝑎𝑙𝑙(𝑧) utilizzata come condizione di bordo sulla parete interna del vaso.
Viene allora utilizzato un modello iperelastico Neo-Hookean quasi incomprimibile con parametri 𝜇 = 1.7 ⋅ 106 Pa, 𝑘 = 34,7 ⋅ 1010 Pa [9] e 𝜌 = 104 kg/m3. Il modello Neo-Hookean è un modello di materiale iperelastico che può essere utilizzato per descrivere la relazione tensione-deformazione per un materiale con comportamento non lineare in elasticità finita. In particolare, tale modello può descrivere con sufficiente accuratezza il comportamento dai vasi sanguigni [10].
Neo Hooke
In elasticità finita il materiale si definisce tramite la densità di energia di deformazione, funzione dello stato elastico elastico tipicamente tale da venire descritta tramite il tensore destro di Cauchy-Green.
Questo permette di calcolare lo stress in coordinate locali tramite il secondo tensore di Piola- kirchoff:
\begin{equation}
\mathbb{S}=2 \frac{\partial W_{s}}{\partial \mathbb{C}}
\end{equation}Che risulta legato agli stress come:
\begin{equation}
\boldsymbol{\sigma}=\boldsymbol{J}^{-1} \mathbb{F} S \mathbb{F}^{T}
\end{equation}Dove:
\begin{equation}
J=\operatorname{det} \mathbb{F}
\end{equation}Oppure in funzione del tensore di Green- Lagrange:
\begin{equation}
\mathbb{E}=\frac{\mathbb{C}-\mathbb{I}}{2}
\end{equation}Dove:
\begin{equation}
\mathbb{E}=\frac{1}{2}\left(\mathbb{F}^{T} \mathbb{F}-\mathbb{I}\right)
\end{equation}Ed F è il gradiente di deformazione:
\begin{equation}
\mathbb{F}=\mathbb{I}+\nabla \mathbf{u}
\end{equation}Energia di deformazione
Nel modello Neo-Hookean [11] viene dunque (10) definita la densità di energia di deformazione in funzione di un contributo volumetrico e un contributo isocoro :
\begin{equation}
W_{\mathrm{s}}=W_{\text {iso }}+W_{\mathrm{vol}}
\end{equation}Dove il contributo volumetrico viene definito in funzione del modulo di Bulk e della deformazione elastica volumetrica 𝐽𝑒𝑙 come:
\begin{equation}
W_{\text {vol }}\left(J_{{el}}\right)=\frac{1}{2} \kappa\left(J_{\mathrm{el}}-\mathbb{I}\right)^{2}
\end{equation}Mentre il contributo isocoro si definisce in funzione del primo invariante isocoro e del secondo parametro di Lamè:
\begin{equation}
W_{\text {iso }}=\frac{1}{2} \mu\left(\overline{\mathbb{I}_{1}}-3\right)
\end{equation}Parete vascolare
Tale modello costitutivo permette di descrivere il guscio rappresentante la parete del vaso. Noto il campo di pressione agente sulla parete interna del vaso, ottenuto dall’analisi fluidodinamica (fig. 5) viene applicato come condizione al bordo nell’analisi strutturale.
Dunque si analizza lo spostamento presente in fig. 7b e lo stress a parete rappresentato in fig. 6. Lo stress viene calcolato lungo le tre direzioni lungo un’ascissa curvilinea passante sul bordo interno della parete. Ne risulta quindi un valore massimo nella direzione tangenziale di 0.2 MPa mentre in direzione radiale è inferiore a 0.5 MPa.
Lo spostamento massimo è inferiore al millimetro ed è concentrato nella zona di raccordo tra la geometria cilindrica e la sfera.
Dalla fig. 7a è possibile osservare la risultante degli stress secondo Von Mises e si osserva come c’è una maggiore concentrazione degli stress sulla zona di raccordo tra circonferenza e parete del vaso. Tale zona è di fatto quella maggiormente sottoposta ad allungamento come è possibile osservare dal diagramma a falsi colori in fig. 7a.
Tali risultati, sebbene frutto di un’analisi disaccoppiata e molto semplificata, mostrano come la geometria rigonfia ha una forte incidenza sull’amplificazione degli stress.
Variabilità geometrica
Dal modello Comsol si passa quindi al modello corrispondente in Matlab tramite il quale è possibile estrarne valori e parametri di interesse. In particolare viene ripetuta l’analisi fluidodinamica a partire da una variazione del diametro del vaso e dell’aneurisma secondo valori fisiopatologici [12, 13].
Viene dunque estratto il campo di pressione sulla parete interna e viene passato ad una simulazione strutturale basata su una geometria variata analogamente. I risultati sono presenti in fig. 9 dove si considera lo stress massimo sulla parete in funzione dei parametri geometrici.
In particolare, tali risultati sono in accordo con la letteratura mostrando come al crescere del diametro dell’aneurisma cresce il valore dello stress massimo. Per lo spessore di parete l’andamento è inverso, ovvero più la parete è sottile più lo stress sarà elevato.
Variazione pressoria
A seguire si implementa una simulazione fluidodinamica non stazionaria, frutto di un tentativo di FSI full-coupled non convergente. In tale simulazione si considerano un inlet e un outlet variabili nel tempo. Per l’inlet si sceglie una curva di velocità e per l’outlet una curva di pressione [14, 2] tipiche per il compartimento arterioso di interesse (indicate in appendice).
Analisi tempo dipendente
Tali curve sono definite in Comsol per punti e poi viene affidata a Comsol l’interpolazione e la trasformazione in una funzione del tempo. Il risultato, in termini di pressione sulla parete, è presente in fig. 10 dove si mostra a diversi instanti temporali. Nonostante le curve sembrino piatte presentano, osservandole nella scala corretta (in Pascal), un andamento simile al caso stazionario.
Si seleziona quella all’instante di tempo dove si raggiunge la pressione massima (0.5 s) e viene nuovamente utilizzata per l’analisi strutturale. I risultati (fig. 11) mostrano come, essendo aumentata la pressione media, aumenta lo stress sulla parete pur presentando un andamento nello spazio molto simile al caso precedente.
Conclusioni
I risultati dell’analisi confermano che la sollecitazione massima aumenta con il diametro dell’aneurisma e al ridursi dello spessore del vaso. Tale sollecitazione inoltre, risulta fortemente legata alla curva di pressione. Le irregolarità geometriche fornisco anche un’amplificazione delle tensioni sulla parete aneurismatica.
Inoltre, l’implementazione di un modello FSI porterebbe a risultati più precisi in quanto permetterebbe di portare in conto l’elevata deformabilità dei vasi che conferisce complianza al sistema cardiocircolatorio.
Riferimenti
- [1] Michal Tendera. “How much does Europe invest in the treatment of cardiovascular diseases?The opinions ex- pressed in this article are not necessarily those of the Editors of the European Heart Journal or of the Eu- ropean Society of Cardiology.”
- [2] Christine M. Scotti et al. “Wall stress and flow dynam- ics in abdominal aortic aneurysms: finite element anal- ysis vs. fluid–structure interaction”. en.
- [3] Nathan M. Wilson, Ana K. Ortiz, and Allison B. Johnson. “The Vascular Model Repository: A Pub- lic Resource of Medical Imaging Data and Blood Flow Simulation Results”.
- [4] Natzi Sakalihasan et al. “Abdominal aortic aneurysms”. en. In: Nature Reviews Disease Primers
- [5] Mamadou Toungara, Gregory Chagnon, and Chris- tian Geindreau. “NUMERICAL ANALYSIS OF THE WALL STRESS IN ABDOMINAL AORTIC ANEURYSM: INFLUENCE OF THE MATERIAL MODEL NEAR-INCOMPRESSIBILITY”.
- [6] C. J. Mills et al. “Pressure-flow relationships and vascular impedance in man”. en. In: Cardiovascular Re- search
- [7] Peter C. Luchsinger et al. “Instantaneous Pressure Dis- tribution Along the Human Aorta”.
- [8] José F. Rodríguez et al. “Mechanical Stresses in Abdominal Aortic Aneurysms: Influence of Diame- ter, Asymmetry, and Material Anisotropy”.
- [9] Mohammed Yahya. “Three dimensional finite-element modeling of blood flow in elastic vessels: effects of arterial geometry and elasticity on aneurysm growth and rupture”.
- [10] Benjamin Owen et al. “Structural modelling of the car- diovascular system”.
- [11] Comsol Multiphysics. “Nearly incompressible hyper- elastic materials”.
- [12] A Thapar et al. “Internal or external wall di- ameter for abdominal aortic aneurysm screening?” [13] Todd E. Rasmussen and John W. Hallett. “Inflam- matory Aortic Aneurysms: A Clinical Review with New Perspectives in Pathogenesis”.
- [14] M.L. Raghavan and David A. Vorp. “Toward a biome- chanical tool to evaluate rupture potential of abdom- inal aortic aneurysm: identification of a finite strain constitutive model and evaluation of its applicabil- ity”.
Appendice
Poiseuille
% load parameters
parameters.length=0.2;
parameters.diam=0.02;
parameters.R=parameters.diam/2;
parameters.viscosity=3.5*1e-3;
parameters.P0=12000;
parameters.PL=11998; %[Pa]
import com.comsol.model.* ;
import com.comsol.model.util.* % load comsol utility
model = ModelUtil.create(’Model’);
path=pwd; % load current directory path
addpath(strcat(path,’\subroutine\’)); % add subrotuine
model.modelPath(path); model.label(’poiseuille.mph’); model=setup_parameters(model,parameters);
model=create_geometry(model);
model=create_physics(model); model=solve_study(model);
[model,speed, radius]=post_process(model);
mphsave(model,’poiseuille_mat.mph’)
pouiseuille_analytic=@(r) (1/(4*parameters.viscosity))*((parameters.PL-parameters.P0)/parameters.length)*(r.^2-(parameters.R*ones(size(r))).^2); % analytic solution: r_span=-parameters.R:0.001:parameters.R;
figure; hold on;
scale=1e3;
plot([0,0],[0,max(max(speed))],’--k’,’LineWidth’,0.5) %plot central axis
ylabel(strcat(’Speed’,’ [m/s]’)) xlabel(strcat(’r’,’ [mm]’))
plot(scale*r_span,pouiseuille_analytic(r_span),’-’,’Color’, ’#A2142F’,’LineWidth’,2)
plot(scale*radius,speed,’o’,’Color’,’#0072BD’,’LineWidth’,1)
legend({’Cyl axis’,’Analytc’,’Numerical’})Post processing
Snippet della subroutine di postprocessing per estrarre i dati di pressione sul bordo interno del vaso (𝚋𝚘𝚞𝚗𝚍𝚊𝚛𝚢 4,5,6,7,8 e 9) e della pressione sulle linee di taglio. Il codice provvede anche a riordinare le ascisse curvilinee aggiornato con i rispettivi indici anche le ordinate.
wall_pressure = mpheval(model,’p’,’selection’,[4 5 6 7 8 9],’edim’,’boundary’);
% extract wall pressure on selected boundary % sort longitudinal value
[line,index]=sort(wall_pressure.p(2,:),’ascend’); pressure_wall=wall_pressure.d1(index);
figure;
plot(line,pressure_wall,’Color’,’#77AC30’,’LineWidth’,1.5) ylabel(strcat(’Pressure’,’ [’,wall_pressure.unit,’]’)); xlabel(’z [m]’); title(’Wall pressure’)
% extract pressure over cutline 1
[U, radius,unit]=mphinterp(model,{’spf.U’,’cln1x’},’dataset’,’cln1’);
clear index
[space,index]=sort(radius’,’ascend’); x_to_plot=space; pressure_cln1=U’;
for j=1:length(pressure_cln1(1,:))
pressure_cln1(:,j)=pressure_cln1(index(:,j),j);
end
% extract pressure over cutline 2
[U1, radius1,unit]=mphinterp(model,{’spf.U’,’cln2x’},’dataset’,’cln2’);
clear index; clear space
[space,index]=sort(radius1’,’ascend’); x_to_plot1=space; pressure_cln2=U1’;
for j=1:length(pressure_cln2(1,:)) pressure_cln2(:,j)=pressure_cln2(index(:,j),j);
end
% compute symmetrization
x_to_plot=[-fliplr(x_to_plot),x_to_plot]; pressure_cln1=[fliplr(pressure_cln1), pressure_cln1]; x_to_plot1=[-fliplr(x_to_plot1),x_to_plot1]; pressure_cln2=[fliplr(pressure_cln2),
pressure_cln2];
figure; scale=1e3; % scale to [mm]
plot(scale*x_to_plot ,pressure_cln1 ,’Color’,’#0072BD’,’LineWidth’,2,’DisplayName’,’cos(x)’);
hold on;
plot(scale*x_to_plot1 ,pressure_cln2 ,’Color’,’#A2142F’,’LineWidth’,2)
plot([0,0],[0,max(max(pressure_cln2))],’--k’)%plot central axis ylabel(strcat(’Speed’,’[’,unit(1),’]’)); xlabel(strcat(’r’,’ [mm]’)); legend({’z=0.5*L’,’’,’z=0.98*L’})Streamline 3D
Linee di flusso tridimensionali per l’analisi fluidodinamica con inlet e outlet stazionari.
Condizioni al bordo tempo dipendenti
Curva di pressione (a) e di velocità (b) utilizzate rispettivamente per l’outlet e l’inlet di pressione facenti riferimento ai parametri fisiologici misurati in Scotti et al. [2].
