1 Introduction

Mechanical ventilation is considered a lifesaving intervention in patients whose lungs prove inadequate to deliver oxygen due to a variety of medical conditions such as acute lung injury and chronic obstructive pulmonary disease [19]. It is a widely used therapeutic technique to deliver medicine in aerosolized form to treat the underlying lung conditions [19]. The effectiveness of ventilation and pulmonary drug delivery during mechanical ventilation is known to depend on several parameters [8] that are grouped as related to: (a) ventilator (ventilation mode, tidal volume, and inspiration waveform); (b) circuit (endotracheal tube—ETT size); (c) device (device type, nature of cycle); (d) drug (dosage, particle size, duration); and (e) patient (age, lifestyle, occupation). Airflow in the human airways has been studied for decades, and a comprehensive review can be found in recent article [16]. Many studies in the past have used the statistical models of Weibel [24] or Horsfield [13] and would typically model four or less generations. In recent years, due to advancements in medical imaging equipment and techniques, studies have begun to use geometries based on computed tomography (CT) scans of patients [5, 25].

Majority of these studies have investigated the influence of various patient-related factors on fluid flow and particle transport in the human airways [25]. There is little or no work on understanding the role of invasive conventional mechanical ventilation (CMV) on pulmonary fluid dynamics and drug delivery [8]. Any systemic study to account for CMV must incorporate ventilation parameters coupled with clinically relevant inlet conditions in order to be of use to clinicians for formulating effective ventilation management techniques. We focus on understanding the effect of ventilator-related parameters on the resulting pulmonary fluid dynamics. Transient simulations were performed on a CT scan-based computer model of the human airways to understand the effect of changes in ventilation waveform on flow structures inside the upper airways; representative waveforms used by commercial CMV in intensive care units are used for the purpose. Our discussion primarily focuses on the flow field in the trachea and in cross sections at main bronchi.

2 Methods

2.1 CT scan-based geometry

Computed tomography (CT) scans were taken from a 57-year-old (male) patient undergoing treatment at the University of Missouri hospital at Kansas City, using a Siemens AS 128 slice CT scanner. The field of view was approximately 34 cm and covered the entire lung. The resolution of each scan was 512 × 512 pixels with a pixel size of about 0.7 × 0.7 mm2 and a slice increment of 3 mm. Images were taken at a pixel depth of 12-bit grayscale. The DICOM files were imported into Mimics® software (Materialize Inc.) for volumetric reconstruction and segmentation. A threshold operation was performed to extract a set of voxels with grayscale values in a selected range of Hounsfield units. A region-growing operation was performed on the set of voxels that had undergone the threshold operation to obtain a contiguous set of voxels. Segmentation was performed on the reconstructed volume to remove under-resolved branches of the lung model. A set of faces was selected for each outlet and was projected onto a plane at the end of the airway perpendicular to the airway the outlet belonged to. To ensure appropriate outlet conditions and prevent backflow, extensions were added to each outlet by extruding the projected faces parallel to the airway and gradually transitioning the airway shape to a circle. A similar process was performed at the inlet of the endotracheal tube to ensure an appropriate inlet condition. The completed geometry is shown in Fig. 1a. The model used in this study includes the trachea and four-to-seven generations. There are a total of 41 outlets with generation numbers ranging from G4 to G7. Due to the nature of the reconstruction technique, the lengths, sizes, branching angles, branch orientations, and shapes of the airways are consistent with the actual geometry of the patient. From the ventilator circuit, the endotracheal tube has been included and the fluid flow from the ventilator into the ETT has been modeled with an inlet condition.

Fig. 1
figure 1

a Upper trachea-bronchial model obtained from CT scans. Model includes endotracheal tube (ETT) inside the trachea (G0). A total of 7 generations and 41 outlets are used. bd Conventional mechanical ventilator waveforms used for study flow transport and mixing in the lungs. e Flowchart detailing steps for inlet velocity generation

2.2 Numerical method and computational mesh

The working fluid for all simulations reported in this study is air at 25 °C (ρ = 1.185 kg/m3; μ = 1.831 × 10−5 kg/m s). At peak inspiratory velocity, the Reynolds number throughout the domain ranges from about 180 to 13,400 indicating that both laminar and turbulent conditions are present in the flow. To model turbulence, large eddy simulation (LES) was used. The filtered conservation of mass and momentum equations for incompressible flow is [21]:

$$\frac{{\partial \bar{u}_{i} }}{{\partial x_{i} }} = 0;\quad \frac{{\partial \bar{u}_{i} }}{\partial t} + \frac{\partial }{{\partial x_{j} }}(\bar{u}_{i} \bar{u}_{j} ) = - \frac{1}{\rho }\frac{{\partial \bar{p}}}{{\partial x_{i} }} + \frac{\partial }{{\partial x_{j} }}\left[ {\nu \left( {\frac{{\partial \bar{u}_{i} }}{{\partial x_{j} }} + \frac{{\partial \bar{u}_{j} }}{{\partial x_{i} }}} \right)} \right] - \frac{{\partial \tau_{ij} }}{{\partial x_{j} }}$$
(1)

where \(\bar{u}_{i}\) is the ith filtered velocity component and x i is the ith position component, ρ is the fluid density, p is the pressure, ν is the fluid kinematic viscosity, and τij are the subgrid-scale (SGS) stresses, which represents the effects of unresolved small-scale fluid motions. The relationship between the subgrid-scale stress and the subgrid-scale eddy viscosity (ν sgs ) is expressed as [22]:

$$\left( {\tau_{ij} - \frac{{\delta_{ij} }}{3}\tau_{kk} } \right) = 2\upsilon_{sgs} \bar{S}_{ij} \;{\text{where,}}\,\bar{S}_{ij} = \frac{1}{2}\left( {\frac{{\partial \bar{u}_{i} }}{{\partial x_{j} }} + \frac{{\partial \bar{u}_{j} }}{{\partial x_{i} }}} \right)$$
(2)

It was chosen for its ability to reproduce the laminar to turbulent transition process. The constant C w was taken to be 0.5 [21]. A second-order central differencing scheme was used for advection, and a second-order backward Euler scheme was used for the transient term. A residual target of 10−4 was set for convergence at each time step than 10−4 [3]. The computational mesh was created using ANSYS® ICEM® meshing software. Vinchurkar and Longest [23] reported that structured hexahedral meshes are prohibitively expensive in terms of time even for simple geometries. Thus, in case of a very complex models such as CT scan-based geometries, unstructured tetrahedral mesh is appropriate and is almost always considered [9]. The geometry of the current study is based of CT scans and is complex for meshing purposes, so a fully tetrahedral mesh was adopted. The mesh near the wall was refined to properly resolve the boundary layer. The average value of y+ normal to the wall ranged from 1.03 to 2.2 for all the four considered waveforms; y+ ranged from 0.5 to 3 for LES has been reported in the literature [1, 7, 10]. Mesh sensitivity was evaluated by calculating the percentage of resolved turbulent kinetic energy (\({\text{TKE}}_{\text{R}}\)) as follows [10, 22]:

$${\text{TKE}}_{\text{R}} \,\% = \frac{{{\text{TKE}}_{\text{R}} }}{{{\text{TKE}}_{\text{R}} + {\text{TKR}}_{\text{UR}} }} \times 100,$$

where TKEUR is the unresolved TKE

A mesh that resolved >80 % of the total energy was considered to be appropriate as it requires moderate computational time [2, 10, 22]. Three mesh sizes considered for the current mesh sensitivity study were—coarse (\(0.5 \times 10^{6}\)), fine (\(1.1 \times 10^{6}\)), and very fine (\(2.5 \times 10^{6}\)). TKE R % values were computed at a vertical plane cutting through the ETT, the trachea, and the main bronchi, where the strong ETT jet release and impingement occur resulting in a highest TKE. The percentage-resolved TKE for the three mesh sizes was 69 % for the coarse mesh, 84 % for the fine mesh, and 92 % for the very fine mesh. Thus, a fine mesh with 1.11 × 106 elements that resolves ~84 % of the TKE for the range of waveforms used in the present study was considered adequate.

The waveforms created for this study were governed by the equations of motion for the lungs given as [19]:

$$P = \frac{V}{C} + QR,\;Q = \frac{\hbox{d}V}{\hbox{d}t}$$
(3)

where V is the volume, P is the airway pressure, and Q is the volumetric flow rate. C and R are lung compliance and airway resistance, respectively, and defined as:

$$C = \frac{\Delta V}{\Delta P},\;R = \frac{{\Delta P}}{{\Delta Q}}$$
(4)

To generate the inspiratory waveforms, one of the three control variables (pressure, volume, or flow rate) was set as a simple shape as found in commercially available ventilators [19] and the other two were solved for numerically using Eq. 3. The tidal volume of the lung was approximated as 650 mL, and the lung properties of compliance and resistance were taken as 0.127 L/cm H2O and 4.33 cm H2O/(L/s), respectively [12], and were kept constant. The ventilator waveforms used for this study are shown in Fig. 1b–d. Two of the waveforms were chosen to replicate a pressure control ventilator; a constant pressure waveform (P-constant), and a sinusoidal (pressure) waveform (P-sine). The third waveform was selected to mimic a volume control ventilator and contained a volume ramp waveform (V-ramp). The fourth waveform was chosen to mimic a ventilator that controls the flow rate into the lungs and is of an ascending ramp type (F-ascending ramp). The sinusoidal pressure waveform had an inspiratory-to-expiratory ratio (I:E) of unity to mimic normal breathing [17]. All other waveforms were set to an I:E ratio of 0.33. Since the expiration is passive during CMV, for all these cases, the expiration is in the shape of an exponential decay, modeling a natural outward breath based on the elastic properties of the lung [19]. Table 1 shows considered cases and the flow parameters considered in current study.

Table 1 CMV parameters for the various waveforms and quantification of ETT effects

2.3 Boundary conditions

The velocity profile at the inlet was specified as a fully developed turbulent velocity profile described by [4]:

$$u(r,t) = 1.224u_{m} (t)\left( {\frac{R - r}{R}} \right)^{{\tfrac{1}{7}}}$$
(5)

where r is the radial distance from the origin, R is the radius of the inlet, and \(u_{m} (t)\) is the time-varying mean velocity which is calculated, based on each waveform shape, prior to CFD simulation as: \(u_{m} (t) = Q(t)/A_{ETT}\), where Q(t) is the flow rate curve obtained by using Eq. (3), see Fig. 1e, and A ETT is the ETT inlet area. The obtained transient mean velocity \(u_{m} (t)\) is then used in Eq. (5) as inlet velocity profile in the CFD run [3]. The walls of the domain were assumed to be smooth, dry, and rigid imposed [16]. A no-slip boundary condition was applied to the walls of the domain. The region surrounding the ETT represents ETT cuff and was set as a wall to prevent flow from entering or escaping the domain [19]. Due to difficulty of obtaining physiologically relevant boundary conditions, the flow at the outlet was specified with mass flow rates defined as fractions of instantaneous fractions of mass flow rates at the inlet of the domain based on date from Horsfield [13]. This is a standard technique that is used in pulmonary CFD studies [3, 25]. Details of outlet boundary conditions specification can be found elsewhere [3].

For validation purposes, a single bifurcation computational domain was created to match the experimental model of Lieber and Zhao [17]. Transient simulations were run with an oscillatory flow pattern that matched the experimental conditions. A Womersley solution to oscillatory flow in a straight tube of diameter D was imposed as a velocity condition at the inlet as [17]:

$$u(r,t) = {\text{Real}}\left[ {\frac{1}{i\rho \omega }\left( {1 - \frac{{J_{0} (\alpha ri^{{\tfrac{3}{2}}} )}}{{J_{0} (\alpha i^{{\tfrac{3}{2}}} )}}} \right)A^{*} e^{i\omega t} } \right]$$
(6)

where ρ is the fluid density, i is \(\sqrt { - 1}\), r is the non-dimensional radial distance, and J 0 is the zeroth Bessel function of the first kind. The parameter A * can be determined by assuming a peak flow rate of \(Q = \text{Re} D\nu /4\) at a cycle time where t/T = 0.25, and T is the period. \(\omega = 2\pi f\) is the angular velocity, where f is the ventilation frequency. Womersley number (\(\alpha = (D/2)\sqrt {\omega /\nu }\)) was 4.3, and the flow rate was 23.43 L/min. This velocity profile was implemented using a FORTRAN subroutine. Recorded velocity profiles from the simulation at time t = T/5 and 7T/10 at stations 2, 10, and 15 [17] (see Fig. 2) show good agreement with the experimental data.

Fig. 2
figure 2

Results from experimental validation using single bifurcation model (inset picture) of Lieber and Zhao [17]. Circles represent velocity measurements from [17]; lines represent computational velocity profiles. “B” indicates that the station is in the bifurcation plane, and “T” indicates it is in the transverse plane. Comparison presented for a t = 7T/10 (inspiration) in the bifurcation plane; b t = 7T/10 (inspiration) in the transverse plane; c t = T/5 (expiration) in the bifurcation plane; and d t = T/5 (expiration) in the transverse plane. Womersley number (\(\alpha = (D/2)\sqrt {\omega /\nu }\)) was 4.3, and the flow rate was 23.43 L/min

3 Results

3.1 Fluid flow structures

The three dominant flow features that were apparent throughout inspiration were: (1) the jet caused by the ETT (referred henceforth as ETT jet), (2) the large rotating structure in the right main bronchus, and (3) swirling flows in the bronchi. These flow structures are shown in Fig. 3a. The swirling flows in the bronchi, at times, resemble Dean vortices [6]; however, some of the swirling is thought to result from the non-circular geometry of the airways as well as flow curvature effects. The Dean number is defined as \(De = \text{Re} \sqrt {D/R_{c} }\), where D is the (local) airway segment diameter and R c is the radius of curvature [6]. De ranged from 715 to 1379 in the right bronchus and from 232 to 466 in the left bronchus for all considered waveforms. The presence and the magnitude of the ETT jet, the large rotating structure in the right main bronchus, and the swirling flows in the bronchi vary throughout the inspiratory phase. Figure 3b shows the early stages of the ETT jet. It is clear that the jet is pointed toward the right main bronchus due to the orientation of the ETT. An ETT jet pointed toward the right main bronchus is common because the path to the right main bronchus from the trachea is more direct for most patients [14]. After the ETT jet impinges on the wall of the right main bronchus, it begins to disperse as shown in Fig. 3c. A low grazing angle and the orientation of the surface where the jet impinges on the wall cause the jet to deflect to the right lung. The maximum velocity of this deflected jet is less than that of the ETT jet. Between t = 0.1–0.2 s, the deflected jet impinges in the lower part of the right main bronchus where the airway turns toward the anterior side of the patient, causing the deflected jet to disperse; however, there is no noticeable deflected (tertiary) jet. The dispersed (deflected) jet begins to flow back upward to the airways of the right upper lobe and to the left main bronchus. It is important to note that a large portion on the flow entering the left main bronchus does not come directly from the ETT but rather passes first through the right main bronchus before recirculating, passing over the carina of the first bifurcation, and entering the left main bronchus, and the dashed black arrow in Fig. 3a indicates the flow direction. The downward flow of the ETT jet and the upward flow of the dispersed jet are what make up the large rotating structure in the right main bronchus (shown in Fig. 3a). The ETT jet, the rotation in the right main bronchus, and the swirling pattern at slice L–L do not change in structure throughout the remainder of inspiration. Similar flow features were observed for all waveforms. The development of these features was similar for all waveforms. However, the rate at which they developed varied based on how quickly the flow waveform accelerated indicating that the flow structures depend upon the flow history as well as the inspiratory flow rate. During expiration, some swirling in the bronchi remained and the ETT caused a constriction in the flow as it leaves the domain, Fig. 3d.

Fig. 3
figure 3

Velocity vectors colored by velocity magnitude for F-ascending ramp waveform: a Fluid flow structures at t/T = 0.1 are labeled as: (1) ETT jet, (2) large rotating structure in the right main bronchus, (3) one example of swirling flows in the bronchi; b early development of the ETT jet at t/T = 0.01875; c impingement of the ETT jet at t/T = 0.025; d velocity contours for P-constant waveform at t/T = 0.275 (peak expiration)

3.2 Velocity distribution

Figure 4 shows the mean axial velocity distribution at different locations in the trachea and the main bronchi during peak inspiration. The mean axial velocity is defined as the time-averaged velocity estimated by averaging the axial velocity at each time step from the initial instant up to peak inspiration [20]. Locations T1 and T2 are located within a distance of 1.75DETT and 3.55DETT downstream from the ETT exit, while locations L and R are located in the left and right bronchus. The radial distance was normalized by the hydraulic diameter of each location; complex axial profiles were observed due to geometry features like ETT, curvatures, and asymmetry. The potential core region of the ETT jet is defined as the zone where the  mean velocity along the centerline of the jet remains constant and equal to the mean speed at the jet exit [18]. It ranged from axial distances of 3.14DETT to 3.46DETT (see Table 1). As observed, the jet forms as a nearly flat, fully developed turbulent velocity at location T1 (see Fig. 4a) and is similar in shape to the velocity waveform assumed at the ETT inlet. The ETT jet is associated with flow reversal in the areas surrounding the jet. However, as described earlier, flow reversal was also observed in the left side of the trachea due to the orientation of the ETT toward the right bronchus. Similar observations of flow reversal on the left side of the trachea with confinement of a laryngeal jet between the large recirculation zones were reported by Xi et al. [25]. Downstream of the potential core region, the jet centerline velocity begins to decay and develop into the Gaussian profile (location T2, see Fig. 4b). Further downstream at location R in the right bronchus (Fig. 4c), the velocity is observed to be skewed toward the inner wall due to the curvature. A large zone of flow reversal exists at the outer wall due to high-speed impingement of the ETT jet at the inner wall. At location L in the left bronchus (Fig. 4d), interestingly, the trend is reversed and the velocity profile is skewed toward the outer wall, while the reversal flow exists at the vicinity or the inner wall. This behavior is due to the observed flow reversal in the outer wall of the right bronchus, which causes a large rotating structure passes over the carina ridge toward the vicinity of the outer wall of the left bronchus (see Fig. 3a). It was noticed that the maximum mean velocity in the trachea was ~2–8 times higher than that in the right bronchus and left bronchus, respectively. This is merely due to the orientation of the ETT toward the right bronchus, which causes the ETT jet to impinge directly this airway segment. The velocity distribution at all four locations has similar profiles for all the waveforms considered in this study. However, the magnitude of the mean velocity varies depending on the flow rate. At all locations, P-sine has the lowest value of the mean velocities due the low flow rate. Although the P-constant waveform has relatively higher flow rate than the V-ramp and F-ascending ramp waveforms, it produces lower mean velocities. This is due to the time averaging over shorter period, since the P-constant reaches the peak inspiration rapidly and then decays exponentially.

Fig. 4
figure 4

Mean axial velocity distribution for the four waveforms (see Fig. 1) at different locations as indicated in the figure. Comparison of velocity profiles at a T1 at 1.75D ETT downstream; b T2 at 3.55D ETT downstream; c R at right bronchus; and d L at left bronchus. [DT1, DT2, DR, and DL are hydraulic diameters at locations T1, T2, R, and L, respectively]. The right y-axis is for the velocity magnitude for the case of no ETT

3.3 Secondary flow

To examine the secondary flow intensity, we use a secondary flow strength (SFS) parameter that is defined as the ratio of the radial velocity (V r ) to the mean velocity at a given location. The SFS was computed at each point on the cross-sectional slice and then averaged at each time step. The locations T at trachea, R at right bronchus, and L at left bronchus are shown in Fig. 1. The presence of the ETT jet dominates the flow, and all waveforms exhibit similar values of SFS; we thus choose to only present the variation in SFS at the four slices for the case of the P-constant waveform in Fig. 5a. Not surprisingly, a strong secondary motion exists at all locations during inspiration phase, due to high-speed ETT jet impinging at the inner wall of the right bronchus and geometry features such as curvatures and successive branching. SFS was found to be >0.6 over the inspiration period (Fig. 5a). At locations T and R, similar SFS behavior and magnitude were observed. At location L (inside left bronchus), SFS was consistently greater than the values at T and R during whole inspiration period with maximum difference of ~25 %. This is because the large portion of the flow in the left bronchus comes from the large swirling flow in the right bronchus, which cause higher radial velocity component compared to the axial velocity. During expiration phase, the secondary flow exists at locations T and R with moderate strength (i.e., 0.3 ≤ SFS ≤ 0.6) [3] and relatively higher at R due to effect of the curvature. In contrary, the secondary flow at location L has similar strength to that during inspiration period (i.e., strong SFS > 0.6) and is attributed to the high curvature and longer airway segment of left bronchus. De ranged from 672 to 1382 in the right bronchus and from 214 to 441 in the left bronchus for all considered waveforms, see Table 1.

Fig. 5
figure 5

a Temporal variation of the secondary flow strength (SFS) for the P-constant waveform. b Contour of TKE for P-constant at peak inspiration (t/T = 0.025). c Contours of vorticity for P-constant at peak inspiration (t/T = 0.025)

3.4 Turbulence

The turbulence inside trachea and main bronchi is characterized in terms of the TKE and vorticity contour at peak inspiration. For brevity, only data for the waveform (P-constant) that has maximum flow rate are shown in Fig. 5b and c. As the ETT jet impinges and deflects, the TKE begins to dissipate. The areas of highest TKE are those where the incoming ETT jet begins to interact with the deflected and dispersed flow coming upward out of the right main bronchus and across the path of the ETT jet. High TKE remains throughout most of the portion of the right lung and through a portion of the left main bronchus. In general, the TKE is stronger with the waveforms that have a higher maximum flow rate. Higher vorticity exists in the right lung due to the ETT jet indicating high turbulence, see Fig. 5c. Areas of high vorticity exist in the middle of the domain around the edges of the ETT jet and at the walls. Higher near-wall vorticity exists where higher near-wall velocities are present. One area of high vorticity exists where the ETT jet impinges and deflects off the wall of the right main bronchus. In the waveforms with higher peak flow rates (i.e., P-constant and F-ascending ramp), the vorticity suddenly drops near the bottom of the wall where the ETT jet is being deflected as jet detaches from the wall. A similar pattern occurs at the first bifurcation as the flow coming from the right lung passes over the carina. Another patch of high vorticity is observed at the top of the left main bronchus. During peak expiration, all waveforms exhibited similar vorticity contours. The P-sine had vorticity that was also similar in structure to other waveforms during inspiration and expiration phases, but had smaller values for vorticity. Vorticity was also present in the left main bronchus where there was a significant swirling motion due to geometry.

3.5 Wall shear stress

The occurrence of high wall shear stress (WSS) during airway ventilation is considered to be a risk to the wall epithelium layer as it erodes the layer resulting in damage and inflammation that leads to ventilator-induced lung injury (VILI). The ventilation through invasive tubes induces a high-speed jet impinging in the first bifurcation segment. Consequently, regions of high wall shear stress exist around the conducting airways. This phenomenon has been observed in many past studies of high-frequency ventilation and during coughing [3, 11]. To examine the influence of the waveforms of the invasive CMV on the wall shear stress, contours of WSS are plotted in Fig. 6 at peak inspiration. The WSS is defined as the force per unit area exerted by the wall on the fluid in a direction on the local tangent plane and is calculated based on \(\tau_{w} = 2\mu S_{ij}\), where \(S_{ij}\) is the strain rate tensor [15, 22]. Data are presented for the case of higher flow rate (P-constant) and lower flow rate (P-sine). When using P-constant waveform, areas of high WSS include the area surrounding the rotating structure in the right main bronchus, the carina of the first bifurcation, and the superior side of the left main bronchus where the flow from the right main bronchus enters the left lung. Not surprisingly, the areas of high WSS are shared with areas of high vorticity near the wall due to the high near-wall velocities. WSS is much higher in the right lung due to the effects of the ETT jet. F-ascending ramp waveform was observed to have relatively lesser patches of maximum WSS due to lower flow rate. Further reduction in the WSS occurred when using V-ramp, and only limited region of low WSS surrounding right bronchus and carina ridge was found. Pressure-controlled sinusoidal waveforms cause merely a minimal WSS in the right bronchus (Fig. 6). This a consequence of the lower flow rate when using sinusoidal waveform shape compared to other waveforms. It was noticed that beyond generation 2, only very low WSS exists and was located mainly around the carinas.

Fig. 6
figure 6

Contours of wall shear stress for the all four waveforms

Figure 7 shows the variation of the WSS over the ventilation cycle at a point P that was subjected to WSS during the ventilation process. The evolution of WSS was observed to follow the flow rate profiles for all cases expect for the V-ramp waveform. Although the V-ramp waveform produces a constant flow rate, the WSS dropped ~40 % after half of the inspiration time. This is attributed to the complex behavior of the flow structure inside the right bronchus segment. The P-constant and P-sine waveforms lead to highest and lowest flow rate, respectively; consequently, the P-constant waveform led to highest WSS, while P-sine waveform induced the lowest with difference of ~75 %. A drastic reduction in WSS was observed during expiration period due to the absence of the ETT jet.

Fig. 7
figure 7

Temporal evolution of wall shear stress (WSS) for all waveforms at point P which has maximum WSS value

4 Discussion

The dominant ETT jet similar to a laryngeal jet reported previously in the literature [4, 18] was observed. TKE in the human airways downstream of a laryngeal jet has been reported by Lin et al. [18]; however, the TKE value is attenuated by factor of 2 at the lower region of the trachea and further diminishes at the lower generations. The ETT jet was released far below the larynx in the trachea and, in this study, was aligned with the opening of the right main bronchus. As a result, the TKE remains strong much lower into the airways when compared with the TKE caused by the laryngeal jet. The eddy viscosity value was found as ~78 % and 45 % of that at the trachea in the first and second generations, respectively, and is found to remain throughout the domain with subsequent reduction in its magnitude relative to that in trachea (~16 % at G3 and ~6.5 % at G6). These flow features have a major role in the particle deposition (the subject of part 2). The ETT jet and the presence of strong turbulence enhanced the deposition and aerosol concentration on the upper generations, particularly, at the jet impingement site at the right main bronchus.

A simulation without the presence of intubation was run to quantify the effects of ETT. Table 1 presents a comparison between the magnitude of different important flow field variables in case of presence and absence of ETT. It was found that De decreased by 28 % to 62 % for all waveforms when ETT was not used. This is attributed to the decrease in the flow rate, which results in lower Re. In addition, the TKE was found to be 4–13 times greater than the case there is no ETT. Furthermore, the wall shear stress ratio (WSSR = WSSETT/WSSnoETT) ranged from 6 to 23 for all used waveforms, indicating a high risk of lung injury due to intubation.

5 Conclusions

LES was used to evaluate fluid flow characteristics for flow in the CT scan-based upper tracheobronchial region under invasive CMV conditions. A range of inspiratory waveforms that typically used to artificially ventilate patients undergoing invasive CMV treatment was generated and used in conjunction with the laws of motion described in Eq. 3. The following main conclusions can be drawn from the present work:

  • The jet caused by the ETT was an important fluid flow feature with effects stronger than the laryngeal jet. These effects include: high flow rates in the trachea, the creation of a zone of recirculation, the initiation of turbulence, and the impaction of particles at the site of the jet’s impingement.

  • Unique velocity profile was observed inside left bronchus airway segment due to presence and orientation of the intubation. The velocity skewed toward the outer wall as opposed to the case of no intubation and all existing observations of flow inside tracheobronchial airways in the literature.

  • The presence and strength of the recirculation zone in the right main bronchus and the swirling flows in the bronchi were found to be dependent upon the flow history as well as the inspiratory flow rate and led to higher deposition at inner wall of the right main bronchus. This implies that these transient fluid flow events may not be captured with steady-state analyses.

  • Strong turbulence was found in the wake of the jet caused by the flow from the endotracheal tube, which dissipated in the lower airways of the domain. The turbulence remains strong much lower into the airways when compared with the turbulence caused by the laryngeal jet. TKE was found to be 4–13 times more than the case when there is no ETT.

  • The ETT jet caused high WSS surrounding the main bronchi, and it was found that using pressure-controlled sinusoidal waveform reduces the wall shear stress by 75 %. Furthermore, the presence of intubation elevated the WSS ~6- to 23-fold, increasing the risk of airways epithelia layer injury.