# Tabu Genetic Cat Swarm Algorithm Analysis of Optimization Arrangement on Mistuned Blades Based on CUDA

Authors: Yi Cai; Junjie Gu; Honggang Pan; Hongyuan Zhang; Tianyu Zhao
Journal: Shock and Vibration
 (2021)

---

## Abstract

Tabu genetic cat swarm optimization algorithm is proposed for optimization arrangement on mistuned blades. Furthermore, it is improved to be an innovative parallel algorithm based on Compute Unified Device Architecture (CUDA), whose performance is analyzed both in continuous and discrete solution space. The lumped parameter model and finite element model of the bladed-disk system are established for dynamics analysis and optimization verification. The complete mistuned parameter identification method and the improved mixed-interface prestressed substructure modal synthesis method are two other highlights. The results indicate that the algorithm in this paper has the advantages of low cost and high efficiency. The vibration localization and amplitude of the mistuned bladed-disk system are both reduced significantly. The optimization analysis method is proved to be applicable to the mistuned bladed-disk system of the aeroengine in engineering practice.


---

## Body

## 1. Introduction

Aeroengine is known as “the crown of mechanical industry.” The vibration failure accounts for over 60% of the total. Furthermore, more than 70% of these are caused by blades. The unbalanced vibration by reason of uneven vibration energy distribution is one of the main causes [1–10]. The bladed-disk system of aeroengine is a circular symmetric structure. The physical properties and geometric conditions of each blade are identical. The modes are uniformly distributed around the circle. There is a small amount of difference in each sector because of the manufacturing error, material property, uneven wear, and so on. This small difference is called mistuning in structural dynamics. It breaks the circulatory symmetry and changes the dynamic characteristics of the bladed-disk system. As a result, the vibration responses are localized and increased. This results in high cycle fatigue failure of blades. The different arrangements on mistuned blades could cause different forced vibration amplitudes. A suitable arrangement is much needed when replacing or installing aeroengine blades. It has been a difficult problem in the production and maintenance engineering of the aeroengine.In recent years, scholars have extensively carried out numerical simulation and theoretical and experimental studies on dynamic characteristics of the mistuned bladed-disk system. Petrov et al. [11] and Petrov and Ewins [12] presented a new method for dynamic analysis of mistuned bladed disks based on the exact relationship between tuned and mistuned systems and solved the optimization problem of determining the worst mistuning patterns. Shahruz et al. [13, 14] defined a sensitivity matrix to determine the occurrence of vibration localization for mistuned periodic structures. Castanier and Pierre [15] and Lim et al. [16] reviewed the literature on reduced-order modeling, simulation, and analysis of the vibration of mistuned bladed disks found in gas-turbine engines. They developed several reduced-order modeling techniques to predict the mistuned forced response efficiently. Yan et al. [17] used component-mode synthesis methods to establish the vibration differential equation of the system for coupling vibration of the mistuned bladed disk in the aeroengine and analyzed the vibration mechanism. Salhi et al. [18] presented two methods for the identification of the vibration characteristics of rotating engine blades from time response measurements. Ganine et al. [19] reviewed the static mode compensation method to predict geometrical mistuning effects on the response of bladed disks and analyzed its limitations. Vargiu et al. [20] extended a preexisting reduction technique suitable for the analysis of mistuned bladed-disk dynamics, the component mode mistuning technique (CMM), in order to allow for the introduction of a sector frequency mistuning pattern. Khemiri et al. [21] studied the forced response of a bladed disk with damping mistuning. Wang et al. [22] presented an efficient method for the nonlinear vibration analysis of mistuned centrifugal impellers with crack damages. Zhou et al. [23] proposed an essentially nonlinear piezoelectric shunt circuit for the practical realization of nonlinear energy sink and then applied to a mistuned bladed disk for blade vibration reduction. Kim [24] investigated the frequency characteristics of natural whirl speeds associated with the in-plain vibration of a rotating mistuned bladed disk mounted on an isotropic support. Capiez-Lernout et al. [25] dealt with the dynamical analysis and uncertainty quantification of a mistuned industrial rotating integrally bladed disk, for which the operating regime under consideration took into account the nonlinear geometrical effects induced by large displacements and deformations. Yuan et al. [26] described an assessment of subset simulation (SubSim) techniques to increase the computational efficiency for the predictions of probabilistic dynamic behavior in mistuned bladed-disc systems. Ma et al. [27–29] developed a new dynamic model of rotor-blade systems. Jamia et al. [30] presented a method of mistuning identification in a bladed disk using the wavelet packet transform. Salas et al. [31] presented a unique way to address the reduced-order model, where each blisk sector was attached as individual substructures with the free-interface approach known as Craig–Chang. Cha [32] examined the performance of friction dampers of a geometric mistuned bladed-disk assembly under random excitations.It can be found from the above literature review that the vibration localization problem of mistuned bladed-disk system has been studied extensively [33–43]. Many research results have been obtained to be theoretical foundation for design and manufacture. However, an effective method that can be really applied directly for optimization arrangement on mistuned blades is lacking.This paper presents tabu genetic cat swarm optimization algorithm for optimization arrangement on mistuned blades. Furthermore, it is improved to be an innovative parallel algorithm based on CUDA. The lumped parameter model and finite element model of the bladed-disk system are established for dynamics analysis. The complete mistuned parameter identification method and the improved mixed-interface prestressed substructure modal synthesis method are proposed for optimization verification. It is of great significance for improving the reliability, reducing the failure rate, and increasing the interchangeability of blades in the aeroengine.


## 2. Dynamics Analysis of Mistuned Bladed-Disk System

### 2.1. Establishment of Lumped Parameter Model

In order to analyze the dynamic characteristics of the bladed-disk system, the lumped parameter model is established, as shown in Figure1. At the same time, the dimensionless parameters are used for being simple and convenient in this paper. As shown in Table 1, n is the number of blades, mb and kb are the quality and stiffness of the single blade, md and kd are the quality and stiffness of the single sector disk, and kt is coupling stiffness.Figure 1
Lumped parameter model of the bladed-disk system.Table 1
Lumped parameter of the bladed-disk system.
ParameternmbmdkbkdktDimension380.11 kg46.86 kg4.5e5 N/m4.95e5 N/m2.22e8 N/mDimensionless38142611.1493Forced vibration equation of the bladed-disk system can be derived:(1)Mq¨+Cq˙+Kq=F,where M, C, and K are mass matrix, damping matrix, and stiffness matrix and q and F are the displacement vector and exciting force vector. The detailed expressions are as follows:(2)M=mb100⋯00md10⋯0000⋯⋮⋮⋮⋮mbn00⋯00mdn,C=cb1−cb100⋅⋅0−cb1cb100⋅⋅000cb2−cb2⋅⋅000−cb2cb2⋅⋅0⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅cbn−cbn0000⋅−cbncbn,K=kb1−kb100⋯0−kb1kb1+kd1+2kt0−kt⋯−kt00kb20⋯⋮⋮⋮⋮⋮kbn−kbn0−kt0⋯−kbnkbn+kdn+2kt,q=x1y1⋯xnynT,F=f1t,0,⋯fnt,0=eiωt+θ0,0,⋯eiωt+θn−1,0,where the expression of the phase angle is(3)θj=2πEjn,j=0,1,…,n−1.ω and E are the frequency and order of exciting force. It is assumed that(4)qt=Aeiωt.Substituting equation (4) into equation (1), it can be obtained that(5)A=−ω2M+iωC+K−1eiθj.


### 2.2. Dynamics Analysis of Mistuned System

The vibration energy of the tuned system should be transferred to the whole structure uniformly in a bladed-disk system. However, mistuning could break this balance. The typical case is mistuning which is caused by small physical differences in stiffness of blades. The mistuning of stiffness is introduced randomly with normal distribution. Three standard deviations of mistuning used in this paper are 1%, 3%, and 5%. They are defined as Mistuned I, Mistuned II, and Mistuned III. The detailed mistuning is shown in Figure2.Figure 2
Stiffness mistuning of blades. (a) Mistuned I. (b) Mistuned II. (c) Mistuned III.
(a)(b)(c)The forced vibration response of the bladed-disk system can be obtained as follows based on equation (4), where the order of exciting force E is equal to the pitch diameter number 6 when system frequency shifts because, in this case, the coupling vibration between the blade and disk is significant. At the same time, the sensitivity of the mistuned parameter is the greatest.It can be found from Figure3 that all the amplitude-frequency characteristic curves are the same. There are two distinct resonance peaks. On the contrary, the amplitude-frequency characteristic curves of the mistuned system are very different. There are many peaks besides two main resonance peaks. In addition, there are certain differences among the maximum amplitudes of Mistuned I, Mistuned II, and Mistuned III. The maximum amplitude-frequency characteristic curves can be obtained by projecting waterfall plot onto the plane.Figure 3
Amplitude-frequency characteristic curve and max-amplitude curve. (a) Tuned system. (b) Mistuned I. (c) Mistuned II. (d) Mistuned III.
(a)(b)(c)(d)It can be seen from Figure4(a) that the resonance peaks increase with the standard deviations of mistuning, while the resonance regions get wider. The resonant frequencies of the mistuned system have been dispersed obviously when the standard deviation of mistuning reaches to be 5%. There are many resonance peaks in this case. It can be considered that the vibration amplitudes increase as a whole. It can be found further from Figure 4(b) that there is a minimum amplitude between the resonance regions. Moreover, the standard deviation of mistuning has almost no effect on it. Therefore, we can make the working frequency of the mistuned bladed-disk system to be near 1.026. From this, the amplitudes and the sensitivities of mistuning can be reduced effectively.Figure 4
Maximum amplitude-frequency characteristic curves of different mistuned cases.
(a)(b)


## 2.1. Establishment of Lumped Parameter Model

In order to analyze the dynamic characteristics of the bladed-disk system, the lumped parameter model is established, as shown in Figure1. At the same time, the dimensionless parameters are used for being simple and convenient in this paper. As shown in Table 1, n is the number of blades, mb and kb are the quality and stiffness of the single blade, md and kd are the quality and stiffness of the single sector disk, and kt is coupling stiffness.Figure 1
Lumped parameter model of the bladed-disk system.Table 1
Lumped parameter of the bladed-disk system.
ParameternmbmdkbkdktDimension380.11 kg46.86 kg4.5e5 N/m4.95e5 N/m2.22e8 N/mDimensionless38142611.1493Forced vibration equation of the bladed-disk system can be derived:(1)Mq¨+Cq˙+Kq=F,where M, C, and K are mass matrix, damping matrix, and stiffness matrix and q and F are the displacement vector and exciting force vector. The detailed expressions are as follows:(2)M=mb100⋯00md10⋯0000⋯⋮⋮⋮⋮mbn00⋯00mdn,C=cb1−cb100⋅⋅0−cb1cb100⋅⋅000cb2−cb2⋅⋅000−cb2cb2⋅⋅0⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅cbn−cbn0000⋅−cbncbn,K=kb1−kb100⋯0−kb1kb1+kd1+2kt0−kt⋯−kt00kb20⋯⋮⋮⋮⋮⋮kbn−kbn0−kt0⋯−kbnkbn+kdn+2kt,q=x1y1⋯xnynT,F=f1t,0,⋯fnt,0=eiωt+θ0,0,⋯eiωt+θn−1,0,where the expression of the phase angle is(3)θj=2πEjn,j=0,1,…,n−1.ω and E are the frequency and order of exciting force. It is assumed that(4)qt=Aeiωt.Substituting equation (4) into equation (1), it can be obtained that(5)A=−ω2M+iωC+K−1eiθj.


## 2.2. Dynamics Analysis of Mistuned System

The vibration energy of the tuned system should be transferred to the whole structure uniformly in a bladed-disk system. However, mistuning could break this balance. The typical case is mistuning which is caused by small physical differences in stiffness of blades. The mistuning of stiffness is introduced randomly with normal distribution. Three standard deviations of mistuning used in this paper are 1%, 3%, and 5%. They are defined as Mistuned I, Mistuned II, and Mistuned III. The detailed mistuning is shown in Figure2.Figure 2
Stiffness mistuning of blades. (a) Mistuned I. (b) Mistuned II. (c) Mistuned III.
(a)(b)(c)The forced vibration response of the bladed-disk system can be obtained as follows based on equation (4), where the order of exciting force E is equal to the pitch diameter number 6 when system frequency shifts because, in this case, the coupling vibration between the blade and disk is significant. At the same time, the sensitivity of the mistuned parameter is the greatest.It can be found from Figure3 that all the amplitude-frequency characteristic curves are the same. There are two distinct resonance peaks. On the contrary, the amplitude-frequency characteristic curves of the mistuned system are very different. There are many peaks besides two main resonance peaks. In addition, there are certain differences among the maximum amplitudes of Mistuned I, Mistuned II, and Mistuned III. The maximum amplitude-frequency characteristic curves can be obtained by projecting waterfall plot onto the plane.Figure 3
Amplitude-frequency characteristic curve and max-amplitude curve. (a) Tuned system. (b) Mistuned I. (c) Mistuned II. (d) Mistuned III.
(a)(b)(c)(d)It can be seen from Figure4(a) that the resonance peaks increase with the standard deviations of mistuning, while the resonance regions get wider. The resonant frequencies of the mistuned system have been dispersed obviously when the standard deviation of mistuning reaches to be 5%. There are many resonance peaks in this case. It can be considered that the vibration amplitudes increase as a whole. It can be found further from Figure 4(b) that there is a minimum amplitude between the resonance regions. Moreover, the standard deviation of mistuning has almost no effect on it. Therefore, we can make the working frequency of the mistuned bladed-disk system to be near 1.026. From this, the amplitudes and the sensitivities of mistuning can be reduced effectively.Figure 4
Maximum amplitude-frequency characteristic curves of different mistuned cases.
(a)(b)


## 3. Analysis of the Intelligent Optimization Algorithm

### 3.1. Tabu Genetic Cat Swarm Optimization Algorithm

The cat swarm optimization algorithm is considered to be a method which can solve those complex optimization problems. It is combining search behavior and tracking behavior of a cat. It has an obvious advantage in computation speed, convergence, and sensitivity to initial conditions. It can be seen that the cat is a feasible solution for the optimization problem which we are going to solve. The behavior of cats can be categorized into two modes. One is called the search mode. It says that the cats are lazy or look around. The other one is called the tracking mode. It says that the cats are tracking a dynamic target. The two modes interact with each other by using mixture ratio (MR). MR indicates the proportion of cats in the whole group which executes the tracking mode.The search mode is used to simulate a cat’s current state. It means resting, looking around, and searching for the next moving location. In view of optimization arrangement on mistuned blades, it converts to express that the new arrangement is generated randomly. The quality of the solution is improved very slowly as the arrangement may be generated repeatedly in the random process. From this, tabu list is introduced to be the short-term memory storage of the algorithm to prevent circular search.The tracking mode is used to simulate the case that a cat is tracking a target. The location of the cat is updated by changing the speed of every single dimension. Furthermore, the speed is changed by adding random disturbance. The heuristic algorithm is needed to avoid the optimization problem getting into local optimal for the typical combinatorial optimization arrangement on mistuned blades. The solutions are improved constantly by introducing the mutation operator and crossover operator of genetic algorithm in this paper so that the global optimal solution can be found in a short time. The change of the solution is combined with the speed and location update process. The update formulas are defined as follows:(6)Vit+1=Vit+r⊗Xbestt⊕xit,xit+1=xit+Vit+1,where Xbest(t) is the location of the cat which has the best fitness value in the current cat swarm and Vi(t) and xi(t) are the speed and location of the ith cat at time t.The speed is to change the location of the cat. It is defined as an exchange list. The position change of the cat is to do exchange processing withx by the speed V. r is a random number between 0 and 1. Its meaning is probability. A random number will be generated for Vi(t) when Vi+1(t) is calculated. The genetic operators are called if the random number is greater than or equal to r, otherwise the speed stays the same. It can make the algorithm to keep good mutant ability by introducing genetic algorithm. The algorithm can avoid premature convergence effectively. The crossover operator and mutation operator play a key role in exciting cat swarm and keeping population diversity.In order to explain the process of crossover and mutation, two arrangements on ten mistuned blades are generated randomly as follows:(7)9,5,3,7,10,1,4,2,6,8,10,8,2,6,1,5,4,3,9,7.It means that No.9 is installed in position 1 and so on. Crossover operator (⊕) is based on the integer crossing method. Two crossover locations are selected. Then, the individual is crossed with the optimal one. It is assumed that 2 and 4 were randomly selected for crossover positions. The specific operation is as follows:(8)Individual9,5¯,3,7¯,10,1,4,2,6,8Optimization10,8¯,2,6¯,1,5,4,3,9,7⟶⊕9,8¯,3,6¯,10,1,4,2,6,8.The new individual has repeat positions obviously. The adjustment method is to replace with the missing blades:(9)9,8,3,6,10,1,4,2,6,8⟶9,8,3,6,10,1,4,2,5,7.Mutation operator (⊗) is based on the method of interchanging two random positions in an individual. Two mutation positions are selected and exchanged. It is assumed that 3 and 6 were selected. The specific operation is as follows:(10)9,8,3¯,6,10,1¯,4,2,5,7⟶9,8,1¯,6,10,3¯,4,2,5,7.The new individuals will be kept only when they are better than the old ones. So, an accurate evaluation criterion is needed. The vibration amplitudes and vibration differences are considered comprehensively, and the fitness function is designed as follows:(11)L=C1meanX+C2maxXC3varX+C4,where X is a vector of the maximum amplitudes of blades, mean(X) and var(X) are average value and variance of the maximum amplitudes of blades, and C1∼C4 are constants which are determined by the mistuned parameter.In addition, there is a great change to be repeated for blade arrangement as the randomness of new positions in the search mode. As a tabu list has short-term memory, the solution that has been searched can be excluded by using it. The search speed can be improved to a certain extent by this means. The size of the array is 25 × 2 in this paper. The memory capacity is 25. One column is the arrangement on mistuned blades. Another column is the fitness value of this arrangement.To sum up, the specific process of optimization arrangement on mistuned blades based on tabu genetic cat swarm algorithm is shown in Figure5.Figure 5
Flowchart of blade arrangement optimized by tabu genetic cat swarm algorithm.


### 3.2. Improved Parallel Algorithm Based on CUDA

CUDA (Compute Unified Device Architecture) is a common parallel computing framework which is launched by NVDIA. It makes that GPU (graphic processing unit) can be as a parallel computing device. From this complex, combinatorial optimization problems can be solved. The memory space of CUDA based on GPU is shown in Figure6.Figure 6
Multistory memory space.The shared memory can be accessed quickly by all threads of one and the same block. The register is a cache of GPU that have fastest access speed. The local memory is private to a single thread. A complete CUDA program is composed of serial program of host and parallel computation kernel function of devices. Kernel function exists as a grid that is composed of several blocks. A block is a basic execution unit. It is composed of several threads that are executed in parallel. Moreover, different blocks cannot communicate with each other so that parallel computation is feasible. Details of the structure are shown in Figure7.Figure 7
The structure of threads.GPU has obvious advantages over CPU in processing capability and memory bandwidth. Floating-point processing capability of mainstream GPU is 10 times more than that of CPU during the same period. Moreover, external memory bandwidth of GPU is 5 times more than that of CPU. GPU allows fine-grained threads to communicate with each other. The detailed frameworks are shown in Figure8.Figure 8
Comparison of frameworks (a) CPU and (b) GPU.
(a)(b)Mass, stiffness, and damping matrices in dynamic equations are all high-dimensional floating-point numbers. Therefore, heterogeneous parallel mode of CPU + GPU can greatly improve the computational efficiency of the program. The complex logical processing and management run in CPU. The intensive parallel computing of large-scale data is executed by GPU. It has significant advantages in computational expense and performance.The exponential cycle computing of high-order matrices takes a lot of time in the optimization process. It cannot meet the needs in the specific application. GPU can be used for parallel optimization as each matrix operation is relatively independent.A lot of data should be concentrated together by matrices which represent complex models. Matrix multiplication is a highly applicable algorithm. The time complexity of the traditional algorithm isO(t3). Striped partition is used for matrix multiplication by GPU. Each thread reads one row and one column in two matrices, as shown in Figure 9. In this way, the work of each thread becomes the product of two one-dimensional vectors. Its time complexity is reduced to O(t2).Figure 9
Matrix multiplication by single thread.The operation is further optimized by matrix grid division. A product of a row and a column is calculated by the threads of one block which can use shared memory. Every block is used for calculating a submatrixCsub. And, every thread of the block is used for calculating an element of the submatrix. Csub is equal to the product of two rectangular matrices. Their dimensions are (wA, Block_size) and (Block_size, wA). They are divided into several square matrices whose dimensions are Block_size. The specific calculation is shown in Figure 10.Figure 10
Call graph of parallel computing.The two corresponding squares are loaded into shared memory from global memory. A thread is used to load one element of each matrix. Then, the product is calculated by each thread. The results of all these products are summarized in a register. When finished, the results are written to global memory. The faster shared memory is used and a large amount of global memory bandwidth is saved. The number of times that A or B is read from global memory are wA/Block_size. In this case, the time complexity isO(t).Thus, it can be seen that the time complexity and cycle indexes of matrix multiplication based on CUDA are reduced significantly. At the same time, the storage speed of shared memory is much faster than read and write speed of local and global memory. In terms of thread scheduling, the threads of GPU are very lightweight and switch within 1 cycle when CPU needs 100 cycles. The computational advantage of GPU is more obvious in case of multithread.The location update of cat swarm is not independent in Section3.1. This calculation method is not suitable for parallel processing. Therefore, it needs to be improved based on CUDA. The basic idea is to initialize multiple populations firstly. Each population updates the location of the cat swarm of tracking mode based on the last optimal solution. Then, the local and global optimal solutions are calculated simultaneously. Large random numbers are required as the planned numbers of cats and populations are numerous. It would have a large increase in the reading and writing of global memory and share memory by GPU that the needed random numbers are generated in time. A random number pool is designed in response to this situation. CuRand function is used to update every time. The improved algorithm is called multipopulation parallel tabu genetic cat swarm optimization algorithm. The detailed process is shown in Figure 11.Figure 11
Flowchart of improved tabu genetic cat swarm parallel algorithm based on CUDA.


### 3.3. Performance Analysis of the Parallel Algorithm

The advantages of the improved algorithm based on CUDA need to be evaluated by comparing the convergence speed and time complexity. The test environment is shown in Table2.Table 2
Test environment configuration.
EquipmentTypeCPUIntel Xeon E5 2600GPUNVIDA Tesla K20Internal storage64G RECCSystem type64-bit systemThe optimization problem of continuous function is tested first. Four classical test functions are used. Their expressions are as follows:(12)Sphere⟶fx=∑j=1nxj2,Rosenbrock⟶fx=∑j=1n100xj+1−xj22+1−xj2,Rastrigrin⟶fx=∑j=1nxj2−10cos2πxj+10,Griewank⟶fx=∑j=1nxj24000−∏j=1ncosxjj+1.The above test functions are inverted to observed optimal solution conveniently. They are shown in Figure12 when the independent variable is two dimension. Furthermore, the spent time is compared in Table 3.Figure 12
Diagram of test function. (a) Sphere. (b) Rosenbrock. (c) Rastrigrin. (d) Griewank.
(a)(b)(c)(d)Table 3
Comparison of computation time.
Test functionTime(s)Speed-up ratioCPUGPUSphere6.85210.119357.4359Rosenbrock7.12590.121258.7946Rastrigrin18.61680.288564.5303Griewank18.80230.294163.9306It can be seen from Table3 that the computation time is less than 0.3 second by using GPU. It has a tremendous increase in efficiency compared with the traditional serial mode of CPU. By this way, the efficiency can be increased dozens of times and the restriction of population size is broken. At the same time, the data structure is simple and the storage space is very small. The graphic memory of Nevida Tesla K20 is 6G. So, the problem of insufficient space can be ignored entirely. Spatial efficiency of the algorithm is exchanged for time efficiency. An effective method is provided for the combinatorial optimization problem that has large solution sets.Traveling salesman problem (TSP) is used for testing performance in the discrete solution space. It is aimed to find the minimum path cost of a single traveler who starts from the starting point and returns to the same point after all the given points are passed. The mathematical expressions are as follows:(13)C=c1,c2,…,cn,L=l11,…,lij,…,lnnG=C,L,,where C is a collection of n cities, L is the intercity distances, and G is a directed graph. The aim is to find the shortest Hamilton loop from G. Three kinds of TSP are selected in this paper. 40 cities, 80 cities, and 120 cities are used for the performance analysis. The results of optimal routes and optimization time are shown in Figure 13 and Table 4.Figure 13
Path planning of different kinds of TSP: (a) 40 cities, (b) 80 cities, and (c) 120 cities.
(a)(b)(c)Table 4
Comparison of computation time.
Number of citiesTime(s)Speed-up ratioCPUGPU4012.95131.41539.15098017.63281.604710.988212025.04312.092611.9675It can be seen from Table4 that using GPU has higher time efficiency. And, the speed-up ratio increases with the number of cities. It has an obvious advantage in solving large-scale high latitude problem. A low cost and high efficiency calculation method is provided for the case that the computation time is strictly limited.


### 3.4. Optimization Results’ Analysis

Mistuned I, Mistuned II, and Mistuned III are optimized based on the above parallel algorithm. The optimization arrangements are as follows:(14)Mistuned I⟶28,37,6,8,38,14,3,27,18,35,31,21,12,20,25,4,30,15,13,19,34,5,29,26,9,22,32,7,23,10,11,33,2,17,36,24,16,1,Mistuned II⟶20,19,15,4,1,10,27,29,9,25,28,16,2,32,21,30,34,36,8,14,33,12,13,31,24,6,26,22,7,35,38,18,5,11,17,3,23,37,Mistuned III⟶4,36,20,23,5,12,6,38,2,35,10,31,13,26,15,8,3,14,34,22,1,30,37,24,33,17,32,19,25,9,18,11,21,7,16,28,29,27.The optimal solution can be obtained by equation (4). It can be seen from Figure 14 and Table 5 that the vibration amplitudes of the bladed-disk system are obviously decreased with optimization arrangement. However, they are always larger than those of the tuned system. In addition, the width of the resonance region varies little. The maximum amplitudes are more concentrated than that before optimization and closer to the upper and lower vibration amplitude of the tuned system. The vibration localization of the mistuned system is reduced significantly by this way.Figure 14
Path planning of different kinds of TSP (a) Mistuned I (b) Mistuned II (c) Mistuned III.
(a)(b)(c)Table 5
Comparison of maximum amplitude and variance of different arrangements.
Mistuned IMistuned IIMistuned IIIMaximumVarianceMaximumVarianceMaximumVarianceTuned system202.55001005.6202.55001005.6202.55001005.6In order342.81252595.0384.26193514.7366.20724829.2Optimization291.88811352.0298.69073035.6287.10584382.8


## 3.1. Tabu Genetic Cat Swarm Optimization Algorithm

The cat swarm optimization algorithm is considered to be a method which can solve those complex optimization problems. It is combining search behavior and tracking behavior of a cat. It has an obvious advantage in computation speed, convergence, and sensitivity to initial conditions. It can be seen that the cat is a feasible solution for the optimization problem which we are going to solve. The behavior of cats can be categorized into two modes. One is called the search mode. It says that the cats are lazy or look around. The other one is called the tracking mode. It says that the cats are tracking a dynamic target. The two modes interact with each other by using mixture ratio (MR). MR indicates the proportion of cats in the whole group which executes the tracking mode.The search mode is used to simulate a cat’s current state. It means resting, looking around, and searching for the next moving location. In view of optimization arrangement on mistuned blades, it converts to express that the new arrangement is generated randomly. The quality of the solution is improved very slowly as the arrangement may be generated repeatedly in the random process. From this, tabu list is introduced to be the short-term memory storage of the algorithm to prevent circular search.The tracking mode is used to simulate the case that a cat is tracking a target. The location of the cat is updated by changing the speed of every single dimension. Furthermore, the speed is changed by adding random disturbance. The heuristic algorithm is needed to avoid the optimization problem getting into local optimal for the typical combinatorial optimization arrangement on mistuned blades. The solutions are improved constantly by introducing the mutation operator and crossover operator of genetic algorithm in this paper so that the global optimal solution can be found in a short time. The change of the solution is combined with the speed and location update process. The update formulas are defined as follows:(6)Vit+1=Vit+r⊗Xbestt⊕xit,xit+1=xit+Vit+1,where Xbest(t) is the location of the cat which has the best fitness value in the current cat swarm and Vi(t) and xi(t) are the speed and location of the ith cat at time t.The speed is to change the location of the cat. It is defined as an exchange list. The position change of the cat is to do exchange processing withx by the speed V. r is a random number between 0 and 1. Its meaning is probability. A random number will be generated for Vi(t) when Vi+1(t) is calculated. The genetic operators are called if the random number is greater than or equal to r, otherwise the speed stays the same. It can make the algorithm to keep good mutant ability by introducing genetic algorithm. The algorithm can avoid premature convergence effectively. The crossover operator and mutation operator play a key role in exciting cat swarm and keeping population diversity.In order to explain the process of crossover and mutation, two arrangements on ten mistuned blades are generated randomly as follows:(7)9,5,3,7,10,1,4,2,6,8,10,8,2,6,1,5,4,3,9,7.It means that No.9 is installed in position 1 and so on. Crossover operator (⊕) is based on the integer crossing method. Two crossover locations are selected. Then, the individual is crossed with the optimal one. It is assumed that 2 and 4 were randomly selected for crossover positions. The specific operation is as follows:(8)Individual9,5¯,3,7¯,10,1,4,2,6,8Optimization10,8¯,2,6¯,1,5,4,3,9,7⟶⊕9,8¯,3,6¯,10,1,4,2,6,8.The new individual has repeat positions obviously. The adjustment method is to replace with the missing blades:(9)9,8,3,6,10,1,4,2,6,8⟶9,8,3,6,10,1,4,2,5,7.Mutation operator (⊗) is based on the method of interchanging two random positions in an individual. Two mutation positions are selected and exchanged. It is assumed that 3 and 6 were selected. The specific operation is as follows:(10)9,8,3¯,6,10,1¯,4,2,5,7⟶9,8,1¯,6,10,3¯,4,2,5,7.The new individuals will be kept only when they are better than the old ones. So, an accurate evaluation criterion is needed. The vibration amplitudes and vibration differences are considered comprehensively, and the fitness function is designed as follows:(11)L=C1meanX+C2maxXC3varX+C4,where X is a vector of the maximum amplitudes of blades, mean(X) and var(X) are average value and variance of the maximum amplitudes of blades, and C1∼C4 are constants which are determined by the mistuned parameter.In addition, there is a great change to be repeated for blade arrangement as the randomness of new positions in the search mode. As a tabu list has short-term memory, the solution that has been searched can be excluded by using it. The search speed can be improved to a certain extent by this means. The size of the array is 25 × 2 in this paper. The memory capacity is 25. One column is the arrangement on mistuned blades. Another column is the fitness value of this arrangement.To sum up, the specific process of optimization arrangement on mistuned blades based on tabu genetic cat swarm algorithm is shown in Figure5.Figure 5
Flowchart of blade arrangement optimized by tabu genetic cat swarm algorithm.


## 3.2. Improved Parallel Algorithm Based on CUDA

CUDA (Compute Unified Device Architecture) is a common parallel computing framework which is launched by NVDIA. It makes that GPU (graphic processing unit) can be as a parallel computing device. From this complex, combinatorial optimization problems can be solved. The memory space of CUDA based on GPU is shown in Figure6.Figure 6
Multistory memory space.The shared memory can be accessed quickly by all threads of one and the same block. The register is a cache of GPU that have fastest access speed. The local memory is private to a single thread. A complete CUDA program is composed of serial program of host and parallel computation kernel function of devices. Kernel function exists as a grid that is composed of several blocks. A block is a basic execution unit. It is composed of several threads that are executed in parallel. Moreover, different blocks cannot communicate with each other so that parallel computation is feasible. Details of the structure are shown in Figure7.Figure 7
The structure of threads.GPU has obvious advantages over CPU in processing capability and memory bandwidth. Floating-point processing capability of mainstream GPU is 10 times more than that of CPU during the same period. Moreover, external memory bandwidth of GPU is 5 times more than that of CPU. GPU allows fine-grained threads to communicate with each other. The detailed frameworks are shown in Figure8.Figure 8
Comparison of frameworks (a) CPU and (b) GPU.
(a)(b)Mass, stiffness, and damping matrices in dynamic equations are all high-dimensional floating-point numbers. Therefore, heterogeneous parallel mode of CPU + GPU can greatly improve the computational efficiency of the program. The complex logical processing and management run in CPU. The intensive parallel computing of large-scale data is executed by GPU. It has significant advantages in computational expense and performance.The exponential cycle computing of high-order matrices takes a lot of time in the optimization process. It cannot meet the needs in the specific application. GPU can be used for parallel optimization as each matrix operation is relatively independent.A lot of data should be concentrated together by matrices which represent complex models. Matrix multiplication is a highly applicable algorithm. The time complexity of the traditional algorithm isO(t3). Striped partition is used for matrix multiplication by GPU. Each thread reads one row and one column in two matrices, as shown in Figure 9. In this way, the work of each thread becomes the product of two one-dimensional vectors. Its time complexity is reduced to O(t2).Figure 9
Matrix multiplication by single thread.The operation is further optimized by matrix grid division. A product of a row and a column is calculated by the threads of one block which can use shared memory. Every block is used for calculating a submatrixCsub. And, every thread of the block is used for calculating an element of the submatrix. Csub is equal to the product of two rectangular matrices. Their dimensions are (wA, Block_size) and (Block_size, wA). They are divided into several square matrices whose dimensions are Block_size. The specific calculation is shown in Figure 10.Figure 10
Call graph of parallel computing.The two corresponding squares are loaded into shared memory from global memory. A thread is used to load one element of each matrix. Then, the product is calculated by each thread. The results of all these products are summarized in a register. When finished, the results are written to global memory. The faster shared memory is used and a large amount of global memory bandwidth is saved. The number of times that A or B is read from global memory are wA/Block_size. In this case, the time complexity isO(t).Thus, it can be seen that the time complexity and cycle indexes of matrix multiplication based on CUDA are reduced significantly. At the same time, the storage speed of shared memory is much faster than read and write speed of local and global memory. In terms of thread scheduling, the threads of GPU are very lightweight and switch within 1 cycle when CPU needs 100 cycles. The computational advantage of GPU is more obvious in case of multithread.The location update of cat swarm is not independent in Section3.1. This calculation method is not suitable for parallel processing. Therefore, it needs to be improved based on CUDA. The basic idea is to initialize multiple populations firstly. Each population updates the location of the cat swarm of tracking mode based on the last optimal solution. Then, the local and global optimal solutions are calculated simultaneously. Large random numbers are required as the planned numbers of cats and populations are numerous. It would have a large increase in the reading and writing of global memory and share memory by GPU that the needed random numbers are generated in time. A random number pool is designed in response to this situation. CuRand function is used to update every time. The improved algorithm is called multipopulation parallel tabu genetic cat swarm optimization algorithm. The detailed process is shown in Figure 11.Figure 11
Flowchart of improved tabu genetic cat swarm parallel algorithm based on CUDA.


## 3.3. Performance Analysis of the Parallel Algorithm

The advantages of the improved algorithm based on CUDA need to be evaluated by comparing the convergence speed and time complexity. The test environment is shown in Table2.Table 2
Test environment configuration.
EquipmentTypeCPUIntel Xeon E5 2600GPUNVIDA Tesla K20Internal storage64G RECCSystem type64-bit systemThe optimization problem of continuous function is tested first. Four classical test functions are used. Their expressions are as follows:(12)Sphere⟶fx=∑j=1nxj2,Rosenbrock⟶fx=∑j=1n100xj+1−xj22+1−xj2,Rastrigrin⟶fx=∑j=1nxj2−10cos2πxj+10,Griewank⟶fx=∑j=1nxj24000−∏j=1ncosxjj+1.The above test functions are inverted to observed optimal solution conveniently. They are shown in Figure12 when the independent variable is two dimension. Furthermore, the spent time is compared in Table 3.Figure 12
Diagram of test function. (a) Sphere. (b) Rosenbrock. (c) Rastrigrin. (d) Griewank.
(a)(b)(c)(d)Table 3
Comparison of computation time.
Test functionTime(s)Speed-up ratioCPUGPUSphere6.85210.119357.4359Rosenbrock7.12590.121258.7946Rastrigrin18.61680.288564.5303Griewank18.80230.294163.9306It can be seen from Table3 that the computation time is less than 0.3 second by using GPU. It has a tremendous increase in efficiency compared with the traditional serial mode of CPU. By this way, the efficiency can be increased dozens of times and the restriction of population size is broken. At the same time, the data structure is simple and the storage space is very small. The graphic memory of Nevida Tesla K20 is 6G. So, the problem of insufficient space can be ignored entirely. Spatial efficiency of the algorithm is exchanged for time efficiency. An effective method is provided for the combinatorial optimization problem that has large solution sets.Traveling salesman problem (TSP) is used for testing performance in the discrete solution space. It is aimed to find the minimum path cost of a single traveler who starts from the starting point and returns to the same point after all the given points are passed. The mathematical expressions are as follows:(13)C=c1,c2,…,cn,L=l11,…,lij,…,lnnG=C,L,,where C is a collection of n cities, L is the intercity distances, and G is a directed graph. The aim is to find the shortest Hamilton loop from G. Three kinds of TSP are selected in this paper. 40 cities, 80 cities, and 120 cities are used for the performance analysis. The results of optimal routes and optimization time are shown in Figure 13 and Table 4.Figure 13
Path planning of different kinds of TSP: (a) 40 cities, (b) 80 cities, and (c) 120 cities.
(a)(b)(c)Table 4
Comparison of computation time.
Number of citiesTime(s)Speed-up ratioCPUGPU4012.95131.41539.15098017.63281.604710.988212025.04312.092611.9675It can be seen from Table4 that using GPU has higher time efficiency. And, the speed-up ratio increases with the number of cities. It has an obvious advantage in solving large-scale high latitude problem. A low cost and high efficiency calculation method is provided for the case that the computation time is strictly limited.


## 3.4. Optimization Results’ Analysis

Mistuned I, Mistuned II, and Mistuned III are optimized based on the above parallel algorithm. The optimization arrangements are as follows:(14)Mistuned I⟶28,37,6,8,38,14,3,27,18,35,31,21,12,20,25,4,30,15,13,19,34,5,29,26,9,22,32,7,23,10,11,33,2,17,36,24,16,1,Mistuned II⟶20,19,15,4,1,10,27,29,9,25,28,16,2,32,21,30,34,36,8,14,33,12,13,31,24,6,26,22,7,35,38,18,5,11,17,3,23,37,Mistuned III⟶4,36,20,23,5,12,6,38,2,35,10,31,13,26,15,8,3,14,34,22,1,30,37,24,33,17,32,19,25,9,18,11,21,7,16,28,29,27.The optimal solution can be obtained by equation (4). It can be seen from Figure 14 and Table 5 that the vibration amplitudes of the bladed-disk system are obviously decreased with optimization arrangement. However, they are always larger than those of the tuned system. In addition, the width of the resonance region varies little. The maximum amplitudes are more concentrated than that before optimization and closer to the upper and lower vibration amplitude of the tuned system. The vibration localization of the mistuned system is reduced significantly by this way.Figure 14
Path planning of different kinds of TSP (a) Mistuned I (b) Mistuned II (c) Mistuned III.
(a)(b)(c)Table 5
Comparison of maximum amplitude and variance of different arrangements.
Mistuned IMistuned IIMistuned IIIMaximumVarianceMaximumVarianceMaximumVarianceTuned system202.55001005.6202.55001005.6202.55001005.6In order342.81252595.0384.26193514.7366.20724829.2Optimization291.88811352.0298.69073035.6287.10584382.8


## 4. Optimization Verification Based on FEM

### 4.1. Mistuned Parameter Identification

The frequencies of blades are easy to be tested in engineering practice. However, they are hard to be used directly by the finite element model and the lumped parameter model. The frequency should be transformed into elastic modulus and stiffness for calculation, respectively.For the lumped parameter model, the mistuned parameterQi is introduced into the stiffness of the ith blade, namely,(15)kbi=kb1+Qi.The relation between frequency and stiffness is(16)2πfbi=kbimb,2πfb=kbmb,where fbiandkbi are the frequency and stiffness of the ith mistuned blade and fb, kb, and mb are the frequency, stiffness, and mass of the tuned blade. So, the mistuned parameter of the lumped parameter model can be derived as follows:(17)Qi=fbifb2−1.For the finite element model, the mistuned parameterPi is introduced into the elastic modulus of the ith blade, namely,(18)Ei=E1+Pi.An innovative method for mistuned parameter identification of the finite element model is proposed in this paper. First, the resonance method is used to test the first-order bending frequency of each blade. The test plan is shown in Figure15.Figure 15
Test plan for obtaining the natural frequencies of blades. (a) Schematic diagram. (b) Testing apparatus.
(a)(b)The finite element model of a single blade is used for modal analysis which is shown in Figure16(b). The linear displacement of the nodes at the contact position between blade and disk is fully constrained. The natural frequencies of first-order bending corresponding to the upper and lower limits of given elastic modules are calculated. The Newton iteration method is used to obtain the corresponding relation of test results and calculation results. The convergence criterion is considered that the error is less than 0.5%. A set of true test results of the aeroengine compressor is introduced in this paper. The tuned frequency is identified as 646 Hz. The results of mistuned parameter identification are shown in Figure 17.Figure 16
Finite element model of the aeroengine bladed disk.
(a)(b)Figure 17
Results of mistuned parameter identification.


### 4.2. Optimization Verification Analysis

The bladed-disk system of the actual aeroengine is used for optimization verification as shown in Figure16(a). There are 38 blades which consist of bodies and tenon. The complete system is formed by connecting the mortise and tenon. So, the nonlinear contact should be considered. At the same time, the system structure is very complex. The number of elements and nodes are 52163 × 38 and 79174 × 38, respectively. From this, an improved mixed-interface prestressed substructure modal synthesis method is proposed. The bladed-disk system is divided into 38 basic sectors. Each one is as a substructure which is shown in Figure 16(b).The dynamic analysis of the bladed-disk system based on FEM is shown in Figure18. The basic idea is as follows. The finite element model of a substructure should be established first. The nodes on common interface of each substructure and two adjacent substructures are selected to be main degrees of freedom which are constrained completely. The exciting force is applied to the nodes on the tip which are selected to be master DOF (degree of freedom) nodes. The prestress at working speed is analyzed and the superelement substructural model is generated. The newly generated file is called in the second place. The using part of the bladed-disk model is generated by coupling the nodes of the interface between adjacent substructures. From this, the condensed solution of master degrees of freedom of superelements could be obtained after the dynamic response analysis. Furthermore, it is expanded to all the elements and nodes of each substructure. The complete solutions are obtained in the end.Figure 18
Flowchart of the substructure modal synthesis method.The above mistuned system in Section4.1 is optimized based on the algorithm in Section 3.2. The final optimization arrangement is obtained as follows:(19)2,10,16,22,29,9,11,8,34,17,24,33,37,6,14,26,21,20,3,27,13,38,28,35,7,1,32,30,15,31,23,5,4,12,18,36,19,25.The dynamic analysis method is based on the above. It can be seen from Figure19 that the maximum amplitudes of order and optimization arrangement are 8.789 and 7.864. The mean variances are 1.1069e − 8 and 7.0690e − 9. The reductions are 10.52% and 36.14%, respectively. This shows that the parallel optimization algorithm presented in this paper is applicable to the mistuned bladed-disk system of aeroengine in engineering practice.Figure 19
Comparison of response amplitudes.


## 4.1. Mistuned Parameter Identification

The frequencies of blades are easy to be tested in engineering practice. However, they are hard to be used directly by the finite element model and the lumped parameter model. The frequency should be transformed into elastic modulus and stiffness for calculation, respectively.For the lumped parameter model, the mistuned parameterQi is introduced into the stiffness of the ith blade, namely,(15)kbi=kb1+Qi.The relation between frequency and stiffness is(16)2πfbi=kbimb,2πfb=kbmb,where fbiandkbi are the frequency and stiffness of the ith mistuned blade and fb, kb, and mb are the frequency, stiffness, and mass of the tuned blade. So, the mistuned parameter of the lumped parameter model can be derived as follows:(17)Qi=fbifb2−1.For the finite element model, the mistuned parameterPi is introduced into the elastic modulus of the ith blade, namely,(18)Ei=E1+Pi.An innovative method for mistuned parameter identification of the finite element model is proposed in this paper. First, the resonance method is used to test the first-order bending frequency of each blade. The test plan is shown in Figure15.Figure 15
Test plan for obtaining the natural frequencies of blades. (a) Schematic diagram. (b) Testing apparatus.
(a)(b)The finite element model of a single blade is used for modal analysis which is shown in Figure16(b). The linear displacement of the nodes at the contact position between blade and disk is fully constrained. The natural frequencies of first-order bending corresponding to the upper and lower limits of given elastic modules are calculated. The Newton iteration method is used to obtain the corresponding relation of test results and calculation results. The convergence criterion is considered that the error is less than 0.5%. A set of true test results of the aeroengine compressor is introduced in this paper. The tuned frequency is identified as 646 Hz. The results of mistuned parameter identification are shown in Figure 17.Figure 16
Finite element model of the aeroengine bladed disk.
(a)(b)Figure 17
Results of mistuned parameter identification.


## 4.2. Optimization Verification Analysis

The bladed-disk system of the actual aeroengine is used for optimization verification as shown in Figure16(a). There are 38 blades which consist of bodies and tenon. The complete system is formed by connecting the mortise and tenon. So, the nonlinear contact should be considered. At the same time, the system structure is very complex. The number of elements and nodes are 52163 × 38 and 79174 × 38, respectively. From this, an improved mixed-interface prestressed substructure modal synthesis method is proposed. The bladed-disk system is divided into 38 basic sectors. Each one is as a substructure which is shown in Figure 16(b).The dynamic analysis of the bladed-disk system based on FEM is shown in Figure18. The basic idea is as follows. The finite element model of a substructure should be established first. The nodes on common interface of each substructure and two adjacent substructures are selected to be main degrees of freedom which are constrained completely. The exciting force is applied to the nodes on the tip which are selected to be master DOF (degree of freedom) nodes. The prestress at working speed is analyzed and the superelement substructural model is generated. The newly generated file is called in the second place. The using part of the bladed-disk model is generated by coupling the nodes of the interface between adjacent substructures. From this, the condensed solution of master degrees of freedom of superelements could be obtained after the dynamic response analysis. Furthermore, it is expanded to all the elements and nodes of each substructure. The complete solutions are obtained in the end.Figure 18
Flowchart of the substructure modal synthesis method.The above mistuned system in Section4.1 is optimized based on the algorithm in Section 3.2. The final optimization arrangement is obtained as follows:(19)2,10,16,22,29,9,11,8,34,17,24,33,37,6,14,26,21,20,3,27,13,38,28,35,7,1,32,30,15,31,23,5,4,12,18,36,19,25.The dynamic analysis method is based on the above. It can be seen from Figure19 that the maximum amplitudes of order and optimization arrangement are 8.789 and 7.864. The mean variances are 1.1069e − 8 and 7.0690e − 9. The reductions are 10.52% and 36.14%, respectively. This shows that the parallel optimization algorithm presented in this paper is applicable to the mistuned bladed-disk system of aeroengine in engineering practice.Figure 19
Comparison of response amplitudes.


## 5. Conclusion

The lumped parameter model of the bladed-disk system is established which has two degrees of freedom in this paper. It can be found that the amplitude-frequency characteristic curves of the tuned system have the two same resonance peaks. However, the amplitude-frequency characteristic curves of the mistuned system are very different. There are many small peaks besides two main resonance peaks. At the same time, the resonance peaks and regions increase with the standard deviations of mistuning. Furthermore, there is a minimum amplitude between the resonance regions. The amplitudes and the sensitivities of mistuning can be reduced effectively when the frequency is near 1.026.Tabu genetic cat swarm optimization algorithm is proposed for optimization arrangement on mistuned blades. Furthermore, it is improved to be a parallel algorithm based on CUDA. Its performance is analyzed in continuous and discrete solution space, respectively. In conclusion, the parallel algorithm is low cost and high efficiency. In result, it can be found that the vibration amplitudes of the bladed-disk system are obviously decreased with optimization arrangement. The vibration localization of the mistuned system is also reduced significantly.The complete mistuned parameter identification method and the improved mixed-interface prestressed substructure modal synthesis method are proposed for optimization verification. It can be found that the maximum amplitudes and mean variance are significantly decreased. Therefore, it can be concluded that the parallel optimization algorithm presented in this paper is applicable to the mistuned bladed-disk system of the aeroengine in engineering practice.



---

*Source: 6653692-2021-06-02.xml*