Enhanced Coarse-Graining of Thermoplastic Polyurethane Elastomer for Multiscale Modeling

To design macroscopic properties of materials, multiscale modeling with molecular simulations has been explored to investigate the structure–property relationship of the materials [1–3]. Multiscale modeling is especially challenged to predict unique properties of polymers and rubbers; e.g., their finite deformation and time-dependent properties. Several researchers have tried to estimate bulk properties of polymers and rubbers using the molecular simulation tools [4–10]. Theodorou and Suter investigated stresses of atoms for uniaxial and shear deformations [4]. They found that the chain topology strongly influences the atomic-level stresses with an active compensation effect of simulation at the atomic level [4]. Fan and Hsu applied a molecular simulation for a uniaxial loading, resulting in a good match of Young's modulus and thermal expansion coefficient but a poor match in Poisson's ratio due to the unavailability to apply Poisson's effect with molecular simulations under uniaxial loading [5].

Other research groups have developed mesoscale/coarse-grained models for polymers/elastomers to predict viscoelastic properties [6–8]. Valavala et al. developed a multiscale tool and used varying initial configurations for full-atomic simulations by setting up upper and lower limits of macroscopic properties [9]. Li et al. developed a constitutive model of a linear polymer using the dynamics of polymer chains at the molecular level, but did not bridge it with a continuum model [10]. A brief literature review of molecular simulation for multiscale modeling of polymers and rubbers is presented in Table 1.

For a computationally viable molecular simulation, coarse-graining methods have been used. For complex structures such as TPE having a long chain with copolymers, we may need a coarse-graining method projecting molecular interaction while keeping computational efficiency.

The objective of this study is to develop an enhanced coarse-graining tool for multiscale modeling of TPU. A coarse-graining method with three beads is suggested with pressure-correction of the force-field. An equivalent continuum model is used to obtain the constitutive relations of TPU by equating the molecular potential energy to the strain energy. By the end of this study, we may answer the following research questions: (1) how sensitive are the iterative procedure and pressure-correction in the enhanced coarse-graining through molecular dynamics (MD) to the effective potentials through inverse Boltzmann method? (2) how useful is the molecular mechanics (MM) for constructing the equivalent energy model for an isochoric loading? and (3) how sensitive are the bulk modulus and Poisson's ratio with the addition of hard segments?

The representative volume element (RVE) of a molecular model of TPU is generated from an equilibrium molecular structure for a force-field. For a computationally viable modeling, we suggest an enhanced coarse-grained model of TPU with three beads constructed from a full-atomic model using the inverse Boltzmann method (IBM). MD simulations are conducted for the coarse-graining of TPU using the large-scale atomic/molecular massively parallel simulator (lammps) software [11] by calibrating density and number average molecular weight of TPU.

The bonded interaction includes valence and cross-coupling. The nonbonded interaction implies van der Walls and electrostatic interaction. For more information on the force-field, the readers are encouraged to read the references [12–16].

To construct an RVE of TPU from a full-atomic model, the chemical structure of TPU is generated using Materials Studio software (Accelrys, San Diego, CA) [17]. TPU is usually synthesized by the reaction between isocyanate and the substances with mobile hydrogen ions, generally polyols, amines, and water [18]. The building block for construction of TPU is illustrated in Fig. 1, where the structure of TPU can be represented by a combination of hard and soft segments. TPU is known as a segmented block copolymer with hard and soft segments where the hard segments have a degree of polymerization in a range of 1–10, and the soft segments have a degree of polymerization in a range of 15–30 [19].

To develop a coarse-grained model that is both computationally efficient and accurate enough to capture the molecular interaction of long polymer chains, we form the segmented block copolymers of TPU with a lower degree of polymerization. The model is constructed with the amorphous module of Materials Studio with a density of 0.85 g/cm³, which is close to the experimentally measured one $(0.97−1.24 g/cm3)$ [20]. An initial length of a cubic simulation box is 7.7 nm. We define beads in the full-atomistic model for coarse-graining. In general, a bead in coarse-graining represents one monomer to describe the interactions among the repeated units, but in this study, we decompose TPU into multiple beads. It is known that the elasticity of TPU is mostly dictated by hard and soft segments, where the hard segments come from diphenylmethanediisocyanate (MDI) and the soft segments come butanediol [18]. Therefore, we consider the segments coming from MDI and butanediol as separate beads, which will enable us to obtain interactions between them. The segments from MDI consist of two phenyl rings connected by a $−CH2$ bond. Due to the long chain of this segment, it may be difficult to capture the internal interactions. Therefore, we divide the hard segment into three beads for an improved coarse-graining of TPU. As illustrated in Fig. 1, the coarse-grained model of TPU has three types of beads: A, B, and C, whereas A + B + C symbolizes the hard segments, B + C represents diisocyanate and only A denotes the soft segment of TPU. Therefore, we construct a full-atomistic chain with a segmented block copolymer structure as $[(ABC)n1−(A)n2−(ABC)n1−(A)n2−(ABC)n1]$⁠, where $n1=n2=4$⁠. This makes 44 beads in a chain, equivalent to 644 individual atoms. Therefore, a total number of atoms of TPU with 50 chains in the full-atomic model is 32,200. To avoid extensive simulation time, we are restricted to the number average molecular weight of the chains as 4684.77 g/mol, which may be too small compared with the physical value $(25,000−100,000 g/mol)$ to capture the phase separation [19]. This will be discussed in the coarse-grained model section for more detail. However, the energy of the original model is minimized using the discover module of Materials Studio for 5000 fs with a 1 fs time step. The model is equilibrated using an NPT ensemble (constant number of atoms, pressure, and temperature) at atmospheric condition (298 K and 1 atm) for 500,000 fs with a 0.1 fs time step. To get an equilibrium model at the molecular level, the full-atomic model is heated to 600 K with a 25 K increment (run for 100,000 fs with a 0.1 fs time step for each 25 K), followed by cooling down to room temperature. At this stage, the density obtained through the simulation is 0.9 g/cm³. Subsequently, an NVT ensemble (constant number of atoms, volume, and temperature) is run for 100,000 fs with a 0.1 fs time step. Twenty trajectories are saved to get the distribution of bond, angle, and radial functions among beads associated with the effective potential functions of the coarse-grained model.

The probability distribution functions among the beads are found from the trajectories of the full-atomic model. The distribution functions for the bond length, $p(l)$ and the angle, $p(θ)$ and the radial distribution of the beads, $g(r)$ are shown in Figs. 2–4, respectively. Most distributions for the bond lengths and angles show a single clear peak, whereas few deviate slightly from this pattern. These single peak distributions of bonds and angles will allow one to develop simplified harmonic potentials on the bond-stretching and angle-bending. The total effective potentials are derived from the bond lengths, angles, and RDFs of atomic simulations with IBM using the following equations [21]:

where $l$ denotes the bond length, the distance between two beads, and $θ$ is the bending angle of three adjacent beads. $P(l)$ and $P(θ)$ are the probability distributions of bond lengths and bending angles, respectively. $g(r)$ is the RDF among the beads apart from the distance $r$⁠. $Ubond(l)$⁠, $Uangle(θ)$⁠, and $Unonb(r)$ are the potential energies for bond-stretching, angle-bending, and pair-wise nonbonded potentials, respectively. $T$ is the temperature in Kelvin (K) and $kB(=0.001987 (kcal/mole))$ is the Boltzmann constant.

where $kstretch$ and $kbend$ are the spring constants for the bond-stretching and the angle-bending, respectively. $l0$ and $θ0$ are the equilibrium bond length and angle, respectively, for the harmonic motion. $ε$ is the depth of the potential well, and $σ0$ is the finite distance at which the interparticle potential is zero. The IBM-derived potentials with the fitted functions for bonds, angles, and pairs are depicted in Figs. 2–4, respectively. The torsional/dihedral potentials of the beads are not presented in this work because they are not the major components for a coarse-grained model [22,23].

The chain of TPU is constructed with a block segmented copolymer as $[(A)n1−(ABC)n2−(A)n1−(BC)n3−(A)n1−(ABC)n2]n$⁠, whereas $n1=14$⁠, $n2=2$⁠, $n3=1$⁠, and $n=5$⁠. Therefore, one chain has 280 beads composed of 230 A-type beads, 25 B-type beads, and 25 C-type beads. Molecular weights of A, B, and C beads are calculated from a full-atomic model; 72.1 g/mol, 135.2 g/mol, and 135.1 g/mol, respectively. One hundred chains are constructed in the model (28,000 beads in total) with the same configuration and the same molecular weight by considering the monodispersity for simplicity. Thus, the number average molecular weight of the coarse-grained model is 23,341 g/mol, where we attempt to keep it as close as 25,000 g/mol. The number average molecular weight of a typical linear polyurethane (PU) elastomer varies from 25,000 to 100,000 g/mol, and its physical properties are known to be rarely affected by a variation of molecular weight within this range [19]. Therefore, we consider this specified coarse-grained model as a reference model of TPU [24].

The IBM-derived potentials are defined for the coarse-grained model. An NVE ensemble (constant number of atoms, volume, and energy) is used for the initial unstable molecular structure. Temperature is controlled at 600 K using the Langevin thermostat. The simulation box is defined slightly bigger than the one required to give some space for relaxing beads. An NPT ensemble is run at 600 K, cooling to 298 K using another NPT ensemble. Following the NPT relaxation at 298 K, the simulation box is deformed to the desired density (0.90 g/cm³) using an NVT ensemble. The model is equilibrated with another NVT at 298 K while keeping the density as close as 0.90 g/cm³. The corresponding pressure developed at this equilibrium model is 482 atm. The deviation appears to be attributed to the fact that IBM with pressure-correction does not fully generate nonbonded effective potentials of the coarse-grained model. The nonbonded effective potentials need pressure-correction for keeping the same thermodynamic states of the two models—the full-atomic and coarse-grained models. Before moving on to the development of pressure-corrected force-fields, we compare the RDFs of the coarse-grained model under the fixed density with those of full-atomic model. It should be noted that three beads (A, B, and C) have six different RDFs which are plotted in Fig. 5. They are in good agreement even though both the full-atomic and coarse-grained models do not possess the same thermodynamic states.

where $A$⁠, a constant, is, in general, $±0.1 kBT$⁠; but the numeric number 0.1 is tunable for coarse or fine matching. As the model generates an enormous pressure to maintain the density, the iteration will lead to more attraction of energy/deeper potential well of pair potentials. After running several iterations, the required pressure comes close to 1 atm. However, the exact matching is difficult due to fluctuation. To confirm the convergence, the model is equilibrated using an NPT ensemble at 1 atm for 1 ns. The final density is found to be $0.86±0.0025 g/cm3$⁠, which is comparable with both the density of TPU $(0.90 g/cm3)$ obtained from a full-atomic simulation and the physical density value of TPU $(0.97−1.24 g/cm3)$ obtained from the experiment [20]. The final RDFs are compared with the target RDFs as shown in Fig. 5. It is expected that the addition of a linear term to the nonbonded potential will compensate the pressure remaining the RDFs the same. However, Fig. 5 shows that there is a marginal increment in the peak of RDFs in some cases: B–B and C–C are distinct but others are in an acceptable range. This may be improved by developing a more realistic pressure-correction technique rather than just adding/subtracting a simple linear term to the nonbonded interactions. The nonbonded potentials with pressure-correction are also compared with the IBM-derived initial potentials in Fig. 5.

To emphasize the significance of pressure-correction for coarse-graining of TPU, we generate a density variation as a function of temperature with and without pressure-correction as illustrated in Fig. 6. The initial structure of TPU compacted to a density of $0.9g/cm3$ at 298 K, heated up to 600 K and cooled down to 298 K with an NPT ensemble with (i) pressure-correction and without (ii) pressure-correction. It should be noted that the force-field without pressure-correction means initial IBM interactions.

As can be seen in Fig. 6, when they are cooled down to 298 K, the structure with force-field without pressure-correction does not come back to the initial density due to weak nonbonded interaction, but it does with the force-field with pressure-correction. The force-field with pressure-correction generates density variation with three phases of TPU—rubbery, transition, and glassy phases (Fig. 7) as the temperature cools down from 600 K. However, the force-field without pressure-correction does not capture the phase transition. This study shows the importance of pressure-correction while building a coarse-grained model of TPU.

Mixing and demixing of polymer chains are strongly dependent on the stereochemical composition and attractive intermolecular interaction of pairs [26,27]. Phase separation of chemically reactive binary mixtures [28] and thermally activated spontaneous and homogeneous nucleation of grains of undercooled iron melt have been investigated with MD simulations [29]. In addition to chemical reaction, the physical arrangement of molecules, types of soft segments and molecular weight, the length of hard segments, thermodynamic incompatibility of segments, hydrogen-bond between the urethane segments, and symmetry of the diisocyanate are known to play an important role to the phase separation of TPU [30]. The discrete hard domains are known to be dispersed in a soft matrix where the hard segments-enriched microphase works as the physical cross-linking of the elastomer [30]. The hard phase enriched region is conspicuously visible in the soft matrix at 1 atm with the pressure-corrected force-field in Fig. 6. A further parametric study is performed in the sections, Case Study of Phase Separation, and Design of Materials at Molecular Level.

The suggested coarse-graining technique is robust. We applied the molecular simulation in the full-atomic models with varying number of monomers and chains during coarse-graining. Bond and angle distributions follow the same path after coarse-graining regardless of its system size. Nonbonded parameters change a little in some cases. However, it becomes negligible after pressure-correction. The overall simulation results are repeatable with a small deviation of ∼1%.

In addition to the reference model of TPU as specified earlier, five additional models for TPU are prepared to observe the effect of hard and soft segments with modified configurations of the segmented copolymer structure as $[(ABC)n1−(A)n2]n$ where $ABC$ and $A$ denote the hard and soft segments, respectively. A total number of chains is kept fixed as 100 for all the models while maintaining a number average molecular weight as close as 25,000 g/mol, as we did earlier in this study. Under this restriction, we can perform a parametric study with $n1$⁠, $n2$⁠, and $n$⁠. The parameters of the five models are reported in Table 2. Note that the reference model in the section, Coarse-Grained Model of TPU, has 34.9 wt.% of hard segments.

All the models are cooled and equilibrated at 1 atm to observe the effect of the hard and soft segments on the phase separation in the equilibrium structures. If the wt.% of the soft segments decreases, the phase mixing increases [30]. Figure 8 shows that the structure having 34.9% hard segment has clear boundaries between the hard and soft segments. The hard segments interpenetrate into each other as the wt.% of the hard segments increase (Fig. 8). The dimensions of the separation domains are known to be the order of 5–100 nm from the experiment [30], whereas Fig. 8(a) shows the size of hard and soft segments as 1.75 nm and 2.63 nm, respectively.

Elastic properties of TPU are obtained using an equivalent continuum model with finite deformation [31], which is suitable for large and amorphous molecular structures. First, an RVE of the coarse-grained model of TPU is chosen to describe the bulk properties of TPU followed by the development of a constitutive equation using the equivalent continuum modeling method [31]. The energies for deformations at the molecular and continuum levels are equated under identical sets of boundary conditions to determine parameters for the constitutive equations.

Note that $c1Ω1$ and $c2Ω2$ denote the strain energy densities corresponding to volumetric and isochoric deformations, respectively. Bulk and shear moduli are functions of $c1$ and $c2$⁠, respectively, but $c1$ and $c2$⁠, respectively, do not directly denote bulk and shear moduli.

where $I$ is the identity tensor. It should be noted that $S$ contains the material parameters $c1$ and $c2$⁠, which need to be determined by equating Eqs. (22) and (23) to the corresponding potential energy densities obtained from molecular simulations.

where $Utotal|iso$ and $Utotal|iso0$ are the potential energies of a molecular model for isochoric deformation modes before and after deformation, respectively. $c2$ is obtained by equating Eqs. (23) and (26) for a set of isochoric strains applied: $c2$ is the slope of the $Ψm|iso$  −  $Ω2$ plot. Green–Lagrangian strain tensor, $E=(1/2)(C−I)$ is used for deformation of an RVE for the molecular model. During the molecular simulations, volumetric strains (⁠$E11=E22=E33)$ of $0.25%$⁠, $0.5%$⁠, $0.75%$⁠, and 1% are applied for a volumetric deformation, and three-dimensional shear strains (⁠$γ23=γ13=γ12, where γij=2Eij when i≠j)$ of $0.25%$⁠, $0.5%$⁠, $0.75%$⁠, and 1% are applied for an isochoric deformation using MM. Before applying the deformation using MM, the model is further equilibrated with an NPT ensemble for 1 ns. Then, the energy of the model is minimized using the conjugate gradient algorithm with stopping criteria: $1.0×10−6$ as a stopping tolerance for energy, $1.0×10−6$ as a stopping tolerance for force, 10,000 as maximum iterations of the minimizer [11]. A fix box/relax option in lammps is used to relax the box and make sure that no external pressure is working on the model. The “fix viscous” command in lammps is also used to dissipate the viscous energy of the beads. After achieving convergence, the model is further NVT relaxed for 100,000 fs. The energy of the model is minimized again with the same stopping criteria of the conjugate gradient algorithm with an NVT ensemble instead of an option with fix box/relax in lammps, which helps further to minimize the energy.

After obtaining the potential energy values, the potential energy densities $Ψm|vol$ and $Ψm|iso$ are obtained from Eqs. (25) and (26), respectively. $Ω1$ and $Ω2$ expressed with invariant forms in Eqs. (22) and (23), respectively, are obtained for the corresponding deformations.

Therefore, $c1$ and $c2$ are ready to get from the slopes of the curves of $Ψm|vol−Ω1$ and $Ψm|iso−Ω2$⁠, respectively. Two different temperatures (1 K and 100 K) are employed for the molecular mechanics simulations to obtain $c1$ and $c2$⁠. We found that a temperature higher than 100 K are not preferable to generate potential energy profiles with MM. The curves on $Ψm|vol−Ω1$ and $Ψm|iso−Ω2$ are depicted in Fig. 9. $c1$⁠, the slope of a linear regression of volumetric deformation in Fig. 9 is obtained to be 145 MPa at 1 K. Interestingly, we received a higher $c1$ of 428 MPa at 100 K, which is thought be due to local minima. $c2$⁠, the slope of isochoric deformation in Fig. 9, is obtained as 5.75 MPa at 1 K and 2.28 MPa at 100 K, which is a reasonable result with an increase in temperature.

Once the total strain energy density function $(Ψc)$ is obtained from Eq. (21), the stress on a macroscale can be estimated for any mechanical loading modes using Eq. (24). Macroscopic mechanical properties under volumetric expansion and isochoric (three-dimensional shear) loading are predicted from this multiscale scheme.

The obtained K, G, E, and $ν$ are reported in Table 3. For a comparison purpose with an experiment, a uniaxial tensile test of TPU is conducted following ASTM D412 [35]. Young's moduli measured from a uniaxial tensile test are 32.14–40.33 MPa for strain rates of 0.068–0.342 s⁻¹. It is reasonable to assume Poisson's ratio of typical TPUs as 0.48–0.5 [36]. The values from the MM simulations and experiment are compared in Table 3. The MM simulation overestimates Young's modulus and shear modulus, but bulk modulus and Poisson's ratio show an excellent agreement with the experimental results. Most amorphous polymers are operated at room temperature, far below their glass transition temperature, which may envision as a “frozen” in a local potential energy equilibrium state [9]. On the other hand, TPU has a glass transition temperature of 200–250 K [20], which is below the room temperature. The low glass temperature indicates the presence of a metastable phase at the room temperature. We will consider the metastable state to predict the macroscopic mechanical properties of TPU for our future study by investigating an equivalent glass transition state at the molecular model with MD simulations.

It is well known that the bulk modulus of elastomers is much greater than the shear modulus. The results from the equivalent continuum approach of the coarse-grained model with molecular simulations exactly project this trend (Table 3).

We may use the multiscale modeling technique as a materials designing tool and investigate a structure–property relationship by varying hard and soft segments of TPU. All of the five models, as prepared with the modified configurations in the Case Study of Phase Separation section, are equilibrated and relaxed with NPT ensembles followed by deformations under volumetric expansion $(E11=E22=E33=1%)$ and isochoric loading $(γ12=γ23=γ31=1%)$ using MM simulations. Using this multiscale scheme, bulk modulus (⁠$K$⁠) and Poisson's ratio (⁠$ν$⁠), which showed a good match with the experimental data in the Comparison with Experimental Results section, are estimated and plotted as a function of density in Fig. 10. The density of TPU apparently increases with an addition of hard segments; an increase of density up to 46% with a rise of wt.% of the hard parts up to 50% (Fig. 10). Figure 10 shows an 8% increase of $K$ with a rise of wt.% of the hard segments up to 50%. We may design the Poisson's ratio with an addition of the hard segments. Figure 10 shows that a change of Poisson's ratio from 0.45 to 0.43 with an increase of wt.% of the hard segments up to 50%.

An MM simulation captures the static behavior of molecules well simply because there is no need to consider the contribution of kinetic energy to stress used for developing potential energy, which is a required step for molecular dynamics. However, due to the simplicity without considering time scale, there is a limitation of MM in implementing rate dependent properties.

We proposed an enhanced coarse-grained model of TPU for the multiscale modeling of copolymers with long chains. The suggested coarse-graining model built with three monomer units of TPU enables us to capture molecular interaction while keeping computational efficiency. The coarse-grained model equipped with pressure-correction of the force-field improves the molecular model to possess similar properties as a bulk model of TPU—varying density with temperature, a close density value of TPU at 1 atm, and the phase separation.

Equating potential energy densities of the coarse-grained model to the strain energy functions of the continuum model at volumetric and isochoric deformation modes, bulk and shear moduli of TPU are directly obtained and used to estimate Young's modulus and Poisson's ratio. Obtained bulk modulus and estimated Poisson's ratio show an excellent match with the experimental results of TPU. The molecular simulation with the coarse-grained model of TPU demonstrates much greater bulk modulus than the shear modulus, which is a typical property of elastomers. Modifying the coarse-grained model of TPU, we successfully demonstrated material design for engineering bulk modulus and Poisson's ratio by varying hard and soft segments of TPU at the molecular level.

The proposed coarse-graining tool will pave a new way to explore the multiscale modeling of copolymers with long chains and can be directly applied to the multiscale modeling of other thermoplastic elastomers; e.g., styrenic block copolymers, thermoplastic olefins, thermoplastic copolyester, etc.

This study was supported by the National Science Foundation (NSF) Industry/University Cooperative Research Center (I/UCRC) Program—Center for Tire Research: Planning Grant (Award No. IIP-1338921) and Research Seed Grant (RSG) of the Vice President for Research and Economic Development (Award No. 63244) at UNT. Computational resources were provided by UNT's High Performance Computing Services.