Molecular Dynamics Simulations of Polymer Grouting Material: Its Mechanical Behavior under Uniaxial Tension, Cyclic Tensile Loading, and Stress Relaxation

At present, polyurethane (PU) has been extensively used as a grouting material in civil engineering. The mechanical properties of PU are the key to achieving the desirable grouting e ﬀ ect. This study presents the research results of the mechanical behavior of PU matrix under tensile, successive cyclic tensile, and stress relaxation at the nanoscale, using the coarse-grained molecular dynamics simulation method. The in ﬂ uences of the number of molecule chains and strain rate on the tensile mechanical properties are discussed, and the tensile deformation mechanism of PU matrix is revealed. The tensile strength of PU matrix is independent of loading path, and after yielding, the strain of PU matrix contains the elastic strain, plastic strain, and viscous strain. In the stress relaxation process, the evolution of the axial stress is mainly caused by the varied van der Waals interactions. The stress relaxation behavior of PU matrix can be described by the viscoelastic model consisting of one elastic element in parallel with one Maxwell element.


Introduction
Polyurethane (PU) [1], the product of polymerization reaction between isocyanate and polyols, is a versatile class of copolymer used in coatings [2], armors [3], and composites [4] for its properties of high durability, high toughness, high capacity to absorb and dissipate energy, and low density. At present, by virtue of the characteristics of the strong permeability resistance [5], fast-setting [6], high expansion force [7], and eco-friendliness [8], the foamed polyurethane has been increasingly adopted as grouting material in civil engineering. As shown in Figure 1, the representative applications of the polymer grouting material include filling soil cavities to lift the sedimentary underground pipelines [9,10] or foundations [11,12] and constructing barriers to prevent water seepage [13,14] or the leakage of pollutants [14,15]. In these applications, the two raw materials, isocyanate and polyols, are simultaneously injected into soil at the constant mass ratio as alternatives to cement-based grouts [16,17], and after polymerization reaction, the PU matrices will be retained in the soil mass to play the reinforcing role [9,10,16]. For example, it has been reported that when the void surrounding the underground pipe was filled through PU grouting, under the same vertical load, the strain and bending moment in the pipe was declined by about 76.9% and 67.9% [10]. The mechanical properties of PU matrices are the key to achieving the desirable grouting effect, and the polymer grouting materials with insufficient strength may cause the refailure of the reinforced engineering structures [9][10][11]18].
To date, numerous experimental studies have been carried out to investigate the tensile, compressive, flexural, and shear strength of polymer grouting materials [19][20][21][22][23]. For instance, the tensile strength and compressive strength of the plain PU composites are about 44.4 MPa and 15.5 MPa, and the appropriate addition of PVA fibers can improve the mechanical properties of PU composites by about 4.5-22.3% and 9.0-12.2% [19]; in Hussain et al.'s study [21], it was found that the experimentally tested flexural strength of PU composites could be about 11.0 MPa; in the study of Zheng et al. [22], it was found that the flexural strength and compressive strength of PU grouting materials with density of 0.1-0.5 g/cm 3 were in the range of about 1.2-18.1 MPa and 1.0-7.5 MPa, respectively. It has been reported that, according to the micro morphology observations [6], PU grouting materials are porous materials containing numerous hollow closed foams, making them perform with obvious discontinuity characteristics. The foams are packed together, and their walls are made of amorphous PU [6], and under the external loading, it is the deformation of the foam walls, i.e., the PU matrix, that cause the failure of polymer grouting materials. By incorporating the beam deformation theory and the membrane stretching theory, the Gibson-Ashby equation that reflects the mechanical behavior of PU grouting materials has been proposed [24]. The Gibson-Ashby equation shows that knowing the mechanical properties of the PU matrix from the molecular scale is the foundation of the theoretical research on the mechanical properties of polymer grouting materials [24], which, however, is rarely seen at present.
In the past decades, the molecular dynamics (MD) simulation has been proven to be a reliable and promising tool to explore the nanoscale mechanical behavior of various materials. For example, studies on the mechanical behavior of amorphous polymer materials [25,26], the microscopic deformation characteristics of inorganic minerals or metals [27,28], the tensile strength of a single carbon nanotube [29,30] or one layer of graphene [31,32], and the nanoscale shear enhancing mechanism of carbon nanotubes in cement [33], all using the MD simulation method, have been reported. A large amount of significant efforts to develop the united atom [34,35] and coarse-grained [36][37][38] have been carried out. When the united atom or coarse-grained MD simulation methods are used, groups of atoms are lumped into single interaction sites with the simplified interaction potentials, and thus, the total number of degrees of freedom of the simulated models can be reduced. The united atom and coarse-grained MD methods require less computation, due to which, the time and length scales of the models can be increased to simulate the bulk behavior. With the united atom or coarse-grained MD simulation methods, the nanoscale mechanical properties and deformation mechanism of amorphous polyethylene under the tensilecompressive load [25] or pure tensile load [35] have been understood. In the study of Park et al. [38], the thermomechanical behavior of the shape-memory polyurethane copolymers was also studied using the coarse-grained MD simulation method. Studies on the nanoscale modulus [39], cohesion strength [40], and self-assembly behavior [41] of PU material with coarse-grained MD simulation method have also been reported. The coarse-grained MD simulation has provided a new approach to understand the mechanical behavior of PU matrix from the molecular perspective.
In this study, using the coarse-grained MD simulation method, the mechanical behavior of the amorphous PU matrix under tensile was firstly studied, and then, both the successive cyclic tensile simulation and stress relaxation simulation were carried out. This study is aimed at addressing the mechanical properties of the amorphous PU matrices at the nanoscale under different types of loadings that are commonly seen in geotechnical engineering [9,11,18,42]. The reason why the tensile simulation was firstly carried out is that under external load, the failure of the PU foams is mainly caused by the tension deformation [6,24,25]. This study is organized as follows: the construction of the coarsegrained model of amorphous PU and the simulation details are included in Section 2. In Section 3, both the influences of the number of the molecular chains and the strain rate on the tensile strength of PU matrix, as well as the tensile deformation mechanism of PU matrix, are discussed. The mechanical behavior of PU matrix under cyclic tensile loading and the stress relaxation characteristics of PU matrix under different tensile strains are included in Section 3. The main conclusions are summarized in Section 4.  2 Geofluids [38], the hard segment of the PU molecular was divided into two different beads, one of which includes two halves of aromatic rings, and the remaining component of the benzenes and the urethane linkage were represented with another bead. As for the coarse-grained PU model proposed by Uddin et al. [39,43], it has the physical and mechanical properties that are a good match with the experimental results; therefore, it was adopted in this study. More details are included as follows. As shown in Figure 2(a), a single PU molecule chain is composed of the hard and soft segments, where the formers come from diphenylmethane diisocyanate (MDI) and butanediol [44]. The hard segments were treated as separate beads to obtain interactions between them. The segments from MDI consist of two phenyl rings connected by a -CH 2 group, and the internal interactions are difficult to capture; therefore, the hard segment was divided into three types of beads [39,43]. The -(CH 2 ) 4 -O group, -(CONH)-C 6 H 6 -CH 2 group, and the -(C 6 H 6 )-NHCOO group are represented by type A bead, type B bead, and type C bead, respectively. The molecular weights of A, B, and C beads are 72.1 g/mol, 135.2 g/mol, and 135.1 g/mol, respectively. The hard segments and soft segments in the PU copolymers often have a degree of polymerization in the range of 1-10 and 15-30, respectively [44]. Based on the commonly seen polymerization degrees of the hard segments and soft segments, the coarse-grained PU model (as shown in Figure 2(b)) was built. Figure 2(b) shows one-fifth of the coarse-grained PU molecule chain, and an entire coarsegrained PU chain contains 280 beads in total, composed of 230 A beads, 25 B beads, and 25 C beads. The polymerization degrees of the hard and soft segments in the coarsegrained model are 2 and 14, respectively, and the molecular weight is 23340 g/mol. The full-atomic polyurethane model mapped from the coarse-grained model is consistent with the typical PU molecule formula presented in literature [44], and thus, it is representative. With this coarse-grained method, for a typical PU matrix with the polymerization of 100, the total number of particles is reduced to 28000 from 389000, making it possible to compute relatively large systems at a physically meaningful simulation scale.

Simulation Models and Methods
The potential energy of the coarse-grained PU model contains three terms, the bond stretching, the bond angle bending, and the van der Waals interactions, whereas the term associated with the dihedral angle torsion is negligible [37]; therefore, the potential energy can be expressed as follows: The function forms of E bond , E angle , and E nonbond in Equation (1) are as follows: where K bond and K angle are the stiffness constants for the bond length and bond angle potential energy, respectively; r 0 and θ 0 are the equilibrium bond length and equilibrium bond angle, respectively; r and θ are the distance between two atoms and the angle between two bonds, respectively; ε is the energy well depth, σ is the distance between two nonbonded atoms at which the potential energy is zero; and r c is the cutoff distance which is defined as 2.5σ [44].
In coarse-grained MD simulations [45,46], Equation (2) and Equation (3) are the most commonly used functions to describe the variations of potential energy associated with bond length and bond angle, and the van der Waals interactions are often expressed by the Lennard-Jones 12-6 function (Equation (4)). Based on the trajectories of the fullatomic polyurethane equilibrium model developed using the polymer consistent force field, the probability distribution functions among the beads can be derived with the inverse Boltzmann (IBM) method [39,43]. The center of mass of each group of atoms is defined as the center of mass of the corresponding bead. Then, the potentials on the bond stretching and the angle bending are fitted with Equation (2) and Equation (3), and the pair potentials are fitted with Equation (4), through which the coarse-grained force field parameters can be obtained. The coarse-grained force field parameters are listed in Tables 1-3.

Constructing the Coarse-Grained Model of Amorphous
PU Matrix. Figure 3 shows the construction procedure of the coarse-grained model. First, an initially large size model was built by replicating the coarse-grained molecule chain, and then, the initial unstable model was relaxed at 600 K using a Langevin thermostat and the NVE ensemble for 10 ns. After that, the model was relaxed for another 10 ns at 600 K and then cooled to 298 K at a rate of 25 K per 50 ps, followed by a steady-state relaxation at 298 K for 10 ns, both with the NPT ensemble and the external pressure of one atmosphere (0.1 MPa). The simulation box was deformed to obtain a model with the desired density, after which, the model was equilibrated for 10 ns while keeping the density unchanged, both using the NVT ensemble at 298 K. The Berendsen thermostat was used to reset the temperature in both the NPT and NVT simulations, and the Berendsen barostat was used to reset the pressure of the system in NPT simulation.
The finally obtained coarse-grained PU models have a density of about 0.86 g/cm 3 , which is comparable with the density obtained from the full-atomic simulation (about 0.90 g/cm 3 ) [39,43] and the experimentally obtained density (0.97-1.24 g/cm 3 ) [47]. In addition, the nanoscaled mechanical properties of the coarse-grained PU model are also comparable with the experimental ones [39,43]. The coarsegrained models are large collections of PU chains that are randomly jumbled together. The cubic simulation cells with the side lengths of 76.9, 102.3, and 162.5 Å were established, and each was composed of 10, 25, and 100 PU chains, respectively ( Figure 3). With these three different models, 3 Geofluids the influence of the number of chains on the tensile strength can be studied.

Simulation Details.
For conducting the tensile simulation, a constant engineering strain rate was applied to the coarse-grained models with the NVT ensemble at the temperature of 298 K. The periodic boundary condition was adopted, and the Berendsen thermostat was used to reset the temperature of the system. According to the symmetric pressure tensor and the geometric size of the models, the tensile stress (σ) and engineering strain (ε) were computed, and the potential energy contributions associated with different terms were tracked as a function of ε. When presenting our results, true stress (σ T ) and true strain (ε T ) were used as alternatives of σ and ε, respectively, and they were calculated as follows: The strain rate of 5 × 10 10 /s was applied to the three models with different side lengths until ε reached 1.0, and the tensile strain was along the x-axis. For the configuration containing 100 PU molecule chains, the uniaxial tension simulation at two other strain rates (2 × 10 11 /s and 1 × 10 10 /s) was also conducted to study the influences of strain rate. The tensile rate adopted here is in the range suggested by previous studies [25,35].
The successive cyclic tension loading simulation was conducted on the 100-chain coarse-grained model. Both the loading and unloading rates of the successive cyclic tensile simulation were 5 × 10 10 /s, and a total of 7 cycles were simulated. For each cycle, the model was unloaded to zero stress at first and then reloaded to a larger strain than the previous cycle. ε T at each unloading point was 0.04, 0.09, 0.18, 0.24, 0.32, 0.45, and 0.60, respectively.
Finally, the stress relaxation simulation was conducted on the 100-chain coarse-grained model. The model was stretched to different true strains (0.05, 0.1, 0.2, 0.4, 0.6, and 0.8) at the rate of 5 × 10 10 /s and then relaxed for a period of time while keeping the strain constant. The evolutions of the axial stress and potential energy of the model with the relaxation time were recorded. Both the model establishment and the tensile simulations were carried out by the LAMMPS package [48], and the open visualization tool OVITO was adopted for molecular visualization [49].

Effects of the Number of PU Chains on the Tensile Mechanical Properties. Figures 4(a)-4(c)
show the variation laws of the stress of the coarse-grained PU models with different numbers of molecular chains under uniaxial tension, from which one can see that when ε T was lower than 0.03, the stress increased almost linearly with the increasing strain, corresponding to the elastic deformation behavior. When ε T was over 0.03, the variation rate of the tensile stress started to decline, implying that plastic deformation had occurred to the model, and at ε T of about 0.20 (yield strain), the tensile stress achieved a plateau. It can be found that the number of molecular chains has little influence on the yield strain or yield strength (both are about 0.20 and 30 MPa). The stress-strain curve becomes smoother with a higher number of PU chains, suggesting that the simulations moved closer to the bulk behavior. Therefore, in all the subsequent simulations, the coarse-grained PU model with 100 molecular chains was used. The snapshots of the deformation process of the configurations containing 10 PU molecule chains are presented in Figure 4(d). It can be seen that the deformation of the model mainly consists of two parts: the slippage between the PU chains and the stretching of PU chains themselves. Before the yielding, the stretching deformation of each independent PU molecule chain was less significant than the slippage between them, whereas it became more and more obvious with the increasing strain after yielding. The slippage deformation occurred mainly along the direction of tensile strain.

Effects of Strain
Rate on the Tensile Mechanical Properties. Figure 5 shows the variations of the stress of the coarsegrained PU matrix containing 100 molecule chains under uniaxial tension at the strain rate of 2 × 10 11 /s and 1 × 10 10 /s. The curves follow the same laws as discussed in Section 3.1 except for the different yield strengths and different yield strains. According to Figures 4(c) and 5, one can see that the yield strength of the model declines with the declined strain rate. At the strain rates of 2 × 10 11 /s, 5 × 10 10 /s, and 1 × 10 10 /s, the yield strength was about 70, 30, and 16 MPa, respectively. On the contrary, with the decreasing strain rate, the yield strain increased gradually. For example, the yield strain was about 0.20 and 0.22 at the strain rate of 5 × 10 10 /s and 1 × 10 10 /s, respectively, which were about double that at the strain rate of 2 × 10 11 /s. Moreover, it can also be seen that when a higher strain rate was adopted, the stress-strain curves became smoother. Figure 6 illustrates the evolution of potential energy (PE total ) of the coarse-grained PU matrix at different tensile strain rates, corresponding to the stress-strain behavior presented in Figures 4(c) and 5. The evolutions of individual components, i.e., the bond potential energy (PE bond ), bond angle potential energy (PE angle ), and van der Waals interactions (PE nonbond ), are also included in Figure 6. E 0 represents the potential energy at initial time.
From Figure 6(a), it can be seen that when the axial strain increased at a rate of 2 × 10 11 /s, both PE total and PE nonbond increased gradually, and the increasing trends are curves with a downward opening, and that the PE total curve and PE nonbond curves almost coincided until ε T reached to about 0.10 (the yield strain). PE total was about 6000 kcal/mol at the yield strain, and after yielding, PE total became lower than P E nonbond ; meanwhile, the PE angle term started to decline in an almost linear way. PE angle was about -6000 kcal/mol at ε T of 0.70. During the tensile procedure, PE bond increased with the increasing strain, but its variation was lower than those of the other two terms (with the maximum increment of lower than 1500 kcal/mol).
The changed PE nonbond , PE angle , and PE bond are primarily attributed to the slippage between different molecules, bond

Geofluids
angle bending behavior, and bond stretching behavior, respectively. Figure 6(a) shows that, at the strain rate of 2 × 10 11 /s, under uniaxial tension, only the slippage between different PU molecules occurred before yielding, but after yielding, in addition to the slippage between different PU molecules that still remained, on each PU molecule, the bond angles started to bend, and the bonds were stretched to synthetically accommodate deformation. The difference between PE nonbond and PE total after yielding was dissipated by the bond angle bending behavior and bond stretching behavior, but primarily via the former. Figures 6(b) and 6(c) show that whether at the strain rate of 5 × 10 10 /s or 1 × 10 10 /s, there was little change in the overall evolution laws of PE total and PE nonbond , but PE angle started to decline from the very beginning of the simulation, and that PE bond was always equal to zero. This indicates that at a lower strain rate, as soon as the tensile load was applied, the bond angles displayed bending behavior, but the bond lengths were little influenced. Figure 6 also shows that the change of total potential energy at yielding is independent of the strain rate. At the three loading rates, both PE total at yielding are about 6000 kcal/mol, but after yielding, the PE total of the model sustaining a lower strain rate has a lower value, which is due to that PE nonbond has significantly declined. With the  Figure 4: (a-c) Stress-strain curves of coarse-grained PU models with different numbers of chains under tensile at a strain rate of 5 × 10 10 /s and (d) the evolution of molecular configuration containing 10 PU chains at different strains. 6 Geofluids decreasing strain rate, the increased yield strain (Figures 4(c) and 5) leads to the increased volume strain of the model. Considering that the yield strength of the model is positively dependent on the potential energy but negatively influenced by the volume [50], when the strain rate declines, PE total is almost unchanged, but the increased volume strain of the coarse-grained PU models at yielding is the main reason for the lower yield strength.     Figures 4(d) and 6 show that the slippage between different molecules mainly dominates the deformation of the PU matrix. The molecules in the original model are tightly compacted together, and the relative slippage between them would occur once the strain was applied (Figure 4(d)), and this would increase the average spacing between different molecules and thus change the van der Waals interactions ( Figure 6). The increased molecular spacing also provides some space for the restoration of the deformed bond angles; therefore, the potential energy associated with the bond angles would decline gradually with the increasing strain. When the model is deformed to a larger size, the average distance between different molecules would be increased, which would result in the less and less pronounced van der Waals interactions (Figure 6). At a lower strain rate, the amount of slippage per unit time between different molecules is smaller, and the time for the variation of bond angles is longer. The bond angles have more time to revert to their equilibrium, and thus, the PE angle begins to decline as soon as the strain appears. However, most of the bond angles are still not in their equilibrium state. This can be illustrated by Figure 7, from which one can see that even at ε T of 0.60, for these 6 types of different bond angles, over 90% of them were smaller than the equilibrium bond angles. Therefore, the PE angle term continuously and steadily declines.
To find out which types of bond angles reverted the most, the evaluations of PE angle of each type of bond angle of the model at a tension rate of 2 × 10 11 /s were calculated, and they are drawn in Figure 8.
From  Figure 6(a). Moreover, even the B-C-A bond angles have the minimum K angle than the others, the total change in PE angle of A-A-A bond angles is more significant than that of the B-C-A bond angles. For example, at ε T of 0.40, the total potential energy of the A-A-A bond angles decreased by about 2400 kcal/mol, which is over 6 times that of B-C-A bond angles (about 380 kcal/mol). This is because the total number of the A-A-A bond angles is the maximum (190 per PU chain); therefore, the total variation of potential energy of the A-A-A bond angles is the most significant. This implies that the change in PE angle depends on both the bond angle stiffness and the total number of bond angles [35]. The A-A-A bond angles are mapped to the bond angles in the soft segment in the all-atom model (Figure 2(a)), which implies that the deformation of the polyurethane molecule mainly occurs in the soft segment.

Mechanical Behavior under Successive Cyclic Tension
Loading. Figure 9 shows the stress-strain curve of the coarse-grained PU matrix under successive tensile loadingunloading cycles. It can be seen that, generally, the contour of the stress-strain curve of the coarse-grained PU sustaining the successive loading-unloading cycles follows that of the tensile stress-strain curve. The strength after each reloading was almost equal to the tensile strength at the same strain. After having passed through the unloading point, the reloading curves precisely follow the original stress-strain curve, which is known as "stress memory" behavior. The "stress memory" phenomenon has also been observed in other polymeric materials like high-density polyethylene [25]. This observation is consistent with the experimental results [51], and it indicates that the strength of polyurethane under tension is independent of the loading path. Figure 9 also shows that, before yielding, the stress curves of unloading and subsequent reloading change linearly, and they are almost coincident, but after yielding, they display the nonlinear variation law. More specifically, the stress curves of reloading increase sharply and linearly within a small strain range, and then, they turn to increase at a much lower rate up to the value of the unloading point. Taking the fourth loop as an example, the slope of the reloading stress curve was about 1850 MPa per strain when ε T increased from 0.22 to 0.24, but it declined to about 90 MPa prestrain in the ε T range of 0.24 to 0.26. The unloading curves vary similarly but in a reverse direction. Due to the nonlinear variation characteristic, after yielding, hysteresis loops are seen on the stress curves of all the unloadingreloading cycles. The hysteresis loops indicate that the strain of the unloading points contains the viscous strain [25]. The viscous strain (ε v ) of polymeric materials, calculated by subtracting the elastic strain (ε e ) and plastic strain (ε p ) from the strain of the unloading point (as shown in Figure 9), originates from the chain dynamics that are determined by the friction between polymer atoms, driven by entropic effects. When the polymer is subjected to relatively small or moderate strain, the mobility of tightly bonded molecules is constrained by the van der Waals interactions that cause friction between polymer atoms. When the deformation increased, the molecules are gradually separated in the space, and the accumulated plastic deformation causes some damage which weakens the van der Waals interactions between the molecules, thus resulting in viscous strain of the polymers. From Figure 9, it can also be read that when the strain of the unloading point increases the viscous strain in PU matrix increases gradually. Figure 10(a) shows the variations of the stress of the coarse-grained PU matrix with the running time at different strains of relaxation, from which one can see that when the stress relaxation begins, the axial stress declines sharply within a small time range, but not to zero. When the strain of relaxation is equal to 0.05 or 0.10, the axial stress declines linearly to a certain value, but when the strain of relaxation is equal to or higher than 0.20, the axial stress declines in a nonlinear way. Besides, when the strain of relaxation is lower than 0.20, after relaxation, the final values of the axial stress increase with the increasing strain of relaxation. For example, when the stress relaxation happened at ε T of 0.20, after relaxation, the axial stress of the model was about 10 MPa, about 33% higher than that relaxed at ε T of 0.10. However, when the strain of relaxation was equal to or higher than 8 Geofluids 0.40, after stress relaxation, all the axial stresses were almost equal to 15 MPa. The stress variation curves of the PU model during the stress relaxation process are similar to the experimental results of the macroscopic polyurethane [52,53], indicating the viscoelastic behavior of the material. Figure 10(b) presents the evolutions of the potential energy of the model during tension and subsequent stress relaxation (the strain of stress relaxation was 0.40). From Figure 10(b), it can be seen that during tension, the variations of the potential energy associated with the van der Waals interactions, bond lengths, and bond angles are similar to those shown in Figure 6. In the stress relaxation process, the PE angle term was maintained at about -2600 kcal/ mol, but the PE nonbond declined gradually from about 14000 kcal/mol to 12000 kcal/mol, which lasted for around 8 ps. The PE nonbond was maintained at 12000 kcal/mol in the following time. Besides, in the stress relaxation process, the variation of PE total was equal to that of PE nonbond . The variation of PE nonbond and PE total corresponds well with that of the axial stress presented in Figure 10(b), indicating that in the stress relaxation, the variation of axial stress of the PU model was mainly dominated by the varied van der Waals interactions. According to Figure 10(a), in the stress relaxation, E t , the relaxation modulus of the coarse-grained PU model, defined as the stress corresponding to unit strain, can be calculated. E t is a parameter of time dependence. Considering that the general Maxwell model illustrated in Figure 11(a) can be used to describe the stress relaxation property of viscoelastic solid materials [53,54], E t can be expressed as the following function:

Mechanical Behavior under the Stress Relaxation.
where E ∞ represents the average of E t when the relaxation time (t) reaches infinity, E i and λ i represent the real-time  In Equation (6), E ∞ can be calculated using the following function: where E 0 represents the relaxation modulus at the initial time. Therefore, Equation (6) can be rewritten as the following: To analyze the relaxation modulus, E t was transformed to the dimensionless form by dividing Equation (8) by E 0 , and the nondimension relaxation modulus can be expressed as the following: The variations of e t of the model in the stress relaxation at different ε T are plotted in Figure 11(b) as the scatter points. The time at the beginning of stress relaxation in Figure 10(a) has been shifted to zero. The variations of e t are fitted with Equation (9), and it was found that the best fit can be obtained for the model at different ε T when N is equal to 1. Therefore, the viscoelastic model for PU matrix during the stress relaxation process is described by the  Figure 9: The stress-strain curve of the coarse-grained PU matrix under successive cyclic tensile loading (the colored one). The stress-strain curve under tensile was also drawn as the black one for comparison. The strain rate is 5 × 10 10 /s.  Geofluids following: The parameters of the viscoelastic model and the goodness of fit (R 2 ) are listed in Table 4, and the theoretical values of e t calculated by Equation (10) have been drawn in Figure 11(b). It can be seen that even though a difference exists, the theoretical values of e t are close to the simulation results. The theoretical e t declines gradually with the increasing relaxation time but not to zero. With the increasing strain of relaxation, the decline of theoretical e t becomes more nonlinear. The theoretical stress value can be obtained by multiplying the theoretical e t with ε T . Since ε T is constant, it can be deduced that the theoretical stress will vary according to the same law as the theoretical e t . Therefore, the viscoelastic model consisting of one elastic element in parallel with one Maxwell element can be used to describe the mechanical behavior of an amorphous polyurethane material under stress relaxation.
In practical engineering, polymer grouting material refers to the product of the reaction between isocyanate and polyols with a constant mass ratio [5,9]. Therefore, in this study, only the influences of different loading types were investigated, and the influences of different raw materials or mass ratios on the mechanical properties of PU matrix will be studied in the further study.

Conclusions
The mechanical behavior of amorphous PU matrix under uniaxial tension, successive cyclic tension, and stress relaxation at nanoscale was investigated with coarse-grained MD simulation method. The conclusions are drawn as follows: (1) The mechanical properties of PU matrix under tensile are little affected by the number of molecules but change significantly when the strain rate varies. With the increasing strain rate, the yield strength increases, and the yield strains decrease. At yield, the total change in the potential energy is independent of the strain rate, and it is about 6000 kcal/ mol. When the strain rate increases, the declined volume strain of the PU matrix is the main reason for the increased yield strength (2) Under tensile, the potential energy associated with the bond stretching varies less significantly than that associated with the bond angle bending or van der Waals interactions, indicating that the deformation of the model is mainly caused by the slippage between different molecule chains and the bending of the bond angles. The bending of bond angles mainly occurs in the soft segments of the molecule chains. After yielding, most of the bond angles are still not in equilibrium (3) The contour of the stress-strain curve of the model under cyclic tension is consistent with that under uniaxial tension, indicating that the model performs with the stress memory behavior, and thus, its tensile  strength is independent of loading path. After yielding, the stress curves of the unloading-reloading cycles contain the hysteresis loops, which imply that the strain of the unloading points contains elastic strain, plastic strain, and viscous strain (4) During the stress relaxation process, the evolution of the axial stress is primarily caused by the varied van der Waals interactions. The mechanical behavior of the amorphous PU matrix under stress relaxation at the nanoscale can be described by the viscoelastic model consisting of one elastic element in parallel with one Maxwell element

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.