Fluid structure interaction simulations of aortic valves : from a 2D towards a 3D case

The heart valve diseases are one of the most common cardiovascular disease (CVD) worldwide, the majority of which involves the aortic and mitral valves. Surgical treatments of aortic valve pathologies include the valve reparation or substitution of the damaged valve with a prosthetic valvular device. In this scenario, the computational modelling plays a fundamental role: heart valve (HV) simulations can provide a relevant insight on the physiological and pathological mechanics of valvular structures, and hemodynamics of the aortic root region. Furthermore, they allow the testing and the optimization process of new valvular prostheses and the virtual testing of new surgical techniques. Despite different approaches can be used to perform HV computational simulations, the most exhaustive and realistic approach is the fluid-structure interaction (FSI) simulation, which accounts for the interaction between the blood and the flexible valvular structures, resulting in a more reliable and complete scenario. This master thesis project wants to investigate the capabilities of Abaqus/CEL (Coupled Eulerian-Lagrangian), one of the FSI package of the commercial software Abaqus (Dassault system), when modelling heart valves. Abaqus/CEL is an extension of the structural solver Abaqus/Explicit, originally realized to solve problems of plasticity: therefore, FSI simulations represent an extreme case where the highly deformable material is a fluid. It results in an explicit, monolithic FSI algorithm, with the eulerian (i.e. IB) description of the fluid domain: the possibility to use a unique solver to perform the FSI analyses, represents a great advantage in computational time, either in the university and industry environment. To the best of our knowledge, this is the first time that this computational environment is used in the biomedical field for heart valve (HV) simulations. From a simplified 2D aortic valve towards a more complete and realistic 3D geometry, the parameters to set-up the FSI simulation have been investigated and the most suitable parameters chosen. Based on these analyses, a final 3D FSI simulation of an aortic valve has been performed and analysed. A 2D-simplified geometry of an aortic valve is used to perform a preliminary test case and the results are compared with those of an analogous model obtained with Tango, an in-house written coupling algorithm for an ALE-based, partitioned and strongly-coupled FSI simulation. The 2D geometry is composed of a straight rigid tube (eulerian domain) with two semi-circular enlargements to mimic the sinuses of Valsava and two flexible thin leaflets (lagrangian domain) placed immediately before the sinuses cavity. Shell elements (S4R) are chosen to discretize the leaflets, while three-dimensional eulerian elements (EC3D8R) are used for the fluid domain. Being three-dimensional elements, a thickness of one element layer is added to the model, resulting in a quasi-2D model. The leaflets are assumed to be linear elastic (Young’s modulus 1MPa, Poisson’s ratio 0.45) with a density of 1100 kg/m3, and rigid aortic walls are assumed via velocity boundary conditions. The blood is modelled as a Newtonian fluid (density 1060 kg/m3, viscosity 0.0035 Pa s), the introduction of a low compressible factor has been necessary to provide the convergence of the simulation. The contact between the leaflets is assumed to be frictionless and no penetration between the solid bodies is ensured by contact algorithm available is Abaqus/CEL. The fluid-solid interaction is managed via contact algorithm which ensures the two domains to not overlap throughout the analysis; moreover, no slip conditions at the fluid-structure interface are simulated. One cardiac cycle (0.8 s) has been simulated, applying the ventricular and aortic curves of pressure to the respective sections. To limit the possible effects introduced by the compressibility and to promote the development of the flow, a preconditioning phase of 0.2 s has been simulated. The simulations have been run on 8 CPUs (3.4 GHz of speed and 5.6 GB of memory). This model has been validated through the comparison with the analogous model realized in Tango. Although the fluid in Abaqus is modelled as slightly compressible, the fluid-dynamics showed to not be deeply influenced by the introduction of this compressibility factor. Despite few differences with the Tango model, due to intrinsic differences of the diverse modellisation, the fluid-dynamics and kinematics are overall comparable: similar velocity flow fields are obtained during the cardiac cycle. Furthermore, in both the analyses, characteristic biofluid-dynamic features are captured in the sinuses region during the closing phase such as: eddy currents. Driven by the difference of pressure between the ventricular and aortic sections, the leaflets achieve the closed configuration, but they are unable to bear the high pressures and, in both cases, they reverse to the ventricular side. For this reason only one cardiac cycle could be simulated. The comparison of the leaflets displacement is also considered satisfactory, although a small time delay affects the movement of the valve in the Abaqus simulation. With respect to the two software performances, the monolithic IB approach adopted by Abaqus is significantly less time consuming: the computational time was of 1h versus the 48h of running time needed for the Tango calculation. Using the quasi-2D geometry as a starting point, its extension to a 3D case is presented. All the set-up properties (eulerian geometry, element type, material model, BCs, loading conditions) have been investigated separately to choose the most suitable parameters for the final 3D HV simulations. A Carpentier-Edwards PERIMOUNT Aortic Heart Valve (Edwards Life-sciences LLC, Irvine, California, size 25 mm) was scanned with a μCT scan. The images were segmented with the commercial software Mimics (Materialise, Leuven, Belgium) to obtain the desired geometry. The reconstructed valve, with a thickness of 6 mm, is placed into a straight rigid tube with three hemispherical enlargements to mimic the sinuses of Valsalva. A cylindrical geometry is chosen for the eulerian domain. The valve is discretized with solid brick elements with reduced integration (C3D8R), while elements EC3D8R are adopted for the eulerian domain. The wall is considered a rigid body. In this 3D model, an hyperelastic material based on a bi-linear approximation of the physiological material properties, is used for the leaflets (density of 1100 kg/m3). The fluid material properties are chosen in accordance with those used in the 2D test case including the choice of the artificial compressibility factor. The contact properties, either solid-solid contact and fluid-structure interactions, are enforced via the Abaqus contact algorithm, as done in the 2D case. After a preconditioning phase of 0.2 s, to promote the flow development, one cardiac cycle (0.8 s) has been simulated by applying the gradient of pressure curve at the ventricular section. To verify the validity of the structural assumptions, a structural simulation on the same HV model, is also performed. All the structural properties are the same used in the FSI analysis. During the systolic phase, the pressure distribution is consistent with the time-dependent load applied at the ventricular section. Starting from a homogeneous distribution, the increasing magnitude in the ventricular portion drives the opening valve movement. As the valve opens, a symmetric flow develops aligned with the centerline of the AR region, achieving a maximum velocity of 1.3 m/s at the sinotubular junction, and once the fully open configuration is reached, eddy currents are observed behind the leaflets and around their tips. The obtained flow field and the eddy currents observed in the sinuses area are highly comparable with previous studies. The analysis of the valve opening movement reveals a symmetric displacement of the three cusps. From the comparison with the equivalent structural simulation, the maximum displacement obtained results lower and presents small time delay. The time delay is also evident if comparing the rapid valve opening time (RVOT) with other FSI simulations presented in literature. One of the reasons may lie in the compressibility introduced in this model: even if it could be considered negligible, the difference could be significant if compared with simulations where the blood is incompressible. From the mechanical point of view, high stresses are observed in the belly region and the maximum value of 28 kPa is measured at the attachment line and below the commissures, where the major flexural movements are experienced. On the contrary, the free edge experiences the lowest stresses throughout the systole. The stress pattern obtained is consistent with data available in literature. In this work, the systole was analysed. In fact, during diastole, some technical issues are introduced, which lead to the inconsistency of the latest portion of the simulation. In the described model, the negative pressure imposed at the inlet generates the suction of the fluid from the domain, and the tube becomes empty. If the boundary conditions are modified (i.e. two positive pressure curves, as done in the 2D case), the leaflets experience high deformation of the elements. A 2D testing case of aortic valve has been performed, and the comparison with an analogous analysis in Tango demonstrates the consistency of the results obtained in Abaqus, in terms of velocity flow field and displacement of the leaflets; it also reveals reduction of the computational time compared to Tango simulation. The 2D model has been expanded to a 3D model of aortic valve: one cardiac cycle has been simulated in Abaqus, accounting for the fluid-structure interaction. In this work only the systolic phase has been realized: both, the fluid-dynamics and mechanics of the leaflets are consistent with data presented in previous works. Unlike most of the published models, which treat the blood as incompressible or nearly compressible, a higher compressibility has been introduced in this model and may result in a time delay that affects the movement of the leaflets in both the 2D and 3D models. Despite the limitations and the need of further investigation on critical issues, this project shows that the FSI simulation of HV using Abaqus/CEL is a promising approach, the numerical results obtained are reasonable and comparable to literature data. Furthermore, the software is commercially available and widely used in both academic and industrial fields. If compared with partitioned approaches, Abaqus provides significant advantages in terms of computational costs. The 3D analysis has been restricted to the systolic phase of the cardiac cycle due to a technical issue, related to the volume of fluid definition, which appears when the curve of pressure becomes negative. The entire cardiac cycle could be simulated with the physiologic ventricular and aortic curves of pressure, as done in the 2D model. In this case a further study on the element type for the leaflets is necessary, in order to circumvent the distortion of elements which are unable to resist the high pressure loads. The influence of the compressibility factor introduced in the Abaqus models has been evaluated through the Mach number and the variation of fluid density. Although the results allow us to consider this influence negligible, a further investigation on this parameter will be performed in a future study. Once all the technical problems are solved, the model can easily be expanded to simulate a more realistic geometry: the distensibility of the aortic walls will introduce the deformation of the vessel which is typical of the physiological case. More complex scenarios can be evaluated by introducing the coronary vessels originating from the sinuses of Valsalva. The interwoven relationship between valvular stenosis and coronary flows can be then evaluated. From the valvular point of view, different devices can be tested by replacing the current prosthetic valve with the device of interest. As a long term result, this will possibly result in a critical comparison between the performances of different biological HV prostheses, allowing for a better evaluation of the device during the pre-operative phase in the clinical practise.

Le valvulopatie, ovvero malfunzionamenti delle valvole cardiache, sono tra le patologie cardiovascolari più diffuse e colpiscono più frequentemente le valvole del “cuore sinistro”. Nella maggior parte dei casi, la terapia prevista è chirurgica: la valvola danneggiata può esser riparata o sostituita con un dispositivo protesico. In ambito biomedico, la modellazione computazionale gioca un ruolo fondamentale: infatti le simulazioni computazionali di valvole cardiache possono fornire informazioni clinicamente rilevanti riguardo il comportamento meccanico delle strutture valvolari e l’emodinamica della regione aortica sia in condizioni fisiologiche che patologiche. Inoltre, consentono di testare virtualmente nuovi dispositivi prostetici, il processo di ottimizzazione, nonché permettono la simulazione e la valutazione di nuove tecniche operatorie. Nonostante la varietà di approcci che possano esser utilizzati per queste simulazioni, l’analisi numerica che consente di simulare l’interazione fluido-struttura (FSI) risulta essere quella più esaustiva e realistica, in grado di simulare uno scenario più completo e affidabile. L’obiettivo di questa tesi è indagare le capacità di Abaqus/CEL (Coupled Eulerian-Lagrangian), per la generazione di modelli di valvole cardiache. L’esistente solver strutturale Abaqus/Explicit, è stato esteso con Abaqus/CEL, includendo la possibilità di simulare problemi di plasticità con una trattazione euleriana del materiale; le simulazioni FSI rappresentano un caso estremo di questa analisi in cui il fluido è considerato un materiale altamente deformabile e viene descritto in modo euleriano. Abaqus/CEL utilizza un algoritmo esplicito per il calcolo della soluzione e risolve l’interazione fluido-struttura in modo monolitico, fatto che consente l’impiego di un unico solver per analisi FSI, limitando così problemi legati al tempo di calcolo tipici delle simulazioni FSI. Da un’analisi della letteratura disponibile, non è stato trovato nessun altro lavoro che utilizzi Abaqus/CEL per la soluzione FSI di problemi relativi alle valvole cardiache. In questa tesi sono stati analizzati parametri necessari per simulazioni di valvole aortiche a partire da una geometria semplificata bidimensionale fino ad un caso più complesso e realistico 3D e, sulla base di queste analisi, è stata eseguita una simulazione finale di interazione fluido-struttura di una valvola aortica 3D. I primi test svolti sono stati realizzati con un modello bidimensionale di valvola aortica, di cui è stato possibile verificare la correttezza grazie al confronto con un caso analogo realizzato con Tango, un algoritmo sviluppato presso l’Università di Ghent, department of heat flow and combustion mechanics, che consente l’interazione tra due solver separati e indipendenti (Fluent e Abaqus) e risolve il problema FSI iterativamente. La geometria 2D è composta da un tubo rigido con due allargamenti semicircolari a rappresentare i seni di Valsava e due leaflets flessibili, posti in corrispondenza delle cavità dei seni. Il dominio fluido è stato discretizzato con elementi euleriani tridimensionali (EC3D8R), mentre elementi shell (S4R) sono stati utilizzati per i foglietti valvolari. Poiché gli elementi euleriani sono tridimensionali, è stato necessario introdurre uno spessore artificiale di uno strato di elementi al modello, rendendo il modello quasi-2D. Il materiale dei leaflets è stato considerato lineare elastico (modulo di Young 1 MPa, coefficiente di Poisson 0.45) con una densità di 1100 kg/m3, mentre la parete aortica è stata modellata come un corpo rigido attraverso condizioni al contorno imposte sulle velocità del fluido in parete. Il sangue è stato considerato come un fluido Newtoniano in regime laminare (densità 1060 kg/m3, viscosità 0.0035 Pa s). Per motivi di convergenza della simulazione, è stato necessario introdurre fattore di comprimibilità nel fluido. Gli algoritmi di contatto disponibili in Abaqus/CEL hanno consentito di gestire sia il contatto tra i due leaflets che l’interazione fluido-struttura, evitando la compenetrazione tra i due corpi solidi e/o la sovrapposizione tra i due domini. In particolare, si è assunta l’assenza di attrito tra i leaflets e su tutte le pareti sono state imposte condizioni di non scivolamento (velocità in parete nulla). Prima di poter applicare il carico desiderato, è stata effettuata una prima fase di pressurizzazione di 0.2 s, necessaria per limitare i possibili effetti della comprimibilità del fluido e favorire lo sviluppo del flusso. La simulazione vera e propria ha previsto l’imposizione di un intero ciclo cardiaco (0.8 s), mediante le due curve fisiologiche di pressione aortica e ventricolare. Tutte le simulazioni sono state eseguite su 8 processori (3.4 GHz di velocità e 5.6 GB di memoria). Nonostante siano presenti alcune differenze tra i modelli realizzati in Abaqus e Tango, dovute alle caratteristiche intrinseche di ciascun approccio, i risultati si sono rivelati confrontabili. In primo luogo, la fluidodinamica in Abaqus non ha risentito eccessivamente del basso fattore di comprimibilità introdottovi. I profili di velocità ottenuti sono confrontabili e, in entrambe le analisi, è stata individuata la presenza di vortici nella regione dei seni, durante la fase di chiusura della valvola. In secondo luogo è stata valutata la cinematica valvolare: dal confronto, il movimento in chiusura dei leaflets, indotto dalla differenza di pressione tra le camere aortica e ventricolare, risulta analogo nelle due simulazioni, nonostante sia evidente un leggero ritardo nel modello in Abaqus. Infine, sono state confrontate le performance dei due software: l’approccio monolitico impiegato da Abaqus è risultato decisamente conveniente in termini di tempo di calcolo: il tempo richiesto da Abaqus, pari a 1h di calcolo, è nettamente inferiore se paragonato alle 48h necessarie per completare l’analisi in Tango. A partire dalla geometria quasi-2D, il modello è stato esteso fino a realizzare una geometria tridimensionale della regione aortica. Data le complessità del nuovo modello, i parametri fondamentali per la simulazione, quali geometria del dominio fluido, tipo di elementi, materiali, condizioni al contorno e condizioni di carico, sono stati analizzati separatamente al fine di facilitare l’identificazione dei parametri più adatti ed analizzare problemi che l’approccio di Abaqus comporta. I dettagli delle singole analisi sono riportati nella sezione 5.3. Combinando i risultati ottenuti dalle singole analisi, è stata eseguita una simulazione finale del modello tridimensionale, di cui verranno discussi i risultati nella seguente sezione. Il modello tridimensionale di regione aortica è costituito da un tubo cilindrico rigido con tre allargamenti emisferici a simulare i seni di Valsava. Una valvola aortica prostetica, spessa 6mm, è stata inserita nel vaso. La geometria della valvola è stata ottenuta a partire da una scansione µCT della valvola protesica biologica Carpentier-Edwards PERIMOUNT (Edwards Life-sciences LLC, Irvine, California, size 25 mm), ricostruita con il software commerciale Mimics (Materialise, Leuven, Belgium). Analogamente al modello 2D, elementi EC3D8R sono stati utilizzati per discretizzare il dominio euleriano, mentre elementi esaedrici C3D8R sono stati scelti per i leaflets, in quanto favoriscono la stabilità della simulazione. In questo modello 3D, per i foglietti valvolari (densità of 1100 kg/m3) è stato scelto un modello di materiale iperelastico, ottenuto dall’approssimazione bi-lineare delle proprietà meccaniche del tessuto fisiologico, mentre per il fluido le proprietà utilizzate sono conformi a quelle utilizzate nel test 2D, fattore di comprimibilità compreso. In egual modo, l’algoritmo di contatto fornito da Abaqus è stato utilizzato per gestire sia il contatto tra i due corpi solidi, sia l’interazione tra fluido e struttura. Dopo una prima fase di precarico (0.2 s), necessario per favorire lo sviluppo del flusso nella regione, la curva del gradiente di pressione fisiologica è stata applicata alla sezione ventricolare, con un ciclo cardiaco di 0.8 s. Infine, per verificare la validità delle ipotesi strutturali, è stata realizzata una simulazione strutturale della sola valvola con proprietà analoghe a quelle utilizzate nell’analisi FSI. Durante la fase di sistole, la distribuzione di pressione nella regione aortica è coerente con il carico applicato. A partire da una distribuzione pressoché omogenea, l’aumento della pressione nella porzione ventricolare determina il movimento di apertura della valvola. Tale movimento è accompagnato dalla formazione di un flusso ad alta velocità, centrale e simmetrico, che raggiunge una velocità massima pari a 1.3 m/s in corrispondenza della giunzione sinotubulare, quando la valvola è in conformazione completamente aperta. In questa conformazione, sono anche state identificate formazioni di vortici nella zona posteriore i leaflets. Tali aspetti biofluidodinamici (flusso, velocità massima del fluido e vorticosità) sono confrontabili con studi precedenti. Una seconda analisi condotta sulla cinematica di apertura della valvola ha evidenziato la simmetria nel movimento delle tre cuspidi; tuttavia, il confronto con la corrispondente analisi strutturale, ha rilevato alcune differenze: nella simulazione FSI, la valvola necessita di un tempo di apertura più elevato per raggiungere un EOA (effective orifice area) massimo inferiore. Il ritardo nell’apertura resta evidente se si confrontano i tempi di apertura rapida della valvola (RVOT) con altre simulazioni FSI disponibili in letteratura, verificando l’ipotesi che il ritardo non sia dovuto all’assenza di fluido (inerzia e resistenza al movimento) nella simulazione strutturale. Una delle cause di tale ritardo potrebbe risiedere nella comprimibilità introdotta artificialmente nel modello FSI in Abaqus che, seppur trascurabile, è significante se confrontata con simulazioni in cui il sangue viene modellato come fluido incomprimibile. Dal punto di vista meccanico, in fase di sistole, gli stress più elevati sono stati osservati in corrispondenza delle zone più flessibili: dapprima in corrispondenza della regione centrale di ciascun leaflet (belly region) e successivamente nell’area più periferica, tra le commissure e il margine di attacco alla parete aortica. In questa regione, nella fase di apertura, sono stati misurati gli stress massimi (28 kPa), mentre i valori più bassi sono stati ottenuti in prossimità di ciascun free edge, risultando in una distribuzione confrontabile con dati presenti in letteratura. In questo lavoro, si è analizzato il periodo di sistole poiché, in fase di diastole, compaiono alcuni problemi tecnici, rendendo i risultati non attendibili. Nel modello descritto, la pressione negativa applicata all’ingresso determina il risucchio del fluido dal dominio, svuotandolo. Al contrario, applicando le due curve fisiologiche di pressione, come nel caso 2D, gli elementi dei leaflets subiscono importanti deformazioni. Grazie al test iniziale con geometria bidimensionale e al suo confronto con una simulazione analoga in Tango è stato possibile verificare la coerenza e affidabilità dei risultati ottenuti in Abaqus/CEL, in termini di flusso di velocità e cinematica dei leaflets, e un elevato vantaggio in termini di costi computazionali. Successivamente il modello 2D è stato esteso a un modello tridimensionale di valvola aortica, sul quale è stata eseguita una analisi di interazione fluido-struttura di un ciclo cardiaco in Abaqus. I risultati fluidodinamici e meccanici, ottenuti in fase di sistole, sono confrontabili con dati presentati in studi precedenti. Tuttavia, a differenza dei modelli pubblicati, che trattano il sangue come fluido incomprimibile o leggermente comprimibile, in questo modello è stato introdotto un più alto fattore di comprimibilità che potrebbe essere parzialmente responsabile del ritardo rilevato nel movimento dei leaflets in entrambi i modelli, 2D e 3D. Nonostante siano stati identificati alcuni limiti e siano necessari ulteriori studi per superare i punti critici della formulazione FSI, questo progetto dimostra come Abaqus/CEL sia un approccio promettente per sviluppare simulazioni FSI di valvole cardiache: i risultati numerici ottenuti sono ragionevoli e confrontabili con dati di letteratura. In più, trattandosi di un software commerciale ampiamente impiegato nell’ambito sia industriale che accademico, il vantaggio in termini di costi computazionali è notevole, se confrontato con approcci di tipo partitioned tradizionalmente utilizzati in questo tipo di simulazioni. Come evidenziato, a causa di problemi tecnici legati alla definizione del volume di fluido, l’analisi del modello 3D è stata limitata alla sola fase di sistole. L’intero ciclo cardiaco potrebbe essere simulato applicando le due curve di pressione fisiologiche, ventricolare e aortica, alle rispettive sezioni, come è stato fatto nel modello 2D. In questo caso, i leaflets subiscono importanti deformazioni agli elementi; risulta necessaria dunque un’ulteriore analisi volta ad identificare il tipo di elemento in grado di resistere a elevati carichi pressori. Gli elementi studiati in questa tesi non hanno portato alla soluzione del problema. L’influenza del fattore di comprimibilità introdotto in questi modelli è stata valutata mediante il numero di Mach e le variazioni di densità del fluido. Nonostante i risultati consentano di considerare tale parametro trascurabile, un’ulteriore approfondimento dovrà essere eseguito in uno studio futuro. Quando tutti i problemi tecnici saranno stati risolti, il modello potrà essere facilmente modificato per ottenere casi più realistici che tengano in considerazione, per esempio, dell’elasticità delle pareti aortiche. Un’ulteriore estensione del modello prevede l’inclusione delle arterie coronarie che si originano dai seni di Valsava. Ciò consentirebbe la valutazione di scenari molto più complessi e patologici: si potranno infatti analizzare i fattori che legano stenosi valvolari al flusso nelle coronarie. Dal punto di vista valvolare, potranno esser testati diversi dispositivi protesici, sostituendo l’attuale protesi con il modello di interesse, adeguatamente scansionato e ricostruito per il modello. Questo approccio potrà, dunque, esser utilizzato per esaminare le prestazioni di differenti bioprotesi valvolari, consentendo una migliore valutazione clinica del dispositivo in fase pre-operatoria.