# Construction of Three-Dimensional Road Surface and Application on Interaction between Vehicle and Road

Authors: Lu Yongjie; Huai Wenqing; Zhang Junning
Journal: Shock and Vibration
 (2018)

---

## Abstract

The quantitative description is given to three-dimensional micro and macro self-similar characteristics of road surface from the perspective of fractal geometry using FBM stochastic midpoint displacement and diamond-square algorithm in conjunction with fractal characteristics and statistical characteristics of standard pavement determined by estimation method of box-counting dimension. The comparative analysis between reconstructed three-dimensional road surface spectrum and theoretical road surface spectrum and correlation coefficient demonstrate the high reconstruction accuracy of fractal reconstructed road spectrum. Furthermore, the bump zone is taken as an example to reconstruct a more arbitrary 3D road model through isomorphism of special road surface with stochastic road surface model. Measurement is taken to assume the tire footprint on road surface to be a rectangle, where the pressure distribution is expressed with mean stiffness, while the contact points in the contact area are replaced with a number of springs. Two-DOF vehicle is used as an example to analyze the difference between three-dimensional multipoint-and-plane contact and traditional point contact model. Three-dimensional road surface spectrum provides a more accurate description of the impact effect of tire on road surface, thereby laying a theoretical basis for studies on the dynamical process of interaction of vehicle-road surface and the road friendliness.


---

## Body

## 1. Introduction

To effectively achieve the active safety control of vehicle, it is essential to attain the status and tire-pavement friction coefficient for a running vehicle, which is based on the consideration of interaction between vehicle and road surface in real time. However, the fact is that the inadequate description of road surface morphology details by existing road surface spectrum models is restricting the study on vehicle-road interaction. In the field of road engineering, road transportation safety is principally studied from perspective of road skid resistance, and road surface topography is considered to be significantly associated with skid resistance of road surface (Kane et al. 2013; Wang et al. 2014) [1, 2]. Since the shape, size, and distribution characteristics of asphalt pavement texture substantially determine the skid resistance of road surface, rational texture feature parameters help to accurately predict the skid resistance of road surface. For the in-depth research of the effect of road surface topography on skid resistance, the PIARC employs the undulating longitudinal wave length of road to describe the road surface topography, having definitely proposed four types of textures, that is, microtexture, macrotexture, megatexture, and roughness. Meanwhile PIARC presented a variety of physical phenomena that correspond to the interaction between running vehicle and road surface, which played a significant role in studying the skid resistance mechanism of road surface, noise, and vehicle safety (PIARC, 1996) [3]. According to a series of researches of PIARC, it is sure that microtexture refers to the morphology with horizontal wavelength of less than 0.5 mm and vertical amplitude ranging from 0 to 0.2 mm, which principally characterizes the surface asperity of aggregate particles and has effect on the actual contact area between tire/road surfaces. Microtexture determines the basic frictional properties of road surface and mainly affects the pavement skid resistance at low speed. Referring to the morphology with horizontal wavelength ranging from 0.5 to 50 mm and vertical amplitude ranging from 0.2 to 10 mm, macrotexture is principally dependent on aggregate shape, size, and distribution and may bring about tire rubber deformation and hysteresis energy loss that result in friction force; macrotexture principally affects the road surface skid resistance in the case of high-speed travel and rainy days (Pulugurtha et al. 2012) [4]. From this point of view, the surface topography of pavement aggregate particle and the distribution of pavement particle protrusion jointly affect the skid resistance of road surface. In recent years, domestic and foreign scholars have carried out a great deal of studies on asphalt pavement morphology information acquisition, characterization, skid resistance evaluation, prediction, and so on and made extremely valuable achievements. Many studies have demonstrated that the friction coefficient is determined by the morphological features between tire and ground, which is also an important factor causing traffic accidents (Kotek and Florková 2014; Qian and Meng 2017) [5, 6]. The UK-based Traffic and Road Research Laboratory (TRRL), as one of the earliest organizations engaged in studies on pavement skid resistance, studied the correlation between the risk of traffic accident on damp road surface and the slippery degree of road surface and developed skid resistance test equipment such as pendulum friction coefficient tester and lateral force coefficient test car SCRIM [7]. Ergun et al. (2005) [8] developed an asphalt pavement microtexture measurement system composed of planar mobile platform, light source system, microscope, CCD camera, image processing system, and so forth of which the horizontal and vertical resolutions are 0.006 mm and 0.01 mm, respectively. Khoudeir et al. (2004) [9] described asphalt pavement’s microtexture characteristics by extracting the statistics (mean value and standard deviation), autocorrelation function deviation, and so on. The gray gradient value of image is analyzed to describe friction property of asphalt pavement, which can be utilized to assess the effect of road surface roughness on pavement wear. Due to the obvious randomness and complexity of the micro and macro morphology of road surface, traditional parameters exhibit instability along with the change in measuring size and range; studies have shown that the morphology illustrates self-similarity and scale independence with the variation of measuring dimension, which means that the road surface topography has indeed fractal characteristics. Kokkalis et al. (2002) [10] proposed a roughness function combining fractal dimension and scale coefficient to describe road surface topography, finding that the skid resistance of asphalt concrete pavement is well correlated with the roughness function value. Zhang et al. (2013) [11] used a laser profiler to test the road surface macroscopic texture profile at 35 test points. It is found that the macroscopic texture profile of asphalt concrete road surface exhibits typical multifractal characteristics, while the multifractal spectrum is rightward hooked and the distribution width of multifractal spectrum reflects the variation range of pavement profile undulation amplitude, which is significantly correlated with the mean profile depth (MPD) of road surface.Researchers in the field of automotive engineering carried out studies on the driving safety, ride comfort, and handling stability of vehicle from the perspective of the mechanical properties of tire, the stress distribution between tire and road surface, and the load transfer. From all perspectives of research it is found that inadequate skid resistance of pavement is a major cause of traffic accident and the road surface topography is significantly associated with vehicle driving performance [12]. The reconstruction of real and accurate road surface roughness has been discussed both in the fields of automotive engineering and road engineering. It has been widely accepted that road surface roughness is a stochastic process [13], which is described as the frequency-domain characteristics of roughness through power spectral density (PSD), power function, and rational function. However, both functions share the same model constitution and provide time-domain model for automobile dynamics study. Domestic and foreign studies on roughness evolved from linear to nonlinear analysis, from frequency-domain to time-domain analysis, and from analytical analysis to simulation calculation. It is particularly that the time-domain model evolved from single point to multipoint model, from sing-track to double-track model, and from 2D to 3D model. Ngwangwa et al. (2014) [14] built a two-dimensional road surface model through artificial neural network, chose some road surface data as training data, and reconstructed the entire road surface using the trained network. Yu et al. (2007) [15] reconstructed three-dimensional road roughness based on multisensor fusion technology. Liu et al. (2014) [16] performed analysis calculation of power spectral density based on the acquired road spectrum data, reconstructed the graded power spectra of various typical road spectra, and analyzed the reconstruction accuracy of road spectrum through correlation coefficient. Road spectrum reconstruction is now principally focused on two-dimensional road surface. However, it fails to take into account the statistical self-similarity of road surface, which results in significant deviation of road spectrum from original road surface spectrum in terms of statistical property. Wullens and Kropp (2004) [17] proposed a three-dimensional contact model for the tyre/road to calculate the dynamic radial contact forces, the local deformation, and normal forced vibrations of the tyre structure. The contact problem is solved using an elastic half-space; the road is assumed as rigid. However, the three-dimensional microroughness of the road surface is also not considered. The fractal theory [18–20] based on fractal interpolation reconstruction method makes up for the shortcomings of traditional refactoring method in respect of road surface roughness self-similarity and employs limited elevation data to reconstruct three-dimensional road surface spectrum similar to the original road surface spectrum while satisfactorily preserving the statistical property and fine structure of road surface data. Wang et al. (2016) [21] built a three-dimensional stochastic road surface model through harmonic superposition, and it was found perfectly consistent with measured road surface according to the comparison in respect of power spectral density. Lu et al. (2014) [22] built a theoretical model of three-dimensional road surface using iterative function method, thereby demonstrating the consistency of reconstructed road surface spectrum with original spectrum. Kogut and Jackson (2006) [23] made a comparison between contact mechanics results obtained with statistical and fractal approaches to characterize surface topography. It is found that differences in the simulated contact area and load can be related solely to the different approach employed for surface characterization. Jiang et al. (2009) [24] proposed that contact stiffness model using fractal geometry topography description is to research the contacts between rough surfaces of machined plane joints. Barnsley (1993) [18] obtained the normal contact stiffness based on fractal theory through a Weierstrass–Mandelbrot function. The results of the theoretical contact stiffness according to the fractal method were also certified by the experimental data. Buczkowski et al. (2014) [25] proposed a modified contact model of fractal rough surface. They concluded that the contact area depends on the contact loading and the contact stiffness increases with increasing contact loading. However, the fractal technique is rarely used to study the problem of contact between tire and ground.In view of the complexity of tire and road surface materials and surface contact characteristics, friction coefficient is far from enough to interpret the tire-pavement contact mechanism. Furthermore, road surface properties are more important for road transportation safety in the event of sudden change. The previous researches on vehicle safety neglected the effect of road surface morphology, and the inadequate description of road surface morphology details in road surface spectra restricts the studies on vehicle-road interaction. For the statistical self-similarity of road surface, finer reconstruction of three-dimensional road surface spectrum could be accomplished through fractal theory. In contrast with two-dimensional road surface spectrum, the three-dimensional road surface spectrum is closer to real road surface. Therefore, this paper employs the fundamental principle [26] of fractal Brownian motion (FBM) and uses the standard deviation of randomly excited roughness of road surface with fractal dimension. The three-dimensional road surface spectrum is reconstructed by diamond-square algorithm and fractal Brownian motion theory. Furthermore, a two-DOF 1/4 car model is built and the tire contact with road surface is modeled as a multipoint-and-plane contact relationship. The coupling between three-dimensional road surface spectrum and three-dimensional tire is analyzed, while a comparative analysis is also done between multipoint-and-plane contact and point contact between tire and rough road surface. It provides a theoretical basis for studies on the dynamic process of interaction between vehicle and road surface, vehicle ride comfort, and road friendliness.


## 2. Time-Domain Model of Random Road Surface Roughness

### 2.1. Power Spectra at Spatial Frequency and Time Frequency

The fitting expression of road surface power spectral densityGq(n) is as follows:(1)Gqn=Gqn0nn0-ω,where n is spatial frequency; m-1 is the reciprocal of wavelength λ; n0 is reference spatial frequency, n0=0.1m-1; ω is frequency index, that is, the diagonal slope on log-log coordinate, which determines the power spectral density frequency structure of road surface; Gq(n0) is the road surface power spectral density at reference spatial frequency n0 called road surface roughness coefficient, m3.When a car travels atu through a road surface with a spatial frequency of n, the equivalent temporal frequency (Hz) is(2)f=un.And the relation between power spectrum densities at temporal frequency and spatial frequency is(3)Gqf=1uGqn,where Gqf is the power spectral density at temporal frequency and Gq(n) is the power spectral density at spatial frequency.


### 2.2. Description of Road Surface Spectrum Time-Domain Model

The PSD of road surface spectrum is statistics corresponding to a certain sort of road surface roughness. However, the reconstructed pavement elevation is not the only one for a given road surface PSD. The road function also corresponds with a sample function in the equivalent pavement elevation of the given road surface spectrum at a certain running velocity.Two prior conditions must be satisfied in order to reconstruct the road time-domain model that satisfies given power spectral density based on a known road surface spectrum: ① The road process is a stationary Gaussian stochastic process. ② The road process is ergodic. The basic process of reconstruction is as follows: Abstract the random fluctuation of road surface elevation into white noise that meets certain conditions, and perform fast Fourier inverse transform to achieve the time-domain model of random roughness of road surface through fitting. Many great methods could be employed to generate road surface elevation time-domain model for stationary Gaussian stochastic process. Principal methods include filtered white noise generation method, stochastic sequence generation method, harmonic superposition method, AR (ARMA) method, and fast Fourier inverse transform generation method. The filtered white noise method is frequently used due to its clear physical significance and easy calculation, as well as the immediate determination of road surface model parameters based on road surface power spectrum values and travel speed. The fitting form of the rational function of power spectral density is as follows:(4)Gqn=2αβ2πα2+n2,where α, β are the constants associated with road surface grade.


### 2.3. Filtered White Noise Generation for Stochastic Road Surface

When a car travels at a constant speedu, the time-domain road surface roughness power spectral density is expressed as follows, due to ω=2πf:(5)Gqω=2π2Gqn0n02uω2.When ω→0, Gq(ω)→∞. Therefore, if the lower cut-off angular frequency is taken into account, the actual power spectral density [12] could be expressed as follows: (6)Gqω=2π2Gqn0n02uω2+ω02,where ω0 is lower cut-off angular frequencyEquation (6) could be taken as the response of white noise excited first-order linear system. According to stochastic vibration theory,(7)Gqω=Hω2Sω,where Hω is frequency response function; Sω is white noise; Wt is power spectral density that takes Sω=1, so(8)Hω=2πn0Gqn0nω0+jω.It can be obtained as (9)q˙t=-2πn0′uqt+2πn0Gqn0uWt,where n0′ is lower cut-off spatial frequency, n0′=0.011m-1; Gq(n0) is road surface roughness coefficient, m3; Wt is white noise Gauss with mean value equal to zero; qt is random elevation displacement of road surface, m.


## 2.1. Power Spectra at Spatial Frequency and Time Frequency

The fitting expression of road surface power spectral densityGq(n) is as follows:(1)Gqn=Gqn0nn0-ω,where n is spatial frequency; m-1 is the reciprocal of wavelength λ; n0 is reference spatial frequency, n0=0.1m-1; ω is frequency index, that is, the diagonal slope on log-log coordinate, which determines the power spectral density frequency structure of road surface; Gq(n0) is the road surface power spectral density at reference spatial frequency n0 called road surface roughness coefficient, m3.When a car travels atu through a road surface with a spatial frequency of n, the equivalent temporal frequency (Hz) is(2)f=un.And the relation between power spectrum densities at temporal frequency and spatial frequency is(3)Gqf=1uGqn,where Gqf is the power spectral density at temporal frequency and Gq(n) is the power spectral density at spatial frequency.


## 2.2. Description of Road Surface Spectrum Time-Domain Model

The PSD of road surface spectrum is statistics corresponding to a certain sort of road surface roughness. However, the reconstructed pavement elevation is not the only one for a given road surface PSD. The road function also corresponds with a sample function in the equivalent pavement elevation of the given road surface spectrum at a certain running velocity.Two prior conditions must be satisfied in order to reconstruct the road time-domain model that satisfies given power spectral density based on a known road surface spectrum: ① The road process is a stationary Gaussian stochastic process. ② The road process is ergodic. The basic process of reconstruction is as follows: Abstract the random fluctuation of road surface elevation into white noise that meets certain conditions, and perform fast Fourier inverse transform to achieve the time-domain model of random roughness of road surface through fitting. Many great methods could be employed to generate road surface elevation time-domain model for stationary Gaussian stochastic process. Principal methods include filtered white noise generation method, stochastic sequence generation method, harmonic superposition method, AR (ARMA) method, and fast Fourier inverse transform generation method. The filtered white noise method is frequently used due to its clear physical significance and easy calculation, as well as the immediate determination of road surface model parameters based on road surface power spectrum values and travel speed. The fitting form of the rational function of power spectral density is as follows:(4)Gqn=2αβ2πα2+n2,where α, β are the constants associated with road surface grade.


## 2.3. Filtered White Noise Generation for Stochastic Road Surface

When a car travels at a constant speedu, the time-domain road surface roughness power spectral density is expressed as follows, due to ω=2πf:(5)Gqω=2π2Gqn0n02uω2.When ω→0, Gq(ω)→∞. Therefore, if the lower cut-off angular frequency is taken into account, the actual power spectral density [12] could be expressed as follows: (6)Gqω=2π2Gqn0n02uω2+ω02,where ω0 is lower cut-off angular frequencyEquation (6) could be taken as the response of white noise excited first-order linear system. According to stochastic vibration theory,(7)Gqω=Hω2Sω,where Hω is frequency response function; Sω is white noise; Wt is power spectral density that takes Sω=1, so(8)Hω=2πn0Gqn0nω0+jω.It can be obtained as (9)q˙t=-2πn0′uqt+2πn0Gqn0uWt,where n0′ is lower cut-off spatial frequency, n0′=0.011m-1; Gq(n0) is road surface roughness coefficient, m3; Wt is white noise Gauss with mean value equal to zero; qt is random elevation displacement of road surface, m.


## 3. Box-Counting Dimensions Method-Based Road Surface Parameter Extraction

The author calculates the fractal dimension of 5,000 m road surface random elevation and makes statistical and fractal analysis of road surface elevation data through box-counting dimensions method-based programming. Fractal dimension could be defined as(10)Lnr=Mnrr=r1-D,where r is measurement scale; Ln(r) is curve length measured with the nth scale; Mn(r) is number of measurement scales; D is fractal dimension of measurement curve.Box-counting dimensions method calculates fractal dimension through covering fractal curve with small boxes with side length ofr. Some boxes are empty, while some cover part of the curve. The number of boxes is counted when the box is not empty. And then the number of nonempty boxes is regarded as Nr. The followed thing is that the size of box gradually reduces, while the Nr will immediately rise. Then the formula is obtained when r→0:(11)D=-limr→0⁡lg⁡Nrlg⁡r.The least square method is used to seek a series ofr and Nr and to fit a straight line in log-log coordinate. The slope of straight line obtained should be the desired fractal dimension. Table 1 shows the fractal dimensions and standard deviations that correspond to road surface spectra at all levels attained from box-counting dimension method.Table 1
Fractal dimension and standard deviation of standard pavement samples.
Road surface grade
A
B
C
D
E
F
G
H
Standard deviation (10−3 mm)
2.2
4.3
8.7
17.4
34.7
69.3
138.7
277.4
Fractal dimension
1.6003
1.6003
1.6003
1.6003
1.6003
1.6003
1.6003
1.6003As shown in Table1, there is no difference under different fractal dimensions with various levels, which indicates obvious self-similarity. The reason is that the inverse Fourier transform assures the coincidence of the straight line of road surface power spectral density obtained through simulation at log-log coordinate, which also reserves the similarity information of road surface roughness.In the case of standard road surface, slight difference is observed between fractal dimensions of road surfaces at various levels obtained by box-counting dimensions method, and the mean value is approximately 1.6, which is taken as the fractal feature index for standard road surface. Since the standard deviation achieved with box-counting dimensions method falls in the international specified range of standard deviation, it could be taken as a statistical indicator for road surface grading. The above-noted two indicators constitute the grading criteria for reconstruction of three-dimensional standard road surface.


## 4. Construction of Three-Dimensional Road Surface Spectrum with Road Surface Morphology Features

The road surface presents random performance and statistical self-similarity, which could be reconstructed through fractal Brownian motion (FBM). The midpoint displacement method is also known as the random midpoint displacement method, which is the simplest and classical method applied in FBM, especially for describing one-dimensional random process. In addition, the diamond-square algorithm is based on midpoint displacement method, which can produce two-dimensional or three-dimensional topography. It can not only simulate three-dimensional pavement, but also get a higher reconstruction accuracy certified in our article later.


### 4.1. Diamond-Square Algorithm Principle

Diamond-square algorithm was originally proposed by Fourniew, Fussell, and Carpenter and is also known as diamond-quadrangle algorithm [27]:(1)
Initialization: the two-dimensional array is initialized and the same elevation value is assigned to the four angles. The size of each dimension is supposed to be 2’s nth power plus 1 (e.g., 33 × 33, 65 × 65, 129 × 129). Figure1 shows the diamond-square algorithm process of a 5 × 5 array, where the elevation value of the 4 angles in Figure 1(a) is initialized and indicated with black spot. Actually, after five-iteration calculation, the road surfaces is separated by a distance of about 3 cm, which can meet our research requirement.(2)
The “diamond” stage: as shown in Figure1(b), a random value is generated using the four points forming a square at the midpoint of such square, that is, the intersection of two diagonals. The midpoint value is equal to the sum of the mean of four corner points’ values and the said random value. Diamond comes into being when a number of squares exist in grid, in which case the center point is indicated with black spot.(3)
The “square” stage: use the four points forming diamond to generate a random value with the same value range as the previous step at the midpoint of that diamond. Also, this midpoint value is equal to the sum of the mean of four corner points’ values and the said random value as shown in Figure1(c), thereby generating a square.(4)
Iterate the process above for the specified number of times.Figure 1
Diamond-square algorithm.
(a)
(b)
(c)
(d)
(e)


#### 4.1.1. Process of Determining One-Dimensional Stochastic Interpolation

The Brownian motion in two-dimensional plane generates three-dimensional morphology of landscape, which means that coordinatesx and y in the plane bring about VH(x,y) as the surface gradient of position (x,y). The change of gradient that occurs during constant speed travel along the straight line path in the XOY plane is fractal Brownian motion. Assuming the travel distance in plane is Δr (where: Δr2=Δx2+Δy2), the curve gradient variation is given by [28, 29](12)ΔV∝ΔrH.The formula of stochastic interpolation algorithm is shown in (13)ZMP=Za+Zb+Zc+Zd4+D′,where D′ represents Gaussian random function and Z means the elevation value of related point. When H≠1/2, (12) leads to [30](14)varVt2-Vt1=t2-t12Hσ2.In (14), V(t) represents the Gaussian distribution when the mean value is equal to zero, and variance is equal to to σ2.AssumingV(t)~N(0,σ2), D~N(0,Δ2), the stochastic interpolation model is(15)Vt2=12Vt-V0+D,where V(0)=0.SolvingΔ2 is the key to achieve the distribution of stochastic increment. The main process is described as follows: The first iteration, V1/2=1/2V1-V0+D1, according to (14)(16)Δ12=σ222H∙1-22H-2.The variance of one-dimensional stochastic increment after n iterations is(17)Δn2=σ22n2H∙1-22H-2.The one-dimensional stochastic increment submits to V(t)~N(0,σ2), Dn~N0,σ2/(2n)2H1-22H-2.


#### 4.1.2. Process of Determining Two-Dimensional Stochastic Interpolation

First, define the matrix range of each point in two-dimensional space as(18)Z=x,y:Xmin≤x≤Xmax,Ymin≤y≤Ymax,∀m,n∈Z+.Define(19)Δx=Xmax-Xminm,Δy=Ymax-Yminn.Hence, each point in the two-dimensional space could be described as(20)xi=Xmin+iΔx,yi=Ymin+jΔy,i=0,1,…,m-1;j=0,1,…,n-1.These discrete points are subjected to the Gaussian distribution; the expected value is 0, while the variance is σ2.Then, the four vertexes are initialized as(21)varX0,1n-X0,0=σ2n2H,varX1m,0-X0,0=σ2m2H,varX1m,1n-X0,0=1m2-1n2Hσ2,where a midpoint X(1/2m,1/2n) is needed; the initial iterative interpolation is expressed as(22)X12m,12n-X0,0=14X0,0+X1m,0+X0,1n+X1m,1n+D1,var12m,12n-X0,0=14H1m2,1n2Hσ2.The first iteration leads to(23)Δ12=σ222H1m2+1n2H-22H-41m2H+1n2H+1m2+1n2H,and the variance of two-dimensional stochastic increment after i iterations is(24)Δi2=σ222iH1m2+1n2H-22H-41m2H+1n2H+1m2+1n2H.Then the stochastic increment of two-dimensional is subjected to X(x,y)~N(0,σ2),Di~N(0,Δi2).


### 4.2. Fractal Brownian Motion Theory-Based Determination of Random Displacement

The random increment in diamond-quadrangle algorithm could be derived from fractal Brownian motion; that is,(25)Δi∝diH,where H is Hurst index, H=2-D2d, D2d=1.6003. di is the segment spacing after i segmentations. The following formula is normally employed for numerical calculation of random displacement:(26)Δi=scale×2-i-H/2×Gauss,or(27)Δi=diH×σ×1-22H-2×Gauss,where “scale” means the scale factor and normally falls within (0,1): Gauss is a random function submitting to standard normal distribution N(0,1); H represents the fractal parameter value of chosen regional terrain; i means the number of iterations for stochastic midpoint displacement; di stands for the segment spacing after i segmentations.The random increment of the first iteration is(28)Δ1=diH1-22H-2σGauss.Accordingly, the random increment after n iterations is (29)Δn=dnH1-22H-2σGauss.In this paper, five iterations are performed (i.e.,n=5), and the sampling interval is 1. The number of iterations is(30)Δ5=1-22H-2σGauss,where H represents fractal parameter and dn means the segment spacing after n iterations. Then the successive number of iterations is(31)zn=z1+z22+Δ5.The number of iterationsi is determined by researchers to satisfy research requirements. Figure 2 shows the flow chart of calculating the elevation of each point in three-dimensional road surface spectrum through fractal Brownian motion and diamond-square algorithm on the premise of five iterations.Figure 2
Algorithm flow chart.


## 4.1. Diamond-Square Algorithm Principle

Diamond-square algorithm was originally proposed by Fourniew, Fussell, and Carpenter and is also known as diamond-quadrangle algorithm [27]:(1)
Initialization: the two-dimensional array is initialized and the same elevation value is assigned to the four angles. The size of each dimension is supposed to be 2’s nth power plus 1 (e.g., 33 × 33, 65 × 65, 129 × 129). Figure1 shows the diamond-square algorithm process of a 5 × 5 array, where the elevation value of the 4 angles in Figure 1(a) is initialized and indicated with black spot. Actually, after five-iteration calculation, the road surfaces is separated by a distance of about 3 cm, which can meet our research requirement.(2)
The “diamond” stage: as shown in Figure1(b), a random value is generated using the four points forming a square at the midpoint of such square, that is, the intersection of two diagonals. The midpoint value is equal to the sum of the mean of four corner points’ values and the said random value. Diamond comes into being when a number of squares exist in grid, in which case the center point is indicated with black spot.(3)
The “square” stage: use the four points forming diamond to generate a random value with the same value range as the previous step at the midpoint of that diamond. Also, this midpoint value is equal to the sum of the mean of four corner points’ values and the said random value as shown in Figure1(c), thereby generating a square.(4)
Iterate the process above for the specified number of times.Figure 1
Diamond-square algorithm.
(a)
(b)
(c)
(d)
(e)


### 4.1.1. Process of Determining One-Dimensional Stochastic Interpolation

The Brownian motion in two-dimensional plane generates three-dimensional morphology of landscape, which means that coordinatesx and y in the plane bring about VH(x,y) as the surface gradient of position (x,y). The change of gradient that occurs during constant speed travel along the straight line path in the XOY plane is fractal Brownian motion. Assuming the travel distance in plane is Δr (where: Δr2=Δx2+Δy2), the curve gradient variation is given by [28, 29](12)ΔV∝ΔrH.The formula of stochastic interpolation algorithm is shown in (13)ZMP=Za+Zb+Zc+Zd4+D′,where D′ represents Gaussian random function and Z means the elevation value of related point. When H≠1/2, (12) leads to [30](14)varVt2-Vt1=t2-t12Hσ2.In (14), V(t) represents the Gaussian distribution when the mean value is equal to zero, and variance is equal to to σ2.AssumingV(t)~N(0,σ2), D~N(0,Δ2), the stochastic interpolation model is(15)Vt2=12Vt-V0+D,where V(0)=0.SolvingΔ2 is the key to achieve the distribution of stochastic increment. The main process is described as follows: The first iteration, V1/2=1/2V1-V0+D1, according to (14)(16)Δ12=σ222H∙1-22H-2.The variance of one-dimensional stochastic increment after n iterations is(17)Δn2=σ22n2H∙1-22H-2.The one-dimensional stochastic increment submits to V(t)~N(0,σ2), Dn~N0,σ2/(2n)2H1-22H-2.


### 4.1.2. Process of Determining Two-Dimensional Stochastic Interpolation

First, define the matrix range of each point in two-dimensional space as(18)Z=x,y:Xmin≤x≤Xmax,Ymin≤y≤Ymax,∀m,n∈Z+.Define(19)Δx=Xmax-Xminm,Δy=Ymax-Yminn.Hence, each point in the two-dimensional space could be described as(20)xi=Xmin+iΔx,yi=Ymin+jΔy,i=0,1,…,m-1;j=0,1,…,n-1.These discrete points are subjected to the Gaussian distribution; the expected value is 0, while the variance is σ2.Then, the four vertexes are initialized as(21)varX0,1n-X0,0=σ2n2H,varX1m,0-X0,0=σ2m2H,varX1m,1n-X0,0=1m2-1n2Hσ2,where a midpoint X(1/2m,1/2n) is needed; the initial iterative interpolation is expressed as(22)X12m,12n-X0,0=14X0,0+X1m,0+X0,1n+X1m,1n+D1,var12m,12n-X0,0=14H1m2,1n2Hσ2.The first iteration leads to(23)Δ12=σ222H1m2+1n2H-22H-41m2H+1n2H+1m2+1n2H,and the variance of two-dimensional stochastic increment after i iterations is(24)Δi2=σ222iH1m2+1n2H-22H-41m2H+1n2H+1m2+1n2H.Then the stochastic increment of two-dimensional is subjected to X(x,y)~N(0,σ2),Di~N(0,Δi2).


## 4.1.1. Process of Determining One-Dimensional Stochastic Interpolation

The Brownian motion in two-dimensional plane generates three-dimensional morphology of landscape, which means that coordinatesx and y in the plane bring about VH(x,y) as the surface gradient of position (x,y). The change of gradient that occurs during constant speed travel along the straight line path in the XOY plane is fractal Brownian motion. Assuming the travel distance in plane is Δr (where: Δr2=Δx2+Δy2), the curve gradient variation is given by [28, 29](12)ΔV∝ΔrH.The formula of stochastic interpolation algorithm is shown in (13)ZMP=Za+Zb+Zc+Zd4+D′,where D′ represents Gaussian random function and Z means the elevation value of related point. When H≠1/2, (12) leads to [30](14)varVt2-Vt1=t2-t12Hσ2.In (14), V(t) represents the Gaussian distribution when the mean value is equal to zero, and variance is equal to to σ2.AssumingV(t)~N(0,σ2), D~N(0,Δ2), the stochastic interpolation model is(15)Vt2=12Vt-V0+D,where V(0)=0.SolvingΔ2 is the key to achieve the distribution of stochastic increment. The main process is described as follows: The first iteration, V1/2=1/2V1-V0+D1, according to (14)(16)Δ12=σ222H∙1-22H-2.The variance of one-dimensional stochastic increment after n iterations is(17)Δn2=σ22n2H∙1-22H-2.The one-dimensional stochastic increment submits to V(t)~N(0,σ2), Dn~N0,σ2/(2n)2H1-22H-2.


## 4.1.2. Process of Determining Two-Dimensional Stochastic Interpolation

First, define the matrix range of each point in two-dimensional space as(18)Z=x,y:Xmin≤x≤Xmax,Ymin≤y≤Ymax,∀m,n∈Z+.Define(19)Δx=Xmax-Xminm,Δy=Ymax-Yminn.Hence, each point in the two-dimensional space could be described as(20)xi=Xmin+iΔx,yi=Ymin+jΔy,i=0,1,…,m-1;j=0,1,…,n-1.These discrete points are subjected to the Gaussian distribution; the expected value is 0, while the variance is σ2.Then, the four vertexes are initialized as(21)varX0,1n-X0,0=σ2n2H,varX1m,0-X0,0=σ2m2H,varX1m,1n-X0,0=1m2-1n2Hσ2,where a midpoint X(1/2m,1/2n) is needed; the initial iterative interpolation is expressed as(22)X12m,12n-X0,0=14X0,0+X1m,0+X0,1n+X1m,1n+D1,var12m,12n-X0,0=14H1m2,1n2Hσ2.The first iteration leads to(23)Δ12=σ222H1m2+1n2H-22H-41m2H+1n2H+1m2+1n2H,and the variance of two-dimensional stochastic increment after i iterations is(24)Δi2=σ222iH1m2+1n2H-22H-41m2H+1n2H+1m2+1n2H.Then the stochastic increment of two-dimensional is subjected to X(x,y)~N(0,σ2),Di~N(0,Δi2).


## 4.2. Fractal Brownian Motion Theory-Based Determination of Random Displacement

The random increment in diamond-quadrangle algorithm could be derived from fractal Brownian motion; that is,(25)Δi∝diH,where H is Hurst index, H=2-D2d, D2d=1.6003. di is the segment spacing after i segmentations. The following formula is normally employed for numerical calculation of random displacement:(26)Δi=scale×2-i-H/2×Gauss,or(27)Δi=diH×σ×1-22H-2×Gauss,where “scale” means the scale factor and normally falls within (0,1): Gauss is a random function submitting to standard normal distribution N(0,1); H represents the fractal parameter value of chosen regional terrain; i means the number of iterations for stochastic midpoint displacement; di stands for the segment spacing after i segmentations.The random increment of the first iteration is(28)Δ1=diH1-22H-2σGauss.Accordingly, the random increment after n iterations is (29)Δn=dnH1-22H-2σGauss.In this paper, five iterations are performed (i.e.,n=5), and the sampling interval is 1. The number of iterations is(30)Δ5=1-22H-2σGauss,where H represents fractal parameter and dn means the segment spacing after n iterations. Then the successive number of iterations is(31)zn=z1+z22+Δ5.The number of iterationsi is determined by researchers to satisfy research requirements. Figure 2 shows the flow chart of calculating the elevation of each point in three-dimensional road surface spectrum through fractal Brownian motion and diamond-square algorithm on the premise of five iterations.Figure 2
Algorithm flow chart.


## 5. Reconstruction and Analysis of Three-Dimensional Rough Road Surface Spectrum

### 5.1. Reconstruction of Three-Dimensional Rough Road Surface Spectrum

The foregoing theory is employed to simulate grade A~H standard highways; Figure3 shows the simulation result after eight-time iterations. The road surface is 600 m in length and 33 m in width. The value is taken from random sample obtained through inverse Fourier transform of standard road surface power spectrum; to simplify the analysis, the X direction and Y direction of road surface and the line segments in any xoy plane share the same fractal characteristics, thus having the same Hurst index in simulation.Figure 3
Reconstruction of three-dimensional rough road.
(a)
A grade pavement
(b)
B grade pavement
(c)
C grade pavement
(d)
D grade pavement
(e)
E grade pavement
(f)
F grade pavement
(g)
G grade pavement
(h)
H grade pavement


### 5.2. Three-Dimensional Rough Road Surface Spectrum Reconstruction Accuracy Analysis

The three-dimensional road surface with 5000 m is reconstructed through fractal interpolation theory and the solution method for power spectral density. Compared with the power spectrum of two-dimensional road surface, the accuracy of road reconstruction proposed in this paper is verified. Due to space limitations, Figure4 only list the comparison diagrams of standard power spectra for grades A~D between theory and simulation results.Figure 4
Power spectrum comparison of road roughness.
(a)
A grade pavement
(b)
B grade pavement
(c)
C grade pavement
(d)
D grade pavementAs shown in Figure4, road surface roughness sequence consistent with statistical property of standardized grade road surface could be achieved by reconstructing road surface roughness using fractal curve. Actually, power spectrum comparison diagram is inadequate for reorganization model accuracy analysis. Since the quantitative analysis is also needed to perform mathematical statistics-based regression analysis by taking road surface roughness power spectrum as regression model. It is necessary to check whether the data of regression model could perfectly fit standard or target data. The power spectral density curve of standard-grade road surface is a straight line, so the regression for simulation of power spectral linear regression analysis should be linear regression, which normally employs residual sum of squares (Se), residual standard deviation (σ^), or correlation index (R2). Any of these three variables could determine the quality of regression equation: Se, σ^, the smaller, the better; R2, the larger, the better, where(32)Se=∑i=1nyi-y^i2,σ^=Sen-2=1n-2∑i=1nyi-y^i2,R2=1-∑i=1nyi-y^i2∑i=1nyi-y¯20<R2<1,where y1,y2,…,yn represents the reconstructed road surface roughness data; y^1,y^2,…,y^n means the original road surface roughness data; n is the number of sampling points; y¯ means the mean of reconstructed road surface roughness data. The reconstructing accuracy for the grade A~D road surface is obtained by the square sum of residuals and correlation coefficient.As shown in Table2, the power spectrum curve attained through the reconstruction of road surface roughness with ideal fractal curve is quite close to standard power spectrum, exhibiting extremely high reconstruction accuracy.Table 2
Regression analysis of standard grade road.
Grade A pavement
Grade B pavement
Grade C pavement
Grade D pavement
Se
0.1542
0.1461
0.5607
1.9363
R2
0.7452
0.9555
0.8751
0.6191


### 5.3. Reconstruction of Three-Dimensional Road Surface Spectrum with Special Road Surface Features

Upon the establishment of three-dimensional stochastic road surface spectrum for each grade of road surface that reflects the road surface topography, some special road surface model can be also simulated to form the more authentic three-dimensional road model. According to transport industry standard JT/T713-2008 “Pavement Rubber Bump” [31], the profile of rubber bump should be approximately trapezoidal, and the bottom width and height of bump should be 300~400 mm and 30~60 mm, respectively. However, no complete standard has been issued for pavement contour curve in terms of specific cross section. The common bump cross-sectional profiles include trapezoid, circular arc, and parabola. Figure 5 shows the parabola-shaped profile.Figure 5
Parabolic cross-sectional profile curve of road bump.The parabolic cross-sectional profile shown in Figure5 is expressed as a mathematical equation:(33)z=H21-cos⁡2πLx+L2.The road bump model file available for simulation is achieved by preparing road bump model file generating program with MATLAB software, developing nodal coordinate matrix through (33), and generating element matrix with MATLAB. Build the parabolic road bump model based on the assumption that the road bump is 300 mm in width and 50 mm in height; identify the installation location of road bump based on the generated three-dimensional road surface; remove the corresponding road surface nodes and cell data from road surface spectrum file based on bump width, and add the node and cell data of bump to appropriate position; realize the isomorphism between road bump and stochastic road surfaces at all levels through newly generated road surface spectrum files. Assuming that the road surface length is 12 m, Figure 6 only presents the road surface model with isomorphism between bump and grade-C road surface.Figure 6
Bump and grade-C random road surface isomorphism simulation result.


## 5.1. Reconstruction of Three-Dimensional Rough Road Surface Spectrum

The foregoing theory is employed to simulate grade A~H standard highways; Figure3 shows the simulation result after eight-time iterations. The road surface is 600 m in length and 33 m in width. The value is taken from random sample obtained through inverse Fourier transform of standard road surface power spectrum; to simplify the analysis, the X direction and Y direction of road surface and the line segments in any xoy plane share the same fractal characteristics, thus having the same Hurst index in simulation.Figure 3
Reconstruction of three-dimensional rough road.
(a)
A grade pavement
(b)
B grade pavement
(c)
C grade pavement
(d)
D grade pavement
(e)
E grade pavement
(f)
F grade pavement
(g)
G grade pavement
(h)
H grade pavement


## 5.2. Three-Dimensional Rough Road Surface Spectrum Reconstruction Accuracy Analysis

The three-dimensional road surface with 5000 m is reconstructed through fractal interpolation theory and the solution method for power spectral density. Compared with the power spectrum of two-dimensional road surface, the accuracy of road reconstruction proposed in this paper is verified. Due to space limitations, Figure4 only list the comparison diagrams of standard power spectra for grades A~D between theory and simulation results.Figure 4
Power spectrum comparison of road roughness.
(a)
A grade pavement
(b)
B grade pavement
(c)
C grade pavement
(d)
D grade pavementAs shown in Figure4, road surface roughness sequence consistent with statistical property of standardized grade road surface could be achieved by reconstructing road surface roughness using fractal curve. Actually, power spectrum comparison diagram is inadequate for reorganization model accuracy analysis. Since the quantitative analysis is also needed to perform mathematical statistics-based regression analysis by taking road surface roughness power spectrum as regression model. It is necessary to check whether the data of regression model could perfectly fit standard or target data. The power spectral density curve of standard-grade road surface is a straight line, so the regression for simulation of power spectral linear regression analysis should be linear regression, which normally employs residual sum of squares (Se), residual standard deviation (σ^), or correlation index (R2). Any of these three variables could determine the quality of regression equation: Se, σ^, the smaller, the better; R2, the larger, the better, where(32)Se=∑i=1nyi-y^i2,σ^=Sen-2=1n-2∑i=1nyi-y^i2,R2=1-∑i=1nyi-y^i2∑i=1nyi-y¯20<R2<1,where y1,y2,…,yn represents the reconstructed road surface roughness data; y^1,y^2,…,y^n means the original road surface roughness data; n is the number of sampling points; y¯ means the mean of reconstructed road surface roughness data. The reconstructing accuracy for the grade A~D road surface is obtained by the square sum of residuals and correlation coefficient.As shown in Table2, the power spectrum curve attained through the reconstruction of road surface roughness with ideal fractal curve is quite close to standard power spectrum, exhibiting extremely high reconstruction accuracy.Table 2
Regression analysis of standard grade road.
Grade A pavement
Grade B pavement
Grade C pavement
Grade D pavement
Se
0.1542
0.1461
0.5607
1.9363
R2
0.7452
0.9555
0.8751
0.6191


## 5.3. Reconstruction of Three-Dimensional Road Surface Spectrum with Special Road Surface Features

Upon the establishment of three-dimensional stochastic road surface spectrum for each grade of road surface that reflects the road surface topography, some special road surface model can be also simulated to form the more authentic three-dimensional road model. According to transport industry standard JT/T713-2008 “Pavement Rubber Bump” [31], the profile of rubber bump should be approximately trapezoidal, and the bottom width and height of bump should be 300~400 mm and 30~60 mm, respectively. However, no complete standard has been issued for pavement contour curve in terms of specific cross section. The common bump cross-sectional profiles include trapezoid, circular arc, and parabola. Figure 5 shows the parabola-shaped profile.Figure 5
Parabolic cross-sectional profile curve of road bump.The parabolic cross-sectional profile shown in Figure5 is expressed as a mathematical equation:(33)z=H21-cos⁡2πLx+L2.The road bump model file available for simulation is achieved by preparing road bump model file generating program with MATLAB software, developing nodal coordinate matrix through (33), and generating element matrix with MATLAB. Build the parabolic road bump model based on the assumption that the road bump is 300 mm in width and 50 mm in height; identify the installation location of road bump based on the generated three-dimensional road surface; remove the corresponding road surface nodes and cell data from road surface spectrum file based on bump width, and add the node and cell data of bump to appropriate position; realize the isomorphism between road bump and stochastic road surfaces at all levels through newly generated road surface spectrum files. Assuming that the road surface length is 12 m, Figure 6 only presents the road surface model with isomorphism between bump and grade-C road surface.Figure 6
Bump and grade-C random road surface isomorphism simulation result.


## 6. Model of Tire Contact with Three-Dimensional Rough Road Surface

The interaction between vehicle and road surface is an extremely complicated dynamic process that involves vehicle dynamics, pavement structural mechanics, and frictional mechanics. When a car travels on road surface, the roughness of road surface is transferred as displacement excitation via tires and suspension to car body and results in the random vibration of car body. The present section builds a 1/4 car body model and further introduces the process of tire contact with random road surface; the load characteristics on tire and tire footprint are integrated through test; the contact surface is considered to be composed of a finite number of points; the contact model is shown in Figure7.Figure 7
Three-dimensional contact model of tire and rough road surface.


### 6.1. Determination of Tire/Road Contact Area

Different pressure is applied to a heavy tire through self-developed actuator system. The test object is heavy-duty radial tire 10.00R20, the standard tire pressure is 830 kPa. The testing site is shown in Figure8; the tire contact distribution is shown in Figure 9.Figure 8
Testing for tire footprint.Figure 9
Tire contact distribution.The contact area is 45,955 mm2 at a standard load of 30,000 N. Test results show that the width of footprint is 200 mm and the contact length is 230 mm. Since the distance between two adjacent points of three-dimensional road surface is 7.8125 mm, there should be 725 contact points under standard tire pressure and load, including 25 points in vertical direction and 29 points in horizontal direction.


### 6.2. Two-DOF Vehicle Model

A Two-DOF 1/4 vehicle suspension model could be expressed by a spring and damper connected in parallel, while the tire could be expressed with a mass block and spring as shown in Figure10, where mt is the unsprung mass, including rim, tire, and axle; ms is the sprung mass, including compartment and load; Kt is tire stiffness coefficient; Ks is suspension system stiffness coefficient; Ct is tire damping constant; Cs is damping constant of shock absorber in suspension system; ζ(t) is ground elevation (road surface roughness), a stochastic process; ys(t) is vertical displacement of sprung part; yt(t) is vertical displacement of unsprung part.Figure 10
1/4 vehicle suspension model.The system motion equation established based on Newton’s Second Law of Motion is(34)mty¨t+Csy˙t-y˙s+Ksyt-ys+Ktyt-ζ=0,msy¨s+Csy˙s-y˙t+Ksys-yt=0,Or(35)MY¨+CY˙+KY=ft,where M=ms00mt represents mass matrix; C=Cs-Cs-CsCs+Ct means damping matrix; K=Ks-Ks-KsKs+Kt is stiffness matrix; Y=[ys,yt]T stands for displacement matrix; f(t) is an excitation matrix; M, C, and K are mass, damping, and stiffness matrixes, respectively, as well as the real symmetric matrixes of n×n, while M is positive definite.


### 6.3. Analysis of Interaction between Tire and Three-Dimensional Rough Road Surface

The subject studied is a heavy-duty automobile of a certain model with tire model 10.00R20; assuming the length and width of tire footprint area remain constant under a standard tire pressure of 830 kPa, the contact area could be determined through test. Select heavy-duty automobile parameters:ms=10109kg, mt=190kg, Ks=75000N/m, Kt=2060000N/m, Cs=30000N/m, Ct=900 N·s/m; tire specifications 10.00R20 R0=526mm, v=20m/s.Take the tire-road surface contact model as a model where there are a finite number of contact points, all of which have the same stiffness; calculate the vehicle-road coupling system under random excitation of road surface using the multipoint-and-plane contact of three-dimensional road surface spectrum and the two-dimensional curve of any profile section, including the vertical acceleration of car body, suspension distortion, and tire force as shown in Figures11 and 12, where the dotted line represents multipoint-and-plane contact model, while the real line represents single point contact model.Figure 11
Vehicle-road response under random three-dimensional road surface spectrum.
(a)
Body vertical acceleration
(b)
Suspension deformation
(c)
Tire forceFigure 12
Vehicle-road response under bump excitation.
(a)
Body vertical acceleration
(b)
Suspension deformation
(c)
Tire force


#### 6.3.1. Response Analysis under Random Three-Dimensional Road Surface Spectrum

Assuming the road surface lengthL0=10m in the case of grade-C random road surface, single point contact model and multipoint-and-plane contact model are used, respectively, to determine the dynamic response of vehicle-road coupling system as shown in Figure 11.As shown in Figure11, the peak values of car body response and tire force of multipoint-and-plane contact model, that is, plane contact model, are much smaller than that of point contact model; the peak value of car body acceleration is smaller by 47%, while the root mean square value is smaller by 57.6%; the peak value of suspension distortion is smaller by 60%, while the root mean square value is smaller by 56%; the peak value of tire force is smaller by 46.7%, while the root mean square value is smaller by 54%. This means the plane contact between tire and ground has buffering and inclusive effect on road surface.


#### 6.3.2. Analysis of Response under Bump Excitation

A parabola-shaped road bump model is built by assuming the road bump width and height to be 300 mm and 50 mm, respectively. As shown in Figure12, single point contact model and multipoint-and-plane contact model are used to achieve the dynamic response of vehicle-road coupling system.According to Figure12, the overall values of car body response and tire force of multipoint-and-plane contact model are smaller than that of point contact model, but the peak value changes slightly when the vehicle is passing through the bump. Furthermore, there is little difference in action time between plane contact model and point contact model.


## 6.1. Determination of Tire/Road Contact Area

Different pressure is applied to a heavy tire through self-developed actuator system. The test object is heavy-duty radial tire 10.00R20, the standard tire pressure is 830 kPa. The testing site is shown in Figure8; the tire contact distribution is shown in Figure 9.Figure 8
Testing for tire footprint.Figure 9
Tire contact distribution.The contact area is 45,955 mm2 at a standard load of 30,000 N. Test results show that the width of footprint is 200 mm and the contact length is 230 mm. Since the distance between two adjacent points of three-dimensional road surface is 7.8125 mm, there should be 725 contact points under standard tire pressure and load, including 25 points in vertical direction and 29 points in horizontal direction.


## 6.2. Two-DOF Vehicle Model

A Two-DOF 1/4 vehicle suspension model could be expressed by a spring and damper connected in parallel, while the tire could be expressed with a mass block and spring as shown in Figure10, where mt is the unsprung mass, including rim, tire, and axle; ms is the sprung mass, including compartment and load; Kt is tire stiffness coefficient; Ks is suspension system stiffness coefficient; Ct is tire damping constant; Cs is damping constant of shock absorber in suspension system; ζ(t) is ground elevation (road surface roughness), a stochastic process; ys(t) is vertical displacement of sprung part; yt(t) is vertical displacement of unsprung part.Figure 10
1/4 vehicle suspension model.The system motion equation established based on Newton’s Second Law of Motion is(34)mty¨t+Csy˙t-y˙s+Ksyt-ys+Ktyt-ζ=0,msy¨s+Csy˙s-y˙t+Ksys-yt=0,Or(35)MY¨+CY˙+KY=ft,where M=ms00mt represents mass matrix; C=Cs-Cs-CsCs+Ct means damping matrix; K=Ks-Ks-KsKs+Kt is stiffness matrix; Y=[ys,yt]T stands for displacement matrix; f(t) is an excitation matrix; M, C, and K are mass, damping, and stiffness matrixes, respectively, as well as the real symmetric matrixes of n×n, while M is positive definite.


## 6.3. Analysis of Interaction between Tire and Three-Dimensional Rough Road Surface

The subject studied is a heavy-duty automobile of a certain model with tire model 10.00R20; assuming the length and width of tire footprint area remain constant under a standard tire pressure of 830 kPa, the contact area could be determined through test. Select heavy-duty automobile parameters:ms=10109kg, mt=190kg, Ks=75000N/m, Kt=2060000N/m, Cs=30000N/m, Ct=900 N·s/m; tire specifications 10.00R20 R0=526mm, v=20m/s.Take the tire-road surface contact model as a model where there are a finite number of contact points, all of which have the same stiffness; calculate the vehicle-road coupling system under random excitation of road surface using the multipoint-and-plane contact of three-dimensional road surface spectrum and the two-dimensional curve of any profile section, including the vertical acceleration of car body, suspension distortion, and tire force as shown in Figures11 and 12, where the dotted line represents multipoint-and-plane contact model, while the real line represents single point contact model.Figure 11
Vehicle-road response under random three-dimensional road surface spectrum.
(a)
Body vertical acceleration
(b)
Suspension deformation
(c)
Tire forceFigure 12
Vehicle-road response under bump excitation.
(a)
Body vertical acceleration
(b)
Suspension deformation
(c)
Tire force


### 6.3.1. Response Analysis under Random Three-Dimensional Road Surface Spectrum

Assuming the road surface lengthL0=10m in the case of grade-C random road surface, single point contact model and multipoint-and-plane contact model are used, respectively, to determine the dynamic response of vehicle-road coupling system as shown in Figure 11.As shown in Figure11, the peak values of car body response and tire force of multipoint-and-plane contact model, that is, plane contact model, are much smaller than that of point contact model; the peak value of car body acceleration is smaller by 47%, while the root mean square value is smaller by 57.6%; the peak value of suspension distortion is smaller by 60%, while the root mean square value is smaller by 56%; the peak value of tire force is smaller by 46.7%, while the root mean square value is smaller by 54%. This means the plane contact between tire and ground has buffering and inclusive effect on road surface.


### 6.3.2. Analysis of Response under Bump Excitation

A parabola-shaped road bump model is built by assuming the road bump width and height to be 300 mm and 50 mm, respectively. As shown in Figure12, single point contact model and multipoint-and-plane contact model are used to achieve the dynamic response of vehicle-road coupling system.According to Figure12, the overall values of car body response and tire force of multipoint-and-plane contact model are smaller than that of point contact model, but the peak value changes slightly when the vehicle is passing through the bump. Furthermore, there is little difference in action time between plane contact model and point contact model.


## 6.3.1. Response Analysis under Random Three-Dimensional Road Surface Spectrum

Assuming the road surface lengthL0=10m in the case of grade-C random road surface, single point contact model and multipoint-and-plane contact model are used, respectively, to determine the dynamic response of vehicle-road coupling system as shown in Figure 11.As shown in Figure11, the peak values of car body response and tire force of multipoint-and-plane contact model, that is, plane contact model, are much smaller than that of point contact model; the peak value of car body acceleration is smaller by 47%, while the root mean square value is smaller by 57.6%; the peak value of suspension distortion is smaller by 60%, while the root mean square value is smaller by 56%; the peak value of tire force is smaller by 46.7%, while the root mean square value is smaller by 54%. This means the plane contact between tire and ground has buffering and inclusive effect on road surface.


## 6.3.2. Analysis of Response under Bump Excitation

A parabola-shaped road bump model is built by assuming the road bump width and height to be 300 mm and 50 mm, respectively. As shown in Figure12, single point contact model and multipoint-and-plane contact model are used to achieve the dynamic response of vehicle-road coupling system.According to Figure12, the overall values of car body response and tire force of multipoint-and-plane contact model are smaller than that of point contact model, but the peak value changes slightly when the vehicle is passing through the bump. Furthermore, there is little difference in action time between plane contact model and point contact model.


## 7. Conclusion

(
1
) There is no difference in fractal dimensions at various levels obtained with box-counting dimensions method, which show obvious fractal characteristics of road surface irregularity.(
2
) The multiscale characterization method is proposed to reconstruct standard-grade road surface morphology through FBM. The contrast analysis between reconstructed spectrum of three-dimensional road surface and theoretical three-dimensional spectrum of two-dimensional road surface in terms of power spectral density demonstrates the higher reconstruction accuracy.(
3
) Two-DOF vehicle model is taken as example to identify the difference between three-dimensional tire-rough road surface multipoint-and-plane contact model and traditional point contact model in terms of the response to car body acceleration, suspension distortion, and tire force. It is included that the multipoint-and-plane contact model reflects the real vertically inclusive characteristics of tire, especially under the excitation of bump.This paper does not take into account such special conditions of road surface as turning, uphill, and downhill, which are to be subjected to more in-depth studies based on practice in future.



---

*Source: 2535409-2018-02-15.xml*