Vibration Response Law of Aircraft Taxiing under Random Roughness Excitation

: Understanding the vibration response pa tt erns of aircraft taxiing under runway roughness excitation is crucial for aircraft design and runway performance evaluation. In this paper, we establish pavement roughness models for both asphalt and concrete surfaces, taking into account their unique structural characteristics. We construct a six-degrees-of-freedom aircraft model using multi-rigid-body system dynamics theory and employ the pseudo-excitation method to examine the in ﬂ uence of pavement roughness types on the steady vibration response of aircraft taxiing at constant speeds. Furthermore, we analyze the non-stationary vibration response pa tt erns of aircraft during takeo ﬀ and landing taxiing using the method of instantaneous frequency response function with space frequency. Lastly, we explore the e ﬀ ect of stochastic structural parameters on aircraft vibration response using the Monte Carlo method. Our ﬁ ndings reveal that the roughness power spectrum di ﬀ ers between asphalt and concrete pavements, and the established roughness models in this paper demonstrate a strong ﬁ t (R2 > 0.95). The type of pavement roughness has a relatively minor impact on the power spectral density distribution of the aircraft vibration response, suggesting that the same roughness model can be used for both asphalt and concrete pavements when high accuracy is not required. The power spectral density distribution of aircraft vibration response varies across di ﬀ erent motion a tt itudes, with the vibration response during landing being signi ﬁ cantly larger than that during takeo ﬀ . Among the aircraft structural parameters, the randomness of the sprung mass has the most substantial e ﬀ ect on the aircraft vibration response, potentially causing the variation coe ﬃ cient of dynamic load on the front landing gear to exceed 0.11. Tire sti ﬀ ness comes next, which can lead to the variation coe ﬃ cient of dynamic load on the main landing gear reaching 0.07. The results have a guiding role in optimizing aircraft structure design and ensuring pavement performance.


Introduction
The airport runway functions as a crucial support platform for aircraft during taxiing, takeoff, and landing. Owing to a combination of factors such as construction variability, external environment, aircraft loading, and other uncertainties, the pavement exhibits random roughness. Pavement roughness is the vertical deviation of the pavement surface relative to the ideal plane. When subjected to random roughness excitation, high-speed taxiing aircraft experience random vibrations, which can negatively impact passenger comfort, accelerate fatigue damage to aircraft components, cause dynamic damage to runway structures, and even compromise flight safety [1]. As a result, comprehending the vibration response patterns of aircraft under random roughness excitation is vital for optimizing aircraft structural design, ensuring runway structure safety, and maintaining operational safety standards.
Numerous studies have demonstrated that road surface or pavement roughness can be characterized as an ergodic, steady Gaussian stochastic process with zero mean, and its statistical properties can be represented by power spectral density (PSD) [2]. In the field of road research, various roughness models have been developed using PSD through the analysis of measured road roughness data [3][4][5][6][7]. However, runway roughness studies are limited by the available data samples, resulting in few investigations into the power spectrum of runway roughness. For instance, Houbolt proposed a PSD expression for airport pavement roughness based on measuring the relative elevation data of longitudinal sections at only two runways, which lacks representativeness [8]. Most researchers have directly adopted road surface roughness models when studying runway roughness [9][10][11][12]. Nevertheless, this approach can introduce significant deviations in the analysis of runway roughness due to the considerable differences between runways and roads concerning structural forms and service object characteristics [13].
In aircraft vibration response research, the majority of scholars utilize measured or simulated pavement roughness samples for time-domain simulations. For instance, Durán and Júnior calculated aircraft acceleration under measured roughness excitation using ProFAA and evaluated pavement roughness accordingly [14]. Dong et al. simulated the pavement roughness curve with the harmonic superposition method and calculated the taxiing dynamic load of a four-degrees-of-freedom aircraft using the Full method [15]. Cheng et al. and Shi et al. employed the filtered white noise method or harmonic superposition method to generate pavement roughness, and determined the dynamic load coefficient of the six-degrees-of-freedom (6-DOF) model at different speeds using Simulink [16,17]. Liu et al. used the cosine function as the roughness input and analyzed the superposition effect of aircraft vibration with automatic dynamic analysis of mechanical systems (ADAMS) [18]. Qian et al. applied PSD to reconstruct runway roughness and investigated the aircraft taxiing dynamic load using ADAMS [19]. However, measured or simulated single roughness curves only represent pavement roughness excitation at a certain probability level during the runway's long-term service. Single deterministic time history analysis based on this does not adequately reflect the random vibration characteristics of the aircraft.
To address this limitation, several researchers have conducted frequency domain analyses based on random vibration theory, establishing a direct correlation between the power spectrum input of pavement roughness and the power spectrum output of aircraft vibration response. For instance, Kirk and Perry used the PSD method to analyze the vibration response of a simplified multiple-degrees-of-freedom aircraft under stationary random excitation [20]. Tung solved the single-degree-of-freedom vibration problem of the aircraft using random vibration theory [21]. Nie and Kortum employed the equivalent linearization method to address the nonlinear parameters of the landing gear and examined random vibration during taxiing using the PSD method [22]. Wang applied the pseudo-excitation method (PEM) to compute the dynamic response of unmanned aerial vehicle landing gear under non-stationary random excitation [23]. Liu et al. obtained the vibration response of a two-degrees-of-freedom aircraft-runway coupling system under roughness excitation using the PEM [24]. However, most of these studies employed simplified aircraft models that were unable to capture the characteristics of a large wheelbase and significant spatial vibration of the aircraft.
In this paper, we address the unique structural characteristics of asphalt and concrete pavements by establishing corresponding pavement roughness models using PSD. We derive the motion balance equations of a 6-DOF aircraft under roughness excitation following D'Alembert's principle. We then use the PEM, the method of instantaneous frequency response function with spatial frequency, and the Monte Carlo method to calculate the PSD distribution of acceleration at the center of gravity (CGA), acceleration at the pilot station (PSA), dynamic load coefficient of the nose landing gear (NGDLC), and dynamic load coefficient of the main landing gear (MGDLC) under various operating conditions. This approach allows us to reveal the vibration response laws of aircraft taxiing under different pavement types, motion attitudes, and stochastic structural parameters. The results can be applied to optimize aircraft design methods and improve the accuracy of pavement roughness evaluation.

Aircraft Nonlinear Ground Dynamics Model
Throughout the taxiing process, pavement roughness excitation is transmitted through the tire to the landing gear, which absorbs some of the energy and transfers the vibration to the fuselage. The aircraft is a typical nonlinear model due to the strong nonlinearity of both the tire and the landing gear. In frequency domain dynamic response analysis, nonlinear elements are typically represented by equivalent linear elements. Therefore, in this section, the nonlinear elements are initially represented by equivalent linear letters (with the equivalent linearization process described in the next section), and the nonlinear aircraft model is established.
Initially, based on the theory of multi-rigid-body system dynamics, the aircraft is simplified into two parts: the sprung mass and unsprung mass. Additionally, the mechanical properties of the landing gear buffer and the tire are simulated as a spring-plus-damping model, and the lift, , is applied to the center of gravity of the aircraft during takeoff and landing on the runway in a straight line. This yields the ground dynamics model of the 6-DOF aircraft [25], as illustrated in Figure 1.  At any given moment, the aircraft is in a state of dynamic equilibrium. Therefore, the motion balance equations for each degree of freedom of the aircraft ground dynamics model under pavement roughness excitation can be derived using the D'Alembert principle.
The vertical vibration balance equation for the sprung mass is shown as Equation (1).
The pitch rotation balance equation for the sprung mass is shown as Equation (2).
The lateral rotation balance equation for the sprung mass is shown as Equation (3).
The vertical vibration balance equations for the unsprung mass of the front, right rear, and left rear landing gear are shown as Equations (4)-(6), respectively.
Combining Equations (1)-(6), the matrix equation for the spatial motion of the 6-DOF aircraft ground dynamics model is established as shown in Equation (7).

[ ][ ] + [ ][ ] + [ ][ ] = [ ][ ] + [ ][ ]
where denotes the aircraft sprung mass; , , and denote the unsprung mass of the front, right rear, and left rear landing gear, respectively; , , and denote the damping coefficients of the front, right rear, and left rear landing gear buffer, respectively; , , and denote the stiffness coefficients of the front, right rear, and left rear landing gear buffer, respectively; , , and denote the damping coefficients of the front, right rear, and left rear landing gear tire, respectively; , , and denote the stiffness coefficients of the front, right rear, and left rear landing gear tire, respectively; and denote the horizontal distances from the rear and front landing gears to the x-axis of the aircraft, respectively; and denote the horizontal distances from the right rear and left rear landing gear to the y-axis of the aircraft, respectively; and denote the rotational inertia of the aircraft around the x-axis and y-axis, respectively; , , and denote the vertical displacement, pitch rotation displacement, and lateral rotation displacement of the sprung mass, respectively; , , and denote the vertical displacements of the front, right rear, and left rear landing gear unsprung mass, respectively; , , and denote the roughness excitation of the front, right rear, and left rear landing gear tire, respectively; is the lift force acting on the aircraft; [ ] is the mass matrix of the aircraft model, which is a matrix of 6 rows and 6 columns.

Statistical Linearization of the Nonlinear System
In the frequency domain, nonlinear systems typically require equivalent linearization, which can be achieved using the statistical linearization method [20,21]. The objective of statistical linearization is to identify an equivalent linear system that has a minimal mean square error when compared to the original nonlinear system. The statistical linearization method can conveniently transform the nonlinear parameters into linear parameters, so the frequency domain analysis of random vibration can be carried out using the existing random vibration method. Therefore, the method is employed to deal with the nonlinear parameters of the aircraft.

Nonlinear Function Expression
(1) Nonlinearity of landing gear elastic element In general, the landing gear spring force increases at an increasing rate as the moving displacement rises, as shown in Figure 2. Multiple regression analysis is used to examine different fitting models. Considering the fitting effect and the complexity of the formula, a cubic polynomial can be applied to describe the nonlinear characteristics of Figure 2. The fitted formula is shown in Equation (8).
where is the spring force; is the moving displacement; and , , and are fitting parameters.
(2) Nonlinearity of landing gear damping element Generally speaking, the damping elements of the landing gear have nonlinear characteristics, which can be expressed by piecewise functions at different velocities, as shown in Equation (9). The relationship curve between the damping force and the moving speed is shown in Figure 3.
where ( ) is the damping force; is the moving speed; and and are fitting parameters.

Calculation of Equivalent Stiffness and Equivalent Damping by Statistical Linearization
In a nonlinear vibration system with n degrees of freedom, the mass of each inertial element is , , . . . , , respectively, and the excitation force ( ) = { ( ), ( ), . . . , ( )} acting on the system is an n-dimensional joint Gaussian process vector with zero mean. Then, in the inertial coordinate system , , . . . , , the differential equation of motion for this system is shown in Equation (10).
where [ ] = [ ] is the inertial matrix of the system in the coordinate system , , . . . , ; ( , ) is the sum of elastic force and damping force acting on ; and ( ) is the exciting force applied to . Suppose the transformation between coordinates , , . . . , and coordinates , , . . . , can be expressed as Equation (11).
Taking Equation (11) into Equation (10), the equation of motion of the system represented by the relative coordinates , , . . . , is obtained, as shown in Equation (12).
where Assume that the statistically equivalent linear system of the initial nonlinear system satisfies the differential equation shown in Equation (13).
where [ ] and [ ] are the equivalent linear stiffness matrix and the equivalent damping matrix, respectively.
Since ( ) is the elastic force of elastic element i, then ( )/ is the stiffness of elastic element i. Therefore, [ ( )/ ] should be the statistical equivalent stiffness coefficient of elastic element i, which can be expressed as Equation (15a).
Similarly, the statistically equivalent damping coefficient, , of each damping element can be derived as shown in Equation (15b).
Consequently, the equivalent stiffness coefficient and the equivalent damping coefficient can be treated as the expected values of Equations (8) and (9) for varying displacement and velocity conditions. In reality, the probability distribution of the vertical vibration displacement of the landing gear and tire during aircraft taxiing becomes relevant. After estimating the probability distribution, Equations (15a) and (15b) can be employed to obtain the statistically equivalent linear coefficient for each element.

Pavement Roughness Models
The Federal Aviation Administration (FAA) gathered and published longitudinal profile elevation data for 37 runway centerlines [26]. Our team also collected pavement roughness data for 20 runways in China using a vehicle-mounted laser profile meter and a Beidou high-precision positioning system, which can obtain the true full-wave runway profile [27]. Figure 4 presents the pavement roughness and the data collection system used by our team. The basic information for the 20 Chinese runways collected by our team is presented in Table 1. Upon analyzing the data, we discovered that asphalt pavement and concrete pavement display distinct digital characteristics in their roughness sequences. Figure 5 illustrates the power spectrum of measured roughness data for the asphalt pavement (a) at Atlantic City International Airport and the concrete pavement (b) at Taizhou Luqiao Airport. In Figure 5, n is the spatial frequency, and its unit is cycle/m. Dividing n by the unit spatial frequency (1 cycle/m), we can obtain the unitless spatial frequency nu. Lgnu is the logarithm of nu to the base 10. There are notable differences in the power spectrum between asphalt pavement and concrete pavement. The PSD distribution curve of asphalt pavement roughness is relatively smooth, while that of concrete pavement decreases significantly faster at high frequencies. This phenomenon can be attributed to the fact that concrete runways consist of slabs as units, with joints and uneven surfaces between two adjacent slabs, while the surface within a single unit is generally flat. In contrast, asphalt runways are flexible pavements and are relatively more continuous in terms of roughness.   The roughness power spectrum of 20 asphalt runways, both domestic and international, is analyzed. After repeated comparison, the asphalt pavement roughness model takes reference from the Sussman model [4], as shown in Equation (16).
where ( ) is the PSD function of pavement roughness; C is the roughness coefficient, and the larger the C value, the rougher the pavement; w is the frequency index, and the larger w, the more significant the long-wave component; and α is the truncation frequency, and the PSD tends to be smooth as n is less than α.
As Figure 5a shows, the asphalt pavement roughness model established in this paper exhibits a strong fitting effect, with R2 values above 0.95.

Cement Pavement Roughness Model
Since the slab length of concrete runways is generally 5 m, there is an inflection point in the power spectrum curve in Figure 5b around n = 0.1 m −1 . By analyzing the roughness power spectrum of 37 concrete runways both domestic and international, this paper utilizes a roughness model with discontinuous segmentation points to represent the concrete pavement, as depicted in Equation (17).
where is the segmentation point, = 0.1 m ; and are the roughness coefficients of the first and second sections, respectively; and and are the frequency indexes of the first and second segments, respectively.
As Figure 5b shows, the cement pavement roughness model established in this paper exhibits a strong fitting effect, with R2 values greater than 0.95.

Pavement Roughness Excitation Input
By utilizing the relationship between temporal frequency f, time angular frequency ω, and spatial frequency n (as described in Equations (18a) and (18b)), the pavement roughness model in the spatial domain can be transformed into pavement roughness excitation input in the time domain.
where v is the aircraft taxiing speed. The asphalt pavement roughness excitation input in the time domain can be expressed as Equations (19a) and (19b).
The cement pavement roughness excitation input in the time domain can be expressed as Equations (20a) and (20b).

Solution Method for Steady Random Vibration of Uniform Taxiing Aircraft
Under pavement roughness excitation, the aircraft taxiing at a constant speed generates steady random vibration. To address this, this paper utilizes the PEM to determine the steady vibration response of the aircraft [28]. This method is highly efficient and theoretically accurate, making it widely used in numerous engineering fields [29].
Firstly, the characteristic decomposition of the PSD matrix [ ( )] of the pavement roughness excitation input is carried out to construct the virtual harmonic excitation, as shown in Equation (21).
where is the eigenvalue of matrix [ ( )] ; { } is the eigenvector corresponding to ; and superscript * represents the conjugate operation.
Secondly, the frequency response function matrix [ ( )] of the aircraft system is obtained by the Fourier transform of Equation (7), as shown in Equation (22). Furthermore, the virtual harmonic response ( , ) of the aircraft can be calculated according to Equation (23).
Finally, the PSD matrix [ ( )] of the actual aircraft response is solved using Equation (24).
where the superscript T represents the transpose operation. Using MATLAB, the vibration response of the aircraft can be solved simply and quickly based on the above steps.

Solution Method for Non-Stationary Random Vibration of Takeoff and Landing Taxiing Aircraft
In the spatial domain, pavement roughness is steady, but the radial motion of the aircraft during takeoff and landing taxiing is a non-uniform process. As a result, the system differential equation becomes a time-varying differential equation, leading to a nonstationary response of the aircraft system. For a non-stationary random vibration system, the vibration response cannot be calculated directly using the frequency response function H(ω). Therefore, this paper employs the method of instantaneous frequency response function with spatial frequency [30] combined with the PEM to solve the transient power spectrum of the vibration response during aircraft takeoff and landing taxiing.
Lin et al. demonstrated through field measurements and theoretical research that aircraft takeoff and landing taxiing on the runway can be regarded as a uniformly accelerated linear motion process. The corresponding velocity v and displacement s can be expressed as described in Equations (25a) and (25b) [31].
where is the initial velocity of uniform variable motion, which is taken as 0 during the aircraft takeoff taxiing, and the landing speed during aircraft landing taxiing; a is the acceleration, which is positive during the aircraft takeoff taxiing and negative during the aircraft landing taxiing.
The time-space frequency relationship of the uniformly varying motion during aircraft takeoff and landing taxiing can be obtained by substituting Equations (25a) and (25b) into Equation (18), as shown in Equation (26).
Based on the characteristic that the frequency response function H(ω) of the aircraft system in non-stationary random vibration is constant at temporal frequency, the instantaneous frequency response function ( , ) at spatial frequency can be derived from the relationship between temporal frequency and spatial frequency. Therefore, by substituting Equation (26) into Equation (22), the instantaneous spatial domain frequency response function matrix of the aircraft system can be expressed as Equation (27).
After obtaining Equation (27), the instantaneous spatial PSD of the aircraft vibration response can be calculated following the PEM in Section 3.1.

Aircraft Parameters
In this paper, the B737-800 is selected as the representative model. Based on the values given in Table 2, the equivalent linear stiffness coefficient for the nose landing gear is 58,027 N/m and the equivalent linear damping coefficient is 131,500 N·s/m. The equivalent linear stiffness coefficient for the main landing gear is 996,531 N/m and the equivalent linear damping coefficient is 572,500 N·s/m.

. Pavement Roughness Parameters
The power spectrum of the measured data from the 57 runways is fitted using the asphalt and cement pavement roughness models developed in this paper. The distributions of the fitting parameters for asphalt and cement pavements are presented in Tables  3 and 4, respectively. The values in brackets in the tables are the parameter values at the maximum probability density, i.e., the typical asphalt pavement and cement pavement roughness parameters. In order to increase the representativeness of the study, typical roughness parameters of domestic runways for both asphalt and cement pavements are employed.

Vibration Response Analysis of Aircraft Taxiing under Different Pavement Types
The PSD distributions of CGA, PSA, NGDLC, and MGDLC under typical roughness excitations of both asphalt pavement and concrete pavement during a 100 km/h taxiing of a B737-800 aircraft are shown in Figure 6a to Figure 6d, respectively. Overall, the PSD distribution of each vibration response is similar for both pavement types, as shown in Figure 6. The red and black lines show a similar trend of rising and falling, and the sensitivity frequencies corresponding to each peak are almost the same. However, the vibration response corresponding to concrete pavement roughness at around 2.6 Hz decreases dramatically. This is attributed to the segmental characterization of the concrete pavement roughness model, which has few broken slabs or cracks within the concrete slab and high flatness.

Vibration Response Analysis of Aircraft Taxiing under Different Motion Attitudes
The runway takeoff distance is set at 2500 m, and the B737-800 aircraft takes off from the ground at a speed of 260 km/h with an acceleration of 1.04 m/s 2 throughout the process. The PSD distributions of CGA, PSA, MGDLC, and NGDLC under typical asphalt pavement roughness excitation during takeoff are shown in Figure 7a-d. The power spectral density distribution of the system response at any instant on the runway is consistent with the uniform steady-state response. As the takeoff distance increases, the taxiing speed of the aircraft increases, and the PSD of the vibration response generally increases, although not monotonically. For example, the PSD of CGA is locally reduced near the airframe resonance frequency within the distance range of 1000-1500 m. The most sensitive frequency of the aircraft remains unchanged in each spatial segment, which is due to the constant natural frequency of the aircraft structure. The landing speed of the B737-800 is set at 220 km/h, the landing distance is set at 1000 m, and the acceleration is −1.86 m/s 2 . The PSD distributions of CGA, PSA, MGDLC, and NGDLC under typical asphalt pavement roughness excitation during landing are shown in Figure 8a-d. In general, the closer to the landing site, the greater the PSD of the vibration response. Thus, during the landing process, the first half of the runway roughness has a stronger impact on the aircraft vibration. Compared to the takeoff process, the PSD of the aircraft vibration response during landing is greater because the aircraft's mass decreases due to fuel consumption during the flight, and the aircraft's acceleration during landing deceleration is larger. Therefore, the roughness of the landing runway should be of greater concern to managers than that of the takeoff runway.

Vibration Response Analysis of Aircraft Taxiing under Stochastic Structural Parameters
Due to the large number of degrees of freedom in the aircraft model, it is complicated to theoretically derive the effect of stochastic structural parameters on the vibration response. The Monte Carlo method can build a mathematical model of the input and output variables, apply the probability distribution of each variable to random sampling, and calculate the output variables using the randomly generated input variables. By repeating the random simulations, the probability distribution function of a certain output variable can be obtained. Additionally, the accuracy of the output variable's probability distribution function increases with the number of simulations. With the aid of the MATLAB (R2019a) software and the Monte Carlo method, this paper analyzes the variation in each vibration response when the probability density function of aircraft structural parameters is modeled as a normal distribution.
The variation coefficients of the sprung mass, rotational inertia, stiffness coefficient of the landing gear buffer, damping coefficient of the landing gear buffer, stiffness coefficient of the tire, and damping coefficient of the tire are all set to 0.05. By using the Monte Carlo method to randomly sample 1000 times at a time, the variation coefficient distribution of the variance of CGA, PSA, NGDLC, and MGDLC can be obtained, as shown in Figure 9. The influence of the sprung mass is the largest, having a significant effect on the four vibration responses, with variation coefficient values exceeding 0.05. Among them, the sprung mass has the greatest impact on NGDLC, reaching 0.113. The rotational inertia mainly affects PSA, with variation coefficients above 0.05, but has little effect on the other three vibration responses. The damping coefficient of the landing gear buffer has a greater impact than the stiffness coefficient of the landing gear buffer. The damping coefficient of the front landing gear affects PSA and NGDLC more, while the damping coefficient of the main landing gear has a greater impact on CGA and MGDLC, as the center of gravity is closer to the main landing gear. The situation with the tire is the opposite of that of the landing gear, with the tire damping coefficient having a much smaller effect than the tire stiffness. The stiffness coefficient of the tire has a great influence on the dynamic load coefficient of the landing gear where it is located, but it has little effect on the other landing gear. At the same time, the stiffness coefficient of the main landing gear tire has a greater impact on the fuselage vibration than the front landing gear, indicating that more attention should be paid to the roughness excitation applied to the main landing gear.

Conclusions
In this paper, pavement roughness models for both asphalt and cement pavements are established based on their structural characteristics. The PEM, the method of instantaneous frequency response function with spatial frequency, and the Monte Carlo method are used to comprehensively analyze the influence of different pavement types, different motion attitudes, and stochastic structural parameters on the vibration response of a 6-DOF aircraft model under random roughness excitation in the frequency domain. The results of the study can be adopted for runway roughness evaluation and aircraft structure design. The main conclusions from this study are as follows: (1) The power spectra of asphalt pavement and cement pavement roughness are different. The corresponding pavement roughness models were developed in this paper and the measured data demonstrates that the fit is excellent, with R2 values greater than 0.95. (2) The pavement type has little influence on the PSD distribution of the aircraft vibration response. When the accuracy requirements are not high, the same pavement roughness model can be used for both cement and asphalt pavements. (3) The vibration response of the aircraft during landing is greater than that during takeoff under the same conditions. Managers should therefore pay more attention to the roughness of the landing runway compared to the takeoff runway. (4) The impact of stochastic structural parameters on the aircraft's vibration response was analyzed using the Monte Carlo method. It was found that the randomness of the sprung mass has the greatest effect on the vibration response and can result in NGDLC reaching 0.113, followed by the tire stiffness and the aircraft's rotational inertia. Funding: This research was funded by the Natural Science Foundation of Shanghai, grant number 23ZR1466300.