Effect of Wrapping Force on the Effective Elastic Behavior of Packed Cylinders

Understanding the cylinder or sphere packing problem is an art as much as a science [1], and this challenging problem has a wide spectrum of real-world applications in multiple industries, including textile production and packing, naval, automobile, and aerospace [2]. Many mathematical models have been developed to explore the closest packing [3–5], inspiring the optimized design with improved performances [1,6]. For example, in the civil engineering field, mechanical behaviors of granular materials, also closely related to sphere packing, play significant roles in many aspects, such as pavement construction [7], combating natural hazards (e.g., landslides) [8], excavation planning [9], etc.

Experimental testing can measure the macroscopic mechanical behavior of granular materials. Although the adoption of the celebrated photo-elastic grain technique can provide many significant discoveries in granular media [10], it is challenging to quantify the force transfer between grains, which is the origin or underlying mechanism of the mechanical behavior of granular materials [11]. These grain scale forces can be well captured by the adoption of discrete element method (DEM) [12,13], which models the contact between two spheres with the Hertz [14] and Mindlin-Deresiewicz [15] theories. DEM models each particle with both translational and rotational degrees-of-freedom. Although it can simulate the particulate behaviors of the granular materials (i.e., shear banding, necking, etc.), the results are sensitive to the scale and parameters, and thus it is computationally expensive to identify appropriate parameters and model scale in order to reach a practical and convergent prediction of the material behavior.

The recently proposed microstructure-based finite element (μFE) model for handling granular medium captures the natural depositional grain scale characteristics of sand (i.e., arbitrary shapes) [16]. Unlike DEM, μFE incorporates deformable grains so that the contact response emerges from the interaction of contacting bodies, thus enabling the modeling of irregular morphologies. However, its key drawback originates from the mesh generation: the surface mesh is a refinement of the constrained Delaunay tetrahedralization [17], and the volumetric mesh filling the grain with tetrahedral elements is bounded by the iso-surfaces. Compared with DEM, the number of degrees-of-freedom of each grain grows from six to hundreds or even thousands, resulting in a tremendous increase in computational cost [18] and hindering its superiority in making a difference in granular material simulation.

Continuum mechanics approaches can circumvent the high computational cost issue associated with DEM and μFE. The pilot attempt in modeling granular material with a continuum approach dates back to the middle of the last century when Duffy and Mindlin published a stress-strain relation for identical spheres packed in the face-centered cubic array [19]. This equation was derived by relating the contact force and displacements of a cubic unit cell through equilibrium and compatibility relations. This work opened the door for the following researchers to generalize this constitutive model to other regular packing patterns, such as simple cubic, tetrahedral, etc. [20–22]. Mindlin’s method, however, cannot be generalized to other packing patterns with non-cubic representative elements [22], but this difficulty was solved by stress homogenization over volume [23,24]. With the use of the energy conservation approach, the secant stiffness tensor can be obtained for all regular packing patterns [25]. All of these works employ field variables of intrinsic macroscopic nature without explicit connections with the underlying discrete material microstructure. These limitations motivated the development of micromechanical theory for elastic granular media with kinematic degrees-of-freedom included [26], but this theory is in linear form, which obviously contradicts the actual behavior at sphere contacts.

Although the force transfer through the contact between particles is through the stress on the contact surface, globally the load transfer through a granular material can be simplified by a lattice network between the center of particles with point-point forces, in which the force is correlated to the center-center distance by a potential function. The recently developed Singum model [27] uses the Wigner–Seitz (WS) cells of a lattice to represent a continuum solid so that the singular point forces can be transformed into the contacting stress between the continuum particle. By applying a displacement variation, from the relationship between the stress and the strain increments, we obtain the elastic constants. This procedure can be applied to general lattice networks and foam materials, which exist in nature or metamaterials, or composites, and the recent work demonstrated its application to a lattice metamaterial with harmonic potential or linear spring bonds [28].

This paper applies the Singum model [27,28] to the assembly of cylinders packing in certain patterns, which are equivalent to a 2D granular materials through the Hertz contacts and can be extended to 3D granular materials in future work. The solution exhibits significance in electric and civil engineering applications. High current electric transmission lines [29] and suspension bridge cables [30] are commonly using hundreds to thousands of wires packed in a certain pattern with wrapping bands. For example, Fig. 1(a) shows a suspension bridge supported by two large cables with banded wires, which can be observed by the cross-section of the cable in Fig. 1(b). Moreover, the cable wires (Fig. 1(b)) sustain a majority of the loads applied onto the deck and play an important role in the capacity and performance of the bridge [30]. These wire bundles are formed in a hexagonal arrangement tightened up by wrapping bands at a certain interval. The effective stiffness of packed cylinders in the cross section changes with the stress in the wrapping bands. It has been an empirical art to tighten the bands for the integrity and safety of the cable. A rigorous relationship between the stiffness and the wrapping force will be very useful for those applications, so that a formulation can be provided for the material and structural design given the cable and wire geometry and elastic constants.

In the following, the problem will be initially proposed. The Singum model [27] is constructed upon a hexagonal packing pattern, which can be generalized to square packing as well. The force-distance relationship between two neighboring cylinders can be formulated by a pairwise potential using the Hertz contact model. The experiments validate the potential function. The constitutive model for both square packing and hexagonal packing of cylinders is developed. The comparison with numerical simulation results proves the capability and accuracy of this model. The application of the Singum model to suspension bridge cable is demonstrated. The research output will enable scientists and engineers to efficiently predict the multi-scale mechanics of granular materials, thus inspiring the design of new metamaterials. Future studies will extend the current framework to study the poly-disperse of particles in the 3D space.

To illustrate the effective elastic behavior of packed cylinders with a wrapping force, numerical confined compression tests in 2D plane strain conditions are performed. Figure 2 shows a bundle of long cylinders confined in a container with four rigid side surfaces, and these cylinders are packed in regular patterns with N_x and N_y units along the x and y directions. For simplicity, we consider smooth identical cylinders with diameter d, elastic modulus E, and Poisson’s ratio ν. No friction is considered between the smooth cylinders. The 2D plane strain problem is assumed by constraining the displacement along the z direction.

All the boundary platens are assumed to be perfectly rigid. The bottom and left platens are hinged together with the bottom fixed on the ground, and the top and right platens are hinged together. The two parts are assembled together by two elastic links to wrap the packed cylinders while keeping the lattice structure. By shortening the two links simultaneously at the same rate, the contact area between cylinders increases while the lattice structure remains the same. Now keeping the wrapping force the same, we apply an infinitesimal displacement variation on the top platen. From the relationship between the displacement variation and external force, we can measure the tangential elastic modulus at the corresponding state of the wrapping force. Both shear and uniaxial loads can be applied to measure the elastic constants in different directions and loading modes.

Since materials have different Poisson ratios, in this research, the effect of a cylinder’s Poisson’s ratio on the overall mechanical properties is studied. Two lattice structures will be considered in this 2D study: square packing and hexagonal packing. In this work, we concentrate on the elastic behavior of 2D lattices through the normal forces of the Hertz contact between cylinders, while the tangential and torsional forces giving rise to irreversible deformations are ignored for the smooth surface. Earlier experimental studies stated that the contributions of shearing and torsional grain contacts are negligible to the volumetric elasticity [31–33], proving the validity of our setting. These two packing patterns show different force transmission mechanisms so the bond length (r) to applied axial strain (⁠$ε$⁠) relation is treated on a case-to-case basis in the subsequent sections. For both cases, we compute the axial stress (σ_a) and confining stress (σ_c) with the Singum model for further analysis. To validate the Singum model, our predictions are compared with the direct numerical simulation results computed with our implemented matlab code, and this development process will be introduced in the following section.

This section briefly introduces the Singum model, followed by the derivations for the inter-particle potential for the Hertz contact problem and general nonlinear pairwise interactions, respectively. The elastic constants are calculated from the inter-particle potential function.

Yin [27] proposed the concept of Singum model to correlate the pairwise interaction with the elastic constants of solids, paving a way for cross-scale modeling. A Singum particle, constructed by Voronoi decomposition, occupies the space of a WS cell with a particle at the center, filling the entire domain without gaps. For example, a hexagonal packing pattern is illustrated in Fig. 3(a) with a unit cell including one cylinder with six neighboring cylinders. The Singum for cylinder 0 can be constructed in Fig. 3(b) by cutting the six bonds with perpendicular lines forming a hexagon. The original radius of each cylinder is $lp0$⁠, and the center-center distance or bond length will change with the interaction force, written as $2lp=2λlp0$ with $λ=lp/lp0$ being the deformation ratio. Under a hydrostatic load, λ for all bonds shall be the same, which is the case this paper investigates.

The Hertz contact theory deals with the mechanics at the contact between non-conforming solids. This theory builds up on the simplification that each body can be regarded as an elastic half-space loaded over a small elliptical region of its plane surface to calculate the local deformation and stress distribution [14]. The Hertz contact theory assumes infinitesimal contact strain and frictionless contact surface to make the aforementioned simplification justifiable.

Although Hertz’s model for the 3D spherical case has been well established [35], the 2D case cylinder contact problem exhibits different forms of force-deformation relations, mainly two categories: implicit and explicit. Johnson [35], Radzimovsky [38], and Goldsmith [39] models mutual approach as an implicit function of contact force in similar forms with logarithmic function included, which requires an iterative process for the solution of P at each given indentation δ, thus limiting their applications in computational programs [40]. In view of this shortcoming, Lankarani and Nikravesh [41] proposed a simplified explicit model considering energy dissipation during the impact process, making it well-suited for implementation, especially for dynamics problems. However, there is a parameter to be determined empirically. In this work, Johnson’s model is selected because its P − δ prediction well fits the results of both the finite element simulations and the experimental testings, which are shown in the  Appendix.

The wrapping force F will produce contact force P between cylinders and change λ from 1 at the undeformed state with the zero wrapping force. Therefore, once the stiffness tensor C is written in terms of λ, the dependence of the elastic modulus on the wrapping force can be obtained. In the following, two packing lattices are considered, respectively.

The Singum model provides the closed form of elasticity, which considers the effects of the wrapping force or contact force. This section will first verify the model by numerical experiments, demonstrate the accuracy of the model, and then apply it to the bridge cable design and analysis.

A matlab program is developed to perform numerical experiments serving as a benchmark for validating the proposed Singum model. The overall flow of this program is as follows: In the initialization process, an array of cylinders with a radius of $lp0$ are automatically generated based on the given lattice and the corresponding inputted N_x and N_y numbers. For example, Fig. 4(a) schematically illustrates N_x = N_y = 5 with 5 × 5 blue cylinders. A list of neighbors for each cylinder is detected and saved for the force computation step. To simulate a hydrostatic loading, which causes all the bonds to shrink by the same ratio, we update both x and y coordinates of all the particles to $x=X(1−ε)$ and $y=Y(1−ε)$⁠, respectively, under any given strain $ε$⁠. This deformation process is clearly illustrated in Fig. 4(b), where the contact forces emerge at the highlighted red overlaps. For each pair of particles in contact, namely x_i and x_j, the magnitude of the contact force is evaluated with Eq. (13), and its direction is accurately defined as the unit vector connecting the centers of two circles. The required hydrostatic wrapping force can be computed from the contact force P based on the assembly of the cylinders, which will be demonstrated subsequently.

The modeling results exhibit a certain variation when N_x and N_y are small due to the boundary effect. However, when N_x and N_y are larger than 20 × 20, the variation becomes negligible. In the following numerical experiments, N_x = 100 and N_y = 101 are used as the default configuration.

The following factors may affect the accuracy in actual applications:

1. The Hertz contact in Eq. (13) assumed that contacted width 2b is equal to the arc of the cylinder. When the contact area is larger, the applicability of this assumption becomes questionable, thus limiting the validity of the current contact model to infinitesimal strain only with the small contact area.

2. The friction between cylinders will play a significant role in actual experiments and applications with a shear load, whereas the present model assumes the perfectly smooth surface of cylinders.

3. The present potential function V(r) is derived from the linear theory of elasticity; whereas nonlinear elastic or inelastic behavior of the materials may produce a considerable discrepancy in actual applications as the contact zone exhibits stress concentration and thus large strain.

4. The displacement δ was calculated by the line integral of the strain along the center-centerline of the two-particle contact; whereas a particle is in contact with more than three particles in the actual applications.

Therefore, although the Hertz contact model has been widely used for granular materials in the literature, the accuracy of Eq. (13) is limited to infinitesimal strain conditions. However, the Singum model can be applied to general potential functions between particles. The present Hertz contact model can be straightforwardly replaced if a P − δ curve for large deformation can be developed numerically or experimentally, from which the potential function can be obtained by the path integral. Particularly, when the balls are hollow, large deformation is expected. Some experimental studies are underway.

In this section, we will use the Singum model to explore the mechanics and physics of packed cylinders with wrapping stresses. Without loss of the generality, we take the Young’s modulus of the cylinder to be E = 210 GPa, and consider a range of Poisson’s ratio, i.e., v = 0.1, 0.3 and 0.5, to check how it affects the effective material properties.

One reason is the packing density of hexagonal packing (0.907) is 15% higher than that of square lattice (0.785). The other reason stems from the difference in force transmission mechanism within the square and hexagonal lattices. The relation between wrapping force and the stretch ratio is the basis of subsequent analysis of the effective elastic modulus.

Figure 8 plots how the effective Young’s modulus (E) varies with σ_m for square and hexagonal lattices. Generally speaking, for both lattices, an increasing trend can be observed as the wrapping force increases, indicating that the effective E of packed cylinders can be manipulated by adjusting the wrapping force. A special point is that a sudden jump in E occurs at the moment when a very small wrapping force is applied compared to the loose condition because of P_,λ = 0 and P_,λλ → ∞ in the neighborhood of λ = 1. This jump indicates the power of wrapping on significantly increasing the effective E.

In addition, given the same wrapping force, the sample with a higher cylinder ν exhibits a higher effective E. This phenomenon may be due to the same reason that E/2(1 − ν²) in Eqs. (16) and (33) increases with ν. Comparing the square lattice with the hexagonal lattice, we can notice that to get a similar elastic modulus, a higher hydrostatic stress is required for the hexagonal lattice, and combined with σ_m − λ relation in Fig. 7, the difference in E at the same stretch is similar.

Figure 9 plots how the effective Poisson’s ratio (ν) varies with wrapping force for hexagonal lattices. Note that ν for square lattices increases from zero at an undeformed state to a small finite number as the sample is compressed. This unphysical phenomenon observed at first glance is actually true because the uniaxial compression will cause the side area parallel to this compression to become smaller, leading to an increase in stress, which gives the nonzero ν value. This is the power of configurational force. However, the Poisson’s ratio for hexagonal lattice is far from zero. When λ → 1, ν → 1/3; whereas ν increases as λ decreases.

Unlike the square lattice, the effective ν for hexagonal lattice is almost 0.35 although a slight increase can be observed by increasing the wrapping force. The reason for this finite ν stems from the lattice geometry: force can be transmitted from vertical direction to lateral direction through the incline bonds.

Overall, the prestress provides significant effects on the effective elasticity of the assembly of the lattice of the cylinders. The simple, explicit form of elasticity enables the design of lattice materials with programmable mechanical properties by adjusting the wrapping force. Although the Hertz contact model has been used for deriving the potential function of Eq. (15), the present model can be applicable to the general form of the potential function V(λ), which can be determined by the experiments directly or by other models.

In the field of constitutive modeling of granular materials, two main streams are kinematics and static approaches, which mainly solve for the movements of particles and the closed-form elastic modulus, respectively [42]. Because of similarities with the Singum model, the static approach is chosen for comparisons. Following Chang’s work [25], we derived the stiffness tensor components for both square and hexagonal lattices as

However, the present Singum model still exhibits its own limits that the established models address specifically, such as plastic deformation of particles and friction between particles, etc., which shall be considered in future work by including moment and shear force at the cutting points at the Singum surface.

In this research, we extended the Singum model by deriving the pairwise potential considering the Hertz contact between two cylinders, enabling the Singum model to efficiently predict the tangent stiffness tensors of particles packed in regular lattices in 2D plane strain conditions. To select an appropriate contact model, we performed experiments and finite element analysis on cylinder samples of different material properties, and Johnson’s model is selected in the Singum model for deriving the inter-cylinder potential. Both square and hexagonal lattices are considered in this research to show the versatility of the Singum model. The dependence of effective elastic modulus on wrapping force predicted by Singum model agrees very well with the numerical verification, regardless of the packing lattices and material properties of cylinders. It is interesting to show the hexagonal lattice exhibits isotropic elasticity while the square lattice anisotropic in the 2D space. The solution can be used in the design of cylinder packs with controllable mechanical properties via adjusting the wrapping force. The significance is demonstrated in our case study on designing cables for suspension bridges. The superiority of the Singum model is demonstrated by performance comparison with common strategies for constitutive modeling of granular materials in literature. In addition to the 2D lattices, the Singum model can be extended to modeling granular materials consisting of spheres packed in 3D lattices, such as face-centered cubic, body-centered cubic, and simple cubic. The research output will shed light on investigating the mechanics of packed cylinders and exploring the optimized design of packing problems.

This work is sponsored by the National Science Foundation IIP #1738802, IIP #1941244, CMMI#1762891, and U.S. Department of Agriculture NIFA #2021-67021-34201, whose support is gratefully acknowledged. Dr. Zadshir and Dr. Yin thank the support from NASA SBIR #80NSSC22PB164. We thank Dr. Liming Li and Mr. James Basirico for their insightful discussion and help with the experiments. We are specially grateful to Professor Raimondo Betti for sharing his photos on bridge cables with us. Dr. Yin conceptualized and supervised the project and paper writing; Cui conducted the numerical studies and drafted the paper; Dr. Zadshir co-advised the project and supervised the experiments; and Teka conducted the experiments and data analysis and wrote the experimental part.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

As a result of the complexity of the cylinders in the contact problem, various types of models have been proposed to describe the relationship between contact force (P) and mutual approach (δ). The most well-known models include Johnson’s model [35], Radzimovsky’s model [38], Goldsmith’s model [39], and Lankarani and Nikravesh’s model [41]. Following these models, the relation between deformed cylinder radius l_p and contact force P is summarized below:

Experimental tests were conducted at the Carleton Laboratory, Columbia University, to investigate the load-displacement (P − δ) relationship. A universal testing machine (UTM) with a maximum capacity of 34 kips (150 kN)was used to apply a compression load. An abrasion-resistant polyurethane rubber rod was acquired from McMaster with Part Number 8695K693 in July 2022 and cut into three specimens with a diameter 54 mm, and varying lengths of 97.8 mm, 103 mm, and 104 mm, respectively.

Using Eq. (A5), the Young’s modulus (E) and Poisson’s ratio (ν) are calculated to be 470 mPa and 0.5, respectively in the linear elastic range.

In order to determine the P − δ relationship of the cylindrical rubber, experimental tests were performed on the three specimens as shown in Fig. 12(b). The specimens were laid in the horizontal direction and held in place using steel plates that were oiled to prevent friction. The specimens were loaded uni-axially using a displacement control of 0.762 mm/min. The time, load, deflection values for all tests were recorded through a data acquisition system. The load-displacement data of all three specimens were averaged and used for comparison as shown in Fig. 13. In addition, the error bars shown in the figure, demonstrate that the load-displacement relationship of the three specimens is very close to each other. Note that although only one cylinder is used in the test, because the steel platens can be considered to be a rigid surface, and with a mirror symmetry it can reproduce the deformation pattern of the contact of two identical cylinders in the finite element method (FEM) simulation of Fig. 12(c).

In order to further verify that an appropriate contact model is selected, we performed finite element analysis with abaqus 2019. As shown in Fig. 12(c), the model geometry strictly follows the experimental configurations, and the material is set to be linearly elastic with Young’s modulus and Poisson’s ratio the same as measured in the experiment. The contact between two cylinders is defined as frictionless and hard contact, which minimizes the penetration of the secondary surface into the primary surface and does not allow the transfer of tensile stress across the interface.

Figure 13 shows the comparison between those aforementioned contact models. Note that the Lankarani and Nikravesh’s (LN) model exhibits a parameter n, which was recommended in the range of [1, 1.5]. Here, the case of n = 1 shows too much off from those from both experiments and finite element analysis; while other cases of n might fit the experiments better. Although the LN model exhibits the advantage of its explicit form solution, to avoid the empirical calibration of the parameter n, we turn to other three models. Johnson’s, Radzimovsky’s, and Goldsmith’s models provide similar predictions, which match reasonably well with the experimental result. Note that rubber exhibits hyperelastic behavior with nonlinear elastic moduli at different levels of strain. Because the single values of E and ν in the linear elastic range are used in the contact models, some deviations from the experimental results are anticipated. Using the linear elastic constants, the finite element results can provide another reference to the contact problem, which agrees very well with Johnson’s prediction. Therefore, Johnson’s model is chosen in this research to derive the potential function for the contact problem between two cylinders.

In actual applications, because the stress at the contacting surface and its neighborhood is much higher than the rest part, the non-linearity of elasticity or inelastic behavior of the material may affect the accuracy of the contact models. Moreover, the friction between particles may change the contact mechanics as well. For multiple contacts between particles, the pairwise contacts may exhibit some loss of accuracy. Therefore, although Johnson’s model provides good agreement with the present FEM and experimental results, its applicability to different materials may change with the load levels and testing geometry or configuration, particularly for finite deformation of many particle systems. More investigation of the applicability of those models is underway. However, as long as a high-fidelity P − δ curve is provided, the present Singum model can straightforwardly use it in the same fashion.