Radial Turbine Thermo-Mechanical Stress Optimization by Multidisciplinary Discrete Adjoint Method

Abstract

This paper addresses the problem of the design optimization of turbomachinery components under thermo-mechanical constraints, with focus on a radial turbine impeller for turbocharger applications. Typically, turbine components operate at high temperatures and are exposed to important thermal gradients, leading to thermal stresses. Dealing with such structural requirements necessitates the optimization algorithms to operate a coupling between fluid and structural solvers that is computationally intensive. To reduce the cost during the optimization, a novel multiphysics gradient-based approach is developed in this work, integrating a Conjugate Heat Transfer procedure by means of a partitioned coupling technique. The discrete adjoint framework allows for the efficient computation of the gradients of the thermo-mechanical constraint with respect to a large number of design variables. The contribution of the thermal strains to the sensitivities of the cost function extends the multidisciplinary outlook of the optimization and the accuracy of its predictions, with the aim of reducing the empirical safety factors applied to the design process. Finally, a turbine impeller is analyzed in a demanding operative condition and the gradient information results in a perturbation of the grid coordinates, reducing the stresses at the rotor back-plate, as a demonstration of the suitability of the presented method.

1. Introduction

Multidisciplinary optimization methods are widely adopted in the development cycle of turbomachinery-related technologies, such as in the energy and mobility businesses [1,2,3]. The increasing complexity of the products in terms of geometry and flow characteristics in design and off-design conditions, along with the fulfillment of requirements of lifetime, cost, and packaging, is reflected in the growth of concurrent design techniques that provide the experts with a comprehensive view on the problem [4,5]. The exploration of the design space is enriched by multiple physical disciplines whose interactions are captured along the optimization, resulting in accurate predictions of the component performance and any possible violation of the structural constraints. This holistic approach to the design problem differs from classical staggered methods, presenting a series of mono-discipline optimizations performed in cascade, leading typically to significantly longer development times. Moreover, such staggered methods may not converge to the same optimum as the holistic one and thus lead to suboptimal designs.

1.1. Background

The recent evolution of the market request towards aggressive reductions in emissions and fuel consumption offers new challenges to the research community. The need for further advances in multidisciplinary design optimization techniques is presented in the NASA 2040 roadmap [6] as one of nine core development areas in the context of integrated multiscale modelling frameworks contributing to the success in efficient design, manufacturing, and certification of future aerospace systems. Similarly, the engineering of advanced turbocharging solutions in the automotive sector demands further increases in efficiency to support novel combustion strategies, with no compromise on the durability targets [7,8].

Enhancements in turbomachinery components performance can be achieved by increasing the detail in the geometrical representation of the domain and reflecting it in a higher number of parameters engaged during the optimization process. This higher resolution allows the exploration of more complex shapes and discloses hidden interactions among features now present in the design space. Additionally, an increased fidelity in the physical description of the multidisciplinary problem is highly regarded and the aim is twofold. First, the optimizer is able to make more reliable predictions of the responses to the geometrical modifications, having access to lumped information accounting for the interactions of more disciplines involved in the computation of the objective functions and constraints. Therefore, the evolution of the optimization benefits from adhering more to the real behavior of the component, finally resulting in a reduced amount of experimental validation testing. Second, a more sophisticated physical description opens the path to adopting a richer set of constraints, and the safety margins usually applied to the design process can be relaxed in favor of more performance-oriented layout choices. Therefore, the optimizer can target more diversified paths towards the achievement of the desired solution, thus reducing the stiffness of the algorithm in performing the shape modifications.

The increased complexity in the models is reflected in the choice of a suitable optimization method between two major classes, notably denoted as gradient-free and gradient-based techniques.

Gradient-free methods [9], also known as “zero-order” methods, solely make use of the function values for the search of the global optimum, neglecting any knowledge of its gradients. The solution is identified after the evaluation of the output values of a large population of candidates with characteristics scattered within the boundaries of the selected design space. The larger the population, the higher the likelihood of identifying the global optimum. In the case of a multidisciplinary problem, each sample in the population must be processed through all the disciplines (FEM and FVM analyses, among others) [10]. Therefore, the number of computations necessary to identify a solution scales up with the number of involved disciplines and the number of design variables necessary to investigate the full domain, resulting quickly in a computationally expensive approach [11,12]. Even if, nowadays, gradient-free methods still represent the state-of-the-art approach to optimizations in industry because of their robustness and the high quality of the solutions, the advent of more stringent requirements and the need of evaluating larger domains are exposing an intrinsic limitation, notably referred to as the “curse of dimensionality” [13].

Gradient-based methods [14,15], with specific reference to “first-order” methods, make use of the gradient information of the function of interest in the search of the optimal solution, and therefore are more efficient than zero-order methods. The evaluation of the gradients may be computationally intensive if performed by Finite Differences or complex step techniques [16], especially in the presence of a large number of design variables and in the case of disciplines involving expensive computations. However, if the number of objective functions is lower than the number of optimization parameters, an efficient calculation of the gradients can be obtained by the adjoint method [17]. This technique, first introduced in [18] for fluid problems, is characterized by a cost for the computation of the sensitivity derivatives of the objective function that is essentially independent of the number of design variables. Therefore, it is particularly suited for problems requiring detailed geometrical descriptions with many degrees of freedom. Moreover, the exploration of the design space through its gradients decreases the number of iterations necessary to achieve the local optimal solution. Therefore, the reduced overhead can be traded for the introduction of more computationally demanding disciplines, returning a higher accuracy in the numerical solutions.

A performance comparison between these two techniques is offered in [19], whose application to the multi-point optimization of a radial turbine impeller shows that the adjoint method is capable of extending the Pareto front identified by a differential evolution (DE) algorithm. In particular, the gradient-based approach delivers a solution with a higher total-to-static efficiency, improving also the moment of inertia figure at the cost of a fraction of the computational time required by the gradient-free method to achieve its global optimum (roughly 6–10%). Other examples of applications of the adjoint methods to the optimization of turbines and compressors are reported in [20,21,22,23,24,25], with focus on the opportunity of increasing the number of variables in the design space and the multidisciplinary aspect of the physical evaluations (including aerodynamics, stresses, and modal analysis).

The higher maturity recently gained by gradient-based approaches and their improved robustness in the analysis of complex flow conditions [26,27] provides evidence of their application to problems of industrial relevance. Therefore, adjoint methods are considered a suitable choice to address the previously mentioned design advancements in turbomachinery-related applications.

1.2. Motivation

The present document focuses on the development of a radial turbine impeller for an automotive turbocharger. In contrast to aeronautical applications, the typical duty cycles for passenger cars span the entire turbine map, from very low mass flow rates and pressure ratios in urban driving situations up to peak power, resulting in requirements stretching the machine design over a wide range of operative conditions. Moreover, these turbines present tight clearances between rotor and housing in the order of few hundreds of microns, whilst the typical thermal operative range spans from ambient temperature up to 1050 °C, without the possibility of any blade-dedicated cooling circuit. Such scenario clearly demonstrates the challenges in performance requirements and durability the designer is confronted with during the product development process.

The structural integrity of the rotor experiencing centrifugal and fluid-induced stresses is a critical aspect in the search for design solutions fulfilling the efficiency and permeability targets. The problem of the Fluid–Structure Interaction (FSI) in radial turbines in relation to the thermal stresses induced by the fluid in contact with the blades is discussed in [28], in an effort to accurately predict the heat exchange between the two media for a reliable estimation of the rotor lifetime in steady state and dynamic operative conditions. The study shows how the incorrect estimation of the thermal stresses could impair the robustness of the component, justifying the inclusion of this additional constraint on top of the computation of the centrifugal stresses.

The introduction of thermal analyses in the framework of an optimization process is a topic of recent interest [29,30,31]. The physical phenomenon is described through the application of a Conjugate Heat Transfer (CHT) procedure, as originally proposed in [32,33], which involves the thermal interaction between a solid body and its surrounding fluid by a coupled solution of the two domains. The challenge in the implementation of this analysis within the landscape of a multidisciplinary optimization of a complex geometry consists of the trade-off between accuracy and significant computational overhead, so far limiting thermal predictions only to the final validation phase in industrial design procedures. In fact, it is common practice during an optimization to replace such intensive computations with reduced-order models or safety margins based on empirical experience. This approach to the evaluation of the thermal stresses would suffer in terms of accuracy or, in the case of large safety margins, would over-constrain the optimization problem.

1.3. Goal of the Paper

This paper addresses the problem of robust optimizations w.r.t. thermo-mechanical stresses by developing a discrete adjoint framework tailored to the implementation of the CHT analysis within the in-house multidisciplinary design and optimization platform “CADO” developed at the von Karman Institute for Fluid Dynamics [34]. Indeed, the selection of a gradient-based method is justified by its computational efficiency in identifying a local minimum, although the computation of the gradients is more convoluted due to the multiple disciplines involved in the performance analysis. The calculation of the thermal stresses in the solid domain and their sensitivities with respect to the grid coordinates is therefore the subject of the present paper, with the aim of including this analysis in a Sequential Quadratic Programming-based optimizer [35].

To the best of our knowledge, only a few studies have considered the CHT analysis in adjoint-based optimizations [36,37], with only a limited number adopting a discrete adjoint formulation [38,39,40] and none treating a test case with a complex three-dimensional flow field or extending the problem to thermo-mechanical evaluations. Therefore, this work is aimed at covering critical gaps in the design of radial turbines with the aim of creating a framework suitable for industrial applications.

In light of this objective, the remaining portion of the document is structured as follows. First, the framework of the primal solver for the execution of the CHT analysis is presented, with extension to the computation of the thermo-mechanical stresses. Then, the implementation of the adjoint workflow for the calculation of the sensitivities of the constraint function with respect to the solid and fluid grid coordinates is discussed, outlining the steps towards the coupling of the domains in reverse mode. Finally, the validation of the model and the results of the application of the procedure to a radial turbine rotor test case are presented and the conclusions are drawn from this study.

2. Primal Solver

The CHT analysis implemented in the current work follows the partitioned coupling approach proposed by [41], in which a fluid and a solid solver are sequentially called in an iterative process with the mutual exchange of boundary conditions till reaching the convergence of the temperature at their interface. The advantage of this strategy is related to the possibility of adopting dedicated meshes and specialized numerical methods for the solution of each domain, increasing the robustness of the convergence of the two analyses. Moreover, in the case of an unsteady computation, the partial decoupling of the two fields, which show very different characteristic time scales, would lead to a faster convergence of the assembly [42,43,44,45].

2.1. Mapping Procedure

The problem under investigation is represented in Figure 1, in which a turbine rotor is analyzed considering a periodic sector. A multi-block structured fluid grid with boundary layer refinement is interfaced to an unstructured solid mesh of second-order tetrahedral elements for the computation of the heat flux between the two domains.

Because the meshes are non-matching at the common interface, the exchange of information is performed through the distance-weighted interpolation (DWI) technique presented in [46]. The procedure implies the mapping of the two regions in order to record the correspondences of each fluid cell center with the closest surrounding solid nodes and vice versa. For each fluid cell, the information of heat flux or temperature detected in the mating solid nodes is processed by a weighting function based on their mutual distances, resulting in the coupling of the two regions. The same routine is applied to the solid grid with respect to its fluid counterpart.

This method exhibits a good level of generality in treating different coupling problems, irrespective of their grids characteristics. For instance, in the turbine test case the exchange of information is effectively performed also at the blade hub fillet, a feature which is present in the solid domain but not detailed in the corresponding fluid region.

An improvement in the accuracy of this coupling procedure derives from the local labelling of the interface. With reference to Figure 2, the surface between the two media is split in sub-regions (e.g., hub, blade pressure side, blade suction side, tip, etc.) and the corresponding solid faces and fluid blocks boundaries are named accordingly. Therefore, the algorithm searching for the fluid cells–solid nodes correspondences is separately executed for each sub-region, instead of attempting the coupling of the two fields through a global search within the entire domain. This arrangement is particularly effective in locations presenting thin walls, such as in proximity of the blade tip, as it avoids any erroneous fluid-solid matching between cells and nodes actually sitting on opposite sides of the blade (i.e., suction side versus pressure side).

Additionally, the original algorithm seeks the nodes–cells correspondences through a “search radius” delimiting for each grid node a spherical volume within which the neighbors are identified. This technique may induce interpolation errors in proximity of the walls, where temperature gradients and fluid grid stretching are significant. In fact, the temperature assigned to a solid node may result from the weighting of temperature values derived from fluid cells located within the boundary layer but not in direct contact with the interface. Moreover, as demonstrated in [46], distinct grid refinements at the solid and fluid side may influence the local accuracy of the wall temperature information passed to the FEM. In order to cope with these issues, the searching criterion is limited to the sole layer of nodes and cells in direct contact with the walls, and the concept of “search radius” is dismissed in favor of a ranking assigned to each node or cell according to their mutual distance. Finally, consistently with the solution proposed in [46] addressing the problem of the aspect ratio of cells on a wall boundary, additional virtual grid points are introduced on the FVM side between two neighboring cell centers. Such virtual points are aimed at improving the uniformity of the local distribution of grid points in regions of significant mesh stretching, and their temperature is interpolated between the closest “master” cells of the original mesh. Indeed, the goal is to increase the density of the attractors in proximity of each solid node, avoiding possible inconsistencies generated from capturing the information from too far fluid cells. Thus, the domain coupling is executed through the distance-weighting interpolation of information belonging to cells and nodes in actual close proximity.

The quality of the interpolation is referenced in Figure 3. The fluid temperature at the walls (Figure 3a) is interpolated by the DWI procedure and passed to the solid nodes laying on the surface (Figure 3b). During a CHT iteration, the solid returns the heat flux to the fluid, as discussed in the next section.

2.2. Partitioned Coupling

The data exchange between the fluid and solid fields is performed through the “heat transfer forward flux back” method (hFFB) [47], whose stability properties are discussed in [48]. In summary, the convergence rate of the conjugate problem depends on the local Biot number, which expresses the ratio between the conductive over the convective thermal resistances. In complex geometries, as in the case of a radial turbine impeller, the Biot number may locally change to values greater or lower than unity, according to the variations in the blade thickness. Hence, [46] demonstrates that methods such as the hFFB or the “heat transfer forward temperature back” (hFTB) method would promote the stabilization of the fluid-solid coupling w.r.t. other techniques, such as the “flux forward temperature back” (FFTB) or the “temperature forward flux back” (TFFB) methods, which strictly require a modulus of the Biot number respectively lower or greater than unity. Instead, the hFTB and hFFB methods exhibit a wider range of convergence. Finally, the hFFB method was selected for the present study since a heat flux boundary condition imposed to the fluid domain would generally improve the convergence stability of the CFD computation more than an imposed wall temperature boundary condition.

Figure 4 reports the partitioned coupling workflow, described as follows.

The CHT analysis is addressed with an initial CFD computation considering adiabatic walls. The fluid temperature $T_{FW}$ and the heat flux normal to the walls $q_{FW}$ are extracted at the boundaries in contact with the solid, and a virtual fluid bulk temperature $T_{fl}$ is calculated as:

$$T_{fl}=T_{FW}-\frac{q_{FW}}{\tilde{h}},$$

(1)

with $\tilde{h}$ indicating a constant user-defined virtual heat transfer coefficient.

Since the initial fluid simulation is adiabatic, Equation (1) results in a virtual fluid bulk temperature equal to the fluid temperature at the walls. The DWI procedure returns the $T_{fl}^{\prime }$ field associated with the solid grid nodes.

Consequently, $\tilde{h}$ and $T_{fl}^{\prime }$ are assigned to the FEM heat transfer model, such that the heat flux to the solid results from the following boundary condition:

$$q_{SW}=\tilde{h}\left(T-T_{fl}^{\prime }\right),$$

(2)

with T as the unknown solid temperature. After solving the FEM heat transfer problem, the temperature and heat flux at this boundary are known.

The heat flux at the interface is calculated from the temperature field T by means of the Fourier’s law applied to all the elements exposing at least one face to the fluid:

$$q_{SW}=-k\frac{dT}{dn},$$

(3)

with $dT/dn$ indicating the temperature gradient normal to the wall and k the thermal conductivity coefficient.

The heat flux $q_{SW}$ is processed by the DWI procedure, returning the heat flux $q_{FW}$ assigned to the fluid cells at the interface. The fluid simulation, now accounting for an external heat flux at the viscous walls, is recomputed and the entire process is re-iterated for several loops till the achievement of the continuity of temperatures at the interface between the two fields.

The user-imposed virtual heat transfer coefficient $\tilde{h}$ influences the predicted wall heat fluxes at any intermediate cycle of the CHT procedure, thus determining the path to convergence of the whole coupling. In the case of the hFFB method, [46] shows that a suitable choice of the virtual coefficient $\tilde{h}<2h$, with h as the physical value of the heat transfer coefficient, would guarantee the convergence of the conjugate problem for any local value of the Biot number. The higher the value of $\tilde{h}$ within the stable region, the faster the convergence of the partitioned coupling method. A discussion about the determination of the $\tilde{h}$ value for a problem of industrial relevance is reported in Section 4.2.

2.3. Solid Heat Transfer Solver

The steady-state energy balance within the solid domain is described by Equation (4), with q as the rate of heat transfer into the system, and is computed through a FEM steady linear solver according to [49].

$$\nabla ·q=0,$$

(4)

The boundary conditions applied to the problem are classified in three families:

$$\begin{array}{c}{T|}_{\partial {\Omega }_1}=const\\ -k\nabla {T|}_{\partial {\Omega }_2}=q\\ -k\nabla {T|}_{\partial {\Omega }_3}=\tilde{h}\left({T|}_{\partial {\Omega }_3}-T_{fl}^{\prime }\right)\end{array},$$

(5)

with the surface of the solid domain defined as $\partial {\Omega }_1{\displaystyle \cup }\partial {\Omega }_2{\displaystyle \cup }\partial {\Omega }_3=\partial \Omega$.

The first equation in (5) implies a fixed temperature imposed at the surface $\partial {\Omega }_1.$ The second type of boundary conditions specifies a heat flux at the surface $\partial {\Omega }_2$, while Fourier’s law at LHS describes the heat flux through the medium. If q = 0, the adiabatic wall boundary condition is defined. The third equation represents the convection condition, whose heat flux through the surface $\partial {\Omega }_3$ is proportional to the temperature difference w.r.t. the surrounding fluid. The virtual heat transfer coefficient $\tilde{h}$ and the interpolated thermal field $T_{fl}^{\prime }$ are reported for the sake of consistency with Equation (2).

Once integrated over the entire domain, Equation (4) is discretized, resulting in the linear system:

$$A\mathit{T}=b,$$

(6)

with the semi-positive definite stiffness matrix A comprising the conductive and convective terms, A = [A_cond + A_conv] $\in R^{n,n}$, and the load vector b $\in R^n$ accounting for the contribution of the boundary conditions. The linear system in (6) is solved by an iterative conjugate gradient method.

2.4. Fluid Solver

The compressible Reynolds-Averaged Navier–Stokes solver with cell-centered spatial discretization over a structured multi-block grid, developed in [19], is invoked in the CHT framework. The solver computes the convective fluxes by Roe’s upwind scheme with MUSCL extrapolation, while a central discretization is applied to the calculation of viscous fluxes. The time marching technique is implemented according to the JT-KIRK scheme proposed in [26] for the improved stabilization of the discrete adjoint solver. The Spalart–Allmaras turbulence model is adopted in this study.

In order to establish the coupling with the solid domain, heat fluxes are imposed at the viscous wall boundaries adopting a thin shear layer approximation. The fluid mesh presents two layers of ghost cells in order to facilitate the computation of the fluxes at the interfaces [50]. The heat fluxes calculated in Equation (3) and processed through the DWI technique are associated with the corresponding fluid block boundaries in contact with the solid walls. With reference to Figure 5, D1 is the first layer of the inner fluid domain; G1 and G2 are the corresponding ghost layers. The fluid temperature in D1 is detected from the previous solver iteration, and the ghost cell temperature in the first layer is updated according to:

$$T_{G1}=T_{D1}+q_{FW}\frac{\Delta n}{k_{fl}},$$

(7)

with $\Delta n$ as the distance between the cell centers and k_fl as the locally computed fluid conductivity. Similarly, the temperature in the G2 layer is updated in cascade, with reference to the newly computed temperature in G1. The temperature values are used to calculate the corresponding cell densities considering a zero-order pressure extrapolation in wall normal direction, while the no-slip condition is assured by reversing the fluid tangential velocity components in the ghost layers. Equation (7) is linearized and accounted in the implicit scheme.

The rate of convergence induced by the Neumann boundary condition is influenced by the actual magnitude of the imposed heat flux in Equation (7). Since the hFFB method is employed in the present development, no under-relaxation factors are invoked for the update of the ghost cells temperatures, as they are replaced by the selection of a suitable virtual heat transfer coefficient $\tilde{h}$, determinant for the stability of the coupling process [46].

2.5. Solid Mechanical Solver

Once the continuity of temperatures and heat flux is achieved at the interface of the two domains, an FEM linear elastic solver with quadratic elements presented in [51] is invoked for the solution of the von Mises stresses.

$$S\mathit{u}=f,$$

(8)

with $S\in R^{n,n}$ as the stiffness matrix, $\mathit{u}\in R^n$ referencing the vector of nodal displacements, and $f\in R^n$ indicating the load.

Within the framework of the present analysis, the solid mechanical solver shares the same mesh as the heat transfer solver, as remarked in Section 3. Therefore, the previously computed nodal temperatures are directly assigned to the new linear system.

The right-hand side vector accounts for the centrifugal loading and the thermal strains, considering the effect of the material contractions and expansions induced by the temperature variations, which are computed as:

$${\epsilon }_{thermal}=\alpha \left(T_i\right)\left(T_i-T_{ref}\right),$$

(9)

with α(T_i) indicating the temperature-dependent thermal expansion coefficient; T_i pointing to the node temperatures resulting from the CHT analysis; and T_ref as the reference temperature at the initial state, at which no thermal stresses are present within the material [52]. In the current study, we refer to T_ref as the ambient temperature. Through the thermal strains, the stresses in the material are dependent on the outcome of the CHT computation, which increases the complexity of the gradient calculation, as presented later on.

The maximum von Mises stress is computed using a p-norm function according to Equation (10), with p = 75.

$${\sigma }_{maxVM}\approx \sqrt[p]{\frac{{{\displaystyle \int }}_V^{}{\sigma }_{VM}^{p}dV}{{\displaystyle \int }dV}}.$$

(10)

Equation (10) fulfills the requirement of global differentiability, as it is a continuous function, and is invoked as a constraint within the optimization problem with the aim of keeping the thermo-mechanical stresses, induced through the Conjugate Heat Transfer, bounded to a prescribed value. The gradient of this constraint is computed by an adjoint methodology, as described in the next section.

3. Adjoint Solver

The structure of the primal solver evaluating the thermo-mechanical stresses is computationally expensive because the CHT procedure requires several loops of CFD-FEM analysis to reach the equilibrium of the temperatures and heat flux at the fluid–solid interface. For this reason, such a framework is less suited for gradient-free optimization methods, as a prohibitive large population of candidates would need to be evaluated using this expensive tool. A gradient-based optimization technique, on the other hand, may offer a faster convergence through a more guided search in the design space. Especially when combined with the adjoint method, the computation of the derivatives is very efficient and almost insensitive to the number of design parameters.

The evaluation of the cost function J in primal mode is a two-step process. The first one considers the values of the design variables $\alpha \in R^n$ and generates the correspondent geometry and finally the grid, whose points coordinates are expressed by the vector $X\in R^m$. The second step involves the flow solver and some post-processing routines, returning the value of J. Following the same approach in reverse mode, also the adjoint evaluation procedure is split in two parts: the first one computing the sensitivities of the cost function J w.r.t. the perturbations of the grid coordinates X, and the second one calculating the sensitivities of the grid coordinates X w.r.t the changes in the design variable α. Therefore, the CAD-based parametrization described in [21] allows transferring the grid sensitivities of the cost function to the design variables, controlling the component shape modifications by the application of the chain rule of differentiation expressed in Equation (11):

$$\frac{dJ}{d\alpha }=\frac{\partial J}{\partial X}\frac{dX}{d\alpha }.$$

(11)

Within the thermo-mechanical analysis process developed in the present work, Equation (11) is invoked twice in order to account for the cost function sensitivities w.r.t. the coordinates of the fluid grid points and the solid grid nodes, respectively. Therefore, the global sensitivity expression can be decomposed as follows:

$${\left(\frac{dJ}{d\alpha }\right)}_{TM}=\frac{\partial J}{\partial X_{fl}}\frac{d{X}_{fl}}{d\alpha }+\frac{\partial J}{\partial X_{sl}}\frac{d{X}_{sl}}{d\alpha }.$$

(12)

The remaining portion of this section is devoted to the calculation of the cost function sensitivities with respect to the fluid and solid grid coordinates, respectively $\partial J/\partial X_{fl}$ and $\partial J/\partial X_{sl}$, operated through the application of a discrete adjoint method. Such sensitivities provide the descent direction for constrained or unconstrained multidisciplinary optimization problems with the aim of minimizing the thermally related response functions.

Consistently with the structure of the primal solver, the adjoint framework treating this interdisciplinary problem is formulated according to a “loose-coupling” approach, different from the “strong-coupling” techniques as described in [53]. In fact, the latter accounts for cross-discipline Jacobian terms implicitly exchanging boundary conditions between the fluid and solid domains, and therefore directly solving the global adjoint system at once. Instead, the selected partitioned-coupling method is less intrusive in the structure of the existing adjoint CFD solver from [19] and therefore more suited to a continuously growing modular multidisciplinary platform.

Figure 6 reports the framework of the thermo-mechanical evaluation, with the adjoint workflow on the RHS. The iterative structure characterizing the solution of the primal problem is maintained in its adjoint counterpart, which walks through the entire chain in reverse mode by the manual differentiation of the full process. Therefore, the adjoint branch is discussed starting from the bottom, till achieving the final grid sensitivities on the top.

3.1. Adjoint Variables

The most demanding terms in the computation of Equation (12) are the sensitivities of the response function w.r.t. the grid coordinates. Since:

$$J\left(X,u\left(X\right)\right),$$

(13)

with u denoting the vector of the state variables, the application of the chain rule of differentiation results in the following decomposition:

$$\frac{dJ}{dX}=\frac{\partial J}{\partial X}+\frac{\partial J}{\partial u}\frac{\partial u}{\partial X}.$$

(14)

Equation (14) holds for both the fluid and the solid grids. The second term at RHS presents the partial derivative of the cost function w.r.t. the state variables $\partial J/\partial u$. Since the current problem exhibits a multidisciplinary footprint, the hierarchical structure of the workflow in Figure 6 implies that the output states contributing to the p-norm calculation in Equation (10) depend on the intermediate ones evaluated by the preceding CFD and FEM solvers. Therefore, according to the principle of the “reverse differentiation” [54] we choose for instance an output variable $u_3$, whose sensitivity $dJ/d{u}_3$ is known, and calculate its sensitivities w.r.t. each intermediate state till the initial one, such that:

$$J\left(u_3\left(u_2\left(u_1\right)\right)\right)\hspace{1em}\hspace{1em}\to \hspace{1em}\hspace{1em}\begin{array}{c}\frac{dJ}{d{u}_2}=\frac{dJ}{d{u}_3}\frac{d{u}_3}{d{u}_2}\\ \frac{dJ}{d{u}_1}=\frac{dJ}{d{u}_2}\frac{d{u}_2}{d{u}_1}\end{array}.$$

(15)

Since the work focuses on the single thermo-mechanical output $J_{TM}$ expressed by Equation (10), herein named y, we can associate with any intermediate state variable $u_i$ a new variable $\overline{u_i}$, called the adjoint variable, defined as:

$$\overline{u_i}=\frac{\partial y}{\partial u_i}.$$

(16)

The upper bar notation follows the convention for adjoint variables reported in [54].

Hence, we apply the chain rule walking throughout the original trace in Figure 6 in backward mode, paying attention to the opposite propagation of the adjoint variables with respect to the physical ones. This technique finally returns the sensitivities of the response function w.r.t. the initial state variables.

The following sections describe the backward propagation process in detail, with the aim of offering a guidance to the development of the adjoint method.

3.2. Adjoint Response Function

The adjoint framework is initiated by seeding the maximum von Mises constraint. Since $y=J_{TM}={\sigma }_{maxVM}$, the input is:

$$\frac{\partial y}{\partial J_{TM}}=\overline{J_{TM}}={\overline{\sigma }}_{\max VM}=1.$$

(17)

The maximum von Mises stress is function of the components of the Cauchy stress tensor. Therefore, in reverse mode the adjoint stress components are calculated as:

$$\frac{dy}{d{\sigma }_{i,j}}={\overline{\sigma }}_{i,j}=\frac{d{J}_{TM}}{d\sigma }\overline{J_{TM}}=\frac{d}{d\sigma }\sqrt[p]{\frac{{{\displaystyle \int }}_V^{}{\sigma }_{VM}^{p}dV}{{\displaystyle \int }dV}}.$$

(18)

On the other hand, Equation (10) shows an explicit dependence over the solid grid coordinates through the two volume integrals appearing at the numerator and denominator. Therefore, the constraint function sensitivities with respect to the solid grid coordinates are accumulated as follows:

$$\frac{dy}{d{X}_{sl}}=\overline{X_{sl}}+=\frac{d{J}_{TM}}{d{X}_{sl}}\overline{J_{TM}}=\frac{d}{d{X}_{sl}}\sqrt[p]{\frac{{{\displaystyle \int }}_V^{}{\sigma }_{VM}^{p}dV}{{\displaystyle \int }dV}}.$$

(19)

In Equation (19), the operator “+=” is reported to emphasize the concept of accumulation of the sensitivities w.r.t. the solid grid coordinates that is propagated backwards throughout the adjoint workflow till the grid generation procedure. Moreover, the $\overline{X}$ notation is adopted for the sake of simplicity, while in reality the sensitivities over the three dimensions (X,Y,Z) are accounted in the development.

Finally, in the remaining portion of the document the notations $dy/d{u}_i$ and $dy/dX$ will be skipped in favor of the more convenient ${\overline{u}}_i$ and $\overline{X}$.

3.3. Adjoint Mechanical Solver

The evaluations through the mechanical solver involve the multi-step approach outlined in Figure 7, where the primal mode is represented along with its adjoint counterpart. Notably, the reverse computation follows the backward propagation of the adjoint variables, closely resembling the primal scheme. The related procedure is discussed in detail in this section.

In primal mode, the stress components ${\sigma }_{i,j}$ are computed from the elasticity matrix of the material under consideration and the strain components. Thus, the differentiation in reverse mode returns the adjoint strains ${\overline{\epsilon }}_{i,j}$.

$${\sigma }_{i,j}=E\left(T\right)\left({\epsilon }_{m\_ij}-{\epsilon }_{th\_ij}\left(T\right)\right)\hspace{1em}\hspace{1em}\to \hspace{1em}\hspace{1em}\{\begin{array}{c}{\overline{\epsilon }}_{m\_ij}={\overline{\sigma }}_{i,j}E\left(T\right)\\ {\overline{\epsilon }}_{th\_ij}\left(T\right)=-{\overline{\sigma }}_{i,j}E\left(T\right)\\ \overline{E}\left(T\right)={\overline{\sigma }}_{i,j}\left({\epsilon }_{m\_ij}-{\epsilon }_{th\_ij}\left(T\right)\right)\end{array}.$$

(20)

Here, the adjoint stresses ${\overline{\sigma }}_{i,j}$ are the once computed by the reverse differentiation discussed in Section 3.2.

The adjoint mechanical strains ${\overline{\epsilon }}_m$ and the adjoint thermal strains ${\overline{\epsilon }}_{th}$ appear in Equation (20). Since the thermal strains exhibit direct dependence from the nodal temperatures through the $\left(T_i-T_{ref}\right)$ term and indirect dependence through the thermal expansion coefficient α(T), the relative contributions to the adjoint solid node temperatures are accumulated and propagated to the next step in the reverse workflow.

$${\epsilon }_{th\_ij}\left(T\right)=\alpha \left(T\right)\left(T-T_{ref}\right)\hspace{1em}\hspace{1em}\to \hspace{1em}\hspace{1em}\{\begin{array}{c}\overline{T}+={\overline{\epsilon }}_{th\_ij}\alpha \left(T\right)\\ \overline{T}+={\overline{\epsilon }}_{th\_ij}\left(T-T_{ref}\right)\frac{d\alpha \left(T\right)}{dT}\end{array}.$$

(21)

Additionally, also the elasticity matrix E depends on the nodal temperature, which leads similarly to contributions for the adjoint solid node temperatures.

$$E\left(T\right)=f\left(T\right)\hspace{1em}\hspace{1em}\to \hspace{1em}\hspace{1em}\overline{T}+=\overline{E}\left(T\right)\frac{df\left(T\right)}{dT}.$$

(22)

The process of accumulation of the contributes to the adjoint temperatures field is facilitated by the commonality of the solid grid coordinates between the mechanical and heat transfer solvers, as no re-interpolations are necessary for the backward propagation of the sensitivities.

Next, recalling that the mechanical strains $\epsilon$ in primal mode are calculated from the vector of nodal displacements u, the adjoint nodal displacements are derived by the reverse differentiation ($\overline{\epsilon }\to \overline{u}$) and provided in input to the adjoint mechanical solver.

The iterative linear system solver adopted in primal mode for the solution of Equation (8) is not directly differentiated but, as discussed in [20], a new linear system is formulated with the adjoint displacements vector appearing at RHS as:

$$S^T\overline{f}=\overline{u}.$$

(23)

The system (23) is solved for the adjoint load vector $\overline{f}$, while the transpose of the stiffness matrix S, previously computed in primal mode, appears at LHS. Therefore, the adjoint stiffness matrix is obtained from the system (24).

$${\overline{S}}_{ij}=-u_j{\overline{f}}_i.$$

(24)

Finally, the system assembly process is algorithmically differentiated, returning the grid sensitivity accumulation:

$${\overline{X}}_{sl}+=\overline{S}\frac{dS}{d{X}_{sl}}+\overline{f}\frac{df}{d{X}_{sl}}.$$

(25)

The system (25) accumulates the grid sensitivities attributed to the mechanical solver. The adjoint temperatures deriving from Equations (21) and (22) are provided as input to the adjoint heat transfer solver at the next section, thus linking the structural solver to the iterative adjoint CHT process.

3.4. Adjoint Heat Transfer Solver

Consistently with the previous section, the iterative linear system solver adopted in primal mode for the solution of Equation (6) is not directly differentiated, whilst the adjoint load vector and the adjoint stiffness matrix are computed through the new linear systems in Equation (26).

$$A^T\overline{b}=\overline{T}.$$

(26a)

$${\overline{A}}_{ij}=-T_j{\overline{b}}_i.$$

(26b)

Equation (26a) shows the transposed stiffness matrix $A^T$ accounting for the conductive and convective terms previously calculated in primal mode, while the adjoint temperature vector $\overline{T}$ from Section 3.3 is imposed at the RHS. The resulting adjoint load vector $\overline{b}$ is assigned to Equation (26b), together with the temperature solution stored at the last CHT loop in primal mode, contributing to the calculation of the adjoint stiffness matrix $\overline{A}$.

Similarly to Equation (25), the system assembly process is algorithmically differentiated and the adjoint contributes to the solid grid coordinates are accumulated. Since the thermal solver and the mechanical solver share the same grid, a unique vector of solid grid sensitivities is accumulated and propagated backwards throughout the CHT workflow.

The convective loading at the RHS of the system (6) accounts for the virtual bulk fluid temperature. Its adjoint counterpart ${\overline{T\prime }}_{fl}$ is computed from $\overline{b}$ (returned by Equation (26a)) and from the user-defined virtual heat transfer coefficient $\tilde{h}$.

$$b_i=\tilde{h}T{\prime }_{f{l}_i}\hspace{1em}\hspace{1em}\to \hspace{1em}\hspace{1em}{\overline{T\prime }}_{f{l}_i}={\overline{b}}_i\tilde{h}.$$

(27)

${\overline{T\prime }}_{fl}$ is not accumulated along the entire adjoint workflow but recomputed at every CHT loop in reverse mode and then passed to the adjoint hFFB routine.

3.5. Adjoint hFFB Procedure (Solid → Fluid)

The adjoint virtual bulk fluid temperature field ${\overline{T\prime }}_{fl}$ associated with the solid nodes is processed by algorithmic differentiation, returning the adjoint ${\overline{T}}_{fl}$ assigned to the neighboring fluid cells centers. Equation (28) reports the DWI procedure in primal mode on the left, with focus on the calculation of the solid temperature at node j, leveraging the information of the cluster of i = (1, ..., n) neighboring fluid cells. On the right, the backward differentiation of the interpolation procedure results in the adjoint temperature field associated to the fluid domain. The operator “+=” stresses the fact each fluid cell may belong to the “neighboring cluster” of different solid nodes j, and therefore would accumulate their adjoint values.

$${T^{\prime }}_{fl}\left(j\right)=\frac{{{\displaystyle \sum }}_{i=1}^n\frac{T_{fl}\left(i\right)}{{\left(dist\left(i\right)\right)}^2}}{{{\displaystyle \sum }}_{i=1}^n\frac{1}{{\left(dist\left(i\right)\right)}^2}}\hspace{1em}\hspace{1em}\to \hspace{1em}\hspace{1em}{\overline{T}}_{fl}\left(i\right)+={\overline{T^{\prime }}}_{fl}\left(j\right)\frac{\frac{1}{{\left(dist\left(i\right)\right)}^2}}{{{\displaystyle \sum }}_{i=1}^n\frac{1}{{\left(dist\left(i\right)\right)}^2}}.$$

(28)

Equation (28) suggests also that additional contributions to the solid grid sensitivities ${\overline{X}}_{sl}$ and to the fluid grid sensitivities ${\overline{X}}_{fl}$ can be accumulated through the adjoint vector of mapped distances between each solid node and the corresponding fluid neighbors, $\overline{dist}\left(i\right)$, as detailed in Section 2.1. Thus, after deriving the adjoint distances $\overline{dist}\left(i\right)$ as ${\overline{T^{\prime }}}_{fl}\left(j\right)·d{T}_{fl}^{\prime }\left(j\right)/d\left(dist\left(i\right)\right)$, the cells-nodes mapping routine is differentiated and the contributions to the two grids are obtained.

Finally, the algorithmic differentiation of Equation (1) returns the adjoint temperature at the fluid walls ${\overline{T}}_{FW}$, which is used as input to the adjoint CFD solver at the next step, and the adjoint wall heat flux ${\overline{q}}_{FW}$, stored for a later accumulation.

$$T_{fl}=T_{FW}-\frac{q_{FW}}{\tilde{h}}\hspace{1em}\hspace{1em}\to \hspace{1em}\hspace{1em}\{\begin{array}{c}{\overline{T}}_{FW}={\overline{T}}_{fl}\\ {\overline{q}}_{FW}={\overline{T}}_{fl}·\left(-\frac{1}{\tilde{h}}\right)\end{array}.$$

(29)

3.6. Adjoint Fluid Solver

The partitioned coupling technique adopted in this work allows for non-intrusive calls to the in-house adjoint CFD solver. The flow field is initialized with the converged flow solution U from the corresponding CHT loop in primal mode, and the viscous wall temperatures in the ghost cells layers are updated again imposing the Neumann boundary condition, as reported in Equation (7).

The adjoint fluid wall temperature ${\overline{T}}_{FW}$, resulting from the reverse hFFB procedure at the previous step, contributes to the linearization of the constraint function with respect to the conservative flow variables U and the fluid grid coordinates $X_{fl}$. The adjoint equation discussed in [19] is repeated here for convenience:

$${\left[\frac{\partial R}{\partial U}\right]}^T\psi ={\left[\frac{\partial J}{\partial U}\right]}^T,$$

(30)

with R as the non-linear residuals of the primal flow solver, ψ indicating the unknown adjoint variables, and J the objective or constraint function of interest. Therefore, the RHS term linearized at the viscous walls exposed to the heat flux from the solid develops by the chain rule as follows:

$$\begin{array}{c}\frac{\partial J}{\partial U}=\frac{\partial J}{\partial V_{bnd}}\frac{\partial V_{bnd}}{\partial V_{domain}}\frac{\partial V_{domain}}{\partial U_{domain}}=\frac{\partial J}{\partial T_{FW}}\frac{\partial T_{FW}}{\partial V_{bnd}}\frac{\partial V_{bnd}}{\partial V_{domain}}\frac{\partial V_{domain}}{\partial U_{domain}}=\\ {\overline{T}}_{FW}\frac{\partial T_{FW}}{\partial V_{bnd}}\frac{\partial V_{bnd}}{\partial V_{domain}}\frac{\partial V_{domain}}{\partial U_{domain}},\end{array}$$

(31)

with V_bnd as the primitive variables computed at the boundary of interest, V_domain as the corresponding variables from the interior of the fluid domain, and the last term $\partial V/\partial U$ indicating the transformation matrix from primitive to conservative variables.

Similarly, the linearization of the cost function w.r.t. the fluid grid coordinates is the following:

$$\frac{dJ}{d{X}_{fl}}=\frac{\partial J}{\partial X_{fl}}-{\psi }^T\frac{\partial R}{\partial X_{fl}}.$$

(32)

Therefore, the first term at RHS linearized at the viscous walls with heat flux results in:

$$\frac{\partial J}{\partial X_{fl}}=\frac{\partial J}{\partial V_{bnd}}\frac{\partial V_{bnd}}{\partial X_{fl}}=\frac{\partial J}{\partial T_{FW}}\frac{\partial T_{FW}}{\partial V_{bnd}}\frac{\partial V_{bnd}}{\partial X_{fl}}={\overline{T}}_{FW}\frac{\partial T_{FW}}{\partial V_{bnd}}\frac{\partial V_{bnd}}{\partial X_{fl}}.$$

(33)

Equations (31) and (33) show the contributions of the adjoint fluid wall temperatures ${\overline{T}}_{FW}$ calculated in Section 3.5 to the fluid grid sensitivities at the wall boundaries exposed to the heat flux.

For the principle of the reverse accumulation, the wall heat flux term, entering in primal mode as the update of the boundary conditions in Equation (7), results in an adjoint fluid wall heat flux ${\overline{q}}_{FW}$ that is computed once the adjoint CFD solver approaches full convergence. This contribution is summed up with the ${\overline{q}}_{FW}$ returned in Section 3.5 by the reversed hFFB procedure.

$$T_{G1}=T_{D1}+q_{FW}\frac{\Delta n}{k_{fl}}\hspace{1em}\hspace{1em}\to \hspace{1em}\hspace{1em}{\overline{q}}_{FW}={\overline{q}}_{F{W}_{\overline{hFFB}}}+{\overline{T}}_{G1}\frac{\Delta n}{k_{fl}}.$$

(34)

A further contribution derives also from the second layer of ghost cells, not reported in Equation (34) for the sake of simplicity.

3.7. Adjoint hFFB Procedure (Fluid → Solid)

The partitioned coupling scheme in reverse mode completes by transferring the adjoint terms from the fluid to the solid domain. ${\overline{q}}_{FW}$ is processed by the differentiated DWI procedure, similarly to Equation (28), and returns the adjoint heat flux ${\overline{q}}_{SW}$ associated with the solid grid nodes. Consistently with the previous call, additional contributions are accumulated to the solid grid sensitivities ${\overline{X}}_{sl}$ and to the fluid grid sensitivities ${\overline{X}}_{fl}$ through the adjoint vector of mapped distances between each fluid cell center and the neighboring solid nodes.

Finally, the routine adopted to calculate the wall heat flux within the solid domain in Equation (3) is algorithmically differentiated, accounting only for the nodes laying on the interface. This process results in the new field of adjoint solid temperatures ${\overline{T}}_{SW}$, passed to the adjoint heat transfer solver described in Section 3.4.

The whole adjoint CHT process is iteratively repeated from Section 3.4, Section 3.5, Section 3.6 and Section 3.7 for as many loops as the ones performed during the primal computation till convergence.

Once walked through the entire workflow in reverse mode, ${\overline{X}}_{fl}$ and ${\overline{X}}_{sl}$, enclosing the contributions from the entire multidisciplinary chain, are introduced in Equation (12), and the constraint function sensitivities w.r.t. the design variables are finally computed. Hence, $dJ/d\alpha$ is evaluated by the SQP-based optimizer and the geometry is updated, opening the path to a new thermo-mechanical evaluation within the history of the shape optimization problem.

4. Validation

4.1. Flat Plate (Primal Mode)

The validation of the CHT process in primal mode is performed with the consideration of the analytic solution offered by [33] for the conjugate problem applied to a flat plate subjected to an incompressible flow. The mesh characteristics and boundary conditions are reported in

Table 1. Flat plate domain characteristics and boundary conditions.

Domain Settings	Value
Fluid domain length/height	0.25 m/0.1 m
Fluid mesh cells count	365,000
Plate thickness/length	0.01 m/0.2 m
Solid mesh nodes count/elements count	35,000
Fluid type	Air
Inlet flow total pressure	1.03 × 105 Pa
Inlet flow temperature	1000 K
Outlet flow static pressure	1.029 × 105 Pa
Plate temperature at bottom face	600 K
Plate thermal conductivity	0.29 W/m K
Virtual heat transfer coefficient $\tilde{h}$	100 W/m2 K

, while the fluid and solid meshes are qualitatively shown in Figure 8. The thermal conductivity adopted for the solid phase is aimed at obtaining an average Biot number around unity.

The convergence history of the maximum temperature difference between two successive fluid–solid iterations $L_{\infty }$ is reported in Figure 9, where a threshold of 1K is selected as the stabilization criterion for the model.

$$L_{\infty }={\left|T_{n+1}^s-T_n^s\right|}_{max},$$

(35)

The distribution of temperature at the interface between the two domains is presented in Figure 10a, comparing Luikov’s differential heat transfer (DHT) solution with the numerical one. Figure 10b shows the temperature variation in a vertical section located at x = 0.05m. In both cases, a sufficient agreement between the numerical results and the analytic solution is obtained within the limitations of the DHT approach, which does not capture the effect of the solid conductivity in the streamwise direction, as explained in [55].

4.2. Rotor Mesh Sensitivity Analysis

A mesh sensitivity study for the CHT analysis is performed on a three-dimensional test case presented in [19], considering a radial turbine impeller with 10 blades and an inlet diameter of 50 mm. A two-dimensional sketch of the domain on the meridional plane of the rotor is reported in Figure 11, along with the locations of the imposed boundary conditions. The method of the characteristics [56] is adopted to establish the number of physical conditions to be assigned at the inflow and outflow boundaries. The total upstream temperature and pressure are specified at the inlet of the domain, along with flow velocity components in the corresponding coordinate system. At the subsonic outflow boundary, a single flow variable is imposed, specifically the downstream static pressure. In conclusion, the boundary conditions applied to the problem are summarized in

Table 2. Turbine rotor boundary conditions.

Boundary Conditions	Value
Inlet total pressure p0	173 kPa
Inlet total temperature T0	1080 K
Inlet flow angle α from radial direction	62 deg
Outlet static pressure ps	101 kPa
Blade rotational speed ω	140,000 RPM

, in which the impeller rotational speed is also included.

The interactions of the rotor fluid and solid meshes are explored over three levels of refinement, as summarized in

Table 3. Fluid and solid mesh refinements, $\tilde{h}$ levels.

Factors Levels	Value
Fluid domain “coarse”—cells count	0.8 M
Fluid domain “mid”—cells count	1.3 M
Fluid domain “fine”—cells count	2.1 M
Solid domain “coarse”—nodes count	105 k
Solid domain “mid”—nodes count	295 k
Solid domain “fine”—nodes count	1.1 M
$\tilde{h}$ low–mid–high (W/m2 K)	800–1000–1300

. Additionally, three values of the virtual heat transfer coefficient $\tilde{h}$ are investigated in order to determine a suitable trade-off between stability and computational time. These factors are explored following the combinations expressed by the L9 Orthogonal Array method [57], returning nine test cases (

Table 4. L9 Orthogonal Array applied to the turbine rotor CHT analysis.

Test Case Number	CFD Mesh	FEM Mesh	$\tilde{\mathit{h}}\left[\mathbf{W}/{\mathbf{m}}^2\mathbf{K}\right]$
1	coarse	coarse	800
2	coarse	mid	1000
3	coarse	fine	1300
4	mid	coarse	1000
5	mid	mid	1300
6	mid	fine	800
7	fine	coarse	1300
8	fine	mid	800
9	fine	fine	1000

).

All the simulations are run till convergence at a relative residual drop of 10⁻⁸ for the computational fluid dynamics (CFD) solver and 10⁻¹⁰ for the FEM solver, while the CHT workflow is stopped for a maximum deviation in wall temperature between two successive fluid-solid iterations L_∞ below 1 K. y+ values below unity are obtained for the “mid” and “fine” levels of the fluid grid refinement, while the coarsest one exhibits a peak in y+ around 2.5 at the blade leading edge. Concerning the maximum number of nodes in the solid mesh, the code was tested up to the finest possible grid before incurring in memory leakage issues with the adopted Conjugate Gradient iterative solver.

The results in

Table 5. Fluid–solid grid sensitivity: summary of CHT computations.

Case#	#CHT Loops to Convergence [-]	Norm. Computational Time for CHT Iteration	Computational Time: CFD—FEM w.r.t. Total	Maximum Solid Temperature [K]	Delta Temperature Integral
1	11	1.0×	97.0–0.5%	1019	0.121
2	10	1.02×	95.2–1.9%	1025	0.053
3	9	1.36×	71.3–25.0%	1027	0.045
4	8	1.76×	98.3–0.3%	1026	0.092
5	10	1.78×	97.2–1.1%	1028	0.038
6	12	2.12×	81.6–16.1%	1028	0.023
7	11	2.4×	98.6–0.2%	1026	0.074
8	11	2.44×	97.2–0.8%	1028	0.019
9	8	2.76×	86.5–12.4%	1028	0.000

demonstrate a low sensitivity of the maximum temperature in the solid domain due to the settings applied to the coupling process, except for the first test case with coarse meshes at both the fluid and solid side. A critical aspect is represented by the quality of the interpolation of the information between the two domains realized by the fine-tuned distribution of virtual grid points on the interface at fluid side, as discussed in Section 2.1. Moreover, the finest solid mesh increases the resolution of the convective loading, resulting in more accurate temperature predictions in the material and, therefore, also more reliable heat fluxes returned to the fluid. Figure 14 shows two cases of solid mesh size, the coarse one (Figure 14a) and the intermediate one (Figure 14b), with the addition of the trailing edge refinement described hereafter. In the case of finer solid grids, the interpolation procedure returns a richer temperature pattern at the interface, closely mirroring the fluid conditions at the walls, whilst the coarsest one approximates the thermal loading distribution with the highest deviations localized in the blade tip region, where large secondary flows and leakages are present, and at the trailing edge. In this respect, a surface temperature analysis at a cross section located at about 25% of the chord was performed in order to evaluate the integral of the temperature deviations between each investigated test case w.r.t. test case 9. The axial position of the planar section was chosen in a region of high convective perturbations due to the presence of a significant flow detachment at the blade leading edge. The results, summarized in Table 5, are supported by the example in Figure 12, showing the temperature contours for case 1 and case 9. It is interesting to note that, in both set-ups, the blade tip returns significant thermal gradients because of the thinner geometry exposed to the flow vortices. However, the coarsest mesh overestimates the temperature drop, resulting in the integral reported in Table 5.

Finally, considerations about the computational time are synthesized in Table 5, with a dominant portion covered by the CFD solver and a significant rise in overhead emerging from the FEM heat transfer computations at the finest grid level. The normalization refers to the total duration of a single CHT iteration, including also the hFFB procedure. The comparison is performed w.r.t. the reference duration of test case 1, presenting the coarsest meshes at both fluid and solid side.

The interactions among the three factors in the Orthogonal Array from Table 4 are studied through the compounded signal S, defined as:

$$S=\omega A+\left(1-\omega \right)B,$$

(36)

with the term A enclosing the normalized difference between the maximum temperature predicted by the test case of interest and test case 9 (showing the finest grids), B representing the normalized computational time, and ω as the weighting coefficient. In the current study, ω = 0.7 in order to bias the objective function towards the accuracy of the coupling process. Indeed, the factors levels finally expected from this study are the ones minimizing the signal S.

The analysis of the Orthogonal Array [57], resulting in Figure 13, reveals the dependence of the signal S from the three variables and their correspondent levels. The chart illustrates the virtual heat transfer coefficient $\tilde{h}$ is selected at its highest value within the stability region of the Biot number discussed in [46], as it influences the B term in Equation (36), determining the rate of convergence of the CHT coupling process. Indeed, the current problem exhibits an optimal $\tilde{h}$ = 1000 W/m² K, since the highest value approaches the limit of the stability region, with incipient oscillations in the value of heat flux exchanged at the interface, before reaching convergence. The solid grid size, instead, characterizes the accuracy of the interpolation of the quantities passed between the two domains through the DWI process and affects the A term in Equation (36). In general, finer solid grids are favored by the current analysis. Remarkably, the fluid mesh size follows an opposite trend because its influence on the A term is less pronounced, as the coarsest level already provides good-quality solutions. On the other hand, the signal S is driven by the increase in computational overhead experienced with the fluid mesh refinement.

Since the weighting coefficient ω privileges the accuracy of the CHT coupling process, the final selection of the factors levels results in the intermediate values for the refinements of both the grids and for the virtual heat transfer coefficient. However, it is recognized that smaller solid element sizes improve the stability of the CHT coupling process because it avoids local poor quality in the discretization of the blade surface in correspondence with thin regions (Figure 14a), potentially inducing drops in the local Biot number and inconsistencies with the selected $\tilde{h}$ value from a stability standpoint. Therefore, the issue is addressed by the generation of a “hybrid” configuration (Figure 14b), envisaging a local solid mesh refinement in correspondence with the blade tip surface and trailing edge, whilst maintaining the intermediate element size in the rest of the domain. This last setup increases the mesh density to a total node count of about 420 k.

Finally, the Orthogonal Array analysis is repeated, assigning the integral of temperature deviations to the A term in Equation (36). The outcome, shown in Figure 15, confirms the previous considerations, except for the fluid mesh size, which promotes the trend towards the coarsest mesh because it still provides sufficiently accurate results. However, since the difference between the coarse and intermediate meshes is moderate, the previous selection of the medium size grid will be pursued for the sake of improved prediction accuracy.

Based on such considerations, the settings returned by the analysis of the Orthogonal Array are considered for the further assessments of the coupling problem.

4.3. Gradients Accuracy Evaluation

The present work accounts for the manual differentiation of the multidisciplinary process shown in Figure 6. The sensitivities are accurately computed by the reverse Algorithmic Differentiation technique [54] for what concerns the FEM solvers and the hFFB process, while the partitioned coupling approach opens the path to the direct integration of the in-house adjoint CFD solver.

The validation of the differentiated FEM solvers is presented herein with reference to the rotor test case from Section 4.2. The adjoint sensitivities $\partial J/\partial \alpha$ of the maximum temperature and of the maximum von Mises stress including the thermal strain w.r.t. selected design parameters (referenced in

Table 6. Design variables adopted in the dJ/dα gradient validation.

Design Variables α
1	Rotor hub meridional contour: Y-coordinate at 20% chord;
2	Rake angle;
3	Back-plate thickness;
4	Turbine shaft diameter;
5	Turbine shaft length;
6	Rotor maximum diameter;
7	Blade height at leading edge;
8	Back-plate/shaft connection axial position;
9	Blade hub thickness;
10	Blade hub fillet radius.

and Figure 16) are compared to the correspondent gradients computed by the Finite Differences (FD) technique. The latter accounts for evaluations by a central differencing scheme, whose optimal step size is searched with the aim of accurate gradients computations. An exemplary case is reported in Figure 17, showing the identification of a suitable step size at the value of 10⁻⁶ m for the seventh design variable in the framework of thermal evaluations by the heat transfer solver. The magnitude of the perturbation must avoid too-small values, resulting in round-off errors, and too large values as well, because they introduce significant truncation errors. This outcome results from the evaluation of the cost function at each step size by solving the linear system obtained after the perturbation of the mesh through the morphing technique described in [20].

The sensitivities of the maximum solid temperature to perturbations of the ten design variables are summarized in Figure 18a, which shows a close agreement between the two methods, both in sign and magnitude. To understand the sensitivity analysis in detail, let us consider the thermal paths in the turbine rotor experiencing the convective loading resulting from the operative condition reported in Table 2, along with a convective condition on the back-plate surface with a uniform fluid temperature of 950 K and a Dirichlet boundary condition of 500 K assigned to the extreme section of the shaft.

Figure 19 reveals that the maximum temperature is detected at the nodes in the region of the hub, in proximity with the leading edge. The computation of such constraint through a p-norm function returns a marked sensitivity from the rotor back-plate thickness, whose enlargement would reduce the influence of the heat sink located at the back-plate outer surface. Similarly, an increase in blade height at the leading edge would extend the thermal path to the colder areas highlighted at the tip, while an elongation of the shaft itself would reduce the gradient field. On the contrary, an increase in the turbine shaft diameter, an advancement of the axial position of the intersection between the rotor back-plate and the shaft, and an enlargement of the blade hub thickness would favor the cooling of the hub upper region. Such sensitivities provide the descent direction for a constrained optimization problem aimed at limiting the maximum rotor temperature for the sake of extended lifetime.

Similarly, the accuracy of the adjoint structural solver sensitivities is assessed by comparison with the same gradients computed by Finite Differences. The results reported in Figure 18b confirm the suitability of the manually differentiated framework.

Moreover, the time for the computation of the adjoint-based gradients for the heat transfer solver and the structural solver was measured. If X is the time required for the computation of the multidisciplinary workflow in primal mode and the problem accounts for n design variables, the calculation of the gradients by Finite Differences would approximately cost n ∙ X. On the contrary, the application of the adjoint method herein described costs about 8.6X for the heat transfer module and 2.3X for the structural analysis, in fair agreement with the original expectations. Such costs mostly depend on the assembly process of the differentiated system (Equation (25)) rather than the solver itself.

Finally, the loose-coupling approach adopted for the development of the CHT workflow allows the direct integration of the in-house adjoint CFD solver, and the accuracy of the computed gradients was already demonstrated by the same author.

5. Results and Discussion

The design work aims at reducing the thermo-mechanical stresses in the turbine impeller presented in Section 4.2. The operative condition in Table 2 is considered, along with a Dirichlet boundary condition with a temperature of 500 K imposed at the extreme cross section of the shaft, emulating the heat sink provided by the oil cooling circuit in the turbocharger bearing housing. The contribution of the thermal strains introduced by the present work improves the accuracy of the computation of the von Mises stresses, whose gradients provide information for the geometry modifications.

According to the grids selection anticipated in Section 4.2, the fluid domain is discretized with a multi-block structured mesh of about 1.3 M cells with boundary layer refinement. The solid domain accounts for an unstructured grid of approximately 420 k nodes and second-order tetrahedral elements.

Figure 20 presents the temperature distribution in the baseline geometry under the selected steady-state operative condition. A significant thermal gradient is identified at the connection between the shaft and the rotor back-plate, with a correspondent local variation of thermal strains according to Equation (9). This is reflected in increased localized stresses, as demonstrated in Figure 21. In fact, Figure 21a shows the mechanical pattern of the baseline geometry, whose structural computation is performed considering only the centrifugal loading. On the other hand, Figure 21b reveals the new stress distribution accounting also for the thermal field, resulting from the convective input and the heat sink at the shaft. Indeed, a local deviation in von Mises stresses of about 20% is identified in the region around the connection with the back-plate.

The suitability of the process described in Section 3 is tested with the computation of the sensitivities of the thermo-mechanical constraint w.r.t. the grid coordinates, as reported in Figure 22, providing the information for a grid perturbation. Therefore, Figure 23 compares the distribution of von Mises stresses in the baseline geometry (left) and in the updated layout (right). A local reduction in the maximum von Mises stress in excess of 35% is achieved with the new configuration. Most of the advancements can be attributed to the increased curvature of the back-plate contour, in particular in the region of the connection with the shaft, originally experiencing the highest concentrated stresses. The new shape reveals a more gradual mechanical pattern at such intersection, supported also by a smoother temperature evolution lowering the influence of the local thermal strains. The new layout returns lower sensitivity values, in line with the goal of satisfying the thermo-structural constraint within an optimization problem.

Other areas of the blade are highlighted in Figure 22, in particular at high blade span and in the upper portion of the back-plate. The activation of the sensitivities in such regions is favored by the contribution of the thermal input. However, these portions of the rotor present a pattern of smaller stresses, not violating the original constraint, and therefore their sensitivity fields are neglected.

The second positive effect of the implementation of the CHT analysis within the design loop is about the more accurate prediction of the actual turbine efficiency, as no more adiabatic walls assumptions are in place. The deviation in the total-to-static efficiency w.r.t. an adiabatic simulation is about 0.2% in the selected operative condition while considering the sole impeller. However, more marked differences are expected in case the CHT analysis is extended to the whole turbine stage, including the volute and the nozzle guide vanes.

6. Conclusions

This paper addresses the problem of the Fluid–Structure Interaction in the design optimization of turbomachinery components under thermo-mechanical constraints. The activity focuses on a radial turbine rotor for turbocharger applications, with the aim of reducing the maximum von Mises stress located in a critical area of the impeller experiencing the influence of the thermal loading.

The problem of the integration of high-fidelity thermal predictions within the framework of an optimization is addressed by the development of a novel multidisciplinary gradient-based approach, suitable for the efficient computation of the gradients of the cost function with respect to a large number of design variables. The introduction of a Conjugate Heat Transfer evaluation procedure, supported by the heat transfer forward flux back method, is undertaken in order to establish the energy coupling of the solid domain with its surrounding fluid, according to a partitioned coupling approach. The gradient calculation is addressed by a discrete adjoint technique, and the resulting sensitivities of the thermo-mechanical constraint w.r.t. the grid coordinates allow for the mechanical optimization of the component, including considerations related to the impact of its thermal operative environment.

The solid–fluid grid coupling process, the analysis procedure, and the multidisciplinary differentiation technique are described. Moreover, the validation of the conjugate problem is discussed with reference to the comparison with an analytic solution and a mesh sensitivity study. The accuracy of the computed gradients is also presented.

The application of the multidisciplinary workflow demonstrates the possibility of controlling the local von Mises stresses within admissible ranges, extending the outreach of the original mechanical evaluations with the inclusion of the thermal stresses. New regions of the impeller, before not affected by such considerations, are now involved, as their relative sensitivities are enhanced by the local thermal field. In particular, in the tested sample an improvement exceeding 35% in local von Mises stresses is achieved by the application of the presented method. Additionally, more accurate performance predictions result from the introduction of the CHT analysis within the design optimization process, no more relying on the adiabatic walls assumption.

This contribution is a new step towards the large-scale optimization of turbomachinery components by gradient-based methods, with the goal of the future introduction of fatigue life evaluations in transient conditions.