Shape optimization for suppressing coherent structure of two-dimensional open cavity flow

This paper presents an optimal design obtained as a shape optimization problem in a domain with a singular point. For shape optimization, the eigenvalue in Snapshot Proper Orthogonal Decomposition (Snapshot POD) is defined as a cost function. The main problems are a Non-stationary Navier–Stokes problem and eigenvalue problem of Snapshot POD. An objective functional is described using Lagrange multipliers and finite element method. Twodimensional open cavity flow is adopted for an initial domain, where the domain includes a singular point. In this paper, two kinds of sensitivities assuming velocity vector in H1 and H2 are used. Using H1 gradient method for domain deformation, all triangles over a mesh are deformed as the cost function decreases. Finally, eigenvalues of Snapshot POD are compared in the initial and optimal domains.


Introduction
This paper presents an optimal design obtained using a shape optimization problem in a domain with a singular point. The author earlier reported specific examination of construction of a shape optimization method with Snapshot Proper Orthogonal Decomposition (snapshot POD; Nakazawa, 2019), in which eigenvalues of Snapshot POD are defined as a cost function. Time average and time fluctuation parts of transient flows are suppressed directly in a two-dimensional cavity flow with an isolated disk, where the isolated disk is used as a design boundary. As described herein, the shape optimization problem suggested in work by Nakazawa (2019) is improved for a domain with a boundary with a singular point, which is used as a design boundary. Thereby, a two-dimensional open cavity flow investigated by Sipp et al. (2007) and Barbagallo et al. (2009) is chosen as an initial domain. Then a numerical calculation is performed. The particular history, background, and procedure of the suggested shape optimization problem are described below.
With the rapid development of computer technology and numerical methods, shape optimization based on computational fluid dynamics (CFD) is playing an important role in fluid mechanics and aerodynamics design. Shape optimization problems in fluid dynamics were first addressed by Pironneau (1973Pironneau ( ,1974 for the respective domains in which the stationary Stokes and Navier-Stokes equations are defined. Subsequently, Haslinger and Makinen (2003), Mohammadi and Pironneau (2001), and Moubachir and Zolesio (2006) constructed fundamental frameworks of flow-field shapeoptimization problems. Recently, many researchers are examining the topic with compressible flow or Non-Newtonian flow such as Plotnikov et al. (2012), Zhang et al. (2015Zhang et al. ( ,2016, Romero (2017).
The original motivation is explained below. Although flow stabilization presents the important challenge of choosing which research field best addresses flow control, few reports of the relevant literature describe flow stabilization by shape 2. Formulation of the Problem

Main problems
For a shape optimization problem considering Snapshot POD, this paper presents the main problems: a non-stationary Navier-Stokes problem and an eigenvalue problem in Snapshot POD. By the way, the eigenvalue problem of Snapshot POD and the non-stationary Navier-Stokes problem are analyzed respectively on Stiefel manifold and Hilbert space. To derive an adjoint problem, we need choose either Stiefel manifold or Hilbert space, and this paper selects Stiefel manifold by evaluating Reynolds Average Navier-Stokes problem with eigenvalues and POD basis obtained from Snapshot POD. Below, the mapping ϕ of the position vector x from the initial domain Ω to the optimal domain ϕ(Ω) is assumed as given. Furthermore, the flow of a viscous incompressible fluid is assumed to occupy a bounded domain Ω.

Non-stationary Navier-Stokes problem
Problem 1 (Non-stationary Navier-Stokes) Find (u, p) and u 0 represents a stationary solution of the stationary Navier-Stokes problem, where the velocity vector represents u = u x , u y T . The definition of Reynolds number is the same as Sipp et al. (2007) and Barbagallo et al. (2009). Letting (w, q) be adjoint variables with respect to the velocity and the pressure, and then by discretizing in the time direction with the finite difference method, a set of necessary variables is found as ζ 1 = {u, p, w, q}, where hereinafter u = {u n } N 2 n=N 1 , p = {p n } N 2 n=N 1 , w = {w n } N 2 n=N 1 , and q = {q n } N 2 n=N 1 . The variational form of the non-stationary Navier-Stokes problem is defined as Nakazawa, Misaka and Poignard, Journal of Fluid Science and Technology, Vol.16, No.1 (2021) by setting m = N 2 − N 1 + 1 with N 1 = T 1 ∆t and N 2 = T 2 ∆t for time step size ∆t, at time t = T 1 , T 2 . The density function G n 1 (x, ζ 1 ) is presented as

Snapshot proper orthogonal decomposition
We define a Snapshot POD analysis from time t = T 1 to T 2 , where a weight function is prepared to extract arbitrary primary components from all primary components.
The correlation coefficient matrix R(N 1 , N 2 ,ũ,ũ) ∈ R m×m is composed as Let eigenvalues and eigenvectors of R be ω ∈ R m andû ∈ R m×m , where R(N 1 , N 2 ,ũ,ũ) is a positive-semidefinite matrix satisfying the eigenvalue 0 ≤ ω, and where Φ i represents the POD basis for the i-th primary component as Using the definitions, we define snapshot POD analysis as described below.
For the optimization problem, let ζ 2 = {ω,û, α, u} ∈ R d×d be the set of necessary variables used in Problem 2, where α is an adjoint variable for the eigenvectorsû.
For the shape optimization problem studied here, Snapshot POD is a main problem. It plays a role as a constraint function. Therefore, the following functional is defined as where Weight function δ j→k is intended to extract from the 1 − m th primary components to the j − k th primary components.

Reynolds Average Navier-Stokes problem
Problem 3 (Reynolds Average Navier-Stokes)Find (ū,p) : Nakazawa, Misaka and Poignard, Journal of Fluid Science and Technology, Vol.16, No.1 (2021) where andū = ū x ,ū y T andp depict the time average velocity vector and time average pressure expressed as Appendix 7.1 presents details of the tensor for the time fluctuation term C 2 . Letting (w,q) be adjoint variables with respect to the velocity and the pressure, then by discretizing in the time direction with the finite difference method, a set of necessary variables is found asζ 1 = {ū,p,w,q}. The variational form of the Reynolds Average Navier-Stokes problem is defined as The density functionḠ 1 x,ζ 1 is

Shape Optimization Problem
The shape optimization problem using Snapshot POD in the open cavity flow is constructed next, with the Lagrange function first defined to deduce the first variation. Next, based on the Kuhn-Tucker condition, the main and adjoint problems are solved to obtain the main and adjoint variables, which are substituted into the first variation to evaluate sensitivity for the shape optimization problem.

Lagrange function
We formulate the following minimization problem of the cost function f as By the way, from defintion of the correlation coefficient matrix R(N 1 , N 2 ,ũ,ũ) ∈ R m×m , the eigenvalue is meaning the L 2 norm of the velocity vectors for the time average part or the fluctuation part, and momentum energy with corresponding modes is decreasing by decreasing the cost function, as a result.
Problem 4 (Shape Optimization)After letting f (ϕ 0 , ω) be defined as Eq. (46), we find ϕ(Ω) such that By application of the Lagrange multiplier method, Lagrange function L for the shape optimization problem in this study can be expressed as The first variation of the above functional is Please see more details about deriving the first variation in Appendix 7.2. Nakazawa, Misaka and Poignard, Journal of Fluid Science and Technology, Vol.16, No.1 (2021) For sensitivity analysis, the adjoint variable method is used to derive a main problem and an adjoint problem by setting After solving the main and the adjoint problems, sensitivity is evaluated by substituting the main and adjoint variables into eq. (50).

Main and adjoint problems
Based on the adjoint variable method, the main problems of Problem 4 are introduced into Problem 3 from L(ϕ 0 , ζ(x)) w ′ ,q ′ = 0 and Problem 2 from L(ϕ 0 , ζ(x)) [α ′ ] = 0. In addition, the adjoint problems of Problem 4 are induced from Such strong forms are given as presented below.
Problem 5 (Adjoint Problem for ω)Given eigenfunctionû of Problem 2, then find α ∈ R m×m such that andû, α are the unitary matrix from Problem 5. Therefore, α is obtained as the inverse matrix or the transposed matrix ofû as Problem 6 (Adjoint Problem forû)Let the solution u of Problem 1 be given. Find α T ∈ R m×m such that In fact, solving Problem 6 is unnecessary because α has already been obtained in Problem 5.

Sensitivity of the shape optimization problem
Here, based on the H 1 gradient method (Azegami, 1986), we evaluate the sensitivity of the shape optimization problem as, with Strain tensor E and the initial domain is deformed as ϕ(Ω) = ϕ 0 (Ω) + ϵφ H 1 (Ω). In fact, the term L(ϕ 0 , ζ(x)) φ depends on the function space forū,w. as shown in Appendix 7.2. Please see details of the domain reshaping method in Appendix 7.3. As mentioned at Appendix 7.2, the domain variation φ should be in W 1,∞ and it is not guranteed that the sensitivity after subsituting all main and adjoint variable, is in W 1,∞ . Somehow it is needed to make the sensitivity to be in W 1,∞ , but, FEM is discritizing by H 1 in space. Thereby, by using the H 1 gradient method, the sensitivity is in H 1 ⊂ W 1,∞ . Many previous studies provides its effectiveness.
The stationary solution (u 0 , p 0 ) is obtained to solve the stationary Navier-Stokes problem using the Newton-Raphson method. The non-stationary solution {(u n , p n )} N n=1 is obtained to solve Problem 1 with the UMFPACK solver presented by Davis (2004) for every time step of n = 1 to N. For the material derivative term, the second-order characteristic curve method is used. Using it, Notsu and Tabata (2014) proved its mathematical proof and numerical availability. After obtaining the non-stationary solution {(u n , p n )} N n=1 , the correlation coefficient matrix R is formed for snapshot POD. The eigenvalue problem for the matrix R is solved in Problem 2 using lapack solver.
Based on the theory of the optimization problem considered herein, the adjoint problem of Problem 7 is solved to obtain (w,q) with UMFPACK solver (Davis, 2004). After evaluating the sensitivity, for domain deformation, the H 1 gradient method is used with UMFPACK solver (Davis, 2004).

Numerical Calculations and Discussion
In this section, numerical calculations and discussion are described. Before presenting numerical results, the author validates the number of nodes N nodes and the number of elements N elements in subsection 5.1. Next, based on finite element meshes (N nodes ,N elements ) decided in 5.1, the suggested shape optimization problem is performed, where sensitivities assuming velocity vector in H 1 and H 2 are called as Volume Integration Type (VIT) and Generalized J Integration Type (GJIT), and both numerical results are compared.

Spatial and temporal discretization, adaptive mesh refinement
Two-dimensional open cavity flow includes some singular points or areas in which stress concentration appears, for example two connections between the Neumann boundary condition and the Dirichlet boundary condition, and more corners of the boundary. Thereby, Adaptive Mesh Refinement (AMR) is applied for increasing numerical accuracy. Details of AMR and comparisons of numerical results with Sipp et al. (2007) and Barbagallo et al. (2009) are presented in Appendix 7.4 and Appendix 7.5. Velocity and pressure are discretized in the spatial direction using finite element method, with respective nodes and elements of (N nodes , N elements ) = (44225, 87286). For discretization in time, the time step size ∆t = 0.001 is used to take time integrations of Problem 1 at Re = 7500. Fig. 1 and the solid line on Fig. 9 show the momentum energy and eigenvalues for Snapshot POD in the initial domain with AMR for (N nodes , N elements )=(44225,87286) at Re = 7500.
Next, the author takes inner products θ i (t) 5 i=2 between the POD basis Φ i and velocity u(t) shown in Fig. 2.
From Fig. 2, long numerical simulation gives us more asymptotic flow field, but it is leading to high caluclation cost. Next, to decrease caluclation cost, a dependency of eigenvalues on sampling period (N 1 , N 2 ) = (9000, 10000) and (19000,20000) is checked in Table. 1. From Table. 1, relative errors between both of sampling period are respectively 0.446%, 0.404%, 4.385%, 5.085%. Therefore this study samples velocity vectors from N 1 = 9000 to N 2 = 10000 for Snapshot POD, for every reshaping step. As described in the next subsection, Volume Integration Type and Generalized J Integration Type are used to evaluate the sensitivity for comparing numerical results. Moreover, three cases are examined as better combinations of eigenvalues in Snapshot POD: only the time average part δ 1→1 (Case 1); the time average and fluctuation parts δ 1→m (Case 2); and only the time fluctuation part δ 2→m (Case 3). In all cases, the step size ϵ = 0.5 and the terminal condition β = 10 −5 mentioned in Section 7.3 are used and the domain is reshaped by the H 1 gradient method.

Numerical results for Volume Integration Type (VIT)
In all cases, the sum of eigenvalues of snapshot POD from i = 2 to m, 4 i=2 ω i are not decreasing, as shown in Fig. 3(a). Fig. 5 depicts sensitivities for Case 1, Case 2, and Case 3 evaluated as Volume Integration Type (VIT) in the initial domain. Nakazawa, Misaka and Poignard, Journal of Fluid Science and Technology, Vol.16, No.1 (2021) Fig. 3(b) shows the sum of eigenvalues of snapshot POD from i = 2 to m, 4 i=2 ω i for Case 1, Case 2, and Case 3. In Case 1, the value is monotonically non-increasing until the eighth step. Thereafter, it is monotonically nondecreasing. Originally, this shape optimization problem does not gurantee that the cost function is convex, because Navier-Stokes problem which is one of a representative nonlinear equation, is used as the constraint function. And more, even if GJIT is used, Case 1 is not abale to obtain the sensitivity evaluating precisely stress concentration near the point Γ d1 ∩ Γ wall , because time fluctuation part (in particular, over the 2nd primary component) involves complex flow behavior and generates stress concentration. In fact, Fig. 4 shows stream function in the initial domain for the 1st to the 4th primary components, and and stream function with the 1st primary component has smooth flow structure, and stream function with over the 2nd primary components involves complex flow structure around the corner. However, this paper does not study mathematical and physical aspects of this cost function's bahevior.

Numerical results for Generalized J Integration Type (GJIT)
For Case 2 and Case 3, the values are monotonically non-increasing until the ninth step and the tenth step, respectively. They become asymptotic at about the twentieth step. Generally, in the case in which one defines a smooth boundary as a design boundary, a cost function is monotonically increasing or decreasing but it seems that such behavior of the value is monotonically non-increasing or non-decreasing. For that reason, huge shear and normal stresses appear near the point Γ d1 ∩ Γ wall . The situation is similarly singular point.
In fact, GJIT relies on the assumption that the velocity gradient and its adjoint variable are in H 2 . Similarly, H 1 is assumed for Volume Integration Type (VIT). Therefore, GJIT is available to assume function regularity C 1 class higher than VIT. Sensitivity near point Γ d1 ∩ Γ wall is evaluated more precisely as a result, even though 4 i=2 ω i is not decreasing monotonically, rather it is monotonically non-increasing or non-decreasing. In fact, Fig. 6 depicts sensitivities for Case 1, Case 2, and Case 3 evaluated using GJIT in the initial domain. Especially, comparing Fig. 5(c) and Fig. 6(c), a sensitivity of the latter one can be evaluated precisely near the point Γ d1 ∩ Γ wall . Fig. 7 and Fig. 8 represent sensitivity for Case 2 and Case 3 evaluated using GJIT in the optimal domains. A corner at Γ d1 ∩ Γ wall in the initial domain becomes a smooth curve in the optimal domain as a result of the shape optimization problem.

Conclusions
As described herein, the author combines AMR with a shape optimization method for controlling the time fluctuation component of a transient flow, effectively based on Snapshot Proper Orthogonal Decomposition (Snapshot POD) and Reynolds Average Navier-Stokes equation (RANS). Particularly, the sum of eigenvalues in Snapshot POD is defined as the cost function. The non-stationary Navier-Stokes problem and the eigenvalue problem in Snapshot POD are used as main problems. The main problems are transformed from strong forms to weak forms with trial functions based on a standard framework of finite element method (FEM). The functional is described using the Lagrange multiplier method with FEM. Next, its material derivative is derived to evaluate sensitivity using adjoint variable method. Sensitivities of two kinds can be evaluated: VIT and GJIT. Respectively, velocity gradient and its adjoint variables are assumed in H 1 class and H 2 class. The initial domain is reshaped iteratively until the cost function satisfies the terminal condition, where Nakazawa, Misaka and Poignard, Journal of Fluid Science and Technology, Vol.16, No.1 (2021)   the H 1 gradient method is used for stable domain deformation. A two-dimensional open cavity flow is adopted for a numerical demonstration because the domain has the singular point on a corner at Γ d1 ∩ Γ wall . It clarifies the difference between VIT and GJIT. Numerical results reveal that in the cases of evaluated sensitivities by VIT, the sums of eigenvalues of snapshot POD from i = 2 to m, 4 i=2 ω i are not decreasing. Only in Case 2 and Case 3 were sensitivities evaluated by GJIT, with results showing that 4 i=2 ω i becoming asymptotic after monotonically non-increasing or non-decreasing. It is considered that such behavior of 4 i=2 ω i throughout reshaping steps results from the singular point. However, this paper is not discussing theoretically about this cost function's behavior, and it is needed to analyze that the time fluctuation part's eigenvalues are decreasing as a result of the shape optimization problem only in Case 2 and Case 3 by GJIT, in the view point of mathematical and physical aspects.

Derivation of the tensor for time fluctuation flow
Regarding the tensor for the time fluctuation term C 2 , the time fluctuation velocity u d (t) is represented as Therein, τ i (t) represents a time periodic function with unity amplitude and with frequency to each POD mode i. Finally, we are able to derive tensor C 2 with Φ ω on a POD basis Φ and eigenvalue ω as where one trigonometric identity engenders 1 2 π 0
In fact, this paper presents consideration of finite element method for numerical calculation. Also, ζ(ϕ(x)) should be in H 1 (ϕ(Ω); R d ). Therefore, Eq. (71) must be transformed as The method of deriving eq. (72) is explained by Kimura (2008). Finally, from eq. (67) and from eq. (68) and eq. (69), the functional in ϕ(Ω) can be rewritten as and in the case of ζ = ζ 1 , ζ 2 , and for H 1 (ϕ(Ω); R d ) and for H 2 (ϕ(Ω); R d ) As a result, the first variation of the functional is For sensitivity analysis, the adjoint variable method is used to derive a main problem and an adjoint problem by setting After solving the main and the adjoint problems, the sensitivity is evaluated by substituting the main and adjoint variables into eq. (76).

Domain Reshaping Method
We obtain the optimal domain numerically using the following iterating scheme, where an integer number k ∈ N denotes the iteration step, and where positive values ϵ, β ∈ R represent arbitrary small parameters, K ∈ N denotes the final step number when the iteration scheme finishes.
Step 5 Operate H 1 gradient method to obtain the smoothed sensitivity φ k H 1 .
Step 7 Judge the convergence: -If terminal condition | f k+1 − f k |/ f 0 < β is satisfied for the cost function f , then stop. -Otherwise, replace k + 1 with k and return to Step 2.
In fact, all the variables, functions, function spaces, and functionals used for this study depend on iteration step k for domain deformation. Hereinafter, for convenience, iteration step k is not described.

Adaptive Mesh Refinement
In this section, AMR distributed in Freefem++ is summarized as explained below.
Letting η(x) be a data function describing any physical state in the domain for finite element method, then a Taylor expansion of the data function η(x) with respect to any interior point x in an element over a mesh can be expressed as for t ∈ [0, 1], where x c represents the position vector at a node of an element in a mesh, and where η h denotes the linear approximation for η. The interpolation error e(x) at a displacement tδx from node x c can be expressed as As described herein, for increasing numerical accuracy on some singular points or areas where a stress concentration appears in two-dimensional open cavity flow, the Hessian matrix ∇∇η = ∇ u 0 T is used, where u 0 is the solution of the stationary Navier-Stokes problem obtained using Newton-Raphson method. Finally, the maximal interpolation error over a mesh is written as The author only works with one variable η for meeting the fixed error tolerance ϵ h that must be equidistributed over the mesh as The procedure of AMR distributed in Freefem++ is the following for i = 1 · · · N nodes , where N nodes ∈ N represents the number of the nodes.
Step 1 Set i = 1 and ϵ h arbitrarily.
Step 2 Let d i stand for the length of the edge a i .
Step 3 Compare ϵ h and d i : -If d i > ϵ h and i ≤ N nodes , then the edge a i must be cut into two edge and return to Step 2.
-If d i ≤ ϵ h and i ≤ N nodes , then replace i with i + 1 and return to Step 2.
-Otherwise stop. Additional details related to the numerical procedure are reported elsewhere in the literature (Castro-Diaz and Hecht, 2006) along with a mathematical proof in continuous and discrete spaces of the domain (Alauzet et al., 2006).

Parameter Validation
The author validates finite element meshes with and without AMR, changing the number of nodes N nodes and the number of elements N elements and performing linear stability analysis, where Fig. 10 represents numerical demonstrations of finite element meshes with and without AMR. For solving eigenvalue problems of linear stability analysis, Arnoldi method distributed in Freefem++ is applied, where 0 + 7.5i is set for the shift parameter. As reported by Barbagallo et al. (2009), 4.66 ± 7.88i is obtained at Re = 7500. Table 2 shows eigenvalues with increasing (N nodes , N elements ). Especially, it seems that the case of the initial domain with AMR approximates the result reported by Barbagallo et al. (2009) 4.66 ± 7.88i for (N nodes , N elements )=(44225,87286).
Finally, by performing Snapshot POD, eigenvalues of Snapshot POD are obtained every (N nodes , N elements ), as shown in Table 4. From them, such eigenvalues are almost converging from (N nodes , N elements )=(44225,87286). Nakazawa, Misaka and Poignard, Journal of Fluid Science and Technology, Vol.16, No.1 (2021)

Linear Stability in the optimal domain
The author investigates linear stability in the optimal domains for Case 2 and Case 3, shown as spectra in Fig. 12. Compared to results presented in Fig. 11(b), the shape optimization problem makes the disturbance more unstable than the initial domain.