Innovatives Supercomputing in Deutschland
inSiDE • Vol. 10 No. 2 • Autumn 2012
current edition
about inSiDE
index  index prev  prev next  next

Simulation of Nasal Cavity Flows for Virtual Surgery Environments

The anatomy and the functionality of the human respiration system is well understood, whereas physical pro- cesses such as the details of the gas transport within the respiration cycle lack insight. The influence of the air flow on physiological functions like the sense of smell and taste is only one example which is of great interest. Furthermore, the impact of pathologies in the human nasal cavity and the lung are still not fully explained and are the topic of intensive scientific research in the fields of medical treatment, engineering and computer science. An impaired respiration capability can lead to a reduction of the respiratory efficacy and can also induce diminished olfac- tory and degustation capabilities. Such pathologies are in general caused by allergies, mal- or deformations of e.g. the nasal septum, where their impact on the flow field in the nose plays a major role in the analysis and under- standing of the patients complaints.

Figure 1: Virtual surgery environment for treatment of individual pathologies in rhinology. A surface of the nasal cavity is extracted from CT-images and is used for flow simulations to analyze the patient's complaints. A virtual surgery is performed and the resulting nasal cavity is analyzed and evaluated until a satisfactory treatment is found and a real surgery is carried out.

In addition, the nasal cavity plays an important role in the protection of the human lung. The large surface of the human nasal cavity allows the moistur- ization and heating of the air to attain optimal conditions for the lower airways. These functions can be impaired by e.g., congenital malformations or rhino- anaplasty. The nasal cavity is also re- sponsible for filtering out particles and aerosols from the inhaled air. Diesel aerosols and respirable dust particles are proofed to cause damage to the human lung and induce the growth of cancer. Particle deposition in the respiration tract is a topic of intensive research to further understand the impact of such pathologies. To analyze the functions of the nasal cavity it is almost impossible to do in-vivo investi- gations, although some parameters like the pressure and the temperature can be investigated with a rhinomanometry and temperature measurements in the pharynx. These aspects only scratch the surface of the knowledge, which is required to fully understand the physics of the human respiration and starve for methods of higher detail.

"Rhinomodel", an interdisci- plinary Research Project

Computer-based methods have been developed over the last decades for the analysis of fluid mechanical problems and are well established in air and space science. They are, however, not very common in medical fields due to the high complexity of human organs and the physical aspects of the respi- ration, i.e. the particle and gas trans- port, temperature and moisturization. Within the scope of several research projects funded by the German Research Foundation (DFG) the Institute of Aerodynamics of RWTH Aachen University (AIA) investigates the human respiration system experimentally with the help of Particle Image Velocimetry methods (PIV), a laser- and photogra- phy-based particle tracking system, and numerically by applying Computa- tional Fluid Dynamics methods (CFD).

Figure 2: Surface extraction pipeline of the human nasal cavity surface from CT-data. The CT-images are preprocessed by a sharpening filter before a segmentation takes place, which splits the cavity into volumes of tissue and air. The surface is extracted and smoothed to obtain an accurate geometry of the nasal cavity.

The flow field in the complex shape of the nasal cavity and the human lung is simulated on High Performance Com- puters (HPC) like the Cray XE6 Hermit System at the Höchstleistungsrechen- zentrum (HLRS) in Stuttgart. These projects are highly interdisciplinary as they also require the expertise of scientists from the medical and com- puter science field. Therefore, a strong cooperation exists with radiologists from the University Hospital in Aachen (UKA), rhinology experts from the Institute of Statistics and Medical Informatics Cologne (IMSIE) and the Virtual Reality (VR) group of the Computing Center of the RWTH Aachen University (RZ RWTH). A key aspect of this research is the relation between the physical properties of the flow and the comfort of the patient, which is raised in a patient survey.

The study cycles through an analysis procedure as shown in Fig. 1, where in a first step it is required to obtain the geometrical representation of the nasal cavity from Computer Tomography (CT) images obtained from the UKA. The CT-images contain density intensities stored in a three-dimensional mesh. The extraction process involves several extraction procedures itself and is depicted in Fig. 2. In a preprocessing step, it is sometimes required to apply an image enhancement filter to the CT- data. Unsharp boundaries between tissue and air complicate the extraction of the surface and hence need to be sharpened with a convolution filter, which applies a filter mask to each voxel involving the surrounding neigh- boring voxels. To detect the volume of interest, a seeded region growing algorithm recursively marks voxels depending on their density value as provided by the CT. A triangular representation is then generated from the segmentation by using the Marching Cubes algorithm, which produces a water-tight surface.

Figure 3: Surface of the nasal cavity based on CT-image data. A saggital slice of the CT-data shows the relative position of the extracted surface. Streamlines, colored by their velocity visualize the flow field.

The segmentation of the CT-data operates on binary data and hence the resulting surface contains unsmooth stairsteps. The surface is then smoothed by applying filter methods such as a windowed sinc function smoothing, which is a low-pass filter reducing the high fre- quent surface variations deflected by the stairsteps. Based on this surface (see Fig. 4) the pipeline enters an opti- mization cycle, where at the beginning a flow simulation with CFD methods is performed to obtain the raw fluid mechanical properties, which are then analyzed by surgeons from the UKA and engineers from the AIA. Based on the gained knowledge, the surgeons perform a virtual surgery in the virtual environment "CAVE" provided by the VR group of the RZ RWTH and edit the extracted surface by virtually ablating or remodeling tissue.

This results in a new triangular representation of the nasal cavity, which again enters the optimization cycle to detect improve- ments to the patient's comfort based on the fluid mechanical properties of the flow. This cycle is repeated until a satisfactory modification of the nasal cavity is obtained which can then be applied in a real surgery. This way, the patient's complaints can be explored and treated individually. The simulation of the complete respiratory cycle can be simulated on the Cray XE6 Hermit system at HLRS, depending on the ac- curacy of the solution, within hours and days, which makes this method being within reach for the application in real world clinical environments.

Figure 4: Surface of the nasal cavity extracted from CT-data. This patient has a warped septum and a swollen inferior turbinate on the right side.

From Billions of Cubes to a Simulation

The flow solver ZFS (Zonal Flow Solver) for the simulation of the respiration cycle is developed at the AIA and works on Cartesian grids. This kind of computational meshes allows a fully automatic grid generation, which is a great advantage over grid generation methods producing boundary fitted meshes. Such meshes have to be set up manually and their configuration for complicated geometries such as the human respiratory system would take several days and is hence not feasible. ZFS combines a parallel grid generator, a Finite Volume (FV) and a Lattice-Boltzmann Method (LBM) solver for the simulation of the fluid flow, a Lagrangian particle tracking and post- processing functionalities. As a base for the flow solvers the space needs to be split into small cubes to discretize the equations of fluid motion, i.e. the Navier-Stokes equations for the FV solver and the Boltzmann equation for the LBM solver, respectively.

Figure 5: Streamlines of the nasal cavity flow colored by the local flow velocity. The bending in front shows the flow characteristics in the region of the inferior turbinate. A mixing region is found in the pharynx, where the left and right cavity merge.

Figure 6: Local boundary refinement in the human nasal cavity. The mesh is refined based on the distance to the boundary. In the center, a frontal slice through the nasal cavity shows the location of the mesh detail.

The grid generator is efficiently parallelized to be executed on thousands of processors on HPC systems, allows a resolution of billions of cells per nasal cavity and takes only seconds to generate on the Cray XE6 Hermit system at HLRS. The parallelism is a major advantage over serial grid generators, which are strongly limited to the available memory on only one system and may take up to several days to generate a grid with a sufficient resolution for the simulations in the respiratory tract. Local grid refinement, as shown in Fig. 6, allows a high resolution in regions, where high gradients are expected, i.e. in wall-bounded shear layers and mixing zones and increase the accuracy of wall-shear stress calculations, which are of great importance in the analysis of nasal cavity flows. Local flow phe- nomena, such as shear layers also need a high resolution for an accurate prediction. Such refinement is obtained by dynamic solution adaptive refine- ment initiated during the simulation. For the simulation, the LBM solver is used since it is more efficient for the low Mach number flow as it appears in the respiration system. In the LBM, the Lattice-BGK equation, which is the discretized BGK-equation, a simplified form of the Boltzmann equation, is solved iteratively for so called particle probability density functions (PPDFs).

Figure 7: A mapping of the wall-shear stresses on the nasal cavity surface. Red regions indicate high shear stresses, which appear in convergent channels, generating higher pressure losses.

In contrast to the Navier-Stokes equa- tions, the Boltzmann equation is a statistical approach for the description of the rate of change of the number of particles in an infinitesimal small volume of fluid. The LBM solver contains two simple steps, the propagation and collision step, and can easily be parallel- ized on thousands of processors. The algorithm computes the new PPDFs in each cell for a certain number of dis- crete directions, i.e. for 14, 18 or 26 neighboring cells in three dimensions. After this collision step the information is exchanged between neighboring cells and is again used in the computation of the new PPDFs. The LBM is valid for quasi-incompressible flows and hence only in the low Mach numbers regime. This leads to a decoupling of the energy equation from the Navier-Stokes equa- tions. To additionally solve for the tem- perature field, a one-directional coupled scalar transport equation is integrated, in which the convection of the scalar is controlled by the velocity. With this procedure all primitive variables, veloc- ity, density, pressure and temperature can be obtained in each time step in every grid cell.

The nasal Cavity, a highly complex and individual Organ

Several aspects play a major role in specifying the quality of a nasal cavity. One of the properties governing the respiration efficacy is the static pres- sure loss resulting from the shape of the cavity. Warped septums and swol- len turbinates are typical diagnoses of genetic endowments and allergies, respectively. Within the DFG-funded project "Rhinomodel", simulations of fluid flows in the nasal cavity have shown to be responsible for an increased pressure loss in such pathologies. The formation of recirculation zones inhibit- ing the direct access to the pharynx is analyzed by considering the character- istics of streamlines as shown in Fig. 3 and Fig. 5. Here, the streamlines are colored by the velocity of the fluid and show strong curvatures in regions of high vorticity.

Figure 8: Simulation of the temperature distribution in the human nasal cavity. The inflow condition is set to 20° C and the wall tem- perature to 37° C body temperature. The inhaled air is heated up depending on the heating capability of the nasal cavity and reaches almost body temperature in the pharynx under optimal conditions.

This form of visualization also supports the investigation of the olfactory capability as fluid not directed along the olfactory organ suggests a diminished sense of smell since the olfactory receptors cannot be reached. The forces acting on the nasal cavity tissue generating high pressure losses are investigated by mapping the wall- shear stresses, which are highest in regions of strong velocity gradients near the wall, onto the geometry.

Fig. 7 shows an example of such a visuali- zation, where the stresses are highest (red) in highly converged channels and in zones of increased vorticity. To accurately compute these stresses, a high grid resolution is required, which is why local grid refinement plays a fun- damental role in the estimation of the wall shear-stresses. The temperature increase is important for attaining optimal conditions for the lower airways and is considered by investigating the passive scalar transport from the nostrils, with a boundary condition of 20°C ambient temperature and heated up by the 37°C tissue temperature by diffusion and convection, to the pharynx. Healthy cavities show an in- crease of the temperature to almost body temperature while pathological cavities often exhibit a reduced heating capability. An example of the tempera- ture distribution in a nasal cavity is shown in Fig. 8.

Another aspect not being verified evidently in literature is whether the flow in the nasal cavity is completely laminar. The type of flow is especially important for investigations in the lower airways, where at the tran- sition from the pharynx to the larynx accurate boundary conditions have to be set. Recent investigations at the AIA show that the flow can indeed become transitional, depending on the local shape of the geometry and the local properties of the flow. Fig. 11 shows a strong production of vortices after the unification of the nasal cavity channels near the pharynx. Such zones are potential candidates of transitional to turbulent flow and are currently under detailed investigation.

Figure 9: Speedup results for a stong scaling test
for a total grid size of 0.403 x 109 cells for the
LBGK D3Q19 algorithm on the Cray XE6 Hermit system at HLRS Stuttgart.
Figure 10: Speedup results for a weak scaling test
for a constant local grid size of 0.537 x 106 cells for
the LBGK D3Q19 algorithm on the Cray XE6 Hermit system at HLRS Stuttgart.

Demand for highly optimized Code and future Requirements

The solution algorithm for the simulation of the nasal cavity flow should produce results within a reasonable time, especially when it is used in a virtual surgery environment. Therefore, it is important to have a highly optimized code, which is implemented on HPC systems. In cooperation with computer science specialists from HLRS Stuttgart and CRAY, ZFS has recently been optimized to yield a higher performance on the High Performance Computing systems of HLRS, with a special focus on Cray XE6 Hermit. It is important to note, that the emphasis is not only placed on the reduction of the compu- tation time, but also on the reduction of the communication time.

Fig. 9 shows the results of a strong scaling experiment on Cray XE6 Hermit, where for a con- stant problem size of 0.403 x 109 cells the number of processes is continu- ously increased from 1,024 to 8,192 and the speedup is measured. The small decrease in speed for 8,192 pro- cesses is due to the increase of the communication overhead. Instead of keeping the overall problem size con- stant, the problem size per process can be kept constant which leads to a weak scaling, for which the results in Fig. 10 show an almost perfect scalability.

For this experiment a maxi- mum number of cells of 109 has been reached. Currently, the scalability to even higher processor numbers up to the complete Cray XE6 system is investigated and first tests show very promising results. Future research at the AIA includes investigations of the aeroelastic problem which occurs dur- ing sleep apnea and also snoring. In this case an interaction of the elastic tissue with the flow field occurs, where the force exerted by the flow leads to a deformation of the tissue which in turn impacts the flow. Such problems require a prediction of unsteady flow fields coupled with a finite element solver for the tissue structure and increases the problem size and the computational effort considerably. These simulations will therefore be ideal can- didates for the HPC systems of HLRS.

Figure 11: The velocity distribution in the right nasal cavity shows a strong production of vorticity where two nasal cavity channels merge and the flow is directed towards the pharynx.


The research for project "Rhinomodel" has been conducted under DFG research grant WE-2186/5. The parallel grid generator has been developed as part of the collaborative research center SFB 686. The financial support by the German Research Foundation (DFG) and the optimization support by the HLRS and CRAY is gratefully acknowledged.

• Andreas Lintermann
Chair of Fluid Mechanics and Institute of Aerodynamics RWTH Aachen University

top  top