Exploration of the optimum ﬁnite element modelling techniques for honeycomb structures for non-pneumatic tyre applications

Honeycomb spokes are among the most commonly used structures in non-pneumatic tyre (NPT) research and development. Even though ﬁnite element (FE) modelling plays a key role in this effort, there is still a lack of knowledge regarding the requirements for accurate FE simulation of the mechanical behaviour of NPT honeycomb structures. The use of an inappropriate FE type can lead to misleading results and act as a barrier for robust research and development. To address this gap in literature


Introduction
A Non-Pneumatic Tyre (NPT) typically consists of flexible solid spokes and does not rely on pressurised air for support. Therefore, it cannot puncture and does not require pressure maintenance. Moreover, the direct link between spoke design and mechanical behaviour also significantly enhances its capacity for customisation [1]. NPTs have been successfully used in automotive, aerospace, and other applications, and ongoing research and development exists to help transform these airless tyres for use in passenger vehicles [2].
The main disadvantage of NPTs is the inability to disperse heat build-up effectively when travelling at high speeds (!50mph) [3,4]. This is due to the spokes rapidly deforming which also causes excessive vibration. This makes them more applicable in low-speed applications such as tractors, mowers, bicycles, and wheelchairs [1,5].
A typical NPT comprises four major components: an outer rubber tread (similar to pneumatic tyres), a polymeric shear-band and/or metallic outer-ring internally attached to the tread, concentric polymeric spokes internally connected to the shear-band/ outer-ring, and an inner metallic hub which sandwiches the spokes inside [1,6]. The tread is in contact with the road and ensures sufficient traction, the shear-band/outer-ring controls and maintains the shape of the tyre, and the metallic hub rigidly connects to the vehicle [7][8][9]. Last but not least, the spokes deform when subjected to a load and return to their original shape when the load is removed to provide support and suspension [7]. The geometry of the spokes can be altered to change its vertical stiffness and its shear stiffness which are directly linked to passenger comfort and rolling efficiency respectively [5]. NPT spokes exhibit nonlinear mechanical behaviour due to buckling which makes Finite Element (FE) modelling a challenging process, and so it is impor- tant to ensure the optimum strategy is used (e.g., element order, type, shape, and size) for accurate results.
Honeycomb structures are among the most commonly used designs of NPT spokes in the literature [6,10,11]. This is because honeycomb spokes can maintain high stiffness whilst being lightweight [12]. Another advantage of honeycomb structures is that their mechanical behaviour can be tuned according to the needs of specific applications or users by altering geometrical features such as spoke angle, height, and thickness ( Fig. 1a) [7,8].
Even though many studies on NPT honeycomb structures include the use of FE modelling, there is still a lack of knowledge regarding the requirements for accurate FE simulation of their mechanical behaviour. This need becomes more pronounced if one considers the range of different applications and therefore the broad range of potential sizes and mechanical requirements of honeycomb NPTs. A design parameter that can strongly affect the accuracy of different modelling techniques is the relative thickness-to-height ratio (T/H, Fig. 1a) of the honeycomb spoke. In this case, relatively thin spokes might be better simulated as 3D shells or relatively thicker spokes using 2D plane FEs.
The use of an inappropriate FE type can lead to misleading results and act as a barrier for robust research and development. In this context, the present study is the first that explores the optimum FE type for different relative dimensions of honeycomb NPT spokes. To this end, the mechanical behaviour of honeycomb NPT spokes will be modelled here using either quadrilateral or triangular 2D plane elements as well as using 3D shells. FE predictions for spokes of different relative sizes will be directly compared against original experimental results to identify the modelling approach achieving the highest accuracy. Methods and results from mechanical testing will be presented first followed by the FE analysis and the comparison between numerical and experimental results.

Methods
Preliminary analyses and literature [5,7] showed that loading in the NPT spokes is at its most intense in the spokes closer to the ground. To emulate this type of loading in the lab, samples of three-spoke clusters (Fig. 1b) were produced, and their mechanical response was studied under compression. These samples were designed assuming they belonged to a hypothetical wheelchair NPT with a diameter of 24 in. and 100 honeycomb spokes in total [13]. The design and dimensions of a reference honeycomb NPT was taken from literature [7]. The top and bottom surfaces of the samples were flat and parallel to one another to enable the application of compression using conventional compression plates. To assess the effect of thickness-to-height ratios, honeycomb designs were produced for plate thickness equal to 1 mm, 1.5 mm, 2 mm and 2.5 mm. All other dimensions were kept constant (Fig. 1a). This resulted in varied thickness-to-height ratios of the spokes, shown in Table 1. At the end, five samples were tested for each one of the four thickness scenarios. Samples were produced via fused deposition modelling (Ultimaker S3) with a soft thermoplastic polyurethane (Ultimaker TPU95A). The behaviour of the 1 mm thickness samples is highlighted in Fig. 1c.
To ensure the reliability of 3D printing, the plate thickness and depth of all samples were measured using a digital calliper. The measured thicknesses are shown in Table 1 and were deemed within an acceptable range.
During testing, the samples were compressed to 50% of the specimens' height using a 10kN loading frame (Insight Electromechanical, MTS Ò Systems) at 2 mm/min. This was to ensure post buckling behaviour was also captured. Force was recorded using a 10kN load cell (MTS Ò Systems). The recorded data was used to plot the force-displacement graphs of each sample (Fig. 1c).
Material characterisation of the 3D printing material was also performed to support FE modelling. To this end, dog-bone samples were designed according to the standardized test method for tensile properties of plastics (ASTM D638). Dog-bone samples were 3D printed using the same material and 3D printers as before. Quasistatic tensile testing (2 mm/min) was carried out for three samples using the same 10kN loading frame to measure their tensile stressstrain behaviour. Strain was measured using an optical extensometer (RTSS Videoextensometer, Limess Messtechnik & Software GmbH). Poisson's ratio was also calculated for one sample using the same extensometer rotated 90°to measure lateral strain. At the end, the measured stress-strain behaviour (averaged over three samples) and Poisson's ratio were used to calculate the mate- rial's hyperelastic coefficients according to the Mooney-Rivlin model (five parameter model) [14]: where W is the strain energy potential; I 1,2 are the deviatoric strain invariants; C10, C01, C20, C11, C02 are material coefficients characterizing the deviatoric deformation of the material; J and d are the compressibility parameters. Material coefficients were calculated using ANSYS Workbench. More specifically, the experimentally measured stress-strain behavior and Poisson's ratio were imported into ANSYS Workbench and the coefficients that achieved the best fit to the experimental data were calculated by the specialized tool within the software.

Results
All honeycomb samples buckled under loading leading to forcedeformation graphs that had a clear maximum value of force (Fig. 1c). More specifically, following bedding error (initial lack of fit between the sample surface and the compression plates) all honeycomb samples initially exhibited an almost linear increase in force with deformation. The start of buckling (%80% of maximum force) [15] caused a gradual continuous reduction in the slope of the force-deformation graph with increased deformation. The force-deformation graphs reached a maximum value beyond which force dropped due to buckling (Fig. 1c).
As expected, the samples became stiffer with increasing thickness. More specifically, increasing thickness from 1 mm to 1.5 mm, 2.0 mm or 2.5 mm led to an increase in maximum force of 147%, 461%, and 731% respectively ( Table 1). The observed substantial increase in maximum force highlights the importance of thickness for determining and optimising the mechanical behaviour of honeycomb spokes.

Design of honeycomb FE models
All FE simulations were conducted using Ansys Mechanical 2021 R2. Two different modelling approaches were followed to simulate the mechanical behaviour of honeycomb spokes. More specifically the compression of the honeycomb spokes was simu-lated either as a 2D plane-stress-with-thickness [14] problem or as a 3D problem using shell elements. For the 2D simulation, higher order quadrilateral (quad) and higher order triangular elements were considered (Plane183) (Fig. 2a). In these models, the honeycomb spokes and supports were meshed using the same element types and sizes. However, in the case of 3D shells, the honeycomb spokes were meshed using higher order 3D shell elements (Shell281), but the supports were meshed using 3D solid elements (Solid186). Appropriate pairs of contact elements (Conta175/ Targe170) were used to ensure the correct transfer of forces and moments between the shell elements of the honeycomb spokes and the 3D solid elements of the two supports (Fig. 2b).
The mechanical behaviour of TPU95 was simulated using the Mooney-Rivlin material coefficients that were previously experimentally calculated. For increased accuracy, the measured average spoke thicknesses and depths were assigned to the FE models (Table 1). For simplicity, each thickness condition is subsequently referred to using the initially assumed values (i.e. 1 mm, 1.5 mm, 2 mm and 2.5 mm).
To simulate compression, the models were assumed in simple contact with friction with rigid compression plates (friction coefficient = 0.62 [16]). A displacement of 5 mm was applied in the direction of compression to a master node connected to the top surface of the upper support via a multipoint constraint contact. The simulated bottom compression plate was fully fixed. Mesh convergence analysis was conducted to decide the most appropriate element size. It was concluded that at least 13,319 and 27,667 Plane183 elements were needed to minimise mesh dependency in the 2D plane for quads and triangles, respectively. 614 contact elements (Conta172, Targe169) were also needed. In the case of 3D Shells, 18,720 Shell281, 6606 Solid186, and 12,617 contact elements (Conta174, Conta175, Targe170) were needed. For an assessment of the relative computational cost of each model, the simulation time was also recorded. All simulations were done using a conventional personal computer system (Intel Core i7, 32 GB RAM).

FE results and comparison against experiments
All simulations led to buckling and the development of a maximum force for a compressive displacement that was smaller than the maximum applied. Based on that it could be said that all models were able to qualitatively simulate the observed mechanical response of the samples.
In quantitative terms, 2D plane models were able to accurately predict maximum force for the two intermediate thicknesses, but not for the thinnest or thickest. Indeed, the absolute difference to the experiments was < 3.5% for thickness 1.5 mm and 2 mm. Absolute error increased to 12% and 15% for 1 mm and 2.5 mm respectively. 3D shells were able to accurately predict maximum force for the thinnest model (1% error). However, error significantly increased with thickness. More specifically, deviation from experiments increased to 15%, 21% and 29% for thicknesses of 1.5 mm, 2 mm and 2.5 mm, respectively.

Discussion
The purpose of this study was to determine the optimum FE method through comparison of experimental and numerical results. The numerical results from the higher order 3D shells have the highest accuracy for the 1 mm thick spokes (1% error). At the same time, 2D plane elements appeared to be most accurate for thicknesses equal to 1.5 mm and 2 mm (<3.5% error). The accuracy of 2D plane models was reduced for the thickest scenario that was investigated (2.5 mm), but they remained significantly more accurate than 3D shells (15% and 29% error respectively).
Shell elements are expected to be accurate for thickness-toheight ratios ranging between 1/10 down to 1/100 [17]. The challenge with regards to honeycomb structures is that defining the relevant height for the calculation of thickness-to-height ratios is not always straight forward. In this study, thickness-to-height ratios were defined using the height of the entire spoke as reference. Based on this assumption, the thickness-to-height ratio for 1 mm thickness was 1/18. Shell elements stopped being accurate for thickness-to-height ratio of 1/12 or larger, indicating that 2D plane elements might be more appropriate for these relative thicknesses.
An alternative way to assess this ratio would be based on the length of the individual plates of the honeycomb spoke. In this case, the thickness-to-length ratio for thicknesses equal to 1 mm and 1.5 mm would be %1/5 and %1/3 respectively (Fig. 1a). These would fall out of the suggested range for using shell elements, however the results for 1 mm thickness are very accurate. Due to the spoke plates being interconnected, it can be suggested that they are acting as one entity rather than single plates, which highlights the entire spoke height as a more representative measurement to inform the selection of appropriate FE types.
The numerical results and computational time of the 2D analyses were very similar for both the higher order triangular and quad elements. As far as this analysis is concerned, quads and triangular elements can be considered as equivalent.
The relative lower accuracy of the 2D models for the 1 mm thickness could suggest that 2D elements are inaccurate for simulating structures with very small thicknesses. It could be further said that these elements give the most accurate results for thickness-to-height ratios of 1/12 and higher. Additional tests may be required to determine the exact threshold for where 2D elements start outperforming shell elements. For the scenarios tested here this threshold is somewhere in the range of 1/18 and 1/12.
These findings with regards to the applicability of shell and 2D plane elements are not limited to NPT applications but can be extended to the study of honeycomb structures in general. A key prerequisite for this is that the applied loading is relevant to the compression simulated here.
A limitation of this study was the exclusion of 3D solid elements within the spokes in the comparative analysis. Solid elements can be more accurate than 2D plane elements, but they are also significantly more computationally expensive. Preliminary simulations using 3D solid elements indicated that > 8 million Solid186 elements were needed to replicate the same number of elements across the thickness of honeycomb plates as in the 2D plane models resulting in a simulation time of > 5 h. Considering that the ultimate goal is to use these FE models for the design optimisation of entire tyres comprising a multitude of spokes (%100 spokes) [5], 3D solid elements were considered to be beyond the scope of this study. Another limitation of this study was that the testing of honeycomb spokes was limited to compression only. This was done because for the intended NPT applications, compression appears to be the most intense type of loading [8]. However, an NPT is subjected to spoke tension at the top of the tyre and shear at the sides when under vertical loading, as well as shear due to rotation. Even though it is highly unlikely that the key findings about the suitability of 2D plane and shell elements would change, the exact thresholds and optimum meshing could indeed change for different loading scenarios. The loading rate of the specimens was kept constant at 2 mm/min, and so the effect of the specimens' viscoelasticity is unknown. This could be an important parameter to investigate when designing NPTs for different applications as loading rates vary. However, same as before, including the viscoelastic behaviour of TPU is likely to change the absolute values of results but it is unlikely to affect the key conclusions of this study. Last but not least, the material coefficients of TPU were calculated using data from standardised tension only. Combining data from tension and compression would improve the accuracy of the simulations.

Conclusions
The selection of an appropriate FE type for the simulation of honeycomb NPT structures depends on the relative dimensions of the spoke and can significantly enhance the accuracy of the analysis. This study is the first to provide specific guidance about the type of element that should be used depending on the thicknessto-height ratios that would exist in different applications. It was found that: 3D shells should be used to simulate honeycomb spokes with thickness-to-height ratios of 1/18 and below. 2D plane elements should be used to simulate honeycomb spokes with thickness-to-height ratios of 1/12 or larger. Triangular and quadrilateral elements appear to be equally accurate. The definition of height is very important. The above suggestions are made assuming height as the height of the entire spoke (H), not the length (L) of individual plates (Fig. 1a).
These conclusions are transferable to applications of honeycomb structures beyond NPTs. A key prerequisite for this is that loading is relevant to the compression studied here.
CRediT authorship contribution statement

Data availability
Data will be made available on request.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.