Optimizing the Honeycomb Spoke Structure of a Non-Pneumatic Wheel to Reduce Rolling Resistance

: Traditional pneumatic tyres are prone to puncture or blowout and other safety hazards. Non-pneumatic tyres use a high-strength, high-toughness support structure to replace the “airbag body” structure of pneumatic tyres, which is made of fibre skeleton materials and rubber laminated layers, thus effectively avoiding the problems of blowout and air leakage. However, discontinuous spokes undergo repeated bending deformation when carrying loads, which leads to energy loss, of which the rolling resistance of non-pneumatic tyres is one of the main sources of energy loss. This paper focuses on the study of gradient honeycomb non-pneumatic tyres. Firstly, a finite element model was established, and the accuracy of the model was verified by numerical simulation and stiffness tests. Secondly, the order of the effect of different spoke thicknesses on rolling resistance was obtained through orthogonal test analysis of four-layer honeycomb spoke thicknesses. Then, four optimized design variables were selected in combination with the spoke angles, and the effects of the design variables on rolling resistance were analyzed in detail by means of the Latin hypercube experimental design. Finally, the response surface model was established, and the non-linear optimization model was solved by the EVOL optimization algorithm considering the tyre stiffness limitations so that the rolling resistance was minimized. The results of the study laid down theoretical and methodological guidance for the design concept and technological innovation of low rolling resistance comfort non-pneumatic tyres.


Introduction
As one of the most important components of a vehicle, tyres are closely related to all aspects of vehicle performance, such as dynamics, fuel economy, braking, handling stability, smoothness, and so on [1].Pneumatic tyres have been widely used in a variety of moving vehicles due to their better vibration damping, load carrying capacity, grip performance, comfort, and other comprehensive performance aspects [2].However, pneumatic tyres have a number of disadvantages, such as the risk of flat tyres, complex manufacturing processes, low recycling rates, and so on [3].According to statistics, the proportion of traffic accidents caused by flat tyres reaches 40%, and the mortality rate of traffic accidents caused by flat tyres at high speeds is close to 100% [4].In contrast, a non-pneumatic tyre (NPT) does not require air for support, which avoids disadvantages such as flat tyres and air leaks, and at the same time has the advantages of being environmentally friendly, highly recyclable, simple to manufacture, and so on [5][6][7].Well-known tyre companies such as Michelin, Goodyear, and Bridgestone have launched corresponding airless tyres, with Michelin's Uptis range of NPTs already in use.
A non-pneumatic tyre usually consists of a rim, elastic spokes, reinforcements, shear bands, and tread, and tyre performance is largely determined by the elastic spokes [8].The main considerations in the design of NPTs are tyre mass, stiffness, ground pressure, and rolling resistance [9].Some spokes derived from bionic designs are able to combine these properties, such as honeycomb spokes, petal spokes, and so on [10].Honeycomb structures have a low mass in addition to high elasticity and stiffness [11].Therefore, the selection of honeycomb structures as spokes for non-pneumatic tyres is appropriate.Many scholars have done a lot of research in these areas.Zang et al. proposed a cellular structure design method based on the tangent method; built hexagonal and circular cellular structures; and investigated the displacement, radial stiffness, and ground pressure patterns of cellular structures with different densities.This laid a foundation for the study of contact performance analysis of non-pneumatic tyres under complex working conditions [12].Based on honeycomb non-pneumatic tyres, Deng et al. investigated the static and dynamic mechanical properties of tyres in good condition and in local damage condition and found that spoke damage has an effect on tyre performance [13].Jin et al. optimized spokes with honeycomb porous structure by the design of experiments and finite element method with the objective of minimum rolling resistance and assessed the sensitivity of the design variables by analysis of variance (ANOVA) [14].Different spoke configurations' effects on tyre contact pressure were compared by Aboul Yazid et al.The findings from the study showed that the honeycomb spoke structure reduces the contact pressure and the stress on the spokes of the tyre [15].In general, if a material has a high stiffness, it has a low elastic strain limit.Kim et al. designed the spokes of NPT with both stiffness and resilience by using an auxiliary honeycomb structure with negative Poisson's ratio [16].
An in-depth study of the honeycomb structure is more helpful in improving the performance of NPT.The regular honeycomb structure is generally a positive hexagonal cell with some energy absorption.It has been shown that irregular structures absorb the energy of load impacts better, resulting in better vibration damping and longer service life [17][18][19].The gradient honeycomb structure is an irrational hexagonal structure, which refers to the loss of periodicity in a certain direction.Sahu et al. investigated the anisotropic gradient honeycomb structure and analyzed its energy absorption capacity by 3D printing prepared samples for static and cyclic compression tests in the in-plane direction, which showed the highest damping capacity with mixed gradient results [20].Xu et al. found that thin-walled structures and porous materials can exhibit more efficient and effective energy absorption performance by introducing gradient properties and that structures and materials with functionally graded properties exhibit better crashworthiness [21].A new in-plane gradient honeycomb was proposed and investigated by Yong et al.The results showed that the compressive strength and energy-absorbing capacity of the positive gradient honeycomb structure were significantly enhanced compared to the non-gradient one [22].It can be seen that the gradient honeycomb structure does have better overall performance than the normal honeycomb structure.This structure can also be applied in the design of the spoke structure of NPT, which can also make NPT have good performance.Wu et al. designed a gradient inverse tetrachiral structure, applied it to NPT, and used the finite element method to study its mechanical properties under compressive load.The results showed that the gradient inverse tetrachiral structure and the NPT based on this structure have some similarity in the deformation mode [23].
The above studies show the good prospect of the gradient structure in the design of NPT spokes.Moreover, at present, NPT mainly adopts spokes with a uniform thickness, and the research mainly focuses on static performance analyses.There are relatively few studies on dynamic performance and rolling resistance, which also proves the novelty of the research in this paper.For this reason, this paper revolves around the rolling resistance of a gradient honeycomb non-pneumatic tyre, which is structured as shown below.In Section 2, a finite element model is developed, and the modelling approach is tested through experiments.In Section 3, an experimental design is used to analyze the sensitivity of different parameters of the spoke structure to the effect of rolling resistance so as to identify the significantly affected components.The optimized design of the significantly affected components is carried out to minimize the rolling resistance and maintain the stiffness performance.

NPT Components and Materials
The 3D models of the NPTs were initially established using the commercial Siemens UG NX12.0 software, while their corresponding finite element models were developed using the Altair HyperMesh software.The elemental composition of honeycomb NPTs is shown in Figure 1 and mainly consists of a tread, two reinforcement layers, a shear band, a number of honeycomb spokes, and a hub.In finite element simulation, the material of the hub is an aluminium alloy (AI 7075-T6), and the material of the two reinforcement layers is high-strength steel (ANSI 4340).Without reinforcements, the edges of the spokes over the contact zone with the ground would buckle and cause an undesirable non-linear effect on the honeycombs.The mechanical properties of AI 7075-T6 and ANSI 4340 are given in Table 1.The honeycomb spokes and shear band are made of polyurethane (PU).The PU composite has both elasticity and stiffness at the same time, and it has a relatively low modulus that allows for large strain with low stress.The tread is made of synthetic rubber.The hyper-elastic and viscous properties of the polyurethane and rubber are defined using the Ogden model and the Prony series, respectively, the material constants of which are identified via experimental tests, as summarized in Table 2.
affected components is carried out to minimize the rolling resistance and maintain th stiffness performance.

NPT Components and Materials
The 3D models of the NPTs were initially established using the commercial Siemen UG NX12.0 software, while their corresponding finite element models were developed using the Altair HyperMesh software.The elemental composition of honeycomb NPTs i shown in Figure 1 and mainly consists of a tread, two reinforcement layers, a shear band a number of honeycomb spokes, and a hub.In finite element simulation, the material o the hub is an aluminium alloy (AI 7075-T6), and the material of the two reinforcemen layers is high-strength steel (ANSI 4340).Without reinforcements, the edges of the spoke over the contact zone with the ground would buckle and cause an undesirable non-linea effect on the honeycombs.The mechanical properties of AI 7075-T6 and ANSI 4340 ar given in Table 1.The honeycomb spokes and shear band are made of polyurethane (PU The PU composite has both elasticity and stiffness at the same time, and it has a relativel low modulus that allows for large strain with low stress.The tread is made of syntheti rubber.The hyper-elastic and viscous properties of the polyurethane and rubber are de fined using the Ogden model and the Prony series, respectively, the material constants o which are identified via experimental tests, as summarized in Table 2.

Geometrical Parameters
The NPT specification was designed with reference to a 205/60 R15 pneumatic tyre The thicknesses of the hub, shear band, reinforcement layers (inner and outer ones), and thread were 5.9, 9.36, 0.64, and 7.9 mm, and the inner radii of the hub and thread were 20 and 315 mm, respectively.The width of the NPT was 205 mm.

Geometrical Parameters
The NPT specification was designed with reference to a 205/60 R15 pneumatic tyre.The thicknesses of the hub, shear band, reinforcement layers (inner and outer ones), and thread were 5.9, 9.36, 0.64, and 7.9 mm, and the inner radii of the hub and thread were 203 and 315 mm, respectively.The width of the NPT was 205 mm.
Honeycomb-shaped spokes are important supporting and energy-absorbing structures for NPTs, and the structural design of honeycomb spokes has a significant impact on various performance aspects of the NPTs.The honeycomb structure of the NPT is shown in Figure 2. It shows a reference cell of a honeycomb spoke structure, which achieves the honeycomb support spoke for NPTs by uniformly distributing the reference cell at 25 equal cycles in the circumferential direction.The length of each edge of the honeycomb spokes is a crucial parameter that changes the structure of the honeycomb and affects the performance of the NPT.The length of each edge can be determined based on its radial projection length, thereby quantitatively determining the specific structure of honeycomb spokes.The specific length parameters of each structure are shown in Table 3.
Appl.Sci.2024, 14, x FOR PEER REVIEW 4 of 17 Honeycomb-shaped spokes are important supporting and energy-absorbing structures for NPTs, and the structural design of honeycomb spokes has a significant impact on various performance aspects of the NPTs.The honeycomb structure of the NPT is shown in Figure 2. It shows a reference cell of a honeycomb spoke structure, which achieves the honeycomb support spoke for NPTs by uniformly distributing the reference cell at 25 equal cycles in the circumferential direction.The length of each edge of the honeycomb spokes is a crucial parameter that changes the structure of the honeycomb and affects the performance of the NPT.The length of each edge can be determined based on its radial projection length, thereby quantitatively determining the specific structure of honeycomb spokes.The specific length parameters of each structure are shown in Table 3.In order to facilitate the study of the effect of the spoke thickness gradient on rolling resistance, the length of each edge was a fixed constant, and the spokes were divided into four layers: Tp, Cp, Lp, and Bp, as shown in Figure 2.Moreover, the average value of the four-layer spokes' thicknesses was 5 mm, although each dimension was allowed to alter.

Meshing of NPT
In order to improve the computation efficiency, the tyre model needed to be appropriately simplified.The 3D model of the honeycomb NPTs was imported into HyperMesh for mesh partitioning.The two-dimensional cross-sectional meshes were first divided, with mixed quadrilateral/triangular meshes for the honeycomb spokes and quadrilateral meshes for the other components.The three-dimensional meshes were subsequently derived using the solid map feature available in the HyperMesh software.Once the meshing was complete, it was then imported into SIMULIA ABAQUS 2020 for simulation and analysis.According to the tyre simulation method [24], some wedge-shaped grids of the spokes used C3D6 elements, and the rest were hexahedron grids that used C3D8RH simplified integration elements.The NPTs of different schemes used the same grid partitioning method, with a total number of nodes of approximately 243,700 and a total number of units of approximately 177,700 for each NPT.The finite element model of the honeycomb NPT is shown in Figure 3.In order to facilitate the study of the effect of the spoke thickness gradient on rolling resistance, the length of each edge was a fixed constant, and the spokes were divided into four layers: Tp, Cp, Lp, and Bp, as shown in Figure 2.Moreover, the average value of the four-layer spokes' thicknesses was 5 mm, although each dimension was allowed to alter.

Meshing of NPT
In order to improve the computation efficiency, the tyre model needed to be appropriately simplified.The 3D model of the honeycomb NPTs was imported into HyperMesh for mesh partitioning.The two-dimensional cross-sectional meshes were first divided, with mixed quadrilateral/triangular meshes for the honeycomb spokes and quadrilateral meshes for the other components.The three-dimensional meshes were subsequently derived using the solid map feature available in the HyperMesh software.Once the meshing was complete, it was then imported into SIMULIA ABAQUS 2020 for simulation and analysis.According to the tyre simulation method [24], some wedge-shaped grids of the spokes used C3D6 elements, and the rest were hexahedron grids that used C3D8RH simplified integration elements.The NPTs of different schemes used the same grid partitioning method, with a total number of nodes of approximately 243,700 and a total number of units of approximately 177,700 for each NPT.The finite element model of the honeycomb NPT is shown in Figure 3.
different sites for different mesh dimensions are shown in Figure 4. Normalization of the results of the numerical simulations reveals that the variation between the peak NPT contact pressure and the maximum stress at different sites did not exceed 20% for differen mesh dimensions and was almost zero at 5 mm.Considering both computational accuracy and cost, the selection of a 5 mm grid size was appropriate.

Interaction and Constraints
A reference point was created at the wheel hub centre of the NPT.The spokes and the reference point of the hub were constrained with kinematic coupling.The RP was se lected as a control point for the coupling.Another reference point was also created on the ground surface and attached to the ground surface body element using a rigid body con straint, as shown in Figure 3.The rigid hub, the spokes, the inner reinforcement, the shear band, the outer reinforcement, and the tread were assembled together with a tie constraint.The contact and interaction between the rigid road surface and the tread was se up, with a 0.3 friction coefficient in the tangential direction and hard contact in the direc tion normal to the terrain.

Loads and Boundary Conditions
The numerical test on NPT rolling analysis consisted of two steps in this study: a vertical deflection and rolling.In the first step, the six degrees of freedom at the centre reference point of the NPT were fully restricted, and a vertical load of 3000 N was applied to the NPT by moving the rigid road surface upwards along the y-axis, which was ramped for a time period of 0.1 s.A vertical deflection for the load was obtained by searching the displacement at the reference point for contact between the tyre and ground.In the second step, the load applied in the first step was maintained while rotating the NPT around the central reference point at an angular velocity of 26 rad/s.The degrees of freedom to rotate The initial size of the grid element was 5 mm.In order to eliminate the influence of grid division on the finite element analysis results, the grid size of the original NPT was gradually refined to 1-10 mm.The NPT peak contact pressure and maximum stresses at different sites for different mesh dimensions are shown in Figure 4. Normalization of the results of the numerical simulations reveals that the variation between the peak NPT contact pressure and the maximum stress at different sites did not exceed 20% for different mesh dimensions and was almost zero at 5 mm.Considering both computational accuracy and cost, the selection of a 5 mm grid size was appropriate.
Appl.Sci.2024, 14, x FOR PEER REVIEW 5 of 17 The initial size of the grid element was 5 mm.In order to eliminate the influence of grid division on the finite element analysis results, the grid size of the original NPT was gradually refined to 1-10 mm.The NPT peak contact pressure and maximum stresses at different sites for different mesh dimensions are shown in Figure 4. Normalization of the results of the numerical simulations reveals that the variation between the peak NPT contact pressure and the maximum stress at different sites did not exceed 20% for different mesh dimensions and was almost zero at 5 mm.Considering both computational accuracy and cost, the selection of a 5 mm grid size was appropriate.

Interaction and Constraints
A reference point was created at the wheel hub centre of the NPT.The spokes and the reference point of the hub were constrained with kinematic coupling.The RP was selected as a control point for the coupling.Another reference point was also created on the ground surface and attached to the ground surface body element using a rigid body constraint, as shown in Figure 3.The rigid hub, the spokes, the inner reinforcement, the shear band, the outer reinforcement, and the tread were assembled together with a tie constraint.The contact and interaction between the rigid road surface and the tread was set up, with a 0.3 friction coefficient in the tangential direction and hard contact in the direction normal to the terrain.

Loads and Boundary Conditions
The numerical test on NPT rolling analysis consisted of two steps in this study: a vertical deflection and rolling.In the first step, the six degrees of freedom at the centre reference point of the NPT were fully restricted, and a vertical load of 3000 N was applied to the NPT by moving the rigid road surface upwards along the y-axis, which was ramped for a time period of 0.1 s.A vertical deflection for the load was obtained by searching the displacement at the reference point for contact between the tyre and ground.In the second step, the load applied in the first step was maintained while rotating the NPT around the central reference point at an angular velocity of 26 rad/s.The degrees of freedom to rotate

Interaction and Constraints
A reference point was created at the wheel hub centre of the NPT.The spokes and the reference point of the hub were constrained with kinematic coupling.The RP was selected as a control point for the coupling.Another reference point was also created on the ground surface and attached to the ground surface body element using a rigid body constraint, as shown in Figure 3.The rigid hub, the spokes, the inner reinforcement, the shear band, the outer reinforcement, and the tread were assembled together with a tie constraint.The contact and interaction between the rigid road surface and the tread was set up, with a 0.3 friction coefficient in the tangential direction and hard contact in the direction normal to the terrain.

Loads and Boundary Conditions
The numerical test on NPT rolling analysis consisted of two steps in this study: a vertical deflection and rolling.In the first step, the six degrees of freedom at the centre reference point of the NPT were fully restricted, and a vertical load of 3000 N was applied to the NPT by moving the rigid road surface upwards along the y-axis, which was ramped for a time period of 0.1 s.A vertical deflection for the load was obtained by searching the displacement at the reference point for contact between the tyre and ground.In the second step, the load applied in the first step was maintained while rotating the NPT around the central reference point at an angular velocity of 26 rad/s.The degrees of freedom to rotate the tyre's centre reference point along the z-axis were released, as shown in Figure 3.The tyre was rotated by pulling away from the ground in the x-axis direction, and the ground was displaced along the x-axis by 1 m.

Experimental Verification
In order to verify the feasibility of the above gradient spoke structure design, tyre stiffness tests were implemented, and qualitative and quantitative analyses were performed.Based on the geometrical parameters in Table 4 and the four-layer spoke thicknesses of 5 mm, 3 mm, 7 mm, and 5 mm, respectively, a physical reproduction was carried out using 3D printing technology.The experimental sample was a 1:2 model made of resin.The tyre stiffness test equipment was used to carry out the stiffness tests, as shown in Figure 5. the tyre's centre reference point along the z-axis were released, as shown in Figure 3.The tyre was rotated by pulling away from the ground in the x-axis direction, and the ground was displaced along the x-axis by 1 m.

Experimental Verification
In order to verify the feasibility of the above gradient spoke structure design, tyre stiffness tests were implemented, and qualitative and quantitative analyses were performed.Based on the geometrical parameters in Table 4 and the four-layer spoke thicknesses of 5 mm, 3 mm, 7 mm, and 5 mm, respectively, a physical reproduction was carried out using 3D printing technology.The experimental sample was a 1:2 model made of resin.The tyre stiffness test equipment was used to carry out the stiffness tests, as shown in Figure 5.The test specimen was mounted on the tester's mounting arm.A level was used to ensure that the tyre pressed evenly against the support platform and was perpendicular to the support platform.The camera was placed in front of the test unit at the same height as the wheel spokes.Finally, a piece of carbon paper was inserted between the test model and the support block.Once the test setup was prepared, the test specimen was gradually squeezed through the mounting arm, and the force transducer was used to collect the vertical force.Since the connection between the rim and the arm was rigid, the vertical The test specimen was mounted on the tester's mounting arm.A level was used to ensure that the tyre pressed evenly against the support platform and was perpendicular to the support platform.The camera was placed in front of the test unit at the same height as the wheel spokes.Finally, a piece of carbon paper was inserted between the test model and the support block.Once the test setup was prepared, the test specimen was gradually squeezed through the mounting arm, and the force transducer was used to collect the vertical force.Since the connection between the rim and the arm was rigid, the vertical displacement of the arm was measured using a micrometre to obtain the tyre deformation.The deformation of the wheel spokes was recorded by means of a video camera.The contact footprint of the test specimen was obtained through carbon paper.
Figure 6 shows the deformation information of the tyre during the simulation and experiment; it can be seen that the spoke deformation situation is very similar by comparison, thus qualitatively verifying the accuracy of the NPT model.Figure 7 demonstrates the change in tyre displacement as the force applied to the tyre gradually increases during the simulation and test.Comparing the two curves, it is not difficult to find that the simulation model basically matched the stiffness curve of the real object, thus quantitatively proving that the NPT modelling method has a high reliability.
experiment; it can be seen that the spoke deformation situation is very similar by comparison, thus qualitatively verifying the accuracy of the NPT model.Figure 7 demonstrates the change in tyre displacement as the force applied to the tyre gradually increases during the simulation and test.Comparing the two curves, it is not difficult to find that the simulation model basically matched the stiffness curve of the real object, thus quantitatively proving that the NPT modelling method has a high reliability.

Rolling Resistance Analysis and Optimization
In this section, the rolling resistance calculation method is first clarified.Then, an orthogonal experimental design was carried out, and the results were analyzed to select the optimal design variables.After that, 15 sets of sample points were extracted using the Latin Hypercubic Design of Experiments method, and the response surface model was developed after analyzing the simulation results.Finally, the model was solved using the EVOL algorithm with the objective of minimizing the rolling resistance.The whole process is shown in Figure 8. ison, thus qualitatively verifying the accuracy of the NPT model.Figure 7 demonstrates the change in tyre displacement as the force applied to the tyre gradually increases during the simulation and test.Comparing the two curves, it is not difficult to find that the simulation model basically matched the stiffness curve of the real object, thus quantitatively proving that the NPT modelling method has a high reliability.

Rolling Resistance Analysis and Optimization
In this section, the rolling resistance calculation method is first clarified.Then, an orthogonal experimental design was carried out, and the results were analyzed to select the optimal design variables.After that, 15 sets of sample points were extracted using the Latin Hypercubic Design of Experiments method, and the response surface model was developed after analyzing the simulation results.Finally, the model was solved using the EVOL algorithm with the objective of minimizing the rolling resistance.The whole process is shown in Figure 8.

Rolling Resistance Analysis and Optimization
In this section, the rolling resistance calculation method is first clarified.Then, an orthogonal experimental design was carried out, and the results were analyzed to select the optimal design variables.After that, 15 sets of sample points were extracted using the Latin Hypercubic Design of Experiments method, and the response surface model was developed after analyzing the simulation results.Finally, the model was solved using the EVOL algorithm with the objective of minimizing the rolling resistance.The whole process is shown in Figure 8.

Rolling Resistance
The rolling resistance of an NPT is calculated from the total viscoelastic energy loss

Rolling Resistance
The rolling resistance of an NPT is calculated from the total viscoelastic energy loss (ALLCD in ABAQUS) per unit of rolling distance [25].The total viscoelastic energy loss was calculated by multiplication of the elemental strain energy, W, and dissipation factor, tan δ: where N is the total number of elements.
A tyre is subject to cyclic loading during rolling, resulting in hysteresis.The hysteretic energy loss of a viscoelastic tyre material is the main contributor to rolling resistance in tyres, which constitutes about 90 to 95% of the total rolling energy loss.Rolling resistance is defined by the energy dissipated per distance rolled, which is given by where F R is the rolling resistance, W d is the energy dissipated, and D is the distance rolled by a tyre.

Design Variables and Design of Experiment (DOE)
When there are many influencing factors, some factors that do not have a significant effect on the response increase the amount of computation and waste computational resources.Therefore, in this work, orthogonal experiments were used to pick out the significant influencing factors as design variables to continue with the optimization design.The orthogonal experiment scheme and simulation results are shown in Table 4. FR represents the rolling resistance, and t1, t2, t3, and t4 represent the thicknesses of the Tp, Cp, Lp, and Bp, respectively, which are shown in Figure 2.
The polar analysis of orthogonal experimental results is shown in Table 5.K represents the sum of the test results for each level of each factor.K-avg represents the average of the test results for each level of each factor, which is equal to K divided by the number of levels of the corresponding factor.R represents the polar value of factors.The magnitude of its value reflects the degree of influence of the factors on the experimental results.The larger the R of the factors, the deeper the influence on the experimental results is approximately.The effects of the factors on the results of the experiment were, in descending order, t1, t2, t3, and t4.The relatively better levels corresponding to each factor were 7, 7, 7, and 5 by using the variance analysis.Despite its low workload and intuitive convenience, the extreme variance analysis method cannot distinguish between fluctuations in test data caused by changes in test conditions during the test and fluctuations in data caused by test errors.Therefore, the experimental data were analyzed by ANOVA, and the results are shown in Table 6.The SS stands for the sum of squared deviations, df stands for degrees of freedom; MS stands for mean square, which is equal to SS divided by df; and F is the statistic of the test, which is equal to the mean square of each item divided by MSD.The symbol '*' means that the factor has a significant effect on the results of the test with a 95 per cent confidence level.The aggregate represents the total sum of the squared deviations.Since the sum of the squares of intergroup deviations at the t4 factor was small, it was considered the sum of the squares of errors to test the significance of the t1, t2, and t3 factors.The degree of influence of the factors on the results of the experiment was consistent with the analysis of extreme variance method.Factors t1 and t2 had a significant effect on the test results, while factors t3 and t4 had no significant effect on the test results.Combining the analysis of polarity and analysis of variance of the test results, the top and middle spoke thicknesses were selected as the optimal design variables in terms of spoke thickness.Considering that the change in the spoke angle affects the length of the spoke, which will have an effect on the rolling resistance, the angles β and γ were comprehensively selected as the optimal design variables in terms of the spoke angle.On the basis of the above four optimal design variables, 15 sets of experimental scenarios were designed using the Latin hypercube sampling method, and the experimental scenarios and simulation results are shown in Table 7, where FR represents rolling resistance, TK represents vertical stiffness, X1 represents angle γ, X2 represents angle β, X3 represents the thickness of Tp, and X4 represents the thickness of Cp.The Latin hypercube sampling results were analyzed using the SIMULIA Isight software 2020.The correlations of the four design variables were first analyzed, and the results are shown in Figure 9. High correlations between multiple independent variables may lead to multicollinearity problems.Severe multicollinearity problems may lead to a loss of significance of variables, failure of the predictive function of the model, unreasonable parameter estimation, and so on.As can be seen from the figure, the correlation coefficients between the four independent variables were all below 0.3, and the correlation between them was low, which proves the reasonableness of the choice of independent variables for the design.The main effects of the four factors on the response are shown in Figure 10.It can be seen that all four factors had a large degree of influence on the rolling resistance.This was a non-linear relationship.A pareto plot was used to reflect the contribution of all factors to the rolling resistance, as shown in Figure 11, with blue indicating positive effects and red indicating negative effects.As can be seen from the figure, X1 2 and X3 2 had the largest contribution to rolling resistance, reaching more than 20 per cent; X4 2 and X2 2 also had a large contribution to rolling resistance; and the X2-X4 cross term, X1-X3 cross term, and X1-X2 cross term had more advanced contributions to rolling resistance.This is consistent with the results of the previous main effects analysis.
are shown in Figure 9. High correlations between multiple independent variables may lead to multicollinearity problems.Severe multicollinearity problems may lead to a loss of significance of variables, failure of the predictive function of the model, unreasonable parameter estimation, and so on.As can be seen from the figure, the correlation coefficients between the four independent variables were all below 0.3, and the correlation between them was low, which proves the reasonableness of the choice of independent variables for the design.The main effects of the four factors on the response are shown in Figure 10.It can be seen that all four factors had a large degree of influence on the rolling resistance.This was a non-linear relationship.A pareto plot was used to reflect the contribution of all factors to the rolling resistance, as shown in Figure 11, with blue indicating positive effects and red indicating negative effects.As can be seen from the figure, X1 2 and X3 2 had the largest contribution to rolling resistance, reaching more than 20 per cent; X4 2 and X2 2 also had a large contribution to rolling resistance; and the X2-X4 cross term, X1-X3 cross term, and X1-X2 cross term had more advanced contributions to rolling resistance.This is consistent with the results of the previous main effects analysis.Since the X2-X4 cross term and X1-X3 cross term had a large effect on the rolling resistance, the interaction between them was specifically analyzed.Figure 12a shows the interaction of X2 and X4 at low levels of X1 and X3.The rolling resistance gradually increases as X2 increases, while it decreases as X4 increases.Figure 12b illustrates the interaction of X2 and X4 at high levels of X1 and X3.It can be seen that the rolling resistance gradually increases as X2 increases, while it decreases as X4 increases.Figure 12c shows the interaction of X1 and X3 on rolling resistance at low levels of X2 and X4.As X1 increases, rolling resistance gradually increases, and as X3 increases, rolling resistance de-  Since the X2-X4 cross term and X1-X3 cross term had a large effect on the rolling resistance, the interaction between them was specifically analyzed.Figure 12a shows the interaction of X2 and X4 at low levels of X1 and X3.The rolling resistance gradually increases as X2 increases, while it decreases as X4 increases.Figure 12b illustrates the interaction of X2 and X4 at high levels of X1 and X3.It can be seen that the rolling resistance gradually increases as X2 increases, while it decreases as X4 increases.Figure 12c shows the interaction of X1 and X3 on rolling resistance at low levels of X2 and X4.As X1 increases, rolling resistance gradually increases, and as X3 increases, rolling resistance decreases and then increases.Figure 12d shows the interaction of X1 and X3 on rolling re- Since the X2-X4 cross term and X1-X3 cross term had a large effect on the rolling resistance, the interaction between them was specifically analyzed.Figure 12a shows the interaction of X2 and X4 at low levels of X1 and X3.The rolling resistance gradually increases as X2 increases, while it decreases as X4 increases.Figure 12b illustrates the interaction of X2 and X4 at high levels of X1 and X3.It can be seen that the rolling resistance gradually increases as X2 increases, while it decreases as X4 increases.Figure 12c shows the interaction of X1 and X3 on rolling resistance at low levels of X2 and X4.As X1 increases, rolling resistance gradually increases, and as X3 increases, rolling resistance decreases and then increases.Figure 12d shows the interaction of X1 and X3 on rolling resistance at high levels of X2 and X4.It can be seen that as X1 increases, rolling resistance gradually decreases, while as X3 increases, rolling resistance also decreases and then increases.

Response Surface Model (RSM)
The response surface model fits the design space using polynomial functions.The functional relationships are approximated relatively accurately on a local scale through fewer trials and are presented in simple algebraic expressions, which are computationally simple and bring great convenience to design optimization [26,27].With the experimental data obtained by DOE, a quadratic response surface model was created for subsequent parameter optimization of tyre rolling resistance.All possible combinations of terms were examined by exhaustive search to select the term that minimized the fitting error.Under the premise of ensuring the fitting accuracy, the expression was simplified to reduce the time of the subsequent optimization process.The simplified functional relationship between tyre rolling resistance and design variables is shown in Equation (3).
where is the rolling resistance, X1 represents angle γ, X2 represents angle β, X 3 represents Tp spoke thickness, and X4 represents Cp spoke thickness.In order to test the accuracy of the response surface model, an error analysis was carried out.The values predicted by the response surface model and the true values are shown in Figure 13a.The results of the error analysis are shown in Table 8.The mean value of the error between the predicted value and the true value was 0.10, the maximum value was 0.22, the root mean square value was 0.12, and the coefficient of determination was 0.86, and they were all within the acceptable range.

Response Surface Model (RSM)
The response surface model fits the design space using polynomial functions.The functional relationships are approximated relatively accurately on a local scale through fewer trials and are presented in simple algebraic expressions, which are computationally simple and bring great convenience to design optimization [26,27].With the experimental data obtained by DOE, a quadratic response surface model was created for subsequent parameter optimization of tyre rolling resistance.All possible combinations of terms were examined by exhaustive search to select the term that minimized the fitting error.Under the premise of ensuring the fitting accuracy, the expression was simplified to reduce the time of the subsequent optimization process.The simplified functional relationship between tyre rolling resistance and design variables is shown in Equation (3).
where F R is the rolling resistance, X1 represents angle γ, X2 represents angle β, X3 represents Tp spoke thickness, and X4 represents Cp spoke thickness.In order to test the accuracy of the response surface model, an error analysis was carried out.The values predicted by the response surface model and the true values are shown in Figure 13a.The results of the error analysis are shown in Table 8.The mean value of the error between the predicted value and the true value was 0.10, the maximum value was 0.22, the root mean square value was 0.12, and the coefficient of determination was 0.86, and they were all within the acceptable range.The constraints should consider the radial stiffness constraints of the tyre in addition to the range of values of the four design variables.The radial stiffness of the tyre affects its load carrying capacity: too small will lead to large deformation of the tyre, and too large will affect the smoothness of the car, so the tyre stiffness needs to be constrained in the optimization procedure.A multivariate quadratic regression equation for radial stiffness was constructed by calculating the radial stiffness of the tyre from the 15 scenarios in the previous Latin hypercube sampling method, which was simplified by retaining the key terms.The expression for the radial stiffness of the tyre is as follows.
where TK is the stiffness, X1 represents angle γ, X2 represents angle β, X3 represents top spoke thickness, and X4 represents middle spoke thickness.The values predicted by the regression equation and the true values are shown in Figure 13b.The results of the error analysis are shown in Table 8.The mean value of the error between the predicted value and the true value was 0.02, the maximum value was 0.06, the root mean square value was 0.03, and the coefficient of determination was 0.99, and they were all within the acceptable range.This proves that the regression equation has an excellent fit.The radial stiffness of the tyre was generally in the range of 150 to 400 kN/m.

Mathematical Model and Optimization Results
On the basis of the above conditions, the corresponding mathematical optimization model was established with the objective of minimizing the rolling resistance of the tyre, as follows:  The constraints should consider the radial stiffness constraints of the tyre in addition to the range of values of the four design variables.The radial stiffness of the tyre affects its load carrying capacity: too small will lead to large deformation of the tyre, and too large will affect the smoothness of the car, so the tyre stiffness needs to be constrained in the optimization procedure.A multivariate quadratic regression equation for radial stiffness was constructed by calculating the radial stiffness of the tyre from the 15 scenarios in the previous Latin hypercube sampling method, which was simplified by retaining the key terms.The expression for the radial stiffness of the tyre is as follows.
where TK is the stiffness, X1 represents angle γ, X2 represents angle β, X3 represents top spoke thickness, and X4 represents middle spoke thickness.The values predicted by the regression equation and the true values are shown in Figure 13b.The results of the error analysis are shown in Table 8.The mean value of the error between the predicted value and the true value was 0.02, the maximum value was 0.06, the root mean square value was 0.03, and the coefficient of determination was 0.99, and they were all within the acceptable range.This proves that the regression equation has an excellent fit.The radial stiffness of the tyre was generally in the range of 150 to 400 kN/m.

Mathematical Model and Optimization Results
On the basis of the above conditions, the corresponding mathematical optimization model was established with the objective of minimizing the rolling resistance of the tyre, as follows: The variables X1, X2, X3, and X4 are called the decision variables, and Equation ( 5) is called the objective function of the constraints of the problem.Obviously, this is a non-linear planning problem.Unlike linear planning, which has the simplex method as a general method, there are no general algorithms for non-linear planning that are suitable for all kinds of problems, and each method has its own specific scope of application.The difficulty in solving this type of problem is finding a globally optimal solution, and it is often easy to fall into the local optimal link in the solution process.Therefore, the optimization problem was solved using the evolutionary algorithm (EVOL).Compared to traditional optimization methods, evolutionary algorithms have the advantages of wide applicability, a high degree of nonlinearity, ease of modification, and parallelism.The final optimum value of the rolling resistance was 9.86 N, corresponding to the values of the optimized design variables of 134 • , 90 • , 6 mm, and 8 mm.In order to verify the optimization results of the response surface model, the optimized design variable values were used for finite element analysis.In order to verify the accuracy of the optimization results of the response surface model, finite element analysis (FEA) was carried out using the optimized values of the design variables, and the numerical simulation of rolling resistance was 10.62 N. The results of the finite element analysis were similar to the results of the response surface model optimization and were within the error allowance, verifying the accuracy of the response surface model's optimization results.

Spoke Stress and Contact Pressure
Because the spoke structure angles of the optimized NPT model were very close to the No. 12 model, the analysis on mechanical characteristics of the two models are selected to explain the differences in the rolling resistance.Figure 14 presents the comparison results of the spoke structures without the reinforcements.It can be seen that the spoke stress of the two models was similar: by increasing the distance from the rim to tread, the magnitude of stress decreased continuously.The position of the maximum stress values mainly appeared at the top spoke near the rim, and the second maximum stress at the middle spoke; however, the maximum mises value of the No. 12 was larger than that of the optimized NPT model.The radial stiffness K is defined as the vertical force divided by the vertical displacement.The stiffness K of the No. 12 model was 199.5 N/mm, and the stiffness K of the optimized model was 307.5 N/mm.This means that the optimized NPT model had a higher stiffness than the No. 12 model.Figure 15 shows the comparison of tyre contact pressure of the two models; it can be seen that the contact area of the optimized model was smaller than that of the No. 12 model because it had a higher stiffness, and the smaller area of the optimized model can lead to the higher contact pressure values under the same vertical force of 3000 N. In order to explain the differences in the contact pressure contour, the partial stress of the spokes within the contact patch are presented in Figure 16.It can be seen that there was a stress difference at the same positions of the spoke A and A': the higher the stiffness K, the more obvious the stress gradient.The higher stress of spoke A means that it will transmit more vertical load, and the two side spokes B and C may withstand relatively less load.The differences in the force transmission of the same position of the spokes are the main reason for the difference in the contact pressure distribution.
means that it will transmit more vertical load, and the two side spokes B and C may withstand relatively less load.The differences in the force transmission of the same position of the spokes are the main reason for the difference in the contact pressure distribution.

Conclusions
Based on the gradient honeycomb NPT, this paper optimizes the spokes of the tyre with the objective of rolling resistance reduction.The conclusions are shown below: (1) A finite element model was established and samples were fabricated by 3D printing for stiffness tests, and it was found that the local deformation of the actual spokes was similar to that of the simulation, which also proved the accuracy of the modelling method.
(2) The role of different factors in influencing the rolling resistance was clarified by the experimental design.The strongest non-linear effect on rolling resistance was found between angle γ and Tp spoke thickness, followed by the interaction between angle β and Cp spoke thickness and the interaction between angle γ and Tp spoke thickness.
(3) The difference in the vertical stiffness before and after the optimized model can make the force transmission of the same position of the spokes undergo change, and this means that it will transmit more vertical load, and the two side spokes B and C may with-stand relatively less load.The differences in the force transmission of the same position of the spokes are the main reason for the difference in the contact pressure distribution.

Conclusions
Based on the gradient honeycomb NPT, this paper optimizes the spokes of the tyre with the objective of rolling resistance reduction.The conclusions are shown below: (1) A finite element model was established and samples were fabricated by 3D printing for stiffness tests, and it was found that the local deformation of the actual spokes was similar to that of the simulation, which also proved the accuracy of the modelling method.
(2) The role of different factors in influencing the rolling resistance was clarified by the experimental design.The strongest non-linear effect on rolling resistance was found between angle γ and Tp spoke thickness, followed by the interaction between angle β and Cp spoke thickness and the interaction between angle γ and Tp spoke thickness.
(3) The difference in the vertical stiffness before and after the optimized model can make the force transmission of the same position of the spokes undergo change, and this

Conclusions
Based on the gradient honeycomb NPT, this paper optimizes the spokes of the tyre with the objective of rolling resistance reduction.The conclusions are shown below: (1) A finite element model was established and samples were fabricated by 3D printing for stiffness tests, and it was found that the local deformation of the actual spokes was similar to that of the simulation, which also proved the accuracy of the modelling method.
(2) The role of different factors in influencing the rolling resistance was clarified by the experimental design.The strongest non-linear effect on rolling resistance was found between angle γ and Tp spoke thickness, followed by the interaction between angle β and Cp spoke thickness and the interaction between angle γ and Tp spoke thickness.
(3) The difference in the vertical stiffness before and after the optimized model can make the force transmission of the same position of the spokes undergo change, and this will result in the difference in the spoke stress and contact pressure.The higher the stiffness, the lower the spoke stress, and the smaller the contact area.
This paper conducted a simulation analysis of rolling resistance and explored the influence of spoke structure parameters on rolling resistance.On the one hand, the result in this paper will provide methodological guidance for the rolling resistance analysis of non-pneumatic tyres with different spoke structures in the future, and on the other hand, it will also lay the foundation for optimizing the spoke parameters.Therefore, the paper provides a good foundation for the structural design and parameter optimization for the lower roll resistance of non-pneumatic tyres.

Figure 2 .
Figure 2. Honeycomb spoke structure parameters and spoke segmentation diagram.

Figure 2 .
Figure 2. Honeycomb spoke structure parameters and spoke segmentation diagram.

Figure 3 .
Figure 3. Modelling and boundary conditions of NPT.

Figure 3 .
Figure 3. Modelling and boundary conditions of NPT.

Figure 3 .
Figure 3. Modelling and boundary conditions of NPT.

Figure 7 .
Figure 7.Comparison of tyre stiffness simulation and experiment.

Figure 7 .
Figure 7.Comparison of tyre stiffness simulation and experiment.

Figure 7 .
Figure 7.Comparison of tyre stiffness simulation and experiment.

Figure 13 .
Figure 13.Comparison between the predicted curve of the response surface model and the actual value.

Figure 13 .
Figure 13.Comparison between the predicted curve of the response surface model and the actual value.
Appl.Sci.2024, 14, x FOR PEER REVIEW 15 of 17 means that it will transmit more vertical load, and the two side spokes B and C may withstand relatively less load.The differences in the force transmission of the same position of the spokes are the main reason for the difference in the contact pressure distribution.

Figure 14 .
Figure 14.Stress comparison of spoke structures.

Figure 14 .
Figure 14.Stress comparison of spoke structures.

Figure 15 .
Figure 15.Comparative analysis of contact pressure.

Figure 16 .
Figure 16.Comparison analysis of partial spoke structure stress.

Table 2 .
Polyurethane and synthetic rubber material parameters.

Table 2 .
Polyurethane and synthetic rubber material parameters.

Table 3 .
Specific parameters of honeycomb structure.

Table 3 .
Specific parameters of honeycomb structure.

Table 4 .
Orthogonal experiment schemes and simulation results.

Table 4 .
Orthogonal experiment schemes and simulation results.

Table 5 .
Extreme variance analysis results.

Table 8 .
Error analysis table.

Table 8 .
Error analysis table.