A Multiaxial Theory of Double Network Hydrogels | Macromolecules

Article Cite This: Macromolecules XXXX, XXX, XXX−XXX

pubs.acs.org/Macromolecules

A Multiaxial Theory of Double Network Hydrogels Vu Ngoc Khiêm,*,† Thanh-Tam Mai,‡ Kenji Urayama,‡ Jian Ping Gong,§ and Mikhail Itskov† †

Department of Continuum Mechanics, RWTH Aachen University, Kackertstr. 9, 52072 Aachen, Germany Department of Macromolecular Science & Engineering, Kyoto Institute of Technology, Sakyo-ku, Kyoto 606-8585, Japan § Faculty of Advanced Life Science, Institute for Chemical Reaction Design and Discovery, Soft Matter GI-CoRE, Hokkaido University, Sapporo, Hokkaido 001-0021, Japan ‡

Downloaded via BUFFALO STATE on August 1, 2019 at 03:15:02 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: The incorporation of a stiff and brittle phase into a soft and flexible polymer network makes double network hydrogels remarkably tougher than the conventional one. It also induces a stress softening in cyclic loading, which is seemingly identical to that of filled rubbers. Therefore, material models proposed for double network hydrogel are mostly based on continuum damage theories of filled elastomers. However, recent experimental studies clearly distinguish double network hydrogels from filled elastomers by the damage cross-effect in multiaxial deformation. In this paper, a micromechanical model for double network hydrogels under multiaxial deformation is proposed on the basis of the analytical network-averaging concept. Accordingly, the directional fracture of the first network results from cumulative damage in all directions. Furthermore, the strong interpenetration of the two networks is taken into account. The proposed analytical model includes very few physically motivated material constants and demonstrates excellent agreement with comprehensive experimental data of double network hydrogels under multiaxial loading conditions. links that can break and then recover after one loading cycle.6 For a thorough review on tough polymer networks interested readers are referred to a recent study by Creton.7 To date, much effort has been made on constitutive modeling of double network hydrogels. Most of the works are based on isotropic continuum damage theories of filled elastomers (see e.g. ref 8), in which the value of the internal variable remains constant when changing the principal loading direction. Wang and Hong9 proposed the first constitutive theory for double network hydrogels using two internal variables. Zhao10 described the damage in double network hydrogels by means of the network alternation theory,8 taking into account different chain lengths in the two networks. This motivated Wang and Gao11 to construct another network alternation theory accounting for chain length distributions in different directions. Bacca et al.12 developed a generalized model for multiple network elastomers on the basis of a neo-Hookean-Gent free energy, taking into account the swelling of the ground network in each synthesis step. The damage in each network is elucidated by an evolution of the limit of chain extensibility resulting from chain scission. Note that due to the assumption of the damage equal in all directions, isotropic continuum damage models cannot capture the strain-induced anisotropy13,14 as well as the

1. INTRODUCTION Hydrogels are cross-linked networks of hydrophilic polymer molecules dispersed in water. They possess a number of promising properties, such as biocompatibility, environmental friendliness, adhesion, antibiofouling behavior, and biodegradability. However, relatively poor mechanical strength of synthetic hydrogels limited their further advancement and applications. A breakthrough in synthesis of hydrogels was made in 2003. By introducing a stiff and brittle network (the first network) into a soft and flexible polymer (the second network), Gong and colleagues1 were able to produce a double network hydrogel whose stiffness and toughness are comparable to natural hydrogels (such as tendons, ligaments, and cartilages). Under deformation, the brittle network breaks into smaller fragments while still remaining connected to the second network.2,3 The double network hydrogel can dissipate large energy and is significantly tougher than the conventional one (single network). The breakage of covalent cross-links in the first network also induces an irreversible stress softening in cyclic loading of double network hydrogels, which is seemingly identical to the Mullins effect in filled rubbers. In recent years, alternative synthesis methods have been developed to introduce the dissipation mechanism into hydrogels, such as by using ionically cross-linked first network4 or by incorporating dipoledipole and hydrogen-bonding interactions.5 To overcome the low fatigue life of double network hydrogels, self-healing hydrogels have been synthesized by including physical cross© XXXX American Chemical Society

Received: May 22, 2019 Revised: July 11, 2019

A

DOI: 10.1021/acs.macromol.9b01044 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 1. Successive fracture of an exemplary part of the first network into several solid fragments. Chains A and B both belong to the second network, but only chain A survives after the breakage of the first network molecules. For the sake of simplicity, only chains (i.e., parts of the polymer strand) are shown.

permanent set in double network hydrogels15 (see ref 16 for detailed discussions on isotropic and anisotropic finite strain continuum damage theories). Liu et al.17 developed an anisotropic microsphere damage model for double network hydrogels with two internal variables in each direction of the unit sphere. This model aimed to describe the necking in uniaxial tension of double network hydrogels by means of a nonconvex constitutive formulation. Khiêm and Itskov18 proposed a constitutive model capturing mechanically induced chemiluminescence in the double network19 as well as filled elastomers.20 The model is based on an anisotropic microvoid growth theory and demonstrated very good agreement with both stress-strain response and light emission intensity measured in uniaxial cyclic tension tests of mechanochemically responsive elastomers. Long et al.21 formulated a constitutive model describing rate-dependent self-healing in dual cross-link hydrogels. Accordingly, reversible stress softening is explained by breakage and reattachment of physical cross-links. Mao et al.22 described the rate-dependent damage in ionically and covalently cross-linked double network hydrogels by decomposing the strain energy into three parts corresponding to elasticity, rate-independent damage, and viscoelasticity. Recently, Vernerey et al.23 constructed an affine damage model for multiple network elastomers.19 In this model, the distribution of chain end-to-end vectors was driven by the fraction of ruptured chain and damage evolution followed the principle of maximum damage dissipation (see ref 16). All of the aforementioned constitutive models were only compared to uniaxial tests of double network hydrogels.24 For example, the damage evolution in the work of Wang and Hong9 was compared with uniaxial tension and compression.25 However, boundary conditions of uniaxial compression and equibiaxial tension are different, so that this validation is still insufficient for characterizing mechanical behavior of double network hydrogels (see e.g. ref 26 and Figure S3 of the Supporting Information). Furthermore, a recent experimental result of double network hydrogel under multiaxial deformation demonstrated a strong damage cross-effect15 which was not revealed by uniaxial tests.9,27 Accordingly, the reductions in stiffness of the material in equibiaxial and biaxial tension are significantly larger than that one in uniaxial tension (and pure shear) at the same level of deformation. Thus, damage in one direction is controlled by the history of deformation not only in this direction but also in other directions. Thus, a proper

constitutive theory for double network hydrogels under multiaxial loading conditions is still missing. In the current work, we present a micromechanical model describing the multiaxial damage in covalently cross-linked double network hydrogels. The damage mechanism in the first network is postulated as a successive fracture of a series of network blocks. To formulate the damage evolution, a novel type of internal variable is introduced on the basis of a recently developed statistical framework (which is termed the analytical network averaging concept18,26,28). Since this approach can account for the spatial organization of polymer strands,18,28 strong interpenetration of the two networks in double network hydrogels can be captured. In other words, the internal variable in one direction depends on deformation history in all directions, which elucidates the cross-effect of multiaxial damage. The paper is organized as follows: in section 2 we present the successive fracture mechanism of block series in the first network, which serves as the foundation of the proposed model (section 3). To describe multiaxial damage in double network hydrogels, the new type of the internal variable is defined in section 3.2. It is then connected to microstructural changes in double network hydrogels via a survival fraction of chains in the first network (section 3.3). On this basis, a micromechanical model is then developed (sections 3.4 and 3.5). Finally, predictive capabilities of the model are illustrated and discussed in section 4 by using comprehensive multiaxial tension experimental data of double network hydrogels. For reader’s convenience, the concept of analytical network averaging is recalled in the Appendices.

2. DAMAGE IN DOUBLE NETWORK HYDROGELS To properly model the mechanical behavior of a double network hydrogel, its microstructural damage, dissipation, and reinforcement mechanisms should be understood. In this study, the network damage is elucidated by successive fracture of the first network. To this end, we model the first network by blocks in series. Each block is composed of short polymer chains (i.e., partial strands, see section 3.1) and is interpenetrated by long chains from the second network. The polymer chains of the second network form entanglements with the first network that are trapped topologically. Under deformation, bonds in the first network whose deformation energy exceeds the critical value break, and smaller blocks are separated from the first network. B


Article

Macromolecules

Figure 2. Difference between microstructural changes in (a) a weakly interpenetrated double network (the second network chains are only entangled with the first network at its boundary) and (b) a strongly interpenetrated double network (the second network chains get entangled at the center of the first network). Note that the broken first network block always entangles with the second network even though some of the entanglements are deactivated. For the sake of simplicity, only chains contributed to the damage are shown.

whose end-to-end vector is in X-direction. Whereas, because of the strong interpenetration of chains in the second network into the first network, the extension in Y-direction causes rupture of bonds in the first network and releases the entanglement of chain B with the first network (Figure 2b). Thus, chains of the second network illustrated by chain B (in X-direction) are deactivated by a deformation in Y-direction, indicating the damage crosseffect.

This damage process can be imagined by considering an exemplary part of the first network containing four blocks in series (Figure 1). They are interpenetrated by two long chains from the second network denoted by A and B. As seen in Figure 1, after two deformation steps, the blocks in the middle connected by chain B lose their anchors, and therefore chain B is disentangled from the first network and becomes deactivated. In other words, only chain A survives after the fracture of the first network. One can observe that during the above damage process the number of Kuhn segments of the second network remains the same, whereas the number of active chains decreases (in this example, it reduces from two chains to one). It sets our approach apart from those based on the network alternation theory9−12 in which the number of Kuhn segments increases.8 Furthermore, because of the existence of the solid phase (the first network) in double network hydrogels, the stretch in the second network chains is larger compared to that in a single network (the socalled strain amplification effect29). The above damage process induces mechanical dissipation and requires more energy for the crack growth that indicates higher fracture toughness of double network hydrogels (compared to the single network gel). As aforementioned, the damage cross-effect results from strong interpenetration of the two networks. It can be demonstrated in two types of hydrogel with weak and strong interpenetration by considering microstructural changes in the network caused by a deformation in Y-direction (Figure 2). In the weakly interpenetrated double network (Figure 2a), the deformation in Y-direction only breaks the first network without any influence on the second network illustrated by chain B

3. THEORY 3.1. Terminology. In this paper, we consider a polymer network by a composition of many linear polymer strands dispersed in different directions n in space (Figure 3).26 A strand is defined by a network path between two permanent cross-links whose functionalities could be three or more. Each polymer strand contains a series of chains. A chain (so-called substrand) is defined as a part of strand between two successive entanglements of quadrafunctional cross-links.

Figure 3. Polymer network, a strand, and a chain. C


Article

Macromolecules

Figure 4. Difference between isotropic and anisotropic continuum damage: (a) deformation history, the sample is subject to a cycle of pure shear; (b) (top) corresponding survival fraction of the first network in the current model (α̅ = 20, n̅ = 15, and L̅ = 1); (b) (bottom) the survival fraction in an isotropic damage model9 (η0 = 0.9, Λ0 = 1.5, and d = 0.5). max

3.2. Internal Kinematic Variable. The concept of internal variable30 has frequently been employed in continuum mechanics to describe solids undergoing irreversible thermodynamic processes. The internal variable is supposed to be associated with microstructural change inside a material (such as bond rupture in filled elastomers8,16,27) and is therefore associated with the deformation history. In this study, we define internal variables on the basis of kinematic measures. As will be shown below, it is convenient for the derivation of damage evolution. As discussed above, the strong interpenetration of the two networks causes damage cross-effect in double network hydrogels. To capture this effect, we make use of the analytical network averaging concept18,26,28 because this approach can take into account the self-organization of molecules in multiple directions.18,28 Accordingly, the first network is considered as a composition of one-dimensional strands dispersed around three principal strain directions. A directional probability distribution function ρ(n) of polymer strands (eq 23) is assumed. All polymer strands dispersed around a mean direction are replaced with a representative strand in this direction subject to an average stretch. By using ρ(n), we can obtain average deformation measures applied on the representative strand by analytical averaging of the macroscopic deformation. In particular, the average stretch Λ̅ i is evaluated as the root-

Λ̅ i

max

Υ̅ i

n

y

1/2

{

1/2 I i y = jjj(1 − wi) 1 + wi Λi 2zzz , 3 k {

i = 1, 2, 3

(2)

=

max Υ̅i(τ ),

τ∈ (−∞ , t ]

i = 1, 2, 3

(3)

where Υ̅ 1, Υ̅ 2, and Υ̅ 3 denote the average area contractions in three principal strain directions. Similarly to eq 1, their values can be obtained from the root-mean-square of the macroscopic n

change in area Υ as16,18,28 i Υ̅i = jjjj k

y

∫S ρ(n)Υ2 dSzzzz n

1/2

{

1/2 I i y = jjj(1 − wi) 2 + wi Υi 2zzz , 3 k {

i = 1, 2, 3

(4)

where I2 represents the second principal invariant of C, Υi2 = cof C: Ei⊗Ei is the square of the change in infinitesimal area with the normal direction Ei, and the cofactor of the right Cauchy−Green tensor is given by cof C = C−T det C (see e.g. ref 31). 3.3. Evolution of Multiaxial Damage in the First Network. In this section, we show explicitly how the internal variables defined in section 3.2 are linked to bond rupture in the first network. Hereafter, variables related to the first network are indicated by an overbar as ●̅ . The survival fraction of chains in the first network (relative to the number of first network chains in the reference state) in a direction i resulting from the successive fracture of a series of polymer blocks can be given by (for more details see Appendix B)

n

∫S ρ(n)Λ2 dSzzzz

max Λ̅ i(τ ),

τ∈ (−∞ , t ]

where Λ̅ 1, Λ̅ 2, and Λ̅ 3 denote the average stretches in three principal strain directions. τ denotes a time point in the past. In the same manner, the maximal average area contraction Υ̅ max previously reached in the deformation history can be i defined as

mean-square of the macroscopic stretch Λ over the unit sphere S (for more details see Appendix A) by i Λ̅ i = jjjj k

=

i = 1, 2, 3

(1)

Ä ÉÑ max ij α ÅÅÅ Ñy 2π Λ̅ i L̅ max Ñz − 2π Λ̅ i L̅ ÑÑÑÑzzzz, Φi = Φ0 expjjjj ̅ ÅÅÅÅn ̅ sin ÑÑÑz j 4π ÅÅÇ n̅ Ö{ k

where I1 represents the first principal invariant of the right Cauchy−Green tensor C = FTF, F is the deformation gradient, and Λi2 = C:Ei⊗Ei is the square of the principal stretch in a principal direction specified by a unit vector Ei. Furthermore, wi = w01/3 represents the directional fraction of oriented first network strands in each principal strain direction i (i = 1, 2, 3). w0 is the effective volume fraction of the swollen first network relative to the total volume of the first and second network in the reference configuration.18 The first internal variable is defined as the maximal average stretch previously reached in the deformation history and denoted as Λ̅ max i . It indicates the maximal stretch ratio that a directional strand has experienced and is given at time t by

i = 1, 2, 3 (5)

where α̅ is the mode likelihood parameter indicating the maximum probability of finding a strand between two cross-links and L̅ is the effective referential end-to-end distance of the first network (by division by the number of chains in the strand; see Appendix B). n̅ is the number of molecular segments in each chain (the smallest possible block) of the first network and Ä É α Å 2π L Ñ Φ0 = exp 4π̅ ÅÅÅÅ2πL̅ − n ̅ sin n ̅ ÑÑÑÑ is the normalization factor. ÅÇ Ö ̅ Ñ

(

D

)


Article

Macromolecules

The mesoscopic stretch Λ̅ i (applied at the strand) and the microscopic stretch λ̅i of the chain in the first network can be related by a stretch amplification exponent q0 as follow:s26

This formulation expresses the remaining percentage of the first network subjected to the deformation Λ̅ max i . Since the internal kinematic variable is stored as a function of deformation, the available fraction of the first network Φi is deformation history dependent. It can be proved that the so-described damage evolution follows the principle of maximal damage dissipation.16 Furthermore, the internal variables are not the same in different directions. This characteristic is illustrated in Figure 4. In this example, a double network hydrogel sample is subject to pure shear deformation up to stretch ΛX = 2.5. As seen in Figure 4b (top), the damage is different in X- and Y-direction, whereas the damage is identical in all directions of isotropic damage models (Figure 4b (bottom)). Furthermore, the cross-effect of damage can also be observed in this example (Figure 4b (top)). Under the tensile loading in X-direction, damage already takes place in Y-direction (0−100 [s]). To describe the statistical mechanics of polymer chains in double network hydrogels, we make use of the recently developed tube model.26 Accordingly, each polymer chain in a directional strand is considered as a freely jointed chain (ideal chain) confined within a tube formed by neighboring chains (Figure 5). Therefore, the free energy Ψ of a virgin single

q

i = 1, 2, 3

λi̅ = Λ̅ i 0,

(8)

The parameter q0 depends on the number of cross-links in the strand.26,28,33 To calculate the microscopic area contraction, we assume directional incompressibility in a single network so that the microscopic area contraction υ̅i can be directly linked to the mesoscopic area contraction Υ̅ i by q

i = 1, 2, 3

υi̅ = Υ̅ i 0,

(9)

The microscopic area contraction υ̅i represents the relation of the reference to the current cross-section area of the tube, whereas the mesoscopic area contraction is defined in eq 4. The total number of active chains in the double network hydrogel decreases when the first network ruptures (see section 2). We further assume gradual fracture of the first network so that every survival polymer block in the first network connects to the same number of second network chains during the whole course of deformation. Therefore, the fraction of all chains available in each principal direction Ni, i = 1, 2, 3, is proportional to the survival fraction of the first network in this direction as Ni =

N0 Φi , 3

i = 1, 2, 3

(10)

Figure 5. Tube model of a single polymer chain. Filled circles represent cross sections of neighboring chains.

In view of eqs 7 and 10, the ideal shear modulus is proportional to the number of available chains. Therefore, the ideal shear moduli of the first network (with the volume fraction w0) can be given for each direction i by μ μci̅ = w0 c0 Φi , i = 1, 2, 3 (11) 3

network can be split into two parts related to the ideal chain and topological constraint induced by neighboring molecules26 as

Furthermore, we take into account that the topological constraint is released by damage in all directions.16 Therefore, the topological shear modulus of the first network (with the volume fraction w0) is given by

3

Ψ=

∑ i=1

μc0 kBT

ψc(λi) +

μt0 kBT

ψt(υi)

μt̅ i = w0

(6)

Therein, ψc (λi) is the free energy of the ideal chain defined as a function of microscopic stretch λi (change in the length of the tube). ψt(υi) is the topological free energy which is formulated in terms of microscopic area contraction υi (change in the crosssection area of the tube). The shear moduli corresponding to the two parts in eq 6 are termed the initial ideal and topological shear moduli. According to the classical statistical theory of rubber elasticity, they are linked to the initial number of chains in the double network hydrogel N0 by the following relations:26,32 μc0 = N0kBT ,

μt0 = N0kBTω

μt0 3

3

∏ ΦΩ i i, i=1

i = 1, 2, 3 (12)

Therein, the released part due to area contraction Ωi is defined in view of eq 28 by Ä ÉÑ max ij α ÅÅÅ Ñyz 2π Υ̅ i L̅ max Ñ Å ̅ j Å j Ωi = Φ0 expjj ÅÅn ̅ sin − 2π Υ̅ i L̅ ÑÑÑÑzzzz, Å ÑÑÖ n̅ k 4π ÅÇ { i = 1, 2, 3 (13) 3.4. Microstructural Change of the Second Network. Similarly to the first network, the second network is also modeled by three representative strands in three principal strain directions. As aforementioned, the fragments of the first network present a solid phase, whereas the second network is considered as a soft phase. Therefore, the microscopic stretch in the second network is larger than the macroscopic draw ratio Λi. This nonaffine deformation can be accounted for by a stretch amplification exponent qi as follows:26

(7)

where kB is the Boltzmann constant, T is the absolute π2 n

temperature, ω = D2 is a geometrical constant of the tube, D is the tube diameter in the reference state, and n denotes the number of Kuhn segments of molecules in the network.26 In the following we shall determine the microscopic stretch λ̅i and microscopic area contraction υ̅ i of the first network when the double network gel is subject to deformation gradient F. We assume that the principal direction of a macrostretch always coincides with that one of a microstretch.

q

λi = Λ i i ,

i = 1, 2, 3

(14)

The change of stretch amplification exponent during the course of deformation can be given by16 E


Article

Macromolecules qi =

q0 1/3

1 − (w0 Φi)

3 N ijj ∂ψ (n ̅ , λi̅ ) ∂λi̅ ∂Ψ − pF−T = w0 ∑ 0 jjjΦi c j ∂F ∂λi̅ ∂F i=1 3 j k 3 3 ∂υi̅ yzzz N0 ijjj ∂ψc(n , λi) ∂λi z jjΦi + − + ∏ ΦΩ ω (1 ) w ∑ z j j 0 j ∂F zz ∂λi ∂F j=1 i=1 3 j { k 3 ∂υi yzzz zz − pF−T + ∏ ΦΩ j jω ∂F zz j=1 (21) {

i = 1, 2, 3

,

P=

(15)

where (w0Φi)1/3 represent the directional fraction of the solid phase. In the same manner with eq 9, the microscopic area contraction υi can be obtained as q

i = 1, 2, 3

υi = Υ i i ,

(16)

The shear moduli corresponding to the second network (with a volume fraction 1 − w0) can be expressed in view of eqs 11 and 12 as μci = (1 − w0)

μc0 3

Φi ,

μti = (1 − w0)

μt0 3

where

3 i=1

i μci̅

zy ψt(υi̅ )zzz z kBT i = 1 k kBT { 3 ij μci yz μti + ∑ jjj ψ (n , λ i ) + ψt(υi)zzz j kBT c z k T B i=1 k { 3 3 i yz N w jj zz z ( ) = ∑ 0 0 jjjΦiψc(n ̅ , λi̅ ) + ∏ ΦΩ ωψ υ j j t i̅ z jj zz 3 j=1 i=1 k { 3 3 i yz N (1 − w0) jjj zz jjΦiψc(n , λi) + ∏ ΦΩ z ( ) +∑ 0 ωψ υ j j t i z j zz 3 j j=1 i=1 k {

∑ jjjjj 3

ψc(n ̅ , λi̅ ) +

μt̅ i

ψc(n , λ) = nkBTκ ln

( )

n sin

(19)

and κ = 9/π2 represents a universal constant. The free energy (eq 19) arises from a closed-form of the exact non-Gaussian statistics.26 The topological free energy takes the form26 ψt(υi) = kBTυi

πλi n

( ) and p denotes a

parameters

physical meaning

μc0 n̅ n q0 μt0 α̅ L̅

initial ideal shear modulus [MPa] average number of Kuhn segments of the first network chains average number of Kuhn segments of the second network chains stretch amplification exponent in a single network gel initial topological shear modulus [MPa] mode likelihood parameter referential distance of strands in the first network

4. RESULTS AND DISCUSSION In this section, we consider a benchmark example to validate the performance of the model presented above. First, the model is fitted to uniaxial tension, pure shear, and equibiaxial tension data of covalently cross-linked double network hydrogels to determine its material constants. Afterward, the same set of material parameters is used to predict stress-strain response of the same double network hydrogel in a biaxial tension test. In Mai et al.,15 multiaxial tension (uniaxial tension, pure shear, equibiaxial tension, and biaxial tension) of covalently crosslinked double network hydrogels was performed at 25 °C. The first network of the double network hydrogel was synthesized by using 2-acrylamido-2-methylpropanesulfonic acid sodium salt (NaAMPS) as the monomer and N,N′-methylenebis(acrylamide) (MBAA) as the cross-linking agent. Synthesis of the second network (polyacrylamide, PAAm) in the presence of the first one results in a double network hydrogel PNaAMPS/ PAAm. The volume fraction of the swollen first network relative to the total volume of the first and second network in the reference configuration is estimated as 35% (see section 1 in the Supporting Information). Since the mechanical response of this double network hydrogel is rate-independent,15 it is suitable for the validation of our quasi-static constitutive theory (eq 18). In the first step, the model (18) is fitted to the first loading cycles of uniaxial tension (up to stretch 1.97), pure shear (up to stretch 1.88), and equibiaxial tension data (up to stretch 1.62) of PNaAMPS/PAAm hydrogels. The obtained material constants

The former term of eq 18 corresponds to the energy contribution of the first network, while the latter describes that of the second network. The free energy of freely jointed chains in eq 18 is given by26 π λ n i πλi n

i

meaning of the material parameters (μc0, n, q0, and μt0) has been demonstrated for single network gels by comparison of our model prediction with molecular dynamic simulations.16

(18)

sin

n

= kBTκ λ − kBTκπ n cot

Table 1. Physical Meaning of the Material Constants (17)

3.5. Constitutive Model. In view of eqs 6, 10, 11, 12, and 17, the free energy of double network hydrogels can be formulated as Ψ=

∂λi

Lagrange multiplier. In view of eq 7, the free energy (eq 18) results in a constitutive model with seven material constants (Table 1). The physical

∏ ΦΩ i i,

i = 1, 2, 3

∂ψc(n , λi)

(20)

Arguments of the free energy functions in eqs 19 and 20 are applied both for the first and the second networks. For the sake of simplicity, we do not explicitly include the decomposition of deformation gradient12 into a swelling and a tension part in eq 18. Because of incompressibility of double network hydrogels,15,25 the first Piola−Kirchhoff stress tensor can be derived as

Table 2. Experimentally Derived Parameter w0 F

35% DOI: 10.1021/acs.macromol.9b01044 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

value L̅ ≈ 6.7 of the normalized distance indicates that in the first network the referential length of the strand is of 67 statistical Kuhn segment lengths. Very good agreement between the model prediction and experimental data with respect to stress softening at other stretch amplitudes in uniaxial tension, pure shear, and equibiaxial tension is automatically obtained (Figure 6). Note that many other models, as for example the Gent model, are even not able to fit simultaneously to a simplified version of the pure shear data in X- and Y-direction15 (in which only unloading curves were considered). We further used the same values of material constants (Table 3) to plot nominal stresses predicted by the proposed model (18) in two directions of the (unequal) biaxial tension test (strain ϵ2 = 0.5ϵ115). As seen in Figure 7, the model demonstrates very good agreement with experimental data as well. The dissipation density D is specified by15

are given in Table 3. Note that for the value of amplification exponent q0 = 2, our model reproduces the tube model by Table 3. Material Constants of the Proposed Model for Covalently Cross-Linked Double Network Hydrogel μc0 [MPa] 0.001 75

n̅ 28.9462

n

q0

μt0 [MPa]

294.8085

2

0.0067

α̅ 0.298

L̅ 6.7068

Heinrich and colleagues33,34 for each single network (w0 = 0). Since material parameters of the Heinrich model have physical and structural meanings (see e.g. ref 35), q0 = 2 is physically relevant. Linearization of eq 18 results in the initial modulus μ E = 3μc0 + 2t0 q2 . By using the initial modulus EUN = 0.4

(

)

[MPa], evaluated at the strain 1%, one can obtain the corresponding ideal and topological shear moduli. The network structural parameters μc0 and μt0 can be obtained from different experimental techniques such as equilibrium swelling or low field nuclear magnetic resonance.36 In view of eqs 11 and 12, the ratio ω = μt0/μc0 = 3.8 suggests that the total cross-link density of the double network hydrogel is relatively low, so that the topological constraint is still prominent (see e.g. refs 16 and 35). Since it is difficult to measure separately the mean values of the chain length n̅ and n in the two networks, these parameters are generally obtained from fitting. Nevertheless, it can be seen that the chain length of the first network is very short compared to the second network. This indicates that the first and second networks are tightly and loosely cross-linked, respectively. If one assumes that each strand contains for example 10 chains, the

2

D=

2

∑∫ i=1

virgin loading

Pi dΛi −

∑∫ i=1

unloading

Pi dΛi

(22)

where Pi are the principal nominal stress in direction i, i = 1, 2 2 and W0 = ∑i = 1 ∫ Pi dΛi is the stored energy density in virgin loading

the primary loading. To further demonstrate the effect of multiaxiality, the stored energy density W0 and the dissipation density D predicted by the model are plotted versus stretch in Xdirection (Λ1) with respect to different states of deformation. The areas are calculated using trapezoidal numerical integration.

Figure 6. Prediction of the proposed model in comparison to multiaxial tension data:15 (a) uniaxial tension, (b) equibiaxial tension, (c) pure shear (stress in the tension direction), and (d) pure shear (stress in the constrained direction). G


Article

Macromolecules

Figure 7. Prediction of the proposed model in comparison to unequal biaxial tension (strain ϵ2 = 0.5ϵ1) data: (a) stress in the first direction of the biaxial tension; (b) stress in the second direction of the biaxial tension. Material constants are given in Table 3.

Figure 8. Effect of multiaxiality in the model and comparison to experimental data:15 (a) stored energy density in the primary loading; (b) dissipation density. Material constants are given in Table 3.

Figure 9. (a) Survival fraction of the first network in the tensile direction Φ1 under multiaxial loading. (b) Reduction the total shear modulus and comparison to experimental data.15 Material constants are given in Table 3.

cross-linked double network hydrogels was reported at Λ ≈ 2, whereas the limit of chain extensibility in the second network was √n ≈ 17 (for the values of material parameters used for the simulations). To further investigate this phenomenon, the number of chains remaining in the first principal direction of the first network is plotted for different deformation states in Figure 9a. One can observe that the survival fraction of first network reduces dramatically at Λ ≈ 2 for both states of deformation, which is in line with the experimental observation. Furthermore, the directional survival fraction of the first network is not the same in uniaxial tension and equibiaxial tension because of the damage cross-effect (Figure 9a). It can also be seen that the

As seen in Figure 8, the model shows very good agreement with experimental data. Thus, the model can accurately capture not only the stress (Figures 6 and 7) but also the energy change (Figure 8) during deformation of double-network hydrogels. Accordingly, the dissipation in equibiaxial tension and biaxial tension is much larger than that in uniaxial tension, when compared at the same value of Λ1, indicating the effect of multiaxiality. It is explained by the strong interpenetration of the two networks described by the analytical network averaging (1). Furthermore, the damage in the presented model is solely caused by bond rupture in the first network but not by failure of the second network chains. Indeed, the yielding of the covalently H


Article

Macromolecules

Figure 10. Permanent deformations predicted by the model under multiaxial cyclic loading: (a) equibiaxial tension and pure shear; (b) uniaxial tension and biaxial tension.

survival fraction of the first network in the constrained direction of pure shear test Φ2 is less than unity (see Figure S2d), which indicates the cross-effect of strain in the tensile direction (Λ1 > 1) on the damage in the constrained direction (Λ2 = 1). Since the presented damage formulation is multiaxial, the reduction in total shear modulus results from damage in three principal directions. As seen in Figure 9b, the model shows qualitatively good agreement with experimental data. The discrepancy in the small strain regime may be induced by the finite extensibility effect (see e.g. ref 26) which is not reflected in the physical shear moduli (eqs 11 and 12). Furthermore, it may also arise from the stiffening effect on polymer chains in the highly swollen first network37 which is not considered in the current model. The presented model belongs to the framework of anisotropic continuum damage theory and can thus capture a permanent set. Residual stretches remaining after cyclic loading are plotted versus maximal stretch amplitude in Figure 10 for different deformation states. The model prediction (Figure 10) slightly overestimates the experimental values; the source of this discrepancy remains unclear, since other predictions show very good agreement with the mechanical response of double network hydrogel at small deformations (see Figure S2 for the model prediction at small strains with the set of material constants given in Table 3). Nevertheless, the proposed model qualitatively captures the evolution of the permanent set in different deformation states. Moreover, as seen in Figure S3, isotropic damage models (such as ref 9) cannot describe the permanent set at all.

network hydrogels has not been obtained before. Since material parameters of the model are physically relevant, it can be useful for designing hydrogels with the desired mechanical response.

■

APPENDIX A The analytical network averaging concept18,26,28 assumes the existence of an orientational distribution function of polymer strands in polymeric materials. For an arbitrary direction n, this function can be represented by a probability density ρ related to a field direction Ei (in case the polymer is subject to mechanical force, the direction Ei coincides with the principal strain direction) as follows:16,38 ρ(n) = ρ(θ , ςi , Ei) =

ςi cosh(ςi cos θ) 4π sinh(ςi)

(23)

where θ denotes the angle between the directional unit vector n and the field direction Ei and ςi is the concentration parameter. According to the analytical network-averaging concept, the average stretch Λ̅ i in direction i is evaluated as the root-meansquare of the macroscopic stretch over the unit sphere by16 n

2

Λ̅ i =

∫S ρ(n)Λ2 dS 2π

= C:

∫0 ∫0

i = C : jjj k

π

2π

∫0 ∫0

5. CONCLUSION In this paper, we proposed a micromechanical model to capture multiaxial damage in covalently cross-linked double network hydrogels. Accordingly, the stress softening and toughening mechanism in double network hydrogels were explained by the successive fracture of a series of network blocks in the first network. This effect was described within the context of thermodynamics of internal variable. To this end, we further explored the concept of analytical network averaging and defined a new type of internal variable. This internal variable is derived by averaging of the deformation measure, taking into account the dispersed distribution of the first network within the second network. Thus, the interpenetration of two polymer networks can be captured, which elucidates the damage crosseffect. To the best of our knowledge, such excellent agreement with multiaxial experiments (including uniaxial tension, pure shear, equibiaxial tension, and biaxial tension) of double

2π

+

∫0 ∫0

+

∫0 ∫0

π

ρ(n)n ⊗ n sin θ dθ dφ π

ρ(n) sin 3 θ cos2 φ dθ dφ H1 ⊗ H1

ρ(n) sin 3 θ sin 2 φ dθ dφ H2 ⊗ H2

y ρ(n) cos2 θ sin θ dθ dφ Ei ⊗ Ei zzz { i y I = C : jjj(1 − wi) + wiEi ⊗ Ei zzz 3 k { I1 = (1 − wi) + wi Λi 2 (24) 3 2π

π

n

where Λ denotes the macroscopic stretch in direction n and θ ∈ [0, π] and φ ∈ [0, 2π] are spherical coordinates of the vector n = sin θ cos φH1 + sin θ sin φH2 + cos θEi with respect to orthonormal vectors (H1, H2, Ei). The coefficient wi is expressed by I


Article

Macromolecules 2π

wi =

∫0 ∫0

dθ dφ =

π

then cof C represents the corresponding multiaxial compression). Thus, Υ̅i = cof C : Ei ⊗ Ei , i = 1, 2, 3, play the role of the principal stretches in the inverse deformation. Therefore, the survival fraction of the first network in direction i caused by this inverse deformation can be given in view of eq 27 by

ρ(θ , ςi , Ei)(cos 2 θ sin θ − sin 3 θ cos2 φ)

ςi2 − 3 coth(ςi)ςi + 3 ςi2

(25) 39,40

wi represent fractions of oriented strands in each direction i (i = 1, 2, 3).18,28 As seen in eq 24, wi control the anisotropy degree in the model, and eq 24 reduces to an isotropic one in case wi = 0 (see refs 26 and 28).

■

■

APPENDIX B As aforementioned, each directional strand of the first network is composed of a series of Mc polymer chains connecting two high functional cross-links. Since the first network is inhomogeneous, its strands can have different number of polymer chains with a maximum number Mmax c . The probability PD (lc) of existence of a polymer strand with the contour length lc (normalized by division by the Kuhn length) connecting two high functional cross-links can be assumed as16 ÄÅ α(l − 1) ÉÑ πl Å Ñ α sin 2 M maxc n expÅÅÅ c2 ÑÑÑ ÅÇ ÑÖ c ̅ PD(lc) = ÄÅ max ÉÑ ÅÅ αMc n̅ Ñ 2π l 2π expÅÅÅ 4π sin M maxc n − sin M max n ÑÑÑÑ c c ̅ ̅ Ñ (26) ÅÇ Ö

( ) (

∫Λ

Mcmax n

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.macromol.9b01044. Brief information about the volume fraction w0 of the swollen first network; Figures S1−S3 (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected].

)

ORCID

Vu Ngoc Khiêm: 0000-0002-8749-364X Thanh-Tam Mai: 0000-0001-6190-1142 Kenji Urayama: 0000-0002-2823-6344 Jian Ping Gong: 0000-0003-2228-2750 Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS M. Itskov thanks the German Research Foundation (DFG) for the financial support of this work (project IT 67/13-1). REFERENCES

(1) Gong, J.; Katsuyama, Y.; Kurokawa, T.; Osada, Y. DoubleNetwork Hydrogels with Extremely High Mechanical Strength. Adv. Mater. 2003, 15, 1155−1158. (2) Na, Y.-h.; Tanaka, Y.; Kawauchi, Y.; Furukawa, H.; Sumiyoshi, T.; Gong, J. P.; Osada, Y. Necking Phenomenon of Double-Network Gels. Macromolecules 2006, 39, 4641−4645. (3) Brown, H. R. A Model of the Fracture of Double Network Gels. Macromolecules 2007, 40, 3815−3818. (4) Sun, J.-Y.; Zhao, X.; Illeperuma, W. R. K.; Chaudhuri, O.; Oh, K. H.; Mooney, D. J.; Vlassak, J. J.; Suo, Z. Highly stretchable and tough hydrogels. Nature 2012, 489, 133−136. (5) Zhang, Y.; Li, Y.; Liu, W. Dipole-dipole and h-bonding interactions significantly enhance the multifaceted mechanical properties of thermoresponsive shape memory hydrogels. Adv. Funct. Mater. 2015, 25, 471−480. (6) Dai, X.; Zhang, Y.; Gao, L.; Bai, T.; Wang, W.; Cui, Y.; Liu, W. A mechanically strong, highly stable, thermoplastic, and self-healable supramolecular polymer hydrogel. Adv. Mater. 2015, 27, 3566−3571. (7) Creton, C. 50th Anniversary Perspective: Networks and Gels: Soft but Dynamic and Tough. Macromolecules 2017, 50, 8297−8316. (8) Marckmann, G.; Verron, E.; Gornet, L.; Chagnon, G.; Charrier, P.; Fort, P. A theory of network alteration for the Mullins effect. J. Mech. Phys. Solids 2002, 50, 2011−2028. (9) Wang, X.; Hong, W. Pseudo-elasticity of a double network gel. Soft Matter 2011, 7, 8576. (10) Zhao, X. A theory for large deformation and damage of interpenetrating polymer networks. J. Mech. Phys. Solids 2012, 60, 319− 332.

Ä ij α ÅÅÅ ij 2π Λ̅ imax L 2π yz Å j − sin max zzz expjjj ÅÅÅMcmax n ̅ jjjsin max j j 4π ÅÅ Mc n ̅ Mc n ̅ z{ k k ÅÇ ÉÑ Ä Ñzy ij α ÅÅÅ 2π max Ñ Ñ + 2π (1 − Λ̅ i L)ÑÑÑzzzz − expjjjj ÅÅÅÅ−Mcmax n ̅ sin max ÅÅÇ π M 4 ÑÑÑÖz{ c n̅ k ÑÉÑy Ñz + 2π (1 − Mcmax n ̅ )ÑÑÑÑzzzz ≈ Φ0 ÑÑÖ { Ä max max É Å Ñ max Λ̅ i L ÑÑÑzyz 2π Λ̅ i L jij αMc ÅÅÅ ÅÅn ̅ sin ÑÑzz = Φ0 − π expjjj 2 ÅÅ Mcmax n ̅ Mcmax ÑÑÑÖz{ k 4π ÅÇ Ä ÉÑ max ij α ÅÅÅ Ñyz 2π Λ̅ i L̅ max Ñ Å ̅ j Å j − 2π Λ̅ i L̅ ÑÑÑÑzzzz expjj ÅÅn ̅ sin Å ÑÑÖ n̅ (27) k 4π ÅÇ { max i L

ASSOCIATED CONTENT

S Supporting Information *

where α is the mode likelihood parameter describing the maximum probability of finding a strand between two crosslinks. Let L be the referential end-to-end distance of the strand (normalized by division by the Kuhn length). The deformed end-to-end distance of the strand oriented in the direction i is Λ̅ max i L. Since polymer strands cannot be extended more than their contour length, the ones with the contour length lc shorter than Λ̅ max i L break. Therefore, the set of survival strands can be max max given by {lc| Λ̅ max i L ≤ lc < Mc n̅ }, where Mc n̅ denotes the maximal (normalized) contour length. Keeping in mind that exp(−Mmax c n̅) ≪ 1, the survival fraction of first network in direction i can be calculated as max Φi = Φ̂(Λ̅ i ) = Φ0

Ä ÉÑ max ij α ÅÅÅ Ñyz 2π Υ̅ i L̅ max max Ñ Å j ̅ ̂ Å j − 2π Υ̅ i L̅ ÑÑÑÑzzzz, Ωi = Φ(Υ̅ i ) = Φ0 expjj ÅÅn ̅ sin ÑÑz j 4π ÅÅÇ n̅ Ö{ k (28) i = 1, 2, 3

PD(lc) dlc = Φ0D

where Φ0 is the normalization factor, α̅ = αMmax c , and the L effective end-to-end distance of the strand L̅ = M max . c

Note that for incompressible single network gels characterized by the condition det C = 1, cof C = C−1 represents an inverse deformation of C (i.e., if C describes a multiaxial tension, J


Article

Macromolecules (11) Wang, Q.; Gao, Z. A constitutive model of nanocomposite hydrogels with nanoparticle crosslinkers. J. Mech. Phys. Solids 2016, 94, 127−147. (12) Bacca, M.; Creton, C.; McMeeking, R. M. A model for the Mullins effect in multi-network elastomers. J. Appl. Mech. 2017, 84, 121009. (13) Itskov, M.; Haberstroh, E.; Ehret, A. E.; Voehringer, M. C. Experimental Observation of the Deformation Induced Anisotropy of the Mullins effect in Rubber. Kautschuk Gummi Kunststoffe 2006, 59, 93−96. (14) Nakajima, T.; Kurokawa, T.; Ahmed, S.; Wu, W.-l.; Gong, J. P. Characterization of internal fracture process of double network hydrogels under uniaxial elongation. Soft Matter 2013, 9, 1955. (15) Mai, T.-T.; Matsuda, T.; Nakajima, T.; Gong, J. P.; Urayama, K. Distinctive Characteristics of Internal Fracture in Tough Double Network Hydrogels Revealed by Various Modes of Stretching. Macromolecules 2018, 51, 5245−5257. (16) Khiêm, V. N.; Itskov, M. An averaging based tube model for deformation induced anisotropic stress softening of filled elastomers. Int. J. Plast. 2017, 90, 96−115. (17) Liu, Y.; Zhang, H.; Zheng, Y. A Micromechanically Based Constitutive Model for the Inelastic and Swelling Behaviors in Double Network Hydrogels. J. Appl. Mech. 2016, 83, 021008. (18) Khiêm, V. N.; Itskov, M. Analytical network-averaging of the tube model: Mechanically induced chemiluminescence in elastomers. Int. J. Plast. 2018, 102, 1−15. (19) Ducrot, E.; Chen, Y.; Bulters, M.; Sijbesma, R. P.; Creton, C. Toughening elastomers with sacrificial bonds and watching them break. Science 2014, 344, 186−189. (20) Clough, J. M.; Creton, C.; Craig, S. L.; Sijbesma, R. P. Covalent Bond Scission in the Mullins Effect of a Filled Elastomer: Real-Time Visualization with Mechanoluminescence. Adv. Funct. Mater. 2016, 26, 9063−9074. (21) Long, R.; Mayumi, K.; Creton, C.; Narita, T.; Hui, C. Y. Time dependent behavior of a dual cross-link self-healing gel: Theory and experiments. Macromolecules 2014, 47, 7243−7250. (22) Mao, Y.; Lin, S.; Zhao, X.; Anand, L. A large deformation viscoelastic model for double-network hydrogels. J. Mech. Phys. Solids 2017, 100, 103−130. (23) Vernerey, F. J.; Brighenti, R.; Long, R.; Shen, T. Statistical Damage Mechanics of Polymer Networks. Macromolecules 2018, 51, 6609−6622. (24) Morovati, V.; Dargazany, R. Micro-mechanical modeling of the stress softening in double-network hydrogels. Int. J. Solids Struct. 2019, 164, 1−11. (25) Webber, R. E.; Creton, C.; Brown, H. R.; Gong, J. P. Large Strain Hysteresis and Mullins Effect of Tough Double-Network Hydrogels. Macromolecules 2007, 40, 2919−2927. (26) Khiêm, V. N.; Itskov, M. Analytical network-averaging of the tube model: Rubber elasticity. J. Mech. Phys. Solids 2016, 95, 254−269. (27) Merckel, Y.; Brieu, M.; Diani, J.; Caillard, J. A Mullins softening criterion for general loading conditions. J. Mech. Phys. Solids 2012, 60, 1257−1264. (28) Khiêm, V. N.; Itskov, M. Analytical network-averaging of the tube model: Strain-induced crystallization in natural rubber. J. Mech. Phys. Solids 2018, 116, 350−369. (29) Bueche, F. Mullins effect and rubber-filler interaction. J. Appl. Polym. Sci. 1961, 5, 271−281. (30) Coleman, B. D.; Gurtin, M. E. Thermodynamics with Internal State Variables. J. Chem. Phys. 1967, 47, 597. (31) Ogden, R. Non-Linear Elastic Deformations; Dover Publications: 1997. (32) Treloar, L. The Physics of Rubber Elasticity; Oxford University Press Inc.: New York, 1975. (33) Heinrich, G.; Straube, E.; Helmis, G. Rubber elasticity of polymer networks: Theories. Adv. Polym. Sci. 1988, 85, 33−87. (34) Klüppel, M.; Schramm, J. A generalized tube model of rubber elasticity and stress softening of filler reinforced elastomer systems. Macromol. Theory Simul. 2000, 9, 742−754.

(35) Marzocca, A. J.; Cerveny, S.; Raimondo, R. B. Analysis of the Variation of Molecular Parameters of NR During Vulcanization in the Frame of the Conformational Tube Model. J. Appl. Polym. Sci. 1997, 66, 1085−1092. (36) Basterra-Beroiz, B.; Rommel, R.; Kayser, F.; Westermann, S.; Valentín, J. L.; Heinrich, G. New Insights Into Rubber Network Structure By a Combination of Experimental Techniques. Rubber Chem. Technol. 2017, 90, 347−366. (37) Morozova, S.; Muthukumar, M. Elasticity at Swelling Equilibrium of Ultrasoft Polyelectrolyte Gels: Comparisons of Theory and Experiments. Macromolecules 2017, 50, 2456−2466. (38) Fisher, R. Dispersion on a Sphere. Proc. R. Soc. London, Ser. A 1953, 217, 295−305. (39) Murthy, N. S.; Bray, R. G.; Correale, S. T.; Moore, R. A. F. Drawing and annealing of nylon-6 fibres: studies of crystal growth, orientation of amorphous and crystalline domains and their influence on properties. Polymer 1995, 36, 3863−3873. (40) Pazur, R. J.; Prud’homme, R. E. X-ray pole figure and small angle scattering measurements on tubular blown low-density poly(ethylene) films. Macromolecules 1996, 29, 119−128.

K


A Multiaxial Theory of Double Network Hydrogels | Macromolecules

Recommend Documents