© M. Edwards et al., Published by EDP Sciences, 2020

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 small-scale 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 small-scale 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 spring-mass-damper system. For a vehicle with a built-in 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 half-length. 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 half-length a. The convention of half-length 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 re-meshed to add extra nodes in the area of contact for particularly small contact area half-lengths. This ensured the contact area had at least four nodes along its length. Depending on the geometry, mesh sizes ranged between 1.0 × 10⁴ and $5.3\times 10^4$ nodes, having degrees of freedom between $2.1\times 10^6$ and $1.1\times 10^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 16-Core 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 half-length 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 half-length values below 0.003 were not able to be reliably meshed.

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.

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 half-length. 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 half-length: 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 half-length $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 two-piece 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 one-in-the-same. 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 one-piece wheel with a small solid axle, (b) an unrealistic one-piece wheel with an impractically large solid axle, and (c) a realistic two-piece 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 horn-shaped 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 non-linearly dependent variables: the axle radius and the contact area half-length (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 half-length $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 half-length. Plotted on a logarithmic z-scale (base 10).

The normalized stiffness was then plotted as a function of the axle/wheel radii ratio and relative contact area half-length. The data is plotted with a logarithmic z-scale (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 one-dimensional polynomials, and a single two-dimensional 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 two-dimensional 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 one-dimensional 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}}=\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 one-dimensional polynomials for both varying axle/wheel radii ratio and contact area half-length. 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 half-length. Plotted on a logarithmic y-scale (base 10).

4.4 Two-dimensional polynomial

The third characterization method developed involves fitting a single two-dimensional 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}}=\exp \left[\sum _{N=0}^5\mathrm{}\left(\sum _{i=1}^N\mathrm{}{C}_{i,N-i}{f}^i{g}^{N-i}\right)\right].$$(7)

For ease of display, the expressions for axle/wheel radii ratio and contact area half-length have been replaced by functions $f$ and $g$, such that

$$f=\ln (r_{\mathrm{A}}/r_{\mathrm{W}}),$$(8)

$$g=a/r_{\mathrm{W}}.$$(9)

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}}=& \exp (95f^5+280f^{4}g-1 319f^3{g}^2+4 004f^2{g}^3-9 965f{g}^4+1 052 708g^5\\ & -249f^4-352f^{3}g+1 109f^2{g}^2+76f{g}^3-305 716g^4+247f^3+123f^{2}g\\ & -454f{g}^2+33 013g^3-117f^2+7\mathrm{fg}-1 625g^2+29f+43g-5)\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 half-length. Plotted on a logarithmic z-scale (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.

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 half-length 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 ₁ = 5811 kN and K ₂ = 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 y-scale (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=& -5 506 860(a/r_{\mathrm{W}}{)}^6+1 766 278(a/r_{\mathrm{W}}{)}^5-222 310(a/r_{\mathrm{W}}{)}^4\\ & +14 130(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=& -13 397 078(a/r_{\mathrm{W}}{)}^6+4 314 984(a/r_{\mathrm{W}}{)}^5-545 898(a/r_{\mathrm{W}}{)}^4\\ & +34 857(a/r_{\mathrm{W}}{)}^3-1 177(a/r_{\mathrm{W}}{)}^2+39(a/r_{\mathrm{W}})+0.3.\end{array}$$(14)

Plugging the relative contact area half-lengths 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 ₁ = 5954 kN and K ₂ = 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 ₁ = 5793 kN and K ₂ = 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.

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 FE-based 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 cross-section 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 half-length. 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 z-scale (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 cross-section 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 cross-sectional area changes along its height, from $2\mathrm{aw}$ on one end to $2r_{\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 half-length. 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 half-lengths. 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 z-scale (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 non-constant 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 pre-calculate 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 half-length will result in a small change the wheel stiffness. Because the contact area half-length 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 time-domain rolling noise models. It is more computationally expensive, but provides greater accuracy in the stiffness estimation.

In the continuous-calculation 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 pre-calculation 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 pre-calculation 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 pre-calculation, users which prioritize short computation times may still choose to implement pre-calculation 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 half-length, 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 half-length above this limit, that in reality it is of no concern. For contact area half-length 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 half-lengths. 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 half-length. 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 one-dimensional polynomials, or a single two-dimensional 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 pre-calculated 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 non-academic 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 (ANR-10-LABX-0060/ANR-11-IDEX-0007).

References

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 half-length $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 one-piece wheel with a small solid axle, (b) an unrealistic one-piece wheel with an impractically large solid axle, and (c) a realistic two-piece 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 half-length. Plotted on a logarithmic z-scale (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 half-length. Plotted on a logarithmic y-scale (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 half-length. Plotted on a logarithmic z-scale (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 y-scale (base 10).
In the text

	Figure 11 (a) Normalized wheel stiffness: FE model results (blue) vs an equivalent beam (red). Plotted on a logarithmic z-scale (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 z-scale (base 10). (b) Percent error between the FE model results and equivalent trapezoidal prism.
In the text