Optimal shape design for a time-dependent Brinkman flow using asymptotic analysis techniques

Abstract

In this paper, we consider the geometric inverse problem of recovering an obstacle ω immersed in a bounded fluid flow Ω governed by the time-dependent Brinkman model. We reformulate the inverse problem into an optimization problem using a least squares functional. We prove the existence of an optimal solution to the optimization problem. Then, we perform the asymptotic expansion of the cost function in a simple way using a penalty method. An important advantage of this method is that it avoids the truncation method used in the literature. To reconstruct the obstacle, we propose a fast algorithm based on the topological derivative. Finally, we present some numerical experiments in two- and three-dimensional cases showing the efficiency of the proposed method.

Keywords

Optimal shape design inverse problem time-dependent Brinkman equation topological sensitivity analysis

1. Introduction

In this paper, we consider the geometric inverse problem of locating an obstacle immersed in a porous medium from internal data. The physical model is described by a time-dependent Brinkman equation. The applications of such problem cover petroleum engineering, biological flows, porous underground prospecting and medical problems.

Geometric inverse problems in fluid mechanics have been studied extensively in the literature and most of the applications considered are modeled by the Stokes equations; see for example [1,4,11,12] and the references therein. However, there are a few papers on geometric inverse problems that consider the stationary Brinkman equation as a physical model. In this context we can cite the approaches in [25]. The authors reconstruct the shape of a penetrable inclusion from boundary measurements using the factorization method. In [35] the authors used the velocity method to perform the shape derivative for the cost functional and presented a gradient descent algorithm to solve the shape inverse problem. The author of [28] solved the inverse problem by using the framework of the topological derivative.

In this work, we extend the sensitivity analysis to the case of a time-dependent Brinkman flow. More precisely, we recast the geometric inverse problem into a minimization problem using a tracking shape functional. Then, we perform the sensitivity analysis using a penalization technique combined with a prior estimate of the state solutions. A fast algorithm based on the derived topological gradient is used to solve the inverse problem numerically.

The method of the topological derivative has been intensively studied in the literature; see for instance [1,9,16,19–21,23,27–32] and the references therein. The performance and stability of the method has been well demonstrated by Ammari and his co-authors [5,6] in the context of anomaly detection.

Let us briefly introduce the method of the topological derivative: It consists in studying the local behavior of a given shape functional j with respect to a singular perturbation of the topology of the domain Ω. A typical singular perturbation of Ω is given by $Ω_{ε} = Ω ∖ ω_{ε} (x_{0})$ , where $ω_{ε} (x_{0}) = x_{0} + ε ω$ , ε is a small parameter and ω is a given shape. It leads to an asymptotic expansion of j of the form: $\begin{matrix} j (Ω_{ε}) = j (Ω) + f (ε) g (x_{0}) + o (f (ε)) . \end{matrix}$ Here $f (ε)$ is a scalar positive function that depends on the chosen perturbation and goes to zero with ε and g is called the topological gradient: $\begin{matrix} (1) & g (x_{0}) = lim_{ε \to 0} \frac{j (Ω_{ε}) - j (Ω)}{f (ε)} . \end{matrix}$ When $Ω_{ε}$ is a perturbation of Ω by a family of transformation $T_{ε} = I + ε V$ , where $V : R^{d} \to R^{d}$ is a Lipschitz vector field and $f (ε) = ε$ , the limit in (1) coincides with the definition of the shape derivative of j with respect to Ω. In view of this observation, the topological gradient can be considered as a generalization of the shape derivative. In some cases, especially when the shape functional is constrained by a an elliptic partial differential equation, the topological gradient can be computed as a singular limit of the shape derivative. When a Dirichlet boundary condition is applied to the boundary $\partial ω_{ε} (x_{0})$ , the adjoint method presented in the literature, which requires a fixed function space, cannot be directly applied to compute the topological gradient. In classical shape optimization, this requirement can be satisfied using a domain parameterization technique [15,24]. This technique involves a fixed domain and a bi-Lipschitz map between that domain and the modified domain. In the context of topology optimization, such a mapping between Ω and $Ω_{ε}$ does not exist. However, a function space independent of ε can be constructed by using a domain truncation method and the asymptotic analysis can be performed.

As mentioned by Masmoudi and his co-authors in [22], asymptotic expansions are calculated using mathematical tools that are complicated in some ways. An averaged adjoint method recently presented in [33] allows the computation of the topological gradient for semi-linear elliptic problems. As the author points out, the case of a Dirichlet boundary condition on the boundary of the inclusion cannot be treated by the proposed method and requires further investigation.

In this paper, we use a simple method based on a penalization technique and a standard prior estimate of the states to compute the topological gradient for our problem. The key idea is to remove the Dirichlet boundary condition on $\partial ω_{ε} (x_{0})$ by introducing a parameter η large enough in $ω_{ε} (x_{0})$ and vanishing in the background $Ω ∖ {\overline{ω}}_{ε} (x_{0})$ .

As far as we know, the topological sensitivity analysis for the time-dependent Brinkman equation has not been studied yet and the aim of this work is to extend the topological gradient method for the nonlinear unsteady Brinkman flow environment.

The rest of the article is as follows. In Section 2, we introduce the forward and the inverse problem. We then introduce the minimization problem and prove the existence of optimal solution. In Section 3, we compute the topological gradient of our problem using the penalization technique. In Section 4, we present some numerical results showing the efficiency of our method. In the last section we give some concluding remarks.

2. Setting of the problem

In this paper we will use the following notations: The Frobenius product of two matrices $Y, Y \in R^{d \times d}$ is defined by $\begin{matrix} X : Y = \sum_{i, j = 1}^{d} X_{i, j} Y_{i, j} . \end{matrix}$ The Frobenius norm $‖ X ‖_{F}$ of a matrix $X \in R^{n \times n}$ is given by $\begin{matrix} ‖ X ‖_{F} = \sqrt{X : X} = \sqrt{\sum_{i, j = 1}^{d} X_{i, j}^{2}} . \end{matrix}$ The product and norm of vectors $U, V \in R^{d}$ , are defined by $\begin{matrix} U \cdot V = \sum_{j = 1}^{d} U_{j} V_{j}, ‖ U ‖ = \sqrt{U \cdot U} = \sqrt{\sum_{j = 1}^{d} U_{j}^{2}} . \end{matrix}$ Let Ω be a bounded Lipschitz domain of $R^{d} (d = 2 or 3)$ containing an incompressible Newtonian fluid with kinematic viscosity $ν > 0$ , density $ρ > 0$ and having inverse permeability $α > 0$ . The boundary $\partial Ω = Γ_{N} \cup Γ_{D}$ with ${\overline{Γ}}_{N} \cap {\overline{Γ}}_{D} = \emptyset$ . Let ω be a bounded Lipschitz domain contained in Ω. For given data $g \in L^{2} (0, T; L^{2} {(Γ_{N})}^{d})$ , the Brinkman system describing the motion of the fluid in Ω in the presence of the obstacle ω is given by $\begin{matrix} (2) & \{\begin{matrix} ρ \partial_{t} u - ν Δ u + α u + \nabla p = 0 & in Ω ∖ \overline{ω} \times [0, T], \\ div u = 0 & in Ω ∖ \overline{ω} \times [0, T], \\ u = 0 & on \partial ω \times [0, T], \\ u = 0 & on Γ_{D} \times [0, T], \\ σ (u, p) n = g & on Γ_{N} \times [0, T], \\ u (\cdot, t = 0) = u^{0} & in Ω ∖ \overline{ω} . \end{matrix} \end{matrix}$ Here n is the outer unit normal vector to the boundary $\partial Ω$ or $\partial ω$ , $u^{0}$ is an initial condition and $σ (u, p)$ is the stress tensor $\begin{matrix} σ (u, p) = 2 ν \nabla^{s} u - p I, \end{matrix}$ where I is the $d \times d$ identity matrix and $\nabla^{s} u$ is the strain tensor given by: $\begin{matrix} \nabla^{s} u = \frac{1}{2} (\nabla u + {(\nabla u)}^{T}) . \end{matrix}$ Introduce the spaces $\begin{array}{l} V (Ω ∖ \overline{ω}) = {v \in H^{1} {(Ω ∖ \overline{ω})}^{d} : div v = 0, v_{∣_{Γ_{D}}} = 0}, \\ L^{2} (0, T; V (Ω ∖ \overline{ω})) = {v : (0, T) \to V (Ω ∖ \overline{ω}), measurable : \int_{0}^{T} {‖ v (\cdot, t) ‖}_{H^{1} {(Ω ∖ \overline{ω})}^{d}}^{2} d t < \infty}, \\ W (Ω ∖ \overline{ω}) = {v \in L^{2} (0, T; V (Ω ∖ \overline{ω})); \partial_{t} v \in L^{2} (0, T; L^{2} {(Ω ∖ \overline{ω})}^{d})}, \end{array}$ and $\begin{matrix} X (Ω ∖ \overline{ω}) = {v \in W; v (\cdot, 0) = u^{0}} . \end{matrix}$ Thus, the variational formulation of the Brinkman model (2) can be written in the following form: $\begin{matrix} (3) & Find u \in X (Ω ∖ \overline{ω}) such that A (u, v) = l (v), \forall v \in L^{2} (0, T, V (Ω ∖ \overline{ω})), \end{matrix}$ where $\begin{array}{l} A (u, v) = \int_{0}^{T} \int_{Ω ∖ \overline{ω}} ρ \partial_{t} u \cdot v d x d t + \int_{0}^{T} \int_{Ω ∖ \overline{ω}} ν \nabla u : \nabla v d x d t + \int_{0}^{T} \int_{Ω ∖ \overline{ω}} α u \cdot v d x d t, \\ l (v) = \int_{0}^{T} \int_{Γ_{N}} g \cdot v d s d t . \end{array}$ Using a standard Galerkin technique [17, Chap. 7], we can prove the following theorem.

Theorem 1.
For $g \in L^{2} (0, T; L^{2} {(Γ_{N})}^{d})$ the problem (3) has a unique solution $u \in W$ . Moreover, there exists a positive constant C such that: $\begin{matrix} (4) & ‖ u ‖_{L^{2} (0, T; H^{1} {(Ω ∖ \overline{ω})}^{d})}^{2} ⩽ C ‖ g ‖_{L^{2} (0, T; L^{2} {(Γ_{N})}^{d})}^{2} . \end{matrix}$

Let $O \subseteq Ω$ be a nonempty open subset. The inverse problem to be considered is the following: $\begin{matrix} (5) & Find the obstacle ω \subset Ω from a given desired velocity u_{d} \in L^{2} (0, T; H^{1} {(O)}^{d}) . \end{matrix}$ To solve the inverse problem (5) numerically, we introduce the following constrained minimization problem: $\begin{matrix} (6) & \{\begin{matrix} minimise J (ω, u) = \int_{0}^{T} \int_{O} ‖ u - u_{d} ‖^{2} d x d t + λ P (ω, Ω), \\ subject to ω \in D_{a d} and u solves the problem (2), \end{matrix} \end{matrix}$ where $\begin{matrix} (7) & D_{a d} = {ω \subset Ω, such that P (ω, Ω) < + \infty} . \end{matrix}$ Here $λ \in] 0, + \infty [$ and $P (ω, Ω)$ denotes the relative perimeter of ω in Ω. For completeness, we give the definition of the relative perimeter.
Definition 1.
The relative perimeter of $ω \subset Ω$ in Ω is defined by $\begin{array}{l} P (ω, Ω) & : = sup {\int_{ω} div Φ d x : Φ \in C_{c}^{1} (Ω; R^{d}) such that | Φ | ⩽ 1} . \end{array}$ If $P (ω, Ω) < \infty$ , we say that ω has a finite perimeter in Ω. In this case $P (ω, Ω)$ coincides with the total variation of the distributional gradient of the characteristic function of ω. $\begin{matrix} P (ω, Ω) = | D χ_{ω} | (Ω) . \end{matrix}$

The penalization term in the definition of the cost function is relevant to prove the existence of an optimal solution to the minimization problem (6), which is discussed in the following theorem. Theorem 2.
The minimization problem (6) admits at least a solution in $D_{a d}$ .
Proof.
Let ${(ω_{k})}_{k}$ be a minimizing sequence of $J$ in $D_{a d}$ : $\begin{matrix} (8) & lim_{k \to \infty} J (ω_{k}) = inf_{ω \in D_{a d}} J (ω) . \end{matrix}$ Since $ω_{k} \in D_{a d}$ , we have $P (ω_{k}, Ω) < \infty$ . Then ${(χ_{ω_{k}})}_{k}$ is bounded in the space of bounded variation $BV (Ω)$ . From the compactness embedding of $BV (Ω)$ in the space $L^{p} (Ω)$ for $1 ⩽ p < \frac{d}{d - 1}$ (see, e.g. [8]), we can extract a subsequence still denoted ${(χ_{ω_{k}})}_{k}$ such that $\begin{matrix} (9) & χ_{ω_{k}} \to χ_{ω^{}} in L^{p} (Ω) and ω^{} \in D_{a d} . \end{matrix}$ Let us consider $u_{k}$ , the solution of the problem (2) associated with $ω_{k}$ . From Theorem 1 we deduce that ${(u_{k})}_{k}$ is bounded in the space $L^{2} (0, T; H^{1} {(Ω)}^{d})$ . Then there exists $u^{} \in L^{2} (0, T; H^{1} {(Ω)}^{d})$ and a subsequence of ${(u_{k})}_{k}$ , still denoted ${(u_{k})}_{k}$ , such that $\begin{array}{l} (10) & u_{k} ⇀ u^{} in L^{2} (0, T; H^{1} {(Ω)}^{d}) as k \to \infty, \\ (11) & u_{k} \to u^{} in L^{2} (0, T; L^{2} {(Ω)}^{d}) as k \to \infty, \\ (12) & \nabla u_{k} ⇀ \nabla u^{} in L^{2} (0, T; L^{2} {(Ω)}^{d}) as k \to \infty, \\ (13) & \frac{d}{d t} u_{k} ⇀ \frac{d}{d t} u^{} in L^{2} (0, T; L^{2} {(Ω)}^{d}) as k \to \infty . \end{array}$ In the next step we show that $u^{} = u (ω^{})$ , where $u (ω^{})$ solves the problem (3) with $ω = ω^{}$ . Taking $u = u_{k}$ and $ω = ω_{k}$ in (2), we get $\begin{matrix} (14) & \int_{0}^{T} \int_{Ω} χ_{Ω ∖ \overline{ω_{k}}} (\partial_{t} u_{k} \cdot v + ν \nabla u_{k} : \nabla v + α u_{k} \cdot v) d x d t = \int_{0}^{T} \int_{Γ_{N}} g \cdot v d s d t, \end{matrix}$ for all $v \in L^{2} (0, T, V (Ω ∖ \overline{ω_{k}}))$ . Using (9)-(13) and passing to the limit in (14), we obtain $\begin{matrix} (15) & \int_{0}^{T} \int_{Ω} χ_{Ω ∖ \overline{ω^{}}} (\partial_{t} u^{} \cdot v + ν \nabla u^{} : \nabla v + α u^{} \cdot v) d x d t = \int_{0}^{T} \int_{Γ_{N}} g \cdot v d s d t, \end{matrix}$ for all $v \in L^{2} (0, T, V (Ω ∖ \overline{ω^{}}))$ . From the uniqueness of the limit we conclude from (15) that $u^{} = u (ω^{})$ . Using the lower semi-continuity of the $L^{2}$ -norm $\begin{array}{l} \int_{0}^{T} \int_{O} {‖ u^{} - u_{d} ‖}^{2} d x d t ⩽ lim_{k \to \infty} inf \int_{0}^{T} \int_{O} ‖ u_{k} - u_{d} ‖^{2} d x d t, \end{array}$ and the lower semi-continuity of the relative perimeter, $\begin{matrix} P (ω^{}, Ω) ⩽ lim_{k \to \infty} inf P (ω_{k}, Ω), \end{matrix}$ we obtain, $\begin{array}{l} J (ω^{}) & ⩽ lim_{k \to \infty} inf J (ω_{k}) ⩽ inf_{ω \in D_{a d}} J (ω) . \end{array}$ Consequently, $ω^{}$ is a minimizer of $J$ . □

3. Topological sensitivity analysis

We define the perturbed domain as $Ω_{ε} = Ω ∖ \overline{D_{ε}}$ , where $D_{ε}$ is a small perturbation of Ω of the form $D_{ε} = x_{0} + ε D$ , characterized by its center $x_{0} \in Ω$ , its size ε and a given domain D containing the origin.

In the presence of the obstacle $D_{ε}$ , the velocity $u_{ε}$ and the pressure $p_{ε}$ solve the following Brinkman system: $\begin{matrix} (16) & \{\begin{matrix} ρ \partial_{t} u_{ε} - ν Δ u_{ε} + α u_{ε} + \nabla p_{ε} = 0 & in Ω_{ε} \times [0, T], \\ div u_{ε} = 0 & in Ω_{ε} \times [0, T], \\ u_{ε} = 0 & on \partial D_{ε} \times [0, T], \\ u_{ε} = 0 & on Γ_{D} \times [0, T], \\ σ (u_{ε}, p_{ε}) n = g & on Γ_{N} \times [0, T], \\ u_{ε} (\cdot, t = 0) = u^{0} & in Ω_{ε} . \end{matrix} \end{matrix}$ To avoid the truncation method in the asymptotic expansion of the cost functional, we solve an equivalent problem in Ω by penalizing the obstacle $D_{ε}$ . $\begin{matrix} (17) & \{\begin{matrix} ρ \partial_{t} u_{ε} - ν Δ u_{ε} + α u_{ε} + η χ_{D_{ε}} u_{ε} + \nabla p_{ε} = 0 & in Ω \times [0, T], \\ div u_{ε} = 0 & in Ω \times [0, T], \\ u_{ε} = 0 & on Γ_{D} \times [0, T], \\ σ (u_{ε}, p_{ε}) n = g & on Γ_{N} \times [0, T], \\ u_{ε} (\cdot, t = 0) = u^{0} & in Ω, \end{matrix} \end{matrix}$ where $η \in R_{+}$ . If η is large enough, the solution of the problem (17) tends to the solution of the problem (16) in the $L^{2}$ -norm; see for example [7,10].

The weak formulation of the problem (17) reads: $\begin{matrix} (18) & A_{ε} (u_{ε}, v) = l_{ε} (v), \forall v \in L^{2} (0, T; V), \end{matrix}$ where $\begin{array}{l} \begin{matrix} A_{ε} (u_{ε}, v) & = \int_{0}^{T} \int_{Ω} ρ \partial_{t} u_{ε} \cdot v d x d t + \int_{0}^{T} \int_{Ω} ν \nabla u_{ε} : \nabla v d x d t \\ + \int_{0}^{T} \int_{Ω} (α + η χ_{D_{ε}}) u_{ε} \cdot v d x d t, \end{matrix} \\ l_{ε} (v) = \int_{0}^{T} \int_{Γ_{N}} g \cdot v d x d t . \end{array}$ Using Galerkin techniques as in Theorem 1, we can prove that the problem (18) has a unique solution $u_{ε} \in X (Ω)$ . The cost function is then defined in the perturbed domain as follows: $\begin{matrix} (19) & j (ε, x_{0}) = J (D_{ε}, u_{ε}) = \int_{0}^{T} J_{ε} (u_{ε} (\cdot, t)) d t, \end{matrix}$ where $\begin{matrix} (20) & J_{ε} (u_{ε} (\cdot, t)) = \int_{O} {‖ u_{ε} (\cdot, t) - u_{d} (\cdot, t) ‖}^{2} d x . \end{matrix}$ Here we neglect the penalization term, since it is not topologically differentiable.

Our goal is to compute the asymptotic expansion of the cost function $J$ to derive the topological derivative and recover the obstacle ω. To compute the asymptotic expansion of the functional $J$ , we need the following results.

Lemma 1.
Let $u_{ε}$ be the solution of the perturbed problem ( 17 ) and $u_{0}$ is the solution of the problem ( 2 ) with $ω = \emptyset$ . Then, we have the following estimate $\begin{array}{l} (21) & \int_{0}^{T} ‖ u_{ε} - u_{0} ‖_{H^{1} {(Ω)}^{d}} d t ⩽ C ε^{\frac{d}{2} + ξ}, \end{array}$ where $0 < ξ ⩽ 1$ and C is a positive constant independent of the parameter ε.
Proof.
For simplicity, we write $V$ for $V (Ω)$ . Subtracting (18) with $ε \neq 0$ from (18) with $ε = 0$ , we get $\begin{array}{l} \int_{0}^{T} \int_{Ω} ρ \partial_{t} (u_{ε} - u_{0}) \cdot v d x d t + \int_{0}^{T} \int_{Ω} ν \nabla (u_{ε} - u_{0}) : \nabla v d x d t \\ + \int_{0}^{T} \int_{Ω} (α + η χ_{D_{ε}}) (u_{ε} - u_{0}) \cdot v d x d t \\ + \int_{0}^{T} \int_{Ω} η χ_{D_{ε}} u_{0} \cdot v d x d t = 0 for all v \in L^{2} (0, T; V) . \end{array}$ If we take $v = u_{ε} - u_{0}$ as test function in the above equation and integrate in time, we get $\begin{array}{l} \frac{ρ}{2} \int_{Ω} {‖ u_{ε} (\cdot, T) - u_{0} (\cdot, T) ‖}^{2} d x + \int_{0}^{T} \int_{Ω} ν {‖ \nabla (u_{ε} - u_{0}) ‖}_{F}^{2} d x d t \\ + \int_{0}^{T} \int_{Ω} α ‖ u_{ε} - u_{0} ‖^{2} d x d t + \int_{0}^{T} \int_{D_{ε}} η ‖ u_{ε} - u_{0} ‖^{2} d x d t \\ + \int_{0}^{T} \int_{D_{ε}} η u_{0} \cdot (u_{ε} - u_{0}) d x d t = 0 . \end{array}$ Considering the Cauchy–Schwarz inequality in space and the boundedness of $u_{0}$ , we obtain $\begin{matrix} (22) & \begin{matrix} \int_{0}^{T} \int_{Ω} ν {‖ \nabla (u_{ε} - u_{0}) ‖}_{F}^{2} d x d t + \int_{0}^{T} \int_{Ω} α ‖ u_{ε} - u_{0} ‖^{2} d x d t \\ ⩽ C_{1} ε^{d / 2} \int_{0}^{T} ‖ u_{ε} - u_{0} ‖_{L^{2} {(D_{ε})}^{d}} d t . \end{matrix} \end{matrix}$ Moreover, if we use the Hölder inequality together with the Sobolev embedding theorem [14, Ch. IV, §8, Section 1.2, pp 140], we have $\begin{array}{l} \int_{0}^{T} ‖ u_{ε} - u_{0} ‖_{L^{2} {(D_{ε})}^{d}} d t & ⩽ C_{2} ε^{\frac{d}{2 q}} \int_{0}^{T} {(\int_{D_{ε}} ‖ u_{ε} - u_{0} ‖^{2 p} d x)}^{\frac{1}{2 p}} d t \\ = C_{2} ε^{\frac{d}{2 q}} \int_{0}^{T} ‖ u_{ε} - u_{0} ‖_{L^{2 p} {(D_{ε})}^{d}} d t \\ ⩽ C_{3} ε^{ξ} \int_{0}^{T} ‖ u_{ε} - u_{0} ‖_{H^{1} {(Ω)}^{d}} d t, \end{array}$ where $1 / p + 1 / q = 1$ with $q ⩾ d / 2$ and $ξ = d / 2 q$ . From (22), we obtain $\begin{array}{l} \int_{0}^{T} \int_{Ω} ν {‖ \nabla (u_{ε} - u_{0}) ‖}_{F}^{2} d x d t + \int_{0}^{T} \int_{Ω} α ‖ u_{ε} - u_{0} ‖^{2} d x d t ⩽ C_{4} ε^{d / 2 + ξ} \int_{0}^{T} ‖ u_{ε} - u_{0} ‖_{H^{1} {(Ω)}^{d}} d t . \end{array}$ This implies that $\begin{matrix} \int_{0}^{T} ‖ u_{ε} - u_{0} ‖_{H^{1} {(Ω)}^{d}} d t ⩽ C_{5} ε^{\frac{d}{2} + ξ}, \end{matrix}$ and the proof is completed. □
Proposition 1.
The functional $J$ is differentiable with respect to $u_{0}$ and $\begin{matrix} D J (\emptyset, u_{0}) φ = 2 \int_{0}^{T} \int_{O} (u_{0} - u_{d}) \cdot φ d x d t, for all φ \in H^{1} {(O)}^{d} . \end{matrix}$ Moreover $\begin{matrix} J (D_{ε}, u_{ε}) - J (\emptyset, u_{0}) = D J_{0} (\emptyset, u_{0}) (u_{ε} - u_{0}) + o (ε^{d}) . \end{matrix}$
Proof.
The variation of the cost function $J$ is given by $\begin{array}{l} J (D_{ε}, u_{ε}) - J (\emptyset, u_{0}) & = \int_{0}^{T} \int_{O} ‖ u_{ε} - u_{d} ‖^{2} d x d t - \int_{0}^{T} \int_{O} ‖ u_{0} - u_{d} ‖^{2} d x d t \\ = \int_{0}^{T} \int_{O} {‖ (u_{ε} - u_{0}) + (u_{0} - u_{d}) ‖}^{2} d x d t - \int_{0}^{T} \int_{O} ‖ u_{0} - u_{d} ‖^{2} d x d t \\ = 2 \int_{0}^{T} \int_{O} (u_{ε} - u_{0}) (u_{0} - u_{d}) d x d t + \int_{0}^{T} \int_{O} ‖ u_{ε} - u_{0} ‖^{2} d x d t . \end{array}$ Using Lemma 1, we obtain $\begin{matrix} J (D_{ε}, u_{ε}) - J (\emptyset, u_{0}) = 2 \int_{0}^{T} \int_{O} (u_{ε} - u_{0}) (u_{0} - u_{d}) d x d t + o (ε^{d}) . \end{matrix}$ Then the functional $J$ is differentiable with respect to $u_{0}$ and $\begin{matrix} D J (\emptyset, u_{0}) φ = 2 \int_{0}^{T} \int_{O} (u_{0} - u_{d}) \cdot φ d x d t . \end{matrix}$ □

Now we give the asymptotic expansion of the function j in the following theorem. Theorem 3.
The function j has the following asymptotic expansion $\begin{matrix} (23) & j (ε, x_{0}) - j (0, x_{0}) = ε^{d} η | D | δ j (x_{0}) + o (ε^{d}), \end{matrix}$ where $\begin{matrix} (24) & δ j (x_{0}) = \int_{0}^{T} u_{0} (x_{0}, t) \cdot v_{0} (x_{0}, t) d t, \end{matrix}$ and $v_{0}$ is the solution of the adjoint equation $\begin{matrix} (25) & A_{0} (v_{0}, w) = - D J (\emptyset, u_{0}) w, \forall w \in L^{2} (0, T; V) . \end{matrix}$
Proof.
Introduce the Lagrangian $\begin{matrix} (26) & L_{ε} (u, v) = j (ε) + A_{ε} (u, v) - l_{ε} (v) \forall u \in X (Ω), v \in L^{2} (0, T; V) . \end{matrix}$ Using the fact that $u_{ε}$ is the solution of (3), we have $\begin{matrix} (27) & L_{ε} (u_{ε}, v) = j (ε), \forall v \in L^{2} (0, T; V) . \end{matrix}$ Then, we get $\begin{array}{l} j (ε, x_{0}) - j (0, x_{0}) & = L_{ε} (u_{ε}, v) - L_{0} (u_{0}, v) \\ = J (D_{ε}, u_{ε}) - J (\emptyset, u_{0}) + A_{ε} (u_{ε}, v) - A_{0} (u_{0}, v) - l_{ε} (v) + l_{0} (v) . \end{array}$ Since the linear form $l_{ε}$ is independent of ε, it follows that $\begin{matrix} (28) & (l_{ε} - l_{0}) (v) = 0, \forall v \in L^{2} (0, T; V) . \end{matrix}$ Therefore $\begin{array}{l} j (ε, x_{0}) - j (0, x_{0}) & = J (D_{ε}, u_{ε}) - J (\emptyset, u_{0}) + A_{ε} (u_{ε}, v) - A_{0} (u_{0}, v) \\ = J (D_{ε}, u_{ε}) - J (\emptyset, u_{0}) + A_{0} (u_{ε} - u_{0}, v) + \int_{0}^{T} \int_{D_{ε}} η u_{ε} \cdot v d x d t . \end{array}$ It follows from Proposition 1 that $\begin{matrix} j (ε, x_{0}) - j (0, x_{0}) = \int_{0}^{T} \int_{D_{ε}} η u_{ε} \cdot v d x d t + A_{0} (u_{ε} - u_{0}, v) + D J (\emptyset, u_{0}) (u_{ε} - u_{0}) + o (ε^{d}) . \end{matrix}$ Taking $v = v_{0}$ in the above equation, where $v_{0}$ is the adjoint solution defined in (25), we obtain $\begin{matrix} (29) & j (ε, x_{0}) - j (0, x_{0}) = \int_{0}^{T} \int_{D_{ε}} η u_{ε} \cdot v_{0} d x d t + o (ε^{d}) . \end{matrix}$ Let us now examine the asymptotic behavior of the first term of right-hand side in (29). We have $\begin{matrix} (30) & \int_{0}^{T} \int_{D_{ε}} η u_{ε} \cdot v_{0} d x d t = \int_{0}^{T} \int_{D_{ε}} η u_{0} \cdot v_{0} d x d t + \int_{0}^{T} \int_{D_{ε}} η (u_{ε} - u_{0}) \cdot v_{0} d x d t . \end{matrix}$ Using change of variable $x = x_{0} + ε y$ , we deduce that $\begin{array}{l} \int_{0}^{T} \int_{D_{ε}} η u_{0} \cdot v_{0} d x d t \\ = η ε^{d} \int_{0}^{T} \int_{D} u_{0} (x_{0}, t) \cdot v_{0} (x_{0}, t) d y d t \\ + η ε^{d} \int_{0}^{T} \int_{D} (u_{0} (x_{0} + ε y, t) \cdot v_{0} (x_{0} + ε y, t) - u_{0} (x_{0}, t) \cdot v_{0} (x_{0}, t)) d y d t \\ = η | D | ε^{d} \int_{0}^{T} u_{0} (x_{0}, t) \cdot v_{0} (x_{0}, t) d t \\ + η ε^{d} \int_{0}^{T} \int_{D} (u_{0} (x_{0} + ε y, t) \cdot v_{0} (x_{0} + ε y, t) - u_{0} (x_{0}, t) \cdot v_{0} (x_{0}, t)) d y d t . \end{array}$ From the Taylor expansion and using the fact that $u_{0}$ and $v_{0}$ are regular near $x_{0}$ , we obtain $\begin{matrix} η ε^{d} \int_{0}^{T} \int_{D} (u_{0} (x_{0} + ε y, t) \cdot v_{0} (x_{0} + ε y, t) - u_{0} (x_{0}, t) \cdot v_{0} (x_{0}, t)) d y d t = o (ε^{d}) . \end{matrix}$ Using Hölder inequality for the second term of the right-hand side in (30), we get $\begin{matrix} \int_{0}^{T} \int_{D_{ε}} η (u_{ε} - u_{0}) \cdot v_{0} d x d t ⩽ η (\int_{0}^{T} ‖ u_{ε} - u_{0} ‖_{L^{2} (D_{ε})} ‖ v_{0} ‖_{L^{2} (D_{ε})} d t) . \end{matrix}$ Applying Cauchy–Schwarz inequality to the right hand-side of the above inequality, we obtain $\begin{matrix} \int_{0}^{T} \int_{D_{ε}} η (u_{ε} - u_{0}) \cdot v_{0} d x d t ⩽ η ‖ u_{ε} - u_{0} ‖_{L^{2} (0, T; L^{2} (D_{ε}))} ‖ v_{0} ‖_{L^{2} (0, T; L^{2} (D_{ε}))} . \end{matrix}$ Since $v_{0}$ is uniformly bounded on $D_{ε}$ and by change of variables, we deduce that $\begin{matrix} \int_{0}^{T} \int_{D_{ε}} η (u_{ε} - u_{0}) \cdot v_{0} d x d t ⩽ C ε^{d / 2} ‖ u_{ε} - u_{0} ‖_{L^{2} (0, T; L^{2} (D_{ε}))} . \end{matrix}$ From Lemma 1, we have $‖ u_{ε} - u_{0} ‖_{L^{2} (0, T; L^{2} (D_{ε}))} = o (ε^{d / 2})$ . We conclude that $\begin{matrix} \int_{0}^{T} \int_{D_{ε}} η (u_{ε} - u_{0}) \cdot v_{0} d x d t = o (ε^{d}), \end{matrix}$ and the proof is completed. □
Remark 1.
Let us recall the topological asymptotic expansion of the penalized problem given in Theorem 3: $\begin{matrix} (31) & j (ε, x_{0}) - j (0, x_{0}) = ε^{d} η | D | δ j (x_{0}) + o (ε^{d}) . \end{matrix}$ For the non-penalized problem, the asymptotic expansion can be determined with complicated mathematical tools and depends on the form of the perturbation. For certain forms, the asymptotic expansion can be computed explicitly and takes the following form: $\begin{matrix} (32) & j (ε, x_{0}) - j (0, x_{0}) = m (ε) δ j (x_{0}) + o (m (ε)) . \end{matrix}$ We would like to emphasize that for a large η the asymptotic expansions (31) and (32) contain the main term $δ j (x_{0})$ , but they do not have the same behavior with respect to ε. From a numerical point of view, the behavior with respect to ε is not so important as pointed out by Masmoudi and his co-authors [22]. The so-called topological gradient $δ j (x_{0})$ is used to locate the obstacles.

Note that the penalization method is a well-known technique often used in the implementation of the finite element approximation method to enforce Dirichlet boundary conditions.

Our approach to computing the topological asymptotic expansion using the penalty method is relatively simple and can easily be extended to other contexts. It has already led to very interesting numerical results for Stokes and Brinkman flows in the stationary case [22,28].

4. Implementation details and numerical results

For the computation of the state and the adjoint state, we use the projection method which was presented in [13,34]. The method is based on the idea of splitting the PDE into a series of simpler equations. The first step is to neglect the incompressibility condition and compute a tentative velocity. The next step is to correct the velocity by performing a projection onto the divergence-free vector fields. The tentative velocity $u^{*}$ at the time step k is computed by solving the equation $\begin{matrix} (33) & \frac{1}{Δ t} (u^{*} - u^{k}) - ν Δ u^{*} + α u^{*} + ζ \nabla p^{k} = 0, \end{matrix}$ where $Δ t$ is the time step size and ζ is an adjustable factor for flexible control of the pressure information. Note here that actually the weak form of (33) is solved, but for readability, we use the strong form to explain the projection method. When ζ is set to zero, we obtain the non-incremental pressure correction scheme, also known as Chorin’s method [13]. When $ζ = 1$ , we obtain the incremental pressure correction scheme, also known as IPCS [18]. The boundary conditions $u^{*}$ are the same as for the original problem. Since $u^{*}$ does not satisfy the incompressibility equation (5.2), we need to correct it. We do this by evaluating the pressure and the velocity in certain terms at the time step $k + 1$ : $\begin{matrix} (34) & \frac{1}{Δ t} (u^{k + 1} - u^{k}) - ν Δ u^{*} + α u^{k + 1} + ζ \nabla p^{k + 1} = 0, \end{matrix}$ Subtracting (33) from (34), we end up with $\begin{matrix} (35) & \frac{1}{Δ t} (u^{k + 1} - u^{*}) + α (u^{k + 1} - u^{*}) + \nabla (p^{k + 1} - ζ p^{k}) = 0 . \end{matrix}$ Now we take the divergence of the equation (35) to get rid of the unknown velocity, using that $\nabla \cdot u^{k + 1} = 0$ , and get a Poisson equation for the pressure $p^{k + 1}$ $\begin{matrix} (36) & Δ (p^{k + 1} - ζ p^{k}) = (\frac{1}{Δ t} + α) \nabla \cdot u^{*} . \end{matrix}$ Solving for the pressure $p^{k + 1}$ , it remains only one unknown in equation (35), and the velocity $u^{k + 1}$ can be computed. The updated velocity is divergence-free, but not exactly satisfy the boundary conditions. It is possible to enforce boundary conditions on the updated velocity.

The pressure correction scheme in weak form can be summarized in three equations to be solved: $\begin{array}{l} (a) \int_{Ω} (\frac{1}{Δ t} + α) u^{*} \cdot v + ν \nabla u^{*} : \nabla v d x = \int_{Ω} \frac{1}{Δ t} u^{k} \cdot v - ζ \nabla p^{k} \cdot v d x, \\ (b) \int_{Ω} \nabla p^{k + 1} \cdot \nabla q d x = \int_{Ω} ζ \nabla p^{k} \cdot \nabla q - (\frac{1}{Δ t} + α) \nabla \cdot u^{*} \cdot q d x, \\ (c) \int_{Ω} (\frac{1}{Δ t} + α) u^{k + 1} \cdot v d x = \int_{Ω} (\frac{1}{Δ t} + α) u^{*} \cdot v - \nabla (p^{k + 1} - ζ p^{k}) d x . \end{array}$ The three steps can be described as $(a)$ tentative velocity, $(b)$ pressure correction and $(c)$ velocity update. It is now clear that the matrices on the left-hand side are independent of the solution at time k and only need to be assembled in the first iteration.

The numerical solutions are obtained using the FEniCS environment [26]. FEniCS is a collection of open-source software components that enable the automatic solution of differential equations. The differential equations can be specified in almost mathematical notation, such as finite element variational problems. The language used to express the weak forms is called UFL (Unified Form Language) [3]. FEniCS can be programmed in both C++ and Python, using a library called DOLFIN, which is the main user interface of FEniCS. A just-in-time compiler called FFC (FEniCS Form Compiler) [2] generates efficient low-level C++ code from the high-level mathematical description (UFL) of the variational problem.

DOLFIN abstracts the spatial discretization problem, but not the temporal discretization problem. Transient DOLFIN-based solvers typically consist of a hand-written temporal loop in which one or more discrete variational (finite element) problems are solved with DOLFIN in each iteration.

To solve the optimization problem, we apply a fast non-iterative algorithm based on the following steps.

Algorithm 1.

Solve the direct (17) and adjoint (25) problem in the domain Ω for an obstacle $D_{ε} = \emptyset$ .

Compute the topological gradient $δ j (x_{0})$ given by (24) for all points $x_{0} \in Ω$ .

Determine the location $x_{0}$ where the topological sensitivity is most negative.

In the following numerical results, the parameters of the problem are chosen as follows: $\begin{matrix} T = 1, Δ t = 0.01, ρ = 1, ν = 0.1 and α = 1 . \end{matrix}$

Fig. 1.
Image of the topological gradient without noise: case of a single obstacle with different positions.
4.1. Numerical simulations in the two dimensional case

In these examples, the computational domain Ω is the unit square $[0, 1] \times [0, 1]$ and $O = Ω$ . The topological gradient is computed on a square grid of $16 \times 16$ points. To generate the data $u_{d}$ , we use the boundary conditions $\begin{array}{l} u = (- sin (π x_{1}), 0) on {1} \times [0, 1], \\ u = (0, 0) on [0, 1] \times {1} \cup [0, 1] \times {0} and \\ u = 0 on \partial ω, \\ p = (0, 0) on {0} \times [0, 1] . \end{array}$ Figure 1 shows the image of the topological gradient for different positions of the obstacle without noise. The red circle is the boundary of the exact obstacle to be reconstructed. The most negative values of the topological gradient are indicated in blue. Figure 2 shows the image of the topological gradient in the case of multiple obstacles, without noise. Figure 3 shows the image of the topological gradient with respect to different relative noise levels. In these two-dimensional numerical results, the quality of the detection of the obstacles with the topological gradient is satisfactory and is stable with respect to a small amount of noise.

Fig. 2.

Image of the topological gradient without noise: case of multiple obstacles.

Fig. 3.

On the left the image of the topological gradient with relative noise level $δ = 10 %$ and on the right the image of the topological gradient with relative noise level $δ = 30 %$ .

4.2. Numerical simulations in the three dimensional case

4.2.1. Reconstruction with real synthetic internal measurement

In these examples, the computational domain Ω is the unit cube ${[0, 1]}^{3}$ , and $O = Ω$ . The topological gradient is computed on the cube grid of $16 \times 16 \times 16$ points. To generate the data $u_{d}$ we use the boundary conditions $\begin{array}{l} u = (- sin (π x_{1}), 0, 0) on {1} \times [0, 1] \times [0, 1], \\ u = (0, 0, 0) on [0, 1] \times [0, 1] \times {1} \cup [0, 1] \times [0, 1] \times {0}, \\ u = (0, 0, 0) on \partial ω, \end{array}$ and $\begin{matrix} p = (0, 0, 0) on {0} \times [0, 1] \times [0, 1] . \end{matrix}$ In the case of single obstacle, the exact shape ω to be recovered is the ball of radius $r = 0.1$ , centered at the point $(0.5, 0.5, 0.5)$ . In the case of multiple obstacle, the exact shape to be recovered is composed of two balls with radius $r = 0.1$ , centered at the points $(0.5, 0.3, 0.2)$ and $(0.5, 0.3, 0.8)$ .

Fig. 4.

On the left the image of the topological gradient on the cube ${[0, 1]}^{3}$ and on the right the image of the topological gradient on the half cube $[0, 0.5] \times [0, 1] \times [0, 1]$ , without noise.

Fig. 5.

On the left the image of the topological gradient on the half cube $[0, 1] \times [0, 0.5] \times [0, 1]$ and on the right the image of the topological gradient on the half cube $[0, 1] \times [0, 1] \times [0, 0.5]$ , without noise.

Fig. 6.

Figures 4–8 show the image of the topological gradient without noise. The most negative values of the topological gradient are indicated in blue. Figures 9–12 show the image of the topological gradient in the case of single and multiple obstacles, with relative noise levels. The quality of the detection of the obstacles with the topological gradient is satisfactory and is stable with respect to the amount of noise levels.

4.2.2. Reconstruction with penalized synthetic internal measurement

In the following example, we use a penalization synthetic measurement $u_{d}$ with $η = 10^{8}$ and $O = Ω$ . The exact obstacle ω to be recovered is the ball of radius $r = 0.1$ centered at the point $0.5, 0.5, 0.5$ . Figures 13–14, show the image of the topological derivative, without noise. In this case, the reconstruction is satisfactory.

Fig. 7.

On the left the image of the topological gradient for the first obstacle on the subdomain $[0, 1] \times [0, 0.3] \times [0, 1]$ and on the right the image of the topological gradient for the second obstacle on the subdomain $[0, 1] \times [0, 0.7] \times [0, 1]$ , without noise: case of multiple obstacles.

Fig. 8.

On the left the image of topological gradient for the first obstacle on the subdomain $[0, 1] \times [0, 1] \times [0, 0.2]$ and on the right the image of the topological gradient for the second obstacle on the subdomain $[0, 1] \times [0, 1] \times [0, 0.8]$ , without noise: case of multiple obstacles.

Fig. 9.

On the left the image of the topological gradient on the cube ${[0, 1]}^{3}$ and on the right the image the topological gradient on the half cube $[0, 0.5] \times [0, 1] \times [0, 1]$ , with relative noise level $δ = 10 %$ .

Fig. 10.

Fig. 11.

On the left the image of the topological gradient on the cube $[0, 1] \times [0, 1] \times [0, 1]$ and on the right the image of the topological gradient on the half cube $[0, 0.5] \times [0, 1] \times [0, 1]$ , with relative noise level $δ = 30 %$ .

Fig. 12.

On the left the image of the topological gradient on the half cube $[0, 1] \times [0, 0.5] \times [0, 1]$ , and on the right the image of the topological gradient on the half cube $[0, 1] \times [0, 1] \times [0, 0.5]$ , with relative noise level $δ = 30 %$ .

Fig. 13.

Fig. 14.

Fig. 15.

On the left the internal data $u_{d}$ at the final time T on the half cube $[0, 0.5] \times [0, 1] \times [0, 1]$ , and on the right the momentum $x \mapsto \int_{0}^{T} ‖ u_{d} (x, t) ‖^{2} d t$ on the half cube $[0, 0.5] \times [0, 1] \times [0, 1]$ , without noise.

Remark 2.

We would like to emphasize that if the measurement $u_{d}$ is available in the whole domain Ω, the inverse problem can be solved by simply displaying the momentum $x \mapsto \int_{0}^{T} ‖ u_{d} (x, t) ‖^{2} d t$ ; see Fig. 15. This direct approach gives good results but cannot be used if the measurement is available in some part of the domain Ω. As shown in Section 4.3, our approach is relatively simple and can be applied to solve the inverse problem from a partial internal measurement.

4.3. Reconstruction using partial internal measurement

In the following example, the measurements $u_{d}$ are only available on the subdomain $O$ given by: $\begin{matrix} O : = Ω \cap {(x, y, z) \in R^{3} : {(x - 0.5)}^{2} + {(y - 0.5)}^{2} + {(z - 0.5)}^{2} ⩾ 0 . 2^{2}} . \end{matrix}$ Figures 16–17 show the image of the topological gradient in the case of partial internal measurement without noise. We observe that the topological gradient gives a satisfactory approximation of the shape of the obstacle.

Fig. 16.

On the left the image of the topological gradient on the cube ${[0, 1]}^{3}$ and on the right image of the topological gradient on the half cube $[0, 0.5] \times [0, 1] \times [0, 1]$ , without noise.

Fig. 17.

5. Conclusion

In this work, we considered a shape detection problem from internal data. The model problem is governed by a time-dependent Brinkman equation. We have reformulated our detection problem into an optimization problem. The topological sensitivity analysis was employed to evaluate the topological gradient of the proposed shape functional. To avoid the truncation method, we have used a simple way based on a penalization technique and a standard prior estimate of the states. A one-shot algorithm based on the topological gradient was used to provide useful information for the positioning of the obstacles.

In a forthcoming work, we intend to combine the topological gradient results with the level-set method to obtain a good approximation of shape problem.

References

A.B.

Abda,

Hassine,

Jaoua and

Masmoudi, Topological sensitivity analysis for the location of small cavities in Stokes flow, SIAM journal on control and optimization 48(5) (2010), 2871–2900. doi:10.1137/070704332.

M.S.

Alnæs,

Logg,

K.-A.

Mardal,

Skavhaug and

H.P.

Langtangen, Unified framework for finite element assembly, International Journal of Computational Science and Engineering 4(4) (2009), 231–244. doi:10.1504/IJCSE.2009.029160.

M.S.

Alnæs,

Logg,

K.B.

Ølgaard,

M.E.

Rognes and

G.N.

Wells, Unified form language: A domain-specific language for weak formulations of partial differential equations, ACM Transactions on Mathematical Software 40(2) (2014).

C.J.

Alves and

A.L.

Silvestre, On the determination of point-forces on a Stokes system, Mathematics and computers in Simulation 66(4–5) (2004), 385–397. doi:10.1016/j.matcom.2004.02.007.

Ammari,

Bretin,

Garnier,

Jing,

Kang and

Wahab, Localization, stability, and resolution of topological derivative based imaging functionals in elasticity, SIAM Journal on Imaging Sciences 6(4) (2013), 2174–2212. doi:10.1137/120899303.

Ammari,

Garnier,

Jugnon and

Kang, Stability and resolution analysis for a topological derivative based imaging functional, SIAM Journal on Control and Optimization 50(1) (2012), 48–76. doi:10.1137/100812501.

Angot,

C.-H.

Bruneau and

Fabrie, A penalization method to take into account obstacles in incompressible viscous flows, Numerische Mathematik 81(4) (1999), 497–520. doi:10.1007/s002110050401.

Bergounioux, Introduction au traitement mathématique des images – méthodes déterministes, Vol. 76, Springer, 2015.

Borrvall and

Petersson, Topology optimization of fluids in Stokes flow, International journal for numerical methods in fluids 41(1) (2003), 77–107. doi:10.1002/fld.426.

10.

Carbou,

Fabrie et al., Boundary layer for a penalization method for viscous incompressible flow, Advances in Differential equations 8(12) (2003), 1453–1480.

11.

Caubet,

Conca and

Godoy, On the detection of several obstacles in 2d Stokes flow: Topological sensitivity and combination with shape derivatives, Inverse Problems & Imaging 10(2) (2016), 327. doi:10.3934/ipi.2016003.

12.

Caubet and

Dambrine, Localization of small obstacles in Stokes flow, Inverse Problems 28(10) (2012), 105007. doi:10.1088/0266-5611/28/10/105007.

13.

A.J.

Chorin, Numerical solution of the Navier–Stokes equations, Mathematics of computation 22(104) (1968), 745–762. doi:10.1090/S0025-5718-1968-0242392-2.

14.

Dautray and

Lions, Functional and Variational Methods, Mathematical Analysis and Numerical Methods for Science and Technology, Vol. 2, 1988.

15.

M.C.

Delfour and

J.-P.

Zolésio, Shapes and Geometries, 2nd edn, Advances in Design and Control, Vol. 22, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2011, Metrics, analysis, differential calculus, and optimization.

16.

H.A.

Eschenauer,

V.V.

Kobelev and

Schumacher, Bubble method for topology and shape optimization of structures, Structural optimization 8(1) (1994), 42–51. doi:10.1007/BF01742933.

17.

Evans, Partial Differential Equations, American Mathematical Society. American Mathematical Society, 2010.

18.

J.-L.

Guermond,

Minev and

Shen, An overview of projection methods for incompressible flows, Computer methods in applied mechanics and engineering 195(44–47) (2006), 6011–6045. doi:10.1016/j.cma.2005.10.010.

19.

J.K.

Guest and

J.H.

Prévost, Topology optimization of creeping fluid flows using a Darcy–Stokes finite element, International Journal for Numerical Methods in Engineering 66(3) (2006), 461–484. doi:10.1002/nme.1560.

20.

Guillaume and

K.S.

Idris, The topological asymptotic expansion for the Dirichlet problem, SIAM journal on Control and Optimization 41(4) (2002), 1042–1072. doi:10.1137/S0363012901384193.

21.

Guillaume and

K.S.

Idris, Topological sensitivity and shape optimization for the Stokes equations, SIAM Journal on Control and Optimization 43(1) (2004), 1–31. doi:10.1137/S0363012902411210.

22.

Hassine,

Jan and

Masmoudi, From differential calculus to 0–1 topological optimization, SIAM journal on control and optimization 45(6) (2007), 1965–1987. doi:10.1137/050631720.

23.

Hassine and

Masmoudi, The topological asymptotic expansion for the quasi-Stokes problem, in: ESAIM: Control, Optimisation and Calculus of Variations, Vol. 10, 2004, pp. 478–504.

24.

Henrot and

Pierre, Variation et optimisation de formes, Mathématiques & Applications (Berlin) [Mathematics & Applications], Vol. 48, Springer, Berlin, 2005, Une analyse géométrique. [A geometric analysis].

25.

Lechleiter and

Rienmüller, Factorization method for the inverse Stokes problem, Inverse Problems & Imaging 7(4) (2013), 1271. doi:10.3934/ipi.2013.7.1271.

26.

Logg,

K.-A.

Mardal and

Wells, Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, Vol. 84, Springer Science & Business Media, 2012.

27.

Mastoid,

Pommier and

Samet, The topological asymptotic expansion for the Maxwell equations and some applications, Inverse Problems 21(2) (2005), 547. doi:10.1088/0266-5611/21/2/008.

28.

Meftahi, Optimal shape design in three-dimensional Brinkman flow using asymptotic analysis techniques, Quarterly of Applied Mathematics 75(3) (2017), 525–537. doi:10.1090/qam/1464.

29.

Pommier and

Samet, The topological asymptotic for the Helmholtz equation with Dirichlet condition on the boundary of an arbitrarily shaped hole, SIAM Journal on Control and Optimization 43(3) (2004), 899–921. doi:10.1137/S036301290241616X.

30.

Samet,

Amstutz and

Masmoudi, The topological asymptotic for the Helmholtz equation, SIAM Journal on Control and Optimization 42(5) (2003), 1523–1544. doi:10.1137/S0363012902406801.

31.

Schumacher, Topologieoptimierung von Bauteilstrukturen unter Verwendung von Lochpositionierungskriterien, PhD thesis, Forschungszentrum für Multidisziplinäre Analysen und Angewandte, 1996.

32.

Sokolowski and

Zochowski, Topological derivative in shape optimization, in: Encyclopedia of Optimization, p. 3908–3918, 2009. doi:10.1007/978-0-387-74759-0_682.

33.

Sturm, Topological sensitivities via a Lagrangian approach for semilinear problems, Nonlinearity 33(9) (2020), 4310. doi:10.1088/1361-6544/ab86cb.

34.

Temam, Sur l’approximation de la solution des équations de Navier–Stokes par la méthode des pas fractionnaires (ii), Archive for rational mechanics and analysis 33(5) (1969), 377–385. doi:10.1007/BF00247696.

35.

Yan,

Liu and

Jing, Shape inverse problem for Stokes–Brinkman equations, Applied Mathematics Letters 88 (2019), 222–229. doi:10.1016/j.aml.2018.09.003.