DPD method shares the common features of coarse-grained MD as the larger length scale and time scale by predigesting the exquisite atomic interactive force in classical MD. A softer interactive force between particles is adopted in DPD simulations and divided into three parts as (Groot and Warren, 1997)

where r_ij is distance between two beads. The FijC is the conservative force. It acts as the repulsive force between two particles within the cutoff R_c, which is linearly related to the particle distance with the repulsive coefficient a_ij. The FijD is the dissipative force which always points to the opposite direction of the relative velocity v⇀ij along the center line of two particles. It represents the effect of viscosity, depending on a coefficient γ and the distance between particles. The FijR is the random force. It also depends on a coefficient σ=2kTγ and the distance between two particles as well as a random variable α with the standard normal distribution. Since the dissipative force and random force act as the energy sink and source, the balance between two forces as a thermostat can be ensured by σ = 3, following the Groot's (Groot and Warren, 1997; Groot and Rabone, 2001) work.

Besides, it is convenient to normalize the mass of particles, the energy and the length of interactive cutoff as non-dimensional m = kT = R_c = 1, where k is the Boltzmann constant and T the temperature. Thus, the Newton's second law governing the particles motion is expressed by d2ri⇀dt2=F⇀i. Meanwhile, the unit temperature T = 1 in normalized DPD refers to absolute temperature 300 K. The time unit in DPD simulations is also normalized as τ=mRc2/kT= 1.

An additional bonding force acting as the bonds between polymer beads is described as

where C is the bond coefficient. A soft bond coefficient 4.0 (Rao et al., 2014; Wei et al., 2017) is used in the crosslinking simulations, while C = 116, 000 is adopted in tension simulations to match the real ratio of C-C bond strength to the heat fluctuation. r₀ is the equilibrium bond length between two polymer beads. We take the mean distance between two beads r₀ = R_c = 1.

Three types of molecules are coarse-grained in our study, i.e., acrylamide(AAm), methylenebisacrylamide (MBAA) and water as shown in Figure 1. Since the mass of all particles in DPD simulations are considered to be the same, one AAm bead refers to two AAm monomers, one MBAA bead refers to one MBAA molecules, and one water bead refers to eight water molecules, to make sure their relative molecular weights are close. The validity of a coarse-grained model is determined by the interactive parameters in force field. According to Groot and Warren (1997) and Liu M. B. et al. (2015), the repulsive coefficient a for the same type of beads can be approximated by

where ρ is the number of beads in the volume of Rc3. Previous researchers chose ρ = 3 to get a better match on the compressibility of the model fluid (Liu M. B. et al., 2015). However, this high number density of beads leads to pure repulsion, because the mean distance between beads is much smaller than the repulsion cutoff R_c in conservative force F^C. It makes the model solution under very high hydrostatic pressure when the number density is higher than 2 for the closest packing beads. Here we choose ρ=2 to simulate the crosslinking process only affected by the viscosity and heat fluctuations, where the repulsion force between closest packing beads just vanishes. For different types of beads, Groot (Groot and Warren, 1997; Groot and Rabone, 2001) suggest the repulsive coefficient a_ij is linearly related to the Flory-Huggins parameters χ_ij with a_ij = a_ii + 3.27kTχ_ij, where the χ_ij characterizes the mixture energy needed to form an equilibrium interaction between two clusters. We assume the monomer beads and crosslinker beads are the same type of beads denoted as polymer beads, and we take the repulsive coefficient χ_ij = 0.57 between polymer beads with water beads (Wei et al., 2017).

The real time scale of a DPD time unit can be estimated by comparing the diffusivity of water beads with the real diffusivity of water as follows (Groot and Rabone, 2001)

where the N_m is the coarse-grained scale representing the water molecule number that a water bead is, D_water the real diffusivity of water 2.43 ×10⁻⁹ m²/s, R_cr the real length of the DPD cutoff which can be obtained by the real volume of one water beads as 6.98 Å. The simulated diffusivity of water D_simu can be obtained from three independent simulations and calculated by Dsimu=〈r2〉6t=0.212±0.002 as shown in Figure 2, where 〈r²〉 is the mean square displacement of all water beads. Therefore, the real time for DPD unit time is about 340.0 ps.

As mentioned in introduction, most current models of polymer network are constructed based on the theoretical hypothesis, such as full crosslinking, crosslinking with certain orientation or certain degree of polymerization etc. However, the real polymer networks in PAAm hydrogels are mostly imperfect. In order to reveal the real conformation of the polymer network, we simulate the crosslinking process of the precursor solution. The crosslinking process of PAAm hydrogel in experiments often lasts hours because it is a growth of polymer chains guided by initiator (Sun and Marshall, 1981), such as tetramethylethylenediamine (TEMED) (Bai et al., 2017; Tang et al., 2017; Zhang et al., 2018; Lei et al., 2019). Although the time scale of DPD simulation is increased by coarse-grain method, it still cannot affords the full simulation of real crosslinking process to capture both the growth of polymer chains and the diffusion of monomers and crosslinkers in precursor solution. However, when the amount of initiator is close to the magnitude of precursor, the growth of polymer chains starts from everywhere in precursor solution, so that it can be simply regarded as the simultaneously bonding process as long as reactant molecules are close enough. Moreover, if we consider unreacted molecules to be identical in the precursor solution, it would not make much difference on the crosslinking polymer network when the diffusion length of unreacted molecules is much larger than their mean distance within our simulation time.

Based on the assumptions above, we conduct the following crosslinking simulations. Random mixture models with water content 80% are built and crosslinked under different temperatures. Temperature in DPD simulations determines the diffusion length of beads. Low temperature with insufficient diffusion length leads to localized polymerization in which it's hard to form an effective polymer network throughout hydrogel. In order to find the proper temperature for crosslinking in DPD simulations, we choose six temperatures, i.e., T = 1, 3, 5, 7, 9, 11, to investigate the effect of diffusion length on the crosslinking process. Meanwhile, the effect of different water contents on the crosslinking process is investigated by simulating the crosslinking process of random mixture models with different water content 99, 98, 97, 96, and 95%. These high water contents are chosen to find the lowest crosslinking threshold of precursor content.

All random mixture models are composed of 125,000 beads with the certain precursor mass ratio AAm:MBAA = 1:0.002. For example, models with water content 80% have 100,000 water beads, 24,954 AAm beads and 46 MBAA beads. The simulation box size is (44.54Rc)3. NVT ensemble is adopted as the thermostat. Time step is set as 0.01 DPD time. The mixture models are relaxed for 10,000 time steps first. Then the crosslinking process is performed by creating bonds between polymer beads every 10 time steps within 500,000 time steps, when the distance between two polymer beads is smaller than one DPD cutoff R_c. The crosslinked models are then cooled down to T = 1 within 100,000 time steps and relaxed for another 100,000 time steps. For each water content and temperature, three independent models are generated for statistical averaging. Figure 3A is one of the models and Figure 3B show its polymer network without water beads, where blue beads are AAm beads, red beads are MBAA beads and green beads are water beads. All the following figures about the model structure are shown without water beads for better view.

Polymer chains in all models are classified into different chain types for structural analysis. Incompressible uniaxial tension tests are simulated in NVT ensemble with T = 1 in three loading rates 0.005, 0.001, and 0.0002 per DPD time unit. It should be noted that the loading rates here are true strain rate for the convenience of the implementation of incompressible uniaxial tension in LAMMPS. Thus, the true strain rates in other two directions are set to be half of the loading rate. Time step for tension tests is set as 0.001 DPD time.

The large deformation of PAAm hydrogel is dominated by the polymer network conformation change in mesoscale, while the fracture of polymer chains is determined by the bond break in atomic scale. In order to obtain the fracture criterion of polymer chains in PAAm hydrogel, we build a full-atom model for PAAm chains with two connected AAm monomers as shown in the inset in Figure 4. This model corresponds to one AAm bead in our DPD simulations. Consistent valence force field (CVFF) (Dauber-Osguthorpe et al., 1988) is adopted to characterize the electrostatic forces, van der Waals interactions, bonds, bond angles, dihedrals and impropers between atoms. Detailed force field parameters can be found in Supplementary Material. Tension is imposed along chain direction in NVT ensemble with the loading rate 0.01/ps and temperature T = 300K. The time step is set as 0.1 fs. The energy-strain curve in Figure 4 shows the fracture of AAm chains occurs when the engineering strain is 0.225. This strain is chosen to be the fracture criterion we use to analyze the fracture properties of PAAm hydrogel models in DPD simulations.

In order to find a proper temperature for crosslinking simulations, six temperatures, i.e., T = 1, 3, 5, 7, 9, 11, is used to control the diffusion length of precursor beads in the crosslinking simulations. Figures 5A–C show part of the PAAm hydrogel models after crosslinking under temperatures T = 1, 5, 9, respectively. To analyze the generated complex structure of PAAm hydrogel models, polymer chains in all models are classified into five types, i.e., effective network chains, branch chains, isolated chains, crosslinking loops and isolated loops, as the schematic diagram shown in Figure 5D. Effective network chains form the crosslinking network throughout the hydrogel model. Branch chains attach to the effective network. It may be a chain with one free end or a loop. Isolated chains, crosslinking loops and isolated loops have no covalent bonds with effective network. The number and length of each type of chains are counted and averaged between the same groups of models. The crosslinking rate can be obtained by N_bond/(N_AAm + 2N_MBAA), where N_bond is the total bond number, N_AAm the number of AAm beads, N_MBAA the MBAA beads, and N_AAm + 2N_MBAA the maximum bonds between these monomers and crosslinkers. The final crosslinking rate for all models is above 99.97%. This indicates the crosslinking process in every model is sufficient. Chain numbers and number fractions for different type of chains are counted for all models. The chain lengths for different type of chains are counted as the bond number in current chain. The chain length fraction is calculated by the current chain length over total chain length.

The diffusion length of polymer beads is the key to the formation of polymer network. Figure 6A shows the diffusion length per DPD time unit in different temperatures. When T = 1, the diffusion length per DPD time unit is 1.13R_c. This is even lower than the mean distance between polymer beads 1.52R_c marked as red dashed line in Figure 6A. It makes the precursor solution fail to form polymer network as shown in Figure 5A. Polymer beads crosslinking with neighbors form large amount of isolated loops with short length, corresponding to the largest chain number as shown in Figure 6B and the highest length fraction of isolated loops. With the increase of temperature from 3 to 11, the diffusion length per DPD time unit is much larger than the mean distance of polymer beads as shown in Figure 6A. Polymer network forms in crosslinking process. Both the number fraction and the length fraction of effective network chains shown in Figures 6B, 7A rise up with the temperature increasing. The length fraction of effective network chains reach to the top when T = 9. Meanwhile, the total chain number and the length fraction of effective chains and isolated loops tends to converge when temperature increases. This convinces us that our crosslinking simulations are sufficient in high temperatures and our PAAm hydrogel models are reliable. We also note that in Figure 7A the length fraction of isolated loops converges to about 7% with the temperature increasing. It indicates that the larger diffusion length, which may results from higher temperature, longer gelation time or more sufficient stirring in experiments, cannot eliminate the formation of isolated loops. The success of forming polymer network proves that our crosslinking simulation is a practical way to construct PAAm hydrogel model.

Further validation of our crosslinking method are conducted by predicting the crosslinking limit of the polymer content. Random mixture models with five water content 99, 98, 97, 96, and 95% are crosslinking with the temperature T = 9. Figure 7B shows the effective chain ratio of models with different water contents. It can be seen that models with water content 99% fail to form polymer network. Meanwhile, one of the three models with water content 97% also fails to form polymer network as well as one of the three models with water content 98%. These simulations give the crosslinking limit of water content 97%, which is very close to the swelling limit of many PAAm hydrogels in experiments (Zhang et al., 2018). The prediction of crosslinking limit also proves the validity of our crosslinking simulations and PAAm hydrogel models.

We choose three PAAm hydrogel models with water content 80% crosslinking at T = 9 for further discussion. Figure 8A shows the number fraction and the length fraction of all types of chains in these models. Error bars for both the number fraction and length fraction indicate that our crosslinking process is reproducible. We denote all chains except effective network chains as ineffective chains. Only 7.1 ± 0.9% of polymer beads form ineffective chains, including branch chains 0.9 ± 0.4% and isolated loops 6.2 ± 0.6%. Because the mean contour length of ineffective chains 8.8 ± 0.7r₀ is much lower than the effective network chains 288.7 ± 10.5r₀, the length fraction of effective chains is 92.9 ± 0.9% though the number fraction is 28.5 ± 1.7%. Such high mean contour length leads to compliant effective network with the chain conformation more winding than other full-atom models (Tönsing and Oldiges, 2001; Oldiges and Tönsing, 2002; Wu et al., 2009; Mathesan et al., 2016; Xu et al., 2016; An et al., 2019; Hou et al., 2019) and coarse-grain models (Nawaz and Carbone, 2014; Wei et al., 2017; Xing et al., 2019) of PAAm hydrogels. This mean contour length is very close to that estimated by N_AAm:2N_MBAA=271 from precursor ratio. The mean end-to-end distance of effective network chains is 47.0 ± 3.6r₀. The distributions of the initial length ratio which is the end-to-end distance over the contour length of each effective network chain in three models are shown in Figure 8B. It can be seen that the initial length ratio of effective chains shows nearly a Maxwell distribution. In order to validate the statistical properties of our models, we also present the distribution of two orientation angles (θ, φ in Figure 9A) and the ratio of the end-to-end distance to the contour length of all effective network chains in Figure 9B. The distribution of the chain orientation angle θ is almost uniform and the distribution of the chain orientation angle φ is almost sine-shaped. It proves that our model is sufficiently large to bridge the molecular and continuum properties of PAAm hydrogel.

Water beads are uniformly distributed around polymer chains. Water molecules in hydrogels are trapped by hydrophilic functional groups in polymer network via hydrogen bonds. The coverage of these hydrogen bonds is only within several water molecules away, meaning that the so-called bound water in hydrogels should be distributed uniformly surrounding polymer network and other chains. The hydrophilia of the polymer beads is embodied in our simulation as the relatively small Flory-Huggins parameter χ. It can be seen that our PAAm hydrogel models give reasonable crosslinking threshold and show sufficiently realistic and large polymer network in PAAm hydrogel. It convinces us to use these models to investigate the structural and mechanical properties of PAAm hydrogel.

In this section, Uniaxial tension tests with constant volume with λ1=λ,λ2=λ3=λ-12, where λ_i is the stretch in i-th direction, are performed in DPD simulations to investigate the large deformation behavior of PAAm hydrogels with water content 80%. The nominal stress of the incompressible uniaxial tension test is obtained from

where σ_i is the Virial stress in each direction.

Three loading rates, i.e., 0.005, 0.001, and 0.0002 per DPD time unit, are imposed on PAAm hydrogel models. Figure 10A shows the nominal stress-stretch curves of three PAAm hydrogel models under three loading rates. It shows that our PAAm hydrogel models capture both the viscoelastic and the hyperelastic behaviors. All three models shows the loading rate-dependent behavior. This is because the viscosity of all beads is embedded in DPD force field F^D. Higher loading rate causes more intensive rebound of bonded beads, leading to the higher dissipative force. Compared to previous experimental results (Lei et al., 2019; Yang et al., 2019), the loading rate-dependent behavior is much more significant since the loading rates we adopt in our DPD simulations are far larger than realistic ones used in experiments. Considering the loading rate-dependent behavior and the computational resources, we choose loading rate 0.001 per DPD time unit for further discussion. The hyperelasticity is shown in the initial part of the stress-stretch curve during uniaxial tension. The hardening stage can also be found where the stretch is large, which is cause by the stretch limit of polymer chains. We also present the bond stretch in effective chains during the uniaxial tension tests as the solid lines shown in Figure 10B. It shows that significant bond stretch occurs in the hardening stage where effective chains are stretched to almost straight. The trend of the nominal stress-stretch curves agrees well with experimental results (Bai et al., 2017; Zhang et al., 2018; Yang et al., 2019). The scattered stretch limit where hardening stage occurs in PAAm hydrogels with the same water content are also found.

Although the number fraction and length fraction of each type of chains in three PAAm hydrogel models only have a <2% standard deviation, their mechanical behaviors in Figure 10A show much differences, especially when the stretch is large. This is not what the current theory of hydrogels can predict. Therefore, we compare our simulation results with current constitutive theory of polymers to reveal the underlying mechanism of the large deformation behavior of PAAm hydrogel. In current constitutive theory (Huang et al., 2020), two key parts, i.e., the free energy of a single chain and the statistical sum of the free energy of all chain, are crucial to bridge the mesoscale chain conformation change to the bulk mechanical response of hydrogels. Considering the polymer chains in our DPD simulations are close to the so-call freely joint chains, the Langevin chain model is adopted to characterize the free energy of a single chain

where kT is the thermodynamics energy unit, N the number of Kuhn segments in current chain, λ_c the stretch of chain with respect to the mean initial length ratio of current chain l0=L0Nbb0, L₀ the mean end-to-end distance of current chain, N_br₀ the mean contour length with the mean bond number N_b and the bond length r₀. Since the inverse Langevin function L-1 is singular when λ_cl₀ = 1, the Langevin model often need additional modification (Mao et al., 2017). Li and Bouklas (2020) proposed the stretching force of a polymer chain caused by the conformation entropy change can also stretch bonds in polymer chains. Thus, the modified chain stretch λ_c turns out to be λ_c/λ_b, where λ_b is the mean stretch of bonds in effective chains. Combining with the bond force F^B in section Governing Equation of DPD Simulations, this modification circumvents the singularity of the inverse Langevin function by giving the relationship of chain stretch λ_c to the mean bond stretch λ_b as follows

In order to statistically sum up the free energy of all effective chains, the mean stretch of all effective chains 〈λc〉=I1/3 and mean number of Kuhn segments 〈N〉 are used to formulate the total free energy of hydrogel model as follows

where n is the effective chain density. I₁ is the first invariant of the deformation gradient of bulk hydrogel model with I1=λ12+λ22+λ32=λ2+2λ-1, which connect the chain conformation change in mesoscale to the bulk deformation in continuum. Thus, the λ_b − λ relationship between mean bond stretch and the stretch of PAAm hydrogel models can be obtained theoretically by solving following equation

Figure 10B compare the theoretically λ_b − λ relationship with that measured in simulations. It is very clear that there are large discrepancies between theoretical predictions and simulation results. Parameter study for the λ_b − λ relation is performed to find the reason for the large discrepancies. Figures 11A,B shows the λ_b − λ relationship with the different bond strength parameter CkT and different mean initial length ratio l₀, respectively. Combining Figures 10B, 11A, we can find that the simulation results show a much softer mean bond stretching process during uniaxial tension of PAAm hydrogel models than theoretical predictions, although such a large bond strength parameter CkT= 116,000 is used in both simulations and theory. From lines in Figures 12A–C, we can find that the softer mean bond stretch of the whole model is moderated by the unsynchronized stretch of different effective chains. Also, the flat stress-stretch curves for hydrogel actually result from the superposition of the unsynchronized stiff bond stretching behavior of large amount of chains. In Figure 11B, we can find that the critical stretch where significant bond stretch occurs depends on the initial length ratio. The discrepancy of critical stretch between theory and simulation is because of the use of mean initial length ratio. As shown in Figure 8B, the initial length ratio of effective chains distributes in a very wide range. Chains with high initial length ratio first stretch to almost straight, so that using the lower mean initial length ratio overestimates the critical stretch in Figure 11B.

The nominal stress-stretch curves and the λ_b − λ relation indicates that current continuum constitutive theory for hydrogel is only valid in small stretch, while the mechanical response under large deformation near fracture needs detailed structural models. The statistical averaging structural properties, such as water content and mean chain length, smears out the wide distribution of the possible conformation of effective chains, while specific structural features are what dominate the large deformation behavior of hydrogel.

The randomness of the polymer network in hydrogel makes it hard to give an accurate fracture criterion. The widely used fracture toughness works poor as the fracture criterion of hydrogels since it often counts the deformation energy far away from the crack. As shown in the PAAm hydrogel models, the stretch limit is the key property that determines the fracture of a polymer chain. A stretch criterion is more general which can be applied in different polymer network systems and more practical for different loading conditions. In this section, we propose a stretch criterion of the fracture of hydrogels.

The critical energy to the fracture of a polymer chains is equal to the fracture energy of one C-C bond

where λbcr is the critical stretch when a C-C bond breaks. This critical free energy corresponds to the critical mean stretch of all bonds λ¯bcr in this chain with

Thus, the fracture criterion for a single chain can be expressed as

Using the λ¯bcr=1.225 when N_b = 2 obtained in Figure 4, the critical stretch of a C-C bond λbcr is 1.318. The critical mean stretch of a C-C bond λbcr is then related to the first invariant of the bulk deformation gradient by substituting the critical mean stretch λ¯bcr into the λ_b − λ relation. As the dots shown in Figure 10B, for three models with the N_b = 600.0, 558.6, 573.6, the theoretical fracture stretch are 11.2, 11.1 and 9.8, and the fracture stretch measured from simulations are 9.1, 8.8, and 7.4, respectively. There are much discrepancies between the theoretical results and simulation results, mainly because the equal-chain-length assumption in current constitutive theory overestimates the mean chain stretch by 〈λc〉=I1/3. Furthermore, the fracture stretch measured from simulations using the mean bond stretch of all effective chains should be the upper limit of the real fracture stretch. Considering the unsynchronized stretch of different effective chains, 50% of effective chains may have been broken when the mean bond stretch of all effective chains reaches fracture stretch. This amount of chain scission is enough for the macroscopic fracture of hydrogel. On the other hand, we can get the lower limit of the fracture stretch by analyzing the fracture stretch of every effective chains. Blue dots in Figures 12A–C are the fracture stretch where the first effective chain (as shown in Figure 12D) is about to break, while red dots are fracture stretch obtained from mean bond stretch of all effective chains. It gives limited ranges of the fracture stretch of our three PAAm hydrogel models as 5.6–9.1, 8.3–8.8, and 5.0–7.4, respectively. These ranges almost cover the fracture stretch of PAAm hydrogels with the water content 78% in previous experiments (Zhang et al., 2018) as 5.5, 8.6 and 6.1.

Our coarse-grain PAAm hydrogel models are the representative models for numerous hydrogels since only the repulsive coefficient a is adjustable for different hydrogel system. Hydrogels share the same large deformation mechanism with the conformation change of the complex polymer network, so that the detailed mesoscale model is such a powerful tool to bridge the molecular movement to bulk deformation.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.