Multi-objective aerodynamic optimization of high-speed train heads based on the PDE parametric modeling

With the increasing running speed, the aerodynamic issues of high-speed trains are being raised and impact the running stability and energy efficiency. The optimization design of the head shape is significantly important in improving the aerodynamic performance of high-speed trains. Existing aerodynamic optimization methods are limited by the parametric modeling methods of train heads which are unable to accurately and completely parameterize both global shapes and local details. Due to this reason, they cannot optimize both global and local shapes of train heads. In order to tackle this problem, we propose a novel multi-objective aerodynamic optimization method of high-speed train heads based on the partial differential equation (PDE) parametric modeling. With this method, the half of a train head is parameterized with 17 PDE surface patches which describe global shapes and local details and keep the surface smooth. We take the aerodynamic drag and lift as optimization objectives; divide the optimization design process into two stages: global optimization and local optimization; and develop global and local optimization methods, respectively. In the first stage, the non-dominated sorting genetic algorithm (NSGA-II) is adopted to obtain the framework of the train head with an optimized global shape. In the second stage, Latin hypercube sampling (LHS) is applied in the local shape optimization of the PDE surface patches determined by the optimized framework to improve local details. The effectiveness of our proposed method is demonstrated by better aerodynamic performance achieved from the optimization solutions in global and local optimization stages in comparison with the original high-speed train head.


Introduction
High-speed trains play an irreplaceable role in modern means of transportation due to their advantages such as high running speed, ride comfort, large transport capacity and low energy consumption. With the increasing running Responsible Editor: Emilio Carlos Nelli Silva Yu Xia yxia@bournemouth.ac.uk 1 et al. 2016). However, this method is not ideal in accurately describing train heads because the local shape in surface patches is uncontrollable and has no ability to deform if the framework remains unchanged. Moreover, the tangential continuity at the joint of two adjacent surface patches cannot be guaranteed because the two adjacent patches are only constrained by the position coordinates of the framework. In order to optimize local shapes instead of the whole shape of a train head and deform the surface smoothly, various shape deformation methods have attracted considerable attention such as the free-form deformation (Li et al. 2016;Zhang et al. 2018) and the arbitrary shape deformation (Sun et al. 2010;Yao et al. 2014). The shape deformation method can focus on improving specified local shapes, but its global shape is non-optimized. Besides, the shape deformation methods usually use one parameter to control the shape change of a local region which limits the degree of freedom of the deformation. For example, the freeform deformation with one design parameter moves control points simultaneously in the same direction and makes the shape change have only one degree of freedom, i.e., sticking out or caving in, as shown in Fig. 1a.
The limitations of above modeling methods are that they cannot create complicated train head shapes with few design variables, strictly keep the surface smooth when the design variables change, and achieve both global and local shape changes. Normally, more design variables can describe more complicated train head shapes and provide more degrees of freedom to control the shapes in the modeling process, but it will be more difficult to maintain the surface smoothness and find an optimal solution in the optimization process. Therefore, a good parametric modeling method able to define detailed train head models is the key to success (Wang et al. 2018). The partial differential equation (PDE)-based modeling method is fully capable in creating a complex shape with three design variables which provides an effective solution to these problems.
PDE-based modeling is to create surfaces by the solution to a PDE with shape control parameters subjected to exact satisfaction of boundary conditions. It has four advantages (Castro et al. 2008). (1) Small data: complicated and detailed shapes can be described with few PDE surface patches. For example, the train head model is constructed by only 17 PDE surface patches in this paper. (2) Efficient shape deformation: applying shape control parameters (3 design variables) can accurately control surface shapes with more degrees of freedom compared with the free-form deformation method, as shown in Fig. 1. (3) Good continuities: any high-order continuity between two adjacent PDE surface patches is readily achieved and naturally maintained.
(4) Physics nature: the PDE method is physics-based which has a potential to create a more realistic appearance. Due to these virtues, PDE-based modeling has the strength in shape designs of high-speed train heads.

Aerodynamic analysis
In the aerodynamic analysis process, the numerical simulation, i.e., the computational fluid dynamics (CFD), is commonly used to investigate the aerodynamic performance of high-speed train heads. By comparing the simulation results with the wind tunnel test of high-speed trains, the CFD is proved to be accurate and acceptable for aerodynamic analysis of high-speed train heads (Cheli et al. 2010;Morden et al. 2015). The approximation method, i.e., the surrogate model, has been applied simultaneously with the CFD simulation to reduce the time of the optimization process (Li et al. 2016). The frequently used surrogate models include the response surface (Ku et al. 2010;Kwon et al. 2001), Kriging model (Li et al. 2016;Zhang et al. 2018;Sun Fig. 1 A comparison between the free-form deformation and PDE-based modeling methods. a The free-form deformation with 9 red control points which are governed by one design variable. b PDE-based deformation by changing the value of 3 control parameters (3 design variables) Yao et al. 2014), support vector machine (SVM) model (Yao et al. 2016;Lee and Kim 2008) and radial basis function network (Muñoz-Paniagua et al. 2014). Although these surrogate models can quickly obtain the optimization results, they all need to be trained by inputting a lot of training sets such as 168 sets in Lee and Kim (2008) before applying it into the optimization process. Especially, when the training sets come from numerical simulations, the computational accuracy will be decreased further since the surrogate model and numerical simulations are both approximation solutions. In addition, there are no reliable rules to estimate the accuracy of a given surrogate model. Depending on the complexity of a nonlinear problem and the quality of training sets, the training process of the surrogate model could easily be non-convergent or locally convergent.
The running environment of a high-speed train is very complex and affected by many aerodynamic characteristics such as the aerodynamic drag force, lift force, side force and micro-pressure wave (Yao et al. 2014). In multi-objective optimization, researchers usually choose some combinations of these aerodynamic characteristics as optimization objectives according to different running environments, and the aerodynamic drag and lift forces are frequently taken into account. Since the aerodynamic drag force of a highspeed train accounts for about 75-80% of the total resistance when the train speed reaches 250-300 km/h (Brockie and Baker 1990;Schetz 2001), a large aerodynamic drag force will significantly increase the energy consumption of the train and seriously affect its running speed. For the aerodynamic lift force, a large positive lift force will reduce the adhesion of the wheel-rail system and a negative lift force will increase the abrasion of the wheel-rail system (Li et al. 2016). Particularly, the unsteady characteristics of the wake flow have a predominant effect on the train which may cause a large fluctuation to the aerodynamic lift of the train tail (Yao et al. 2014). Therefore, the aerodynamic drag of the whole train and aerodynamic lift of the train tail are closely related to the safety and reliability of a running high-speed train.
As discussed above, the aerodynamic drag and lift forces are very important in influencing the aerodynamic performance of high-speed trains. We take them as two optimization objectives in the following aerodynamic analysis process.

Our work and contributions
In this paper, we proposed a novel multi-objective optimization method of high-speed train heads based on the PDE parametric modeling. We parameterize the half of a high-speed train head with 17 PDE surface patches by solving a vector-valued fourth-order PDE and apply the CFD simulation to analyze the aerodynamic performance, i.e., the aerodynamic drag and lift forces of the train. The optimization process of our method consists of two stages, i.e., the global and local optimizations. In the global optimization stage, we use ten design parameters as design variables to control the framework of the train head and the non-dominated sorting genetic algorithm (NSGA-II) to obtain the Pareto front of the globally optimized solutions. In the local optimization stage, we apply the three shape control parameters in PDE (1) below as design variables to change the local shapes of all PDE surface patches, and the Latin hypercube sampling (LHS) method is employed to generate the sample sets of three shape control parameters for finding the locally optimized solution. The detailed optimization process is explained in Section 3. Our main contributions can be summarized as follows: -A new multi-objective aerodynamic optimization method of high-speed train heads based on the PDE parametric modeling is proposed to optimize the head shape. -PDE surface patches are introduced to describe more detailed surface shapes with fewer design variables and guarantee the smoothness between two adjacent surface patches. -The global and local optimization methods in the multi-objective aerodynamic optimization process are developed to optimize both global and local shapes.
The remaining parts of this paper are organized as follows.
The related work on the optimization of high-speed train heads and the PDE-based modeling method are briefly reviewed in Section 2. An overview of the proposed multiobjective aerodynamic optimization process is given in Section 3. Our method of PDE-based modeling is described in Section 4. The optimization algorithm and the CFD simulation method are given in Section 5 and Section 6, respectively. The optimization results and discussions are presented in Section 7, and finally, the conclusion is drawn in Section 8.

Related work
Our proposed method is related to the optimization of highspeed train heads and the PDE-based modeling method. In this section, we briefly review the most related work in the two fields.

Optimization of high-speed train heads
In early stages, researchers usually focused on studying the aerodynamic performance of a given high-speed train using experimental methods such as the wind tunnel test and the moving model test (Schetz 2001). Maeda et al. (1989) presented a method to estimate the aerodynamic drag of trains and evaluated the accuracy of the results by measuring the total resistances in the open air and in a tunnel. De Wolf and Demmenie (1997) developed a train tunnel test facility which can launch models up to 500 km/h to measure the interacting pressure waves and their reduction in tunnels for high-speed trains. Cheli et al. (2010) presented an aerodynamic analysis on two different versions of the high-speed train EMUV250 to study the cross wind behavior by a combined use of wind tunnel investigations and numerical CFD simulations. However, the experiment method, greatly depending on engineering experience of researchers to seek the optimal solution, is limited by some disadvantages such as expensive testbed and high cost of time.
With the development of computer technology and increase of the computational power, numerical simulations have been widely applied in the optimization design of high-speed train heads and usually combined with surrogate models to accelerate the optimization process. Kwon et al. (2001) used the response surface methodology and the axisymmetric compressible Euler equations to optimize the nose shape and introduced the Hicks-Henne shape functions to define the design space. Ku et al. (2010) employed the vehicle modeling function for multi-objective optimization of the high-speed train nose and performed a multi-step design optimization using the BFGS algorithm with a response surface model. Lee and Kim (2008) adopted the Hicks-Henne shape functions to parameterize the high-speed train nose and presented an approximate optimization method to reduce the micro-pressure wave by using an SVM-based metamodel and sequential quadratic programming. Sun et al. (2010) proposed an optimization approach to improve the aerodynamic drag of a CRH3 high-speed train head by combining CFD analysis with the genetic algorithm and introduced an arbitrary shape deformation technology. These methods, however, focus on studying single-objective optimization problems which are inadequate to find optimal solutions.
The multi-objective optimization of high-speed train heads has been much investigated in the recent decade. The task of multi-objective optimization is to find a set of solutions, i.e., Pareto-optimal solutions, which represent a trade-off among objective functions (Deb 1999). Suzuki and Nakade (2013) developed a multi-objective optimization method of high-speed train heads using an evolutionary algorithm to estimate the aerodynamic drag and pressure variation, and parameterized the train head shape by Bspline curves and Coons patches. Li et al. (2016) optimized the aerodynamic drag and lift forces of a CRH2 high-speed train head using NSGA-II based on a Kriging model and applied five design variables to control the local shapes of the train head with a free-form deformation method. By using the similar modeling and analysis methods in Li et al. (2016), Zhang et al. (2018) studied the aerodynamic drag, lift and side forces of a high-speed train running on an embankment under crosswinds. Muñoz-Paniagua et al. (2014) defined the geometrical parameterization of the nose shape of a high-speed train by three design variables, and adopted the genetic algorithm with a radial basis function network to minimize the maximum pressure gradient and aerodynamic drag of the high-speed train entering a tunnel. Yao et al. (2014) constructed a multi-objective optimization process of a high-speed train head using the modeling method of arbitrary shape deformation and the optimization method of NSGA-II based on a Kriging model. They also studied a multi-objective particle swarm algorithm with a SVM regression model for the aerodynamic optimization of high-speed train heads (Yao et al. 2016).

PDE-based modeling
PDE-based geometric modeling was pioneered by Bloor and Wilson (1990) where they used a vector-valued fourthorder PDE with one shape control parameter to generate free-form surfaces and create different shapes like propeller blades and phone handsets. Subsequently, the PDE-based modeling technology has been widely used in industrial geometric designs such as the submarine (Bloor and Wilson 1994) and aircraft design (Athanasopoulos et al. 2009). In order to interactively create surfaces in real time, Ugail et al. (1999) developed an interactive design tool based on PDE-based modeling. Since one shape control parameter in a PDE was not enough to diversify the PDE surface shapes, a fourth-order PDE with three shape control parameters was proposed by Zhang and You (2002) to generate various shapes of 3D vase models. They also developed a sixth-order PDE with four shape control parameters which provides enough degrees of freedom to satisfy C 2 continuous boundary conditions and offers more shape control parameters to control surface shapes . Since PDE surfaces have not become an industrial standard in computer-aided design and manufacturing systems, Wang et al. (2019) developed a method of achieving optimal conversion of PDE surfaces into NURBS surfaces for the shape design of high-speed train heads.
PDE-based modeling methods include analytical and numerical ones. Since accurate analytical methods are very difficult to obtain and can only deal with simple modeling problems like blending between two primary surfaces (Bloor and Wilson 1989), various approximate analytical methods have been developed such as the Fourier series method (Bloor and Wilson 1990;Athanasopoulos et al. 2009;Ugail et al. 1999) and weighted residual method You et al. 2004). However, these methods are applicable to specific boundary conditions only and cannot describe complicated surface shapes. In contrast, numerical methods are most effective in solving PDEs subjected to arbitrary boundary conditions. Although numerical methods require heavy computations, they are more powerful than accurate and approximate analytical methods in shape description of complicated 3D models. The popular numerical methods are the finite difference method (Wang et al. 2019;Du and Qin 2005) and finite element method (Brown et al. 1998).
In summary, the existing optimization methods of highspeed train heads perform unsatisfactorily in optimizing both global and local shapes. Especially, the parametric modeling algorithms used in the existing optimization methods cannot accurately and completely describe global shape and local details of high-speed train heads. In order to address this problem, we will develop a multi-objective aerodynamic optimization design method of high-speed train heads by taking full advantage of the PDE-based modeling technology.

Overview of the design flow
As shown in Fig. 2, the design flow of our method can be basically divided into three parts: (1) PDE-based parametric modeling, (2) global optimization, and (3) local optimization.
In the first part, i.e., PDE-based parametric modeling, the half of a high-speed train head is first decomposed into 17 simple parts and the boundaries of each part are represented by 10 contour lines. After initializing the design variables which control the shape of contour lines, the framework of the train head is obtained as shown in Fig. 4a. Then, each part is described by a PDE surface patch generated from the finite difference solution of a vector-valued fourth-order PDE given in (1). Since any two adjacent PDE surface patches share the same boundary conditions, all the PDE surface patches are automatically and smoothly stitched together to represent the PDE surface model of the high-speed train head.
The second part is the global optimization process which contains the first part. By inputting the PDE surface model of the train head into the CFD simulation, we can obtain the aerodynamic drag coefficient of the whole train (C d ) and lift coefficient of the train tail (C l ) which are two optimization objectives in NSGA-II. The objective function of NSGA-II is to minimize the C d and |C l |. There are 10 design variables which control the shapes of contour lines of the train head. The criterion of the global optimization convergence is that the mean error of the non-dominated front F 1 between two successive generations is less than or equal to a threshold value ε. If the global optimization does not converge, the values of the design variables will be updated and inputted into the first part to create a new framework of the train head. Otherwise, we will obtain the globally optimized solution, i.e., the optimized framework of the train head.
Based on the optimized framework of the train head from the second part, the third part is to optimize local shapes of all PDE surface patches in the CFD simulation. In the local optimization process, only one PDE surface patch is selected and optimized each time. Therefore, only three shape control parameters of the PDE surface patch are taken to be the design variables and set to some initial values. By using the LHS method, 125 sample sets are created to generate various shapes of the PDE surface patch. Then, the corresponding train heads with these shapes of the Fig. 2 Overall design flow PDE surface patch are inputted into the CFD simulation. If the aerodynamic performance of C d and C l are improved, the PDE surface patch will be replaced by the new shape. Otherwise, the original shape of the PDE surface patch will be kept.
Since there are 17 PDE surface patches of the train head, the above calculations are repeated 17 times until all the PDE surface patches have been optimized, and then we obtain the locally optimized solution, i.e., the final optimized shape of the train head.
The global and local optimizations can be conducted simultaneously. However, the design variables of the simultaneous optimization will include 10 parameters of the contour lines and 3 control parameters of 17 PDE surface patches. The total number of design variables is 61, which cannot be simultaneously computed due to our limited computational resources, and it is also more difficult to find optimal results. Due to this reason, we conduct global and local optimizations respectively in this paper.

PDE-based parametric modeling method
The PDE-based parametric modeling method provides an efficient solution for the optimization of complicated geometry. Since the high-speed train head is a symmetrical structure, there is no need to represent the whole shape of the train head with PDE surface patches. Thus, we can represent the half of the train head with PDE surface patches, and then obtain the whole shape of the PDE surface-represented train head by reflecting the half of the train head through its symmetry plane.

The fourth-order PDE and the finite difference method
The mathematical model of PDE-based parametric modeling can take several forms due to different orders, i.e., the second, fourth and sixth orders (Castro et al. 2008). Since a second-order PDE cannot guarantee the tangent continuity between two adjacent PDE surface patches, and a sixthorder PDE is difficult to be solved and has low calculation efficiency, we select a fourth-order PDE as the mathematical model of surface patches of the train head which well balances performance and efficiency (Wang et al. 2019). The vector-valued fourth-order PDE is defined as: T is a vectorvalued position function which represents the generated parametric surface, a 1 , a 2 and a 3 are three shape control parameters, and u and v are the parametric variables of describing the parametric location of a point (x, y, z) on a PDE surface patch which are defined by u ∈ [0, 1] and v ∈ [0, 1]. The methods of solving a PDE can be analytical or numerical. Analytical methods represent a PDE surface patch by a vector-valued continuous function which can be calculated rapidly and precisely. However, analytical methods are usually applicable to a low-order or simple PDE and difficult to solve the PDE for the surface patches defined by four boundaries. Compared with analytical methods, numerical methods are more powerful in solving a high-order PDE and can deal with various complicated surface modeling problems although numerical methods are computationally more complicated and expensive. Therefore, numerical methods are more suitable in solving the PDE shown in (1). One popular numerical method is the finite difference method due to its simplicity, intuitiveness, and high efficiency in solving a vector-valued fourth-order PDE to create surface patches with four boundaries (Wang et al. 2019).
The finite difference method is to discretize the parametric region of a surface into I ×J regularly and uniformly distributed nodes as shown in Fig. 3. The small dots, squares and triangles represent the unknown inner nodes, the known boundary nodes and the virtual nodes beyond the boundaries of a PDE surface patch, respectively. The virtual nodes will be involved in the finite difference equations of the first partial derivatives on the boundaries and used to guarantee the boundary tangent continuity.
In Fig. 3, f i,j = f(u j , v i ) (i = 2, 3, . . . , I − 1, j = 2, 3, . . . , J − 1) represents an arbitrary inner node (i, j ) on the finite difference grid. Based on the Taylor series expansion of function f(u, v), the central difference approximation of ∂f i,j ∂u and ∂f i,j ∂v can be derived as follows: where h denotes the grid interval. The central difference approximations of the fourth partial derivatives derived from (2) and (3) can be expressed as: After substituting (4), (5) and (6) into (1), we obtain the following linear algebra equation at the inner node (i, j ): According to the central difference approximations defined by (4), (5) and (6), the finite difference equations at the inner nodes next to the four boundaries of the grid involve the boundary nodes and the virtual nodes. All the boundary nodes are known. The virtual nodes can be determined from the known boundary tangents, i.e., first partial derivatives of the function f(u, v) with respect to the parametric variable u or v on the four boundaries of the PDE surface patch. After merging these boundary conditions into the finite difference equations, Equation (7) can be written in the following matrix form: By using the PDE-based parametric modeling method, we build an original simplified high-speed train head model according to some practical constraints, such as the appropriate slenderness ratio, the appropriate space in the driving cab, and the good driver's perspective. According to the shape changes, we divide half of the high-speed train head into 17 parts as shown in Fig. 4a. Each part is represented by one PDE surface patch as shown in Fig. 4b. All PDE surface patches are generated by (8) and controlled by boundary curves, boundary tangents and three shape control parameters while the numerical solution (8) of (1) gives the approximated values of the grid nodes for every PDE surface patch.

Global deformation controlled by ten design parameters
We first optimize the whole shape of the high-speed train head by proposing a global deformation method. In order to achieve high efficiency of shape optimization of the highspeed train head, we choose ten contour lines, i.e., A 1 A 6 , B 1 B 6 , C 1 C 6 , D 1 D 6 , A 2 D 2 , A 3 D 3 , A 4 D 4 , A 5 D 5 , C 1 D 1 and A 1 E 1 C 1 to define half of the high-speed train head Fig. 4 The half of the high-speed train head model. a Contour lines. b PDE surface patches as shown in Fig. 4a and reduce the design variables. Each of the contour lines consists of some line segments and each line segment is between two adjacent joint vertices. For example, the contour line A 4 D 4 consists of three line segments A 4 B 4 , B 4 C 4 and C 4 D 4 . Since the two contour lines A 1 C 1 and A 6 D 6 decide the key position information of the train nose cone and coach, respectively, we fix their shapes to ensure the smooth transition between different train parts.
In order to further reduce the number of the total design variables in the optimization process, we use one design parameter to describe the shape of one contour line through the following equation: where the subscript i indicates the ith contour line, N i indicates the number of the total points on the ith contour line,L in and L in are the vector-valued position of the nth point on the ith contour line before and after the deformation, and δ i is the design parameter which control the shape of the ith contour line. By applying (9) to deform a contour line, the deformation becomes bigger and bigger when moving from two end vertices to the center of the contour line, such as the red contour line A 1 A 6 shown in Fig. 5a. Since each of the contour lines of the train head model is connected to some other line segments, the deformation of each of the contour lines will impact the shapes of these line segments. For example, when A 1 A 6 in Fig. 5a is deformed, the line segments A 2 B 2 , A 3 B 3 , A 4 B 4 and A 5 B 5 which are connected to A 1 A 6 are deformed accordingly. We employ the following deformation algorithm to describe the deformation of each connected line segment: We take the contour lines A 1 A 6 and A 1 E 1 C 1 as an example to show the effect of the contour line deformation on the shape of the high-speed train head in Fig. 5. In Fig. 5a, the contour lines in black and red present the original and deformed shapes, respectively. In Fig. 5b, the blue meshes represent the PDE surface model of the original half train head and the transparent surface represents the deformed shape after adjusting A 1 A 6 and A 1 E 1 C 1 by using (9) and (10). Table 1 shows ten design parameters, the corresponding contour lines and deformation directions. By adjusting the values of ten design parameters, the purpose of controlling the global deformation of the whole high-speed train head model with fewer design parameters is achieved.
The process of the global deformation consists of five steps. First, ten design parameters δ 1 , δ 2 , . . . , δ 10 are generated randomly within their design spaces. Second, one contour line is randomly selected and its deformation is determined by introducing the design parameter of the selected contour line into (9). Third, the deformations of all the line segments connected to the deformed contour line are determined by (10). Fourth, the shapes of all contour lines are updated with the deformations obtained in the second and third steps. Fifth, another contour line is selected and the second, third and fifth steps are repeated until all the ten contour lines are deformed.
Using the ten design parameters to control the global deformation of the high-speed train head greatly improves the efficiency in the first stage of the multi-objective Fig. 5 The effect of the contour line deformation on the shape of the half train head. a The deformation of the contour lines A 1 A 6 and A 1 E 1 C 1 and their connected contour lines. b The PDE surface model before and after adjusting the two contour lines optimization process of the high-speed train head. After optimizing the two optimization objectives, i.e., aerodynamic drag and aerodynamic lift, the optimized values of the ten design parameters are found and the globally optimized framework of the high-speed train head is obtained.

Local deformation controlled by three shape control parameters
After obtaining the optimized framework of the high-speed train head, its surface model can be generated using (8) by filling 17 PDE surface patches into the framework as shown in Fig. 4b. For acquiring a better hydro-mechanical property, we will further adjust the shape of the high-speed train head by deforming the local shapes of PDE surface patches in the second stage. As explained previously, the shape of each PDE surface patch is controlled by the three shape control parameters in PDE (1), i.e., a 1 , a 2 and a 3 .
With different values of the three shape control parameters, we obtain different surface shapes. For example, we set the three shape control parameters a 1 , a 2 and a 3 for the PDE surface patch A 3 A 4 B 4 B 3 to two groups of different values: Table 1 The design parameters and deformation directions of ten contour lines

Design parameter
Contour line Deformation direction 49, a 2 = −2.49, a 3 = −1.68 and a 1 = −0.13, a 2 = −0.13, a 3 = −1.68, and we obtain different surface shapes shown in Fig. 6. Using the three shape control parameters in PDE (1) to control the deformation of all PDE surface patches, we can further adjust the local shape of the high-speed train head which improves the surface quality and optimization results in the second stage. After the local shape optimization, the ultimate optimized shape of the high-speed train head is obtained.

NSGA-II
In the first stage of the multi-objective optimization design process, NSGA-II is applied to obtain the optimized framework of the high-speed train head. NSGA-II is an evolutionary multi-objective optimization algorithm and it is proposed by Deb et al. (2002). This algorithm is suitable for solving complex multi-objective optimization problems and has fast and accurate search performance. Fast nondominated sorting approach with elite strategy is used in NSGA-II, which greatly improves the sorting speed. Moreover, the use of the elite strategy ensures that the good solution will not be discarded. The implementation process of NSGA-II is shown in Fig. 7.
The NSGA-II algorithm includes five steps. First, the initial population P t with a size N is randomly generated and non-dominated sorting is performed. Second, the selection, crossover, and mutation operations are applied on population P t to generate an offspring population Q t with the same size N. Third, the two populations P t and Q t are combined to form the new population R t with a size 2N. Then, the combined population R t is sorted based on the non-dominated sorting approach to get non-dominated front F i (i = 1, 2, 3, . . .). Meanwhile, the crowding distance of each individual in F i is calculated. Fourth, according to the Fig. 6 The influence of the three shape control parameters on the shape deformation of the PDE surface patch A 3 A 4 B 4 B 3 . a The PDE surface patch obtained from a 1 = −2.49, a 2 = −2.49, a 3 = −1.68; b The PDE surface patch obtained from a 1 = −0.13, a 2 = −0.13, a 3 = −1.68; c The shape comparison between (a) and (b) order of i from small to large, F i is added into the next generation populationP t . When the addition of a certain F i causes the size ofP t to exceed the population size N, individuals in the F i will be added intoP t according to the crowding distance in descending order instead of adding the whole F i intoP t . Fifth, if the termination condition, i.e, the number of iterations, reaches its maximum, the procedure ends. Otherwise, theP t is set as the initial population and restart the first step.

Latin hypercube sampling
In the second stage of the multi-objective aerodynamic optimization process, we look for the local optimization Fig. 7 The implementation process of NSGA-II by adjusting the three shape control parameters in PDE (1) to improve the local shapes of PDE surface patches of the high-speed train head. In order to obtain the applicable values of the three shape control parameters, LHS method is applied to initialize them in the defined design space. LHS is a stratified sampling technology to approximate the random sampling from multiple parameter distributions (McKay et al. 2000), and it is the generalization of Latin square to multi-dimensions and each axis-aligned hyperplane contains only one sample which ensures all portions of the sample space are sampled and improves the sampling accuracy.
The three shape control parameters in PDE (1) are taken as three different input variables in the LHS process. In order to avoid an overlarge deformation of the PDE surface patch, we define the design space of each variable as [−3, 0.1] after a dozen experiments. The range of each variable is divided into 5 intervals with equal marginal probability, and a sample is randomly selected from each interval. Since each variable generates 5 samples, we can obtain totally 125 combinations of the three variables. Each combination including a sample set of a 1 , a 2 and a 3 can construct a different shape of the PDE surface patch which results in a different high-speed train head. Therefore, 125 high-speed train models are automatically generated and applied to the CFD simulation in the local optimization process. The LHS result of the three shape control parameters is shown in Fig. 8.

CFD simulation method
In this paper, we conduct the CFD simulation by using two software products: ICEM and FLUENT. ICEM is used to divide the whole high-speed train and the computational domain into aerodynamic meshes as shown in Fig. 10, and FLUENT is applied to carry out the subsequent fluid dynamics analysis for the high-speed train. To achieve automation of the optimization design, the process of the mesh generation and the aerodynamic calculation of the high-speed train are executed automatically by invoking the script files of ICEM and FLUENT  a 1 , a 2 and a 3 respectively. These script files can be performed by batch commands.

The model and the computational domain
In order to evaluate the practical aerodynamic performance of the high-speed train head, we construct a whole train model by adding a middle coach and a tail (the same as the head) to the train head. Since the real high-speed train has a complex shape which may greatly increase the time consumption of the parametric modeling, mesh generation and CFD simulation in the optimization process, and we mainly study the relationship between the shape of the train head and the aerodynamic drag and lift forces, the constructed high-speed train model is simplified by ignoring bogies and other auxiliary structures, as shown in Fig. 9a. In this paper, the generated high-speed train model is on the same scale of a real high-speed train which runs at the speed of 300km/h in the open air without a crosswind.
We first construct a computational domain to simulate the flow field around the high-speed train. We indicate the train length with the symbol L (L = 78 m) and take it as a characteristic length. The height and width of the train is 0.04435L and 0.04470L, respectively. Then, we determine the size of the computational domain according to the characteristic length L. As shown in Fig. 9b, the distance between the entrance of the computational domain and the nose cone of the train head is 1L, and the distance between the nose cone of the train tail and the exit of the computational domain is 1.5L. The distance from the ground to the top boundary of the computational domain is 0.5L, and the distance from the train center to the boundary on both sides of the computational domain is 0.5L. The distance between the train wheel and the ground is 0.00235L which represents the height of rail tracks.

CFD simulation
The software FLUENT is adopted to carry out the CFD simulation after obtaining the mesh file generated from the ICEM, and calculate the coefficients of aerodynamic drag and lift of the high-speed train model. Note that we have not yet conducted the wind tunnel or real vehicle experiments and have no available experimental data. In order to ensure that our simulation results are reliable, our strategy of the simulation setup makes reference to previous similar work (Li et al. 2016;Zhang et al. 2018;Morden et al. 2015;Yao et al. 2014;Cheli et al. 2010), which has indicated that the CFD computation is able to predict the flow correctly.
Governing equations of fluid flows The Navier-Stokes (N-S) equations are the governing equations of fluid flows and have different forms for incompressible and compressible flows (Cebeci et al. 2005). In this paper, the flow around the high-speed train is considered to be an incompressible flow in the CFD simulation based on the following reasons.
(1) The speed of the high-speed train is 300 km/h and the resultant Mach number is 0.245. (2) When the highspeed train is in the open air without passing other trains or going through a tunnel, the impact of the air density on the flow can be ignored. (3) The numerical simulation using an incompressible flow solver is accurate compared with experimental data (Cheli et al. 2010;Morden et al. 2015).
(4) The incompressible flow solver is commonly used in the studies of the aerodynamic optimization of high-speed train heads (Sun et al. 2010;Li et al. 2016;Zhang et al. 2018). The incompressible N-S equations take the following tensor forms: where u i is the velocity components in the i direction, ρ is the density, p is the pressure, and v is the fluid kinematic viscosity.

Approximate approaches for N-S equations
The high-speed train has a large Reynolds number so that the flow of the high-speed train is highly turbulent (Wang et al. 2017). Since the high turbulent flow has fluctuations of pressure, temperature and velocity over a wide range of frequencies, solving the incompressible N-S equations is a formidable challenge (Cebeci et al. 2005). The commonly used approximate computational approaches include direct numerical simulation (DNS), large eddy simulation (LES) and Reynolds-averaged Navier-Stokes (RANS) equations. Compared with DNS and LES, RANS has lower computational cost. Besides, RANS is also good at accurately predicting the pressure distribution and the frictional resistance (Bensow et al. 2006). Therefore, since we focus on studying the drag and lift of the high-speed train and have less interest in the flow details, we apply the RANS approach in this paper. The RANS equations apply the Reynolds decomposition on the instantaneous incompressible N-S equations (11) and (12), which splits the flow into its mean and fluctuating components (Morden et al. 2015). The incompressible RANS equations are defined as: where τ ij = u i u j is the Reynolds stress tensor which cannot be formally expressed in terms of mean flow variables and thus a turbulent model is used to close the equations.

Turbulent model
The frequently used turbulent models include k − ε, k − ω, and shear stress transport (SST) k − ω models. We select the SST k − ω model in this paper due to the following reasons.
(1) The SST k − ω model blends the advantages of classical k − ω and k − ε models, and can better model flows on the smooth surfaces of a high-speed train (Wang et al. 2017). (2) The SST k − ω model is recommended as the optimal RANS model based on experimental verifications (Morden et al. 2015). (3) The SST k − ω model is commonly selected as the turbulence model in the studies of aerodynamic optimization of highspeed train heads (Yao et al. 2014;Yao et al. 2016;Zhang et al. 2018) Wall function Since the gradients of velocity near the wall is steep, we need to use a large number of thin meshes to accurately capture the gradients. However, these thin meshes will result in poor mesh quality and high computational cost. In order to solve this problem, we can use a single large mesh instead of many thin meshes plus a nonlinear function called the wall function to simulate the gradient variation. In this paper, we adopt the standard wall function because it works reasonably well for a broad range of wall-bounded flows (ANSYS I 2009) and has been most widely used in the CFD simulation of high-speed train heads (Yao et al. 2014;Muñoz-Paniagua et al. 2014;Yao et al. 2015;Li et al. 2016;Zhang et al. 2018). The standard wall function is defined as (ANSYS I 2009): U * = y * , y * < 11.225 1 κ ln(Ey * ), y * > 11.225 (15) where U * is the dimensionless velocity, y * is the dimensionless wall distance, κ = 0.4187 is the von Kármán constant, and E = 9.793 is the empirical constant.
Other CFD setup strategies We use the pressure-based segregated solver, and the SIMPLE scheme is used to couple the pressure and velocity. In the spatial discretization scheme, we apply the least squares cell based gradient with the second-order interpolation such as the secondorder pressure and the second-order upwind momentum. As shown in Fig. 9b, the left and right of the computational domain are set as the velocity inlet and pressure outlet boundaries, respectively, and the two sides and the top are set as the symmetric boundary. The high-speed train surface is set as the non-slip wall boundary conditions, and the ground is set as the slip wall boundary conditions to simulate the ground effect and the slip velocity is equal to the speed of the train. In addition, we define that the farfield pressure is 1 atm, the temperature is 288 K, and the reference area is the maximum cross-sectional area of the train which is 12.0638 m 2 .
In previous aerodynamic research work of high-speed trains, Cheli et al. (2010) and Morden et al. (2015) have conducted wind tunnel tests to verify their CFD simulation results. For example, Cheli et al. (2010) studied the static aerodynamic coefficients of a high-speed train with different wind angles, and showed the lateral and vertical force coefficients from the simulation and experiment are very close when the wind angle is below 20 • . In the work of  Morden et al. (2015), the C d and C l from simulation and experiment results are 0.14 against 0.13 and 0.17 against 0.23, respectively, which show a good agreement. The CFD setup strategies of their verified simulation models give us a solid reference. Therefore, our CFD setup strategy is similar to their work as shown in Table 2, and it ensures that our results are as reliable as possible, though we have no experiment data support.

Mesh generation
When conducting a specific division of spatial meshes in ICEM, the quality and quantity of meshes have a significant influence on the computational efficiency, astringency and precision of the CFD simulation results. We apply Hexahedral meshes to divide the whole computational domain and distribute prism meshes in the area around the body of the high-speed train. Since the train head bears most of the aerodynamic drag in simulation process, the mesh division around the train head is refined for improving the computational accuracy. In addition, in order to reduce computational cost, the mesh size of the middle coach is slightly larger than that of the train head and tail. The area around the body of the whole train utilizes five layers of fine prism meshes to accurately simulate the flow field around the train body. The y + is non-dimensional wall distance of the first cell from the wall based on the fluid local velocity. The range of the y + values obtained around the train body is from 5.99 to 77.91. Since it is difficult to get all surfaces of the train model to have the desired y + value exactly, we have made the majority of y + within the recommended ranges. Note that the y * in (15) plays the same role as y + in measuring the dimensionless distance. The difference is that the velocity scale of y * is based on the turbulent kinetic energy and the velocity scale of y + is based on the wall shear stress (ANSYS I 2009). The number of cell volumes within the boundary layer region around the train body is about 2.1 million. The mesh details are shown in Fig. 10.
Since too many meshes will increase the simulation time and too few meshes will reduce the computational accuracy in CFD simulation process, it is important to generate a proper number of meshes (Yao et al. 2016). We build four sets of meshes with different mesh quantities, i.e., 7.2 million, 9.8 million, 14.3 million, and 22.2 million, to evaluate the influence of different meshes on the aerodynamic drag and lift. Table 3 shows the results of C d and |C l | of the four sets of meshes. The values of C d and  |C l | obtained from the second set of meshes are 10.04%, 9.49% larger and 5.99%, 9.25% smaller than those obtained from the third and fourth sets of meshes, respectively. Due to the high number of simulations needed to be carried out for the optimization algorithm, which require a large amount of computational cost, we make a compromise between the accuracy and the computational cost, and accept a discrepancy of approximated 10%. Therefore, the second set of meshes with 9.8 million meshes is adopted for all flow field calculations in this paper.
We also carry out a small scale computation for 22.2 million meshes following the same optimization flow of 9.8 million meshes. The optimized train head with 22.2 million meshes has a smaller drag and a larger lift compared with 9.8 million meshes, which is in accordance with Table 3, and the error of optimized C d between 9.8 and 22.2 million meshes is within 10%, which indicate that the acceptable error of 10% covers the error caused by 9.8 million meshes.

Results and discussion
The multi-objective optimization process of the highspeed train head includes the global optimization and local optimization stages. We first get the optimized framework of the train head in the global optimization stage (Section 7.1), then optimize the local shape of PDE surface patches on the train head with the framework and finally obtain the optimized shape of the train head in the local optimization stage (Section 7.2). The aerodynamic performance of the train head in different stages will be compared and analyzed in Section 7.3.

Global optimization
In NSGA-II, C d and C l are set as the two optimization objectives. The population size is set to be 40 (N = 40) and the number of generations is set to be 10. Moreover, the crossover and mutation probability are set as 0.9 and 0.1, respectively. The threshold value of the convergence criterion is set as ε = 0.01. The ranges of the ten design variables controlling the deformation of the train head framework are shown in Table 4. In order to avoid the distortion of the train head shape, the suitable lower and By using NSGA-II, we obtain 400 solutions in which there are 13 Pareto-optimal solutions constructing a Paretooptimal front as shown in Fig. 11. In order to find the most satisfactory solution among Pareto-optimal solutions and inspired by the minimum distance algorithm discussed in Li et al. (2016), we develop a cost function which assigns a proper weight to each objective and aggregates all of the objectives together, defined by: where i represents the ith Pareto-optimal solution, S is the total number of Pareto-optimal solutions, f C d (i) and f C l (i) are the ith Pareto-optimal solutions of C d and C l , respectively, min f C d = min f C d (1), f C d (2), ..., f C d (S) , min |f C l | = min |f C l (1)|, |f C l (2)|, ..., |f C l (S)| , and ω c is the user-defined weight. Fig. 11 The Pareto-optimal front in the global optimization stage By applying different value of ω c , we can obtain different global optimized results from the 13 Pareto-optimal solutions. Since C d plays a more important role in reducing the air resistance to the forward motion of a train compared with C l , we set ω c = 0.001 so that C d is the main contributor based on previous work (Brockie and Baker 1990;Schetz 2001). The final solution corresponding to the minimum D is selected as the global optimization result which is shown as the red star marker in Fig. 11. The ten design parameters of the global solution construct the optimized framework of the high-speed train head which will be used in the next local optimization stage.

Local optimization
In order to demonstrate the process of the local optimization, we take the patch A 3 A 4 B 4 B 3 as an example. The design variables are the three shape control parameters of the patch A 3 A 4 B 4 B 3 . Like the global optimization, the optimization objectives for the local optimization are C d and C l . We apply LHS to sample the design variables into 125 sample sets and use a dominated sorting to select the acceptable solutions whose C d and C l are not both dominated by the results of the optimized framework of the high-speed train head. In addition, (16) is employed to determine the final solution from the sample sets. Since C d is the main contributor in the global optimization stage, we enhance the influence of C l in the local optimization stage by setting ω c = 0.5. The final solution is shown as the red star marker in Fig. 12.

Discussion
Since our proposed multi-objective optimization method based on the PDE parametric modeling includes global and local optimization stages, we compare the results from the two stages to discuss the advantages of our method. Figure 13 shows the original (a), the globally optimized (b) and the locally optimized (c) high-speed train head models. After the global optimization, there are distinct deformations on the train head compared (b) with (a). For example, the height of the cab decreases and the nose cone and the front spoiler move forward, as shown in (d). The  local optimization further optimizes PDE surface patches and the locally optimized train head is shown in (c). It can be seen clearly from (d) that the locally optimized shape of the patch A 3 A 4 B 4 B 3 representing the cab window on the train head becomes concave.
The results of the optimization objectives C d and C l of the original, the globally optimized and the locally optimized high-speed trains are shown in Table 5. Compared with the original train, the C d and |C l | of the global optimized train is reduced by 7.58% and 16.56%, respectively. After the local optimization of the globally optimized train head, the two optimization objectives are further reduced. Compared with the original train, the C d and |C l | of the locally optimized train is reduced by 7.90% and 38.85%, respectively.
Aerodynamic drag and lift forces are mainly caused by the pressure force which mainly exists on the surface of the train head and tail. The shape deformation of the train head has a direct impact on the pressure distributions (Li et al. 2016). In order to discuss the aerodynamic performance of the train before and after the multi-objective optimization, the pressure distributions of the train head and tail are presented in Fig. 14. Since the scales of the pressure distribution of the train head and tail are different, we use two color bars in different scales for the train head and tail to clearly indicate the changes before and after the optimization. The left column in Fig. 14 shows the pressure distributions of the original, globally optimized and locally optimized train heads. There are mainly three high-pressure zones near the nose cone, cab window and roof, respectively, which are indicated by the circle of red dashed lines. After the global optimization, the pressure near the cab window significantly decreases and the pressure variation near the roof, i.e., a negative pressure followed by a positive pressure, almost disappears, but the pressure near the nose cone is a little larger than that of the original train head. After the local optimization, the pressure near the nose cone is reduced in addition to further decrease of the pressure near the cab window. Similarly, the right column in Fig. 14 shows the pressure distributions of train tails. After the global optimization, the high-pressure zone in front of the nose cone is larger than that of the original train, which gives the train tail a forward and an upward push. Moreover, the pressure near the cab window is decreased in both global and local optimizations, and the region with negative pressure near the roof is significantly reduced after the global optimization and slight improved after the local optimization. Through the multi-objective optimization process with the global and local optimization stages, the final optimized shape of the high-speed train head is obtained which has an improved pressure distribution and small aerodynamic drag and lift forces. Figure 15 shows the streamlines around the high-speed train before and after the multi-objective optimization. From the overall views and the closer views of the train heads and tails, the streamlines are smooth near the train heads, and a flow separation occurs near the front spoiler of the train tail, which produce a trailing wake vortex as shown in (a). After global and local optimizations, the wake vortex is reduced and the streamlines around the train tail become smoother, which enlarge the high-pressure zone in front of the nose cone and give a forward push to reduce the drag, as shown in Fig. 14. The shape of the high-speed train head has a significant effect on the drag. By optimizing the shape, the air flow around the train can be smoother and the drag can be reduced.
In order to demonstrate the advantage of our method, we compare our method with current shape optimization methods of high-speed train heads in terms of optimizing objectives C d and C l , as shown in Table 6. Note that our results in Table 6 base on the choice of ω c in (16), i.e., ω c = 0.001 in the global optimization and ω c = 0.5 in the local optimization, and different choice of ω c will produce different values of C d and C l as well as their reduction ratios. The high-speed trains in different studies are under the same running conditions, i.e., the high-speed train is in the open air without passing each other or going through a tunnel and the running speed is 300 km/h. The data in Table 6 indicate that although our original model already has good aerodynamic performance, i.e., the smallest values of C d and |C l | in comparison with existing optimization studies, our method can still achieve Fig. 15 Comparison of the streamlines around the original (a), globally optimized (b) and locally optimized (c) train heads (the left column), whole trains (the middle column) and train tails (the right column)

Table 6
Comparison of our method with other optimization methods Our method  optimal shape; RE, reduce (%); -, no data available the maximum reduction rates 7.9% and 38.85% of C d and |C l |, respectively, among all the optimization methods and improve the drag and lift further. Therefore, our method is more effective in improving the aerodynamic performance of high-speed train heads.

Conclusion
In this paper, a novel multi-objective aerodynamic optimization design process of a high-speed train head is proposed. The PDE-based parametric modeling method is applied to construct the parametric model of the high-speed train head, which can describe the complicated shape in detail with few design variables and keep the surface smooth. NSGA-II is adopted to obtain Pareto-optimal solutions in the global optimization stage of the high-speed train head and take the aerodynamic drag of the whole train and the aerodynamic lift of the train tail as the optimization objectives. Then, an optimized framework of the high-speed train head is selected from the Pareto-optimal solutions using an improved minimum distance algorithm. Based on the obtained optimized framework, LHS is introduced into the local optimization stage to obtain the final optimized shape of the train head by generating various sample sets of the three shape control parameters of PDE surface patches and optimizing the shape of each patch. We have demonstrated our method by analyzing the aerodynamic characteristics, pressure distributions and streamlines of the optimization solutions in both global and local optimization stages compared with the original high-speed train head, and evaluating the drag and lift coefficients compared with other optimization methods. The analysis results indicate that our method is able to optimize both global and local shapes and significantly improve the aerodynamic performance of the high-speed train head.
In the future, we wish to research multi-objective aerodynamic optimization of high-speed train heads under more complicated running environments such as crosswinds and tunnels. We also plan to develop an analytical PDEbased modeling method which can create more accurate and complex shape deformations of train heads more efficiently.
Funding This research is supported by the PDE-GIR project which has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 778035, and the National Natural Science Foundation of China No. 51475394.

Conflict of interest
The authors declare that they have no conflict of interest.

Replication of results
The paper contains sufficient details that is needed to replicate the results with limitations due to the heuristic nature of the applied evolutionary algorithms. The code and data of this study will be made available by request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.