The Influence of Carbon Black Colloidal Properties on the Parameters of the Kraus Model

The Payne Effect (also known as the Fletcher–Gent Effect) has a fundamental impact on the behavior of filled rubber composites and therefore must be considered during their design. This study investigates the influence of carbon black (CB) surface area and structure on the observed Payne Effect and builds on the existing models of Kraus and Ulmer to explain this phenomenon. Dynamic strain sweeps were carried out on natural rubber (NR) compounds containing eight different grades of CB at equivalent volume fractions. The loss and storage moduli were modeled according to the Kraus and Ulmer equations, using a curve optimization tool in SciPy. Subsequent regression analysis provided strong correlations between the fitting parameters and the CB structure and surface area. Using this regression analysis, this work provides further insight into the physical meaning behind the Kraus and Ulmer models, which are phenomenological in nature.


Introduction
When subjected to dynamic strains beyond 0.1%, filled elastomer composites exhibit a sigmoidal decrease in storage modulus (G ) and a corresponding peak in loss modulus (G ). This phenomenon is known as the Payne Effect or occasionally referred to as the Fletcher-Gent Effect [1][2][3][4]. The Payne Effect is mainly attributed to the strain-dependent breakdown and reformation of particle-particle contacts, held together by weak van der Waal's forces [4][5][6]. This breakdown results in irreversible energy dissipation, which has fundamental consequences on the performance of an elastomer composite. For certain applications, including tires, belts, and anti-vibration devices, good fatigue resistance is required, and therefore low levels of energy dissipation are preferred to avoid excess heat build-up. In such instances, a small Payne Effect is desirable. For fracture applications, however, energy dissipation allows for the removal of energy from a propagating crack, thus improving the tear resistance of the elastomer [7]. In these cases, a larger Payne Effect would allow for higher levels of energy dissipation. Most engineering applications require a trade-off between these phenomena; hence, there is motivation to model the Payne Effect.
Several models have been proposed to describe the Payne Effect, including the van de Waale, Tricot, and Gerspacher models, the network junction model, the chain-slippage model, and the links-nodes-blobs model [8][9][10][11][12]. Kraus proposed the first quantitative model based on the breakdown and reformation of assumed van der Waals forces between filler aggregates [13]. This model assumes that the amount of network breakdown per cycle, R b , is proportional to the number of existing CB contacts, N, as a function of strain amplitude, f b : where k b is a constant. The network reformation rate, R m is assumed to be proportional to the number of broken contacts as a function of strain amplitude, f m : (2) where N 0 is the number of intact contacts in the network and k m is a constant. At equilibrium, when R b = R m , N is given as: The excess in storage modulus over that at infinite strain, G (γ) − G ∞ , is assumed to be proportional to the number of existing particle contacts. Using power laws, f b and f m can be rewritten as: where γ is half the peak-to-peak strain amplitude and m is a constant. With this in mind, the Kraus model for G can be rewritten as: where G 0 is G at zero strain amplitude and γ c is at characteristic strain equal to (k m /k b ) 1/2m . The G 0 − G ∞ term, commonly denoted as ∆G , equates to the magnitude of the Payne Effect. According to Kraus, the loss modulus at a given strain, G (γ), in excess over the loss modulus at infinite strain, G ∞ , arises from the energy dissipation associated with the breakdown and reformation of individual contacts. In other terms: where C 1 is a constant. With f b = γ m , Equation (6) can be differentiated to obtain an expression for the maximum G , i.e., G m . This can be rearranged to derive the Kraus expression for loss modulus: where γ c is the strain at which G (γ) is a maximum. The Kraus model of the storage modulus, given by Equation (5), was found to be applicable at various loadings of CB under different loading conditions, including tension, torsion, and shear [14,15]. The Kraus model for the loss modulus, Equation (7), demonstrated lower success in a study carried out by Ulmer, which reported significant deviations in the low strain regions of the model [16]. Kraus predicted that G reduces to G ∞ at zero strain because the extent of the broken network available for reformation also reduces to zero, which is not the case experimentally. Consequently, Ulmer proposed an additional contribution which is assumed to be independent of the CB network, but rather the result of an additional loss contribution associated with the process of polymer network breakdown and reformation. This additional term is proportional to the number of particle contacts: Applying the Kraus strain dependencies of f b and f m , the C 1 k b N f b term becomes the Kraus expression for G (γ) − G ∞ obtained from Equation (7). The additional term, C 2 N, is evaluated by replacing N with Equation (3). The second term is a simple, empirical choice of function that decreases exponentially with increasing strain: where G m − G ∞ of the Kraus model is called ∆G k , and C 2 N 0 is denoted as ∆G 2 . Equation (9) gives a considerably better description of the loss modulus than the Kraus derived alternative (Equation (7)). There are some disagreements between the estimates of m and γ c , obtained from the Kraus G (γ) model (Equation (5)) and the modified G (γ) (Equation (9)) but this is usually less than ±5%. Ulmer demonstrated that Equations (5) and (9) can successfully model the Payne Effect for various types of rubber including NR, styrene butadiene rubber, nitrile butadiene rubber, and butyl rubber. Furthermore, the model accurately predicted the behavior to within ±2% error for compounds containing between 25 and 50 parts per hundred rubber (phr) of high abrasion furnace CB. The accuracy of these models slightly reduces beyond CB loadings of 80 phr, where the error reported increased to ±5%. The Kraus and modified Kraus equations are limited because they are phenomenological in nature and do not provide any information on how to predict the Payne Effect. Observations dating back to the 1960s have reported that the Payne Effect is heavily influenced by the nature of the reinforcing particles [17,18]. Despite these observations, the current models for predicting this behavior do not consider the basic properties of the filler. CB represents the most popular filler system for rubber-based applications, its reinforcing potential being largely governed by aggregate morphology and surface area. CB exists as aggregates, formed by the fusion of para-crystalline primary particles. The spatial arrangement of these primary particles is defined as the level of 'structure'. The structure of CB can be quantified using oil absorption measurements [19,20]. Primary particle size can range between 5 nm and 200 nm and is inversely proportional to the surface area. Surface area is typically measured using gas adsorption techniques [21]. CB manufacturers are able to precisely control the structure and surface area of CB, giving rise to a wide selection of commercially available CBs. In this paper, the term 'colloidal properties' refers to the structure and surface area of CB. These parameters, alongside surface activity, define the reinforcing potential of a CB filler [22]. Surface activity is related to the reactivity of the chemical groups present on the CB surface as well as the surface energy distribution [23]. It relates to how well a filler is able to interact with the polymer matrix. This particular feature will not be discussed in this paper because the selected CBs all display broadly similar surface chemistries. Many studies have demonstrated the effect of surface area on polymer-filler interface, inter-aggregate distances, and filler networking [24][25][26]. Meanwhile, the increased structure gives rise to occluded rubber which is screened from globally applied strains, effectively increasing the filler volume fraction [27][28][29]. More recently, researchers performed micro-structural studies on the role of structure and surface area of CB to explain observed phenomenology [30]. Particular studies demonstrated the effects of CB structure and surface area on the Payne Effect, concluding that structure and surface area dominate the low strain behavior, with the effects of surface area diminishing to zero as strain increases [31,32].
This work aims to develop an understanding of how the colloidal properties influence the terms in the existing Kraus and the modified Kraus models (Equations (5) and (9)). The goal of this research is to introduce physical meaning to these phenomenological equations and provide CB selection criteria based on their colloidal properties, to achieve the desired Payne Effect.

Materials
The data modeled in this study represents a small subset of data collected by Birla Carbon (Marietta, GA, USA). Kyei-Manu et al. proposed a broader study into the effects of CB colloidal properties on the mechanical and dynamic properties of the same materials [31]. Eight grades of CB were selected for this study, based on their relative positions on the colloidal plot, as shown in Figure 1, to cover a broad range of surface areas and structures. Figure 1 also shows some commercially available CBs (N772, N660, N330, N347, N220, and N115). These grades are included in the figure to provide additional context but are not evaluated in this work. The eight selected grades of CB were added to analogous NR-based compounds at equal filler loading of 50 phr. Details of the formulations can be found in Table 1. A naming convention was adopted in this paper to immediately identify the type of CB based on its structure and surface area; details of this are provided in Table 2. The structure (given as compressed oil absorption number, COAN) and surface area (given as statistical thickness surface area, STSA) are denoted as superscript and subscript, respectively. For example, N550 is referred to as CB 84 37 because it has a structure value of 84 mL·100g −1 and surface area value of 37 m 2 ·g −1 . The corresponding rubber compound produced using N550 is referred to using the same naming convention.
Compounding was carried out by Birla Carbon (Marietta, GA, USA) using a 1.6 L capacity Banbury Mixer. A three-stage mixing process was performed, to promote uniform dispersion of the CB. Compounds were vulcanized into sheets measuring 11 mm × 11 mm × ∼2 mm via compression molding. This was carried out at 150 • C for a time of T90 + 5 min, where T90 represents the time taken to reach 90% of the maximum torque on an Alpha Technologies moving die rheometer (MDR) located in Hudson, OH, USA. Interferometric microscopy (IFM) was used to analyze the compound dispersion indices (2-100 µm), following method D in the ASTM standard D2663 [33]. This testing was carried out at Birla Carbon (Marietta, GA, USA).

Dynamic Strain Sweep Characterization
Dynamic strain sweeps, between 0.1% and 39.5% single strain amplitude, were performed at 10 Hz with zero mean strain using an ARES G2 torsional rheometer from TA Instruments, located in New Castle, DE, USA. Testing was carried out at 60 • C. The specimen geometries were disks of 8 mm diameter and approximately 2 mm thickness. Specimens were adhered to the rheometer parallel plate geometry using Loctite 480 adhesive from Henkel, Hemel Hempstead, UK. A compressive pre-force of 100 g was applied to the disks for the duration of the tests. It is noted that with this particular specimen geometry, the strain amplitude varies across the specimen radius; the reported values are therefore the maximum values for the extremity of the disk radius. The samples were pre-conditioned six times at the specified dynamic strain amplitude before collecting the torque-time data. This process was repeated at each strain segment from low to high. Each compound was tested twice.

Curve Fitting and Optimization
The G (γ) data were fitted according to Kraus (Equation (7)) by using a curve optimization tool in SciPy which minimizes the sum of the squared residuals. A full breakdown of the code with added comments can be found in the supplementary information (S1 and S2). An initial optimization was run with a non-negative constraint, meaning that all parameters must take on values equal to or greater than zero. Using the same curve optimization tool, an additional non-negative constraint fitting of the G (γ) data was carried out as per Equation (9). Further constrained optimizations were performed, details of which are discussed in the subsequent Results and Discussion (Sections 3.4 and 3.5).

IFM Microscopy Results
The Payne Effect is sensitive to the state of filler dispersion, so this feature was important to consider during this work. Compounds with lower levels of filler dispersion are shown to exhibit increased values of G 0 and decreased values in the median strain required to cause a 50% reduction in ∆G [34]. The dispersion indices, as determined by IFM, are given in Table 2. Most of the compounds display good levels of dispersion with dispersion indices >98. CB 62 161 is noted to have a dispersion index of 81.5, which is significantly lower than the other compounds. At a fixed surface area, a decrease in structure leads to an increase in attractive contacts per unit volume of CB; the force necessary to separate aggregates in a pellet, therefore, increases [35]. As a result, grades of CB with a comparatively low structure-to-surface area ratio, for example, CB 62 161 , become increasingly challenging to disperse.

Strain Dependence of G and Non-Negative Curve Optimization
The raw data for the strain dependence of the storage modulus, G , is given in Figure 2. From the first observation, it is clear that the surface area and structure of CB greatly influence the behavior. An initial non-negative curve optimization was performed as per Equation (5), results from which are shown in Table 3. Throughout this paper, fitting parameters are presented as averages across the two repeated measurements; details about individual fitting parameters can be found in the supplementary information (S3-S9). Results across the two independent fittings are in good agreement, proving that there is good agreement between the sets of experimental data. The individual fittings all have R 2 values greater than 0.997, which indicates that they are all strong fits.

Strain Dependence of G and Non-Negative Curve Optimization
The raw data, showing the strain dependence of the loss modulus, G , are given in Figure 3. As with the G data, there are obvious differences between the various grades of CB upon initial inspection. The results from the initial unconstrained curve optimizations are given in Table 4. The non-negative constraint meant that values of G ∞ across all the fittings were 0, so this information is not included in the table. The R 2 values for these optimizations are high, with all individual fittings obtaining an R 2 value of at least 0.996. The repeated data sets are all in good agreement with one another. It is worth noting that values for m and γ c are different from those found in the previous optimization of G (γ), given in Table 3. By definition, γ c should be the strain at which the peak in the loss modulus is located, since this corresponds to the strain at which maximum cluster-cluster breakdown occurs. Unfortunately, the value of γ c cannot be taken directly from torsional experiments because this particular geometry shifts the G (γ) peak from its uniform strain position, therefore slight differences in this parameter are expected [16]. The m values are shifted to much higher values, the same observation was previously made by Ulmer.  An additional optimization was attempted whereby the G and G optimizations were run simultaneously using γ c and m as common fitting parameters. The R 2 coefficients for these curve optimizations are lower than the aforementioned fittings which are particularly reflected in the loss modulus; examples of these fittings are given in Figure 4. A previous evaluation of the Kraus model suggests that the universal value of m lies between 0.5 and 0.6, independent of specific filler type [36]. In Table 3, the average m value obtained is 0.44 ± 0.4. While these results may only appear to slightly deviate from one another, it was realized that the γ c and m terms effectively correct for one another, which could add further error to the fittings. For this reason, these initial optimizations were deemed to be inappropriate.

Selection of m
In theory, there should be agreement between the m and γ c values across the models for storage and loss modulus [16]. In practice, there are slight variations between the individual optimizations, as reported in Tables 3 and 4; this was also recognized by Ulmer when he proposed the modified Kraus model for the loss modulus. The m term exhibits a notable strain-dependency that has not been accounted for by either model, which adds a significant complication to its assignment. During torsional shear measurements, specimens do not undergo uniform strain, which adds further complexity to quantifying m. In the previous unconstrained optimizations, it was found that m and γ c are interdependent. This means that any uncertainty in assigning m will have direct consequences for γ c . Therefore, to better understand the effects of CB colloidal properties on γ c , a simple approach was taken to constrain m to a fixed value. The simple empirical choice of m = 0.55 was selected because this was the average value obtained across all the unconstrained fittings. This value coincides with the literature, which finds the universal range of m to be between 0.5 and 0.6, regardless of filler type [36].

Constrained Curve Optimizations of G (γ) and G (γ)
Results from the constrained curve optimizations of the storage and loss modulus are given in Tables 5 and 6. It is worth noting that only the γ c parameter was affected by the constraint. Applying the constrained value of m reduced the average error in the value of γ c across the optimizations for storage and loss modulus from 33.0% to 24.0%. The R 2 values for the individual optimizations remain strong for the G (γ) data but are slightly reduced for the G (γ) optimizations. In particular, the selected m value does not appear to be appropriate for CB 62 161 . The R 2 value for this CB has been significantly reduced by the m constraint. In the initial fitting given in Table 4, this species of CB has the highest m value of 0.82. This corresponds to a sharper peak in the loss modulus, which could relate to the lower dispersion index observed in this compound. Table 5. Fitting of G (γ) according to Equation (5) where the value of m has been assigned as 0.55.  Table 6. Fitting of G (γ) according to Equation (9) where the value of m has been assigned as 0.55. The γ 2 term appears to be independent of filler type, showing low variation across all compounds, as seen in Table 6. This parameter was devised by Ulmer as part of an additional term, independent of the filler network, which decreases exponentially with increasing strain. It is, therefore, assumed that this parameter is independent of CB grade and that further optimization of the other fitting parameters could be achieved by constraining this value. For these reasons, an additional curve optimization was performed, assigning γ 2 as 0.23, the average value obtained from Table 6. Results from this optimization are given in Table 7. Table 7. Fitting of G (γ) according to Equation (9) where γ 2 has been assigned as 0.23. This additional constraint has caused minor changes to the other fitting parameters but R 2 values are largely similar. The optimization of CB 62 161 results in a far lower R 2 value compared to the other compounds. This may be related to the lower dispersion index obtained for CB 62 161 , as previously determined in Table 1. The level of filler dispersion is known to affect various aspects of the Payne Effect [34]. Therefore, this particular grade of CB has been excluded from further regression analysis.

Regression Analysis
The Kraus and the modified Kraus equations use a phenomenological approach to successfully model the Payne Effect. The aim of this work is to bring physical meaning to these phenomenological equations. To achieve this, multiple linear regression analysis was performed on the fitting parameters obtained in Tables 5 and 6. This analysis was carried out in Origin 2019, which bases calculations on Equation (10).
where Y is the dependent property being examined (in this case the fitting parameters), C is the intercept, β St is the coefficient of structure, β SA is the coefficient of surface area and is the error. Full regression results are provided in Tables 8 and 9. Only those values with a corresponding p-value < 0.05 are statistically significant. Values that do not meet this requirement have been highlighted in red bold font. From the storage modulus regression results given in Table 8, a reversal in dominating colloidal property is observed as the strain amplitude is increased. At low strains, indicated at G 0 , both structure and surface area contribute towards the observed modulus value. The magnitude of β SA is larger than β St , which indicates that surface area has a larger effect on G 0 . At higher strains in the region of G ∞ , the surface area becomes statistically non-relevant and the structure dominates the behavior. At small strains, rigid filler networks augment the storage modulus, with the extent of networking being attributed to the number of aggregates per unit volume, which is related to their surface area. At high strains, these particle networks were broken down and the primary stiffening mechanisms are caused by strain amplification, a product of occluded rubber, which is determined by structural effects, among other factors. Despite subtle differences between γ c values determined from the storage and loss modulus optimizations, in both cases, the results are defined by surface area with structure playing a statistically insignificant role. As CB surface area increases, γ c values decrease. The negative correlation being observed relates back to the average number of aggregates per unit volume being higher for higher surface area CBs. A subsequent lowering in inter-aggregate distance requires smaller strains to reach critical breakdown. Results from the regression analysis of G are given in Table 9. The G m parameter refers to the height of the peak in G (γ) which is affected by the surface area and structure of CB, with the surface area having a much larger influence. G m defines the peak in the loss modulus, a product of energy dissipative mechanisms including filler-filler breakdown and polymer-filler slippage phenomena. Statistically, these events have a much higher probability of occurring as the surface area is increased. The ∆G 2 term is an additional polymer network term that controls the height of the small strain augmentation in loss modulus. Despite this contribution being polymer-related in definition, the surface area of CB appears to have a strong positive effect on the outcome. During their work, Ulmer investigated several different types of rubber, demonstrating that ∆G 2 is dependent on the polymer matrix. As a result, it was concluded that the small strain augmentation is a result of polymer networking. In this work, the polymer remains constant across all eight compounds. The strong correlation between CB surface area and ∆G 2 indicates that an additional contribution to the height of the small strain modulus may arise from intact filler networking. As previously discussed, the level of filler networking is dependent on CB surface area, which explains the positive correlation being observed in ∆G 2 .

Conclusions
During this work, the Kraus and the modified Kraus equations were able to successfully model the Payne Effect for a variety of NR-based compounds filled with various CBs which differed in their colloidal properties. An additional constrained optimization was carried out whereby the value of m was set to 0.55, to overcome issues faced with the interdependence of the m and γ c parameters. During this optimization, it was noted that CB dispersion affects the sharpness of the peak in the loss modulus and therefore the value of m. This meant that the selected m value was inappropriate for the CB 62 161 , which displayed the lowest compound dispersion as determined by IFM. The R 2 values obtained across all other optimizations were only slightly reduced following the addition of the m constraint. A final constraint on γ 2 was justified because this parameter showed low variation across all compounds. This constraint had a minimal effect on the R 2 values.
To bring physical meaning to these models, the roles of CB surface area and structure on the various fitting parameters were assessed using multiple regression analysis. As a result of this analysis, the following conclusions can be made: • At low strains, indicated at G 0 , there is a small structural contribution to the storage modulus arising from polymer occlusion, but the larger influence is that of the surface area which dictates the level of filler networking. • At higher strain regions of G ∞ , the surface area becomes statistically insignificant as the filler network is broken down, leaving behind only the effects of polymer occlusion, which is dictated by CB structure. • As surface area increases, γ c , the strain at which maximum cluster-cluster breakdown occurs decreases (when m is assumed to be constant). This relates to the reduction in the inter-aggregate distance as filler surface area is increased. • The peak in the loss modulus, G m , is affected by surface area and structure but the surface area has the dominating contribution. Energy dissipation arising from fillerfiller breakdown and polymer-filler slippage phenomena has a higher probability of occurring as the surface area is increased. • The additional polymer network term, ∆G 2 , correlates positively with filler surface area. This can be linked to higher levels of filler networking observed in high surface area CBs which causes augmentation to the small strain modulus.
As such, it is possible to tailor the magnitude of the Payne Effect. For most fatigue applications seeking low heat build-up and therefore a small Payne Effect, low surface area and highly structured CBs are recommended. For fracture applications that require high energy dissipation and a larger Payne Effect, high surface area and low structure CBs would be preferred. It is noted that these recommendations are only suitable for NR-based compounds with 50 phr CB, although similar observations have previously been made for NBR containing 60 phr CB [32]. Future work should be carried out to determine the effects of changing the filler volume fraction and polymer matrix to build upon the results found herein.