A port-Hamiltonian formulation of coupled heat transfer

ABSTRACT Heat transfer and cooling solutions play an important role in the design of gas turbine blades. However, the underlying mathematical coupling structures have not been thoroughly investigated. In this work, the port-Hamiltonian formalism is applied to the conjugate heat transfer problem in gas turbine blades. A mathematical model based on common engineering simplifications is constructed and further simplified to reduce complexity and focus on the coupling structures of interest. The model is then cast as a port-Hamiltonian system and examined for stability and well posedness.


Introduction
With the German government's plans to increase the share of renewable energy in the power grid to over 60% by 2050, gas turbines will most likely gain in importance and take on new roles and tasks in the power grid. Due to their short start-up times and high efficiency, they are particularly well suited as reserve power plants and can cushion power drops and demand peaks. They are also relevant in the 100% renewable energy source scenario, when hydrogen and methane are used as energy storage. These changing application roles go hand in hand with changing design requirements, especially in terms of efficiency, reliability and flexibility of operation.
In order to accurately incorporate these requirements into the design process, the use of high-level multiphysics simulations is required, combining fluid dynamics, structural mechanics, heat conduction, convective heat transport and 1D flow networks, among others. The GivEn project [1] aims to integrate multiphysics simulation with a multicriteria shape optimization process.
To achieve this goal and obtain useful results, other requirements must first be met. One of them is to ensure a suitable coupling structure between the different parts of the multiphysics simulation. In this work, we want to model and investigate a special coupling structure -the heat transfer at the walls of the cooling channels of the turbine blade. It couples the heat conduction in the metal of the turbine blade with the transport of the cooling fluid flowing through the cooling channels.
To improve the thermal efficiency, modern gas turbines are operated at very high temperatures (1200°C to 1500°C). Since these temperatures greatly exceed the range in which the turbine blade's metal can be used safely, active cooling of the turbine blades is necessary for them to withstand these temperatures. One of the techniques used for this purpose is convection cooling, i.e. cooling channels are installed within the blade. These internal cooling channels are small ducts within the turbine blade that are filled with a stream of cooling fluid, usually (comparatively cool) air extracted from the compressor. Figure 1 gives a graphical representation of the cooling channels within a turbine blade, as well as other cooling techniques not discussed here.
This approach leads to a so-called conjugate heat transfer (CHT) problem (strong thermal interactions between solids and fluids) [3], which is the heat transfer between the turbine blade, the internal flow in the cooling channel and the hot external flow. In this work, we will focus on the role of the cooling channels, where the heating of the blade by the external flow is discussed in [3]. In order to maximize heat transfer between the blade and the fluid, the flow within the cooling channels is kept deliberately turbulentfor example, by use of so-called rib turbulators, periodic protrusions and recessions in the channel walls. For a more in-depth review of turbine blade cooling, see, e.g [2,4]. The latter work also contains additional graphical illustrations of the turbine blade and its cooling channels. A discussion of a more realistic model for the flow within the internal cooling channels themselves can be found in [5,6]. Although the resulting model is a bit dated from today's point of view, the underlying considerations are still valid. Numerical simulations for the coupled system of blade and cooling channels can be found in [7,8] although both works focus more on modelling each subsystem and only remark on the coupling in passing. We formulate a so-called port-Hamiltonian system (pHs) modelling this conjugate heat transfer and study in detail the resulting coupling structure. Since it is closely related to the Hamiltonian formalism originally developed in theoretical physics, this port-Hamiltonian framework is a natural fit for modelling physical systems and, in particular, their interconnections since two port-Hamiltonian systems connected by a suitable coupling structure in turn form a port-Hamiltonian system. The formalism also makes conservation laws, a fundamental property of virtually any physical system, explicit. Moreover, a suitable port-Hamiltonian formulation makes the process of discretizing a continuous system for numerical simulation relatively simple while ensuring that conservation laws still hold in the discretized system [9]. Moreover, several desirable properties such as stability or controllability are either inherent to port-Hamiltonian systems or can be guaranteed by some easily checked additional conditions.
The paper is structured as follows. First, in Section 2 we motivate and introduce a mathematical model of the coupled system to be investigated, including a rescaling of the variables. Next, in Section 3 we rewrite each subsystem of our model in the port-Hamiltonian framework, after which we combine these port-Hamiltonian systems and study the properties of the resulting coupled system. Finally, we summarize and interpret the results in Section 5 and give some concluding remarks.

The model system
Our model system is based on a highly simplified model of heat transfer within the blade of a gas turbine since we are mainly interested in the mathematical coupling structure.
In the design process of gas turbine blades, much work is done with twodimensional slices, which are then stacked and interpolated to form the final threedimensional blade. This is done because simulation and optimization is much easier and faster that way than with a full three-dimensional model and still gives 'good enough' results. The cooling channel, which is more or less perpendicular to the slices, then has exactly one point of contact with each slice (neglecting 180 � hairpin curves) where heat transfer can occur.
Of course, the simplifications made in this work mean that the model system we consider is not suitable for any kind of realistic engineering work or quantitative simulation of the real system. This is purely meant as a way to qualitatively investigate the coupling structure. A system relevant for engineering purposes would, at the very least, need to replace the simple transport equations modelling the cooling channel by a more complex model for compressible flows as discussed in, e.g. [5]. It would also need to model the heat equation in the more complex geometries of the two-dimensional turbine blade slices. A full system modelling the entire turbine blade including the actual three-dimensional flow within the cooling channels and its interaction with the external gas is still infeasible as far as we are aware.
Since we are primarily interested in the thermal coupling structure between the turbine blade and the cooling channel, and a 2D-slice introduces additional complexity into the model, we decided to further simplify the system and use a one-dimensional rod instead. Thus, our model consists of a thermally conductive rod (a � x m � b) and a cooling channel (i � x c � o) through which a fluid is pumped. The left end of the rod at x m ¼ a is in contact with a thermal reservoir at a given temperature T ext ðtÞ. The right end of the rod at x m ¼ b is in contact with the wall of the cooling channel. Figure 2 gives a graphical representation of this arrangement.
The temperature of the metal rod T m ¼ T m ðx m ; tÞ is modelled by a heat equation supplied with Robin boundary conditions at x m ¼ a and x m ¼ b. The temperatures T in ¼ T in ðx c ; tÞ and T out ¼ T out ðx c ; tÞ of the inflowing and outflowing parts of the cooling channel are described by simple transport equations. This is, again, a simplification as it assumes that the convective heat transport dominates, and we neglect the diffusion in the cooling channel medium. For the usual flow rates and cooling fluids used in gas turbines, such as air or water vapour, this is a valid assumption as they have a very low thermal conductivity compared to their heat capacity. The coupling at the point x c ¼ c is such that the temperature of the outflowing fluid at the point c is determined by the temperature of the inflowing at this point plus the temperature change due to heat flowing out of the metal rod according to the boundary condition at x m ¼ b.

The system of equations
The above setting leads to the following equations for the temperatures in the metal rod and in the cooling channel: with the following boundary conditions: Hence, to summarize, the temperature field in the metal rod is described by the heat Equation (1a) supplied with the two Robin boundary conditions (1d), (1e) modelling a Fourier-type heat transfer. Written more expansively, the heat equation can be considered as a combination of Fourier's law q ¼ À k @T @x , stating that the heat flux q is proportional to the negative temperature gradient, together with a Dulong-Petit-like law, linking the change of internal energy Q to the temperature change cρ @T @t ¼ @Q @t ¼ À @q @x . The boundary conditions take the form of Newton's law of cooling. Next, in the inflow part of the cooling channel (i < x c < c) a given temperature profile T inflow ðtÞ at the left boundary x c ¼ i is convected with the speed v, see (1b), (1f). Finally, in Equation (1g) the coupling of the two systems at the point x m ¼ b equals x c ¼ c is described the heat flux from the metal rod to the cooling channel, depending on the inflowing temperature T in at x c ¼ c.
We remark that, due to Equation (1e), the last boundary condition (1g) can alternatively be written as c c vðT out ðcÞ À T in ðcÞÞ ¼ À k @T m @x m ðbÞ: (1h)

The rescaled system
A rescaling of system (1) such that the temperature is defined on a unit interval ½0; 1� allows us to write the system (1) in a more compact way, and we see how geometric dimensions enter the system. For this purpose, we introduce new space variables and the rescaled temperature functions # j ð� j ðx m=c Þ; tÞ ¼ T j ðx m=c ; tÞ for each of the indices j 2 fm; in; outg respectively. Restating the system in three scaled spatial variables � j then results in with the following boundary conditions: Now we are prepared to rewrite this system (2) as a port-Hamiltonian system in the next section.

The port-Hamiltonian formulation of infinite-dimensional systems
As mentioned earlier, we now want to formulate our model system in the port-Hamiltonian framework as this makes it easier to check for certain properties, especially stability.
Since the state variables # m , # in , # out in the model system (2) are continuous in space, we cannot apply the finite-dimensional port-Hamiltonian framework as described in, for example [10,11]. Instead, we use a generalization for distributed parameter systems as presented in, for example [12][13][14][15][16]. Since the system (2) only contains first-order derivatives w.r.t. time, we also restrict ourselves to linear first-order systems for simplicity. However, the second-order derivative in space representing thermal diffusion, in combination with our choice of Hamiltonian, will cause our system to be dissipative as seen later on. We therefore need to use a formalism that accounts for dissipation.
Port-Hamiltonian systems (pHs) can be viewed as a combination of a Dirac structure (or Stokes-Dirac structure in the case of infinite dimensional systems) and a Hamiltonian. We start with a definition of the underlying Dirac structure. Definition 3.1 ((Stokes-)Dirac structure [12,14,16,17], cf. also [15]).
Let F be a linear space, E its dual and h�; �i : E � F ! R their dual product. Further Note that E and F here do not just contain the storage ports but can also contain dissipative ports and external ports -which in the case of Stokes-Dirac structures contain distributed and boundary ports, in turn.

Remark 1.
In the context of port-Hamiltonian systems, we mainly encounter the following special case of a (Stokes-)Dirac structure: This results in the following definition of a pHs, as given in [13, Definition 7.1.2]: Definition 3.2 (port-Hamiltonian system [13]). Let P 1 2 R n�n invertible and selfadjoint, P 0 2 R n�n skew-adjoint, i.e. P T 0 ¼ À P 0 and H 2 R n�n symmetric such that mI � H � MI with constants m; M > 0. Further, let X ¼ L 2 ð½a; b�; R n Þ be a Hilbert space with the inner product gð�Þ T Hð�Þf ð�Þd�: Then the differential equation is a linear first-order port-Hamiltonian system (pHs) on a one-dimensional space with the associated Hamiltonian Θð�; tÞ T Hð�ÞΘð�; tÞd�:

Remark 2.
With the usual choice of f ¼ @Θ @t and e ¼ HΘ, we can see the connection between the port-Hamiltonian system of Definition 3.2 and the Stokes-Dirac structure of Definition 3.1, since the operator The most important difference between a Dirac structure and a Stokes-Dirac structure is the presence of a boundary port that governs the power flow across the boundary and takes the place of boundary conditions in 'regular' PDEs. For a port-Hamiltonian system, as in Definition 3.2, the boundary port takes the following form, cf. [ where f @ and e @ can be considered the boundary flow and boundary effort, respectively, and R 0 the transformation matrix mapping the effort variables at the boundaries to the boundary port. While this is essentially a simple variable substitution, it makes the power flow across the boundary obvious, since now dH dt ¼ e T @ f @ holds (in the absence of other external ports).
We equip the port-Hamiltonian system with boundary conditions of the form where uðtÞ is a time-dependent input function, W B 2 R n�2n has full rank, W B �W T B � 0 and � ¼ 0 We note that the above property guarantees that the port-  (5) and (6), we see that Many port-Hamiltonian systems do not fit directly into the formalism above, especially if the order of space and time derivatives does not match. However, the formalism can be extended to dissipative systems, cf. [14,Chapter 6]: The operator G � R is the formally adjoint operator of G R (i.e. the adjoint of G R neglecting boundary conditions) and S is a coercive operator on L 2 ð½0; 1�; RÞ. Since we are only considering linear first-order systems, the operators take the following form: As shown by Villegas [14,Chapter 6.3], the operator J À G R SG � R in Equation (8) is equivalent to the expanded skew-symmetric operator J e together with the closure relation e r ¼ Sf r in the expanded system In this formulation, the Dirac structure induced by J e has an additional port, called the resistive port. This port is then terminated by the resistive closure relation (9b). Obviously, this also means that the boundary conditions can also depend on e r : Since W B only gives us n conditions, this leaves the other n open for use as outputs, so we can define the input u and output y as To ensure that we do not use quantities such as output that are already set by the input, and we require that W B W C � � is of full rank. If W B �W T B � 0 then the port-Hamiltonian system with dissipation is a boundary control system [14, Theorem 6.11] and thus the system possesses (mild and classical) solutions. However, there are no exponential stability results available for this class of port-Hamiltonian systems.

The PHS-formulation of the model system
In the following, we will first model each part of our model system as a port-Hamiltonian system, and then look at the combined system in Section 4.3.

The heat equation
The heat Equation (2a) cannot be written as a port-Hamiltonian system of the form given in Definition 3.2. However, it can be written as a system with dissipation of the form given in Equation (8). In this work, we will derive a port-Hamiltonian formulation of the heat equation by choosing a suitable Hamiltonian and then comparing the heat equation with the structure given in Definition 3.2 and Equation (8) in order to find the correct operators J , G R , etc. As we will see, this approach can work quite well and even recover some of the physics underlying the heat equation. The opposite approach, starting with the physical laws and showing that they form a Stokes-Dirac structure is also possible as seen, for example, in Section 2.3.3 of [19].
By choosing the associated quadratic Hamiltonian as well as the usual choice of flow variable f ¼ @# m @t and effort variable e ¼ H# m , we find that with l m @ � and thus we recover the heat Equation (2a) from Equation (8).
Note that in this case the Hamiltonian (12) is not the physical energy so the dissipation present in this system does not automatically violate the law of conservation of energy. This is by design, since choosing this Hamiltonian allows us to use both the temperature and heat flux in the boundary port. It would, of course, be possible to use the physical energy as the Hamiltonian instead. However, this would create several issues. For one, we would have to introduce the thermodynamic entropy into the system formulation, either replacing the state variable, or as a function of the state variable. In this case, the irreversibility of the underlying physical system would need to be handled explicitly, instead of hiding it in the dissipative term. The choice of using a Hamiltonian other than the physical energy is therefore a practical choice for analysis since it avoids the need for a framework for irreversible port-Hamiltonian systems, cf [17,20]. On top of that, this change would also change the quantities involved in the boundary portfor example, replacing the heat flux with the entropy flux. While the issues caused with regard to the coupling conditions would be surmountable, they would require replacing the transport equations with some more complex models, as, to our knowledge, the transport equations are incompatible with this choice of Hamiltonian (and non-quadratic Hamiltonians in general).
By putting the heat equation into the form of (9), we obtain the following equations This formulation makes it clear how spatial derivatives of # can occur in the boundary conditions. Remembering that f ¼ @# m @t and e ¼ H# m , it follows that e r ¼ k c m l m @# m @� . We show that this port-Hamiltonian formulation organically gives rise to some physically meaningful quantities. For example, f r is a rescaled temperature gradient, causing the heat flux and e r is a rescaled heat flux, driving temperature change over time so that Equation (13b) becomes Fourier's law of thermal conduction.
The corresponding boundary conditions (2d), (2e) can now be rewritten in a formulation with inputs: With and therefore We can see that W B is of full rank. As h a , h b and l m are all physical constants of the system and thus positive, we find that W B �W T B > 0 holds. Therefore, the heat equation with the chosen boundary conditions has unique (classical and mild) solutions, which are nonincreasing in norm, cf. [14,Theorem 6.9].
As outputs, we choose the temperatures at the boundaries (scaled by some constants for convenience), so we get Then the combined matrix e w B e W C � � has full rank, which means we are not measuring quantities we already set as input.

The transport equations
We divide the cooling channel into two parts -an incoming and an outgoing channel, denoted by the indices in and out respectively -and model each of these sections separately as a transport equation. Since the transport equation does not have any dissipative terms, we can write each one directly as a linear first-order port-Hamiltonian system as defined in Definition 3.2. Since the resulting port-Hamiltonian systems are quite simple, we have combined them into a single system for brevity of notation.
@# in @t @# out @t i.e. with the choices Next, rewriting the boundary conditions (2f) and (2g) in a formulation with inputs, we obtain |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } As outputs, we choose the temperature at the ends of each cooling channel part, i.e.  (7), we obtain (21) From the last line we can see that this system is not stable for all variable choices since the stability condition W B �W T B � 0 is not always satisfied. Positive semi-definiteness is only given for In the special case of l in ¼ l out , this condition can be simplified to 2c This condition links the fluid's capacity to transport heat and the boundary's heat transfer. This is most likely connected to the (unphysical) temperature discontinuity at the coupling point contained in our model. The magnitude of that temperature jump depends on h b , c c and v. As the real physical system would not have a temperature discontinuity, we assume that this condition puts a bound on how far we can deviate from physics before the model starts to break down completely. Mathematically, violating this condition means we can no longer guarantee existence and uniqueness of the solution.

The combined system
Since we now have a port-Hamiltonian model for all three sub-systems, we can now consider the interconnected system. Remembering Equations (1) and (2), there are two coupling conditions. Equation (1e) couples the heat equation to the incoming cooling channel, relating heat flux and temperature difference similar to Newton's law of cooling. As a second condition, we have either Equation (1g) or (1h). Both equations couple the outflowing cooling channel to the other two sub-systems in a way that ensures physical energy conservation. Their difference lies in the choice of representing the energy flowing out of the heat equation. Note that while the coupling conserves physical energy, it is not necessarily conservative with regard to the Hamiltonians chosen in the previous sections. It is therefore not clear, a priori whether the combined system is a port-Hamiltonian system. In this section, we investigate whether the coupled system is port-Hamiltonian. When we combine the three port-Hamiltonian sub-systems discussed before, we get the following system of equations. As we can see in Equation (26), two boundary conditions are not connected to an input anymore, but are instead set to zero. These are the aforementioned coupling conditions between the two subsystems. Accordingly, we also have only two outputs. Specifically, the system has the form given by Equations (9) resp. (8), with the choices of Inserting these choices into Equation (9) we get the following form for the combined system: We can see from Equation (29) that both eigenvalues are non-negative if the following two conditions hold: Since the boundary condition (1h) is equivalent to (1g) due to (1e), we can use that one to develop our coupling structure instead. In doing so, we find: (30) Here it becomes immediately clear that the matrix for l in ¼ l out can never be positive definite, but only semi-definite, if 2h b l m À m v � 0 and À c c h b l in c m þ c m l m h b ¼ 0 holds, which yields 2c c v � h b , the same restriction as seen in Section 4.2. Thus, it is more restrictive in this respect than the previous coupling.
Summarizing, the combined system possesses unique (classical and mild) solutions on ½0; 1Þ and these solution are non-increasing in the energy norm. We remark, that it is an open question whether exponential stability of the extended port-Hamiltonian system implies exponential stability of the combined system, that is, whether the condition W B �W T B > 0 guarantees exponential stability of the combined system.

Conclusion
While the heat equation as described in Section 4.1 is exponentially stable for all (physically meaningful, i.e. positive) values of the constants involved, this is not the case for the transport equations of the cooling channel described in Section 4.2 nor for the combined system of Section 4.3. For both systems, it is possible to find values of the constants that are physically reasonable but do not satisfy the stability criteria, regardless of which formulation is chosen for the coupling conditions. It is also noteworthy that the two formulations of the coupling conditions studied, while technically equivalent, imply different regions of stability. It is still an open question as to why this is the case. Although some of the stability conditions have a clear, physical motivation, others -particularly those related to the coupling condition -appear to be nonsensical. The most obvious example of the latter would be that the ratio between the lengths of the cooling channel parts can determine the stability of the system. As there is no feedback mechanism from the end of the outflowing cooling channel to the rest of the system, there also cannot be a physically sound explanation for this behaviour.
Therefore, the most likely explanation or interpretation for this is that our model system is oversimplified and does not properly capture the properties of the real system. The fact that all heat transfer occurs at a single point as opposed to an extended region with a non-zero physical dimension, causes a discontinuity in the temperature distribution within the cooling channel, making it a likely candidate for the source of the above problems. In the future work, we will investigate whether the observed stability problems persist if we consider instead the coupling of a twodimensional heat equation with heat transfer along the entire length of the cooling channel.
In any case, the port-Hamiltonian formalism has proven to be a very useful tool to study the properties of a model system describing several interconnected physical processes and to find problems and limitations of the model. With its help, it has been shown that simplifications -even those widely used in industry -are not always useful from a mathematical point of view and must be carefully evaluated before use.
Our future work will focus on investigating new discretization techniques based on these coupling strategies and studying a more realistic two-dimensional setting.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Nomenclature/Notation
Th naming of variables denoting physical quantities follows the conventions set in Table 1. Variables that do not directly referencing physical quantities are generally defined when introduced.
The index m is used for variables that refer to the metal subsystem, i.e. the heat conducting rod. The indices in and out are used for the inflowing and outflowing parts of the cooling channel, respectively. Finally, the index c is used for variables that refer to the entirety of the cooling channel.