**Preface to "Elastomers: From Theory to Applications"**

This Special Issue "Elastomers: From Theory to Applications" focuses on the current state of the art of elastomers, both in modern applications and from a theoretical perspective. The main characteristic of elastomer materials is the high elongation and (entropy) elasticity of these materials, and the ability to swell multiple times in a suitable solvent. The use of filled elastomers, especially of new kinds of elastomer nanocomposites, is of high interest for rubber technologies. Nanostructured adaptable gels allow for the tailoring of responsive materials or filtering systems. These materials enable widespread applications in engineering fields, ranging from modern tire technologies to medical applications and consumer goods. Elastomers are also useful in various biomaterial applications. Bio-elastomers are widely available in nature and have been shown to have specific properties that are often far superior to their synthetic counterparts. Furthermore, nowadays, additional functions play a major role in elastomeric composites, e.g., in the use of dielectric, electrorheological, or magnetorheological elastomers for smart applications in soft robotics, etc.

All elastomers share typical features, such as entropy-driven elasticity, the presence of entanglements, and topological constraints of network chain conformations. These features still offer fascinating scientific challenges in the synthesis, characterization and application, as well as for the theory of polymer networks and the cross-scale modeling of elastomeric materials, their polymeric constituents and, finally, the corresponding solids under real-life conditions.

The series of papers within this Special Issue offer the latest research in several specific sub-areas of elastomer research or review selected fields. The scope of the Special Issue encompasses the leading scientific contributions in the synthesis, characterization, and theory of elastomers. Of particular interest are new structures and functionalities incorporated into elastomers, leading to enhanced properties of crosslinked elastomeric materials, and/or to a better understanding of the structure–property relationships and application behaviour.

Thus, the present volume provides a comprehensive overview of the recent developments in the field of elastomers and will be of interest to both academic researchers and industrial professionals.

G. H. acknowledges the DFG (German Research Foundation) Project 380321452/GRK2430 for financial support.

> **Gert Heinrich and Michael Lang** *Editors*

## *Review* **Theory of Flexible Polymer Networks: Elasticity and Heterogeneities**

#### **Sergey Panyukov**

P. N. Lebedev Physics Institute, Russian Academy of Sciences, Moscow 117924, Russia; panyukov@lpi.ru Received: 28 February 2020; Accepted: 21 March 2020; Published: 1 April 2020

**Abstract:** A review of the main elasticity models of flexible polymer networks is presented. Classical models of phantom networks suggest that the networks have a tree-like structure. The conformations of their strands are described by the model of a combined chain, which consists of the network strand and two virtual chains attached to its ends. The distribution of lengths of virtual chains in real polydisperse networks is calculated using the results of the presented replica model of polymer networks. This model describes actual networks having strongly overlapping and interconnected loops of finite sizes. The conformations of their strands are characterized by the generalized combined chain model. The model of a sliding tube is represented, which describes the general anisotropic deformations of an entangled network in the melt. I propose a generalization of this model to describe the crossover between the entangled and phantom regimes of a swollen network. The obtained dependence of the Mooney-Rivlin parameters *C*<sup>1</sup> and *C*<sup>2</sup> on the polymer volume fraction is in agreement with experiments. The main results of the theory of heterogeneities in polymer networks are also discussed.

**Keywords:** elastomers; polymer networks; elastic modulus; loops; heterogeneities

#### **1. Introduction**

Polymer networks and gels belong to a unique class of materials that have the properties of both solids and liquids. Such "soft solids" are elastically deformed at macroscopic scales, while at short distances and short times the network strands experience liquid fluctuations, which are responsible for the exceptional properties of polymer networks, such as the reversibility of huge deformations (when stretched up to 3000%, see [1]).

In this work, I review the main microscopic (molecular) approaches to the description of the elasticity of flexible polymer networks and establish relationships between these approaches. Phantom networks are considered in Section 2, and networks with topological entanglements are studied in Section 3. Spatial heterogeneities developed in swollen and deformed polymer networks are investigated in Section 4.

The classical theory proposed by Flory more than half a century ago, Ref. [2] suggests that polymer networks have a tree-like structure. Although this theory takes into account the presence of loops in the network, it is assumed that they all are infinite in size. Such an approach is insufficient to describe the thermodynamics of polymer networks with loops of finite size, and explicit consideration of these loops is required. The mathematical model of such networks is based on the replica method. This approach takes into consideration that the properties of the network depend not only on the "current" conditions in which the network is placed but also on the conditions of its preparation since it is these conditions that determine the molecular structure of the formed network. The replica approach takes into account the excluded volume interaction of linear chains both at preparation and experiment conditions by generalizing the theory of the "*<sup>n</sup>* → 0-component" field *<sup>ϕ</sup>*<sup>4</sup> proposed by de Gennes. The main properties of the order parameter *ϕ* introduced for the description of such "soft solids" are discussed in Section 2.

The conformations of the chains in the network are described by the combined chain model [3]. This model is based on the concept of virtual chains, which determine the effective elasticity of the tree-like structures in the network. In polydisperse networks, the distribution of virtual chains is also polydisperse. I show that the order parameter in the replica network model stores information about this distribution. The found solution of the replica model is used to calculate the distribution function of the lengths of virtual chains in polydisperse networks.

The structure of real networks differs significantly from the classical picture of ideal trees. Typical loops of such networks have finite dimensions and strongly overlap with each other. The impact of such loops on the conformations of the network chain is described by the extended combined chain model [4]. I discuss the connection of this model with the predictions of the replica theory, which offers an analytical description of such mutually overlapping loops. Along with typical loops, the network may contain topological defects, which are also taken into account by the replica method. The condition for the loss of elasticity by polymer networks, both due to the presence of topological defects and near the gel point is discussed.

With an increase in the length of polymer chains, the effects of topological entanglements become significant. The properties of entangled polymer liquids are significantly different from the properties of entangled soft solids—polymer networks. While non-concatenated rings in the melt are compacted into fractal loopy globules [5], internal stresses in stretched polymer networks prevent such a collapse. Entangled networks are described by the non-affine tube model that generalizes the concept of virtual chains to the case of entangled polymer systems [6]. I discuss the physics of deformation of network strands in this model in Section 3.

Anisotropic deformations of the entangled network lead to slippage of chains along the contour of the entangled tube, which is described by the slip tube model [3]. This model suggests the additivity of the phantom and entangled contributions to the free energy of the network. Although this additivity approximation describes well the uniaxial deformations of the network in the melt, it cannot be used to describe the crossover of the Mooney Rivlin dependence between entangled and phantom deformation regimes observed during swelling of the polymer network. I propose the generalization of the slip tube model to account for this crossover and show that the obtained concentration dependences of the Mooney-Rivlin parameters *C*<sup>1</sup> and *C*<sup>2</sup> are in agreement with the experimental data.

The irregularity of the molecular structure of polymer networks leads to the presence of spatial inhomogeneities. The main results of the theory of heterogeneities in polymer networks are discussed in Section 4.

#### **2. Phantom Networks**

The polymer network is obtained by crosslinking linear chains. In models of phantom polymer networks, it is assumed that in the process of thermodynamic fluctuations, polymer chains can freely "pass" through each other. The phantom approximation works well only for not too long network strands; otherwise, topological restrictions related to the mutual impermeability of polymer chains should be taken into account. In this section, I discuss the main approaches to the theoretical description of phantom networks. The elastic free energy of a network of Gaussian chains is quadratic in the deformation ratios *λα* = *Lα*/*L*0*<sup>α</sup>* of the network dimensions along the principal axes of deformation *α* = *x*, *y*, *z* in the deformed state (*Lα*) and under conditions of its preparation (*L*0*α*):

$$\mathcal{F}\_{ph} = G \sum\_{\alpha} \frac{\lambda\_{\alpha}^{2} - 1}{2} \tag{1}$$

The interaction of strand monomers is usually taken into account by imposing the incompressibility condition

$$
\lambda\_x \lambda\_y \lambda\_z = 1/\phi. \tag{2}
$$

where *φ* is the polymer volume fraction. In the case of uniaxial deformation we have

$$
\lambda\_{\mathbf{x}} = \phi^{-1/3} \lambda\_{\mathbf{r}} \quad \lambda\_{\mathbf{y}} = \lambda\_{\mathbf{z}} = \phi^{-1/3} \lambda^{-1/2} . \tag{3}
$$

In the affine network model, it is assumed that the ends of each network strand are pinned to a "non-fluctuating elastic background" which deforms affinely with the network surface [3]. The elastic modulus of such a network is proportional to the density *ν* of its strands,

$$G = kT\nu\tag{4}$$

and *kT* is the thermal energy. Equation (4) is attractive for its simplicity. However, in real networks, the strand ends are not fixed in space—they are connected to neighboring strands through cross-links. The elastic modulus *G* depends on the molecular structure of the network and the polydispersity of its strands. At first glance, the "bookkeeping" of the order of connection of all network strands with each other and the numbers of monomers of each of these strands seems impossible for networks of macroscopic sizes. However, this is precisely what the replica method does.

#### *2.1. Replica Method*

The replica method was proposed a long time ago [7] and was used later to calculate correlation functions of density fluctuations in polymer networks [8]. The main goal of this approach is to build a mathematical model that takes into account most of effects inherent to real polymer networks. Unlike widespread numerical simulations, which can also take these effects into account, the replica method is an analytical theory, which is more amenable for understanding network elasticity as a function of many different molecular features.

When the polymer network is deformed, its structure remains unchanged, formed under the conditions of its preparation. The main idea of the replica approach is to consider the polymer network in an expanded space with coordinates **xˆ** = (**x**(0), **x**(1), ... , **x**(*m*)) of dimension 3(1 + *m*). The replica space consists of an "initial system" under conditions of the network preparation (with coordinates **x**(0)) and *m* replicas of the "final system" under experimental conditions (with coordinates **x**(*k*), *k* = 1, ... , *m*), see Figure 1.

**Figure 1.** The replica system consists of the initial system (*k* = 0)—the polymer network at preparation conditions (blue cube) and *m* identical replicas of the final system (*k* = 1, ... , *m*)—the network, deformed by factors *λα* along the principal axes *α* = *x*, *y*, *z* of deformation (yellow cuboids). The polymer network in the replica space of dimension 3(1 + *m*) is mainly localized inside the (green) cylinder with the diameter *<sup>R</sup> aN*¯ 1/2 directed along unit vectors {**eˆ** *<sup>α</sup>*}, corresponding to affine deformation of the polymer. The networks (and all their monomers and strands) in the initial system and any of the *m* replicas of the final system are the projections of the network in the replica space onto the corresponding subspaces *k* = 0, . . . , *m*.

In each of these systems, the polymer network with all of its monomers and bonds is the projection of the network in the replica space onto the corresponding subspaces. Therefore, in all these systems, the network structure is exactly the same, while the conformations of network strands can vary significantly. The free energy *F* of the final system (a deformed network) is expressed through the analytic continuation to *m* = 0 of the free energy *Fm* of the replica system [8,9]:

$$F = \lim\_{m \to 0} \frac{dF\_m}{dm}.\tag{5}$$

#### *2.2. Liquid-Solid Order Parameter*

According to de Gennes, polymer networks are "soft solids", which emphasizes the presence of liquid and solid-state degrees of freedom interacting with each other. In the Landau approach, solids are described by an order parameter—a number (in general case, a tensor), equal to zero for a liquid and nonzero only in the solid phase. Unlike ordinary low molecular weight solids, the properties of polymer networks depend on the characteristics of both the initial and final systems. Therefore, the order parameter for polymer networks is not a number, but a function of the difference of coordinates *x* (*k*) *<sup>α</sup>* − *λαx* (0) *<sup>α</sup>* of final and initial systems. This dependence reflects the dual nature of polymer networks: At scales larger than the characteristic size of the network cycles, the network deforms affinely as an ordinary solid, see Figure 1. On a smaller spatial scale, the network chains fluctuate in space and have conformations similar to the chains in a polymer liquid. The unique properties of polymer networks are determined by the interaction of their solid-state and liquid degrees of freedom [10].

In the liquid phase, the order parameter is independent of coordinates. In the case of Gaussian polymer networks, the order parameter *ϕ* (*ς*) is a function of a single variable

$$\boldsymbol{\xi} = \frac{1}{2} \left[ \hat{\mathbf{x}}^2 - \sum\_{\boldsymbol{\alpha}} \left( \hat{\mathbf{e}}\_{\boldsymbol{\alpha}} \hat{\mathbf{x}} \right)^2 \right],\tag{6}$$

where **eˆ** *α* is a unit vector along the direction *x* (*k*) *<sup>α</sup>* = *λαx* (0) *<sup>α</sup>* of the affine deformation in the replica space. At *ς* = 0 (that is, at *x* (*k*) *<sup>α</sup>* = *λαx* (0) *<sup>α</sup>* ), the order parameter *ϕ* (0) is determined only by the conditions of the network preparation and does not depend on the experimental conditions. This last dependence is encrypted in the dimensionless order parameter *χ* (*t*) = *ϕ* (*ς*) /*ϕ* (0) which is the function of the dimensionless variable *t* = *ς*/*R*. In the case of well-developed networks obtained by crosslinking polymer chains with *f*-functional monomers far beyond the gel point, the characteristic radius is *<sup>R</sup> bN*¯ 1/2, see Figure 1. Here *<sup>b</sup>* is the monomer size and *<sup>N</sup>*¯ is the average number of network strand monomers.

In the case of monodisperse networks, the function *χ* (*t*) was calculated by Edwards (see Equation (3.19) of Ref. [11]):

$$\chi\left(t\right) = e^{-\frac{f-2}{f-1}t} \tag{7}$$

In the case of polydisperse networks, this function *χ* (*t*) is determined by the differential equation (see Appendix A)

$$t\chi'''(t) = \chi\left(t\right) - \chi^{f-1}\left(t\right) \tag{8}$$

For large *t* 1, the function *χ* (*t*) decreases more slowly compared to the monodisperse case (7), as a stretched exponent

$$
\chi\left(t\right) \sim t^{1/4} e^{-2t^{1/2}}\tag{9}
$$

For small *t*, the solution of Equation (8) can be approximated by a power function

$$\chi\left(t\right) \simeq \left(1 + ct\right)^{-k} \tag{10}$$

Substituting it in Equation (8) and equating to zero the first two expansion coefficients (*t* and *t* 2) in powers of *t*, we find

$$k = \frac{3}{f - 2}, \quad c = \frac{(f - 2)^2}{f + 1} \tag{11}$$

The replica method is also called "replica trick" because of the hidden meaning of its mathematical constructions. As will be shown below, Equation (8) describes the relationship between distribution functions of the molecular trees that characterize the structure of the polymer network. Note that Equation (8) is obtained by estimating the Hamiltonian (A2) in the Appendix A using the steepest descent method. Therefore, it corresponds to the mean-field approximation, which works due to the presence of a large number of overlapping network strands.

#### *2.3. Overlap Parameter*

Like in liquid polymer systems, in a typical network there are many overlapping network strands. An important characteristic of the network is the overlap parameter, which is defined as the number of network strands within the volume *R*<sup>3</sup> pervaded by one network strand. In the case of a solution of chains with monomer density *<sup>ρ</sup>*(0), *<sup>P</sup>*(0) *<sup>ρ</sup>*(0)*R*3/*N*. In Gaussian networks, the size of the network strand with Kuhn length *<sup>b</sup>* is *<sup>R</sup> bN*1/2, and

$$P^{(0)} \simeq \rho^{(0)} b^3 N^{1/2} \tag{12}$$

with a large parameter *P*(0), polymer networks obtained far from the gel point can be considered as *P*(0) overlapping yet independent "elementary" networks, the strands of which consist of *N* monomers. There is on the order of one network strand per volume *R*<sup>3</sup> in the "elementary" network, and its elastic modulus is estimated as *Gelem kT*/*R*3, *kT* is the thermal energy. Since the real network consists of *P*(0) overlapping elementary networks, its elastic modulus is *P*(0)-times larger [12]:

$$G = P^{(0)} G\_{elem} \simeq kT \rho^{(0)} / N \tag{13}$$

The polymer network can also be represented as *P*(0) polymer trees (or "layers") overlapping with each other. The perfect network model assumes an infinite number of layers, and also that these trees have an infinite number of generations. In real networks *P*(0) < *Pe* is finite and not too large (the entanglement overlap parameter *Pe* ∼ 20–30, see Section 3 below), and the different root trees are interconnected by loops. The number of network strands forming a closed loop of minimum length, which binds together different layers of the network, depends logarithmically on the overlap parameter [4],

$$d \approx \frac{1}{f - 1} \ln P^{(0)} \tag{14}$$

The case *<sup>P</sup>*(0) <sup>≈</sup> 1 with the loop size *<sup>l</sup>* <sup>≈</sup> 1 describes the "elementary networks" in Equation (13). In accordance with the results of numerical simulations [13], the loop length *l* decreases with dilution and with increasing functionality *f* of cross-links. Note that the minimal size loops do not directly determine the network elasticity since each elastically effective chain is simultaneously part of a large number of loops. The cumulative effect of these typical loops of the network is self-averaged, and can be described in the effective mean-field approximation, see Equation (37) below.

The overlap parameter also determines the density fluctuation of the network obtained in the case of incomplete conversion above the gel point, at *p* > *pc*. In the mean-field approximation, the density of gel monomers is

$$
\rho\_{\mathcal{S}}^{(0)} \simeq (p - p\_{\mathcal{C}}) \,\rho^{(0)}\,,\tag{15}
$$

and its connectivity radius is estimated as *<sup>ξ</sup> <sup>b</sup>* [*N*/ (*<sup>p</sup>* <sup>−</sup> *pc*)]1/2. Fluctuations in the number of monomers of the polymer network on the scale of the connectivity radius occur by attaching or detaching sol molecules with a characteristic number of monomers

$$L \simeq \mathcal{N} / \left(p - p\_c\right)^2\tag{16}$$

Therefore, relative fluctuations in the gel density are estimated as

$$\frac{\delta \rho\_{\mathcal{S}}^{(0)}}{\rho\_{\mathcal{S}}^{(0)}} \simeq \frac{L}{R\_{bb}^3 \rho\_{\mathcal{S}}^{(0)}} \simeq \frac{1}{(p - p\_c)^{3/2} \, P^{(0)}},\tag{17}$$

In the case of a large overlap parameter, these fluctuations are small outside the narrow Ginzburg region,

$$p - p\_{\mathcal{E}} \gg \left( P^{(0)} \right)^{-2/3} \tag{18}$$

#### *2.4. Combined Chain Model: Monodisperse Networks*

The elasticity of the polymer network has an entropic origin and is determined by the change in strand conformations during the network deformation. The conformations of a strand depend on its interaction with the elastic environment, which can be described by virtual chains attached to the ends of this strand. The network strand together with the two virtual chains is called a combined chain. The ends of this combined chain are attached to an affine deformable non-fluctuating background, see Figure 2b.

**Figure 2.** (**a**) The model of perfect network with tree-like connection of network strands (solid curves) and infinite size loops (dotted lines). (**b**) A network strand with *N* monomers can be considered as a part of the combined chain, connected through two effective chains with *n* monomers (dashed green curves) to the non-fluctuating background. (**c**) Actual network with finite size loops. (**d**) The combined chain consists of three effective chains. The filled circles and crosses indicate cross-links and nonfluctuating background, respectively. The end-to-end vector of the combined chain **X** deforms affinely with network deformation.

In the perfect network model, it is assumed that the network has a tree structure, whereas all its cycles are of infinite size, see Figure 2a. Each of the trees attached to the ends of the network strand is modeled by a virtual chain. The effective number *n* of its monomers is related to the numbers *ni* of virtual chain monomers on the next generation of the tree *i* = 1, ... , *f* − 1 and the corresponding numbers of strand monomers *Ni* as

$$\frac{1}{n} = \sum\_{i=1}^{f-1} \frac{1}{m\_i}, \qquad m\_i = N\_i + n\_i \tag{19}$$

In the case of a monodisperse network with a fixed number *Ni* = *N* of monomers of its strands, all *ni* = *n* and the solution of Equation (19) determines the number of monomers of the virtual chain [14]

$$m = N / \left(f - \mathfrak{Z}\right) \tag{20}$$

The average vector **R** between the ends of the network strand is determined from the force balance condition. In the case of a monodisperse network,

$$\mathbf{R} = \frac{\mathbf{x}}{1 + 2n/N} \tag{21}$$

where *n* is the number of monomers of the virtual chain, Equation (20). The mean square fluctuation of the vector **R** is determined by the parallel connection of the real chain and two serially-connected virtual chains

$$\left< \left( \Delta \mathbf{R} \right)^{2} \right> = \frac{b^{2}}{1/N + 1/\left< 2n \right>} = \frac{2}{f} b^{2} N \tag{22}$$

The vector **X** between the ends of the combined chain is deformed affinely with macroscopic network deformation. The elastic modulus of the perfect network with monomer density *ρ*(0) and chain concentration *ρ*(0)/*N* is

$$G = \frac{kT\rho^{(0)}/N}{1 + 2n\_1/N} = kT\nu \left(1 - \frac{2}{f}\right) \tag{23}$$

Here *ν* = *ρ*(0)/*N* is the concentration of strands and (1 + 2*n*1/*N*) <sup>−</sup><sup>1</sup> is the strand fraction of the combined chain. A comparison of expressions (23) and (4) shows that the elastic modulus of the perfect network is less than the modulus of the affine network by a factor of 1 − 2/ *f* . It is generally accepted that the 2/ *f* factor is associated with the presence of fluctuations in the strand size, which are also proportional to this factor (see Equation (22)). In fact, fluctuations have nothing to do with it [15]: an exact elastic modulus of Gaussian phantom network is different from that of the perfect network model and coincides with the elastic modulus of the classical non-fluctuating grid, in which each strand is replaced by a corresponding elastic thread with the same elastic stiffness coefficient [16]. This coincidence is due to affine deformation of the average distances between the cross-links in such networks.

Equation (23) can be recast in more general form with the network modulus proportional to the difference of the number densities of network strands *ν* and cross-links *μ* = 2*ν*/ *f* , since there are *f* /2 strands per crosslink [14]:

$$G = kT\left(\mu - \nu\right) \tag{24}$$

This expression is quite universal; with general densities *ν* and *μ*, it describes polydisperse networks, as well as networks with incomplete conversion, prepared above the gel point far from the Ginzburg region defined by Equation (18). In general, *μ* − *ν* makes a sense of the cyclic rank of the network of unit volume. To define it, mentally cut one of the strands so that the network does not break into two disconnected parts. The cyclic rank is equal to the maximum number of such cuts, and it makes sense the number of independent network loops. An important limitation of expression (24) is the assumption of a model of perfect networks about the infinite sizes of all of the network cycles.

In real networks, not all loops transmit stress in the network. For example, primary loops, connected to the network at only one crosslink, are not capable of permanently supporting a stress. To exclude the contribution of such elastically ineffective loops and dangling chain ends, it was proposed to account in Equation (24) only for the elastically effective strands and crosslinks. Elastically effective strands deform and store elastic energy upon network deformation [17]. Elastically effective crosslinks connect at least two elastically effective strands [18]. The elastic modulus of real networks can significantly differ from such a modified expression (24), since each of finite size loops is

characterized by its "elastic effectiveness", which depends on the type of the loop. The replica method provides a universal tool for accounting for such effects, see Equation (37) below.

#### *2.5. Combined Chain Model: Polydisperse Networks*

Real networks are polydisperse, and therefore, the numbers of monomers *n* of virtual chains are random variables. Information on the distribution of the lengths of virtual chains is important for studying the distribution of tension in real chains of the network. Below I calculate the distribution function *p* (*n*) of virtual chain lengths. The distribution of inverse variables *s* = 1/*n* is described by the function

$$q\left(s\right) = \frac{1}{s^2} p\left(\frac{1}{s}\right), \quad s = \frac{1}{n} \tag{25}$$

Appendix B provides an algorithm for calculating this distribution for an arbitrary distribution *P* (*N*) of the number of monomers *N* of the real network strands. The problem is reduced to solving a system of nonlinear integral equations. In the most interesting case of the exponential distribution *P* (*N*) = *e*−*N*/*N*¯ /*N*¯ , the asymptotic behavior *<sup>n</sup> <sup>N</sup>*¯ of the solution of these equations can be found

$$p\left(n\right) \sim e^{-\left(f-1\right)^{2}n/N} \tag{26}$$

Of course, knowledge of only asymptotic is not enough to describe conformations of a chain in the polydisperse network. The replica method offers a much more constructive approach, allowing us to calculate this function over the entire range of *n* values. One can show, generalyzing calculations of Ref. [9] that the distribution function *q* (*s*) in Equation (25) is determined by inverse Laplace transform of the function *χf*−<sup>1</sup> (*t*) with the dimensionless order parameter *χ* (*t*) = *ϕ* (*ς*) /*ϕ* (0),

$$\int q\left(s\right)e^{-s\hat{N}t}ds = \chi^{f-1}\left(t\right)\tag{27}$$

Due to the boundary condition *χ* (0) = 1 for the function *χ* (*t*), the distribution function *q* (*s*) is normalized to unity. To better understand the meaning of Equation (27), we substitute the exponential function *χ* (*t*) for the monodisperse networks, Equation (7), into this equation. The result is a *δ*-functional distribution,

$$q\left(\mathbf{s}\right) = \delta\left[\mathbf{s} - \left(f - \mathbf{2}\right) / N\right],\tag{28}$$

in accordance with expression (20) for the number of monomers *n* = 1/*s* of virtual chains.

In the case of polydisperse networks, this correspondence, Equation (27), establishes a connection between the standard method of distribution functions and the replica approach. Since *χf*−<sup>1</sup> (*t*) is the Laplace transform of the convolution function *χ*˜(*f*−1) (*t*) (see Appendix B), the function *χ* (*t*) is the Laplace transform of the distribution function *χ*˜ (*s* ) defined by Equation (A7). Equation (8) for this function corresponds to Equation (A8) of the replica approach for the exponential distribution of the numbers of monomers of real strands.

Substituting the asymptotic expression (9) for the function *χ* (*t*), into Equation (27), we find a more accurate asymptotic expression for the distribution function

$$p\left(n\right) \simeq \frac{1}{\bar{N}} \left(\frac{n}{\bar{N}}\right)^{f/2-1} e^{-\left(f-1\right)^2 n/\bar{N}}, \quad n \gg \frac{\bar{N}}{\left(f-2\right)^2} \tag{29}$$

Using the function *χ* (*t*) from Equation (10), we get

$$p\left(n\right) \simeq \frac{c}{\bar{N}\Gamma\left[\left(f-1\right)k\right]} \left(\frac{\bar{N}}{cn}\right)^{k\left(f-1\right)+1} e^{-\frac{\bar{N}}{cn}}\tag{30}$$

Thus, the distribution *p* (*n*) of virtual chains vanishes with all its derivatives as *n* → 0 and decreases exponentially on the scale *<sup>N</sup>*¯ / (*<sup>f</sup>* <sup>−</sup> <sup>2</sup>) <sup>2</sup> for large *n*. Calculating the average number of monomers of the virtual chain with the distribution (30), we find

$$
\pi \cong \frac{f+1}{2f-1} \frac{N}{f-2} \tag{31}
$$

Note that the obtained distribution of the lengths of the virtual chains is narrower than the initial distribution of the lengths of the network strands.

#### *2.6. Generalized Combined Chain Model*

A perfect network assumes a tree-like structure on all spatial scales. In real networks, due to the excluded volume effect, trees cannot "grow" to infinity (the Malthus effect), since too large trees cannot fit in real 3D space. Therefore, in real networks, the size of a typical tree is finite, and the network consists of a large number of loops of finite length, see Figure 2c. They are strongly overlapped and interconnected with each other. The impact of such typical loops, which are responsible for the solid-state elasticity of the network, can be described by the generalized combined chain model.

In this model, in addition to two virtual chains of *n*<sup>1</sup> monomers, the combined chain includes an additional virtual chain of *n*<sup>2</sup> monomers, see Figure 2d. This virtual chain represents an effective elastic of all loops of finite size in the network, shunting the real strand of *N* monomers. The average end-to-end vector *R* of the network strand is related to the end-to-end vector *X* of the combined chain via the force balance condition, generalizing Equation (21):

$$\mathbf{R} = \frac{\mathbf{X}}{1 + 2n\_1 \left(1/N + 1/n\_2\right)}.\tag{32}$$

In this model, the mean square fluctuation of the vector **R** is

$$\left\langle \left( \Delta \mathbf{R} \right)^{2} \right\rangle = \frac{b^{2}}{1/N + 1/\left(2n\_{1}\right) + 1/n\_{2}} = \frac{b^{2}}{1/N + 1/\left(2n\right)}\tag{33}$$

The second equality can be treated as the contribution of two effective virtual chains connecting the network strand to the non-fluctuating background. They have renormalized number of monomers

$$2n = \frac{2n\_1 n\_2}{2n\_1 + n\_2} \,\text{\,\,}\tag{34}$$

corresponding to parallel connection of two effective chains: with 2*n*<sup>1</sup> and *n*<sup>2</sup> monomers.

The elastic modulus of this model is

$$G = \frac{kT\rho^{(0)}/N}{1 + 2n\left(1/N + 1/n\_2\right)} \frac{1}{1 + N/n\_2} \tag{35}$$

The factor (1 + *N*/*n*2) <sup>−</sup><sup>1</sup> is the fraction of energy related to the network strand when it is connected in parallel with the effective chain of *n*<sup>2</sup> monomers, see Figure 2d. This result can be rewritten in the form

$$G = \frac{kT\rho^{(0)}}{N + 2n} - \frac{kT\rho^{(0)}}{N + n\_2}.\tag{36}$$

where 2*n* is the number of monomers of the effective virtual chains, Equation (34). The first term in this expression has the same meaning as for the perfect network in Equation (23). The negative second contribution in Equation (36) is due to the finite size of typical loops in the network. The number of monomers of the shunt virtual chain, *n*2, can be calculated using the replica method.

#### *2.7. Finite-Size Loops of Real Networks*

In real networks, there are always cyclic defects of finite dimensions. Sparse structural defects can be taken into account within the framework of the perfect network model, assuming that they are connected with the affinely deformed non-fluctuating background through the root trees. These trees can be modeled by the virtual chains [3]. Since part of the network strands is "spent" on creating the cyclic fragment itself, the length of the effective chain linking this fragment to the elastic non-fluctuating background is different from Equation (20). In polydisperse networks, the distribution function of such virtual chains is determined by Equation (27), in which *f* − 1 has the meaning of the number of branches on the first generation of this tree. In the ideal defect gas approximation, the contribution of the structural defects to the elastic modulus of a perfect network was calculated in works [19–22]. Note that typical network loops are ignored in this approximation, which takes into account only explicitly treated loops of small concentration.

In real networks, there are both loops, binding different network layers (see Equation (14)) and topological defects, which are not sparsely distributed. The replica method allows calculating the effect of both typical loops and cyclic defects of arbitrary concentration on the elasticity of the network. The elastic modulus of the network obtained by random end-linking polydisperse chains by cross-links with functionality *f* and concentration *ρ* (0) *<sup>f</sup>* and cyclic fragments with functionality { *fi*} and arbitrary concentrations *xiρ* (0) *f* is [4]

$$\mathbf{G} = kT \rho\_f^{(0)} \frac{f/2 - 1 + \sum\_{i} \left( f\_i/2 - 1 \right) \mathbf{x}\_i}{1 + \sum\_{i} i \mathbf{x}\_i} - kT v^{(0)} \rho^{(0)} \frac{Q(0)}{2} \tag{37}$$

The functionality *fi* of the cyclic defect defined as the number of "external" network strands joined to the cyclic fragment. For the ring fragment with *i* "internal" strands and cross-links, *fi* = *i*(*f* − 2), since two of the functional groups of each cross-link are involved in creating the ring. In the case of incomplete conversion, *p* < 1 (but still far from the gelation threshold, |1 − *p*| 1), in expression (37), instead of functionalities *f* and *fi*, their average values should be substituted, *p f* and *p fi*. Such a dependence of the network elastic modulus on the conversion is in agreement with numerical simulations of Gaussian networks, performed in Ref. [16].

The fractions *xi* of cyclic defects (small fragments of the network with *i* - 1 "internal" cross-links) depend on the type of reactions leading to the formation of the network [23,24]. All *xi* universally depend on one dimensionless parameter *x*<sup>1</sup> characterizing the conditions for network preparation [25]. The primary loops with concentration *x*<sup>1</sup> are tied to the network at only one crosslink. The denominator of the first term in expression (37) describes an increase in the effective number of monomers between cross-links due to primary loops, *<sup>N</sup>*¯ <sup>→</sup> *<sup>N</sup>*¯ (<sup>1</sup> <sup>+</sup> *<sup>x</sup>*1). Therefore, although primary loops cannot bear shear stress in the network, they renormalize the length of the elastically effective chains, which deform and store elastic energy upon network deformation [7,9].

For general *f* > 2, the first term in Equation (37) can be interpreted as the contribution of elastically effective network strands with renormalized number of monomers due to the presence of topological defects in the network—primary loops and cyclic fragments of finite size. Only this contribution is predicted in the classical model of phantom networks. Its expansion in a series in the parameters *xi* with accuracy up to first order terms reproduces the result of the so-called "network theory with strand pre-strain" [21]. At network preparation conditions, the strands of a cyclic fragment are contracted and the surrounding strands are stretched (pre-strained) with respect to Gaussian sizes of these chains [22]. These effects are "automatically" taken into account in the framework of the replica approach in Equation (37).

The last term in Equation (37) is always negative and gives an impact of a large number of typical loops of the network. The factor *Q* (0) in the second term of expression (37) determines the probability of the formation of typical closed loops. These loops shunt elastically effective chains of the network, an effect that is taken into account by an additional effective chain in the model of the generalized

combined chain in Figure 2d. Unlike perfect networks, the monomers of real networks mutually repel each other, "crowding out" the network loops on small length scales, see Figure 2. Therefore, the elasticity of a polymer network, Equation (37), explicitly depends on excluded volume parameter *v*(0), characterizing monomer interaction at the preparation conditions. This effect is ignored by the classical theory of polymer networks.

Note that the elastic modulus of the polymer network *G* in Equation (37) vanishes at a finite fraction *<sup>x</sup>*<sup>1</sup> of cyclic fragments. Such a network with overlap parameter *<sup>P</sup>*(0) 1 has tree-like connections of overlapped cyclic fragments with a finite monomer density. To better understand this result, we compare it with the more familiar case of a disordered low-molecular-weight solid. In this case, the network fragments cannot overlap since the overlap parameter *P*(0) = 1. Therefore, the elastic modulus turns to zero only at the point of percolation transition, at which the network density also vanishes.

A similar effect was observed in numerical simulations of networks obtained near the gel point, which is an analog of the percolation transition in low-molecular-weight solids. It is shown that the elastic modulus of such networks vanishes at a finite density of the network monomers [26]. One can estimate the corresponding shift in critical conversion *pc<sup>μ</sup>* at which *G* = 0, compared with the gel point at the conversion *pcγ*, at which the average degree of polymerization of soluble molecules diverges. For end-linked networks, the looping probability is *<sup>Q</sup>* (0) (*bN*1/2)−<sup>3</sup> <sup>∼</sup> 1/(*NP*(0)). Dropping all the numerical factors, we find from Equation (37) the mean-field estimate for the elastic modulus

$$G \simeq kT \frac{\rho\_{\mathcal{S}}^{(0)}}{L} - kT \frac{\rho\_{\mathcal{S}}^{(0)}}{NP^{(0)}} \tag{38}$$

Here *L N*/ (*p* − *pcγ*) <sup>2</sup> and *ρ* (0) *<sup>g</sup>* (*<sup>p</sup>* <sup>−</sup> *pcγ*) *<sup>ρ</sup>*(0) is the density of gel monomers, see Equations (16) and (15). Therefore, *G* = 0 at *pc<sup>μ</sup>* − *pc<sup>γ</sup>* ∼ 1/ √ *P*(0), in accordance with the result of numerical simulations [26]. A more rigorous description of this effect requires the inclusion of fluctuation corrections. Inside the strongly fluctuating Ginzburg region, Equation (18), the overlap parameter of elastically effective chains is small (which is why strong density fluctuations develop in this region), and they can be described by critical exponents of percolation theory, see Reference [26].

#### **3. Entangled Networks**

The physics of topological entanglements is perhaps one of the most "entangled" sections of physics of polymer networks since it is impossible to give a sufficiently rigorous and constructive formalism that adequately describes entanglements for networks of macroscopic size. Therefore, the consideration of topological entanglements is usually carried out using some uncontrolled assumptions, allowing to switch from substantially "multi-chain" problem to considering the conformations of one network strand. The transfer of elastic shear stresses from individual strands of a network to its solid-state degrees of freedom (represented by the affinely deformed non-fluctuating background) is described by virtual chains. In phantom networks, such a transfer is performed only through cross-links at the strand ends. In entangled networks, the stress is carried by all entangled segments of the strand due to its topological interaction with the loops formed by adjacent network strands.

Entangled segments are characterized by the large overlap parameter

$$P\_{\varepsilon} \simeq \rho^{(0)} b^3 N\_{\varepsilon}^{1/2} \tag{39}$$

where *Ne* is the number of monomers of the entangled strand. In the case of flexible polymers, the parameter *Pe* 20–30. The number of monomers in an entanglement strand *Ne* and the overlap parameter *Pe* depend on the polymer chemistry [27].

As shown in [5], in the melt of non-concatenated rings the condition *Pe* = *const* is satisfied at all scales down to the distance *<sup>a</sup>*<sup>0</sup> *bN*1/2 *<sup>e</sup>* between adjacent entanglements. Entangled loops of different sizes overlap with similar size loops at the same overlap parameter *Pe*. As a result, ring molecules are packaged into the self-similar structures called fractal loopy globules. Unlike ordinary globules, in which chains with Gaussian statistics experience random reflections from the surface of the globule, the ring chains in the fractal loopy globules are not Gaussian and form loops at all scales, starting from a size

$$a\_0 \simeq b \mathcal{N}\_{\mathcal{e}}^{1/2} \tag{40}$$

of the entangled segment. Ring sections with the number of monomers *g* < *Ne* smaller than entanglement strands are still Gaussian with size *<sup>r</sup> bg*1/2. Larger subsections of a ring are compressed with fractal dimension *df* = 3 and size

$$r \simeq a\_0 \left(\text{g/N}\_{\varepsilon}\right)^{1/3}, \quad a\_0 < r < R \tag{41}$$

up to the size *R* of a loopy globule in a melt

$$R \simeq bN^{1/3}N\_{\varepsilon}^{1/6} \tag{42}$$

Large typical values of the parameter *Ne* 1 allow developing the mean-field theory of entanglements in polymer networks, since the total effect of large number *Pe* 1 of loops interacting with the entanglement segment is effectively averaged. The strands of the entangled polymer networks are confined in an effective entangled tube with a diameter *a*, which, under the conditions of the network preparation, is equal to the size *a*<sup>0</sup> of the entangled segment. The tube rotates by a random angle of about 90◦ on the scale *a*0, see Figure 3a.

**Figure 3.** Nonaffine tube model. (**a**) In an undeformed state, the chain fluctuates in entangled tube of diameter *<sup>a</sup>*<sup>0</sup> (**b**) The tube diameter *az* in the stretching direction is less than the affine length *<sup>R</sup>aff <sup>z</sup>* . The dashed lines show the trajectory of the primitive path, which is obtained by smoothing the affine deformed coordinates of the original Gaussian chain on the tube diameter scale *a*.

In the case of networks obtained by crosslinking a melt of linear chains, the presence of a large number of overlapping entangled chains allows using the mean-field approximation. In this approximation, topological interactions are described by virtual chains attached to all monomers (*s* = 1, . . . , *N*) of the chain with coordinates **x**(*s*). The elastic energy of Gaussian virtual chains is

$$\sum\_{s} \frac{\Im \left( \mathbf{r} \left( s \right) - \mathbf{X} \left( s \right) \right)^{2}}{2b^{2} N\_{\varepsilon}^{2}} \tag{43}$$

Here *r<sup>α</sup>* (*s*) = *x<sup>α</sup>* (*s*) /*λα* and *X<sup>α</sup>* (*s*) (*α* = *x*, *y*, *z*) are undeformed coordinates of two ends of the *s*-th virtual chain. As shown in [3], with this choice of the potential of virtual chains, they do not contribute to the elastic stress of the polymer network, which is determined only by the contribution of all real polymer chains.

#### *3.1. Physics of Entangled Network Deformation*

Consider the case of a highly entangled network, *N*/*Ne* 1, with many entanglements per network chain. As shown in [28], the network deforms affinely like a solid only at spatial scales exceeding the affine length *Raff <sup>α</sup>* , see Figure 3. On shorter scales, chains take on liquid-like conformations. In the case of anisotropic deformation of the gel by the ratios *λα* along the axes *α* = *x*, *y*, *z*, the affine length is also anisotropic and depends on the direction *α*:

$$R\_{\alpha}^{aff} \simeq \lambda\_{\alpha} a\_0 \tag{44}$$

We define an affine strand of the size *Raff <sup>α</sup>* directed along the axis *<sup>α</sup>*, which consists of *<sup>N</sup>aff α* monomers and deforms by the stretching ratio *λα*. The number of chain entanglements does not change upon network deformation; therefore, the diameter of the deformed tube *a<sup>α</sup>* is equal to the size of the entangled segment of *Ne* monomers. In the stretched network this segment is a section of a longer stretched affine strand of size *Raff <sup>α</sup>* > *aα*:

$$a\_{\alpha} \simeq R\_{\alpha}^{aff} \left( \mathcal{N}\_{\varepsilon} / \mathcal{N}\_{\alpha}^{aff} \right) \tag{45}$$

Fluctuations of the affine strand *b Naff α* 1/2 are limited by entanglements and therefore equal to the tube diameter

$$b\left(N\_a^{aff}\right)^{1/2} \simeq a\_a\tag{46}$$

The solution of Equations (45) and (46) is

$$N\_{\mathfrak{A}}^{aff} = \lambda\_{\mathfrak{A}} N\_{\mathfrak{E}} \tag{47}$$

Thus, the entanglement efficiency decreases with a network stretching, and the entanglements are not effective at *Naff <sup>α</sup>* > *N*; such an entangled network behaves like a phantom one. The free energy of the entangled network is estimated as the elastic energy of all its stretched entangled segments,

$$\frac{\mathcal{F}\_{\mathcal{E}}}{V} \simeq \frac{\rho kT}{N\_{\mathcal{E}}} \sum\_{\alpha} \frac{\left(R\_{\alpha}^{aff}\right)^{2}}{a\_{\alpha}^{2}} \simeq \frac{\rho kT}{N\_{\mathcal{E}}} \sum\_{\alpha} \lambda\_{\alpha} \tag{48}$$

A more rigorous estimate for the free energy of the non-affine tube model with the number of monomers of the chain *N Ne* gives [29]

$$\frac{\mathcal{F}\_{\varepsilon}}{V} = \frac{\rho kT}{2N\_{\varepsilon}} \sum\_{\mathfrak{a}} \left( \lambda\_{\mathfrak{a}} + \frac{1}{\lambda\_{\mathfrak{a}}} \right)\_{\prime} \tag{49}$$

This expression includes an additional to Equation (48) contribution ∼1/*λα* describing the compression of the chain in the tube.

#### *3.2. Slip-Tube Model*

This model generalizes the nonaffine tube model by taking into account the slippage of chains along the tube contour. In the case of anisotropic deformation of the network, its chains are more stretched in the direction of maximum extension relative to other directions. The increased chain tension in this direction draws the chain from the tube segments along these directions. The redistribution of the chain length between different sections of the tube is taken into account by renormalization of parameters of the non-affine tube model for each direction of deformation

$$
\lambda\_{\mathfrak{u}} \to \lambda\_{\mathfrak{u}} / \mathfrak{g}\_{\mathfrak{u}}^{1/2} \, \prime \qquad \mathcal{N} \to \mathcal{N} / \mathfrak{g}\_{\mathfrak{u}} \, \tag{50}
$$

The parameters *gα* are normalized by the condition of preserving the full length of the chain,

$$\sum\_{\mathfrak{a}} \mathbf{g}\_{\mathfrak{a}} = \mathbf{3} \tag{51}$$

and in the case of uniaxial deformation, Equation (3), there is only one independent redistribution parameter *gz*:

$$\mathbf{g}\_{\mathbf{x}} = (\mathbf{3} - \mathbf{g}\_{\mathbf{z}}) / \mathbf{2} \tag{52}$$

In addition to the above renormalization of Equation (49), the free energy of this model has an entropic contribution taking into account the chain slippage along the tube,

$$\frac{\mathcal{F}\_{\varepsilon}}{V} = \frac{\rho kT}{2N\_{\varepsilon}} \sum\_{a} \left( \frac{\lambda\_{a}}{g\_{a}^{1/2}} + \frac{g\_{a}^{1/2}}{\lambda\_{a}} \right) - \frac{\rho kT}{3N\_{\varepsilon}} \sum\_{a} \ln g\_{a\prime} \tag{53}$$

The parameters *g<sup>α</sup>* are found by minimizing the free energy (62), and the minimization equation reads as

$$h\_z\left(\mathcal{g}\_z\right) = h\_x\left(\mathcal{g}\_x\right), \quad h\_\mathbb{il}\left(\mathcal{g}\_\mathbb{il}\right) = \frac{\partial}{\partial \mathcal{g}\_\mathbb{il}} \frac{\mathcal{F}\_\mathcal{e}}{\nu kT} \tag{54}$$

The elastic stress of the slip-tube model is

$$\frac{\sigma\_{\alpha\alpha}^{\varepsilon}}{\nu kT} = \frac{\lambda\_{\alpha}}{\nu kT} \frac{d\mathcal{F}\_{\varepsilon}}{d\lambda\_{\alpha}}.\tag{55}$$

The effect of entanglements substantially depends on the degree of network stretching, and the chains of strongly stretched networks become effectively phantom. Such a crossover from entangled to affine behavior is usually described using the assumption of additivity of the phantom and affine contributions to the free energy:

$$\mathcal{F} = \mathcal{F}\_{\text{ph}} + \mathcal{F}\_{\text{\textquotedblleft}\text{\textquotedblright}} \tag{56}$$

This assumption also leads to the additivity of the corresponding contributions to the elastic tensor,

$$
\sigma\_{\mathfrak{uu}} = \sigma\_{\mathfrak{uu}}^{\text{ph}} + \sigma\_{\mathfrak{uu}}^{\varepsilon} \tag{57}
$$

In the case of uniaxial network deformation, Equation (3), a convenient representation of stress is the Mooney ratio of the total stress *<sup>σ</sup>zz* − *<sup>σ</sup>xx* to the functional dependence *<sup>λ</sup>*<sup>2</sup> − *<sup>λ</sup>*−<sup>1</sup> predicted by phantom network models:

$$f^\*\left(\lambda^{-1}\right) = \frac{\sigma\_{zz} - \sigma\_{xx}}{\lambda^2 - \lambda^{-1}}\tag{58}$$

As was shown in Ref. [3], the Mooney plot of the experimental data is well described by the slip tube model.

#### *3.3. Crossover from Unentangled to Entangled Regimes*

The Mooney-Rivlin dependence (58) is usually fitted by the phenomenological equation

$$f^\*\left(\lambda^{-1}\right) \simeq 2\mathbb{C}\_1 + 2\mathbb{C}\_2/\lambda\tag{59}$$

It is usually assumed that the parameter 2*C*<sup>1</sup> is only due to chemical cross-links, while the parameter 2*C*<sup>2</sup> includes all of the entanglement contributions. According to the non-affine tube model, part of the entanglement contribution is also included in the parameter 2*C*1[6].

In this section, we study the concentration dependence of the Mooney-Rivlin coefficients. Although the free energy additivity approximation (56) can be used to describe the polymer network in the melt, a refined theory should be developed to find the Mooney dependence in the entire concentration range. The importance of such a study is emphasized by the fact that the networks prepared near the crossover from unentangled to entangled regimes are a fairly typical case.

Consider a network, consisting of *ν* Gaussian chains with *N* monomers whose endpoints are displaced affinely with the global deformation of the network. We assume that each chain of this network is confined in an effective "non-affinely deformed" tube, described by the virtual chain's potential (43). The free energy of this model is calculated as the sum of contributions of all its strands and is found in Appendix C:

$$\frac{\mathcal{F}}{\nu kT} = \sum\_{\alpha} \left( \frac{\lambda\_{\alpha}^{2} - 1}{2} \frac{k\_{\alpha}}{\tanh k\_{\alpha}} + \ln \frac{\sinh k\_{\alpha}}{k\_{\alpha}} \right),\tag{60}$$

where *<sup>k</sup><sup>α</sup>* <sup>=</sup> *<sup>N</sup>*/ (*Neλα*). In the "phantom" limit *<sup>N</sup> <sup>N</sup>aff <sup>α</sup>* = *λαNe* Equation (60) turns into the free energy of affine network model, Equation (4). In the limit of strongly entangled networks, *<sup>N</sup> <sup>N</sup>aff <sup>α</sup>* , this expression reproduces the free energy of the non-affine tube model, Equation (49). The first term in Equation (60) can be rewritten as

$$\left(\frac{\lambda\_a R\_0}{R\_{fl,a}}\right)^2\tag{61}$$

where *R*<sup>0</sup> = *bN*1/2 is the chain size under network preparation conditions and *Rf l*,*<sup>α</sup>* is the fluctuation size in the *α*-direction, equal to *R*<sup>0</sup> in the phantom regime and the tube diameter *a<sup>α</sup>* in the entangled regime. The logarithmic term in this expression describes the entropy of the chain compression in the tube.

By renormalizing the parameters of this model (50) to take into account the effect of chain redistribution in the entangled tube, and adding the contribution of entropy of slippage along the tube (53), we arrive at the expression for free energy

$$\frac{\mathcal{F}}{\nu kT} = \sum\_{\mathfrak{a}} \left( \frac{\lambda\_{\mathfrak{a}}^2 - g\_{\mathfrak{a}}}{2} \frac{k\_{\mathfrak{a}}}{\tanh k\_{\mathfrak{a}}} + g\_{\mathfrak{a}} \ln \frac{\sinh k\_{\mathfrak{a}}}{k\_{\mathfrak{a}}} - \frac{N}{3N\_{\mathfrak{e}}} \ln g\_{\mathfrak{a}} \right), \tag{62}$$

Solving equations (54) for this model, we find the stress (55) in the polymer network, which determines the Mooney-Rivlin dependence (58). The parameter *C*<sup>2</sup> (*φ*) of this dependence is determined by the expression

$$\text{2C}\_2\left(\phi\right) = \frac{df^\*\left(\lambda^{-1}\right)}{d\lambda^{-1}}.\tag{63}$$

We plot the dependences *C*<sup>1</sup> (*φ*) and *C*<sup>2</sup> (*φ*) (63) for *λ*−<sup>1</sup> = 0.7 and *N*/*Ne* = 2.4 in Figure 4. The parameter *C*<sup>1</sup> of the Mooney-Rivlin dependence (59) weakly depends on the polymer volume fraction *φ*. In accordance with the experimental data [30], the dependence of the parameter *C*<sup>2</sup> on *φ* is almost linear, and *C*<sup>2</sup> vanishes at *φ* 0.2. The vanishing of *C*<sup>2</sup> at finite *φ* is related to the transition between phantom and entangled regimes. When the polymer swells, its chains move away from each other, which leads to a weakening of the effective topological potential that holds the chains in the tube. Chain fluctuations in this potential increase with the swelling and can reach the fluctuations of the phantom chain. This concentration corresponds to the disappearance of the contribution of entanglements to the network elasticity and the vanishing of the Mooney-Rivlin coefficient, *C*<sup>2</sup> = 0.

**Figure 4.** Dependence of the Mooney-Rivlin coefficients *C*<sup>1</sup> (*φ*) (upper curve) and *C*<sup>2</sup> (*φ*) (lower curve) on the volume fraction of polymer *φ*: theory (solid curves) and experimental data (points, Kg/cm2) for swollen rubber [30].

#### **4. Heterogeneities in Polymer Networks**

Network inhomogeneity is a common feature of polymer networks and gels due to the randomness of the crosslinking process and the presence of additional topological defects such as dangling chain ends, cross-linker shortcuts, and chains forming loops. The origin of nanostructural inhomogeneities and their characterization by light, neutron, and X-ray scattering as well as by NMR spectroscopy and optical, electron, and X-ray microscopies is reviewed in Ref. [31], and the main methods of their study are outlined in the review [32]. The first attempts to take such heterogeneities into account were made by using the model of randomly cross-linked networks containing fractal regions, such as percolation clusters [33]. Networks with fractal heterogeneities can be swollen and deformed by unfolding the fractal regions without significant elastic entropy penalty [34]. It is shown that strong heterogeneities lower the shear modulus of the network if a part of strands is so short to be considered rigid and not able to deform. Networks with a large number of such very short strands have higher breaking energies.

The effect of entanglements on heterogeneities in polymer networks is studied in Ref. [35] using molecular dynamics simulations of polymer networks made by either end-linking or randomly crosslinking a melt of linear precursor chains. The end-linking leads to nearly ideal monodisperse networks, while random cross-linking produces strongly polydisperse networks. It is shown that the microscopic strain response, the diameter of the entanglement tube, and stress–strain relation weakly depend on the linking process by which the networks were made.

The replica method was used in Ref. [28] to calculate the characteristic size and amplitude of the spatial nonuniformities of the network due to defects of its structure and topological restrictions. Using this method it is shown that inhomogeneities can arise as consequences of a stretching of polymer networks [36]. Although the replica theory [8] provides a complete solution of the statistical mechanics of polymer gels, it uses replica trick which is unfamiliar to the majority of people in the polymer community. A more intuitive phenomenological approach capturing all the main physical ingredients of the complete theory is developed in Ref. [37].

#### *4.1. Theory of Heterogeneities in Polymer Networks*

In this section, we review the main results of the theory of random spatial heterogeneities developed in swollen and deformed polymer networks [38]. The non-triviality of this problem stems from the fact that information about network structure is "encrypted" in the pattern of cross-links joining polymer chains, which represent a very small fraction of the network volume. The initial crosslinking pattern is reproduced only partially due to thermal fluctuations arising in the new equilibrium state after crosslinking. Conformations of polymer strands in such networks with fixed topological structure can be varied in a wide range depending on experimantal conditions.

The density profile of monomers in the polymer network can be recovered from the Fourier component of the deviation of the density from its average value *ρ*:

$$
\rho(\mathbf{x}) = \rho + \int \tilde{\rho}\left(\mathbf{q}\right) e^{i\mathbf{q}\mathbf{x}} \frac{d\mathbf{q}}{(2\pi)^3} = \rho\_{\ell\eta}(\mathbf{x}) + \delta\rho(\mathbf{x})\tag{64}
$$

Here, *δρ*(**x**) are random deviations (due to thermal fluctuations) of the density from the equilibrium density profile *ρeq*(**x**) describing spatial inhomogeneities in polymer networks. The fluctuation contribution to the free energy of the network with a given distribution of monomer density can be represented as the sum of the entropy contribution and the (osmotic) contribution of volume interactions ∼ *v*:

$$\mathcal{F}\left[\vec{\rho}\right] = \frac{kT}{2} \int \left[\frac{|\vec{\rho}\left(\mathbf{q}\right) - \tilde{n}\left(\mathbf{q}\right)|^2}{\tilde{\mathbf{g}}\left(\mathbf{q}\right)} + v\left|\vec{\rho}\left(\mathbf{q}\right)|^2\right] \frac{d\mathbf{q}}{(2\pi)^3} \tag{65}$$

Here *<sup>n</sup>* (**q**) is the density profile in the "elastic reference state" [37], maximizing the entropy of the polymer network. The density *<sup>n</sup>* (**q**) is determined by the molecular structure of the network and it vanishes in the short-wavelength limit *<sup>q</sup> <sup>R</sup>*<sup>−</sup>1, since solid-state degrees of freedom are determined only at length scales exceeding the monomer fluctuation radius *R*. On a smaller scale, the gel has liquid degrees of freedom, contributing to the temperature structural factor

$$\mathfrak{F}(\mathbf{q}) = \rho N \left[ \frac{1}{Q^2/2 + (4Q^2)^{-1} + 1} + \frac{2Q^2}{(1 + Q^2)^2 (\boldsymbol{\lambda} \cdot \mathbf{Q})^2} \right] \tag{66}$$

The first term in square brackets determines the contribution of the liquid degrees of freedom of the polymer network. The dimensionless wavevector **Q** = *R***q** is normalized by the monomer fluctuating radius, and the vector *<sup>λ</sup>* · **<sup>Q</sup>** has components *λxQx*, *λyQy*, *λzQz* . At large *Q* 1 the term *Q*2/2 in the denominator of the first term of the right hand side of Equation (66) gives the usual Lifshitz entropy of polymer solutions, *g*˜ (**q**) = 2*ρN*/*Q*<sup>2</sup> [39]. The term (4*Q*2)−1, first introduced by de Gennes for heteropolymer networks [40], describes the suppression of density fluctuations on length scales larger than the monomer fluctuation radius *R*. The second term in square brackets in Equation (66) determines the contribution of solid-state degrees of freedom of the polymer network. It remains finite in the long-wavelength limit and retains its angular dependence on the anisotropic deformation even in the limit *q* → 0.

Equilibrium monomer density profile is found by minimizing the free energy (65)

$$\tilde{\rho}\_{\epsilon\eta} \left( \mathbf{q} \right) = \frac{\tilde{n} \left( \mathbf{q} \right)}{1 + v \tilde{\mathbf{y}} \left( \mathbf{q} \right)} \tag{67}$$

As can be seen from Equation (67), in good solvent with positive excluded volume *v* > 0, the equilibrium density profile is more homogeneous than that of the corresponding elastic reference state. The excluded volume interaction suppresses not only thermal fluctuations in the polymer network, but also frozen-in spatial inhomogeneities. Although the static spatial fluctuations are permanent inhomogeneities, they can reversibly increase and diverge at the spinodal line, at which *vg*˜ (0) = −1. This observation is experimentally confirmed in Ref. [41].

In (ergodic) heterogeneous systems, there are two types of averages, i.e., the *thermal* or *time averages* and *ensemble* or *space averages*, denoted by *X* and *X*, respectively. The Fourier component of the thermal correlator of density fluctuations *δρ* (**x**) = *ρ* (**x**) − *ρeq* (**x**) is found by averaging with the Gibbs weight *e*−F[*ρ*˜]/*kT*, Equation (65):

$$\tilde{D}\left(\mathbf{q}\right) = \left\langle \left| \delta\tilde{\rho}\left(\mathbf{q}\right) \right|^2 \right\rangle = \frac{\tilde{\mathbf{g}}\left(\mathbf{q}\right)}{1 + v\tilde{\mathbf{g}}\left(\mathbf{q}\right)}\tag{68}$$

Therefore, in addition to permanent spatial heterogeneities, a polymer network undergoes thermal dynamical density fluctuations which diverge at the spinodal line.

All information about the heterogeneity of the molecular structure of the network is contained in the monomer density *<sup>n</sup>* (**q**) in the "elastic reference state", which can be considered as a random variable characterized by a correlator

$$\nu\left(\mathbf{q}\right) = \overline{\tilde{n}\left(\mathbf{q}\right)\tilde{n}\left(-\mathbf{q}\right)} = \frac{1}{(1+Q^2)^2} \left[6\rho N + \frac{9}{\lambda\_x \lambda\_y \lambda\_z} \tilde{S}^{(0)}\left(\lambda \cdot \mathbf{q}\right)\right] \tag{69}$$

Therefore, it is important to understand the physical meaning of the different terms of this expression:

Each of two factors (1 + *Q*2)−<sup>1</sup> describes the thermal "smearing" of the density response to random stresses in the polymer network. The first term in square brackets comes from the correlator of frozen random stresses, which exist even in a homogeneous polymer network. Such stresses appear due to heterogeneities in the distribution of cross-links. The second term in brackets in Equation (69), which is proportional to the correlator of the affine deformed density pattern in the reference state *<sup>S</sup>*(0) (*<sup>λ</sup>* · **<sup>q</sup>**), strongly depends on network preparation conditions.

(1) In the case of cross-linking in a melt, density fluctuations are suppressed and *S*˜(0) (**q**) = 0. The fact that *ν* (**q**) does not vanish upon swelling of such a network, despite the absence of density fluctuations at the preparation conditions, means that there are still inhomogeneities in the crosslink density, which manifest themselves upon swelling.

(2) In the case of instant cross-linking of a semi-dilute polymer solution, *S*˜(0) (**q**) is given by the Ornstein-Zernicke expression

$$\left|\tilde{S}^{(0)}\left(\mathbf{q}\right) = \left\langle \left| \delta \tilde{\rho}^{(0)}\left(\mathbf{q}\right) \right|^{2} \right\rangle \right\rangle = \frac{\rho^{(0)}N}{v^{(0)}\rho^{(0)}N + Q^{2}/2},\tag{70}$$

which is finite for any second virial coefficient *v*<sup>0</sup> at network preparation conditions.

(3) In case of equilibtium chemical cross-linking, the total structure factor of the gel in the state of preparation is determined by the expression for the polymer liquid, Equation (70), with renormalized excluded volume parameter, *<sup>v</sup>*(0) <sup>→</sup> *<sup>v</sup>* (0) *ren* <sup>=</sup> *<sup>v</sup>*(0) <sup>−</sup> (*ρ*(0)*N*)−1. The decrease in *<sup>v</sup>*(0) is due to the fact that the monomers forming cross-links give an additional attractive contribution to the effective second virial coefficient *v* (0) *ren*. The amplitude of heterogeneities increases when approaching the *cross-link saturation threshold* at *<sup>v</sup>*(0)*ρ*(0)*<sup>N</sup>* <sup>=</sup> 1 at which both the structure factor *<sup>S</sup>*(0) (**<sup>q</sup>** <sup>→</sup> <sup>0</sup>) and the characteristic size of heterogeneities at preparation conditions

$$\xi \cong R/\sqrt{v^{(0)}\rho^{(0)}N - 1}$$

diverge.

The theory also predicts the appearance of a maximum degree of spatial gel inhomogeneity at a critical polymer network concentration, as confirmed by experiments on PAAm gels [42,43].

#### *4.2. Scattering Intensity*

The scattering intensity on wavevector **q** is proportional to the structure factor, which is given by the sum

$$\tilde{S}\left(\mathbf{q}\right) = \overline{\left\langle \left| \tilde{\rho}\left(\mathbf{q}\right) \right|^{2} \right\rangle} = \tilde{D}\left(\mathbf{q}\right) + \tilde{C}\left(\mathbf{q}\right) \tag{71}$$

of contributions of the thermal fluctuations, Equation (68), and the inhomogeneous equilibrium density variations due to defects of the topological structure of the network,

$$\tilde{\mathcal{C}}\left(\mathbf{q}\right) = \overline{\left|\tilde{\rho}\_{\epsilon\eta}\left(\mathbf{q}\right)\right|^{2}} = \frac{\nu\left(\mathbf{q}\right)}{\left[1 + \nu\tilde{\mathbf{y}}\left(\mathbf{q}\right)\right]^{2}}\tag{72}$$

where bar means ensemble average.

The experimental data can be visualized using contour plots of neutron scattering from random inhomogeneities of network structure in anisotropic-deformed swollen gels [44]. Under uniaxial gel stretching, an increase in the scattered intensity on small wave vectors *q* in the stretching direction is observed, which is enveloped by elliptical patterns at larger *q* values with a maximum oriented normal to this axis. This behavior is opposite of what is expected in theories assuming only thermal fluctuations and is called "abnormal butterfly patterns". The elliptical patterns at large wavevector *q* originate from the correlator of static inhomogeneities, Equations (72) and (69), which contains the term *S*<sup>0</sup> (*λ* · **q**) that "remembers" the affinely deformed inhomogeneous structure of the network. The butterfly patterns along the stretching axis in the small *q* range is due to strong angular dependence of the thermal structure factor, (68). This function comes into numerator of the thermal correlator in Equation (68), describing "normal butterfly patterns". At high cross-link concentrations the correlator of static inhomogeneities gives the main contribution to the scattering intensity. Since the function *g*˜ (**q**) is included in the denominator of Equation (72), this expression describes the so named "abnormal butterfly patterns". As a result, under uniaxial extension, a crossover from normal to abnormal butterfly scattering patterns occurs with an increase of the strength of inhomogeneity or the swelling ratio [45]. In addition to such butterfly scattering patterns, the theory also predicted "Lozenge" patterns if only part of all network chains, i.e., their deuterated fragments, could scatter neutrons, in accordance with scattering experiments [46].

The static inhomogeneity in poly (N-isopropyl acrylamide) gel (PNIPA) has been investigated in Ref. [47] by the methods of small-angle neutron scattering (SANS) and neutron spin echo. The obtained SANS scattering amplitude *S*˜ (*q*) was successfully decomposed into the thermal and static components, respectively, *D*˜ (*q*) and *C*˜ (*q*) in Equation (71). It was revealed that *C*˜ (*q*) becomes dominant in the *q*-region, where the abnormal butterfly scattering under stretching is observed. As the temperature increases toward the temperature for volume phase transition, *C*˜ (*q*) approximated by the square of the Lorentzian shape increases more drastically than *D*˜ (*q*) of the Lorentzian shape. These experimental findings are also well described in the theoretical framework of this section.

In Ref. [48] gels prepared by two cross-linking methods were studied using SANS technique. One is chemical cross-linking with BIS and the other is gamma-ray cross-linking of a PNIPA solution. It is shown that the degree of the inhomogeneity is much larger in chemically cross-linked gels than in the gamma-ray cross-linked gels. The experimental data are in quantitative agreement with predictions of the theory of scattering intensity *S*˜ (**q**) on chemically and instantaneously cross-linked gels, respectively.

#### *4.3. Heterogeneities in Charged Gels*

A study of the structure factor of weakly charged polyelectrolyte gels under uniaxial stretching was performed by Mendes et al. [49], who observed after introducing ions to the gel, the disappearance of the "butterfly pattern" in the small angle scattering intensity, as well as an increase in the scattering intensity in the direction perpendicular to the gel stretching. The origin of this maximum has been elucidated in SANS experiment by Shibayama et al. who studied the deformed state of weakly charged polymer gels PNIPA/AAc [50]. In this experiment, an anisotropic scattering maximum is observed, which indicates that the spatial distribution of the charged groups changes as a result of gel deformation and therefore is strongly coupled with the static inhomogeneities.

All the observed patterns of SANS intensity were well reproduced using a generalization of the above theory to the case of charged polymer networks [51]. The only effect of electrostatic interactions is to replace the excluded volume parameters *v*(0) and *v* with effective virial coefficients for both the final state and for the state of preparation:

$$\begin{aligned} v &\rightarrow v\left(\mathbf{q}\right) \equiv v + 1/s\_{DH}\left(\mathbf{q}\right), \\ v^{(0)} &\rightarrow v^{(0)}\left(\mathbf{q}\right) \equiv v^{(0)} + 1/s\_{DH}^{(0)}\left(\mathbf{q}\right) \end{aligned}$$

where

$$\frac{1}{s\_{DH}\left(\mathbf{q}\right)} = \frac{4\pi l\_B a^2}{\kappa^2 + q^2}$$

is the inverse of the Debye-Hückel structure factor, *κ*−<sup>1</sup> is the Debye screening length, *α* is the degree of ionization (fraction of charged monomers) and *lB* is the Bjerrum length. *s* (0) *DH* (**q**) is obtained by substituting *κ* = *κ*(0) and *α* = *α*(0) into the above equation. As shown in Ref. [52], the above theoretical prediction for *S* (**q**) well reproduces the observed scattering intensity functions of wekly charged PNIPA/AAc gels.

In general, an introduction of cross-links to a polymer solution leads to an increase in the scattering intensity due to static inhomogeneities. However, a reverse phenomenon, called the "inflection" in scattering intensity, was predicted by the above theory [53] and observed in weakly charged gels and polymer solutions [54]. While the gel becomes more inhomogeneous with increasing the degree of cross-linking in a good solvent, the inhomogeneities can be suppressed in a poor solvent, although in a relatively small regions of cross-link concentration. However, this phenomenon is interesting due to its physical significance.

#### *4.4. Good Solvent: Scaling Approach*

Earlier in this paper, we considered the networks with Gaussian chains obtained in the melt or the *θ*-solution of linear chains. In general, polymer networks can also be obtained by crosslinking semi-diluted polymer solutions. The networks can be placed in a good solvent, in which the polymer chains swell compared to the case of theta solvent. This section explains the basic ideas of the scaling approach to describing networks prepared/swollen in a good solvent.

The mean-field approach can be adapted to descripbe the elasticity of such gels [55] using the well-known de Gennes blob picture of semi-dilute solutions [56]. The key idea is the spatial scale separation: while static density inhomogeneities exist only on scales comparable to or larger than the monomer fluctuation radius *R*, thermal density fluctuations are dominated by smaller scales and are quite similar to those in semi-dilute polymer solutions.

Consider a gel prepared by random cross-linking of chains in a semi-dilute polymer solution in a good solvent at the monomer density *ρ*(0) that is swollen to density *ρ* < *ρ*(0). The monomer density fluctuations should be taken into account both in the initial (where the gel was prepared) and in the final (where it is being studied) states of the gel. On length scales shorter than the correlation lengths in these states,

$$\xi^{(0)} = b^{-5/4} \left(\rho^{(0)}\right)^{-3/4} \qquad \text{and} \qquad \xi = b^{-5/4} \rho^{-3/4} \,. \tag{73}$$

density fluctuations are large, and the gel behaves like a polymer solution ("liquid-like" regime). On scales exceeding the blob sizes (the corresponding wave vectors *q*(0) = 1/*ξ*(0) and *q* = 1/*ξ*), density fluctuations are small, and the mean-field description with appropriately renormalized parameters can be used [8]

$$\begin{aligned} b^{(0)} \rightarrow b\_{ren}^{(0)} &= b \left( \rho^{(0)} b^3 \right)^{-1/8}, \\ b \rightarrow b\_{ren} &= b \left( \rho b^3 \right)^{-1/8}, \quad \lambda\_\alpha \rightarrow \lambda\_\alpha b\_{ren}^{(0)} / b\_{ren} \end{aligned}$$

where the subscript *ren* stands for renormalized values that differ from the "bare" ones. The renormalized second virial coefficients *v* (0) *ren* and *vren* are:

$$
v^{(0)} \rightarrow v\_{\rm ren}^{(0)} = b^3 \left(\rho^{(0)} b^3\right)^{1/4}, \qquad v \rightarrow v\_{\rm ren} = b^3 \left(\rho b^3\right)^{1/4} \tag{74}$$

Equations (73) and (74) complete the renormalization of the mean-field theory: in order to describe a gel in a good solvent on length scales larger than the thermal correlation length, one only have to replace the bare parameters in the previously derived expressions for the free energy, correlation functions, etc., by their renormalized values.

#### *4.5. Amplification of Cross-Linking Density Pattern*

Is it possible to "write" some "useful" information in a polymer network cross-link pattern, and under what conditions it can be "read" back? The answer to this question is given in Ref. [57]. After the formation of a homogeneous (on length scales large compared to its "mesh" size *R*) network, large-scale patterns can be generated in it by further cross-linking followed by swelling (and possibly stretching) of the network, which leads to a nonuniformly swollen gel. This additional cross-linking can be done by adding light-sensitive cross-links to the transparent network. By focusing the laser beam in the regions inside the gel, it is possible to "write" information into the gel structure in the form of 2*D* or 3*D* patterns of cross-link density.

It was shown that although such information is hidden at the preparation conditions, it can be reversibly "developed" by gel swelling, since unobservable variations of the cross-link density in the melt are transformed into observable variations of monomer density in the swollen gel. The gel regions with increased cross-link concentration can be considered as inclusions with enhanced elastic modulus. It has been shown that in swollen gels that stretch isotropically upon absorption of the solvent, the observed monomer density pattern is an affinity stretched initial cross-linking pattern. Such gels can serve as a magnifying glass that enlarges the initially written pattern without distorting its shape. The corresponding magnification factor can be very large in the case of super-elastic networks.

A strong enhancement of the contrast between the high and the low monomer density regions can be obtained by placing the gel in a poor solvent with a negative second virial coefficient *v* < 0. The observable image becomes distorted, especially near the edges and corners of the pattern. Gel boundaries should be fixed in order to avoid its contraction. Gel contraction can also be prevented by focusing a laser beam only on a part of the localized pattern and heating it, resulting in a local change of the quality of solvent.

#### **5. Discussion**

Despite its more than half-century history, even now the theory of polymer network elasticity asks more questions than it gives answers. In this work, we tried to bridge the gap between the main developed approaches, which allowed us to give answers to some of these questions.

In real networks, there are both topological defects and typical loops, which are not sparsely distributed. Using the replica method, the cumulative effect of a large number of strongly overlapping typical loops of finite size can be described in the effective mean-field approximation. This approximation works for polymer networks due to large overlap parameter of network strands, *<sup>P</sup>*(0) 1, see Equation (12).

In the replica approach, the properties of polymer networks are described by the liquid-solid (sol-gel) order parameter *ϕ*<sup>1</sup> (*ς*), which is actually not a parameter, but a function of the variable *ς*, and is determined by a complex integro-differential equation. The only, but very important exception is the elastic modulus, which is expressed through the value of the order parameter at *ς* = 0, and the corresponding equations for *ϕ*<sup>1</sup> (0) become algebraic. This greatly simplifies the calculations of the elastic modulus of the network with an arbitrary number of cyclic fragments of finite size. In contrast to the classical theory of phantom networks, the mean field of loops, *ϕ*<sup>1</sup> (0) explicitly depends on the

excluded volume parameter *v*(0) and the density of monomers *ρ*(0) at network preparation conditions. This dependence takes into account the limitations imposed by the packing of highly overlapping typical loops in real 3*D* space on the molecular structure of the network being formed.

**T**he resulting elastic modulus of real networks has two main contributions, see Equation (37). The contribution of elastically effective network strands is taken into account in the classical model of phantom networks. The classical expression for this contribution is renormalized due to the presence of topological defects in the network—primary loops and cyclic fragments of finite size. The second contribution to the elastic modulus is always negative and gives an impact of a large number of typical loops of the network. It describes the effect of shunting of the elastically effective chains by finite size polymer loops and depends on the interaction of monomers at the preparation conditions of the network.

Both contributions to the network elastic modulus can be represented by the generalized combined chain model with an additional effective chain (see Figure 2b) that accounts the cumulative effect of finite size loops in real networks. The virtual chains transmit local stresses from the network fragments to the solid-state degrees of freedom of such soft solids. We established the connection of the distribution function of the lengths of virtual chains with the order parameter of the replica network model and calculated this distribution function for the exponential distribution of the lengths of the network strands.

We also proposed a generalization of the slip tube model of entangled polymer networks, which allows us to describe the crossover between the entangled and phantom regimes of network swelling. The dependence of the Mooney-Rivlin parameters *C*<sup>1</sup> and *C*<sup>2</sup> on the polymer concentration calculated for this model is in agreement with the experiments.

Although from the very beginning the heterogeneities were recognized as one of the most essential features of gels [58], it took long time to formulate their theoretical description due to complexity of this problem. Up to today, the understanding of the gel structure has greatly improved owing to both theoretical developments and a large number of experimental studies. Effect of cross-links on heterogeneous structure of polymer gels, abnormal butterfly patterns, microphase separation, and so on, are well understood with the aid of the theretical approaches discussed in Section 4.

The list of remaining questions is much longer:

Polymer networks can be obtained by crosslinking semidiluted polymer solutions. How does a strong nonlinear deformation of such networks occur in the phantom regime? How do such networks fracture when they are strongly stretched? The molecular theory (Lake-Thomas model) of the fracture of dry networks is constructed in Ref. [59]. What happens when compressing such networks? What is deformation mechanism of heterogeneous super-tough networks, kinetics of crack growth and phase transitions in such networks? Stress-induced microphase separation in multicomponent networks was studied in [60]. What is the effect of microphase separation on the elasticity of such networks?

These and many other questions are still awaiting answers.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Replica Model of the Network**

Consider a network prepared by random end-linking a solution of polydisperse linear chains by *f*-functional cross-links. Following the idea of *ϕ*<sup>4</sup> formulation of the excluded volume problem [56], we introduce *n* → 0 component field, *ϕ* (**xˆ**) in the replica space with components *ϕ<sup>i</sup>* (*i* = 1, ... , *n*), which is related to the monomer density in this space as

$$\rho(\mathfrak{k}) = \frac{1}{2}\varrho^2\left(\mathfrak{k}\right) = \frac{1}{2}\sum\_{i=1}^n \varphi\_i^2\left(\mathfrak{k}\right). \tag{A1}$$

Unlike the main text, in this Appendix the symbol *n* is used only to indicate the number of components of the field *ϕ*. The free energy *Fm* of the network in the replica system is described by Hamiltonian

$$H\left[\boldsymbol{\varrho}\right] = \int d\mathbf{\hat{x}} \left[\frac{1}{2}\mu\boldsymbol{\varrho}^2 + \frac{b^2}{2}\left(\boldsymbol{\widehat{\nabla}\boldsymbol{\varrho}}\right)^2 - \frac{z\_f}{f!}\boldsymbol{\varrho}\_1^f\right]$$

$$+ \sum\_k \int d\mathbf{x}^{(k)} \frac{\boldsymbol{\upsilon}^{(k)}}{8} \left[\prod\_{l\neq k} \int d\mathbf{x}^{(l)} \boldsymbol{\varrho}^2\right]^2. \tag{A2}$$

Here, ∇ is the gradient of the variables **xˆ**, *μ* is monomer chemical potential and *z <sup>f</sup>* is the fugacity of *f*-functional cross-links. The last sum in Equation (A2) describes the two-body monomer interaction in the initial system characterized by the excluded volume parameter *v*(0) and in final systems (*k* = 1, . . . , *m*) with the excluded volume parameters *v*(*k*) = *v*.

In the mean-field approximation, the steepest descent value of the field *ϕ*<sup>1</sup> is found from the extremum of Hamiltonian in Equation (A13). The field *ϕ*1(**xˆ**) has the meaning of the order parameter of the liquid-solid (sol-gel) phase transition. The steepest descent solution *ϕ*<sup>1</sup> (*ς*) > 0 depends on a single variable *ς*, Equation (6). At *ς* = 0, the differential Equation (8) for this function becomes algebraic, an explicit solution of the minimum conditions is found

$$\varphi\_1\left(0\right) = \left[\frac{(f-1)!}{\bar{N}z\_f}\right]^{1/(f-2)}\text{ \AA} \tag{A3}$$

In the limit *m* → 0, the dimensionless order parameter is determined by differential Equation (8). The free energy of the network depends only on *ϕ*<sup>1</sup> (0). Substituting Equation (A2) into Equation (5), we find the free energy of the network[8]:

$$F = V^{(0)} G \left[ \sum\_{\underline{a}} \frac{\lambda\_{\underline{a}}^2}{2} + \ln \left( a N^{1/2} \right) \right] \tag{A4}$$

with the elastic modulus

$$\frac{G}{kT} = \rho\_f^{(0)}\left(\frac{f}{2} - 1\right),\tag{A5}$$

where *ρ* (0) *<sup>f</sup>* is the cross-link concentration.

#### **Appendix B. Distribution of Virtual Chains**

According to Equation (19) the distribution function *q* (*s*) = *χ*˜(*f*−1) (*s*) of the inverse number of monomers *s* = 1/*n* can be presented as a convolution of *f* − 1 distribution functions *χ*˜ (*s* ) of random variables *s* = 1/*m*. The multiple convolution is defined recurrently as

$$\tilde{\chi}^{(k)}\left(\mathbf{s}\right) = \left(\tilde{\chi}^{(k-1)} \* \tilde{\chi}\right)\left(\mathbf{s}\right) \equiv \int\_0^s \tilde{\chi}^{(k-1)}\left(\mathbf{s} - \mathbf{s}'\right) \tilde{\chi}\left(\mathbf{s}'\right) d\mathbf{s}' \tag{A6}$$

The distribution function *χ*˜ (*s* ) of an individual branch of the tree is related to the distribution function *π* (*m*) of the inverse values *m* = 1/*s* by an equation similar to Equation (25):

$$
\tilde{\chi}\left(\mathbf{s}'\right) = \frac{1}{\left(\mathbf{s}'\right)^2} \pi\left(\frac{1}{\mathbf{s}'}\right) \tag{A7}
$$

Finally, the distribution function *π* (*m*) of the monomer numbers *m* = *N* + *n* of one branch is a convolution of distribution functions of the monomer numbers of the network strands *P* (*N*) and virtual chains *p* (*n*),

$$
\pi \left( m \right) = \left( P \ast p \right) \left( m \right) = \int\_0^m P \left( m - n \right) p \left( n \right) dn \tag{A8}
$$

Equations (25) and (A6)–(A8) completely determine the distribution function *p* (*n*) of the number of monomers *n* of virtual chains for a given distribution function *P* (*N*) of the number of monomers *N* of real network strands. The only problem is solving this system of nonlinear integral equations. Here the replica method comes to the rescue, see Appendix A.

#### **Appendix C. Free Energy of Entangled Networks**

The network free energy for given positions of attachment points **X**(*s*) of virtual chains to the affinely deformed background is

$$F\left[\mathbf{X}\right] = -\nu kT \sum\_{a} \ln \int D\mathbf{x}\_{a}(\mathbf{s}) e^{-H\_{a}\left[\mathbf{x}\_{a}\right]/kT} \,, \tag{A9}$$

$$\frac{H\_{\rm a} \left[ \chi\_{\rm a} \right]}{kT} = \int\_{0}^{N} ds \left\{ \frac{3\dot{\chi}\_{\rm a}^{2}}{2b^{2}} + \frac{3[\chi\_{\rm a}(s)/\lambda\_{\rm a} - X\_{\rm a}(s)]^{2}}{2b^{2}N\_{\rm c}^{2}} \right\}.\tag{A10}$$

The simplest way to calculate the above integral is to use the normal mode expansion of the chain trajectories, which satisfy the boundary conditions *xα*(*N*) = *λαXα*(*N*) and *xα*(0) = *λαXα*(0) (**X** (*s*) are the positions of attachment points of virtual chains in the undeformed network):

$$\mathbf{x}\_{\mathfrak{a}}(\mathbf{s}) = \lambda\_{\mathfrak{a}} R\_{\mathfrak{a}} \frac{\mathbf{s}}{N} + \sum\_{\omega} \tilde{\mathbf{x}}\_{\mathfrak{a}}(\omega) \,\boldsymbol{e}^{i\omega \mathbf{s}}, \qquad \omega = \frac{2\pi k}{N}, \quad k = \pm 1, \cdots \tag{A11}$$

where *R<sup>α</sup>* = *Xα*(*N*) − *X<sup>α</sup>* (0) and

$$\mathbf{X}\_{\mathfrak{d}}\left(\mathbf{s}\right) = \lambda\_{\mathfrak{d}}\mathbf{R}\_{\mathfrak{d}}\frac{\mathbf{s}}{N} + \sum\_{\omega} \tilde{\mathbf{X}}\_{\mathfrak{d}}\left(\omega\right) \mathbf{e}^{i\omega\mathbf{s}}, \qquad \tilde{\mathbf{x}}\_{\mathfrak{d}}\left(-\omega\right) = \tilde{\mathbf{x}}\_{\mathfrak{d}}^{\*}\left(\omega\right), \quad \tilde{\mathbf{X}}\_{\mathfrak{d}}\left(-\omega\right) = \tilde{\mathbf{X}}\_{\mathfrak{d}}^{\*}\left(\omega\right). \tag{A12}$$

The mode with *ω* = 0 is excluded from the summation because it corresponds to the translational displacement of the chain, which is prohibited in the network. Substituting the above expressions into Equation (A10) we get

$$\frac{H\_{\rm a} \left[ \mathbf{x}\_{\rm a} \right]}{kT} = \frac{3R\_{\rm a}^{2}}{2b^{2}N} \lambda\_{\rm a}^{2} + \frac{3}{2b^{2}} \sum\_{\omega \neq 0} \left[ \omega^{2} |\tilde{\mathbf{x}}\_{\rm a} \left( \omega \right)|^{2} + \frac{|\tilde{\mathbf{x}}\_{\rm a} \left( \omega \right) - \lambda\_{\rm a} \tilde{\mathbf{x}}\_{\rm a} \left( \omega \right)|^{2}}{N\_{\rm e}^{2} \lambda\_{\rm a}^{2}} \right] \tag{A13}$$

Changing the integration variables, *x*˜*<sup>i</sup>* (*ω*) = *u*˜ (*ω*) + *iv*˜(*ω*), we find

$$\frac{F[\mathbf{X}]}{\nu kT} = \sum\_{\underline{a}} \left[ \frac{3R\_{\underline{a}}^2}{2b^2 N} \lambda\_{\underline{a}}^2 + \frac{3\lambda\_{\underline{a}}^2}{2b^2} \sum\_{\omega > 0} \frac{\omega^2 |\breve{X}\_{\underline{a}}(\omega)|^2}{1 + \omega^2 N\_{\varepsilon}^2 \lambda\_{\underline{a}}^2} \right.$$

$$-\sum\_{\omega > 0} \ln \frac{2b^2 N\_{\varepsilon}^2 \lambda\_{\underline{a}}^2}{3N(1 + \omega^2 N\_{\varepsilon}^2 \lambda\_{\underline{a}}^2)} \right].\tag{A14}$$

Here sums are taken only over positive *ω* since the only positive frequency components *u*˜ (*ω*) and *v*˜(*ω*) can be taken as independent variables: the condition for the reality of the function *x*(*s*), Equation (A12), imposes restrictions *u*˜ (−*ω*) = *u*˜ (*ω*) and *v*˜(−*ω*) = −*v*˜(*ω*) on the real and imaginary parts of *x*˜ (*ω*).

This free energy should be averaged over positions of attachment points, **X**(*s*), with correlation functions

$$\overline{\tilde{X}\_{\mathfrak{A}}(\omega)}\,\tilde{X}\_{\mathfrak{F}}(-\omega) = \frac{b^2}{3} \left(\frac{1}{\omega^2} + N\_e^2\right)\delta\_{a\mathfrak{F}}\tag{A15}$$

which corresponds to randomly placed brash of chain with two virtual chains of *N*<sup>2</sup> *<sup>e</sup>* monomers

$$\overline{\left[\mathbf{X}\left(\mathbf{s}\right) - \mathbf{X}\left(\mathbf{s'}\right)\right]^2} = b^2 \left|\mathbf{s} - \mathbf{s'}\right| + 2N\_\ell^2 \tag{A16}$$

Averaging the free energy (A14) with the correlation function, (A15) and changing the integration variables *<sup>X</sup>*˜(*ω*) <sup>→</sup> (*R*, *<sup>p</sup>*˜ (*ω*), *<sup>q</sup>*˜(*ω*)) when integrating over **<sup>X</sup>** (*ω*) we find

$$\begin{split} \frac{\mathcal{F}}{\omega kT} &= \sum\_{\underline{a}} \left[ \frac{\lambda\_{\underline{a}}^2}{2} + \sum\_{\omega > 0} \frac{\lambda\_{\underline{a}}^2 - 1}{1 + \omega^2 N\_e^2 \lambda\_{\underline{a}}^2} \\ &+ \sum\_{\omega > 0} \ln \left( 1 + \frac{1}{N\_e^2 \lambda\_{\underline{a}}^2 \omega^2} \right) \right] + const. \end{split} \tag{A17}$$

Here all the terms which do not depend on the extension coefficients *λα* are absorbed in *Const*. The sum and the product can be expressed in terms of hyperbolic functions, see Equation (60).

#### **References**


© 2020 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **E**ff**ect of Fillers on the Recovery of Rubber Foam: From Theory to Applications**

#### **Thridsawan Prasopdee <sup>1</sup> and Wirasak Smitthipong 1,2,3,\***


Received: 24 October 2020; Accepted: 17 November 2020; Published: 19 November 2020

**Abstract:** Natural rubber foam (NRF) can be prepared from concentrated natural latex, providing specific characteristics such as density, compression strength, compression set, and so on, suitable for making shape-memory products. However, many customers require NRF products with a low compression set. This study aims to develop and prepare NRF to investigate its recoverability and other related characteristics by the addition of charcoal and silica fillers. The results showed that increasing filler loading increases physical and mechanical properties. The recoverability of NRF improves as silica increases, contrary to charcoal loading, due to the higher specific surface area of silica. Thermodynamic aspects showed that increasing filler loading increases the compression force (*F*) as well as the proportion of internal energy to the compression force (*F*u/*F*). The entropy (*S*) also increases with increasing filler loading, which is favorable for thermodynamic systems. The activation enthalpy (Δ*H*a) of the NRF with silica is higher than the control NRF, which is due to rubber–filler interactions created within the NRF. A thermodynamic concept of crosslinked rubber foam with filler is proposed. From theory to application, in this study, the NRF has better recoverability with silica loading.

**Keywords:** rubber foam; filler; charcoal; silica; compression set; recovery; thermodynamics

#### **1. Introduction**

Natural rubber (NR) is a biobased polymer usually applied in the rubber industry to manufacture products such as gloves, pillows, mattresses, tires, and so on. It is derived from the *Hevea brasiliensis* tree as colloidal latex, which presents other substances or nonrubber components, such as proteins, fatty acids, inorganic matters, and so on. This natural latex can be stabilized by the base chemical to maintain prolonged storage [1–4].

Many products, such as pillows and mattresses, are made from rubber foam prepared from concentrated natural latex, which provides specific characteristics. Rubber foams are generally porous and elastic with ventilated surfaces. These foams are made into lightweight products for comfort applications such as pillows and mattresses [5,6]. From the perspective of mechanical properties, rubber foam can either be soft or firm depending on the formula of the compounded latex. Concentrated natural latex is mixed with chemical agents, consisting of a blowing agent, vulcanizing agent, accelerators, an activator, antioxidants, a gelling agent, and so on, to form compounded latex [1,7]. All the chemical agents have to be approximately ground into the micron scale because they can be

easily mixed with the concentrated natural latex since the rubber particle is also present in the micron scale [8].

Rubber foam typically presents shape-memory polymer characteristics, which can be altered by useful properties from the external stimuli [9]. This ability enables responsive materials with form-fitting, actuation, and sensing characteristics in industries such as furniture, biomedicine, aerospace, and so on [10,11]. The viscoelastic properties of rubber foam play an important role in the shape-memory effect. However, many customers require natural rubber foam products that return to their initial form in minimal time or have better recoverability. One can enhance the elasticity of rubber foam to improve the recovery performance of rubber foam products.

In 2015, Bashir et al. [12] used eggshell powder (ESP) as a filler in NRF. The contribution of ESP to an increase in mass increases the density of NRF with increasing ESP loading. At low ESP loading, the tensile strength of the ESP-filled NRF initially drops because the filler is not enough to reinforce the NRF. The tensile strength increases after adding more filler loadings due to the reinforcement effect of ESP. The increment of ESP loading causes the ESP-filled NRF to lose recoverability, caused by nonelastic deformation, indicating the deformation of the hard phase from the increase of ESP loading. Another filler in NRF is rice husk powder (RHP). Ramasamy et al. [13] found that the recovery percentage decreases with the increase of RHP loading, whereas the compression strength is increased with increasing RHP loading.

Many previous works have shown that the recoverability of NRF is decreased with increasing filler loading. The main objective of this study is to enhance the recoverability of NRF; thus, we investigated a crosslinked NRF filled with commercial filler consisting of activated charcoal and silica, which are widely used in the rubber industry. This is a useful channel to better understand the relationship of the mechanical properties and thermodynamic theory of NRF, which can be applied for the recovery of rubber foam products.

#### **2. Materials and Methods**

#### *2.1. Materials*

The following materials were used: high-ammonia-concentrated natural latex (Num Rubber and Latex Co., Ltd., Trang, Thailand), potassium oleate (KO: 20% aqueous dispersion prepared by Thanodom Technology Co., Ltd., Bangkok, Thailand), sulfur (S: 50% aqueous dispersion prepared by Thanodom Technology Co., Ltd., Bangkok, Thailand), zinc diethyldithiocarbamate (ZDEC: 50% aqueous dispersion prepared by Thanodom Technology Co., Ltd., Bangkok, Thailand), zinc-2-mercaptobenzothiazole (ZMBT: 50% aqueous dispersion prepared by Thanodom Technology Co., Ltd., Bangkok, Thailand), Wingstay L (WingL: 50% aqueous dispersion prepared by Thanodom Technology Co., Ltd., Bangkok, Thailand), zinc oxide (ZnO: 28% aqueous dispersion prepared by Thanodom Technology Co., Ltd., Bangkok, Thailand), diphenylguanidine (DPG: 28% aqueous dispersion prepared by Thanodom Technology Co., Ltd., Bangkok, Thailand), sodium silicofluoride (SSF: 23% aqueous dispersion prepared by Thanodom Technology Co., Ltd., Bangkok, Thailand), toluene (AR grade, RCI Labscan Co., Ltd., Bangkok, Thailand), bamboo charcoal (BET 2–4 m2/g for a 40 μm charcoal particle size, CharcoalHome Co., Ltd., Bangkok, Thailand), precipitated silica (BET 170–190 m2/g for a 10–20 nm silica particle size, Extol Technology Co., Ltd., Dongguan, China).

#### *2.2. Preparation of Rubber Foams*

Rubber foams were prepared in nine types, as presented in Table 1. For each sample, high-ammonia-concentrated natural latex was weighed, and the ammonia was eliminated using a blender at 80 rpm for 1 min. Filler was later added and mixed for 1 min. The speed of the blender was increased to 160 rpm, and KO was subsequently added and mixed for 10 min until the volume of foam increased by approximately three times. The blending speed was reduced to 80 rpm. A group of chemicals (S, ZDEC, and ZMBT as vulcanizing chemicals and WingL) was added, and mixing continued

for 1 min. After that, another group of chemicals (ZnO and DPG) was added, and mixing continued for 1 min. SSF was added at last and mixed for 3 min. During the mixing of SSF, gel-forming was repeatedly tested until the gel point was almost reached. The foam was then poured into molds, and the lids were closed afterward. The samples were left at room temperature for 45 min. Finally, the samples were vulcanized in a hot air oven at 90 ◦C for 2 h, removed from the molds, washed, and dried at 70 ◦C.


**Table 1.** Composition of raw materials for the natural rubber foam (NRF) used in this study.

<sup>1</sup> Parts per hundred of rubber.

#### *2.3. Characterization*

Density was measured from the weight of the rubber foam sample and the measured volume of the rubber foam sample, calculated as follows:

$$\text{Density} = \frac{M}{V} \tag{1}$$

where *M* is the weight of the rubber foam sample (kg), and *V* is the volume of the rubber foam sample (m3).

The crosslinking density of the NRF was determined by the swelling method. Small pieces of the NRF were immersed into toluene to reach the equilibrium of swelling. The crosslinking density of rubber (ν) can be calculated according to the Flory–Rehner equation [14–16] as follows:

$$\nu = -\frac{1}{V\_s} \left[ \frac{\ln(1 - V\_\mathbf{r}) + V\_\mathbf{r} + \chi V\_\mathbf{r}^2}{V\_\mathbf{r}^\frac{1}{3} - \frac{V\_\mathbf{r}}{2}} \right] \tag{2}$$

where *V*<sup>s</sup> is the molar volume of toluene (106.9 cm3/mol), *V*<sup>r</sup> is the volume fraction of rubber in the swollen network, and χ is the parameter between the rubber and solvent interaction (0.43 + 0.05 *V*r). The NRF sample in toluene was kept in the dark at an ambient temperature. It was removed from the toluene and weighed every day until the equilibrium was reached (approximately 7 days). Finally, the sample was dried in an oven at 60 ◦C for 24 h and weighed again.

The functional group of the NRF sample was determined by attenuated total reflection–Fourier transform infrared spectroscopy (ATR–FTIR VERTEX 70, Bruker, Billerica, MA, USA). The NRF sample was put on a Ge crystal probe at 500–4000 cm<sup>−</sup>1.

Compression stress of the NRF was determined in 3 replications by a texture analyzer (TA.XT Plus, Stable Micro Systems, Godalming, Surrey, UK) using a platen probe with a diameter of 100 mm and a speed of 0.1 mm/s, adapted from ISO 3386. The NRF sample was 45 mm width × 45 mm length × 21.5 mm thickness, compressed at 75% from the foam surface at room temperature.

The compression set was conducted according to ISO 1856 Method C by measuring the height of the sample (*d*o). The sample was then compressed at 75 ± 4%height for 72 h at an ambient temperature. Finally, the sample was released from the compression for 30 min, the height of the sample was measured again (*d*r), and %compression set was calculated as follows:

$$\% \text{compression set} = \left[ \frac{d\_{\rm o} - d\_{\rm r}}{d\_{\rm o}} \times 100 \right] \tag{3}$$

Furthermore, %recovery of rubber foam can be calculated as follows:

$$\% \text{recovery} = 100 - \% \text{compression set} \tag{4}$$

The morphological properties of the NRF sample were characterized by scanning electron microscopy (SEM) analysis (FEI, Quanta 450 FEI, Eindhoven, Netherlands). The NRF sample was cut into small pieces and coated with gold in a sputter coater (Polaron Range SC7620, Quorum Technologies Ltd., Kent, UK).

For the thermodynamic aspects, the NRF was determined in 3 replications by a texture analyzer (TA.XT Plus, Stable Micro Systems, Godalming, Surrey, UK) using a platen probe with a diameter of 100 mm and speed of 0.1 mm/s. The sample was 45 mm width × 45 mm length × 21.5 mm thickness and prepared at the test temperature 10 min before starting the experiment. The sample was then compressed from 20% strain to 70% strain from the original foam shape at various temperatures (25, 35, 45, 55, and 65 ◦C). After that, the force–temperature relationship graph was plotted, and the obtained variables were used to calculate the compression force (*F*) due to the changes of internal energy (*F*u) and entropy (*F*s) associated with the deformation process.

Other thermodynamic aspects were studied by a temperature sweep of the NRF sample using dynamic mechanical analysis (DMA1, Mettler Toledo, Columbus, OH, USA). The NRF sample was cut into a sample of 7 mm width × 7 mm length × 7 mm thickness and tested at temperatures from 80 to 80 ◦C. The changes of Gibbs free energy (Δ*G*) and entropy (Δ*S*) were calculated. The activation enthalpy of the transition process ((Δ*H*a)avg) for the relaxation of the backbone motion of rubber chains is related to the area under the tan δ peak, which can be calculated by [17]:

$$(\Delta H\_{\text{a}})\_{\text{avg}} = \frac{(\ln E\_{\text{g}} - \ln E\_{\text{r}}) \pi RT\_{\text{g}}^2}{t\_{\text{A}}} \tag{5}$$

where *t*<sup>A</sup> is the area under the tan δ peak, *E*<sup>g</sup> is the storage modulus at the glassy state, *E*<sup>r</sup> is the storage modulus at the rubbery state, *R* is the gas constant (8.3145 J/mol·K), and *T*<sup>g</sup> is the glass transition temperature (K) of the NRF.

#### **3. Results and Discussion**

#### *3.1. Physical and Chemical Properties*

The density of the NRF samples increases with increasing filler loading, as presented in Figure 1. This is due to increasing filler loading, which causes an increase in the mass of the NRF with filler. Kudori and Ismail [18] found that foam density increases as the filler size decreases. The smaller filler size hinders pore formation and increases the continuous matrix amount. In the present study, silica filler (nanometer) is smaller than charcoal (micrometer). Therefore, NRF with silica loading exhibits a higher density than NRF with charcoal loading at a given filler concentration. Increasing charcoal loading barely affects the density, which may be because charcoal acts as a nucleating agent during the process of foam growth [19–21].

**Figure 1.** Effect of filler loading on the density of the NRF samples.

As presented in Figure 2, the crosslinking density of the NRF samples increases with the presence of filler, which may be due to the additional carbon–sulfur linkages formed by the chemical reaction between the rubber and filler [22]. Another reason is that the amplification of the deformation of rubber chains in the NRF with filler loading is more than the control NRF. Fillers in NRF extend rubber chains due to the interaction of rubber chains at the filler surface, i.e., some rubber chains may be occluded in the voids of the filler, causing the extension of rubber chains and leading to an increased crosslinking density. However, at the same loading of vulcanizing chemicals, increasing filler loading causes fewer differences in the crosslinking density.

**Figure 2.** Effect of filler loading on the crosslinking density of the NRF samples.

The chemical compositions of the control NRF and NRF with fillers were analyzed by ATR–FTIR (Figure 3). There is no significant difference in the functional groups of NRF [7,23]. There is almost no difference for the NRF with charcoal loading due to carbonization at high temperatures, which causes the charcoal powder to exhibit a hydrophobic nature [24]. However, there is a band growing at 1100 cm−<sup>1</sup> for the NRF filled with silica. This band corresponds to the vibration absorption of the silane group (Si–O–C) [25], present in the rubber network, which usually exhibits within the ranges 800–850 and 1100–1200 cm<sup>−</sup>1.

**Figure 3.** ATR–FTIR spectra of the NRF samples: (**a**,**b**) NRF with charcoal loading; (**c**,**d**) NRF with silica loading.

#### *3.2. Mechanical and Morphological Properties*

The control NRF and NRF with fillers were compressed up to 75% (Figure 4). The compression strength at maximum compression shows that the NRF with silica loading has higher compression strength than the NRF with charcoal loading at a given filler concentration. There are two different regions in the compression stress–strain curves of foam materials: elasticity at the low-strain region and solidity at the high-strain region [7]. Increasing filler loading increases the solidity or stiffness of the NRF at high strain, where the foam cells with each other, leading to the immediate increase of compression stress. There is also the stress-induced crystallization of rubber chains that affects the increase of compression stress at high strain. The addition of more than 2 phr of charcoal causes the foam to be sticky, explaining the unfavorable interaction within the foam structure. Although the crosslinking density of the NRF with various silica loadings is almost identical, the compression strength is significantly different. The better interaction within the foam structure is due to the smaller silica filler size, which has more specific surface areas compared to the charcoal filler. We found a linear relationship between compression strength and filler loading in both types of fillers. This relationship depends on the rubber–filler interaction, presented as:

$$\text{C}\_{\text{Ch}} = 0.285 \cdot \text{[Ch]} + 5.304 \tag{6}$$

$$\mathbf{C}\_{\rm Si} = 1.721 \mathbf{\dot{\omega}} \mathbf{\dot{\omega}} \mathbf{\dot{\omega}} + 4.649 \tag{7}$$

where *C*Ch is the compression strength of the NRF with charcoal loading (kPa), *C*Si is the compression strength of the NRF with silica loading (kPa), [Ch] is the concentration of charcoal filler (phr), and [Si] is the concentration of the silica filler (phr).

**Figure 4.** Compression stress of the NRFs: (**a**) NRF with charcoal loading; (**b**) NRF with silica loading.

The compression set describes the elastic behavior of the NRF, which relates to the material's recovery percentage. Figure 5 shows that increasing charcoal loading increases the compression set percentage and decreases the recovery percentage. As mentioned above, the addition of more than 2 phr of charcoal causes the foam to be sticky. Thus, when the NRF with more than 2 phr of charcoal loading is compressed at 75%height of its thickness for a long period (72 h), the ability to return to its original shape is decreased. On the contrary, increasing silica loading decreases the compression set percentage and increases the recovery percentage. Decreasing the compression set percentage indicates higher elasticity. Hence, the NRF with silica loading possesses higher elasticity than the NRF with charcoal loading. Microsized charcoal has been shown to behave like eggshell powder and rice husk powder in previous works [12,13], which decreased the percentage of NRF recovery when filler loading is increased and vice versa with NRF-filled nanosized silica. Therefore, we can propose the relationship between recoverability of NRF and filler concentration as the following polynomial equation:

$$\text{V}\% \text{R}\_{\text{ch}} = -0.3057 \cdot \text{[Ch]}^2 + 1.5206 \cdot \text{[Ch]} + 94.697 \tag{8}$$

$$\text{V}\,\%R\_{\text{Si}} = -0.1165\,\text{[Si]}^2 + 1.5602\,\text{[Si]} + 94.789\,\tag{9}$$

where %*R*Ch is the recoverability of the NRF with charcoal loading (%), %*R*Si is the recoverability of the NRF with silica loading (%), [Ch] is the concentration of charcoal filler (phr), and [Si] is the concentration of silica filler (phr).

**Figure 5.** Compression recovery properties of the NRFs: (**a**) %compression set of the NRF with different filler loadings; (**b**) %recovery of the NRF with different filler loadings.

The morphological properties of NRF were investigated by SEM. The micrographs indicated that all types of NRF contain a cellular structure that exhibits an interconnected network of open cells, as presented in Figure 6. Porosity analysis was determined by the ImageJ software by adjusting the threshold of the images. The white region corresponds to the pore shape, whereas the dark region corresponds to the open holes or pores (Figure 6). The interconnected porosity of these NRFs is an important parameter that affects the mechanical properties [7]. This result can be explained by the cell density value. The cell density (*d*cell) is calculated as follows [26]:

$$d\_{\rm cell} = \frac{3}{4\pi r^3} \left( 1 - \frac{\rho}{\rho\_s} \right) \tag{10}$$

where ρ is the foam density, ρ<sup>s</sup> is the solid phase density (NR 0.93 g/cm3), and *r* is the average cell radius.

**Figure 6.** SEM images of: (**a1**) control NRF; (**b1**) NRF/2 Ch; (**c1**) NRF/4 Ch; (**d1**) NRF/6 Ch; (**e1**) NRF/8 Ch; (**f1**) NRF/2 Si; (**g1**) NRF/4 Si; (**h1**) NRF/6 Si; (**i1**) NRF/8 Si. The SEM images with ImageJ analysis of: (**a2**) control NRF; (**b2**) NRF/2 Ch; (**c2**) NRF/4 Ch; (**d2**) NRF/6 Ch; (**e2**) NRF/8 Ch; (**f2**) NRF/2 Si; (**g2**) NRF/4 Si; (**h2**) NRF/6 Si; (**i2**) NRF/8 Si.

As presented in Table 2, increasing charcoal loading increases the average NRF pore size, decreasing the porosity and cell density. On the other hand, increasing silica loading decreases the average pore size, increasing the porosity and cell density. Although the density of the NRF increases with increasing filler loading, the cell density is more complicated. For charcoal loading, the density slightly increases or remains almost constant with the addition of more than 4 phr of charcoal. The cell density of the NRF with charcoal loading decreases and becomes almost constant with the addition of more than 4 phr of charcoal. Since charcoal can act as the nucleating agent, which can promote

foam growth, excess charcoal may stop acting as a filler and become a nucleating agent, resulting in an almost constant density with a larger average pore size and smaller porosity and cell density. For silica loading, the density increases as increasing silica loading, indicating the decreasing of the average pore size while the cell density increases.


**Table 2.** Average pore size, porosity, and cell density of the control NRF and NRF with filler loading.

#### *3.3. Thermodynamic Aspects*

Thermodynamic studies of the deformation process in uncrosslinked rubbers have already been discussed [27–29]. Most of the works showed the results of the temperature dependence of the stress in the extended state. In the present work, the mechanical compression properties of the crosslinked NRF samples are remarkable, especially for the %compression set and %recovery; thus, it is interesting to investigate the thermodynamic aspects. From the perspective of thermodynamics, the elasticity of rubber attributes to the changes in the conformations of rubber molecules from the unstrained molecules to the strained molecules. Such changes are related to the changes of internal energy and entropy associated with the deformation process as the following relationship [27,30,31]:

$$F\_{\parallel} = \left(\frac{\partial \mathcal{U}}{\partial \mathcal{L}}\right) - T\left(\frac{\partial \mathcal{S}}{\partial \mathcal{L}}\right) = \left.F\_{\mathbf{u}}\right| + \left.F\_{\mathbf{s}}\right|\tag{11}$$

$$F\_{\rm u} = \left(\frac{\partial \mathcal{U}}{\partial L}\right) \tag{12}$$

$$F\_8 = -T \left(\frac{\partial \mathcal{S}}{\partial L}\right) \tag{13}$$

where *F* is the compression force causing changes in the length of NRF (*L*), *U* is the internal energy of NRF, *T* is the temperature used in the experiment, and *S* is the entropy of NRF. When plotting the compression force graph as a function of the conducted temperature, the interception at 0 K is equal to *F*u, and the slope multiplied by the temperature is equal to *F*s.

To investigate the relationships between force (compression mode) and temperature, the NRF samples were compressed up to 20%strain, 30%strain, 40%strain, 50%strain, 60%strain, and 70%strain in the temperature controller chamber (at 298.15, 308.15, 318.15, 328.15, and 338.15 K). The relationships between compression force and temperature of the control NRF, NRF with 8 phr of charcoal, and NRF with 8 phr of silica are presented in Figures 7–9, respectively. The graphs of the other samples are presented in Figures S1–S6. The results reveal that the compression force to the sample increases with increasing %strain. At a given strain, the compression force decreases with increasing temperature. Moreover, the slope decreases at a higher strain due to a decrease of the rubber chains' degree of freedom in the NRF, which is unfavorable for entropy.

**Figure 7.** Force at a constant strain as a function of the temperature of the control NRF with a minimum of R<sup>2</sup> = 0.9 in each strain.

**Figure 8.** Force at a constant strain as a function of the temperature of the NRF with 8 phr of charcoal (NRF/8 Ch) with a minimum of R2 = 0.9 in each strain.

**Figure 9.** Force at a constant strain as a function of the temperature of the NRF with 8 phr of silica (NRF/8 Si) with a minimum of R2 = 0.9 in each strain.

Figures 7–9 show the values of the compression force (*F*) and relative internal energy contributing to the compression force (*F*u/*F*) at 298.15 K can be calculated as indicated in Table 3. The results of the other samples are presented in Table S1. The values of *F*u and *F* increase with increasing compression strain, whereas the values of *F*u/*F* decrease. Since the internal energy (*U*) is varied by the compression force (*F*), the internal energy increases with increasing compression force. The entropy (*S*) can be varied by the length of the NRF, indicating the degree of freedom of rubber chains during the compression process. The compression causes a reduction in the length of the NRF, leading to a decrease in the rubber chains' degree of freedom. Thus, the entropy of compressed NRF is also reduced.


**Table 3.** Compression strain, compression limit, *F*u, *F*, and *F*u/*F* values of the control NRF and NRF with various fillers at 298.15 K.

The *F*u/*F* values of uncrosslinked rubber in the extension mode are typically in the range of 0.1–0.2 [27]. In this work, the *F*u/*F* values of the crosslinked NRF in the compression mode are in the range of 0.6–0.8, which are approximately three times higher than those of the literature review. The difference in *F*u/*F* values could come from the material structures (uncrosslinked rubber vs. crosslinked rubbers) and the test methods (extension mode vs. compression mode). The NRFs with fillers have higher *F*u/*F* values than the control NRF at a given strain level, possibly explained by the interaction of rubber and filler, which promotes changes of entropy in the deformation process. The slope direction of the *F*u/*F* values of the NRF with charcoal and NRF with silica are different (Figure 10). This is due to the control NRF and NRF with charcoal loading possess different degrees of freedom of rubber chains at different compression limits, and lower compression limit leads to lower *F*u/*F* values. However, the NRFs with silica loading possess a similar degree of freedom of rubber chains at different compression limits. This indicates that the stability of the degree of freedom of rubber chains at different compression strains or compression limits (λ) is related to the high mechanical property of the NRF with silica loading. Although the compression strength of the NRF with filler increased with the increment of filler, the ratio of *F*u/*F* indicates that the addition of filler affects the mechanical properties in the aspect of thermodynamics.

**Figure 10.** Relative internal energy contribution to the compression force (*F*u/*F*): (**a**) NRF with charcoal loading compared to control NRF; (**b**) NRF with silica loading compared to control NRF.

We also investigated the change in Gibbs free energy (Δ*G*) and entropy (Δ*S*) in the NRF with fillers compared to the control NRF. These thermodynamic parameters can be calculated by the Flory–Huggins equation and statistical theory of rubber elasticity as follows [17,32]:

$$
\Delta G = \left. RT \right[ \ln(1 - V\_{\rm r}) + V\_{\rm r} + V\_{\rm r}^2 \chi] \tag{14}
$$

$$
\Delta S = -\frac{\Delta G}{T} \tag{15}
$$

where *R* is the gas constant (8.3145 J/mol·K), and *T* is the test temperature (300.15 K).

From the perspective of the crosslinking density, Δ*G* and Δ*S* are shown in Table 4. The volume fraction of rubber (*V*r) with fillers is higher than the control NRF. The swelling behavior of the NRF with various filler loadings decreases with increasing filler loading, indicating that the filler enhances the rubber swelling resistance against the penetration of the solvent. Since the filler is the hard phase, which is impermeable to solvent molecules, there must be a higher interaction between phases and more rubber chains attached to the filler surface. Hence, the swelling ability of NR is reduced while increasing the volume fraction of rubber [17].

**Table 4.** Calculated parameters from crosslinking density results for the control NRF and NRF with various fillers.


Table 4 shows that all samples present a negative Δ*G*, which decreases with increasing filler and is a favorable spontaneous system. This is due to the restriction in the ability of the rubber chain motion in the presence of filler, resulting in a decrease in the Gibbs free energy, which can be attributed to good compatibility between polymer and filler [17,32]. Δ*S* increases with increasing filler loading, which is favorable in thermodynamics. This result is in good agreement with the result of the *F*u/*F* value.

Based on the dynamic mechanical properties of the sample, the storage modulus (*E'*) and tan δ results of the NRF with filler loading were determined by temperature sweep using dynamic mechanical analysis or DMA. The results are presented in Figure 11. The addition of filler, both charcoal and silica, affects the storage modulus and tan δ.

**Figure 11.** Storage modulus (*E'*) and tan δ as a function of the temperature of the NRF with various fillers: (**a**) storage modulus of the NRF with charcoal loading compared to the control NRF; (**b**) tan δ of the NRF with charcoal loading compared to the control NRF; (**c**) storage modulus of the NRF with silica loading compared to the control NRF; (**d**) tan δ of the NRF with silica loading compared to the control NRF.

The DMA results presented in Table 5 reveal that the storage modulus in both the glassy state and rubbery state of the NRF with filler loading is higher than the control NRF. The addition of filler decreases the free volume within the foam, which causes more rigidity, resulting in a higher storage modulus in the glassy state [33]. The storage modulus in the rubbery state of the control NRF is lower than the NRF with filler loading. This is due to the effect of the filler on the relaxation time of rubber chains. Increasing filler loading increases the volume fraction of filler (Δ*V*f), which causes a higher stress relaxation rate of rubber molecules where they require more time to unload the applied force [17,33]. This affects the degree of freedom of the rubber molecules to be more pronounced, i.e., when there is a greater number of interactions between the rubber chains and filler, the stress relaxation rate is increased, resulting in an increase in entropy [33,34]. Therefore, we can propose a model of the control NRF (Figure 12a) compared to the NRF with filler loading (Figure 12b). The thermodynamic meaning of this work can be explained as follows: the change in entropy (Δ*S*) of the NRF with filler loading is more pronounced compared to the control NRF due to the stress relaxation rate of rubber

chains from the amplification of chain deformation between the rubber and filler interaction, as shown in Equation (16).

$$
\Delta S = \frac{\Delta V\_{\text{f}} \cdot \Delta S \text{t}}{T} \tag{16}
$$

where Δ*S* is the change of entropy, Δ*V*<sup>f</sup> is the change of volume fraction of filler, Δ*S*t is the change of stress relaxation rate of rubber molecules, and *T* is the temperature.

**Table 5.** Obtained parameters by the dynamic mechanical analysis (DMA) test for the control NRF and NRF with various fillers.


**Figure 12.** Model of the NRF with rubber chains (blue spot) and fillers (black spot): (**a**) control NRF; (**b**) NRF with filler, which increases the volume fraction of filler (Δ*V*f) induces more stress relaxation (Δ*S*t), resulting in a more pronounced entropy (Δ*S*).

Table 5 also presents the *T*<sup>g</sup> value or peak of tan δ of the NRF with various fillers. The addition of filler causes this value to shift toward higher temperatures when compared to the control NRF. The shift of the *T*<sup>g</sup> value toward the higher temperatures indicates ionic and hydrogen bonding interactions between the rubber chains and filler [17]. However, the nonpolar charcoal might not disperse well in the concentrated natural latex or agglomerate and, instead, form filler–filler networks within the foam. This may cause a synergy effect where the filler–filler networks might defeat the movement of the free chains of rubber. Therefore, the addition of charcoal affects higher hysteresis with increasing tan δ max and *t*A, resulting in lower activation enthalpy (Δ*H*a) than the control NRF.

At the same time, NRF is well-reinforced with silica. The rubber chains within the NRF with silica are hindered to freely move, and there are interactions between rubber–filler within the NRF. Thus, it has a higher activation enthalpy than the control NRF. The tan δ max of the NRF with filler has a higher value than control NRF refers to more dissipation energy of the NRF in the existence of filler. The values of Δ*H*<sup>a</sup> in this work are in the same order as in the works of Sadeghi Ghari and Jalali-Arani [17].

#### **4. Conclusions**

In this work, natural rubber foam (NRF) was prepared in two conditions: NRF with charcoal loading and NRF with silica loading. The results showed that increasing filler loading increases the density and mechanical properties of rubber foam. Since rubber chains may be occluded in the voids of filler causing the expansion of rubber chains, leading to increasing crosslinking density, somehow, at the same loading of vulcanizing chemicals, increasing filler loading is affected less in the crosslinking density.

Increasing filler loading increases the compression stress of NRF. The compression strength of the NRF with silica loading is higher than NFR with charcoal loading due to the better interaction within the foam structure caused by the smaller silica size, which presents a more specific surface area compared to charcoal filler. Since charcoal can act as a nucleating agent and promote foam growth, the excess of charcoal might change from being filler and become the nucleating agent, resulting in a larger average pore size and smaller porosity and cell density. Increased silica loading results in a decrease in the average pore size while the cell density increases. Thermodynamic aspects showed the relationships between the force and temperature of the NRF samples. The compression force (*F*) and internal energy force (*F*u) values of the NRF samples increase with increasing compression strain. The NRF with various fillers has higher *F*u/*F* values than the control NRF. The slope direction of the *F*u/*F* value of the NRF with charcoal and NRF with silica are different, which comes from the different degrees of freedom of rubber chains in the NRF samples.

The swelling behavior of the NRF with filler loading decreases with increasing filler loading compared to the control NRF. The Δ*G* decreases while Δ*S* increases with increasing filler loading, which demonstrates a favorable thermodynamic system. The Δ*H*<sup>a</sup> of the control NRF is lower than NRF with silica due to the movement limitation of rubber chains, whereas NRF filled with charcoal is more complicated. Increasing filler loading increases the volume fraction of filler, causing a higher stress relaxation rate due to the attempt to relax itself of molecules; therefore, the entropy is more pronounced. All the results from the theory to applications indicate that NRF, which normally behaves like a shape-memory material, can be developed into an anti-shape-memory material by the addition of silica loading, which is favorable in thermodynamics.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4360/12/11/2745/s1, Figure S1. Force at a constant strain as a function of the temperature of the NRF with 2 phr of charcoal (NRF/2 Ch) with a minimum of R<sup>2</sup> = 0.9 in each strain. Figure S2. Force at a constant strain as a function of the temperature of the NRF with 4 phr of charcoal (NRF/4 Ch) with a minimum of R<sup>2</sup> = 0.9 in each strain. Figure S3. Force at a constant strain as a function of the temperature of the NRF with 6 phr of charcoal (NRF/6 Ch) with a minimum of R<sup>2</sup> = 0.9 in each strain. Figure S4. Force at a constant strain as a function of the temperature of the NRF with 2 phr of silica (NRF/2 Si) with a minimum of R<sup>2</sup> = 0.9 in each strain. Figure S5. Force at a constant strain as a function of the temperature of the NRF with 4 phr of silica (NRF/4 Si) with a minimum of R<sup>2</sup> = 0.9 in each strain. Figure S6. Force at a constant strain as a function of the temperature of the NRF with 6 hr of silica (NRF/6 Si) with a minimum of R2 = 0.9 in each strain. Table S1. Compression strain, compression limit, *F*u, *F*, and *F*u/*F* values of the control NRF and NRF with various fillers at 298.15 K.

**Author Contributions:** Conceptualization, W.S.; methodology, W.S.; software, T.P.; validation, W.S.; formal analysis, T.P.; investigation, W.S. and T.P.; resources, W.S.; data curation, W.S. and T.P.; writing—original draft preparation, T.P.; writing—review and editing, W.S.; visualization, W.S. and T.P.; supervision, W.S.; project administration, W.S.; funding acquisition, W.S. Both authors have read and agreed to the published version of the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** This research was supported by a graduate study development scholarship from the National Research Council of Thailand as of the 2020 fiscal year. This research was also supported by the Specialized Center of Rubber and Polymer Materials in Agriculture and Industry (RPM), Faculty of Science, Kasetsart University, Bangkok, Thailand.

**Conflicts of Interest:** The authors declare no conflict of interest in this research. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Adherence Kinetics of a PDMS Gripper with Inherent Surface Tackiness**

#### **Umut D. Çakmak 1,\*, Michael Fischlschweiger 2, Ingrid Graz <sup>3</sup> and Zoltán Major <sup>1</sup>**


Received: 18 September 2020; Accepted: 19 October 2020; Published: 22 October 2020

**Abstract:** Damage and fiber misalignment of woven fabrics during discontinuous polymer processing remain challenging. To overcome these obstacles, a promising switchable elastomeric adherence gripper is introduced here. The inherent surface tackiness is utilized for picking and placing large sheets. Due to the elastomer's viscoelastic material behavior, the surface properties depend on loading speed and temperature. Different peeling speeds result in different adherence strength of an interface between the gripper and the substrate. This feature was studied in a carefully designed experimental test set-up including dynamic thermomechanical, as well as dynamic mechanical compression analyses, and adherence tests. Special emphases were given to the analyses of the applicability as well as the limitation of the viscoelastic gripper and the empirically modeling of the gripper's pulling speed-dependent adherence characteristic. Two formulations of poly(dimethylsiloxane) (PDMS) with different hardnesses were prepared and analyzed in terms of their applicability as gripper. The main insights of the analyses are that the frequency dependency of the loss factor tanδ is of particular importance for the application along with the inherent surface tackiness and the low sensitivity of the storage modulus to pulling speed variations. The PDMS-soft material formulation exhibits the ideal material behavior for an adhesive gripper. Its tanδ varies within the application relevant loading speeds between 0.1 and 0.55; while the PDMS-hard formulation reveals a narrower tanδ range between 0.09 and 0.19. Furthermore, an empirical model of the pulling speed-dependent strain energy release rate G(v) was derived based on the experimental data of the viscoelastic characterizations and the probe tack tests. The proposed model can be utilized to predict the maximum mass (weight-force) of an object that can be lifted by the gripper

**Keywords:** silicon rubber (PDMS); dynamic thermomechanical analyses; storage modulus; mechanical loss factor; viscoelastic gripper

#### **1. Introduction**

The picking and placing of limp solids (textile, woven fabrics, soft, and flexible electronics) during assembling in manufacturing processes remain challenging. The transfer to the desired position has to be achieved without damaging or misaligning the limp solid (i.e., fiber unravel, pull out, contamination by the gripper, etc.). A variety of technologies are available to overcome this challenge including vacuum suction, ingressive and adhesive grippers. The first two methods directly interact with the substrate's surface to be transferred and deform it to some extent: vacuum suction results in a local lift up of the limp substrate and ingressive picking in penetrating as well as interlocking

with fibers. An adhesive gripper, on the other hand, adheres on the surface, and the substrate itself is neither deformed nor misaligned. The adhesive part is usually an elastomer (e.g., copolymers and plasticized elastomers) with permanently tacky surface [1,2]. A major drawback is that the adherence cannot be switched off and so the release of a substrate is critical for a textile or soft electronic part. The modification of surface properties and specifically the enhancement of the adsorption abilities (e.g., for catalysts) [3–7] and can help to achieve switchable grippers. Alternatively, transfer printing is utilized in the assembly of micro/nanofabrication, where a silicon elastomer stamp with a viscoelastic surface behavior is used to deliver a part from a donor to a receiver substrate (see, e.g., Cheng et al. (2012) [8]). At the interface of two elastic materials the adhesion force is rather a constant than tunable [9]. However, the adhesion can depend on the peeling speed (high → lifting and low → release), the mechanical loading protocol (directional shearing induced delamination), temperature (laser pulse) or could be controlled by the surface relief structure as it was pointed out by Cheng et al. (2012) [8], Li et al. (2012) [10], and Chen et al. (2013) [9]. Extensive research efforts were presented to exploit the switchable adhesion in kinetically controlled [9,11,12], shear-assisted [8,13], direction-controlled [14], laser-driven non-contact [10,15], and microstructure enabled transfer printing [16–19]. Furthermore, in robotics various designs for gripper are presented in order to pick up and handle arbitrary objects (see, e.g., Brown et al. (2010) [20]). This influence of external physical loadings on the elastomeric gripper is exploited in order to achieve a controlled picking and placing of a substrate. The pronounced viscoelastic behavior of the gripper with a high loading rate sensitivity of the mechanical properties is therefore important.

However, up to now no detailed study on the viscoelastic characteristics of an elastomeric adhesive gripper is presented. A throughout understanding of the interplay between loading rate sensitivity of the gripper's bulk property and the inherent surface tackiness is of particular necessity to optimize the performance of an adhesive gripper. This paper presents an approach to gain insights to the temperature and loading rate dependent storage and dissipative (loss) properties by dynamic thermomechanical analyses of adhesive grippers as well as probe tack tests. The grippers under investigation are made of two different formulations of poly(dimethylsiloxane). Two woven fabrics (limp substrate) and two woven reinforced thermoplastic matrix composites (stiff substrate) are used to examine the gripper capabilities by the probe tack test. In "Material and Specimens", the material selection is briefly discussed and the specimen geometries for the mechanical characterization ("Experimental") presented. The loading rate dependent bulk and interface properties are modeled by simple linear viscoelastic theory assumptions and linear elastic fracture mechanics approach ("Experimental"). Based on the comprehensive experimental characterization the main mechanical features required for a viscoelastic adhesive gripper is summarized and presented in "Results and Discussion".

#### **2. Materials and Specimens**

The silicone rubbers (poly(dimethylsiloxane); PDMS) under investigation were mixed and cured with two different base polymer:catalyst ratios (PDMS-soft → 20:1 and PDMS-hard → 10:1). The volume fraction of catalyst, curing temperature, and time to achieve a solid material are thereby critical, whose surface exhibits tackiness [21] high enough to pick up fabrics from a stack; the higher the tackiness the higher is the risk of contamination in a repetitive assembly process. Details about the fabrication of the material can be found in Cakmak et al. (2014a) [22].

Two categories of specimens were prepared. One category was dedicated for the general mechanical characterization of PDMS by dynamic thermomechanical analysis (DTMA); the other was for component tests under an application-related measurement protocol as well as dynamic mechanical analysis (DMA). DTMA was performed with the "barbell" specimen and the component tests were investigated with a cylindrical stamp. Figure 1 shows the specimen geometries and the respective dimensions.

**Figure 1.** Specimen geometries: (**a**) barbell specimen and (**b**) gripper stamp specimen. Dimensions are in mm; R: radius.

In the following section, the experimental procedure ("Experimental") for the characterization of the employed PDMS is presented. Furthermore, the design and the application-related experiments of the gripper stamp are reported there.

#### **3. Experimental**

After the initiation of the idea to adopt the rate and delamination loading-dependent (switchable) adherence from transfer printing to a large-scale application of limp solids (fabrics), some promising preliminary experiments were performed. Motivated by these preliminary results, an experimental scheme limited to mechanical characterizations only was defined. All experiments were performed with an electro-mechanic/dynamic actuator (TestBench, Bose Corp., ElectroForce Systems Group, MN, USA).

#### *3.1. Dynamic (Thermo-) Mechanical Analysis (D(T)MA)*

In order to gain more information concerning the viscoelastic nature of the PDMS, D(T)MA experiments were performed with the barbell specimen as well as the gripper stamp. A classical D(T)MA test procedure was defined including frequency dependent experiments from 1 Hz to 16 Hz under isothermal condition at seven different temperatures (−30 to +30 ◦C). The loading mode was uniaxial tensile with a sinusoidal excitation of a mean strain level of 1% and a dynamic amplitude of 0.5%. A similar test procedure was performed with the gripper stamp in order to determine the viscoelastic behavior of the PDMS-soft material under component relevant condition. Component relevant condition means, thereby, a uniaxial compression dynamic mechanical analysis (DMA) just before delamination from the substrate. Compression DMA was examined at room temperature (22 ◦C) with a mean level of −0.5 mm and amplitude of 0.5 mm excitation from 0.1 Hz to 47 Hz. The data of the D(T)MA experiments were analyzed in WinTest software (Bose Corp., ElectroForce Systems Group, MN, USA) and the storage (*E'*) as well as the transient (*E"* and tanδ = *E"*/*E'*) mechanical material properties were exported for further analyses. The complex modulus E\* is given by Equation (1) and can be modeled by the well-known Prony series in the frequency (ω = 2πf) domain [23]. In the series

*E*<sup>0</sup> refer to the instantaneous modulus, τ*<sup>i</sup>* is the relaxation time and *gi* = (*Ei*−1−*Ei*)/*E*<sup>0</sup> the Prony series coefficient corresponding to the relaxation time.

$$E^\*(\omega) = E'(\omega) + j \cdot E'^\prime(\omega) = E\_0 \cdot \left(1 - \sum\_{i} \frac{\mathcal{G}i}{1 + \left(\tau\_i \cdot a\_T \omega\right)^2} + j \cdot \sum\_{i} \frac{\mathcal{G}i \tau\_i \mu\_T \omega}{1 + \left(\tau\_i \cdot a\_T \omega\right)^2}\right) \tag{1}$$

If the material's thermorheological behavior is simple, then the time–temperature superposition can be applied and resulting in master curve at a reference temperature [22,24]. The temperature dependent shift factor aT can be given by the well-known WLF function [25] at a certain temperature *T* where *T*g is the glass transition temperature:

$$\left| \log\_{10}(a\_T) \right| = 17.44 \cdot \frac{T - T\_{\mathcal{S}}}{51.6 + T - T\_{\mathcal{S}}} \tag{2}$$

By inverse Laplace transformation of *E*\* into the time-domain, the relaxation modulus in time is obtained: 

$$E(t) = E\_0 \cdot \left(1 - \sum\_{i} g\_{i\cdot} \left(1 - e^{-t/\tau\_i} \right) \right) = E\_{\infty} \cdot \frac{1 - \sum\_{i} g\_{i\cdot} \left(1 - e^{-t/\tau\_i} \right)}{1 - \sum\_{i} g\_{i}} \tag{3}$$

where *E*<sup>∞</sup> = *E*<sup>0</sup> (1 − i gi ) is the Young's modulus after infinite long time.

For further analyses of the rate-dependent adhesion, Equation (2) will be more convenient. From an experimental point of view, the DMA tests are preferable, as the storage and the transient viscoelastic behaviors are determined at once.

#### *3.2. Probe Tack Test*

The gripper stamp made of PDMS-soft was investigated in order to determine the rate dependent adherence properties for different substrate surfaces. Basically, the adherence test was a tensile delamination test (see Figure 2) at various pulling speeds (0.1 mm/s up to 100 mm/s) in accordance to the surface tack measurement set-up of Çakmak et al. (2011) [21]. The stamp was pressed onto the substrate, held for 10 s at −10 N, followed by the pulling and delamination from the substrate surface, while the force F was measured. In a real application, the stamp will lift up the substrate, but here the adherence was of particular interest and therefore the substrate was fixed. The substrates (100 <sup>×</sup> 100) mm2 were attached to the support by means of double-side adhesive tape. To avoid inertia effects at high pulling speeds, the force was measured at the fixed side of the test setup. The maximum force corresponds to the fully delaminated stamp and is the limiting factor in-service for the considered pulling speed. If components with higher weight forces are applied, the stamp is not capable to pick them up. Application-relevant substrates were investigated and comprised dry carbon and glass fiber woven fabrics as well as poly(amide) matrix organo sheets reinforced with carbon and glass fiber weave.

The delamination process is assumed to be driven mainly by the crack propagation from the edge of the stamp when a tractive force *F* was applied. For the sake of completeness, it should be mentioned that the other mechanisms would be the crack propagation starting from inside at some interfacial defects and decohesion at theoretical contact strength [26]. The real delamination process will be rather a combination of the aforementioned mechanisms. However, visual observations of the experiments revealed that the crack initiation and propagation at the edge was more dominant. From fracture mechanical point of view, the stress singularity at the crack tip can be sufficiently described by the stress intensity factor *K*I. When the stamp is fully propagated along the interface, the crack length is equal to the radius *R* of the stamp and *K*<sup>I</sup> is given by

$$K\_I = \frac{F}{\pi \cdot R^2} \cdot \sqrt{\pi \cdot R} \cdot f \tag{4}$$

where *F* is the tractive force (*F*/π*R*<sup>2</sup> is the tractive stress σ) and *f* is the geometry coefficient and for a pillar shaped geometry *f* = <sup>1</sup> <sup>2</sup> . The energy needed to create an area surface during a plain strain loading is connected to the stress intensity factor KI by

$$G = \frac{K\_I^2 \cdot \left(1 - \nu^2\right)}{2 \cdot E} \tag{5}$$

where ν is the Poisson's ratio and *E* is the Young's modulus. Since the material under investigation can be considered as incompressible (ν = 0.5) and by combining Equation (4) with (5), the following relation can be found.

**Figure 2.** Component test set-up.

$$\mathbf{G}\_{\mathbf{c}}(\mathbf{v}, t\_{\mathbf{s}}(\mathbf{v})) = \frac{\mathbf{3} \cdot \sigma(\mathbf{v})^2 \cdot \pi \cdot \mathbf{R}}{\mathbf{32} \cdot \mathbf{E}(t\_{\mathbf{s}})} \tag{6}$$

Equation (6) shows the critical energy release rate per area depending on the pulling speed v and the separation time *t*s. The separation time is in turn a function of the pulling speed, meaning that the higher v, the lower *t*s. This dependency is mainly driven by the viscoelastic nature of the stamp material and is also examined experimentally. To determine the actual value of the Young's Modulus *E*(*t*s), knowledge of the separation time is crucial. The pulling velocity-dependent tractive stress σ(v) has to be obtained by the component tests and is mainly determined by the surface energies of each solid in contact and the surface interaction energy of them in direct contact (cf. Duprè energy of adhesion). It is expected that the investigated substrates have a different trend in terms of measured tractive stresses. However, the velocity dependency of σ(v) is a respond given by the stamp's viscoelastic nature and can be described by the following square root function sufficiently enough.

$$
\sigma(\mathbf{v}) = \sigma\_0 \cdot \sqrt{1 + \left(\frac{\mathbf{v}}{\mathbf{v}\_0}\right)^n} \tag{7}
$$

When the square root function is determined by simply fitting the measured data and combining Equations (3) and (7) with (6) lead to the pulling speed-dependent critical energy release rate:

$$\mathbf{G}\_{\mathbf{c}}(\mathbf{v}, t\_{\mathbf{s}}(\mathbf{v})) = \mathbf{G}\_{0} \cdot \frac{1 + \left(\frac{\mathbf{v}}{\mathbf{v}\_{0}}\right)^{\mathrm{n}} \cdot \left(1 - \sum\_{i}^{3} \mathbf{g}\_{i}\right)}{1 - \sum\_{i}^{3} \mathbf{g}\_{i} \cdot (1 - e^{-t\_{\mathbf{s}}/\tau\_{r}})} \tag{8}$$

where *G*<sup>0</sup> is the critical energy release rate near zero pulling speed v and infinite separation time *ts*, *v*<sup>0</sup> is the reference pulling speed and n is the scaling exponent. Equation (8) is in accordance with the empirical form shown and utilized by Shull (2002) [27] and Feng et al. (2007) [11]. The modification is related to the time dependent Young's modulus and so the Prony series is incorporated. Parameter identification of the function in Equation (8) is carried out based on the experimental data of the component tests as well as DMA. The model will be used to predict the critical energy release rate for the desired in-service pulling speed.

#### **4. Results and Discussion**

As was mentioned in the "Materials and Specimens" section, two different formulations were considered for the application of the gripper stamp. These two materials will be herein after referred to as PDMS-soft and PDMS-hard. Basically, only the PDMS-soft material was suitable for the use as gripper; however, basic viscoelastic characterization (D(T)MA) was performed for both of the materials in order to show the range of possible mechanical behavior. Only PDMS-soft was investigated with the component testing procedure.

#### *4.1. Dynamic (Thermo-) Mechanical Analysis (D(T)MA)*

Figure 3 shows the DTMA results of the investigated materials. Each color corresponds to the results of isothermal testing at various frequencies, which are indicated as single points.

**Figure 3.** Loss factor tanδ over storage modulus *E'* characteristic of the investigated materials at various temperatures.

This diagram is very convenient for analyzing the trend of the inherent respond time's temperature dependency as a result of the external mechanical loading. If all examined data reveal that tanδ is a unique function of *E'*, then the respond times (here the relaxation times) are equally temperature-dependent and time–temperature superposition (TTSP) is applicable (cf. Çakmak and Major (2014b) [28] and Tschoegl et al. (2002) [24]). Besides that, PDMS-soft formulation's characteristic is located within the diagram rather at the top right than the trend of PDMS-hard as expected. Important is that for both cases the TTSP can be applied, shift factors determined to construct the master curves at reference temperatures and these factors could be modeled with WLF function, among others.This is illustrated in Figure 4 for PDMS-soft (the candidate material for the application).

**Figure 4.** Panels (**a**,**b**) show the temperature and frequency dependent storage modulus and loss factor, respectively. In panel (**c**), the obtained time-temperature shift factors aT are presented as points and the WLF function as dashed line.

If the shift factor function (WLF) is applicable, then the temperature influence on the contact mechanical behavior can be easily predicted (e.g., Feng et al. (2007) [11]). This is of particular interest when the environmental conditions in-service of the gripper (i.e., during manufacturing of products) are fluctuating.

Figure 5 demonstrates the compression DMA results as before illustrated for the DTMA results. Now the gripper stamps made of the PDMS-formulations are investigated under contact and mechanically excited before delamination occurs. As mentioned earlier, the surface tackiness is an important requirement for the gripper stamp. However, the analyses of the viscoelastic behavior reveal that the loss factor characteristics are equally significant. It is believed that a high frequency (loading rate) dependency is essential to exploit the rate dependent adhesion. If the material has almost no rate-dependent tanδ, which is the case for PDMS-hard formulation, then the material will not be

suitable as gripper with rate dependent adhesion. The storage modulus and its rate dependency are of minor importance. It should be low and remain low with increasing loading rate.

**Figure 5.** Loss factor tanδ over storage modulus *E'* characteristic of the investigated formulation at 20 ◦C under compression DMA loading procedure.

The storage modulus of the PDMS-soft material is modeled with the 3rd order Prony series (cf. Equation (1)) and the respective Prony parameters are listed in Table 1. Figure 6 shows the experimental data of the frequency dependent storage modulus and the mentioned fit curve.


**Table 1.** Prony parameters of PDMS-soft formulation.

**Figure 6.** Linear plot of the frequency (ω = 2π)-dependent storage modulus *E'* of PDMS-soft formulation (black points) with the 3rd order Prony series fit of the data (cf. Table 1)**.**

With the Prony series, the time dependent storage modulus needed for computing and modeling the critical energy release rate is obtained.

#### *4.2. Probe Tack Test*

From the tack tests the tractive stresses are evaluated, which are measured as the maximum stresses before delamination at various pulling speeds and are shown in Figure 7. The experimental data confirm the assumption of the square root relation between the stresses and the pulling speeds shown in Equation (7). This function's fit parameters for the investigated substrates are given in Table 2.

**Figure 7.** Measured tractive stresses σ<sup>t</sup> between the PDMS-soft gripper stamp and the investigated substrates at various pulling speeds (points) with the fit function (Equation (7)) as dashed lines.

**Table 2.** Fit parameter of the tractive stress and the corresponding calculated G0.


With *E*<sup>∞</sup> and σ0, the critical energy release rate G0 near zero pulling speed and infinite separation time can be calculated. Table 2 shows the computed values and it can be observed, that the fabrics have a similar G0 as well as the organo sheets. The surfaces of the fabrics are comprised by the fibers and the empty space in between and so the surface adhesion is determined by this meso-structure; on the other hand, the organo sheets' surfaces are dominated by the thermoplastic PA6.

Moreover, the time to separate the gripper stamp from the fixed substrates is evaluated (see Figure 8) in order to determine the actual storage modulus according to Equation (3). Figure 8 reveals that the separation time is with a good approximation independent of the substrate in contact with the gripper stamp. The red line shows the exponential decay trend over the pulling speed.

**Figure 8.** Pulling speed related time to separate the gripper stamp from the substrates; red line shows the trend of all measurement points.

From the data presented above, the critical energy release rates G are computed according to Equation (6) and are shown to be pulling speed-dependent in Figure 9. Also the modeled energy release rates are shown with good agreement with the experimental data. The model seems to be appropriate to model and predict the critical energy release rate for in-service use. From this, the pulling speed needed to lift up/release an object with a certain mass could be calculated as well.

**Figure 9.** Critical energy release rate G versus the pulling speed v including the model of Equation (8) for each substrate.

#### **5. Conclusions**

The presented methodology to evaluate the viscoelastic as well as surface tack characteristics revealed that the gripper made of the soft PDMS formulation is capable for grasping, holding, transporting, and finally assembling (draping) of fabrics as well as semi-finished thermoplastic parts

(organo sheets with PA6 matrix). Based on the observations and the experimental data, the essential viscoelastic as well as adherence kinetic behaviors are determined. A high loading speed dependency of the tanδ characteristic is needed (i.e., with higher loading speed, higher tanδ) to exploit the rate-dependent adherence. The tanδ of PDMS-soft varies between 0.1 (@0.1 Hz) and 0.55 (@46.5 Hz), while the PDMS-hard formulation reveals a narrower tanδ range between 0.09 and 0.19. If the material has almost no rate-dependent tanδ, which is the case for PDMS-hard formulation (see Figure 5), then the material is inapplicable as a gripper. The storage modulus and its loading speed dependency are of minor importance (PDMS-soft: 0.25 MPa (@0.1 Hz) to 0.5MPa (@46.5 Hz); PDMS-hard: 2.4 MPa (@0.1 Hz) to 4.5 MPa (@46.5 Hz)). It should be low and remain low with increasing loading speed.

Grasping and assembling are achieved by the utilization of the rate-dependent adhesion. For the application, the knowledge of the gripper's pulling speed dependent adherence characteristic is essential. Here, a model (see Equation (8)) is derived based on the experimental results and describes sufficiently the rate dependency of the gripper. It can be applied for the prediction of the critical energy release rate and, more importantly, the maximum mass (weight force) to lift and hold. We believe that the proposed model is of particular interest for the application and operation of a robot handling system with adhesive grippers. In the presented work, the aging, and therefore the altering of the gripper material's mechanical inherent viscoelastic properties, are omitted. For the prediction of the service (replacement) interval times, the aging behavior has to be known and has to be verified in future.

**Author Contributions:** Conceptualization, U.D.Ç.; Investigation, U.D.Ç. and M.F.; Methodology, U.D.Ç. and M.F.; Project administration, U.D.Ç.; Software, U.D.Ç.; Supervision, Z.M.; Validation, U.D.Ç. I.G. and Z.M.; Writing–original draft, U.D.Ç.; Writing–review & editing, U.D.Ç., M.F., I.G. and Z.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Austrian Research Promotion Agency (FFG) in the framework of the Competence Headquarter Program operated by Engel Austria GmbH and by the LIT "ADAPT" grant number LIT2016-2-SEE-008.

**Acknowledgments:** Open Access Funding by the University of Linz.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Elastic Properties of Polychloroprene Rubbers in Tension and Compression during Ageing**

#### **Rami Bouaziz 1,2, Laurianne Tru**ff**ault 2, Rouslan Borisov 3, Cristian Ovalle 1, Lucien Laiarinandrasana 1, Guillaume Miquelard-Garnier 2,\* and Bruno Fayolle 2,\***


Received: 18 September 2020; Accepted: 11 October 2020; Published: 14 October 2020

**Abstract:** Being able to predict the lifetime of elastomers is fundamental for many industrial applications. The evolution of both tensile and compression behavior of unfilled and filled neoprene rubbers was studied over time for different ageing conditions (70 ◦C, 80 ◦C and 90 ◦C). While Young's modulus increased with ageing, the bulk modulus remained almost constant, leading to a slight decrease in the Poisson's ratio with ageing, especially for the filled rubbers. This evolution of Poisson's ratio with ageing is often neglected in the literature where a constant value of 0.5 is almost always assumed. Moreover, the elongation at break decreased, all these phenomena having a similar activation energy (~80 kJ/mol) assuming an Arrhenius or pseudo-Arrhenius behavior. Using simple scaling arguments from rubber elasticity theory, it is possible to relate quantitatively Young's modulus and elongation at break for all ageing conditions, while an empirical relation can correlate Young's modulus and hardness shore A. This suggests the crosslink density evolution during ageing is the main factor that drives the mechanical properties. It is then possible to predict the lifetime of elastomers usually based on an elongation at break criterion with a simple hardness shore measurement.

**Keywords:** polychloroprene; ageing; mechanical properties; rubber elasticity theory; crosslinking

#### **1. Introduction**

Rubbers (also called elastomers) are loosely crosslinked networks of amorphous polymers having a very low glass transition temperature [1]. Due to their macromolecular structure, these materials are relatively soft (Young's modulus on the order of 1–10 MPa) but display a large deformability (up to about 10 times their original length), mostly reversible, due to an elasticity that is entropic in nature [2]. Amongst this class of materials, polydiene-based elastomers such as polyisoprene or polychloroprene (neoprene) have remarkable mechanical properties (especially a large elongation at break) along with a good stability towards chemicals [3]. However, their use is often limited by the fact that they harden and become brittle over the long term. These changes in mechanical properties are most often caused by an oxidation phenomenon [4]. Indeed, the oxidation process leads to the production of free alkyl/peroxyde radicals along the chain which can then react with the double bonds by addition [5]. These inter-macromolecular addition reactions then contribute to an increase in the crosslinking density [6].

The classical way to predict the lifetime of elastomers is to consider that the time to reach an end-life criterion in terms of elongation at break (such as 50% of initial elongation at break) follows a pseudo-Arrhenius behavior. Based on accelerated ageing time test, an extrapolation at room temperature can then be done to predict lifetime. However, this kind of extrapolation is highly questionable since there are many experimental results showing that end life does not always follow such an empirical law [7].

To avoid such extrapolation, a possible way is to develop an approach in two steps. The first step is to build a kinetic model taking into account all the chemical mechanisms involved during ageing: according to this model, it is then possible to simulate the chemical state and the chain scission/crosslinking events whatever the exposure conditions in terms of time or temperature of exposure (see for instance for the polychloroprene [8]). The second step is to relate the chemical state and the chain scission/crosslinking events to the relevant property such as elongation at break or hardness. This kind of structure–property relationship has to be valid whatever the accelerated ageing conditions to be used in real conditions.

The prediction of the lifetime of elastomer parts first requires reliable constitutive models to be established in 3D, i.e., taking the hydrostatic pressure into account. The basic Mooney–Rivlin deformation energy *W* commonly used to assess the lifetime of engineering elastomeric structures is expressed as Equation (1) [9]:

$$\mathcal{W} = \mathcal{C}\_{10}(I\_1 - 3) + \mathcal{C}\_{01}(I\_2 - 3) + 0.5\mathcal{K}(I - 1)^2\tag{1}$$

where *C*10, *C*<sup>01</sup> and *K* are the material parameters, *K* being in particular the bulk modulus of the elastomer; *I*<sup>1</sup> and *I*<sup>2</sup> the first and second invariants of the deformation gradient tensor; *J* the Jacobian which can be written as the volume ratio *V*/*V*0, with *J* = 1 for perfectly incompressible materials.

Deriving the uniaxial Cauchy stress σ from Equation (1) leads to:

$$
\sigma = \frac{1}{J} \Big( 2C\_{10} + \frac{2C\_{01}}{\lambda} \Big) \Big( \lambda^2 - \frac{1}{\lambda} \Big) + K(J-1) \tag{2}
$$

with λ the elongation, defined as λ = l/l0 where l0 is the initial sample length.

The relationship between hyperelastic (Mooney–Rivlin) material parameters (*C*10, *C*01), bulk modulus (*K*) and elastic properties at small strain (Young's modulus *E* and Poisson's ratio ν) can then be summarized as follows:

$$E = \Im(2\mathcal{C}\_{10} + 2\mathcal{C}\_{01})\tag{3}$$

$$K = \frac{E}{3(1 - 2v)}\tag{4}$$

In this context, relating the evolution of both hyperelastic and elastic material properties to the progress of the crosslinking process is one of the objectives of this study, since it is required for lifetime prediction of the elastomers.

In the case of unfilled elastomers, the value of the Poisson's ratio is often considered to be equal to 0.5, corresponding to an incompressible state of the material, i.e., a deformation that occurs at constant volume [10,11]. However, elastomers are actually "ever so slightly" compressible and it is noteworthy that measuring the Poisson's ratio is not straightforward. In the case of filled elastomers, the value of the Poisson's ratio has been typically measured between 0.46 and 0.49 [12–14]. By using Equation (4) and considering a typical elastomeric Young's modulus of 10 MPa, such variation of the Poisson's ratio leads to an increase of the bulk modulus *K* from 42 MPa to 167 MPa. These values are in contradiction with the supposed incompressibility assumption letting the bulk modulus *K* tends to infinity, and a small uncertainty in the Poisson's ratio measured value can easily lead to an uncertainty of as much as an order of magnitude for *K*.

To the best of our knowledge, there are no results available allowing us to characterize such a pair of hyperelastic and elastic parameters during an oxidative crosslinking, which is then required to perform finite elements simulations of industrial parts. In the literature, the changes in mechanical behavior during the crosslinking process induced by oxidation are characterized mainly by tensile tests to assess *E* [15] or toughness in mode I [16]. However, although different loading conditions are investigated through these experiments, none of them gives access simultaneously to two parameters as described above. Therefore, a relevant mechanical characterization is required for a proper assessment of *E* and ν (or *K*), which can monitor not only the changes in terms of Young's modulus during the crosslinking process due to ageing, but also the possible changes in terms of compressibility modulus or Poisson's ratio.

In this study, the objective is to characterize the mechanical behavior of filled and unfilled crosslinked polychloroprene through tensile and oedometric compression tests to assess both the change of *E* and ν (or *K*) over time at different oxidative conditions. The choice of several oxidative conditions has been driven to check possible relationships between these values independently of oxidation rate, i.e., temperature of exposure. A specific attention is paid to promote spatial homogeneous oxidation for the samples considered here in order to assess only intrinsic material relationships.

Finally, the evolution of the elastic properties *E* and *K* during ageing is compared to elongation at break or hardness which are classically used to follow oxidative ageing of elastomers. The second objective of the study is then to propose correlations between all these mechanical properties, which shall be valid not only during the whole oxidation process but also for both unfilled and filled elastomers. As a result, a critical value for a specific property (easy to measure, such as hardness shore A) could be used as a proxy for critical values of properties that are much harder to assess (such as elongation at break), for lifetime prediction purposes.

#### **2. Materials and Methods**

Thin rectangular plates (thickness = 2 mm) of unfilled neoprene rubber and of neoprene filled with 40 phr (~30 wt %) of carbon black were provided by Nuvia (Villeurbanne, France). To study the impact of the thermal oxidative ageing on the mechanical properties, the plates were aged at three different temperatures (70 ◦C, 80 ◦C and 90 ◦C) during specific ageing times. It was verified with microindentation tests that ageing was homogeneously distributed throughout the samples (see Appendix A). All samples were thoroughly characterized after each ageing time. Specifically, differential scanning calorimetry (DSC) was used to follow the consumption of antioxidants/stabilizers under oxidative ageing. Mechanical tests, namely tensile, hardness, and oedometric tests, were developed to follow the evolution of the mechanical properties, namely Young's modulus, elongation at break and bulk modulus, with the ageing time. All mechanical tests were carried out at ambient temperature.

#### *2.1. DSC*/*OIT*

The antioxidants/stabilizers consumption was analyzed through the determination of the oxidation induction time (OIT) measured with a Q20 DSC from TA Instruments (New Castle, DE, USA). Polychloroprene samples with a mass ranging between 7.0 and 8.0 mg were placed in aluminum standard pans (TA Instruments) without a lid. Reference consisted of an empty aluminum pan, also without a lid. After equilibrating the temperature at 50 ◦C and an isothermal of five minutes, the samples were heated to 180 ◦C at 10 ◦C/min under nitrogen flow at 50 mL/min. Temperature was then maintained at 180 ◦C for five minutes before substituting nitrogen for oxygen for 250 min. OIT corresponds to the intersection of the tangent of the baseline at the time of atmosphere change with the tangent at the beginning of the oxidation peak.

#### *2.2. Tensile Tests*

The tensile tests were conducted on H3 dog-bone specimens (4 mm wide, 2 mm thick and 20 mm long obtained with a cutting-die) with a 5966 Instron tensile machine (Instron, Élancourt, France). Two samples were tested per ageing condition. Forces were monitored using a load cell of 10 kN. The specimens were stretched up to failure with a crosshead speed set to 20 mm/min. Due to the ability of the samples to undergo very large strains, pneumatic clamping jaws were used to avoid slippage.

Young's modulus was determined by fitting the stress–elongation curves in the λ = [1; 1.45] region with a Mooney–Rivlin incompressible model in uniaxial tension, i.e., following Equations (2) (adapted for engineering stress) and (3).

#### *2.3. Hardness Tests (Shore A)*

The hardness of the material was measured using a shore A durometer (model LX-A-Y from RS PRO, Corby, UK). The measured values indicate the resistance to indentation of the tested material on a scale between 0 (depth of indentation equal to the sample thickness) and 100 (no indentation).

#### *2.4. Oedometric Compression Test*

The oedometric compression tests were conducted on discs (diameter = 25 mm and thickness = 2 mm) obtained from the rubber plates with a cutting die. A stack of three discs was inserted in a lab-made apparatus consisting in a steel die with a compression applied by means of a rigid stamp. A small gap (less than 0.5 mm) between the samples and the steel die was initially set to allow the specimens to be positioned inside the die. The compression load was applied with a tensile/compression Instron machine (model 5982, Instron France) with a load cell of 100 kN. The specimens were compressed to about 95 kN with a displacement speed set to 1 mm/min.

#### **3. Results**

#### *3.1. Chemical Stability: OIT Changes*

Figure 1a reports the changes of OIT as a function of ageing time for filled and unfilled polymers at the three temperatures of exposure. The initial OIT gap between unfilled and filled samples can be attribute d to different processing conditions leading to a higher antioxidants loss in the case of the filled samples, such as a supplementary blending step at high temperature, using an extruder or an internal mixer.

**Figure 1.** (**a**) Oxidation induction time (OIT) and (**b**) relative OIT (OITr) evolution as a function of ageing time.

As expected, OIT decreases with ageing time, indicating a consumption of the antioxidants/stabilizers during the samples exposure. OIT reaches a plateau value close to 0 for the filled rubbers indicating an almost total consumption of the antioxidants. This plateau is reached at 2000–3000 h for samples aged at 90 ◦C and after about 6000 h for the filled samples aged at 70 ◦C, while it was not obtained for the unfilled samples aged at the same temperature during the total time of the experiment (8000 h).

In order to take into account the fact that unfilled sample shows a higher initial OIT value, we plot here in Figure 1b "relative OIT" (OITr) which is defined equal to the actual OIT for the filled samples, and for the unfilled values as OIT − 0.6 (which is the OIT difference at t = 0 h between unfilled (1.7) and

filled (1.1) samples). Thus, Figure 1b shows the relative OIT for the three temperatures of exposure and allows to compare both elastomers in terms of OIT decrease. It clearly appears that the relative OIT decrease rate (i.e., the slope of the "linear" region shown as dotted lines in Figure 1b) follows the same trend for both elastomers whatever the exposure temperature. As a result, filled and unfilled elastomers show the same ageing kinetic and differ only by their initial state.

A classical way to quantify the relative OIT change during ageing as a function of temperature is to consider the slope *x* of the linear region follows a pseudo-Arrhenius behavior with temperature [7], i.e., *x* = *x*0*e Ea* <sup>R</sup>*<sup>T</sup>* with *Ea* the activation energy (J/mol), R the gas constant = 8.314 J/mol/K and *T* the temperature (in kelvin). *Ea* is then considered as a characteristic of the process and its activation by temperature. Here, the temperature of exposure activates the OIT drop rate with an activation energy close to 78 kJ/mol considering such behavior.

#### *3.2. Tensile Tests*

Figure 2 displays representative tensile curves for unfilled (a) and filled (b) samples aged at 90 ◦C and tested at different ageing times. The plots show the engineering stress (the applied load divided by the unstretched cross-sectional area) as a function of the elongation λ.

In Figure 2a, it can be observed that the unfilled rubber undergoes a very large strain, leading to a pronounced hardening before failure. For polychloroprene rubbers, it has been reported that this hardening is due to strain-induced crystallization (SIC) [17]. As proposed in [18], SIC-related stress hardening during tensile test appears amplified by the presence of fillers (see Figure 2b). This, along with the filler properties, leads to the samples initially stiffer but with a smaller elongation at break than the unfilled rubbers. Both SIC onset associated to the stress–elongation curve upturn and elongation at break decrease monotonously with ageing. Furthermore, the upturn phenomenon is only clearly visible for the nonaged sample and disappears for the aged samples. However, all tensile curves show a similar trend with ageing: stiffening and elongation at break decrease, both trends being attributed to extra crosslinking induced by ageing.

**Figure 2.** (**a**) Tensile curves to failure of samples aged at 90 ◦C for unfilled and (**b**) filled rubbers.

Moduli changes can be extracted from the stress–elongation curves for both rubbers at all exposure conditions. Figure 3 shows the evolution of the Young's modulus (*E*) with ageing time for unfilled (a) and filled (b) rubbers. This confirms that *E* continuously increases with ageing and that this increase is more pronounced at higher temperatures of exposure for a given ageing time. Similar trends have been obtained on a carbon black filled butadiene in [19]. Since modulus is driven by the crosslink density, it is thus possible to follow chemical modification through this quantity. Activation energy for the modulus (calculated using the same procedure as the one described for the OIT assuming

here an Arrhenius behavior) increase rate is close to 90 kJ/mol for the unfilled and 84 kJ/mol for the filled rubber.

**Figure 3.** (**a**) Evolution of Young modulus (*E)* with ageing time for unfilled and (**b**) filled rubbers at 70 ◦C, 80 ◦C and 90 ◦C.

In order to study the fillers' effect on the modulus evolution, the ratio of the filled and unfilled moduli is plotted as a function of ageing time at the three temperatures in Figure 4. It appears this ratio is almost constant for all temperatures and ageing time, with a value of about 4.5. It is worth noting that this value is higher than the Guth and Gold equation [20,21] given by Equation (5):

$$\frac{E\_{\text{fulled}}}{E\_{\text{unified}}} = \, 1 + 2.5\varphi + 14.1\varphi^2 \,\text{.}\tag{5}$$

where ϕ is the effective volume fraction of filler.

Assuming a density of 2 g/cm<sup>3</sup> for carbon black and 1.2 g/cm3 for the elastomer [22], the volume fraction of fillers here is about 0.2, which yields a value of 2 for the modulus ratio using Equation (5).

A modified equation by Guth [20] takes into account the shape factor *f* of the carbon black (see Equation (6)):

$$\frac{E\_{\text{filled}}}{E\_{\text{unified}}} = 1 + 0.67f\phi + 1.62f^2\phi^2\tag{6}$$

With a shape factor *f* = 6.5 as determined by Mullins and Tobin [23], Equation (6) gives a value of 4.6 for modulus ratio, consistent with the experimental data, as can be seen in Figure 4.

**Figure 4.** Ratio of the filled and unfilled moduli as a function of ageing time for three ageing temperatures.

Finally, if we consider modulus changes as a tracer of the chemical degradation, a correlation between chemical stability (i.e., OIT) and the evolution of normalized Young's modulus (defined as *E*(*t*)/*E*(*t* = 0)) shall be evidenced. Figure 5 presents normalized modulus as a function of normalized OIT (OITn = OIT(*t*)/OIT(*t* = 0)) for both filled and unfilled rubbers at the three ageing temperatures. All the data display a similar trend and fall within the same envelope (see dotted lines in Figure 5): starting from high normalized OIT values (unaged rubbers), normalized Young's modulus increases as normalized OIT decreases until a limit value of OIT which corresponds to a total consumption of stabilizers (OITn reaching values close 0).

**Figure 5.** Normalized modulus En as a function of normalized oxidation induction time OITn.

From Figure 2, it is also possible to extract the elongation at break evolution, which is a relevant parameter often used as a criterion to define end-life duration in industrial applications of such materials [24]. As expected, Figure 6 shows the elongation at break decreases with ageing, again in a much more pronounced way at higher temperatures.

This embrittlement process is necessarily associated with the crosslink density increase already witnessed through modulus changes. From a kinetic point of view, elongation at break decrease rate shows an activation energy close to 86 kJ/mol for the unfilled and 74 kJ/mol for the filled rubber, consistent with the values obtained for the moduli.

**Figure 6.** Evolution of the elongation at break with ageing time.

#### *3.3. Hardness Tests (Shore A)*

Although values obtained from hardness tests are a priori too qualitative to identify the elastic parameters, they are easy to obtain and as such often used in the industry. Figure 7 shows the evolution of the materials hardness with ageing time. It appears that the hardness of the rubbers increases with ageing similarly to the Young's modulus.

**Figure 7.** Hardness Shore A evolution with ageing time.

#### *3.4. Bulk Modulus*

The small-strain elastic mechanical behavior of isotropic materials is fully characterized by a couple of parameters, for example Young's modulus (*E*) and bulk modulus (*K*). If the evolution of the Young's modulus with ageing is often reported in the literature [15], no data concerning such evolution of the bulk modulus of elastomers are available according to our knowledge.

Figure 8 shows such evolution for unfilled (see Figure 8a) and filled (b) rubbers. The values obtained are in the GPa range, similar to previous reports for elastomers [25]. Contrary to what has been observed for Young's modulus (see Figure 3), the effect of ageing on the bulk modulus is small both in terms of time and temperature dependency (i.e., less than 5% variations). It can also be noticed that the bulk moduli of filled rubbers are slightly higher than the ones of unfilled rubbers (from about 1100 MPa to 1200 MPa).

Let us recall Young's modulus can be related to the crosslink density, whereas bulk modulus is essentially controlled by van der Waals forces [26,27]. Chemical degradation yielding a crosslinking process does not modify significantly the van der Waals' energy.

**Figure 8.** (**a**) Evolution of the bulk modulus (*K*) with ageing time for unfilled and (**b**) filled rubbers at 70 ◦C, 80 ◦C and 90 ◦C.

#### **4. Discussion**

#### *4.1. Structure-Properties Relationships*

The Young's modulus of rubbers can be related to the crosslink density following Equation (7):

$$E = \text{\textquotedblleft N}\rho \text{RT.} \tag{7}$$

with N the crosslink density in mol/kg and ρ the elastomer density (kg/m3) [28]. As a consequence, the crosslink density modifications which can be deduced from *E* measurements may be linked to OIT, elongation at break and hardness changes. Although FTIR measurements only slightly evidence oxidation products in the timeframe of the ageing experiments (see infrared spectra in Appendix B), it is well known that oxidation leads to crosslinking in the case of polychloroprene rubbers during long-term exposure in air, due to addition reactions between alkyl/peroxyl radicals and double bonds [6,8]. This is consistent with the fact that OIT decreases when *E* increases.

To confirm this, the activation energies estimated in the previous section for all these properties are listed and can be compared in Table 1: it appears that all activation energies are close to 80 kJ/mol, suggesting that all these properties changes have the same chemical origins and are driven by the same mechanisms. In [4], similar activation energy characterizing the degradation was found (~90 kJ/mol), the values for the unfilled neoprene also being slightly higher than those of the filled material. Finally, note that the hardness shore A activation energy value is a bit smaller for the unfilled samples but that is probably due to experimental uncertainty as the variations of hardness over ageing time are very small for these samples.

**Table 1.** Activation energy for all parameters understudied.


Let us also examine direct correlation between *E* and hardness. Many relationships between *E* and hardness shore A (*SA*) are proposed in the literature [29,30]. In [30], a finite element simulation study leads to a linear relation between the logarithm of the elastic modulus and the hardness values (see Equation (8)):

$$\log\_{10} \mathbf{E} = \mathbf{c}\_1 \; \mathbf{S}\_A - \mathbf{c}\_2 \tag{8}$$

with *c*<sup>1</sup> and *c*<sup>2</sup> empirical constants determined from the best fit as *c*<sup>1</sup> = 0.0235 and *c*<sup>2</sup> = 0.6403 for 20 < *SA* < 80 [30].

In Figure 9, *E* is plotted as a function of *SA*. According to the previous approach, our experimental results lead to *c*<sup>1</sup> = 0.0245 and *c*<sup>2</sup> = 1.075 for both filled and unfilled samples, independently of ageing, in the 40 < *SA* < 85 range, which is in good agreement with the values reported previously [30].

**Figure 9.** Young's modulus (E) versus hardness shore A for all samples (filled, unfilled) and temperatures.

Finally, we would like now to study possible correlations between λ*<sup>b</sup>* and *E*. As a first approximation, the elongation at break can be related to the chain dimensions [1]: assuming an initial Gaussian conformation *<sup>R</sup>***<sup>0</sup>** <sup>=</sup> <sup>√</sup> *nl* with *R0* the chain end-to-end distance, *<sup>n</sup>* the average number of repeating units in the chain and *l* the monomer typical size, and fracture occurring when the chain extension reaches a maximal value *Rmax* <sup>=</sup> *<sup>n</sup>***l**, we can write Equation (9):

$$
\lambda\_b \propto \frac{R\_{\text{max}}}{R\_0} = \sqrt{n} \tag{9}
$$

Similarly, the Young's modulus of rubbers can be related to the molar mass between crosslinks, according to Equation (7) rewritten under the form *E* <sup>=</sup> **<sup>3</sup>** <sup>ρ</sup>**R***<sup>T</sup> Mc* with *Mc* the molar mass between crosslinks (in kg/mol), which is proportional to *<sup>n</sup>*. Hence, <sup>λ</sup>*b* should scale as **<sup>1</sup>**/ √ *E* which is indeed observed in Figure 10. Interestingly, this simple scaling remains valid whatever the ageing conditions and for both filled and unfilled materials.

**Figure 10.** <sup>λ</sup>*b* versus **<sup>1</sup>**/ √ *E* for all samples (filled, unfilled) and temperatures.

The linear fit gives the following value with E in MPa: <sup>λ</sup>*b* <sup>=</sup> **20.7**/ √ *E*−**4**.

Assuming again a density of 1.2 g/cm3, and with a molar mass *M0* of 88 g/mol for the polychloroprene monomer, the calculated prefactor is <sup>3</sup>ρR*<sup>T</sup> <sup>M</sup>*<sup>0</sup> <sup>≈</sup> 10 MPa−1/2. This is in very good agreement with the value obtained from the fit. Note the calculated value is obtained by assuming the typical molar mass between crosslinks is the molar mass of the whole chain, which explains the slightly underestimated value obtained.

To conclude this part, correlations between *E*, λ*<sup>b</sup>* and hardness shore A show a direct relation according to simple physical relations from rubber elasticity theory. Interestingly, these correlations are valid for the filled and unfilled rubbers during ageing. From a practical purpose, the choice of a critical value for one of these three properties can be simply translated in terms of crosslink density changes during ageing.

#### *4.2. Poisson's Ratio*

All mechanical properties determined here have to allow the description of the behavior law which can be, for example, later used in a finite element code: *E* and *v* or the shear modulus *G* and *K* for instance. The final objective would be to perform simulations to predict strain and stress in an elastomer part under loading and after ageing.

Figure 11 shows the evolution of the Poisson's ratio (*v*) with ageing time of both unfilled (a) and filled (b) rubbers aged at three different temperatures. Values of *v* were obtained using the following relationship that links *<sup>E</sup>* and *<sup>K</sup>* (see Equation (4), rewritten as follows): *v* <sup>=</sup> **0.5 (1**−*E*/**3***K*).

**Figure 11.** (**a**) Poisson's ratio versus ageing time for unfilled and (**b**) filled rubbers at 70 ◦C, 80 ◦C and 90 ◦C.

It can be noted again that the often used approximation for rubbers *v* = 0.5 is an upper limit (see Figure 11) obtained when *K* tends to infinity. For this reason, this limit value cannot be used in finite elements simulations.

We can observe a small linear decrease of *v* from 0.4996 to 0.4986 at 90 ◦C and 0.4988 at 80 ◦C in the case of filled rubbers. However, in the case of unfilled rubbers (and of filled rubbers at 70 ◦C), the decrease of *v* may be considered as negligible, with a value of 0.4995 ± 0.0004. To the best of our knowledge, this is one of the first measurements of the evolution of the Poisson's ratio with ageing for elastomers.

#### **5. Conclusions**

In this article, we described the evolution of Young's modulus, elongation at break and hardness shore A as a function of ageing conditions. In this study, we also developed an experimental protocol allowing to estimate the Poisson's ratio evolution for the same sample. If the Poisson's ratio is only marginally affected by ageing for the unfilled rubbers, a slight decrease is witnessed in the case of the filled rubbers. This point, previously unseen in the literature, has to be taken into account for 3D element finite simulations of aged rubber industrial parts.

It is shown that Young's modulus and elongation at break can be quantitatively linked to each other using simple scaling arguments coming from the rubber elasticity theory, whatever the samples (filled or unfilled) and ageing time and temperature. This strongly suggests that the mechanical properties evolution is mainly driven by the extra crosslinking of the neoprene matrix, induced by ageing.

Since Young's modulus and hardness are also related quantitatively to each other, it can then be possible to predict the lifetime of elastomers based on an elongation at break criterion with a simple hardness shore measurement. This result may have numerous applications for industries where predicting such lifetime is of fundamental importance. This precise measurement of the Poisson's ratio could also open the way to more quantitative mechanical models.

**Author Contributions:** Conceptualization, L.L., G.M.-G. and B.F.; methodology, R.B. (Rami Bouaziz), L.T., G.M.-G.; C.O., R.B. (Rouslan Borisov), L.L., B.F. validation, R.B. (Rami Bouaziz), L.T., G.M.-G., L.L. and B.F.; formal analysis, R.B. (Rami Bouaziz), L.T., G.M.-G. and B.F.; investigation, R.B. (Rami Bouaziz), L.T. and C.O.; resources, L.L. and B.F.; data curation, R.B. (Rami Bouaziz), L.T. and R.B. (Rouslan Borisov); writing—original draft preparation, R.B. (Rami Bouaziz), L.T., G.M.-G. and B.F.; writing—review and editing, R.B. (Rami Bouaziz), C.O., L.L., R.B. (Rouslan Borisov) G.M.-G. and B.F.; visualization, R.B. (Rami Bouaziz) and L.T.; supervision, L.L., C.O., G.M.-G. and B.F.; project administration, R.B. (Rouslan Borisov); funding acquisition, L.L., G.M.-G., B.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by EDF and Nuvia provided the rubber samples used in the study.

**Acknowledgments:** The authors would like to thank M. Kuntz (EDF), B. Basile (Freyssinet) and C. Guyonnet (Nuvia) for fruitful discussions. L. Truffault and R. Bouaziz also thank S. Baiz and C. Gorny from the PIMM lab for their help with the samples preparation and microindentation measurements.

**Conflicts of Interest:** The authors declare no conflict of interest. EDF participated in the design of the study and in the interpretation of the obtained data. EDF agreed to the publication of the results presented in the study.

#### **Appendix A**

Small pieces of neoprenes (about 10 mm long, 5 mm large and 2 mm thick) were cut and embedded in epoxy resin. The epoxy embedded neoprenes were then polished in order to obtain flat areas. Local moduli were measured each 200 μm along the neoprenes thickness with a microhardness tester from CSM + Instruments equipped with a conospherical indenter. Hertz's model was used to determinate the modulus from the loading region of the strength–indentation depth curves.

The modulus profiles along the thickness measured by microindentation for three filled neoprenes aged at 90 ◦C during different times are presented in Figure A1. For the most aged sample (3000 h), the highest modulus increase between the center of the sample and its edge was around 17%. This value was slightly higher (33%) for the unaged sample. Our results demonstrate the absence of oxidation profile within the samples thickness. Indeed, these increases are both very small compared to the increases reported by Wise et al. [31] for neoprene aged at 100 ◦C (>100% increase) and by Le Gac et al. who reported a modulus approaching 1 GPa in the 100 μm superficial zone for neoprene aged at 120 ◦C [32].

Furthermore, it is interesting to note that the ratio between the average moduli measured by microindentation and the average moduli measured by traction is similar for the different exposition conditions (~1.5–2).

**Figure A1.** Evolution/profile of modulus through the thickness measured by microindentation of a 2 mm-thick sample aged at 90 ◦C for different durations.

#### **Appendix B**

Unfilled neoprene plates were analyzed by Fourier transform infrared spectroscopy (FTIR) in attenuated total reflectance (ATR) mode using a Frontier FTIR spectrometer from Perkin Elmer. Spectra were obtained in the 4000–650 cm−<sup>1</sup> region with a resolution of 4 cm−<sup>1</sup> and 16 scans.

The spectra of neoprene aged at 90 ◦C for different times are presented in Figure A2. The spectral range between 2000 and 1000 cm−<sup>1</sup> is of particular interest when studying polychloroprene degradation due to the presence of bands attributed to oxidation products [4,5,32]. It should be noted that unfilled neoprene has been aged up to 6000 h for this analysis, while in the rest of the study the longest ageing time was 3000 h. The intensity of the band centered at 1590 cm−<sup>1</sup> clearly increases with ageing, but the change becomes only significant after 4900 h. This band can be attributed to carboxylate species. Their presence results from the conversion of acid chlorides and carboxylic acids during ageing in the presence of a small amount of metal oxides such as ZnO, classically used as an activator for the sulfur vulcanization of rubbers [33].

**Figure A2.** FTIR spectra of neoprene aged at 90 ◦C for different durations.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
