Enhancement of flow measurements using fluid-dynamic constraints

Journal of Computational Physics 344 (2017) 558–574

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Enhancement of flow measurements using fluid-dynamic constraints H. Egger a,∗ , T. Seitz a , C. Tropea b a b

Institute for Numerical Analysis and Scientific Computing, Department of Mathematics, TU Darmstadt, Germany Institute for Fluid Mechanics and Aerodynamics, Center of Smart Interfaces, TU Darmstadt, Germany

a r t i c l e

i n f o

Article history: Received 22 March 2016 Received in revised form 19 January 2017 Accepted 30 April 2017 Available online 15 May 2017 Keywords: Velocity measurements Denoising Optimal control with pdes Fluid dynamics Navier–Stokes equations Inverse problems Regularization

a b s t r a c t Novel experimental modalities acquire spatially resolved velocity measurements for steady state and transient flows which are of interest for engineering and biological applications. One of the drawbacks of such high resolution velocity data is their susceptibility to measurement errors. In this paper, we propose a novel filtering strategy that allows enhancement of the noisy measurements to obtain reconstruction of smooth divergence free velocity and corresponding pressure fields which together approximately comply to a prescribed flow model. The main step in our approach consists of the appropriate use of the velocity measurements in the design of a linearized flow model which can be shown to be well-posed and consistent with the true velocity and pressure fields up to measurement and modeling errors. The reconstruction procedure is then formulated as an optimal control problem for this linearized flow model. The resulting filter has analyzable smoothing and approximation properties. We briefly discuss the discretization of the approach by finite element methods and comment on the efficient solution by iterative methods. The capability of the proposed filter to significantly reduce data noise is demonstrated by numerical tests including the application to experimental data. In addition, we compare with other methods like smoothing and solenoidal filtering. © 2017 Elsevier Inc. All rights reserved.

1. Introduction Since numerous years the visualization of flow fields has had a significant impact on the systematic understanding and development of fluid-dynamic models as well as on the calibration and verification of computational methods. While traditional experimental techniques were able to provide only partial information about the flow field, novel measurement modalities such as particle tracking, tomographic particle imaging, or magnetic resonance velocimetry deliver spatially resolved three-dimensional velocity measurements [1–4]. In principle, these new methods therefore allow to image complex flow patterns in a wide range of engineering applications and even in biological in-vivo studies. Distributed flow measurements provide valuable information about simple and complex flows, but they are typically contaminated by measurement errors which limit their usability in practice to some extent. In order to make the flow measurements more suitable for further analysis, e.g., for model discrimination or for the assessment of derived quantities like pressure drop or wall shear stress, some sort of data post-processing is required [5]. A widely used technique in this context is data smoothing which can be accomplished, for instance, by Tikhonov regularization [6,7] given by

*

Corresponding author. E-mail address: [email protected] (H. Egger).

http://dx.doi.org/10.1016/j.jcp.2017.04.080 0021-9991/© 2017 Elsevier Inc. All rights reserved.

H. Egger et al. / Journal of Computational Physics 344 (2017) 558–574

min u − uδ 2 + α ∇ u2 .

559

(1)

u

Here and below uδ denotes the flow measurements and the minimizer will be the enhanced velocity field. Note that the penalization of the gradient term leads to a smoothed velocity reconstruction. The choice of the regularization parameter α allows a certain trade-off between smoothness and fit to the data. The underlying quadratic minimization problem can be solved efficiently by Fourier transform or multigrid iterative solvers which makes this filter very efficient in practice. Note that the above procedure and also various other image processing methods [8] successfully reduce high frequency components in the noisy measurements but do not utilize any information about the underlying physics. In many applications, the flow under consideration is incompressible and one might want to incorporate such prior knowledge into the reconstruction process. Requiring the improved velocity field to be divergence free and using a smoothing procedure similar to above, we obtain a constrained minimization problem

minu − uδ 2 + α ∇ u2 u

s.t.

∇ · u = 0.

(2a) (2b)

This quadratic programming problem can again be solved efficiently by iterative methods. Various computational strategies leading to related divergence free reconstructions have been investigated recently in the literature under the name divergence-free or solenoidal filtering; see e.g. [9–14]. Let us note that, although some noise reduction has been observed even for the case α = 0, the divergence constraint alone does formally not guarantee smoothness of the reconstruction. This can be seen from the Helmholtz decomposition of vector fields [15] and will be illustrated by numerical tests below. A natural extension of the solenoidal filtering approach, which takes into account only the incompressibility, is to incorporate into the reconstruction process also a model for the momentum balance. Since distributed measurement techniques typically acquire time averaged data, it seems reasonable to assume steady flow conditions and to consider, as a first step, the stationary Navier–Stokes equations as the governing physical model. The reconstruction could then be defined via

min u − uδ 2 + α f2 f,u,p

− ν u + u · ∇ u + ∇ p = f,

s.t.

∇ · u = 0.

(3a) (3b)

Here and below ν > 0 denotes the viscosity parameter which is assumed to be constant here. Additional boundary conditions have to be prescribed for a complete specification of the flow model. The residual f in the momentum equation measures the deviation from an idealized flow model due to unmodeled effects like time dependence or non-Newtonian behavior. Since a prescribed flow model is satisfied by the reconstructed fields, one may call such an approach a fluiddynamically consistent filter. Due to the presence of the differential terms in the flow model (3b), the reconstruction will be smooth automatically and no additional penalization of the velocity gradients is required. Moreover, some information about the pressure p is obtained. The system (3a)–(3b) has the form of an optimal control problem governed by the Navier–Stokes equations. Such problems have been investigated extensively in the literature; see e.g. [16–19] for steady and [20–24] for unsteady flow. Note that the nonlinearity in the momentum equation poses severe challenges, both, for the analysis and for the numerical solution. It is well-known, for instance, that the Navier–Stokes system admits a unique solution only for sufficiently small data [25,26]. Moreover, due to the nonlinear constraints, the optimization problem (3a)–(3b) is non-convex and may in general have many local minima. Both aspects make the computational solution demanding or even infeasible. In this paper, we therefore propose a strategy that allows us to take advantage of the benefits and at the same time to overcome the drawbacks of the previous approach. The basic step here is to use the distributed velocity measurements uδ in order to replace the nonlinear term in the momentum equation by some linearization; one may think of uδ · ∇ u as an approximation for u · ∇ u, although such a simple choice would not yield a well-posed problem in general due to lack of smoothness in the data. A proper linearization of the convective term u · ∇ u will however allow us to replace the nonlinear problem (3a)–(3b) by a quadratic optimal control problem with linearized constraints which has a unique minimizer that can be computed efficiently. The way in which we use the measurements in the governing equations is closely related to the equation error method, which is well-established in the context of parameter estimation [27,28]. In summary we thus obtain a well-posed and analyzable reconstruction method that produces a smooth divergence free velocity field and a corresponding pressure distribution which agree well with the measurements and at the same time approximately satisfy the prescribed fluid-dynamic model. A proper choice of the regularization parameter α will allow us to find a good balance between data fit and model errors. The remainder of the manuscript is organized as follows: In Section 2, we introduce the linearized fluid flow model underlying our reconstruction approach and we formulate appropriate boundary conditions. We then establish the wellposedness of the linearized flow problem and derive error estimates for the linearization step. In Section 3, we introduce the optimal control problem which serves as mathematical formulation of our reconstruction procedure. We prove existence and uniqueness of minimizers, provide some error estimates, and highlight a direct connection to the solenoidal filtering. Our approach is formulated in infinite dimensions and some discretization strategy is required in order to obtain implementable algorithms. In Section 4, we therefore outline the discretization by finite element methods and we briefly discuss

560

H. Egger et al. / Journal of Computational Physics 344 (2017) 558–574

efficient strategies for the numerical solution of the discretized optimal control problem. The validity of our theoretical results is illustrated in Section 5 by some computational tests for a simple model problem in which we also compare with the smoothing and the solenoidal filtering approaches outlined above. In order to demonstrate the applicability to problems of practical interest, we report in Section 6 about flow reconstructions for experimental data obtained with magnetic resonance velocimetry. The presentation concludes with a short summary and a discussion of open problems and possible directions for further research. 2. The linearized flow model Let us first introduce the linearized flow model that is later on used as a constraint in the reconstruction process and establish its well-posedness. Utilizing the fact that the true flow field satisfies a model of similar structure, we then derive some perturbation error estimates for the linearization step. For illustration of the results and later reference, we then discuss in some detail the Poiseuille flow between two parallel plates. 2.1. Geometric setting We start with fixing the geometric setting we have in mind. Let  ⊂ Rd , d = 2, 3 be some bounded Lipschitz domain. We assume that the boundary ∂ is piecewise smooth and can be split into three distinct parts ∂ wall , ∂in , and ∂out such that ∂in ∩ ∂out = ∅. One may think of a channel where ∂in is the inflow, ∂out the outflow, and ∂ wall the wall of the channel. 2.2. The linearized flow model The key step in our approach is to replace the nonlinear convective term u · ∇ u in (3b) by an appropriate linearization. For this, we make use of the simple identity

w · ∇ u = 12 w · ∇ u + 12 div(u ⊗ w) − 12 (∇ · w)u, (4)   where (w · ∇ u)i = j w j ∂ j ui and div(u ⊗ w)i = j ∂ j (ui w j ) following the usual notation. We use bold symbols throughout to denote vector valued functions and spaces of such functions. For an incompressible fluid, we have ∇ · u = 0 and can thus write the convective term u · ∇ u equivalently as 1 div(u ⊗ u). 2

1 u · ∇u + 2

Such a form of the convective term is sometimes employed in the design and analysis of numerical methods for incompressible flow due to its improved stability properties. As we will see below, the convective term written in this form gives an anti-symmetric contribution in the weak formulation of the flow problem and thus has no negative effect on the a-priori estimates for the solution, even if the convective velocity is not divergence-free. Using the velocity measurements uδ in part of the quadratic terms yields an approximation u · ∇ u ≈ 12 uδ · ∇ u + 12 div(u ⊗ uδ ) which can be used to replace the nonlinear convective term in the momentum equation (3b). This leads to the linearized flow model

−ν u + 12 uδ · ∇ u + 12 div(u ⊗ uδ ) + ∇ p = f

in ,

(5a)

∇ ·u=0

in .

(5b)

This is a generalized Oseen problem with convective velocity uδ that will neither be smooth nor divergence free in general. To complete the description of the flow model, we impose the following boundary conditions

(−ν ∇ u +

1 u ⊗ uδ 2

u=g

on ∂in ,

(5c)

u=0

on ∂ wall ,

(5d)

on ∂out .

(5e)

+ pI) · n = h

The full velocity field is thus prescribed at the inflow boundary and we use a Neumann-type boundary condition at the outflow. No-slip conditions are used at the walls of the channel but other types of the boundary conditions could be incorporated as well. The functions f, g, and h arising as data in the flow model will later enter the reconstruction process as additional parameters which are to be determined. 2.3. Well-posedness of the linearized flow model Since the data uδ are obtained by measurements, one can in general not require their spatial smoothness. It is therefore not clear a-priori, if the model (5a)–(5e) is meaningful from a mathematical point of view. As a first step, we thus want to establish the well-posedness of the linearized flow model.

H. Egger et al. / Journal of Computational Physics 344 (2017) 558–574

561

Theorem 2.1. Let uδ ∈ L3 (). Then for any data f ∈ L2 (), g ∈ H10 (∂in ), and h ∈ L2 (∂out ), the problem (5a)–(5e) has a unique weak solution u ∈ H1 () and p ∈ L2 (). Moreover

  uH1 () + pL2 () ≤ C fL2 () + gH1 (∂in ) + hL2 (∂out ) with C depending only on uδ L3 () , on the parameter ν , and on the geometry. In practice the measurements of the flow velocity will always be bounded in magnitude and the technical assumption uδ ∈ L 3 () therefore does not pose any real restriction. Proof. The result follows with standard arguments for the analysis of the stationary flow equations. Since the momentum equation is a bit non-standard, we sketch the main steps of the proof. The weak solution of problem (5a)–(5e) is characterized by the following mixed variational problem: Find u ∈ H1 () and p ∈ L2 () with u = g on ∂in and u = 0 on ∂ wall such that

a(u, v) + c (uδ ; u, v) + b(v, p) = (f, v) + (h, v)∂out b(u, q) = 0

(6) (7)

for all v ∈ H1 () and q ∈ L 2 () with v = 0 on ∂ wall ∪ ∂in . Here a(u, v) = ν (∇ u, ∇ v) is the bilinear form for the viscous term, c (uδ ; u, v) = 12 (uδ · ∇ u, v) − 12 (u, uδ · ∇ v) represents the convective term, and b(u, q) = −(∇ · u, q) the weak form for the divergence operator. Without loss of generality, we may assume that g = 0 in the sequel. Due to the special form of the convective terms, the bilinear form c (uδ ; ·, ·) is anti-symmetric, which particularly implies c (uδ ; u, u) = 0. The assumption uδ ∈ L3 () further implies that c (uδ ; u, v) is bounded for u, v ∈ H1 (). Standard arguments used for the analysis of the Oseen problem then yield the assertions; see [26, Ch. II] for details. 2 As can be seen from the arguments used in the proof, the special form of the convective term and of the outflow boundary condition were essential here to obtain the well-posedness and the energy estimate under minimal regularity assumptions on the measured flow field. 2.4. Estimates for the linearization error The above linearization procedure introduces some perturbations which we would like to quantify next. To be able to do so, we set up a flow model of similar structure which describes the true flow. Let u† and p† denote the true velocity and pressure fields which are assumed to be sufficiently smooth. Then

−ν u† + 12 u† · ∇ u† + 12 div(u† ⊗ u† ) + ∇ p† = f† , †

∇ · u = 0,

(8a) (8b)

for some appropriate function f† which can just be defined by the left hand side of the first equation. In a similar way, we can define functions g† and h† such that

u† = g† †

u =0 †

(−ν ∇ u +

1 † u 2

†

†

⊗ u + p I) · n = h

†

on ∂in ,

(8c)

on ∂ wall ,

(8d)

on ∂out .

(8e)

This system has the same form as (5a)–(5e) but with data f, g, h and convective velocity field uδ replaced by the functions f† , g† , h† , and u† , respectively. This allows us to estimate the difference between the solutions of (8a)–(8e) and the linearized flow model (5a)–(5e) in the following systematic way. Theorem 2.2. Let uδ ∈ L3 () and let (u, p) and (u† , p† ) be defined as above. Then

u − u† H1 () + p − p† L2 ()   ≤ C f − f† L2 () + g − g† H1 (∂in ) + h − h† L2 (∂out ) + u† − uδ L3 () with C depending only on the bounds for the data, the parameter ν , and the geometry. Proof. Let ( u,  p) be the solution of problem (5a)–(5e) with uδ replaced by u† . Then we may write

u − u† = w + z

and

p − p† = π + ψ,

562

H. Egger et al. / Journal of Computational Physics 344 (2017) 558–574

with functions w = u −  u, π = p − p, z =  u − u† , and ψ =  p − p† that can be estimated separately. From the definition of w and π , we observe that w = 0 on ∂in ∪ ∂ wall and

a(w, v) + c (uδ ; w, v) + b(v, π ) = c (u† − uδ ; u, v) b(w, q) = 0 for all v ∈ H1 () and q ∈ L 2 () with v = 0 on ∂ wall ∪ ∂in . Choosing v = w and q = −π as test functions and applying the Poincaré–Friedrichs inequality yields

c w2H1 () ≤ a(w, w) = a(w, w) + c (uδ ; w, w) + b(w, π )

= c (u† − uδ ; u, w) ≤ C u† − uδ L3 ()  uH1 () wH1 () . Since  uH1 () can be bounded uniformly by Theorem 2.1, we further obtain wH1 () ≤ C u† − uδ L3 () and π  L 2 () can be bounded by wH1 () as well with the usual arguments. Next observe that z = g − g† on ∂in and z = 0 on ∂ wall , and that

a(z, v) + c (u† ; z, v) + b(v, ψ) = (f − f† , v) + (h − h† , v)∂out b(z, q) = 0 for all test functions q ∈ L 2 () and v ∈ H1 () with v = 0 on ∂ wall ∪ ∂in . By Theorem 2.1 we thus obtain zH1 () +   ψ L 2 () ≤ C

f − f† L2 () +g − g† H1 (∂in ) +h − h† L2 (∂out ) . The assertion of Theorem 2.2 then follows by a combination of these estimates. 2 Remark 2.3. Let us note that, if the flow model is a reasonable approximation for the physical conditions, we may assume that f† , g† , and h† are known up to small perturbations. For illustration, let us discuss a particular example which will also serve as our test problem later on. 2.5. Poisseuille flow In simple geometries, the stationary Navier–Stokes equations can be solved analytically. Consider for instant the laminar flow between two parallel plates with distance d along a path of length L. Then the resulting flow field is given by

p† (x, y ) = p0 +

p L − p0 L



x

and

u† (x, y ) =

p L − p0 2ν L

 (dy − y 2 ), 0 ,

(9)

where p0 , p L denote the pressures at position x = 0 and x = L, respectively. Similar formulas are available for channels with various geometries [29]. Note that the solution (u† , p† ) given by the Poisseuille law (9) satisfies the system (8a)–(8e) with f† = 0 and functions g† and h† that can be computed from (9). The estimates of the previous theorem then read

u − u† H1 () + p − p† L2 ()   ≤ C f − f† L2 () + g − g† H1 (∂in ) + h − h† L2 (∂out ) + u† − uδ L3 () . The total error in the solution thus results from misspecifications of the model data f − f† , g − g† , h − h† on the one hand, and from perturbations uδ − u† in the flow measurements on the other. This observation will be the guideline for the formulation of our reconstruction method in the next section. 3. The reconstruction method For the systematic enhancement of the distributed velocity measurements uδ obtained with one of the imaging modalities mentioned in the introduction, we now consider the optimal control problem



min u − uδ 2L2 () + α f − f∗ 2L2 () + g − g∗ 2H1 (∂ ) + h − h∗ 2L2 (∂ ) out in

f,g,h,u,p

s.t. (5a)–(5e),



(10a) (10b)

which can be interpreted as a linearization of problem (3a)–(3b). The choice of function spaces over which is minimized is guided by the theorems of the previous section and the norms in the regularization terms. The functions f∗ , g∗ , and h∗ serve as approximations for the unknown problem data f† , g† , and h† in the theoretically exact fluid-dynamic model (8a)–(8e) and enter as additional model parameters. The reconstructed field thus minimizes a weighted sum of deviations from the velocity data and the prescribed flow model and model data. The choice of the regularization parameter α allows us to balance the two error contributions.

H. Egger et al. / Journal of Computational Physics 344 (2017) 558–574

563

3.1. Existence of a unique minimizer Due to the well-posedness of the linearized flow problem (5a)–(5e), we can express u = u(f, g, h) and p = p(f, g, h) in terms of the functions f, g, and h. This allows us to eliminate the fields u and p from (10a)–(10b) and to obtain the following equivalent quadratic minimization problem

min J α (f, g, h)

(11)

f,g,h

with reduced cost functional J α defined by





J α (f, g, h) = u(f, g, h) − uδ 2L2 () + α f − f∗ 2L2 () + g − g∗ 2H1 (∂ ) + h − h∗ 2L2 (∂ ) . out in The existence of a unique minimizer for (11) follows with standard arguments in convex analysis. By the equivalence with the problem (10a)–(10b), we also obtain the well-posedness of the original formulation. Theorem 3.1. Let uδ ∈ L3 () and let f∗ ∈ L2 (), g∗ ∈ H10 (∂in ) and h∗ ∈ L2 (∂out ). Then for any regularization parameter α > 0, the reduced minimization problem (11) has a unique solution (fα , gα , hα ) with components fα ∈ L2 (), gα ∈ H10 (∂in ), and hα ∈ L2 (∂out ). Together with uα = u(fα , gα , hα ) and pα = p(fα , gα , hα ) this yields the unique solution of problem (10a)–(10b). Proof. The mapping (f, g, h) → (u(f, g, h), p(f, g, h)) is affine linear and continuous. As a consequence, the functional J α is quadratic, bounded from below, strictly convex, lower semi-continuous, and coercive. This implies existence of a unique minimizer. 2 As we will see below, both formulations (10a)–(10b) as well as (11) are both well suited as a starting point for the design of efficient numerical solution procedures. 3.2. Estimates for the reconstruction error As a theoretical justification for the proposed method let us next present some quantitative estimates for the reconstruction error which illustrate what kind of numerical results can be expected and which allow us to draw some conclusions about the proper choice of the regularization parameter. Theorem 3.2. Let (uα , pα ) denote the velocity and pressure components of the unique solution of problem (10a)–(10b) and assume that u† − uδ L3 () ≤ δ . Then the following estimates hold true:





(i) uα − u† 2L2 () ≤ C δ 2 + α f† − f∗ 2L2 () + g† − g∗ 2H1 (∂ ) + h† − h∗ 2L2 (∂ ) . out in





(ii) uα − u† 2H1 () ≤ C δ 2 + δ 2 /α + f† − f∗ 2L2 () + g† − g∗ 2L2 (∂ ) + h† − h∗ 2L2 (∂ ) . out in The second bound also holds for the error pα − p† 2L2 () in the pressure. u,  p) denote the solution of (6)–(7) with f, g, h replaced by f† , g† , and h† , respectively. Then  u − u† H1 () + Proof. Let (

 p − p†  L 2 () ≤ C u† − uδ L3 () , which follows like the estimate for the function w in the proof of Theorem 2.2. By definition of uα as a minimizer, we further have

uα − uδ 2L2 ()

  ≤ uα − uδ 2L2 () + α fα − f∗ 2L2 () + gα − g∗ 2H1 (∂ ) + hα − h∗ 2L2 (∂ ) out in   ≤  u − uδ 2L2 () + α f† − f∗ 2L2 () + g† − g∗ 2H1 (∂ ) + h† − h∗ 2L2 (∂ ) . in

out

The first assertion now follows by combining the two estimates and using the assumption on the data error together with the continuous embedding of L 3 () into L 2 () on bounded domains. For the second estimate, we use the triangle inequality to obtain

uα − u† H1 () + pα − p† L2 () ≤ ( u − u† H1 () +  p − p† L2 () ) + (uα −  uH1 () + pα − pL2 () ). The first term can be estimated by  u − u† H1 () + pα − p L 2 () ≤ C δ as above. Proceeding as in the proof of Theorem 2.1, the remaining term can be bounded by

564

H. Egger et al. / Journal of Computational Physics 344 (2017) 558–574

uα −  uH1 (∂) + pα − pL2 ()   † ≤ C fα − f L2 () + gα − g† H1 (∂in ) + hα − h† L2 (∂out ) . Let us consider the first term on the right hand side in more detail: Using the triangle inequality, we can see that fα − f† L2 () ≤ fα − f∗ L2 () + f† − f∗ L2 () ; note that the second term already appears in the final result. By the definition of the minimizers, the first term can be further estimated by

fα − f∗ 2L2 () ≤ α −1 uα − uδ 2L2 () + fα − f∗ 2L2 () + gα − g∗ 2H1 (∂ ≤ α −1  u − uδ 2L2 () + f† − f∗ 2L2 () + g† − g∗ 2H1 (∂

in

in )

+ hα − h∗ 2L2 (∂

out )

+ h† − h∗ 2L2 (∂ )

out )

.

Together with the estimates for  u − u† and the bound on the data error, this yields

 fα − f∗ 2L2 () ≤ C δ 2 /α + f† − f∗ 2L2 () + g† − g∗ 2H1 (∂

in )

+ h† − h∗ 2L2 (∂

out )



.

The same arguments lead to estimates for gα − g∗ 2H1 (∂ ) and hα − h∗ L2 (∂out ) which completes the proof of the second in assertion. 2 Remark 3.3. The estimates in Theorem 3.2 show that a proper choice of the regularization parameter α allows to obtain a balance between data fit and model error. In particular, a good fit to the measurements can always be obtained by choosing α sufficiently small. If the proposed flow model is a good description of the physical conditions, i.e., if the model errors f† − f∗ , g† − g∗ , and h† − h∗  are sufficiently small, one can actually choose the regularization parameter in the order of one and still obtain a good fit for the velocity field in the stronger norm and also a good reconstruction of the pressure. The approach is therefore robust with respect to the choice of the regularization parameter; see also the numerical results below. 3.3. Poisseuille flow Let us again return to the model problem discussed in Section 2.5. In this case, we may choose the model data as f∗ = f† , g∗ = g† , and h∗ = h† with f† = 0 and known functions g† , h† that can be computed from the Poisseuille law (9). The estimates of the previous theorem then simplify to

uα − u† L2 () ≤ C δ

and

√ uα − u† H1 () + pα − p† L 2 () ≤ C δ(1 + 1/ α ).

The reconstruction errors will therefore be in the order of the measurement errors if we choose the regularization parameter α in the order of one! In particular, any fixed choice of the regularization parameter will lead to O (δ) convergence for the velocity errors in the L 2 - and the H 1 -norm. A more sophisticated choice of the regularization parameter is however required if the underlying flow model does not describe the physical situation sufficiently well, i.e., if the model data f∗ , g∗ , and h∗ are not chosen appropriately. 3.4. Further properties As a final step of our theoretical considerations, let us comment on two further important properties of the reduced cost functional J α and the velocity component uα of the minimizer. Theorem 3.4. Let uα = u(fα , gα , hα ) be defined as above. Then (i) minf,g,h J α (f, g, h) ≤ minf,g,h J β (f, g, h) whenever α ≤ β . L2

(ii) uα → u S F with α → 0, where u S F is the solution of (2a)–(2b) with α = 0. Proof. By definition of uα as component of the minimizer, we have

min J α (f, g, h) = J α (fα , gα , hα ) ≤ J α (fβ , gβ , hβ ) ≤ J β (fβ , gβ , hβ ), f,g,h

where in the last step we used the condition α ≤ β and the positivity of the regularization term. To show the second assertion, let us note that for any ε > 0 one can find a smooth divergence free approximation  uS F SF SF SF S F for the velocity field u that is obtained by (2a)–(2b) with α = 0 and such that u − u L2 () ≤ ε . By plugging  u into (5a)–(5e), we obtain corresponding residuals  f,  g, and  h. Using the definition of uα , we further get

H. Egger et al. / Journal of Computational Physics 344 (2017) 558–574

565

uα − uδ 2L2 ()

  ≤ uα − uδ 2L2 () + α fα − f∗ 2L2 () + gα − g∗ 2H1 (∂ ) + hα − g∗ 2L2 (∂ ) out in   ≤  u S F − uδ 2L2 () + α  f − f∗ 2L2 () +  g − g∗ 2H1 (∂ ) +  h − h∗ 2L2 (∂ ) . out

in

This shows that

lim sup uα − uδ L2 () ≤  u S F − uδ L2 () ≤ u S F − uδ L2 () + ε α →0

for any

ε > 0. We therefore conclude that lim supα →0 uα − uδ L2 () ≤ u S F − uδ L2 () . Since uα L 2 () ≤ C for all α > 0,

¯ ∈ L2 (), and by the lower semi-continuity of the norm, we we can select a subsequence {uα } converging weakly to some u can deduce

u¯ − uδ L2 () ≤ lim inf uα − uδ L2 () ≤ lim sup uα − uδ L2 () ≤ u S F − uδ L2 ()

α →0

Moreover, since ∇ · uα = 0 for all ¯ = uS F . 2 conclude that u

α →0

α > 0, we also have ∇ · u¯ = 0. But since u S F is the unique minimizer of (2a)–(2b), we may

Remark 3.5. For α → 0, not only the minimal values of the cost functional J α but also the data residuals uα − uδ L2 () can be shown to decrease monotonically. This allows us to choose the regularization parameter α via a discrepancy principle. The fact that the reconstructed velocity fields uα converges to the solution u S F of the solenoidal filtering problem with α → 0 indicates that our reconstruction approach yields an enhancement of and is at least as stable as the solenoidal filtering (2a)–(2b) with α = 0. 4. Numerical realization In order to obtain implementable algorithms that may be used for actual computations, we still have to discretize the reconstruction method proposed in this paper. Since the numerical approximation of optimal control problems is wellunderstood, we only sketch the main ideas here. 4.1. Discretization of the fluid-dynamic model For discretization of the state system (5a)–(5e), we use a standard Galerkin method with an inf-sup stable pair of finite element spaces [30,15]. This leads to an algebraic system of the form

Au + B p = Mf + ε1 REg + Nh,

(12a)

Bu = 0,

(12b)

where A = ν K + C(uδ ) + ε1 R, the matrix K represents the vector Laplacian, M the mass matrix, C(uδ ) is the anti-symmetric

convective term, and B is the discrete divergence operator. The Dirichlet boundary conditions are incorporated here by a penalty approach with ε being a small parameter. The matrix R represents the corresponding integrals over the boundary ∂in ∪ ∂ wall , and E realizes an extension of boundary values into the domain by zero. Similarly, the matrix N accounts for integrals over ∂out . The vector uδ denotes the finite element representation of the velocity measurements. When using a finite element discretization, all matrices will be sparse. For details on the implementation, we refer to standard textbooks [15,26]. 4.2. Discretization of the optimal control problem Using the above notation, the discretized optimal control problem can be written as



min u − uδ 2M + α f − f∗ 2M + g − g∗ 2G + h − h∗ 2H

u,p,f,g,h

s.t. (12a)–(12b).



(13a) (13b)

Here r2M = r Mr denotes the norm induced by the positive definite matrix M. The Gramian matrices G and H induce the corresponding norms on the boundary. The first order optimality conditions for (13a)–(13b) are obtained by differentiating the associated Lagrangian. Since the problem under investigation is quadratic and strictly convex, the first order conditions are necessary and sufficient. Due to the use of finite elements, the discrete optimality system is sparse symmetric and indefinite and can be solved efficiently by appropriate preconditioned iterative methods [31–33].

566

H. Egger et al. / Journal of Computational Physics 344 (2017) 558–574

The discrete state system (12a)–(12b) can again be used to express u = u(f, g, h) and p = p(f, g, h) as functions of the unknown model data. This yields the corresponding reduced problem



min u(f, g, h) − uδ 2M + α f − f∗ 2M + g − g∗ 2G + h − h∗ 2H



f,g,h

(14)

which is quadratic and strictly convex and can be solved by a preconditioned conjugate gradient method. The computation of the gradients can be realized efficiently via adjoint problems which have a similar structure as (12a)–(12b); see e.g. [34] for more details on the implementation of a related problem. 4.3. Notes on other filtering approaches For later reference, let us also briefly discuss the implementation of the smoothing and the solenoidal filtering approaches outlined in the introduction. Using the same notation as above, the discrete version of the smoothing filter (1) can be expressed as

min u − uδ 2M + α u2K

(15)

u

and the minimizer is characterized by the regularized normal equations

(M + α K)u = Muδ .

(16)

This system is symmetric and positive definite and can be solved efficiently by the conjugate gradient method. Parameter robust multigrid preconditioning [35] may be applied to obtain an algorithm of optimal complexity. Also the discrete version of the solenoidal filtering approach (2a)–(2b) can be written in a similar manner. Using the above notation, we obtain

min u − uδ 2M + α u2K u

s.t. Bu = 0,

(17)

and the optimality system for this constrained minimization problem reads

(M + α K)u + B p = Muδ

(18)

Bu = 0.

(19)

The structure of this system is similar to that of the state system (12a)–(12b), which is why we denote the Lagrange multiplier for the divergence constraint again by p here. Existence and uniqueness of a solution is guaranteed for all α ≥ 0, if inf-sup stable finite elements are used for the discretization of u and p. Again, iterative methods with multigrid preconditioning may be applied for the efficient solution [36,37]. 4.4. Summary As can be seen from the discussion above, all linear filtering approaches considered in this paper can be discretized systematically and in a uniform framework by finite element methods. The resulting optimality systems can then always be solved efficiently by iterative solvers and appropriate preconditioning techniques. All resulting filters can therefore be considered to be algorithms of optimal complexity. 5. Computational results In order to illustrate the properties of the proposed reconstruction approach, we now present some preliminary computational results. For ease of presentation, we only consider a simple two-dimensional test problem. In all simulations, we use a regular triangulation of the computational domain and we assume that the measured velocity field uδ is given at each vertex of the mesh. For the discretization of the flow equations, we use here the Mini element [30,38]. Other inf-sup stable finite elements, in particular, such leading to exactly divergence free discrete velocity fields, could however be used as well. 5.1. A test problem for channel flow Let us consider the steady laminar flow between two parallel plates discussed in Section 2.5. As computational domain, we choose  = (0, L ) × (0, H ) with boundaries ∂in = {0} × (0, H ), ∂out = { L } × (0, H ), and ∂ wall = (0, L ) × {0, H }. For our simulations, we set H = 1 and L = 5. The model parameters are set to ν = 0.01, p 0 = 1, and p L = 0. The exact solution of (8a)–(8b) given by the Poisseuille law (9) then reads

p† (x, y ) = 1 − x/5

and

u† (x, y ) = (10 y (1 − y ), 0),

and the corresponding right hand side and boundary data are given by

f† (x, y ) = (0, 0),

g† (0, y ) = (10 y (1 − y ), 0) ,

and

h† (5, y ) = (−50 y 2 (1 − y )2 , 0).

These functions will serve as the reference solution and as the model data for our computational tests.

H. Egger et al. / Journal of Computational Physics 344 (2017) 558–574

567

Fig. 1. Linearization errors u − u† H1 () (left) and p − p†  L 2 () (right) for different values of the noise level δ and various mesh sizes h. We obtain O (δ) convergence for the error in the velocity and in the pressure, as stated in Theorem 2.2. Some saturation due to discretization errors can be observed for large mesh size h.

5.2. Linearization error With our first test, we would like to illustrate the estimates for the linearization error given in Theorem 2.2. To do so, we construct perturbed data uδ by adding random noise to u† such that uδ − u† L3 () = δ , and then compute the solution

(u, p) of the linearized problem (5a)–(5e) with data f = 0, g = g† , and h = h† by the finite element method outlined above. The resulting errors are displayed in Fig. 1. The estimate of Theorem 2.2 reads u − u†  H 1 + p − p†  L 2 ≤ C (δ + data errors). Note that on coarse meshes, the discretization error also contributes to the data error and we therefore observe a certain saturation phenomenon as δ goes to zero. A similar behavior would be obtained in the presence of model errors other than the discretization error. On refined meshes, the numerical results reveal the expected O (δ) convergence of the total linearization error as predicted by our theory. 5.3. Verification of the estimates for the reconstruction error Let us next illustrate the two estimates of Theorem 3.2. To do so, we compute approximations for the minimizers of (10a)–(10b) by solving the discretized optimal control problem (13a)–(13b). The tests are repeated for different values of the noise level δ . Following the remarks in Section 3.3, we should obtain

uα − u† L2 () ≤ C δ

and

√ uα − u† H1 () ≤ C (δ + δ/ α )

when setting f∗ = f† , g∗ = g† , and h∗ = h† . The first estimate does not depend on α while the second predicts a blow-up with α → 0, which is a manifestation of the ill-posedness of the underlying data smoothing problem. In the absence of model errors, we should however obtain errors in the size O (δ) in both, the L 2 - and the H 1 -norm, even if α is chosen in the order of one. To avoid the influence of discretization errors, we choose a rather fine mesh with mesh size h = 0.032 and nt = 17920 triangles for all computations. To further evaluate the influence of errors in the model data, the tests are repeated with parameters f∗ = f† chosen randomly and such that f∗ − f† L2 () = 5. The corresponding results are summarized in Fig. 2. As predicted by our analysis, the data residuals uα − uδ L2 () decrease monotonically in α and δ , while the reconstruction errors uα − u† H1 () + pα − p†  L 2 () show a semi-convergence behavior. These observations are in good agreement with the results of Theorem 3.2 and the first assertion of Theorem 3.4. 5.4. Parameter choice by the discrepancy principle The monotonicity of the Tikhonov functional stated in Theorem 3.4 allows to show that also the data residuals uα − uδ 2L2 () are monotonically decreasing with α → 0, which can also be seen in the two plots in the first line of Fig. 2. This motivates the use of a discrepancy principle for the choice of the regularization parameter in order to achieve an automatic balance of the error contributions due to data noise and model errors. Using the standard procedure [6], we define

αdis (δ) = max{α = α0 2−k : uα − uδ L2 () ≤ τ δ, k ∈ N}, where α0 > 0 and τ > 1 are given. For our computations, we choose α0 = 1 and τ = 1.2. In Table 1, we list the parameters selected by this strategy and the corresponding data residuals and reconstruction errors that are obtained for the test case

568

H. Egger et al. / Journal of Computational Physics 344 (2017) 558–574

Fig. 2. Convergence of the residuals uα − uδ L2 () and the reconstruction errors uα − u† H1 () + pα − p†  L 2 () with regularization parameter Left: exact model data f† − f∗ L2 () = 0; right: perturbed model data f† − f∗ L2 () = 5.

α → 0.

Table 1 Data residuals uα − uδ L2 () and total reconstruction errors uα − u† H1 () + pα − p†  L 2 () for regularization parameter αdis chosen by the discrepancy principle and the choice αopt minimizing the total reconstruction error.

δ

αdis

Residual

Error

αopt

Residual

Error

0.1000 0.0500 0.0250 0.0125 0.0063

1.56e-02 3.91e-03 9.77e-04 2.44e-04 3.81e-06

0.1143 0.0585 0.0292 0.0145 0.0074

2.7189 1.9909 1.4758 1.1311 0.7341

1.95e-03 4.88e-04 1.22e-04 1.53e-05 3.81e-06

0.0864 0.0398 0.0194 0.0105 0.0074

2.2246 1.6092 1.2144 0.9159 0.7341

corresponding to the right plot in Fig. 2. Note that the model data were improperly specified here, i.e., f∗ = f† in this example. For comparison, we also display the results for the values αopt of the regularization parameter for which the reconstruction errors are minimized. The discrepancy principle typically selects the regularization parameter αdis somewhat larger than the optimal value αopt . Repeating the tests with exact model data f∗ = f† , g∗ = g† , and h∗ = h† , we obtain αdis ≈ constant independent of the noise level δ . A similar behavior is observed for αopt . This is in agreement with the first estimate in Theorem 3.2 and can be seen in the plots on the left side of Fig. 2. 5.5. Comparison with other filters To further evaluate the overall performance of the proposed filter, let us make also a short comparison with the other filtering approaches discussed in the introduction. As outlined in Section 4.1, these can be implemented in a similar manner as the method presented in this paper which allows a fair comparison.

H. Egger et al. / Journal of Computational Physics 344 (2017) 558–574

569

Table 2 Reconstruction errors for velocity and pressure and discrete divergence of the reconstructed velocity field for different filter strategies at noise level δ = 0.1 on a mesh 9153 vertices, 17920 elements and mesh size h = 0.032. Filtering method

u − u†  L 2

u − u†  H 1

p − p†  L 2

∇h · u L 2

Smoothing Solenoidal (α = 0) Solenoidal w. smoothing Fluid-dyn. consistent

0.07960 0.08110 0.07581 0.05704

3.5893 20.8378 3.4003 3.0351

*.**** 1.2901 1.2901 0.1251

0.9986 0.0000 0.0000 0.0000

Fig. 3. Top view (left) with computational domain (red box) and front view (right) of the geometry. It consists of 8 arrays of displaced columns with diameter 10 mm, height 20 mm and 25 mm distance to each other. Measurements from the rectangular area [3, 148] × [50, 100] mm2 in red are later used for the reconstruction. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

In the following tests, we compare the smoothing filter (1), the solenoidal filtering (2a)–(2b) with and without smoothing, and the fluid-dynamically consistent filter (10a)–(10b) proposed in this paper. Whenever required, the regularization parameter α is selected via the discrepancy principle with τ = 1.2. In Table 2, we list different measures for the reconstruction error for the various choices of the filter. The discrete divergence ∇h · u here corresponds to the projection of ∇ · u onto the discrete pressure space. All filters yield small errors for the velocity in the L 2 -norm which is an immediate consequence of their construction. The H 1 -norm of the error in the velocities is comparable for the filters involving some sort of smoothing. For the solenoidal filter without smoothing, the H 1 -norm errors increase with decreasing mesh size due to the ill-posedness of the reconstruction problem. The smoothing filter yields a velocity reconstruction which is not discrete divergence free and no information about the pressure is obtained. The three filters involving a divergence constraint yield some reconstruction of the pressure. Those obtained with the fluid-dynamically consistent filter are however at least an order of magnitude better than those obtained with the solenoidal filters. In summary, the fluid-dynamically consistent filter proposed in this paper yields the best reconstruction of the flow fields with respect to all error measures listed in the table. 6. Application to experimental data To demonstrate the performance of the proposed filtering strategy in a realistic setting, we now apply our algorithms to experimental data. We consider the flow of water through a pin-fin array, i.e., a rectangular channel which is penetrated by cylindrical obstacles. A cross section of the geometry is depicted in Fig. 3. This test case has been formulated as a benchmark problem at the 15th ERCOFTAC SIG15 Workshop on Turbulence Modeling on October 2011 in EDF at Chatou in France [39,40]. For our computations, we use flow measurements obtained by magnetic resonance velocimetry (MRV) at the University Medical Center in Freiburg [41]. The inflow velocity in the experiments was approximately U = 360 mm/s. With a kinematic viscosity of ν = 0.6 mm2 /s and a length scale of L = 20 mm, corresponding to the height of the channel, this experimental setup amounts to a Reynolds number of about Re = 12.000. Since the measurement process takes some time, the recorded flow data actually correspond to a time averaged velocity field. The full measurements consist of 329 × 211 × 31 voxels of size 0.6534 × 0.6534 × 0.70 mm3 for each of the three velocity components. A brief inspection of the flow data revealed that the flow is almost two dimensional, i.e., the vertical flow component u z of the velocity is almost zero and the horizontal components u x and u y are almost independent of the height z in the channel, except in close vicinity of the top and bottom walls. For illustration, we depict in Fig. 4 the data for the velocity u x in flow direction at the inflow boundary for different levels of the vertical position z. As can be seen from the plot, the velocity measurements of u x are rather noisy but almost constant across the inflow boundary. It thus seems reasonable to assume a two dimensional flow away from the walls of the channel.

570

H. Egger et al. / Journal of Computational Physics 344 (2017) 558–574

Fig. 4. Inflow velocity u x for different values of height z averaged over ten voxel layers in z-direction. Inflow velocity can bee seen to be almost constant within the computational area, but still rather noisy even after averaging.

Fig. 5. Velocity data and computational mesh used for all computations. The colored boundaries represent areas of different boundary conditions, i.e. inflow (green), outflow (magenta), no-slip (red) and periodic (blue). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

6.1. Data preprocessing and computational domain For our numerical tests, we use a two dimensional flow model and only utilize data of the horizontal velocity components u x and u y from a single slice at the mid plane of the channel at z = 10 mm and from the area [3, 148] ×[50, 100] mm2 depicted as red rectangle in Fig. 3. The geometry and mesh used for our computations as well as the flow data uδ after interpolation to the finite element mesh are shown in Fig. 5. 6.2. Choice of boundary conditions and model data No slip conditions were prescribed for the boundaries of all obstacles and periodic boundary conditions were used at the top and bottom boundary of the rectangular area serving as computational domain. The in- and outflow boundaries at the left and right end were treated similarly as in the examples presented in the previous section. Due to lack of better knowledge, we chose f∗ = 0 and h∗ = 0, and we utilized the H1 -seminorm ∂ y g2L2 (∂ ) for regularization of the inflow data in

thus promoting a constant inflow profile. The viscosity parameter ν was finally chosen as follows: Since the data correspond to time averaged velocity fields, the stationary flow problem should actually be considered as a Reynolds averaged Navier–Stokes (RANS) model [42]. This means that the viscosity ν in the model should have the form ν = ν0 + νt with ν0 denoting the physical and νt the artificial turbulent viscosity. Note that in practice νt is often much larger than ν0 . To estimate the total viscosity parameter ν , we proceeded as follows: We computed the residuals u − uδ  by solving the forward problem (5a)–(5e) for different values of ν and also computed the residuals when applying our filter for different values of the viscosity ν and the regularization parameter α . From the results depicted in Fig. 6 one can see that a choice ν = 80 mm2 /s leads to small residuals for all values of α , and we decided to use this value for all further computations. Let us emphasize that the total viscosity is about hundred times larger then the kinematic viscosity ν0 = 0.6 mm2 /s of water here. A good estimate for turbulent viscosity νt thus is an important ingredient for the fluid dynamic consistent filter.

H. Egger et al. / Journal of Computational Physics 344 (2017) 558–574

Fig. 6. Residual for different values of

571

ν in mm2 /s with αν being constant.

Fig. 7. Velocity reconstruction of the fluid dynamic consistent filter with

ν = 80 mm2 /s for α = 1, 10−5 , and 10−12 .

6.3. Results for the fluid-dynamically consistent filter We are now ready to apply the proposed filter to the experimental data. In Fig. 7, we display the velocity reconstructions obtained with the fluid dynamic consistent filter for different choices of the parameter α . Let us note that the reconstructions are very stable with respect to the choice of the regularization parameter α , which is in perfect agreement with the a-priori estimates derived in Theorem 3.2. Even if α is chosen to be one, the reconstruction shows the relevant flow features, e.g., the recirculation zones behind the obstacles. Only for very small α , the data noise becomes clearly visible in the reconstructions. The reconstructions are therefore very robust with respect to the choice of the regularization parameter; see Theorem 3.2. Also recall that the no-slip boundary condition at the obstacles and the discrete divergence free condition are realized exactly for all values of the regularization parameter α > 0. 6.4. Comparison to the solenoidal filter and pressure reconstruction In order to compare the performance of the proposed filter, we also computed reconstructions using the solenoidal filter with smoothing outlined in the introduction. The corresponding results are depicted in Fig. 8. Note that the results are much more sensitive to the choice of the regularization parameter here. For large α we obtain over smoothed reconstructions,

572

H. Egger et al. / Journal of Computational Physics 344 (2017) 558–574

Fig. 8. Reconstruction with divergence-free (solenoidal) filter for

Fig. 9. Pressure distribution for optimal

α = 30, α = 0.5 and α = 0.001.

α for the fluid-dynamically consistent filter and the divergence free filter with smoothing.

and small values of the regularization parameter result in noisy reconstructions. A rather accurate choice of α is required to obtain a comparable reconstruction quality as with the fluid dynamic consistent filter. We next compare the two filters also with respect to another quality measure: Both, the fluid dynamic consistent and the solenoidal filter, automatically generate an estimate for the pressure p. In Fig. 9, we depict the corresponding reconstructions. Note that the solenoidal filter does not provide any reasonable information about the physical pressure, while the pressure reconstruction obtained with the fluid dynamic consistent filter has the expected behavior: a linear pressure drop along the flow and stagnation respectively depletion zones before and behind each obstacle. As mentioned previously, the measurements and the reconstructions obtained with the stationary flow model amount to Reynolds averaged flow fields. The reconstructed pressure therefore also contains contributions of the turbulent kinetic energy [42]. A further post-processing would thus be required to gain quantitative information about the actual pressure drop along the channel. 6.5. Reconstruction in three dimensions From our arguments at the beginning of the section, we would expect to obtain similar reconstructions when using three dimensional data and a corresponding flow model. To demonstrate this, we repeated the computations for the solenoidal filter in three dimensions, now taking into account the measurements of all three velocity components in the box [3, 148] × [50, 100] × [5, 15] mm3 . In Fig. 10, we depict the reconstructed horizontal velocity components at the mid plane z = 10 mm corresponding to the area of the two dimensional reconstructions. A quick comparison with the plots of Fig. 8 shows that the results are in fact almost indistinguishable from those obtained in the two dimensional case. The third velocity component was two orders of magnitude smaller than the other components and is therefore not shown here. These results clearly

H. Egger et al. / Journal of Computational Physics 344 (2017) 558–574

573

Fig. 10. Velocity reconstructions for the horizontal velocity components at the central slice z = 10 mm obtained with the solenoidal filter (2a)–(2b) with regularization parameters α = 30, α = 0.5 and α = 0.001.

demonstrate that the assumption of an essentially two dimensional flow field was appropriate for the data acquired in this experiment. 7. Discussion In this paper we considered the reconstruction of the velocity and pressure fields of an incompressible fluid from distributed measurements of the flow velocities. For the stable solution of this inverse problem, we proposed a novel filter which minimizes a weighted sum of the data residual and the mismatch of a specified flow model. This strategy was formulated as an optimal control problem constrained by a flow model. In order to guarantee the well-posedness of our approach, we utilized a linearized flow model which directly incorporates the measured velocity field. This allowed us to show the existence and uniqueness of minimizers and to derive estimates for the reconstruction errors in various norms. The theoretical results were illustrated by numerical tests including a comparison to other filters discussed in the literature. The applicability to realistic problems was demonstrated by application of the proposed method to experimental data obtained by magnetic resonance velocimetry for a benchmark problem. The strategy of using a linearized flow model as a constraint in the reconstruction process could be generalized in various ways: While we directly used the velocity measurements in order to specify our linearized flow model here, some pre-filtered velocity field could be used as well. Our analysis also covers this case and could possibly be refined leading to sharper estimates. Repeating the argument, one could also define an incremental reconstruction approach. Preliminary numerical tests for such multi-step algorithms showed a further significant improvement of the reconstructed flow fields. A full analysis would however exceed the scope of the current presentation. In order to handle more general flow regimes, some sort of turbulence model should be incorporated as a next step and a refined modeling of the constitutive equations and the boundary conditions should be considered. Both aspects are subject of current research. Acknowledgements The authors would like to gratefully acknowledge the support by the German Research Foundation (DFG) via grants IRTG 1529, GSC 233, and Eg-331/1-1. Part of the work of the second author was carried out during a research stay at Waseda University, Tokyo. The hospitality and kind support of Waseda university is gratefully acknowledged.

574

H. Egger et al. / Journal of Computational Physics 344 (2017) 558–574

References [1] C.J. Elkins, M. Markl, N. Pelc, J.K. Eaton, 4D Magnetic resonance velocimetry for mean velocity measurements in complex turbulent flows, Exp. Fluids 34 (2003) 494–503. [2] C.J. Elkins, M.T. Alley, Magnetic resonance velocimetry: applications for magnetic resonance imaging in the measurement of fluid motion, Exp. Fluids 43 (2007) 823–858. [3] S.F. Herrmann, K. Hinsch, Holographic particle image velocimetry, Meas. Sci. Technol. 15 (2004) 1–9. [4] F. Pereira, M. Gharib, D. Dabiri, D. Modarress, Defocussing digital particle image velocimetry: a 3-component 3-dimensional DPIV measurement technique. Application to bubbly flows, Exp. Fluids 29 (2000) 78–84. [5] J.A. Liburdy, E.F. Young, Processing of three-dimensional particle tracking velocimetry data, Opt. Lasers Eng. 17 (1992) 209–227. [6] H.W. Engl, M. Hanke, A. Neubauer, Regularization of Inverse Problems, Kluwer Academic Publishers, Dordrecht, 1996. [7] A.N. Tikhonov, Solution of incorrectly formulated problems and the regularization method, Sov. Math. Dokl. 5 (1963) 1035–1038. [8] O. Scherzer, M. Grasmair, H. Grossauer, M. Haltmeier, F. Lenzen, Variational Methods in Imaging, Applied Mathematical Sciences, vol. 167, Springer, New York, 2009. [9] J. Busch, D. Giese, L. Wissmann, S. Kozerke, Reconstruction of divergence-free velocity fields from cine 3D phase-contrast flow measurements, Magn. Reson. Med. 69 (2013) 200–210. [10] C.M. de Silva, J. Philip, I. Marusic, Minimization of divergence error in volumetric velocity measurements and implications for turbulence statistics, Exp. Fluids 54 (2013) 1557. [11] I. Macedo, R. Castro, Learning Divergence-Free and Curl-Free Vector Fields with Matrix-Valued Kernels, Technical report, Instituto Nacional de Matematica Pura e Aplicada, 2008. [12] C.P. McNally, Divergence-free interpolation of vector fields from point values – exact ∇ · b = 0 in numerical simulations, Mon. Not. R. Astron. Soc. 413 (2011) L76–L80. [13] F. Ong, M. Uecker, U. Tariq, A. Hsiao, M.T. Alley, S.S. Vasanawala, M. Lustig, Robust 4D flow denoising using divergence-free wavelet transform, Magn. Reson. Med. 73 (2015) 828–842. [14] D. Schiavazzi, F. Coletti, G. Iaccarino, J.K. Eaton, A matching pursuit approach to solenoidal filtering of three-dimensional velocity measurements, J. Comput. Phys. 263 (2014) 206–221. [15] V. Girault, R.-A. Raviart, Finite Element Approximation of the Navier–Stokes Equations, Lecture Notes in Mathematics, vol. 749, Springer-Verlag, Berlin, 1979. [16] M. Desai, K. Ito, Optimal controls of Navier–Stokes equations, SIAM J. Control Optim. 32 (1994) 1428–1446. [17] M. Gunzburger, L. Hou, T. Svobodny, Analysis and finite element approximation of optimal control problems for the stationary Navier–Stokes equations with Dirichlet controls, Modél. Math. Anal. Numér. 25 (1991) 711–748. [18] R. Herzog, K. Kunisch, Algorithms for pde-constrained optimization, GAMM-Mitt. 33 (2010) 163–176. [19] L. Hou, T. Svobodny, Optimization problems for the Navier–Stokes equations with regular boundary control, J. Math. Anal. Appl. 177 (1993) 342–367. [20] F. Abergel, R. Temam, On some control problems in fluid mechanics, Theor. Comput. Fluid Dyn. (1990) 303–325. [21] H. Fattorini, S. Sritharan, Existence of optimal controls for viscous flow problems, Proc. R. Soc. Lond. Ser. A 439 (1992) 81–102. [22] A. Fursikov, Control problems and theorems concerning the unique solvability of a mixed boundary value problem for the three-dimensional Navier– Stokes and Euler equations, Math. USSR Sb. 43 (1982) 281–307. [23] M. Hinze, K. Kunisch, Second order methods for optimal control of time-dependent fluid flow, SIAM J. Control Optim. 40 (2001) 925–946. [24] M. Gunzburger, L. Hou, T. Svobodny, Boundary velocity control of incompressible flow with an application to viscous drag reduction, SIAM J. Control Optim. 30 (1992) 167–181. [25] C. Foias, R. Temam, Structure of the set of stationary solutions of the Navier–Stokes equations, Commun. Pure Appl. Math. 30 (1977) 149–164. [26] R. Temam, Navier–Stokes Equations. Theory and Numerical Analysis, third edition, Studies in Mathematics and Its Applications, vol. 2, North-Holland Publishing Co., Amsterdam, 1984. [27] H.T. Banks, K. Kunisch, Estimation Techniques for Distributed Parameter Systems, Systems & Control: Foundations & Applications, vol. 1, Birkhäuser, 1989. [28] M. Hanke, O. Scherzer, Error analysis of an equation error method for the identification of the diffusion coefficient in a quasi-linear parabolic differential equation, SIAM J. Appl. Math. 59 (1999) 1012–1027. [29] G.K. Batchelor, An Introduction to Fluid Dynamics, Cambridge University Press, Cambridge, 1967. [30] D. Boffi, F. Brezzi, L.F. Demkowicz, R.G. Durán, R.S. Falk, M. Fortin, Mixed Finite Elements, Compatibility Conditions, and Applications, Lecture Notes in Mathematics, vol. 1939, Springer-Verlag, Berlin, 2008. [31] A. Draganescu, A.M. Soane, Multigrid solution of a distributed optimal control problem constrained by the Stokes equations, Appl. Math. Comput. 219 (2013) 5622–5634. [32] J. Schöberl, W. Zulehner, Symmetric indefinite preconditioners for saddle point problems with applications to PDE-constrained optimization problems, SIAM J. Matrix Anal. Appl. 29 (2007) 752–773. [33] S. Takacs, A robust all-at-once multigrid method for the Stokes control problem, Numer. Math. 130 (2015) 517–540. [34] H. Egger, M. Schlottbom, Efficient reliable image reconstruction schemes for diffuse optical tomography, Inverse Probl. Sci. Eng. 19 (2011) 155–180. [35] J. Schöberl, Multigrid methods for a parameter dependent problem in primal variables, Numer. Math. 84 (1999) 97–119. [36] J. Schöberl, Robust multigrid preconditioning for parameter-dependent problems. I. The Stokes-type case, in: Multigrid Methods V, Stuttgart, 1996, in: Lect. Notes Comput. Sci. Eng., vol. 3, Springer, Berlin, 1998, pp. 260–275. [37] W. Zulehner, Analysis of iterative methods for saddle point problems: a unified approach, Math. Comput. 71 (2002) 479–505. [38] D. Boffi, F. Brezzi, M. Fortin, Mixed Finite Element Methods and Applications, Springer Series in Computational Mathematics, vol. 44, Springer, Heidelberg, 2013. [39] F.E. Ames, L.A. Dvorak, M.J. Morrow, Turbulent augmentation of internal convection over pins in staggered-pin fin arrays, J. Turbomach. 127 (1) (2005) 183–190. [40] F.E. Ames, L.A. Dvorak, Turbulent transport in pin fin arrays: experimental data and predictions, J. Turbomach. 128 (1) (2005) 71–81. [41] F. Wassermann, S. Grundmann, Flow visualization of a staggered pin fin array using magnetic resonance velocimetry, in: 15th ERCOFTAC-SIG15/IAHR Workshop, October 2011. [42] G.K. Batchelor, The Theory of Homogeneous Turbulence, Cambridge University Press, 1953.