# Integration of second-order sensitivity method and cokriging surrogate model

Zebin Zhang, Martin Buisson, Pascal Ferrand, Manuel Henner

► To cite this version:

Zebin Zhang, Martin Buisson, Pascal Ferrand, Manuel Henner. Integration of second-order sensitivity method and cokriging surrogate model. Mathematics, 2021, 9 (4), pp.401. 10.3390/math9040401. hal-03573259

HAL Id: hal-03573259

https://hal.science/hal-03573259

Submitted on 4 Apr 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Integration of Second-Order Sensitivity Method and CoKriging Surrogate Model †

Zebin Zhang 1,*, Martin Buisson 2, Pascal Ferrand 2 and Manuel Henner 3

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Zhang, Z.; Buisson, M.; Ferrand, P.; Henner, F. Integration of Second-Order Sensitivity Method and CoKriging Surrogate Model.

*Mathematics* **2021**, *9*, 401. https://doi.org/10.3390/math9040401

Academic Editor: Giovanni Angiulli

Received: 26 January 2021

Accepted: 13 February 2021

Published: 18 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

**Abstract:** The global exploring feature of the surrogate model makes it a useful intermedia for design optimization. The accuracy of the surrogate model is closely related with the efficiency of optima-search. The cokriging approach described in present studies can significantly improve the surrogate model accuracy and cut down the turnaround time spent on the modeling process. Compared to the universal Kriging method, the cokriging method interpolates not only the sampling data, but also on their associated derivatives. However, the derivatives, especially high order ones, are too computationally costly to be easily affordable, forming a bottleneck for the application of derivative enhanced methods. Based on the sensitivity analysis of Navier-Stokes equations, current study introduces a low-cost method to compute the high-order derivatives, making high order derivatives enhanced cokriging modeling practically achievable. For a methodological illustration, second-order derivatives of regression model and correlation models are proposed. A second-order derivative enhanced cokriging model-based optimization tool was developed and tested on the optimal design of an automotive engine cooling fan. This approach improves the modern optimal design efficiency and proposes a novel direction for the large scale optimization problems.

**Keywords:** sensitivity method; surrogate model; cokriging; optimization; CFD

## 1. Introduction

With its inherent benefit of uncertainty prediction, Kriging and cokriging methods have drawn much attention in recent decades. An overwhelming number of approaches emerge posteriori to the design and analysis of computer experiments (DACE) approach [1,2], including surrogate-based optimization (SBO), surrogate-assisted evolutionary algorithms [3]. Some recent approaches show impressive performance for high dimensional problems with 30–200 design variables [4,5]. However, the improvement on model accuracy remains to be further studied. Studies on model reliability, either from local or global points of view [6,7], follows two main paths. First, by applying different training strategies [8], one tries to improve the local precision of the surrogate model, especially that of the optima neighborhoods, so-called “design space exploitation”. This is particularly useful, and at times indispensable for large-scale optimization. Recent studies also deal with different strategies to find the balance between exploration and exploitation by using multi-fidelity [9] and an ensemble of different surrogate methods [10]. Second, by introducing assistant information such as derivatives to the original dataset, the surrogate model can be improved. The most typical development concerns the integration of the adjoint method and the gradient-enhanced cokriging model in the aerodynamic shape

<footnote n="1">¹ School of Mechanical and Power Engineering, Zhengzhou University, Zhengzhou 450001, China</footnote>

<footnote n="2">² Laboratoire de Mécanique des Fluides et d'Acoustique (LMFA), Ecole Centrale de Lyon, 69134 Ecully, France; martin.buisson@ec-lyon.fr (M.B.); pascal.ferrand@ec-lyon.fr (P.F.)</footnote>

<footnote n="3">³ Simulation and Reliability Metier, Valeo Thermal Systems, 78320 La Verrière, France; manuel.henner@valeo.com</footnote>

<footnote>* Correspondence: zebin.zhang@zzu.edu.cn</footnote>

<footnote>† This paper is an extended version of our paper published in the World Congress on Global Optimization (WCGO 2019), the 6th international conference on Optimization of Complex Systems: Theory, Models, Algorithms and Applications, Metz, France, 8–10 July 2019.</footnote>
optimization [11]. What is worth noticing here, the co-Kriging method [9], which uses the correlations between several models of various fidelities, is not the same with the cokriging method in this study which interpolates the samples and the associated derivatives at the same time, generating one-single surrogate model. One remarkable advantage of adjoint method lies in the independence of the number of design variables, making this approach an appropriate choice for high-dimensional problems. However, the gradients obtained from this method depend on some given objective functions because the reversal mode differentiation is used, consequently, adding new objectives requires a new adjoint-state evaluation. Furthermore, for the modeling of multimodal functions, the higher is the order of derivatives used, the better is the model accuracy [12,13].

Yet the computational cost, being the main obstacle, limits the application of the high-order approach. With a classical differentiation method such as finite differences, depending on the discretization scheme, the number of evaluations required multiply by a factor of *n* (the number of parameters) for the differentiation of every higher degree of order. Although the higher accuracy brought by the derivatives, the total computational effort remains quasi-identical with a classical derivative-free method. The integration of derivatives in a surrogate-modeling process is not efficient unless the derivatives can be obtained in a cheap way.

This paper introduces the integration of the sensitivity analysis method and the derivative enhanced cokriging surrogate model. The former allows us to pursue the high-order derivatives in a cost effective way. The latter can make the best of the auxiliary information and generate a far-more accurate surrogate model.

The sensitivity analysis is to obtain the sensitivity of the output to the input by using their differential relations [14]. In the optimal design, the sensitivity analysis is often used for the quantitative study of the variation of the objective functions in response to some perturbations of the parameters. It can also be used to evaluate the robustness of the optimal design, performed at the post-processing phase. However, the implementation of sensitivity associated methods to the modeling phase is also feasible. The sensitivity information at the design site can be used to build the surrogate model. This is particularly useful for the problems featuring important nonlinear characteristics.

For the aerodynamic shape optimization, there are several ways to compute the sensitivities of design variables: finite difference method, complex step method [15,16], automatic differentiation and sensitivity equation method. E.Colin, Etienne, Pelletier and Borggaard brought forth the idea of continuous sensitivity equation method, illustrating how to obtain a fast solution of a flow field and its sensitivities with the finite element method [17]. The sensitivity equations are formed by implicitly differentiating the Navier-Stokes equations. Mahieu et al. proposed the generalized formulation of this approach which allows computing the first and second order derivatives of any parameters [18]. Please note that some key problems remain unsolved due to the high nonlinearity of Navier-Stokes equations and the derivability difficulty of turbulence model used in Reynolds-averaged Navier-Stokes equations.

A high-order extension of sensitivity equation method has been studied by Aubert, Ferrand et al. [19,20]. Instead of differentiating the governing equations directly, derivatives are computed *via* a direct mode automatic differentiation [21] applied to the discretized reference solutions [22,23]. As with the adjoint-state method, this approach applies its reverse mode to get high-order derivatives. Some non-commercialized tools have been developed, namely Turb'Flow, Turb'opty and Turb'post, which are flow solver, sensitivity solver and flow extrapolation tool respectively. Based on these tools, Soulat has accomplished an aerodynamic shape optimization for a fan system used in the automotive engine cooling module [24]. High order derivatives were used to build a polynomial surrogate model, which can be viewed as an aerodynamic database extrapolated from one single reference design. Taylor expansion was employed for the extrapolation. Optimization based on this model brought a significant reduction of aerodynamic tip loss due to secondary flow. However, only one single reference design has been considered in his study, so the truncation
tion error presented by nature, and the effective region covered by the model is limited to a certain radius. By introducing derivative-assisted cokriging method, it is possible to enlarge the covered region.

In summary, to make the best of their respective advantages, the high order sensitivity method and cokriging surrogate model are integrated in current studies. The high-order sensitivity approach is firstly introduced, then its integration with Hessian-enhanced cokriging model is presented. This methodology is employed for the aerodynamic shape optimization of an automotive engine cooling fan.

## 2. Sensitivity Approach in CFD Simulation Based Aerodynamic Shape Design

At stationary equilibrium state, the Navier–Stokes equations can be expressed in the following form:

$$F(q(P), P) = 0 \qquad (1)$$

With $F$ the flux vectors: the mass, the momentum and the energy in terms of con- servative variables $q(\rho, \rho U, \rho E)$, and the transport of turbulent variables such as $k$ and $\omega$ used in Reynolds-Averaged Navier–Stokes equations (RANS) [20]. $P$ stands for the design parameter, which can be either geometric or physic parameter. Differentiation of this equation with respect to $P$ yields:

$$\frac{\partial F}{\partial q}|_{P}(q, P) \cdot \frac{dq}{dP} + \frac{\partial F}{\partial P}|_{q}(q, P) = 0 \qquad (2)$$

It can also be written:

$$\frac{\partial F}{\partial q}|_{P}(q, P) \cdot q^{(1)} = -\frac{\partial F}{\partial P}|_{q}(q, P) \cdot dP \qquad (3)$$

in which $q^{(1)}$ is the first order variation of $q$ by varying $P$ by $dP$.

The Jacobian matrix $G(q, P) = \frac{\partial F}{\partial q}|_{P}$ allows writing the $n^{th}$ order differentiation as $\frac{d^n F}{dP^n} = \frac{\partial F}{\partial q} \cdot q^{(n)} = G \cdot q^{(n)}$. The right part $R = -\frac{\partial F}{\partial P} \cdot \Delta P$ is the flux vector variation in function of $\Delta P$. Hence the equation can be represented in this form:

$$G \cdot q^{(1)} = R \qquad (4)$$

The $n^{th}$ order differentiation $q^{(n)}$ can be constructed by recursively differentiating Equation (4).

$$\begin{aligned} G \cdot q^{(2)} &= R^{(1)} - G^{(1)} \cdot q^{(1)} \\ G \cdot q^{(3)} &= R^{(2)} - G^{(2)} \cdot q^{(1)} - 2G^{(1)} \cdot q^{(2)} \\ \dots \\ G \cdot q^{(n)} &= R^{(n-1)} - \sum_{i=1}^{n-1} C_{n-1}^{i} G^{(i)} \cdot q^{(n-i)} \end{aligned} \qquad (5)$$

Please note that the linear system is based on the Jacobian matrix $G$, which is usually available in a classical implicit Computational Fluid Dynamics (CFD) solver. The high-order derivatives $q^{(n)}$ can then be used to reconstruct the flow fields corresponding to parameter variations $\Delta P$, as shown in Figure 1.

The implementation of this methodology requires a coherence between the algorithm used for the discretization of the Navier–Stokes equations and that of the resolution of sensitivity equations. However, the derivability difficulty presented in the discrete RANS equations when using SST-K-$\omega$ turbulence model with Roe discretization scheme. Negative dissipation term has been found during the resolution of sensitivity equations, leading to a negative turbulent viscosity, which caused the crash of the solver. Positivity preservation study has been carried out but there is not yet any constructive conclusion.
<figure><img src="image_4.png" /><figcaption>Figure 1. Diagram of parametrization process (Ref: [19]).</figcaption></figure>

For a methodological illustration, this study adopts an alternative approach for the
derivative computation by using the finite difference method. This does not deprive the
originality of the methodology. To respect the original idea of the sensitivity method, a mesh
deformation method is applied [25]. The latter relies on a morphing technique which calcu-
lates the mesh nodes displacements with a Radial Basis Function (RBF) approach [26,27].
The mesh deformation is a RBF type propagation out of the displacements of some user-
defined control points which are usually positioned on the boundaries of computational
domain. A such type deformation is illustrated on a 2-D surface mesh in Figure 2, where
the dots refer to the control points. Displacements are attributed to the inner control points,
the outer ones are fixed deliberately. The displacements of mesh nodes are calculated
following a radial basis function defined *a priori*. This technique is employed to obtain
deformed mesh of an automotive engine cooling fan blade (Figure 3).

<figure><img src="image_5.png" /><figcaption>Figure 2. Mesh deformation driven by control points (Ref: [26]).</figcaption></figure>

One single reference mesh is used for one geometrical configuration. Derivatives can
be deduced from the results of some new CFD simulations based on deformed meshes.
Compared with the finite difference discretization scheme which parametrizes the geometri-
cal configuration directly, followed by a re-meshing process for each new configuration, the
mesh deformation approach deals with one single reference mesh, which conserves the total
mesh element number and assures the similar discretization of the computational domains.

Figure 3 compares the original surface mesh (on the left) with a deformed mesh (on the right) following a variation of sweep angle. As it is shown from the global view of the blade and local zoomed view to the trailing edge, the refinement levels are preserved, the total mesh element number remains unchanged after the deformation. The deformation of the blade surface mesh is propagated into the computational domain in a RBF pattern as presented before. The derivative calculation takes the same method presented in the reference [28].
<figure><img src="image_6.png" /><figcaption>Figure 3. Deformation example: sweep angle variation and the zoomed view to the trailing edge.</figcaption></figure>

### 3. Derivative-Enhanced Cokriging Method

Surrogate models of three objectives in function of three parameters (see for details in section of “Case study”), can be built out of the database constituted of the CFD and sensitivity analysis results. Different approaches are available, being either parametric or non-parametric [29]. A polynomial model, based on one single reference design, is taken in a previous study [28]. The cokriging method which interpolates several reference designs and their derivatives, can effectively reduce the truncation error and enlarge the covered region.

By assuming a stochastic process of all the interpolations upon the given sample points, the Kriging method models an unknown function through a mean function and a covariance function. It is often the best linear unbiased prediction from some given evaluations, i.e., the prediction error is minimized. The regression model can be considered to be the mean function of all the possible interpolation functions subject to the existing evaluations. Constant average functions are frequently used as the regression model [30,31]. However, once differentiated, high-order derivative-enhancing effect will be erased. On the contrary, a too high order regression model may provoke some unrealistic local oscillation similar to that of Runge phenomenon while equidistant samples are used for the interpolation [32,33]. Experience shows the second order polynomial function is good compromise and it is, therefore chosen for current study which uses databases up to second order. The regression matrix for a cokriging model consists of a regression part and a differentiated part which takes the following form:

$$ \mathcal{F}(x) = \begin{bmatrix} \mathcal{F}_0 \\ d(\mathcal{F}) \\ d^2(\mathcal{F}) \end{bmatrix} \qquad (6) $$

where $d(\mathcal{F})$ and $d^2(\mathcal{F})$ are respectively the first order and second order derivatives of the regression matrix $\mathcal{F}$.

#### Correlation Model

The correlation model reflects the spatial correlations between the points of the design site. The prediction error of a Kriging model is assumed to follow a random process, a careful choice of law of probability for this process is essential. Generally, the law of probability is unknown. Assumption has to be made in order to characterize the covariogram in the random field. The Gaussian model is the most-frequently used, which assumes a priori the second-order stationary conditions of the design site, making Kriging modeling a Gaussian process. The interpolation properties primarily depend on the local behavior of the random field. Near to the origin, the Gaussian model shows a quadratic
behavior. For the sake of derivability, it is chosen for current study, the correlation between two samples can be expressed as:

$$R(\theta, s_k^j, s_l^j) = \prod_{j=1}^{n} \exp(-\theta_j d_j^2) \qquad (7)$$

where $d_j = s_k^j - s_l^j$ is the Euclidean distance between two samples $s_k$ and $s_l$ at $j^{th}$ dimension. Where $j \in [1, n]$ denotes the dimension index, $k, l \in [1, m]$ are the index of the sample point.

The correlation model can be established by evaluating a hyper-parameter namely $\theta$, which can be obtained by the minimizing of the following function:

$$\varphi(\sigma, R) = \det(R)^{\frac{1}{m}} * \sigma(\theta)^2 \qquad (8)$$

where $m$ is the number of sample points, $\sigma$ is the standard deviation of the stochastic process of the samples, $R$ is the correlation matrix.

The correlation matrix can be written:

$$R = \begin{bmatrix} \operatorname{cov}(S, S) & \operatorname{cov}(S, dS) & \operatorname{cov}(S, d^2S) \\ \operatorname{cov}(dS, S) & \operatorname{cov}(dS, dS) & \operatorname{cov}(dS, d^2S) \\ \operatorname{cov}(d^2S, S) & \operatorname{cov}(d^2S, dS) & \operatorname{cov}(d^2S, d^2S) \end{bmatrix} \qquad (9)$$

where $cov$ stands for the covariance, $S$ represents the normalized design site, interpreted by the position of sampling points in a unit hyper-cube, $dS$ and $d^2S$ are the differentiated forms of $S$.

The covariance between any reference point $s_l$ and all the rest reference points is given by:

$$\operatorname{cov}(s_l^i, s_k^j) = \sigma^2 R(s_l^i, s_k^j) \qquad (10)$$

where $i, j \in [1, n]$ denote the dimension indices, $k \in [1, m]$, $k \neq l$ denote the sample indices, $\sigma^2$ is the variance of the sample.

An example of the correlation between the first order derivative of $i^{th}$ dimension for $k^{th}$ point and the first order derivative of $j^{th}$ dimension for $l^{th}$ point is given:

$$\operatorname{cov}(\frac{\partial S}{\partial s_k^i}, \frac{\partial S}{\partial s_l^j}) = \sigma^2 \frac{\partial^2 R(s_k^i, s_l^j)}{\partial s_k^i \partial s_l^j} \qquad (11)$$

Once the regression matrix $\mathcal{F}$ and the covariance matrix $R$ are calculated, the cokriging model $\hat{y}(x)$ can be deduced:

$$\hat{y}(x) = f(x)^T \beta + r(x)^T R^{-1} (Y - \mathcal{F} \beta), \qquad (12)$$

where $r(x)$ is the correlation vector between any design to be predicted and existing reference points, $Y$ is the response matrix, regression coefficient vector $\beta$ can be calculated by a least square estimation which gives:

$$\beta^* = (\mathcal{F}^T R^{-1} \mathcal{F})^{-1} \mathcal{F}^T R^{-1} Y \qquad (13)$$

Variance $\sigma^2$ is given by:

$$\sigma^2 = \frac{1}{n} (Y - \mathcal{F} \beta)^T R^{-1} (Y - \mathcal{F} \beta), \qquad (14)$$

The derivatives calculated by the sensitivity approach is integrated with the above described derivative-assisted cokriging method. Models can be built for all the objective
functions. Optimization can be performed based on these models. This methodology is applied to the optimal design of an automotive engine cooling fan.

**4. Case Study**

An automotive engine cooling fan, as shown in Figure 4, has been optimized based on the presented methodology. A reference fan design has been adapted for the particular surrounding of the vehicle underhood, where the downstream flow is radially deviated from its axis by the engine.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Geometrical parameter definition.</figcaption></figure>

The pressure rise $\Delta p$ from upstream to downstream of the fan wheel, the torque $T$ acting on the fan and the aerodynamic efficiency $\eta$ are considered to be objectives to be optimized. The first study has been performed with only one database, details can be found in the work of Zhang et al. [28] and references therein. The three parameters considered are: variations of stagger angles at mid-span and tip as shown in Figure 4, flow rate $Q$. For simplicity, the variations are illustrated with $\gamma_m$ and $\gamma_t$ for mid-span stagger angle variation and tip stagger angle variation respectively. Derivatives of the three objectives with respect to three parameters are computed based on CFD simulations performed on the deformed mesh.

Two geometrical configurations, functioning at two operating points, generates four reference designs. Based on the derivative-assisted reference data, surrogated models are built and are used for optimizations.

For one single operating point with first and second order derivatives, at most it describes an evolution of order 2. Yet, the characteristic curve in Figure 5, which describe two objectives: $\Delta p$ and efficiency at different flow rates, shows an inflection at the operating point (flow rate 2500 m³/h) and a linear evolution at the flow rates higher than 2500 m³/h, a behavior which can only be approximated by a curve of order 3 or higher. With one single flow rate, one cannot reproduce this characteristic for such an extended range. By associating two points, however, it is possible to depict a characteristic of order 5, hence enlarge the range of extrapolation.

Figure 6 illustrates the schema of coupling of different databases. The reference values are 2300 m³/h and 2600 m³/h for the flow rate, 0° and 3° for the stagger angle at mid-span, the stagger angle at the tip takes its reference relative value 0°. The circles represent the extrapolation from one single reference point, i.e., they are polynomial models. The rectangulars show the coupling of two or four databases. The characteristic curves of the two fans are shown in Figure 5.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Reference fan characteristic curves: Δp-Q and η-Q. T00M00 and T00Mp3 represent two different geometric configurations. The stagger angle at mid-span of T00Mp3 is 3° more than that of T00M00 (Reference configuration).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Metamodeling schema.</figcaption></figure>

With a processor of type Inter(R) Core(TM) i7-2640M, 8G physical RAM, the building of models C and G have consumed about 6 minutes of CPU time, the construction of model D takes about 23 minutes because of the inverse operation of a larger size of covariance matrix R. Time is considerably short comparing with a CFD run, which typically takes more than 1000 CPU hours.

Figure 7 (left) shows the prediction on Δp by the model C, where the stagger angle is fixed at its reference value. The two stars shows the Δp of two reference points. The error bar (red) is given by the surrogate model. CFD evaluations are illustrated with crosses. Notice that the blue line in the model does not pass by all the evaluated points because of all the CFD results, except for the reference ones, are used to calculate the derivatives but not to build the model directly. The error bar shows clearly the error is minimized as approaching the reference points.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Comparisons among models C, G and D: $\Delta p$-Q. $\gamma_m = 0$.</figcaption></figure>

Figure 7 (right) shows the comparison of three surrogate models: C, G and D. For the flow rate, the curve of model C and D passes by two reference points $Q = 2300 \text{ m}^3/\text{h}$ and $Q = 2600 \text{ m}^3/\text{h}$. The curve of model G does not pass by those points, which is originally from the error of second order cross derivative $\frac{\partial^2 \Delta p}{\partial(\gamma_m Q)}$. At high flow rate, model D is closer to model G, but not model C, this fact shows the influence from the dimension of $\gamma_m$. In fact, $\gamma_m$ has a similar effect with the flow rate. By increasing the $\gamma_m$, the air incidence is reduced, the blade is discharged, a phenomenon that can also be observed by increasing the flow rate. Model D, which couples four different databases, allows us to obtain more reliable results on a larger range of parameter variation. Hence this model is used to exploit the optima according to different criteria.

### 4.1. Mono-Objective Optimization Based on Model D

Several optimizations have been done, one of them is the maximization of efficiency by varying all the parameters at their full ranges. The goal of this optimization is firstly to have an idea about the maximum efficiency that can be found in the current design space, secondly to testify this approach of uniting the optimization algorithm and the surrogate model.

The genetic algorithm [34] is used with an initial population of 5000, evaluated on 100 generations. Thanks to the surrogate model, 500,000 evaluations have been done within a few hours, a task which is not achievable by CFD means. What is worth mentioning, the time-cost can be further reduced by using the long term memory assistance approach for evolutionary algorithms where the same solutions are not evaluated again [35]. This will be highly recommended when the number of parameters is higher, which should be considered in the future studies.

The result is shown in Table 1.

<table><caption>**Table 1.** Comparison between the reference design and the optimum on efficiency.</caption><thead><tr><th>Case</th><th>γ<sub>t</sub> (°)</th><th>γ<sub>m</sub> (°)</th><th>Q (m<sup>3</sup>/h)</th><th>Δp (Pa)</th><th>T (Nm)</th><th>η (%)</th></tr></thead><tbody><tr><td>Reference CFD</td><td>0</td><td>0</td><td>2300</td><td>226.4</td><td>0.9073</td><td>54.36</td></tr><tr><td>OptimA model</td><td>-1.30</td><td>2.25</td><td>2324.4</td><td>209.6</td><td>0.8265</td><td>55.90</td></tr><tr><td>OptimA CFD</td><td>-</td><td>-</td><td>-</td><td>209.4</td><td>0.8237</td><td>55.94</td></tr></tbody></table>

To assess the validity of the prediction, the results are compared with a CFD run. The differences are no more than 1% for all the three objectives, an error within the numerical uncertainty of the CFD simulation. Hence this approach is validated for the maximum efficiency design.

The optimization proposes a set of parameters which modifies the reference geometry and eventually leads to the required performance. The analysis of physics helps us to understand the local modification on the flow field by the optimizer.
The pressure distributions on six different radii have been drawn to have a global view of the fan blade charges (Figure 8).

<figure><img src="image_11.png" /><figcaption>Figure 8. The pressure distribution at six constant radii.</figcaption></figure>

This graph highlights the strong acceleration at the leading edge of the blade, which gradually amplifies at larger radius. The most important outcome of the optimizer is a discharged blade particularly at the region near the tip.

To better understand the flow structure modification, the radial profile of the axial velocity at about 4 mm downstream of the trailing edge, circumferentially averaged, is mapped out (Figure 9), where the reference fan “Reference” is compared with the efficiency-optimized fan “optimA”.

<figure><img src="image_12.png" /><figcaption>Figure 9. Axial velocity at the downstream of the fan wheel.</figcaption></figure>

Globally the radial evolution of the axial velocity of the optimum design is similar with that of the reference design. Both of them show a forced vortex characteristic, because only two stagger angles at mid-span and tip are modified. This also explains the similarity between two profiles near the hub. The axial velocity is small near the hub and linearly increases until a radius of 190 mm, beyond which the secondary flow and the wall effect of the shroud predominate. The optimum configuration maintains a null axial velocity at the tip, which allows the blade to work more efficiently at higher radius. This analysis shows the complexity of the flow near the shroud: tip vortex, secondary flow, brutal deviation of inlet velocity from axial to radial, high charges at higher radius. This study is carried out on only two geometrical parameters. According to this analysis, it will be useful to add more parameters such as the stagger angle at the hub section of the blade for the future study.

## 4.2. Multi-Objective Optimization Based on Model D

Based on model D, a multi-objective optimization has been performed to pursue an optimal design in term of performance at two different operating conditions: $Q_n = 2300 \text{ m}^3/\text{h}$
and $Q_i = 2800 \text{ m}^3\text{/h}$. At nominal operating point $2300 \text{ m}^3\text{/h}$, a higher efficiency is wanted, for $2800 \text{ m}^3\text{/h}$, a higher pressure rise $\Delta p$ is preferred in order to have an extended range of flow rate. The flow rate corresponding to $\Delta p = 0$ is called the "transparent" point, beyond which the fan module generates negative pressure rise, bringing extra drag for the vehicle. A transparent point with higher $Q$ value allows vehicle running at higher speed without generating extra drag. According to the fan characteristic curve shown in Figure 5, it is more likely to have a higher flow rate at transparent point if the $\Delta p$ at $2800 \text{ m}^3\text{/h}$ is higher.

Two objectives: efficiency at $Q_n = 2300 \text{ m}^3\text{/h}$, namely $\eta_n$ and pressure rise at $Q_i = 2800 \text{ m}^3\text{/h}$, namely $\Delta p_i$, have been considered. The algorithm NSGA-II [36] has been employed with 5000 individuals and 100 generations due to the inexpensive model-based evaluations. For each point, two performance evaluations have been done thanks to the surrogate model, one evaluation at $2300 \text{ m}^3\text{/h}$ and the other at $2800 \text{ m}^3\text{/h}$, hence 10,000 evaluations have been done for each generation. The bi-objective Pareto front is illustrated (Figure 10).

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Pareto front of optimization: maximization of $\eta_n$ (2300 m³/h) and $\Delta p_i$ (2800 m³/h).</figcaption></figure>

In Figure 10, the initial individuals, marked as red points, form a 2-dimensional projection on the objective surface $\eta_n - \Delta p_i$. The frontiers are clearly depicted, where the Pareto front can be seen on the top-right part, marked with black dots, formed of 502 survived individuals.

In order to illustrate the possible exploitation with the surrogate model, two optima have been adopted on the Pareto front, one favors the transparent point (optimB) and the other valorize the efficiency of nominal condition (optimC). Optimization results are compared with the reference configuration "Ref" in Table 2.
<table><caption>**Table 2.** Result of multi-objective optimization at two operating conditions.</caption><thead><tr><th>Geometry</th><th>γ<sub>t</sub> (°)</th><th>γ<sub>m</sub> (°)</th><th>Q (m<sup>3</sup>/h)</th><th>Δp (Pa)</th><th>T (Nm)</th><th>η (%)</th></tr></thead><tbody><tr><td>Ref</td><td>0</td><td>0</td><td>2300</td><td>226.4</td><td>0.9073</td><td><strong>54.36</strong></td></tr><tr><td>Ref</td><td>0</td><td>0</td><td>2800</td><td><strong>184.8</strong></td><td>0.9019</td><td>54.23</td></tr><tr><td>OptimB</td><td>-2.88</td><td>0.04</td><td>2300</td><td>242.0</td><td>0.9689</td><td><strong>54.36</strong></td></tr><tr><td>OptimB</td><td>-2.88</td><td>0.04</td><td>2800</td><td><strong>199.8</strong></td><td>0.9716</td><td>54.39</td></tr><tr><td>OptimC</td><td>-2.65</td><td>0.99</td><td>2300</td><td>232.1</td><td>0.9163</td><td><strong>55.20</strong></td></tr><tr><td>OptimC CFD</td><td>-</td><td>-</td><td>2300</td><td>235.1</td><td>0.9268</td><td><strong>55.28</strong></td></tr><tr><td>OptimC</td><td>-2.65</td><td>0.99</td><td>2800</td><td><strong>184.8</strong></td><td>0.9074</td><td>53.87</td></tr><tr><td>OptimC CFD</td><td>-</td><td>-</td><td>2800</td><td><strong>183.0</strong></td><td>0.9091</td><td>53.38</td></tr></tbody></table>

The numbers in italics are those values of objectives concerned in this optimization. By keeping the same efficiency with the reference geometry, OptimB manage to raise the pressure rise at 2800 m³/h by 15 Pa, or 8.1%. In addition, if the pressure rise at 2800 m³/h is kept unchanged, the efficiency at 2300 m³/h can be improved by 0.84, or 1.5% higher than the reference configuration. These improvements are analyzed according to the flow structure modifications.

The variation of stagger angles changes the radial distribution of fan load. Comparing to the reference geometry, the optimB shows an important variation of tip stagger angle, which charges the tip with a high incidence, so that a higher flow rate can be accepted. By keeping almost the same stagger angle at the mid-span, the efficiency stays unchanged comparing with the reference geometry. For the optimC, the tip stagger angle has little variation comparing with optimB, but a more important variation at mid-span. This reduces the incidence and eventually augments the ηn (2300 m³/h) by keeping the same level of Δpᵢ (2800 m³/h).

To better understand the optimization result shown on the Pareto front, a Self-Organizing Map (SOM) [37] visualization was created with all the individuals on the Pareto front. SOM can be used to make the connection between design-space and objective-space effectively. It is a type of neural network designed to understand high-dimensional data with help of its low-dimensional representations. Each variable (parameter or objective) is depicted as a two-dimensional map which preserves topological properties of the Pareto optimal solutions. One individual is always found at the same 2-D position from one map to another, distinguished with the other individuals according to their colors, which represent the magnitude of the variable.

In Figure 11, there are four maps which represent the values of two parameters and two objectives on the Pareto fronts. Only negative values are observed on the tip stagger angle plot. All the solutions are kept away from their reference value. By contrast, large variation (5.36°) is shown for the mid-span stagger angle, the positive values correspond to the individuals which give better efficiency but less pressure rise Δpᵢ, the negative values correspond to the individuals which give a less efficiency but a better Δpᵢ. A clear correlation can be drawn among three variables: mid-span stagger angle, efficiency ηn and pressure rise Δpᵢ. CFD simulations have been performed to validate the optimA and optimC, the maximum relative error is 1.3%, found for the Δp of optimC at 2300 m³/h, an error within the numerical uncertainty. The approach is again validated for this multi-objective optimization.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** SOM visualization of Pareto front, optimA (green diamond point) and optimC (blue square point).</figcaption></figure>

**5. Conclusions**

The integration of the second-order sensitivity method and the cokriging surrogate model was studied in this paper. The derivative-assisted information significantly improve the model accuracy. Optimizations are performed based on established surrogate model. One of the optimization results shows obvious improvements on both the aerodynamic efficiency and the torque for an engine cooling fan. One multi-objective optimum succeeds in enlarging the operating range of the fan, the other manages to keep the same range and improve the efficiency at nominal condition.

The contributions of this paper are twofold:

1. The cokriging method can make the best use of the derivatives and significantly improves the surrogate model accuracy;

2. The implementation of parameterization tool (Turb'Opty) has been pre-illustrated. The cycle of optimal design can be effectively shortened through low-cost derivatives. This study proves the validity and efficiency of this methodology.

**Author Contributions:** Conceptualization, Z.Z.; methodology, Z.Z.; software, Z.Z. and M.B.; validation, P.F. and M.H.; writing—original draft preparation, Z.Z.; writing—review and editing, M.B. and P.F.; visualization, Z.Z.; supervision, P.F. and M.H.; project administration, P.F.; funding acquisition, M.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by "Association Nationale de la Recherche et de la Technologie" (ANRT) in France.

**Acknowledgments:** The authors acknowledge the supports from CIRT and Fluorem S.A.S.

**Conflicts of Interest:** The authors declare no conflict of interest.
## Abbreviations

The following abbreviations and symbols are used in this manuscript:

<table><tr><td>CFD</td><td>= Computational Fluid Dynamics</td></tr><tr><td>RANS</td><td>= Reynolds-Averaged Navier-Stokes equations</td></tr><tr><td>RBF</td><td>= Radial Basis Function</td></tr><tr><td>F</td><td>= Flux over any arbitrary grid cell</td></tr><tr><td>q</td><td>= Flow variables</td></tr><tr><td>P</td><td>= Design parameter</td></tr><tr><td>G</td><td>= Jacobian matrix</td></tr><tr><td>m</td><td>= Number of sample points</td></tr><tr><td>n</td><td>= Number of dimensions</td></tr><tr><td>R</td><td>= Kriging correlation matrix</td></tr><tr><td>F</td><td>= Kriging regression matrix</td></tr><tr><td>θ</td><td>= Kriging correlation parameter</td></tr><tr><td>d<sub>j</sub></td><td>= Distance between two samples</td></tr><tr><td>S</td><td>= Set of a sampling</td></tr><tr><td>s<sup>i</sup><sub>k</sub></td><td>= i<sup>th</sup> dimension of k<sup>th</sup> sample</td></tr><tr><td>ŷ</td><td>= cokriging predictor</td></tr><tr><td>f x</td><td>= Regression function for prediction</td></tr><tr><td>β</td><td>= Coefficient vector of regression function</td></tr><tr><td>β*</td><td>= Generalized least squares solution of regression problem</td></tr><tr><td>r</td><td>= Correlation vector of a design to be predicted</td></tr><tr><td>Y</td><td>= Kriging objective matrix</td></tr><tr><td>σ</td><td>= Standard deviation</td></tr><tr><td>Q</td><td>= Volumic flow rate</td></tr><tr><td>Δp</td><td>= Pressure variation through fan system</td></tr><tr><td>η</td><td>= Static efficiency</td></tr><tr><td>T</td><td>= Torque of electric motor</td></tr><tr><td>γ<sub>m</sub></td><td>= Stagger angle at mid-span</td></tr><tr><td>γ<sub>t</sub></td><td>= Stagger angle at tip</td></tr><tr><td>C<sub>p</sub></td><td>= Pressure coefficient</td></tr></table>

## References

1. Sondergaard, L.; Nielsen, H. *Aspects of the Matlab Toolbox DACE*; Technical Report; Technical University of Denmark: Kongens Lyngby, Denmark, 2002.

2. Bartz-Beielstein, T.; Zaefferer, M. Model-based methods for continuous and discrete global optimization. *Appl. Soft Comput.* **2017**, *55*, 154–167. [CrossRef]

3. Jin, Y. Surrogate-assisted evolutionary computation: Recent advances and future challenges. *Swarm Evol. Comput.* **2011**, *1*, 61–70. [CrossRef]

4. Dong, H.; Wang, P.; Yu, X.; Song, B. Surrogate-assisted teaching-learning-based optimization for high-dimensional and computationally expensive problems. *Appl. Soft Comput.* **2021**, *99*, 106934. [CrossRef]

5. Cai, X.; Gao, L.; Li, X.; Qiu, H. Surrogate-guided differential evolution algorithm for high dimensional expensive problems. *Swarm Evol. Comput.* **2019**, *48*, 288–311. [CrossRef]

6. Forrester, A.I.J.; Keane, A.J.; Bressloff, N.W. Design and Analysis of “Noisy” Computer Experiments. *AIAA J.* **2012**, *44*, 2331–2339. [CrossRef]

7. Leifsson, L.; Koziel, S.; Tesfahunegn, Y.A. Multi-objective Aerodynamic Optimization by Variable-Fidelity Models and Response Surface Surrogates. *AIAA J.* **2016**, *54*, 531–541. [CrossRef]

8. Li, Y.; Shen, J.; Cai, Z.; Wu, Y.; Wang, S. A Kriging-Assisted Multi-Objective Constrained Global Optimization Method for Expensive Black-Box Functions. *Mathematics* **2021**, *9*, 149. [CrossRef]

9. Koziel, S.; Bekasiewicz, A.; Couckuyt, I.; Dhaene, T. Efficient multi-objective simulation-driven antenna design using co-Kriging. *IEEE Trans. Antennas Propag.* **2014**, *62*, 5900–5905. [CrossRef]

10. Ye, P.; Pan, G. Selecting the Best Quantity and Variety of Surrogates for an Ensemble Model. *Mathematics* **2020**, *8*, 1721. [CrossRef]

11. Han, Z.-H.; Gortz, S. Alternative Cokriging Model for Variable-Fidelity Surrogate Modeling. *AIAA J.* **2012**, *50*, 1205–1210. [CrossRef]

12. Rumpfkeil, M.P. Optimizations Under Uncertainty Using Gradients, Hessians, and Surrogate Models. *AIAA J.* **2013**, *51*, 444–451. [CrossRef]
13. Laurenceau, J.; Meaux, M.; Sagaut, P. Comparison of Gradient-Based and Gradient-Enhanced Response-Surface-Based Optimizers. *AIAA J.* **2010**, *48*, 981–994. [CrossRef]

14. Saltelli, A. *Global Sensitivity Analysis: The Primer*; John Wiley: Hoboken, NJ, USA, 2008.

15. Lyness, J.N.; Moler, C.B. Numerical Differentiation of Analytic Functions. *Siam J. Numer. Anal.* **1967**, *4*, 202–210. [CrossRef]

16. Lyness, J.N. Numerical Algorithms Based on the Theory of Complex Variable. 1967. Available online: https://dl.acm.org/doi/abs/10.1145/800196.805983 (accessed on 26 January 2021).

17. Colin, E.; Etienne, S.; Pelletier, D.; Borggaard, J. Application of a sensitivity equation method to turbulent flows with heat transfer. *Int. J. Therm. Sci.* **2005**, *44*, 1024–1038. [CrossRef]

18. Mahieu, J.-N.; Etienne, S.; Pelletier, D.; Borggaard, J. A Second-order sensitivity equation method for laminar flow. *Int. J. Comput. Fluid Dyn.* **2005**, *19*, 143–157. [CrossRef]

19. Aubert, S.; Ferrand, P.; Pacull, F.; Buisson, M. Fast CFD For Shape and Flow Parameterization With Meta-models built on High-order Derivatives. Application to Fast Design. In Proceedings of the ICAS 2010, 27th Congress of the International Council of the Aeronautical Sciences, Nice, France, 19–24 September 2010.

20. Soulat, L.; Ferrand, P.; Moreau, S.; Aubert, S.; Buisson, M. Efficient Optimisation for Design Problems in Fluid Mechanics. *J. Comput. Fluids* **2013**, *82*, 73–86 [CrossRef]

21. Dwight, R.; Brezillon, J.; Vollmer, D. Efficient Algorithms for Solution of the Adjoint Compressible Navier-Stokes Equations with Applications Odas. *DLR* **2006**, *60*, 365–389.

22. Bischof, C.H.; Bucker, H.M.; Rasch, A.; Slusanschi, E.; Lang, B. Automatic Differentiation of the General-Purpose Computational Fluid Dynamics Package FLUENT. *J. Fluids Eng.* **2007**, *129*, 652–658. [CrossRef]

23. Bischof, C.; Corliss, G.; Green, L.; Griewank, A.; Haigler, K.; Newman, P. Automatic Differentiation of Advanced CFD Codes for Multidisciplinary Design. *Comput. Syst. Eng.* **1992**, *3*, 625–637. [CrossRef]

24. Soulat, L. Definition, Analyse et Optimisation Aerodynamique d'un Nouveau Concept de Traitement de Carter au Moyen D'outils Numeriques. Application aux Turbomachines a Basse Vitesse. Ph.D. Thesis, Ecole Centrale de Lyon, Écully, France, 2010.

25. Zhang, Z.; Buisson, M.; Ferrand, P.; Henner, M. Databases Coupling for Morphed-Mesh Simulations and Application on Fan Optimal Design. Optimization of Complex Systems: Theory, Models, Algorithms and Applications. WCGO 2019. In *Advances in Intelligent Systems and Computing*; Springer: Cham, Switzerland, 2019; Volume 991, pp. 981–990.

26. Rozenberg, Y.; Benefice, G.; Aubert, S. Fluid Structure Interaction Problems In Turbomachinery Using RBF Interpolation And Greedy Algorithm. In Proceedings of the ASME Turbo Expo 2014: Turbine Technical Conference and Exposition, Düsseldorf, Germany, 16–20 June 2014; Volume 2B. [CrossRef]

27. Rendall, T.C.S.; Allen, C.B. Unified fluid- structure interpolation and mesh motion using radial basis functions. *Int. J. Numer. Methods Eng.* **2008**, *74*, 1519–1559. [CrossRef]

28. Zhang, Z.; Buisson, M.; Ferrand, P.; Henner, M.; Gillot, F. Meta-model Based Optimization of a Large Diameter Semi-radial Conical Hub Engine Cooling Fan. *Mech. Lndustry* **2015**, *15*, 102. [CrossRef]

29. Wang, G.; Shan, S. Review of Metamodeling Techniques in Support of Engineering Design Optimization. *J. Mech. Des.* **2007**, *129*, 370–380. [CrossRef]

30. March, A.; Willcox, K. Provably Convergent Multifidelity Optimization Algorithm Not Requiring High-Fidelity Derivatives. *AIAA J.* **2012**, *50*, 1079–1089. [CrossRef]

31. Yamazaki, W.; Mavriplis, D.J. Derivative-Enhanced Variable Fidelity Surrogate Modeling for Aerodynamic Functions. *AIAA J.* **2013**, *51*, 126–137. [CrossRef]

32. James, E. On the Runge example. *Amer. Math. Mon.* **2018**, *94*, 329–341. [CrossRef]

33. De Marchi, S.; Marchetti, F.; Perracchione, E.; Poggiali, D. Polynomial interpolation via mapped bases without resampling. *J. Comput. Appl. Math.* **2019**, *364*. [CrossRef]

34. David E. Goldberg. In *Genetic Algorithms in Search, Optimization and Machine Learning*, 1st ed.; Addison-Wesley Longman Publishing Co., Inc.: Boston, MA, USA, 1989.

35. Črepinšek, M.; Liu, S.-H.; Mernik, M.; Ravber, M. Long Term Memory Assistance for Evolutionary Algorithms. *Mathematics* **2019**, *7*, 1129. [CrossRef]

36. Deb, K.; Agrawal, S.; Pratap, A.; Meyarivan, T. A Fast Elitist Non-dominated Sorting Genetic Algorithm for Multi-objective Optimization: NSGA-II. In *Parallel Problem Solving from Nature PPSN VI*; Schoenauer, M., Ed.; PPSN 2000; Lecture Notes in Computer Science; Springer: Berlin/Heidelberg, Germany, 2000; Volume 1917.

37. Vesanto, J.; Himberg, J.; Alhomieni, E.; Parhankangas, J. SOM Toolbox for Matlab5 Report A57; Helsinki University of Technology: Espoo, Finland, 2000.