A Finite Element-Finite Volume Code Coupling for Optimal Control Problems in Fluid Heat Transfer

1. Introduction

One of the most important aspects in the design of an engineering system is foreseeing the physical behavior of the variables considered. Knowing the range in which a specific quantity can evolve is a crucial issue in checking the correct operation of a given apparatus. In this context, the numerical simulation of the partial differential equations that describe physical phenomena turns out to be essential. Specifically, the standard approach for a numerical simulation consists of building a parameter set, such as boundary conditions, source terms, and material physical properties that, given as an input to the conservation equations, provides a numerical output that describes the evolution of the physical system. On the other hand, this approach leads to a try-and-fail approach since we do not have direct control over the output results, and thus it can become expensive from a computational point of view. However, numerical optimization is proposed as a viable alternative to the former approach because the desired outputs turn into input parameters that can be selected in the model. In fact, the idea is to minimize a cost function that represents the distance between a variable and its desired target value under specific physical constraints that are the governing equations [1]. We can therefore introduce the objective term that represents a physical quantity to be obtained, such as a velocity or temperature matching, but also the enhancement or the reduction of a specific variable. The control term is a physical parameter that is used to perform the minimization, and it can be distributed over the domain, or it can act only on the boundaries. Possibly, the control can modify the geometry of the problem in the context of shape optimization [2]. Lastly, objective and control terms must follow the governing equations that have the role of constraints describing the physics under investigation.

In this work, we present a numerical optimization algorithm that considers generic heat transfer problems. In particular, we present two examples that are representative of physical phenomena commonly present in several engineering applications. The first one involves the solution of the Navier-Stokes and energy equations (under the Boussinesq approximation) in the context of a natural convection problem considering the Rayleigh-Bénard setting as described in [3]. The second one is described by the Conjugate Heat Transfer (CHT) model, which describes the heat transfer between a solid and a fluid region as seen in [4]. Specifically, we adopt an adjoint-based optimization, which entails the solution of a coupled system given by the state, adjoint, and control equations [5].

Regarding the first application, several results have already been shown in the literature. One of the first works appears in [6], where a natural convection optimal control problem has been considered for a two-dimensional square cavity and for a temperature matching problem. In [7], an optimal control problem for the Boussinesq system has been presented in the context of liquid metals in a differentially heated cavity, while in [8] the same approach is described with the addition of a dynamic turbulence model. Several papers have shown the existence, uniqueness, and regularity of the solution for this kind of problem. The interested reader can consult [9,10,11,12,13,14,15,16,17,18,19,20], where different settings are taken into account, including the use of different control approaches. For the specific problem of the boundary control problem considering the stationary Rayleigh-Bénard-Marangoni setting, one can consult [21]. Regarding the optimal control of the heat transfer between a fluid and solid region, or more in general heat transfer problems, we refer to [4,22,23].

Nevertheless, validated and state-of-the-art codes already exist nowadays, at least for the state equations (i.e. the constraints of the optimal control system). Therefore, the aim of this paper is to leverage these codes to build the optimal control system. Here we extend the work presented in [24,25] with the aim of coupling two CFD codes for the implementation of a heat transfer optimal control problem. The idea is to use the finite volume (FV) code OpenFOAM [26] for the solution of the state equation and the in-house finite element (FE) code FEMuS [27] for the adjoint and control equations. In fact, the latter equations arise naturally from a variational formulation of the considered problem, justifying the use of a finite element approach. In this framework, the latter equations should be solved because their derivation requires the concept of variation of the unknowns, with the possibility of also including the gradient of the variations. It appears clear from a mathematical point of view that the finite element approach lends itself to these representations which are often missing or cumbersome to evaluate in standard finite volume frameworks. The code coupling, i.e. the data transfer between the different codes, is implemented by using the external library MEDCOUPLING [28].

Additionally, two main approaches can be followed when dealing with optimal control problems. The first one consists of discretizing the optimality systems after minimizing the functional considering an infinite-dimensional framework (Optimize then discretize). Differently, the second one discretizes the Lagrangian equation, and then it derives the different equations to be solved (Discretize then optimize). Under generic assumptions, these two approaches do not produce the same set of equations [1,5]. The numerical code coupling by using finite volume and finite element discretization can implement only the first approach, since with the second one, the discretized equations, such as the adjoint and the control, follow directly from the discretization used for the Lagrangian function.

The paper is organized as follows: in the first section, we derive the optimality systems for the heat transfer problems, one considering the Boussinesq equations, and a second one for the Conjugate Heat Transfer between solid and fluid. After that, we provide a brief description of the numerical algorithm for the code coupling involving the two CFD codes via the external library. Lastly, we provide two numerical examples in order to validate the proposed numerical approach considering numerical benchmarks taken from the literature.

2. Optimal Control Problems

In this section, we present and describe the mathematical formulation of the system of equations considered for both numerical problems. We start by considering the state equations that indeed describe the different methods to transfer heat, i.e., the Boussinesq system and the Conjugate Heat Transfer. After that, since the examples involve different objective and control terms, we derive each optimality system considering the specific application setup.

We first introduce some mathematical notation from functional analysis [29]. Considering a generic domain $\mathrm{\Omega }\in {\mathbb{R}}^N$, $H^s\left(\mathrm{\Omega }\right)$ denotes the standard Sobolev space $W^{s,2}\left(\mathrm{\Omega }\right)$ of order $s\in \mathbb{R}$, and clearly the space of square-integrable functions $L^2\left(\mathrm{\Omega }\right)$ refers to the zero-order Sobolev space $H^0\left(\mathrm{\Omega }\right)$,

$$L^2\left(\mathrm{\Omega }\right)=\left\{f:\mathrm{\Omega }\to \mathbb{R}\left|{\int }_{\mathrm{\Omega }}{\left|f\right|}^{2}d\mathrm{\Omega }<+\infty \right.\right\}.$$

Of particular interest is the space $H^1\left(\mathrm{\Omega }\right)$ defined as

$$H^1\left(\mathrm{\Omega }\right)=\left\{f\in L^2\left(\mathrm{\Omega }\right)\left|{\displaystyle \frac{\partial f}{\partial x_k}}\in L^2\left(\mathrm{\Omega }\right),k=1,2,\dots ,N\right.\right\},$$

and we also introduce the subspace of functions that vanish on the boundary $\mathrm{\Gamma }\in \partial \mathrm{\Omega }$,

$$H_{\mathrm{\Gamma }}^1\left(\mathrm{\Omega }\right)=\left\{f\in H^1\left(\mathrm{\Omega }\right)\left|f=0\mathrm{on}\mathrm{\Gamma }\right.\right\}.$$

For a generic element f of the Sobolev spaces, the norms can be written as

$${\left|\right|f\left|\right|}_{L^2\left(\mathrm{\Omega }\right)}={\int }_{\mathrm{\Omega }}{\left|f\right|}^{2}d\mathrm{\Omega },$$

and

$${\left|\right|f\left|\right|}_{H^1\left(\mathrm{\Omega }\right)}={\int }_{\mathrm{\Omega }}{\left(\right|f|}^2+{|\nabla f|}^2)d\mathrm{\Omega },$$

where ${\displaystyle {\left|\nabla f\right|}^2=\sum _{k=1}^N{\left|\frac{\partial f}{\partial x_k}\right|}^2}$. For functions defined only on the boundary of the domain $\mathrm{\Gamma }$, then its norm is defined as

$${\left|\right|f\left|\right|}_{H^1\left(\mathrm{\Gamma }\right)}={\int }_{\mathrm{\Gamma }}{\left(\right|f|}^2+|{\nabla }_{s}f{|}^2)d\mathrm{\Gamma },$$

where the surface gradient is ${\nabla }_{s}f=\nabla f-(\nabla f·\mathbf{n})\mathbf{n}$. Lastly, for $(f,g)\in L^2$ such that $\left(fg\right)\in L^1\left(\mathrm{\Omega }\right)$, where $L^1$ is the Sobolev space of integrable functions $W^{0,1}$, their inner product is defined as

$$(f,g)={\int }_{\mathrm{\Omega }}fgd\mathrm{\Omega }.$$

In the case of vector-valued spaces, we indicate them with boldface letters ${\mathbf{L}}^2,{\mathbf{H}}^1$. Their norm is defined in a completely analogous manner to the scalar case. Taking $\mathbf{u},\mathbf{v}\in {\mathbf{L}}^2\left(\mathrm{\Omega }\right)$, the inner product is defined as

$$(\mathbf{u},\mathbf{v})={\int }_{\mathrm{\Omega }}\mathbf{u}·\mathbf{v}d\mathrm{\Omega }.$$

Next, we introduce the bilinear forms

$$\begin{array}{cccc}\hfill a(f,g)& ={\int }_{\mathrm{\Omega }}\nabla f·\nabla gd\mathrm{\Omega }\hfill & \hfill & \forall f,g\in H^1\left(\mathrm{\Omega }\right),\hfill \\ \hfill a(\mathbf{u},\mathbf{v})& ={\displaystyle \frac{1}{2}}{\int }_{\mathrm{\Omega }}\left(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T\right):\left(\nabla \mathbf{v}+\nabla {\mathbf{v}}^T\right)d\mathrm{\Omega }\hfill & \hfill & \forall \mathbf{u},\mathbf{v}\in {\mathbf{H}}^1\left(\mathrm{\Omega }\right),\hfill \\ \hfill b(\mathbf{u},f)& =-{\int }_{\mathrm{\Omega }}f\nabla ·\mathbf{u}d\mathrm{\Omega }\hfill & \hfill & \forall \mathbf{u}\in {\mathbf{H}}^1\left(\mathrm{\Omega }\right),\forall f\in L^2\left(\mathrm{\Omega }\right),\hfill \end{array}$$

and the trilinear forms

$$\begin{array}{cccc}\hfill c(\mathbf{u},\mathbf{v},\mathbf{w})& ={\int }_{\mathrm{\Omega }}\left(\left(\mathbf{u}·\nabla \right)\mathbf{v}\right)·\mathbf{w}d\mathrm{\Omega }\hfill & \hfill & \forall \mathbf{u},\mathbf{v},\mathbf{w}\in {\mathbf{H}}^1\left(\mathrm{\Omega }\right),\hfill \\ \hfill c(\mathbf{u},f,g)& ={\int }_{\mathrm{\Omega }}\left(\mathbf{u}·\nabla f\right)gd\mathrm{\Omega }\hfill & \hfill & \forall \mathbf{u}\in {\mathbf{H}}^1\left(\mathrm{\Omega }\right),\forall f\in H^1\left(\mathrm{\Omega }\right),\forall g\in L^2\left(\mathrm{\Omega }\right),\hfill \end{array}$$

that will be useful in the description of the optimality systems.

We introduce the functional $\mathcal{F}$ that we define as

$$\mathcal{F}\left(\mathbf{u},T,q\right)=\mathcal{J}\left(\mathbf{u},T\right)+{\left|\right|q\left|\right|}_{H^s},$$

where $\mathcal{J}$ represents the generic objective functional, i.e. $\left|\right|\nabla \times {\mathbf{u}\left|\right|}_{L^2\left(\mathrm{\Omega }\right)}$ or $\left|\right|T-T_d{\left|\right|}_{L^2\left(\mathrm{\Omega }\right)}$, and ${\left|\right|q\left|\right|}_{H^s}$ the regularization term with $s=0$ for distributed control and $s=1$ or $s=\pm {\displaystyle \frac{1}{2}}$ for boundary control. In this paper, we will address Dirichlet boundary control problems with $H^1\left(\mathrm{\Gamma }\right)$ regularization while Neumann boundary control problems with $L^2\left(\mathrm{\Gamma }\right)$ regularization. In both cases we are restricting the space where we are looking for the solution with respect to its native space, respectively $H^{1/2}$ and $H^{-1/2}$. For more details regarding this aspect and the different numerical approaches to address this issue, the interested reader can consult [30,31,32,33].

2.1. Rayleigh-Bénard Problem

In this section, we derive the optimality system for the Rayleigh-Bénard problem considering a square domain as reported in Figure 1. The Rayleigh-Bénard system exhibits buoyancy-driven flow, where the natural convection develops due to the temperature gradient imposed between the two horizontal walls ${\mathrm{\Gamma }}_{top}$ and ${\mathrm{\Gamma }}_1$. The objective of this optimization case is to stabilize the flow inside the domain, thereby minimizing the vorticity. This will be achieved by introducing appropriate cavity heating, initially using a volumetric heat source within the domain, and then applying a heat flux condition on the boundary.

The coupled Navier-Stokes system and energy equation form the Boussinesq system of equations that model the Rayleigh-Bénard problem. We write the stationary form of the Navier-Stokes equations as

$$\left\{\begin{array}{cc}\nabla ·\mathbf{u}=0\hfill & \mathrm{on}\mathrm{\Omega },\hfill \\ -2\nu \nabla ·\mathcal{D}\left(\mathbf{u}\right)+(\mathbf{u}·\nabla )\mathbf{u}=-\nabla p+\beta \mathbf{g}(T-T_{ref})\hfill & \mathrm{on}\mathrm{\Omega },\hfill \\ \mathbf{u}=0\hfill & \mathrm{on}\partial \mathrm{\Omega },\hfill \end{array}\right.$$

(1)

where $\mathbf{u}$ is the velocity field, p is the pressure field, $\nu$ is the kinematic viscosity, $\beta$ is the thermal expansion coefficient, $\mathbf{g}$ is a vector representing the gravitational acceleration in the negative y direction, and $\mathcal{D}\left(\mathbf{u}\right)={\displaystyle \frac{1}{2}}(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T)$ is the rate of deformation tensor of the flow. Next, we write the stationary energy equation in the temperature formulation as

$$\left\{\begin{array}{cc}-\alpha \Delta T+\mathbf{u}·\nabla T=Q\hfill & \mathrm{on}\mathrm{\Omega },\hfill \\ T=g_D\hfill & \mathrm{on}{\mathrm{\Gamma }}_{top}\cup {\mathrm{\Gamma }}_1,\hfill \\ {\displaystyle \frac{\partial T}{\partial \mathbf{n}}}=g_N\hfill & \mathrm{on}{\mathrm{\Gamma }}_N,\hfill \\ {\displaystyle \frac{\partial T}{\partial \mathbf{n}}}=0\hfill & \mathrm{on}{\mathrm{\Gamma }}_l\cup {\mathrm{\Gamma }}_r,\hfill \end{array}\right.$$

(2)

where T is the temperature field, Q is a generic volumetric heat source, $g_N$ is a generic wall heat flux, $g_D$ is a given temperature profile ad the boundary and $\alpha$ is the thermal diffusivity coefficient. The temperature gradient inside the cavity is generated by fixing the temperature along the top wall ${\mathrm{\Gamma }}_{top}$ and part of the bottom wall ${\mathrm{\Gamma }}_1$ according to $g_D$. As mentioned earlier, the vorticity minimization is achieved by controlling the volumetric heat source Q and the wall heat flux $g_N$ separately, however, the generic derivation for both types of control is reported below.

The weak formulation of the state constraints (1) and (2) is written using the linear, bilinear, and trilinear forms introduced earlier

$$\left\{\begin{array}{cc}b(\mathbf{u},f)=0\hfill & \forall f\in L^2\left(\mathrm{\Omega }\right),\hfill \\ \nu a(\mathbf{u},\psi )+c(\mathbf{u},\mathbf{u},\psi )=-b(\psi ,p)+(\beta \mathbf{g}(T-T_{ref}),\psi )\hfill & \forall \psi \in {\mathbf{H}}_{\partial \mathrm{\Omega }}^1\left(\mathrm{\Omega }\right),\hfill \\ \alpha a(T,\varphi )+c(\mathbf{u},T,\varphi )=(Q,\varphi )+{(g_N,\varphi )}_{{\mathrm{\Gamma }}_N}\hfill & \forall \varphi \in H_{{\mathrm{\Gamma }}_{top}\cup {\mathrm{\Gamma }}_1}^1\left(\mathrm{\Omega }\right).\hfill \end{array}\right.$$

(3)

Then we specialize the generic functional to this case, i.e.

$$\mathcal{F}(\mathbf{u},Q,g_N)={\displaystyle \frac{1}{2}}{\int }_{\mathrm{\Omega }}{|\nabla \times \mathbf{u}|}^{2}d\mathrm{\Omega }+{\displaystyle \frac{{\lambda }_Q}{2}}{\int }_{\mathrm{\Omega }}{\left|Q\right|}^{2}d\mathrm{\Omega }+{\displaystyle \frac{{\lambda }_N}{2}}{\int }_{{\mathrm{\Gamma }}_N}{\left|g_N\right|}^{2}d\mathrm{\Gamma },$$

(4)

where ${\lambda }_Q,{\lambda }_N$ are the regularization coefficients for the distributed and the Neumann boundary control cases. Note that, in the case of boundary control, since we are dealing with the temperature derivative, the regularization is intended to be in $L^2$. We also introduce the Lagrange multipliers $(\xi ,\mu ,\theta )$, referred to as the adjoint variables, corresponding to the adjoint velocity, adjoint pressure, and adjoint temperature, respectively. Finally, we write the Lagrangian as

$$\begin{array}{cc}\hfill \mathcal{L}(\mathbf{u},p,T,\xi ,\mu ,\theta )=& {\displaystyle \frac{1}{2}}{(\nabla \times \mathbf{u},\nabla \times \mathbf{u})}_{\mathrm{\Omega }}+{\displaystyle \frac{{\lambda }_Q}{2}}{(Q,Q)}_{\mathrm{\Omega }}+{\displaystyle \frac{{\lambda }_N}{2}}{(g_N,g_N)}_{{\mathrm{\Gamma }}_N}+\hfill \\ \hfill -& b(\mathbf{u},\mu )+\hfill \\ \hfill +& \nu a(\mathbf{u},\xi )+c(\mathbf{u},\mathbf{u},\xi )+b(\xi ,p)-(\beta \mathbf{g}(T-T_{ref}),\xi )+\hfill \\ \hfill +& \alpha a(T,\theta )+c(\mathbf{u},T,\theta )-(Q,\theta )-{(g_N,\theta )}_{{\mathrm{\Gamma }}_N}.\hfill \end{array}$$

(5)

Taking the variational derivative of the Lagrangian in the adjoint directions, we re-obtain the state equation system (3). Then we take the variation of the Lagrangian in the state directions $(\tilde{\mathbf{u}},\tilde{p},\tilde{T})$, obtaining the adjoint equation system

$$\left\{\begin{array}{cc}b(\xi ,\tilde{p})=0\hfill & \forall \tilde{p}\in L^2\left(\mathrm{\Omega }\right),\hfill \\ \nu a(\tilde{\mathbf{u}},\xi )+c(\tilde{\mathbf{u}},\mathbf{u},\xi )+c(\mathbf{u},\tilde{\mathbf{u}},\xi )-b(\tilde{\mathbf{u}},\mu )=(\nabla \times \mathbf{u},\nabla \times \tilde{\mathbf{u}})-c(\tilde{\mathbf{u}},T,\theta )\hfill & \forall \tilde{\mathbf{u}}\in {\mathbf{H}}_{\partial \mathrm{\Omega }}^1\left(\mathrm{\Omega }\right),\hfill \\ \alpha a(\tilde{T},\theta )+c(\mathbf{u},\tilde{T},\theta )=(\beta \mathbf{g}\tilde{T},\xi )\hfill & \forall \tilde{T}\in H_{{\mathrm{\Gamma }}_{top}\cup {\mathrm{\Gamma }}_1}^1\left(\mathrm{\Omega }\right).\hfill \end{array}\right.$$

(6)

Finally, we differentiate the Lagrangian in the two control variables direction $(Q,g_N)$

$$\begin{array}{cccc}\hfill {\lambda }_Q(Q,\tilde{Q})& =(\theta ,\tilde{Q})\hfill & \hfill & \forall \tilde{Q}\in L^2\left(\mathrm{\Omega }\right),\hfill \\ \hfill {\lambda }_N{(g_N,\widetilde{g_N})}_{{\mathrm{\Gamma }}_N}& ={(\theta ,\widetilde{g_N})}_{{\mathrm{\Gamma }}_N}\hfill & \hfill & \forall \widetilde{g_N}\in L^2\left({\mathrm{\Gamma }}_N\right),\hfill \end{array}$$

(7)

clearly taking ${\lambda }_N=0$ we solve the distributed control case, while taking ${\lambda }_Q=0$ we solve the boundary control case. From these two equations, we can derive the functional gradient used in the steepest descent method as

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \mathcal{L}}{\partial Q}}& =\theta +{\lambda }_{Q}Q,\hfill \\ \hfill {\displaystyle \frac{\partial \mathcal{L}}{\partial g_N}}& =\theta +{\lambda }_N{g}_N,\hfill \end{array}$$

(8)

respectively.

Finally, we write the strong form of the optimality system. We must isolate the velocity variation $\tilde{\mathbf{u}}$ to achieve this. To do so, we perform the following integrations by parts

$$\begin{array}{cc}\hfill (\nabla \times \mathbf{u},\nabla \times \tilde{\mathbf{u}})& =-{\int }_{\mathrm{\Omega }}(\nabla \times \nabla \times \mathbf{u})·\tilde{\mathbf{u}}d\mathrm{\Omega }+{\int }_{\mathrm{\Gamma }}\left(\left(\nabla \times \mathbf{u}\right)\times \mathbf{n}\right)·\tilde{\mathbf{u}}d\mathrm{\Gamma }\hfill \\ & =-{\int }_{\mathrm{\Omega }}\tilde{\mathbf{u}}\nabla (\nabla ·\mathbf{u})d\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}\tilde{\mathbf{u}}\Delta \mathbf{u}d\mathrm{\Omega }+{\int }_{\mathrm{\Gamma }}\left(\left(\nabla \times \mathbf{u}\right)\times \mathbf{n}\right)·\tilde{\mathbf{u}}d\mathrm{\Gamma },\hfill \\ \hfill c(\mathbf{u},\tilde{\mathbf{u}},\xi )& ={\int }_{\mathrm{\Omega }}\nabla \tilde{\mathbf{u}}:\left(\mathbf{u}\otimes \xi \right)d\mathrm{\Omega }=-{\int }_{\mathrm{\Omega }}\tilde{\mathbf{u}}·\left(\nabla ·\left(\mathbf{u}\otimes \xi \right)\right)d\mathrm{\Omega }+{\int }_{\mathrm{\Gamma }}\left(\tilde{\mathbf{u}}·\mathbf{n}\right)·\left(\mathbf{u}\otimes \xi \right)d\mathrm{\Gamma }\hfill \\ & =-{\int }_{\mathrm{\Omega }}\tilde{\mathbf{u}}·\left((\nabla ·\mathbf{u})\xi +\mathbf{u}·\nabla \xi \right)d\mathrm{\Omega }+{\int }_{\mathrm{\Gamma }}\tilde{\mathbf{u}}·\left(\mathbf{u}·\mathbf{n}\right)·\xi d\mathrm{\Gamma }.\hfill \end{array}$$

(9)

Considering both the divergence of the velocity and the adjoint velocity on the boundary to be zero, we obtain the following strong form for the adjoint velocity system of equations, which corresponds with the steady formulation of the Oseen equation  [34]

$$\left\{\begin{array}{cc}\nabla ·\xi =0\hfill & \mathrm{on}\mathrm{\Omega },\hfill \\ -2\nu \nabla ·\mathcal{D}\left(\xi \right)+\xi {(\nabla \mathbf{u})}^T-(\mathbf{u}·\nabla )\xi -\nabla \mu +\theta \nabla T=\Delta \mathbf{u}\hfill & \mathrm{on}\mathrm{\Omega },\hfill \\ \xi =0\hfill & \mathrm{on}\partial \mathrm{\Omega }.\hfill \end{array}\right.$$

(10)

Isolating the temperature variation $\tilde{T}$ in the adjoint system  (6) after some manipulation

$$\begin{array}{cc}\hfill a(\tilde{T},\theta )& =-{\int }_{\mathrm{\Omega }}\tilde{T}\Delta \theta d\mathrm{\Omega }+{\int }_{\mathrm{\Gamma }}\tilde{T}(\nabla \theta ·\mathbf{n})d\mathrm{\Gamma },\hfill \\ \hfill c(\mathbf{u},\tilde{T},\theta )& =-{\int }_{\mathrm{\Omega }}\mathbf{u}\tilde{T}·\nabla \theta +{\int }_{\mathrm{\Gamma }}(\mathbf{u}·\mathbf{n})\tilde{T}\theta d\mathrm{\Gamma },\hfill \end{array}$$

(11)

where in the second one we used the divergence-free property of $\mathbf{u}$. Since the temperature variation $\tilde{T}$ vanishes on the boundary ${\mathrm{\Gamma }}_i$ due to the Dirichlet condition on T, we obtain the following strong form system for the adjoint temperature $\theta$

$$\left\{\begin{array}{cc}-\alpha \Delta \theta -\mathbf{u}\nabla \theta =\beta \mathbf{g}·\xi \hfill & \mathrm{on}\mathrm{\Omega },\hfill \\ \theta =0\hfill & \mathrm{on}{\mathrm{\Gamma }}_i,\hfill \\ \nabla \theta ·\mathbf{n}=0\hfill & \mathrm{on}\partial \mathrm{\Omega }\setminus {\mathrm{\Gamma }}_i.\hfill \end{array}\right.$$

(12)

Lastly, the strong form of the control equations, both for the distributed and the boundary control case, is obtained directly from equations (7) leading to

$$\begin{array}{cc}\hfill Q& =-{\displaystyle \frac{\theta }{{\lambda }_Q}},\hfill \\ \hfill g_N& =-{\displaystyle \frac{\theta }{{\lambda }_N}}.\hfill \end{array}$$

(13)

2.2. Solid/Fluid Heat Transfer Problem

In this section, we derive the optimality system for the solid-fluid domain depicted in Figure 2. Moreover, during the following derivation, we will include some additional steps compared to the previous case for the sake of clarity. Unlike the previous case, the system here is considered to operate in a forced convection regime. As a result, the buoyancy term in the Navier-Stokes equations can be neglected. Furthermore, we impose a controlled inlet flow $\mathbf{w}$ on the boundary ${\mathrm{\Gamma }}_i$ and an outflow condition on ${\mathrm{\Gamma }}_o$. The corresponding system of equations reads

$$\left\{\begin{array}{cc}\nabla ·\mathbf{u}=0\hfill & \mathrm{on}{\mathrm{\Omega }}_f,\hfill \\ -2\nu \nabla ·\mathcal{D}\left(\mathbf{u}\right)+(\mathbf{u}·\nabla )\mathbf{u}+\nabla p=0\hfill & \mathrm{on}{\mathrm{\Omega }}_f,\hfill \\ \mathbf{u}=0\hfill & \mathrm{on}{\mathrm{\Gamma }}_c\cup {\mathrm{\Gamma }}_4,\hfill \\ \mathbf{u}=\mathbf{w}\hfill & \mathrm{on}{\mathrm{\Gamma }}_i,\hfill \\ \sigma (\mathbf{u},p)·\mathbf{n}=0\hfill & \mathrm{on}{\mathrm{\Gamma }}_o,\hfill \end{array}\right.$$

(14)

where the stress tensor is defined as $\sigma (\mathbf{u},p)=-p\mathbf{I}+2\nu \mathcal{D}\left(\mathbf{u}\right)$, and $\nu$ is the kinematic viscosity for the fluid. The temperature is solved both in the solid and the fluid domain, and on the interface, a conjugate heat transfer condition is imposed. In the following, we define a CHT interface condition as the condition of continuity of temperature as well as normal heat flux. We fix the temperature on ${\mathrm{\Gamma }}_i$ while we consider the other walls $\partial \mathrm{\Omega }\setminus {\mathrm{\Gamma }}_i$ to be insulated. Furthermore, in the solid region, a volumetric heat source Q is enforced. We write separately the system of equations for the solid region

$$\left\{\begin{array}{cc}-{\alpha }_s\Delta T_s=Q\hfill & \mathrm{on}{\mathrm{\Omega }}_s,\hfill \\ T_s=T_f\hfill & \mathrm{on}{\mathrm{\Gamma }}_c,\hfill \\ {\alpha }_s{\displaystyle \frac{\partial T_s}{\partial \mathbf{n}}}={\alpha }_f{\displaystyle \frac{\partial T_f}{\partial \mathbf{n}}}\hfill & \mathrm{on}{\mathrm{\Gamma }}_c,\hfill \\ {\displaystyle \frac{\partial T_s}{\partial \mathbf{n}}}=0\hfill & \mathrm{on}\partial {\mathrm{\Omega }}_s\setminus {\mathrm{\Gamma }}_c,\hfill \end{array}\right.$$

(15)

and for the fluid region

$$\left\{\begin{array}{cc}-{\alpha }_f\Delta T_f+\mathbf{u}·\nabla T_f=0\hfill & \mathrm{on}{\mathrm{\Omega }}_f,\hfill \\ T_f=1\hfill & \mathrm{on}{\mathrm{\Gamma }}_i,\hfill \\ T_f=T_s\hfill & \mathrm{on}{\mathrm{\Gamma }}_c,\hfill \\ {\alpha }_f{\displaystyle \frac{\partial T_f}{\partial \mathbf{n}}}={\alpha }_s{\displaystyle \frac{\partial T_s}{\partial \mathbf{n}}}\hfill & \mathrm{on}{\mathrm{\Gamma }}_c,\hfill \\ {\displaystyle \frac{\partial T_f}{\partial \mathbf{n}}}=0\hfill & \mathrm{on}{\mathrm{\Gamma }}_o\cup {\mathrm{\Gamma }}_4,\hfill \end{array}\right.$$

(16)

where ${\alpha }_s$ and ${\alpha }_f$ are the thermal diffusivity coefficients for the solid and fluid regions, respectively.

Now, we derive the weak form of the Navier-Stokes system after performing the following steps:

$$\begin{array}{cc}\hfill {\int }_{{\mathrm{\Omega }}_f}\psi \nabla ·\mathcal{D}\left(\mathbf{u}\right)d\mathrm{\Omega }& =-{\int }_{{\mathrm{\Omega }}_f}\mathcal{D}\left(\mathbf{u}\right):\mathcal{D}\left(\psi \right)d\mathrm{\Omega }+{\int }_{\mathrm{\Gamma }}\psi ·\mathcal{D}\left(\mathbf{u}\right)\mathbf{n}d\mathrm{\Gamma },\hfill \\ \hfill {\int }_{{\mathrm{\Omega }}_f}\nabla p·\psi d\mathrm{\Omega }& =-{\int }_{{\mathrm{\Omega }}_f}p\left(\nabla ·\psi \right)d\mathrm{\Omega }+{\int }_{\mathrm{\Gamma }}p\left(\psi ·\mathbf{n}\right)d\mathrm{\Gamma }.\hfill \end{array}$$

While in the Boussinesq case both boundary terms were zero due to the choice of the test function $\psi$, in this application they give a contribution on ${\mathrm{\Gamma }}_i$. Then, we can write the weak formulation for the Navier-Stokes equations system as follows

$$\left\{\begin{array}{cc}b(\mathbf{u},f)=0\hfill & \forall f\in L^2\left({\mathrm{\Omega }}_f\right),\hfill \\ \nu a(\mathbf{u},\psi )+c(\mathbf{u},\mathbf{u},\psi )+b(\psi ,p)-{\left(\sigma \mathbf{n},\psi \right)}_{{\mathrm{\Gamma }}_i}=0\hfill & \forall \psi \in {\mathbf{H}}_{{\mathrm{\Gamma }}_c\cup {\mathrm{\Gamma }}_4}^1\left({\mathrm{\Omega }}_f\right),\hfill \\ {(\mathbf{u}-\mathbf{w},\mathbf{s})}_{{\mathrm{\Gamma }}_i}=0\hfill & \forall \mathbf{s}\in {\mathbf{H}}^{-1/2}\left({\mathrm{\Gamma }}_i\right),\hfill \end{array}\right.$$

(17)

where the Dirichlet boundary condition at the inlet is also formulated weakly in order to obtain the control in the desired function space. For the temperature system, we first integrate the Laplacian term in the solid and the fluid systems to get

$$\begin{array}{cc}\hfill {\int }_{{\mathrm{\Omega }}_s}-{\alpha }_s\Delta T_s\varphi d\mathrm{\Omega }& =-{\int }_{{\mathrm{\Omega }}_s}{\alpha }_s\nabla T_s·\nabla \varphi d\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_c}{\alpha }_s(\nabla T_s·\mathbf{n})\varphi d\mathrm{\Gamma },\hfill \\ \hfill {\int }_{{\mathrm{\Omega }}_f}-{\alpha }_f\Delta T_f\varphi d\mathrm{\Omega }& =-{\int }_{{\mathrm{\Omega }}_f}{\alpha }_f\nabla T_f·\nabla \varphi d\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_c}{\alpha }_f(\nabla T_f·\mathbf{n})\varphi d\mathrm{\Gamma }.\hfill \end{array}$$

After summation of the above equations, the two boundary terms on ${\mathrm{\Gamma }}_c$ cancel each other due to the heat flux continuity. Then, we can write the temperature equation for the whole domain in weak form as

$$\alpha a(T,\varphi )+c{(\mathbf{u},T,\varphi )}_{{\mathrm{\Omega }}_f}={(Q,\varphi )}_{{\mathrm{\Omega }}_s}\forall \varphi \in L^2\left({\mathrm{\Omega }}_s\right),$$

(18)

where of course the diffusivity constant is $\alpha ={\alpha }_{s}H\left({\mathrm{\Omega }}_s\right)+{\alpha }_{f}H\left({\mathrm{\Omega }}_f\right)$ being H the Heaviside function.

Now we introduce the cost functional

$$\mathcal{F}(T,\mathbf{w})={\displaystyle \frac{1}{2}}{\int }_{{\mathrm{\Gamma }}_d}|T-T_d{|}^{2}d\mathrm{\Gamma }+{\displaystyle \frac{{\lambda }_1}{2}}{\int }_{{\mathrm{\Gamma }}_i}{\left|\mathbf{w}\right|}^{2}d\mathrm{\Gamma }+{\displaystyle \frac{{\lambda }_2}{2}}{\int }_{{\mathrm{\Gamma }}_i}{\left|{\nabla }_s\mathbf{w}\right|}^{2}d\mathrm{\Gamma },$$

(19)

where ${\mathrm{\Gamma }}_d\subset (\partial {\mathrm{\Omega }}_f\cup \partial {\mathrm{\Omega }}_s)$ is the observation boundary and an $H^1$ regularization is chosen since we deal with a Dirichlet boundary control [35]. Even if we use two different weights, ${\lambda }_1$ and ${\lambda }_2$, the existence of the solution for a Dirichlet boundary optimal control necessarily requires the regularization of the gradient term if we consider a CHT problem [4,36], i.e. ${\lambda }_2>0$. Using again the Lagrange multipliers $(\xi ,\mu ,\tau ,\theta )$, we write the Lagrangian $\mathcal{L}$ as

$$\begin{array}{cc}\hfill \mathcal{L}(\mathbf{u},p,T,\xi ,\mu ,\tau ,\theta )=& {\displaystyle \frac{1}{2}}{\left(T-T_d,T-T_d\right)}_{{\mathrm{\Gamma }}_d}+{\displaystyle \frac{{\lambda }_1}{2}}{\left(\mathbf{w},\mathbf{w}\right)}_{{\mathrm{\Gamma }}_i}+{\displaystyle \frac{{\lambda }_2}{2}}{\left({\nabla }_s\mathbf{w},{\nabla }_s\mathbf{w}\right)}_{{\mathrm{\Gamma }}_i}+\hfill \\ \hfill -& b(\mathbf{u},\mu )+\hfill \\ \hfill +& \nu a(\mathbf{u},\xi )+c(\mathbf{u},\mathbf{u},\xi )+b(\xi ,p)-{\left(\sigma \mathbf{n},\xi \right)}_{{\mathrm{\Gamma }}_i}+{(\mathbf{u}-\mathbf{w},\tau )}_{{\mathrm{\Gamma }}_i}+\hfill \\ \hfill +& \alpha a(T,\theta )+c{(\mathbf{u},T,\theta )}_{{\mathrm{\Omega }}_f}-{(Q,\theta )}_{{\mathrm{\Omega }}_s}.\hfill \end{array}$$

(20)

Taking the variational derivative of the Lagrangian in the adjoint direction, we re-obtain the state equations (17) and (18). In order to derive the adjoint velocity system, we differentiate in the velocity $\mathbf{u}$ and pressure p direction to get

$$\left\{\begin{array}{cc}b(\xi ,\tilde{p})=0\hfill & \forall \tilde{p}\in L^2\left({\mathrm{\Omega }}_f\right),\hfill \\ \nu a(\tilde{\mathbf{u}},\xi )+c(\tilde{\mathbf{u}},\mathbf{u},\xi )+c(\mathbf{u},\tilde{\mathbf{u}},\xi )-b(\tilde{\mathbf{u}},\mu )+{(\tilde{\mathbf{u}},\tau )}_{{\mathrm{\Gamma }}_i}+c(\tilde{\mathbf{u}},T,\theta )=0\hfill & \forall \tilde{\mathbf{u}}\in {\mathbf{H}}_{{\mathrm{\Gamma }}_c\cup {\mathrm{\Gamma }}_4}^1\left({\mathrm{\Omega }}_f\right),\hfill \\ {(\tilde{\sigma }\mathbf{n},\xi )}_{{\mathrm{\Gamma }}_i}=0\hfill & \forall \tilde{\sigma }\mathbf{n}\in {\mathbf{H}}^{-1/2}\left({\mathrm{\Gamma }}_i\right).\hfill \end{array}\right.$$

(21)

Similarly, by taking the variation of the Lagrangian in the temperature T direction, we obtain the weak form of the adjoint temperature equation

$$\alpha a(\tilde{T},\theta )+c{(\mathbf{u},\tilde{T},\theta )}_{{\mathrm{\Omega }}_f}=-{(T-T_d,\tilde{T})}_{{\mathrm{\Gamma }}_d}\forall \tilde{T}\in H_{{\mathrm{\Gamma }}_i}^1\left(\mathrm{\Omega }\right).$$

(22)

Lastly, differentiating the Lagrangian with respect to the control variable $\mathbf{w}$ we obtain the Euler equation in weak form

$${\displaystyle \frac{{\lambda }_1}{2}}{(\mathbf{w},\tilde{\mathbf{w}})}_{{\mathrm{\Gamma }}_i}+{\displaystyle \frac{{\lambda }_2}{2}}{({\nabla }_s\mathbf{w},{\nabla }_s\tilde{\mathbf{w}})}_{{\mathrm{\Gamma }}_i}={(\tau ,\tilde{\mathbf{w}})}_{{\mathrm{\Gamma }}_i}\forall \tilde{\mathbf{w}}\in {\mathbf{H}}^{-1/2}\left({\mathrm{\Gamma }}_i\right).$$

(23)

In order to iteratively solve this system through the gradient descent method, we need the explicit formulation for the Lagrangian derivative in the control direction, which clearly results in

$${\displaystyle \frac{\delta \mathcal{L}}{\delta \mathbf{w}}}={\lambda }_1\mathbf{w}-{\lambda }_2{\Delta }_s\mathbf{w}-\tau .$$

(24)

Once again, we transition from the weak form of the adjoint velocity system to the strong form. Specifically, we start isolating the velocity variation $\tilde{\mathbf{u}}$. To do so, we perform the same steps reported in (9), i.e. we perform some manipulations to obtain

$$\begin{array}{cc}\hfill a(\tilde{\mathbf{u}},\xi )& =-{\int }_{\mathrm{\Omega }}\tilde{\mathbf{u}}\nabla ·\mathcal{D}\left(\xi \right)d\mathrm{\Omega }+{\int }_{\mathrm{\Gamma }}\mathbf{u}·\mathcal{D}\left(\xi \right)\mathbf{n}d\mathrm{\Gamma },\hfill \\ \hfill b(\tilde{\mathbf{u}},\mu )& ={\int }_{\mathrm{\Omega }}\tilde{\mathbf{u}}·\nabla \mu d\mathrm{\Omega }-{\int }_{\mathrm{\Gamma }}(\tilde{\mathbf{u}}·\mathbf{n})\mu d\mathrm{\Gamma }.\hfill \end{array}$$

Taking the velocity variation $\tilde{\mathbf{u}}$ to vanish on ${\mathrm{\Gamma }}_c$ and ${\mathrm{\Gamma }}_4$ due to Dirichlet conditions, and the adjoint velocity $\xi$ to vanish on ${\mathrm{\Gamma }}_i$, ${\mathrm{\Gamma }}_c$, ${\mathrm{\Gamma }}_4$ as shown in the weak formulation of the system, we obtain the following strong form for the adjoint velocity system of equations

$$\left\{\begin{array}{cc}\nabla ·\xi =0\hfill & \mathrm{on}{\mathrm{\Omega }}_f,\hfill \\ -2\nu \nabla ·\mathcal{D}\left(\xi \right)+\xi {(\nabla \mathbf{u})}^T-(\mathbf{u}·\nabla )\xi -\nabla \mu +\theta \nabla T=0\hfill & \mathrm{on}{\mathrm{\Omega }}_f,\hfill \\ \xi =0\hfill & \mathrm{on}{\mathrm{\Gamma }}_i\cup {\mathrm{\Gamma }}_c\cup {\mathrm{\Gamma }}_4,\hfill \\ -\mu \mathbf{n}+2\nu \mathcal{D}\left(\xi \right)\mathbf{n}+(\mathbf{u}·\mathbf{n})\xi =0\hfill & \mathrm{on}{\mathrm{\Gamma }}_o,\hfill \\ \tau =-\mu \mathbf{n}+2\nu \mathcal{D}\left(\xi \right)\mathbf{n}\hfill & \mathrm{on}{\mathrm{\Gamma }}_i.\hfill \end{array}\right.$$

(25)

Performing the same steps as in (11), we write the strong form equation system for the adjoint temperature $\theta$ as

$$\left\{\begin{array}{cc}-\alpha \Delta \theta -\mathbf{u}\nabla \theta =(T-T_d)H_{{\mathrm{\Gamma }}_d}\hfill & \mathrm{on}\mathrm{\Omega },\hfill \\ \theta =0\hfill & \mathrm{on}{\mathrm{\Gamma }}_i,\hfill \\ \nabla \theta ·\mathbf{n}=0\hfill & \mathrm{on}\partial \mathrm{\Omega }\setminus {\mathrm{\Gamma }}_i.\hfill \end{array}\right.$$

(26)

Lastly, the strong form of the control equation is easy to obtain directly from equation (23)

$${\lambda }_1\mathbf{w}-{\lambda }_2{\Delta }_s\mathbf{w}=\tau \mathrm{on}{\mathrm{\Gamma }}_i.$$

(27)

2.3. Minimization Method

The optimality system consists of the state equations for $(T,\mathbf{u},p)$, the adjoint equations for $(\theta ,\xi ,\mu )$, and the control equation, or optimality condition, for $Q,g_N$, and $\mathbf{w}$. These three systems of equations are clearly coupled but are solved in a segregated manner, as mentioned, to reduce computational effort. The state and adjoint systems are solved using a simple fixed-point iteration method, while for the control variables, to ensure convergence to the optimum, we employ a basic gradient method with a backtracking line search algorithm, as shown in Algorithm 1. Here we refer to the control variable with the generic symbol q. It can be demonstrated that, for a sufficiently small step size that the backtracking line search algorithm always guarantees, the gradient method yields descent direction and is globally convergent [37]. As mentioned, the control update in line 8 is performed using a gradient method. Firstly, the control variation $\delta q$ in the three different cases is calculated as

$$\begin{array}{cc}\hfill & \delta Q=-\rho (\theta +{\lambda }_{Q}Q),\hfill \\ \hfill & \delta g_N=-\rho (\theta +{\lambda }_N{g}_N),\hfill \\ \hfill & {\lambda }_1\delta \mathbf{w}-{\lambda }_2\Delta \delta \mathbf{w}=-\rho (-\tau +{\lambda }_1\mathbf{w}-{\lambda }_2{\Delta }_s\mathbf{w}),\hfill \end{array}$$

and then the update is performed

$$q_{k+1}=q_k+\delta q.$$

Algorithm 1:Steepest gradient method

1: Set $k=0$ 2: Compute state $T_k,{\mathbf{u}}_k,p_k$ 3: Compute functional ${\mathcal{F}}_k$ 4: Compute adjoint ${\theta }_k,{\xi }_k,{\mu }_k$ 5: while true do 6:     Set initial step size $\rho ={\rho }_0$ 7:     while ${\mathcal{F}}_{k+1}>{\mathcal{F}}_k\mathrm{and}\rho >ϵ$ do 8:         Update control $q_{k+1}$ 9:         Compute state $T_{k+1},{\mathbf{u}}_{k+1},p_{k+1}$ 10:         Compute functional ${\mathcal{F}}_{k+1}$ 11:         Update step size $\rho =\rho /2$ 12:     end while 13:     if $\left|{\displaystyle \frac{{\mathcal{F}}_{k+1}+{\mathcal{F}}_k}{{\mathcal{F}}_{k+1}}}\right|<\tau$ then 14:         break 15:     end if 16:     Compute adjoint ${\theta }_{k+1},{\xi }_{k+1},{\mu }_{k+1}$ 17: end while

3. Coupling Algorithm

In this section, we present the numerical algorithm for the code coupling that solves the state equations separately from the adjoint and control equations. The previous section has given a detailed explanation on the reasoning to adopt a finite element framework for solving the adjoint and control equations. Moreover, the coupling application leverages the well-validated OpenFOAM solvers, making it a reliable tool for solving the state equations and eliminating the need for re-implementing traditional conservation equations.

Data exchange between the two codes has been made possible by implementing three C++ classes. The first class, med_coupling, manages the data transfer routines with the external library MEDCOUPLING, specifically for data manipulation and interpolation. Furthermore, we have two specific wrapper classes for the CFD codes. The first one, referred to as of_interface, is responsible for the communication between the routines of OpenFOAM and the MED data structure. The corresponding class for the FEMuS code is labeled as femus_interface. More details regarding the specific implemented classes and the coupling strategy can be found in [24]. With the previously described code architecture, we describe the specific algorithms for the Rayleigh-Bénard problem and the solid/fluid heat transfer problem in the following sections.

3.1. Rayleigh-Bénard Problem

For the first application, we want to couple FEMuS and OpenFOAM in order to solve a Rayleigh-Bénard problem describing the natural convection inside a cavity. The state equations to be solved are the Navier-Stokes equations in the incompressible formulation together with the temperature equation, as described in the previous section. This system of equations is then solved by adding the adjoint and the control equations, which can be both a volumetric source or a boundary condition for the temperature equation. In this work, we present these two different data type transfers, i.e. a volumetric and a boundary approach. For more details regarding the data transfer with the external library MEDCOUPLING, the interested reader can consult [24], where each routine used in this work is explained.

The data exchange from OpenFOAM to FEMuS involves transferring the state variables: the velocity $\mathbf{u}$ and the temperature T fields. Thus, the get and set functions extract these volumetric fields from OpenFOAM and write them into FEMuS for the solution of the adjoint equations. Note that this data exchange implies a transfer from a finite volume representation to a finite element data structure. The same approach applies to the control variable q, where the get and set functions transfer this variable from FEMuS to OpenFOAM. This data exchange can be used in both volumetric and boundary-based methods, with data interpolation now originating from a finite element grid and converting data to a finite volume grid. For the volumetric data transfer, q is the volumetric source of the temperature equation, while for boundary control, q is transferred to the Neumann boundary condition of the same equation. A schematic representation of this algorithm is reported in Figure 3.

3.2. Solid/Fluid Heat Transfer Problem

For this test, we exploit already validated numerical tools available in OpenFOAM for the solution of the state equations. In particular, we refer to the multi-region technique, which provides tools to exchange data between two regions, such as the conjugate heat transfer approach.

The schematic representation for the application is reported in Figure 4. The data transfer between the CFD codes follows the same volumetric and boundary data exchange described in the previous application. However, the of_interface class has been improved to support the multi-region solver of OpenFOAM, allowing for the retrieval of field values within different domains. In this application, the temperature field is transferred from OpenFOAM to FEMuS not only across the entire fluid domain but also at the interface between the two regions. Moreover, OpenFOAM supplies the volumetric velocity field to FEMuS, while FEMuS provides the boundary control q to the fluid domain in OpenFOAM.

It is worth noting that OpenFOAM exploits a segregated approach in the implementation of the CHT boundary condition between the solid and the fluid, i.e. the two regions are solved separately, and then the boundary conditions are enforced on ${\mathrm{\Gamma }}_c$, taking into account the continuity of the temperature field and the conservation of the heat fluxes. On the other hand, FEMuS solves the adjoint temperature in a monolithic way on a single domain which includes both the solid and the fluid region. Therefore, we do not enforce additional boundary conditions on the common interface ${\mathrm{\Gamma }}_c$ since the assembling of a unique matrix allows us to satisfy the continuity conditions on temperature and heat flux in a natural way. Of course, different material properties and source terms are taken into account depending on the specific domain where the local element of the mesh is assembled.

Lastly, we remark that OpenFOAM solves the Navier-Stokes system only by a projection method, and we selected the PIMPLE algorithm. On the other hand, a monolithic approach is implemented in FEMuS for the adjoint Navier-Stokes, which indeed couples the adjoint velocity and pressure by assembling and solving them simultaneously.

3.3. Numerical Discretization

We briefly discuss the discretization framework of the equations. We already know that the state equations are solved with a finite volume approach while the adjoint and control equations use a finite element discretization. In the introduction, we also mentioned the use of the optimize then discretize concept to build our optimal control problem that supports different discretization spaces for the approximation of the state and adjoint equations. We introduce the finite dimensional spaces $X_h\subset X=H_{{\mathrm{\Gamma }}_D}^1\left(\mathrm{\Omega }\right)$, $Q_h\subset L^2\left(\mathrm{\Omega }\right)$ and $S_h\subset L^2\left(\mathrm{\Omega }\right)$. Specifically, given a domain partition ${\mathcal{T}}_h$, the element K, and the space ${\mathcal{P}}_r$ of polynomials of degree at most r, we set

$$\begin{array}{cc}\hfill & X_h=\{v_h\in C^0\left(\overline{\mathrm{\Omega }}\right):v_h{|}_K\in {\mathcal{P}}_2\forall K\in {\mathcal{T}}_h\},\hfill \\ \hfill & Q_h=\{v_h\in C^0\left(\overline{\mathrm{\Omega }}\right):v_h{|}_K\in {\mathcal{P}}_1\forall K\in {\mathcal{T}}_h\},\hfill \\ \hfill & S_h=\{v_h\in L^2\left(\mathrm{\Omega }\right):v_h{|}_K\in {\mathcal{P}}_0\forall K\in {\mathcal{T}}_h\}.\hfill \end{array}$$

We refer to the vector-valued version of the same spaces with a bold symbol. Regarding the discretized adjoint variables, we consider a standard Taylor-Hood pair for the adjoint velocity and pressure [38,39], thus we have $({\xi }_h,{\mu }_h)\in {\mathbf{X}}_h\times Q_h$. For the adjoint temperature and the control, we consider biquadratic finite elements, i.e. ${\theta }_h,q_h\in X_h$.

The finite volume approximation of the state equations requires only one type of discretized space. Therefore, we consider the piecewise constant approximation for velocity, pressure and temperature, i.e. ${\mathbf{u}}_h\in {\mathbf{S}}_h$ and $p_h,T_h\in S_h$.

These discretization spaces interact through the numerical interpolation of the physical variables from a finite volume to a finite element grid and vice-versa. A comprehensive analysis of the interpolation techniques used can be found in [40,41].

4. Numerical Results

In this section, we present two numerical examples in order to validate the effectiveness of the above-presented algorithm. As a first numerical test, a Rayleigh-Bénard natural convection problem is described considering a square cavity. Later, we present a multi-region problem with a solid and a fluid domain that interact via a standard conjugate heat transfer model. For each application, we recall the objective functional, the control term in the state equations, and the corresponding boundary conditions. Numerical results are then compared with already published works with the same problem setting.

4.1. Rayleigh-Bénard Problem

For the first numerical test, we refer to the example reported by Lee [3], where a vorticity minimization has been studied on a square domain. Specifically, the natural convection obtained from the Rayleigh-Bénard configuration has been simulated with both a distributed and a Neumann boundary control for the temperature equation. Among the several cases reported in [3], we selected two of them that present different behaviors of the control and the optimal solution. Therefore, after a comparison of the uncontrolled solution, the cases with ${\lambda }_Q={10}^{-3}$ and ${\lambda }_Q={10}^{-5}$ will be discussed (same values have also been considered for the Neumann case, i.e. for ${\lambda }_N$).

In Figure 5, the domain $\mathrm{\Omega }={[0,1]}^2$ is shown, with the partition of the bottom boundary into two sections, ${\mathrm{\Gamma }}_1$ and ${\mathrm{\Gamma }}_N$. The boundary conditions for the temperature field are as follows: a homogeneous Neumann is enforced on ${\mathrm{\Gamma }}_l$ and ${\mathrm{\Gamma }}_r$; on ${\mathrm{\Gamma }}_{top}$ we have a homogeneous Dirichlet boundary condition, while on ${\mathrm{\Gamma }}_1$ we have a fixed temperature value $T=10$. Lastly, on ${\mathrm{\Gamma }}_N$ we have the nonhomogeneous Neumann boundary condition that acts as control.

The grid used for the simulation consists of $20\times 20$ elements for the FE code and $40\times 40$ elements for the FV code, both equally spaced. This discretization setup leads to a very similar amount of degrees of freedom for each code (the finite element grid is indeed equipped with biquadratic elements, therefore implying 1681 dofs).

In the case of the distributed control, we fix ${\lambda }_N=0$, while for the Neumann boundary control we set ${\lambda }_Q=0$. Regarding the material properties of the fluid, we set $\nu =\beta =\alpha =1$.

We recall the optimality system in strong form for this problem, given by (28), (29) and (34), along with the corresponding boundary conditions reported in Table 1.

$$State:\left\{\begin{array}{cc}\nabla ·\mathbf{u}=0\hfill & \mathrm{on}\mathrm{\Omega },\hfill \\ -2\nu \nabla ·\mathcal{D}\left(\mathbf{u}\right)+(\mathbf{u}·\nabla )\mathbf{u}=-\nabla p+\beta \mathbf{g}(T-T_{ref})\hfill & \mathrm{on}\mathrm{\Omega },\hfill \\ -\alpha \Delta T+(\mathbf{u}·\nabla )T=Q\hfill & \mathrm{on}\mathrm{\Omega },\hfill \end{array}\right.$$

(28)

$$Adjoint:\left\{\begin{array}{cc}\nabla ·\xi =0\hfill & \mathrm{on}\mathrm{\Omega },\hfill \\ -2\nu \nabla ·\mathcal{D}\left(\xi \right)+\xi {(\nabla \mathbf{u})}^T-(\mathbf{u}·\nabla )\xi =\nabla \mu -\theta \nabla T+\Delta \mathbf{u}\hfill & \mathrm{on}\mathrm{\Omega },\hfill \\ -\alpha \Delta \theta -(\mathbf{u}·\nabla )\theta =\beta \mathbf{g}·\xi \hfill & \mathrm{on}\mathrm{\Omega },\hfill \end{array}\right.$$

(29)

$$\begin{array}{cccccccc}\hfill Control:& \hfill & \hfill Q& =-{\displaystyle \frac{\theta }{{\lambda }_Q}}\mathrm{on}\mathrm{\Omega },\hfill & \hfill \mathrm{or}& \hfill & \hfill g_N& =-{\displaystyle \frac{\theta }{{\lambda }_N}}\mathrm{on}{\mathrm{\Gamma }}_N.\hfill \end{array}$$

(30)

4.1.1. Uncontrolled Case

The results for the uncontrolled case (namely, the control variable is set to zero) are shown in Figure 6. The contour plot of the temperature field is shown on the left side, while the velocity field is shown on the right.

4.1.2. Distributed Control

Figure 7 shows the results for the temperature fields in the distributed control case, compared to the same contour plot from [3] which is shown with a dashed line. We observe that both temperature distributions get closer and closer to a purely vertical temperature gradient. It appears that the reduction in the horizontal temperature gradient leads to a temperature profile and hence a buoyancy force that is nearly constant along a single horizontal line. This evens out the horizontal differences in buoyancy and contributes to the minimization of the vorticity. A good agreement with respect to the reference paper can be observed for the isolines of the temperature field, in particular for the lowest ${\lambda }_Q$ value.

The contour plot of the control Q is shown in Figure 8, on the left for the case with ${\lambda }_Q={10}^{-3}$ and on the right for ${\lambda }_Q={10}^{-5}$. We recall that the control Q acts as the volumetric source of the temperature equation, and we can see the different distributions for the two values of ${\lambda }_Q$. The lower the regularization factor ${\lambda }_Q$, the more the Q field concentrates in magnitude at the bottom of the domain.

Figure 9 shows the adjoint velocity vector glyphs for each controlled case.

Similarly, Figure 10 shows the velocity vector glyphs for each controlled case. Note that in this case we have rescaled the velocity profiles by a common scalar, in order to visualize the overall reduction in magnitude.

Specifically, the computed values of ${\left|\mathbf{u}\right|}_{max}$ are reported in Table 2, where it can be noticed that the velocity magnitude decreases by increasing the control contribution, i.e. with a lower ${\lambda }_Q$ value.

In Table 3, we report the value of the objective term, the regularization term, and the total functional for the distributed case, with a comparison to the same results from [3]. In the last column, we have reported the relative error of the total functional, computed as

$${\epsilon }_r={\displaystyle \frac{|{\mathcal{F}}_C-{\mathcal{F}}_R|}{{\mathcal{F}}_R}},$$

(31)

where the subscripts C and R stand for Coupled and Reference, respectively. For the same reason, the results presented in this work are labeled with C in the table. As we can see, the proposed algorithm is able to reproduce the same values of the cost functional terms with an error below $2\%$, validating the accuracy of the numerical code coupling. A good agreement is evident in each of the above cases. Slight differences appear for the uncontrolled case, while the results of the code coupling algorithm become closer and closer to the literature ones by decreasing the control parameter ${\lambda }_Q$.

4.1.3. Neumann Boundary Control

Regarding the Neumann boundary control, in Figure 11 we plot the control $g_N$ along the boundary ${\mathrm{\Gamma }}_N$ for each controlled case (${\lambda }_N={10}^{-3}$ and ${\lambda }_N={10}^{-5}$). In this case, this control represents the Neumann boundary condition for the temperature equation, which can be understood as the wall heat flux provided in order to minimize the vorticity. For the highest ${\lambda }_N$ case the control $g_N$ seems to grow monotonically, while for ${\lambda }_N={10}^{-5}$ we can notice a peak in the left part of ${\mathrm{\Gamma }}_N$, which is attached to ${\mathrm{\Gamma }}_1$ where we have a fixed temperature boundary condition. We remark that on the point between ${\mathrm{\Gamma }}_1$ and ${\mathrm{\Gamma }}_N$ the adjoint temperature $\theta$ is zero. Hence, $\theta$ cannot help pushing the control at this point and so the control $g_N$ tends to have a local maximum nearby in order to compensate for the lack of push.

In Figure 12, the contour plot for the temperature is reported for the Neumann boundary control case for each ${\lambda }_N$, together with the reference case [3] represented by dashed lines.

Some discrepancies appear between the solid isolines and the corresponding dotted lines as reported in [3], especially for the case ${\lambda }_N={10}^{-3}$ in the bottom part of the domain, where indeed we have the control $g_N$. This situation improves with the lower ${\lambda }_N$.

The velocity fields are reported in Figure 13 for both controlled cases. Moreover, in Table 4 we report a comparison of the maximum value of the velocity magnitude for the computed cases, and we notice a difference of one order of magnitude between the results of the two ${\lambda }_N$ values.

In Table 3 we report the value of the objective term, the regularization term, and the total functional for the Neumann boundary control case, with a comparison of the same results from [3].

Differently from the distributed control case, the values computed in Table 5 for the Neumann control are sometimes significantly distant from the literature. Anyhow, the total value of the functional $\mathcal{F}$ with the proposed algorithm is smaller than the one obtained in [3].

4.2. Solid/Fluid Heat Transfer Problem

In the second test, we reproduce the numerical examples proposed in [4], which we use as a benchmark. Here a fluid domain exchanges heat with a neighboring solid region through a common boundary. This configuration is typically referred to as a conjugate heat transfer problem. The geometry is described in Figure 14, which shows two neighboring regions, ${\mathrm{\Omega }}_s$ for the solid and ${\mathrm{\Omega }}_f$ for the fluid. The boundary walls of the fluid domain are given by ${\mathrm{\Gamma }}_4$ and ${\mathrm{\Gamma }}_c$, which is the common interface with the solid. On the walls a no-slip condition is enforced, i.e., a homogeneous Dirichlet boundary condition for the velocity field. On the other hand, ${\mathrm{\Gamma }}_i$ and ${\mathrm{\Gamma }}_o$ are the inlet and the outlet boundaries of the channel through which the fluid flows. We impose a Dirichlet boundary condition at the inlet and a stress-free condition at the outlet.

The boundaries of the solid domain are ${\mathrm{\Gamma }}_1$, ${\mathrm{\Gamma }}_2$, ${\mathrm{\Gamma }}_3$ and ${\mathrm{\Gamma }}_c$. The boundary conditions for the solid region are given by a Neumann boundary condition for each external edge, while on the common interface ${\mathrm{\Gamma }}_c$ we enforce the continuity of the temperature and of the heat flux between the two domains. The aforementioned benchmark case has been obtained within a finite element environment. This presents some challenges when transitioning to the finite volume method.

In the FE formulation, the point at the junction ${\mathrm{\Gamma }}_i\cap {\mathrm{\Gamma }}_1$ belongs to boundary sections that require a temperature Dirichlet condition on one side and a Neumann condition on the other. In the FE context with Lagrange shape functions, the point itself is a degree of freedom carrier, and the Dirichlet condition can be enforced in essential way.

On the other hand, in the FV framework, the degrees of freedom are given by the cell centers, and a Dirichlet boundary condition involves the face of an element. In this situation, the intersection point will not behave in the FV code in the same way as in the FE code. Therefore, the FE versus the FV solutions show a rather different behavior at this intersection point. We have observed with auxiliary numerical tests that the results in [4] have been obtained with a mesh size that yields a FE solution that is quite far from grid convergence. In spite of that, in order to compare with the reference paper values, we decided to choose an FV mesh with approximately the same dof count. At such a resolution, the Dirichlet value computed by the FV code on the lower-left corner of the solid domain turns out to be quite far from the FE value. We were able to fix this difference by simply enforcing a Dirichlet condition also on the first solid element above ${\mathrm{\Gamma }}_c$.

We use a $40\times 40$ element grid for the FE code and a similar-sized (from the d.o.f. point of view) $80\times 80$ cell grid for the FV code to replicate the results proposed in [4]. The latter grid is divided as a $60\times 80$ grid in the fluid region and a $20\times 80$ grid in the solid region, refined towards the solid/fluid interface ${\mathrm{\Gamma }}_c$. Concerning the thermal diffusivity, we have ${\alpha }_s=1$ and ${\alpha }_f=0.1$, while the kinematic viscosity for the momentum equations is $\nu =1/1000$. Lastly, the volumetric heat source in the solid region is set $Q=6$.

We recall the optimality system in strong form for this problem, given by Eqs. (32), (33) and (34), along with the boundary conditions in Table 6:

$$State:\left\{\begin{array}{cc}\nabla ·\mathbf{u}=0\hfill & \mathrm{on}{\mathrm{\Omega }}_f,\hfill \\ -2\nu \nabla ·\mathcal{D}\left(\mathbf{u}\right)+(\mathbf{u}·\nabla )\mathbf{u}=-\nabla p\hfill & \mathrm{on}{\mathrm{\Omega }}_f,\hfill \\ -{\alpha }_f\Delta T_f+(\mathbf{u}·\nabla )T_f=0\hfill & \mathrm{on}{\mathrm{\Omega }}_f,\hfill \\ -{\alpha }_s\Delta T_s=Q\hfill & \mathrm{on}{\mathrm{\Omega }}_s,\hfill \\ T_s=T_f\hfill & \mathrm{on}{\mathrm{\Gamma }}_c,\hfill \\ \mathbf{n}·(k_s\nabla T_s-k_f\nabla T_f)=0\hfill & \mathrm{on}{\mathrm{\Gamma }}_c,\hfill \end{array}\right.$$

(32)

$$Adjoint:\left\{\begin{array}{cc}\nabla ·\xi =0\hfill & \mathrm{on}{\mathrm{\Omega }}_f,\hfill \\ -2\nu \nabla ·\mathcal{D}\left(\xi \right)+\xi {(\nabla \mathbf{u})}^T-(\mathbf{u}·\nabla )\xi =\nabla \mu -\theta \nabla T+\Delta \mathbf{u}\hfill & \mathrm{on}{\mathrm{\Omega }}_f,\hfill \\ \tau =-\mu \mathbf{n}+2\nu \mathcal{D}\left(\xi \right)\mathbf{n}\hfill & \mathrm{on}{\mathrm{\Gamma }}_i,\hfill \\ -\alpha \Delta \theta -(\mathbf{u}·\nabla )\theta =\beta \mathbf{g}·\xi \hfill & \mathrm{on}\mathrm{\Omega },\hfill \end{array}\right.$$

(33)

$$Control:{\lambda }_1\mathbf{w}-{\lambda }_2{\Delta }_s\mathbf{w}=\tau \mathrm{on}{\mathrm{\Gamma }}_i.$$

(34)

Notice that the state equations are reported explicitly over the two regions since OpenFOAM solves in a segregated way, while the adjoint uses the monolithic formulation. We also remark that, as in [4], the Dirichlet condition on the inlet velocity ${\mathbf{u}|}_{{\mathrm{\Gamma }}_i}$ is enforced only in the normal direction.

4.2.1. Definition of the Target Temperature

Following [4], the target of the control problem is the temperature profile $T_d$ on ${\mathrm{\Gamma }}_d\subset {\mathrm{\Gamma }}_c$ that is obtained with an analytical input velocity ${\mathbf{u}}_0{|}_{{\mathrm{\Gamma }}_i}={(2(1.5y-2y^2),0)}^T$, i.e. a Poiseuille flow. In the following, we refer to this as the uncontrolled solution. The target $T_d$ is then fed to the optimal control system. For the solution of the optimality system, the first control in the steepest descent method is initialized with a different inlet boundary condition for the velocity field, i.e. ${\mathbf{u}|}_{{\mathrm{\Gamma }}_i}={((1.5y-2y^2),0)}^T$, and therefore the control acts on the interface boundary to reach a velocity field that matches the temperature distribution $T_d$ as close as possible. It appears clear that the optimal solution should provide a control on the boundary that reproduces the original boundary inlet condition for the velocity field ${\mathbf{u}}_0$. In this way, it is easy to check the quality of the optimal control solution. Of course, due to the presence of a regularization term, i.e. a parameter ${\lambda }_i>0$, the control is not able to reproduce exactly that original velocity profile, but it provides the closest solution according to the value of ${\lambda }_i$.

In Figure 15, we report the target temperature solution. On the left we show the overall behavior as a contour map and on the right we plot the trend on a portion of the interface ${\mathrm{\Gamma }}_c$, which corresponds to the target temperature $T_d$. We recall that the objective functional was defined only on a portion of the solid/fluid interface, i.e. the region ${\mathrm{\Gamma }}_d=\{(x,y)\in {\mathbb{R}}^2:x\in [0.075,1],y=0.75\}$. The isolines in the left figure confirm the interface conditions of the CHT model, i.e. the continuity of the temperature and the change in slope due to the continuity of the heat flux across different conductivities.

4.2.2. Dirichlet Boundary Control

In Figure 16, the Dirichlet velocity control on ${\mathrm{\Gamma }}_i$ is reported as a solid line for two controlled cases. In particular in Case 1 we set the pair $({\lambda }_1,{\lambda }_2)$ to $(5\times {10}^{-2},5\times {10}^{-3})$, while in Case 2 we set $({\lambda }_1,{\lambda }_2)=(8\times {10}^{-3},8\times {10}^{-4})$. The comparison with the reference data of [4] is shown with circular markers.

We notice a good agreement between our results and the reference data for each controlled case. The inlet velocity control that we obtain is not symmetric with respect to the midpoint $(0,0.375)$ of ${\mathrm{\Gamma }}_i$. In fact, since we want to match the temperature on a portion of ${\mathrm{\Gamma }}_c$, and since the solid has a volumetric heat source, we expect to have an advection contribution of the fluid that cools down the solid, i.e., a velocity field which is able to remove the heat generated. Therefore, the peak of the velocity profile is shifted towards the interface boundary. We also notice, as expected, that a lower combination of the values of ${\lambda }_i$ yields a velocity profile closer to the target one represented by the dashed line.

In Figure 17, we show the temperature profiles on ${\mathrm{\Gamma }}_d$ and the corresponding relative percentage errors with respect to the target value for both controlled cases. The temperature profile obtained with the initial guess of the inlet velocity ${\mathbf{u}|}_{{\mathrm{\Gamma }}_i}$ is also reported. Regarding the relative error with respect to $T_d$, we define

$${\epsilon }_r={\displaystyle \frac{|T_i-T_d|}{T_d}}·100,$$

(35)

where $T_i$ corresponds to the two controlled cases.

As expected, the numerical temperature field does not coincide with the target temperature $T_d$ since we weight the control contribution with the ${\lambda }_i$ parameters. As a result, the case with the lowest ${\lambda }_i$ values provides a temperature distribution on ${\mathrm{\Gamma }}_c$ that is closer to $T_d$. In fact, Case 2 presents a maximum error below $0.2\%$, almost three times less than Case 1, which reaches a value close to $0.6\%$. This aspect is confirmed later when we provide the values of the cost functional terms, specifically the objective term.

The adjoint velocity and the adjoint temperature have been reported in Figure 18 for both controlled cases. More precisely, we report the vector glyph for the adjoint velocity on the left and the contour map for the adjoint temperature on the right. With regards to this latter plot, the adjoint temperature is multiplied by ${10}^3$, since its absolute values are relatively low.

As we know from the literature [4,42], the adjoint velocity of a fluid channel configuration with an inlet and an outlet behaves like a flow moving in the opposite direction. Indeed, the adjoint inlet is now ${\mathrm{\Gamma }}_o$, while the other boundaries behave as a wall for the adjoint velocity, since velocity Dirichlet boundary conditions are enforced on them.

In analogy with the previous test, in Table 7 we report the values of the cost functional terms, both for the uncontrolled case and for the two controlled cases. Unlike the boundary control of the Rayleigh-Bénard application, in this case the boundary control is in ${\mathbf{H}}^1\left({\mathrm{\Gamma }}_i\right)$, hence we take into account two contributions for the regularization part, which are given in separate columns.

The table confirms the expected behavior of the functional value in the last column, showing a lower number for the better controlled case, i.e. the one with the lowest ${\lambda }_i$. In fact, with the latter values of ${\lambda }_i$, we obtain a total functional almost six time lower than the other controlled case. This reduction is mainly due to the objective term, which has a decrease of one order of magnitude, despite the regularization terms being similar.

4.3. Steepest Descent Convergence Order

In order to verify the correctness of our numerical implementation, we have computed the order of convergence of the steepest descent algorithm for all the proposed cases. We know from the theory that the steepest descent method has a linear order of convergence [43]. This test allows us to check the reliability of the coupling routines and the magnitude of the error made in the coupling process. We recall the formula for the convergence order p as

$$p=\underset{k\to \infty }{lim}{\displaystyle \frac{log\left(\right||{\mathcal{F}}_{k+1}-{\mathcal{F}}_{min}|\left|\right)}{log\left(\right||{\mathcal{F}}_k-{\mathcal{F}}_{min}|\left|\right)}},$$

(36)

where k represents the numerical iterations and ${\mathcal{F}}_{min}$ is the minimum.

Table 8 and Table 9 present the computed convergence order p for all the analyzed cases. This value is obtained by averaging the convergence order values of the steepest descent iterations computed using formula (36). The results show good agreement with the theoretical value of the gradient descent algorithm, which is equal to one. Therefore, we can confidently affirm that the use of segregated codes with different discretization techniques for the solution of optimal control systems does not affect the convergence order of the steepest descent method.

5. Conclusions

This paper presents an extension of a code coupling technique for optimal control problems already published by the authors. Specifically, we have considered optimal control problems for two fluid heat transfer applications, the first one described by the Rayleigh-Bénard system for a natural convection cavity configuration while the second one considers the conjugate heat transfer between a solid and a fluid region. The presented approach is based on the solution of the optimality system employing two different CFD codes and the data exchange between them using the external MEDCOUPLING library. In fact, the state equations are solved in the finite volume code OpenFOAM while the adjoint and the control system are discretized with the finite element code FEMuS.

Two numerical examples have been reported with a comparison of the same application setup taken from the literature. We achieve good agreement with reference literature data obtained with a monolithic approach, validating the goodness of the numerical approach.

The computational approach described in this paper can be a valuable tool for spreading the optimal control features into widely-used CFD codes, enabling an engineering design approach based on physical target values to be reached.

References yes