Shape optimization of an elliptic cylinder located in adiabatic flows

Abstract

The purpose of this paper is to present a shape optimization method for the drag force reduction problem of an elliptic cylinder located in fluid flows. The fluid is assumed to be in an adiabatic state. The conservation equations of mass and momentum are used with the density and momentum employed as field variables. For the discretization, the finite element method is applied using SUPG stabilization based on a linear interpolation function is used. The performance function consists of the square sum of the drag forces, integrated between the starting and final times. The coordinates of the surface of the cylinder are updated to obtain the minimal performance function. For this minimization procedure, we employ the weighted gradient method using the solution of the adjoint equation. A numerical computation is executed to demonstrate that the proposed method is effective for the drag reduction problem.

Keywords:

• shape optimization
• elliptic cylinder
• adiabatic fluid flow
• mass and momentum formulation
• adjoint equation method
• finite element method

Public Interest Statement

The characteristic feature of this paper is to present a numerical method to obtain the minimum drag force subjected to a cylindrical structure located in adiabatic flows, which is the flows with compressibility and without temperature change. For the computation of velocity and density, which is a function of pressure, the finite element method with the SUPG stabilization is effectively used. For the minimization procedure, the adjoint equation method based on the weighted gradient method is employed. To ensure the numerical algorithm, the drag force minimization problem has been implemented. Applications of the method presented in this paper is widely spread to the drag minimization in all engineering field.

1. Introduction

A number of published studies regarding the shape optimization of a cylinder located in fluid flows have considered the flows as incompressible, adiabatic, and compressible using the finite element method (Azegami, 2006; Ishiyama & Kawahara, 2008; Katamine, Azegami, & Ito, 2005; Mohammadi & Pironneau, 2009; Nakajima & Kawahara, 2008; Nakajima & Kawahara, 2010; Ogawa & Kawahara, 2003; Okumura & Kawahara, 2000; Okumura, Hikino & Kawahara, 2013; Sakamoto & Kawahara, 2010; Yagi & Kawahara, 2005; Yagi & Kawahara, 2007; Yoshida & Kawahara, 2008). Fluid flow problems are often solved under the assumption of incompressible flows. On the other hand, there are many problems in which we should consider density changes, for example, in mechanical and wind engineering. Essentially, density changes induce changes in temperature. However, there are many problems in which it is possible to formulate governing equations that incorporate density changes, but only include temperature changes implicitly. In this paper, we call these adiabatic flows.

The purpose of this paper is to present a shape optimization method based on the idea that the fluid flows are assumed to be in an adiabatic state, and are analyzed by employing the density and momentum as field variables. The previous studies by the authors (Kawahara, 2016; Nasu, Nojima & Kawahara, 2013; etc.) have detected out that for the evaluation of the fluid force applied to the cylinder located in fluid flows, we must consider the compressibility, but, need not consider the temperature change especially in case of low Reynolds number flows. Several methods have been presented for solving adiabatic flows using finite element methods, such as the two step method (Kawahara & Hirano, 1983, 1983; Kawahara & Miwa, 1984), bubble function method (Nakajima & Kawahara, 2008, 2010; Okumura & Kawahara, 2000; Okumura et al., 2013; Sakamoto & Kawahara, 2010; Yagi & Kawahara, 2005, 2007, 2009, 2008), and SUPG method (Kawahara, 2016; Maruoka, Uchiyama & Kawahara, 2017; Nasu et al., 2013). The method employed in this paper is the SUPG method based on the primitive variables, the formulation of which is entirely innovative. The shape optimization problem concerning an elliptic cylinder located in adiabatic flows is discussed, considering the density and momentum as field variables. The basic equations constitute the conservation equations of mass and momentum, and the Poisson law for the constitutive equations. The finite element method with SUPG stabilization with primitive variables is employed for the discretization.

In this paper, the optimum shape of a cylinder refers the shape on which the acting fluid forces are minimized. To measure the fluid forces, the performance function is introduced, which is a square sum of the fluid forces integrated from the initial to the final time. To determine the optimum shape, the coordinates of the cylinder surface are regulated as the performance function tends toward the minimum. By calculating the first variation of the performance function with respect to the surface coordinates and equating this with zero, the adjoint equation and the gradient can be derived. The simple weighted gradient method is employed for the minimization procedures. We have carried out many studies (e.g. Yamazaki & Kawahara, 2009) comparing the efficiency between by the single-order adjoint method and by the second-order adjoint method. We have the conclusion that the single-order adjoint method is simple and efficient. The shape of the cylinder is updated according to the gradient of the performance function with respect to the surface coordinates. The volumeof the cylinder should be kept constant. A mesh refinement is also introduced, based on the Delaunay triangulation medthod.

A drag force minimization method for an elliptic cylinder located in adiabatic fluid flows is presented. The streamline shape is finally obtained. The comparison between the computed results and experimental results Kawahara (2016) shows good agreement. Thus, we demonstrate that the proposed method is suitable for application to the shape optimization of a cylinder located in fluid flows.

2. Basic equation

The indicial notation and the summation convention are introduced for expressing equations. The conservations of mass and momentum are employed as the basic equations, and can be expressed using the density $\rho$ and velocity $u_i(i=1,2)$ in computational domain $\mathrm{\Omega }$ as follows:(1) $$\begin{array}{ccc}\hfill \dot{\rho }+{(\rho u_i)}_{,i}& =0\hfill & \hfill \text{in}\mathrm{\Omega },\end{array}$$(1) (2) $$\begin{array}{ccc}\hfill \rho ({\dot{u}}_i+u_j{u}_{i,j})+p_{,i}-{\tau }_{ij,j}& =0\hfill & \hfill \text{in}\mathrm{\Omega },\end{array}$$(2)

where p denotes pressure and ${\tau }_{\mathit{ij}}$ are viscous stress given by:(3) $$\begin{array}{c}\hfill {\tau }_{\mathit{ij}}=\mu (u_{i,j}+u_{j,i})+\lambda u_{k,k}{\delta }_{\mathit{ij}}.\end{array}$$(3)

In Equation (3), $\mu$ is viscosity, and $\lambda$ is bulk viscosity, and ${\delta }_{\mathit{ij}}$ represents the Kronecker delta, where:(4) $$\begin{array}{cc}\hfill {\delta }_{\mathit{ij}}& =\left\{\begin{array}{cc}1\hfill & (i=j)\hfill \\ 0\hfill & (i\ne j)\hfill \end{array}\right.,\hfill \end{array}$$(4) (5) $$\begin{array}{cc}\hfill \lambda & =-\frac{2}{3}\mu .\hfill \end{array}$$(5)

Under the assumption of an adiabatic state, the state equation can be expressed as follows:(6) $$\begin{array}{c}\hfill \left(\frac{p}{p_0}\right)={\left(\frac{\rho }{{\rho }_0}\right)}^{\gamma },\end{array}$$(6)

where p is pressure, $\rho$ is density, and $\gamma$ is ratio of specific heat. The reference values of the pressure and density are denoted by $p_0$ and ${\rho }_0$, respectively. The acoustic velocity c and Mach number $M_a$ are given by(7) $$\begin{array}{c}\hfill c=\sqrt{\frac{\gamma p}{\rho }},\end{array}$$(7) (8) $$\begin{array}{c}\hfill M_a=\frac{U}{c}.\end{array}$$(8)

where U is the reference velocity. Precise derivations of Equations (1)–(8) can be found in Kawahara (2016).

3. Discretization

The finite element method with SUPG stabilization is employed for the spatial discretization. Linear functions based on triangular finite elements are employed for the interpolation functions, which are expressed as ${\mathrm{\Phi }}_{\alpha }(\alpha =1-3)$. These interpolation functions are used for both the density and momentum. The pressure p is expressed as(9) $$\begin{array}{c}\hfill p_{,i}=\gamma \frac{p_0}{{\rho }_0}{\left(\frac{\rho }{p_0}\right)}^{\gamma -1}{\rho }_{,i}=\kappa {\rho }_{,i}\end{array}$$(9)

where $\kappa$ is assumed to be constant. Using weighting functions ${\rho }^{\ast }$ and $u_i^{\ast }$, the weighted residual equations of Equations (1) and (2) are derived as follows:(10) $$\begin{array}{cc}\hfill {\int }_{\mathrm{\Omega }}{\rho }^{\ast }\left\{\dot{\rho }+{(\rho u_j)}_{,j}\right\}d\mathrm{\Omega }& +\hfill \\ \hfill \sum _{e=1}^N{\int }_{{\mathrm{\Omega }}_e}{\tau }_s{\rho }_{,i}^{\ast }{a}_i\left\{\dot{\rho }+{(\rho u_j)}_{,j}\right\}d{\mathrm{\Omega }}_e& +\hfill \\ \hfill \sum _{e=1}^N{\int }_{{\mathrm{\Omega }}_e}{\tau }_s{\rho }_{,i}^{\ast }{a}_{\mathit{ij}}\left\{\rho ({\dot{u}}_j+u_k{u}_{j,k}+p_{,j})+p_{,j}\right\}d{\mathrm{\Omega }}_e& =0,\hfill \end{array}$$(10) (11) $$\begin{array}{cc}\hfill {\int }_{\mathrm{\Omega }}{u}_i^{\ast }\left\{\rho ({\dot{u}}_i+u_j{u}_{i,j})+p_{,i}-{\tau }_{ij,j}\right\}d\mathrm{\Omega }& +\hfill \\ \hfill \sum _{e=1}^N{\int }_{{\mathrm{\Omega }}_e}{\tau }_s{u}_{i,k}^{\ast }{b}_{\mathit{il}}^k\left\{\rho (\dot{u_l}+u_j{u}_{l,j})+p_{,l}\right\}d{\mathrm{\Omega }}_e& +\hfill \\ \hfill \sum _{e=1}^N{\int }_{{\mathrm{\Omega }}_e}{\tau }_s{u}_{i,j}^{\ast }{b}_i^j\left\{\dot{\rho }+{(\rho u_k)}_{,k}\right\}d{\mathrm{\Omega }}_e& =0,\hfill \end{array}$$(11)

where ${\tau }_s$ is the stabilization parameter for the SUPG method. In Equations (10) and (11), the terms a and b are derived from the advection terms in Equations (1) and (2), and can be written as:

where $u=u_1$, $v=u_2$, and $i,j,k=1,2$. We can formulate the SUPG stabilization method following the similar procedures described precisely by Kawahara (2016). The final discretized equations based on Equations (10) and (11) can be derived in the following forms:(14) $$\begin{array}{cc}\hfill M_{\alpha \beta }{\dot{\rho }}_{\beta }+B_{\alpha \beta i}{m}_{\beta i}+{\tau }_s{a}_i(D_{\alpha i\beta }{\dot{\rho }}_{\beta }+C_{\alpha i\beta j}{m}_{\beta j})& +\hfill \\ \hfill {\tau }_s{a}_{\mathit{ij}}(D_{\alpha i\beta }{\dot{m}}_{\beta j}+E_{\alpha i\beta \gamma k}{u}_{\beta k}{m}_{\gamma j}+\kappa C_{\alpha i\beta j}{\rho }_{\beta })& =0,\hfill \end{array}$$(14) (15) $$\begin{array}{cc}\hfill M_{\alpha \beta }{\dot{m}}_{\beta i}+H_{\alpha \beta j\gamma }{u}_{\beta j}{m}_{\gamma i}-\kappa B_{\beta \alpha i}{\rho }_{\beta }+S_{\alpha i\beta j}{m}_{\beta j}& +\hfill \\ \hfill {\tau }_s{b}_{\mathit{il}}^k(D_{\alpha k\beta }{\dot{m}}_{\beta l}+E_{\alpha k\beta \gamma j}{u}_{\beta j}{m}_{\gamma l}+\kappa C_{\alpha k\beta l}{\rho }_{\beta })& +\hfill \\ \hfill {\tau }_s{b}_i^j(D_{\alpha j\beta }{\dot{\rho }}_{\beta }+C_{\alpha j\beta k}{m}_{\beta k})& =F_{\alpha i},\hfill \end{array}$$(15)

where$$\begin{array}{cc}& M_{\alpha \beta }={\int }_{\mathrm{\Omega }}({\mathrm{\Phi }}_{\alpha }{\mathrm{\Phi }}_{\beta })d\mathrm{\Omega },B_{\alpha \beta i}={\int }_{\mathrm{\Omega }}({\mathrm{\Phi }}_{\alpha }{\mathrm{\Phi }}_{\beta ,i})d\mathrm{\Omega },\hfill \\ \hfill & C_{\alpha i\beta j}={\int }_{\mathrm{\Omega }}({\mathrm{\Phi }}_{\alpha ,i}{\mathrm{\Phi }}_{\beta ,j})d\mathrm{\Omega },D_{\alpha i\beta }={\int }_{\mathrm{\Omega }}({\mathrm{\Phi }}_{\alpha ,i}{\mathrm{\Phi }}_{\beta })d\mathrm{\Omega },\hfill \\ \hfill & E_{\alpha i\beta \gamma k}={\int }_{\mathrm{\Omega }}({\mathrm{\Phi }}_{\alpha ,i}{\mathrm{\Phi }}_{\beta }{\mathrm{\Phi }}_{\gamma ,k})d\mathrm{\Omega },\hfill \\ \hfill & H_{\alpha \beta j\gamma }={\int }_{\mathrm{\Omega }}({\mathrm{\Phi }}_{\alpha }{\mathrm{\Phi }}_{\beta ,j}{\mathrm{\Phi }}_{\gamma })d\mathrm{\Omega },\hfill \\ \hfill & S_{\alpha i\beta j}=\nu \left\{{\int }_{\mathrm{\Omega }}({\mathrm{\Phi }}_{\alpha ,k}{\mathrm{\Phi }}_{\beta ,k})d\mathrm{\Omega }{\delta }_{i,j}\right.\hfill \\ \hfill & \left.+{\int }_{\mathrm{\Omega }}({\mathrm{\Phi }}_{\alpha ,j}{\mathrm{\Phi }}_{\beta ,i})d\mathrm{\Omega }-\frac{2}{3}{\int }_{\mathrm{\Omega }}({\mathrm{\Phi }}_{\alpha ,i}{\mathrm{\Phi }}_{\beta ,j})d\mathrm{\Omega }\right\}.\hfill \\ \hfill & F_{\alpha i}={\int }_{\mathrm{\Gamma }}{\mathrm{\Phi }}_{\alpha }(-p{\delta }_{i,j}+{\tau }_{i,j})n_{j}d\mathrm{\Gamma }.\hfill \end{array}$$

For the temporal discretization, the implicit scheme and forward difference method are adopted. Denoting the time point by n, the following approximation is introduced:(16) $$\begin{array}{ccc}\hfill {\rho }_{\beta }& \approx \hfill & \hfill {\rho }_{\beta }^{n+1},\end{array}$$(16) (17) $$\begin{array}{ccc}\hfill m_{\beta i}& \approx \hfill & \hfill m_{\beta i}^{n+1},\end{array}$$(17) (18) $$\begin{array}{ccc}\hfill {\dot{\rho }}_{\beta }& \approx \hfill & \hfill \frac{1}{\mathrm{\Delta }t}\left({\rho }_{\beta }^{n+1}-{\rho }_{\beta }^n\right),\end{array}$$(18) (19) $$\begin{array}{ccc}\hfill {\dot{m}}_{\beta i}& \approx \hfill & \hfill \frac{1}{\mathrm{\Delta }t}\left(m_{\beta i}^{n+1}-m_{\beta i}^n\right).\end{array}$$(19)

Figure 1. Computational domain.

4. Shape optimization

Consider an elliptic cylinder located in fluid flows as shown in Figure 1. The cylinder is subjected to fluid forces, consisting of drag and lift forces. The problem studied in this paper concerns which type of shape for the cylinder minimizes the drag force. As an index of the magnitude of drag force, a performance function is employed given by the square sum of the difference between the computed and objective forces. The coordinates on the surface of the cylinder are employed as the control variables. Thus, the problem is to find the surface coordinates of the cylinder that minimize the performance function under the constraints of the basic equations of fluid flows. In order to obtain the surface coordinates, the adjoint equation method is employed.

In this study, the minimization problem concerning fluid forces acting on a cylinder is considered. The performance function J is defined by the square sum of the fluid forces as follows:(20) $$\begin{array}{c}\hfill J=\frac{1}{2}{\int }_{t_0}^{t_f}\left(F_{\alpha i}-{\overline{F}}_{\alpha i}\right)Q_{\alpha i\beta j}\left(F_{\beta j}-{\overline{F}}_{\beta j}\right)dt,\end{array}$$(20)

where $Q_{\alpha i\beta j}$ denote weighting parameters, $F_{\alpha i}$ represent drag and lift forces, and ${\overline{F}}_{\alpha i}$ are the objective fluid forces. For instance, in the case of the minimum drag force problem, the objective forces are given by ${\overline{F}}_{\alpha i}=0$ and $Q_{\alpha 2\beta 2}=1$, with the others equal to zero. In the case of minimum lift force problem, ${\overline{F}}_{\alpha i}$ is taken to be symmetric on the upper and lower surfaces. In order to obtain the optimal shape, the performance function defined by the fluid forces should be minimized, satisfying constraint conditions given by the basic equations in Equations (14) and (15). The adjoint variables for the basic equation are defined as ${\eta }_{\alpha }$ and ${\xi }_{\alpha i}$. The performance function is extended by adding the inner product between the adjoint variables and Equations (14) and (15). The extended performance function $J^{\ast }$ is expressed as follows:(21) $$\begin{array}{cc}\hfill J^{\ast }& =\frac{1}{2}{\int }_{t_0}^{t_f}\left(F_{\alpha i}-{\overline{F}}_{\alpha i}\right)Q_{\alpha i\beta j}\left(F_{\beta j}-{\overline{F}}_{\beta j}\right)dt\hfill \\ \hfill & +{\int }_{t_0}^{t_f}{\eta }_{\alpha }\left(M_{\alpha \beta }{\dot{\rho }}_{\beta }+B_{\alpha \beta i}{m}_{\beta i}\right)dt\hfill \\ \hfill & +{\int }_{t_0}^{t_f}{\xi }_{\alpha i}\left(M_{\alpha \beta }{\dot{m}}_{\beta i}+H_{\alpha \beta j\gamma }{u}_{\beta j}{m}_{\gamma i}\right.\hfill \\ \hfill & \left.-\kappa B_{\beta \alpha i}{\rho }_{\beta }+S_{\alpha i\beta j}{m}_{\beta j}-F_{\alpha i}\right)dt.\hfill \end{array}$$(21)

Solving the minimization problem with the constraint conditions results in solving the stationary conditions of the extended performance function $J^{\ast }$. The first-order adjoint equation is expressed as follows:(22) $$\begin{array}{cc}\hfill \delta J^{\ast }& ={\int }_{t_0}^{t_f}\delta {\eta }_{\alpha }\left(M_{\alpha \beta }{\dot{\rho }}_{\beta }+B_{\alpha \beta i}{m}_{\beta i}\right)dt\hfill \\ \hfill & +{\int }_{t_0}^{t_f}\delta {\xi }_{\alpha i}\left(M_{\alpha \beta }{\dot{m}}_{\beta i}+H_{\alpha \beta j\gamma }{u}_{\beta j}{m}_{\gamma i}\right.\hfill \\ \hfill & \left.-\kappa B_{\beta \alpha i}{\rho }_{\beta }+S_{\alpha i\beta j}{m}_{\beta j}-F_{\alpha i}\right)dt\hfill \\ \hfill & +{\int }_{t_0}^{t_f}\delta {\rho }_{\beta }\left(-M_{\alpha \beta }{\dot{\eta }}_{\alpha }-\kappa B_{\beta \alpha i}{\xi }_{\alpha i}\right)dt\hfill \\ \hfill & +{\int }_{t_0}^{t_f}\delta m_{\beta i}(-M_{\alpha \beta }{\dot{\xi }}_{\alpha i}+H_{\alpha \beta j\gamma }{u}_{\beta j}{\xi }_{\alpha i}\hfill \\ \hfill & +S_{\alpha k\beta i}{\xi }_{\alpha k}+B_{\alpha \beta i}{\eta }_{\alpha })dt\hfill \\ \hfill & +{\int }_{t_0}^{t_f}\delta F_{\alpha i}\left\{{\xi }_{\alpha i}-Q_{\alpha i\beta j}\left(F_{\beta j}-{\overline{F}}_{\beta j}\right)\right\}dt\hfill \\ \hfill & +G_{\delta k}\delta X_{\delta k}.\hfill \end{array}$$(22)

where $X_{\delta k}$ denote the coordinates of the surface of cylinder on $\delta$ node in the k direction, and $G_{\delta k}$ are the derivatives of the performance function corresponding to $X_{\delta k}$.

The stationary condition means that the first variation of the extended performance function is equal to zero. Therefore, the first-order adjoint equations are obtained as follows:(23) $$\begin{array}{cc}\hfill M_{\alpha \beta }{\dot{\rho }}_{\beta }+B_{\alpha \beta i}{m}_{\beta i}& =0\text{in}\mathrm{\Omega },\hfill \end{array}$$(23) (24) $$\begin{array}{cc}\hfill M_{\alpha \beta }{\dot{m}}_{\beta i}+H_{\alpha \beta j\gamma }{u}_{\beta j}{m}_{\gamma i}-& \hfill \\ \hfill \kappa B_{\beta \alpha i}{\rho }_{\beta }+S_{\alpha i\beta j}{m}_{\beta j}& =F_{\alpha i}\text{in}\mathrm{\Omega },\hfill \end{array}$$(24) (25) $$\begin{array}{cc}\hfill -M_{\alpha \beta }{\dot{\eta }}_{\alpha }-\kappa B_{\beta \alpha i}{\xi }_{\alpha i}& =0\text{in}\mathrm{\Omega },\hfill \end{array}$$(25) (26) $$\begin{array}{cc}\hfill -M_{\alpha \beta }{\dot{\xi }}_{\alpha \beta }+H_{\alpha \beta j\gamma }{u}_{\beta j}{\xi }_{\alpha i}+& \hfill \\ \hfill S_{\alpha k\beta i}{\xi }_{\alpha k}+B_{\alpha \beta i}{\eta }_{\alpha }& =0\text{in}\mathrm{\Omega },\hfill \end{array}$$(26) (27) $$\begin{array}{cc}\hfill {\xi }_{\alpha i}-Q_{\alpha i\beta j}\left(F_{\beta j}-{\overline{F}}_{\beta j}\right)& =0\text{on}{\mathrm{\Gamma }}_4,\hfill \end{array}$$(27) (28) $$\begin{array}{cc}\hfill G_{\delta k}\delta X_{\delta k}& =0\text{in}\mathrm{\Omega },\hfill \end{array}$$(28)

Then, the gradient used to update the surface coordinate can be derived as(29) $$\begin{array}{c}\hfill \text{grad}{(J^{\ast })}_{\delta k}=G_{\delta k},\end{array}$$(29)

where $G_{\delta k}$ is expressed as follows:(30) $$\begin{array}{cc}\hfill G_{\delta k}& {\displaystyle ={\int }_{t_0}^{t_f}{\eta }_{\alpha }\left(\frac{\partial M_{\alpha \beta }}{\partial X_{\delta k}}{\dot{\rho }}_{\beta }+\frac{\partial B_{\alpha \beta j}}{\partial X_{\delta k}}{m}_{\beta j}\right)dt}\hfill \\ \hfill & {\displaystyle +{\int }_{t_0}^{t_f}{\xi }_{\alpha i}\left(\frac{\partial M_{\alpha \beta }}{\partial X_{\delta \beta }}{\dot{m}}_{\beta i}+\frac{\partial H_{\alpha \beta j\gamma }}{\partial X_{\delta k}}{u}_{\beta j}{m}_{\gamma i}\right.}\hfill \\ \hfill & \left.+\frac{\partial B_{\alpha \beta j}}{\partial X_{\delta k}}{\rho }_{\beta }+\frac{\partial S_{\alpha i\beta j}}{\partial X_{\delta k}}{m}_{\beta j}\right)dt.\hfill \end{array}$$(30)

5. Minimization

It is necessary to preserve the volumeof the cylinder at the lth iteration cycle. A constant volumeis preserved, maintaining the volumeof the whole domain as follows:(31) $$\begin{array}{c}\hfill \sum _{e=1}^N\left(A_e\left(X_{\delta k}^{(l)}\right)\right)-A_0=0\text{in}\mathrm{\Omega },\end{array}$$(31)

where $A_e(X_{\delta k}^{(l)})$ represents the volumeof each element, $A_0$ denotes the volumeof the initial domain, and $X_{\delta k}^{(l)}$ are the coordinates on the boundary of the cylinder at the lth iteration cycle.

A modified performance function $K^{(l)}$ can be introduced by adding a penalty term to the extended performance function, which is expressed as follows:(32) $$\begin{array}{cc}\hfill K^{(l)}& =J^{\ast (l)}+\frac{1}{2}{\int }_{{\mathrm{\Gamma }}_4}\left(X_{\xi i}^{(l+1)}-X_{\xi i}^{(i)}\right)\hfill \\ \hfill & ·W_{\xi i\delta k}^{(l)}\left(X_{\delta k}^{(l+1)}-X_{\delta k}^{(l)}\right)d\mathrm{\Gamma }\hfill \\ \hfill & +\zeta \left\{\sum _{e=1}^N\left(A_e\left(X_{\delta k}^{(l)}\right)\right)-A_0=0\right\},\hfill \end{array}$$(32)

where l and $W_{\xi i\delta k}^{(l)}$ are the iteration count and stabilizing weight, respectively, and $\zeta$ is an adjoint variable for the volumeconstraint. In the case that the modified performance function $K^{(l)}$ converges to a minimum, the penalty term can be equal to zero. Minimizing the modified performance function $K^{(l)}$ is equivalent to minimizing the extended performance function $J^{\ast }$. Thus,(33) $$\begin{array}{c}\hfill \delta K^{(l)}=0.\end{array}$$(33)

By virtue of the stationary condition, the update of the coordinates around the cylinder is carried out by applying the following update equation:(34) $$\begin{array}{c}\hfill W_{\xi i\delta k}^{(l)}{X}_{\delta k}^{(l+1)}=W_{\xi i\delta k}^{(l)}{X}_{\delta k}^{(l)}-\text{grad}{(J^{\ast })}_{\xi i}^{(l)},\end{array}$$(34)

where(35) $$\begin{array}{c}\hfill \text{grad}{(K^{(l)})}_{\xi i}=G_{\xi i}^{(l)}-\zeta \left(\sum _{e=1}^N\frac{\partial A_e\left(X_{\delta k}^{(l)}\right)}{\partial X_{\xi i}}\right).\end{array}$$(35)

It is necessary to determine the adjoint variable $\zeta$ in order to preserve the volumeof the cylinder at each iteration. The variable $\zeta$ is obtained by setting the difference between the initial volumeand the present volumeto be zero:(36) $$\begin{array}{c}\hfill A_0-\sum _{e=1}^N\left(A_e\left(X_{\delta k}^{(l+1)}\right)\right)=0\end{array}$$(36)

By substituting Equation (34) into Equation (36), an implicit equation for $\zeta$ can be obtained.

Figure 2. Computational domain.

Figure 3. Boundary conditions.

6. Numerical study

Here, a drag force minimization problem is computed for an elliptic cylinder located in an fluid flow in an adiabatic state. The Reynolds and Mach numbers are given by $Re=1$ and $M_a=0.5$. The ratio of specific heat is $\gamma =1.4$. The computational domain and the boundary conditions are shown in Figures 2 and 3, respectively. Constant flows with $u=1$ and $v=0$ are imposed on the inflow boundary. On the outflow boundary, the pressure is prescribed as $p=0$. On the boundaries at both sides, slip conditions are imposed, i.e., $v=0$. The velocities on the surface of the cylinder are assumed to be zero. The density is taken to be ${\rho }_0=1.0$. The finite element mesh is illustrated in Figures 4 and 5. The total numbers of nodes and elements are 6,188 and 12,128, respectively. The number of surface nodes of the cylinder is 128. An initial ratio of 1:2 for the vertical and horizontal lengths is employed. The adopted time increment is $\mathrm{\Delta }t=0.001$.

Figure 4. Finite element mesh.

Figure 5. Finite element mesh around the cylinder.

Figure 6. Performance function.

Figure 7. Time history of drag force.

Figure 8. Finite element mesh of initial shape.

Figure 9. Finite element mesh of optimal shape.

The variation of the performance function is illustrated in Figure 6. It can be seen that the ratio of J is reduced to 0.5881, which constitutes a $41.2\%$ reduction. The drag forces acting on the cylinder at the initial and final stages are illustrated in Figure 7. Figures 8 and 9 present the finite element mesh at the initial and final stages, respectively. A vertical ratio of 1:4.75 is attained. The shape at the final stage is the streamline shape. Figures 10 and 11 illustrate the pressure distribution. At the final stage, we observe that the pressure at the front part of the cylinder is significantly reduced.

Figure 10. Pressure distribution of initial shape.

Figure 11. Pressure distribution of optimal shape.

7. Conclusion

In this paper, a shape optimization problem was presented regarding an elliptic cylinder located in fluid flows in an adiabatic state, and solved using the finite element method. The basic equations are formulated using the density and momentum as field variables. The finite element method using SUPG stabilization based on linear interpolation functions was effectively applied to solve the basic equations. The square sum of the drag forces acting on the cylinder, integrated over the working time, was introduced as the performance function. The weighted gradient method based on the adjoint equation was employed for the optimization. Precise forms of the SUPG finite element method based on the primitive variables, the adjoint equations, and the gradient of the performance function were concretely formulated. Starting with an elliptic shape, a streamline shape was obtained by performing a numerical computation using the proposed method. Thus, it has been demonstrated that the shape optimization method proposed in this paper is effective for practical engineering problems.

8. Funding

The authors received no direct funding for this research.

Additional information

Notes on contributors

Yuta Mochizuki

Yuta Mochizuki is a student of Department of Civil Engineering, Chuo University, Tokyo, Japan. He is now working for Tokyo Metropolitan Railway Co. Ltd.

Mutsuto Kawahara

Mutsuto Kawahara is professor of Emeritus, Department of Civil Engineering, Chuo University, Tokyo, Japan. His main research subject was the development of the numerical method, especially finite element method in fluid flows and shallow water flows. The optimization problems including force reduction and optimal estimation, etc. are recent research subject.