Shape optimization of unsteady natural convection fields

Abstract

This paper presents a numerical method for two shape optimization problems, namely, prescribing the temperature history distribution on sub-boundaries and maximizing the discharged heat on sub-boundaries of unsteady natural convection fields. The square error integral between the actual temperature distribution and the target temperature distribution on the sub-boundaries during a specified period of time was used as the objective functional for the prescribed temperature history distribution. The shape gradients of these shape determination problems were derived theoretically using the Lagrange multiplier method, adjoint variable method, and the material derivative formulae. Reshaping was performed by the traction method, which was proposed as an approach for solving shape optimization problems. Numerical programs for the shape determination problems are developed based on FreeFem++ in order to verify the proposed method.

Keywords:

• Shape optimization
• shape identification
• optimum design
• inverse problem
• natural convection field
• adjoint method

AMS Subject Classifications:

• 65: Numerical analysis
• 65L09: Inverse problems

1. Introduction

Research on the shape optimization of flow fields began with Pironneau [1,2], who proposed a method for updating shapes by treating displacements in the normal direction of the boundary as design variables in the distribution system and evaluating the shape gradient functions (sensitivities) by selecting the node coordinates on the boundary of a finite-element model as design variables. When the finite element method (FEM) was used for numerical analysis, the mesh within the region had to be updated after the boundary was updated [3]. An unstable phenomenon involving the undulation of boundary shapes has been known to occur when shape gradient functions related to the displacement of such design variables were evaluated, and the resulting sensitivities were used to move nodes [4]. Jameson [5] proposed a numerical method for updating the boundary shapes using an elliptic equation defined on the boundary. This method replaced the actual sensitivities with smooth sensitivities in order to facilitate the use of the second-order differentiation conducted during the sensitivity evaluation in updating shapes in flow fields. Mohammadi et al. [6] also proposed a similar method.

The traction method [7,8] was proposed for solving equations governing domain variables for linear elastic fields. Theoretical consideration for retaining the boundary smoothness following changes is possible with this method, and the undulating phenomenon did not occur [9]. Furthermore, remeshing is not required for ordinary shape update analyses because the domain variation of all nodes within a region can be treated as design variables in the FEM while still retaining smoothness. This method enables sensitivity calculations of large-scale design variables for complex coupled fields based on the adjoint variable method [10].

Shape determinations in natural convection fields have been studied by Momose et al. [11], Park et al. [12,13], and Aounallah et al. [14]. Momose et al. [11] have proposed a shape sensitivity analysis method seeking to maximize the velocity in the sub-domain in natural convection fields; however, this method was limited to steady-state problems. Park et al. [12,13] used the adjoint variable method for shape identification problems in unsteady natural convection fields. They proposed to limit the design variables that express boundary shapes to the minimum number of variables required for the steady-state problem. Aounallah et al. [14] in two-dimensions expressed the equations governing unsteady natural convection fields using flow functions and vorticity. They conducted an analysis using Bezier curves, which facilitate expressions of smooth boundary shapes with a minimal number of design variables, similar to Park et al. [12,13]. Topology shape optimization analyses have recently been conducted on heat convection fields by Yaji et al. [15] and Alexandersen et al. [16]. Yaji et al. [15] optimized the topology of forced convection fields, while Alexandersen et al. [16] used the adjoint variable method to optimize the topology of natural convection fields. However, in all these cases, the analysis was limited to steady-state problems.

This study considers two shape determination problems for unsteady natural convection fields. The first problem involves identification of shapes when prescribing the temperature distribution history on the sub-boundaries in natural convection fields. The problem of minimizing the integrated squared error between the actual temperature distribution and the target temperature distribution on sub-boundary is formulated in Section 4.1, and the shape gradient is derived using either the Lagrange multiplier method or the adjoint variable method and the material derivative function in Section 4.2. The second problem involves the determination of shapes for maximizing the amount of heat discharged on a sub-boundary. The problem is formulated, and the shape gradient density function is derived in Section 5 in a similar manner to that of the first problem. The traction method, Section 6, is applied to examples in Section 7 based on the shape gradient density functions derived for each problem.

2. Domain variation

Before formulating the shape optimization problem, a method of representing domain variation using the speed method will be discussed briefly. A more detailed explanation is presented in [17].

Assume that a bounded domain $\mathrm{\Omega }$ of ${\mathbb{R}}^d(d=2,3)$ with boundary $\mathrm{\Gamma }$ is variable. One approach to describing the domain variation is to use a one-parameter family of one-to-one mappings ${\overrightarrow{T}}_s(\mathrm{\Omega }):\mathrm{\Omega }\mapsto {\mathrm{\Omega }}_s$ or its inverse ${\overrightarrow{T}}_s^{-1}({\mathrm{\Omega }}_s):{\mathrm{\Omega }}_s\mapsto \mathrm{\Omega }$, where s denotes the domain variation history. When a domain functional $J_{\mathrm{\Omega }}$ and a boundary functional $J_{\mathrm{\Gamma }}$ of a distributed function $\mathit{\psi }$ are considered, their derivatives ${\dot{J}}_{\mathrm{\Omega }}$ and ${\dot{J}}_{\mathrm{\Gamma }}$ with respect to s at $s=0$ are given by the formulae of the material derivative:(1) $$\begin{array}{cc}& J_{\mathrm{\Omega }}={\int }_{\mathrm{\Omega }}\mathit{\psi }\mathrm{d}x,{\dot{J}}_{\mathrm{\Omega }}={\int }_{\mathrm{\Omega }}{\mathit{\psi }}^{\prime }\mathrm{d}x+{\int }_{\mathrm{\Gamma }}\mathit{\psi }\overrightarrow{n}·\overrightarrow{V}\mathrm{d}x,\hfill \\ \hfill & J_{\mathrm{\Gamma }}={\int }_{\mathrm{\Gamma }}\mathit{\psi }\mathrm{d}\mathrm{\Gamma },{\dot{J}}_{\mathrm{\Gamma }}={\int }_{\mathrm{\Gamma }}\left\{{\mathit{\psi }}^{\prime }+\left({\mathrm{\nabla }}_n\mathit{\psi }+\mathit{\psi }\mathit{\kappa }\right)\overrightarrow{n}·\overrightarrow{V}\right\}\mathrm{d}\mathrm{\Gamma },\hfill \end{array}$$(1)

where $\overrightarrow{n}$ is an outward unit normal vector to the boundary, ${\mathrm{\nabla }}_n(·)\equiv \mathrm{\nabla }(·)·\overrightarrow{n}$, and $\mathit{\kappa }$ denotes the quantity ($d-1$) times the average boundary curvature. The shape derivative ${\mathit{\psi }}^{\prime }$ of the distributed function $\mathit{\psi }$ indicates that the derivatives are fixed in spatial coordinates. The derivative $\overrightarrow{V}({\mathrm{\Omega }}_s)$ of ${\overrightarrow{T}}_s(\mathrm{\Omega })$ with respect to s given by(2) $$\begin{array}{cc}& \overrightarrow{V}({\mathrm{\Omega }}_s)=\frac{\mathit{\partial }{\overrightarrow{T}}_s\left(\mathrm{\Omega }\right)}{\mathit{\partial }s}=\frac{\mathit{\partial }{\overrightarrow{T}}_s}{\mathit{\partial }s}\left({\overrightarrow{T}}_s^{-1}({\mathrm{\Omega }}_s)\right)\hfill \end{array}$$(2)

is referred to as the velocity because of the analogy between s and time.

3. Governing equations for unsteady natural convection fields

Consider the unsteady natural convection field in the region $\mathrm{\Omega }$ of ${\mathbb{R}}^d(d=2,3)$ in time interval [0, T]. Consider determining the flow velocity $u(\overrightarrow{x},t)={(u_i(\overrightarrow{x},t))}_{i=\overline{1,d}}$, pressure $p(\overrightarrow{x},t)$ and temperature $\mathit{\theta }(\overrightarrow{x},t)$ in space $\overrightarrow{x}\in \mathrm{\Omega }$ at time $t\in [0,T]$. The Boussinesq-approximated Navier–Stokes equation, continuity equation, and energy equation are the governing equations for unsteady natural convection fields. They can be expressed as follows:(3) $$\begin{array}{cc}& \mathit{\rho }\frac{\mathit{\partial }{u}_i}{\mathit{\partial }t}+\mathit{\rho }{u}_j{u}_{i,j}=-p_{,i}+\mathit{\mu }{u}_{i,jj}-e_i\mathit{\rho }\mathit{\beta }g(\mathit{\theta }-{\mathit{\theta }}_0),(\overrightarrow{x},t)\in \mathrm{\Omega }\times [0,T],\hfill \end{array}$$(3) (4) $$\begin{array}{cc}& u_{i,i}=0,(\overrightarrow{x},t)\in \mathrm{\Omega }\times [0,T],\hfill \end{array}$$(4) (5) $$\begin{array}{cc}& \frac{\mathit{\partial }\mathit{\theta }}{\mathit{\partial }t}+u_j{\mathit{\theta }}_{,j}=\mathit{\alpha }{\mathit{\theta }}_{,jj},(\overrightarrow{x},t)\in \mathrm{\Omega }\times [0,T],\hfill \end{array}$$(5)

where the boundary $\mathrm{\Gamma }=\mathit{\partial }\mathrm{\Omega }={\mathrm{\Gamma }}_u\cup {\mathrm{\Gamma }}_{\mathit{\sigma }}={\mathrm{\Gamma }}_{\mathit{\theta }}\cup {\mathrm{\Gamma }}_q\cup {\mathrm{\Gamma }}_h$, see Figure 1.

Figure 1. Natural convection field.

The boundary conditions and initial conditions are:(6) $$\begin{array}{cc}& u_i(\overrightarrow{x},t)={\widehat{u}}_i(\overrightarrow{x},t),t\in [0,T],\overrightarrow{x}\in {\mathrm{\Gamma }}_u,\hfill \end{array}$$(6) (7) $$\begin{array}{cc}& {\mathit{\sigma }}_i(\overrightarrow{x},t)={\widehat{\mathit{\sigma }}}_i(\overrightarrow{x},t)=(-p{\mathit{\delta }}_{ij}+\mathit{\mu }{u}_{i,j})n_j=0,t\in [0,T],\overrightarrow{x}\in {\mathrm{\Gamma }}_{\mathit{\sigma }},\hfill \end{array}$$(7) (8) $$\begin{array}{cc}& \mathit{\theta }(\overrightarrow{x},t)=\widehat{\mathit{\theta }}(\overrightarrow{x},t),t\in [0,T],\overrightarrow{x}\in {\mathrm{\Gamma }}_{\mathit{\theta }},\hfill \end{array}$$(8) (9) $$\begin{array}{cc}& -\mathit{\alpha }\mathit{\theta }{(\overrightarrow{x},t)}_{,j}{n}_j=\widehat{q}(\overrightarrow{x},t),t\in [0,T],\overrightarrow{x}\in {\mathrm{\Gamma }}_q,\hfill \end{array}$$(9) (10) $$\begin{array}{cc}& -\mathit{\alpha }\mathit{\theta }{(\overrightarrow{x},t)}_{,j}{n}_j=\widehat{h}(\mathit{\theta }(\overrightarrow{x},t)-{\widehat{\mathit{\theta }}}_f),t\in [0,T],\overrightarrow{x}\in {\mathrm{\Gamma }}_h,\hfill \end{array}$$(10) (11) $$\begin{array}{cc}& u_i(\overrightarrow{x},0)=u_{i_{ini}}(\overrightarrow{x}),\overrightarrow{x}\in \mathrm{\Omega },\hfill \end{array}$$(11) (12) $$\begin{array}{cc}& p(\overrightarrow{x},0)=p_{ini}(\overrightarrow{x}),\overrightarrow{x}\in \mathrm{\Omega },\hfill \end{array}$$(12) (13) $$\begin{array}{cc}& \mathit{\theta }(\overrightarrow{x},0)={\mathit{\theta }}_{ini}(\overrightarrow{x}),\overrightarrow{x}\in \mathrm{\Omega }.\hfill \end{array}$$(13)

Here, $\mathit{\rho }$ represents the density, ${\mathit{\theta }}_0$ represents the base temperature, $\mathit{\mu }$ represents the dynamic viscosity, $e_1=e_3=0,e_2=-1$, and g represents the acceleration of gravity. Moreover, $\mathit{\beta }$ represents the coefficient of thermal expansion, $\mathit{\alpha }$ represents the thermal diffusivity, $\widehat{q}$ and $\widehat{h}$ represent the heat flux and the heat transfer coefficient normalized by heat capacity, ${\widehat{\mathit{\theta }}}_f$ represents the external temperature, $u_{i_{}\text{ini}}$ represents the initial flow velocity, $p_{ini}$ represents the initial pressure, ${\mathit{\theta }}_{}\text{ini}$ represents the initial temperature, ${\mathit{\delta }}_{ij}$ represents the Kronecker delta. Finally, $(\widehat{·})$ represents the known function on the boundary.

The weak forms of the respective governing Equations (4(4) $$\begin{array}{cc}& u_{i,i}=0,(\overrightarrow{x},t)\in \mathrm{\Omega }\times [0,T],\hfill \end{array}$$(4) )–(6(6) $$\begin{array}{cc}& u_i(\overrightarrow{x},t)={\widehat{u}}_i(\overrightarrow{x},t),t\in [0,T],\overrightarrow{x}\in {\mathrm{\Gamma }}_u,\hfill \end{array}$$(6) ) can be expressed with adjoint flow velocity $w(\overrightarrow{x},t)={(w_i(\overrightarrow{x},t))}_{i=\overline{1,d}}$, adjoint pressure $q(\overrightarrow{x},t)$ and adjoint temperature $\mathit{\xi }(\overrightarrow{x},t)$ as follows:(14) $$\begin{array}{cc}& {\int }_0^T\{t^V(u_{,t},w)+a^V(u,w)+b^V(u,u,w)+c(w,p)+f^V(w,\mathit{\theta })-l(w)\}\mathrm{d}t=0,\hfill \\ \hfill & \forall w\in W,\hfill \end{array}$$(14) (15) $$\begin{array}{cc}& {\int }_0^T\{c(u,q)\}\mathrm{d}t=0,\forall q\in Q,\hfill \end{array}$$(15) (16) $$\begin{array}{cc}& {\int }_0^T\{t^H({\mathit{\theta }}_{,t},\mathit{\xi })+a^H(\mathit{\theta },\mathit{\xi })+b^H(u,\mathit{\theta },\mathit{\xi })+f_q^H(\mathit{\xi })+f_h^H(\mathit{\theta },\mathit{\xi })-f_{hf}^H(\mathit{\xi })\}\mathrm{d}t=0,\hfill \\ \hfill & \forall \mathit{\xi }\in \mathrm{\Xi }.\hfill \end{array}$$(16)

Furthermore, $t^V(u_{,t},w)$, $a^V(u,w)$, $b^V(v,u,w)$, c(w, p), $f^V(w,\mathit{\theta })$, l(w), $t^H({\mathit{\theta }}_{,t},\mathit{\xi })$, $a^H(\mathit{\theta },\mathit{\xi })$, $b^H(u,\mathit{\theta },\mathit{\xi })$, $f_q(\mathit{\xi })$, $f_h^H(\mathit{\theta },\mathit{\xi })$, and $f_{hf}^H(\mathit{\xi })$ are defined as follows:(17) $$\begin{array}{cc}\hfill t^V(u_{,t},w)& ={\int }_{\mathrm{\Omega }}\mathit{\rho }{w}_i\frac{\mathit{\partial }{u}_i}{\mathit{\partial }t}\mathrm{d}x,a^V(u,w)={\int }_{\mathrm{\Omega }}\mathit{\mu }{w}_{i,j}{u}_{i,j}\mathrm{d}x,\hfill \\ \hfill b^V(v,u,w)& ={\int }_{\mathrm{\Omega }}\mathit{\rho }{w}_i{v}_j{u}_{i,j}\mathrm{d}x,c(w,p)=-{\int }_{\mathrm{\Omega }}{w}_{i,i}p\mathrm{d}x,\hfill \\ \hfill f^V(w,\mathit{\theta })& ={\int }_{\mathrm{\Omega }}{w}_i{e}_i\mathit{\rho }\mathit{\beta }g(\mathit{\theta }-{\mathit{\theta }}_0)\mathrm{d}x,l(w)={\int }_{{\mathrm{\Gamma }}_{\mathit{\sigma }}}{w}_i{\widehat{\mathit{\sigma }}}_i\mathrm{d}\mathrm{\Gamma },\hfill \\ \hfill t^H({\mathit{\theta }}_{,t},\mathit{\xi })& ={\int }_{\mathrm{\Omega }}\mathit{\xi }\frac{\mathit{\partial }\mathit{\theta }}{\mathit{\partial }t}\mathrm{d}x,a^H(\mathit{\theta },\mathit{\xi })={\int }_{\mathrm{\Omega }}\mathit{\alpha }{\mathit{\xi }}_{,i}{\mathit{\theta }}_{,i}\mathrm{d}x,b^H(u,\mathit{\theta },\mathit{\xi })={\int }_{\mathrm{\Omega }}\mathit{\xi }{u}_j{\mathit{\theta }}_{,j}\mathrm{d}x,\hfill \\ \hfill f_q^H(\mathit{\xi })& ={\int }_{{\mathrm{\Gamma }}_q}\mathit{\xi }\widehat{q}\mathrm{d}\mathrm{\Gamma },f_h^H(\mathit{\theta },\mathit{\xi })={\int }_{{\mathrm{\Gamma }}_h}\widehat{h}\mathit{\xi }\mathit{\theta }\mathrm{d}\mathrm{\Gamma },f_{hf}^H(\mathit{\xi })={\int }_{{\mathrm{\Gamma }}_h}\widehat{h}\mathit{\xi }{\widehat{\mathit{\theta }}}_f\mathrm{d}\mathrm{\Gamma }.\hfill \end{array}$$(17)

Here, the superscripts V and H denote terms about the viscosity and heat, the subscripts q and h, hf denote terms about the heat flux and the heat transfer. Tensors described in this paper use the Einstein summation convention and differentiation ${(·)}_{,i}=\mathit{\partial }(·)/\mathit{\partial }{x}_i$, while ${(·)}_{,t}$ expresses the time derivative of the function.

The flow velocity u, its adjoint w, and the other variables are considered to be elements of the following functional spaces:(18) $$\begin{array}{cc}\hfill U& =\{u(\overrightarrow{x},t)\in H^1{(\mathrm{\Omega }\times [0,T])}^d|u\text{satisfies}\text{(7) and}\text{(12)}\},\hfill \end{array}$$(18) (19) $$\begin{array}{cc}\hfill Q& =\{q(\overrightarrow{x},t)\in L^2(\mathrm{\Omega }\times [0,T])|{\int }_{\mathrm{\Omega }}q\mathrm{d}x=0\text{if measure}({\mathrm{\Gamma }}_{\mathit{\sigma }})=0\},\hfill \end{array}$$(19) (20) $$\begin{array}{cc}\hfill \mathrm{\Theta }& =\{\mathit{\theta }(\overrightarrow{x},t)\in H^1(\mathrm{\Omega }\times [0,T])|\mathit{\theta }\text{satisfies}\text{(9)--(11),}\text{(14)}\},\hfill \end{array}$$(20) (21) $$\begin{array}{cc}\hfill W& =\{w(\overrightarrow{x},t)\in H^1{(\mathrm{\Omega }\times [0,T])}^d|w(\overrightarrow{x},t)=0,t\in [0,T],\overrightarrow{x}\in {\mathrm{\Gamma }}_u,\hfill \\ \hfill & w(\overrightarrow{x},T)=0,\overrightarrow{x}\in \mathrm{\Omega }\},\hfill \end{array}$$(21) (22) $$\begin{array}{cc}\hfill \mathrm{\Xi }& =\{\mathit{\xi }(\overrightarrow{x},t)\in H^1(\mathrm{\Omega }\times [0,T])|\mathit{\xi }(\overrightarrow{x},t)=0,t\in [0,T],\overrightarrow{x}\in {\mathrm{\Gamma }}_{\mathit{\theta }},\hfill \\ \hfill & \mathit{\xi }(\overrightarrow{x},T)=0,\overrightarrow{x}\in \mathrm{\Omega }\}.\hfill \end{array}$$(22)

4. Prescribing temperature on a sub-boundary

4.1. Problem formulation

In this section, the problem of minimizing the square integration errors between the actual temperature ${\mathit{\theta }|}_{{\mathrm{\Gamma }}_D\times [t_1,t_2]}$ from time $t=t_1\in [0,T]$ to $t=t_2\in [0,T]$ and the target temperature ${\mathit{\theta }}_D{|}_{{\mathrm{\Gamma }}_D\times [t_1,t_2]}$ on sub-boundary ${\mathrm{\Gamma }}_D\subset \mathrm{\Gamma }$ is formulated. We assume $t_1<t_2$. The domain variation of this natural convection field region $\mathrm{\Omega }$ is denoted by $\overrightarrow{T_s}$, and the domain $\mathrm{\Omega }$ is assumed to vary to reach ${\mathrm{\Omega }}_s=\overrightarrow{T_s}(\mathrm{\Omega })$. For simplicity, we assume that the sub-boundaries ${\mathrm{\Gamma }}_D$, ${\mathrm{\Gamma }}_{\mathit{\sigma }}$ and ${\mathrm{\Gamma }}_h$ are invariable, that is $\overrightarrow{T_s}({\mathrm{\Gamma }}_D)={\mathrm{\Gamma }}_D$, $\overrightarrow{T_s}({\mathrm{\Gamma }}_{\mathit{\sigma }})={\mathrm{\Gamma }}_{\mathit{\sigma }}$ and $\overrightarrow{T_s}({\mathrm{\Gamma }}_h)={\mathrm{\Gamma }}_h$ for domain variation. The square integration error problem for temperature distribution from time $t=t_1$ to $t=t_2$ is formulated as follows:(23) $$\begin{array}{cc}& \text{Given}\mathrm{\Omega }\text{find}{\mathrm{\Omega }}_s\text{that minimizes}{\int }_{t_1}^{t_2}{E}_{{\mathrm{\Gamma }}_D}(\mathit{\theta })\mathrm{d}t\text{subject to}\text{(15)--(17)}\text{and}\hfill \\ \hfill & {\int }_{\mathrm{\Omega }}\mathrm{d}x\le {\mathit{\beta }}_{V}M,\hfill \end{array}$$(23)

where ${\mathit{\beta }}_V$ is a coefficient related to the initial domain measure M, and(24) $$\begin{array}{c}\hfill E_{{\mathrm{\Gamma }}_D}(\mathit{\theta })={\int }_{{\mathrm{\Gamma }}_D}{(\mathit{\theta }-{\mathit{\theta }}_D)}^2\mathrm{d}\mathrm{\Gamma }.\end{array}$$(24)

4.2. Shape gradient function

The square integration error problem can be rewritten as a retention problem without any restrictions from the Lagrange multiplier method. In this instance, the Lagrange function $L(u_i,p,\mathit{\theta },w_i,q,\mathit{\xi },\mathrm{\Lambda })$ is given as follows:(25) $$\begin{array}{cc}\hfill L=& {\int }_{t_1}^{t_2}{E}_{{\mathrm{\Gamma }}_D}(\mathit{\theta })\mathrm{d}t\hfill \\ \hfill & -{\int }_0^T\{t^V(u_{,t},w)+a^V(u,w)+b^V(u,u,w)+c(w,p)+f^V(w,\mathit{\theta })-l(w)\}\mathrm{d}t\hfill \\ \hfill & -{\int }_0^T\{c(u,q)\}\mathrm{d}t\hfill \\ \hfill & -{\int }_0^T\{t^H({\mathit{\theta }}_{,t},\mathit{\xi })+a^H(\mathit{\theta },\mathit{\xi })+b^H(u,\mathit{\theta },\mathit{\xi })+f_q^H(\mathit{\xi })+f_h^H(\mathit{\theta },\mathit{\xi })-f_{hf}^H(\mathit{\xi })\}\mathrm{d}t\hfill \\ \hfill & +\mathrm{\Lambda }\left({\int }_{\mathrm{\Omega }}\mathrm{d}x-{\mathit{\beta }}_{V}M\right).\hfill \end{array}$$(25)

where $w\in W$, $q\in Q$, and $\mathit{\xi }\in \mathrm{\Xi }$ were introduced as Lagrange multiplier functions or the adjoint functions, with respect to the weak forms. The non-negative real constant number $\mathrm{\Lambda }$ is Lagrange multiplier with respect to the volume constraint. The derivative of L with respect to domain variation is derived using the velocity field (3(3) $$\begin{array}{cc}& \mathit{\rho }\frac{\mathit{\partial }{u}_i}{\mathit{\partial }t}+\mathit{\rho }{u}_j{u}_{i,j}=-p_{,i}+\mathit{\mu }{u}_{i,jj}-e_i\mathit{\rho }\mathit{\beta }g(\mathit{\theta }-{\mathit{\theta }}_0),(\overrightarrow{x},t)\in \mathrm{\Omega }\times [0,T],\hfill \end{array}$$(3) ), as follows, [7,8],(26) $$\begin{array}{cc}& \dot{L}=\hfill \\ \hfill & -{\int }_0^T\{t^V(u_{,t},w^{\prime })+a^V(u,w^{\prime })+b^V(u,u,w^{\prime })+c(w^{\prime },p)+f^V(w^{\prime },\mathit{\theta })-l(w^{\prime })\}\mathrm{d}t\hfill \\ \hfill & -{\int }_0^T\{c(u,q^{\prime })\}\mathrm{d}t\hfill \\ \hfill & -{\int }_0^T\{t^H({\mathit{\theta }}_{,t},{\mathit{\xi }}^{\prime })+a^H(\mathit{\theta },{\mathit{\xi }}^{\prime })+b^H(u,\mathit{\theta },{\mathit{\xi }}^{\prime })+f_q^H({\mathit{\xi }}^{\prime })+f_h^H(\mathit{\theta },{\mathit{\xi }}^{\prime })-f_{hf}^H({\mathit{\xi }}^{\prime })\}\mathrm{d}t\hfill \\ \hfill & -{\int }_0^T\{t^V(u_{,t}^{\prime },w)+a^V(u^{\prime },w)+b^V(u^{\prime },u,w)+b^V(u,u^{\prime },w)+c(u^{\prime },q)\hfill \\ \hfill & +b^H(u^{\prime },\mathit{\theta },\mathit{\xi })\}\mathrm{d}t\hfill \\ \hfill & -{\int }_0^T\{c(w,p^{\prime })\}\mathrm{d}t\hfill \\ \hfill & -{\int }_0^T\{f^V(w,{\mathit{\theta }}^{\prime })+t^H({\mathit{\theta }}_{,t}^{\prime },\mathit{\xi })+a^H({\mathit{\theta }}^{\prime },\mathit{\xi })+b^H(u,{\mathit{\theta }}^{\prime },\mathit{\xi })+f_h^H({\mathit{\theta }}^{\prime },\mathit{\xi })\}\mathrm{d}t\hfill \\ \hfill & +{\int }_{t_1}^{t_2}{E}_{{\mathrm{\Gamma }}_D}({\mathit{\theta }}^{\prime })\mathrm{d}t\hfill \\ \hfill & +\dot{\mathrm{\Lambda }}\left({\int }_{\mathrm{\Omega }}\mathrm{d}x-{\mathit{\beta }}_{V}M\right)+l_G(\overrightarrow{V}).\hfill \end{array}$$(26)

Here, ${(·)}^{\prime }$ represents the derivative with respect to domain variation of the function fixed on the spatial coordinates, and(27) $$\begin{array}{cc}& l_G(\overrightarrow{V})={\int }_{{\mathrm{\Gamma }}_u}G\overrightarrow{n}·\overrightarrow{V}\mathrm{d}\mathrm{\Gamma },\hfill \end{array}$$(27)

where(28) $$\begin{array}{cc}& G={\int }_0^T\left\{-\mathit{\mu }{w}_{i,j}{u}_{i,j}+w_{i,i}p-\frac{\mathit{\partial }\mathit{\theta }}{\mathit{\partial }t}\mathit{\xi }-\mathit{\alpha }{\mathit{\xi }}_{,i}{\mathit{\theta }}_{,i}-\mathit{\xi }{u}_j{\mathit{\theta }}_{,j}-{\mathrm{\nabla }}_n(\mathit{\xi }\widehat{q})-(\mathit{\xi }\widehat{q})\mathit{\kappa }\right\}\mathrm{d}t+\mathrm{\Lambda },\hfill \\ \hfill \end{array}$$(28)

and $u_i,p,\mathit{\theta },w_i,q,\mathit{\xi }$, and $\mathrm{\Lambda }$ are determined by the following conditions:(29) $$\begin{array}{cc}& {\int }_0^T\{t^V(u_{,t},w^{\prime })+a^V(u,w^{\prime })+b^V(u,u,w^{\prime })+c(w^{\prime },p)+f^V(w^{\prime },\mathit{\theta })-l(w^{\prime })\hfill \\ \hfill & +c(u,q^{\prime })\}\mathrm{d}t=0,\forall w^{\prime }\in W,\forall q^{\prime }\in Q,\hfill \end{array}$$(29) (30) $$\begin{array}{cc}& {\int }_0^T\{t^H({\mathit{\theta }}_{,t},{\mathit{\xi }}^{\prime })+a^H(\mathit{\theta },{\mathit{\xi }}^{\prime })+b^H(u,\mathit{\theta },{\mathit{\xi }}^{\prime })+f_q^H({\mathit{\xi }}^{\prime })+f_h^H(\mathit{\theta },{\mathit{\xi }}^{\prime })\hfill \\ \hfill & -f_{hf}^H({\mathit{\xi }}^{\prime })\}\mathrm{d}t=0,\forall {\mathit{\xi }}^{\prime }\in \mathrm{\Xi },\hfill \end{array}$$(30) (31) $$\begin{array}{cc}& {\int }_0^T\{t^V(u_{,t}^{\prime },w)+a^V(u^{\prime },w)+b^V(u^{\prime },u,w)+b^V(u,u^{\prime },w)+c(u^{\prime },q)+b^H(u^{\prime },\mathit{\theta },\mathit{\xi })\hfill \\ \hfill & +c(w,p^{\prime })\}\mathrm{d}t=0,\forall u^{\prime }\in U,\forall p^{\prime }\in Q,\hfill \end{array}$$(31) (32) $$\begin{array}{cc}& {\int }_0^T\{f^V(w,{\mathit{\theta }}^{\prime })+t^H({\mathit{\theta }}_{,t}^{\prime },\mathit{\xi })+a^H({\mathit{\theta }}^{\prime },\mathit{\xi })+b^H(u,{\mathit{\theta }}^{\prime },\mathit{\xi })+f_h^H({\mathit{\theta }}^{\prime },\mathit{\xi })\}\mathrm{d}t\hfill \\ \hfill & -{\int }_{t_1}^{t_2}{E}_{{\mathrm{\Gamma }}_D}({\mathit{\theta }}^{\prime })\mathrm{d}t=0,\forall {\mathit{\theta }}^{\prime }\in \mathrm{\Theta },\hfill \end{array}$$(32) (33) $$\begin{array}{cc}& \mathrm{\Lambda }\ge 0,{\int }_{\mathrm{\Omega }}\mathrm{d}x\le {\mathit{\beta }}_{V}M,\mathrm{\Lambda }\left({\int }_{\mathrm{\Omega }}\mathrm{d}x-{\mathit{\beta }}_{V}M\right)=0.\hfill \end{array}$$(33)

The derivative of the Lagrange function agrees with the derivative of the evaluation function, establishing the following relationship:(34) $$\begin{array}{c}\hfill \dot{L}{{|}_{u_i,p,\mathit{\theta },w_i,q,\mathit{\xi },\mathrm{\Lambda }}={\dot{E}}_{{\mathrm{\Gamma }}_D}|}_{u_i,p,\mathit{\theta },w_i,q,\mathit{\xi },\mathrm{\Lambda }}=l_G(\overrightarrow{V}).\end{array}$$(34)

Since $G\overrightarrow{n}$ in Equation (28(28) $$\begin{array}{cc}& G={\int }_0^T\left\{-\mathit{\mu }{w}_{i,j}{u}_{i,j}+w_{i,i}p-\frac{\mathit{\partial }\mathit{\theta }}{\mathit{\partial }t}\mathit{\xi }-\mathit{\alpha }{\mathit{\xi }}_{,i}{\mathit{\theta }}_{,i}-\mathit{\xi }{u}_j{\mathit{\theta }}_{,j}-{\mathrm{\nabla }}_n(\mathit{\xi }\widehat{q})-(\mathit{\xi }\widehat{q})\mathit{\kappa }\right\}\mathrm{d}t+\mathrm{\Lambda },\hfill \\ \hfill \end{array}$$(28) ) is a coefficient function of the velocity field $\overrightarrow{V}$ that provides minute variations in the domain, $G\overrightarrow{n}$ is referred to as a sensitivity function or shape gradient function. Furthermore, the scalar function G is referred to as the shape gradient density function.

Equation (30(30) $$\begin{array}{cc}& {\int }_0^T\{t^H({\mathit{\theta }}_{,t},{\mathit{\xi }}^{\prime })+a^H(\mathit{\theta },{\mathit{\xi }}^{\prime })+b^H(u,\mathit{\theta },{\mathit{\xi }}^{\prime })+f_q^H({\mathit{\xi }}^{\prime })+f_h^H(\mathit{\theta },{\mathit{\xi }}^{\prime })\hfill \\ \hfill & -f_{hf}^H({\mathit{\xi }}^{\prime })\}\mathrm{d}t=0,\forall {\mathit{\xi }}^{\prime }\in \mathrm{\Xi },\hfill \end{array}$$(30) ) is a weak form of the Navier–Stokes equation and the continuity equation. Equation (31(31) $$\begin{array}{cc}& {\int }_0^T\{t^V(u_{,t}^{\prime },w)+a^V(u^{\prime },w)+b^V(u^{\prime },u,w)+b^V(u,u^{\prime },w)+c(u^{\prime },q)+b^H(u^{\prime },\mathit{\theta },\mathit{\xi })\hfill \\ \hfill & +c(w,p^{\prime })\}\mathrm{d}t=0,\forall u^{\prime }\in U,\forall p^{\prime }\in Q,\hfill \end{array}$$(31) ) is a weak form of the energy equation in the state equation. Equation (32(32) $$\begin{array}{cc}& {\int }_0^T\{f^V(w,{\mathit{\theta }}^{\prime })+t^H({\mathit{\theta }}_{,t}^{\prime },\mathit{\xi })+a^H({\mathit{\theta }}^{\prime },\mathit{\xi })+b^H(u,{\mathit{\theta }}^{\prime },\mathit{\xi })+f_h^H({\mathit{\theta }}^{\prime },\mathit{\xi })\}\mathrm{d}t\hfill \\ \hfill & -{\int }_{t_1}^{t_2}{E}_{{\mathrm{\Gamma }}_D}({\mathit{\theta }}^{\prime })\mathrm{d}t=0,\forall {\mathit{\theta }}^{\prime }\in \mathrm{\Theta },\hfill \end{array}$$(32) ) is a weak form of the Navier–Stokes equation and continuous state equation for the adjoint problem, Equation (33(33) $$\begin{array}{cc}& \mathrm{\Lambda }\ge 0,{\int }_{\mathrm{\Omega }}\mathrm{d}x\le {\mathit{\beta }}_{V}M,\mathrm{\Lambda }\left({\int }_{\mathrm{\Omega }}\mathrm{d}x-{\mathit{\beta }}_{V}M\right)=0.\hfill \end{array}$$(33) ) is a weak form of the energy equation in the state equation for the adjoint problem, and (34(34) $$\begin{array}{c}\hfill \dot{L}{{|}_{u_i,p,\mathit{\theta },w_i,q,\mathit{\xi },\mathrm{\Lambda }}={\dot{E}}_{{\mathrm{\Gamma }}_D}|}_{u_i,p,\mathit{\theta },w_i,q,\mathit{\xi },\mathrm{\Lambda }}=l_G(\overrightarrow{V}).\end{array}$$(34) ) is a constraint equation related to the Lagrange multiplier $\mathrm{\Lambda }$. The traction method can be applied if the shape gradient function can be evaluated by analyzing $u_i,p,\mathit{\theta },w_i,q,\mathit{\xi }$, and $\mathrm{\Lambda }$ based on these equations.

Furthermore, considering the continuity equation for the adjoint equation, and assuming that the flow velocity on the design boundaries satisfies $u_i=0$, the shape gradient density function G can be expressed from (29(29) $$\begin{array}{cc}& {\int }_0^T\{t^V(u_{,t},w^{\prime })+a^V(u,w^{\prime })+b^V(u,u,w^{\prime })+c(w^{\prime },p)+f^V(w^{\prime },\mathit{\theta })-l(w^{\prime })\hfill \\ \hfill & +c(u,q^{\prime })\}\mathrm{d}t=0,\forall w^{\prime }\in W,\forall q^{\prime }\in Q,\hfill \end{array}$$(29) ) as(35) $$\begin{array}{cc}& G=G_0+G_1,\hfill \\ \hfill & G_0={\int }_0^T\left\{-\mathit{\mu }{w}_{i,j}{u}_{i,j}-\frac{\mathit{\partial }\mathit{\theta }}{\mathit{\partial }t}\mathit{\xi }-\mathit{\alpha }{\mathit{\xi }}_{,i}{\mathit{\theta }}_{,i}-{\mathrm{\nabla }}_n(\mathit{\xi }\widehat{q})-(\mathit{\xi }\widehat{q})\mathit{\kappa }\right\}\mathrm{d}t,\hfill \\ \hfill & G_1=\mathrm{\Lambda },\hfill \end{array}$$(35)

where $G_0$ is the shape gradient density function for the square integration errors and $G_1$ is the shape gradient density function for the constraint conditions.

The same procedure can also be used to derive the adjoint equation and shape gradient function when prescribing the temperature ${\mathit{\theta }|}_{{\mathrm{\Omega }}_D\times [t_1,t_2]}$ in sub-domain ${\mathrm{\Omega }}_D$.

The adjoint equations for the shape identification problem can be expressed as follows:(36) $$\begin{array}{cc}& {\int }_0^T\{t^V(u_{,t}^{\prime },w)+a^V(u^{\prime },w)+b^V(u^{\prime },u,w)+b^V(u,u^{\prime },w)+c(u^{\prime },q)+b^H(u^{\prime },\mathit{\theta },\mathit{\xi })\hfill \\ \hfill & +c(w,p^{\prime })\}\mathrm{d}t=0,\forall u^{\prime }\in U,\forall p^{\prime }\in Q,\hfill \end{array}$$(36) (37) $$\begin{array}{cc}& {\int }_0^T\{f^V(w,{\mathit{\theta }}^{\prime })+t^H({\mathit{\theta }}_{,t}^{\prime },\mathit{\xi })+a^H({\mathit{\theta }}^{\prime },\mathit{\xi })+b^H(u,{\mathit{\theta }}^{\prime },\mathit{\xi })+f_h^H({\mathit{\theta }}^{\prime },\mathit{\xi })\}\mathrm{d}t\hfill \\ \hfill & -{\int }_{t_1}^{t_2}{E}_{{\mathrm{\Omega }}_D}({\mathit{\theta }}^{\prime })\mathrm{d}t=0,\forall {\mathit{\theta }}^{\prime }\in \mathrm{\Theta },\hfill \end{array}$$(37) (38) $$\begin{array}{cc}& \mathrm{\Lambda }\ge 0,{\int }_{\mathrm{\Omega }}\mathrm{d}x\le {\mathit{\beta }}_{V}M,\mathrm{\Lambda }\left({\int }_{\mathrm{\Omega }}\mathrm{d}x-{\mathit{\beta }}_{V}M\right)=0.\hfill \end{array}$$(38)

5. Maximizing heat discharge on a sub-boundary

5.1. Problem formulation

The problem of maximizing the discharged heat from time $t=t_1$ to $t=t_2$ on a sub-boundary ${\mathrm{\Gamma }}_D\subset {\mathrm{\Gamma }}_h$ is formulated in the same way as in Section 4.1. For simplicity, we assume that the sub-boundaries ${\mathrm{\Gamma }}_{\mathit{\sigma }}$, ${\mathrm{\Gamma }}_{\mathit{\theta }}$ and ${\mathrm{\Gamma }}_h$ are invariable, that is $\overrightarrow{T_s}({\mathrm{\Gamma }}_{\mathit{\sigma }})={\mathrm{\Gamma }}_{\mathit{\sigma }}$, $\overrightarrow{T_s}({\mathrm{\Gamma }}_{\mathit{\theta }})={\mathrm{\Gamma }}_{\mathit{\theta }}$ and $\overrightarrow{T_s}({\mathrm{\Gamma }}_h)={\mathrm{\Gamma }}_h$ for domain variation. The heat discharge maximization problem from time $t=t_1$ to $t=t_2$ is formulated as follows:(39) $$\begin{array}{c}\hfill \text{Given}\mathrm{\Omega }\text{find}{\mathrm{\Omega }}_s\text{that maximizes}{\int }_{t_1}^{t_2}{J}_{{\mathrm{\Gamma }}_D}(\mathit{\theta })\mathrm{d}t\text{subject to}\text{(15)--(17)}\text{and}\text{(24)},\end{array}$$(39)

where(40) $$\begin{array}{c}\hfill J_{{\mathrm{\Gamma }}_D}(\mathit{\theta })={\int }_{{\mathrm{\Gamma }}_D\subset {\mathrm{\Gamma }}_h}\widehat{h}(\mathit{\theta }-{\widehat{\mathit{\theta }}}_f)\mathrm{d}\mathrm{\Gamma }.\end{array}$$(40)

5.2. Shape gradient function

The Lagrange function $L(u_i,p,\mathit{\theta },w_i,q,\mathit{\xi },\mathrm{\Lambda })$ for this problem is given by the following formula:(41) $$\begin{array}{cc}\hfill L=& -{\int }_{t_1}^{t_2}{J}_{{\mathrm{\Gamma }}_D}(\mathit{\theta })\mathrm{d}t\hfill \\ \hfill & -{\int }_0^T\{t^V(u_{,t},w)+a^V(u,w)+b^V(u,u,w)+c(w,p)+f^V(w,\mathit{\theta })-l(w)\}\mathrm{d}t\hfill \\ \hfill & -{\int }_0^T\{c(u,q)\}\mathrm{d}t\hfill \\ \hfill & -{\int }_0^T\{t^H({\mathit{\theta }}_{,t},\mathit{\xi })+a^H(\mathit{\theta },\mathit{\xi })+b^H(u,\mathit{\theta },\mathit{\xi })+f_q^H(\mathit{\xi })+f_h^H(\mathit{\theta },\mathit{\xi })-f_{hf}^H(\mathit{\xi })\}\mathrm{d}t\hfill \\ \hfill & +\mathrm{\Lambda }\left({\int }_{\mathrm{\Omega }}\mathrm{d}x-{\mathit{\beta }}_{V}M\right).\hfill \end{array}$$(41)

Assuming that the flow velocity satisfies $u_i=0$ at the design boundary, the shape gradient density function for the problem is given by (36(36) $$\begin{array}{cc}& {\int }_0^T\{t^V(u_{,t}^{\prime },w)+a^V(u^{\prime },w)+b^V(u^{\prime },u,w)+b^V(u,u^{\prime },w)+c(u^{\prime },q)+b^H(u^{\prime },\mathit{\theta },\mathit{\xi })\hfill \\ \hfill & +c(w,p^{\prime })\}\mathrm{d}t=0,\forall u^{\prime }\in U,\forall p^{\prime }\in Q,\hfill \end{array}$$(36) ).

The adjoint flow velocity $w_i$, adjoint pressure q, adjoint temperature $\mathit{\xi }$, and Lagrange multiplier $\mathrm{\Lambda }$ related to the heat discharge maximization problem are determined by the following adjoint equations:(42) $$\begin{array}{cc}& {\int }_0^T\{t^V(u_{,t}^{\prime },w)+a^V(u^{\prime },w)+b^V(u^{\prime },u,w)+b^V(u,u^{\prime },w)+c(u^{\prime },q)\hfill \\ \hfill & +b^H(u^{\prime },\mathit{\theta },\mathit{\xi })+c(w,p^{\prime })\}\mathrm{d}t=0,\forall u^{\prime }\in U,\forall p^{\prime }\in Q,\hfill \end{array}$$(42) (43) $$\begin{array}{cc}& {\int }_0^T\{f^V(w,{\mathit{\theta }}^{\prime })+t^H({\mathit{\theta }}_{,t}^{\prime },\mathit{\xi })+a^H({\mathit{\theta }}^{\prime },\mathit{\xi })+b^H(u,{\mathit{\theta }}^{\prime },\mathit{\xi })+f_h^H({\mathit{\theta }}^{\prime },\mathit{\xi })\}\mathrm{d}t\hfill \\ \hfill & +{\int }_{t_1}^{t_2}{J}_{{\mathrm{\Gamma }}_D}({\mathit{\theta }}^{\prime })\mathrm{d}t=0,\forall {\mathit{\theta }}^{\prime }\in \mathrm{\Theta },\hfill \\ \hfill & \text{and}(39).\hfill \end{array}$$(43)

6. Solution

6.1. Traction method

The traction method has been proposed for solving the velocity field $\overrightarrow{V}\in D$ that indicates domain variation. Domain variation is based on the following governing equation [7]:(44) $$\begin{array}{c}\hfill a^E(\overrightarrow{V},y)=-<G\overrightarrow{n},y>,\forall y\in D,\end{array}$$(44)

and(45) $$\begin{array}{c}\hfill D=\{\overrightarrow{V}\in {(H^1(\mathit{\Omega }))}^d|\text{constraints of domain variation}\}.\end{array}$$(45)

Here, $a^E(\overrightarrow{V},y)$ is a bilinear form indicating the strain energy on a linear elastic body, which is defined by the following equation with respect to the displacement vector function $u_i$, $v_i$:(46) $$\begin{array}{c}\hfill a^E(u,v)={\int }_{\mathrm{\Omega }}{A}_{ijkl}{u}_{k,l}{v}_{i,j}\mathrm{d}x,\end{array}$$(46)

where $A_{ijkl}$ is a constant elastic tensor that is positive definite. Equation (44(44) $$\begin{array}{c}\hfill a^E(\overrightarrow{V},y)=-<G\overrightarrow{n},y>,\forall y\in D,\end{array}$$(44) ) indicates that the velocity field $\overrightarrow{V}$ is analyzed as a displacement field when the negative shape gradient function $-G{n}_i$ is applied as an external force. The domain variation resulting from the traction method can be derived as a displacement field when the shape gradient function is applied to the pseudo-elasticity problem as an external force.

Furthermore, this analysis is a heat convection field analysis. For simplicity, the elastic tensor $A_{ijkl}$ is given by(47) $$\begin{array}{c}\hfill A_{ijkl}={\mathit{\delta }}_{ik}{\mathit{\delta }}_{jl}.\end{array}$$(47)

Taking account of a bilinear Equation (44(44) $$\begin{array}{c}\hfill a^E(\overrightarrow{V},y)=-<G\overrightarrow{n},y>,\forall y\in D,\end{array}$$(44) ) and $G=G_0+G_1$, the domain variation is expressed as $\overrightarrow{V}={\overrightarrow{V}}_0+\mathrm{\Lambda }{\overrightarrow{V}}_1$. We obtain(48) $$\begin{array}{cc}& a^E({\overrightarrow{V}}_0,y)=-<G_0\overrightarrow{n},y>,\forall y\in D.\hfill \\ \hfill & a^E(\mathrm{\Lambda }{\overrightarrow{V}}_1,y)=-<G_1\overrightarrow{n},y>,\forall y\in D.\hfill \end{array}$$(48)

The Lagrange multiplier $\mathrm{\Lambda }$ is determined from(49) $$\begin{array}{c}\hfill {\int }_{\mathrm{\Gamma }}\overrightarrow{n}·{\overrightarrow{V}}_{0}d\mathrm{\Gamma }+\mathrm{\Lambda }{\int }_{\mathrm{\Gamma }}\overrightarrow{n}·{\overrightarrow{V}}_{1}d\mathrm{\Gamma }=0.\end{array}$$(49)

Finally, because ${\overrightarrow{V}}_0,{\overrightarrow{V}}_1,\mathrm{\Lambda }$ are given, the domain $\overrightarrow{x}\in \mathrm{\Omega }$ is moved into $\overrightarrow{x}+\overrightarrow{V}\in {\mathrm{\Omega }}_s$.

6.2. Analysis procedure

The proposed shape determination analysis can be performed by repeating the following steps:

Step 1.	Give initial shape.
Step 2.	Analyze the flow velocity $u(\overrightarrow{x},t)$, pressure $p(\overrightarrow{x},t)$ and temperature $\mathit{\theta }(\overrightarrow{x},t)$, based on the state Equation (30(30) $$\begin{array}{cc}& {\int }_0^T\{t^H({\mathit{\theta }}_{,t},{\mathit{\xi }}^{\prime })+a^H(\mathit{\theta },{\mathit{\xi }}^{\prime })+b^H(u,\mathit{\theta },{\mathit{\xi }}^{\prime })+f_q^H({\mathit{\xi }}^{\prime })+f_h^H(\mathit{\theta },{\mathit{\xi }}^{\prime })\hfill \\ \hfill & -f_{hf}^H({\mathit{\xi }}^{\prime })\}\mathrm{d}t=0,\forall {\mathit{\xi }}^{\prime }\in \mathrm{\Xi },\hfill \end{array}$$(30) ) and (31(31) $$\begin{array}{cc}& {\int }_0^T\{t^V(u_{,t}^{\prime },w)+a^V(u^{\prime },w)+b^V(u^{\prime },u,w)+b^V(u,u^{\prime },w)+c(u^{\prime },q)+b^H(u^{\prime },\mathit{\theta },\mathit{\xi })\hfill \\ \hfill & +c(w,p^{\prime })\}\mathrm{d}t=0,\forall u^{\prime }\in U,\forall p^{\prime }\in Q,\hfill \end{array}$$(31) ) using initial condition.
Step 3.	Stop when the objective functional has converged.
Step 4.	Analyze the adjoint flow velocity $w(\overrightarrow{x},t)$, adjoint pressure $q(\overrightarrow{x},t)$ and adjoint temperature $\mathit{\xi }(\overrightarrow{x},t)$, using the given flow velocity $u(\overrightarrow{x},t)$ and temperature $\mathit{\theta }(\overrightarrow{x},t)$, based on the adjoint Equations (32(32) $$\begin{array}{cc}& {\int }_0^T\{f^V(w,{\mathit{\theta }}^{\prime })+t^H({\mathit{\theta }}_{,t}^{\prime },\mathit{\xi })+a^H({\mathit{\theta }}^{\prime },\mathit{\xi })+b^H(u,{\mathit{\theta }}^{\prime },\mathit{\xi })+f_h^H({\mathit{\theta }}^{\prime },\mathit{\xi })\}\mathrm{d}t\hfill \\ \hfill & -{\int }_{t_1}^{t_2}{E}_{{\mathrm{\Gamma }}_D}({\mathit{\theta }}^{\prime })\mathrm{d}t=0,\forall {\mathit{\theta }}^{\prime }\in \mathrm{\Theta },\hfill \end{array}$$(32) ) and (33(33) $$\begin{array}{cc}& \mathrm{\Lambda }\ge 0,{\int }_{\mathrm{\Omega }}\mathrm{d}x\le {\mathit{\beta }}_{V}M,\mathrm{\Lambda }\left({\int }_{\mathrm{\Omega }}\mathrm{d}x-{\mathit{\beta }}_{V}M\right)=0.\hfill \end{array}$$(33) ) or Equations (42(42) $$\begin{array}{cc}& {\int }_0^T\{t^V(u_{,t}^{\prime },w)+a^V(u^{\prime },w)+b^V(u^{\prime },u,w)+b^V(u,u^{\prime },w)+c(u^{\prime },q)\hfill \\ \hfill & +b^H(u^{\prime },\mathit{\theta },\mathit{\xi })+c(w,p^{\prime })\}\mathrm{d}t=0,\forall u^{\prime }\in U,\forall p^{\prime }\in Q,\hfill \end{array}$$(42) ) and (43(43) $$\begin{array}{cc}& {\int }_0^T\{f^V(w,{\mathit{\theta }}^{\prime })+t^H({\mathit{\theta }}_{,t}^{\prime },\mathit{\xi })+a^H({\mathit{\theta }}^{\prime },\mathit{\xi })+b^H(u,{\mathit{\theta }}^{\prime },\mathit{\xi })+f_h^H({\mathit{\theta }}^{\prime },\mathit{\xi })\}\mathrm{d}t\hfill \\ \hfill & +{\int }_{t_1}^{t_2}{J}_{{\mathrm{\Gamma }}_D}({\mathit{\theta }}^{\prime })\mathrm{d}t=0,\forall {\mathit{\theta }}^{\prime }\in \mathrm{\Theta },\hfill \\ \hfill & \text{and}(39).\hfill \end{array}$$(43) ), using the conditions $w_i(\overrightarrow{x},T)=0$ and $\mathit{\xi }(\overrightarrow{x},T)=0$ in (22(22) $$\begin{array}{cc}\hfill \mathrm{\Xi }& =\{\mathit{\xi }(\overrightarrow{x},t)\in H^1(\mathrm{\Omega }\times [0,T])|\mathit{\xi }(\overrightarrow{x},t)=0,t\in [0,T],\overrightarrow{x}\in {\mathrm{\Gamma }}_{\mathit{\theta }},\hfill \\ \hfill & \mathit{\xi }(\overrightarrow{x},T)=0,\overrightarrow{x}\in \mathrm{\Omega }\}.\hfill \end{array}$$(22) ) and (23(23) $$\begin{array}{cc}& \text{Given}\mathrm{\Omega }\text{find}{\mathrm{\Omega }}_s\text{that minimizes}{\int }_{t_1}^{t_2}{E}_{{\mathrm{\Gamma }}_D}(\mathit{\theta })\mathrm{d}t\text{subject to}\text{(15)--(17)}\text{and}\hfill \\ \hfill & {\int }_{\mathrm{\Omega }}\mathrm{d}x\le {\mathit{\beta }}_{V}M,\hfill \end{array}$$(23) ).
Step 5.	Calculate the shape gradient density function G based on these results using Equation (36(36) $$\begin{array}{cc}& {\int }_0^T\{t^V(u_{,t}^{\prime },w)+a^V(u^{\prime },w)+b^V(u^{\prime },u,w)+b^V(u,u^{\prime },w)+c(u^{\prime },q)+b^H(u^{\prime },\mathit{\theta },\mathit{\xi })\hfill \\ \hfill & +c(w,p^{\prime })\}\mathrm{d}t=0,\forall u^{\prime }\in U,\forall p^{\prime }\in Q,\hfill \end{array}$$(36) ).
Step 6.	Calculate the velocity field $\overrightarrow{V}$ based on the traction method using Equation (44(44) $$\begin{array}{c}\hfill a^E(\overrightarrow{V},y)=-<G\overrightarrow{n},y>,\forall y\in D,\end{array}$$(44) ) to update the shape. Return to Step 2.

7. Numerical experiments

To verify the proposed method, numerical programs for the above two shape determination problems were developed based on FreeFem++. This section presents simple numerical experiments analyzed using the derived shape gradient function and the traction method to solve the problem of prescribing temperature history at a sub-boundary and the problem of maximizing heat discharge.

7.1. Prescribing temperature on a sub-boundary

An analysis model was considered similarly to the model presented by Park et al. [13] and shown in Figure 2. In Figure 2, the target shape is represented by dotted lines, while the initial shape is represented by solid lines. The boundary conditions of the temperature field were assumed as follows. The lower surface boundary AB was assumed to be the high temperature base boundary ${\mathrm{\Gamma }}_{\mathit{\theta }}$ and was set to $\widehat{\mathit{\theta }}=1$ K. Surfaces DA and BC on the left and right, respectively, were assumed to be low temperature base boundaries ${\mathrm{\Gamma }}_{\mathit{\theta }}$ and were set to $\widehat{\mathit{\theta }}$=0 K. Surface DC at ${\mathrm{\Gamma }}_q$ was assumed to be a natural boundary, and the thermal insulation condition was set to $\widehat{q}=0$. The temperature for the target shape on surface DC was assumed to be the target temperature ${\mathit{\theta }}_D(\overrightarrow{x},t)$ in Equation (25(25) $$\begin{array}{cc}\hfill L=& {\int }_{t_1}^{t_2}{E}_{{\mathrm{\Gamma }}_D}(\mathit{\theta })\mathrm{d}t\hfill \\ \hfill & -{\int }_0^T\{t^V(u_{,t},w)+a^V(u,w)+b^V(u,u,w)+c(w,p)+f^V(w,\mathit{\theta })-l(w)\}\mathrm{d}t\hfill \\ \hfill & -{\int }_0^T\{c(u,q)\}\mathrm{d}t\hfill \\ \hfill & -{\int }_0^T\{t^H({\mathit{\theta }}_{,t},\mathit{\xi })+a^H(\mathit{\theta },\mathit{\xi })+b^H(u,\mathit{\theta },\mathit{\xi })+f_q^H(\mathit{\xi })+f_h^H(\mathit{\theta },\mathit{\xi })-f_{hf}^H(\mathit{\xi })\}\mathrm{d}t\hfill \\ \hfill & +\mathrm{\Lambda }\left({\int }_{\mathrm{\Omega }}\mathrm{d}x-{\mathit{\beta }}_{V}M\right).\hfill \end{array}$$(25) ). The boundary conditions for the flow field assumed that all surrounding boundaries were wall boundaries set to $u_1=u_2=0$.

Figure 2. Numerical model and boundary conditions for Park’s problem.

The initial conditions were set to ${\mathit{\theta }}_{}\text{ini}=0$ and $u_{i_{}\text{ini}}=0$ for the entire domain. Internally generated heat was ignored, and the pressure was considered to be unique in order to achieve an average value of 0. The time was set to $t_1=0$ and $t_2=T$, while time integration was performed from $t=0$ s to $t=T=600$ s with a time increment of $\mathrm{\Delta }t=0.8$ s.

The prescribed temperature on boundary ${\mathrm{\Gamma }}_D$ was considered for the upper surface DC. The design boundary ${\mathrm{\Gamma }}_{}\text{design}$ was considered to be the high temperature base boundary of lower surface AB. The remaining boundaries were set using constraints on the shape variation, specifically low temperature base boundaries DA and BC and thermal insulation boundary DC.

The temperature field $\mathit{\theta }$, flow field velocity $\overrightarrow{u}$ and pressure p, adjoint temperature $\mathit{\xi }$, adjoint flow velocity $\overrightarrow{w}$ and adjoint pressure q, and shape updating analysis (velocity field $\overrightarrow{V}$) were all performed using FreeFem++ [18,19]. A triangular secondary element was used for $\mathit{\theta },\overrightarrow{u},\mathit{\xi },\overrightarrow{w}$, and $\overrightarrow{V}$, while a triangular primary element was used for p, q.

Based on the FEM, an analysis was performed using the function ‘convect’ of FreeFem++ on the state equation in the time direction. The function ‘adaptmesh’ of FreeFem++ was used to perform remeshing for shape updates; ‘adaptmesh’ is a function that generates an adaptive mesh for improving solution accuracy of the state variables derived by finite-element analysis. Adaptmesh was used to improve the accuracy of the shape gradient density function G and velocity field $\overrightarrow{V}$. The following physical values were used based on the book [18]: density $\mathit{\rho }$ = 1.0 kg/m${}^3$, dynamic viscosity $\mathit{\mu }=1.0\times {10}^{-3}$ Pa s, Prandtl number $Pr=7.01$, thermal diffusivity $\mathit{\alpha }=\mathit{\mu }/(\mathit{\rho }Pr)$, gravitational acceleration $g=9.8$ m/s${}^2$ downward, base temperature ${\mathit{\theta }}_0=0$ K, and coefficient of thermal expansion $\mathit{\beta }=2.06\times {10}^{-4}$ 1/K. The dimensions of the shape were set to $h=0.5$ m, $b=2$ m and $a=0.125$ m.

The value of the objective functional at iteration k was denoted by $J^{[k]}$, and convergence was reached when(50) $$\begin{array}{c}\hfill \frac{|J^{[k-1]}-J^{[k]}|}{J^{[k]}}<\mathit{\epsilon },\end{array}$$(50)

where $\mathit{\epsilon }>0$ is a small number.

For $\mathit{\epsilon }={10}^{-5}$, a comparison of the initial, identified and target shapes is shown in Figure 3. The temperature $\mathit{\theta }$, flow velocity $\overrightarrow{u}$, and streamlines at the end time $t=T=600$ s for the respective shapes are shown in Figure 4. From Figure 4, it was confirmed that the temperature distribution, flow velocity distributions, and streamlines in the identified shape agreed with those in the target shape.

Figure 3. Shapes with finite-element meshes for Park’s problem.

Figure 4. Temperature distributions, flow velocity distributions, and streamlines at time $t=T=600$ for Park’s problem.

A comparison of the target shape (dotted lines) and the outline of the obtained shape (solid lines) is shown in Figure 5. The target and identifying shapes are practically identical. Figure 6 shows the iterative history for the objective functional, normalized with initial value.

Figure 5. Comparison of target shape and identified shape for Park problem.

Figure 6. Numerical results: Iterative history for Park’s problem.

Figure 7. Numerical model and boundary condition for heat discharge maximization problem.

Figure 8. Mesh and distribution of temperature and streamlines at time $t=T=$ 600 for initial shape for heat discharge maximization problem.

Figure 9. Numerical results: Mesh and distributions of temperature and streamlines at time $t=T=$ 600 for optimum shape for heat discharge maximization problem.

Figure 10. Iterative histories for heat discharge maximization problem.

These results confirm that the objective functional, which is reduced to practically zero, and the identifying shapes, which are practically identical to the target shapes, can be obtained using the proposed method.

7.2. Maximizing heat discharge on a sub-boundary

A disc with five holes, as depicted in Figure 7, was used as an analysis model for the problem of the maximizing heat discharge for an unsteady natural convection field.

The boundary conditions considered hole A on the lower side to be the known temperature boundary ${\mathrm{\Gamma }}_{\mathit{\theta }}$ and was set to $\widehat{\mathit{\theta }}=1$ K. Holes B, C, D and E were considered to be thermal insulation boundaries ${\mathrm{\Gamma }}_q$ and were set to $\widehat{q}=0$ mK/s. The external boundaries considered $\widehat{h}=1$ m/s and ${\widehat{\mathit{\theta }}}_f=0$ K as heat transfer boundaries ${\mathrm{\Gamma }}_h$. The boundary conditions for the flow field assumed that all surrounding boundaries were wall boundaries set to $u_1=u_2$= 0. The initial conditions of the entire domain were set to ${\mathit{\theta }}_{}\text{ini}=0$ and $u_{i_{}\text{ini}}=0$. Internally generated heat was ignored, and the pressure was uniquely set to achieve an average of 0. The time was set to $t_1=0$ and $t_2=T$, and time integration was performed from $t=0$ s to $t=T=600$ s with $\mathrm{\Delta }t=0.8$ s time increment. The four hole boundaries B, C, D, and E comprising thermal insulation boundary ${\mathrm{\Gamma }}_q$ were considered to be design boundaries ${\mathrm{\Gamma }}_{}\text{design}$. Other boundaries were constrained with respect to domain variation in order to analyze the problem of maximizing the heat discharge at the external heat transfer boundary ${\mathrm{\Gamma }}_h$.

The physical values were the same as those used in the example of Section 7.1, while the dimensions were set to $d_0=1.2$ m, $d_i=0.3$ m, and $a=0.7$ m. The coefficient ${\mathit{\beta }}_V$, which restrains the size of the domain, was varied for three values: in Case A, ${\mathit{\beta }}_V=1$; in Case B, ${\mathit{\beta }}_V=0.9$; and in Case C, ${\mathit{\beta }}_V=1.1$. For $\mathit{\epsilon }={10}^{-5}$, the obtained shape for optimization resulted in a sharpened shape along the horizontal direction for hole boundaries B and E, even when the FreeFem++ adaptmesh function was used after each shape update. It was therefore difficult to converge to the heat discharge value from the given initial guess. As a result, the value of $\mathit{\epsilon }$ was set to $\mathit{\epsilon }={10}^{-3}$.

The mesh, temperature, and streamlines of the initial shape at end time $t=T=600$ s are shown in Figure 8, with the corresponding numerical results shown for Cases A, B, and C in Figure 9. The heat discharge, which is the objective functional being sought by shape updates, and domain size changes are shown in Figure 10. In Figure 10, the objective functional and the volume are normalized with initial value. The shape of hole boundaries B and E was updated in the horizontal direction in all cases based on the results shown in Figure 9. The temperature was elevated in the heat transfer boundary near hole boundaries B and E. As a result, the obtained shape had increasing heat discharge on the external heat transfer boundary ${\mathrm{\Gamma }}_h$. A review of the Case B results reveals that the shape was modified for vertical extension on hole boundaries C and D, where the influence of the heat discharge is relatively low. This resulted in a shape that satisfies the restraining condition ${\mathit{\beta }}_V=0.9$ to reduce the size of the domain. Similarly, the shape of boundaries C and D was vertically modified in Case C as well, resulting in a shape that had an increased domain size. The results shown in Figure 10 confirm that, in all cases, the heat discharge increased while satisfying the constraining conditions of the domain size. An improvement of about 5% occurred with respect to the initial shape.

8. Conclusion

This paper proposed a numerical method for two problems: the problem of shape identification for prescribing the temperature history for a sub-boundary in unsteady natural convection fields, and the problem of shape optimization for maximizing heat discharge on a sub-boundary. These problems were formulated, and the shape gradient density functions were derived. The verification of the proposed method was accomplished using simple examples of two-dimensional problems by applying the traction method to the shape gradient distribution functions.

In addition, in this study, we did not consider noise in the input measurement and a corresponding stopping criterion for iterations. Because these are important matters about shape determination of inverse problems, we are going to examine them as future research themes.

Notes

No potential conflict of interest was reported by the authors.