A finite element solver and energy stable coupling for 3D and 1D fluid models ^{†}^{†}thanks: This work has been supported in part by RFBR grants 120100283, 110100971, 110100767, 120133084
Abstract
The paper develops a solver based on a conforming finite element method for a 3D–1D coupled incompressible flow problem. New coupling conditions are introduced to ensure a suitable bound for the cumulative energy of the model. We study the stability and accuracy of the discretization method, and the performance of some stateoftheart linear algebraic solvers for such flow configurations. Motivated by the simulation of the flow over inferior vena cava (IVC) filter, we consider the coupling of a 1D fluid model and a 3D fluid model posed in a domain with anisotropic inclusions. The relevance of our approach to realistic cardiovascular simulations is demonstrated by computing a blood flow over a model IVC filter.

Keywords: geometrical multiscale modeling, 3D1D coupling, fluid flows, cardiovascular simulations, finite element method, iterative methods
1 Introduction
Coupling a 3D fluid flow model and a system of hyperbolic equations posed on a 1D graph is a well established approach for numerical simulations of blood flows in a system of vessels [32]. Such a geometric multiscale strategy is particularly efficient, when the attention to local flow details and the qualitative assessment of global flow statistics are both important. The relevance to cardiovascular simulations and challenging mathematical problems of coupling parabolic 3D and hyperbolic 1D equations put 3D–1D flow problems in the focus of intensive research. Thus, the coupling of a 3D fluid/structure interaction problem with a reduced 1D model merged to outflow boundary, which acts as an absorbing device, was studied in [15]. The coupling of a 3D fluid problem with multiple downstream 1D models in the context of a finite element method was considered in [38]. In [33], a system of a 3D fluid/structure interaction problem and a 1D finite element method model of the whole arterial tree was implemented to model the carotid artery blood flow; and in [6] a unified variational formulation for multidimensional models was introduced. A splitting method, extending the pressurecorrection scheme to 3D–1D coupled systems, was studied in [26].
In most of these studies, the 3D model was a generic fluidelasticity or rigid fluid model, while numerical validations were commonly done for cylindric type 3D domains (with rigid or elastic walls); several authors considered geometries with bifurcation [38, 33] or constrained geometries (modeling a stenosed artery) [6]. More complicated geometries occur in simulations of blood flows, if one is interested in modeling the effect of endovascular implants, such as inferior vena cava (IVC) filters. In numerical simulations, a part of a vessel with an intravenous filter leads to the computational 3D domain with strongly anisotropic inclusions. A downstream flow behind the implant may exhibit a complex structure with traveling vortices, swirls, and recirculation regions (the latter may occur if plaque is captured by the filter). Moreover, the IVC flow is strongly influenced by the contraction of the heart, and both forward (towards the heart) and reverse (from the heart) flows occur within one cardiac cycle. Downstream coupling conditions for such flows may be a delicate issue. Thus, the flow over an IVC filter is an interesting and challenging problem for a 3D1D flow numerical solver.
The coupling conditions of 3D and 1D fluid models and their properties were studied by several authors. Coupling conditions and algorithms based on subdomain iterations were introduced in [15], and the stability properties of each subproblem were analyzed separately. The first analysis of two models together was done in [16]. In that paper, it was noted that if the NavierStokes equations are taken in the rotation form and suitably coupled with a 1D downstream flow model, then one can show a bound for the joint energy of the system. It is, however, well known that using a finite element method for the rotation form of the NavierStokes equations needs special care [22], and setting appropriate outflow boundary conditions can be an issue. In the present paper, we introduce an energy consistent coupling with a 1D model for the convection form of the NavierStokes equations. The joint energy of a coupled 3D1D model is appropriately balanced and dissipates for viscous flows.
Handling highly anisotropic structures is a wellknown challenge in numerical flow simulations and analysis. There are only a few computational studies addressing the dynamics of blood flows in vessels with implanted filters. Recently, Vassilevski et al. [34] numerically approached the problem of intravenous filter optimization using a finitedifference method on octree cartesian meshes to resolve the geometry of implants. In that paper, it was also discussed how the effect of an implant can be accounted in a 1D model through a modification of a vessel wall state equation (see also [36, 37] for the development of this method for atherosclerotic blood vessels). In the present paper, we take another approach and locally resolve the full 3D model, while keeping the state equation unchanged. We report on a finite element method for modeling a 3D1D coupled fluid problem, when the 3D domain has anisotropic inclusions. Naturally, this leads to meshes containing possibly anisotropic tetrahedra. We study the performance of the finite element method both by considering the accuracy of solutions and by monitoring the convergence of one stateoftheart linear algebra solver for the systems of linear algebraic equations to be solved on every time step of the method. We are interested in the ability of the solver to predict such important statistics as the drag force experienced by an intravenous implant.
The remainder of the paper is organized as follows. In section 2, we review 3D and 1D fluid models and discuss coupling conditions. The stability properties of the coupled model are also addressed in section 2. Section 3 presents a timestepping numerical scheme and an algebraic solver. In section 4, we validate the 3D finite element solver and the coupled method by considering the benchmark problem of a flow past a 3D cylinder and a problem with an analytical solution. The application of the method to simulate a blood flow over a model IVC filter is given in section 5. Numerical experiments were performed using the Ani3D finite element package [40], which was used to generate tetrahedra subdivisions of 3D domains, to build stiffness matrices, and to implement the linear algebra solvers described in section 3.
2 The 1D3D coupled model
This section reviews 3D and 1D fluid models and describes the coupling of the models. In this study, the 3D model is assumed ‘rigid’.
2.1 The 3D model
Consider a flow of a viscous incompressible Newtonian fluid in a bounded domain . We shall distinguish between the inflow part of the boundary, , the noslip and nopenetration part (rigid walls), , and the outflow part of the boundary, . On the inflow part we assume a given velocity profile. The outflow boundary conditions are defined by setting the normal stress tensor equal to a given vector function . Thus, the 3D model is the classical NavierStokes equations in pressurevelocity variables:
(1) 
Here is the outward normal vector to . The system is also supplemented with initial condition () for in .
We remark that the notion of ‘inflow’ and ‘outflow’ boundary is used here and further in the text conventionally, since the inequalities or are not necessarily pointwise satisfied on or , respectively. In applications we consider, the mean flux, averaged in space and in time, is expected to be negative at and positive at . However, for certain , the flux may take positive values at (negative at ).
If the solution to (1) is sufficiently smooth and the inflow boundary conditions are homogeneous, the following energy balance holds:
Here and in the rest of the paper, denotes the norm. If one assumes
(2) 
then solutions to (1) satisfy the a priori energy inequality and this opens possibilities for showing partial wellposedness results. Even though, the assumption (2) was used for the purpose of analysis in the literature on 3D1D blood flow models (see [15, 16]), it is hard to verify (2) for practical flows. Moreover, the inequality (2) no longer holds if reverse flows occur, as, for example, happens in IVC [27, 39].
2.2 The 1D model
A onedimensional model can be derived from the NavierStokes equations posed in a long axisymmetric elastic pipe by integrating over cross section, making some simplifying assumptions and considering integral average quantities as unknowns, see, e.g., [1, 32]. Let be the cross section of the pipe normal to , is the area of and is the axial velocity. Introduce the averaged variables: the mean axial velocity and the mean pressure:
We consider the model given by the following system of equations for unknowns :
(3) 
For initial conditions, the mean velocity and the cross section area are prescribed, , . Here is the external pressure, is a function modelling the source or sink of the fluid, as may be required in hemodynamic simulations, if a blood loss happens in a vessel. Further, we assume and , so from now has the meaning of the difference between the fluid pressure and the external pressure. The term accounts for external forces, such as gravity or friction. Following [1, 28], we set
(4) 
Here is the viscosity coefficient, is the pipe diameter, is the reference area (in the hemodynamic applications is the cross section area of a vessel at rest) and
The last equation in (3) relates the pressure to the cross section area. The function is defined by the elasticity model of the pipe walls, is the elasticity model parameter. We use the one from [28]:
(5) 
Other algebraic defining relations linking the mean pressure and the cross section area are known from the literature, see, e.g., [32]. They are equally well suited for the purpose of this paper.
In [15], the authors derived the energy equality for the onedimensional fluid model in variables, where is the flux. Up to possibly different choices of and , the formulation in [15] is equivalent to (3). Written in terms of , and , the energy equality for (3) is (recall that we assumed ):
(6) 
with the energy functional
For given in (5), the second term in the definition of is always positive, making positive for all . The choice of in (4) ensures that the second term on the lefthand side of (6) is positive as well. Thus, for the homogenous boundary conditions the energy of the 1D model dissipates: .
System (3) is hyperbolic and can be integrated along characteristics. To see this, we write (3) in the divergence form:
with , , . Denote by and (i=1,2) the left eigenvectors and eigenvalues of the Jacobian . One finds
and
Thus, system (3) can be written in the characteristic form:
(7) 
Under physiological conditions in hemodynamics, it holds , implying that the eigenvalues have opposite signs. Therefore, is the incoming characteristic from point and is the incoming characteristic from point . Two boundary conditions, one at and one at , are enough to close the system.
2.3 The coupling of 1D and 3D models
We now consider two domains and , as shown in Figure 1. In we pose the 1D fluid model and in we pose the full threedimensional NavierStokes equations. The onedimensional domain is coupled to the downstream boundary of .
There are several options to define the coupling conditions of 1D and 3D models, as discussed, for example, in [15]. One can ask for the continuity of the mean pressure, the mean axial velocity, the flux, the normal cross section area, the averaged normal stress, or the entering characteristic. In general, the continuity of all these quantities cannot be satisfied simultaneously. A choice has to be made.
One common choice, see, e.g., [33, 6], is to impose the continuity of the normal stress and the flux:
(8)  
(9) 
This choice is, however, known to be energy inconsistent in the following sense. Assume the homogeneous boundary conditions on and the downstream end of . Then the cumulative energy balance of the coupled model is
(10) 
with . is a positive coefficient defined from (4). Easy to see that for coupling conditions (8)–(9) the right hand side of (10) reduces to
In general, it is not clear if this quantity is nonpositive and thus if the cumulative energy is properly dissipated, as it holds for the full 3D NavierStokes equations with the Dirichlet homogeneous boundary conditions. To circumvent this inconsistency, it was suggested in [16] to replace condition (8) by the continuity of socalled total stress:
(11) 
Condition (11) together with the continuity of the flux, i.e. condition (9), makes the right hand side of (10) equal to zero. Hence, the cumulative energy dissipates. However, setting the total stress equal to a constant is not a consistent outflow condition for the simplest Poiseuille flow. Moreover, (11) is the natural boundary condition for the rotation form of the NavierStokes equations. Using it with the common convection or conservation forms leads to nonlinear coupling conditions and often requires iterative treatment [7]. Although the rotation form is an interesting alternative to the standard convection form, it is still not a standard option in the existing software and its use requires certain care [22, 23].
The condition for the normal stress, as in (8), is the natural boundary condition for the commonly used convection and conservation forms of the NavierStokes equations. Such a condition has been shown to be surprisingly useful as outflow boundary condition [19]. Thus, instead of keeping (9) and changing (8) to the total stress condition, we retain (8) and instead of (9) assume the continuity of the linear combination of the fluid flux and the energy flux:
(12) 
One easily checks that the combination of (8) and (12) ensures the righthand side of (10) equal to zero and thus the correct energy balance and energy inequality are valid as formulated in the following theorem.
Theorem 2.1
Remark 2.1
Since has the meaning of the difference between fluid and external pressures, it can be negative. In this case, more than one value of may satisfy (12). To ensure that the coupled model is not defective, one has to prescribe a particular rule for choosing the root of the cubic equation (12). In our numerical experiments, we take which is the closest to .
Remark 2.2
The boundary condition (12) is the combination of fluid and energy fluxes and so it does not guarantee to conserve the ‘mass’ of the entire coupled system. Although, we do not observe any perceptible generation or loss of mass in our numerical experiments, it does not necessarily mean that for all problems this effect should be negligible. Actually, one may consider any other linear combination of fluid and energy fluxes coupling on 3D1D boundary to compromise between energy stability and mass conservation.
In practice, one may also be interested in coupling the 1D fluid model to the upstream boundary of the 3D domain. Hence, we now consider three domains , , and , as shown in Figure 2. In and the simplified 1D model is posed and in the full threedimensional NavierStokes equations are solved. The domain is coupled to the inflow (upstream) boundary of and is coupled to the outflow (downstream) boundary of . The downstream coupling is described above. In the literature, it is common not to distinguish between upstream and downstream coupling boundary conditions. For example, in [6, 33] conditions (8),(9) are assumed both on upstream and downstream boundaries. Following this paradigm, one may consider (8),(9) or energy consistent conditions (8), (12) as the coupling conditions on between 1D model posed in and 3D model posed in . Note that in entirely 3D fluid flows simulations, inflow and outflow boundary conditions usually differ. If a numerical approach to 1D3D problem is based on subdomains splitting (see the next section for an example), then it is appropriate to distinguish between upstream and downstream coupling conditions. Thus, we impose the upstream coupling conditions in such a way that the 3D problem is supplied with the Dirichlet inflow boundary conditions. This is a standard choice for incompressible viscous fluid flows solvers and is especially convenient if third parties or legacy codes are separately used to compute 3D and 1D solutions, and they communicate only through coupling conditions.
For the upstream boundary, , we introduce a reference velocity profile , , such that . Then the boundary condition on is Dirichlet, given by
(14) 
Setting ensures the continuity of the flux (9) on . If is found to satisfy the equation
then the coupling condition (12) is valid on . Two more scalar boundary conditions are required for the 1D model in . We assume that or are given in and in an absorbing condition is prescribed: in computations we set in ; another reasonable absorbing condition would be setting the incoming characteristic equals zero. On the downstream end of , we also assume an absorbing boundary condition.
We summarize the properties for the 3D1D coupling introduced in this section:
3 Discretization and algebraic solver
In this section, we introduce a splitting numerical timeintegration algorithm based on subdomain splitting. Further, we consider a fully discrete problem and review one stateoftheart algebraic solver.
Denote by , , , and approximations to the corresponding unknown variables at time . Given these approximations, we compute , , , and for () in three steps:

Step 1. Integrate (7) for , with given on the upstream end of the interval and the absorbing downstream condition at .

Step 2. Set according to (14), with , and compute and as the linear extrapolations of and from times and . Solve the linearized NavierStokes problem in for :

Step 3. Find from
Now, using for boundary condition in and the absorbing boundary condition in , we integrate (7) for to find , in .
For the numerical integration of the 1D model equations, we use a first order monotone finite difference scheme applied to the characteristic form (7), see [35]. To handle the 3D model, one has to solve on every time step the linearized NavierStokes equations, also known as the Oseen problem:
(15) 
where , , the body forces term and the advection velocity field depend on previous time velocity approximations, , , and . Here and in the remainder of this section, we dropped out the timestep dependence () index for unknown velocity and pressure.
To discretize the Oseen problem (15), we consider a conforming finite element method. Denote the finite element velocity and pressure spaces by and , respectively. Let be the subspace of of all FE velocity functions vanishing at .
The finite element problem reads: Find , , and satisfying
(16) 
with
Let for . We assume the ellipticity, the continuity, and the stability conditions:
(17) 
(18) 
(19) 
with positive meshindependent constants , , , and . Condition (18) is wellknown as the LBB or infsup stability condition [9].
Let and be bases of and , respectively. Define the following matrices:
The linear algebraic system corresponding to (16) (the discrete Oseen system) takes the form
(20) 
The right hand side accounts for body forces and inhomogeneous velocity boundary conditions. To solve (20), we consider a Krylov subspace iterative method, with the block triangular preconditioner [12, 21]:
(21) 
The matrix is a preconditioner for the matrix , such that may be considered as an inexact solver for linear systems involving . The matrix is a preconditioner for the pressure Schur complement of (20), . In the algorithm, one needs the actions of and on subvectors, rather than the matrices , explicitly. Once good preconditioners for and are available, a preconditioned Krylov subspace method, such as GMRES or BiCGstab, is the efficient solver. In the literature, one can find geometric or algebraic multigrid (see, e.g., [12] and references therein) or domain decomposition [17, 31] algorithms which provide effective preconditioners for a range of and various meshes. We use one Vcycle of the algebraic multigrid method [30] to define .
Defining an appropriate pressure Schur complement preconditioner is more challenging. In this paper, we follow the approach of Kay et al. [20]. First, we define the pressure mass and velocity mass matrices:
The original pressure convectiondiffusion (PCD) preconditioner, proposed in [20], is defined through its inverse:
(22) 
Here denotes an approximate solve with the pressure mass matrix. Matrices and are approximations to convectiondiffusion and Laplacian operators in , respectively. Both and (explicitly or implicitly) assume some pressure boundary conditions to be prescribed.
If defines continuous pressure approximations, one can use the conforming discretization of the pressure Poisson problem with Neumann boundary conditions:
Likewise, Neumann boundary conditions are conventionally used to define the pressure convectiondiffusion problem on . However, the optimal boundary conditions setup both for and depends on the type of the boundary and flow regime, see [14, 25].
We use a modified PCD preconditioner defined below. This modification partially obviates the issue of setting pressure boundary conditions and is consistent with the Cahouet–Chabard preconditioner [10], if the inertia terms are neglected. The Cahouet–Chabard preconditioner is the standard choice for the timedependent Stokes problem and enjoys the solid mathematical analysis in this case [24]. To define the preconditioner, we introduce the discrete advection matrix for continuous pressure approximations as
Then the modified pressure convectiondiffusion preconditioner (mPCD) is (compare to (22) ):
where is a diagonal approximation to the velocity mass matrix.
Regarding the numerical analysis of the algebraic solver used here, we note the following. The eigenvalues bounds of the preconditioned Schur complement:
(23) 
were proved for and the LBB stable finite elements in [13] and for a more general case in [25]. The constants are independent of the meshsize , but may depend on the ellipticity, continuity and stability constants in (17)–(19), and thus may depend on the problem parameters. In particular, the pressure stability constant , and so from (23), depends on the geometry of the domain [11] (tending to zero for long or narrow domains) and for certain FE pairs depends on the anisotropy ratio of a triangulation [2]. Both of this dependencies require certain care in using the approach for computing flows in 3D elongated domains with thin and anisotropic inclusions (prototypical for simulating a flow over IVC filter).
Characterizing the rate of convergence of nonsymmetric preconditioned iterations is a difficult task. In particular, eigenvalue information alone may not be sufficient to give meaningful estimates of the convergence rate of a method like preconditioned GMRES [18]. Nevertheless, experience shows that for many linear systems arising in practice, a wellclustered spectrum (away from zero) usually results in rapid convergence of the preconditioned iteration. This said, we should mention that a rigorous proof of the GMRES convergence applied to (20), with blocktriangular preconditioner (21), is not available in the literature (except the special case, when is symmetric [21, 5]). Thus, the numerical assessment of the approach is of practical interest.
4 Accuracy and efficiency assessment
In this section, we validate the accuracy and stability of the solver for the 3D1D coupled fluid model by (i) comparing the computed discrete solutions against an analytical solution for a problem with simple geometry; (ii) computing the drag coefficient and the pressure drop value for the flow around the 3D circular cylinder. The TaylorHood P2P1 elements were used for the velocitypressure approximation. The resulting linear algebraic systems (20) were solved by the preconditioned BiCGstab method. The initial guess in the BiCGstab method was zero on the first time step and equal to the for the subsequent time steps. The stopping criteria was the decrease of the Euclidean norm of the residual.
4.1 Test with analytical solution
First we consider an example with analytical solution. The 3D domain is . The circular cross sections are the inflow and outflow boundaries. Domains and are two intervals of length 5. The analytical solution is given by
(24) 
with , , , . This solution satisfies the continuity of flux condition on the coupling boundaries. The righthand sides , and were set accordingly. In this test, the 3D domain was triangulated using the global refinement of an initial mesh, resulting in the sequence of meshes (further denoted by mesh 1, mesh 2, mesh 3), with the number of tetrahedra , respectively. Since we use the first order scheme for the 1D problem, the mesh size in and was divided by 4 on each level of refinement: . The corresponding time step was halved for every spacial refinement, so we use for mesh 1, mesh 2, and mesh 3, respectively.
mesh 1  0.37  0.30  0.30  8.10E2 

mesh 2  8.60E2 (2.11)  6.13E2 (2.29)  5.97E2 (2.33)  4.04E2 (1.00) 
mesh 3  2.66E2 (1.69)  1.68E2 (1.87)  1.62E2 (1.88)  2.02E2 (1.00) 
in  in  

mesh 1  0.25  8.84E4  0.19  7.02E4 
mesh 2  5.78E2 (2.11)  2.03E4 (2.12)  3.67E2 (2.37)  1.00E4 (2.81) 
mesh 3  1.43E2 (2.02)  5.01E5 (2.02)  5.64E3 (2.70)  2.95E5 (1.76) 
Based on the energy balance (13), the natural norms for measuring error in are and for velocity and for pressure. These norms and, additionally, for velocity error are shown in Table 1. The error norms in the 1D domains coupled to the inflow and the outflow boundaries of are shown in Table 2. We observe the expected second order of convergence for all variables except for pressure in 3D. We remark that the integral in time error norms were computed approximately using the quadrature rule:
Other integral norms were computed in the same way.
mesh 1  10.7  10.7  13.4  14.5 

mesh 2  5.59  6.69  8.09  8.99 
mesh 3  4.42  4.42  6.50  7.88 
1  0.5  0.1  0.05  0.01  0.005  0.001  

0.1  11.78  10.78  14.00  *  *  *  * 
0.05  7.16  6.84  7.21  12.00  *  *  * 
0.01  4.28  4.45  5.65  6.27  6.30  6.63  7.14 
The average number of iterations of the BiCGstab method for solving (20), with blocktriangular preconditioner (21), are shown in Tables 3 and 4. Table 3 shows that the convergence of the preconditioned BiCGstab method depends only slightly on the viscosity parameter and improves when the grid is refined. The latter observation is consistent with independent eigenvalue bounds in (23) and with numerical results reported in [12] for steady problems. Such robust behavior with respect to is observed only for sufficiently small values of the time step . The results in Table 4 show that for small s the preconditioned BiCGstab method fails to converge unless .
4.2 Flow around circular cylinder
Motivated by the simulation of blood flows over an intravenous filter, we experiment with flows in a 3D domain having an inclusion and coupled with 1D model at the outflow boundary. Interesting statistics for such applications are the drag force acting on inclusions and the pressure drop. To validate the ability of the 3D solver to predict these statistics, we consider two benchmark problems of channel flows past a 3D cylinder with circular cross sections [29, 8]. The 3D flow domain with a cylinder is shown in Figure 3.
The noslip and nopenetration boundary condition is prescribed on the channel walls and the cylinder surface. The inlet velocity is given by
here is the width of the channel. The kinematic viscosity of the fluid in this test is and its density is . The Reynolds number, , is defined based on the cylinder width and . We consider the following two benchmark problems from [29]:

Problem P1: Steady flow with ();

Problem P2: Unsteady flow with varying Reynolds number for , .
The benchmarks setups do not specify outflow boundary conditions. Hence, on the outflow boundary we apply the 3D1D coupling using the new conditions (8), (12) so that numerical performance of the coupling can be verified.
The statistics of interest are the following:

The difference between the pressure values in points and .

The drag coefficient given by an integral over the surface of the cylinder :
(25) Here is the normal vector to the cylinder surface pointing to and is a tangent vector.
For problem P2, the reference velocity in (25) is .
mesh  % err  % err  

coarse  6.149  0.58%  0.1679  1.81%  11.5 
fine  6.196  0.17%  0.1678  1.87%  10.5 
Schäfer & Turek  [6.05, 6.25]  [0.165, 0.175]  
Braack & Richter  6.185  0.1710 
mesh  % err  

coarse  3.273  0.76%  0.115  11.7 
fine  3.311  0.39%  0.107  10.6 
Schäfer & Turek  [3.2,3.3]  [0.11, 0.09]  
Bayraktar et al.  3.298  – 
For these benchmark problems, the paper [29] collects several DNS results based on various finite element, finite volume discretizations of the NavierStokes equations and the Lattice Boltzmann method. In [29], the authors provided reference intervals, where the statistics are expected to converge. Using a higher order finite element method and locally refined adaptive meshes, more accurate reference values of and were found in [8] for the steady state solution (problem P1) and in [3] for unsteady problem P2. For the computations we use two meshes: a ‘coarse’ and a ‘fine’ ones, both adaptively refined towards cylinder. The coarser mesh is build of 35803 tetrahedra, which results in 53061 velocity d.o.f. and 8767 pressure d.o.f. for the TaylorHood P2P1 element. The finer mesh consists of 51634 tetrahedra, which results in 73635 velocity d.o.f. and 12321 pressure d.o.f. Both coarse and fine mesh consist of regular tetrahedra. The refinement ratio is about 20 and 60 for the coarse and the fine meshes, respectively. We remark that the fine mesh has four times as many tetrahedra touching the cylinder as the coarse mesh. The time steps are and for the coarse and the fine meshes, respectively.
We first show in Tables 5 and 6 results for problems P1 and P2 obtained with the coarse and the fine meshes. For all settings, the computed values are within “reference intervals” from [29] (except for problem P2, but in this case the upper reference bound appears to be tough). The computed drag coefficients were well within 1% of reference values and pressure drop within 2%. This is a good result for the number of the degrees of freedom involved. Indeed, the results shown in [8, 3, 29] for meshes with about the same number of degrees of freedom show comparable or worse accuracy. In Figure 4, we show the computed evolution of the drag coefficient for problem P2 and compare it to the reference results. The computed drag coefficients match the reference curve very well. We conclude that the conforming finite element method with the coupling outflow conditions is a reliable and stable approach for the simulation of such flow problems.
5 Simulations of a flow over a model IVC filter
The development of endovascular devices is the challenging problem of cardiovascular medicine. One example is the design of vascular filters implanted in inferior vena cava (IVC) to prevent a blockage of the main artery of the lung or one of its branches by a substance that has traveled from elsewhere in the body through the bloodstream. The filter is typically made of thin rigid metal wires as illustrated in Figure 5 (left). Numerical simulation is an important tool that helps in finding an optimal filter design. Thin and anisotropic construction of a IVC filter requires adaptive grid refinement and makes computations of flows in such domains not an easy task. In this section, we demonstrate the ability of the numerical method to treat such problems in a stable way. One statistic of interest here is the drag force experienced by a filter. We recall that in this paper we do not account for the elastic properties of the vessel walls, which are otherwise important in practice.
We consider a segment ( long) of IVC with elliptic cross section . The filter is placed on the distance from inflow, it is long and the diameter of its 12 wire legs is . Blood is assumed to be incompressible fluid with dynamic viscosity equal to and density equal to .
A blood flow in IVC is strongly influenced by the contraction of the heart. The IVC have pulsatile waveforms with two peaks and reverse flow [39] occurring on every cardiac cycle. We consider the Doppler blood flow waveforms of IVC reported in [39] and approximate them by a smooth periodic function plotted in Figure 5 (right). Note that the presence of significant reverse flows in IVC differs this problem from computing arteria flows, where such phenomenon does not typically occur.
On the inflow and outflow, the 3D vessel is coupled to 1D models as described in section 2.3. Each 1D model consists of equations (3)–(5) posed on intervals of length. Periodic velocity with waveform as shown in Figure 5 is prescribed on the upstream part of the 1D model coupled to . The maximum 1D model velocity of yields the maximum inlet velocity in of about . This agrees with the measurements in [39]. The coupling conditions are the same regardless of the mean flow direction.
The mesh was adapted towards the filter, so the ratio of largest and smallest element diameters was about , the maximum elements anisotropy ratio was about . The resulting mesh is illustrated in Figure 6. The time step in 3D model was set equal to . The BiCGstab iterative method, with preconditioner (21) was used to solve discrete Oseen subproblems. The stopping criterion was the reduction of the residual by the factor of . The average number of linear iterations on every time step was about 35. We found that choosing time step larger for this problem, leads to the significant increase of the linear iteration counts and makes ‘long time’ computations nonfeasible.
We visualize the computed solutions in Figure 7 by showing the values of the component of the velocity in several cutplanes orthogonal to axis. Behind the filter the velocity component eventually has negative values, indicating the occurrence of circulation zones and ‘returning’ flows. Note that the solution behind the filter is no longer axialsymmetric: a perturbation to solution induced by nonsymmetric tetrahedral grid is sufficient for the von Karman type flow instability to develop behind the filter.
Figure 8 (left) shows the time evolution of the drag force experienced by the filter. After the instantaneous start, the flow needs few cycles to obtain the periodic regime. In general, the drag force follows the pattern of the inflow waveform. In particular, the filter experiences forces both in downstream and upstream directions at different periods of the cardiac cycle. The right plot in Figure 8 shows the mean axial velocity in the middle point of the 1D model before and after the 3D domain with cava filter. It is remarkable that after few cycles, when the flow is periodic, the waveforms in the 1D domains coupled to upstream and downstream boundaries are very close. This suggests that the coupling conditions are efficient in conserving averaged flow quantities such as mean flux.
6 Conclusions
We reviewed the 3D and 1D models of fluid flows and some existing coupling conditions for these models. New coupling conditions were introduced and shown to ensure a suitable bound for the cumulative energy of the model. The conditions were found to perform stable in several numerical tests with analytical and benchmark solutions. For the example of the flow around IVC filter, the coupled numerical model was found to capture the periodic flow regime and correct 1D waveforms before and after 3D domain. The model was able to handle ‘opposite direction’ flow, i.e. the flow where the ‘upstream’ boundary (boundary with Dirichlet boundary conditions) becomes the outflow boundary for a period of time. The preconditioned BiCGstab method with one stateoftheart preconditioner applied to the linearized finite element NavierStokes problem performs well. However, often the time step should be taken small enough to make the linear solver converge sufficiently fast. Overall, the coupled 3D1D model together with the conforming finite element method and preconditioned iterative strategy was demonstrated as a reliable tool for the simulation of such biological flows as the flow over an inferior vena cava filter.
Acknowledgements
The authors are grateful to A. Danilov (INM RAS, Moscow) for his help in building tetrahedra meshes in ANI3D, S. Simakov (MIPT, Moscow) for providing us with 1D fluid solver code, and Yu.Vassilevski (INM RAS, Moscow) for the implementation of the algebraic solver for Problem (20). We would like to thank A. Quaini and S. Canic (UH, Houston, TX) for fruitful discussions and pointing to papers [27, 39].
References
 [1] M.V. Abakumov, K.V. Gavrilyuk, N.B. Esikova, A.V. Lukshin, S.I. Mukhin, N.V. Sosnin, V.F. Tishkin, A.P. Favorskij, Mathematical model of the hemodynamics of the cardiovascular system, Differ. Equations. 33 (7) (1997) 895–900.
 [2] Th. Apel, H.M. Randrianarivony, Stability of discretizations of the Stokes problem on anisotropic meshes, Mathematics and Computers in Simulation. 61 (2003) 437–447.
 [3] E. Bayraktar, O. Mierka, S. Turek, Benchmark computations of 3D laminar flow around a cylinder with CFX, OpenFOAM and FeatFlow, International Journal of Computational Science and Engineering. 7 (2012), 253–266.
 [4] M. Benzi, M.A. Olshanskii, An augmented Lagrangianbased approach to the Oseen problem, SIAM J. Sci. Comput. 28 (2006) 2095–2113.
 [5] M. Benzi, M.A. Olshanskii, Fieldofvalues convergence analysis of Augmented Lagrangian preconditioners for the linearized NavierStokes problem, SIAM J. Numer. Anal. 49 (2011) 770–788.
 [6] P.J. Blanco, R.A. Feijoro, S.A. Urquiza, A unified variational approach for coupling 3D–1D models and its blood flow applications, Comput. Methods Appl. Mech. Engrg. 196 (2007) 4391–4410.
 [7] P.J.Blanco, S. Deparis, A.C.I. Malossi, On the continuity of mean total normal stress in geometrical multiscale cardiovascular problems, EPFLARTICLE182892, 2012
 [8] M. Braack, T. Richter, Solutions of 3D NavierStokes benchmark problems with adaptive finite elements, Computers & Fluids. 35 (2006) 372–392.
 [9] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer, 1991.
 [10] J. Cahouet, J. P. Chabard, Some fast D finite element solvers for the generalized Stokes problem, Internat. J. Numer. Methods Fluids. 8 (1988) 869–895.
 [11] E.V. Chizhonkov, M.A. Olshanskii, On the domain geometry dependence of the LBB condition, Math. Modelling Numer. Anal.: . 34 (2000) 935–951.
 [12] H.C. Elman, D.J. Silvester, A.J. Wathen, Finite Elements and Fast Iterative Solvers: With Applications in Incompressible Fluid Dynamics, Numer. Math. Sci. Comput., Oxford University Press, New York, 2005.
 [13] H.C. Elman, D. Loghin, A. J. Wathen, Preconditioning techniques for Newton’s method for the incompressible Navier–Stokes equations, BIT. 43 (2003) 961–974.
 [14] H.C. Elman, R.S. Tuminaro, Boundary conditions in approximate commutator preconditioners for the NavierStokes equations, Electronic Transactions on Numerical Analysis. 35 (2009) 257–280.
 [15] L. Formaggia, J.F. Gerbeau, F. Nobile, A. Quarteroni, On the coupling of 3D and 1D NavierStokes equations for flow problems in compliant vessels, Computer Methods in Applied Mechanics and Engineering. 191 (2001) 561–582.
 [16] L. Formaggia, A. Moura, F. Nobile, On the stability of the coupling of 3D and 1D fluidstructure interaction models for blood flow simulations, ESAIM: Mathematical Modelling and Numerical Analysis. 41 (4) (2007) 743–769.
 [17] M. Garbey, Yu. A. Kuznetsov, Yu. V. Vassilevski, Parallel Schwarz method for a convectiondiffusion problem, SIAM J. Sci. Comput. 22 (2000) 891–916.
 [18] A. Greenbaum, V. Pták, Z. Strakos̆, Any nonincreasing convergence curve is possible for GMRES, SIAM J. Matrix Anal. Appl. 17 (1996) 465–469.
 [19] J.G. Heywood, R. Rannacher, S. Turek, Artificial boundaries and flux and pressure conditions for the incompressible NavierStokes equations, International Journal for Numerical Methods in Fluids. 22 (1996) 325–352.
 [20] D. Kay, D. Loghin, A.J. Wathen, A preconditioner for the steadystate Navier–Stokes equations, SIAM J. Sci. Comput. 24 (2002) 237–256.
 [21] A. Klawonn, G. Starke, Block triangular preconditioners for nonsymmetric saddle point problem, Numer. Math. 81 (1999) 577–594.
 [22] W. Layton, C.C. Manica, M. Neda, M.A. Olshanskii, L.G. Rebholz, On the accuracy of the rotation form in simulations of the NavierStokes equations, Journal of Computational Physics. 228 (2009) 3433–3447.
 [23] M.A. Olshanskii, A low order Galerkin finite element method for the NavierStokes equations of steady incompressible flow: A stabilization issue and iterative methods, Comp. Meth. Appl. Mech. Eng. 191 (2002) 5515–5536.
 [24] M.A. Olshanskii, J. Peters, A. Reusken, Uniform preconditioners for a parameter dependent saddle point problem with application to generalized Stokes interface equations, Numerische Mathematik. 105 (2006) 159–191.
 [25] M.A. Olshanskii, Yu.V. Vassilevski, Pressure Schur complement preconditioners for the discrete Oseen problem, SIAM J.Sci.Comp. 29 (2007) 2686–2704.
 [26] G. Papadakis, Coupling 3D and 1D fluidï¿½structureinteraction models for wave propagation in flexible vessels using a finite volume pressurecorrection scheme, Commun. Numer. Meth. Engng. 25 (2009) 533–551.
 [27] W. T. Pua, T. Ishiwatac, A. L. Juraszeka, Q. Mac, S. Izumo, GATA4 is a dosagesensitive regulator of cardiac morphogenesis, Developmental Biology. 275 (2004) 235–244.
 [28] S.S. Simakov, A.S. Kholodov, Computational study of oxygen concentration in human blood under low frequency disturbances, Mathematical models and computer simulations. 1 (2) 2009 283–295.
 [29] M. Schäfer, S. Turek, The benchmark problem “Flow around a cylinder”. In Hirschel EH (ed.), Flow Simulation with HighPerformance Computers II, vol. 52. Notes on Numerical Fluid Mechanics, Vieweg. 1996; 547–566
 [30] J.W. Ruge, K. Stüben, Algebraic multigrid, in Multigrid Methods, edited by S. F. McCormick, p. 73 (SIAM, Philadelphia, PA, 1987).
 [31] A. Quarteroni, A. Valli, Domain Decomposition Methods for Partial Differential Equations, Oxford University Press, Oxford, UK, 1999.
 [32] A. Quarteroni, L. Formaggiaa, A. Veneziani, Cardiovascular Mathematics: Modeling and Simulation of the Circulatory System, SpringerVerlag Itali, Milano 2009
 [33] S.A. Urquiza, P.J. Blanco, M.J. Vernere, R.A. Feijoro, Multidimensional modelling for the carotid artery blood flow, Comput. Methods Appl. Mech. Engrg. 195 (2006) 4002–4017.
 [34] Y.V. Vassilevski, S.S. Simakov, S.A. Kapranov, A multimodel approach to intravenous filter optimization, International Journal for Numerical Methods in Biomedical Engineering. 26 (2010) 915–925.
 [35] Yu. Vassilevskii, S. Simakov, V. Salamatova, Yu. Ivanov, T. Dobroserdova, Numerical issues of modelling blood flow in networks of vessels with pathologies, Russian Journal of Numerical Analysis and Mathematical Modelling. 26 (6) (2011) 605–622.
 [36] Y. Vassilevski, S. Simakov, V. Salamatova, Y. Ivanov, T. Dobroserdova, Vessel wall models for simulation of atherosclerotic vascular networks, Mathematical Modelling of Natural Phenomena. 6 (7) (2011) 82–99.
 [37] Y. Vassilevski, S. Simakov, V. Salamatova, Y. Ivanov, T. Dobroserdova, Blood flow simulation in atherosclerotic vascular network using fiberspring representation of diseased wall, Mathematical Modelling of Natural Phenomena. 6 (5) (2011) 333–349.
 [38] I.E. VignonClementel, C.A. Figueroa, K.E. Jansen, C.A. Taylor, Outflow boundary conditions for threedimensional finite element modeling of blood flow and pressure in arteries, Comput. Methods Appl. Mech. Engrg. 195 (2006) 3776–3796.
 [39] D. Zhang, T. Kanzaki, Doppler waveforms: the relation between ductus venosus and inferior vena cava, Ultrasound in Med. & Biol. 31 (9) (2005) 1173–1176.
 [40] ANI3D: Advanced Numerical Instruments. http://sourceforge.net/projects/ani3d