1 Introduction

Cardiovascular diseases (CVDs) are accounted for 45% of all deaths in Europe and for 37% of all deaths in the European Union [36, 51]. Among these casualties, ascending thoracic aortic aneurysm (aTAA) is the 19th common cause of human death [30].

Recently, the capacity to simulate in-silico patient-specific models, namely a digital twin on a computer, has pointed out the opportunity to give a structured, reproducible and predictive framework for interpreting and integrating clinical data, thus paving the way for developing personalized and preventive management strategies for CVDs. The medical digital twin is becoming a powerful numerical means to get such a result, mainly because of the important progress in techniques dealing with acquisition, reconstruction and visualization [12, 26], as well as the availability of always more efficient numerical algorithms [47]. Despite these enhancements, however, the direct application of the digital twin in a medical environment is still limited due to the high computational cost demanded to perform high fidelity simulations when computational fluid dynamics (CFD) techniques and tools are used.

Health-care provision can conceptually be simplified into four main processes: acquisition of clinical data, diagnosis, therapy planning and intervention. Ideally, numerical simulations ought to be carried out in a very fast way, to rapidly provide a quantitative output on each new patient-specific geometrical configuration and, therefore, a useful response from a practical point of view. With particular attention to aneurysmatic pathology, to overcome the aforementioned limitation while supplying responses in a time acceptable for clinical timing, lots of stochastic approaches have been proposed. Such techniques could account material properties and geometrical variations [21, 23, 25] as well as sets of uncertainties [15, 16] that can affect finite element (FE) simulations. Unfortunately, despite the advantages characterizing these approaches, none of them is able to provide, and show, the full field of the variables of interest in real-time, making their application unsatisfying in a clinical environment. The continuously growing availability of computational power and the simultaneous algorithmic improvements make possible nowadays the high-fidelity numerical resolution of complex problems via standard discretization procedures. However, these schemes remain prohibitevely expensive in many-query and real-time contexts, both in terms of CPU time and memory demand. This follows from the large amount of degrees of freedom (DOFs) they imply, resulting from the (fine) spatial discretization such as in case of hemodynamics studies. In light of this, recent advances in numerical tools assisting computational analyses, such as mesh morphing techniques [17] and reduced order models (ROM) [28, 33] can help in deepening the knowledge of complex phenomena pertaining clinical scenarios. The potential of the mesh morphing technique applied to cardiovascular field was introduced in [9] where Biancolini and co-workers implemented a pipeline based on radial basis functions (RBF) to investigate the effects of morphological alterations of the carotid bifurcation on the fluid-dynamic field. More recently, in Capellini et al. [19, 20] a new approach able to cope with the hemodynamic alterations in the ascending thoracic aorta during a bulge formation was proposed. On the other hand, reduced order modeling methods have received a significant attention in the last decades. ROM have been used, and successfully applied, in several fields of classical mechanics for systems with high complexity, in view of replacing the full model with a low DOF one so as to drastically reduce the overall cost of computing. Despite the ROM technique has been widely used in classical mechanics [2, 43], their implementation in the field of biomechanics is very recent and limited to few works [41, 44]. More recently Caiazzo et al [18] and Luboz et al. [40] explored the use of ROM for biomechanical applications. However, it is well known that fluid systems are difficult to reduce efficiently due to several non-linearities. Model reduction of the Navier-Stokes equations is a challenging task because their solutions tend to exhibit complex phenomena at multiple temporal and spatial scales, which means they are difficult to reduce to low-dimensional models without losing at least some of the scales [38].

In this paper we present a combined use of ROM and RBF technique for mesh morphing, with the aim to enable a very powerful and viable medical digital twin. A simplified CFD model of aortic aneurysms is presented by considering the hemodynamic effects that occur spanning a broad family of ascending aortic configurations, where RBF mesh morphing is introduced to fast and accurately control the geometrical parametrization of the vessel. The article is structured as follows: the next section describes the theoretical foundation of the approach (2). In Sect. 3 we present the digital twin data modelling, while in Sect. 4 we show and discuss the results gained through the proposed approach. Finally, we report the conclusions in Sect. 5.

2 Theoretical background

2.1 RBF

RBF are mathematical functions able to interpolate, on a distance basis, scalar information known only at discrete points, typically referred to as source points [12]. The quality and the behaviour of the interpolation depends on both the function and the kind of the chosen radial function. RBF can be roughly classified depending on the kind of support they guarantee (local or global), meaning the domain in which the radial function is not zero valued [31]. Some of the most common functions are gathered in Table 1.

Table 1 Common RBF with global and local support
Full size table

RBF can be defined in an n dimensional space, and are function of the distance which, in the case of morphing, can be assumed as the Euclidean norm of the distance between two points in the space. A linear system of order equal to the number of points used [17] shall be solved in order to evaluate the system sought coefficients. Once such coefficients are found, the displacement of a given node of the mesh, either inside (interpolation) or outside the morphing domain (extrapolation), can be calculated as the superposition of the radial contribution of each source point. In this way, it is then possible to define at known points the displacement in the space and, thus, to retrieve the value at mesh nodes obtaining a mesh deformation that leaves the grid topology unaltered [7]. The interpolation function is composed by the basis \(\phi \) and the polynomial term h with a degree that depends on the kind of the chosen basis. This latter contribution, in particular, is added to assure uniqueness of the problem and polynomial precision, allowing to recover exactly rigid body motions. If N is the total number of source points \(\mathbf{x} _{{k}_i}\), the displacement \(\textit{s}\) at a given point in the space \({{\varvec{x}}}\) can be determined according to the following relation:

$$\begin{aligned} s(\mathbf{x} )=\sum _{i=1}^N \gamma _i \phi ( \Vert \mathbf{x} - \mathbf{x} _{{k}_i} \Vert ) + h(\mathbf{x} ) \end{aligned}$$
(1)

An interpolation exists if coefficients \(\gamma _i\) and weights of the linear polynomial \(\beta _i\) can be found such that the given value at source points \(\mathbf{x} _{{k}_i}\) can be retrieved exactly. At source points the function passage conditions are imposed together with the orthogonality conditions:

$$\begin{aligned} s(\mathbf{x} _{{k}_i}) = g_i, 1 \le i \le N \qquad \mathrm {and} \qquad \sum _{i=1}^N \gamma _i p(\mathbf{x} _{{k}_i}) = 0 \end{aligned}$$
(2)

for all the polynomials p of degree less or equal to polynomial h being \(g_i\) the assigned displacement vector at source point \(\mathbf{x} _{{k}_i}\). A single interpolant exists if the basis is conditionally positive definite [42]. If the degree is \(m \le 2\) [7] a linear polynomial can be used:

$$\begin{aligned} h(\mathbf{x} ) = \beta _1 + \beta _2 x_1 + \beta _3 x_2 + \cdots + \beta _{n+1} x_n \qquad \text {in} \qquad \mathrm{I\!R}^{\text {n}} \end{aligned}$$
(3)

The system built to calculate coefficients and weights can be written in matrix form for an easy implementation:

$$\begin{aligned} \left[ \begin{array}{cc} \mathbf{M} &{} \mathbf{P} \\ \mathbf{P} ^T &{} \mathbf{0} \end{array}\right] \left\{ \begin{array}{c} \varvec{\gamma }\\ \varvec{\beta } \end{array} \right\} = \left\{ \begin{array}{c} \mathbf{g} \\ \mathbf{0} \end{array} \right\} \end{aligned}$$
(4)

where g is is the vector of known displacements for each source point and M is the interpolation matrix with the radial distances between source points:

$$\begin{aligned} M_{ij} = \phi (\Vert \mathbf{x} _{{k}_i} - \mathbf{x} _{{k}_j} \Vert ), 1 \le i \le N, 1 \le j \le N \end{aligned}$$
(5)

P is the constraint matrix containing the coordinates of source points in the space:

$$\begin{aligned} \mathbf{P} = \begin{pmatrix} 1 &{} x_{{k}_1} &{} y_{{k}_1} &{} z_{{k}_1} \\ 1 &{} x_{{k}_2} &{} y_{{k}_2} &{} z_{{k}_2} \\ \vdots &{} \vdots &{} \vdots \\ 1 &{} x_{{k}_N} &{} y_{{k}_N} &{} z_{{k}_N} \\ \end{pmatrix} \end{aligned}$$
(6)

Once the weights and and coefficients of the system are obtained, the displacement values can be retrieved for the three directions at a given x point as [11] :

$$\begin{aligned} {\left\{ \begin{array}{ll} S_x(\mathbf{x} ) = \sum _{i=1}^N \gamma _i^x \phi ( \Vert \mathbf{x} - \mathbf{x} _{{k}_i} \Vert ) + h_x(\mathbf{x} )\\ S_y(\mathbf{x} ) = \sum _{i=1}^N \gamma _i^y \phi ( \Vert \mathbf{x} - \mathbf{x} _{{k}_i} \Vert ) + h_y(\mathbf{x} )\\ S_z(\mathbf{x} ) = \sum _{i=1}^N \gamma _i^z \phi ( \Vert \mathbf{x} - \mathbf{x} _{{k}_i} \Vert ) + h_z(\mathbf{x} ) \end{array}\right. } \end{aligned}$$
(7)

In computer-aided engineering (CAE) applications, RBF can be used to drive the morphing of the computational model discretized domain nodes \(\mathbf{x} _i\) by applying predefined displacement fields (\({g}_i\)) to a set of purposely generated points (\(\mathbf{x} _{{k}_i}\)), called RBF source points or centres, and, then, interpolating them in order to update the nodal positions of the mesh (target points)[11] adopting a specific radial function \(\phi \). The main characteristics of such an approach are the meshless property, the preservation of mesh consistency and the low disk usage. Main advantages include the exact control of nodes during smoothing, the prevention of remeshing noise, the possibility to be fully integrated in the computing process and the high performance in handling large models [27].

Notable examples of RBF applications are fluid-dynamic and structural shape optimizations [14, 35], fluid structure interaction [3, 8], and motorsport [45], nautical [13, 50], road infrastructures [48] and medical [9] analyses.

Fig. 1
figure 1

ROM workflow implemented in ANSYS Release 19.1 with indication of the main blocks and the estimation of the time for each step

Full size image

2.2 ROM

Reduced-order models for fluidynamic applications operate on parametrically generated data, the so called snapshots, represented by either surface quantities (e.g. surface pressure and shear stress) or volume quantities. ROM have been applied during years in several technical disciplines including multiphysics [34] and medical [6, 28] to handle applications concerning design [52], optimization [37] and control [5]. In this paper the solution recently implemented in the ANSYS software is described. Adopting the terminology used in the ANSYS training material and working context, the typical ROM workflow is divided in two main blocks:

  1. 1.

    ROM Creation including the use of ROM Builder and 3D solvers;

  2. 2.

    ROM Consumption including the use of 3D solvers and Stand Alone Systems.

While the first step can require hours or days to be set up, the second step demands a minimal intervention of the user to provide solutions almost in real-time. Figure 1 depicts through a block diagram the ROM workflow in the ANSYS framework environment. The key concept of the ROM environment is its integration with the response surface (RS) technique. By using RS, the values of output parameters are obtained as a function of input parameters on the design points defined in the setup of the design of experiments. Finally, interpolation methods are used for the building up of the entire model. In a fluid-dynamic context, the ROM is applied on the outputs of the CFD solver ANSYS\(^{\textregistered }\) Fluent\(^{\textregistered }\)- Release 19.1 (Fluent) which computes a set of calculations, called “snapshots”, that are representative of the parameter space. This set of frames is successively compressed into a small number of modes using a proper orthogonal decomposition (POD) [29]. This technique allows each solution frame to be linearly described through a small set of scalar coefficients exploiting the RS technique which, in turn, is based on the Kriging interpolation available in ANSYS\(^{\textregistered }\) DesignXplorer\(^{\mathrm{TM}}\)-Release 19.1 (DesignXplorer). Once built, the ROM interactively gives an accurate approximation of the solver solution for a new set of input parameters applied to the model nearly in real-time. The accuracy of ROM results depends on the number and quality of snapshots, whose suitability makes the ROM able to entirely represent the physical behaviour encountered in the parameter space, just exploiting the number of modes extracted using the POD algorithm.

3 Materials and methods

3.1 Geometry

The geometrical models of both healthy and pathological patients were obtained according to what was previously described in Capellini et al. [19][20]. Briefly, patient image data were retrospectively segmented and analyzed with VMTK (www.vmtk.org) which was integrated with custom Python scripts. Firstly, for each patient the 3D surface model was reconstructed and, then, the global thoracic aorta centerline \(\varsigma \) was computed and automatically subdivided to extract the ascending tract [4]. Finally, the cross sections along \(\varsigma \) at regular intervals of 1 mm were traced, and the area was computed for each cross section. The minimum diameter (\(D_{min}\)) and the maximum diameter (\(D_{max}\)) were computed for the section with the maximum area, while the tortuosity values have been calculated according to Piccinelli et al. [46]. A statistically significant shape geometry for both healthy (\(SS_H\)) and aneurysmatic (\(SS_A\)) subjects was thus defined. Due to the interest in the ascending region, the same descending aorta of the \(SS_H\) was also used for the \(SS_A\) model. The final 3D shape geometries of healthy and pathological patients are reported in Fig. 2.

Fig. 2
figure 2

\(SS_H\) (a) and \(SS_A\) (b) geometries with indication of the \(\varsigma \) of the ascending aorta (point-dot line) and of the starting and final points of the \(\varsigma \) (\(SP_c\) and \(FP_c\), respectively). In b the cross-section of maximum diameter of the aTAA is indicated (\(D_{max}\)) and the bulge extension reported along the \(\varsigma \) from the start bulge cross-section (\(SB_e\)) to the final bulge cross-section (\(FB_e\))

Full size image

The extension of the bulge was defined as the portion of \(\varsigma \) whose points correspond to a diameter of the bulge above the 95\(\%\) of the maximum diameter. For what just defined, the \(\varsigma \) of the ascending aorta starts at \(SP_c\), that is the sinu-tubular junction, and finishes at \(FP_c\) that is located in the end of the ascending portion of aorta (see Fig. 2). The variation of the parameters has been assumed considering the result of the morphological analysis presented in Capellini et al. [20] where a total of 45 geometries of pathological patients selected from the most representative have been processed. In particular we focused on the statistical analysis performed for the morphological parameters such as the \(D_{max}\), the tortuosity of the ascending aorta and the bulge extension. The mean value of \(D_{max}\) resulted to be 50.82 mm with a standard deviation of 5.13 mm, the tortuosity mean value was equal to 0.113 with a standard deviation of 0.0354 and the bulge extension showed an average value of 29.81 mm with a standard deviation of 11.88 mm.

3.2 RBF set-up

The RBF morphing procedures concerning the set-ups of the RBF solutions enabling the wanted shape variations, have been accomplished with the RBF Morph software (www.rbf-morph.com). The mesh morphing has been applied through a set of source points of the ascending thoracic aortic region. The localization of these generated points is reported in Fig. 3a,b from two perspectives.

Fig. 3
figure 3

Position of the source points generated for the aorta model and two cylinder-shaped delimiting volumes (a). Nodes of the inferior boundary source-points region (\(ISP_{RBF}\)); nodes of the maximum diameter boundary (\(DSP_{RBF}\)) and nodes of the superior boundary source-points region (\(SSP_{RBF}\)) (b). Local coordinate reference system (c) and five different morphing effects (d–h)

Full size image

In particular, Fig. 3a depicts two cylinders assigned to define the volume delimiting the effects due to the application of mesh modifications: such a volume corresponds to the boolean union of the two cylinders. Three sets of RBF source points have been suitably created at:

  1. 1.

    the inferior boundary (\(ISP_{RBF}\));

  2. 2.

    the boundary where the maximum diameter is expected (\(DSP_{RBF}\));

  3. 3.

    the superior boundary (\(SSP_{RBF}\)).

The actual modifications of the geometrical model, qualitatively and qualitatively expressed in terms of maximum diameter, extension and tortuosity, have been enabled by applying two scaling transformations and three rigid translations. A total of five RBF shapes have been generated according to the predefined field of the geometries. The set-ups of such variations share the same rationale and differ only in the type of movement assigned to a specific subset of the generated RBF source points. In more detail, during the morphing process the \(SSP_{RBF}\) and \(ISP_{RBF}\) of the ascending aorta were kept fixed, while the \(DSP_{RBF}\) were used to control the type of the morphing action and the corresponding extent. As just introduced, the RBF solutions are characterized by a specific rigid movement or scaling transformations imposed to the \(DSP_{RBF}\). Figure 3c depicts the local coordinate reference system (CRS) specified for accomplishing the set-up of the RBF mesh morphing solutions. The origin of CRS is the point c identified by the intersection between the area over which \(DSP_{RBF}\) lay and the centerline \(\varsigma \) of the aorta. The x-axis is aligned with the direction toward which the maximum bulge occurs (anterior area of the aorta), while the versor of the z-axis is tangent to centerline \(\varsigma \) in c. The procedure and the final effect caused the shape modifications is summarized as follows:

  1. 1.

    scaling in x+ direction: the bulge maximum diameter is controlled applying a scaling of an extent of 20\(\%\) with respect to an axis parallel to the z axis and positioned in the x- axis (close to the posterior aortic surface). With this set-up, the vessel is mostly inflated towards the x+ direction (Fig. 3.d). This shape parameter is called \(\text {BW}_1\) (bulge-width);

  2. 2.

    scaling in y direction: the bulge maximum diameter is controlled applying a scaling of an extent of 20\(\%\) with respect to an axis parallel to the z axis and positioned near the center c of the CRS. With this set-up, the vessel is inflated towards the y+/- direction (Fig. 3e). This shape parameter is called \(\text {BW}_2\);

  3. 3.

    rigid motion in x+ direction: the \(DSP_{RBF}\) are translated in the x+ axis direction. With this set-up, an increase of the vessel projected tortuosity in the x-z plane is implemented (Fig. 3f). This shape parameter is called \(\text {BO}_1\) (bulge-offset);

  4. 4.

    rigid motion in y+ direction: the \(DSP_{RBF}\) are translated in the y+ axis direction. With this set-up, an increase of the vessel projected tortuosity in the y-z plane is implemented (Fig. 3g). This shape parameter is called \(\text {BO}_2\);

  5. 5.

    rigid motion in z+ direction: the \(DSP_{RBF}\) are translated in the z+ axis direction. This set-up allows the definition of a new maximum diameter location along the centerline \(\varsigma \) (Fig. 3h). This shape parameter is called \(\text {BO}_3\).

To appreciate the induced variations, each image shall be compared to the one dealing with the baseline configuration (Fig. 3a,b). Relating to the amplification values that are employed to enable morphing, the first three RBF solutions range between 0 and 3, whilst the last two between -1 and +1. Such values have been set taking into account the bounds dictated by the real modifications, so that the combinations of the shape modifications are such to cover the real geometrical variations.

Fig. 4
figure 4

Example of RBF Morph GUI for interactive shape update. In red the effect of the morphing procedure on the ascending aorta surface. (color figure online)

Full size image

3.3 Mesh and CFD set-up

The CFD field in steady state conditions was determined by means of the finite-volume commercial solver Fluent. The mesh model was generated according to what was already reported in [20]. In particular, a measuring extension equal to three diameters was added to both the ascending aortic inlet and all outlets in a perpendicular direction to the vessel cross-section. An unstructured mesh of about 1.5 million tetrahedral elements was generated in order to discretize the fluid volume. Moreover, 6 inflation layers, with a growth-rate of 1.2 and a total thickness of 1.8 mm were created. The blood flow was assumed to behave as incompressible and Newtonian in laminar regime, and to have a density of 1060 kg/m3 and a kinematic viscosity of 3.5 \(\times \) 10 -3 Pa s. The assumption of the Newtonian fluid for blood with a constant viscosity is acceptable in large vessels [1]. The boundary conditions were defined starting from velocity inlet and outlet pressure profiles described in the previous work [20]. Specifically, the velocity value at systolic peak time (1.115 m/s) was prescribed as velocity inlet condition at ascending aortic inlet, while a constant pressure was imposed at all outlets considering the pressure values at systolic peak time previously determined for each of them (BCA and LSA 14,398 Pa, LCCA 14,665 Pa, descending aorta 13,865 Pa). The CFD simulations were performed on a Workstation equipped with 2 Intel Xeon 2.80 GHz processors with 10 cores each, and 128 GB of shared RAM using 20 parallel processes. The total time required for the computation of each model would have been 8 h approximately. However, since RBF morphing maintains the mesh topology unchanged, the CFD solution for each geometrical variation was strongly accelerated by initializing the flow field using the converged solution obtained for the baseline configuration. By employing this strategy, each design point calculation took about 30 min only.

3.4 ROM set-up

As introduced in Sect. 2.2, the ROM implementation provided in Fluent requires that several 3D simulations are performed in view of providing the ROM algorithm with a proper set of snapshots. In this work, to take into account the huge variation in terms of geometry and flow that are undeniably enabled by the 5 shape variations, 40 snapshots were evaluated employing the Kriging method. In such a manner, the design points in the variable space were suitably seeded. After that, the POD was calculated extracting the first 5 orthogonal modes from the gained snapshots. It is worth to mention that up to 10 modes were taken into account in the preliminary stage of the study, but only 5 were finally retained since the maximum local difference in terms of pressure was evaluated to be below 1%.

4 Results

The results presented in this section are not intended to get an exhaustive knowledge with respect to aTAA numerical modelling but, instead, to assess the real capabilities of the inspection tool conceived for handling a medical scenario generated by a workflow making use of RBF mesh morphing to parameterize shapes, and accelerated by ROM. It is worth to underline that the same approach could be quite straightforwardly replicated for similar in silico simulation scenarios where the digital twin, based on CFD simulations, could be a key enabler.

4.1 RBF verification

The ability of getting a parametric shape has been first verified by exploring at which extent the five shape parameters can be combined keeping the quality of the CFD mesh at an acceptable level. Table 2 shows that, despite the large extent of shape modifications explored, the mesh morphing action allows to preserve a valid mesh with a reasonable increment of the skewness.

Table 2 Skewness values for four combinations of transformations by modifying the two bulge-width (BW) and three bulge-offset (BO) values shown in Sect. 3.2
Full size table
Fig. 5
figure 5

Pressure contours comparison between ROM (a–c) and CFD (d–f) evaluation for three different morphologies (a–d), (b–e) and (c–f)

Full size image

The parametric shape can be inspected in an interactive session of Fluent or used in batch by binding the shape amplifications to ANSYS® WorkbenchTM parameters. The interactive usage can be tweaked acting on the Multi-Sol panel of RBF Morph that is shown in Fig. 4. As far as computing performances are concerned, the time required to morph the full CFD mesh is 43s on a i7 laptop with 16 GB RAM, while the time required to prepare a preview for surface mesh of the whole model takes less than 1s. Further shape previews happen almost in real time.

Fig. 6
figure 6

Maximum errors between ROM and CFD for the pressures on the surface (a,b) and in the volume (c,d)

Full size image

4.2 ROM verification

In order to verify the ROM strategy in this new hemodynamics context, reduced and non-reduced approaches, namely the full differential equations system describing the CFD case, were compared. Such a comparison was accomplished using a single geometrical parameter that gradually morphs the healthy aorta into an aneurysmatic one, and monitoring the difference in the main quantities of interest during morphing. In Fig. 5 ROM and CFD simulations were compared at three stages of accretion ranging from healthy, (a) for ROM and (d) for CFD, to fully developed aneurysm, (c) for ROM and (f) for CFD. Images (b) and (e) refer respectively to ROM and CFD evaluations for the medium developed shape.

To provide a fair term of comparison between ROM and CFD simulations, the same RBF shape parameter was employed to morph the numerical grid of both cases. To further guarantee a proper simulation accuracy, the maximum skewness and aspect ratio were constantly monitored during calculation. Pressure maps compare satisfactorily, with slight differences between distribution as depicted in Fig. 5b,e.

Fig. 7
figure 7

Example of CFD-ROM results for different morphologies obtained with RBF mesh morphing techniques considering the scaling effect of two parameters (x/y) at different ratios: 0/0 (a); 0/2 (b); 3/0 (c) and 3/3 (d). For x and y direction see Fig. 3c

Full size image

To quantify these differences, the maximum error between CFD and ROM pressures was calculated for this latter configuration, namely the worst one. The percentage error distributions on aorta surfaces and in the aorta volume are shown in Fig. 6a,b and c,d, respectively. The former presents a maximum value of 1.5%, while the latter is characterized by the highest discrepancy between ROM and CFD results that is around 2.5%. This value, that can be improved by increasing the number of modes and snapshots employed in the ROM, sounds reasonable considering that the ROM result can be visualised almost in real-time whilst the CFD computation completion requires at least 30 min. In Fig. 6a,b discrepancies over the aorta surfaces show a maximum error of 1.5%. The highest discrepancy is around 2.5% and it is found in the volume as shown in Fig. 6c,d. In Fig. 7 four different morphologies are reported as example of pressure maps.

4.3 ROM interaction

Having successfully verified the ROM accuracy employing a single geometrical variation, the workflow was also tested in a more challenging scenario including the full set of parameters described in Sect. 3.2. In this study we have experimented how to consume the ROM model generated by the ROM builder of ANSYS. It is worth to comment that all ROM data are distilled into a single file with extension .romz intended for an offline consumption that could be a custom application. We have explored two possible options offered by ANSYS. The first one is based on a specific post-processing software called ROM viewer developed by ANSYS, whilst the second one foresees to use the ROM results to compute a new full computed CFD model (solution and results) that can be reintroduced in a standard CAE workflow. These solutions are respectively deepened in the following sections.

4.3.1 ROM consumption with ROM viewer

By using the standalone ROM viewer software, it is possible to consume a previously calculated .romz file. Since the .romz file is a standalone digital object, it includes not only the ROM of the fluid system, but also the numerical grid thus enabling a standalone visualization of the ROM results. In the right side of Fig. 8 the aTAA case loaded in the ROM viewer interface is shown. As visible, the five sliders positioned in the top left part of the GUI allow the user to tweak interactively the amplification values of the shape parameters.

Fig. 8
figure 8

Example of ROM viewer launcher interface with aTAA case (a) and ROM result for verifying the geometry from Fig. 5b

Full size image

The 1.5 GB .romz file is loaded in 30 min but, once the ROM information is in memory, each result can be evaluated in almost real-time by interactively moving the sliders. In the right side of Fig. 8 the ROM post-processing using the same parametrization employed in Fig. 5 is shown employing the ROM viewer.

4.3.2 ROM consumption with Fluent

As previously introduced, another viable option to post-process the ROM file is importing the .romz file in Fluent. Following this strategy the ROM result can be computed directly into Fluent by changing the parameters value from DesignXplorer. In this way it is possible not only to accurately post-process the updated flow, but also to insert the ROM evaluation in a classical CAE workflow by simply substituting the original CFD simulation. The time required to load and make available ROM results depends on the mesh dimension and is equal to the one necessary for loading the computational grid in the Fluent environment.

5 Discussions and conclusions

The basic concept of a digital twin is that a computer model can predict specific functioning of an individual’s body. It is well recognized that patient-specific computational models can provide additional information extrapolating data not available at clinical level. In this context, CFD simulations can aid the identification of patients prone to adverse outcomes by providing detailed information about haemodynamic factors. Moreover, numerical models may support clinicians by virtually simulating different interventional strategies. Despite all these advantages, numerical simulations are not currently applied in clinical field for the decision-making process around a disease. One of the main limitations regarding this transferability of knowledge and tools from the engineering scenario to clinical one, is the computational time required for the simulations that do not match with the health-care system. In literature several studies have tried to answer to this need by providing an overall solution in terms of probability and sensitivity simulations [22, 24, 32, 39]. However, even if these studies defined the formulation of a “virtual-patient library”, they were unable to provide the results of the entire field of investigation and, consequently, they are no directly applicable for a clinical environments. In this study we have demonstrated how the coupling of ROM approach and RBF mesh morphing allows to explore full field CFD results in an interactive and almost real time way. Geometrical models of the thoracic aorta district for an averaged healthy and a pathological patient are here adopted as the reference of patient specific data and a CFD model of the healthy one was considered as the starting one. Mesh morphing actions can bring this model onto the pathological one, as previously presented in [19, 20], or tweaked to represent a new specific pathological shape according to five parameters defined to control the bulge position and shape. The effectiveness of the new shape parameterization was verified showing how the mesh morphing of the surface can blend parameters to get reasonable shapes and that a valid volume mesh of the CFD domain is achieved even for the most aggressive combination of shape parameters. After a verification of the ROM approach, achieved by demonstrating how the morphing from healthy to pathological shape can be predicted using ROM, instead of a full CFD run, getting negligible errors, we have explored the ROM consumption according to two approaches. The first one relies on the ROM viewer tool of ANSYS which is capable to read the .romz file and work in a standalone fashion, whilst the second approach foresees to read into Fluent a new result that becomes available for standard post-processing as the ones computed by the full CFD solution. It is worth to comment that the proposed ROM approach, once integrated with proper custom visualization tools, could be coupled to an haptic device extending what demonstrated in [10, 49] to inspect the shape effect on the flow in real time by interactively manipulating the digital twin representation. The interest of the approach presented in this work goes beyond the preliminary application here addressed. Indeed, reduced model can be applied to a broad variety of problems for which real-time simulation in complex parameterized geometries is the key issue. On a medical perspective, we can say that morphing of subject-specific models of aorta segments from CT data is a stepping stone to explore the definition of shape degeneration during bulge formation on a single subject as well as on a complete population. Moreover, it could be a promising tool to: (i) fast re-mesh when conducting sensitivity studies (e.g. on device design or positioning) (ii) easily compare results sets from two or more meshes (since morphing generates isotopological meshes, that share the same node numbering and connectivity) (iii) improve the speed and automation of subject-specific FE model generation while keeping surface regularity. The soundness of the ROM concept for the aTAA interactive modelling has been here explored for a simplified flow condition, i.e. steady state flow at systolic peak conditions, but could be consistently extended to full transient CFD simulations to account for pulsatility of flow.