Rubber elasticity of realizable ideal networks

Cross-linked rubbers have been valuable materials from an industrial point of view since the discovery of vulcanization of uncross-linked natural rubber with sulphur by Goodyear in 1839.¹ They are used in a wide range of areas, such as tires, paints, and oil seals.² In terms of fundamental science, it is important to understand the properties of rubber materials, such as their entropic elasticity. The entropic elasticity is qualitatively different from the energetic elasticity that is usually observed in non-rubber materials such as metals and glasses and it is one of the most important properties on which a group of soft materials is based. Thus, to more effectively utilize cross-linked rubbers in industrial areas and further develop soft material science, it is crucial to deeply understand the properties of cross-linked rubbers.

It is well known that normal vulcanized rubbers under small uniaxial deformation have the following approximate relation:¹ 

where σ_n is the nominal stress, i.e., the stress obtained by dividing the tensile force exerted in the elongational direction by the cross-section in the unstrained state, G is the shear modulus, and λ is the extension ratio in the elongational direction. Several prominent scientists have proposed promising theories of the rubber elasticity to give a molecular interpretation to Eq. (1), such as Kuhn,³ Wall,⁴ Flory and Rehner,⁵ and James and Guth.⁶ The affine network theory proposed by Kuhn is still one of the most basic and useful molecular models. He assumed that the rubber networks are composed of only elastically effective bridged chains and have no structural defects, such as dangling chains, loops, and trapped entanglements,^1,7,8 which are contained in actual cross-linked rubbers. Here, we call such defectless networks “ideal” networks. Using his own ideal networks, Kuhn³ derived G in Eq. (1) from a molecular standpoint:

where ν is the number density of bridged chains, k_B is the Boltzmann constant, and T is temperature.

In actual cross-linked materials, there are a significant number of structural defects. In addition, partial chains with various lengths and cross-linking points are heterogeneously distributed. The partial chains have their own excluded volume (EV) and interact with each other by the repulsive interactions. These factors are not taken into account in Kuhn’s model. Regardless of the remarkable discrepancy between the assumptions in Kuhn’s theory and the actual situation, Eq. (2) has thus far been useful for experimental estimation of the degree of cross-linking from the shear modulus.¹ 

As a result of the unexpected success of Kuhn’s theory, many researchers are satisfied with the current state of understanding of cross-linked rubbers and have not considered the ideal networks underlying Kuhn’s theory. Recently, Sakai^9–11 and Shibayama^12,13 have performed a series of experimental studies on ideal networks with tetra-polyethylene glycol (PEG) gels. Tetra-PEG gels consist of partial chains with the same length and are much more homogeneous both spatially and topologically than normal cross-linked networks.^14–16 They have reported that the structural defects such as dangling chains and loops are actually negligible in Tetra–PEG gels. Nevertheless, tetra-PEG gels are expected to contain a significant number of trapped entanglements and grain boundaries, which can make tetra-PEG gels isotropic elastic bodies.

The effect of the EV of cross-linked rubbers on the elastic properties has not been theoretically nor experimentally clarified, although attracting attention of many polymer scientists.¹⁷ In practice, the EV interaction between any pair of constituent atoms cannot be turned off. In theory, this is a type of many-body problem and it is very difficult to obtain an exact analytic solution. Many studies^8,18–24 have succeeded in avoiding the difficulty by applying some type of approximate technique, such as the mean-field approximation. The outcomes are fascinating, but they are only approximations and cannot uncover the distinct difference in the elastic behavior owing to EV interactions.

Differing from experimental and theoretical approaches, coarse-grained molecular dynamics (CGMD) simulations are a promising tool to investigate defectless ideal networks. In CGMD simulations, the number and length of the partial chains, as well as the number and functionality of the cross-linking points, can be precisely controlled. Many coarse-grained simulations (including CGMD simulations) of cross-linked networks have been performed to uncover their mechanical properties.²⁵ However, most of the target structures are either randomly cross-linked networks^26,27 or randomly end-linked networks,²⁸ because one of the main aims is to clarify the elastic behavior of actual vulcanized rubbers with random structures. Everaers and Kremer^29–31 have investigated interpenetrating diamond networks, which do not have dangling chains and loops at all, by CGMD simulations. Their results are both interesting and instructive in thinking over the ideal networks. However, we should note that the interpenetrating networks contain many trapped entanglements and are much complicated than the ideal networks.

In CGMD simulations, it is easy to obtain hypothetical networks without EVs, which can be achieved by only turning off the non-bonding interactions between constituent particles in silico. Gao and Weiner^32–34 investigated the EV effects on the force–length relation of single chains and the stress–extension relation of cross-linked networks by coarse-grained simulations. They discussed the elasticity of the polymeric systems in terms of the configuration space available to each partial chain.³² However, they did not thoroughly investigate the difference in the shear modulus of cross-linked networks owing to EV interactions. In the following, we call hypothetical networks without EV interactions “phantom” networks and ordinary networks with the repulsive interactions “real” networks.

The aim of the present study is to verify the validity of Kuhn’s affine network theory of rubber elasticity by CGMD simulations. The elastic behavior of phantom ideal networks is clarified by both CGMD simulations and a modified version of Kuhn’s theory. In addition, we investigated how the EV interactions between individual partial chains influence the elastic properties of the network. The remainder of this paper is organized as follows. In Section II, the coarse-grained model of the ideal networks is explained. In Section III, the stress–extension curve and shear modulus obtained from a series of simulations are discussed. Finally, we give a brief conclusion in Section IV.

First, we will review coarse-graining of cross-linked networks. The bead–spring model proposed by Kremer and Grest^35,36 successfully reproduces several phenomena of polymer materials, both static and dynamic. We used or modified their model to represent the ideal networks. As shown in Fig. 1, the partial chains and cross-linking points, which are the main elements of networks, are modeled with particles of the same size connected by springs. The ideal networks are constructed with two types of interactions: a bonding potential U_b(r) and a non-bonding potential U_nb(r), where r is the interparticle distance.

The bonding potential is designed to represent the bond between two connected particles and can be split into two terms:

The first term U_LJ(r) is the purely repulsive Lennard-Jones (LJ) potential:^37,38

where σ and ϵ are the diameter and interaction strength of the LJ particles, respectively, and r_c is the cutoff length set to 2^1/6σ. The second term U_FENE(r) is the finitely extensible nonlinear elastic (FENE) potential:³⁹ 

where k is the elastic coefficient and R₀ is the fully extended length of the bond.

The non-bonding potential U_nb(r) is chosen according to whether the network is phantom or real. As explained in Section I, the former are hypothetical networks without EV interactions, while the latter are realistic networks with EV interactions. It is technically easy to implement the two types of networks. When we want to use the phantom networks, we only have to turn off the non-bonding potential, i.e., U_nb = 0. In contrast, when we use the real networks, the non-bonding potential is defined as the repulsive LJ potential [Eq. (4)], i.e., U_nb(r) = U_LJ(r).

In the present study, we consider four types of regularly ordered networks. The unit cells are schematically shown in Fig. 2. From a structural point of view, these crystalline-like networks are characterized by three variables: the length of the bridged chains l, their number density ν, and the functionality of the cross-linking points f. For simplicity, we assume that the length of the partial chains is mono-dispersed and defined by the number of particles they contain (except for the particles belonging to the cross-linking points). The networks are classified by the f value: “gyroid” (f = 3), “diamond” (f = 4), “cube” (f = 6), and “octopus” (f = 8).

We will now discuss the feasibility of the ideal networks. The main features of our networks are as follows:

1. All of the partial chains are elastically effective bridged chains with the same length.

2. All of the cross-linking points have the same functionality and are homogeneously allocated.

3. There is no trapped entanglement between any partial chain, unless special procedures are added.

Strictly speaking, our networks are orthotropic rather than isotropic bodies because the alignment of each constituent chain is isotropic but extremely discrete. They are qualitatively different from actual vulcanized rubbers. As a result, it is expected that their elastic properties, such as their nominal stress and shear modulus, strongly depend on the elongational direction. Nonetheless, all of the features listed above are compatible with the assumptions for the ideal network on which Kuhn’s theory is based. Thus, we treat the crystalline-like networks in Fig. 2 as a type of ideal network, even at the expense of the discrete orientation of the partial chains. It is desirable that we should extract the universal properties of the networks independent of the elongational direction. However, this is beyond the scope of the present paper and will be discussed in the future.

Here, we discuss the relationship among the three quantities l, ν, and f. In general, these three quantities are closely related. For the ideal networks shown in Fig. 2, the following relation is established:

where ρ is the number density of the constituent particles. When l ≫ 1, Eq. (6) can be approximately represented by ν ≃ ρ/l, which means that once l and ρ are chosen, ν is automatically decided and almost independent of f. In the present study, the approximation partially holds.

To obtain the stress–extension curves of the ideal networks, we performed CGMD simulations of uniaxial elongation of the networks. It is assumed that each particle has the same mass m and obeys usual Langevin dynamics under the NVT ensemble.^25,36 Hereafter, we set the quantities σ, ϵ, and m to the units of length, energy, and mass, respectively. Time, temperature, and stress (including pressure and shear modulus) are measured in units of $τ≡(mσ2/ϵ)1/2$⁠, ϵ/k_B, and ϵ/σ³, respectively. In the next sections, for simplicity, the units are omitted, except in the figures. In the simulations, the usual periodic boundary conditions are applied to a simulation box in the three axis directions. The main parameters are adjusted as follows: the total number of constituent particles of the ideal networks N is in the range 17, 496–31, 104. The repetition number of unit cells in each axis direction M is the same and set to greater than or equal to 3. It has been confirmed that networks composed of such repeated unit cells do not show significant finite size effects.^29,31 Therefore, we conclude that the wide range of N matters nothing. The details of the molecular information about the ideal networks are given in supplementary material. The particle number density is set to ρ = 0.85 σ⁻³, which corresponds to that of actual polymer melts.³⁶ The temperature is kept constant at T = 1.0 ϵ/k_B. The time increment used to solve the set of Langevin equations is fixed at 0.01 τ. For the two parameters included in the potential U_FENE(r) [Eq. (5)], we set k = 30.0 ϵ/σ² and R₀ = 1.5 σ to prevent artificial crossability between partial chains.^35,36 The CGMD calculations were performed with COGNAC⁴⁰ in OCTA⁴¹ and LAMMPS.⁴² 

We will now explain how to calculate the uniaxial elongation process of the ideal networks. Before deformation, we first have to prepare the networks equilibrated at ρ = 0.85 σ⁻³. For the initial configuration of the networks, we arrange the particles belonging to each chain at equal intervals of about 1 σ along the straight line bridging two nearest-neighboring cross-linking points, which results in the crystalline-like networks shown in Fig. 2. At this stage, the particle densities are much lower than ρ = 0.85 σ⁻³. Thus, isotropic compression of the networks is performed up to the desired density slowly enough not to cause artificial intersection of chains. Equilibration calculations are then performed for at least 10⁴ τ. During uniaxial deformation, the elongation rate is kept at 10⁻³ σ/τ, the inverse of which is much longer than the terminal relaxation time of each partial chain. The networks are also deformed with their volume conserved (Poisson’s ratio = 0.5), which is realized by equivalently compressing the networks in the two directions perpendicular to the elongational direction. When the network is elongated in the z-axis direction, the normal true stress of the network σ_tr is defined with the diagonal elements of the true stress tensor σ_αα (α = x, y, z) in the following way:

where subtraction of the second term is performed to remove the influence of the hydrostatic pressure determined by the external conditions.^1,43 Moreover, σ_tr is converted into its nominal version σ_n by the relation σ_n = σ_tr/λ. The converted stress corresponds to the experimentally obtained stress. Finally, we obtain the shear modulus G of the ideal network by fitting the approximate relation [Eq. (1)] to our result for the stress–extension relation.

First, the stress–extension curves, which are also called stress–strain (SS) curves, are discussed. Figure 3 shows the overall SS curves for the four ideal networks with the same partial chain length (l = 13). The solid lines correspond to the phantom networks without EV interactions, while the dashed lines correspond to the real networks with EV interactions. The SS curves for the randomly cross-linked networks composed of partial chains with number-averaged length l_n = 13.8 are also shown for comparison. The disordered networks are prepared in the following way. The initial system is composed of a single long chain with length 20, 000. The particle density is kept constant at ρ = 0.85 during preparation. After sufficient equilibration of the chain, we replace particles along the chain with reactive particles at an equal interval of 13. The number of allocated reactive particles is 1, 400. By connecting these reactive particles pairwise, a randomly cross-linked network with practically no dangling chains is obtained with a cross-linking yield of 98.1 %. There are no significant differences between the ideally and randomly cross-linked networks, except that each cross-linking point of the former consists of one particle while that of the latter consists of two particles.

From Fig. 3, the extensibility of the networks is strongly dependent on their type (phantom or real). Let us consider the situation where the networks have uniaxially elongated to their limit. Using the number of extended partial chains spanned in the elongational direction, the edge length of the undeformed simulation box L (see supplementary material for specific numerical values), and the extension limit of the bond in the FENE potential R₀ = 1.5 [Eq. (5)], we can estimate the extension ratio limits for the ideal networks λ_max. For l = 13, the limit values are λ_max(gyroid) = 14.5, λ_max(diamond) = 13.3, λ_max(cube) = 5.8, and λ_max(octopus) = 8.4. By sorting the values, the following order is obtained: λ_max(cube) < λ_max(octopus) < λ_max(diamond) < λ_max(gyroid). This order exactly agrees with the actual order of the extensibility of the networks.

Here, we will discuss the SS curve for the randomly cross-linked network and compare it with that for the diamond network. Although they have the same functionality (f = 4), the latter extends much longer than the former. This reflects the difference in the distributions of cross-linking points: the distribution of cross-linking points in the random cross-linked network is heterogeneous, while that in the diamond network is both homogeneous and regular (Fig. 2). The difference in the extensibility can also be explained by the minimum path in the network in the extension direction. The minimum path in the randomly cross-linked network with l_n = 13.8 is composed of 117 particles. Conversely, that in the diamond network with l = 13 is composed of 280 particles. In addition, the edge lengths of the cubic simulation boxes in the undeformed state are 28.7 for the former and 31.7 for the latter. Similar to the preceding discussion about the four types of ideal networks, the estimated λ_max values for the two networks are λ_max(random) = 6.1 for the randomly cross-linked network and λ_max(diamond) = 13.2 for the diamond network. Figure 3 shows that the latter is twice or more as extensible as the former for both the phantom and real networks.

Finally, we will discuss the effect of the EV interactions inside the networks on their elastic properties. From Fig. 3, the nominal stress of the phantom ideal network σ_n is larger than that of the corresponding real network over a wide range of extension ratios λ. Similar behavior is observed for the randomly cross-linked network, although the range is much narrower. The four types of ideal networks contain no structural defects, such as trapped entanglements, dangling chains, and loops.^1,7 Even in the randomly cross-linked network, the effect of trapped entanglements on the elasticity of the network is probably negligible because the number-averaged partial chain length (l_n = 13.8) is shorter than the entanglement length (35–70) estimated from previous CGMD simulations of linear polymer melts.^36,44 Thus, the higher values of σ_n for the phantom networks are mainly because of the lack of EV interactions inside the networks. It is surprising that the EV interactions affect the mechanical properties of the networks over a wide range of λ. From Eq. (1), it is clear that the shear modulus of the phantom network is also larger than that of the corresponding real network, which has not been previously reported. In Section III D, the effect of the EV interactions on G will be discussed in more detail.

As mentioned in Section II A, the partial chains in our ideal networks are oriented in a finite number of discrete directions. In addition, from the radial distribution functions of the partial chains in undeformed states (not shown in this paper), they are always non-Gaussian. The latter fact indicates that of the several assumptions in Kuhn’s theory,^1,3 the Gaussian assumption that each partial chain with no external force obeys Gaussian statistics is not satisfied for our networks. Without the Gaussian assumption, Tobolsky^45,46 and Flory⁴⁷ derived a generalized version of Eq. (2):

where ⟨R²⟩ is the mean square end-to-end distance of the mono-dispersed bridged chains in the undeformed network and $⟨R2⟩0$ is that of free chains of the same length with no EV. The ratio of $⟨R2⟩/⟨R2⟩0$ is often referred to as the “front factor” and it is sensitive to the process of cross-linking. When each partial chain satisfies the Gaussian assumption, the factor becomes unity and Eq. (8) reduces to Eq. (2). The front factors $⟨R2⟩/⟨R2⟩0$ for the four types of ideal networks can be expressed as

where b_s is the size of the particles constituting the partial chains. This is an adjustable parameter and it is set to b_s = r_min, where r_min is the interparticle distance that gives the minimum value of U_b(r) [Eq. (3)]. Derivation of the front factors $⟨R2⟩/⟨R2⟩0$ is given in supplementary material.

The dependence of the shear moduli G of the four types of phantom ideal networks on the partial chain density is shown in Fig. 4. From Fig. 4, a set of data extracted from the SS curves agree well with the prediction [Eq. (8)] by the modified theory, regardless of the type of ideal network. The compatibility with the modified theory is far better than that with the original Kuhn theory. For networks with high functionality, such as cube (f = 6) and octopus (f = 8) networks, a few data in a high-density region of bridged chains show slight deviation above the theoretical curve. The origin of the deviation will be discussed in Section III D.

Here, we will discuss the relation between G and f for the phantom ideal networks. Figure 4 indicates that G for our networks monotonically increases with increasing f under the condition of the same ν. We also expect that G is related not only to f but also to the type of network. However, a more detailed relationship remains unclear.

Here, we compare the shear moduli G of the real ideal networks with those of the corresponding phantom networks and then discuss the effect of EV interactions on G. Figure 5 shows the dependence of G on the partial chain density ν for the four types of networks. Clearly, G of the real network, as well as that of the corresponding phantom network, increases with increasing ν, but the former is always smaller than the latter except in the high-density region of bridged chains in the networks with high functionality, i.e., cube (f = 6) and octopus (f = 8).

We will interpret the magnitude relation of G from the viewpoint of the microscopic origin of the elasticity of the cross-linked networks. First, our ideal networks, both phantom and real, have no trapped entanglement. In addition, the elasticity that all the ideal networks show under small uniaxial elongation is purely entropic, except for the networks composed of extremely short partial chains and cross-linking junctions with high functionality (see supplementary material for details). The tension under uniaxial deformation along the z-axis f_z can be expressed as^1,7,8

where S is the entropy of the whole system in the deformed state, V is the volume of the network, and L_z is the edge length of the deformed network in the elongational direction. The nominal stress σ_n corresponding to f_z is given by

where L is the edge length of the network before deformation.

Here, we consider the configuration space available to each partial chain.^32,34 To distinguish the quantities associated with the phantom and real networks, we will hereafter use the subscripts “ph” and “real”, respectively. Let us assume that both ends of the chains are located with their relative positions kept constant and with no fluctuation, which is the same condition as the original Kuhn theory.^1,3 If all the constituent particles have their own EV, the space where each chain can move freely becomes smaller than that of the corresponding chain without EV interactions.³² A smaller accessible space of each real chain finally leads to smaller conformational entropy of the whole network, which leads to S_ph > S_real. In addition to the relation, the smoothly monotonic decrease in S as a function of L_z and approximately the same extension limit L_z,max in the two type of networks probably results in ∂S_ph/∂L_z > ∂S_real/∂L_z as schematically shown in Fig. 6. From Eqs. (13) and (14), it is concluded that σ_n of the real network under small uniaxial deformation becomes smaller than that of the corresponding phantom network (σ_n,ph > σ_n,real). From Eq. (1), G of the former also becomes smaller than that of the latter (G_ph > G_real), which is consistent with most of the results in Fig. 5. We expect that simulations of coarse-grained networks without EV interactions will be very effective not only to confirm the validity of various molecular theories for rubber elasticity, but also to directly investigate the effect of EV interactions on the elasticity of networks with any structure.

Finally, we will discuss networks composed of extremely short partial chains. The small number of conformations available to the short chains leaves each chain less Gaussian and then the elasticity of the networks less entropic. In addition, the EV interactions between the short chains in the real networks probably increase the energetic contribution to the elasticity. The inversion phenomena of G in Fig. 5 can be mainly attributed to the qualitative change in the origin of the elasticity. Equations (3)–(5) also indicate that for both the phantom and real networks, EV interactions exist even between bonded particles, which can cause the networks with extremely short chain lengths to show somewhat artificial energetic elasticity. The deviation of G above the theoretical prediction, which is discussed in Section III C, will be partly driven by the artificial elasticity of the phantom networks.

We have performed CGMD simulations for four types of realizable ideal networks (gyroid, diamond, cube, and octopus networks) to verify the validity of Kuhn’s affine network theory. We clarified the relation between the structural and elastic properties of both the phantom and real versions of the networks. The main conclusions are as follows:

1. The four types of ideal networks, whether phantom or real, show extremely high extensibility, which increases in the following order: cube < octopus < diamond < gyroid. The order can be explained by the edge length of the unstrained simulation box and both the length and number of partial chains spanned in the elongational direction. For the same functionality and a similar chain length, the ideal network shows higher elongation than the randomly cross-linked network. This is because of the homogeneous and regular arrangement of the cross-linking points in the ideal networks.

2. The dependence of the shear modulus of the phantom ideal networks on the partial chain density is in good agreement with the prediction from a generalized network theory without the Gaussian assumption. This means that the crystalline-like ideal networks are incompatible with the Gaussian assumption. Thus, to compare various types of cross-linked networks, the generalized Kuhn theory should be used as the standard tool.

3. The shear modulus for the real ideal network is smaller than that for the corresponding phantom network. This is because of the difference in the conformational entropy of each partial chain inside the network. From the size of the space available to the chain, the conformational entropy of each chain in the real network is much smaller than that in the phantom network. The decrease in the conformational entropy results in decreases in both the nominal stress and shear modulus of the network.

We believe that CGMD simulations of networks without EV is a promising method not only to verify the validity of molecular theories for rubber elasticity, but also to directly investigate the EV effect on the elasticity of networks.

See supplementary material for the details of the molecular information about the four types of ideal networks and the temperature dependence of the shear moduli G and the derivation of the front factors $⟨R2⟩/⟨R2⟩0$ of these networks.

This study was partly performed at the rubber and elastomer consortium hosted by AIST. We appreciate fruitful discussions with all of the members. We would also like to thank Prof. T. Kawakatsu, Prof. T. Indei, and Prof. K. Nakajima for valuable comments.