Shape optimization problem for transient Non-Newtonian fluid in hybridized discontinuous Galerkin method

This paper presents a shape optimazation method for transient Non-Newtonian fluid which is playing important roles of calculating blood flow, oil flow and so on. So far, the author constructed a shape optimization problem for suppressing transient Newtonian fluid by using Snapshot POD, and extends it toward to Non-Newtonian fluid, here. For such the suggested shape optimization, the eigenvalue in Snapshot POD is defined as a cost function, where the constraint functions are the Oldroyd-B model, and an eigenvalue equation of Snapshot POD. For numerical calculations, a two–dimensional cavity flow with a disk-shaped isolated body is adopted for an initial domain. To descritize the Oldroyd-B model spatially, Galerkin Method (GM) and Hybridized Discontinuous Galerkin Method (HDGM) are used to compare numerical accuracies. As a result, it is considered that HDGM is able to obtain better solutions than GM during numerical validations. Finally, eigenvalues of Snapshot POD are compared in the initial and optimal domains obtained by HDGM.


Introduction
Non-Newtonian fluids play important roles in widely diverse engineering and industry phenomena such as blood flow in a vein inside a living body and oil flow in a pipe inside the engine of a vehicle. Such non-Newtonian fluids have properties of fluids and solids. Their very complex behavior in spatial and temporal spaces might appear even at a low Reynolds number. Below is an overview presenting aspects of Non-Newtonian fluids.
The Carreau model, which is studied mainly as a representative blood flow model, shows Newtonian fluid behavior in a low shear rate as a power-law fluid at an intermediate shear rate, and as a Newtonian fluid again at a high shear rate. Such flow behavior is interesting for researchers who are studying fluid dynamics from the perspective of linear stability. Giannetti et al. (2011), Haque et al. (2012) and Citro et al. (2015) and many engineers have examined linear stability for a Carreau model. Zhang et al. (2013) addressed the Oldroyd-B and FENE-P models used as polymer viscoelastic fluid flow models. In fact, in engineering and industry especially, optimal designs are expected to be necessary to minimize arbitrary cost functions to facilitate the practical use of non-Newtonian fluids. For mathematical aspects of shape optimization problems related to Non-Newtonian fluids, Sokolowski and Stebel (2011) mathematically analyzed shape sensitivity. Zhang et al. (2015Zhang et al. ( ,2016 and Romero et al. (2017) solved topology optimization problems for Carreau models and Power law models. Nevertheless, few numerical simulations for shape optimization problems have been undertaken to examine time-dependent non-Newtonian fluids. In a nod to the challenge above, this paper presents a solution for shape optimization problems for the transient Oldroyd-B model.
The author constructed a shape optimization method using Snapshot POD (Nakazawa, 2019(Nakazawa, ,2021, for which the eigenvalue in a Snapshot POD is defined as the cost function. Such a Snapshot POD has the same numerical procedure as primary component analysis (PCA), and plays a role as a Reduced Order Model (ROM)-based method as well as Dynamic Mode Decomposition (DMD). Using it, the number of modes of Snapshot POD (if contribution rates are higher than arbitrary large numbers as 0.99) is generally less than the number sampled from each time step. From the perspective of fluid dynamics, the time-dependent flow can be decomposed into a time average part (the first primary component) and a fluctuation part (over the second primary components). The cost function can be defined to distinguish both parts. Moreover, under the definition of Snapshot POD given in this paper, the eigenvalue represents the L 2 norm of the velocity vectors for the time average part or the fluctuation part. The momentum energy with corresponding modes is decreasing by decreasing the cost function, as a result. A brief summary of the shape optimization problem is presented below.
The sum of eigenvalues in Snapshot POD is defined as the cost function. For this study, the non-stationary Navier-Stokes problem and the eigenvalue problem in Snapshot POD are used as the main problems. The main problems are transformed from strong forms to weak forms with trial functions based on a standard application of finite element method (FEM). The functional is described using Lagrange multipliers with FEM. Next, its first variation, which is the same as the material derivative, is derived to evaluate sensitivity using adjoint variable method. At that time, it is assumed that the velocity vector is in H 2 for avoiding singularity of sensitivity to the greatest extent possible based on the Generalized J integral described by Ohtsuka (1981Ohtsuka ( , 1985Ohtsuka ( , 1986, by Khludnev (2000, 2018) and by Kimura (2008). After evaluating the sensitivity, an initial domain is reshaped iteratively to obtain an optimal domain. Then the H 1 gradient method (Azegami, 1996) is used for stable domain deformation.
As described above, this paper introduces the Oldroyd-B model to the author's earlier research and extends it to non-Newtonian fluids, where the constitutive equation of Oldroyd-B model is added to a constraint function, and for which the two-dimensional open cavity flow with an disk-shaped isolated body is adopted, as described by Nakazawa (2019). By the way, this paper is focusing on how to control effectively time fluctuation part of velocity. Therefore, the sensitivity is based on the only Reynolds average Navier-Stokes equation obtained by subsituting velocity and pressure after calculating the nonstationary Navier-Stokes equation and the constituteive equation of the Oldroyd-B model. It is meaning that the adjoint problem of the Orldoyd-B model is not solved here. The reason to choice the Oldroyd-B model is that Lukácová-Medvidiová et al. (2017) derived mathematically an approximate expression with the characteristic curve method for the material derivative term in the Peterlin model, and more the Oldroyd-B model is involving the material derivative term and simpler than Peterlin model. FreeFEM++ (Hecht, 2012) is used for all numerical calculations. Discretizing schemes of two kinds are used to compare calculation results: Galerkin Method (GM) and Discontinuous Galerkin Method (DGM). DGM was constructed for the stable calculation of a compressible fluid, a non-Newtonian fluid, or any another complex flow by Dolejši and Feistauer (2015), Cockburn et al. (2000), and Pietro and Ern (2012). In recent years, Oikawa (2013) reported Hybridized Discontinuous Galerkin Method (HDGM) for convection-diffusion problems with mixed boundary conditions. A feature is that it can greatly reduce the number of globally coupled degrees of freedom compared with classical DGM. Higher coercivity of a convective part is achieved by adding an upwinding term. Furthermore, it was shown that the approximate solution given by the scheme is close to the solution of a purely convective problem when the viscosity coefficient is small. By solving the shape optimization problems in DG and HDGM, the eigenvalues of Snapshot POD can be compared in both optimal domains.

Initial domain and domain variation
Let Ω 0 be a fixed bounded Lipschitz domain in R d (d ∈ N). Additionally, let Ω be an open subset of Ω 0 , where a position vector is denoted as x ∈ R d for spatial dimension d = 2. The boundary is represented as ∂Ω = Γ D ∪ Γ N for the Dirichlet and the Neumann boundary Γ D ∪ Γ N . For mesh deformation, Γ Design ⊂ Γ D is used as the design boundary as the cost function is decreasing.
For a shape optimization problem, the mapping ϕ of the position vector x from the initial domain Ω to the optimal domain ϕ(Ω) = ϕ 0 (Ω) + ϵφ(Ω) with a small parameter ϵ < 1 is assumed, where ϕ 0 (Ω) and φ(Ω) depict an identity map and a domain variation. Letting ζ be a scholar value function, an objective functional can be represented as where G(x, ζ(x)) is a density function. By taking the material derivative of the objective functional, the first variation is defined as where (·) ′ represents the Fréchet derivative with respect to ζ. Using adjoint variable method, we take L(ϕ 0 , ζ(x)) [ ζ ′ ] = 0 to solve the main and adjoint problems. Next, by substituting the main and the adjoint variables into L(ϕ 0 , ζ(x)) [ φ ] , the sensitivity is evaluated.

Governing equation
The nondimensional velocity u, pressure p, and tensor B for the constitutive equation of Oldroyd-B model are assumed to be satisfied in this domain Ω. A Reynolds number Re is defined with reference length and the reference speed. The Weissenberg number is written as Wi and α = 1 Wi (B − I) is prepraed for convenience sake, ϵ 1 represets an artificial viscosity parameter, I represents the identity matrix. The initial boundary value problem for Oldroyd-B model discussed in this paper is defined as shown below.
Therein, t ∈ (0, T ) and ν respectively represent the time and the outward variable on the boundary, and u D denotes a boundary condition on the Dirichlet boundary Γ D , and u 0 and B 0 respectively represent the stationary solutions of the stationary Stokes equation and the constitutive equation of Oldroyd-B model. Obtaining velocity vectors for every time step after solving Oldroyd-B model, Snapshot POD is performed to obtain eigenvalue ω and POD basis Φ for a sampling number m, Finally, we prepare the Reynolds Average Navier-Stokes equation using POD basis as shown below. where and the time average velocityū and the time fluctuation velocityũ are represented as p is time-averaged pressure for the time-averaged velocityū.

Shape optimization problem
Concrete techniques of the shape optimization problem constructed by Nakazawa et al. (2021) are used here. Particularly after solving the governing equation (3)-(10) and sampling velocity vectors from t = T 1 to T 2 , an eigenvalue problem of Snapshot POD is solved to obtain eigenvalue and POD basis with every primary component. Subsequently the Reynolds Average Navier-Stokes equation (12)-(15) is derived to be used as the constraint functions. Regarding to the cost function, the sum of eigenvalues is defined as where a weight function δ j→k is intended to be extracted from the 1 − m th primary components to the j − k th primary components. Particulary, j and k denote the number of a primary component, and Based on a standard framework of finite element method and Lagrange Multiplier method, the objective functional L(x,ζ(x)) with a set of necessary variablesζ(x) = {ū,p,w,q} (x) is prepared with adjoint variablesw,q to a solutionū,p of the Reynolds Average Navier-Stokes equation. Next, the material derivative is taken with respect toζ(x) for the objective functional. The main problem (12)-(15) and the adjoint problem (20) The main and adjoint variablesζ(x) are substituted into L(ϕ 0 , ζ(x)) [ φ ] to evaluate the sensitivity with the Generalized J integration type (GJI-type) constructed in Nakazawa et al. (2021).

Numerical Schemes
For all numerical calculations, FreeFEM++ (Hecht, 2012) is used. Also, an UMFPACK solver (Davis, 2004) is used to solve the entire stiffness matrix obtained after discretization in space. We use GM and DGM as a spatial discretization scheme for solving the initial boundary value problem for Oldroyd-B. Then we compare both calculation results.
For the former, the numerical scheme suggested by Lukácová-Medvidiová et al. (2017) is used for Peterlin model. Particularly for the nonstationary Navier-Stokes equation, the first-order characteristic curve method (Notsu and Tabata, 2015) is introduced to discretize the material derivative term spatially and temporally. Brezzi-Pitkäranta pressure stabilization (Brezzi-Pitkäranta,1984) is used for spatial discretization of the pressure term. For the constitutive equation, the first-order characteristic curve method (Notsu and Tabata, 2015) is used for the material derivative term, too. For the latter, Hybridized DGM (HDGM; Oikawa, 2014) is adopted to discretize the initial boundary value problem spatially, where the upwind scheme for the advection term is implemented automatically. For temporal discretization, the implicit method is introduced.
To solve an eigenvalue problem of Snapshot POD and the adjoint problem, and to perform H 1 gradient method, the same numerical schemes as those presented by Nakazawa et al. (2021) are adopted.

Numerical Calculations and Discussion
For numerical calculations, the author adopts the two-dimensional cavity flow with a disk-shaped isolated body shown as where the disk boundary ∂Ω m is used as the design boundary. For the Dirichlet boundary condition Γ D , we defined The Neumann boundary condition Γ N is not used here.

Validations of spatial discretization, sampling time, artificial viscosity and Weissenberg number
For validations of spatial discretization, sampling time and artificial viscosity, the time step size ∆t = 0.001, is fixed to take time integrations of the initial-boundary problem at Re = 100 below.

Numerical calculation for shape optimization problems
Based on the result presented in the subsection above, this paper presents a solution for the suggested shape optimization problem and presents comparison of results obtained using both GM and HDGM under ∆t = 0.001, Re = 100, NBVX=5477, ϵ 1 = 10 −7 and (T 1 , T 2 ) = (3, 6). The author changes α from 0.1∆t −1 to 0.025∆t −1 , where Case 1, Case 2, and Case 3 respectively represent δ 1→1 with the time average part, δ 1→m with the time average and time fluctuation parts, and δ 2→m with the time fluctuation part.
First, GM is used as spatial discretizing scheme. The cost functions are converging at 9 and 25 step only for Case 3 at α = 0.1∆t −1 , 0.075∆t −1 as shown in Fig. 2. The optimal shapes and the sensitivities are visualized in Fig. 4 of Section 6.1. However, the cost functions are not converging for all cases at α = 0.05∆t −1 , 0.025∆t −1 .
Next, HDGM is used as a spatial discretizing scheme. The cost function is converging for all cases. Comparison with Fig. 4, where all partial differential equations are discretized in GM, reveals that the behavior of the sum of eigenvalues with the time fluctuation parts in Fig. 3 is not smoothly, but monotonically non-decreasing. Figures 5, 6, 7, and 8 of Section 6.2 portray optimal shapes of Case 1, Case 2, and Case 3 at α = 0.1∆t −1 , α = 0.075∆t −1 , α = 0.05∆t −1 and α = 0.025∆t −1 .  Table 6 Eigenvalues of Snapshot POD for different Weissenberg numbers α = 1 Wi in the initial domain at the first and the second primary components in the GM

Conclusion
As described in this paper, the author compares GM and HDGM for spatial discretization to demonstrate the suggested shape optimization problem for a complex flow. For numerical demonstrations, the two-dimensional cavity flow with the disk-shape isolated body is adopted as an initial domain. Oldroyd-B model is chosen as a represented complex flow. The Oldroyd-B model is adopted because Lukácová-Medvidiová et al. (2017) provided a numerical scheme in GM that is greatly beneficial for comparing numerical results with one reported in HDGM by Oikawa (2013 ). As shown in Fig. 1, HDGM, with which the upwind scheme for evaluating a convection term stably is involved, is more useful for calculation in the case of increasing Weissenberg Number than GM. In fact, by performing the shape optimization problem, HDGM can suppress the time fluctuation part more than GM. This paper is focusing on constructing a numerical scheme of a shape optimzation problem for suppressing transient Non-Newtonian fluid, and does not mention about a physical meaning of an optimal shape. 6. Appendix 6.1. Optimal shapes and sensitivity in GM