A Constitutive Model for Alginate-Based Double Network Hydrogels Cross-Linked by Mono-, Di-, and Trivalent Cations

In this contribution, a micro-mechanically based constitutive model is proposed to describe the nonlinear inelastic rubber-like features of alginate-based double network hydrogel cross-linked via various counterions. To this end, the lengthening of the polysaccharide polymer chain after a fully stretched state is characterized. A polymer chain is firstly considered behaving entropically up to the fully stretched state. Then, enthalpic behavior is accounted for concerning the following lengthening. To calculate enthalpic behavior, the macroscopic material properties, such as elastic modulus, are integrated into the proposed model. Thus, a new energy concept for a polymer chain is proposed. The model is constituted by the proposed energy concept, the network decomposition model, the Arruda–Boyce eight chain model and the network alteration theory. The model is compared against the cyclic tensile test data of alginate-based double network hydrogels cross-linked via mono-, di-, and trivalent cations. Good agreement between the model and experiments is obtained.


Introduction
In the field of bioengineering, biocompatible materials with particular and controllable structures are frequently required [1]. Such materials are utilized for applications, such as drug delivery [2], bulking agent [3] or extracellular matrix [4]. In this regard, hydrogels are one of the most important candidates in bioengineering.
Hydrogels are a class of three-dimensional polymer networks consisting of hydrophilic polymer chains and containing water up to 99.9% of their dry weights [5]. Due to their biocompatibility, water absorption capacity and environment-sensitive structures, they have attracted attention since the beginning of the 1950s in various fields ranging from cosmetics to the food industry and electronics to bioengineering applications [6]. Particularly, polysaccharide hydrogels are used as biomaterials in various applications due to their mimic ability to extracellular matrix role of living tissues [7]. However, the mechanical weakness of conventional hydrogels has obstructed the usage of these materials in many implementations requiring a load-bearing material [8].
Since the beginning of the last decade, many attempts have been carried out to strengthen and toughen hydrogels [9]. Nano-composite gels [10], slide ring gels [11], double network hydrogels (DNH) [12] are some of which proposed in the literature and show very promising mechanical properties for load-bearing applications.
Among tough hydrogels, a double network concept, which has so far provided one of the best mechanical properties for hydrogels, was proposed by Gong et al. [12]. A DNH is made up mainly of highly cross-linked polyelectrolyte and sparsely cross-linked neutral polymer networks. The polyelectrolyte network acts as a skeleton of the polymer matrix, whereas the neutral network keeps the structure together. To achieve tough hydrogels within this concept, some additional conditions with regard to cross-linked networks, cross-linking density and molar ratio of the networks should be satisfied [13]. However, these hydrogels also demonstrate inelastic features such as stress softening, necking, yielding and do not show self-healing after deformation, which is very critical for biological applications [14,15].
Polysaccharides are abundant materials in nature. Cellulose, which is the most predominant material found in nature, is also a class of polymeric polysaccharides. Most of the polysaccharides can be in the form of hydrogel due to their natural properties [16]. Alginate, agarose, chitosan, hyaluronic acid are some of the polysaccharides forming hydrogel. Due to their biological adaptation to the body as well as structural similarities to living tissues, polysaccharide hydrogels are important for biological implementations. In literature, numerous studies have been carried out about polysaccharide hydrogels for different fields [7,[17][18][19].
To overcome the mechanical weakness of polysaccharide hydrogels, toughening methods, some of which are mentioned above, are used [20]. Sun et al. [21] synthesized an alginate/PAAm DNH showing time and temperature-dependent self-healing after deformation. One of the most important inelastic features of natural tough hydrogels mimicking living tissues is self-healing. In this respect, the cross-linking procedure of tough hydrogels is critical. Generally, chemical procedures for self-healing tough hydrogels were based on changing the cross-linkages in the highly cross-linked polyelectrolyte network from covalent bonds that break in an irreversible way to weaker interactions such as ionic metal complexes which break in a reversible way. For example, the polyelectrolyte alginate network is formed by complex formation between calcium ions and the carboxylic acids of the alginate [21]. This alginate network was used in alginate/PAAm DNH with time and temperature-dependent self-healing behavior. Later, alginate/PAAm DNH is improved by using ion-sensitive segments to control the dimension and shape fidelity of the material [22]. To achieve better material response under deformation, different divalent and trivalent cations are utilized instead of Ca +2 as well [23][24][25]. The result of these studies shows that cations having a larger ion radius, such as Fe +3 , provide better stability in the biological environment and better mechanical features of the material.
Under large deformation DNH generally exhibits J-and S-shaped nonlinear stress behavior along with inelastic features similar to the rubber-like materials. To characterize inelastic features, such as the Mullins effect and hysteresis, background information for rubber-like materials is a useful tool. Constitutive models within the perspective of rubber elasticity are widely used to model complicated behaviors of rubber-like materials [26][27][28][29][30]. Rubber-like behaviors of DNH are also described using different concepts of the rubber elasticity. First models assumed the damage in DNH as a chain breakage [14,31]. Later, Zhao [32] decomposed the polymer matrix into short and long polymer chains and used Arruda-Boyce eight chain model [33] and the network alteration theory [34] to characterize the damage in DN hydrogel and the stress softening. Külcü [35] proposed a model for the stress softening and self-healing behavior of alginate-based DNH by using the network decomposition concept [29] and thermal effects in DNH put into a thermal chamber for selfhealing. Recently, Morovati and Dargazany represented a micro-mechanical constitutive model for the stress softening in DNH. Their model is based on the network evolution concept proposed by Dargazany and Itskov [27] and shows good agreement with experimental data. Mao et al. [36] constituted a viscoelasticity model for the complicated behavior of DNH under large deformation. Their model successfully describes the viscoelastic features as well as the stress softening in DN hydrogels.
In this contribution, a micro-mechanically based model is proposed to describe the enthalpic behavior of a polysaccharide polymer chain as well as an entropic one and the effect of various types of ionic cross-linkages in polysaccharide-based DNH. To constitute the present model, a new form of an energy function of a single polymer chain is introduced. Also, the Arruda-Boyce eight chain model [33], the network decomposition approach [27,29] and the network alteration concept [37] are pursued in the model.

Results and Discussion
The proposed model consists of eight material parameters for inelastic and hyperelastic polymer networks. Six material parameters belong to the polymer network exhibiting inelasticity, of which five of them are utilized for the entropic behavior while one parameter is taken for the enthalpic behavior. Two material parameters are accounted for the hyperelastic part of the model. Table 1 shows the material parameters along with their description. In the verification of the model with experimental data, the Levenberg-Marquardt algorithm has been used to find the values of the material parameters. Table 1. List of material parameters and their description.

Material Parameter Description
Inelastic (  Lastly, to depict the effect of the new contribution, the enthalpic part is removed from the model. Then, the sum of the stress contributions of the entropic and hyperelastic parts are compared with the proposed model (see Figure 3). The modeling results for Fe +3 , Al +3 and Ba +2 are utilized, as algiate/PAAm DNHs cross-linked via these counterions exhibit relatively more complicated behavior. The influence of the enthalpic part is clearly seen in hardening part of the material showing S-shaped behavior. It should be considered here that materials are subjected to tensile deformation, which is more than ten times of initial length. Although such huge deformation may result in a dramatic change in the material behavior and material structure, the proposed model is able to achieve a considerable match with the experimental data. The error margin of the model against experimental data is rather low. However, a yielding-like behavior is shown in the experimental data of DNHs cross-linked via trivalent cations during the deformation (see [23] and Figure 2d,e). To physically describe this behavior in the modeling, approaches of elastoplasticity might be useful to have a better characterization of this types of DNH.

Conclusions
In this study, a micromechanical constitutive model is proposed for polysaccharidebased DNH cross-linked via mono-, di-and trivalent cations. Due to the lengthening of the polysaccharide chain after a fully-stretched state, the free energy of a single chain is considered to involve enthalpic behavior as well as entropic one in the model. As a modeling tool, the Arruda-Boyce eight chain model [33], the network decomposition concept [38], the network alteration theory [37] and Hooke's law are used. First, the DNH matrix is decomposed into polyelectrolyte and neutral networks, where the latter one is assumed to be a purely elastic material. Then, the polyelectrolyte network is decomposed into entropic and enthalpic parts. For entropic chain behavior, the Langevin chain is used. To consider enthalpic behavior after the fully-stretched state, a straight polymer chain is assumed to be a cantilever and relatedly Hooke's law is utilized to describe stress contribution of the polymer chain due to its lengthening. Reorganization of the DNH matrix is taken into account by the network alteration theory. Finally, the model is compared against various DNH types. Good agreement between the model and experimental data is obtained.

Constitutive Modeling
A DNH consists of polyelectrolyte and neutral polymer networks [12]. While the latter is made up of loosely cross-linked long polymer chains, highly cross-linked polymer chains form the polyelectrolyte network. Alginate-PAAm hydrogel (APH) considered in the present manuscript is synthesized using the DNH concept [21,23]. Alginate of APH is a class of polysaccharides and acts as a polyelectrolyte network cross-linked via various types of cations in the DNH matrix. Polyacrylamide of APH is the neutral network and mainly ensures the stability of the polymer matrix. Within the perspective of the constitutive modeling, the strain energy (Ψ) of a DNH can be investigated by decomposing the DNH matrix into polyelectrolyte and neutral polymer networks as where superscripts • M , • P and • N denote the DNH matrix, polyelectrolyte and neutral polymer networks, respectively. In studies dealing with constitutive modeling, the energy of a polymer chain in microscale is mainly dominated by entropy [26,30]. The Langevin approximation is one of the most common approaches to model the entropic energy of a polymer chain [27,33]. However, in case of polysaccharides, during the deformation each monomer experiences a lengthening by about 10% [39] (see Figure 4). This fact necessitates a modification in the treatment of the modeling of a polymer chain, which is to consider the enthalpic behavior of a polymer chain in addition to the entropic one. A modified freely jointed chain (m-FJC) concept has been proposed to involve enthalpic and entropic energies simultaneously [40,41]. However, the derivation of m-FJC is mathematically complex and requires an approximation. This fact suspends the modeling approach from the physical motivation in the micro-scale. In the current paper, a micro-mechanically motivated model considering enthalpic and as well as entropic behaviors of a polymer chain is proposed. To this end, Ψ P is decomposed into entropic (• S ) and enthalpic (• H ) and Ψ M is rewritten as Assuming the polyacrylamide network is purely elastic [15], a hyperelastic model can be used to produce its mechanical behavior. Ψ N is thus given by [42] Ψ N = µ f (I 1 , α) where I 1 and I 2 are the first and second invariants of the right Cauchy-Green deformation tensor, respectively, µ is the shear modulus, α is the scalar material parameter and Introducing the number of polymer chains in the alginate network N P , Ψ P is expressed as To model the entropic behavior of the alginate network, the Arruda-Boyce eight chain model [33] is used due to its relatively simple and efficient mathematical expression. In the Arruda-Boyce model, a unit cube from a polymer network has eight chains along the half diagonals of the cube. Deformation of the unit cubes results in the stretching of the polymer chains in a network by the same ratio. The entropic energy of a single chain is captured by the Langevin statistics in the Arruda-Boyce eight chain model and can be expressed as where k is the Boltzmann constant, T is the absolute temperature, n is the segment number,r is the normalized average contour length and β is the inverse Langevin function.
As the stretch ratio is assumed to be the same for all chains, a fully stretched polymer chain can be considered as a cantilever to describe the enthalpic behavior (see Figure 5). Therefore, Hooke's law can be integrated into the enthalpic energy of a single chain as where E is the elastic modulus of the polymer matrix and taken from the experimental data [23], is the strain, which a polymer chain undergoes. Figure 6 shows the comparison between the entropic energy of a Langevin chain and the proposed energy function for ψ P . A probability density function is accounted for a network exhibiting inelastic behavior to define the network evolution under deformation (see Figure 7). To this end, a set of available chain D is firstly specified as where n min and n max denote minimum and maximum segment numbers in a network, respectively. n max is taken as a material parameter, where is n min is given by where ν is a material parameter utilized to prevent singularity at the fully stretched state of a Langevin chain. Also, note that the affine deformation is taken into account. Thus, macro-stretch Λ is equal to λ. D indicates the chain stretch and entropic behavior. However, fully stretched chain is deformed further due to deformation of polymer segments. Thus, a new set of available chain for entropic behavior D H is written as where n H min is the minimum segment number of a polymer chains showing enthalpic behavior. Considering 10% strain value for a fully-stretched polymer chain [39], n H min is represented as The Gaussian distribution gives a sufficient approximation if a polymer network is polymerized randomly. Assuming the damage in a network under loading takes place with respect to the chain length from the shortest chain to longer ones [29] and considering crosslinking density in the probability distribution [27,28,35], the probability density function (PDF) is written as, P(n,r, κ) = κ 6 where κ is the material parameter describing the cross-linking density of a polymer network and To normalize the PDF, the network alteration theory, which is to assume constant active segment number in a network, [37] is taken into account as Φ = 1 D P(n,r, κ)ndn .
The total energy of the alginate network is represented by Finally, the first Piola-Kirchhoff stress tensor P of the DNH matrix based on the incompressibility condition is given by where F is the deformation gradient.
Funding: This research received no external funding.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author, upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.