Study Mullins effect of polyurethane reinforcement with halloysite nanotube by molecular dynamics simulation

Abstract

Molecular dynamics simulation was applied to study the irreversible strain through loading and unloading cyclic tests of polyurethane (PU) reinforced with halloysite nanotube (HNT). The influences of the stretching cycle rate, different temperatures, the volume of halloysite nanotube and the density rate of the hard and soft domain of PU were studied on the permanent set. The results illustrate that the residual strain was increasing when the stretching loading is increasing, for example, with increasing strain load to 250% the residual strain increased to 55%. In contrast, the increasing volume fraction of HNT and hard part content of PU lead to lower residual strain. The recovery of the permanent set is achievable by increasing temperature from 1 K to 200 K residual strain is decreased to 52%. An Ogden constitutive and the theory of pseudo-elasticity were adopted to simulate this composite in the ABAQUS software. This model has proposed a reasonable prediction of plastic deformation.

Introduction

Elastomer materials show substantial stress softening when subjected to stretching. As presented in Figure 1, in a variety of materials such as polymer,^1,2 rubber^3–8 and living tissues,^9,10 the residual strain was recognized. Mullins and Tobin (1957) ¹¹ and the first time, Mullins (1969) ⁸ proposed expression for modelling the hyperplastic behaviour for this type of material. The majority of stress softening happens at the first cycle and sequenced cycles of stretching provoke the same or less stress softening. Thereby, the residual strain has remained after the first cycle. This residual strain is so-called permanent set. ¹²

Figure 1.

Schematic of the Mullins effect.

Most studies on the Mullins effect concentrate on filler;^3,11 however, stress softening is also observed only on matrix during stress loading (instead of strain). Therefore, both methods indicate a universal principle of these materials. Different types of physical interpretations for this phenomenon have been proposed to date. Generally, we can group them into three categories; ¹³ (1) interface between chains in the matrix, (2) interface between matrix and filler surface and (3) bond strength within the reinforcement component. ¹² Experimental observation of swelling analyses determined that the fracture in the matrix does not only depend on non-bond interaction between chains or chain-breaking during the axial test; therefore, bond breaking is not the only factor to explain the permanent set. ¹³ Technically, the controversial aspect of the study of large deformation of rubber material is scanning microstructure changes include the interface between polymer chains and polymer chain and filler during the mechanical stretching. However, Ma et al. ¹² studied the permanent set of a styrene-butadiene–rubber reinforced with silica particles; they employed a coarse-grained (CG), molecular model, without considering different segmental of the polymer composition.

Polyurethane/halloysite is composed of more than 95% (by weight) content of polyurethane polymer and less than 5% hard mineral ingredient in the form of halloysite structural. Natural and synthetic polymers reinforced with halloysite nanotubes (HNTs) play an important role in designing new composites. Carbon nanotubes are the most commonly studied nanomaterial; however, toxicity and high costs make them less attractive in industrial applications. The halloysite nanotube is extensively available around the world at low cost and so has received more attention from researchers. HNT has converted to an interesting reinforcement particle material for developing the functionality of the polymer. ¹⁴

Most research on halloysite has concentrated on determining the mechanical properties of HNT under different loading conditions. ¹⁵ The elastic modulus of the halloysite layer was ascertained by Lu et al. ¹⁶ in the range 130 ± 24 GPa. Lvov et al. ¹⁷ (2008) calculated Young’s modulus of a roll HNT; they obtained Young’s modulus in the range 230–340 GPa. The elastic modulus of different arrangements of HNT was calculated by Guimaraes et al. ¹⁸ The result has shown that the elastic modulus of the zigzag arrangement is in the range of 230–300 and 300–340 GPa for the armchair arrangement. ¹⁸ Halloysite forms tubes with a length of 600–1000 nm and an internal radius of 12–30 nm and an external radius of 30–100 nm .^18,19 Halloysite is a natural alumina-silicate (Al₂Si₂O₅(OH)₄.2H₂O) composed of a two-layer structure, where the outer face consists of siloxane structure, and the inner face is composed of alumina structure as described in Figure 2(a). ¹⁵

Figure 2.

(a) Structure of halloysite nanotube ¹⁵ and (b) schematic of Polyurethane. ²⁰

Polyurethane (PU) is a copolymer formed by combining a di-isocyanate (OCN–R–NCO) and a polyol (OH–R^′–OH), where R is aromatic or aliphatic and R^′ is a linear hydrocarbon chain. Also, PU is composited of soft segments (SS) which is from the polyol and hard segments (HS) create of the di-isocyanate of crystalline. ²⁰ PU has a variety of applications in a range of products. ²¹ Figure 2(b) shows the composition of a monomer polyurethane chain employs in this research. The chain is composed of five repeating monomers (n = 5) and 15 repeating polyol parts (m = 15) and one di-isocyanate part.

Molecular Mechanics is basically a method that uses potentials that are functions of interatomic distance. These potentials are determined from experiments or DFT calculations. The potential energy function is derived from the positions of the atoms that form a molecule. Some conventional potential energy functions used are Dreiding, UFF and PCFF.^22–25 Generally, the selection of an appropriate potential is attributed to the properties of the arrangement. The total energy U of a molecular structure is obtained with the following equation

$\begin{array}{l}U={\displaystyle \sum _{allbonds}\frac{1}{2}{K}_b{\left(b-b_0\right)}^2}+{\displaystyle \sum _{allangles}\frac{1}{2}{K}_{\theta }{\left(\theta -{\theta }_0\right)}^2}+{\displaystyle \sum _{alltorsions}{K}_{\phi }\left[1-\cos \left(n\phi \right)\right]}+{\displaystyle \sum _{\mathrm{\chi }}{K}_{\mathrm{\chi }}{\mathrm{\chi }}^2}+\\ {\displaystyle \sum _{i,jnonbonded}{\epsilon }_{ij}\left[2{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^9-3{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^6\right]}+{\displaystyle \sum _{i,jnonbonded}\frac{q_i{q}_j}{4\pi {\epsilon }_0\epsilon r_{ij}}}\end{array}$

(1)

The first four terms in equation (1) introduce the stretch bonds energy (b), bend angles energy (θ), torsion angles energy (φ) and distort energy with respect to a reference plane (χ). The final two functions define the nonbonded interactions as a sum of van der Waals and electrostatic interactions. All terms depend on the distance (r _ij ) between atoms. Each force field defines a different type of coefficient for a variety of chemical bonding. Because of computational accuracy, the polymer consistent forcefield (PCFF) with the interface extension (PCFF-IFF) has been used for this study.^25,26

In inhomogeneous material, different segments present different distributions of stiffness properties. Therefore, the role of each component on the behaviour has to be considered in the nanocomposite. Because of the complexity of elastomer systems which is combined with a filler and interfaces between components, studying the Mullins effect and understanding microstructure movement in the nano-scale within the stretching become more challenging. Although a few types of research have been studied at a molecular scale,^13,27 however, studying the molecular mechanism of the Mullins effect remains controversial. In this section, we develop an MD simulation of a reinforcement PU with HNT to study the atomic scale of the Mullins aspects such as permanent set, stress softening and recovery. Finally, an Ogden-type constitutive model was derived to model composite with FEM analysis.

Molecular dynamic model

The targeted system is a polyurethane (PU) reinforced with 1% wt halloysite nanotube (HNT). This percentage is according to the desirable range of ref. ²⁸ The models in this section were created by the amorphous polymer builder starting with 780 chains of PU and single HNT with 1803 atoms. ²⁹ The HNT is randomly embedded in the cell and the PU chain was randomly placed around it. The total number of atoms is 447,739. All the MD simulations were carried out with a ‘large scale atomic/molecular massively parallel simulator (LAMMPS)’ ³⁰ with the PCFF-IFF force field.^25,26

According to Dewapriya et al. ³¹ (2020), an annealed proceed was applied upon nanocomposite arrangement, then three steps equilibrium under an NPT ensemble at zero pressure and sequence of 300 K, 500 K and 300 K temperatures. During this process, we adopted the PPPM summation method ³² with the Lennard-Jones cutoff radius of 18 Å. Also, the periodic boundary conditions (PBC) were applied in states and the time step was 0.25 fs and the total time step is 1200 ps, the final density of the composite was 1.1071 g/cm³(Figure 3).

Figure 3.

Density of models after annealing cycle and relaxation process for PU (m = 15) with 1% wt HNT.

In MD simulation, the following approach is applied to study mechanical stretching tests. The PU and HNT systems were stretched by increasing the PBC length with constant speed (strain rate is 5*109 s-1) in the z-direction while the other box length directions were unchanged. When we reach a desirable strain, strain rate deformation is removed. During stretching, the NVT ensemble (with the Langevin thermostat) ³³ is used to keep the temperature at 1 K. During deformation, we adopted the PPPM summation method ³² with a Lennard-Jones cutoff radius of 18 Å. Also, the total time was 50 ps.

Results and discussion

In this section, MD simulations observation was used to study the permanent set, recovery and microstructure properties to ascertain the root of the Mullins effect.

Cyclic loading

Permanent set is defined according to incapacity to return to the original state. The first detection of the permanent set is observed when the PU/HNT composite was stretched in a cyclic to a maximum of 250% strains and returned to the zero-stress state. It can be noticed from Figure 4(a) that, after release load, the PU/HNT composite cannot back to the original state. Instead, a residual strain remains. The residual strain is 43% when the strain load is 250%.

Figure 4.

Stress–strain response of the PU/HNT composite (a) for three sequence cycles with strain amplitude 2.5, (b) for three sequence cycles with three different strain amplitudes from 1 to 1.5 to 2.5 is 2.5.

As displayed in Figure 4(a), the first stretch generates the majority of stress softening. The stress is reduced considerably between the first stretch cycle. Also, during the subsequent stretching cycles, the softening rate occurs. The area surrounded by cyclic curves was applied to determine the stress softening which identifies the dissipated energy density. ¹² The value in the first cycle is at 10.82 J cm³, and 8.96 J cm³ for the second and 4.036 J cm³ for the third cycle.

As presented in Figure 4(b), we have recognized residual strain was expanded with increasing strain load for example with increasing strain load to 250% the residual strain increase to 25%.

Figure 4(b) shows three loading-unloading applied strain cycles extending from 100% to 150%, to 250%. The stress softening was also observed when the strain loading is surpassed by the previous strain. Also, the stress softening decreases at the subsequent stretch cycles. The dissipated energy density is 1.51 J cm³, 2.92 J cm³ and 11.36 J cm³ for first, the second and third cycles, respectively. Figure 5 shows the simulation model. In all three models, deformation can be detected by two components, the HNT and the matrix. However, the deformation of the polymer chain matrix is more noticeable due to the PU having a much lower yield strength.

Figure 5.

Loading-unloading applied strain cycles, (a) PU/HNT, (b) strain = 1, (c) strain = 1.5, (d) strain = 2.5, (e) halloysite nanotube.

Temperature evaluates

Generally, the Mullins effect is dependent on the temperature and time, in which the temperature impact exponentially and time linearly controls the recovery.^34,35 Also, recovery of the residual set is achievable even at room temperature if relaxed for a long time. ³⁶ The permanent set of nanocomposites is recovered by increasing temperature. It has been proved by various experiments, ¹² as pictured in Figure 6, by increasing temperature from 1 K to 200 K residual strain is decreased from 0.37 to 0.23.

Figure 6.

Recovery of residual strain at different temperatures.

The original simulation is performed at T = 1 K, without external heat loading. As illustrated in Figure 6, the blue dashed curve refers to the higher temperature (T = 200 K) red curve refers to the stretching nanocomposite at a medium temperature (T = 100 K) stretching was under the NVT ensemble. As depicted in Figure 6, stress softening is recovered at a higher temperature system.

We also studied the weight of the hard and soft segment ratio PU on the mechanical response with changing repeating soft segment (SS) density through changing m to 20 and 10. The dispersion energy is declined with an increase in soft segment density as shown in Figure 7(b).

Figure 7.

The stress–strain response for (a) different HNT densities in the composite, (b) different the hard and soft segment ratio densities in HNT/PU.

Figure 7(b) displays the mechanical response of HNT/PU composite versus a fraction of hard segment content. The hard domains influence crosslinks and restrict the soft segment mobility. Furthermore, the soft segment controls chain mobility. Hence, increasing hard segment content increased elastic modulus and tensile strength due to reducing segment mobility. In contrast, increases in hard segment content reduce stress softening. However, increasing hard segment content reduces the soft segment mobility to diffusion which is unfavourable.

Additionally, we examine the influence of HNT percentage in nanocomposite on the starching cycles. The HNT volume fractions are increased from 1.0% to 3.0%. The residual strain is declined with an increase in HNT density as observed in Figure 7(a), which is compatible with experimental observations. ²⁸ Stress softening was observed in all three cases.

Generally, the larger the HNT density leads to a lower percentage of the fraction of monomers. This can be interpreted that the HNT provides strong constraints on the polymer chains. Increasing HNT density makes more restriction of some disentanglements because of the strong chemical bonding and difficulty to break them. Overall, increased reinforcement causes enhanced entanglements, as well as stronger localized stress concentration. However, the disentanglement is dependent on cross-link breakage among chains and between the chain and reinforcement filler. Thereby, there is an optimum value for HNT density and beyond entanglement would be decreased. Previous experimental research presented this percentage is around 0.8–1% wt. ²⁸

In order to evaluate the microstructure changes and consistent regeneration which are key characteristics of the Mullins effect, we look at the chain deformation in our MD simulation. Different macroscopic stresses are shown in Figure 8. It can be detected that the folding and unfolding behaviour is compatible with the loading and unloading process.

Figure 8.

Expansion of atomic stress for different macroscopic strains.

Moreover, it can be observed that the atomic stress during loading tends to be higher than the release process (Figure 8(c) vs Figure 8(e)). The inverse correlation between atomic stress is a clear obvious parameter to identify disentanglement.

Chain entanglement evolution

The methodology of Hossain et al. ³⁷ was selected to quantify the entanglement of the system. According to their definition, the entanglement parameter was defined as the total fraction of flexion nodes. First, calculate the angle between two vectors. These vectors are created from each atom to neighbouring atoms that are separated by 10 atoms on the same chain. For example, one vector is generated of atom i with the atom (i+10) and another vector is created between the atom i with the atom (i-10). Second, using a threshold angle of 90 to attribute that the atom is classified as entangled or not. Figure 9(a) displays the distribution of the angle of atoms for the initial and stretched states (strain = 2.5). Figure 9(b) displays the distribution of the entanglement parameter. Five curves contain three different strain rates for PU with repeating chain sequence m15 and an additional two curves for different repeating chain sequences.

Figure 9.

(a) The distribution of the angle of atoms for two states, (b) the variation of the entanglement parameter during a stretch–release process.

When there is no external strain, the distribution of the angle of atoms approximately approaches the Gaussian distribution. When the strain increases, the angle between two vectors increases, consequently the entanglement parameter is decreased monotonically. In contrast, during the release, the entanglement parameter raises. However, it is not able to reach the initial state. This unrecovered parameter is associated with the permanent set and chain length. Conclusively, the chain entanglement is associated with the Mullins effect.

Ogden Pseudo-Elasticity Model

There are enormous efforts that have been dedicated to formulating and explaining the non-linear behaviours of elastomer materials. Ogden and Roxburgh^38,39 (1999) introduced the idealized Mullins effect in elastomer materials using a softening damage variable. It is successful estimates to model the hyperplastic stress behaviour with proper expression. The benefit of this model is comparatively easy to obtain a proper expression for hyperplastic behaviour. It can be directly derived from experimental or molecular dynamic data. These pseudo-elastic models are extensively employed due to adjustability and simplicity.^40,41 In this section, we consider the basic equations of pseudo-elastic material models for HNT/PU composite.

The analytical expression of the strain energy for hyperplastic materials od Ogden^39,42 (1972, 1982) method is given as follows

$W_{Ogden}={\displaystyle \sum _{i=1}^3\frac{2{\mathrm{\mu }}_i}{{\mathrm{\alpha }}_i^2}}\left({\mathrm{\lambda }}_1^{{\mathrm{\alpha }}_i}+{\mathrm{\lambda }}_2^{{\mathrm{\alpha }}_i}+{\mathrm{\lambda }}_3^{{\mathrm{\alpha }}_i}-3\right)$

(2)

The uniaxial stretch is given as follows

${\mathrm{\lambda }}_1=\mathrm{\lambda },\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathrm{\lambda }}_2={\mathrm{\lambda }}_3={\mathrm{\lambda }}^{-1/2}$

(3)

For uniaxial tensile stress, the Ogden (N = 3) strain energy is given at the following equation ⁴³

$W_{Ogden}={\displaystyle \sum _{i=1}^3\frac{2{\mathrm{\mu }}_i}{{\mathrm{\alpha }}_i^2}}\left({\mathrm{\lambda }}^{{\mathrm{\alpha }}_i}+2{\mathrm{\lambda }}^{-{\mathrm{\alpha }}_i/2}-3\right)$

(4)

The stress tensor for anisotropic material is obtained as follows

${\sigma }_1=\sigma ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\sigma }_1={\sigma }_2=0$

(5)

The Mullins effect is estimated by applying a strain energy term and a damage function as defined by equation (6)

$W\left(\mathrm{\lambda },{\mathrm{\lambda }}^{-1/2},\mathrm{\eta }\right)=\mathrm{\eta }{W}_0\left(\mathrm{\lambda },{\mathrm{\lambda }}^{-1/2}\right)+\mathrm{\varphi }\left(\mathrm{\eta }\right)$

(6)

where W0 (λ, λ−1/2) indicate the strain energy in the stretching pathway and W (λ, λ^−1/2, η) for releasing path, and Φ (η) indicate the damage function of the material during the stretch and releasing states. ⁴³ This function determines the energy dissipated during one cyclic loading and it is controlled by variable η which is given by equation (7)

$\eta =1-\frac{1}{r}erf\left(\frac{W^m-W}{m+\mathrm{\beta }{W}^m}\right)$

(7)

where W ^m is the maximum value of W; r, β and m are Mulling effect parameters; and erf(x) is the Gaussian error function obtained as following ⁴³

$erf\left(x\right)=\frac{2}{\sqrt{\pi }}{\displaystyle \underset{0}{\overset{x}{\int }}\exp \left(-y^2\right)}dy$

(8)

In the loading process, η is equal 1, and in the unloading process, η was declined in a range of 0 < η ≤ 1. Finally, with the above introduction, the stresses are given by

$\sigma =\mathrm{\lambda }\frac{\partial W\left(\mathrm{\lambda },{\mathrm{\lambda }}^{-1/2},\mathrm{\eta }\right)}{\partial \mathrm{\lambda }}=\mathrm{\lambda \eta }\frac{\partial W_0}{\partial \mathrm{\lambda }}$

(9)

The MATLAB package was applied to fit the Ogden constitutive parameters through the Levenberg–Marquardt nonlinear least-squares technique. The parameters of the model were fitted to the MD result from the cyclic loading tests with different strain amplitudes of the PU-HNT composite. Table 1 provides our result.

Table 1.

Ogden constitutive parameters for loading and unloading curves of PU/HNT composite.

<Parameters
Ogden constitutive model						Mullins parameters
μ1/MPa	α1	μ2/MPa	α2	μ3/MPa	α3	r	m/MPa	β
0.1829	−9.8825	1.5516	4.9834	3.500	0.4972	1.18	1.2246	0.52

The hyperplastic behaviour and Mullins effect were examined in the ABAQUS software with a rectangular planer model under the strain cycle. As depicted in Figure 10, the simulation curves for loading and cyclic loading are basically in agreement with the MD curves. This model has proposed a reasonable prediction of plastic deformation.

Figure 10.

(a) indicates f loading, and (b) shows loading and unloading MD and Ogden constitutive models of the simulation.

Conclusion

To implement HNT in different applications, it is necessary to consider materials exposed to different loading conditions. HNT has extraordinary properties, in specific, stiffness and strength, which make it attractive as the reinforcing element in a composite material. To achieve the full benefit of these properties, it is necessary to understand the combination with a matrix material. The reason for applying an HNT as reinforced is cross-link between matrix and influence on the polymer chain entanglement and disentanglement. The characteristics of the combination composite material are controlled by the features of the HNT and PU as well as the geometry of the interface. We are determined, MD simulation can demonstrate essential aspects of the Mullin effect and microstructure disentanglement of PU chains which is an effective parameter to identify the residual strain. The analysis of the mechanical response of HNT–PU nanocomposites revealed that the PU chain entanglement is correlated with increasing strain rate. After the first stretch, the disentanglement is not fully reversible, and the permanent set and stress softening begin from the irreversible disentanglement. Generally, the entanglement originates from the interaction interface between HNT and PU and PU chains together. Significant recovery is achieved when the structure is exposed to a higher temperature. The achieve maximum recovery of the residual strain can be explained bond breakage such as between chains, HNT/PU cross-link, and HNT bond break during stretching and recovery them with increasing temperature. Note that the maximum recovery is available when the temperature does not overpass the glass transition temperature. Also, our simulation predicts that with an increase in density of HNT, dismantlement decreases because of increasing chemical bonding density and difficulty to break them; thus, increased reinforcement causes enhanced entanglements. Also, increasing the hard segment content of PU chains increased elastic modulus and tensile strength due to reducing segment mobility. In contrast, increases in hard segment content reduce stress softening. An Ogden constitutive and the theory of pseudo-elasticity have proposed a reasonable prediction of plastic deformation of the HNT in the nanocomposite.

Supplemental Material

Supplemental Material