Issue 
Acta Acust.
Volume 4, Number 2, 2020



Article Number  4  
Number of page(s)  11  
Section  Computational and Numerical Acoustics  
DOI  https://doi.org/10.1051/aacus/2020003  
Published online  16 April 2020 
Scientific Article
Polynomial relations for cylindrical wheel stiffness characterization for use in a rolling noise prediction model
^{1}
Matelys Research Lab, 7 Rue des Maraîchers, Bât B, 69120 VaulxenVelin, France
^{2}
INSA–Lyon, 20 Avenue Albert Einstein, 69621 Villeurbanne cedex, France
^{*} Corresponding author: matthew.edwards@matelys.com
Received:
6
January
2020
Accepted:
10
March
2020
In vehicle tire/road contact modeling, dynamic models are typically used which incorporate the vehicle’s suspension in their estimation: thus relying on a known stiffness to determine the movement of the wheel in response to roughness excitation. For the case of a wheeled device rolling on a floor (such as a delivery trolley moving merchandise around inside a commercial building), there is often no suspension, yet the wheel is still too soft to able to be considered mechanically rigid (as is the case in train/rail contact). A model which is aimed at incorporating the dynamic effects of the trolley in predicting the sound generated by rolling needs to provide a robust way of estimating the wheel’s effective stiffness. This work presents an original technique for estimating the stiffness of a solid cylindrical wheel. A parametric study was conducted in order to identify the dependence of the wheel stiffness on each of the relevant variables: including the wheel’s radius, axle size, width, applied load, and material properties. The methodology may be used to estimate the stiffness of new wheel types (i.e. different geometries and materials) without needing to solve a finite element model each time. Such a methodology has application beyond the field of acoustics, as the characterization of shapes with nonconstant cross sections may be useful in the wider field of materials science.
© M. Edwards et al., Published by EDP Sciences, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
In the field of acoustics, prediction models which estimate the sound produced in a rolling event are widespread. Mostly existing for automotive tire/road and train wheel/rail contact, these models calculate the interaction forces between the two bodies (the wheel and the surface upon which the wheel rolls), which are produced by the smallscale relative roughness in the contact area [1, 2]. These forces induce motion in the two bodies, and their resulting structural vibrations are what radiates sound to the surrounding area.
These models may typically be described as consisting of three main subsystems: a contact model, a dynamic model, and a propagation model. The contact model calculates the contact force between the two bodies as a function of (primarily) the smallscale relative roughness between them. The dynamic model calculates the motion of the two bodies as a function of the previously calculated contact force. Finally, the propagation model calculates the radiated sound as a function of the previously calculated motion of the two bodies.
Recently, work has been done to develop a prediction model applicable to indoor rolling scenarios, such as that of a delivery trolley rolling on a floor inside a commercial space [3]. Such a situation presents unique challenges which are not present in automotive tire/road or train wheel/rail rolling contact. One of these differences is that of the relative stiffnesses of the two contacting bodies. In the dynamic models of train wheel/rail rolling contact, the wheel and rail themselves are generally considered mechanically rigid, due to their metallic construction and thus extremely high stiffnesses [4]. Instead, the train bogie suspension and railroad sleepers are modeled as dynamic systems which respond to the contact force as their input [5–8]. In the dynamic models used in automotive tire/road rolling contact, the road is generally considered mechanically rigid, its vibration contributing minimally to the propagation of sound [9–12]. Calculation of the complex motion of the tire is instead given priority. In a dynamic model used for indoor trolley wheel/floor rolling contact, the wheel and floor cannot always be considered mechanically rigid. They may both move and deform on a macro level in response to the forces in the area of contact between them, due to their often relatively low elasticities.
An indoor rolling trolley may be represented in the dynamic model as some form of a springmassdamper system. For a vehicle with a builtin suspension system, the stiffness is straightforward to calculate: it is simply the stiffness of the suspension spring. However, most trolleys do not have suspension systems: they are simply made up of a wheel turning about a rigid axle. In such a scenario, the stiffness used in the dynamic model is instead the stiffness of the wheel itself. (This is not to be confused with the local contact stiffness between the wheel and the floor, which represents the small scale deformation due to the interpenetration of the two bodies.) Considering the vast range of sizes and materials that trolley wheels come in, this presents a challenge for incorporation into a sound prediction model that wishes to be widely applicable for a range of indoor rolling scenarios.
One option may be to build a finite element (FE) model for a given wheel, but this is costly, as it requires a new FE model to be run for every single unique geometry/material composition. Additionally, a further problem exists in the fact that the geometry of the wheel in a rolling event is dependent on the motion of the wheel itself. According to Hertzian contact theory, the size of the contact area formed between a cylindrical wheel and a flat floor is a function of the wheel radius, wheel width, the compressional load on the wheel, and the elastic properties (Young’s modulus and Poisson’s ratio) of the wheel and floor [13]. In a static scenario, the load is constant, and the size of the contact area may be easily calculated. In a dynamic scenario, the load changes continuously as the motion of the wheel and floor respond to their relative interaction. Consequently, the “flat spot” formed on the bottom of the wheel by the contact area (and thus the effective geometry of the wheel itself) also changes continuously. If FE modeling were to be used in the calculation of each unique wheel’s stiffness, either the dynamic phenomenon of the wheel geometry would need to be neglected, or a FE model would need to be run at each individual time step in the dynamic model. Both options are neither desirable nor practical.
An alternative option would be to run a series of FE models in the form of a parametric study, extract trends in the data, and develop a system for estimating the stiffness of any geometry or material composition, without the need to solve a new FE model each time. In this paper, a straightforward method is proposed to characterize the compressional stiffness of a cylindrical wheel in contact with a flat surface, including how it changes with time throughout the course of the rolling event. This is inspired by the methodology used by Sim and Kim to estimate the elastic properties of viscoelastic materials [14]. Polynomial relations are derived from high order FE models of a cylindrical wheel under static vertical compression, linking the wheel/axle radii ratio and contact area halflength. A series of one dimensional polynomials or a single two dimensional polynomial may be found which completely describe how a cylindrical wheel’s stiffness changes as a function of its changing geometry. Such polynomial functions, or alternatively a classical lookup table may be incorporated into an indoor rolling noise model, allowing the wheel stiffness to be calculated throughout the rolling event. This allows for a more precise estimation of the movement of the wheel, leading to increased accuracy in sound predictions.
2 Parametric study
A parametric study was conducted in order to identify the degree of influence that each parameter has on the wheel stiffness. A series of FE models were run for cylindrical wheels of various geometries and material compositions, from which trends were identified. For the parametric study, four geometric variables, shown in Figure 1, were chosen to be investigated: wheel radius r _{W}, axle radius ${r}_{\mathrm{A}}$, wheel width w, and contact area halflength a. The convention of halflength is adopted from Hertzian contact theory, and simply refers to half the length of the contact area in the direction of rolling (the full contact area having length $2a$). Additionally, Young’s modulus E and Poisson’s ratio ν of the wheel material, as well as the input displacement ${z}_{\mathrm{in}}$, were also included.
Figure 1 Description of the wheel geometry. 
2.1 Meshing
The wheel geometries were constructed using pythonOCC [15]. Meshing was done using Tetgen [16]. The mesh densities of the FE models were first defined based on the width of the wheel, having a maximum edge length of w/7. They were then remeshed to add extra nodes in the area of contact for particularly small contact area halflengths. This ensured the contact area had at least four nodes along its length. Depending on the geometry, mesh sizes ranged between 1.0 × 10^{4} and $5.3\times 1{0}^{4}$ nodes, having degrees of freedom between $2.1\times 1{0}^{6}$ and $1.1\times 1{0}^{7}$.
To solve the FE model, the flat contact surface on the bottom of the wheel was fixed in place. A uniform displacement in the downward vertical direction was imposed on the axle surface, simulating the mass of the trolley borne by the wheel under load. A glued condition between the axle and wheel was assumed for the purpose of FE calculations.
The stiffness of each FE model was calculated by dividing the force on the axle surface by the displacement of the axle surface:
$$K={\left(\frac{{F}_{z}}{{U}_{z}}\right)}_{\mathrm{axle}\mathrm{surface}},$$(1)where the vertical normal stress ${\sigma}_{\mathrm{zz}}$ and vertical displacement ${u}_{z}$ are integrated across the whole axle surface to find it’s overall force and displacement
$${F}_{z}=\int \mathrm{}{\sigma}_{\mathrm{zz}}\mathrm{dA},$$(2)
$${U}_{z}=\frac{1}{A}\int \mathrm{}{u}_{z}\mathrm{dA}.$$(3)
A Windows 10 PC with an AMD Ryden Threadripper 1950X 16Core 3.40 GHz processor and 64 GB of RAM was used for all calculations. The FE models were solved in FreeFem++ [17] with the UMFPACK direct solver.
2.1.1 Investigated variables
For each of the variables chosen to be investigated, a range was selected which encompasses the full scope of any trolley wheel which may be reasonably expected to be found in the real world. Table 1 summarizes the range of values that were chosen for each variable. The axle radius, wheel width, and contact area halflength are expressed as ratios: relative to the wheel radius. This nomenclature allows for ease of comparison due to the use of dimensionless quantities. Due to limitations with the meshing software, relative contact area halflength values below 0.003 were not able to be reliably meshed.
Range of values used for each parameter.
In lieu of running a calculation for every possible permutation of the values of all six parameters (which would have resulted in over four billion individual runs), each parameter was given its own investigation. Thus for a given parameter, all the others were fixed at a constant value, and the stiffness calculated for each value of the parameter in question. For example, for the runs investigating the effect of the wheel width, all other parameters remained constant, and only the wheel width was varied from one run to another. The only exception to this was Poisson’s ratio, which was investigated for all eleven of its values across all parameter investigations. Table 2 shows the constant values chosen for the parametric study. For a given parameter with N possible values, 11N runs were conducted. This resulted in 1342 runs overall.
Values used for the fixed variables.
3 Parametric study results
After all the desired conditions were solved, the results were analyzed, and the calculated stiffness of the wheel plotted against each parameter in question. The results of the FE simulations reveal patterns in how the wheel stiffness is influenced by each of the various parameters: some expected, some not. For each plot, fixed parameters (those which are identical for all data points) are shown in the bottom right hand corner. Points having the same horizontal position are for varying values of Poisson’s ratio (0–0.495). Schematics of the extreme wheel geometries on each plot are shown for reference: drawn to scale with respect to one another.
The parameters having a linear influence on the wheel stiffness, shown in Figure 2, are Young’s modulus and the wheel width. For a given FE model, a change in either of these values results in a proportional change to the wheel stiffness. Since Young’s modulus of a material is essentially a measure of that material’s intrinsic stiffness independent of shape, logic follows that it would be linearly proportional to the actual wheel stiffness. Similarly, increasing the width of the wheel essentially acts like adding springs in parallel, thus increasing the stiffness in a linear fashion.
Figure 2 Results of the parametric study: (a) Young’s modulus $E$ and (b) relative wheel width $w/{r}_{\mathrm{W}}$. Points having the same x coordinate (but a different y coordinate) are for varying values of Poisson’s ratio $\nu $. Schematics of the extreme wheel geometries are shown for reference (drawn to scale with respect to one another). 
The two parameters exhibiting unique relationships with the wheel stiffness, shown in Figure 3, are the axle radius and the contact area halflength. Neither may be described by a simple expression, and will thus become the topic of further investigation in Section 4. For the axle radius: as the ratio between the axle radius and the wheel radius approaches unity, the stiffness approaches infinity. Since the axle is considered rigid in the FE models, as the volume of wheel material approaches zero, the compressed wheel gets closer and closer to exhibiting that of a rigid connection, thus having infinite stiffness. For the contact area halflength: as the size of the contact area approaches zero, the stiffness begins to drop off exponentially. While not able to be directly calculated, it is believed that with even smaller contact areas, the stiffness would continue to fall towards zero at an ever increasing rate.
Figure 3 Results of the parametric study: (a) relative axle radius ${r}_{\mathrm{A}}/{r}_{\mathrm{W}}$ and (b) relative contact area halflength $a/{r}_{\mathrm{W}}$. Points having the same x coordinate (but a different y coordinate) are for varying values of Poisson’s ratio $\nu $. Schematics of the extreme wheel geometries are shown for reference (drawn to scale with respect to one another). 
Looking on the right side of Figure 3a, it may not appear at first glance to be a reasonable geometry for a trolley wheel. After all, an axle that is makes up 96% of the volume of the wheel/axle assembly may hardly be considered realistic. Indeed, this is not the type of scenario being exemplified by these high axle/wheel radii ratios. They are instead exemplifying the scenario of a twopiece construction wheel, where a hard wheel core is wrapped in a softer outer layer. Such a situation is relatively common for indoor trolley wheels. Here the core is sufficiently stiffer than the softer outer layer that it may be considered mechanically rigid. Thus, for the purpose of the FE models, the inner wheel core and axle may be considered oneinthesame. Schematics of these kinds of wheels are shown in Figure 4. Here, wheels (a) and (c) are being exhibited by the FE models, not wheel (b).
Figure 4 Three example wheels: (a) A realistic onepiece wheel with a small solid axle, (b) an unrealistic onepiece wheel with an impractically large solid axle, and (c) a realistic twopiece wheel with a soft outer wheel layer, a hard inner wheel core, and a solid axle. In the scenario of a high wheel/axle ratio, wheel (c) is being exhibited, not wheel (b). 
Figure 5 shows the three parameters which have no influence on the wheel stiffness: the wheel radius and input displacement. The first comes with a caveat, as the ratio between the axle radius and the wheel radius does indeed have a strong effect on the wheel stiffness. However, when this ratio is held constant, increasing the wheel radius does not result in any change in wheel stiffness. So perhaps it is better to say that the wheel radius does have an effect on the stiffness, but that this effect can be fully described by the changing axle/wheel radii ratio. In regards to input displacement, the lack of dependence is in accordance with the use of a linear material in the FE models, where increasing compression does not result in increased stiffness.
Figure 5 Results of the parametric study: (a) wheel radius ${r}_{\mathrm{W}}$ and (b) input displacement ${z}_{\mathrm{in}}$. Points having the same x coordinate (but a different y coordinate) are for varying values of Poisson’s ratio $\nu $. Schematics of the extreme wheel geometries are shown for reference (drawn to scale with respect to one another). 
Finally, there is one more parameter which deserves discussion: Poisson’s ratio. An elastic solid of uniform cylindrical or rectangular prism shape, compressed in the axial direction, is known to have a high dependency on Poisson’s ratio [14]. However, looking across all six plots in Figures 2–5, a change in Poisson’s ratio results in essentially no change in stiffness (all else held constant). The points which exhibit the largest deviation are all for the highest values of Poisson’s ratio (above 0.49), and they do not occur in a consistent fashion (i.e. not always an increase and not always a decrease). In essence, the stiffness is remaining nearly unchanged for all values of Poisson’s ratio all the way up until around 0.49, where only then the modeling instability that exists when approaching ν = 0.5 starts to have a minor effect. This near lack of dependence on Poisson’s ratio was certainly unexpected. In fact, when loaded in the vertical direction it is the effect of the wheel shape which dominates. The cylindrical form allows the wheel material to “move out of the way” under compression into the hornshaped free spaces between the wheel and the floor. A visualization of this is show in Figure 6. The shape of the geometry itself is linking the movement of the two transverse directions to that of the longitudinal direction in a way that negates any influence that the Poisson effect may have.
Figure 6 Diagram of the wheel compression under an arbitrary load $Q$. The movement of the wheel into the areas formed between the wheel and the floor overrides any effect a change in Poisson’s ratio may have on the stiffness. 
4 Generation of the abacus
To generate the abacus, a second series of FE models were run for the two chosen nonlinearly dependent variables: the axle radius and the contact area halflength (both normalized by the wheel radius). Nineteen values between 0.05 and 0.98 were chosen for the axle/wheel radii ratio ${r}_{\mathrm{A}}/{r}_{\mathrm{W}}$. Eighteen values between 0.003 and 0.1 were chosen for the relative contact area halflength $a/{r}_{\mathrm{W}}$. Constant values of $w/{r}_{\mathrm{W}}=0.4$, $E=1$ Pa, ν = 0.3, and ${z}_{\mathrm{in}}=0.1$ mm were chosen for the wheel width, Young’s modulus, and Poisson’s ratio, respectively. Though it should be noted that, due to their relationship with the wheel stiffness (linear for $w$ and $E$, and independent for ν and ${z}_{\mathrm{in}}$), their choice is arbitrary.
4.1 Normalized wheel stiffness
Figure 7 shows the results of the second series of FE calculations. Equation (4) was first used to find the normalized wheel stiffness by dividing by Young’s modulus and the wheel width:
$${K}_{\mathrm{norm}}=\frac{K}{\mathrm{wE}}.$$(4)
Figure 7 Normalized wheel stiffness as a function of axle/wheel radii ratio and relative contact area halflength. Plotted on a logarithmic zscale (base 10). 
The normalized stiffness was then plotted as a function of the axle/wheel radii ratio and relative contact area halflength. The data is plotted with a logarithmic zscale (base 10) in order to visualize how the stiffness changes for both low and high values of ${r}_{\mathrm{A}}/{r}_{\mathrm{W}}$. As the axle/wheel radii ratio grows, so does the influence of the size of the contact area on the stiffness. For ratios 0–0.8, there is a 1.5x factor in the change in normalized stiffness. Above 0.8 the stiffness increases exponentially, up to a factor of nearly 16x.
Three methods have been developed for implementing the data from the abacus into a rolling noise model: a lookup table, a series of onedimensional polynomials, and a single twodimensional polynomial. They are described in detail in the following sections.
4.2 Lookup table
A straightforward method of implementing the wheel stiffness abacus is to use the data points in Figure 7 as a lookup table. The normalized wheel stiffness for a wheel of any reasonable size (i.e. for $0.05\le {r}_{\mathrm{A}}/{r}_{\mathrm{W}}\le 1$ and $0.003\le a/{r}_{\mathrm{W}}\le 0.1$) may be calculated via twodimensional linear interpolation of the data points. This encompasses nearly every type of cylindrical solid wheel which would be expected to be found on an indoor trolley.
4.3 Series of onedimensional polynomials
In lieu of using a lookup table, polynomial relations may instead be generated to characterize the relationship between the influencing parameters and the wheel stiffness. This has the benefit of simplifying the implementation procedure into a rolling noise model. To do so, univariate 9th order exponential polynomials following the form of equation (5) were first fitted to the data in Figure 7 with ${r}_{\mathrm{A}}/{r}_{\mathrm{W}}$ as the dependent variable: giving one polynomial for each relative contact area half length:
$${P}_{a/{r}_{\mathrm{W}}}({r}_{\mathrm{A}}/{r}_{\mathrm{W}})=\frac{K}{\mathrm{wE}}=\mathrm{exp}\left(\sum _{i=1}^{9}\mathrm{}{C}_{i}^{a/{r}_{\mathrm{W}}}({r}_{\mathrm{A}}/{r}_{\mathrm{W}}{)}^{i}\right).$$(5)
This was chosen as it is the highest order polynomial for which a close fit is observed. Equation (5) produces polynomials having coefficient of determination ${R}^{2}>0.998$ for all $a/{r}_{\mathrm{W}}$. Higher order polynomials began to exhibit unstable behavior (large peaks and troughs started to emerge between the known data points) and no longer accurately mapped to the data sets. These polynomials characterize how the normalized stiffness changes with varying axle/wheel radii ratio for a given contact area half length.
A second series of univariate 6th order polynomials were then fitted to data with $a/{r}_{\mathrm{W}}$ as the dependent variable: giving one polynomial for each axle/wheel radii ratio. These polynomials, following the form of equation (6), characterize how the normalized stiffness changes with varying contact area half length for a given axle/wheel radii ratio:
$${P}_{{r}_{\mathrm{A}}/{r}_{\mathrm{W}}}(a/{r}_{\mathrm{W}})=\frac{K}{\mathrm{wE}}=\sum _{i=1}^{6}\mathrm{}{C}_{i}^{{r}_{\mathrm{A}}/{r}_{\mathrm{W}}}(a/{r}_{\mathrm{W}}{)}^{i}.$$(6)
Figure 8 shows the generated onedimensional polynomials for both varying axle/wheel radii ratio and contact area halflength. Together, they may be used to calculate the stiffness of a wheel of nearly any geometry. The procedure for how this is done in practice is shown in Section 5.
Figure 8 (a) Nineth order 1D polynomials for normalized wheel stiffness as a function of axle/wheel radii ratio. (b) Sixth order 1D polynomials for normalized wheel stiffness as a function of contact area halflength. Plotted on a logarithmic yscale (base 10). 
4.4 Twodimensional polynomial
The third characterization method developed involves fitting a single twodimensional polynomial to the entire data set. This greatly simplifies the implementation process by reducing the entire calculation procedure to a single step. A homogeneous bivariate exponential polynomial of the form shown in equation (7) was fit to the data points:
$$P\left(f,g\right)=\frac{K}{\mathrm{wE}}=\mathrm{exp}\left[\sum _{N=0}^{5}\mathrm{}\left(\sum _{i=1}^{N}\mathrm{}{C}_{i,Ni}{f}^{i}{g}^{Ni}\right)\right].$$(7)
For ease of display, the expressions for axle/wheel radii ratio and contact area halflength have been replaced by functions $f$ and $g$, such that
$$f=\mathrm{ln}({r}_{\mathrm{A}}/{r}_{\mathrm{W}}),$$(8)
The resulting 2D polynomial is given by equation (10). Note that the coefficients have been rounded here for ease of display:
$$\begin{array}{ll}P(f,g)=\frac{K}{\mathrm{wE}}=& \mathrm{exp}(95{f}^{5}+280{f}^{4}g1319{f}^{3}{g}^{2}+4004{f}^{2}{g}^{3}9965f{g}^{4}+1052708{g}^{5}\\ & 249{f}^{4}352{f}^{3}g+1109{f}^{2}{g}^{2}+76f{g}^{3}305716{g}^{4}+247{f}^{3}+123{f}^{2}g\\ & 454f{g}^{2}+33013{g}^{3}117{f}^{2}+7\mathrm{fg}1625{g}^{2}+29f+43g5)\end{array}.$$(10)
The 2D polynomial is plotted with the FE model results in Figure 9. It has an extremely good fit, with a coefficient of determination of ${R}^{2}=0.996$. This polynomial fully characterizes how the normalized stiffness changes with nearly any geometry.
Figure 9 Fifth order homogeneous bivariate polynomial for normalized wheel stiffness as a function of both axle/wheel radii ratio and relative contact area halflength. Plotted on a logarithmic zscale (base 10). 
5 Comparison of the three methods
In order to compare the three methods, as well as demonstrate how each may be used to calculate the stiffness of a given wheel, an example characterization of performed. Let us take two wheels with properties given in Table 3.
Example wheel parameters.
For a cylindrical wheel, Hertzian contact theory states the size of the contact area will be
$$a=\sqrt{\frac{4Q{r}_{\mathrm{W}}}{\pi \mathrm{wE}\mathrm{\prime}}},$$(11)where Q is the applied load and $E\mathrm{\prime}$ is the apparent Young’s modulus between the wheel and floor.
$$E\mathrm{\prime}={\left(\frac{1{\nu}_{\mathrm{wheel}}^{2}}{{E}_{\mathrm{wheel}}}+\frac{1{\nu}_{\mathrm{floor}}^{2}}{{E}_{\mathrm{floor}}}\right)}^{1}.$$(12)
For the example, we will use a concrete floor (${E}_{\mathrm{floor}}=3.3$ GPa, ${\nu}_{\mathrm{floor}}=0.2$). If the wheel is put under a static load of $Q=200$ N, Equations (11) and (12) estimate a contact area halflength of ${a}_{1}=0.80$ mm and ${a}_{2}=1.83$ mm for each wheel. Thus the dimensionless geometric quantities necessary to calculate the stiffnesses are ${r}_{\mathrm{A}}/{{r}_{\mathrm{W}}}_{1}=0.30$ and $a/{{r}_{\mathrm{W}}}_{1}=0.02$ for wheel 1 and ${r}_{\mathrm{A}}/{{r}_{\mathrm{W}}}_{2}=0.897$ and $a/{{r}_{\mathrm{W}}}_{2}=0.03$ for wheel 2.
Using the lookup table, the procedure is straightforward. Plugging the dimensionless geometric quantities into the lookup table and interpolating yields normalized wheel stiffnesses of K _{norm,1} = 0.31 and K _{norm,2} = 1.05. Using equation (4), we obtain the true estimated wheel stiffnesses: K _{1} = 5811 kN and K _{2} = 4960 kN.
Figure 10 demonstrates the process of building the 1D polynomials which are specific to our two example wheels. Each of the first series of 9th order 1D polynomials are evaluated for the given axle/wheel radii ratio. This is represented by the vertical line drawn at the given ${r}_{\mathrm{A}}/{r}_{\mathrm{W}}$ values. This provides a new set of data points, to which the second 6th order polynomials can now be fit.
Figure 10 Example of the 1D polynomial characterization procedure for ${r}_{\mathrm{A}}/{{r}_{\mathrm{W}}}_{1}=0.30$ and ${r}_{\mathrm{A}}/{{r}_{\mathrm{W}}}_{2}=0.897$. Plotted on a logarithmic yscale (base 10). 
Equations (13) and (14) show the final 1D polynomials for each example wheel. Note that the coefficients have been rounded here for ease of display:
$$\begin{array}{ll}{P}_{1}(a/{{r}_{\mathrm{W}}}_{1})={\left(\frac{K}{\mathrm{wE}}\right)}_{1}=& 5506860(a/{r}_{\mathrm{W}}{)}^{6}+1766278(a/{r}_{\mathrm{W}}{)}^{5}222310(a/{r}_{\mathrm{W}}{)}^{4}\\ & +14130(a/{r}_{\mathrm{W}}{)}^{3}498(a/{r}_{\mathrm{W}}{)}^{2}+12(a/{r}_{\mathrm{W}})+0.2,\end{array}$$(13)
$$\begin{array}{ll}{P}_{2}(a/{{r}_{\mathrm{W}}}_{2})={\left(\frac{K}{\mathrm{wE}}\right)}_{2}=& 13397078(a/{r}_{\mathrm{W}}{)}^{6}+4314984(a/{r}_{\mathrm{W}}{)}^{5}545898(a/{r}_{\mathrm{W}}{)}^{4}\\ & +34857(a/{r}_{\mathrm{W}}{)}^{3}1177(a/{r}_{\mathrm{W}}{)}^{2}+39(a/{r}_{\mathrm{W}})+0.3.\end{array}$$(14)
Plugging the relative contact area halflengths into the above equations yields normalized wheel stiffnesses of ${K}_{\mathrm{norm},1}$ = 0.31 and K _{norm,2} = 1.03. Using equation (4), we obtain the true estimated wheel stiffnesses: K _{1} = 5954 kN and K _{2} = 4867 kN.
The process for using the 2D polynomial is also very straightforward. All that is needed is to solve equation (7) for the given dimensionless geometric quantities. Doing so yields normalized wheel stiffnesses of K _{norm,1} = 0.30 and K _{norm,2} = 1.09. Using equation (4), we obtain the true estimated wheel stiffnesses: K _{1} = 5793 kN and K _{2} = 5146 kN.
Table 4 summarizes the results of the three methods. For wheel 1, there is a $2.4\%$ difference between the results of the lookup table and 1D polynomials, a 0.3% difference between the lookup table and 2D polynomial, and a 2.7% difference between the 1D polynomials and 2D polynomial. For wheel 2, there is a 1.9% difference between the results of the lookup table and 1D polynomials, a 3.7% difference between the lookup table and 2D polynomial, and a 5.6% difference between the 1D polynomials and 2D polynomial. All three methods provide reasonably close results to one another.
Summary of example results.
The major difference in methods becomes apparent when comparing their calculation times. As an example, a simple MATLAB script was written to perform a loop of 10 000 iterations using each calculation method. Using the same computer that was used for the parametric study (described in Sect. 2.1), the 1D polynomials and the 2D polynomial methods each finished in less than 0.1 s. The lookup table method, however, took more than 9 s to complete. When running a Python version of the same script, the difference is even more pronounced (0.1, 0.1, and 40.4 s for the 1D polynomials, 2D polynomial, and lookup methods, respectively). In a time domain model implementation where the wheel stiffness is recalculated at each moment in time throughout the course of the simulation, implementation of one of the two polynomial methods would provide significant improvement over the lookup method in terms of computational efficiency. This is discussed further in Section 7.1.
6 Comparison with analytical methods
The results of the previous FEbased estimations were compared with two analytical methods to identify whether these simpler analytical solutions could be used in their place as an accurate representation of the wheel stiffness. Two analytical methods were investigated: an equivalent beam and an equivalent trapezoidal prism.
6.1 Equivalent beam
In the rolling noise model presented in [3], the stiffness used in the dynamic model is taken as that of an equivalent beam, who’s crosssection is defined by the contact area, and who’s height is defined by the distance between the floor and the axle (the presence of an axle is ignored in [3], but is included here for comparison). This may be given by
$${K}_{\mathrm{beam}}=\frac{2\mathrm{awE}}{\sqrt{{{r}_{\mathrm{W}}}^{2}{a}^{2}}{r}_{\mathrm{A}}}.$$(15)
This may be rearranged to have the same form as equation (4) in order to achieve the normalized equivalent beam stiffness:
$${K}_{\mathrm{beam},\mathrm{norm}}=\frac{{K}_{\mathrm{beam}}}{\mathrm{wE}}=\frac{2(a/{r}_{\mathrm{W}})}{\sqrt{1(a/{r}_{\mathrm{W}}{)}^{2}}({r}_{\mathrm{A}}/{r}_{\mathrm{W}})}.$$(16)
Figure 11 shows the estimated wheel stiffness found using equations (10) and (16), as well as the percent error between the two as a function of axle/wheel radii ratio and contact area halflength. When compared to the wheel abacus, the approximation given by the equivalent beam is not very accurate. Calculating the equivalent beam stiffnesses for the wheel geometries used in this work yields a dataset with relatively poor fit to the wheel abacus: having an average percent error of 62% (median 67%), with some values as high as 247% for extremely small axle/wheel radii ratios. The standard deviation of the percent error is 10%.
Figure 11 (a) Normalized wheel stiffness: FE model results (blue) vs an equivalent beam (red). Plotted on a logarithmic zscale (base 10). (b) Percent error between the FE model results and equivalent beam. 
6.2 Equivalent trapezoidal prism
The issue with the simple beam is that it has a constant crosssection throughout its height. The cylindrical wheel, with its two load points being the contact patch from below and the axle surface from above, is perhaps more akin to a trapezoidal prism. To that effect, an equivalent prism whose crosssectional area changes along its height, from $2\mathrm{aw}$ on one end to $2{r}_{\mathrm{A}}w$ on the other, may yield a closer approximation of the equivalent stiffness. Such a prism would have a stiffness given by
$${K}_{\mathrm{trap},\mathrm{norm}}=\frac{{K}_{\mathrm{trap}}}{\mathrm{wE}}={\left({\int}_{0}^{h}\mathrm{}\frac{z}{2x\left(z\right)}\mathrm{d}z\right)}^{1},$$(17)where $h=\sqrt{{{r}_{\mathrm{W}}}^{2}{a}^{2}}{r}_{\mathrm{A}}$ is the height of the prism in the z direction: the distance between the flat spot and the bottom of the axle. The length of the prism is given by $x$, which is now a function of $z$. The function has been inverted inside the integral (and then again outside after integration) to reflect the fact that when integrating along the height of the prism, one is essentially adding a number of springs in series, and thus requires an inverse summation.
Figure 12 shows the estimated wheel stiffness found using equations (10) and (17), as well as the percent error between the two as a function of axle/wheel radii ratio and contact area halflength. This formulation, on top of being quite a bit more complicated, unfortunately yields results which are still not satisfactorily close to the those given by the wheel abacus. The equivalent trapezoidal prism results have an average percent error of 67% (median 55%) with respect to the wheel abacus, with some values as high as 197% for extremely small contact area halflengths. The standard deviation of the percent error is 9%.
Figure 12 (a) Normalized wheel stiffness: FE model results (blue) vs. an equivalent trapezoidal prism (red). Plotted on a logarithmic zscale (base 10). (b) Percent error between the FE model results and equivalent trapezoidal prism. 
While the trapezoidal prism does account for the changing cross sectional area of the equivalent shape, it still does not account for the decoupling of the lateral and vertical deformations (i.e. the absence of Poisson’s effect). It is believed that the curvature of the wheel is the source of this phenomenon, and thus cannot be captured by a trapezoidal prism.
Neither an equivalent beam nor an equivalent trapezoidal prism is sufficiently accurate to be used reliably as a replacement for the wheel abacus. On the other hand, this shows that the implementation of the wheel abacus into an indoor rolling noise model would indeed be a beneficial improvement.
7 Discussion
On a larger scale, the results shown in this work (particularly in regards to the discovery of the lack of influence of the Poisson effect) demonstrate that a similar modeling technique may be used to characterize the properties of other shapes as well. Characterizations have been performed for constant circular cross sections [14, 18] (for which the Poisson effect plays a large role), and characterizations of materials with constant square cross sections yield results which are nearly identical to those of a constant circular cross section below Poisson’s ratios of about 0.47. However, analysis of other nonconstant cross sections, or even more complex shapes, could potentially provide beneficial insights in applications and industries beyond rolling contact modeling.
7.1 Implementation into a rolling noise model
There are two ways in which the stiffness estimation procedures presented in this work may be implemented into a rolling noise model: pre or continuous calculation.
One option is to precalculate the stiffness of the wheel at the start of the model computation process. This may be used in models which operate in either the time or frequency domain. Here, Hertzian equations are used to compute the size of the contact area under static load in the absence of roughness, and the stiffness is then estimated based on this wheel geometry. Thus the wheel stiffness remains constant throughout the entire rolling model computation process.
In reality, however, the wheel stiffness is not constant. It depends on the size of the contact area, and any change in contact area halflength will result in a small change the wheel stiffness. Because the contact area halflength changes throughout the rolling event as a new roughness profile continuously “moves into” the area between the wheel and the floor, the wheel stiffness itself changes continuously as well. To account for this phenomenon, the second option is to calculate the wheel stiffness at each instant throughout the rolling event, providing a unique stiffness value for every discrete moment in time. Consequently, this may only be used in timedomain rolling noise models. It is more computationally expensive, but provides greater accuracy in the stiffness estimation.
In the continuouscalculation method, the presence of roughness in the contact area is accounted for in calculating the size of the contact area for the purpose of estimating the wheel stiffness (This is not to be confused with the local contact stiffness between the wheel and the floor, which represents the small scale deformation due to the interpenetration of the two bodies). In the precalculation method, a simplifying assumption is implicitly made that the presence of roughness will change the value of $a/{r}_{\mathrm{W}}$ from its original static Hertzian estimation by a small enough amount that it may be ignored for the purpose of calculating the wheel stiffness. In both cases, the value of $a/{r}_{\mathrm{W}}$ is always calculated prior to estimating the wheel stiffness.
Which method should be used depends on the methodology of the model (time or frequency domain), the magnitude of the roughness profile, and the priorities of the user. A roughness profile which is relatively smooth may not result in a contact area which changes greatly throughout the course of the rolling event, and thus precalculation may be more practical. Rougher profiles however, particularly those containing large discontinuities such as floor joints or wheel flats, may benefit from continuous calculation. Finally, as continuous calculation will take longer to complete than precalculation, users which prioritize short computation times may still choose to implement precalculation instead.
7.2 Scope of the method’s applicability
In both phases of the parametric study, care was taken to run FE models for a wide range of parameter values. Thus, for the most part, any kind of trolley wheel which may be reasonably expected to be found in the real world can have this method applied. Due to their linear dependence, for the wheel width, radius, and Young’s modulus, any value may be used. Axle radii of 5%–98% of the wheel radius are valid, which for all intents and purposes is comprehensive, as a wheel outside this range would be of no practical use. For the contact area halflength, the method is assumed to be valid for any value below 10% of the wheel radius. Again, the wheel softness and applied load would need to be so high to achieve a contact area halflength above this limit, that in reality it is of no concern. For contact area halflength values below 0.3% of the wheel radius, the wheel stiffness is assumed to continue to follow the polynomial profile. This likely results in an overestimation of the stiffness for these extremely small contact area halflengths. However, this is considered acceptable in order to allow the model to allow for loss of contact between the wheel and the floor: a possibility in the presence of wheel flats and/or floor joints. In this type of situation, the size of the contact area tends to zero before disappearing at the moment of loss of contact.
8 Conclusion
This work presents an original technique for estimating the stiffness of a cylindrical indoor trolley wheel. A parametric study was conducted in order to identify the dependence of the wheel stiffness on each of the relevant variables. The stiffness is linearly dependent on the wheel width and Young’s modulus, and largely independent from Poisson’s ratio. A unique dependency exists for the axle/wheel radii ratio and contact area halflength. Using the information from these relationships, an abacus was created for estimating the stiffness of virtually any cylindrical wheel. Three methods were presented for estimating the wheel stiffness: a lookup table, a series of onedimensional polynomials, or a single twodimensional polynomial. The polynomials have extremely good fit to the data, and all three methods provide estimates that are reasonably close to one another. These methods may be implemented into a rolling noise model to provide either a precalculated or continuously updating estimation of the wheel stiffness.
Acknowledgments
This work was done as part of Acoutect: an innovative training network composed of five academic and seven nonacademic participants. This consortium comprises various disciplines and sectors within building acoustics and beyond, promoting intersectoral, interdisciplinary and innovative training and mobility of the researchers within the project. This project has received funding from the European Union’s Horizon 2020 research and innovation program under grant agreement No. 721536. This work was also performed within the framework of the Labex CeLyA of Université de Lyon, operated by the French National Research Agency (ANR10LABX0060/ANR11IDEX0007).
References
 P.J. Remington: Wheel/Rail Rolling Noise: What do we know? What don’t we know? Where do we go from here? Journal of Sound and Vibration 120 (1988) 203–226. https://doi.org/10.1016/0022460X(88)904300. [Google Scholar]
 A. Kuijpers, G. Blokland: Tyre/Road Noise Models in the Last two Decades: A Critical Evaluation. The Hague, Holland, 2001. [Google Scholar]
 M. Edwards, F. Chevillotte, L. Jaouen, F.X. Bécot, N. Totaro: Rolling Noise Modeling in Buildings, Chicago, IL USA, 2018. [Google Scholar]
 P. Remington, J. Webb: Estimation of wheel/rail interaction forces in the contact area due to roughness. Journal of Sound and Vibration 193 (1996) 83–102. https://doi.org/10.1006/jsvi.1996.0249. [Google Scholar]
 D.J. Thompson, B. Hemsworth, N. Vincent: Experimental validation of the TWINS prediction program for rolling noise, Part 1: Description of the model and method, Journal of Sound and Vibration 193 (1996) 123–135. https://doi.org/10.1006/jsvi.1996.0252. [Google Scholar]
 T.X. Wu, D.J. Thompson: A double Timoshenko beam model for vertical vibration analysis of railway track at high frequencies, Journal of Sound and Vibration 224 (1999) 329–348. https://doi.org/10.1006/jsvi.1999.2171. [Google Scholar]
 T. Mazilu: Green’s functions for analysis of dynamic response of wheel/rail to vertical excitation, Journal of Sound and Vibration 306 (2007) 31–58. https://doi.org/10.1016/j.jsv.2007.05.037. [Google Scholar]
 A. Pieringer, W. Kropp, D.J. Thompson: Investigation of the dynamic contact filter effect in vertical wheel/rail interaction using a 2D and a 3D nonHertzian contact model, Wear 271 (2011) 328–338. https://doi.org/10.1016/j.wear.2010.10.029. [Google Scholar]
 T. Clapp, A. Eberhardt, C. Kelley: Development and validation of a method for approximating road surface textureinduced contact pressure in tirepavement interaction, Tire Science and Technology 16 (1988) 2–17. [CrossRef] [Google Scholar]
 F. Wullens, W. Kropp: A threedimensional contact model for tyre/road interaction in rolling conditions, Acta Acustica United with Acustica 90 (2004) 702–711. [Google Scholar]
 E. Rustighi, S.J. Elliott: Stochastic road excitation and control feasibility in a 2D linear tyre model, Journal of Sound and Vibration 300 (2007) 490–501. https://doi.org/10.1016/j.jsv.2006.06.076. [Google Scholar]
 W. Kropp, P. Sabiniarz, H. Brick, T. Beckenbauer: On the sound radiation of a rolling tyre, Journal of Sound and Vibration 331 (2012) 1789–1805. https://doi.org/10.1016/j.jsv.2011.11.031. [Google Scholar]
 H. Hertz, D.E. Jones, G.A. Schott: Miscellaneous Papers, London: Macmillan. Macmillan and co., New York, 1896. [Google Scholar]
 S. Sim, K.J. Kim: A method to determine the complex modulus and poisson’s ratio of viscoelastic materials for FEM applications, Journal of Sound and Vibration 141 (1990) 71–82. https://doi.org/10.1016/0022460X(90)90513Y. [Google Scholar]
 T. Paviot: pythonOCC (Jun 2007). https://www.pythonocc.org. [Google Scholar]
 H. Si: Tetgen (2013). http://www.tetgen.org. [Google Scholar]
 F. Hecht, O. Pironneau, A. Le Hyaric: FreeFem++ (July 2017) https://freefem.org. [Google Scholar]
 ISO 184375:2011: Mechanical vibration and shock – Characterization of the dynamic mechanical properties of viscoelastic materials – Part 5: Poisson ratio based on comparison between measurements and finite element analysis, 2014. https://www.iso.org/standard/50428.html. [Google Scholar]
Cite this article as: Edwards M, Chevillotte F, Bécot FX, Jaouen L & Totaro N. 2020. Polynomial relations for cylindrical wheel stiffness characterization for use in a rolling noise prediction model. Acta Acustica, 4, 4.
All Tables
All Figures
Figure 1 Description of the wheel geometry. 

In the text 
Figure 2 Results of the parametric study: (a) Young’s modulus $E$ and (b) relative wheel width $w/{r}_{\mathrm{W}}$. Points having the same x coordinate (but a different y coordinate) are for varying values of Poisson’s ratio $\nu $. Schematics of the extreme wheel geometries are shown for reference (drawn to scale with respect to one another). 

In the text 
Figure 3 Results of the parametric study: (a) relative axle radius ${r}_{\mathrm{A}}/{r}_{\mathrm{W}}$ and (b) relative contact area halflength $a/{r}_{\mathrm{W}}$. Points having the same x coordinate (but a different y coordinate) are for varying values of Poisson’s ratio $\nu $. Schematics of the extreme wheel geometries are shown for reference (drawn to scale with respect to one another). 

In the text 
Figure 4 Three example wheels: (a) A realistic onepiece wheel with a small solid axle, (b) an unrealistic onepiece wheel with an impractically large solid axle, and (c) a realistic twopiece wheel with a soft outer wheel layer, a hard inner wheel core, and a solid axle. In the scenario of a high wheel/axle ratio, wheel (c) is being exhibited, not wheel (b). 

In the text 
Figure 5 Results of the parametric study: (a) wheel radius ${r}_{\mathrm{W}}$ and (b) input displacement ${z}_{\mathrm{in}}$. Points having the same x coordinate (but a different y coordinate) are for varying values of Poisson’s ratio $\nu $. Schematics of the extreme wheel geometries are shown for reference (drawn to scale with respect to one another). 

In the text 
Figure 6 Diagram of the wheel compression under an arbitrary load $Q$. The movement of the wheel into the areas formed between the wheel and the floor overrides any effect a change in Poisson’s ratio may have on the stiffness. 

In the text 
Figure 7 Normalized wheel stiffness as a function of axle/wheel radii ratio and relative contact area halflength. Plotted on a logarithmic zscale (base 10). 

In the text 
Figure 8 (a) Nineth order 1D polynomials for normalized wheel stiffness as a function of axle/wheel radii ratio. (b) Sixth order 1D polynomials for normalized wheel stiffness as a function of contact area halflength. Plotted on a logarithmic yscale (base 10). 

In the text 
Figure 9 Fifth order homogeneous bivariate polynomial for normalized wheel stiffness as a function of both axle/wheel radii ratio and relative contact area halflength. Plotted on a logarithmic zscale (base 10). 

In the text 
Figure 10 Example of the 1D polynomial characterization procedure for ${r}_{\mathrm{A}}/{{r}_{\mathrm{W}}}_{1}=0.30$ and ${r}_{\mathrm{A}}/{{r}_{\mathrm{W}}}_{2}=0.897$. Plotted on a logarithmic yscale (base 10). 

In the text 
Figure 11 (a) Normalized wheel stiffness: FE model results (blue) vs an equivalent beam (red). Plotted on a logarithmic zscale (base 10). (b) Percent error between the FE model results and equivalent beam. 

In the text 
Figure 12 (a) Normalized wheel stiffness: FE model results (blue) vs. an equivalent trapezoidal prism (red). Plotted on a logarithmic zscale (base 10). (b) Percent error between the FE model results and equivalent trapezoidal prism. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.