Computational fluid dynamics benchmark dataset of airflow in tracheas

Computational Fluid Dynamics (CFD) is fast becoming a useful tool to aid clinicians in pre-surgical planning through the ability to provide information that could otherwise be extremely difficult if not impossible to obtain. However, in order to provide clinically relevant metrics, the accuracy of the computational method must be sufficiently high. There are many alternative methods employed in the process of performing CFD simulations within the airways, including different segmentation and meshing strategies, as well as alternative approaches to solving the Navier–Stokes equations. However, as in vivo validation of the simulated flow patterns within the airways is not possible, little exists in the way of validation of the various simulation techniques. The data presented here consists of very highly resolved flow data. The degree of resolution is compared to the highest necessary resolutions of the Kolmogorov length and time scales. Therefore this data is ideally suited to act as a benchmark case to which cheaper computational methods may be compared. A dataset and solution setup for one such more efficient method, large eddy simulation (LES), is also presented.


Value of Data
The tracheal geometry can act as a benchmark geometry for assessment of computational methods for airflow in the tracheal region.
The extracted data (locations of extraction given) and associated refinement statistics have been provided from the DNS simulation. Researchers can compare their methods precisely with the provided case at the same locations.
A detailed description of methods for calculating turbulent refinement statistics are outlined for future investigations.
UA provides a simple measure for the degree of non-uniformity in the flow and is particularly useful for flow in constricted and curved tubular geometries. Benchmark values of this metric are provided in this paper.

Data
The dataset includes: the tracheal geometry in STL format (Supplementary 1) and a CSV file of points that represent the centerline (Supplementary 2); extracted fluid mechanical metrics and turbulent statistics (including locations of extraction) from a simulation approaching the level of DNS. Additionally, all methods to calculate fluid mechanical metrics are outlined, including turbulent statistics (Section 4.2) and UA (Section 4.3). UA data has been provided for a number of idealised flow scenarios.

Reference and LES validation data
Two different tracheal airflow datasets are presented: a reference solution and a large eddy simulation. The reference solution is near DNS and is computationally expensive as it tries to resolve all scales, while LES is cheaper as it introduces a model that represents the effect of the unresolved scales on the resolved scales. As the reference solution is near DNS it can be used as benchmark data to assess the validity of a turbulence model, such as LES, and also to assess the degree of turbulence within the flow.
Analysis of convergence data was performed on one geometry (case B in Bates et al. [1] and shown in Fig. 1). The data was extracted from this geometry as it demonstrates large curvature and constriction. Hence represents a challenging benchmark from a flow computation point of view. A highly resolved simulation, with 9.2 million elements and a time step of 0.01 ms was performed as a reference case and was found to be approaching the resolution required for direct numerical simulation, as shown in Table 1. The Kolmogorov length and time scales were calculated as described in Section 4.2 and the dissipation values did not change significantly with refinement level, suggesting these values had converged. Therefore, the estimates of Kolmogorov scales were deemed sufficiently accurate. The cube root of cell volume, 3 ffiffiffi ffi V p and the time step Δt were used as characteristic length and time scales respectively. The ratio 3 ffiffiffi ffi V p =η represents the mesh refinement level divided by the smallest length scale in the flow. The parameter Δt=τ η is the ratio between the time step and the smallest time scale it is necessary to resolve, neglecting convection, while UΔt=η shows how far such a feature can be convected during one time step as a ratio to the smallest length scale in the flow.
The highly resolved data was used as a reference case to which reduced resolution data could be compared. Simulation meshes were generated with numbers of elements ranging from 945,000 to 9.2 million and quasi-steady simulations were performed on each mesh with time steps of 1, 0.1 and 0.01 ms. These coarser simulations were run both without further turbulence modelling and with the LES model.
The parameters found to offer the best compromise between matching the benchmark data and computational expense were a mesh of 2.2 million elements, with five prism layers starting at 0.15 mm high and increasing with a geometric progression of 1.3, a time step of 0.1 ms and with LES turbulence modelling. This mesh is shown in Fig. 2. Several criteria were used to judge the data: convergence, such as mean and fluctuating (see definitions in Section 4.2); velocity along lines normal and tangential to the flow path, as shown in Fig. 3; overall mean pressure drop from the glottis to one diameter above the apex of the bifurcation at the carina; as well as the amplitude and frequency of pressure oscillations. Table 2 shows these values to be close to identical for the reference and LES simulations. Each LES simulation took between 2800 and 4000 core hours to model the required period of inhalation. This time makes running large data sets at different flow rates prohibitively expensive. Fig. 4 shows spectra of fluctuating velocity data at the same point in both the Reference and LES simulations. The spectra are calculated using the method described in Varghese et al. [6], but using just two Hann windows, due to the limited sample size. The spectra are then normalised as described   for the LES simulations was defined as where C w ¼ 0:544 and V is the cell volume.

Turbulent statistics
This section provides data on important turbulent statistics provided with the dataset.
The Kolmogorov length and time scales are both calculated from the dissipation, ε which is defined by Delafosse et al. [2] as where u 0 1 , u 0 2 and u 0 3 represent the components of fluctuating velocity and νthe kinematic viscosity. The Kolmogorov length scale, η is then calculated by and the Kolmogorov time scale, τ η by  To interpret turbulent quantities, averaging operations are utilised. For a specific flow variable (f ) the time average is defined by, where T is the time period over which the variable is sampled. Numerically this procedure is defined by f ðx; y; z; tÞ where N indicates the number of samples taken over the period T. Random turbulent fluctuations (f 0 ) are represented by the instantaneous deviation from this mean value defined by Using Eq. (7) the variance of these fluctuations is defined by Numerically Eq. (8) can be expressed by which represents the square of the root mean square value (f 0 rms ), which is commonly quoted in literature.
The turbulent kinetic energy is related to the velocity fluctuations by As an example of typical averaging windows, in Bates et al. [1] all computations were simulated for 0.25 s of steady inhalation. Of this period, the first 0.1 s was ignored, to allow for transients caused by the impulsive start to die away. This period was determined by analysis of measures such as overall pressure losses, point velocities and velocity variances and determining when a plateau has been achieved. All turbulent statistics data were calculated over the remaining 0.15 s, a period required for mean values to converge, sampled at every time step. Table 3 shows UA data for several types of fully developed flow which have analytical solutions, all calculated for circular ducts.