Content area
This paper presents an open-source pipeline for simulating flow and flow-related processes in (embedded) tubular structures. Addressing a gap in computational fluid dynamics (CFD) and simulation sciences, it facilitates the transition from raw three-dimensional imaging, graph networks or computer aided design (CAD) models of tubular objects to refined, simulation-ready meshes. This transition, traditionally labourintensive, is streamlined through a series of innovative steps that include surface mesh processing, centre-line construction, anisotropic mesh generation and volumetric meshing, leading to finite element method (FEM) simulations. The pipeline leverages a range of open-source software and libraries, notably GIBBON, FEniCS and Paraview, to provide flexibility and broad applicability across different simulation scenarios, ranging from biomedical to industrial applications. We demonstrate the versatility of our approach through five applications, including the mesh generation for soil-root systems, lung airways, microcirculation networks and portal vein networks, each originating from a different data source. Moreover, for several of these cases, we incorporate CFD simulations and strategies for 3D-1D coupling between the embedding domain and the embedded structures. Finally, we outline some future perspectives aimed at enhancing accuracy, reducing computational time and incorporating advanced modelling and boundary condition strategies to further refine the framework's capabilities.
Keywords:
computational fluid dynamics, 3D-1D mixed-dimension, embedded networks, vascular networks, conforming mesh, lung airways
This paper presents an open-source pipeline for simulating flow and flow-related processes in (embedded) tubular structures. Addressing a gap in computational fluid dynamics (CFD) and simulation sciences, it facilitates the transition from raw three-dimensional imaging, graph networks or computer aided design (CAD) models of tubular objects to refined, simulation-ready meshes. This transition, traditionally labourintensive, is streamlined through a series of innovative steps that include surface mesh processing, centre-line construction, anisotropic mesh generation and volumetric meshing, leading to finite element method (FEM) simulations. The pipeline leverages a range of open-source software and libraries, notably GIBBON, FEniCS and Paraview, to provide flexibility and broad applicability across different simulation scenarios, ranging from biomedical to industrial applications. We demonstrate the versatility of our approach through five applications, including the mesh generation for soil-root systems, lung airways, microcirculation networks and portal vein networks, each originating from a different data source. Moreover, for several of these cases, we incorporate CFD simulations and strategies for 3D-1D coupling between the embedding domain and the embedded structures. Finally, we outline some future perspectives aimed at enhancing accuracy, reducing computational time and incorporating advanced modelling and boundary condition strategies to further refine the framework's capabilities.
(ProQuest: ... denotes formulae omitted.)
1. Introduction
In the rapidly evolving field of computational fluid dynamics (CFD) and simulation sciences, the precise and accurate representation of tubular structures, including those embedded within complex environments, emerges as a crucial element across a broad spectrum of applications. These applications span diverse fields, from detailed biophysical interactions within root-soil systems, as evidenced by previous studies [1-6], to the exploration and biomimicry of ant nest architectures and their tunnelling strategies [7,8]. They extend further to the biomedical realm from the modelling of flows in blood vessels [9-11] and within lung airways [12,13] to understanding tissue perfusion dynamics [14-20]. The relevance of (embedded) tubular structures also extends to industrial applications, for example the design and optimization of heat exchangers [21], the efficiency of oil extraction processes at the scale of wells and reservoirs [22,23] and safety considerations in tunnel ventilation systems [24] through to the engineering complexity involved in ensuring the stability and integrity of pipeline-soil systems [25,26].
Recent developments in imaging techniques (such as synchrotron or laboratory X-ray computed tomography (CT) [27-29], neutron CT [30] and multifluorescence high-resolution episcopic microscopy [31]), coupled with novel segmentation and image processing methods [32,33], have provided unprecedented access to realistic structures comprising tubular objects down to the finest length scales. Alternatively, when the required imaging resolution is lacking, it is possible to algorithmically generate biomimetic and morphologically realistic synthetic tubular objects algorithmically [34-36]. This availability of structural data from imaging or synthetic sources presents both challenges and opportunities. Combining such data with mathematical models of physical processes provides the opportunity to interrogate the link between structural and functional relationships [37]. However, this relies on robust tools to extract the structural features from the raw imaging data [33].
Despite this opportunity, the pathway from obtaining three-dimensional imaging, synthetic data or computer-aided design (CAD) models of (embedded) tubular objects to generating meshes that are ready for numerical simulations presents significant hurdles. This process is often labour-intensive and fraught with challenges, particularly when it requires the transition from raw data to refined, simulation-ready formats. Additionally, when these models are combined with FEM to perform relevant simulations such as CFD and advection-diffusion-reaction simulations, followed by the need for careful post-processing, the whole process becomes much more complex. This issue becomes even more challenging when trying to complete all these steps within a single open-source platform. In fact, the absence of seamless integration between various stages can make workflows more complex, particularly because of the format incompatibilities among different tools. Moreover, implementing an open-source platform often comes at the cost of using different libraries/software and ensuring input/output compatibility. However, opensource platforms increase the transparency and reproducibility of the research and, thereby, accelerate innovation within the scientific community.
Consequently, it is crucial to develop a comprehensive highly automated open-source pipeline that integrates the various steps of the process into a single framework. This unified approach would streamline workflows and promote the use of advanced simulation techniques across a wide range of applications. Such a workflow would mitigate the labour-intensive work associated with processing geometries and enabling the automated analysis of large datasets. This is particularly pivotal in fields such as the biomedical sector, where the capability to efficiently handle pre-clinical data (e.g. from animal models) and/or patient-specific data can significantly enhance diagnostic and therapeutic strategies on large multi-object datasets.
This paper is structured as follows: §2 describes our methods, starting with a comprehensive review of the software choices that underpin our pipeline and highlighting the interoperability of software packages and the automation strategies employed to streamline these processes. The different compartments are then described including the surface mesh processing from network-based and image-based inputs, the skeletonization of tubular structures, the boundary condition automation, the mesh generation and conversion, the FEM for the mathematical models integrated in the framework and finally the post-processing and visualization of the simulation output. We showcase various examples in §3 to demonstrate the applicability and efficacy of our framework across a range of scenarios, underscoring its potential to significantly advance the field of (embedded) tubular structure simulations. Finally, a list of potential enhancements is outlined in §4, delineating the path for future advancements before concluding the paper in §5.
2. Methods
2.1. Software/package choices
While other open-source frameworks such as Sim Vascular [9] and CRIMSON [38] facilitate physics-based simulations of tubular structures, they are primarily tailored for cardiovascular applications and have limited applicability in scenarios involving different physical processes, such as root-soil interactions or embedded tubular systems. In developing our pipeline, we have strategically selected a combination of powerful general-purpose open-source software packages to ensure efficiency and flexibility across various stages of our workflow. It is worth mentioning that all the software and packages were invoked through a unified MATLAB script to maximize automation capabilities.
The first part of the pipeline generates and processes surface meshes from different inputs (threedimensional imaging, CAD models, synthetic/skeletonized networks), using an open-source MATLABbased [39] toolbox called GIBBON (the Geometry and Image-Based Bioengineering add-ON) [40]. GIBBON includes a wide range of CAD, geometry and image processing tools. This toolbox also enables interfacing with free open-source software such as TetGen [41] and Geogram [42] that both possess high-quality meshing capabilities.
The second part of our pipeline involves extracting the centre-lines of the tubular structures. We integrate two versatile skeletonization algorithms, VesselVio [43] and TreeSkel [44], that have been used respectively on vascular and airways datasets.
For the automatic labelling of surface and volumetric mesh generation and conversion in the third and fourth segments of the pipeline, we rely on the features offered by GIBBON and meshio [45,46]. The meshio Python library is designed for handling mesh data. It provides functionalities for reading, writing and converting mesh files in various formats commonly used in computational simulations and visualization.
There are numerous open-source finite element software options which have been developed to solve free fluid flow and Darcy-scale porous-medium flow problems as well as coupled partial differential equation (PDE)-based transport models, such as Dumu· [47,48] and FEBio [49,50] (a detailed list can be found in [51]). We opt for FEniCS [52-54], which is a general-purpose open-source software framework designed to streamline the process of discretizing PDE using the FEM. FEniCS leverages the capabilities of the unified form language (UFL) [55] and the FEniCS form compiler (FFC) [56] to automatically produce optimized low-level C++ code for evaluating equations expressed in finite element variational forms. Additionally, it provides interfaces to PETSc and Trilinos for linear algebra computations. To use FEniCS, the user is required to provide the high-level variational form of the differential equation and does not have to perform any coding for the discretized form at the level of the cell and element scale. FEniCS is also well-suited for advanced solid mechanics simulations, particularly when integrated with MFront [57-59], allowing therefore a coupling between flow and mechanics.
Finally, for post-processing purposes, ParaView [60,61] (built on the visualization toolkit (VTK) library [62]) stands out as a comprehensive visualization tool, facilitating intuitive exploration and analysis of simulation outcomes. Its broad array of features and compatibility with multiple data formats render it an essential resource for extracting valuable information from simulation results.
Figure 1 illustrates the interaction of the different software/packages used in the Tube2FEM pipeline. Next, we describe the components of this workflow in turn.
2.2. Generating a watertight surface mesh
This workflow is intentionally flexible and can handle different inputs, including segmented threedimensional images issued from different imaging tools (X-ray computed tomography (XTC), optical coherence tomography (OCT), magnetic resonance imaging (MRI), etc.) and synthetic/skeletonized graph networks.
For the first type of input, namely, the segmented three-dimensional image, the creation of the surface mesh can be efficiently accomplished by employing the computer graphics marching cubes algorithm [63]. We use the implementation of the Marching Cube from the MATLAB file exchange [64] for this purpose. It is also possible to use the GIBBON functions im2patch and patch2tri to get a watertight, i.e. closed and manifold, surface mesh.
For the alternative input form, namely skeletonized or synthetic graph networks, continuous surface models are obtained by using level-set images. The procedure for creating level-sets typically includes three steps. The initial step involves embedding the spatial data within an image domain. The second step entails computing a distance function between the spatial data nodes and the image's voxel grid. The last step involves using the distance function to produce a (signed) level-set image. From there, continuous surface geometries can be obtained through the calculation of isosurfaces. A detailed description of this method can be found in [10] and an implementation of it already exists in GIBBON. Based on the graph network's size and the desired level of precision, calculating the distance function between the network nodes and the image voxel grid can become significantly time-consuming. To address this, we enhance the efficiency of the function by creating a binary image from the graph network, enabling the distance map calculation to occur only between a dilated version of this binary image and the nodes of the graph network.
2.3. Surface mesh processing
The watertight surface mesh, outlined in the previous section, exhibits voxel artefacts. A remeshing step is therefore required to eliminate these artefacts. We use Geogram which is an external library in GIBBON to remesh the triangulation defined by the mesh faces (F) and the vertices (V) [42]. In particular, the code Vorpalite is used and optimizes the Voronoi diagram from the point of view of sampling regularity. The resulting mesh consists of a near-isotropic distribution of triangles that transforms the stair-shaped surface patterns to a smoother appearance. We follow that by using the patchSmooth function available in GIBBON to complete the smoothing operation. In particular, we use the option of 'Humphrey's Class' that preserves best the overall volume of the watertight surface mesh.
However, after the smoothing step, small holes might appear in the surface mesh, leading to the loss of its watertightness. To repair it, we use PyMeshFix [65] from the PyVista library [66]. For compatibility purposes, we use the trimesh library [67] to read and write the defective and repaired mesh, respectively.
2.4. Skeletonization
In order to run physics-based simulations, it is necessary to assign boundary conditions on the edges of the network (for example to impose fluid pressure, velocity or solute concentration conditions). This requires creating flat surfaces at the terminal branches of the network and labelling them. While this task can be done manually for small networks, it becomes cumbersome and time consuming for medium and large networks. Thus, automating this process is crucial to make this pipeline viable for a broader spectrum of applications. One approach to automation involves baselining the cutting and labelling process around the terminal branches of a network's centre-line volume. Here we describe the step in the pipeline which creates the centre-line for the tubular volumes.
We assume that the centre-line of the vessels is equivalent to the skeleton obtained by skeletonization algorithms in the context of image-processing. An easily accessible open-source tool for threedimensional tubular-dataset analysis is therefore needed. Moreover, to facilitate automation, this tool must be designed to operate as a standalone executable application, negating the need for graphical user interface (GUI) interaction. VesselVio [43] meets these requirements and is compatible with both Windows and MacOS. It also provides high-level Python libraries, just-in-time compilation and parallel processing for rapid and detailed feature extraction techniques. This latest feature allows the transformation of the skeleton image to a graph format. For Linux users, the TreeSkel package [44], employing the minimum cost paths algorithm [68], offers an alternative solution for skeletonization.
The smoothed surface mesh was converted back into a binary image, due to the potential alterations introduced by the remeshing and smoothing operations on the original binary image. Vessel Vio, or alternatively TreeSkel, was then used to derive the skeleton of the network. Furthermore, plugins were created to import the graph format output into MATLAB. However, when Vessel Vio is used, overlaying the skeleton on the surface mesh reveals a noticeable 'retraction' effect, where the skeleton's edge points do not extend fully to the surface mesh. This is expected since VesselVio employs a custom implementation of a widely used medial axis parallel thinning algorithm [69]. This retraction poses a challenge as the edge points are crucial for guiding the slicing and labelling of surfaces to create boundary conditions. To address this issue, a ray-tracing algorithm [70] is employed to elongate the edge sub-segments until they intersect with the surface mesh. These intersection points are then integrated into the skeleton, serving as the new edge points.
2.5. Anisotropic mesh to enhance computational accuracy and efficiency
Before using the skeleton edge points and terminal branches to guide the cutting and labelling of the surface mesh, it is essential to perform an important step, especially for networks composed of vessels with a wide range of diameters. In fact, users may encounter a challenge in selecting an appropriate element size for the isotropic mesh. Choosing a small element size better preserves the morphology of the smaller tubular structures, yet results in excessive meshing of the larger ones, leading to a significant increase in the computational time required for the physics-based simulations. On the other hand, opting for a larger element size can lead to inaccuracies in representing the morphology of smaller tubes, along with a deficiency in the number of elements on the edge surfaces where boundary conditions are applied. To address this issue, employing an anisotropic meshing approach, where the degree of anisotropy is guided by the tubes' radii, is required. A minimum distance is thus calculated between the nodes of the surface mesh and the skeleton nodes using the minDist function found in GIBBON. This data is then fed to TetGen [41], an external library interfaced by GIBBON, which generates a volumetric mesh and from which we extract the anisotropic surface mesh. In this resulting mesh, smaller vessels are finely detailed, whereas larger vessels are meshed more coarsely.
It is important to note that the mesh refinement presented in this section is based solely on morphological considerations. However, depending on the underlying physics of the problem, further refinement may be required to adequately resolve element sizes, capture sharp gradients and minimize numerical errors-particularly in cases involving high Reynolds numbers and strong gradients at the interfaces.
2.6. Slicing and labelling to enable boundary condition assignment
To assign boundary conditions, it is essential to generate terminal cross-section surfaces through the tubular structure. For broader applicability, these surfaces should be flat, particularly to provide sufficient flexibility for use of vectorial entities, and they must contain a sufficient number of mesh elements to avoid numerical instabilities. This can be accomplished by slicing the surface mesh's terminal branches along their edges and subsequently remeshing to ensure the entire mesh remains watertight. Thus, accurate slicing operations are required at each terminal branch.
A slicing plane is defined by a point and a normal vector. For each slicing operation performed on the tubular geometry, the point of origin is an edge point of the skeleton, and the normal vector corresponds to terminal sub-segment of the skeleton. However, since a slicing plane is infinite, it might indiscriminately affect other sections of the mesh. To enhance the accuracy of this operation, a Boolean constraint is introduced, ensuring that only interconnected elements nearest to the skeleton edge point are impacted. This is achieved by using the minDist function to apply the Boolean condition, followed by the triSurfSlice function for executing precise slicing on the targeted branch. Both functions are part of the GIBBON toolbox.
Following that, the now open surface mesh is re-sealed by remeshing the terminal branches using a series of functions, the most important of which is regionTriMesh3D (available in GIBBON). Each newly formed surface is assigned a distinct label, and an automatically generated text file, containing boundary conditions applied to those surfaces, is created. This text file can be directly communicated to FEniCS and is essential in the automation of the finite element simulations.
2.7. Volumetric meshing
The labelled surface mesh is now prepared for a second processing by TetGen [41] to create a volumetric mesh, ready for physics-based simulations. We input specific parameters into TetGen, including meshing options, faces, nodes, face labels and the number of regions. The runTetGen command in GIBBON is then used to produce a tetrahedral mesh. The output includes the tetrahedral elements, element IDs, faces, face IDs and the coordinates of the nodes. It is also possible to convert the tetrahedral elements to hexahedral ones using the tet2hex function. GIBBON also provides the capability to extrude the triangular surface elements into pentahedra, making the mesh suitable for applications like fluid-structure interaction.
2.8. Mesh exporting and conversion
Different finite element (FE) software packages require various mesh formats. The FE software used in our pipeline is FEniCS, which supports a limited number of formats, including .xml and .xdmf. To accommodate this, we have developed a function capable of converting the output from a TetGen mesh into a version 2 ASCII Gmsh file (.msh format). Subsequently, a Python script employing the meshio library is used to convert the .msh file into a .xdmf format, ensuring compatibility with FEniCS.
2.9. Mathematical models
We integrate physics-based mathematical models into this workflow, which are commonly utilized in simulating processes that involve (embedded) tubular objects. The first model is the incompressible Navier-Stokes model to simulate the fluid dynamics in the tubular objects. The second model is to simulate the advection-diffusion-reaction of specie(s) (solutes, oxygen, drugs, nanoparticles, etc.) in the embedded or embedding domain. While these are specific examples, the broad framework can be easily adapted to other physical models (e.g. solid mechanics).
2.9.1. Computational fluid dynamics
We assume that the fluids of interest are non-compressible Newtonian fluids and so fluid flow can be described via the incompressible Navier-Stokes equations, written in the form
... (2.1)
... (2.2)
where f is a given force per unit of volume, u is the velocity of the fluid, p its pressure, p its density and o(u,p) is the stress tensor. For Newtonian fluids, the stress tensor is given by
... (2.3)
where I is the identity matrix and e(u) is the strain-rate tensor,
... (2.4)
The Navier-Stokes equations exhibit nonlinearity, transience and a non-trivial coupling between pressure and velocity. While coupled methods exist to solve the system of equations, it is more popular to decompose the problem into several more manageable equations. The initial splitting scheme was proposed by Chorin [71] and Témam [72]. This method was further refined by Goda [73], who introduced a velocity correction step, and his approach became widely recognized in the literature as the incremental pressure correction scheme (IPCS). Despite the existence of alternative splitting schemes, a thorough comparison presented in [74] indicates that, in terms of both efficiency and accuracy, the IPCS generally outperforms other methods. Moreover, the IPCS method is more straightforward to implement when contrasted with alternative splitting schemes. It consists of solving three equations:
- the tentative velocity step (a convection-diffusion-reaction equation),
- the pressure-correction step (a Poisson equation),
- the velocity update step (a projection).
We customize the implementation of this method provided in the FEniCS tutorial [75].
2.9.2. Species transport in the embedded or embedding domains
We include general-purpose transport of chemical species, described using an advection-diffusion- reaction equation,
... (2.5)
Here the temporal evolution of a chemical species of concentration c is linked to its spatial evolution through diffusion, advection by a vector field u and/or a reaction within the domain. It is worth noting that this equation is usually coupled (weakly) with the aformentioned Navier-Stokes through the velocity field u especially when the advection effects are important.
2.9.3. Mixed-dimension models: 30-10 coupling
Our framework is generic and suitable for flow systems with high, intermediate or low Reynolds numbers. However, there are numerous situations where it is practical to consider steady-state flows at low Reynolds numbers, such as blood flow in networks of small blood vessels and plant root systems. In these cases, it is reasonable to disregard changes in flow velocity and inertial effects. This is reflected mathematically by eliminating both the time-dependent and inertial terms in the Navier-Stokes equations, leading to the use of the simpler Stokes flow equations. Moreover, in such situations, a one-dimensional model of the tubular structure is often adequate to capture the essential aspects of flow physics and is often used in the literature [14,15,17]. Nevertheless, when the focus is on how the one-dimensional tubular objects interact with their three-dimensional surrounding environment, a significant challenge emerges due to the mixed-dimensional (3D-1D) nature of the problem. In fact, various numerical methods vary in their approach to bridging the dimensional incompatibility between network and embedding domain by distributing source terms. The contribution of the source term within the encompassing bulk medium can be represented through line source terms [76,77], surface source terms [78,79] or volume source terms [4,18]. In many of these methods, it is assumed that the diameters of the network branches are significantly smaller than the dimensions of the surrounding domain. This facilitates a simpler mesh creation process that results in an overestimation of the bulk volume by covering the space occupied by the tubular structures. Additionally, due to the elimination of their physical structure, the tubes no longer offer a resistance to flow in the embedding domain. Finally, as these mixed-dimension models typically represent the network by a series of finite cylinders, inaccuracies arise at junctions where source terms should be distributed.
Most of these problems in 3D-1D coupling settings can be avoided if the exact network surface mesh is known. Our framework allows for the building of such continuous surface meshes from networks of onedimensional segments, as shown in 82.2. Consequently, simulations conducted on the network of onedimensional segments can be projected onto the network surface mesh using an minimal distance-based approach between the nodes of the one-dimensional segments and those of the surface mesh. Then, using the TetGen [41] interface with GIBBON, a conforming mesh that integrates the embedded system and its surrounding environment can be created. The interface conformity between the two domains allows the direct use of the distributed source terms on the network surface without the necessity of extrapolation.
2.10. Scientific visualization and post-processing
Geometries and meshes are visualized in MATLAB using the GIBBON open-source toolbox, whereas the simulation results are exported in .vtk format and are visualized and post-processed using ParaView [60,61] (built on the VTK library [62]). For automation purposes, the use of the ParaView GUI must be avoided. For that, generalized macros were written and called from MATLAB through a Command-Line Interface (CLI).
3. Results and discussion
This section highlights five case studies that comprehensively illustrate the capabilities of the Tube2FEM framework. Particularly, we present instances of tubular geometries corresponding to three distinct applications (lung airways, root-soil systems, micro/macro circulation), each emphasizing various aspects of the framework, including surface mesh reconstruction, meshing, CFD simulations, 3D-1D approaches, visualization, etc. We note that all the case studies presented in the upcoming section are automated, with each having all of its commands centralized within a single MATLAB script. Table 1 summarizes all the key user-defined parameters, method specifications and some of the output results for all of the five case studies. It is also worth noting that all computations were carried out on a personal laptop with the following specifications: 11th Gen Intel(R) Core(TM) 17-1185С7 processor @ 3.00 GHz, 1805 Mhz, four cores, eight logical processors and 32 GB of RAM. Finally, with the exception of the cases featuring patient-specific data, all the codes are made available on GitHub at: https://github.com/CheikhSleiman/Tube2 FEM.
3.1. Case study 1: radius-dependent anisotropic mesh generation of lung airways
For our initial example, we employ Tube2FEM to generate simulation-ready meshes from segmented lung airways acquired using medical CT scans. The primary difficulty with this type of tree-like tubular geometry is its topological complexity, since it contains more than 10 branching generations, and there is significant variation in diameter between the smallest and largest airways. Indeed, the challenge stems from choosing a mesh element size that maintains the airways" morphology without leading to increased computational time for the simulations. Tube2FEM addresses this issue by employing a radius-dependent anisotropic meshing approach. The process begins with the segmented lung airways volume, as depicted in figure 2a, from which the tubular structure's centre-line is determined using skeletonization algorithms, as detailed in §2.4. Additionally, the surface mesh is generated from the three-dimensional image using the marching cubes algorithm or similar techniques, as described in §2.2. However, this mesh initially exhibits voxel artefacts, necessitating further processing steps such as remeshing and smoothing, discussed in §2.3. Subsequently, the processed surface mesh is overlaid with the skeleton, as shown in figure 2c. The next phase involves projecting the centre-line radii field onto the surface mesh, based on the minimal distance criterion outlined in §2.5. This projection creates a radius-dependent weighting map that guides the mesh's anisotropy, shown in figure 2d. The final comparison between the preprocessed surface mesh, the isotropically remeshed and smoothed version, and the anisotropic mesh is presented in figure 2e-g respectively. Finally, figure 2h-j respectively illustrates the clipped terminal branches, sealed outlets and labelled boundary regions, demonstrating the readiness of the mesh for numerical simulation.
3.2. Case study 2: multi-domain mesh generation of a root-soil system
In this case study, Tube2FEM is used to generate a multi-domain, simulation-ready mesh of a root network embedded within a soil medium. The segmented root network data, previously examined in [4], was accessible online only in Gmsh (.msh format) as referenced in [80]. Since the root-soil system was originally imaged through X-ray tomography before being segmented and processed into a meshing format, we converted the dataset back into a standard imaging format (.tiff) and reprocessed it from the beginning using our framework. We then use this segmented (.tiff) dataset illustrated in figure За as the input for our procedure. As in the previous case study, we extract the surface mesh using the marching cubes algorithm or equivalent techniques, as outlined in §2.2. Remeshing and smoothing processes are then carried out to eliminate voxel artefacts and achieve a desired mesh density, as illustrated in figure 3c,d. Indeed, in order to produce reliable results in hydraulic root-soil simulations, the mesh needs to be carefully locally refined at the interface. This is particularly important as drying soil around the roots creates large and highly localized pressure gradients. The effective way to capture these gradients is through local mesh refinement at the interface [4]. To represent the plant's transpiration activity above the soil, a designated surface was labelled to impose a transpiration flux boundary condition, as illustrated in figure 3b. A soil medium is then built around the root mesh using CIBBON's triBox function, and both media are inputted into TetGen, as detailed in §2.7, to create a multi-domain tetrahedral mesh. Figure 3e shows both domains, featuring a mesh density gradient from the interface to the boundaries of the soil domain.
3.3. Case study 3: generating surface and volumetric meshes of a microcirculation blood network
This case study demonstrates the ability of Tube2FEM to generate simulation-ready labelled surface and volumetric meshes from networks of one-dimensional tubular objects. These networks can either be synthetic (for example, via generative algorithms) or result from skeletonization processes applied to imaging datasets. While the procedure is general, we apply it to a widely adopted microvascular network dataset used in the literature to explore microcirculatory blood flow and transport processes, studied initially in [81].
The level-set method, detailed in §2.2, is used here to create a watertight surface mesh. The starting point is the vascular network centre-line over which the vessel radii field is presented, as illustrated in figure 4a. The second step entails embedding the spatial graph within an image domain and then computing a distance function between the centre-line nodes and the image's voxel grid. This method can become significantly time consuming if the image domain is not carefully chosen. To address this, we create a binary image from the graph network as illustrated in figure 4b. Several image dilation operations are then executed to establish the required image domain. The signed level-set image, shown is figure 4c can be interpreted as the normalized minimal distance map between the centre-line nodes and the dilated binary image. The isosurface, corresponding to level-set intensity = 1, can be extracted from the level-set image. The resulting surface mesh, which has been remeshed and smoothed to eliminate artefacts, is displayed in figure 4d. The following step involves applying the slicing and labelling rithm, in §2.5, to assign boundary conditions. This is achieved by overlaying the centre-line network with the processed surface mesh and identifying the edge nodes and the directionality of the edge segments. Finally, having the watertight labelled surface mesh illustrated in figure 4e, it is possible now to generate the tetrahedral mesh using the TetGen interface in GIBBON. The volumetric mesh in shown in figure 4f.
3.4. Case study 4: 3D-1D coupling strategy to simulate tissue perfusion
This case study aims to demonstrate the ease of 3D-1D coupling between an embedding tubular network and its surrounding domain using Tube2FEM. For convenience, we apply this procedure on the same graph network used in case study 3 and the same reconstructed surface mesh which was created using the level-set method. The final goal is to simulate the extra-vascular solute transport driven by a heterogeneous distribution of solutes within the vascular network.
In this specific scenario, a one-dimensional representation of each vessel suffices to capture the flow physics. This renders the computational time cost of the simulated intravascular processes very cheap in comparison with a three-dimensional case scenario. Consequently, the initial step in this process volves solute transport across a network of one-dimensional vessels. This simulation can be conducted using methods such as finite element, finite volume among others [14,15]. However, when tention to the interactions between the one-dimensional tubular objects and their three-dimensional surroundings, a considerable challenge arises from the mixed-dimensional (30-10) nature of the boundary. As discussed in §2.9, different numerical methods adopt various strategies to overcome the dimensional mismatch between the network and the embedding domain by distributing source terms (line, surface or volume sources terms). For this case study, we distribute the source terms across the surface, surface, since Tube2FEM is capable of reconstructing the exact surface mesh of the vascular network, as demonstrated in case study 3. To achieve this, we proceed to overlap the simulation results on the centre-line with the reconstructed surface mesh as illustrated in figure 5b. Then the solute field on the network of one-dimensional vessels is projected on the mesh using a minimal distance computation between the centre-line nodes and the mesh ones. This results in the surface field shown in figure 5c. Following this, we create the surrounding tissue volumetric mesh using TetGen-GIBBON interface. In the input options for TetGen, we ensure that the volumetric tissue mesh should be hollowed out from the vascular domain, which is represented by the reconstructed surface mesh. We also ensure that the mesh density of the vascular surface mesh is sufficiently fine to more effectively resolve the volumetric mesh density gradient near the vascular-tissue interface. The result of these operations can be visualized in figure 5d. Additionally, we note that the internal surface of the volumetric mesh, created by the removal of the vascular domain, conforms to the reconstructed vascular surface mesh. Such conformity is essential for facilitating consistent interactions between the different domains. The volumetric mesh produced is first exported in the (.msh) format and subsequently converted to the (.xdmf) format using the meshio package. This latter format ensures compatibility with the FEniCS software, which is employed to carry out the finite element analysis. We also export the solute field evaluated on the vascular surface mesh in a tabular format (xy, YN,ZN,Cs), Where xn, YN and zy are the vascular surface mesh nodes coordinates and с, is the concentration of the solute on the node М. This table is then communicated to FEniCS and is used in assigning boundary conditions. Finally, we use the Poisson equation to model a simple steady-state solute diffusion process in the tissue as illustrated in figure 5e. We note that this procedure can be used for much more complex problems, and can be further developed, for instance to account for time-dependent intravascular flow variations and two-way coupling.
3.5. Case study 5: detailed computed-tomography-to-simulation procedure for tubular geometries: application to portal vein network
This final case study is designed to showcase a comprehensive framework for conducting physics-based simulations on tubular geometries obtained via imaging. While the framework is broadly applicable, we demonstrate its use through a portal vein network derived from a patient's CT scan. This segmented network is illustrated in figure 6. Similarly to previous case studies, the portal vein surface mesh can be extracted using the marching cubes algorithm or an equivalent method detailed in §2.2. The surface is then remeshed to a desired density and smoothed to get rid of geometrical artefacts. Following that, we call VesselVio software in MATLAB and skeletonize the binary image. The overlap of the resulting centreline and the surface mesh can be seen in figure 6b. Subsequently, the skeleton edge nodes are identified, and when combined with the direction of the skeleton terminal sub-segments, they determine the slicing planes. The surface mesh after the slicing operations is visualized in figure 6c. The open geometry is then sealed by remeshing and labelling the open surfaces near the skeleton terminal branches. The detailed slicing, sealing and labelling procedures are detailed in §2.6. In total, figure 6d shows 47 labelled surfaces, 46 of which are created by sealing the sliced geometry and one being the vessel's wall. Those labelled surfaces are crucial to assign boundary conditions for the physics-based models. For the Navier-Stokes (NS) simulation, we apply a time-dependent velocity to the inlet, zero pressure to the outlets and zero velocity to the vessel walls (no-slip boundary condition). Figure 6e illustrates the velocity streamlines at a specific time point, resulting from solving the NS problem. We then save all the velocity fields at the different time steps and use them to inform the advection of the chemical specie in the advection- diffusion-reaction problem. Figure 6f shows the resulting concentration of the chemical field within the portal vein network. It is important to highlight that those simulation results are demonstrative of the capacity of the pipeline to solve complex flow and flow-related problem. The set of parameters and associated boundary conditions used in those simulations do not necessarily represent patient-specific physiological ranges.
4. Perspective
Anticipating future advancements of our pipeline, we highlight critical improvements targeted at enhancing the accuracy and reducing the computational time required for simulating (embedded) tubular structures. The outlined perspectives below may be explored in future work to enhance the framework's efficiency and technical capabilities.
- Use of high-performance finite element Navier-Stokes solver. Oasis is a Python-based highperformance finite element Navier-Stokes solver [82]. It utilizes foundational components from the FEniCS project and is designed for addressing large-scale applications involving intricate geometries, particularly on massively parallel computing clusters.
- Achieve faster periodic convergence in case of cyclic loading. The computational expense of threedimensional biomedical fluid dynamics simulations is often attributed to the need for computing multiple cycles (cardiac, respiratory, etc.) before achieving a periodic solution. The work of Pfaller et al. [83] showed the possibility of fastening the periodic convergence by generating appropriate initial conditions using the simulation results of reduced-order one-dimensional models.
- Automating a comprehensive verification of CFD simulation quality. A few decades ago, multiple fluid mechanics editorial policy boards stated that 'there is a need for higher standards on the control of numerical accuracy" and that 'a single calculation of a fixed grid is not acceptable" [84]. However, recent research indicates that even grid refinement is insufficient for evaluating simulation quality. Instead, a comprehensive CFD investigation involving solver numerics, mesh and time-step refinement is essential. For instance, it has been shown in [85] that robust and minimally dissipative CFD solvers can tolerate surprisingly coarse resolutions, whereas solvers using low-order and/or stabilization schemes may require much higher resolutions to detect relevant flow patterns.
- Include lumped parameter network (LPN) models for boundary conditions. While ensuring an accurate geometric representation of the computational domain is paramount in simulating fluid flow, it is equally essential to emphasize the importance of realistic boundary conditions. In general, the distribution of flow and pressure field within the simulated domain is often unknown and challenging to specify at the inflow/outflow boundaries. In fact, the physiology of the flow at the inlet/outlet is often too complex to be represented by a fixed-value Dirichlet of Neumann-type boundary condition. An alternative strategy involves linking the solution at the inflow/outflow boundaries of the modelled domain with simplified lumped parameter (zerodimensional models) or one-dimensional models representing the missing domain [86,87]. This often requires solving an ODE to describe the physics of the flow at the boundary condition level.
- Implementing two-way coupling for embedded systems. Regarding the embedded systems, it is necessary to establish a two-way coupling between the embedded and the embedding domain. This coupling is necessary to ensure the conservation of mass (e.g. for mass transfer problems in vascular/tissue, root-soil and well-soil systems), momentum (e.g. to account for mechanical deformations in expanding or contracting airway/lung systems) and energy (e.g. for heat transfer problem in heat exchangers).
5. Conclusion
In this paper, we developed Tube2FEM, an open-source highly automated framework for flow and flow-related processes simulation in (embedded) tubular structures. This system is capable of handling various input forms, including segmented tomographies and synthetic or skeletonized networks. It offers a range of features that include creating surface models, automated slicing and labelling for the assignment of boundary conditions. Furthermore, it supports surface and volumetric mesh generation and allows performing simulations using the FEM. It also offers capabilities for post-processing and visualizing the results. The framework utilizes a selection of carefully curated open-source softwares and libraries, chosen to guarantee a high degree of versatility for various simulation scenarios. Among the most significant of these are GIBBON, which is essentially employed for image and/or geometry processing, FEniCS for conducting FE simulations, and Paraview, which is utilized for the post-processing tasks.
The case studies presented in this work demonstrate the comprehensive capabilities of the Tube2FEM framework across a range of applications of (embedded) tubular geometries, from lung airways and root-soil interactions to macro/microcirculation and tissue perfusion. We showcase in all case studies advanced meshing capabilities including surface mesh creation and extraction, isotropic/anisotropic remeshing, smoothing, multi-domains and labelling for automated boundary conditions. Additionally, we showcase in case studies 4 and 5 a range of physics-based models applied to tubular and embedded tubular systems. Particularly, a framework for 3D-1D mixed dimensions simulation for embedded tubular system was created and was used to simulate the solute transport from a microcirculation network to its surrounding biological tissue. In addition, a CT-to-simulation framework was also developed to conduct CFD and advection-diffusion-reaction simulations on image-based geometries.
Finally, future work will focus on enhancing the automation and efficiency of the workflow, expanding the framework's capability to handle complex boundary conditions and couplings, and on integrating high-performance solvers to achieve faster computational time while maintaining the results precision.
Ethics. Lung data: approval for this study of clinical and CT data was obtained from the local research ethics committees and Leeds East Research Ethics Committee: 20/YH/0120. Liver data: all data collected was with patient consent. The ethics approval for anonymized laparoscopic data collection was titled Image guidance with image overlay for laparoscopic surgery" and was approved by the NRES Committee London- Stanmore. Reference no.: 14/L0/1264 on 16 November 2017.
Data accessibility. Data and relevant code for this research work are stored in GitHub: [88] and have been archived within the Zenodo repository: [89].
Declaration of Al use. We have used Al-assisted technologies in creating this article.
Authors' contributions. H.C.S.: conceptualization, data curation, formal analysis, methodology, software, visualization, writing - original draft; K.M.M.: formal analysis, investigation, software, writing - review and editing; D.C.D.O.: formal analysis, investigation, validation, writing -review and editing; J.J.: data curation, formal analysis, methodology, writing -review and editing; N.M.: data curation, investigation, writing - review and editing; B.R.D.: data curation, formal analysis, investigation, writing-review and editing; S.W.-S.: formal analysis, funding acquisition, investigation, project administration, supervision, writing-review and editing; R.J.S.: conceptualization, formal analysis, investigation, project administration, resources, supervision, validation, writing - review and editing.
All authors gave final approval for publication and agreed to be held accountable for the work performed therein.
Conflict of interest declaration. Joseph Jacob (JJ) reports fees from Boehringer Ingelheim, Roche, NHSX, Takeda and GlaxoSmithKline and grants from Gilead Sciences and Microsoft Research unrelated to the submitted work.
Funding. The authors are grateful for funding from CRUK (C44767/A29458 and C23017/A27935; H.C.S., S.W.-S., R.J.S.) and from the EPSRC (EP/W007096/1; H.C.S., S.W.-S., R.J.S. and EP/V034537/1; D.C.D.O.). J.J. was supported by Wellcome Trust Clinical Research Career Development Fellowship 209553/Z/17/Z, Wellcome Trust Career Development Fellowship 227835/7/23/7 and the NIHR Biomedical Research Centre at University College London. The liver data of case study 5 were collected as part of an NIHR Invention for Innovation (i4i) Programme grant (II-LA-1116-20005) awarded to Prof. Clarkson and Prof. Davidson. Views expressed are those of the authors and not necessarily those of the NIHR or the Department of Health and Social Care.
References
1. RooseT, Fowler AC. 2004 A mathematical model for water and nutrient uptake by plant root systems. J. Theor. Biol. 228, 173-184. (doi:10.1016/ jjtbi.2003.12.013)
2. Rudolph-Mohr N, Bereswill 5, Totzke С, Kardjilov № Oswald SE. 2021 Neutron computed laminography yields 3D root system architecture and complements investigations of spatiotemporal rhizosphere patterns. Plant Soil 469, 489-501. (doi:10.1007/s11104-021-05120-7)
3. Anselmucci Е, Floriana A, Ando E, Viggiani ©, Lenoir № Peyroux В, Sibille L. 2021 The use of x-ray tomography to investigate soil deformation around growing roots. Géotechnique Lett. 11, 1-19. (doi:10.1680/jgele.20.00114)
4. Koch T.2022 Projection-based resolved interface 1D-3D mixed-dimension method for embedded tubular network systems. Comput. Math. Appl. 109, 15-29. (doi:10.1016/j.camwa.2022.01.021)
5. Koch T, Wu H, Schneider М. 2022 Nonlinear mixed-dimension model for embedded tubular networks with application to root water uptake. J. Comput. Phys. 450, 110823. (doi:10.1016/j.jcp.2021.110823)
6. Bereswill 5, Gatz-Miller H, Su D, Tôtzke С, Kardjilov N, Oswald SE, Mayer KU. 2023 Coupling non-invasive imaging and reactive transport modeling to investigate water and oxygen dynamics in the root zone. Vadose Zone J. 22, e20268. (doi:10.1002/vzj2.20268)
7. Espinoza DN, Santamarina JC. 2010 Ant tunneling- a granular media perspective. Granul. Matter 12, 607-616. (doi:10.1007/s10035-010-0202y)
8. Belachew M, Yamamoto К, Nichols E, Zhang D, Frost JD, Arson С. 2024 Ant nest geometry, stability, and excavation-inspiration for tunneling. Acta Geotech. 19, 1295-1313. (doi:10.1007/s11440-024-02232-7)
9. Updegrove A, Wilson NM, Merkow J, Lan H, Marsden AL, Shadden SC. 2017 SimVascular: an open source pipeline for cardiovascular simulation. Ann. Biomed. Eng. 45, 525-541. (doi:10.1007/s10439-016-1762-8)
10. Moerman KM et al. 2022 Development of a patient-specific cerebral vasculature fluid-structure-interaction model. J. Biomech. 133, 110896. (doi: 10.1016/j.jbiomech.2021.110896)
11. Leoni M, Szasz J, Meier J, Gerardo-Giorda L. 2022 Blood flow but not cannula positioning influences the efficacy of Veno-Venous ECMO therapy. Sci. Rep. 12, 20950. (doi:10.1038/s41598-022-23159-7)
12. Кай |, Pichelin M, Montesantos S, Murdock A, Fromont S, Venegas J, Caillibotte G. 2017 The influence of lung volume during imaging on CFD within realistic airway models. Aerosol Sci. Technol. 51, 214-223. (doi:10.1080/02786826.2016.1254721)
13. Sadafi H, Monshi Tousi N, De Backer W, De Backer J. 2024 Validation of computational fluid dynamics models for airway deposition with SPECT data of the same population. Sci. Rep. 14, 5492. (doi:10.1038/s41598-024-56033-1)
14. d'Esposito A et al. 2018 Computational fluid dynamics with imaging of cleared tissue and of in vivo perfusion predicts drug uptake and treatment responses in tumours. Nat. Biomed. Eng. 2, 173-187. (doi:10.1038/s41551-018-0306-y)
15. Berg M, Davit Y, Quintard M, Lorthois S. 2020 Modelling solute transport in the brain microcirculation: is it really well mixed inside the blood vessels? J. Fluid Mech. 884. (doi:10.1017/jfm.2019.866)
16. KremhellerJ, Vuong A, Schrefler BA, Wall WA. 2019 An approach for vascular tumor growth based on a hybrid embedded/homogenized treatment of the vasculature within a multiphase porous medium model. Int. J. Numer. Methods Biomed. Eng. 35, e3253. (doi:10.1002/cnm.3253)
17. Sweeney PW, d'Esposito A, Walker-Samuel S, Shipley RJ. 2019 Modelling the transport of fluid through heterogeneous, whole tumours in silico. PLoS Comput. Biol. 15, 1006751. (doi:10.1371/journal.pchi. 1006751)
18. Koch T, Schneider M, Helmig R, Jenny P. 2020 Modeling tissue perfusion in terms of 1d-3d embedded mixed-dimension coupled problems with distributed sources. J. Comput. Phys. 410, 109370. (doi:10.1016/j.jcp.2020.109370)
19. Koppl T, Vidotto E, Wohlmuth ?. 2020 A 3D-1D coupled blood flow and oxygen transport model to generate microvascular networks. Int. J. Numer. Methods Biomed. Eng. 36, e3386. (doi:10.1002/cnm.3386)
20. Fritz M, Kóppl T, Oden JT, Wagner A, Wohlmuth ?, Wu ?. 2022 A 1D-0D-3D coupled model for simulating blood flow and transport processes in breast tissue. Int. J. Numer. Method. Biomed. Eng. 38, e3612. (doi:10.1002/cnm.3612)
21. Aslam Bhutta MM, Hayat ? Bashir MH, Khan AR, Ahmad KN, Khan $. 2012 CFD applications in various heat exchangers design: a review. Appl. Therm. Eng. 32, 1-12. (doi:10.1016/j.applthermaleng.2011.09.001)
22. Jackson GT, Balhoff MT, Huh C, Delshad M. 2011 CFD-based representation of non-Newtonian polymer injectivity for a horizontal well with coupled formation-wellbore hydraulics. J. Pet. Sci. Eng. 78, 86-95. (doi:10.1016/j.petrol.2011.05.007)
23. Jafari A, Hasani M, Hosseini M, Gharibshahi R. 2020 Application of CFD technique to simulate enhanced oil recovery processes: current status and future opportunities. Pet. Sci. 17, 434-456. (doi:10.1007/s12182-019-00363-7)
24. Levoni P, Angeli D, Stalio ?, Agnani E, Barozzi GS, Cipollone M. 2015 Fluid-dynamic characterisation of the Mont Blanc tunnel by multi-point airflow measurements. Junn. Undergr. Space Technol. 48, 110-122. (doi:10.1016/j.tust.2015.03.006)
25. O'Rourke TD, Jung JK, Argyrou ?. 2016 Underground pipeline response to earthquake-induced ground deformation. Soil Dyn. Earthg. Eng. 91, 272-283. (doi:10.1016/).soildyn.2016.09.008)
26. Choudhury D, Chaudhuri CH. 2023 A critical review on performance of buried pipeline subjected to pipe bursting and earthquake induced permanent ground deformation. Soil Dyn. Earthq. Eng. 173, 108152. (doi:10.1016/j.soildyn.2023.108152)
27. Viggiani G, Tengattini A. 2019 Recent developments in laboratory testing of geomaterials with emphasis on imaging. In Proc. of the XVII ECSMGE2019: Geotechnical Engineering Foundation of the Future, pp. 5267-5293. The Icelandic Geotechnical Society. (doi:10.32075/17ECSMGE-20191112)
28. Walsh CL et al. 2021 Imaging intact human organs with local resolution of cellular structures using hierarchical phase-contrast tomography. Nat. Methods 18, 1532-1541. (doi:10.1038/s41592-021-01317-x)
29. Xian RP et al. 2022 A multiscale X-ray phase-contrast tomography dataset of a whole human left lung. Sci. Data 9, 264. (doi:10.1038/s41597-02201353-y)
30. Tengattini A, Lenoir N, Ando E, Giroud B, Atkins D, Beaucour J, Viggiani G. 2020 NeXT-Grenoble, the Neutron and X-ray tomograph in Grenoble. Nucl. Instruments Methods Phys. Res. Sect. 968, 163939. (doi:10.1016/j.nima.2020.163939)
31. Walsh ?, Holroyd NA, Finnerty E, Ryan SG, Sweeney PW, Shipley RJ, Walker-Samuel $. 2021 Multifluorescence high-resolution episcopic microscopy for 3D imaging of adult murine organs. Adv. Phot. Res. 2, 2100110. (doi:10.1002/adpr.202100110)
32. Holroyd NA, Li Z, Walsh C, Browng E, Shipley RJ, Walker-Samuel S. 2023 tUbe net: a generalisable deep learning tool for 3D vessel segmentation. Microbiology. (doi:10.1101/2023.07.24.550334)
33. Walsh CL, Berg M, West H, Holroyd NA, Walker-Samuel $, Shipley RJ. 2024 Reconstructing microvascular network skeletons from 3D images: what is the ground truth? Comput. Biol. Med. 171, 108140. (doi:10.1016/j.compbiomed.2024.108140)
34. Guy AA, Justin AW, Aguilar-Garza DM, Markaki AE. 2020 3D printable vascular networks generated by accelerated constrained constructive optimization for tissue engineering. IEEE Trans. Biomed. Eng. 67, 1650-1663. (doi:10.1109/tbme.2019.2942313)
35. Brown E, GuyA, Holroyd N, Sweeney P. 2023 Physics-informed deep generative learning for quantitative assessment of the retina. See http://dx. doi.org/10.1101/2023.07.10.548427.
36. Cruz de Oliveira DM. 2024 Feto-placental vascular generator: first release (version 1). Zenodo. See https://zenodo.org/doi/10.5281/zenodo. 10557280.
37. Berg M, Holroyd N, Walsh ?, West H, Walker-Samuel 5, Shipley ?. 2022 Challenges and opportunities of integrating imaging and mathematical modelling to interrogate biological processes. Int. J. Biochem. Cell Biol. 146, 106195. (doi:10.1016/j.biocel.2022.106195)
38. Arthurs CJ et al. 2021 CRIMSON: An open-source software framework for cardiovascular integrated modelling and simulation. PLoS Comput. Biol. 17, e1008881. (doi:10.1371/journal.pchi. 1008881)
39. MATLAB. 2022 MATLAB. 9.12.0 (R2022a). Natick, MA, USA: The MathWorks Inc.
40. Moerman KM. 2018 GIBBON: The geometry and image-based bioengineering add-on. J. Open Source Softw. 3, 506. (doi:10.21105/j0ss.00506)
41. Si H. 2015 TetGen, a Delaunay-based quality tetrahedral mesh generator. ACM Trans. Math. Softw. 41, 1-36. (doi:10.1145/2629697)
42. Lévy B, Bonneel N. 2013 Variational anisotropic surface meshing with Voronoi parallel linear enumeration. In Proc. of the 21st Int. Meshing Roundtable, pp. 349-366. Berlin, Germany: Springer Berlin Heidelberg. (doi:10.1007/978-3-642-33573-0_21)
43. Bumgarner JR, Nelson RJ. 2022 Open-source analysis and visualization of segmented vasculature datasets with VesselVio. Cell Rep. Methods 2, 100189. (doi:10.1016/j.crmeth.2022.100189)
44. Pakzad A, SzmulA. 2022 Open curvilinear skeletonization based on path cost for airway trees. See https://github.com/ashkanpakzad/TreeSkel.
45. GitHub. 2023 GitHub-nschloe/meshio. See https://github.com/nschloe/meshio (accessed 26 October 2023).
46. Schlômer N et al. 2018 Nschloe/Meshio V1.11.7. Zenodo. https://zenodo.org/record/1173116
47. Flemisch ? et al. 2011 DuMu·: DUNE for multi-phase, component scale physics, flow and transport in porous media. Adv. Water Resour. 34, 1102-1112. (doi:10.1016/j.advwatres.2011.03.007)
48. Koch Tet al. 2021 DuMu· 3 - an open-source simulator for solving flow and transport problems in porous media with a focus on model coupling. Comput. Math. Appl. 81, 423-443. (doi:10.1016/j.camwa.2020.02.012)
49. FEBio website. FEBio Software Suitedate. See https://febio.org (accessed 1 March 2025).
50. Maas SA, Ellis BJ, Ateshian GA, Weiss JA. 2012 FEBio: finite elements for biomechanics. J. Biomech. Eng. 134, 011005. (doi:10.1115/1.4005694)
51. Bilke L, Flemisch ?, Kalbacher T, Kolditz O, Helmig ?, Nagel T. 2019 Development of open-source porous media simulators: principles and experiences. Transp. Porous Media 130, 337-361. (doi:10.1007/s11242-019-01310-1)
52. Logg A, Wells G. 2010 DOLFIN: Automated finite element computing. ACM Trans. Math. Softw. 3, 37. (doi:10.1145/1731022.1731030)
53. Logg A, Wells G, Mardal KA. 2012 Automated solution of differential equations by the finite element method. In The FEniCS book (eds A Logg, KA Mardal, G Wells), vol. 84. Berlin, Germany: Springer Berlin Heidelberg. (doi:10.1007/978-3-642-23099-8)
54. Alnaes M et al. 2015 The FEniCS project version 1.5. Arch. Num. Softw. 3. (doi:10.11588/ans.2015.100.20553)
55. Alnaes MS, Logg A, Olgaard KB, Rognes ME, Wells GN. 2014 Unified form language. ACM Trans. Math. Softw. 40, 1-37. (doi:10.1145/2566630)
56. Kirby RC, Logg A. 2006 A compiler for variational forms. ACM Trans. Math. Softw. 32, 417-444. (doi:10.1145/1163641.1163644)
57. Helfer T, Michel B, Proix JM, Salvo M, Sercombe J, Casella M. 2015 Introducing the open-source mfront code generator: application to mechanical behaviours and material knowledge management within the PLEIADES fuel element modelling platform. Comput. Math. Appl. 70, 994-1023. (doi:10.1016/j.camwa.2015.06.027)
58. Helfer T, Bleyer J, Frondelius T, Yashchuk I, Nagel T, Naumov D. 2020 The MFrontGenericInterfaceSupport project. J. Open Source Softw. 5, 2003. (doi:10.21105/j0ss.02003)
59. Helfer T, Bleyer J.2019 FEniCS and MFront for complex non linear solid mechanics simulation. See http://rgdoi.net/10.13140/RG.2.2.35501.54247.
60. Ahrens J, Geveci B, Law C. 2005 ParaView: an end-user tool for large-data visualization. In Visualization handbook (eds CD Hansen, CR Johnson), ??. 717-731. Burlington, MA: Elsevier. (doi:10.1016/B978-012387582-2/50038-1)
61. Ayachit U, Geveci ?, Moreland ?, Patchett J, Ahrens J. 2012 The paraview visualization application. Boca Raton, FL: Chapman and Hall/CRC. (doi: 10.1201/b12985-23)
62. VTK website. VTK - The Visualization Toolkitdata. See https://vtk.org (accessed 1 March 2025).
63. Lorensen WE, Cline HE. 1998 Marching cubes: a high resolution 3D surface construction algorithm. In Seminal graphics: pioneering efforts that shaped the field, pp. 347-353, vol. 1. New York, NY: Association for Computing Machinery. (doi:10.1145/280811.281026)
64. Hammer P. 2023 Marching cubes. See https://www.mathworks.com/matlabcentral/fileexchange/32506-marching-cubes.
65. Attene ?. 2010 A lightweight approach to repairing digitized polygon meshes. Vis. Comput. 26, 1393-1406. (doi:10.1007/s00371-010-0416-3)
66. Sullivan C, Kaszynski A. 2019 PyVista: 3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (VTK). J. Open Source Softw. 4, 1450. (doi:10.21105/j0ss.01450)
67. Dawson ? et al. trimesh. See https://trimesh.org/.
68. Jin D, lyer KS, Chen ?, Hoffman EA, Saha PK. 2016 A robust and efficient curve skeletonization algorithm for tree-like objects using minimum cost paths. Pattern Recognit. Lett. 76, 32-40. (doi:10.1016/j.patrec.2015.04.002)
69. Lee TC, Kashyap RL, Chu CN. 1994 Building skeleton models via 3-D medial surface axis thinning algorithms. CVGIP 56, 462-478. (doi:10.1006/ cgip. 1994.1042)
70. Tuszynski J. 2023 Triangle/Ray Intersection. See https://www.mathworks.com/matlabcentral/fileexchange/33073-triangle-ray-intersection.
71. Chorin AJ. 1968 Numerical solution of the Navier-Stokes equations. Math. Comput. 22, 745-762. (doi:10.1090/s0025-5718-1968-0242392-2)
72. Témam ?. 1969 Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires (II). Arch. Ration. Mech. Anal. 33, 377-385. (doi:10.1007/BF00247696)
73. Goda K. 1979 A multistep technique with implicit difference schemes for calculating two- or three-dimensional cavity flows. J. Comput. Phys. 30, 76-95. (doi:10.1016/0021-9991(79)90088-3)
74. Valen-Sendstad K, Logg A, Mardal KA, Narayanan H, Mortensen M. 2012 A comparison of finite element schemes for the incompressible Navier- Stokes equations, pp. 409-430, vol. 84. Berlin, Germany: Springer Berlin Heidelberg. (doi:10.1007/978-3-642-23099-8_21)
75. Langtangen HP, Logg A. 2016 Solving PDEs in Python. Cham, Switzerland: Springer International Publishing. (doi:10.1007/978-3-319-52462-7)
76. Hsu ?, Secomb TW. 1989 A Green's function method for analysis of oxygen delivery to tissue by microvascular networks. Math. Biosci. 96, 61-78. (doi:10.1016/0025-5564(89)90083-7)
77. Cattaneo L, Zunino P. 2014 Computational models for fluid exchange between microcirculation and tissue interstitium. Networks Heterog. Media 9, 135-159. (doi:10.3934/nhm.2014.9.135)
78. Koppl T, Vidotto E, Wohlmuth ?, Zunino P. 2018 Mathematical modeling, analysis and numerical approximation of second-order elliptic problems with inclusions. Math. Model. Methods Appl. Sci. 28, 953-978. (doi:10.1142/s0218202518500252)
79. Laurino F, Zunino P. 2019 Derivation and analysis of coupled PDEs on manifolds with high dimensionality gap arising from topological model reduction. ESAIM 53, 2047-2080. (doi:10.1051/m2an/2019042)
80. Koch T. 2020 Benchmark (1.2 - numerical results reference solution. DaRUS. See https://darus.uni-stuttgart.de/citation?persistentld=doi:10. 18419/darus-471.
81. Secomb TW, Hsu ?, Park EYH, Dewhirst MW. 2004 Green's function methods for analysis of oxygen delivery to tissue by microvascular networks. Ann. Biomed. Eng. 32, 1519-1529. (doi:10.1114/b:abme.0000049036.08817.44)
82. Mortensen M, Valen-Sendstad ?. 2015 Oasis: a high-level/high-performance open source Navier-Stokes solver. Comput. Phys. Commun. 188, 177-188. (doi:10.1016/j.cpc.2014.10.026)
83. Pfaller MR, Pham J, Wilson NM, Parker DW, Marsden AL. 2021 On the periodicity of cardiovascular fluid dynamics simulations. Ann. Biomed. Eng. 49, 3574-3592. (doi:10.1007/s10439-021-02796-x)
84. Roache PJ, Ghia KN, White FM. 1986 Editorial policy statement on the control of numerical accuracy. J. Fluids Eng. 108, 2. (doi:10.1115/1.3242537)
85. Khan MO, Valen-Sendstad K, Steinman DA. 2015 Narrowing the expertise gap for predicting intracranial aneurysm hemodynamics: impact of solver numerics versus mesh and time-step resolution. Am. J. Neuroradiol. 36, 1310-1316. (doi:10.3174/ajnr.a4263)
86. Vignon-Clementel IE, Alberto Figueroa ?, Jansen KE, Taylor CA. 2006 Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries. Comput. Methods Appl. Mech. Eng. 195, 3776-3796. (doi:10.1016/j.cma.2005.04.014)
87. Kuchumov AG et al. 2022 Patient-specific 0D-3D modeling of blood flow in newborns to predict risks of complications after surgery. Health Risk Analysis 159-167. (doi:10.21668/health.risk/2022.4.15.eng)
88. CheikhSleiman. 2025 Tube2FEM. GitHub. https://github.com/CheikhSleiman/Tube2FEM
89. Sleiman HC. 2025 CheikhSleiman/Tube2FEM: v1.0.0. Zenodo. (doi:10.5281/zenodo.15441489)
© 2025. This work is published under https://creativecommons.org/licenses/by/4.0/ (the "License"). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.