INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN BIOMEDICAL ENGINEERING Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 Published online 15 February 2012 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/cnm.1490
A novel formulation for Neumann inflow boundary conditions in biomechanics Volker Gravemeier 1,2, * ,† , Andrew Comerford 2 , Lena Yoshihara 2 , Mahmoud Ismail 2 and Wolfgang A. Wall 2 1 Emmy
Noether Research Group on “Computational Multiscale Methods for Turbulent Combustion”, Technische Universität München, Boltzmannstr. 15, D-85747 Garching, Germany 2 Institute for Computational Mechanics, Technische Universität München, Boltzmannstr. 15, D-85747 Garching, Germany
SUMMARY Neumann boundary conditions prescribing the total momentum flux at inflow boundaries of biomechanical problems are proposed in this study. This approach enables the simultaneous application of velocity/flow rate and pressure curves at inflow boundaries. As the basic numerical method, a residual-based variational multiscale (or stabilized) finite element method is presented. The focus of the numerical examples in this work is on respiratory flows with complete flow reversals. However, the proposed formulation is just as well suited for cardiovascular flow problems with partial retrograde flow. Instabilities, which were reported for such problems in the literature, are resolved by the present approach without requiring the additional consideration of a Lagrange multiplier technique. The suitability of the approach is demonstrated for two respiratory flow examples, a rather simple tube and complex tracheobronchial airways (up to the fourth generation, segmented from end-expiratory CT images). For the latter example, the boundary conditions are generated from mechanical ventilation data obtained from an intensive care unit patient suffering from acute lung injury. For the tube, analytical pressure profiles can be replicated, and for the tracheobronchial airways, a correct distribution of the prescribed total momentum flux at the inflow boundary into velocity and pressure part is observed. Copyright © 2012 John Wiley & Sons, Ltd. Received 7 June 2011; Revised 7 September 2011; Accepted 9 November 2011 KEY WORDS:
respiratory mechanics; cardiovascular mechanics; Neumann inflow boundary conditions; residual-based variational multiscale finite element method
1. INTRODUCTION Computational fluid dynamics has become an essential tool in understanding various flow phenomena in biological systems, for example, in cardiovascular and respiratory biomechanics. Recent applications include, among others, the improvement of pulmonary drug delivery [1, 2], the formulation of protective ventilation strategies [3], the design of medical devices [4, 5], and surgical planning [6]. According to a reasonable modeling strategy or because of limited resources, respectively, the computational domain is usually restricted to a small part of the biological system, for example, a section of the bronchial tree or a carotid bifurcation. As a consequence, artificial inflow and outflow boundaries are introduced in the computational model. This poses a problem for the definition of physiologically reasonable flow and/or pressure conditions at these boundary sections. At the outlets, various types of boundary conditions have been applied in the literature, such as zero-traction or (non)uniform pressure boundary conditions [7–13], velocity or flow rate boundary
*Correspondence to: Volker Gravemeier, Emmy Noether Research Group on “Computational Multiscale Methods for Turbulent Combustion”, Technische Universität München, Boltzmannstr. 15, D-85747 Garching, Germany. † E-mail: [email protected] Copyright © 2012 John Wiley & Sons, Ltd.
NEUMANN INFLOW BOUNDARY CONDITIONS IN BIOMECHANICS
561
conditions (e.g., [9,14]), or coupling with reduced-dimensional models ( e.g., [15–19]). At the inlets, the velocity or flow rate has usually been prescribed in the form of Dirichlet boundary conditions (e.g., [8–14,18]). Alternatively, Neumann inlet conditions may be prescribed as, for example, in [20] (The difference of that formulation compared with the one proposed in this work will be detailed in the following text). However, this choice potentially leads to an unstable solution, as shown in [21]. The same problem is encountered at the outlets whenever flow partially or completely reverses, that is, the outlets become (partial) inlets. Partial retrograde flow usually arises because of flow recirculation induced by the bifurcating nature of the geometry. Complete reversal of the flow direction occurs, for example, in airflow simulations during the expiratory phase. In this case, Neumann boundary conditions are applied on complete inflow boundaries over an intentionally prescribed amount of time. And in view of the aforementioned potential instability, this can be considered the worst-case scenario in biomechanical flow problems. Despite the potential existence of this instability, many examples can be found in the literature using Neumann boundary conditions where the flow enters the domain (e.g., [8, 11, 13]). However, stable solutions sometimes simply have to be considered as ‘lucky strikes’, or their stability is due to undesirable artificial numerical effects, and it is unlikely that such a boundary condition setup works properly for all scenarios. Recently, a number of methodologies have been proposed to deal with the aforementioned problem in blood flow simulations. For instance, Formaggia et al. introduced a total pressure boundary condition for hemodynamic applications in [17]. That approach required a reformulation of the Navier–Stokes equations to have the normal total stress as natural boundary condition. The flux at the outlet was enforced using a Lagrange multiplier formulation. However, no sophisticated numerical examples were presented. Therefore, it is unclear whether all numerical difficulties associated with complex flows can be addressed in that way. An alternative method for large-scale blood flow simulations was proposed in [21]. In that approach, a conventional Neumann boundary condition was applied at the outlet boundary. In addition, the velocity profile was constrained by means of a Lagrange multiplier technique. That technique was also applied in combination with so-called defective boundary conditions (i.e., boundary conditions based on averaged values) in [22]. Recently, this idea of applying defective boundary conditions was revisited in conjunction with using Nitsche’s method in [23]. So far, none of the aforementioned methods have been utilized in respiratory flow simulations. In this study, a different type of Neumann inflow boundary conditions is proposed that (i) allows the simultaneous application of velocity/flow rate and pressure curves, (ii) resolves the instability at the outlet, and (iii) does not require the introduction of additional techniques such as Lagrange multipliers or Nitsche’s method. All subsequent developments are based on a formulation of boundary conditions proposed, for example, in [24] for convection–diffusion equations and later in [25] for the incompressible Navier–Stokes equations. Therein, the total momentum flux is prescribed at inflow boundaries and only the traction at outflow boundaries. If Neumann-type boundary conditions are applied, additional boundary terms will appear in a finite element formulation as used here, which usually do not need to be considered explicitly. These additional boundary terms appear in a different form (and at different boundaries), if either a conservative or a convective form of the convective part of the Navier–Stokes equations is considered. Hence, both forms are addressed in the following text. The inevitable necessity of including these additional boundary terms was recently observed for turbulent flow problems, in particular (e.g., [26–28]). The flow configurations considered in those studies, a diffuser, a backward-facing step, and a flow past a square-section cylinder, all feature (artificial) outflow boundaries for turbulent flows. Thus, it may likely occur that vortices leaving the computational domain provoke reverse flow in local sections of the outflow boundary. Again, the outflow boundary turns out to be an inflow boundary, temporarily and partially. If such a change is not adequately accounted for via the explicit addition of the aforementioned boundary terms, the computation will become unstable, eventually. A residual-based variational multiscale (or stabilized) finite element method is used in this study. However, there appears to be no reason that prevents the proposed Neumann inflow boundary conditions from also being used in combination with, for instance, inf–sup stable finite elements. A rather standard variational multiscale approach featuring stabilization terms such as streamline upwind Petrov–Galerkin, pressure-stabilizing Petrov–Galerkin, and grad–div terms (see, e.g., [29, 30] for Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
562
V. GRAVEMEIER ET AL.
details) is pursued here, because turbulence effects will not be particularly investigated in this study. In fact, the maximum Reynolds numbers of the flow configurations exemplarily investigated here are approximately 3000. For wall-bounded flows at such Reynolds numbers, significant low-Reynolds-number effects still have to be expected, at least in channel configurations (e.g., [31]). Thus, turbulence effects in the more complex flow domains considered in this work, if present at all, might also be assumed to be of rather weak impact. Moreover, the present (and similar) residual-based variational multiscale methods have been shown to be valid methods for the numerical simulation of laminar-to-turbulent flow, anyway; see, for example, [26] for turbulent incompressible flow, in particular, and [32] for an investigation of laminar, transitional, and turbulent variable-density flow at low Mach numbers. However, the reader is also referred to [33] for a recent investigation of turbulent flow in a respiratory flow configuration using an algebraic variational multiscale-multigrid method (AVM3 ) specifically developed for the numerical simulation of turbulent flow in complex geometries; see, for example, [34, 35] for details of that method. The outline of the remainder of this article is as follows. The numerical method will be presented in Section 3, after governing equations and boundary conditions will have been introduced in Section 2. Numerical results obtained from two example configurations will be provided in Section 4. Finally, conclusions from this study will be drawn in Section 5. 2. GOVERNING EQUATIONS AND BOUNDARY CONDITIONS Fluid motion in a domain in the form of incompressible flow, as assumed here, is mathematically described by the incompressible Navier–Stokes equations. They consist of the mass-conservation (or continuity) equation, r u D 0,
(1)
and the momentum-conservation equation, which may either be written in conservative form with respect to the convective term as @u C r .u ˝ u/ C rp 2r " .u/ D f @t
(2)
or in convective (or conventional) form as @u C u ru C rp 2r " .u/ D f. @t
(3)
In Equations (1)–(3), the velocity is denoted by u, the (kinematic) pressure by p (i.e., dynamic pressure pdyn divided by density ), and the (kinematic) viscosity by (i.e., dynamic viscosity divided by density ), which is assumed constant in the following. Furthermore, the rate-of-deformation tensor is given by " .u/ D
1 ru C ruT . 2
An initial condition for the velocity field u in is prescribed as follows: u .x, t D 0/ D u0 .x/ . The boundary to is distinguished in two respects. First, an inflow as well as an outflow part of the domain are defined. For this purpose, the normal component of the velocity vector on the boundary, un D u n, where n denotes the outward-pointing unit normal vector on the boundary, is used. With un , the inflow boundary is defined as follows: in .t / D fx 2 j un .x, t / 0g. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
NEUMANN INFLOW BOUNDARY CONDITIONS IN BIOMECHANICS
563
For a unique decomposition of to be achieved, the outflow boundary is defined as out .t / D in .t / , accordingly. Note that inflow and outflow boundaries depend on time as a result of their dependence on the velocity field. However, for simplicity of notation, this dependence on time will be omitted in the following. Second, a Dirichlet part D and a Neumann part N of the boundary are distinguished such that D \ N D ; and D [ N D . Both Dirichlet and Neumann boundaries may encompass inflow as well as outflow sections of the boundary: in=out D D D \ in=out ,
in=out N D N \ in=out .
Dirichlet boundary conditions are prescribed on D as u D uD . In the numerical examples in the following text, Dirichlet boundary conditions are exclusively used for representing (nonmoving) walls, such that the usual no-slip condition uD D 0 is applied. Neumann boundary conditions are prescribed differently on inflow and outflow boundaries: the in total momentum flux is prescribed on N , un u pn C 2" .u/ n D hin ,
(4)
pn C 2" .u/ n D hout .
(5)
out , but only the traction on N
In a unified form, the Neumann boundary conditions on the complete Neumann boundary N may be written as follows: uin n u pn C 2" .u/ n D h,
(6)
where hD
hin hout
in on N out . on N
In Equation (6), a normal velocity uin n D
un jun j 2
in is introduced, which is equal to un on N and zero elsewhere. Accordingly, a normal velocity
uout n D
un C jun j 2
out and zero elsewhere. may be defined, which is equal to un on N
3. NUMERICAL METHOD With a discretization of the domain via nel finite elements e with characteristic element length h, where e D 1, ..., nel , a residual-based variational multiscale finite element method is formulated as follows: find p h 2 Sph and uh 2 Suh such that
h
q ,r u
h
nel X C rq h , Mp RhM
Copyright © 2012 John Wiley & Sons, Ltd.
eD1
e
D 0 8 q h 2 Vph
(7)
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
564
V. GRAVEMEIER ET AL.
and either in conservative form ! h h @u rvh , uh ˝ uh r vh , p h C " vh , 2" uh v , @t
C v
h
, uhn uh
nel X h h C , R r v C C out
N
eD1
nel X C uh rvh , Mu RhM
e
eD1
e
D vh , f C vh , h
N
8 vh 2 Vuh
(8)
or in convective form ! h h @u C vh , uh ruh r vh , p h C " vh , 2" uh v , @t
vh , uhn uh
C in
N
nel X r vh , C RhC
D vh , f C vh , h
e
eD1 N
C
nel X uh rvh , Mu RhM eD1
e
(9)
8 vh 2 Vuh .
In Equations (7)–(9), the weighting functions q h and vh are introduced for pressure and velocity, respectively. Corresponding weighting and solution function spaces are denoted by Vph , Vuh , Sph , and Suh . Furthermore, L2 -inner products ., / over various domains and boundaries, respectively, constitute the individual terms in Equations (7)–(9). For general considerations on finite element methods for flow problems, the reader is referred to [36]. The crucial term, which needs to be added to account for the (normal)velocity part at inflow or outflow boundaries, appears as the first term in the second line of Equations (8) and (9), respectively. Note that the word ‘added’ is actually rather inappropriate in this context; it is used here merely because this term is usually not included, particularly in the convective formulation due to in the fact that there is usually no Neumann part of the inflow boundary N . It will be briefly demonstrated in the following that the respective term essentially arises as a result of integration-by-parts procedures and the choice of boundary conditions according to Equation (6). The conservative formulation, Equation (8), is integrated by parts with respect to the convective, pressure, and viscous terms. As usual, the flux vector h is introduced for the Neumann boundary terms arising in the course of this procedure, and it is shifted to the right-hand side. Consequently, a correct balance at the inflow boundary is achieved. However, at the outflow boundary, where only the traction is prescribed, the convective part of the Neumann boundary term is not prescribed. Hence, the aforementioned ‘additional’ term, which carries a positive sign and is only active at an outflow boundary, remains on the left-hand side. On the contrary, the convective formulation, Equation (9), is integrated by parts only with respect to the pressure and viscous term. Introducing the flux vector in h in this case results in an ‘additional’ term in the form of the convective part on N , which then remains on the left-hand side; this term carries a negative sign and is only active at an inflow boundary accordingly. It is noted that, for these additional terms, the normal component of the velocity field is evaluated at each Gaussian integration point of the respective boundary element. Hence, even elements exhibiting simultaneous inflow and outflow can be taken into account appropriately. At this point, the difference of the present formulation compared with the one in [20] should be emphasized. In [20], the traction subject to Equation (5) instead of the total momentum flux according to Equation (4) is prescribed, whereas the aforementioned term accounting for the (normal) velocity part at inflow or outflow boundaries does not appear in the formulation in [20]. However, because the finite element formulation of the continuity equation, given here in Equation (7), was integrated by parts in [20], similar boundary terms arise. Hence, as it is required for the present approach to prescribe a velocity profile to represent the velocity part in Equation (4), the authors of [20] also (need to) prescribe a velocity profile. However, they use an augmented Lagrangian method Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
565
NEUMANN INFLOW BOUNDARY CONDITIONS IN BIOMECHANICS
to weakly enforce a shape of the inlet velocity profile as previously proposed for cardiovascular retrograde outlet flows in [21]; such a constraint is not required for the present method. Besides the standard Galerkin term, a pressure-stabilizing Petrov–Galerkin term appears in Equation (7) (the second term on the left-hand side). Moreover, both a grad–div (or bulkviscosity or least-squares on incompressibility constraint) and a streamline upwind Petrov–Galerkin term are included in Equations (8) and (9), respectively (the last two terms on the left-hand side). For details on such residual-based multiscale/stabilization terms, the reader is referred to [29]. The multiscale/stabilization terms are based on the discrete residuals of continuity and momentum-conservation equation, RhC D r uh , @uh C uh ruh C rp h 2r " uh f. @t Note that, as usual, the discrete residual of the momentum conservation is formulated in convective form, irrespective of whether the conservative or the convective form is used for the standard Galerkin terms. Also as usual, the second derivatives contained in the viscous term within RhM are not taken into account when using linearly interpolated elements, as will be carried out in the numerical examples in Section 4. However, it is noted that those second derivatives might be reconstructed, for example, according to a procedure proposed in [37]. For the following definitions of the stabilization parameters, both a ‘static’ characteristic element length, 13 RhM D
hp D
6Ve
p 3 based on the element volume Ve , and a ‘dynamic’ characteristic element length, hu D
2 nen X i D1
rNi
uh kuh k
based on the shape functions Ni for all nen element nodes i, representing the stream length within the respective element, are introduced. The stabilization parameter C used in this study is a slightly modified version of the definition proposed in [38]: C D where
kuh khp C , 2
! mk kuh khp C D min ,1 . 2
The stabilization parameters Mp and Mu are based on a combined definition of the ones proposed in [39, 40]: Mp D where
1 fT t M1p
1 C
4 M2p mk h2 p
,
! mk kuh khp ,1 , M1p D max 2 ! 4fT t ,1 , M2p D max mk h2p
Copyright © 2012 John Wiley & Sons, Ltd.
Mu D
1 fT t M1u
1 C
4 M2u mk h2 u
,
! mk kuh khu M1u D max ,1 , 2 4fT t M2u D max , 1 . mk h2u Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
566
V. GRAVEMEIER ET AL.
The parameter mk is usually assumed to be 1/3 for linearly interpolated elements. Note that merely the definition of the characteristic element length used differentiates Mp and Mu . The parameter fT is related to the time-integration scheme. In this study, a generalized-˛ method as originally proposed for problems of fluid mechanics in [41] is applied. The particular variant of this time-integration scheme used here was recently presented in [28]. The parameter 1 is chosen as in [28, 41], that is, 1 D 1=2, resulting in ˛m D 5=6, ˛f D 2=3, and D 2=3. For this particular scheme, the parameter fT is obtained as fT D ˛f =˛m D 24=45. Alternatively, for a generalized trapezoidal rule (or one-step- time-integration scheme), for instance, this parameter would be computed as fT D . 4. NUMERICAL EXAMPLES Because we would like to demonstrate the suitability of the presented approach for complete flow reversal, we concentrate on numerical examples related to respiratory mechanics in this section. However, the methodology can also be used for partial retrograde flow, which commonly occurs in blood flow simulations. The following applies to both numerical examples: A convective formulation as in Equation (9) is implemented: the additional boundary term is
only active at inflow boundaries, as a result. Linearly interpolated finite elements are used exclusively. It could indeed be confirmed for both numerical examples that, without the inclusion of the cru-
cial term to account for the (normal) velocity part at the respective inflow boundaries (the first term in the second line of Equation (9)), as addressed in Section 3, all computations become unstable, eventually. 4.1. Tube The problem domain is represented by a tube with cross section according to x12 C x22 r 2 and radius r D 0.0075 m, resulting in a cross-sectional area A D 177.0 106 m2 . The length in x3 -direction is l3 D 0.15 m (i.e., from x3 D 0.0 m to x3 D 0.15 m). A constant temperature of the air (297 K) is assumed in , which results in a density D 1.173 kg/m3 and a kinematic viscosity D 15.7 106 m2 /s. A zero velocity field is prescribed as the initial condition: u0 D 0. No-slip boundary conditions are assumed at the hull of the tube. During the inspiration phase, the inflow boundary is located at x3 D 0.0 m and the outflow boundary at x3 D 0.15 m; the boundaries are reversed during the expiration phase. At the respective outflow boundaries, zero-traction conditions are prescribed: hout D 0. At the respective inflow boundaries, hin is prescribed according to Equation (4). A usual parabolic velocity profile in a tube is prescribed for the velocity part in hin subject to u3 D u3,c
x2 C x2 1 1 2 2 r
,
where the centerline velocity u3,c is computed as u3,c D 2Q=A on the basis of the prescribed flow rate Q. The maximum value of the prescribed flow rate is 0.6035 l/s, yielding a corresponding centerline velocity u3,c D 6.83 m/s and a mean velocity u3,m D 3.415 m/s. The Reynolds number Rem based on mean velocity and tube diameter is 3263, which is equivalent to the Reynolds number Rec based on centerline velocity and tube radius. It is emphasized that the velocity profile described previously is not imposed as a Dirichlet boundary condition but as one part of the Neumann boundary condition according to Equation (4), besides the pressure part, which will be described in the following paragraph. Furthermore, it is emphasized that the form of the velocity profile (parabolic in the present case) represents the only assumption required a priori for being able to prescribe the proposed Neumann boundary conditions. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
NEUMANN INFLOW BOUNDARY CONDITIONS IN BIOMECHANICS
567
For the pressure part, the analytical value for the (kinematic) pressure at the inflow boundary p in is prescribed. The momentum balance in a stationary state formulated in x3 -direction results in @p D @x3
@ 2 u3 @ 2 u3 C @2 x 1 @2 x 2
.
(10)
Thus, integrated along the tube length, a value p in D 4=r 2 u3,c l3 D 1.1438 m2 =s2 is obtained. The viscous part is omitted because of its minor significance in the present case. The tube is spatially discretized by 80 elements in x3 -direction. Furthermore, 24 elements are arranged along the diameter of the cross section, refined towards the hull (Figure 1). This discretization results in an overall number of 34,560 elements. A constant time-step length t D 0.004 s is used. In a first test case, relatively fast sinusoidal start-up and slow-down phases of 0.2 s, respectively, with a ‘holding’ phase of 1.1 s in between, where the maximum inflow value is kept constant, are prescribed both for inspiration and expiration. As a second test case, two sine cycles are computed, with the length of a sine cycle covering one inspiration and one expiration phase of equivalent length, respectively, assumed to be T D 3 s. The prescribed flow rate curves for both cases are depicted in Figure 2. The first test case is used for investigating whether the expected velocity and pressure profiles in a stationary state can be obtained via the present definition of boundary conditions. It is sufficient to consider the pressure profiles along the tube; they are extracted during the ‘holding’ phases both for inspiration and expiration. The pressure profiles for the inspiration and expiration phases, respectively, are depicted in Figure 3, along with the respective analytical profiles. It may be observed that, except for slight deviations at the respective ends of the tube, the analytical profiles are matched by the numerically obtained profiles. For the second test case, Figure 4 shows velocity results at
Figure 1. Spatial discretization of tube. 0.8 Q [l/s] 0.6 0.4 0.2 0.0
0
1
2
3
4
5
t [s] 6
-0.2 -0.4 -0.6 -0.8
Figure 2. Prescribed flow rate curves for tube; dashed line: curve for test case 1, solid line: curve for test case 2. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
568
V. GRAVEMEIER ET AL.
0.00
0
20
40
60
80
100 p [10-2m2/s2] 120
0.05
0.10
x3 [m] 0.15
Figure 3. Pressure profiles for tube obtained during ‘holding’ phases in test case 1; solid (orange) line, inspiration phase; dashed (orange) line, expiration phase; and dotted (black) lines, respective analytical profiles.
Figure 4. Velocity vectors on colored velocity magnitude distribution (red color indicates high velocity and blue color low velocity); from left to right, t = 0.75, 1.30, 1.50, 1.70, 2.25, 2.80, 3.00, and 3.20 s.
various instances in time. The depiction series starts with the peak of the inspiration phase, passes the flow reversal when turning from inspiration to expiration for the first time and the peak of the expiration phase, and concludes with an impression of the velocity field at 3.2 s within the second inspiration phase. 4.2. Tracheobronchial airways The second example considers flow and pressure in a realistic tracheobronchial geometry undergoing mechanical ventilation. With this example, we demonstrate the applicability of the boundary condition when flow and pressure curves are given, which are different in form. The geometry used in the simulations is segmented from end-expiratory CT images (slice thickness and pixel size 0.7344 mm) of a healthy patient using the commercially available segmentation software MIMICS (Materialise NV, Leuven, Belgium). The computational geometry is shown in Figure 5. The geometry begins at the level of the trachea and progresses up to the fourth generation. The domain Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
569
NEUMANN INFLOW BOUNDARY CONDITIONS IN BIOMECHANICS
Figure 5. Tracheobronchial geometry used in the simulations.
is discretized in space utilizing mixed (predominantly hexahedral) element types, giving a total of 785,033 elements and 623,478 nodes (2.5M degrees of freedom). In the near wall regions, boundary layer meshes are fitted to capture the relevant flow physics in these locations. A constant time-step length t D 0.001 s is used. As boundary conditions for the simulation, mechanical ventilation data obtained from an intensive care unit patient suffering from acute lung injury are applied. Figure 6 shows the pressure and flow at the level of the trachea for this patient. Note that in the present simulation, the pressure component of the traction is scaled down from 103 to 101 because of the extremely high traction the former pressure level would create on the boundary surfaces as a result of the low air density. For such a simulation, a very long start-up time would be required to arrive at this high pressure level. Remedies to this issue will be the subject of future research. It is emphasized that the absolute pressure drop between the inflow and outflow boundaries remains unchanged during scaling. For the flow and pressure at the five peripheral boundaries to be determined, a zero-dimensional (0D) airway model is generated using the same CT images as the ones used for the three-dimensional (3D) tracheobronchial geometry (cf. Figure 7). The 0D model was recently shown to successfully mimic the underlying physiology of a lung in [42], thus generating the correct pressure and flow rate values within the lung. The flow and pressure determination utilizes the following steps. The mechanical ventilation data are first applied to the 0D model. Flow and pressure data are then extracted from the 0D model at the same locations where the 3D model is trimmed. Finally, extracted flow and pressure values are applied to the 3D model as inflow and outflow boundaries according to Equations (4) and (5), respectively. For simplicity, constant velocity profiles are chosen at all inflow boundaries for both inspiration and expiration phases.
0.6
3000
0.4
2500
0.2
2000
p (Pa)
Q(l/s)
0 −0.2
1500 1000
−0.4 500
−0.6
0
−0.8 −1
0
0.5
1
1.5
2
2.5
time(s)
3
3.5
4
−500
0
0.5
1
1.5
2
2.5
3
3.5
4
time (s)
Figure 6. Ventilator-measured flow (left) and pressure (right) curves at the level of the trachea used in the simulations (note that a positive flow represents flow into the domain). Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
570
V. GRAVEMEIER ET AL.
Figure 8 shows pressure and volume-rendered velocity distributions in the tracheobronchial model. The results demonstrate clearly the difference in behavior of the pressure and velocity over the course of the cycle. In particular, we see that between t D 0.8 and t D 1.4 s, the velocity distribution is essentially unchanged, but the pressure continues to rise; if the pressure and velocity curves had not been simultaneously applied, the pressure would have remained constant. At t D 1.84 s
Figure 7. Patient-specific zero-dimensional (0D) airway model.
Figure 8. Pressure (top) and velocity magnitude (bottom) in the tracheobronchial region. From left to right: t D 0.8, 1.4 (maximum pressure), 1.6 and 1.84 s (maximum expiration). Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
571
NEUMANN INFLOW BOUNDARY CONDITIONS IN BIOMECHANICS 30
3 sim input
2
sim input
25
Pressure (Pa)
velocity (m/s)
1 0 −1 −2
15 10 5
−3
0
−4 −5
20
0
0.5
1
1.5
2
2.5
3
3.5
−5
4
0
0.5
1
1.5
time(s)
2.5
3
3.5
4
45
6 sim input
5 4
35
3
30
2 1 0 −1
sim input
40
Pressure (Pa)
velocity (m/s)
2
time(s)
25 20 15 10 5
−2 −3
0
−4
−5
0
0.5
1
1.5
2
2.5
3
3.5
4
0
0.5
time(s)
1
1.5
2
2.5
3
3.5
4
time(s)
Figure 9. Comparison of the simulated (sim) and ventilator derived (input) velocity and pressure. Mean velocity (left, top) and pressure (right, top) at the trachea inlet. Mean velocity (left, bottom) and pressure (right, bottom) at the right bottom outlet. Negative velocity represents flow out of the domain.
(maximum expiration), we observe reversal of the high and low pressure regions and the formation of a dominant flow, which comes from the left superior lobe. The nonsymmetrical nature of the flow distribution is due to the realistic boundary conditions obtained from the 0D lung model. Figure 9 details a comparison of the ventilator-derived pressure and velocity used to determine the prescribed Neumann boundary condition (cf. Equation (4) for inflow and Equation (5) for outflow) with the pressure and velocity obtained from the simulation results, referred to as ‘input’ and ‘sim’, respectively. Overall, we see excellent agreement, however clearly there are some discrepancies, particularly when the slew rates are high (cf. the difference in pressure and velocity at t1.5 s). This mismatch during these rapid acceleration phases is due to transient fluid inertia. In particular, the velocity cannot keep up with the rapidly changing pressure gradient; hence, distribution of the total momentum flux between the velocity and pressure terms is affected. 5. CONCLUSIONS A particular formulation of Neumann boundary conditions that prescribes the total momentum flux at inflow boundaries of biomechanical problems has been proposed. This approach enables the simultaneous application of velocity/flow rate and pressure curves, which may be required in problems of respiratory mechanics. One of the most prominent features of such respiratory configurations are complete flow reversals, that is, inflow boundaries usually become outflow boundaries (and vice versa) multiple times during the simulation. The suitability of the proposed boundary condition formulation has been demonstrated for two respiratory flow examples, a rather simple tube and tracheobronchial airways (up to the fourth generation, segmented from end-expiratory CT images). A residual-based variational multiscale (or stabilized) finite element method has been used for all numerical simulations. For the second complex example, the boundary conditions have been generated from mechanical ventilation data obtained from an intensive care unit patient suffering Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
572
V. GRAVEMEIER ET AL.
from acute lung injury. Among other things, it has been possible to replicate analytical pressure profiles in the first (tube) example as well as a correct distribution of the prescribed total momentum flux at the inflow boundary into velocity and pressure part in the second (tracheobronchial airways) example. Although the focus of the numerical examples in this study has been on respiratory flows with complete flow reversals, the proposed formulation is suited for cardiovascular flow problems with partial retrograde flow as well. Instabilities observed for such problems (e.g., [21]) are resolved by the present approach without requiring the additional consideration of a Lagrange multiplier technique. The appropriateness of the proposed method for situations where the flow partially reverses, was already demonstrated in [26–28] for various turbulent flow problems; vortices may leave the computational domain and provoke reverse flow in local sections of the respective outflow boundaries in such situations. In the future, these Neumann boundary conditions prescribing total momentum flux will be applied to cardiovascular and respiratory flows which are fully coupled to 0D models of the downstream vessels (see, for example, [15] and [19]). In those approaches, the flow is passed from the 3D domain to the 0D domain and the pressure is returned to the 3D domain during the subsequent time step. The approach presented in this paper allows the same method to be applied; however, the total traction rather than the pressure will be returned to the 3D domain.
ACKNOWLEDGEMENTS
Support from the Emmy Noether Program (first author) and the Protective Artificial Respiration (PAR) Program through projects WA 1521/6-2, WA 1521/9-2, and WA 1521/13-1 (second and third author) of the Deutsche Forschungsgemeinschaft (DFG) is gratefully acknowledged. The authors would like to thank Professor J. Guttmann, Department of Anesthesiology and Intensive Care Medicine, University Hospital Freiburg, Germany, for kindly providing the mechanical ventilation data.
REFERENCES 1. Kleinstreuer C, Zhang Z, Li Z, Roberts WL, Rojas C. A new methodology for targeting drug-aerosols in the human respiratory system. International Journal of Heat and Mass Transfer 2008; 51:5528–5589. 2. Comerford A, Bauer G, Wall WA. Nanoparticle transport in a realistic model of the tracheobronchial region. International Journal for Numerical Methods in Biomedical Engineering 2010; 26:904–914. 3. Wall WA, Wiechert L, Comerford A, Rausch S. Towards a comprehensive computational model for the respiratory system. International Journal for Numerical Methods in Biomedical Engineering 2010; 26:807–827. 4. Figueroa CA, Taylor CA, Yeh V, Chiou AJ, Zarins CK. Effect of curvature on displacement forces acting on aortic endografts: a three dimensional computational analysis. Journal of Endovascular Therapy 2009; 16:284–294. 5. Probst M, Lülfesmann M, Nicolai M, Bücker M, Behr M, Bischof CH. Sensitivity of optimal shapes of artificial grafts with respect to flow parameters. Computer Methods in Applied Mechanics and Engineering 2010; 199:997–1005. 6. Taylor CA, Draney MT, Ku JP, Parker D, Steele BN, Wang K, Zarins CK. Predictive medicine: computational techniques in therapeutic decision-making. Computer Aided Surgery 1999; 4:231–247. 7. Stroud JS, Berger SA, Saloner D. Numerical analysis of flow through a severely stenotic carotid artery bifurcation. Journal of Biomechanical Engineering 2002; 124:9–20. 8. Zhang Z, Kleinstreuer C. Transient airflow structures and particle transport in a sequentially branching lung airway model. Physics of Fluids 2002; 14:862–880. 9. Nowak N, Kakade PP, Annapragada AV. Computational fluid dynamics simulation of airflow and aerosol deposition in human lungs. Annals of Biomedical Engineering 2003; 31:374–390. 10. Green AS. Modelling peak-flow wall shear stress in major airways of the lung. Journal of Biomechanics 2004; 37:661–667. 11. Kabilan S, Lin C-L, Hoffman EA. Characteristics of airflow in a CT-based ovine lung: a numerical study. Journal of Applied Physiology 2006; 102:1469–1482. 12. Freitas RK, Schröder W. Numerical investigation of the three-dimensional flow in a human lung model. Journal of Biomechanics 2008; 41:2446–2457. 13. Xia G, Tawhai MH, Hoffman EA, Lin C-L. Airway wall stiffening increases peak wall shear stress: a fluid-structure interaction study in rigid and compliant airways. Annals of Biomedical Engineering 2010; 38:1836–1853. 14. Calay RK, Kurujareon J, Holdö AE. Numerical simulation of respiratory flow patterns within human lung. Respiration Physiology 2002; 130:201–221. 15. Vignon-Clementel IE, Figueroa CA, Jansen KE, Taylor CA. Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries. Computer Methods in Applied Mechanics and Engineering 2006; 195:3776–3796. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
NEUMANN INFLOW BOUNDARY CONDITIONS IN BIOMECHANICS
573
16. Blanco PJ, Feijòo RA, Urquiza SA. A unified variational approach for coupling 3D–1D models and its blood flow applications. Computer Methods in Applied Mechanics and Engineering 2007; 196:4391–4410. 17. Formaggia L, Moura A, Nobile F. On the stability of the coupling of 3D and 1D fluid-structure interaction models for blood flow simulations. ESAIM Mathematical Modelling and Numerical Analysis 2007; 41:743–769. 18. Ma B, Lutchen KR. An anatomically based hybrid computational model of the human lung and its application to low frequency oscillatory mechanics. Annals of Biomedical Engineering 2006; 34:1691–1704. 19. Comerford A, Förster C, Wall WA. Structured tree impedance outflow boundary conditions for 3D lung simulations. Journal of Biomechanical Engineering 2010; 132:081002. 20. Kim HJ, Vignon-Clementel IE, Figueroa CA, LaDisa JF, Jansen KE, Feinstein JA, Taylor CA. On coupling a lumped parameter heart model and a three-dimensional finite element aorta model. Annals of Biomedical Engineering 2009; 37:2153–2169. 21. Kim HJ, Figueroa CA, Hughes TJR, Jansen KE, Taylor CA. Augmented Lagrangian method for constraining the shape of velocity profiles at outlet boundaries for three-dimensional finite element simulations of blood flow. Computer Methods in Applied Mechanics and Engineering 2009; 198:3551–3566. 22. Formaggia L, Gerbeau J-F, Nobile F, Quarteroni A. Numerical treatment of defective boundary conditions for the Navier–Stokes equations. SIAM Journal on Numerical Analysis 2002; 40:376–401. 23. Vergara C. Nitsche’s method for defective boundary value problems in incompressible fluid-dynamics. Journal of Scientific Computing 2011; 46:100–123. 24. Hughes TJR, Franca LP, Hulbert M. A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective–diffusive equations. Computer Methods in Applied Mechanics and Engineering 1989; 73:173–189. 25. Hughes TJR, Wells GN. Conservation properties for the Galerkin and stabilised forms of the advection–diffusion and incompressible Navier–Stokes equations. Computer Methods in Applied Mechanics and Engineering 2005; 194:1141–1159. 26. Bazilevs Y, Michler C, Calo VM, Hughes TJR. Isogeometric variational multiscale modeling of wall-bounded turbulent flows with weakly-enforced boundary conditions on unstretched meshes. Computer Methods in Applied Mechanics and Engineering 2010; 199:780–790. 27. Gravemeier V, Wall WA. An algebraic variational multiscale-multigrid method for large-eddy simulation of turbulent variable-density flow at low Mach number. Journal of Computational Physics 2010; 229:6047–6070. 28. Gravemeier V, Kronbichler M, Gee MW, Wall WA. An algebraic variational multiscale-multigrid method for large eddy simulation: generalized-˛ time integration, Fourier analysis and application to turbulent flow past a square-section cylinder. Computational Mechanics 2011; 47:217–233. 29. Hughes TJR, Scovazzi G, Franca LP. Multiscale and stabilized methods. In Encyclopedia of Computational Mechanics, Stein E, de Borst R, Hughes TJR (eds). Wiley: Chichester, 2004. 30. Braack M, Burman E, John V, Lube G. Stabilized finite element methods for the generalized Oseen problem. Computer Methods in Applied Mechanics and Engineering 2007; 196:853–866. 31. Moser RD, Kim J, Mansour NN. Direct numerical simulation of turbulent channel flow up to Re = 590. Physics of Fluids 1999; 11:943–945. 32. Gravemeier V, Wall WA. Residual-based variational multiscale methods for laminar, transitional and turbulent variable-density flow at low Mach number. International Journal for Numerical Methods in Fluids 2011; 65:1260–1278. 33. Comerford A, Gravemeier V, Wall WA. An algebraic variational multiscale-multigrid method for large-eddy simulation of turbulent pulsatile flows in complex geometries with detailed insight into pulmonary airway flow, 2011. Preprint, submitted for publication in International Journal for Numerical Methods in Fluids. 34. Gravemeier V, Gee MW, Wall WA. An algebraic variational multiscale-multigrid method based on plain aggregation for convection–diffusion problems. Computer Methods in Applied Mechanics and Engineering 2009; 198:3821–3835. 35. Gravemeier V, Gee MW, Kronbichler M, Wall WA. An algebraic variational multiscale-multigrid method for large eddy simulation of turbulent flow. Computer Methods in Applied Mechanics and Engineering 2010; 199:853–864. 36. Donea J, Huerta A. Finite Element Methods for Flow Problems. Wiley: Chichester, 2003. 37. Jansen KE, Collis SS, Whiting CH, Shakib F. A better consistency for low-order stabilized finite element methods. Computer Methods in Applied Mechanics and Engineering 1999; 174:153–170. 38. Franca LP, Frey SL. Stabilized finite element methods: II. The incompressible Navier–Stokes equations. Computer Methods in Applied Mechanics and Engineering 1992; 99:209–233. 39. Franca LP, Valentin F. On an improved unusual stabilized finite element method for the advective–reactive–diffusive equation. Computer Methods in Applied Mechanics and Engineering 2000; 190:1785–1800. 40. Barrenechea G, Valentin F. An unusual stabilized finite element method for a generalized Stokes problem. Numerische Mathematik 2002; 92:652–677. 41. Jansen KE, Whiting CH, Hulbert GM. A generalized-˛ method for integrating the filtered Navier–Stokes equations with a stabilized finite element method. Computer Methods in Applied Mechanics and Engineering 2000; 190:305–319. 42. Ismail M, Comerford A, Wall WA. Coupled and reduced dimensional modelling of respiratory mechanics during spontaneous breathing, 2011. Preprint, submitted for publication in Annals of Biomedical Engineering.
Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
A novel formulation for Neumann inflow boundary conditions in biomechanics Volker Gravemeier 1,2, * ,† , Andrew Comerford 2 , Lena Yoshihara 2 , Mahmoud Ismail 2 and Wolfgang A. Wall 2 1 Emmy
Noether Research Group on “Computational Multiscale Methods for Turbulent Combustion”, Technische Universität München, Boltzmannstr. 15, D-85747 Garching, Germany 2 Institute for Computational Mechanics, Technische Universität München, Boltzmannstr. 15, D-85747 Garching, Germany
SUMMARY Neumann boundary conditions prescribing the total momentum flux at inflow boundaries of biomechanical problems are proposed in this study. This approach enables the simultaneous application of velocity/flow rate and pressure curves at inflow boundaries. As the basic numerical method, a residual-based variational multiscale (or stabilized) finite element method is presented. The focus of the numerical examples in this work is on respiratory flows with complete flow reversals. However, the proposed formulation is just as well suited for cardiovascular flow problems with partial retrograde flow. Instabilities, which were reported for such problems in the literature, are resolved by the present approach without requiring the additional consideration of a Lagrange multiplier technique. The suitability of the approach is demonstrated for two respiratory flow examples, a rather simple tube and complex tracheobronchial airways (up to the fourth generation, segmented from end-expiratory CT images). For the latter example, the boundary conditions are generated from mechanical ventilation data obtained from an intensive care unit patient suffering from acute lung injury. For the tube, analytical pressure profiles can be replicated, and for the tracheobronchial airways, a correct distribution of the prescribed total momentum flux at the inflow boundary into velocity and pressure part is observed. Copyright © 2012 John Wiley & Sons, Ltd. Received 7 June 2011; Revised 7 September 2011; Accepted 9 November 2011 KEY WORDS:
respiratory mechanics; cardiovascular mechanics; Neumann inflow boundary conditions; residual-based variational multiscale finite element method
1. INTRODUCTION Computational fluid dynamics has become an essential tool in understanding various flow phenomena in biological systems, for example, in cardiovascular and respiratory biomechanics. Recent applications include, among others, the improvement of pulmonary drug delivery [1, 2], the formulation of protective ventilation strategies [3], the design of medical devices [4, 5], and surgical planning [6]. According to a reasonable modeling strategy or because of limited resources, respectively, the computational domain is usually restricted to a small part of the biological system, for example, a section of the bronchial tree or a carotid bifurcation. As a consequence, artificial inflow and outflow boundaries are introduced in the computational model. This poses a problem for the definition of physiologically reasonable flow and/or pressure conditions at these boundary sections. At the outlets, various types of boundary conditions have been applied in the literature, such as zero-traction or (non)uniform pressure boundary conditions [7–13], velocity or flow rate boundary
*Correspondence to: Volker Gravemeier, Emmy Noether Research Group on “Computational Multiscale Methods for Turbulent Combustion”, Technische Universität München, Boltzmannstr. 15, D-85747 Garching, Germany. † E-mail: [email protected] Copyright © 2012 John Wiley & Sons, Ltd.
NEUMANN INFLOW BOUNDARY CONDITIONS IN BIOMECHANICS
561
conditions (e.g., [9,14]), or coupling with reduced-dimensional models ( e.g., [15–19]). At the inlets, the velocity or flow rate has usually been prescribed in the form of Dirichlet boundary conditions (e.g., [8–14,18]). Alternatively, Neumann inlet conditions may be prescribed as, for example, in [20] (The difference of that formulation compared with the one proposed in this work will be detailed in the following text). However, this choice potentially leads to an unstable solution, as shown in [21]. The same problem is encountered at the outlets whenever flow partially or completely reverses, that is, the outlets become (partial) inlets. Partial retrograde flow usually arises because of flow recirculation induced by the bifurcating nature of the geometry. Complete reversal of the flow direction occurs, for example, in airflow simulations during the expiratory phase. In this case, Neumann boundary conditions are applied on complete inflow boundaries over an intentionally prescribed amount of time. And in view of the aforementioned potential instability, this can be considered the worst-case scenario in biomechanical flow problems. Despite the potential existence of this instability, many examples can be found in the literature using Neumann boundary conditions where the flow enters the domain (e.g., [8, 11, 13]). However, stable solutions sometimes simply have to be considered as ‘lucky strikes’, or their stability is due to undesirable artificial numerical effects, and it is unlikely that such a boundary condition setup works properly for all scenarios. Recently, a number of methodologies have been proposed to deal with the aforementioned problem in blood flow simulations. For instance, Formaggia et al. introduced a total pressure boundary condition for hemodynamic applications in [17]. That approach required a reformulation of the Navier–Stokes equations to have the normal total stress as natural boundary condition. The flux at the outlet was enforced using a Lagrange multiplier formulation. However, no sophisticated numerical examples were presented. Therefore, it is unclear whether all numerical difficulties associated with complex flows can be addressed in that way. An alternative method for large-scale blood flow simulations was proposed in [21]. In that approach, a conventional Neumann boundary condition was applied at the outlet boundary. In addition, the velocity profile was constrained by means of a Lagrange multiplier technique. That technique was also applied in combination with so-called defective boundary conditions (i.e., boundary conditions based on averaged values) in [22]. Recently, this idea of applying defective boundary conditions was revisited in conjunction with using Nitsche’s method in [23]. So far, none of the aforementioned methods have been utilized in respiratory flow simulations. In this study, a different type of Neumann inflow boundary conditions is proposed that (i) allows the simultaneous application of velocity/flow rate and pressure curves, (ii) resolves the instability at the outlet, and (iii) does not require the introduction of additional techniques such as Lagrange multipliers or Nitsche’s method. All subsequent developments are based on a formulation of boundary conditions proposed, for example, in [24] for convection–diffusion equations and later in [25] for the incompressible Navier–Stokes equations. Therein, the total momentum flux is prescribed at inflow boundaries and only the traction at outflow boundaries. If Neumann-type boundary conditions are applied, additional boundary terms will appear in a finite element formulation as used here, which usually do not need to be considered explicitly. These additional boundary terms appear in a different form (and at different boundaries), if either a conservative or a convective form of the convective part of the Navier–Stokes equations is considered. Hence, both forms are addressed in the following text. The inevitable necessity of including these additional boundary terms was recently observed for turbulent flow problems, in particular (e.g., [26–28]). The flow configurations considered in those studies, a diffuser, a backward-facing step, and a flow past a square-section cylinder, all feature (artificial) outflow boundaries for turbulent flows. Thus, it may likely occur that vortices leaving the computational domain provoke reverse flow in local sections of the outflow boundary. Again, the outflow boundary turns out to be an inflow boundary, temporarily and partially. If such a change is not adequately accounted for via the explicit addition of the aforementioned boundary terms, the computation will become unstable, eventually. A residual-based variational multiscale (or stabilized) finite element method is used in this study. However, there appears to be no reason that prevents the proposed Neumann inflow boundary conditions from also being used in combination with, for instance, inf–sup stable finite elements. A rather standard variational multiscale approach featuring stabilization terms such as streamline upwind Petrov–Galerkin, pressure-stabilizing Petrov–Galerkin, and grad–div terms (see, e.g., [29, 30] for Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
562
V. GRAVEMEIER ET AL.
details) is pursued here, because turbulence effects will not be particularly investigated in this study. In fact, the maximum Reynolds numbers of the flow configurations exemplarily investigated here are approximately 3000. For wall-bounded flows at such Reynolds numbers, significant low-Reynolds-number effects still have to be expected, at least in channel configurations (e.g., [31]). Thus, turbulence effects in the more complex flow domains considered in this work, if present at all, might also be assumed to be of rather weak impact. Moreover, the present (and similar) residual-based variational multiscale methods have been shown to be valid methods for the numerical simulation of laminar-to-turbulent flow, anyway; see, for example, [26] for turbulent incompressible flow, in particular, and [32] for an investigation of laminar, transitional, and turbulent variable-density flow at low Mach numbers. However, the reader is also referred to [33] for a recent investigation of turbulent flow in a respiratory flow configuration using an algebraic variational multiscale-multigrid method (AVM3 ) specifically developed for the numerical simulation of turbulent flow in complex geometries; see, for example, [34, 35] for details of that method. The outline of the remainder of this article is as follows. The numerical method will be presented in Section 3, after governing equations and boundary conditions will have been introduced in Section 2. Numerical results obtained from two example configurations will be provided in Section 4. Finally, conclusions from this study will be drawn in Section 5. 2. GOVERNING EQUATIONS AND BOUNDARY CONDITIONS Fluid motion in a domain in the form of incompressible flow, as assumed here, is mathematically described by the incompressible Navier–Stokes equations. They consist of the mass-conservation (or continuity) equation, r u D 0,
(1)
and the momentum-conservation equation, which may either be written in conservative form with respect to the convective term as @u C r .u ˝ u/ C rp 2r " .u/ D f @t
(2)
or in convective (or conventional) form as @u C u ru C rp 2r " .u/ D f. @t
(3)
In Equations (1)–(3), the velocity is denoted by u, the (kinematic) pressure by p (i.e., dynamic pressure pdyn divided by density ), and the (kinematic) viscosity by (i.e., dynamic viscosity divided by density ), which is assumed constant in the following. Furthermore, the rate-of-deformation tensor is given by " .u/ D
1 ru C ruT . 2
An initial condition for the velocity field u in is prescribed as follows: u .x, t D 0/ D u0 .x/ . The boundary to is distinguished in two respects. First, an inflow as well as an outflow part of the domain are defined. For this purpose, the normal component of the velocity vector on the boundary, un D u n, where n denotes the outward-pointing unit normal vector on the boundary, is used. With un , the inflow boundary is defined as follows: in .t / D fx 2 j un .x, t / 0g. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
NEUMANN INFLOW BOUNDARY CONDITIONS IN BIOMECHANICS
563
For a unique decomposition of to be achieved, the outflow boundary is defined as out .t / D in .t / , accordingly. Note that inflow and outflow boundaries depend on time as a result of their dependence on the velocity field. However, for simplicity of notation, this dependence on time will be omitted in the following. Second, a Dirichlet part D and a Neumann part N of the boundary are distinguished such that D \ N D ; and D [ N D . Both Dirichlet and Neumann boundaries may encompass inflow as well as outflow sections of the boundary: in=out D D D \ in=out ,
in=out N D N \ in=out .
Dirichlet boundary conditions are prescribed on D as u D uD . In the numerical examples in the following text, Dirichlet boundary conditions are exclusively used for representing (nonmoving) walls, such that the usual no-slip condition uD D 0 is applied. Neumann boundary conditions are prescribed differently on inflow and outflow boundaries: the in total momentum flux is prescribed on N , un u pn C 2" .u/ n D hin ,
(4)
pn C 2" .u/ n D hout .
(5)
out , but only the traction on N
In a unified form, the Neumann boundary conditions on the complete Neumann boundary N may be written as follows: uin n u pn C 2" .u/ n D h,
(6)
where hD
hin hout
in on N out . on N
In Equation (6), a normal velocity uin n D
un jun j 2
in is introduced, which is equal to un on N and zero elsewhere. Accordingly, a normal velocity
uout n D
un C jun j 2
out and zero elsewhere. may be defined, which is equal to un on N
3. NUMERICAL METHOD With a discretization of the domain via nel finite elements e with characteristic element length h, where e D 1, ..., nel , a residual-based variational multiscale finite element method is formulated as follows: find p h 2 Sph and uh 2 Suh such that
h
q ,r u
h
nel X C rq h , Mp RhM
Copyright © 2012 John Wiley & Sons, Ltd.
eD1
e
D 0 8 q h 2 Vph
(7)
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
564
V. GRAVEMEIER ET AL.
and either in conservative form ! h h @u rvh , uh ˝ uh r vh , p h C " vh , 2" uh v , @t
C v
h
, uhn uh
nel X h h C , R r v C C out
N
eD1
nel X C uh rvh , Mu RhM
e
eD1
e
D vh , f C vh , h
N
8 vh 2 Vuh
(8)
or in convective form ! h h @u C vh , uh ruh r vh , p h C " vh , 2" uh v , @t
vh , uhn uh
C in
N
nel X r vh , C RhC
D vh , f C vh , h
e
eD1 N
C
nel X uh rvh , Mu RhM eD1
e
(9)
8 vh 2 Vuh .
In Equations (7)–(9), the weighting functions q h and vh are introduced for pressure and velocity, respectively. Corresponding weighting and solution function spaces are denoted by Vph , Vuh , Sph , and Suh . Furthermore, L2 -inner products ., / over various domains and boundaries, respectively, constitute the individual terms in Equations (7)–(9). For general considerations on finite element methods for flow problems, the reader is referred to [36]. The crucial term, which needs to be added to account for the (normal)velocity part at inflow or outflow boundaries, appears as the first term in the second line of Equations (8) and (9), respectively. Note that the word ‘added’ is actually rather inappropriate in this context; it is used here merely because this term is usually not included, particularly in the convective formulation due to in the fact that there is usually no Neumann part of the inflow boundary N . It will be briefly demonstrated in the following that the respective term essentially arises as a result of integration-by-parts procedures and the choice of boundary conditions according to Equation (6). The conservative formulation, Equation (8), is integrated by parts with respect to the convective, pressure, and viscous terms. As usual, the flux vector h is introduced for the Neumann boundary terms arising in the course of this procedure, and it is shifted to the right-hand side. Consequently, a correct balance at the inflow boundary is achieved. However, at the outflow boundary, where only the traction is prescribed, the convective part of the Neumann boundary term is not prescribed. Hence, the aforementioned ‘additional’ term, which carries a positive sign and is only active at an outflow boundary, remains on the left-hand side. On the contrary, the convective formulation, Equation (9), is integrated by parts only with respect to the pressure and viscous term. Introducing the flux vector in h in this case results in an ‘additional’ term in the form of the convective part on N , which then remains on the left-hand side; this term carries a negative sign and is only active at an inflow boundary accordingly. It is noted that, for these additional terms, the normal component of the velocity field is evaluated at each Gaussian integration point of the respective boundary element. Hence, even elements exhibiting simultaneous inflow and outflow can be taken into account appropriately. At this point, the difference of the present formulation compared with the one in [20] should be emphasized. In [20], the traction subject to Equation (5) instead of the total momentum flux according to Equation (4) is prescribed, whereas the aforementioned term accounting for the (normal) velocity part at inflow or outflow boundaries does not appear in the formulation in [20]. However, because the finite element formulation of the continuity equation, given here in Equation (7), was integrated by parts in [20], similar boundary terms arise. Hence, as it is required for the present approach to prescribe a velocity profile to represent the velocity part in Equation (4), the authors of [20] also (need to) prescribe a velocity profile. However, they use an augmented Lagrangian method Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
565
NEUMANN INFLOW BOUNDARY CONDITIONS IN BIOMECHANICS
to weakly enforce a shape of the inlet velocity profile as previously proposed for cardiovascular retrograde outlet flows in [21]; such a constraint is not required for the present method. Besides the standard Galerkin term, a pressure-stabilizing Petrov–Galerkin term appears in Equation (7) (the second term on the left-hand side). Moreover, both a grad–div (or bulkviscosity or least-squares on incompressibility constraint) and a streamline upwind Petrov–Galerkin term are included in Equations (8) and (9), respectively (the last two terms on the left-hand side). For details on such residual-based multiscale/stabilization terms, the reader is referred to [29]. The multiscale/stabilization terms are based on the discrete residuals of continuity and momentum-conservation equation, RhC D r uh , @uh C uh ruh C rp h 2r " uh f. @t Note that, as usual, the discrete residual of the momentum conservation is formulated in convective form, irrespective of whether the conservative or the convective form is used for the standard Galerkin terms. Also as usual, the second derivatives contained in the viscous term within RhM are not taken into account when using linearly interpolated elements, as will be carried out in the numerical examples in Section 4. However, it is noted that those second derivatives might be reconstructed, for example, according to a procedure proposed in [37]. For the following definitions of the stabilization parameters, both a ‘static’ characteristic element length, 13 RhM D
hp D
6Ve
p 3 based on the element volume Ve , and a ‘dynamic’ characteristic element length, hu D
2 nen X i D1
rNi
uh kuh k
based on the shape functions Ni for all nen element nodes i, representing the stream length within the respective element, are introduced. The stabilization parameter C used in this study is a slightly modified version of the definition proposed in [38]: C D where
kuh khp C , 2
! mk kuh khp C D min ,1 . 2
The stabilization parameters Mp and Mu are based on a combined definition of the ones proposed in [39, 40]: Mp D where
1 fT t M1p
1 C
4 M2p mk h2 p
,
! mk kuh khp ,1 , M1p D max 2 ! 4fT t ,1 , M2p D max mk h2p
Copyright © 2012 John Wiley & Sons, Ltd.
Mu D
1 fT t M1u
1 C
4 M2u mk h2 u
,
! mk kuh khu M1u D max ,1 , 2 4fT t M2u D max , 1 . mk h2u Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
566
V. GRAVEMEIER ET AL.
The parameter mk is usually assumed to be 1/3 for linearly interpolated elements. Note that merely the definition of the characteristic element length used differentiates Mp and Mu . The parameter fT is related to the time-integration scheme. In this study, a generalized-˛ method as originally proposed for problems of fluid mechanics in [41] is applied. The particular variant of this time-integration scheme used here was recently presented in [28]. The parameter 1 is chosen as in [28, 41], that is, 1 D 1=2, resulting in ˛m D 5=6, ˛f D 2=3, and D 2=3. For this particular scheme, the parameter fT is obtained as fT D ˛f =˛m D 24=45. Alternatively, for a generalized trapezoidal rule (or one-step- time-integration scheme), for instance, this parameter would be computed as fT D . 4. NUMERICAL EXAMPLES Because we would like to demonstrate the suitability of the presented approach for complete flow reversal, we concentrate on numerical examples related to respiratory mechanics in this section. However, the methodology can also be used for partial retrograde flow, which commonly occurs in blood flow simulations. The following applies to both numerical examples: A convective formulation as in Equation (9) is implemented: the additional boundary term is
only active at inflow boundaries, as a result. Linearly interpolated finite elements are used exclusively. It could indeed be confirmed for both numerical examples that, without the inclusion of the cru-
cial term to account for the (normal) velocity part at the respective inflow boundaries (the first term in the second line of Equation (9)), as addressed in Section 3, all computations become unstable, eventually. 4.1. Tube The problem domain is represented by a tube with cross section according to x12 C x22 r 2 and radius r D 0.0075 m, resulting in a cross-sectional area A D 177.0 106 m2 . The length in x3 -direction is l3 D 0.15 m (i.e., from x3 D 0.0 m to x3 D 0.15 m). A constant temperature of the air (297 K) is assumed in , which results in a density D 1.173 kg/m3 and a kinematic viscosity D 15.7 106 m2 /s. A zero velocity field is prescribed as the initial condition: u0 D 0. No-slip boundary conditions are assumed at the hull of the tube. During the inspiration phase, the inflow boundary is located at x3 D 0.0 m and the outflow boundary at x3 D 0.15 m; the boundaries are reversed during the expiration phase. At the respective outflow boundaries, zero-traction conditions are prescribed: hout D 0. At the respective inflow boundaries, hin is prescribed according to Equation (4). A usual parabolic velocity profile in a tube is prescribed for the velocity part in hin subject to u3 D u3,c
x2 C x2 1 1 2 2 r
,
where the centerline velocity u3,c is computed as u3,c D 2Q=A on the basis of the prescribed flow rate Q. The maximum value of the prescribed flow rate is 0.6035 l/s, yielding a corresponding centerline velocity u3,c D 6.83 m/s and a mean velocity u3,m D 3.415 m/s. The Reynolds number Rem based on mean velocity and tube diameter is 3263, which is equivalent to the Reynolds number Rec based on centerline velocity and tube radius. It is emphasized that the velocity profile described previously is not imposed as a Dirichlet boundary condition but as one part of the Neumann boundary condition according to Equation (4), besides the pressure part, which will be described in the following paragraph. Furthermore, it is emphasized that the form of the velocity profile (parabolic in the present case) represents the only assumption required a priori for being able to prescribe the proposed Neumann boundary conditions. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
NEUMANN INFLOW BOUNDARY CONDITIONS IN BIOMECHANICS
567
For the pressure part, the analytical value for the (kinematic) pressure at the inflow boundary p in is prescribed. The momentum balance in a stationary state formulated in x3 -direction results in @p D @x3
@ 2 u3 @ 2 u3 C @2 x 1 @2 x 2
.
(10)
Thus, integrated along the tube length, a value p in D 4=r 2 u3,c l3 D 1.1438 m2 =s2 is obtained. The viscous part is omitted because of its minor significance in the present case. The tube is spatially discretized by 80 elements in x3 -direction. Furthermore, 24 elements are arranged along the diameter of the cross section, refined towards the hull (Figure 1). This discretization results in an overall number of 34,560 elements. A constant time-step length t D 0.004 s is used. In a first test case, relatively fast sinusoidal start-up and slow-down phases of 0.2 s, respectively, with a ‘holding’ phase of 1.1 s in between, where the maximum inflow value is kept constant, are prescribed both for inspiration and expiration. As a second test case, two sine cycles are computed, with the length of a sine cycle covering one inspiration and one expiration phase of equivalent length, respectively, assumed to be T D 3 s. The prescribed flow rate curves for both cases are depicted in Figure 2. The first test case is used for investigating whether the expected velocity and pressure profiles in a stationary state can be obtained via the present definition of boundary conditions. It is sufficient to consider the pressure profiles along the tube; they are extracted during the ‘holding’ phases both for inspiration and expiration. The pressure profiles for the inspiration and expiration phases, respectively, are depicted in Figure 3, along with the respective analytical profiles. It may be observed that, except for slight deviations at the respective ends of the tube, the analytical profiles are matched by the numerically obtained profiles. For the second test case, Figure 4 shows velocity results at
Figure 1. Spatial discretization of tube. 0.8 Q [l/s] 0.6 0.4 0.2 0.0
0
1
2
3
4
5
t [s] 6
-0.2 -0.4 -0.6 -0.8
Figure 2. Prescribed flow rate curves for tube; dashed line: curve for test case 1, solid line: curve for test case 2. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
568
V. GRAVEMEIER ET AL.
0.00
0
20
40
60
80
100 p [10-2m2/s2] 120
0.05
0.10
x3 [m] 0.15
Figure 3. Pressure profiles for tube obtained during ‘holding’ phases in test case 1; solid (orange) line, inspiration phase; dashed (orange) line, expiration phase; and dotted (black) lines, respective analytical profiles.
Figure 4. Velocity vectors on colored velocity magnitude distribution (red color indicates high velocity and blue color low velocity); from left to right, t = 0.75, 1.30, 1.50, 1.70, 2.25, 2.80, 3.00, and 3.20 s.
various instances in time. The depiction series starts with the peak of the inspiration phase, passes the flow reversal when turning from inspiration to expiration for the first time and the peak of the expiration phase, and concludes with an impression of the velocity field at 3.2 s within the second inspiration phase. 4.2. Tracheobronchial airways The second example considers flow and pressure in a realistic tracheobronchial geometry undergoing mechanical ventilation. With this example, we demonstrate the applicability of the boundary condition when flow and pressure curves are given, which are different in form. The geometry used in the simulations is segmented from end-expiratory CT images (slice thickness and pixel size 0.7344 mm) of a healthy patient using the commercially available segmentation software MIMICS (Materialise NV, Leuven, Belgium). The computational geometry is shown in Figure 5. The geometry begins at the level of the trachea and progresses up to the fourth generation. The domain Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
569
NEUMANN INFLOW BOUNDARY CONDITIONS IN BIOMECHANICS
Figure 5. Tracheobronchial geometry used in the simulations.
is discretized in space utilizing mixed (predominantly hexahedral) element types, giving a total of 785,033 elements and 623,478 nodes (2.5M degrees of freedom). In the near wall regions, boundary layer meshes are fitted to capture the relevant flow physics in these locations. A constant time-step length t D 0.001 s is used. As boundary conditions for the simulation, mechanical ventilation data obtained from an intensive care unit patient suffering from acute lung injury are applied. Figure 6 shows the pressure and flow at the level of the trachea for this patient. Note that in the present simulation, the pressure component of the traction is scaled down from 103 to 101 because of the extremely high traction the former pressure level would create on the boundary surfaces as a result of the low air density. For such a simulation, a very long start-up time would be required to arrive at this high pressure level. Remedies to this issue will be the subject of future research. It is emphasized that the absolute pressure drop between the inflow and outflow boundaries remains unchanged during scaling. For the flow and pressure at the five peripheral boundaries to be determined, a zero-dimensional (0D) airway model is generated using the same CT images as the ones used for the three-dimensional (3D) tracheobronchial geometry (cf. Figure 7). The 0D model was recently shown to successfully mimic the underlying physiology of a lung in [42], thus generating the correct pressure and flow rate values within the lung. The flow and pressure determination utilizes the following steps. The mechanical ventilation data are first applied to the 0D model. Flow and pressure data are then extracted from the 0D model at the same locations where the 3D model is trimmed. Finally, extracted flow and pressure values are applied to the 3D model as inflow and outflow boundaries according to Equations (4) and (5), respectively. For simplicity, constant velocity profiles are chosen at all inflow boundaries for both inspiration and expiration phases.
0.6
3000
0.4
2500
0.2
2000
p (Pa)
Q(l/s)
0 −0.2
1500 1000
−0.4 500
−0.6
0
−0.8 −1
0
0.5
1
1.5
2
2.5
time(s)
3
3.5
4
−500
0
0.5
1
1.5
2
2.5
3
3.5
4
time (s)
Figure 6. Ventilator-measured flow (left) and pressure (right) curves at the level of the trachea used in the simulations (note that a positive flow represents flow into the domain). Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
570
V. GRAVEMEIER ET AL.
Figure 8 shows pressure and volume-rendered velocity distributions in the tracheobronchial model. The results demonstrate clearly the difference in behavior of the pressure and velocity over the course of the cycle. In particular, we see that between t D 0.8 and t D 1.4 s, the velocity distribution is essentially unchanged, but the pressure continues to rise; if the pressure and velocity curves had not been simultaneously applied, the pressure would have remained constant. At t D 1.84 s
Figure 7. Patient-specific zero-dimensional (0D) airway model.
Figure 8. Pressure (top) and velocity magnitude (bottom) in the tracheobronchial region. From left to right: t D 0.8, 1.4 (maximum pressure), 1.6 and 1.84 s (maximum expiration). Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
571
NEUMANN INFLOW BOUNDARY CONDITIONS IN BIOMECHANICS 30
3 sim input
2
sim input
25
Pressure (Pa)
velocity (m/s)
1 0 −1 −2
15 10 5
−3
0
−4 −5
20
0
0.5
1
1.5
2
2.5
3
3.5
−5
4
0
0.5
1
1.5
time(s)
2.5
3
3.5
4
45
6 sim input
5 4
35
3
30
2 1 0 −1
sim input
40
Pressure (Pa)
velocity (m/s)
2
time(s)
25 20 15 10 5
−2 −3
0
−4
−5
0
0.5
1
1.5
2
2.5
3
3.5
4
0
0.5
time(s)
1
1.5
2
2.5
3
3.5
4
time(s)
Figure 9. Comparison of the simulated (sim) and ventilator derived (input) velocity and pressure. Mean velocity (left, top) and pressure (right, top) at the trachea inlet. Mean velocity (left, bottom) and pressure (right, bottom) at the right bottom outlet. Negative velocity represents flow out of the domain.
(maximum expiration), we observe reversal of the high and low pressure regions and the formation of a dominant flow, which comes from the left superior lobe. The nonsymmetrical nature of the flow distribution is due to the realistic boundary conditions obtained from the 0D lung model. Figure 9 details a comparison of the ventilator-derived pressure and velocity used to determine the prescribed Neumann boundary condition (cf. Equation (4) for inflow and Equation (5) for outflow) with the pressure and velocity obtained from the simulation results, referred to as ‘input’ and ‘sim’, respectively. Overall, we see excellent agreement, however clearly there are some discrepancies, particularly when the slew rates are high (cf. the difference in pressure and velocity at t1.5 s). This mismatch during these rapid acceleration phases is due to transient fluid inertia. In particular, the velocity cannot keep up with the rapidly changing pressure gradient; hence, distribution of the total momentum flux between the velocity and pressure terms is affected. 5. CONCLUSIONS A particular formulation of Neumann boundary conditions that prescribes the total momentum flux at inflow boundaries of biomechanical problems has been proposed. This approach enables the simultaneous application of velocity/flow rate and pressure curves, which may be required in problems of respiratory mechanics. One of the most prominent features of such respiratory configurations are complete flow reversals, that is, inflow boundaries usually become outflow boundaries (and vice versa) multiple times during the simulation. The suitability of the proposed boundary condition formulation has been demonstrated for two respiratory flow examples, a rather simple tube and tracheobronchial airways (up to the fourth generation, segmented from end-expiratory CT images). A residual-based variational multiscale (or stabilized) finite element method has been used for all numerical simulations. For the second complex example, the boundary conditions have been generated from mechanical ventilation data obtained from an intensive care unit patient suffering Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
572
V. GRAVEMEIER ET AL.
from acute lung injury. Among other things, it has been possible to replicate analytical pressure profiles in the first (tube) example as well as a correct distribution of the prescribed total momentum flux at the inflow boundary into velocity and pressure part in the second (tracheobronchial airways) example. Although the focus of the numerical examples in this study has been on respiratory flows with complete flow reversals, the proposed formulation is suited for cardiovascular flow problems with partial retrograde flow as well. Instabilities observed for such problems (e.g., [21]) are resolved by the present approach without requiring the additional consideration of a Lagrange multiplier technique. The appropriateness of the proposed method for situations where the flow partially reverses, was already demonstrated in [26–28] for various turbulent flow problems; vortices may leave the computational domain and provoke reverse flow in local sections of the respective outflow boundaries in such situations. In the future, these Neumann boundary conditions prescribing total momentum flux will be applied to cardiovascular and respiratory flows which are fully coupled to 0D models of the downstream vessels (see, for example, [15] and [19]). In those approaches, the flow is passed from the 3D domain to the 0D domain and the pressure is returned to the 3D domain during the subsequent time step. The approach presented in this paper allows the same method to be applied; however, the total traction rather than the pressure will be returned to the 3D domain.
ACKNOWLEDGEMENTS
Support from the Emmy Noether Program (first author) and the Protective Artificial Respiration (PAR) Program through projects WA 1521/6-2, WA 1521/9-2, and WA 1521/13-1 (second and third author) of the Deutsche Forschungsgemeinschaft (DFG) is gratefully acknowledged. The authors would like to thank Professor J. Guttmann, Department of Anesthesiology and Intensive Care Medicine, University Hospital Freiburg, Germany, for kindly providing the mechanical ventilation data.
REFERENCES 1. Kleinstreuer C, Zhang Z, Li Z, Roberts WL, Rojas C. A new methodology for targeting drug-aerosols in the human respiratory system. International Journal of Heat and Mass Transfer 2008; 51:5528–5589. 2. Comerford A, Bauer G, Wall WA. Nanoparticle transport in a realistic model of the tracheobronchial region. International Journal for Numerical Methods in Biomedical Engineering 2010; 26:904–914. 3. Wall WA, Wiechert L, Comerford A, Rausch S. Towards a comprehensive computational model for the respiratory system. International Journal for Numerical Methods in Biomedical Engineering 2010; 26:807–827. 4. Figueroa CA, Taylor CA, Yeh V, Chiou AJ, Zarins CK. Effect of curvature on displacement forces acting on aortic endografts: a three dimensional computational analysis. Journal of Endovascular Therapy 2009; 16:284–294. 5. Probst M, Lülfesmann M, Nicolai M, Bücker M, Behr M, Bischof CH. Sensitivity of optimal shapes of artificial grafts with respect to flow parameters. Computer Methods in Applied Mechanics and Engineering 2010; 199:997–1005. 6. Taylor CA, Draney MT, Ku JP, Parker D, Steele BN, Wang K, Zarins CK. Predictive medicine: computational techniques in therapeutic decision-making. Computer Aided Surgery 1999; 4:231–247. 7. Stroud JS, Berger SA, Saloner D. Numerical analysis of flow through a severely stenotic carotid artery bifurcation. Journal of Biomechanical Engineering 2002; 124:9–20. 8. Zhang Z, Kleinstreuer C. Transient airflow structures and particle transport in a sequentially branching lung airway model. Physics of Fluids 2002; 14:862–880. 9. Nowak N, Kakade PP, Annapragada AV. Computational fluid dynamics simulation of airflow and aerosol deposition in human lungs. Annals of Biomedical Engineering 2003; 31:374–390. 10. Green AS. Modelling peak-flow wall shear stress in major airways of the lung. Journal of Biomechanics 2004; 37:661–667. 11. Kabilan S, Lin C-L, Hoffman EA. Characteristics of airflow in a CT-based ovine lung: a numerical study. Journal of Applied Physiology 2006; 102:1469–1482. 12. Freitas RK, Schröder W. Numerical investigation of the three-dimensional flow in a human lung model. Journal of Biomechanics 2008; 41:2446–2457. 13. Xia G, Tawhai MH, Hoffman EA, Lin C-L. Airway wall stiffening increases peak wall shear stress: a fluid-structure interaction study in rigid and compliant airways. Annals of Biomedical Engineering 2010; 38:1836–1853. 14. Calay RK, Kurujareon J, Holdö AE. Numerical simulation of respiratory flow patterns within human lung. Respiration Physiology 2002; 130:201–221. 15. Vignon-Clementel IE, Figueroa CA, Jansen KE, Taylor CA. Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries. Computer Methods in Applied Mechanics and Engineering 2006; 195:3776–3796. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
NEUMANN INFLOW BOUNDARY CONDITIONS IN BIOMECHANICS
573
16. Blanco PJ, Feijòo RA, Urquiza SA. A unified variational approach for coupling 3D–1D models and its blood flow applications. Computer Methods in Applied Mechanics and Engineering 2007; 196:4391–4410. 17. Formaggia L, Moura A, Nobile F. On the stability of the coupling of 3D and 1D fluid-structure interaction models for blood flow simulations. ESAIM Mathematical Modelling and Numerical Analysis 2007; 41:743–769. 18. Ma B, Lutchen KR. An anatomically based hybrid computational model of the human lung and its application to low frequency oscillatory mechanics. Annals of Biomedical Engineering 2006; 34:1691–1704. 19. Comerford A, Förster C, Wall WA. Structured tree impedance outflow boundary conditions for 3D lung simulations. Journal of Biomechanical Engineering 2010; 132:081002. 20. Kim HJ, Vignon-Clementel IE, Figueroa CA, LaDisa JF, Jansen KE, Feinstein JA, Taylor CA. On coupling a lumped parameter heart model and a three-dimensional finite element aorta model. Annals of Biomedical Engineering 2009; 37:2153–2169. 21. Kim HJ, Figueroa CA, Hughes TJR, Jansen KE, Taylor CA. Augmented Lagrangian method for constraining the shape of velocity profiles at outlet boundaries for three-dimensional finite element simulations of blood flow. Computer Methods in Applied Mechanics and Engineering 2009; 198:3551–3566. 22. Formaggia L, Gerbeau J-F, Nobile F, Quarteroni A. Numerical treatment of defective boundary conditions for the Navier–Stokes equations. SIAM Journal on Numerical Analysis 2002; 40:376–401. 23. Vergara C. Nitsche’s method for defective boundary value problems in incompressible fluid-dynamics. Journal of Scientific Computing 2011; 46:100–123. 24. Hughes TJR, Franca LP, Hulbert M. A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective–diffusive equations. Computer Methods in Applied Mechanics and Engineering 1989; 73:173–189. 25. Hughes TJR, Wells GN. Conservation properties for the Galerkin and stabilised forms of the advection–diffusion and incompressible Navier–Stokes equations. Computer Methods in Applied Mechanics and Engineering 2005; 194:1141–1159. 26. Bazilevs Y, Michler C, Calo VM, Hughes TJR. Isogeometric variational multiscale modeling of wall-bounded turbulent flows with weakly-enforced boundary conditions on unstretched meshes. Computer Methods in Applied Mechanics and Engineering 2010; 199:780–790. 27. Gravemeier V, Wall WA. An algebraic variational multiscale-multigrid method for large-eddy simulation of turbulent variable-density flow at low Mach number. Journal of Computational Physics 2010; 229:6047–6070. 28. Gravemeier V, Kronbichler M, Gee MW, Wall WA. An algebraic variational multiscale-multigrid method for large eddy simulation: generalized-˛ time integration, Fourier analysis and application to turbulent flow past a square-section cylinder. Computational Mechanics 2011; 47:217–233. 29. Hughes TJR, Scovazzi G, Franca LP. Multiscale and stabilized methods. In Encyclopedia of Computational Mechanics, Stein E, de Borst R, Hughes TJR (eds). Wiley: Chichester, 2004. 30. Braack M, Burman E, John V, Lube G. Stabilized finite element methods for the generalized Oseen problem. Computer Methods in Applied Mechanics and Engineering 2007; 196:853–866. 31. Moser RD, Kim J, Mansour NN. Direct numerical simulation of turbulent channel flow up to Re = 590. Physics of Fluids 1999; 11:943–945. 32. Gravemeier V, Wall WA. Residual-based variational multiscale methods for laminar, transitional and turbulent variable-density flow at low Mach number. International Journal for Numerical Methods in Fluids 2011; 65:1260–1278. 33. Comerford A, Gravemeier V, Wall WA. An algebraic variational multiscale-multigrid method for large-eddy simulation of turbulent pulsatile flows in complex geometries with detailed insight into pulmonary airway flow, 2011. Preprint, submitted for publication in International Journal for Numerical Methods in Fluids. 34. Gravemeier V, Gee MW, Wall WA. An algebraic variational multiscale-multigrid method based on plain aggregation for convection–diffusion problems. Computer Methods in Applied Mechanics and Engineering 2009; 198:3821–3835. 35. Gravemeier V, Gee MW, Kronbichler M, Wall WA. An algebraic variational multiscale-multigrid method for large eddy simulation of turbulent flow. Computer Methods in Applied Mechanics and Engineering 2010; 199:853–864. 36. Donea J, Huerta A. Finite Element Methods for Flow Problems. Wiley: Chichester, 2003. 37. Jansen KE, Collis SS, Whiting CH, Shakib F. A better consistency for low-order stabilized finite element methods. Computer Methods in Applied Mechanics and Engineering 1999; 174:153–170. 38. Franca LP, Frey SL. Stabilized finite element methods: II. The incompressible Navier–Stokes equations. Computer Methods in Applied Mechanics and Engineering 1992; 99:209–233. 39. Franca LP, Valentin F. On an improved unusual stabilized finite element method for the advective–reactive–diffusive equation. Computer Methods in Applied Mechanics and Engineering 2000; 190:1785–1800. 40. Barrenechea G, Valentin F. An unusual stabilized finite element method for a generalized Stokes problem. Numerische Mathematik 2002; 92:652–677. 41. Jansen KE, Whiting CH, Hulbert GM. A generalized-˛ method for integrating the filtered Navier–Stokes equations with a stabilized finite element method. Computer Methods in Applied Mechanics and Engineering 2000; 190:305–319. 42. Ismail M, Comerford A, Wall WA. Coupled and reduced dimensional modelling of respiratory mechanics during spontaneous breathing, 2011. Preprint, submitted for publication in Annals of Biomedical Engineering.
Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. 2012; 28:560–573 DOI: 10.1002/cnm
