Riding comfort optimization of railway trains based on pseudo-excitation method and symplectic method

Journal of Sound and Vibration 332 (2013) 5255–5270

Contents lists available at SciVerse ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

Riding comfort optimization of railway trains based on pseudo-excitation method and symplectic method You-Wei Zhang a, Yan Zhao a,n, Ya-Hui Zhang a, Jia-Hao Lin a, Xing-Wen He b a State Key Laboratory of Structural Analysis for Industrial Equipment, Faculty of Vehicle Engineering and Mechanics, Dalian University of Technology, Dalian 116023, PR China b Division of Engineering and Policy for Sustainable Environment, Faculty of Engineering, Hokkaido University, Sapporo 060-8628, Japan

a r t i c l e i n f o

abstract

Article history: Received 5 February 2012 Received in revised form 6 May 2013 Accepted 7 May 2013 Handling Editor: H. Ouyang Available online 19 June 2013

This research is intended to develop a FEM-based riding comfort optimization approach to the railway trains considering the coupling effect of vehicle–track system. To obtain its accurate dynamic response, the car body is modeled with finite elements, while the bogie frames and wheel-sets are idealized as rigid bodies. The differential equations of motion of the dynamic vehicle–track system are derived considering wheel–track interaction, in which the pseudo-excitation method and the symplectic mathematical method are effectively applied to simplify the calculation. Then, the min–max optimization approach is utilized to improve the train riding comfort with related parameters of the suspension structure adopted as design variables, in which 54 design points on the car floor are chosen as estimation locations. The K–S function is applied to fit the objective function to make it smooth, differentiable and have superior integrity. Analytical sensitivities of the K–S function are then derived to solve the optimization problem. Finally, the effectiveness of the proposed approach is demonstrated through numerical examples and some useful discussions are made. & 2013 Elsevier Ltd. All rights reserved.

1. Introduction Track irregularities are the most important sources of random excitation of coupled vehicle–track systems. Therefore, riding comfort optimization of the running train subjected to track irregularities is a significant issue for railway train designs. The suspension system is the key mechanical structure affecting the riding comfort of trains. Enormous efforts have been devoted by many researchers to optimize the suspension system to improve the train riding comfort, taking advantages of nearly all the available optimization methods. Generally, there are two kinds of approaches mainly used in the optimization problems: gradient-based approximation techniques and intelligent optimization methods. Researches on the gradient-based approaches started much earlier and enormous work has been done [1–8]. The drawback of this technique is that it may be trapped in local minimum points, thus the optimized result greatly depends on the initial starting points. Besides, the calculation of gradient information is complicated and costly, especially when many design variables are considered. The intelligent optimization methods were introduced into this field merely within the last decades, including genetic algorithms [9–11], neural networks [12,13], simulated annealing [14], AGOP methods [15], and so on. In order to achieve feasible global solutions, these methods generally require a large starting population. In general, no matter which method is used, iterative dynamic response analyses of the structures are inevitable. If considerably large

n

Corresponding author. Tel.: +86 411 84706337; fax: +86 411 84708400. E-mail address: [email protected] (Y. Zhao).

0022-460X/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jsv.2013.05.018

5256

Y.-W. Zhang et al. / Journal of Sound and Vibration 332 (2013) 5255–5270

computational capacity is required at each step of iteration, the optimization process will be extremely time-consuming. As a result, only very simple vehicle models have been employed for the optimization problems in the previous researches, and very few studies have focused on the optimization of complicated vehicle–track interaction system subjected to random track irregularities owing to the extremely high computational demands. The major difficulties of the random vibration analyses of FEM-based vehicle–track systems are the huge degrees of freedom (dof) of the coupled system and the low efficiency of the conventional random vibration analysis methods. When using the finite element method (FEM) to establish the equation of motion of the track structure, a sufficient length of track section, including dozens of sleeper spacings before and after the vehicle, must be considered in the numerical model to ensure reasonable boundary conditions. Therefore, even if a low-dof spring-mass-damper vehicle model is used, the dofs of the coupled system are still considerably large. For such a coupled system, the conventional random vibration methods are usually very time consuming, because of calculating the frequency response function matrix and performing large matrix multiplications. If the Monte-Carlo simulation is used instead, enormous track irregularity samples, very short integral time steps, and a relatively long time history for each sample to reach the steady-state are required [16–18]. On the other hand, a real car body is actually an elastic structure and has very complicated dynamic characteristics. If the spring-mass-damper vehicle model is used considering the car body as rigid mass, it can lead to numerical errors and the random vibrations at different parts of the car body cannot be properly expressed. In order to properly evaluate the riding comfort of all passengers at different locations, an elaborate finite element (FE) model of the car body is necessary. However, the computational capacity needed by such complicated models is even more unacceptable for the optimization problems. In order to overcome the above difficulties, by firstly considering the track as an infinite length of substructural chains, Lu et al. [19] developed an efficient method to analyze the vertical random responses of the coupled system using a fixedvehicle model, in which the pseudo-excitation method (PEM) [20,21] and the symplectic mathematical method [22,23] are effectively applied. Xu et al. [24] derived the sensitivity equations of the coupled vehicle–road and vehicle–bridge system for the optimization of vehicle suspension systems to improve the riding comfort using PEM. Based on these work, a FEM-based riding comfort optimization approach of the head car of a high-speed train using dynamic analyses of the coupled vehicle–track system subjected to random track irregularities is developed in this paper. Both PEM and symplectic mathematical methods are used to calculate the pseudo responses and analytical sensitivities, which can achieve an improvement of the computational efficiency by several orders of magnitude compared with the past methods. Therefore, the optimization problems with an elaborate finite element vehicle model can be solved on a standard personal computer. The commercial software ANSYS is used to establish the FE model of the car body and calculated its free mode, and so the equation of motion of the car body is derived using mode superposition method. 54 design points in 3 rows corresponding to the passenger seat positions at the left-side floor of the vehicle are chosen for riding comfort estimation and optimization. The bogie frames and wheel-sets are treated as rigid bodies, connected by primary- and secondary-suspension structures. The track is regarded as an infinitely long substructural chain consisting of three layers, i.e. the rails, sleepers and ballast. The rails are modeled as space beam finite elements and the sleepers and ballast are regarded as rigid bodies. They are connected by springs and dampers between them. Considering that the major research target in this paper is the car body, the low time-consuming fixed-vehicle model is used, which is sufficient to ensure the numerical accuracy [19]. The detailed execution steps of the developed optimization approach and its academic features are described as follows: (1) The PEM [20,21] is applied to convert complicated random vibration analyses into simple pseudo response analyses by transforming three types of rail irregularities, i.e. the longitudinal level irregularity, alignment irregularity and the cross level irregularity, into deterministic generalized single-point harmonic excitations. Then, the response power spectral densities (PSDs) of the coupled system can be obtained conveniently by means of harmonic analysis. This method will make the calculation process of the FEM-based random vibration analysis and riding comfort optimization of the coupled system highly simple and efficient. (2) The symplectic mathematical method [22,23] is then used to establish the low-dof equations of motion of the 3-dementional track structure. Based on this method, only the track elements currently in contact with the wheels need to be included in the computation, which decreases the computational dofs of the track model by 1–2 order of magnitude compared with conventional FE track model. Since the periodic continuous conditions are used, the accurate theoretical solution of the infinitely long track is obtained. (3) Based on the above two improvements, the equations of motion of the coupled system using the fixed-vehicle model subjected to pseudo-excitations is derived and the random vibration analyses are performed. The acceleration PSDs of the design points can then be obtained accordingly, using which the weighted root mean square (rms) values of the acceleration response defined in the international standard ISO2631 [25,26] are calculated to evaluate the riding comfort of these points. (4) In order to improve the riding comfort of all passengers at different locations in the car, a min–max optimization problem are proposed, defining total 8 parameters of the stiffness and damping coefficients of the primary- and secondary-suspension structures in both vertical and lateral directions as design variables. The K–S function [27,28] is then applied to fit the curve of the original objective function. Because the K–S function is smooth, differentiable and has superior integrity, it is possible to derive the analytical sensitivity of the new objective function. (5) By differentiating the equations of motion of the coupled system, the first- and second-order stationary sensitivity equations can be derived, which are subjected to harmonic pseudo-excitations and have low dofs due to the application

Y.-W. Zhang et al. / Journal of Sound and Vibration 332 (2013) 5255–5270

5257

of PEM and symplectic mathematical method. Then, the accurate analytical sensitivities of the pseudo-responses can be calculated directly by solving these equations. The corresponding analytical sensitivities of the response PSDs, the weighted rms values of acceleration responses and the K–S function can be obtained accordingly by some simple mathematical operations. Since there are no other approximations used in the calculation, the proposed method for sensitivity analysis is highly accurate and retains the efficient features of PEM and symplectic mathematical method. (6) Finally, by inputting the above objective function and its first- and second-order analytical sensitivities into the ‘fmincon’ function in the optimization toolbox of MATLAB [29], which is based on the sequential quadratic programming (SQP), the optimized design variables can be obtained and the riding comfort optimization of the train is realized.

The effectiveness of the proposed optimization approach is then demonstrated through several numerical examples. Also, the optimization results using the developed FE model are compared with those by using a rigid body vehicle model. At last, the influences of vehicle velocities and class of track irregularity on optimization results are also discussed.

Fig. 1. Analytical model of coupled vehicle–track system. (a) Front view, (b) top view, (c) left view.

5258

Y.-W. Zhang et al. / Journal of Sound and Vibration 332 (2013) 5255–5270

2. Analytical models of vehicle–track interaction system The analytical models of the coupled vehicle–track system are shown in Figs. 1–3. The vehicle is composed of a body, two bogie frames, four wheel-sets, which are connected by the primary- and secondary-suspension structures between them. Fig. 1 shows the numerical model of the vehicle, in which the body is expressed as a rigid mass. In order to accurately calculate its detailed responses of different parts, an elaborate FE model of the car body as shown in Fig. 2 is developed by ANSYS in this research. Considering that the car body is symmetric, 54 design points in 3 rows corresponding to the passenger seats positions at the left-side floor of the car are chosen for riding comfort estimation and optimization. The bogie frames and wheel-sets are treated as rigid bodies formulized by rigid dynamics. Each bogie frame has 5 dofs, i.e. vertical, pitching, lateral, yawing, and roll displacements, while each wheel has 3 dofs, i.e. vertical, lateral and roll displacements, respectively. The track is regarded as an infinitely long substructural chain consisting of three layers, i.e. the rails, sleepers and ballast, as shown in Fig. 3. Each substructure of the track consists of a pair of rail segments between two adjacent sleepers, and the corresponding ballast segment. The rail segments are modeled as space beam elements with 5 dofs on each node ignoring the axial dof. The sleepers and ballast are both regarded as rigid bodies. A sleeper has 3 dofs, i.e. the vertical, lateral, and roll displacements, while a ballast element has only the vertical displacement. Assuming that the vehicle is moving with a uniform velocity and the wheels are closely attached to the rails without relative sliding, then the relationship between the displacements of wheel-sets and rails is as follows: 9 2 8 0 z > = 6 < wj > 1 uwj ¼ ywj ¼ 6 42 > ; :φ > 0 wj

1 2

0

0

1 2

1 2b

0

8 n 9 3> 9 > 8 > yrjL > > > > > > r vj > > n > = 7< zrjL = < r 0 7 5> yn > þ > aj > > 1 1 > : 2b > rjR > r cj ; − 2b > > > ; : znrjR > 1 2

(1)

where zwj, ywj and φwj are the vertical, lateral and roll displacements of the jth wheel-set; ynrjL , znrjL , ynrjR and znrjR are the vertical and lateral displacements of the track at the wheel-rail contact points on the left and right sides, which can be obtained from the shape function matrix and the displacements of the rails; rvj, raj and rcj are the longitudinal level, alignment and cross level track irregularities at the wheel-rail contact points, respectively. The parameter b is half of the distance between the two rails.

Fig. 2. FE model of car body.

Fig. 3. Analytical model of track structures.

Y.-W. Zhang et al. / Journal of Sound and Vibration 332 (2013) 5255–5270

5259

3. PEM applied to random vibration analysis of vehicle–track system The PEM [20,21] is an accurate and highly efficient algorithm for stationary/non-stationary random structural response analysis particularly for large scale structures with a huge number of dofs. This method can either convert a stationary random vibration analysis into a deterministic harmonic analysis or convert a non-stationary random vibration analysis into an ordinary time-history analysis. It is not only simple and efficient, but also theoretically accurate. The random responses of structures with millions of dofs and hundreds of excitation points can be calculated in a remarkably short time on a standard personal computer by using PEM. Considering that the major research object in the present paper is the car body, the fixedvehicle model is used [19], i.e. the vehicle remains stationary on the track while the track irregularities are moving backwards with the vehicle velocity. Hence, the pseudo-excitation of the coupled system with a fixed-vehicle model will be harmonic. In this section, the application of PEM to the coupled vehicle–track system is outlined. The three types of track irregularities are assumed as mutually incoherent zero-mean value stationary random processes with PSDs of Sv(ω), Sa(ω) and Sc(ω). Considering the four wheel-sets of the vehicle are excited by the same track, the four excitations are fully-coherent and the vectors of track irregularities at the wheel-rail contact points are given as follows: rj ðtÞ ¼ fr j ðt−t 1 Þ r j ðt−t 2 Þ r j ðt−t 3 Þ r j ðt−t 4 ÞgT ;

ðj ¼ v; a; cÞ

(2)

where t1 ¼0, t2 ¼2l1/v, t3 ¼2l2/v, t4 ¼2(l1+l2)/v; 2l1 and 2l2 are the distances between the wheel pair and between the bogie frames, respectively. The input PSD matrix of rj(t) is then expressed as below. 2

1 6 eiωðt2 −t1 Þ 6 j Srr ðωÞ ¼ 6 iωðt −t Þ 4e 3 1

eiωðt1 −t2 Þ

eiωðt4 −t1 Þ

eiωðt4 −t2 Þ

1

eiωðt1 −t 3 Þ

eiωðt1 −t 4 Þ

iωðt 2 −t 3 Þ

iωðt 2 −t 4 Þ

e

iωðt 3 −t 2 Þ

e

1

3

7 e 7 7S ðωÞ iωðt 3 −t 4 Þ 5 j e

eiωðt4 −t 3 Þ

(3)

1

Defining the response of the coupled system as u(t), the PSD of its response subjected to single type track irregularity can be given by the following equation, based on the conventional random vibration theory: Sjuu ðωÞ ¼ Hn ðωÞSjrr ðωÞHT ðωÞ

(4)

n

where H(ω) is the frequency response function matrix and superscripts and T represent the complex conjugate and the transpose, respectively. Then the response PSD of the coupled system subjected to total three types of track irregularities can be expressed as the superposition of those of single type ones. Suu ðωÞ ¼

∑ Sjuu ðωÞ

(5)

j ¼ v;a;c

As indicated by Eqs. (4) and (5), it will be extremely time consuming to directly compute the frequency response function matrix and perform matrix multiplications, particularly for large scale structures. Therefore in this research, the PEM is effectively applied, assuming that the structure is subjected to the following pseudo-excitations. qffiffiffiffiffiffiffiffiffiffiffi r~ j ðω; tÞ ¼ fe−iωt1 e−iωt 2 e−iωt3 e−iωt 4 gT Sj ðωÞeiωt ; ðj ¼ v; a; cÞ (6) It is obvious that these are generalized single-point harmonic excitations and Eqs. (4) and (5) can be rewritten as follows: Suu ðωÞ ¼

T n ∑ u~ j ðω; tÞu~ j ðω; tÞ

(7a)

j ¼ v;a;c

u~ j ðω; tÞ ¼ HðωÞr~ j ðω; tÞ

(7b)

It is clear that u~ j ðω; tÞ is the stationary pseudo-response of the structure when it is subjected to the pseudo-excitation of r~ j ðω; tÞ. Then, the response acceleration PSD can be obtained as follows: Su€ u€ ðωÞ ¼ ω4 Suu ðωÞ

(8)

It is clear that Eq. (7) is much easier to compute than Eqs. (4) and (5). Therefore, PEM leads to a considerable reduction in computational effort by transforming the three types of rail irregularities into corresponding deterministic generalized single-point harmonic excitations, which makes the FEM-based random vibration analysis and riding comfort optimization of the large scale coupled system possible.

4. Equations of motion of coupled vehicle–track system subjected to pseudo-excitations In this section, the equations of motion of the coupled vehicle–track system subjected to pseudo-excitations are established using mode superposition method and symplectic mathematical method.

5260

Y.-W. Zhang et al. / Journal of Sound and Vibration 332 (2013) 5255–5270

4.1. Equations of motion of the vehicle The equation of motion of the car body is expressed as follows: € c þ Cc u _ c þ Kc uc ¼ CCct u _ t þ CKct ut Mc u

(9)

where uc and ut are the displacement vectors of the car body and the bogie frames; Mc, Cc and Kc are the mass, damping and stiffness matrices; and CCct and CKct are the load coefficient matrices, respectively. By applying the mode superposition method, the mode equation of motion of the car body can be obtained as below. ct

ct

_ t þ Ka ut Mc q€ c þ Cc q_ c þ Kc qc ¼ Ca u

(10)

where qc is the mode displacement vector of the car body; Mc , Cc and Kc are respectively the mode mass, damping and stiffness matrices; and the corresponding load coefficient matrices are obtained as below, where Ψc is the mode matrix. ct

ct

Ca ¼ ΨTc CCct ; Ka ¼ ΨTc CKct

(11)

It is noticed that in the computation, Mc , Kc and Ψc are obtained directly by ANSYS and Cc is synthesized by proportional damping. The bogie frames are assumed to be rigid bodies, which can be modeled by basic dynamic theory and the equation of motion is as below. tc

tc

€ t þ Ct u _ t þ Kt ut ¼ Ca q_ c þ Ka qc þ Q Ctw u _ w þ Q Ktw uw Mt u

(12)

where uw ¼ fuTw1 uTw2 uTw3 uTw4 gT is the displacement vector of the wheels; Mt, Ct and Kt are the mass, damping and stiffness tc ct tc ct matrices; Ca and Ka are the transpose of matrix Ca and Ka in Eq. (11) to couple the car body and the bogie frames; and Q Ctw K and Q tw are the load coefficient matrices to generate the loads acting on the bogie frames subjected to the wheel displacements. Considering that the displacements of the wheels are non-independent and can be expressed by the displacements of the rails as shown in Eq. (1), the total equations of motion of the vehicle including the body and the bogie frames can be integrated as follows: 3 3 " #( ) 2 ) 2 ) " # " # ct ( ct ( 0 0 Cc −Ca Kc −Ka qc q_ c q€ c Mc 5 5 _ þ (13) þ 4 tc þ4 ¼ u u C K w tc _t €t Q tw w Q tw ut u u Mt −Ka Kt −Ca Ct The above equations are then redefined as below. € v þ Cv u _ v þ Kv uv ¼ Q Cvt u _ w þ Q Kvt uw Mv u

(14)

The response of the vehicle will be stationary due to the fixed-vehicle model subjected to harmonic excitations. In this case, Eq. (14) can be rewritten as follows: ðKv þ iωCv −ω2 Mv Þuv ¼ Kdv uv ¼ ðiωQ Cvt þ Q Kvt Þuw

(15)

4.2. Equations of motion of the track In this paper, the track is regarded as an infinitely long substructural chain as shown in Fig. 4. In the substructural chain, the unloaded identical substructures are denoted as sub, while the loaded ones subjected to the load fe(t) are indicated as substructure subn, respectively. The substructure of the track is shown in Fig. 3. The equation of motion of the loaded substructure is as Eq. (16). € e þ Ce u _ e þ Ke ue ¼ f e ðtÞ þ f b Me u

Fig. 4. Infinitely long track substructural chain.

(16)

Y.-W. Zhang et al. / Journal of Sound and Vibration 332 (2013) 5255–5270

5261

where ue is the displacement vector of the substructure; Me, Ce and Ke are the mass, damping and stiffness matrices; and fb is the external nodal force vector, respectively. When the substructure is subjected to a harmonic load whose circular frequency is ω, Eq. (16) can be rewritten as below. ð−ω2 Me þ iωCe þ Ke Þue ¼ Kde ue ¼ f e ðtÞ þ f b Its block form can be expressed as follows: 2 0 Kaa K0ab 6 0 6 Kba K0bb 4 K0ia K0ib

38 9 8 9 8 9 K0ai > ua > > f ae > > pa > < = < = < = 7 K0bi 7 5> ub > ¼ > f be > þ > −pb > :u ; :f ; : ; 0 i ie K0ii

(17)

(18)

where ua and ub are the external displacement vectors and ui is the internal displacement vector; fae, fbe and fie are the external load vectors; and Pa and Pb are the corresponding external nodal force vectors, respectively. Pa and Pb can be calculated by the external stiffness matrices Pα and Pβ as follows: pa ¼ −Pβ ua ; pb ¼ pα ub

(19)

For an unloaded identical substructure, the following relation is proved [23]: ( ) " #( ) ( ) ua ub ua Saa Sab ¼ ¼ SðωÞ pb p pa Sba Sbb a

(20)

where, Saa ¼ −ðKab Þ−1 Kaa ; Sab ¼ ðKab Þ−1 Sba ¼ −Kba þ Kbb ðKab Þ−1 Kaa ; Sbb ¼ K0bb −K0bi ðK0ii Þ−1 K0ib Kaa ¼ K0aa −K0ai ðK0ii Þ−1 K0ia ; Kab ¼ K0ab −K0ai ðK0ii Þ−1 K0ib Kba ¼ ðKab ÞT ; Kbb ¼ K0bb −K0bi ðK0ii Þ−1 K0ib herein, S(ω) is the symplectic transfer matrix [22,23] that satisfies the below equation. " # 0 In ; JTn ¼ J−1 ST Jn S ¼ Jn ; Jn ¼ n ¼ −Jn −In 0

(21)

(22)

where In is the n-dimensional unit matrix; and if μ is an eigenvalue of S(ω), then so is 1/μ. These eigenvalues are called wave propagation constants. Hence, the corresponding eigenvectors can be used to constitute the following matrix: " # Xa Xb (23) Φ ¼ ½φ1 φ2 ⋯ φ2n  ¼ Na Nb here, φj(j¼1,2,…n) is the eigenvector corresponding to μj(|μj|≤1). Then the external stiffness matrices can be obtained by −1 Pα ¼ Na X−1 a ; Pβ ¼ −Nb Xb

(24)

For the loaded substructure, substituting Eq. (19) into Eq. (18) gives 2 0 38 9 8 9 K0ab Kaa þ Pβ K0ai > ua > > f ae > = 6 7< = < 0 0 0 6 Kia 7 ui ¼ f be Kii Kib 4 5> > > :u ; : f > ; b ie K0ba K0bi K0bb þ Pα Then any required response of the kth substructure can be obtained by the following equation: 82 3 > Xb μk X−1 > b > > 6 7 > > 6 7ua k o0 > Xb μk1 X−1 > b 4 5 > > > 0 −1 0 0 −1 0 k −1 k1 −1 > < −ðKii Þ Kia Xb μ Xb −ðKii Þ Kib Xb μ Xb 3 uek ¼ 2 > Xa μk1 X−1 > a > >6 7 > > 6 7u b k 4 0 > Xa μk X−1 > a 4 5 > > > 0 0 −1 0 0 −1 −1 −1 k1 k > −ðK Þ K Xa μ X −ðK Þ K Xa μ X : ii

ia

a

ii

ib

(25)

(26)

a

where, k 40 corresponds to the substructures on the front side of the loaded substructure and k o0 to the substructures on the rear side. Here μ is a diagonal matrix of order n with its jth diagonal element equal to μj(|μj|≤1).

5262

Y.-W. Zhang et al. / Journal of Sound and Vibration 332 (2013) 5255–5270

The four wheel-sets of the studied vehicle lead to four corresponding loaded track substructure. The low-dofs equation of motion of the track can be obtained as Eq. (27), by coupling the equations of motion of these substructures based on Eq. (26) with their influence to each other taken into account. C _ K C _ K € Kdtr utr ¼ Q M tt uw þ Q tt uw þ Q tt uw þ Q tv uv þ Q tv uv

(27) QM tt ,

Kdtr

Q Ctt ,

Q Ktt ,

Q Ctv

Q Ktv

where, utr is the track displacement vector; is the dynamic stiffness matrix; and and are the corresponding load coefficient matrices, respectively. If each rail in the track substructure is divided into 3 FE elements, the computational dofs of Eq. (27) are 96 using symplectic method and about 3500 using the conventional finite long track model. It can be seen that the symplectic method not only considerably reduces the dofs of the track, but also avoids the error causes by cutting the infinitely long track as limited track sections. 4.3. Equations of motion of coupled vehicle–track system It is known from the above two sections that the wheel displacement vector uw is in the load vectors of both equations of motion of the vehicle and the track. Thus it can be expressed using the rail-wheel relationship in Eq. (1) as below. uw ¼ Twt utr þ

∑ Γj r j

(28)

j ¼ v;a;c

where, Twt is the transform matrix and Γj(j¼v,a,c) are the coefficient matrices corresponding to the track irregularities. Substituting Eqs. (6) and (28) into Eqs. (15) and (27), the equation of motion of the coupled vehicle–track system subjected to pseudo-excitations can be obtained as follows: 2 3( ) 2 j 3 u~ vj Rv qffiffiffiffi Kdv −iωRCvt −RKvt 4 5 ¼ 4 j 5 Sj eiωt ; ðj ¼ v; a; cÞ (29) C K d ~ u trj −iωR tv −R tv Ktr Rtr The above equations are then rewritten as below. Kd u~ j ¼ Rj

qffiffiffiffi Sj eiωt

(30)

where u~ vj and u~ trj are the pseudo-responses of the vehicle and the track when the coupled system is subjected to the jth pseudo track irregularity; R Cvt , RKvt , RCtv and RKtv are the coupling matrices; and Rjv and Rjtr are the load coefficient matrices of the jth track irregularity, respectively. It is obvious that Eq. (30) can be solved as linear equations, and then the displacement and acceleration PSDs of the coupled system can be obtained by Eqs. (7) and (8). 5. Riding comfort analysis using international standard ISO2631 In this section, the international standard ISO2631 is used to evaluate the riding comfort of the design points. This international standard is formulated by International Standardization Organization and has been in effect since 1974. A riding comfort index of weighted rms acceleration for mechanical vibration and shock evaluation of human exposure to whole body vibration is proposed in this standard, which is usually taken as the objective function in optimization problems. According to this international standard ISO2631, the weighted rms acceleration of design point n at the directions of x, y and z in Fig. 2 are defined as below. !1=2 Z anq w ¼

80

0:5

w2q ðf ÞSnq ðf Þdf u€ u€

;

ðn ¼ 1; 2; :::; 54; q ¼ x; y; zÞ

(31)

where wq is the weighted factor. Considering that the passenger seats are not included in the model of car body, the following parameters for standing persons are adopted: ( 1 ð0:5 Hz of o 2 HzÞ wx;y ðf Þ ¼ 2 ð2 Hz of o 80 HzÞ (32a) f

wz ðf Þ ¼

8 0:5 > > > f > < 4

1 > > > > : 12:5 f

ð0:5 Hz o f o2 HzÞ ð2 Hz o f o4 HzÞ ð4 Hz o f o12:5 HzÞ

(32b)

ð12:5 Hz o f o80 HzÞ

where, Snq ðf Þ is the frequency acceleration PSD of point n at the direction of q with the following relationship: u€ u€ ðf Þ ¼ 2πSnq ðωÞ Snq u€ u€ u€ u€

(33)

Y.-W. Zhang et al. / Journal of Sound and Vibration 332 (2013) 5255–5270

5263

herein, Snq ðf Þ can be obtained by Eq. (8). The total value of weighted rms acceleration is as follows: u€ u€ !1=2 awn ¼

2

2 ∑ kq ðanq w Þ

(34)

q ¼ x;y;z

here in above kx ¼ky ¼kz ¼1 is adopted in this paper. 6. Riding comfort optimization of coupled vehicle–track system The vehicle suspension systems, as depicted in Fig. 1, are the most important mechanical structures affecting the riding comfort of passengers. In this paper, the vertical and lateral stiffness and damping coefficients of the vehicle primary- and secondary-suspensions, i.e. total 8 parameters shown in Eq. (35), are considered as design variables for riding comfort optimization. α ¼ fC z2 C y2 C z1 C y1 K z2 K y2 K z1 K y1 gT

(35)

In order to reduce the weighted rms acceleration responses of all the 54 design points, the following max–min optimization problem is proposed. 8 α > < Find Min maxðawn Þ ðn ¼ 1; 2; …; 54Þ (36) > : S:t: α ≤α ≤α l

u

in which awn is the weighted rms acceleration of point n, and αl and αu are the lower and upper limit of the design variable, respectively. 6.1. Application of K–S function The K–S function is firstly proposed by Kreisselmeier and Steinhauser [27] in 1979. It has been widely used in the optimization field nowadays, because of its remarkably global, smooth and differentiable features. It is usually employed to fit the curves of the multiple functions. Since the objective function of the defined min–max optimization problem is not differentiable and the differences of the values of weighted rms acceleration responses at different design points are small, the K–S function is applied to fit the objective function in Eq. (36) as below.  54  1 (37) KSðαÞ ¼ ln ∑ eρawn ðαÞ ρ n¼1 It is proven that the following property of K–S function exists if ρ-∞ [28]: maxðawn Þ ≤KS ≤maxðawn Þ þ Then Eq. (36) can be rewritten as follows:

8 > < Find Min > : S:t:

1 lnð54Þ ρ

(38)

α KSðαÞ

(39)

αl ≤α ≤αu

The new objective function in Eq. (39) is differentiable, thus it is possible to derive its analytical sensitivity. 6.2. Sensitivity equations of vehicle–track system based on PEM The first- and second-order sensitivity equations of the coupled vehicle–track system subjected to pseudo-excitations are derived in this section. Denoting the kth design variable of vector α in Eq. (35) as αk, the first-order sensitivity equation can be obtained as below by differentiating Eq. (30) with respect to αk. ∂u~ j ∂R j qffiffiffiffi iωt ∂Kd u~ Kd Sj e − ¼ (40) ∂αk ∂αk ∂αk j Then, the second-order sensitivity equation is derived by differentiating the first-order sensitivity equation with respect to design variable α1 as follows:  2  ∂2 u~ j ∂2 Rj qffiffiffiffi iωt ∂ Kd ∂K ∂u~ j ∂Kd ∂u~ j u~ j þ d Sj e − ¼ þ (41) Kd ∂αk ∂αl ∂αk ∂αl ∂αk ∂αl ∂αk ∂αl ∂αl ∂αk It is obvious that the load vectors of Eqs. (40) and (41) are harmonic because of the application of PEM, and the number of dofs of these two equations are remarkably small because of the application of symplectic mathematic method, hence the

5264

Y.-W. Zhang et al. / Journal of Sound and Vibration 332 (2013) 5255–5270

computational effort is considerably reduced. The first- and second-order sensitivities of u~ j can be obtained easily by computing the stationary solution of Eqs. (40) and (41) in turn, and higher-order sensitivities can also be obtained similarly. 6.3. Sensitivity analysis of K–S function According to the sensitivities of pseudo-responses obtained in the previous section, the first- and second-order sensitivities of the K–S function, i.e. the objective function of the optimization problem, can be derived step-by-step. Firstly the sensitivities of the displacement PSD of the coupled system is obtained by differentiating Eq. (7a) with respect to design variable αk and then αl as follows: ! T n ∂u~ j ðω; tÞ T ∂u~ j ðω; tÞ ∂Suu ðωÞ n ~ ~ u j ðω; tÞ þ u j ðω; tÞ ¼ ∑ (42a) ∂αk ∂αk ∂αk j ¼ v;a;c T n n ∂2 u~ j ðω; tÞ T ∂u~ j ðω; tÞ ∂u~ j ðω; tÞ ∂2 Suu ðωÞ u~ j ðω; tÞ þ ¼ ∑ ∂αk ∂αl ∂αk ∂αl ∂αk ∂αl j ¼ v;a;c ! T T n ∂u~ j ðω; tÞ ∂u~ j ðω; tÞ ∂2 u~ j ðω; tÞ n þ u~ j ðω; tÞ þ ∂αl ∂αk ∂αk ∂αl

(42b)

According to Eq. (8), the sensitivities of the acceleration PSD are as below. ∂Su€ u€ ðωÞ ∂Suu ðωÞ ∂2 Su€ u€ ðωÞ ∂2 Suu ðωÞ ¼ ω4 ; ¼ ω4 ∂αk ∂αk ∂αk ∂αl ∂αk ∂αl

(43)

Then the following equations are derived by differentiating Eq. (31): Z 80 Z 80 2 2 ∂Snq ðf Þ ∂2 Snq ðf Þ ∂ðanq ∂2 ðanq wÞ w Þ u€ u€ ¼ w2q ðf Þ u€ u€ df ; ¼ w2q ðf Þ df ∂αk ∂α ∂α ∂α ∂α ∂α k k l k l 0:5 0:5 ðn ¼ 1; 2; :::; 54;

q ¼ x; y; zÞ

(44)

Furthermore, the sensitivities of the weighted rms acceleration of point n can be obtained similarly from Eq. (34) as follows: ! nq 2 ∂awn 2 ∂ðaw Þ ¼ ∑ kq (45a) =ð2awn Þ ∂αk ∂αk q ¼ x;y;z ∂2 awn ¼ ∂αk ∂αl

∑

q ¼ x;y;z

2 kq

! 2 ∂2 ðanq ∂awn ∂awn w Þ −2 =ð2awn Þ ∂αk ∂αl ∂αl ∂αk

(45b)

Finally, the following equation is obtained by rewriting Eq. (37): 54

eρKSðαÞ ¼ ∑ eρawn ðαÞ

(46)

n¼1

Differentiating Eq. (46) gives the first- and second-order sensitivities of the K–S function as follows:  54  ∂KS ∂awn ¼ ∑ eρawn =eρKS ∂αk ∂αk n¼1 ∂2 KS ¼ ∂αk ∂αl



   ∂awn ∂awn ∂2 awn ∂KS ∂KS ρeρawn þ eρawn −ρeρKS =eρKS ∂αl ∂αk ∂αl ∂αk ∂αk ∂αl n¼1

(47.a)

54

∑

(47b)

It can be seen that there is no additional approximation is used in the calculation, thus the proposed method for analytical sensitivity analysis is accurate and retains the efficient features of random vibration analysis based on PEM and symplectic mathematical method. 6.4. Optimization of the coupled system Since the objective function is given in Eq. (37) and its first- and second-order sensitivities are derived in Eq. (47), the optimization problem (39) can be solved by any existing means. In the present paper, the ‘fmincon’ function in the optimization toolbox of MATLAB [29] is used, which is based on SQP. In each step of optimization iteration, the value of K–S function, the gradient vector consisting of the first-order sensitivities and the Hessian matrix consisting of the second-order sensitivities are updated. Using the proposed method, the optimization program can be easily operated on a personal computer.

Y.-W. Zhang et al. / Journal of Sound and Vibration 332 (2013) 5255–5270

5265

Table 1 Vehicle properties. Symbol

Definition

Value

Unit

Mt Jtθ Jtψ Jtφ Mw Jwφ Cx2 Cx1 Kx2 Kx1 h2 h3 2b1 2b2 2l1 2l2

Bogie frame mass Bogie frame rotational inertia around y-axis Bogie frame rotational inertia around z-axis Bogie frame rotational inertia around x-axis Wheel-set mass Wheel-set rotational inertia around x-axis Secondary-suspension longitudinal damping Primary-suspension Longitudinal damping Secondary-suspension Longitudinal stiffness Primary-suspension Longitudinal stiffness Distance between secondary-suspension lower surface and bogie frame center Distance between Primary-suspension lower surface and bogie frame center Primary-suspension Lateral span Secondary-suspension Lateral span Wheel-set spacing Bogie frame spacing

2400 1249 2280 1847 1900 685 4.9  106 0 3.8  105 2.74  107 0.34 0.118 0.2 2.46 2.5 17.5

kg kg m2 kg m2 kg m2 kg kg m2 N s/m N s/m N/m N/m m m m m m m

Table 2 Track properties. Symbol E Iy Iz G J ρ A Ms Js Mb Czp Cyp Czb Cyb Czf Cw Kzp Kyp Kzb Kyb Kzf Kw hr1 hr2 2br 2b l

Definition

Value

Unit 11

2.06  10 3.22  10−5 5.24  10−6 7.92  1010 3.74  10−5 7.86  103 7.72  10−3 237 123 683 5  104 5.2  104 5.88  104 4  104 3.1  104 8  104 7.8  107 2.94  107 2.4  108 5  107 6.5  107 7.8  107 9.45  10−2 8.15  10−2 0.15 1.51 0.545

Rail bending stiffness Rail moment of inertia around y-axis Rail moment of inertia around z-axis Rail shear modulus Rail polar moment of inertia Rail density Rail cross-sectional area Sleeper mass Sleeper rotational inertia around x-axis Ballast mass Rail-pad damping Rail fastening lateral damping Ballast vertical damping Ballast lateral damping Sub-grade damping Ballast shear damping Rail-pad vertical stiffness Rail fastening lateral stiffness Ballast vertical stiffness Ballast lateral stiffness Sub-grade stiffness Ballast shear stiffness Distance between rail neutral axle and top Distance between rail neutral axle and bottom Rail bottom width Distance between the two rails Sleeper spacing

Pa m4 m4 Pa m4 kg/m3 m2 kg kg m2 kg N s/m N s/m N s/m N s/m N s/m N s/m N/m N/m N/m N/m N/m N/m m m m m m

Table 3 Initial values and limits of design variables. Symbol

Definition

Initial value

Lower limit

Upper limit

Cz2 Cy2 Cz1 Cy1 Kz2 Ky2 Kz1 Ky1

Secondary-suspension vertical damping Secondary-suspension lateral damping Primary-suspension vertical damping Primary-suspension lateral damping Secondary-suspension vertical stiffness Secondary-suspension lateral stiffness Primary-suspension vertical stiffness Primary-suspension lateral stiffness

4  104 1.18  105 4  104 0 4.8  105 3.8  105 2.34  106 1.1  107

1  103 1  104 1  103 1  103 1  104 1  104 1  105 1  106

1  105 1  106 1  105 1  105 1  106 1  106 1  107 1  108

The units are N s/m for the dampings and N/m for the stiffnesses.

5266

Y.-W. Zhang et al. / Journal of Sound and Vibration 332 (2013) 5255–5270

Table 4 The 7–35th natural frequencies of the car body. Order Frequency Order Frequency Order Frequency

7 16.06 17 36.32 27 48.66

8 19.13 18 38.10 28 49.34

9 22.91 19 39.71 29 50.14

10 24.65 20 40.13 30 50.79

11 25.19 21 41.08 31 51.24

12 27.11 22 42.20 32 53.27

13 31.46 23 45.28 33 54.16

14 32,79 24 46.32 34 54.87

15 34.29 25 46.57 35 55.33

16 34.97 26 47.95

The units are Hz for the natural frequencies.

Table 5 Parameters for the American track spectrum. Parameters

Class 6 Class 5 Class 4

k Av (cm2/rad/m) Aa (cm2/rad/m) Ωc (rad/m) Ωs (rad/m)

0.25 0.0339 0.0339 0.8245 0.4380

0.25 0.2095 0.0762 0.8245 0.8209

0.25 0.5376 0.3027 0.8245 1.1312

7. Numerical examples To demonstrate the effectiveness and efficiency of the proposed optimization approach, case studies using several numerical examples are carried out in this part. The parameters of the vehicle and the track are shown in Tables 1 and 2, respectively. The initial values and the lower and upper limits of the design variables are shown in Table 3. The first 35 vibration modes of the car body are taken into account, in which the first 6 of them are free modes whose corresponding natural frequencies are zeros and the 7–35th natural frequencies are shown in Table 4. The American track spectrum, i.e. function (a) of Eq. (60), will be adopted when the track wavelength is longer than 1 m; otherwise function (b) of Eq. (60) proposed by the Chinese Railway Science Academy will be used [19]. Parameters k, Av, Aa, Ωc and Ωs take different values for different classes of track, as listed in Table 5. 8 kA Ω2c 2 > Sv ðΩÞ ¼ Ω2 ðΩv2 þΩ > 2 ðcm =rad=mÞ > > cÞ > < kA Ω2c 2 −3:15 Sa ðΩÞ ¼ Ω2 ðΩa2 þΩ (48) ðaÞ ðmm2 =cycle=mÞ ðbÞSv ðf Þ ¼ 0:036f 2 ðcm =rad=mÞ cÞ > > > 2 > 4kAv Ωc > ðcm2 =rad=mÞ : Sc ðΩÞ ¼ ðΩ2 þΩ2 ÞðΩ 2 þΩ2 Þ c

s

The American track spectrum class 6 and vehicle velocity of 200 km/h are adopted in the calculations unless specified otherwise. 7.1. Optimization results of the coupled vehicle–track system The results of the riding-comfort optimization of the coupled vehicle–track system are discussed in this section. Fig. 5 gives the original and optimized weighted rms accelerations of all the design points in three rows. In general, the weighted rms accelerations of points near the windows are larger than those of points away from the windows. The global peak value appears at point A that is near the center of row 3, while the local peak values of rows 1 and 2 appear respectively near the center or at the end of the rows. The weighted rms accelerations of all the design points are greatly reduced after optimization, which means the riding-comfort of the passengers at different locations is overall improved. Besides, the differences among the 3 rows are also significantly reduced. The original and optimized peak values of weighted rms acceleration of the 3 rows are shown in Table 6. It can be seen that the points whose values are larger in the original design are greater improved. At peak point A, the weighted rms acceleration reduced by 58.34 percent after optimization, which demonstrates the effectiveness of the proposed method. As a further investigation, the acceleration PSDs of the design points in the most sensitive frequency range of human body, i.e. 4.0–12.5 Hz in the vertical direction and 0.5–2.0 Hz in the lateral direction, are calculated. Only the results of design points in row 3, whose weighted rms accelerations are relatively larger, are shown in Fig. 6. It can be seen in Fig. 6(a) that in the original design, great differences exist among the vertical accelerations of the design points. In the frequency range less than 5.0 Hz, the vertical acceleration PSDs at the ends of the car body are significantly higher than those at the center, while the reverse is true in the frequency range more than 5.0 Hz. After the optimization as shown in Fig. 6(b), the differences of the results are greatly reduced, and the overall values also greatly decreased. Fig. 6(c) and (d) give the lateral acceleration PSDs of the design points before and after the optimization. The values at the ends of the car body are always larger than those at the center in the concerned frequency range, and the results reduced less after optimization compared

Y.-W. Zhang et al. / Journal of Sound and Vibration 332 (2013) 5255–5270

Fig. 5. Weighted rms accelerations of the design points. (a) Original, (b) optimized. Table 6 Comparison of original and optimized results of 3 rows.

Original (m/s2) Optimized (m/s2) Reduction (%)

Row 1

Row 2

Row 3

1.0172 0.7211 29.11

1.2391 0.7298 41.10

1.8462 0.7691 58.34

Fig. 6. Acceleration PSDs of row 3. (a, c) Original, (b, d) optimized.

5267

5268

Y.-W. Zhang et al. / Journal of Sound and Vibration 332 (2013) 5255–5270

with the vertical direction. The original and optimized peak values of point B are respectively 0.0863 m2/s4/Hz and 0.0726 m2/s4/Hz, only being reduced by 15.87 percent. In contrast, the reduction of vertical acceleration PSD in Fig. 6(a) and (b) are more than 50 percent. 7.2. Comparison with results based on a rigid body model In order to compare the optimization effect using the proposed FE model with that using the rigid body model, a 3-dementional coupled vehicle–track analysis model is established with the car body regarded as a rigid mass as shown in Fig. 1. The parameters, including the mass and rotational inertia around the three axes of the car body and the distance between body center and secondary-suspension upper surface, are calculated by ANSYS. Then, the parameters of its primary- and secondary-suspensions are optimized to minimize the weighted rms acceleration of the car body, using the same approach as in the FE model. Finally, the optimized parameters of the primary- and secondary-suspensions are input into the developed FEM-based model to calculate the weighted rms accelerations of the design points for comparison, whose results are shown in Fig. 7. It can be seen that the optimization effect of the rigid body model is relatively not obvious. The peak value of point A in Fig. 7 is 1.5595 m/s2, which is only a 15.53 percent reduction than the original design. In comparison, the peak value reduced by 58.34 percent when the FEM-based model is used, as shown in Table 6. It is clear

Fig. 7. Optimized results using a spring-mass-damper model.

Fig. 8. Optimized weighted rms accelerations of row 3. (a) Under different velocities, (b) under different classes of track irregularity.

Y.-W. Zhang et al. / Journal of Sound and Vibration 332 (2013) 5255–5270

5269

Fig. 9. Optimized dimensionless parameters of vehicle suspension. (a) Under different velocities, (b) under different classes of track irregularity.

that the elasticity of the car body has a great effect on the optimization results and the results obtained by a rigid body model are unreliable.

7.3. Influence of vehicle running conditions The vehicle velocity and the class of track irregularity are important factors for the response of the coupled system, whose influences on the optimization are discussed in this section. Fig. 8 gives the optimized weighted rms accelerations of design points in row 3 at different vehicle velocities and classes of track irregularity. It can be seen that the higher the vehicle velocity is or the lower the class of track irregularity is, the stronger the responses of the vehicle will occur. The same results are indicated in Ref. [19] and the following conclusion can be drawn: in very bad vehicle running conditions, the riding comfort cannot be completely improved only by means of optimizing the passive vehicle suspension system, thus more effective means, e.g. active control, should be considered to improve the vehicle riding comfort. Fig. 9 gives the optimized dimensionless parameters of the vehicle suspension system under different running conditions, i.e. the optimized parameters divided by the mean of them, in which the design variable vector in Eq. (35) is used as abscissa. It can be seen that the optimized values of these parameters are irregular for either the vehicle velocity or the class of track irregularity. Therefore, in the design of the riding comfort of the vehicle, the vehicle running conditions need to be critically taken into account.

8. Concluding remarks In this research, an innovative optimization method for riding comfort of coupled vehicle–track system is proposed based on the PEM and the symplectic mathematical method, in which the elaborate FEM-based head car model of a high-speed train is developed. The stationary equations of motion of the coupled system with remarkably small number of dofs and their first- and second-order sensitivity equations are derived, which can not only considerably reduce the computational effort, but also obtain higher accuracy of the solution. The proposed method can meet the needs of practical engineering design and the optimization program can operate on a personal computer. In particular, since an accurate and highly efficient method is developed for the FEM-based random vibration analysis of the coupled vehicle–track system, the optimization problem proposed in the present paper can be solved using not only the gradient-based methods but also the intelligent optimization approaches.

5270

Y.-W. Zhang et al. / Journal of Sound and Vibration 332 (2013) 5255–5270

Acknowledgments Supports for this work from the National Natural Science Foundation of China (No. 11172056) and the National 973 Plan Project (No. 2010CB832704). References [1] T. Dahlberg, Parametric optimisation of a 1-dof vehicle travelling on a randomly profiled road, Journal of Sound and Vibration 55 (1977) 245–253. [2] T. Dahlberg, An optimised speed controlled suspension of a 2-dof vehicle travelling on a randomly profiled road, Journal of Sound and Vibration 64 (1979) 541–546. [3] P. Eberhard, D. Bestle, U. Piram, Optimisation of damping characteristics in nonlinear dynamic systems, First World Congress of Structural and Multidisciplinary Optimisation (1995) 863–870. [4] T.C. Sun, Albert C.J. Luo, H.R. Hamidzadeh, Dynamic response and optimization for suspension systems with nonlinear viscous damping, Journal of Multibody Dynamics 214 (2000) 181–187. [5] L.F.P. Etman, R.C.N. Vermeulen, J.G.A.M. van Heck, A.J.G. Schoofs, D.H. van Campen, Design of a stroke dependent damper for the front axle suspension of a truck using multibody system dynamics and numerical optimization, Vehicle System Dynamics 38 (2002) 85–101. [6] P.S. Els, P.E. Uys, Investigation of the applicability of the dynamic-Q optimisation algorithm to vehicle suspension design, Mathematical and Computer Modelling 37 (2003) 1029–1046. [7] A.F. Naude, J.A. Snyman, Optimization of road vehicle passive suspension systems-part 1: optimisation algorithm and vehicle model, Applied Mathematical Modelling 27 (2003) 249–261. [8] M.J. Thoresson, P.E. Uys, P.S. Els, J.A. Snyman, Efficient optimisation of a vehicle suspension system, using a gradient-based approximation method, Part 1: mathematical modelling, Mathematical and Computer Modelling 50 (2009) 1421–1436. [9] A.E Baumal, J.J McPhee, P.H Calamai, Application of genetic algorithms to the design optimization of an active vehicle suspension system, Computer Methods in Applied Mechanics and Engineering 163 (1998) 87–94. [10] R. Alkhatib, G.N. Jazar, M.F. Golnaraghi, Optimal design of passive linear suspension using genetic algorithm, Journal of Sound and Vibration 275 (2004) 665–691. [11] T.J. Yuen, R. Raml, Comparison of computational efficiency of MOEA/D and NSGA-II for passive vehicle suspension optimization, Proceedings of the 24th European Conference on Modelling and Simulation, 2010, pp. 219–225. [12] M. Gobbi, G. Mastinu, C. Doniselli, L. Guglielmetto, E. Pisino, Optimal and robust design of a road vehicle suspension system, Vehicle System Dynamics (Supplement 33) (2000) 3–22. [13] Y.G. Kim, C.K. Park, T.W. Park, Design optimization for suspension system of high speed train using neural network, JSME International Journal Series C 46 (2003) 727–735. [14] I. Rajendran, S. Vijayarangan, Simulated annealing approach to the optimal design of automotive suspension systems, International Journal of Vehicle Design 43 (2007) 11–30. [15] A. Kuznetsov, M. Mammadov, I. Sulta, E. Hajilarov, Optimization of a quarter-car suspension model coupled with the driver biomechanical effects, Journal of Sound and Vibration 330 (2011) 2937–2946. [16] C.H. Li, G.F. Wan, Random vibration of track structure, Journal of the China Railway Society 20 (1998) 97–101 (in Chinese). [17] X. Lei, N.A. Noda, Analyses of dynamic response of vehicle and track coupling system with random irregularity of track vertical profile, Journal of Sound and Vibration 258 (2002) 147–165. [18] D. Dumitriu, D. Baldovin, T. Sireteanu, Numerical simulation of the vertical interaction between railway vehicle and track, Proceedings of the Annual Symposium of the Institute of Solid Mechanics, Bucharest (CD-ROM), 2006. [19] F. Lu, D. Kennedy, F.W. Williams, J.H. Lin, Symplectic analysis of vertical random vibration for coupled vehicle–track systems, Journal of Sound and Vibration 317 (2008) 236–249. [20] J.H. Lin, Y.H. Zhang, Pseudo Excitation Method in Random Vibration, Science Press, Beijing, 2004 (in Chinese). [21] J.H. Lin, Y.H. Zhang, Vibration and Shock Handbook, Chapter 30:Seismic Random Vibration of Long-span Structures, CRC Press, Boca Raton, FL, 2005. [22] W.X. Zhong, F.W. Williams, Wave problems for repetitive structures and symplectic mathematics, Proceedings of the Institution of Mechanical Engineers, Part C 206, 1992, pp. 371–379. [23] J.H. Lin, Y. Fan, P.N. Bennett, F.W. Williams, Propagation of stationary random waves along substructural chains, Journal of Sound and Vibration 180 (1995) 757–767. [24] W.T. Xu, J.H. Lin, Y.H. Zhang, D. Kennedy, F.W. Williams, Pseudo-excitation-method -based sensitivity analysis and optimization for vehicle ride comfort, Engineering Optimization 41 (2009) 699–711. [25] International Organization for Standardization ISO 2631-1, Mechanical vibration and shock - evaluation of human exposure to whole-body vibration part 1: general requirements, 1997. [26] Z.S. Yu, Vehicle theory, China Machine Press, Beijing, 2000 (in Chinese). [27] G. Kreisselmeier, R. Steinhauser, Systematic control design by optimizing a vector performance index, Proceedings of IFAC Symposium on Computer Aided Design of Control Systems, Zürich (1979) 113-117. [28] J. Barthelemy, M.F. Riley, Improved multilevel optimization approach for the design of complex engineering systems, AIAA Journal 26 (1988) 353–360. [29] http://www.mathworks.cn/help/toolbox/optim/index.html.