## Abstract

Polyelectrolyte (PE) gels consist of crosslinked polymer networks that are grafted with ionizable groups and ionic solution. Many stimuli-responsive gels, including pH-responsive, electric-responsive, and light-responsive ones, are PE gels. Most soft biological components are also PE gels. Due to the increasing scientific interests and applications of PE gels, a comprehensive model is needed. In PE gels, not only solvent, but also ions and other small molecules all diffuse inside, and the flows of the different components are coupled. This phenomenon is called cross-diffusion, meaning the flow of one species is not only driven by its own chemical potential gradient, but also influenced by the flow of other species. In this work, we develop a rigorous nonequilibrium thermodynamics framework to study the coupled deformation and diffusion of the PE gels where cross-diffusion is emphasized and quantified. Specific forms of free energy and kinetic laws are proposed. A finite element method is developed and implemented into abaqus through a user element subroutine. The model is used to simulate the deformation of biological axon and PE gels.The numerical results are compared with experimental data. It is shown that cross-diffusion generates anomalous effects not only on the flux but also on the deformation of PE gels.

## 1 Introduction

Gels are crosslinked polymer networks with solvent molecules that are imbibed into the network through predominantly entropic force. When the polymer chains are incorporated with ionizable groups, they become polyelectrolyte (PE) gels. Under certain conditions, these groups can be dissociated into two oppositely charged parts with one part attached to the main chain as fixed charges and the other part detached to become free ions in the solution inside the gel. These free ions disturb the local osmotic balance, which drives the solvent to flow in or out, leading to macroscopic swelling or shrinking of the gels. Many stimuli-responsive gels, including pH-responsive [1–3], electric-responsive [4,5], and light-responsive ones [6–8], are PE gels. They have been widely applied in many applications such as cell culture scaffold [9–11], drug delivery [12–14], microfluidics [15,16], sensors and actuators [17,18], soft robots [19], fuel cell membrane [20], loud speaker [21], artificial axon [22], wearable electronics [23], and many others. Beyond synthetic gels, most soft biological components are also PE gels [24–28], where proteins and polysaccharides aggregate into networks and expand themselves in a liquid environment loaded with various ions and macromolecules. These bio-gels play indispensable roles in every level of living organisms from cells, tissues, to organs by selectively allowing solvent and small molecules to diffuse in and out. As a result, nutrients are supplied, wastes are expelled, signals are transmitted, and various biological functions are realized.

In PE gels, not only solvent, but also ions and other neutral molecules all diffuse inside, and the flows of the different components are coupled. This phenomenon is called cross-diffusion, meaning the flow of one species is not only driven by its own chemical potential gradient but also influenced by the flow of other species. One important mechanism responsible for cross-diffusion is the frictional force between the different molecules [29,30], for which people used Maxwell–Stefan approach [29,31,32] to substitute the Fickian method—the flow of one species is proportional to the chemical potential gradient of all the species in different proportions. Other possible mechanisms of cross-diffusion include the electrostatic interactions, excluded volume effects, and other complex interactions between different mobile species and solid components [33]. In the existence of cross-diffusion, anomalous effects could happen. For example, when ionic solutions flow across charged porous membranes, the solvent could diffuse from low concentration solutions to high concentration solutions or the diffusion rate of certain species could decrease when the concentration difference of that species increases [34–37]. For another example, in ternary systems, the diffusivity could be zero, infinity, or negative corresponding to the osmotic diffusion, diffusion barrier, and reverse diffusion, respectively [29].

Due to the increasing scientific interests and applications of PE gels, a comprehensive model is needed. For biological tissues, Lai et al. [30] and Gu et al. [37,38] developed a biphasic model considering the material consists of a solid phase and a liquid phase that interact through frictional forces. The mixture theory was later extended to include the ions in the solution as the third phase, and they formulated the kinetic law as the flow of one species is driven by its chemical potential gradient and also dragged by the frictional force from interacting with other species. Although the theory has successfully modeled many different phenomena in biological tissues, it is phenomenological. The parameters in the model can not be directly related to the molecular properties of the material. Later based on the Flory-Huggins theory, Hong et al. [39] built a physics-based model for gels and extended it to include the fixed charges and ion effects for PE gels [40]. However, in the kinetics part of the modeling, only the flow of solvent is considered, and it is only driven by its own chemical potential gradient. The flow of ions and the cross-diffusion between solvent and solutes are ignored, which could be important for many gel applications.

In this paper, we follow the nonequilibrium thermodynamic approach for PE gels and formulate the kinetics in a more general form. The plan of this paper is organized as follows. Section 2 provides the nonequilibrium thermodynamics framework of the gel. Section 3 presents the kinetic law with cross-diffusion considered. In Sec. 4, a specific free energy function is given for PE gels and explicit forms of the constitutive and kinetic relations are provided. Finally, in Sec. 5, the computation results based on the new model are compared with the experimental results found in literature. The effects of cross-diffusion on gel behavior are also discussed.

## 2 Thermodynamic Framework

### 2.1 Kinematics.

*B*

_{0}(Fig. 1). Each material point is denoted by

**X**in

*B*

_{0}. After deformation, the material point

**X**moves to a new position with coordinate

**x**in the current configuration

*B*

_{t}. The current coordinate is described by the function

**x**= Γ(

**X**,

*t*). The deformation gradient

**F**is denoted as

### 2.2 Mass Conservation.

*C*

_{m}(

**X**,

*t*). Here, we formulate the case that there is only one type of fixed charged species in the network and use

*m*= 0 to denote it. we use

*m*= 1 to represent the solvent and

*m*> 1 for other soluble species. The mass change of each species in an arbitrary control volume $B0$ is attributed to the flux of that species through the boundary. The mass conservation of each species within the control volume $B0$ gives

**J**

_{m}is the nominal flux of species

*m*and

*m*

_{d}+ 1 is the total number of species in the system. Since the fixed charge is always attached to the polymer network,

**J**

_{0}= 0. Here and in the rest of the paper, we use Div and Grad to represent the divergence and gradient with respect to reference coordinates, and div and grad to represent the divergence and gradient with respect to the current coordinates.

### 2.3 Linear and Angular Momentum Balance.

**P**(

**X**,

*t*) is the first Piola-Kirchhoff stress tensor and

**B**(

**X**,

*t*) is the body force per reference volume. The Cauchy stress

**is related to**

*σ***P**by

**=**

*σ**J*

^{−1}

**PF**

^{T}. The angular momentum balance gives

### 2.4 Electrostatic Effects.

*a*)

*b*)

*ϕ*(

**X**,

*t*) is the electric potential. Here,

*Q*

_{e}is the net charge density per reference volume, which include both free ions and fixed charges. Therefore,

*z*

_{m}is the valence of

*m*th species and

*e*is the elementary charge.

### 2.5 First Law of Thermodynamics.

*U*(

**X**,

*t*) as the internal energy per unit reference volume. The right-hand side terms denote the different mechanisms to change the energy of the control volume. The first and second terms denote the mechanical energy input through body force

**B**(

**X**,

*t*) and traction

**T**(

**X**,

*t*) =

**P**·

**N**, respectively. Here,

**N**is the unit normal of referential surface $\u2202B0$. The third term accounts for the energy change due to the change of electric field. The fourth term is the energy carried by the chemical species transporting across the boundary with

*h*

_{m}being the enthalpy per particle of the

*m*th species. The expression works for any arbitrary volume $B0$. By using Eqs. (4), (6), and (7) and applying the divergence theorem, we can obtain the strong form of the first law of thermodynamics,

### 2.6 Second Law of Thermodynamics.

**X**,

*t*) is the entropy per unit reference volume of the dry gel. The left-hand side is the entropy variation of the control volume and the right-hand side is the entropy carried by the chemical species transporting across the boundary with

*θ*

_{m}being the entropy per particle of the

*m*th species.

### 2.7 Consequences of First and Second Laws.

*W*=

*U*−

*T*Θ, we obtain an inequality for Helmholtz free energy

*W*,

*m*th mobile species is defined as

*μ*

_{m}=

*h*

_{m}+

*ez*

_{m}

*ϕ*−

*Tθ*

_{m}. Accordingly, we postulate that the Helmholtz free energy is a function of deformation gradient

**F**, nominal electric displacement $D~$, and nominal concentration of different chemical species

*C*

_{m}, as state variables of the system,

*C*

_{1}is the nominal concentration of solvent.

**X**,

*t*) is the Lagrange multiplier that is introduced to enforce the incompressible constraint (Eq. (14)), and its physical meaning is the osmotic pressure. Here, we assume the solvent to be neutral, that is

*z*

_{1}= 0.

*a*)

*b*)

*c*)

*d*)

**f**

_{m}= −Grad

*μ*

_{m}is called thermodynamic force, while

**J**

_{m}is the corresponding thermodynamic flow [47]. If the system is not far from equilibrium, a linear relation between the flux and thermodynamics force can be written as

*L*

_{mn}must be a non-negative definite tensor. If we further assume the equality in the dissipative process (Eq. (19)) only holds in the equilibrium state, i.e.,

**f**

_{n}=

**0**, then

*L*

_{mn}must be positive definite. Furthermore, from Onsager reciprocal relations [48],

*L*

_{mn}is a symmetric tensor. Here,

*L*

_{mn}is a phenomenological parameter. Its physical meanings depend on specific kinetic models, which will be discussed in Sec. 3.

## 3 Kinetic Law

*L*

_{mn}being a diagonal matrix only and the chemical potential is only related to species concentration but not other driving forces. The physical meaning is that the flux of one species only dependents on its own concentration gradient. However, as has been evidenced by experimental observations, the flow in general should be driven by chemical potential gradient rather than concentration gradient, and importantly, the cross-diffusion often happens in multi-component solutions and mixtures, where the flow of one species is not only driven by its own chemical potential gradient, but also influenced by the flow of other species due to the interaction between different species. In order to model the cross-diffusion and meanwhile introduce the minimum number of nonzero parameters in

*L*

_{mn}matrix, we take the Maxwell–Stefan approach [38]. The governing equations for the flow of the chemical species are derived from the linear momentum conservation—the thermodynamic force, namely, the chemical potential gradient, of one species is balanced by the frictional forces between this species and all the other species when relative velocities are present between them. Since the flow in gels is usually slow, the inertia is neglected. The linear momentum conservation gives [38],

*a*)

*b*)

*s*to denote the solid polymer network, 1 to denote solvent,

*m*and

*n*to denote other solutes. The variable

*c*denotes the current concentration and is related to the reference concentration by the relation

*c*=

*C*/

*J*. The parameter

*f*

_{mn}is the frictional coefficients between components

*m*and

*n*, and

**v**

_{m}is the velocity of the

*m*th component.

*f*

_{mn}matrix is nonzero. The flow of any species can influence the flow of any other species. For simplicity, we consider the case that the concentrations of solutes are small compared to the solvent, which is also usually the case in most practically used PE gels. In this case, the frictional forces between different solutes and between solutes and polymer networks are negligible. Consequently, the nonzero friction coefficients are

*f*

_{1n}≠ 0,

*n*=

*s*, 2, …,

*m*

_{d}describing the friction between solvent and other species. If we consider the flux of each mobile species as the flow relative to the polymeric matrix, that is

**j**

_{m}=

*c*

_{m}(

**v**

_{m}−

**v**

_{s}), then we have

*a*)

*b*)

It is shown from Eq. (21*a*) that the flux of solvent depends not only on the chemical potential gradient of the solvent but also on the chemical potential gradients of all the solute species. The cross-diffusion comes from the frictional force between solutes and the solvent. Even when the solvent chemical potential is a constant in a certain region, if the chemical potential of any solute is nonzero, the flow of that solute will drag the solvent to flow with it. It is also found that the flux of solvent is proportional to the concentration of the moving solute species. Higher solute concentration provides bigger frictional force and higher flux of solvent. Comparing Eq. (21*b*) with Eq. (21*a*), we can see that the flux of a solute species is a combination of the convection with the solvent and the diffusion relative to the solvent. The cross diffusion between two different solute species is generated through the solvent flux.

*m*th solute as

*k*is the Boltzmann constant and

*T*is the temperature. These kinetic relations (Eq. (21)) can be rewritten in the reference configuration as

*a*)

*b*)

## 4 A Specific Material Model

### 4.1 Explicit Form of Free Energy Function for Polyelectrolyte Gels.

*W*

_{el}is due to the stretch of the polymer network. Here, we adopt the simplest chain model, the Neo-Hookean model, that is

*N*is the number of chains per reference volume, which is related to cross-link density.

*χ*

_{h}is the Flory–Huggins interaction parameter which represents the enthalpy of mixing between polymer and solvent, and $\mu 10$ is the chemical potential of solvent at a standard condition.

*m*th component at the standard condition.

*τ*. The amount of ionized group that has dissociated is denoted by

*C*

_{0}and the unionized amount is

*τ*−

*C*

_{0}. The free energy of dissociation consists of both the entropic and enthalpic parts,

The first term accounts for the entropy for mixing the dissociated groups and associated groups. The second and third terms together account for the enthalpy due to dissociation with $\mu f0$ and $\mu nf0$ being the chemical potential of dissociated groups and undissociated groups at the standard condition.

*ɛ*

_{e}is the permittivity of the gel.

### 4.2 Constitutive Relations.

*a*)

*b*)

*c*)

*d*)

**I**is the second-order identity tensor and ⊗ is the dyadic product.

### 4.3 Electrostatic Double Layer and Electroneutrality.

*r*

_{D}is the Debye length and is derived as

*m*th species in the external solution. For typical electrolyte gels, the value of

*r*

_{D}is in the order of nanometers. The Debye length is the characteristic length of the double layer near the interface where the electroneutrality is not satisfied. For most of the gels in practical applications, the size of the gel is much bigger than the size of the double layer. Therefore, both inside the gel and in the external solution far away from the boundary layer, the electroneutrality is satisfied. In the calculations for bulk gels, rather than solving the Poisson equation (Eq. (33)), the electric potential only needs to be solved by applying the electroneutrality condition as

## 5 Results and Discussion

To summarize the theory, the governing equations include the mass conservation equations (Eq. (3)), momentum conservation equations (Eqs. (4) and (5)), and Maxwell equations (Eq. (6)) or simply the electroneutrality equations (Eq. (35)). The field equations also include the constitutive relations (Eq. (31)) and kinetic relations (Eq. (24)). Together with mechanical and chemical boundary conditions, as well as the initial conditions, we can solve boundary value problems. To solve complicated geometry and inhomogeneous fields, we develop a user element subroutine in abaqus. In Sec. 5.1, we will first solve a 2D problem of axon swelling using finite difference method and compare it with the finite element result to verify the user element subroutine. The effect of cross-diffusion is discussed. We will then calibrate the model by calculating a 3D problem of a PE gel undergoing inhomogeneous cyclic swelling and shrinking and comparing the finite element results with the existing experimental results in literature. Finally, we discuss how the cross-diffusion would lead to abnormal deformations through two examples: 1D constrained swelling and transient bending.

### 5.1 Two-Dimensional Swelling and Verification of User Element Subroutine.

In the first example, we model the free swelling of an infinitely long neural axon. Typically, the length of an axon is more than 1000 times larger than the diameter of axon [53], so we regard it as a 2D plane strain problem. We compare the theoretical and finite element results to verify the user element code we developed in abaqus. Axonal cytoskeletons are composed of microtubules, neurofilaments, and actin filaments. They are weakly crosslinked networks [54] with charged groups. In biological processes, different species can diffuse in and out of the axon across the axon membrane both actively through ion channels and passively through diffusion driven by osmotic imbalance [55]. In this example, we focus on modeling the swelling and shrinking of the axon due to the passive diffusion. Specifically, the axon is first equilibrated with the surrounding solutions containing ions and neutral impermeable solutes (IS) such as proteins. When the axon is exposed to an external solution that contains the same amount of ions but without the impermeable solutes, an osmotic imbalance is developed between inside and outside of the axon, which drives the solvent to migrate into the axon (Fig. 2).

^{+}and K

^{+}, which play important roles in axonal signal transduction. The negative ion is Cl

^{−}. To simplify the problem, we only consider these three principle ions and omit all other ions. The concentrations of Na

^{+}, K

^{+}, and Cl

^{−}are kept at 10 mM, 100 mM, and 110 mM in the external solution. The polymer backbone of the axon is negatively charged [56]. The reference concentration of the fixed negative charges is taken to be

*C*

_{0}= 200 mM. Initially, the axon equilibrates in a solution that contains Na

^{+}, K

^{+}, Cl

^{−}, and an IS of concentration $C\xafIS=402.7mM$. The stretching ratios of the axon in the radial and angular directions are the same. Denoting

*λ*

_{0}as the initial stretching ratio, we have the deformation gradient as

*a*)) and exploiting the stress free condition, we have the force balance equation as

*a*) is the Maxwell stress and is neglected considering the macroscopic size of the axon.

*b*), that

^{+}, K

^{+}, and Cl

^{−}in the external solution. The right-hand side is the chemical potential of ions in the external solvent. Under the chemical equilibrium condition, the chemical potential of each ionic species inside the axon also equals the chemical potential outside. We have, from Eq. (31

*c*), that

*t*= 0

^{+}, the axon is transferred to a new solution that has the same ion concentration as in previous solution but without the ISs. Since the chemical potential of the external solution increases, the solvent will diffuse into the gel until a new equilibrium state is reached. For this axisymmetric problem, all the variables are functions of the reference radial coordinate

*R*only. The coordinate of the material element in the current configuration is

*r*(

*R*). The deformation gradient in the polar coordinate is

**P**is also a diagonal tensor where the radial and angular components are regarded as

*P*

_{R}and $P\Theta $. According to the constitutive relation (Eq. (31

*a*)), we have

*a*)

*b*)

*a*)

*b*)

*J*

_{m}represents the corresponding flux in the radial direction. The angular component of

**J**is zero. The kinetic relations are

*a*)

*b*)

Finally, the electroneutrality condition (Eq. (40)) is always satisfied.

The above set of equations is solved by using both the finite difference (FD) method and the finite element method (FEM) by using the abaqus uel subroutine we developed. In the calculation, the initial radius of the cylinder in the dry state is *R*_{0}, and the following material parameters are used: *N*Ω = 0.01, *χ*_{h} = 0.1, *D*_{1} = 5 × 10^{−8}m^{2/s}, and $DNa+=DK+=DCl\u2212=5\xd710\u221210m2/s$.

The stretching ratio of the axon is plotted as a function of the normalized time in Fig. 3. It can be seen that the FD and FEM results agree very well. In the short time period, the difference between the FD and FEM is relatively bigger, but still less than $1.43%$. The difference is mainly due to the different initial steps used in FD and FEM. In FD calculation, we start with a relatively larger time-step to avoid divergence, which is $D1\Delta t/R02=0.03$, while the initial step in FEM is $D1\Delta t/R02=10\u22126$.

Although the overall swelling of the axon increases monotonically, the local diffusion of the mobile species does not. To illustrate the effect of cross-diffusion, we plot the nominal concentrations of Na^{+} and Cl^{−} at the center of the axon in Fig. 4(a). If there is no cross-diffusion, the concentration of each mobile species should increase or decrease monotonically from its initial value to the equilibrium value. However, as shown in Fig. 4(a), the concentration $CNa+$ increases initially, but then decreases to the new equilibrium value. According to the chemical potential gradient of Na^{+}, it should have the overall tendency to diffuse into the gel. In this case, the solvent is also driven to diffuse into the gel and as the solvent migrates into the gel, it drags more Na^{+} ions to flow with it. As a result, the concentration of Na^{+} increases initially, but later, the ions have to flow out of the gel to compensate the overflow initially dragged by the solvent. Mathematically, it is the first term in Eq. (47*b*) that dominates the flow of Na^{+} in the first stage. As time goes, the second term in Eq. (47*b*) becomes dominant and the concentration of Na^{+} decreases until reaches the final equilibrium concentration. In Fig. 4(b), we also plot the evolution of the concentration of Cl^{−} ions in the center of the gel. Similar as Na^{+} ions, although the $CCl\u2212$ increases in general, there is an overshoot due to the water flux carrying Cl^{−} with it, which leads to the out flow of Cl^{−} in the later stage.

### 5.2 Three-Dimensional Swelling and Model Calibration.

Here, we use the experimental results in Ref. [57] to calibrate the model. In this experiment, the acrylamide (AAm) and methyl chloride quaternized *N*, *N*-dimethylamino ethylacrylate (DMAEA-Q) are copolymerized to form a hydrogel. In water, the DMAEA-Q monomers can be ionized to form fixed positive charges on the main chain and free ions Cl^{−}. In one set of experiment, the gel is put into pH 7 NaCl solution of different ionic strengths. Under this condition, the amount of hydrogen ions and hydroxide ions dissociated from water is much smaller than sodium ions and chloride ions. Therefore, in the simulation, we neglect the hydrogen and hydroxide ions and only consider Na^{+} and Cl^{−} as mobile ions.

*C*

_{0}is the concentration of the fixed positive charges on the main chain.

Solving Eqs. (49)–(52), we can obtain the equilibrium volume swelling ratio of $J0=\lambda 03$ as a function of ionic strengths of external solution. As shown in Fig. 5, we fit this theoretical curve with the experimental data. The gel swells in low ionic strength solution and shrinks in high ionic strength solution. The fitting parameters obtained through this process include *N*, *χ*_{h}, and *C*_{0}. Their values are listed in Table 1. From Fig. 5, we can see that the theoretical results fit well with the experimental data, with relative error no more than 8%. It is worth mentioning that in the experiment, the swelling ratio is measured with respect to the state of the gel as prepared, while in the theory, the swelling ratio is defined with respect to dry gel. In order to convert the two different reference states, we use the information provided in the Ref. [57] and calculate the water content of the gel as prepared to be 68.9%. So, the swelling ratio of the gel as prepared with respect to dry gel is 3.215. In Fig. 5, all the swelling ratios are calculated with respect to dry gel.

Value | Unit | Description | |
---|---|---|---|

N | 1.75 × 10^{25} | m^{−3} | Number of chains per reference volume |

Ω | 1 × 10^{−28} | m^{3} | Volume of one solvent molecule |

C_{0} | 9.67 × 10^{26} | m^{−3} | Number of fixed charge per reference volume |

χ_{h} | 0.206 | Dimensionless | Flory–Huggins interaction parameter |

D_{1} | 4.69 × 10^{−8} | m^{2}/s | Diffusivity of solvent |

$DNa+$ | 5 × 10^{−10} | m^{2}/s | Diffusivity of sodium ions |

$DCl\u2212$ | 5 × 10^{−10} | m^{2}/s | Diffusivity of chloride ions |

T | 300 | K | Absolute temperature |

Value | Unit | Description | |
---|---|---|---|

N | 1.75 × 10^{25} | m^{−3} | Number of chains per reference volume |

Ω | 1 × 10^{−28} | m^{3} | Volume of one solvent molecule |

C_{0} | 9.67 × 10^{26} | m^{−3} | Number of fixed charge per reference volume |

χ_{h} | 0.206 | Dimensionless | Flory–Huggins interaction parameter |

D_{1} | 4.69 × 10^{−8} | m^{2}/s | Diffusivity of solvent |

$DNa+$ | 5 × 10^{−10} | m^{2}/s | Diffusivity of sodium ions |

$DCl\u2212$ | 5 × 10^{−10} | m^{2}/s | Diffusivity of chloride ions |

T | 300 | K | Absolute temperature |

We then simulate the kinetic process and fit the simulation result with the experimental data to extract the kinetic properties of the PE gel. In experiment, the gel is periodically alternated between two pH 7 NaCl solutions: one with low ionic strength *I* = 0.05M and the other with high ionic strength *I* = 0.20M. The cyclic swelling and shrinking of the gel is recorded as a function of time. The gel undergoes inhomogeneous deformation with the edge swells (or shrinks) faster than the inner part of the gel. To simulate this process, we use FEM to solve the problem.

The simulation is done in abaqus through a user element subroutine. Continuum 3D brick elements with 20 nodes and reduced integration have been employed to discretize the model. The computational results are not sensitive to further discretization of the elements. For the boundary conditions, since the cylinder is symmetric about *X*, *Y*, and *Z* surface, only 1/8 of the cylinder is simulated. The chemical potentials of solvent and ions are kept as constants on the surface of the cylinder and are computed by the ionic strength of the external solution.

As shown in Fig. 6, the gel is originally in equilibrium with the external solution of ionic strength *I* = 0.05M. After it is put into the solution of higher ionic strength *I* = 0.2M, the gel shrinks with the edge shrinking faster than the inner part which develops an inhomogenous field of deformation. Overtime, the inner part of the gel also shrinks and the deformation of the gel becomes nearly homogeneous. Similarly, when the gel is then put back into the solution of lower ionic strength *I* = 0.05M again, the gel starts to swell with the edges swelling faster, but eventually every part of the gel swells almost equally. We calculate the overall volume change of the gel as a function of time as it alternates between the two solutions. We fit the numerical results with the experimental data. The result is shown in Fig. 7. Through this process, we obtain the additional material parameters to characterize the kinetic properties of the PE gel, including the diffusivity of solvent *D*_{1}, Na^{+} ions $DNa+$, and Cl^{−} ions $DCl\u2212$. The values are listed in Table 1.

### 5.3 One-Dimensional Swelling and Anomalous Displacement.

In previous examples, although the cross-diffusion happens for the ion transport, the overall swelling (or shrinking) of the gel or axon is still monotonic. In this section, we will show that in certain conditions, cross-diffusion may also induce abnormal swelling (or shrinking) of the gel.

*L*

_{0}. The gel is put into a container as shown in Fig. 8 and submerged in a NaCl solution of concentration $C\xafNa+=C\xafCl\u2212=C\xafion$. The gel swells against the container wall in the

*X*and

*Y*directions. After the gel reaches the wall, the stretching ratio of the gel along

*X*and

*Y*axes is

*λ*

_{f}= 2. This value is chosen so that the gel is always under compression and remains in contact with the

*X*and

*Y*surfaces of the container throughout all the following deformation processes even later when the gel is put into higher salt concentration solutions. After the gel reaches the container wall, it can only slide in the

*Z*-direction with the

*Z*= 0 surface fixed in space and the solvent and ions can only diffuse into or out of the gel through the positive

*Z*surface that is still in contact with the external solution. After the gel is in equilibrium with the initial NaCl solution of concentration $C\xafion$, the initial stretching ratio of the gel in the

*Z*direction

*λ*

_{z0}needs to be calculated. The deformation gradient of the gel is

The chemical equilibrium of solvent, the chemical equilibrium of the mobile ions, and the electroneutrality relation are the same as Eqs. (50), (51), and (52). Through these relations, the initial swelling ratio of the gel in the *Z*-direction *λ*_{z0} can be obtained.

*Z*surface of the gel. We calculate the transient swelling process of the gel under this condition. In this geometry, the

*X*and

*Y*coordinates of each material element remain the same, and only the position in

*Z*-direction changes, which is a function of

*Z*only and denoted as

*z*(

*Z*). The deformation gradient is

*Z*-direction equals zero:

*a*)

*b*)

*J*

_{m}represents the corresponding flux in the

*Z*-direction. The fluxes in the

*X*- and

*Y*- directions are zero. The kinetic relations are

*a*)

*b*)

The chemical potential of the solvent is the same as Eq. (43) and the chemical potential of the mobile ion species is the same as Eq. (44) with *m* = Na^{+}, Cl^{−}. Finally, the electroneutrality condition (Eq. (52)) is always satisfied. With the above-listed relations and the chemical boundary condition, that is, a constant value of chemical potential on the positive *Z* surface corresponding to the external solution of three times original ion concentration, and no flux on any other surface, we can fully solve the transient boundary value problem.

The parameters used in solving the problem are listed in the Table 1. Here, we vary the values of ion diffusivity and initial ion concentration. We calculate the transient swelling of the gel after the external ion concentration is raised to three times of the original values. Figure 9 plots the displacement of the free surface of the gel that is in contact with the external solution as a function of normalized time for two initial ion concentrations. Because of the sudden increase of the external salt concentration, water tends to flow out of the gel and ions tend to flow into the gel. As suggested by Eq. (58*a*), the direction of solvent flux *J*_{1} depends on three chemical potential gradients. As shown in Fig. 9(a), when $C\xafion$ is small, the solvent flux *J*_{1} is dominated by the solvent chemical potential gradient term d*μ*_{1}/d*Z* that drives the solvent to flow out of the gel. As a result, the gel shrinks through the whole kinetic process. Comparatively, if $C\xafion$ is large, the solvent flux *J*_{1} becomes dominated by the ion chemical potential gradient term d*μ*_{m}/d*Z*. The flow of ions drags the solvent to flow into the gel with them. As a result, the gel swells rather than shrinks in the beginning. Later, after the ion flow dies down, the value of d*μ*_{m}/d*Z* becomes smaller and the flux of solvent becomes dominant again by the solvent chemical potential gradient term that drives the solvent to flow out of the gel. Eventually, the gel shrinks. Therefore, for higher initial salt concentration of the external solution, it is more likely to develop non-monotonic swelling behavior of the gel.

Besides ion concentration, the diffusivities of ions and solvent are also important factors that determine whether the gel deforms monotonically or not. As shown in Fig. 9(b), if *D*_{ion} is slightly smaller than *D*_{1}, the chemical potential gradient of ions quickly diminishes. The flow of solvent is primarily dominated by its own chemical potential gradient and the gel shrinks monotonically. Comparatively, if *D*_{ion} is much smaller than *D*_{1}, the chemical potential gradient of ions remains large for a long period of time. The flow of solvent is influenced by both ion flows and its own chemical potential gradient. As a result, the gel shows non-monotonic deformation. In Fig. 10, we summarize the gel behavior whether it exhibits monotonic shrinking or non-monotonic deformation under different values of the initial salt concentration of the external solution and different values of the ratio between the ion diffusivity and solvent diffusivity in the form of a phase diagram. The triangles represent the cases when the gel swells first and then shrinks, while the circles represent the cases when the gel shrinks monotonically.

### 5.4 Bending With Abnormal Directions Due to Cross-Diffusion.

Beyond 1D problems, in this section, we will demonstrate the capability to simulate 3D geometries using our developed uel in abaqus. Since the beam bending geometry has been widely used in soft robotic grippers, here we simulate this problem. As shown in Fig. 11, We simulate a gel. Its dimensions in the dry state are *L*, *W*, and *H* in the *X*, *Y*, and *Z* directions respectively. The surfaces of the beam are covered with a soft impermeable film except the positive *Y*-surface. The ions and solvent can only migrate in and out of the gel through this surface. The same as in previous example, we use the calibrated material parameters as listed in Table 1 in the simulation. Here, we vary the diffusivity values and ion concentrations to explore the parametric dependence of cross-diffusion.

Initially, the gel is equilibrated in the 0.2M NaCl solution. At time *t* = 0^{+}, we transfer the gel to a 0.24M NaCl solution. In this stage, the chemical potential of solvent in the external solution is smaller than that inside the gel, while the chemical potential of Na^{+} ions and Cl^{−} ions in the external solution is higher than that inside the gel. Therefore, the solvent has the overall tendency to flow out of the gel and the ions have the tendency to flow into the gel. If there is no cross-diffusion, the gel beam will bend to the right, and then back to straight. The reason is that the chemical potential of solvent is lower in the external solution, so the solvent inside the gel tends to diffuse out of the gel with the solvent closer to the positive *Y* surface diffusing out first, which causes the gel beam to temporarily bend to the right, but as time goes, the solvent in the deep left side of the gel also gets enough time to migrate out of the gel, and eventually the gel shrinks the same everywhere, which causes the gel to bend back and become straight again. On the contrary, if cross-diffusion occurs, the beam may bend to the left first, then to the right, and eventually back to straight. The reason is the following. As illustrated in Fig. 11, in stage 1, although the solvent has an overall tendency to migrate out of the gel, the ions have a tendency to migrate into the gel, and as the ions flow into the gel, they drag solvent molecules with them into the gel. As a result, the right part of the gel swells, which causes the gel to bend to the left. In stage 2, as the solvent migrates into the gel, the chemical potential difference of the solvent inside and outside of the gel becomes even bigger. Eventually, the solvent has to come out. As the solvent comes out of the gel driven by the chemical potential gradient, the right part of the gel shrinks, which causes the gel beam to bend to the right. In final stage 3, the gel shrinks the same everywhere, which causes it to bend back to the left and become straight again. To quantify the bending, we denote the horizontal displacement of the left top corner of the beam as *u*_{2}. As shown in Fig. 12(a), the normalized displacement *u*_{2}/*W* is plotted as a function of time for several different values of the ion diffusivities. If the diffusivity of ions is small (i.e., the dash-dot curve for *D*_{ion}/*D*_{1} = 10^{−3} and dash curve *D*_{ion}/*D*_{1} = 10^{−2} in Fig. 12(a)), the gel bends to the left first and then to the right before it backs to straight. If the diffusivity of ions is large (i.e., the solid curve for *D*_{ion}/*D*_{1} = 10^{−1} in Fig. 12(a)), the chemical potential gradients of ions quickly diminish; so, the initial influx of solvent grabbed by ions is too small to cause the gel to bend to the left. Instead, the beam only bends to the right and then back to straight. For a similar reason, the gel with higher ion concentration first bends to the left then to the right and for the gel with lower ion concentration only bends to the right (Fig. 12(b)).

## 6 Conclusion

In this work, a large deformation theory coupled with a general kinetic law is developed for polyelectrolyte gels. A general nonequilibrium thermodynamic framework is presented. The different mobile species including ions, solvent, and neutral species are assumed to interact through frictional forces. A specific kinetic law is derived based on linear momentum balance of the mobile molecules. A specific free energy function is proposed, based on which specific forms of governing equations, constitutive relations and kinetic relations are derived. It is a fully coupled model of diffusion and deformation. On one hand, the deformation of the gel influences the flow of solvent and ions. On the other hand, the flow of solvent and ions also causes deformation. Due to cross-diffusion, the flow of solvent may induce the flow of ions and the flow of ions may also cause the flow of solvent. A finite element method is developed and implemented into abaqus through a user element subroutine. The model has been used to simulate the deformation of biological axon and PE gels. The numerical results are compared with experimental data. In the end, it is shown that the cross-diffusion will not only generate anomalous effects on the flux but also on the mechanical deformation.

## Acknowledgment

The materials are based upon work supported by the Air Force Office of Scientific Research under award no. FA9550-19-1-0395 (Dr. B.-L. “Les” Lee, Program Manager). Y.H. also acknowledges the funding support from National Science Foundation (Grant no. 1554326).

### Appendix A: Formulas in Current Configuration

*J*= 1 + Ω

*C*

_{1}by adding a large bulk modulus to the elastic energy.

*κ*is the bulk modulus of the gel and is much larger than the shear modulus

*NkT*. The constitutive relations are

*a*)

*b*)

*c*)

Here, the Maxwell stress in Eq. (31*a*) is ignored.

### Appendix B: Dimensionless Equations

*kT*/Ω, all lengths by the characteristic length of the system

*L*

_{0}, time by $L02/D1$, all concentrations by 1/Ω, electric potential by

*kT*/

*e*, all diffusivities by

*D*

_{1}, and all fluxes by

*D*

_{1}/Ω

*L*

_{0}and define dimensionless chemical potential as

*a*)

*b*)

*c*)

*d*)

*e*)

*f*)

*g*)

*h*)

### Appendix C: Numerical Solution Procedure

*S*

_{μ}and

*S*

_{j}are the corresponding subsurfaces of the boundary where prescribed chemical potential $\mu \xafm$ and prescribed flow $j\xafm$ are applied.

*S*

_{u}and

*S*

_{t}are the corresponding subsurfaces of the boundary where prescribed displacement $u\xaf$ and prescribed traction force $t\xaf$ are applied.

**w**and

*s*are applied to get the weak forms

*A*= 1, 2, … and

*D*= 1, 2, … denote the nodes of the element. $uiA$ and $\mu mD$ denote the nodal displacement and nodal value of chemical potentials.

*N*

^{A}and

*M*

^{D}are shape functions. As suggested in Ref. [58], different shape functions for displacement and chemical potentials are necessary to avoid numerical oscillation. Here, we use the second order shape function for displacement and first order shape function for chemical potentials.

**w**and

*s*.

*a*)

*b*)

*c*)

*t*denotes the values at time

*t*+ Δ

*t*.