Abstract
Meso-scale models for hydrogels are crucial to bridge the conformation change of polymer chains in micro-scale to the bulk deformation of hydrogel in macro-scale. In this study, we construct coarse-grain bead-spring models for polyacrylamide (PAAm) hydrogel and investigate the large deformation and fracture behavior by using Dissipative Particle Dynamics (DPD) to simulate the crosslinking process. The crosslinking simulations show that sufficiently large diffusion length of polymer beads is necessary for the formation of effective polymer. The constructed models show the reproducible realistic structure of PAAm hydrogel network, predict the reasonable crosslinking limit of water content and prove to be sufficiently large for statistical averaging. Incompressible uniaxial tension tests are performed in three different loading rates. From the nominal stress-stretch curves, it demonstrated that both the hyperelasticity and the viscoelasticity in our PAAm hydrogel models are reflected. The scattered large deformation behaviors of three PAAm hydrogel models with the same water content indicate that the mesoscale conformation of polymer network dominates the mechanical behavior in large stretch. This is because the effective chains with different initial length ratio stretch to straight at different time. We further propose a stretch criterion to measure the fracture stretch of PAAm hydrogel using the fracture stretch of C-C bonds. Using the stretch criterion, specific upper and lower limits of the fracture stretch are given for each PAAm hydrogel model. These ranges of fracture stretch agree quite well with experimental results. The study shows that our coarse-grain PAAm hydrogel models can be applied to numerous single network hydrogel systems.
Introduction
A hydrogel is a network of polymer chains swollen in water. Synthetic hydrogels have developed rapidly since the landmark research by Wichterle and Lím (1960). Due to the large water content, hydrogels can be bio-compatible, highly compliant and exhibit low friction, making them widely used for personal care and medical applications, such as superabsorbent diapers (Masuda, 1994), contact lenses (Caló and Khutoryanskiy, 2015), drug delivery (Li and Mooney, 2016; Liu et al., 2017), tissue engineering (Lee and Mooney, 2001; Haque et al., 2012; Lienemann et al., 2012; Xu et al., 2018), and wound dressing (Li et al., 2017). However, the large water content also makes hydrogels have very low resistance to deformation and fracture (Zheng et al., 2018; Xu and Liu, 2019) and hard to be used as load-bearing structures. Since Gong's (Gong et al., 2003) work, diverse efforts have been made to design tough hydrogels by building dissipations into the networks (Zhao, 2014). Examples include double-network hydrogels (Gong et al., 2003; Henderson et al., 2010), poly(vinyl alcohol) hydrogels with crystalline domains (Peppas and Merrill, 1976; Stauffer and Peppast, 1992), hydrogels with hybrid physical and chemical crosslinkers (Kong et al., 2003; Shull, 2012; Sun et al., 2012), and hydrogels with transformable domains (Brown et al., 2009). It is well-known that the hyper-elastic nature of hydrogel originates from its crosslinking polymer network (Liu Z. et al., 2015). Yet the crosslinking polymer network in hydrogels still stays in the realm of hypothesis since the dynamic experimental observations and determinations are inconvenient at such a micro level. In order to reveal the true nature of the crosslinking polymer network during the deformation and fracture of hydrogels, we have to bring ourselves down to the mesoscale or even molecular scale, where molecular dynamics (MD) simulations and Monte-Carlo simulations could be an effective approach.
Researchers have proposed all-atom models (Tönsing and Oldiges, 2001; Oldiges and Tönsing, 2002; Wu et al., 2009; Mathesan et al., 2016; Xu et al., 2016; Hou et al., 2019) to investigate the structural and physical properties of hydrogels, such as the hydrogen-bond configuration and thermal conductivity. These models usually contain only several short polymer chains whose lengths are in the same scale of the persistent length of polymer chains. This makes them hard to depict the compliant polymer network in hydrogels. Besides, most researches focus on the equilibrium properties of hydrogels without bearing loads. An et al. (2019) investigated the mechanical properties of hydrogels via MD simulations, but the strain is small. Jang et al. (2007) obtained the stress-strain curve under large strain, but the stress level is several orders of magnitude larger (as well as An et al.'s work) than real situations. Although many general force fields, including consistent valence force field (CVFF) (Dauber-Osguthorpe et al., 1988), optimized potentials for liquid simulations (OPLS) force field (Kaminski et al., 1994), GROMACS force field (Berendsen et al., 1995), DREIDING force field (Mayo et al., 1990), etc., have been proposed and can be applied to the polymer-solvent system, the sophisticated atomic models with too many structural details hinders the extension of length scale and time scale. Since the large deformation behavior of hydrogels originates from the conformation change of the polymer network in mesoscale, a coarse-grain model of hydrogels balancing the structural features and computational time is in great need.
Dissipative particle dynamics (DPD) (Hoogerbrugge and Koelman, 1992; Koelman and Hoogerbrugge, 1993) is a mesoscale particle method that bridges the gap between macroscopic and microscopic simulations. Combining with coarse-grain methods, it is very suitable for simulating gaseous or fluid systems. DPD has been successfully applied to diverse areas of interests, especially to simulate the equilibrium and dynamical properties of polymers in solution as well as polymer gels (Spenley, 2000; Maiti and McGrother, 2004; Symeonidis et al., 2005; Zhao, 2014). When applying DPD, or other coarse grained MD simulations to study polymers or polymer gels, a physical crosslinking process is crucial to construct real polymer network. However, many previous works (Wu et al., 2009; Nawaz and Carbone, 2014; Mathesan et al., 2016; Jin et al., 2018; Hou et al., 2019; Xing et al., 2019) construct polymer models based on hypothetical cross-linked network structures which inevitably loses some structural components, such as branch chains, polymer loops, unreacted monomers, short segments et al. For instance, Xing et al. (2019) constructed 3D cross-linked networks for DNA hydrogels, while the chain length between cross-linking points was set to be constant. Jin et al. (2018) built randomly cross-linked polymer networks with the real cross-linking densities, while all the polymer chains are forced to crosslink.
In this study, the polyacrylamide hydrogel is taken to be the representative material. We build the bead-spring models of PAAm hydrogel by using DPD to simulate the crosslinking process from the mixture of monomers, cross-linkers and water. Details on the model construction are discussed in section Materials and Methods. The large deformation and fracture behavior of our PAAm hydrogel models are discussed in section Results.
Materials and Methods
In this section, a series of bead-spring models for PAAm hydrogel are constructed and tested using DPD simulations. In addition, a full-atom model for PAAm chain is also constructed and tested using classical MD simulations to provide the fracture criterion of the polymer chains in DPD simulations. All these simulations are performed in large-scale atomic/molecular massive parallel simulator (LAMMPS) (Plimpton, 1995).
Governing Equation of DPD Simulations
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 rij is distance between two beads. The is the conservative force. It acts as the repulsive force between two particles within the cutoff Rc, which is linearly related to the particle distance with the repulsive coefficient aij. The is the dissipative force which always points to the opposite direction of the relative velocity along the center line of two particles. It represents the effect of viscosity, depending on a coefficient γ and the distance between particles. The is the random force. It also depends on a coefficient 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 = Rc = 1, where k is the Boltzmann constant and T the temperature. Thus, the Newton's second law governing the particles motion is expressed by . 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 .
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. r0 is the equilibrium bond length between two polymer beads. We take the mean distance between two beads r0 = Rc = 1.
Coarse-Graining
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 . 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 Rc in conservative force FC. It makes the model solution under very high hydrostatic pressure when the number density is higher than for the closest packing beads. Here we choose 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 aij is linearly related to the Flory-Huggins parameters χij with aij = aii + 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).
Figure 1
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 Nm is the coarse-grained scale representing the water molecule number that a water bead is, Dwater the real diffusivity of water 2.43 ×10−9m2/s, Rcr 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 Dsimu can be obtained from three independent simulations and calculated by as shown in Figure 2, where 〈r2〉 is the mean square displacement of all water beads. Therefore, the real time for DPD unit time is about 340.0 ps.
Figure 2
Modeling and Testing
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 . 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 Rc. 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.
Figure 3
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.
Fracture Criterion of Polyacrylamide Chain
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.
Figure 4
Results
Model Validation
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 Nbond/(NAAm + 2NMBAA), where Nbond is the total bond number, NAAm the number of AAm beads, NMBAA the MBAA beads, and NAAm + 2NMBAA 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.
Figure 5
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.13Rc. This is even lower than the mean distance between polymer beads 1.52Rc 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.
Figure 6
Figure 7
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.7r0 is much lower than the effective network chains 288.7 ± 10.5r0, 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 NAAm:2NMBAA=271 from precursor ratio. The mean end-to-end distance of effective network chains is 47.0 ± 3.6r0. 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.
Figure 8
Figure 9
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.
Large Deformation Behavior
In this section, Uniaxial tension tests with constant volume with , 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 FD. 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.
Figure 10
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 the mean end-to-end distance of current chain, Nbr0 the mean contour length with the mean bond number Nb and the bond length r0. Since the inverse Langevin function is singular when λcl0 = 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 FB 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 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. I1 is the first invariant of the deformation gradient of bulk hydrogel model with , 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 and different mean initial length ratio l0, 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 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.
Figure 11
Figure 12
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.
Fracture Criterion
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 is the critical stretch when a C-C bond breaks. This critical free energy corresponds to the critical mean stretch of all bonds in this chain with
Thus, the fracture criterion for a single chain can be expressed as
Using the 1.225 when Nb = 2 obtained in Figure 4, the critical stretch of a C-C bond is 1.318. The critical mean stretch of a C-C bond is then related to the first invariant of the bulk deformation gradient by substituting the critical mean stretch into the λb − λ relation. As the dots shown in Figure 10B, for three models with the Nb = 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 . 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.
Conclusions
In this study, we propose a method to construct the mesoscale PAAm hydrogel models and investigate the large deformation and fracture mechanism of PAAm hydrogel using DPD simulations. The coarse-grain PAAm hydrogel models are constructed by simulating the crosslinking process in experiments. Different temperatures are tested for achieving sufficient diffusion length of polymer beads for sufficient crosslinking. It shows that the formation of polymer network in crosslinking process only occurs when the diffusion length of polymer beads is much larger than the mean distance of polymer beads. However, the increasing of diffusion length, which may caused by the increasing of temperature, gelation time or stirring in experiments, cannot eliminate the formation of isolated loops. Our PAAm hydrogel models have realistic structure of polymer network, including the compliant effective network, branch chains, isolated chains, crosslinking loops and isolated loops. The degree of polymerization of the effective network chains is almost the same to the theoretical estimation. Our crosslinking simulations also show the upper limit of the water content for forming polymer network is about 97% which is close to the swelling limit of PAAm hydrogel in experiments. The scale of our PAAm hydrogel models are proved to be sufficiently large by the uniformly distributed orientation of effective chains. Incompressible uniaxial tension tests are performed in three different loading rates using our PAAm hydrogel models with the water content 80%. The nominal stress-stretch curves reflect both the hyperelasticity and the viscoelasticity of our PAAm hydrogel models. However, the scattered large deformation behaviors of PAAm hydrogel models with the same water content indicate that the mesoscale conformation of polymer network have great impact on the mechanical behavior in large stretch, because the effective chains with a wide range of initial length ratio stretch to straight at different time during deformation. Furthermore, we propose a stretch criterion of the fracture of PAAm hydrogel using the fracture stretch of C-C bonds. By analyzing our PAAm hydrogel models, specific upper and lower limit of the fracture stretch are given for each PAAm hydrogel models. These ranges agree quite well with experimental results. Our coarse-grain PAAm hydrogel models can be applied to numerous single network hydrogel systems.
Statements
Data availability statement
All datasets generated for this study are included in the article/Supplementary Material.
Author contributions
JL, SX, and ZLi: investigation. ZLiu: resources, supervision, and funding acquisition. JL and SX: writing—original draft preparation. ZLi and ZLiu: writing—review and editing.
Funding
This work was supported by the National Natural Science Foundation of China through grant numbers 11820101001, 11811530287, and 11572236.
Conflict of interest
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.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2020.00115/full#supplementary-material
References
1
AnM.DemirB.WanX.MengH.YangN.WalshT. R. (2019). Predictions of thermo-mechanical properties of cross-linked polyacrylamide hydrogels using molecular simulations. Adv. Theory Simulations2:1800153. 10.1002/adts.201800153
2
BaiR.YangQ.TangJ.MorelleX. P.VlassakJ. J.SuoZ. (2017). Fatigue fracture of tough hydrogels. Extreme Mech. Lett. 15, 91–96. 10.1016/j.eml.2017.07.002
3
BerendsenH. J.van der SpoelD.van DrunenR. (1995). GROMACS: a message-passing parallel molecular dynamics implementation. Comput. Phys. Commun.91, 43–56. 10.1016/0010-4655(95)00042-E
4
BrownA. E.LitvinovR. I.DischerD. E.PurohitP. K.WeiselJ. W. (2009). Multiscale mechanics of fibrin polymer: gel stretching with protein unfolding and loss of water. Science325, 741–744. 10.1126/science.1172484
5
CalóE.KhutoryanskiyV. V. (2015). Biomedical applications of hydrogels: a review of patents and commercial products. Euro. J. Polym. 65, 252–267. 10.1016/j.eurpolymj.2014.11.024
6
Dauber-OsguthorpeP.RobertsV. A.OsguthorpeD. J.WolffJ.GenestM.HaglerA. T. (1988). Structure and energetics of ligand binding to proteins: Escherichia coli dihydrofolate reductase-trimethoprim, a drug-receptor system. Proteins Struct. Funct. Bioinformatics4, 31–47. 10.1002/prot.340040106
7
GongJ. P.KatsuyamaYKurokawaTOsadaY. (2003). Double-network hydrogels with extremely high mechanical strength. Adv Mater. 15, 1155–1158. 10.1002/adma.200304907
8
GrootR. D.RaboneK. (2001). Mesoscopic simulation of cell membrane damage, morphology change and rupture by nonionic surfactants. Biophys. J. 81, 725–736. 10.1016/S0006-3495(01)75737-2
9
GrootR. D.WarrenP. B. (1997). Dissipative particle dynamics: bridging the gap between atomistic and mesoscopic simulation. J. Chem. Phys.107, 4423–4435. 10.1063/1.474784
10
HaqueM. A.KurokawaT.GongJ. P. (2012). Super tough double network hydrogels and their application as biomaterials. Polymer53, 1805–1822. 10.1016/j.polymer.2012.03.013
11
HendersonK. J.ZhouT. C.OtimK. J.ShullK. R. (2010). Ionically cross-linked triblock copolymer hydrogels with high strength. Macromolecules43, 6193–6201. 10.1021/ma100963m
12
HoogerbruggeP.KoelmanJ. (1992). Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics. Europhys. Lett. 19:155. 10.1209/0295-5075/19/3/001
13
HouD.XuJ.ZhangY.SunG. (2019). Insights into the molecular structure and reinforcement mechanism of the hydrogel-cement nanocomposite: an experimental and molecular dynamics study. Compos. B Eng. 177:107421. 10.1016/j.compositesb.2019.107421
14
HuangR.ZhengS. J.LiuZ. S. (2020). Recent advances of the constitutive models of smart materials- hydrogels and shape memory polymers, Int. J. Appl. Mech. 12:2050014. 10.1142/S1758825120500143
15
JangS. S.GoddardW. A.KalaniM. Y. S. (2007). Mechanical and transport properties of the poly (ethylene oxide)–poly (acrylic acid) double network hydrogel from molecular dynamic simulations. J. Phys. Chem. B111, 1729–1737. 10.1021/jp0656330
16
JinK.BarreiroD. L.Martin-MartinezF. J.QinZ.HammM.PaulC. W.et al. (2018). Improving the performance of pressure sensitive adhesives by tuning the crosslinking density and locations. Polymer154, 164–171. 10.1016/j.polymer.2018.08.065
17
KaminskiG.DuffyE. M.MatsuiT.JorgensenW. L. (1994). Free energies of hydration and pure liquid properties of hydrocarbons from the OPLS all-atom model. J. Phys. Chem.98, 13077–13082. 10.1021/j100100a043
18
KoelmanJ.HoogerbruggeP. (1993). Dynamic simulations of hard-sphere suspensions under steady shear. Europhys. Lett. 21:363. 10.1209/0295-5075/21/3/018
19
KongH. J.WongE.MooneyD. J. (2003). Independent control of rigidity and toughness of polymeric hydrogels. Macromolecules36, 4582–4588. 10.1021/ma034137w
20
LeeK. Y.MooneyD. J. (2001). Hydrogels for tissue engineering. Chem. Rev.101, 1869–1880. 10.1021/cr000108x
21
LeiJ.ZhouZ.LiuZ. (2019). Side chains and the insufficient lubrication of water in polyacrylamide hydrogel—a new insight. Polymers11:1845. 10.3390/polym11111845
22
LiB.BouklasN. (2020). A variational phase-field model for brittle fracture in polydisperse elastomer networks. Int. J. Solids Struct.182, 193–204. 10.1016/j.ijsolstr.2019.08.012
23
LiJ.CelizA. D.YangJ.YangQ.WamalaI.WhyteW.et al. (2017). Tough adhesives for diverse wet surfaces. Science357, 378–381. 10.1126/science.aah6362
24
LiJ.MooneyD. J. (2016). Designing hydrogels for controlled drug delivery. Nat. Rev. Mater. 1:16071. 10.1038/natrevmats.2016.71
25
LienemannP. S.LutolfM. P.EhrbarM. (2012). Biomimetic hydrogels for controlled biomolecule delivery to augment bone regeneration. Adv. Drug Deliv. Rev.64, 1078–1089. 10.1016/j.addr.2012.03.010
26
LiuJ.PangY.ZhangS.ClevelandC.YinX.BoothL.et al. (2017). Triggerable tough hydrogels for gastric resident dosage forms. Nat. Commun.8:124. 10.1038/s41467-017-00144-z
27
LiuM. B.LiuG. R.ZhouL. W.ChangJ. Z. (2015). Dissipative particle dynamics (DPD): an overview and recent developments. Arch. Comput. Methods Eng. 22, 529–556. 10.1007/s11831-014-9124-x
28
LiuZ.TohW.NgT. Y. (2015). Advances in mechanics of soft materials: a review of large deformation behavior of hydrogels. Int. J. Appl. Mech. 07:1530001. 10.1142/S1758825115300011
29
MaitiA.McGrotherS. (2004). Bead–bead interaction parameters in dissipative particle dynamics: relation to bead-size, solubility parameter, surface tension. J. Chem. Phys.120, 1594–1601. 10.1063/1.1630294
30
MaoY.TalaminiB.AnandL. (2017). Rupture of polymers by chain scission. Extreme Mech. Lett. 13, 17–24. 10.1016/j.eml.2017.01.003
31
MasudaF. (1994). Trends in the development of superabsorbent polymers for diapers, in superabsorbent polymers. Am. Chem. Soc. 573, 88–98. 10.1021/bk-1994-0573.ch007
32
MathesanS.RathA.GhoshP. (2016). Molecular mechanisms in deformation of cross-linked hydrogel nanocomposite. Mater. Sci. Eng. C59, 157–167. 10.1016/j.msec.2015.09.087
33
MayoS. L.OlafsonB. D.GoddardW. A. (1990). DREIDING: a generic force field for molecular simulations. J. Phys. Chem.94, 8897–8909. 10.1021/j100389a010
34
NawazS.CarboneP. (2014). Coarse-graining poly(ethylene oxide)–poly(propylene oxide)–poly(ethylene oxide) (PEO–PPO–PEO) block copolymers using the MARTINI force field. J. Phys. Chem. B118, 1648–1659. 10.1021/jp4092249
35
OldigesC.TönsingT. (2002). Molecular dynamic simulation of structural, mobility effects between dilute aqueous CH 3 CN solution and crosslinked PAA Part 1. Struct. Phys. Chem. Chem. Phys. 4, 1628–1636. 10.1039/b110238a
36
PeppasN. A.MerrillE. W. (1976). Poly(vinyl alcohol) hydrogels: reinforcement of radiation-crosslinked networks by crystallization. J. Polym. Sci. Polym. Chem. Ed. 14, 441–457. 10.1002/pol.1976.170140215
37
PlimptonS. (1995). Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys.117, 1–19. 10.1006/jcph.1995.1039
38
RaoZ.HuoY.LiuX. (2014). Dissipative particle dynamics and experimental study of alkane-based nanoencapsulated phase change material for thermal energy storage. RSC Adv.4, 20797–20803. 10.1039/C4RA02699C
39
ShullK. R. (2012). Materials science: a hard concept in soft matter. Nature489, 36–38. 10.1038/489036a
40
SpenleyN. (2000). Scaling laws for polymers in dissipative particle dynamics. Europhys. Lett. 49:534. 10.1209/epl/i2000-00183-2
41
StaufferS. R.PeppastN. A. (1992). Poly (vinyl alcohol) hydrogels prepared by freezing-thawing cyclic processing. Polymer33, 3932–3936. 10.1016/0032-3861(92)90385-A
42
SunJ.-Y.ZhaoX.IlleperumaW. R.ChaudhuriO.OhK. H.MooneyD. J.et al. (2012). Highly stretchable and tough hydrogels. Nature489:133. 10.1038/nature11409
43
SunM.MarshallE. (1981). Gels. Sci. Am.244, 124–138. 10.1038/scientificamerican0181-124
44
SymeonidisV.KarniadakisG. E.CaswellB. (2005). Dissipative particle dynamics simulations of polymer chains: scaling laws and shearing response compared to DNA experiments. Phys. Rev. Lett. 95:076001. 10.1103/PhysRevLett.95.076001
45
TangJ.LiJ.VlassakJ. J.SuoZ. (2017). Fatigue fracture of hydrogels. Extreme Mech. Lett. 10, 24–31. 10.1016/j.eml.2016.09.010
46
TönsingT.OldigesC. (2001). Molecular dynamic simulation study on structure of water in crosslinked poly (N-isopropylacrylamide) hydrogels. Phys. Chem. Chem. Phys.3, 5542–5549. 10.1039/b109281m
47
WeiQ.WangY.ZhangY.ChenX. (2017). Aggregation behavior of nano-silica in polyvinyl alcohol/polyacrylamide hydrogels based on dissipative particle dynamics. Polymers9:611. 10.3390/polym9110611
48
WichterleO.LímD. (1960). Hydrophilic gels for biological use. Nature185, 117–118. 10.1038/185117a0
49
WuY.JosephS.AluruN. R. (2009). Effect of cross-linking on the diffusion of water, ions, and small molecules in hydrogels. J. Phys. Chem. B113, 3512–3520. 10.1021/jp808145x
50
XingZ.NessC.FrenkelD.EiserE. (2019). Structural and linear elastic properties of DNA hydrogels by coarse-grained simulation. Macromolecules52, 504–512. 10.1021/acs.macromol.8b01948
51
XuL.ZhaoX.XuC.KotovN. A. (2018). Water-rich biomimetic composites with abiotic self-organizing nanofiber network. Adv. Mater. 30:1703343. 10.1002/adma.201703343
52
XuS.LiuZ. (2019). A nonequilibrium thermodynamics approach to the transient properties of hydrogels. J. Mech. Phys. Solids127, 94–110. 10.1016/j.jmps.2019.03.008
53
XuS.WangY.HuJ.LiuZ. (2016). Atomic understanding of the swelling and phase transition of polyacrylamide hydrogel. Int. J. Appl. Mech.8:1640002. 10.1142/S1758825116400020
54
YangC.YinT.SuoZ. (2019). Polyacrylamide hydrogels. I. network imperfection. J. Mech. Phys. Solids131, 43–55. 10.1016/j.jmps.2019.06.018
55
ZhangE.BaiR.MorelleX. P.SuoZ. (2018). Fatigue fracture of nearly elastic hydrogels. Soft Matter14, 3563–3571. 10.1039/C8SM00460A
56
ZhaoX. (2014). Multi-scale multi-mechanism design of tough hydrogels: building dissipation into stretchy networks. Soft Matter10, 672–687. 10.1039/C3SM52272E
57
ZhengS.LiZ.LiuZ. (2018). The fast homogeneous diffusion of hydrogel under different stimuli. Int. J. Mech. Sci. 137, 263–270. 10.1016/j.ijmecsci.2018.01.029
Summary
Keywords
polyacrylamide hydrogel, dissipative particle dynamics, large deformation behavior, effective network, fracture criterion
Citation
Lei J, Xu S, Li Z and Liu Z (2020) Study on Large Deformation Behavior of Polyacrylamide Hydrogel Using Dissipative Particle Dynamics. Front. Chem. 8:115. doi: 10.3389/fchem.2020.00115
Received
04 December 2019
Accepted
07 February 2020
Published
25 February 2020
Volume
8 - 2020
Edited by
Kerstin G. Blank, Max Planck Institute of Colloids and Interfaces, Germany
Reviewed by
Akash Arora, Massachusetts Institute of Technology, United States; Armand Soldera, Université de Sherbrooke, Canada
Updates
Copyright
© 2020 Lei, Xu, Li and Liu.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zishun Liu zishunliu@mail.xjtu.edu.cn
This article was submitted to Polymer Chemistry, a section of the journal Frontiers in Chemistry
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.