Discrete adjoint for coupled conjugate heat transfer

The typical method to solve multi-physics problems such as Conjugate Heat Transfer (CHT) is the partitioned approach which couples separate solvers through boundary conditions. Effective gradient-based optimisation of partitioned CHT problems requires the adjoint of the coupling to maintain the efficiency of the original multi-physics coupling, which is a significant challenge. The use of automatic differentiation (AD) has the potential to ease this burden and leads to generic gradient computation methods. In this paper, we present how to automate the generation of adjoint fluid and solid solvers for a CHT adjoint using Automatic Differentiation (AD). The derivation of the adjoint of the loose coupling algorithms is shown for three fixed-point coupling algorithms. Application of the coupled adjoint algorithm is shown to two CHT optimisation benchmark cases for inverse design and shape optimisation. The results demonstrate that Robin-based coupling algorithms have faster runtimes and are an attractive option for CHT optimisation problems.


Introduction
Conjugate heat transfer (CHT) describes the process of heat transfer between a fluid and solid and is ubiquitous in engineering applications such as turbine blade cooling, modelling of heat exchangers and cooling of electronics.
CHT problems may be solved using a monolithic approach in which both fluid and solid equations are solved simultaneously by a single solver. However, the partitioned or segregated approach is often adopted where separate solvers for the fluid and structure are loosely coupled through boundary conditions. These conditions need to be updated iteratively until the temperature and heat flux are continuous between the two domains (Ganine et al. 2013;Verstraete and Scholl 2016).
Recent work in CHT has evolved beyond merely solving the CHT problem to an increased interest in shape optimisation (Mousavi and Nadarajah 2011;Zeinalpour et al. 2016;Ferlauto 2014;Gkaragkounis et al. 2018).This has led to the need for efficient optimisation methods such as gradient-based approaches. Although gradient-based approaches are only guaranteed to converge to local minima, they are preferred because they typically require less function evaluations. This is advantageous in applications like CHT where the cost of each function evaluation can be high. Adjoint methods have been shown to be highly efficient as the cost of obtaining gradients can be made almost independent of the number of design variables. However, for partitioned coupling approaches, the flexibility of using different solvers for both domains results in an increased level of complexity with regard to obtaining the required gradients.
Adjoint methods can be grouped into continuous or discrete methods. In the continuous method the adjoint equations are analytically derived before being discretised whilst the discrete method discretises the state equations before formulating the discrete adjoint equations (Giles et al. 2003). Arguments in favour of either approach revolve around stability, accuracy and computational effort. The discrete adjoint is considered to have the advantage with regards to stability and accuracy but not with computational effort. However, the use of automatic differentiation (AD) can significantly reduce the implementation effort associated with producing the adjoint code.
The majority of the work related to CHT optimisation has favoured the continuous adjoint formulation (Mousavi and Nadarajah 2011;Mousavi 2012;Zeinalpour et al. 2016;Gkaragkounis et al. 2018) whilst Burghardt and Gauger (2020) use the discrete adjoint but with a monolithic solver. In the present study, we demonstrate a partitioned/segregated coupling methodology in which the discrete adjoint formulation is achieved through AD. The use of AD has the potential to lead to very generic and less labour intensive implementations of adjoint methods for CHT. We focus on coupling boundary conditions and highlight the advantages of Robin boundary conditions in the fluid domain. Furthermore, the gradient exchange required by partitioned coupling approaches is demonstrated using three fixed-point coupling algorithms. This paper is organised as follows: we first describe the CHT problem, governing equations and the coupling procedure. We then discuss the adjoint procedure for the partitioned approach and gradient verification. Two CHT optimisation problems are then solved using the adjoint method. Finally, a summary is given in the conclusion.

Primal problem
Let us consider flow over a flat plate with finite thickness. The free stream flow temperature is T ∞ , whilst the bottom of the plate is maintained at a lower temperature T b . Consequently, heat is transferred at the interface between the solid plate and the fluid. The primal problem is to accurately compute the wall temperature at the interface between the fluid and solid, which is unknown a priori and can only be computed by considering the coupled problem (see Fig 1).
In order to solve the CHT problem using a partitioned approach, the fluid and solid governing equations must be coupled leading to a coupled system of equations where U denotes the fluid state variables, W the solid state variables and i the coupling iteration. F represents the Reynolds Averaged Navier-Stokes equations: where t denotes pseudo time, Q the source term. The state vectors U, and the inviscid and viscous flux vectors ⃗ f c and ⃗ f v are defined as where, , p, and ⃗ v are the fluid density, pressure and velocity vector, respectively, ⃗ n the surface normal vector, e the internal energy per unit mass, is the stress tensor for Newtonian fluids and is the fluid thermal conductivity, ̃ is the modified eddy viscosity which is obtained using the Spalart-Allmaras turbulence model (Allmaras and Johnson 2012) and S in Eq. (2) refers to the governing equation of the solid domain Ω s , that is the steady state heat conduction equation where s is the conductivity of the solid.
Each domain depends on the other through the boundary conditions specified at the fluid-solid interface. The solid state W, which is the solid's temperature, affects the fluid state U through the viscous flux component of the energy equation and through the heat flux at the non-adiabatic fluid-solid interface which must be specified using a boundary condition. Similarly, the fluid state variables U affect the state of the solid W through the interface boundary conditions specified whilst solving the heat equation.
Consequently, separate stand alone solvers for the fluid and solid, which are loosely coupled through interface boundary conditions, can be used to solve the primal problem. These boundary conditions are updated iteratively until the temperature and heat flux are continuous between the two domains. That is until where T fw is the interface temperature in the fluid domain Ω f , T sw the interface temperature in the solid domain Ω s , f and s the thermal conductivity of the fluid and solid, respectively, n the surface normal, and q fw and q sw the interface heat flux in the fluid and solid domains, respectively (see Fig. 1).
In this work, the fluid equations are solved using the in-house mgOpt flow solver, a vertex centred, finite volume solver, which solves the 3-D compressible RANS equation using unstructured grids (Xu et al. 2015;Christakopoulos 2013). The equations for the solid are solved using the finite element method with the opensource solver, CalculiX, developed by Dhondt (2004). These separate solvers are loosely coupled to solve CHT problems.

Coupling algorithms
In the present work, an implicit coupling method is used to solve the primal problem with a partitioned approach. Different coupling algorithms exist depending on which type of boundary conditions are exchanged between both domains. In both domains, boundary conditions could be specified as Dirichlet, Neumann, or Robin, leading to a total of 7 different types of coupling algorithms (Scholl 2017). In this work, we use the three different fixed-point coupling algorithms described in the following sections.

Temperature forward flux back (TFFB)
In the TFFB method (Divo et al. 2003), the solid interface heat flux distribution, q i sw , where i is the current coupling iteration, is imposed as a boundary condition to the fluid domain. The fluid solver F solves the flow equations resulting in a fluid interface temperature distribution, T i fw . This temperature is then imposed as boundary condition for the solid domain and the solid conduction solver, S , provides an updated heat flux distribution q i+1 sw . This loop is continued until there is no change in the boundary conditions exchanged by both solvers (see Fig. 2),

Temperature forward solid coefficient back (TFRB)
The TFRB coupling algorithm imposes a Dirichlet boundary condition in the solid domain and uses a Robin boundary condition in the fluid domain. The method also requires a virtual conductivity R to be specified, with R =L , ̃ being the virtual solid conductivity, and L a solid length scale. The values of ̃ and L are non-physical quantities chosen by the user which affect the speed of convergence and stability of the method (Scholl 2017;Scholl et al. 2018). Assuming the algorithm begins in the solid domain, an initial guess for the wall temperature T i sw is imposed on the solid and used to obtain the heat flux q i sw . Next, the virtual solid sink temperature T i s is calculated using Eqn. (10). The virtual conductivity and solid sink temperature are used for a Robin boundary condition in the fluid domain (Eqn. (12)). The flow solver then returns an update of the interface temperature which is given to the solid as a Dirichlet boundary condition (see Fig. 2).

Heat transfer coefficient forward solid coefficient back (hFRB)
The hFRB method uses Robin boundary conditions in both domains and is the same as TFRB on the fluid side (see Fig. 2). Similar to the TFRB method, the algorithm can start with an initial guess of the interface temperature on the solid side. Next, the virtual solid temperature T s is then calculated and passed to the flow solver along with the virtual conductivity R .
The outputs from the flow solver are the interface temperature and heat flux ( T fw , q fw ). These are used in Eqn. (17) to calculate the ambient fluid temperature ( T i sink ) for a user specified value of the virtual heat transfer coefficient ( h ). The solid solver then returns an update on the interface temperature and heat flux and the exchange is continued until convergence (Scholl 2017;Scholl et al. 2018).

Discrete adjoint
For gradient-based CHT optimisation, we require the gradient of the objective function I w.r.t the design variables where U and W represent the fluid and solid state variables, respectively. The term g T u is expensive to solve hence it is advantageous to use the adjoint formulation. The diagonal terms are the Jacobians of each discipline whilst the off-diagonal terms show how the states of one discipline affect the state of the other e.g. how the fluid temperature affects the solid flux and vice-versa. Therefore the sensitivity of the cost function, Eqn. (19), for the coupled problem can be written as The cost of calculating the gradient can be reduced through the adjoint method. The adjoint variables are the solution to the adjoint equation represents the fluid and solid adjoint variables, respectively, and Therefore, after only one solve of Eqn. (26) for the adjoint variable v, the sensitivity can be calculated as In the partitioned coupling approach, the Jacobian of the coupled system in Eqn. (26) is not calculated. Instead, the adjoint solution is obtained through an iterative approach, similar to the primal coupling, to compensate for the missing off-diagonals. The adjoint system of each discipline is solved using the adjoint solution the other discipline from the previous iteration (Martins et al. 2002) where i is the current coupling iteration. The partial derivatives ( S U , F W ) in equations (31) and (32) depend on the type of coupling algorithm used (see Table 1) and are obtained by differentiating the solvers w.r.t. the coupling boundary conditions and coordinates. Reverse mode automatic differentiation is applied to the fluid and heat conduction solvers (Imam-Lawal 2020) using Tapenade (Hascoët and Pascual 2013), a source-transformation AD tool. The use of AD significantly reduces the effort required to obtain the adjoint solvers and in the case of the fluid solver, the procedure was fully automated. This allows for easy implementation of new coupling boundary conditions and optimisation objective functions.
The iterative approach shown in equations (31) and (32) results in a reversed/inverted version of each coupling algorithm. The total sensitivity in Eqn. (30) is then obtained by accumulating the gradients output after each solid adjoint iteration.

Gradient calculation and verification
We demonstrate accuracy of the partitioned adjoint methodology using a flat plate test case. The plate has a fixed temperature T b at the bottom and comes in contact with fluid of a different temperature (see Fig 3). A 3D CHT simulation is performed to obtain the interface temperature and heat flux. A perturbation in the temperature of a red node at the bottom results in a change in the heat flux into the fluid domain. This perturbation travels through the coupling and results in a new interface temperature T w at the blue node. Similarly, a coordinate perturbation (x, y, z) changes the volume of the plate and leads to a change in the heat flux at the fluid-solid interface and consequently alters the solution of the coupled problem. Therefore, the effect of these perturbations on the interface temperature is described by two gradients y, z]. The partitioned adjoint approach to solving equations (31) and (32) is now described.

Temperature forward flux back (TFFB)
The adjoint run of the TFFB algorithm, starts with seeding a vector T i fw which is set to 1 for the blue node in Fig. 3 whilst all others are set to 0. T fw is used by the adjoint flow solver ( ) to obtain the adjoint heat flux q sw . This is then passed directly to the adjoint solid ( ) solver to obtain an update of the adjoint temperature. The reverse loop is performed for the same number of iterations as for the primal solution and at the end of each coupling iteration, the gradient T b is accumulated. This single adjoint solve obtains the gradient of the interface temperature at the blue node w.r.t to all design variables (red nodes). A block structure representation of the reverse differentiated coupling algorithm is shown in Fig. 4. In fact, both dT w dT b and dT w d⃗ x are obtained in one reverse solve highlighting the advantage and cost savings of the adjoint method. Table 2 shows gradients for whilst Table 3 shows the gradients for dT w d⃗ x . Both Tables show good agreement between the tangent, adjoint and central difference (CD) methods. Bold numbers in Tables 2-7 highlight the differences in values between the approaches.

Temperature forward solid coefficient back (TFRB)
The adjoint TFRB run is started by seeding a vector of the adjoint fluid temperature T fw to the flow solver to obtain T sb . The reverse differentiated routine for Eqn. (10) is also required to obtain the adjoint heat flux and temperature as shown in equations (36). The solid solver then uses the adjoint heat flux  and dT w d⃗ x and is shown in Tables 4 and 5. Good agreement is seen between the tangent, adjoint and central difference (CD) methods. (35)

Heat transfer coefficient forward solid coefficient back (hFRB)
For the hFRB adjoint, the adjoint wall temperature T fw is used by the flow solver to obtain the adjoint solid sink temperature T s . This is then converted into an adjoint heat flux q sw and adjoint temperature T sw using the differentiated routine which calculates the Robin parameters T sink and h for the solid domain. The solid solver then returns an adjoint sink temperature T sink and adjoint virtual heat transfer coefficient h . The differentiated Robin preprocessing step uses the virtual heat transfer coefficient h to calculate the adjoint wall temperature T fw whilst the adjoint sink temperature T sink is assigned to the first off wall node, as shown in Eqn. (42).

Advantage of fluid Robin-based coupling algorithms
All the gradients obtained in Tables 2 -7 were obtained by fully converging the fluid and solid primal and adjoint equations as well as performing an equal number of primal and reverse coupling iterations. However, for all coupling algorithms, significant time savings can be obtained by only partially converging the fluid primal and adjoint equations each time without significantly impacting the accuracy of the gradients (Imam-Lawal 2020). Furthermore, rather than performing an equal number of primal and reverse coupling iterations, fewer reverse coupling iterations could be performed, at the cost of a slight reduction in gradient accuracy. On the flat plate validation case in Sec. 3.1, for Robin-based coupling algorithms, reducing the number of reverse coupling iterations from 19 to 3 only resulted in a difference of approximately 0.5% in the gradients. This reduction in the number of reverse coupling iterations can also be combined with partial convergence of the fluid adjoint to significantly reduce runtime without great loss of gradient accuracy. However the TFFB algorithm always required the same number of reverse iterations as the primal to obtain accurate gradients.
Consequently, Robin-based coupling algorithms reduce the wall clock time required to obtain reasonably accurate gradients (Imam-Lawal 2020). The Robin-based coupling algorithms have also been shown by Verstraete and Scholl (2016), Scholl (2017) to require less coupling iterations to converge the primal. Therefore, the time saving in both the primal and adjoint coupling runs will significantly speed up the optimisation process.

Inverse design optimisation
To demonstrate the efficacy of the partitioned adjoint methodology, we start with a simple inverse problem for which the final solution is known. Inverse problems are solved by providing a desired solution and adjusting design variables in order to achieve the desired target. In this problem, we  seek the bottom temperature ( T b ) which results in the best match for a given the interface wall temperature ( T target ) as shown in Fig. 5. Table 8 and Fig. 1 show the remaining parameters.
The inverse problem is solved by formulating an optimisation problem, which allows the use of classical direct methods to solve the physics involved. An objective function (J) is defined as the difference between the desired interface temperature ( T target ) and the obtained interface temperature ( T w ) for an estimated bottom temperature ( T b ). and minimising Eqn. (44) results in an interface temperature which matches the target temperature. The objective function depends implicitly on the bottom temperature T b through the solution of the coupled problem. Each mesh node at the bottom of the plate has an independent value of T b specified as a boundary condition, and is used in this work as a design variable ( ) that needs to be changed to drive J to zero. The selected mesh following a mesh dependence study results in 226 design variables for the present work. The objective is hence to minimise Eqn. (44) subject to the constraints of satisfying both the state equations of both domains and maintain continuity of state variables ( T w , q w ) across the interface.
A gradient-based method is used to reduce the deviation of the current interface wall temperature with the desired one. The gradient of the objective function (sensitivity) w.r.t the control variables, , is given as The gradients of the temperature w.r.t the design variables, i.e. the temperature specified on the bottom of the flat plate, are computed using the adjoint approach for calculating described in Sec. 3. The target temperature is obtained by solving the primal problem with a bottom temperature of 600K, hence it is guaranteed that a solution to the problem exists. T b is initially taken as 400K and refers to the estimated bottom temperature that should yield T target .

Results
The objective is minimised using the BFGS optimisation algorithm from the SciPy library (Jones et al. 2014). The difference between the temperatures obtained from the direct and inverse solution is defined as Figure 6a shows that the reduction in the objective function for all three coupling algorithms. The results for the TFRB algorithm in Fig. 6b show that the inverse solution is significantly closer to the target than the initial guess. Similar results were obtained for the TFFB and hFRB coupling algorithms (omitted for clarity). The successful solution of the inverse problem shows that the partitioned discrete adjoint methodology described is effective for solving CHT optimisation problems. Table 9 shows the time taken for the gradient norm to fall below 0.1 for all coupling algorithms. The Robin-based algorithms are seen to be faster and this is attributed to needing fewer primal and reverse coupling iterations. Despite the loss of gradient accuracy due to partially converging the fluid adjoint and performing ≈ 80% less reverse coupling iterations than primal coupling iterations, the Robin-based algorithms solve the inverse problem. This shows that the time savings are worth the sacrifice in gradient accuracy.

MarkII turbine blade
Modern gas turbines are equipped with internal cooling channels which cool the internal structure and prevent damage from high turbine inlet temperatures. Hence, an optimisation problem is formulated to minimise the maximum temperature in a turbine blade by changing only the location of each cooling channel.
The vane geometry is modelled after the Mark II turbine vane which has been investigated experimentally by Hylton et al. (1983). The blade is convectively cooled by ten cooling channels and a 2D simulation of the problem is carried out. Matching meshes are used for both domains as shown in Figure 7with 49,532 nodes in the fluid domain and 5,714 nodes in the solid. A near wall spacing y + of less than 1 is used for the fluid.
In the fluid domain, boundary conditions for the freestream pressure and temperature, as well as the Mach number at inlet are specified, whilst a gauge pressure of (46) Error = T Target − T Inverse . zero is specified at the outlet. In the solid domain, a constant temperature is imposed in each of the cooling channels (Table 10). The blade is made of ASTM 310 stainless steel and the thermal conductivity is a function of temperature taken as where T is the temperature at a point in the solid. The density and heat capacity are taken as 7900 kgm −3 and 586.5 J kg −1 K −1 respectively. We define an objective function (J) as the maximum temperature in the solid domain where p is a user-defined integer taken as 10 and T ref is a user-defined constant. The objective function depends on the solid temperature field which is obtained by solving the CHT problem. The CHT problem was solved using the hFRB algorithm due to stability reasons.
(47) s = 6.811 + 0.020176 ⋅ T,  The coordinates (x, y) of each of the cooling channels are taken as control variables ( ) that need to be changed to drive J to zero. As a result, the control variable in the discrete problem is an array of size N, where N is the number of channels (10), times the x&y coordinates leading to a total of 20. The gradients of the objective function w.r.t the design variables are computed using the partitioned adjoint approach described in Sec. 3. Figures 8 and 9 show the temperature, pressure, and heat transfer coefficient distributions on the baseline geometry. The computational surface temperature distribution deviates from the experiment near the leading edge but this discrepancy reduces closer to the trailing edge. The primary cause of the discrepancy on the suction surface of the leading edge is the due to the limitation of the Spalart-Allmaras turbulence model. This model is not well suited for modelling transitional flow and transitional models can be used to better match the experiment (Hongjun et al. 2013;Hao et al. 2014;Lin et al. 2014).

Baseline analysis
At the trailing edge where we see good agreement between the experiment and CHT results, the presence of the cooling channels leads to oscillations in the interface temperature and heat transfer coefficient distributions. We also see that the solid temperature is highest at the trailing edge.

Optimisation results
An unconstrained optimisation is performed using the BFGS optimisation algorithm from the SciPy library (Jones et al. 2014). The Inverse Distance Weighted (IDW) interpolation method is used to propagate the displacement of the cooling channels to the internal solid mesh nodes. The interface boundary nodes are kept fixed to maintain the match with the fluid domain. The IDW algorithm is also reverse differentiated to obtain fully accurate adjoint gradients of the entire design chain.
We achieve approximately 0.9% reduction in the objective function which corresponds to a 3.14 K drop in the maximum temperature as shown in Fig. 10. The irregular shape of the curves in Fig. 10 are a result of the SciPy line search. Figure 11a shows the displacement of the cooling channels. The black lines represent the channels of the optimised blades whilst the solid surface is the initial geometry and temperature distribution. The optimisation is terminated prematurely as the first cooling channel (channel 1) moves out of the solid domain in the next optimisation step. The channels at the leading edge (channels 1-3) are displaced upwards and outwards towards the higher temperature regions, the trailing edge channels (channels 7-10) are displaced downwards, whilst the middle channels (channels 4-6) are displaced towards the suction surface. Figure 11b shows the node based change in the blade temperature distribution plotted on the initial geometry, where Δ T is the difference between the optimised and baseline temperature. The greatest reduction in temperature is seen at the leading edge near the new location of channel 1.

Conclusion
A partitioned discrete adjoint method for the partitioned approach to Conjugate Heat Transfer modelling has been presented. Loose coupling of the fluid and solid domains with three fixed-point coupling algorithms was presented and their adjoints were derived. The discrete adjoint solvers were developed using reverse mode automatic differentiation (AD) which reduces the effort of development and in particular maintenance. The process of developing the fluid adjoint solver was fully automated which provides a strong case for adopting the discrete adjoint over the continuous adjoint for multidisciplinary optimisation problems. The exchange of sensitivity information between solvers was also described and the gradients verified. The presented methodology was then applied to two CHT

Declarations
Conflict of interest On behalf of all authors, the corresponding author confirms that there are no conflict of interest relating to this work. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.