## Abstract

When cylinders are packed and wrapped by the bands around the surface, the effective elastic behavior in the cross section of the assembly, which is of significance to its stability and integrity, can be controlled by the wrapping force in the band. The wrapping force is transferred to the cylinders through the Hertz contact between each pair of neighboring cylinders, which is validated by the experiments. The Singum model is introduced to study the mechanical behaviors of the packed cylinders with two-dimensional (2D) packing lattices, in which an inner cylinder is simulated by a continuum particle of Singum and the inter-cylinder force is governed by the Hertz contact model so as to derive the effective stress-strain relationship. The wrapping force will produce configurational forces given a displacement variation, which significantly changes the effective stiffness of the packed cylinders. The hexagonal packing exhibits isotropic elasticity whereas the square packing is anisotropic. The efficacy of our model is demonstrated by comparing the closed form elasticity against the numerical simulation and the previous models. The explicit form of elasticity can be used for packing design and quality control of cable construction and installation.

## 1 Introduction

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.

## 2 Problem Statement

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 ($\epsilon $) 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.

## 3 Formulations

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.

### 3.1 The Singum Construction and Modeling.

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\lambda lp0$ with $\lambda =lp/lp0$ being the deformation ratio. Under a hydrostatic load, *λ* for all bonds shall be the same, which is the case this paper investigates.

*C*

_{ijkl}is the stiffness tensor. However, the stiffness of a Singum particle is elastic but not linear. For a nonlinear elastic continuum, the tangential stiffness tensor at the spatial coordinate can be defined in the same fashion by the variations of the stress and strain at the current stress state, which will be illustrated subsequently.

*F*

^{I}at

**x**

^{I}at the boundary of each edge (

*I*= 1, 2, …, 6) of the Singum 0, and because of the equilibrium in the absence of body force, the boundary condition is written as

*δ*(

**x**) is the Dirac Delta function;

*σ*

_{ij}and

*n*

_{i}are the Cauchy stress tensor and surface out-normal vector of a continuum particle, respectively. The stress integral with a Singum particle

*S*

_{ij}can be written as follows [34]:

*V*

^{I}[27] as

*x*[36]

*δd*

_{ij}= 1/2(

*δu*

_{i,j}+

*δu*

_{j,i}) represents a linear displacement gradient tensor, which is related to the variation of the Eulerian strain at the current configuration of a stretch ratio

*λ*as [36]

*I*0 can be disregarded because each pair of the bond share the same center-center distance. It is shown in Fig. 3 that the total number of neighbors changes for different packing patterns and here is 6. The summation in Eq. (9) is reduced to the summation of $niInjI$ and $niInjInkInlI$, which can be written in the following identities for the hexagonal lattice

*V*(

*λ*) can be obtained by the experiments or Hertz’s model, which will be introduced in the next section. It is interesting that the hexagonal lattice exhibits an isotropic elasticity in the cross section, which has been observed in the graphene lattice as well [27,28].

### 3.2 Singum Potential for the Hertz Contact Problem.

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.

*P*per unit length, the half-width of the rectangular contact area is given by [35] as

*P*(

*λ*) can be obtained implicitly. The potential function

*V*(

*λ*) can be written as

*P*(

*λ*) can be numerically solved by Eq. (14) given

*λ*; and the integral is only the half of the bond so that a multiplier of two is applied. Alternatively, the formulation can be simplified by replacing $P(1\u2212\nu 2)/\pi Elp0$ with a single dimensionless variable for a more concise form. The derivatives of

*V*

_{,λ}(

*λ*) and

*V*

_{,λλ}(

*λ*) are written as follows:

*P*can be further derived from Eq. (14) as

*λ*and

*P*is given by Eq. (14), one can use the above equation to explicitly obtain the derivatives of

*V*and stiffness tensor

**C**given

*P*or

*λ*. Note that when

*λ*= 1,

*P*= 0, Eq. (17) provides

*P*

_{,λ}= 0 and

*P*

_{,λλ}→ ∞, which causes small but significantly changing elasticity when

*λ*is close to one.

### 3.3 Elastic Constants of the 2D Lattices Varying With the Wrapping Force.

*V*

_{,λ}(

*λ*) and

*V*

_{,λλ}(

*λ*) into the Singum model Eq. (9), one can compute the stiffness tensor of the packed cylinders given the 2D packing lattice structure, and the corresponding compliance matrix in the Voigt notation is written as

*E*

_{1}=

*E*

_{2},

*v*

_{1}=

*v*

_{2}and

*G*=

*E*/2(1 +

*v*), when the cylinders are distributed in other patterns, the elastic tensor can be anisotropic. For example, the square packing pattern will not satisfy

*G*=

*E*/2(1 +

*v*), which will be shown subsequently. Oblique packing or other patterns may not satisfy

*E*

_{1}=

*E*

_{2},

*v*

_{1}=

*v*

_{2}either.

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.

**C**and contact force

*P*can be obtained. Using the relation between

*P*,

*λ*, and wrapping force, one can design and control the stiffness, which will be demonstrated, subsequently.

*C*of square lattice and pairwise potential can be written as

## 4 Results and Discussion

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.

### 4.1 The Setup of the Numerical Experiments.

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\u2212\epsilon )$ and $y=Y(1\u2212\epsilon )$, respectively, under any given strain $\epsilon $. 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:

The Hertz contact in Eq. (13) assumed that contacted width 2

*b*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.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.

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.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.

### 4.2 The Compressibility of 2D Lattices Varying With *λ*.

*σ*

_{m}, which is generated by the wrapping force, will lead to an incremental mean stress as

*σ*

_{m}is related to

*P*as

*λ*of its original length for square and hexagonal lattices. The Singum prediction agrees very well with the numerical experiments, justifying the correctness of the Singum model. However, the computational resources consumed by Singum approach are much less than those taken by the numerical experiments because no particle generation or neighbor detection is needed for explicit solutions. For both lattices, as the cylinder’s Poisson’s ratio gets larger, a higher compressible load is required to compress the sample by the same amount. It is because the effective elastic modulus

*E*/2(1 −

*ν*

^{2}) in Eq. (13) increases with

*ν*. It is noted that at each

*λ*, the ratio between $\sigma mHex$ of hexagonal lattice and $\sigma mHex$ of the square lattice is

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.

### 4.3 Young’s Modulus and Poisson’s Ratio of 2D Lattices Varying With *σ*_{m}.

*λ*can be noticed from Eqs. (22)–(24) and Eqs. (27)–(29). Note that although both lattices show Young’s modulus and Poisson’s ratio the same in both directions, the hexagonal lattice is isotropic whereas the square lattice is not because it does not satisfy

*G*=

*E*/2(1 +

*v*). Putting them together with Eqs. (19)–(21), the effective Young’s modulus and Poisson’s ratio for square lattice can be derived as

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 − *ν*^{2}) 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.

### 4.4 Shear Modulus of 2D Lattices Varying With *σ*_{m}.

*E*and

*ν*, the effective shear modulus

*G*for square lattices is

*G*for hexagonal lattices is

*G*) varies with wrapping force for square and hexagonal lattices. The square lattice exhibits a negative shear resistance provided in Eq. (37), which is not physical but shows the instability of the square lattice under shear loading. In reality, with any shearing disturbance, the square lattice under a hydrostatic load will collapse and start transforming into the hexagonal lattice. Similarly to the effective

*E*, the

*G*of the hexagonal lattice increases with wrapping force. Also, compared to the loose case, a sudden jump in shear modulus can also be observed by wrapping the cylinders with only a small amount of force. A higher shear modulus can be observed for cylinders with higher

*ν*in Fig. 10(b).

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.

### 4.5 Comparisons With the Existing Models in the Literature.

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

*ν*and

*G*predicted by this model are all zero, which is different from the numerical experiments and the Singum prediction when a prestress exists or the lattice is subjected to a hydrostatic stress. The issue is caused by the effect of configurational change. The similar physics has been investigated by Eshelby of the existing force’s effect on the configuration change by crack propagation [43]. The concept of configurational force or material force has been applied to multiphysical problems [37]. The configurational stress in Eq. (8) captures the non zero Poisson’s ratio for square lattices, which highlights the physical rigor of the Singum model. On the other hand, when

*λ*= 1, for configurational stress is indeed zero with

*V*

_{,λ}= 0, the above equations can be recovered from the Singum model as well. In addition, unlike the tangent stiffness, the secant stiffness is not convenient to be adopted for stress updates in numerical simulations with incrementally increased strain. Additionally, the volume to which stress is homogenized was referred to as the undeformed volume, limiting the applicability to the infinitesimal strain range, while the Singum model can be straightforwardly extended to finite deformation with a high-fidelity interaction potential function.

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.

### 4.6 Case Study of a Suspension Bridge Cable.

*σ*

_{m}and the stretch ratio

*λ*at each wire contact point can be computed with Eq. (30), and here

*λ*=

*R*

_{d}/

*R*

_{0}relates the deformed and undeformed radius of the overall cable. The hoop stress

*σ*

_{h}in the band is related to

*σ*

_{m}as

*R*

_{d}and

*t*

_{h}are the deformed radius and thickness of a band. The wrapping force can be straightforwardly computed as

## 5 Conclusions

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.

## Acknowledgment

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.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

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

### Appendix: Numerical Verification and Experimental Validation of Johnson’s Model

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.

*E*and Poisson’s ratio

*ν*are calculated using the axial strain data from the experiment as follows.

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.

## References

*μ*Fe Model and DEM for an Assembly of Spheres Under Triaxial Compression