Cryo-STEM Tomography of Intact Vitrified Fibroblasts

1.   Introduction

Optimal control problems (OCPs) for partial differential equations is a research area that features challenging theoretical questions for mathematical analysis, the need of devising efficient numerical algorithms, with a broad range of meaningful applications in engineering, like, e.g., aerodynamics ^[29,30] and biofluid dynamics ^[33,40]. In this paper we address one of the most classical problems in this field, namely the boundary control for the Navier-Stokes equations. Much has been done for this problem, both on the mathematical analysis and on its numerical approximation; there are, however, some issues that have not yet been addressed – the case of the state problem with mixed Dirichlet and Neumann conditions is one of those. At the best of our knowledge, a comprehensive framework to address both ($i$) the well posedness analysis for the OCP related to this system and ($ii$) the derivation of a system of optimality conditions, is still missing. Setting this framework, and exploiting it in view of the numerical approximation of the OCP, is the scope of our paper.

More precisely, we study the mixed boundary value problem for the incompressible steady Navier-Stokes equations, both in two and three dimensions, where we prescribe homogeneous Neumann conditions (on the normal stress) on a portion of the boundary, homogeneous Dirichlet conditions on the velocity field on another boundary portion, and a velocity control on the remaining portion of the domain boundary^*. For the sake of application, we consider the minimization of quadratic cost functionals involving either the square of the $ L^2 $ norm of the vorticity, or a velocity tracking-type term, on the whole computational domain. For the problem at hand, we introduce an ad-hoc saddle point formulation of the state problem which, besides velocity and pressure fields, requires a third variable representing the normal Cauchy stress and a Lagrange multiplier penalizing the Dirichlet data (this latter involving the control function, too). This saddle point formulation is more apt for numerical approximation than the one with a single field $ u $ that is obtained after projecting the Navier-Stokes equations onto a divergence free space. For the control problem, the saddle point formulation enables us to derive a system of first-order optimality conditions, and to carry out a rigorous mathematical analysis of the OCP.

Splitting the boundary $ \partial \Omega $ into portions on which different boundary conditions are set is rather common in several fluid dynamics applications. Dirichlet conditions are used for the problem at hand to impose a velocity profile on the inflow boundary and on the control boundary, and a no-slip condition at the lumen boundary, since it is assumed that there fluid particles perfectly adhere to the vessel wall. Neumann conditions prescribe a force per unit area as the normal component of the stress tensor; if the force is vanishing, the Neumann boundary portion is a free outflow – this is indeed the condition we assign on the outflow boundary. Note also that when Dirichlet boundary conditions are specified on the entire boundary $ \partial \Omega $, the pressure is not uniquely defined.

Similar problems have been studied by several authors. What follows is an account of the most relevant works on the subject. In a classical paper by Gunzburger ^[21], the same control problem is addressed with however full Dirichlet boundary conditions (and a different cost functional), by a saddle-point formulation; a similar approach has been applied to viscous drag reduction in ^[19]. In a series of other papers the problem is faced working on the subspace of divergence free velocity fields. In ^[12], the unsteady problem is studied, with Dirichlet conditions on the whole boundary; in ^[13], an exterior problem in two dimensions for the steady Navier-Stokes equations with boundary Neumann control is considered for the minimization of a drag functional; an inequality state constraint is also enforced on the velocity field. Distributed and Neumann controls are also considered, e.g., in ^[20]. Boundary control problems for vorticity minimization have been analyzed in ^[32], in the case of steady Stokes equations, and in ^[31], considering instead unsteady Navier-Stokes equations; in both cases, however, a pure Dirichlet problem is considered. In ^[17] the analysis of existence of optimal boundary control for the Navier-Stokes equations with mixed boundary conditions is carried out in three dimensions, without deriving a system of optimality conditions. Regarding the case of unsteady Navier-Stokes equations, a boundary OCP similar to the one studied in our paper has been analyzed in ^[24], however dealing with a pure Dirichlet problem. The analysis and the approximation of OCPs involving distributed controls and in the case of unsteady Navier-Stokes flows has been considered in ^[22,23], again dealing with full Dirichlet boundary conditions. Mixed Dirichlet-Neumann boundary conditions have been instead considered in ^[4], however without carrying out the well-posedness analysis. We believe that our presentation is more comprehensive, and allows the reader to have a better picture of the whole mathematical analysis for both the state and the control problem and the derivation of a system of optimality conditions.

The structure of the paper is as follows. In Section 2 we introduce a ad-hoc saddle point formulation of the state problem with a Lagrange multiplier that represents the normal Cauchy stress on the whole boundary and carry out its well-posedness analysis in Section 3; some technical results related to the Stokes solver are reported in Section 4; In Section 5 we set up the ground for the analysis of the control problem; we then investigate the differentiability of the control-to-solution map in Section 6, and derive first order optimality conditions for the OCP in Section 7; In Section 8 we illustrate a possible way to approximate our control problem, by relying on a descent iterative method for the minimization of the cost functional, and using the Galerkin finite element method to discretize both the state and the adjoint problem; Numerical results shown in Section 9 concern the solution of an interior flow problem, involving a simplified two-dimensional tract of blood vessel where a bypass is inserted. The goal is to regularize the flow by minimizing a tracking-type quadratic functional, expressing the distributed observation of the (square of) the $ L^2 $ norm of the difference $ {\bf v} - {\bf z}_d $ between the fluid velocity $ {\bf v} $ and the target velocity $ {\bf z}_d $, this latter being provided by the solution of the corresponding Stokes problem over the whole domain. Conclusions will then follow in Section 10.

2.   Problem formulation

In this paper we consider a state system given by the steady Navier-Stokes equations, in velocity and pressure formulation, for steady, incompressible and viscous fluid flows. For this state system, we analyze a minimization problem for quadratic functionals of the velocity field through a velocity control operated on a portion of the inflow boundary. For a concrete example we denote by $ \Omega \subset \mathbb{R}^{d} $, $ d = 2, 3 $ the domain in Figure 1 below, occupied by the fluid and by $ \Gamma = \partial \Omega $ its boundary. The inflow boundary $ \Gamma_{1} = \Gamma _{c}\cup \Gamma _{in} $ is composed by $ \Gamma _{c} $, the control boundary, and by $ \Gamma _{in}, $ where a Dirichlet data is assigned; $ \Gamma _{out} $ and $ \Gamma _{w} $ are the outflow and the wall boundary portions, where zero stress and zero velocity are assigned, respectively.

Let

denote the state space, endowed with the norm $ \Vert \mathbf{v}\Vert _{X} = \Vert \nabla \mathbf{v}\Vert _{0} $ where

Our goal is to minimize either the vorticity functional

or a velocity tracking-type cost functional

where $ \tau \geq 0 $, $ \mathbf{u}\in \mathcal{U} = H_{0}^{1}\left(\Gamma _{c}\right) ^{d} $ is a control variable and $ (\mathbf{v}, \pi) = (\mathbf{v}(\mathbf{u}), \pi (\mathbf{u}))\in X\times L^{2}\left(\Omega \right) $ is a weak solution of the following system:

where $ \mathbf{v} $ and $ \pi $ denote the fluid velocity and pressure, respectively, $ \nu $ is a given kynematic viscosity, and $ {\bf n} $ is the unit outward normal direction. In the case of $ \tilde{J}_2 $, a flow matching problem is solved, aiming at driving the velocity field toward a target velocity profile $ {\bf z}_d \in X $; this latter can be given, e.g., by the Stokes flow in the same domain, whenever the goal is to minimize recirculation and vortex generation for the sake of regularization of flows.

If $ \tau > 0 $, the second term in (2.1) or (2.2), called penalization term, besides helping in the well-posedness analysis, prevents using "too large" controls in the minimization of the cost functional.

We handle the Dirichlet boundary conditions on $ \Gamma _{1} $ (involving the control function) by a Lagrange multiplier approach. Let $ e $ and $ b $ be the following bilinear forms

and $ c $ be the trilinear form

Using Hölder inequality, and the Sobolev inequality

we have that

Let us now define

assuming that $ \mathbf{v}_{in}\in H_{00}^{1/2}\left(\Gamma _{in}\right) ^{d} $, we have that $ \mathbf{g}\in H_{00}^{1/2}(\Gamma _{1})^{d} $, since $ \mathbf{u}\in H_{0}^{1}\left(\Gamma _{c}\right) ^{d}\subset \mathcal{T}_c $ with continuous injection, that is, for some positive constant $ c_{0} $,

For the definition of these functional spaces, see, e.g., [7,Vol. Ⅲ,Chap. Ⅶ]. Observe that if $ \mathbf{v}\in X $, then $ \mathbf{v}|_{\Gamma _{1}}\in H_{00}^{1/2}(\Gamma _{1})^{d} $.

For ease of notation, hereon we denote by $ \mathcal{T}_1 = H_{00}^{1/2}(\Gamma _{1})^{d} $ and by $ \left\langle \cdot, \cdot \right\rangle _{\Gamma _{1}} $ the duality pairing between $ \mathcal{T}_1 $ and its dual $ \mathcal{T}_1^{\ast} $. In the same way, we will denote by $ \mathcal{T}_c = H_{00}^{1/2}(\Gamma _{c})^{d} $, and by $ \left\langle \cdot, \cdot \right\rangle _{\Gamma _{c}} $ the duality pairing between $ \mathcal{T}_c $ and its dual $ \mathcal{T}_c^{\ast} $.

Let $ V = X\times L^{2}\left(\Omega \right) \times \mathcal{T}_1^{\ast} $. The weak formulation of (2.3) that we consider reads as follows: find $ (\mathbf{v}, \pi, \mathbf{t)}\in V $ such that

Note that a formal integration by parts in the first equation gives for the boundary stress $ \pi \mathbf{n}-\nu \frac{\partial \mathbf{v}}{ \partial {n}} $,

and

Thus, the first equation in (2.7) enforces in a natural way the homogeneous Neumann condition on $ \Gamma _{out} $.

By grouping together the last two equations in (2.7), we obtain the following equivalent formulation: find $ \mathbf{v}\in X $, $ (\pi, \mathbf{t)}\in Q = L^{2}(\Omega)\times \mathcal{T}_1^{\ast} $ such that

Upon defining the bilinear form $ \tilde{b}:X\times Q\rightarrow \mathbb{R} $ as

and the norm

we reformulate problem (2.9) as: find $ \mathbf{v}\in X, \left(\pi, \mathbf{t}\right) \in Q $ such that

with $ G\in Q^{\ast } $ given by

3.   Analysis of the state problem

Despite their ubiquitous use when modeling viscous incompressible flows, Navier-Stokes equations endowed with mixed boundary conditions have been seldom analyzed from a theoretical standpoint. If homogeneous Neumann conditions are set on a boundary portion, energy inequalities or a priori estimates of a weak solution to the Navier-Stokes system cannot be obtained because there is no control over the energy flux on this boundary portion. Some results can be found, e.g., in ^[41] regarding the existence and uniqueness of solution in weighted $ L^p $ Sobolev spaces and Hölder spaces for problems set over polyhedral domains; see also, e.g., ^[3,35,36] for the time-dependent case and ^[2] for a well-posedness analysis in divergence-free spaces. Additional conditions that estimate backward flows over the Neumann boundary have been considered, e.g., in ^[34].

In this paper we derive a well-posedness result for the steady Navier-Stokes system in presence of mixed Dirichlet-Neumann boundary conditions including the pressure field among the state variables, in non-weighted Sobolev spaces. We show that if $ \left\Vert \mathbf{g}\right\Vert _{\mathcal{T}_1} $ is small enough, then (2.10) is well posed. We first need to introduce the following continuous operators, related to the bilinear form $ \tilde{b} $: $ B:X\rightarrow Q^{\ast } $ defined by

and $ B^{\ast }:Q\mapsto X^{\ast } $ defined by

We denote by $ N_{B} $ the kernel of $ B, $ by $ N_{B}^{\bot } $ its orthogonal (in $ X $) and by $ N_{B}^{0} $ the polar of $ N_{B} $, that is

The so called inf-sup condition holds for the bilinear form $ \tilde{b} $:

Lemma 3.1. There exists a constant $ \tilde{\beta} > 0 $ such that

Proof. Given $ (\tilde{\phi}, \tilde{\boldsymbol{\lambda}})\in Q $, according to Riesz Theorem there exists a unique $ \boldsymbol{\lambda}^{\ast }\in \mathcal{T}_1 $ such that

with $ \Vert \boldsymbol{\lambda}^{\ast }\Vert _{\mathcal{T}_1} = \Vert \tilde{\boldsymbol{\lambda}}\Vert _{ \mathcal{T}_1^{\ast}} $. Then we look for $ \tilde{\mathbf{v}}\in X $ that satisfies

with

Such $ \tilde{\mathbf{v}} $ can be obtained as follows. Let $ \mathbf{r}\in X $ be the lifting of

where

with $ \int_{\Gamma _{out}}\mathbf{g}_{0}\cdot \mathbf{n} = 1 $. Let us now set $ \tilde{\mathbf{v}} = \mathbf{w}+\mathbf{r} $, with $ \mathbf{w}\in H_{0}^{1}\left(\Omega \right) ^{d} $ to be determined. We have

that is, $ \mbox{div}\mathbf{w} = \tilde{\phi}-\mbox{div}\mathbf{r} $. Moreover,

It follows (see, e.g., ^[15]) that there exists $ \mathbf{w}\in H_{0}^{1}\left(\Omega \right) ^{d} $ with $ \mbox{div}\mathbf{w} = \tilde{\phi}- \mbox{div}\mathbf{r} $ so that $ \tilde{\mathbf{v}} = \mathbf{w}+\mathbf{r} $ meets the required assumptions. Thus we can write

from which $ \left(3.1\right) $ follows with $ \tilde{\beta} = \frac{1 }{2\tilde{C}} $, by taking the supremum over $ \tilde{\mathbf{v}} $ and the infimum over $ (\tilde{\phi}, \tilde{\boldsymbol{\lambda}}) $.

The following lemma is well known (see [15,Chap. 1,Thm. 4.1]).

Lemma 3.2. Condition (3.1) is equivalent to each one of the following properties:

(ⅰ) the operator $ B $ is a continuous isomorphism from $ N_{B}^{\bot } $ onto $ Q^{\ast } $ and

(ⅱ) the operator $ B^{\ast} $ is a continuous isomorphism from $ Q $ onto $ N_{B}^{0} $ and

We are ready to prove the following result.

Theorem 3.3. There exists $ \eta > 0, $ depending only on $ \nu $, $ \tilde{\beta}, C_{S} $, such that if $ \Vert \mathbf{g(u)}\Vert _{\mathcal{T}_1}\leq \eta $, then problem (2.10) has a unique solution $ \left(\mathbf{v, }\pi, \mathbf{t}\right) = \left(\mathbf{v} (\mathbf{u}), \pi (\mathbf{u}), \mathbf{t} (\mathbf{u})\right) {\in X \times L^{2}(\Omega) \times \mathcal{T}_1^{\ast}} $. Moreover,

and

where $ C_{1} = 5\tilde{\beta}^{-1} $ and $ C_{2} $ depend only on $ \nu $, $ \tilde{\beta}, C_{S} $.

Remark 3.4. A possible choice of $ \eta $ is

Proof. Write $ \mathbf{v = z+v}_{0} $ with $ \mathbf{v}_{0}\in N_{B}^{\bot } $, $ \mathbf{z}\in N_{B} $. Define

To solve problem (2.10), we first determine $ \mathbf{v}_{0} $ by solving

where $ G\in Q^{\ast } $ is defined in (2.11). Then we determine $ \mathbf{z} \in N_{B} $ from the equation

where $ F_{0}\in N_{B}^{\ast } $ is the restriction to $ N_{B} $ of $ F\in X^{\ast }, $ given by^†

^†Rigorously, $ F_{0} = F\circ P $ where $ P\in \mathcal{L(}X^{\ast }, N_{B}^{\ast }) $ is defined as $ \left\langle P\mathbf{f}, \boldsymbol{\varphi}\right\rangle _{N_{B}^{\ast }, N_{B}} = \left\langle \mathbf{f}, \boldsymbol{\varphi}\right\rangle _{X^{\ast}, X}\; \; \forall \mathbf{f}\in X^{\ast }, \boldsymbol{\varphi}\in N_{B} $. Clearly, $ \left\Vert F_{0}\right\Vert _{N_{B}{}^{\ast }}\leq \left\Vert F\right\Vert_{X{}^{\ast }} $.

Finally, with $ \mathbf{v = z+v}_{0} $, we find $ (\pi, \mathbf{t})\in Q $ from

Note that

and that $ \left\langle F, \boldsymbol{\psi}\right\rangle -\tilde{e}\left(\mathbf{z, z}, \boldsymbol{\psi}\right) = -e(\mathbf{v}, \boldsymbol{\psi})-c\left(\mathbf{ v, v}, \boldsymbol{\psi }\right), $ $ \boldsymbol{\psi}\in X $. We divide the proof in four steps.

Step 1. Equation (3.6) is equivalent to

which, by Lemma 3.2, ($i$), has a unique solution $ \mathbf{v}_{0}\in N_{B}^{\bot } $. Moreover, from (3.2) and (2.11),

Step 2. Now we reduce (3.7) to a fixed point equation. Let $ \mathbf{z}\in N_{B} $ be given and consider the linear equation

For each $ \mathbf{z} $, the bilinear form $ \tilde{e}\left(\mathbf{z; \cdot, \cdot }\right) $ is continuous, by (2.5). Moreover,

Thus, if

then

Therefore, assuming

with $ \eta $ as in (3.5), we have

and $ \tilde{e}\left(\mathbf{z; \cdot, \cdot }\right) $ is $ N_{B}- $coercive. Let us define

and, for $ m $ integer, define

By the Lax-Milgram Theorem, for each $ \mathbf{z}\in S_{4\delta _{\mathbf{u}}} $, we can define a solution map operator $ T (\mathbf{z}) \in \mathcal{L} (N_{B}^{\ast }, N_{B}) $ that to each $ L\in N_{B}^{\ast } $ associates the unique solution $ \mathbf{w = }T(\mathbf{z}) L $ of the equation

Moreover,

from which

If $ \mathbf{z}^{\ast } $ is a fixed point for the operator $ \mathbf{z}\mapsto T\left(\mathbf{z}\right) F_{0}, $ then

which is nothing but equation (3.7).

Step 3. We show that, if $ \left\Vert L\right\Vert _{N_{B}^{\ast }} $ is small enough, $ \mathbf{z}\mapsto T\left(\mathbf{z}\right) L $ is a strict contraction from $ S_{4\delta _{\mathbf{u}}} $ into itself. From (3.15), if

then $ \mathbf{z}\mapsto T\left(\mathbf{z}\right) L $ maps $ S_{4\delta _{\mathbf{u}}} $ into itself. We now estimate $ \left\Vert (T\left(\mathbf{z}_{1}\right) -T\left(\mathbf{z }_{2}\right))L\right\Vert _{N_{B}}. $ Thanks to the continuity of $ \tilde{e} $, we can rewrite equation (3.15) as

or

where $ A\left(\mathbf{z}\right) \in \mathcal{L} (N_{B}, N_{B}^{\ast }). $ Clearly $ T\left(\mathbf{z}\right) = A\left(\mathbf{z}\right) ^{-1} $ and therefore we have the identity

from which

Now

whence

and

We deduce that if

then $ \mathbf{z}\mapsto T\left(\mathbf{z}\right) L $ is a strict contraction and therefore has a unique fixed point in $ S_{\delta }. $

If $ L = F_{0} $, we have ($ \mathbf{g = g}\left(\mathbf{u}\right) $)

since

Moreover,

since

Thus, $ T\left(\mathbf{z}\right) F_{0} $: $ S_{4\delta _{\mathbf{u} }}\rightarrow S_{4\delta _{\mathbf{u}}} $ is a strict contraction and has a unique fixed point $ \mathbf{z}^{\ast }\in N_{B} $ with

Step 4. It remains to solve equation (3.8). Observe that equation (3.7) means that the operator $ E\left(\mathbf{v}\right) \in X^{\ast } $ given by

where $ \mathbf{v} = \mathbf{z}^{\ast }+\mathbf{v}_{0}, $ belongs to $ N_{B}^{0} $. Therefore, by Lemma 3.2, $ \left(ii\right) $, there exists a unique $ (\pi, \mathbf{t})\in Q $ such that

with

and the proof is complete with $ C_{1} = 5\tilde{\beta}^{-1} $ and $ C_{2} $ depending only on the constants $ \nu $, $ \tilde{\beta}, C_{S} $.

Remark 3.5. Theorem 3.3 holds with a right hand side $ \mathbf{f}\in X^{\ast } $ with small enough $ \left\Vert \mathbf{f} \right\Vert _{X^{\ast }} $ norm. Then the unique solution $ \left(\mathbf{v, }\pi, \mathbf{t}\right) $ satisfies the estimates

and

where $ C_1 $ and $ C_2 $ only depend on the constants $ \nu $, $ \tilde{\beta}, C_{S} $.

Under the hypotheses of Theorem 3.3, stressing the dependence on $ \mathbf{u} $ of the solution to problem (2.10), we define the solution map $ \mathcal{S}:\mathbf{u\mapsto }\left(\mathbf{v}\left(\mathbf{u}\right), \pi \left(\mathbf{u}\right), \mathbf{t}\left(\mathbf{u}\right) \right). $ We will need the following result, concerning the map $ \mathbf{u\mapsto v} \left(\mathbf{u}\right) $, since only this component is involved in the cost functional.

Proposition 3.6. Under the hypotheses of Theorem 3.3, possibly with a smaller $ \eta $ depending only on $ \nu $, $ \tilde{\beta}, C_{S} $, there holds

with $ C_{L} $ depending only on $ \nu $, $ \tilde{\beta}, C_{S} $.

Proof. Write $ \mathbf{v} (\mathbf{u}_{j}) = \mathbf{z} (\mathbf{u}_{j}) +\mathbf{v}_{0} (\mathbf{u}_{j}) = \mathbf{z}_{j}+\mathbf{v}_{0j} $, with $ \mathbf{z}_{j}\in N_{B}, $ $ \mathbf{v}_{0}{}_{j}\in N_{B}^{\bot }, $ $ j = 1, 2. $ Recall that

From (3.10), since the equation for $ \mathbf{v}_{0}{}_{j} $ is linear and $ \mathbf{v}_{2}-\mathbf{v}_{1} = 0 $ on $ \Gamma _{in} $, we get

The equation for $ \mathbf{z = z}_{2}-\mathbf{z}_{1} $ reads, after routine calculations,

for all $ \boldsymbol{\varphi }\in X $. Letting $ \mathbf{\varphi = z} $, we have, using $ \Vert \mathbf{z}_{j} \Vert _{X}\leq 5\tilde{\beta} ^{-1}\eta $

Similarly, using (3.21) and $ \Vert \mathbf{v}_{0j}\Vert _{X}\leq 5\tilde{\beta}^{-1}\eta $,

and

Similar calculations give

Collecting all the above inequalities, and noting that

we deduce

Thus, if

then (3.19) follows with $ C_{L} $ depending only on $ \nu, \tilde{ \beta}, C_{S}. $

4.   The Stokes solver

The same technique used in the proof of Theorem 3.3 provides existence and uniqueness of the weak solution $ (\mathbf{v, }\pi, \mathbf{t})\in V $ of the Stokes problem

for any $ \left(\mathbf{F, }h\mathbf{, G}\right) \in V^{\ast } = X^{\ast }\times L^{2}(\Omega)\times \mathcal{T}_1 $. Moreover $ (\mathbf{v, }\pi, \mathbf{t)} $ satisfies the estimate

where $ C $ only depends on the constants $ \nu, C_{S} $. Let us introduce the Stokes solver

defined as follows: given $ \left(\mathbf{F, }h\mathbf{, G}\right) \in V^{\ast } $,

is the unique solution of (4.1). By (4.2), $ R $ is a (linear) continuous isomorphism between $ V^{\ast } $ and $ V $. Moreover, the following result holds.

Lemma 4.1. The Stokes solver $ R $ is self-adjoint.

Proof. Let $ (\mathbf{\tilde{F}, }\tilde{h}\mathbf{, \tilde{G})}\in V^{\ast } $ and $ R(\mathbf{\tilde{F}, }\tilde{h}\mathbf{, \tilde{G})} = (\mathbf{\tilde{v}}_{s}, \tilde{\pi}_{s}, \mathbf{\tilde{t}}_{s}) $. We have to show that

Indeed, we have

Substituting $ (\boldsymbol{\varphi, }\phi, \boldsymbol{\lambda)} = \left(\mathbf{ \tilde{v}}_{s}, \tilde{\pi}_{s}, \mathbf{\tilde{t}}_{s}\right) $ into (4.1) and $ (\mathbf{\varphi, }\phi, \mathbf{\lambda)} = \left(\mathbf{v} _{s}, \pi _{s}, \mathbf{t}_{s}\right) $ into (4.4) we get

and

Thus,

or, using the other equations,

which easily gives (4.3).

5.   Existence of an optimal control

We are now ready for the well-posedness analysis of the optimal control problem

Here $ J $ is given by either (2.1) or (2.2) and

is a closed, convex set in $ H_{0}^{1}(\Gamma _{c})^{d} $, where $ \eta $ denotes the constant defined in (3.5). Note that the definition of $ \mathcal{U}_{ad} $ depends on the positive constant $ c_{0} $. If we assume that $ \Vert \mathbf{v}_{in}\Vert _{\mathcal{T}_1}\leq \eta /2 $ then

and there exists a unique solution $ \mathbf{v}\left(\mathbf{u}\right) \in X $, $ \left(\pi \left(\mathbf{u}\right), \mathbf{t} \left(\mathbf{u}\right) \right) \in Q $, with $ \Vert \mathbf{v} ({\bf u})\Vert _{X}\leq C_1 \Vert \mathbf{g}(\mathbf{u})\Vert _{\mathcal{T}_1} \leq C \eta $. The following result holds.

Theorem 5.1. If $ \Vert \mathbf{v}_{in}\Vert_{\mathcal{T}_1}\leq \eta /2, $ there exists an optimal control $ \mathbf{\hat{u}}\in \mathcal{U}_{ad} $ and a corresponding optimal state $ \mathbf{\hat{v}}\in X $, $ (\hat{\pi}, \mathbf{\hat{t}})\in Q $.

Proof. Define the set of feasible points by

where $ S_{\delta } = \left\{ \mathbf{v}\in X:\left\Vert \mathbf{v}\right\Vert _{X}\leq \delta \right\} $. The proof follows from classical arguments (see, e.g., [27,Sect. 1.5]), after proving the following four properties:

$ (1) $ $ {J}(\mathbf{v, u})\geq 0 $ and $ \mathcal{F} $ is nonempty by Theorem 3.3.

$ (2) $ Let $ \left\{ \mathbf{v}_{k}, \pi _{k}, \mathbf{t}_{k} \mathbf{, u}_{k}\right\} \in \mathcal{F} $ be a minimizing sequence. Then $ \left\{ \mathbf{v}_{k}, \pi _{k}, \mathbf{t}_{k}\mathbf{, u}_{k}\right\} $ is bounded. Indeed $ \left\{ \mathbf{v}_{k}, \mathbf{u}_{k}\right\} \subset S_{\delta }\times \mathcal{U}_{ad} $ and therefore is bounded in $ X\times H_{0}^{1}(\Gamma _{1})^{d} $. From (3.4),

and hence $ \left\{ \pi _{k}, \mathbf{t}_{k}\right\} $ is bounded in $ Q $.

$ (3) $ $ \mathcal{F} $ is weakly sequentially closed. We only need to check that if $ \left\{ \mathbf{v}_{k}, \pi _{k}, \mathbf{t}_{k}\mathbf{, u} _{k}\right\} $ is a solution to

and $ \left(\mathbf{v}_{k}, \pi _{k}, \mathbf{t}_{k}\right) \rightharpoonup \left(\mathbf{v}, \pi, \mathbf{t}\right) $ in $ V $, $ \mathbf{u} _{k}\rightharpoonup \mathbf{u} $ in $ {H_{0}^{1}\left(\Gamma _{c}\right) ^{d}} $, then $ \left(\mathbf{v}, \pi, \mathbf{t, u}\right) \ $is a solution of (2.7). Since $ \mathbf{v}_{k}\rightharpoonup \mathbf{v} $ in $ X $ and the trace operator $ \mathbf{v\mapsto v}_{\left\vert \Gamma _{1}\right. } $ is continuous from $ X $ into $ \mathcal{T}_1 $, it follows that, for every $ \boldsymbol{\varphi}\in X, $ $ \phi \in L^{2}(\Omega) $ and $ \mathbf{ \lambda }\in \mathcal{T}_1^{\ast}, $

Also, $ \pi _{k}\rightharpoonup \pi $ in $ L^{2}(\Omega)\ $implies that $ b(\boldsymbol{\varphi}, \pi _{k})\rightarrow b(\boldsymbol{\varphi}, \pi). $ Moreover, the embedding of $ {H_{0}^{1}(\Gamma _{1})^{d}}\ $into $ \mathcal{T}_c $ is continuous and hence $ \mathbf{ u}_{k}\rightharpoonup \mathbf{u} $ in $ {H_{0}^{1}\left(\Gamma _{c}\right) ^{d}} $ implies that $ \langle \boldsymbol{\lambda}\boldsymbol{, }\mathbf{g}\left(\mathbf{u} _{k}\right) \rangle _{\Gamma _{1}}\rightarrow \langle \boldsymbol{\lambda} \boldsymbol{, }\mathbf{g}\left(\mathbf{u}\right) \rangle _{\Gamma _{1}} $ and $ \langle \mathbf{t}_{k}, \boldsymbol{\varphi}\rangle _{\Gamma _{1}}\rightarrow \langle \mathbf{t}, \boldsymbol{\varphi}\rangle _{\Gamma _{1}} $.

Let us write

and recall that $ X $ is compactly embedded into $ L^{4}\left(\Omega \right) ^{d} $, so that $ \mathbf{v}_{k}\rightharpoonup \mathbf{v} $ in $ X $ implies that $ \mathbf{v}_{k}\rightarrow \mathbf{v} $ in $ L^{4}\left(\Omega \right) ^{d}. $ Thus

Finally, thanks to the weak convergence of $ \left\{ \mathbf{v}_{k}\right\} $ in $ X $, we have

Passing to the limit for $ k\rightarrow \infty $ into (5.3), we deduce that $ \left(\mathbf{v}, \pi, \mathbf{t, u}\right) $ is a solution of (2.7).

$ (4) $ $ J $ is continuous and convex and therefore also sequentially weakly lower semicontinuous.

6.   Differentiability of the solution map

In order to derive first order optimality conditions, we need to show that the solution map $ \mathcal{S}:\mathbf{u\mapsto }\left(\mathbf{v}\left(\mathbf{u}\right), \pi \left(\mathbf{u}\right), \mathbf{t}\left(\mathbf{u} \right) \right) $ is Fréchet differentiable in $ \mathcal{U}_{ad}. $ Precisely, we show the following result.

Theorem 6.1. Let $ \mathbf{\bar{u}}\in \mathcal{U}_{ad}, \mathbf{k}\in H_{0}^{1}(\Gamma _{c})^{d} $ such that $ \mathbf{\bar{u}}+\mathbf{k}\in \mathcal{U}_{ad}. $ The solution map $ \mathcal{S} $ is Fréchet-differentiable at $ \mathbf{\bar{u}} $ and its derivative $ \mathcal{S}^{\prime }\left(\mathbf{\bar{u}}\right) $ is given by

where, setting $ \mathbf{\bar{v} = v}\left(\mathbf{\bar{u}}\right) $, $ \left(\mathbf{w}, q, \mathbf{s}\right) \in V $ is the solution to

Remark 6.2. Since $ \left\Vert \mathbf{\bar{u}} \right\Vert _{X} < \eta $, the bilinear form

is $ X- $coercive hence the linear problem (6.1) is well posed.

Proof. Let $ \mathbf{\bar{u}}\in \mathcal{U}_{ad} $, $ \mathbf{k}\in H_{0}^{1}(\Gamma _{c})^{d} $ such that $ \mathbf{\bar{u}}+\mathbf{k} \in \mathcal{U}_{ad} $. Denote by $ \mathcal{S}\left(\mathbf{\bar{u}}+\mathbf{k}\right) = (\mathbf{v}^{_{\mathbf{k}}}, \pi ^{_{\mathbf{k}}}, \mathbf{t}^{_{ \mathbf{k}}}) $ and $ \mathcal{S} (\mathbf{\bar{u}}) = \left(\mathbf{\bar{v}}, \bar{\pi}, \mathbf{\bar{t}}\right) $ the corresponding solutions in $ V $ to (2.7). Write

Then the difference $ (\mathbf{d, }p\mathbf{, r)} = (\mathbf{v^{_{\mathbf{k} }}-\bar{v}}, \pi ^{_{\mathbf{k}}}-\bar{\pi}, \mathbf{t^{_{\mathbf{k}}}-\bar{t}}) = \mathcal{S}\left(\mathbf{\bar{u}}+\mathbf{k}\right) -\mathcal{S} \left(\mathbf{\bar{u}}\right) $ satisfies the system

for all $ (\boldsymbol{\varphi}, \phi, \mathbf{\lambda)}\in V $. Therefore,

is a solution to the problem

for all $ (\boldsymbol{\varphi}, \phi, \mathbf{\lambda)}\in V, $ where $ \mathbf{ z = d-w} $. From estimates (3.17), (3.18) we obtain

and

where $ C $ depends only on $ \nu $, $ \tilde{\beta}, c_{0}, C_{S} $ and $ \mathbf{f } $ is the functional $ \boldsymbol{\varphi} \longmapsto -c\left(\mathbf{ d, d }, \boldsymbol{\varphi}\right) $.

Using (3.19) we obtain

from which

From (6.6) and (6.7) we have

where $ C_{3} $ depends only on the constants $ \nu $, $ \tilde{\beta}, C_{S} $.

Remark 6.3. From the previous proof, it is clear that $ \mathcal{S} $ is differentiable also with respect to the $ \mathcal{T}_c $-norm.

7.   First order optimality conditions

We are now ready to write the first order optimality conditions for our OCP (5.1). Precisely, we can show the following result.

Theorem 7.1. Let $ \mathbf{\hat{u}} $ be an optimal control and $ (\mathbf{\hat{v}}, \hat{\pi}, \mathbf{\hat{t}}) $ the corresponding optimal state. There exists a unique multiplier $ \left(\mathbf{\hat{z}, } \hat{q}, \mathbf{\hat{s}}\right) \in V $ such that $ \mathbf{\hat{u}} $ and $ \left(\mathbf{\hat{z}, }\hat{q}, \mathbf{\hat{s}}\right) $ satisfy the following optimality system:

● state problem (2.7);

● adjoint problem: for every $ (\boldsymbol{\psi, }\phi, \boldsymbol{\lambda)}\in V, $

● variational inequality: for every $ \mathbf{k}\in \mathcal{U}_{ad} $,

where $ J({\bf u}) = \tilde{J}({\bf v}({\bf u}), {\bf u}) $ is the reduced cost functional.

Proof. Define the nonlinear operator $ N:V\times \mathcal{T}_c\rightarrow V^{\ast } $ by

The nonlinear component of $ N $ is a continuous bilinear operator and therefore $ N $ is Fréchet differentiable^‡, with

^‡Since $ H_{0}^{1}\left(\Gamma _{c}\right) $ is continuously embedded into $ H_{00}^{1/2}\left(\Gamma _{c}\right) $ it follows that $ N $ is differentiable also as an operator from $ H_{0}^{1}\left(\Gamma _{c}\right) ^{d} $ into $ V^{\ast }\, $ and the two differentials coincide.

and

where

Consider now the Stokes solver $ R:V^{\ast }\rightarrow V $ defined in Section 4. In terms of $ R $ and $ N $, the state problem can be written as

with $ \mathcal{G}:V\times \mathcal{T}_1 \rightarrow V^{\ast } $.

In Lemma 7.3 below we prove that:

$ \left(a\right) $ the adjoint operator $ N_{(\mathbf{v}, \pi, \mathbf{t})}^{^{\prime }}\left(\mathbf{v}, \pi, \mathbf{t, u}\right) ^{\ast :}:V\rightarrow V^{\ast } $ is given by

where $ \mathbf{n} $ is the exterior unit normal to $ \Gamma _{c} $;

$ \left(b\right) $ the adjoint operator $ N_{\mathbf{u}}^{^{\prime }}\left(\mathbf{v}, \pi, \mathbf{t, u}\right) ^{\ast }:V\rightarrow {\mathcal{T}_c^{\ast}} $ is given by

By a well known theorem (see, e.g., [27,Thm. 1.45]) the existence of a unique multiplier $ \left(\mathbf{\hat{z}}, \hat{q}, \mathbf{\hat{s}}\right) \in V $ follows from the properties $ \left(i\right), \left(ii\right), \left(iii\right) $ listed below.

$ \left(i\right) $ $ \tilde{J} $ and $ \mathcal{G} $ are Fréchet differentiable, with

and $ \tilde{J}_{\mathbf{\pi }}^{\prime }\left(\mathbf{v, u}\right) = 0, $ $ \tilde{J}_{\mathbf{t}}^{\prime }\left(\mathbf{v, u}\right) = \mathbf{0}. $ Moreover, since $ R $ is self-adjoint,

$ \left(ii\right) $ From Theorem 6.1 we know that, for every $ \mathbf{u}\in \mathcal{U}_{ad} $, the state equation $ \left(\text{7.5}\right) $ defines a unique solution map $ \mathcal{S}\left(\mathbf{u}\right) $ which is Fréchet-differentiable.

$ \left(iii\right) $ For every $ \mathbf{\bar{u}}\in \mathcal{U}_{ad} $, the partial derivative $ \mathcal{G}_{\left(\mathbf{v, }q\mathbf{, s}\right) }^{\prime } (\mathbf{\bar{v}}, \bar{\pi}, \mathbf{\bar{t}, \bar{u}}) $, with $ (\mathbf{\bar{v}}, \bar{\pi}, \mathbf{\bar{t}}) = ({\bf v} (\mathbf{\bar{u}}), \pi(\mathbf{\bar{u}}), \mathbf{t} (\mathbf{\bar{u}})) $ is a continuous isomorphism between $ V $ and $ V^{\ast } $. Indeed, given $ (\mathbf{F, }h, \mathbf{G}) \in V^{\ast } $, the equation

in weak form reads:

Since $ \mathbf{\bar{u}}\in \mathcal{U}_{ad} $, the bilinear form

is $ X- $coercive, therefore problem (7.9) is well posed and $ \mathcal{G}_{\left(\mathbf{v, }q\mathbf{, s}\right) }^{\prime }\left(\mathbf{ \bar{v}}, \bar{\pi}, \mathbf{\bar{t}, \bar{u}}\right) $ is a continuous isomorphism between $ V $ and $ V^{\ast } $.

Introducing the Lagrangian functional

where $ \left(\mathbf{z, }q, \mathbf{s}\right) $ is a multiplier, the first order optimality conditions are obtained by differentiating $ \mathcal{L}. $ Precisely, the adjoint problem reads $ (\mathcal{\hat{L}}_{\left(\cdot \right) }^{\prime } = \mathcal{L}_{\left(\cdot \right) }^{\prime } ((\mathbf{\hat{v}}, \hat{\pi}, \mathbf{\hat{t}}) \mathbf{, \hat{u}, } \left(\mathbf{\hat{z}}, \hat{q}, \mathbf{\hat{s}}\right)) \, ) : $

or, equivalently, upon introducing the adjoint operator $ \mathcal{G}_{\left(\mathbf{v}, \pi, \mathbf{t}\right) }^{\prime } (\mathbf{\hat{v}}, \hat{\pi }, \mathbf{\hat{t}, \hat{u}}) ^{\ast }, $

for every $ \left(\boldsymbol{\psi}, \eta, \mathbf{r}\right) \in V. $

Since $ R $ is self-adjoint, we have (see the proof of Lemma 7.3)

and therefore (7.10) corresponds to the following system:

for every $ \left(\boldsymbol{\psi}, \eta, \mathbf{r}\right) \in V $, which is (7.1). The variational inequality reads, recalling (7.7),

which is (7.2).

As in Remark 6.2, since $ \mathbf{\hat{u}}\in \mathcal{U}_{ad} $, the bilinear form (adjoint of the bilinear form (6.2))

is $ X- $coercive and the linear problem (7.1) is well posed. To write the adjoint system in strong form, note that

and

since div $ \mathbf{\hat{v}} = 0, $ $ \mathbf{\hat{z}} = \mathbf{0} $ on $ \partial \Omega \backslash \Gamma _{\text{out}} $.

Moreover, in the case of the vorticity functional $ \tilde{J}_1 $,

Therefore, the adjoint system in strong form reads, setting the rotation tensor $ \nabla \mathbf{\hat{v}-(}\nabla \mathbf{\hat{v})}^{\bot } = \mathbf{\hat{W}}, $

We can also deduce that, regarding the multiplier,

In the case of the velocity tracking-type cost functional $ \tilde{J}_2 $, the strong form of the adjoint problem is straightforward to obtain, and reads as follows:

In this case, we find instead that

Remark 7.2. From (7.2) we deduce that the derivative of the reduced functional $ J\left(\mathbf{u}\right) = \tilde{J}\left(\mathbf{v}\left(\mathbf{u}\right), \mathbf{u}\right) $ at $ \mathbf{u} = \mathbf{\hat{u}}\in \mathcal{U}_{ad} $ is given by

Then, introducing the embedding operator $ B:\mathcal{T}_c^{\ast} \rightarrow H^{-1}\left(\Gamma _{c}\right) ^{d} $ defined by

we find that

where $ \mathbf{\hat{w}}\in H_{0}^{1}\left(\Gamma _{c}\right) ^{d} $ is the unique weak solution to

The variational inequality (7.2) can be written in the form

If $ B\mathbf{\hat{s}\neq 0} $ and $ \tau = 0 $, (7.19) gives

If $ B\mathbf{\hat{s}\neq 0} $ and $ \tau > 0 $, we have $ \mathbf{\hat{u} = -\hat{w}/ }\tau $ if $ \mathbf{\hat{w}/}\tau \in \mathcal{U}_{ad} $, otherwise $ \mathbf{\hat{u}} $ is given by (7.20).

Lemma 7.3. The following expressions hold:

$\left(a\right) $ $ N_{(\mathbf{v}, \pi, \mathbf{t})}^{^{\prime }}\left(\mathbf{v}, \pi, \mathbf{t, u}\right) ^{\ast}:V\rightarrow V^{\ast } $ is given by

where $ \mathbf{n} $ is the exterior unit normal to $ \Gamma _{c} $.

$\left(b\right) $ $ N_{\mathbf{u}}^{^{\prime }}\left(\mathbf{v}, \pi, \mathbf{t, u}\right) ^{\ast }:V\rightarrow \mathcal{T}_c^{\ast} $ is given by

Proof. $ (a) $ We check (7.21). Using (7.12), (7.13) and recalling (7.3), we can write

for every $ \left(\mathbf{z, }q, \mathbf{s}\right) \in V $ and every $ \left(\boldsymbol{\psi}, \eta, \mathbf{r}\right) \in V^{\ast }. $

$ (b) $ We check (7.22). We have, using (7.4),

for every $ \left(\mathbf{z, }q, \mathbf{s}\right) \in V $ and $ \mathbf{k}\in {\mathcal{T}_c} $.

8.   Numerical approximation

For the numerical solution of the OCP (5.1) we use an iterative method for the minimization of the cost functional, and we exploit the Galerkin finite element method for the numerical discretization of the state equation, the adjoint equation and the optimality equation. In this way we can take advantage of the system of first-order (necessary) optimality conditions we have derived in the previous sections. Several alternatives have been proposed in literature; see, e.g., ^{[14,25,28,43]} for the case of sequential quadratic programming methods applied to boundary control of steady Navier-Stokes flows, ^[9] for an augmented Lagrangian method, and ^[5] for inexact Newton methods. A review of numerical strategies for OCPs can be found, e.g., in ^[1,18,26].

8.1.   Numerical discretization of state, adjoint and optimality equations

The state system consists of the steady Navier-Stokes equations, while the adjoint system is given by a steady Oseen problem. Following the saddle-point approach described so far, for both these problems Dirichlet boundary conditions are enforced in weak form, by means of the Lagrange multipliers $ {\bf t} $ and $ {\bf s} $, respectively. These two latter variables are indeed related with the boundary stress over the Dirichlet boundary, see, e.g., equations (2.8) and (7.15)–(7.17). As a matter of fact, for the numerical approximation, we can rewrite the state and the adjoint equations in a conventional manner, by introducing appropriate functional spaces and suitable lifting terms to take into account non-homogeneous Dirichlet conditions. This option allows us to avoid the explicit computation of the Lagrange multipliers $ {\bf t} $ and $ {\bf s} $.

For the numerical approximation of both these problems, let us introduce two subspaces $ X_h \subset X $ and $ Q_h \subset L^2(\Omega) $, of dimension $ N_V, N_Q < + \infty $, respectively, and let us denote by $ {\bf v}_h \in X_h $ and $ p_h \in Q_h $ the FE approximations for the velocity and the pressure fields ^[44]. Similarly, $ {\bf z}_h \in X_h $ and $ q_h \in Q_h $ will denote the FE approximations for the adjoint velocity and pressure fields, respectively. We also denote by $ {\bf u}_h \in \mathcal{U}_{h} \cap \mathcal{U}_{ad} $ the FE approximation of the control variable, where $ \mathcal{U}_{h} \subset \mathcal{U} $ is a subspace of the control space $ \mathcal{U} = H_{0}^{1}\left(\Gamma_{c}\right) ^{d} $ of dimension $ N_U < + \infty $.

The Galerkin-FE approximation of the state system (2.7) thus reads as follows: we seek for $ ({\bf v}_h, p_h) \in X_h \times Q_h $ such that

where $ {\bf v}_h^D = {\bf v}_h^D({\bf u}_h) \in X_h $ is a discrete (lifting) function interpolating the Dirichlet data,

All the forms are continuous over the discrete spaces $ X_h $ and $ Q_h $. Conversely, the approximation stability is ensured by imposing that the coercivity and inf-sup conditions are still valid at the discrete level. In particular, the coercivity of $ e $ over $ X_h $ is inherited from that over $ X $; on the other hand, we require that $ b $ is inf-sup stable over $ X_h \times Q_h $, so that the following discrete Brezzi inf-sup condition, holds ^[6]:

This last property is ensured e.g., by choosing $ X_h \times Q_h $ as the space of Taylor-Hood $ \mathbb{P}_2 - \mathbb{P}_1 $ finite elements for the velocity and the pressure, respectively (in both the state and the adjoint problem); however, this choice is not restrictive --- the whole construction keeps holding for other spaces combinations as well ^[16].

We can now derive the matrix formulation corresponding to the Galerkin-FE approximation (8.1). Let us denote by $ \{\boldsymbol{\phi}_i^{{\bf v}}\}_{i = 1}^{N_v} $ and $ \{\phi_i^{p}\}_{i = 1}^{N_p} $ the Lagrangian basis of the FE spaces $ X_h $ and $ Q_h $, respectively, so that we can express the FE velocity and pressure as

We remark that the solution to (8.1) is vanishing on the whole Dirichlet boundary, so that the corresponding velocity approximation fulfilling the boundary conditions is given by $ {\bf v}_h + {\bf v}_h^D $. Thus, by denoting as $ \underline{\bf v}_h \in \mathbb{R}^{N_v} $ and $ \underline{\bf p}_h \in \mathbb{R}^{N_p} $ the vectors of the degrees of freedom appearing in (8.3), problem (8.1) can be rewritten as:

where, for $ 1 \leq i, j \leq N_v $ and $ 1 \leq k \leq N_p $,

We solve the nonlinear saddle-point problem (8.4) by means of a fixed-point method (also referred to as Picard iteration, see, e.g., [11,Sect. 8.2] it has a larger ball of convergence than Newton's method (see e.g., ^[11], Chapter 7.2 and references therein). Moreover, if a small data condition (like, e.g., (3.3) for the case at hand) is satisfied, the fixed-point method is globally convergent. Thus, starting from an initial guess $ (\underline{\bf v}^{(0)}_h, \underline{\bf p}^{(0)}_h) $, for $ k\geq 1 $ we solve

to obtain $ (\underline{\bf v}^{(k)}_h, \underline{\bf p}^{(k)}_h) $, until $ \|\underline{\bf v}^{(k)}_h - \underline{\bf v}^{(k-1)}_h\|_V \leq \varepsilon_{ \rm{tol}}^{NS} $, given a small tolerance $ \varepsilon_{ \rm{tol}}^{NS} > 0 $. As initial guess, we take the Stokes solution of (8.4). Each Oseen system (8.6) is solved by means of a sparse LU factorization.

Regarding instead the numerical approximation of the adjoint variables, the adjoint system (7.11) is linear with respect to both the adjoint velocity and pressure. Moreover, it depends on the state variables only through the term $ c\left(\mathbf{\hat{v}, \psi, \hat{z}}\right) +c\left(\mathbf{\psi, \hat{v}, \hat{z}}\right) $ and the right-hand side of (7.11); hence, the FE arrays corresponding to these two terms must be reassembled at each step during the optimization procedure. Given the approximation $ {\bf v}_h = {\bf v}_h ({\bf u}_h) $ of the (state) velocity and $ {\bf u}_h $ of the control function, the Galerkin-FE approximation of the adjoint system (7.11) thus reads as follows: we seek for $ ({\bf z}_h, q_h) \in X_h \times Q_h $ such that

where

Problem (8.7) can be rewritten algebraically as:

where, for $ 1 \leq i, j \leq N_v $ and $ 1 \leq k \leq N_p $,

8.2.   Numerical optimization

We now describe the numerical optimization method exploited for the solution of our OCP. We rely on a gradient-based algorithm exploiting the gradient $ \nabla J\left(\mathbf{{u}}_h\right) $ to iteratively update the control until a suitable convergence criterion is fulfilled; a similar approach can be found, e.g., in ^[8]. Notable instances are descent methods, such as the gradient, (non-linear) conjugate gradient, quasi-Newton or Newton methods. In the simplest case of a gradient (or steepest descent) method, starting from an initial guess $ {\bf u}^{(0)} $, we iteratively generate a sequence

where $ \lambda^{(k)} > 0 $ is a step size, until, e.g., a suitable stopping criterion is fulfilled, and $ P_{\mathcal{U}_{ad}} $ denotes the projection onto $ \mathcal{U}_{h} \cap \mathcal{U}_{ad} $, to be evaluated similarly to what we have done in (7.20).

As stopping criterion, we might require that either

or

for a given tolerance $ \varepsilon > 0 $.

At each step of the gradient method, we solve the state Navier-Stokes system (8.1), from which we obtain $ {\bf v}^{(k)} $, and the value of the cost functional $ \tilde{J}({\bf v}^{(k)}, {\bf u}^{(k)}) $. Then, we solve the adjoint system (8.7) ($ {\bf v}^{(k)} $ and $ {\bf u}^{(k)} $ being given) and obtain $ ({\bf z}^{(k)}, q^{(k)}) $. From these variable, we determine the boundary stress

and approximate (by the Galerkin finite element method) the solution of the following problem,

in order to evaluate the gradient

We highlight that at each step of the steepest descent method, we take as initial guess of the Newton method used to solve the state Navier-Stokes system the state velocity computed at the previous step of the optimization procedure, instead of computing it by solving a Stokes system, for the sake of computational efficiency.

9.   Numerical results

In this section we present the numerical results related with an OCP for the optimal control of blood flows – described through a steady Navier-Stokes model – through a bypass graft. This latter usually provides blood flow through an alternative bridging path in order to overcome critically occluded arteries; one of the most dangerous cases is related to coronary arteries, which supply the oxygen-rich blood perfusion to the heart muscle. The design of the connection between the graft vessel and the host arterial vessel is a critical factor in avoiding post-operative recurrence of the occlusion, since fluid dynamic phenomena such as recirculation, oscillating or untypical high/low shear rates, and stagnation areas, can cause the growth of another stenosis downstream the arterial-graft connection. Hence, a typical design of a bypass graft aims at minimizing some cost functionals related to haemodynamic quantities ^[10,42]. However, a rigorous model for blood circulation should take into account (ⅰ) the flow unsteadiness, (ⅱ) the arterial wall deformability, and possibly (ⅲ) complex rheological model to characterize the agammaegated nature of blood. Optimal control (and optimal design) problem require the repeated simulation of these flow equations (and the evaluation of the cost functional to be minimized); to reduce the computational burden, we adopt steady incompressible Navier-Stokes equations for laminar Newtonian flows.

We focus on the minimization of the tracking-type functional (2.2) in order to drive the blood velocity towards a specified velocity target state $ {\bf z}_d $, featuring a regular pattern; this latter is given by the solution of the Stokes system in the same domain. Other possible cost functionals for the problem at hand have been employed, e.g., in ^[37,38]; the effect of uncertainty (e.g., affecting either the residual flow across $ \Gamma_{in} $, or the geometrical configuration of the bypass) has been explored in ^[38,39,45].

Here a two-dimensional geometrical configuration is considered, although the whole theoretical analysis and the computational pipeline are perfectly valid for three-dimensional problems as well. For the case at hand, we consider a fluid flow in a tract of blood vessel where a bypass is inserted. As a matter of fact, the whole bypass graft is not simulated; its action is represented via a velocity control $ {\bf u} $ acting on the boundary $ \Gamma_c \subset \partial \Omega $, the interface where the final portion of the graft and the host vessel meet (see Figure 1). We consider two different scenarios: ($i$) a total occlusion of the host artery (setting $ {\bf v}_{in} = {\bf 0} $ on $ \Gamma_{in} $, and ($ii$) the presence of a residual blood flow ($ {\bf v}_{in} \neq {\bf 0} $) on that boundary. We choose $ \rho = 1 $, $ \bar{v} = 25 cm \, s^{-1} $, $ \bar{d} = |\Gamma_c| \approx 0.75 cm $ from which, defining the Reynolds number $ Re = {\rho \, \bar{v} \, \bar{d}} / {\mu} $, which is set equal to $ Re = 200 $, the dynamic viscosity coefficient $ \mu $ can be deduced. We first compute the solution to the Stokes system, to determine the target velocity profile $ {\bf z}_d $; we set the relaxation parameter $ \lambda^{(k)} = \tau = 0.1 $ and the regularization parameter $ \tau = 0.05 $. We consider the stopping criterion based on the difference of two successive iterates (8.10), for which we select $ \varepsilon = 0.01 $; the tolerance for the stopping criterion of the Newton method is instead given by $ \varepsilon_{ \rm{tol}}^{NS} = 10^{-5} $. For our computations we use a mesh with $ 27,233 $ triangular elements and $ 14,048 $ vertices.

We first consider the case of total occlusion of the host arterial vessel, that is, $ {\bf v}_{in} = {\bf 0} $. In Figure 2 we report the velocity (both its magnitude and the streamlines) and pressure fields corresponding to the initial control function $ {\bf u}_h^{(0)} $ and the optimal control $ \hat{\bf u}_h $ obtained when the steepest descent algorithm stops, after $ N_{it} = 157 $ iterations. In this case, the cost functional is reduced of about $ 37% $, decreasing from $ J({\bf u}_h^{(0)}) = 288.241 $ to $ J(\hat{\bf u}_h) = 181.445 $. A reduction in the magnitude of the adjoint velocity field, which can be remarked by comparing the plots displayed in Figure 3, also shows that the derivative of the cost functional evaluated at $ {\bf u}_h = \hat{\bf u}_h $ is smaller than the corresponding quantity calculated when $ {\bf u}_h = {\bf u}_h^{(0)} $.

The magnitude of the vorticity is reported in Figure 4: a global decrease of this quantity is indeed obtained also by minimizing the tracking-type cost functional, as it results by comparing the vorticity magnitude for the initial and the optimal control.

Similar results are obtained when a residual flow across $ \Gamma_{in} $ is present; in this case, we set

with $ a = 50 $ and $ b = 10^{-2} $, to model a (limited) blood flow across the occluded vessel. Its presence actually reduces the recirculation around the graft-artery junction, (compare Figure 5 and Figure 2). Also in this case a remarkable reduction of the value of the cost functional is obtained: the steepest descent algorithm stops after $ N_{it} = 151 $ iterations, and the cost functional is reduced of about $ 37% $, decreasing from $ J({\bf u}_h^{(0)}) = 289.796 $ to $ J(\hat{\bf u}_h) = 183.565 $. Similar considerations can be obtained by comparing the adjoint velocity and the vorticity fields calculated for the initial and the optimal controls (Figures 6 and 7).

Conclusions

In this paper we have proposed a saddle point approach to deal with boundary control problems for the steady Navier-Stokes equations in presence of mixed Dirichlet-Neumann boundary conditions. Both the well posedness analysis for the optimal control problem and the derivation of a system of first-order optimality conditions have been carried out, taking advantage of a suitable treatment of boundary Dirichlet controls realized by means of Lagrange multipliers. The minimization of quadratic cost functionals of the velocity has been considered for the sake of flow regularization. A steepest descent method has then been applied for numerical optimization, exploiting the Galerkin finite element method for the discretization of the derived system of optimality conditions. Numerical results have been shown for simplified two-dimensional fluid flows in a tract of blood vessel where a bypass is inserted, proving the efficacy of our approach.

Conflict of interest

The authors declare no conflict of interest.