Kinetics of Polyelectrolyte Gels

Zhang, Haohui; Dehghany, Mohammad; Hu, Yuhang

doi:10.1115/1.4046737

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.

Gels can absorb a large amount of solvent, especially polyelectrolyte gels. Their volume change through swelling and shrinking can be as large as 1000 times [41–44]. Here, we formulate a large deformation theory. We take the dry network as reference state and denote it as B₀ (Fig. 1). Each material point is denoted by X in B₀. 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

F = \frac{\partial x}{\partial X}

(1)

and the volume change is determined by

J = \det F

(2)

Fig. 1

View large Download slide

The reference configuration and current configuration of gels

2.2 Mass Conservation.

Gels contain solvent and different solute species. We denote the number of each species per unit reference volume as 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

B_{0}

is attributed to the flux of that species through the boundary. The mass conservation of each species within the control volume

B_{0}

gives

{\dot{C}}_{m} + Div J_{m} = 0, m = 0, 1, 2, \dots, m_{d}

(3)

where 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. 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.

Since the motion of the gel is usually slow, the inertia effect is neglected. The linear momentum balance gives

Div P^{T} + B = 0

(4)

where 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⁻¹PF^T. The angular momentum balance gives

P F^{T} = F P^{T}

(5)

2.4 Electrostatic Effects.

The PE gels contain a large amount of free ions and fixed charges. The Maxwell equations concerning electric field need to be considered, that is

Div \tilde{D} = Q_{e}

(6a)

\tilde{E} = - Grad ϕ

(6b)

where

\tilde{D} (X, t)

and

\tilde{E} (X, t)

are the nominal electric displacement and electric field, and ϕ(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,

Q_{e} = \sum_{m = 0}^{m_{d}} e z_{m} C_{m}

(7)

where z_m is the valence of mth species and e is the elementary charge.

2.5 First Law of Thermodynamics.

For an arbitrary control volume

B_{0}

of the material, energy can flow across the boundary

\partial B_{0}

in the form of mechanical and electrical works, and also through the flow of chemical species. The first law of thermodynamics requires the variation of the energy satisfies [45],

\begin{aligned} δ \int_{B_{0}} U d V_{0} & = \int_{B_{0}} B \cdot δ x d V_{0} + \int_{\partial B_{0}} T \cdot δ x d A_{0} - \int_{\partial B_{0}} ϕ δ \tilde{D} \cdot N d A_{0} \\ - \sum_{m = 1}^{m_{d}} \int_{\partial B_{0}} (h_{m} + e z_{m} ϕ) J_{m} δ t \cdot N d A_{0} \end{aligned}

(8)

Here, we assume the system undergoes isothermal processes. The left-hand side is the energy variation of a control volume of the material with 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

\partial B_{0}

⁠. 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 mth species. The expression works for any arbitrary volume

B_{0}

⁠. By using Eqs. , , and and applying the divergence theorem, we can obtain the strong form of the first law of thermodynamics,

\begin{aligned} δ U + \sum_{m = 0}^{m_{d}} e z_{m} ϕ δ C_{m} - P : δ F - \tilde{E} \cdot δ \tilde{D} \\ + Div (\sum_{m = 1}^{m_{d}} (h_{m} + e z_{m} ϕ) J_{m}) δ t = 0 \end{aligned}

(9)

2.6 Second Law of Thermodynamics.

The second law of thermodynamics requires that the entropy increment in a control volume should always equal or greater than the entropy imported into the control volume. The inequality relation can be expressed in the reference configuration as

δ \int_{B_{0}} Θ d V_{0} \geq - \sum_{m = 1}^{m_{d}} \int_{\partial B_{0}} θ_{m} J_{m} δ t \cdot N d A_{0}

(10)

where Θ (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 mth species.

Applying the divergence theorem, the differential form of the second law of thermodynamics becomes

δ Θ + Div (\sum_{m = 1}^{m_{d}} θ_{m} J_{m}) δ t \geq 0

(11)

2.7 Consequences of First and Second Laws.

Combining the first and second laws of thermodynamics (Eqs. and ) by introducing the Helmholtz free energy W = U − TΘ, we obtain an inequality for Helmholtz free energy W,

δ W - P : δ F - \tilde{E} \cdot δ \tilde{D} + \sum_{m = 1}^{m_{d}} (e z_{m} ϕ - μ_{m}) δ C_{m} + \sum_{m = 1}^{m_{d}} J_{m} \cdot Grad μ_{m} \leq 0

(12)

When deriving this inequality, the mass conservation equation is used. The chemical potential of the mth 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

\tilde{D}

⁠, and nominal concentration of different chemical species C_m, as state variables of the system,

W = W (F, \tilde{D}, C_{m})

(13)

Since most gels are soft and cannot sustain large stress, we can safely assume that the individual polymer chains and solvent molecules are incompressible. Additionally, the solute amount is usually much smaller than solvent, so their contribution to the total gel volume can be neglected. Therefore, the volume change of the gel is entirely due to the change of solvent concentration,

J = 1 + Ω C_{1}

(14)

where Ω is the volume of one solvent molecule and C₁ is the nominal concentration of solvent.

Substituting Eq. into Eq. and applying the constraint of Eq. , we have,

\begin{aligned} (\frac{\partial W}{\partial F} - P - Π J F^{- T}) : δ F + (\frac{\partial W}{\partial \tilde{D}} - \tilde{E}) \cdot δ \tilde{D} \\ + (\frac{\partial W}{\partial C_{1}} - μ_{1} + Π Ω) δ C_{1} + \sum_{m = 2}^{m_{d}} (\frac{\partial W}{\partial C_{m}} - μ_{m} + e z_{m} ϕ) δ C_{m} \\ + \sum_{m = 1}^{m_{d}} J_{m} δ t \cdot Grad μ_{m} \leq 0 \end{aligned}

(15)

where Π(X, t) is the Lagrange multiplier that is introduced to enforce the incompressible constraint (Eq. ), and its physical meaning is the osmotic pressure. Here, we assume the solvent to be neutral, that is z₁ = 0.

As the Coleman–Noll procedure [46] suggests, Eq. must be valid for all possible admissible thermodynamic processes, so the coefficient before the variation of each state variable must equal to zero, which derives the following constitutive relations:

P = \frac{\partial W}{\partial F} - Π J F^{- T}

(16a)

μ_{1} = \frac{\partial W}{\partial C_{1}} + Π Ω

(16b)

μ_{m} = \frac{\partial W}{\partial C_{m}} + e z_{m} ϕ, m = 2, 3, \dots, m_{d}

(16c)

\tilde{E} = \frac{\partial W}{\partial \tilde{D}}

(16d)

and the inequality term that defines the dissipative processes becomes

\sum_{m = 1}^{m_{d}} J_{m} \cdot Grad μ_{m} \leq 0

(17)

Here, the dissipative process is only attributed to diffusion. The gradient term 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

J_{m} = \sum_{n = 1}^{m_{d}} L_{m n} f_{n}, m = 1, 2, 3, \dots, m_{d}

(18)

Substituting this relation (Eq. ) into the inequality that defines the dissipation process (Eq. ), we have

\sum_{m = 1}^{m_{d}} \sum_{n = 1}^{m_{d}} L_{m n} f_{m} \cdot f_{n} \geq 0

(19)

To satisfy this non-negative condition, L_mn must be a non-negative definite tensor. If we further assume the equality in the dissipative process (Eq. ) 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

For the diffusion of chemical species, the simplest model is Fick's first law, which is a special case for the relation in Eq. with 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],

c_{1} grad μ_{1} = f_{s 1} (v_{s} - v_{1}) + \sum_{n = 2}^{m_{d}} f_{1 n} (v_{n} - v_{1})

(20a)

c_{m} grad μ_{m} = f_{s m} (v_{s} - v_{m}) + \sum_{n = 1, n \neq m}^{m_{d}} f_{m n} (v_{n} - v_{m}), m = 2, 3, \dots, m_{d} (20 b)

(20b)

Here, we use 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 mth component.

In the most general case, every component in the 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

j_{1} = - \frac{c_{1}}{f_{1 s}} \sum_{n = 1}^{m_{d}} c_{n} grad μ_{n}

(21a)

j_{m} = - \frac{c_{m}}{f_{1 s}} \sum_{n = 1}^{m_{d}} (c_{n} grad μ_{n}) - \frac{c_{m}^{2}}{f_{1 m}} grad μ_{m}, m = 2, 3, \dots, m_{d}

(21b)

It is shown from Eq. (21a) 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. (21b) with Eq. (21a), 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.

We define the diffusivity of solvent as,

D_{1} = \frac{k T c_{1}}{f_{1 s}}

(22)

and the diffusivity of the mth solute as

D_{m} = \frac{k T c_{m}}{f_{1 m}}, m = 2, 3, \dots, m_{d}

(23)

where k is the Boltzmann constant and T is the temperature. These kinetic relations (Eq. ) can be rewritten in the reference configuration as

J_{1} = - \frac{D_{1}}{k T} \sum_{n = 1}^{m_{d}} C_{n} Grad μ_{n}

(24a)

J_{m} = - \frac{D_{1} C_{m}}{k T C_{1}} \sum_{n = 1}^{m_{d}} (C_{n} Grad μ_{n}) - \frac{D_{m} C_{m}}{k T} Grad μ_{m}, m = 2, 3, \dots, m_{d} (24 b)

(24b)

4 A Specific Material Model

4.1 Explicit Form of Free Energy Function for Polyelectrolyte Gels.

Now we introduce a specific free energy function for PE gels. Following Flory-Huggin's theory [49–52], we express the free energy of the gel as a summation of the stretching energy due to stretching of the polymer chains and the mixing energy due to the mixing among polymer, solvent and solutes. In addition, the dissociation of the ionizable groups and electric field inside the gel also contribute to the total free energy [40,45]. The total free energy form is

W = W_{el} + W_{mix} + W_{ion} + W_{fix} + W_{pol}

(25)

The elastic energy W_el is due to the stretch of the polymer network. Here, we adopt the simplest chain model, the Neo-Hookean model, that is

W_{el} = \frac{N k T}{2} (F : F - 3 - 2 \ln J)

(26)

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

The mixing energy of polymer and solvent takes the form

W_{mix} = k T C_{1} (\ln (\frac{Ω C_{1}}{1 + Ω C_{1}}) + \frac{χ_{h}}{1 + Ω C_{1}}) + C_{1} μ_{1}^{0}

(27)

where χ_h is the Flory–Huggins interaction parameter which represents the enthalpy of mixing between polymer and solvent, and

μ_{1}^{0}

is the chemical potential of solvent at a standard condition.

The concentration of mobile solutes is typically low, so the mixing between solvent and solutes is mostly entropic,

W_{ion} = k T \sum_{m = 2}^{m_{d}} C_{m} [\ln (\frac{C_{m}}{J C_{m}^{ref}}) - 1] + \sum_{m = 2}^{m_{d}} C_{m} μ_{m}^{0}

(28)

where

C_{m}^{ref}

is the concentration of each solute species in a reference condition and

μ_{m}^{0}

is the chemical potential of the mth component at the standard condition.

The ionizable groups on the polymer chains can dissociate and become fixed charges on the main chains and the mobile ions dissolved in the solution. The degree of dissociation depends on the surrounding chemical conditions. The concentration of the total ionizable groups in the reference configuration is a constant τ. The amount of ionized group that has dissociated is denoted by C₀ and the unionized amount is τ − C₀. The free energy of dissociation consists of both the entropic and enthalpic parts,

\begin{aligned} W_{dis} & = k T [C_{0} \ln (\frac{C_{0}}{τ}) + (τ - C_{0}) \ln (\frac{τ - C_{0}}{τ})] + μ_{f}^{0} C_{0} \\ + μ_{n f}^{0} (τ - C_{0}) \end{aligned}

(29)

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 $μ_{f}^{0}$ and $μ_{n f}^{0}$ being the chemical potential of dissociated groups and undissociated groups at the standard condition.

Finally, the polarization energy is

W_{pol} = \frac{1}{2 ε_{e} J} (F \cdot \tilde{D}) \cdot (F \cdot \tilde{D})

(30)

where ɛ_e is the permittivity of the gel.

4.2 Constitutive Relations.

Substituting the free energy function (Eq. ) into Eq. , we can get the explicit forms of the constitutive relations of the PE gels,

\begin{aligned} P & = N k T (F - F^{- T}) - k T (\sum_{m = 2}^{m_{d}} C_{m}) F^{- T} - J Π F^{- T} \\ + \frac{1}{2 ε_{e} J} (2 (F \cdot \tilde{D}) \otimes (F \cdot \tilde{D}) - (F \cdot \tilde{D}) \cdot (F \cdot \tilde{D}) I) \cdot F^{- T} \end{aligned}

(31a)

μ_{1} = μ_{1}^{0} + k T (\ln (1 - \frac{1}{J}) + \frac{1}{J} + \frac{χ_{h}}{J^{2}}) + Π Ω

(31b)

μ_{m} = μ_{m}^{0} + k T \ln \frac{C_{m}}{J C_{m}^{ref}} + e z_{m} ϕ, m = 2, 3, \dots, m_{d}

(31c)

\tilde{E} = \frac{1}{ε_{e} J} F^{T} \cdot F \cdot \tilde{D}

(31d)

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

4.3 Electrostatic Double Layer and Electroneutrality.

The combination of Maxwell equations (Eq. ) and the constitutive relation (Eq. ) gives the Poisson equation of electric potential,

Div (J ε_{e} F^{- T} \cdot F^{- 1} \cdot Grad ϕ) = - \sum_{m = 0}^{m_{d}} e z_{m} C_{m}

(32)

This equation can be nondimensionalized as

r_{D}^{2} Div (J F^{- T} \cdot F^{- 1} \cdot Grad \frac{e ϕ}{k T}) = - \frac{\sum_{m = 0}^{m_{d}} z_{m} C_{m}}{\sum_{m = 1}^{m_{d}} z_{m}^{2} {\bar{C}}_{m}}

(33)

where r_D is the Debye length and is derived as

r_{D} = \sqrt{\frac{ε_{e} k T}{2 e^{2} I}}

(34)

where

I = \frac{1}{2} \sum_{m = 2}^{m_{d}} z_{m}^{2} {\bar{C}}_{m}

is the ionic strength with

{\bar{C}}_{m}

being the concentration of the mth 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. ), the electric potential only needs to be solved by applying the electroneutrality condition as

\sum_{m = 0}^{m_{d}} e z_{m} C_{m} = 0

(35)

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).

Fig. 2

View large Download slide

Schematic of the axon swelling process. An infinitely long axon is equilibrated with the external solution that contains Na⁺, K⁺, Cl⁻ and impermeable solute (IS). When the axon is changed to a solution contains only Na⁺, K⁺, and Cl⁻ of the same concentration, but without IS, the axon swells.

The system contains Na⁺ 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₀ = 200 mM. Initially, the axon equilibrates in a solution that contains Na⁺, K⁺, Cl⁻, and an IS of concentration

{\bar{C}}_{IS} = 402.7 mM

⁠. The stretching ratios of the axon in the radial and angular directions are the same. Denoting λ₀ as the initial stretching ratio, we have the deformation gradient as

F = [\begin{matrix} λ_{0} & 0 & 0 \\ 0 & λ_{0} & 0 \\ 0 & 0 & 1 \end{matrix}], J_{0} = λ_{0}^{2}

(36)

Substituting this deformation gradient into the constitutive relation (Eq. ) and exploiting the stress free condition, we have the force balance equation as

N k T (λ_{0}^{2} - 1) - k T (C_{{Na}^{+}} + C_{K^{+}} + C_{{Cl}^{-}}) - Π λ_{0}^{2} = 0

(37)

Here, the last term in Eq. is the Maxwell stress and is neglected considering the macroscopic size of the axon.

Not only the axon is in force equilibrium, but also chemical equilibrium. The chemical potential of the solvent inside the axon equals that of the external solvent. Setting the chemical potential of pure solvent in the standard state as zero, we have, from Eq. , that

\begin{aligned} k T (\ln (1 - \frac{1}{λ_{0}^{2}}) + \frac{1}{λ_{0}^{2}} + \frac{χ_{h}}{λ_{0}^{4}}) + Π Ω = & - k T Ω ({\bar{C}}_{{Na}^{+}} + {\bar{C}}_{K^{+}} \\ + {\bar{C}}_{{Cl}^{-}} + {\bar{C}}_{IS}) \end{aligned}

(38)

where

{\bar{C}}_{{Na}^{+}}, {\bar{C}}_{K^{+}},

and

{\bar{C}}_{{Cl}^{-}}

denote the concentration of Na⁺, 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. , that

k T \ln \frac{C_{N a^{+}}}{λ_{0}^{2}} + e ϕ = k T \ln {\bar{C}}_{N a^{+}}, k T \ln \frac{C_{K^{+}}}{λ_{0}^{2}} + e ϕ = k T \ln {\bar{C}}_{K^{+}}, k T \ln \frac{C_{C l^{-}}}{λ_{0}^{2}} - e ϕ = k T \ln {\bar{C}}_{C l^{-}}

(39)

where the right-hand side is the chemical potential of the ions in the external solution.

The electroneutrality requires

C_{{Na}^{+}} + C_{K^{+}} - C_{{Cl}^{-}} - C_{0} = 0

(40)

Solving Eqs. (37)–(40), we can get the initial stretching ratio of the axon, that is λ₀ = 2.0.

We next solve the transient swelling of the axon. At 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

F = [\begin{matrix} \frac{d r}{d R} & 0 & 0 \\ 0 & \frac{r}{R} & 0 \\ 0 & 0 & 1 \end{matrix}], J = \frac{d r}{d R} \frac{r}{R}

(41)

The stress P is also a diagonal tensor where the radial and angular components are regarded as P_R and

P_{Θ}

⁠. According to the constitutive relation (Eq. ), we have

P_{R} = N k T (\frac{d r}{d R} - \frac{d R}{d r}) - k T (C_{N a^{+}} + C_{K^{+}} + C_{C l^{-}}) \frac{d R}{d r} - Π \frac{r}{R}

(42a)

P_{Θ} = N k T (\frac{r}{R} - \frac{R}{r}) - k T (C_{{Na}^{+}} + C_{K^{+}} + C_{{Cl}^{-}}) \frac{R}{r} - Π \frac{d r}{d R}

(42b)

The chemical potential of the solvent in Eq. can be expressed as

μ_{1} = k T (\ln (1 - \frac{1}{J}) + \frac{1}{J} + \frac{χ_{h}}{J^{2}}) + Π Ω

(43)

and the chemical potentials of the mobile ion species in Eq. become

μ_{m} = μ_{m}^{0} + k T \ln \frac{C_{m}}{J C_{m}^{ref}} + e z_{m} ϕ, m = N a^{+}, K^{+}, C l^{-}

(44)

The evolutionary process satisfies the linear momentum equation,

\frac{d P_{R}}{d R} + \frac{P_{R} - P_{Θ}}{R} = 0

(45)

and mass conservation equations,

\frac{d C_{1}}{d t} + \frac{d J_{1}}{d R} + \frac{J_{1}}{R} = 0

(46a)

\frac{d C_{m}}{d t} + \frac{d J_{m}}{d R} + \frac{J_{m}}{R} = 0, m = {Na}^{+}, K^{+}, {Cl}^{-}

(46b)

Here, J_m represents the corresponding flux in the radial direction. The angular component of J is zero. The kinetic relations are

J_{1} = - \frac{D_{1}}{k T} (C_{1} \frac{d μ_{1}}{d R} + C_{{Na}^{+}} \frac{d μ_{N a^{+}}}{d R} + C_{K^{+}} \frac{d μ_{K^{+}}}{d R} + C_{C l^{-}} \frac{d μ_{C l^{-}}}{d R})

(47a)

J_{m} = \frac{C_{m}}{C_{1}} J_{1} - \frac{D_{m} C_{m}}{k T} \frac{d μ_{m}}{d R}, m = N a^{+}, K^{+}, C l^{-}

(47b)

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₀, and the following material parameters are used: NΩ = 0.01, χ_h = 0.1, D₁ = 5 × 10⁻⁸m^2/s, and $D_{{Na}^{+}} = D_{K^{+}} = D_{{Cl}^{-}} = 5 \times 10^{- 10} m^{2} / 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 $D_{1} Δ t / R_{0}^{2} = 0.03$ ⁠, while the initial step in FEM is $D_{1} Δ t / R_{0}^{2} = 10^{- 6}$ ⁠.

Fig. 3

View large Download slide

The stretching ratio of the axon changes as a function of nondimensionalized time $D_{1} t / R_{0}^{2}$ ⁠. The line is the result from finite difference method and the dots are the result from finite element method by using a user element subroutine.

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 $C_{{Na}^{+}}$ 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. (47b) that dominates the flow of Na⁺ in the first stage. As time goes, the second term in Eq. (47b) 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 $C_{{Cl}^{-}}$ 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.

Fig. 4

View large Download slide

The concentration of the (a) Na⁺ and (b) Cl⁻ in the center of the axon is plotted as a function of the normalized time $D_{1} t / R_{0}^{2}$

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.

In this example, we simulate a 3D cylinder with a diameter of 5 mm and a length of 10 mm as prepared, which is the same as in experiment. The first step is to compare the equilibrium swelling ratio at different ionic strengths at pH = 7. When it reaches equilibrium, the swelling ratio is the same in all three directions, so the deformation gradient is

F = [\begin{matrix} λ_{0} & 0 & 0 \\ 0 & λ_{0} & 0 \\ 0 & 0 & λ_{0} \end{matrix}], J_{0} = λ_{0}^{3}

(48)

Substituting Eq. into Eq. , we have the force balance equation in stress free state as

N k T (λ_{0}^{2} - 1) - k T (C_{{Na}^{+}} + C_{{Cl}^{-}}) - Π J_{0} = 0

(49)

the chemical equilibrium of solvent as

k T (\ln (1 - \frac{1}{J_{0}}) + \frac{1}{J_{0}} + \frac{χ_{h}}{J_{0}^{2}}) + Π Ω = - k T Ω ({\bar{C}}_{{Na}^{+}} + {\bar{C}}_{{Cl}^{-}})

(50)

and the chemical equilibrium of the mobile ions as

k T \ln \frac{C_{N a^{+}}}{J_{0}} + e ϕ = k T \ln {\bar{C}}_{N a^{+}}, k T \ln \frac{C_{C l^{-}}}{J_{0}} - e ϕ = k T \ln {\bar{C}}_{C l^{-}}

(51)

where

{\bar{C}}_{{Na}^{+}}

and

{\bar{C}}_{{Cl}^{-}}

are the sodium ion and chloride ion concentrations in the external solution. The right-hand sides of Eqs. and are the chemical potentials of ions in the external solution. Besides, we also have the electroneutrality relation (Eq. ), which can be expressed as

C_{{Na}^{+}} - C_{{Cl}^{-}} + C_{0} = 0

(52)

where C₀ 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 $J_{0} = λ_{0}^{3}$ 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₀. 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.

Fig. 5

View large Download slide

The experimental and theoretical equilibrium swelling ratio of AAm-DMAEA-Q gel at different ionic strengths

Table 1

Parameters calibrated from fitting with the experimental results

	Value	Unit	Description
N	1.75 × 10²⁵	m⁻³	Number of chains per reference volume
Ω	1 × 10⁻²⁸	m³	Volume of one solvent molecule
C₀	9.67 × 10²⁶	m⁻³	Number of fixed charge per reference volume
χ_h	0.206	Dimensionless	Flory–Huggins interaction parameter
D₁	4.69 × 10⁻⁸	m²/s	Diffusivity of solvent
$D_{N a^{+}}$	5 × 10⁻¹⁰	m²/s	Diffusivity of sodium ions
$D_{C l^{-}}$	5 × 10⁻¹⁰	m²/s	Diffusivity of chloride ions
T	300	K	Absolute temperature

	Value	Unit	Description
N	1.75 × 10²⁵	m⁻³	Number of chains per reference volume
Ω	1 × 10⁻²⁸	m³	Volume of one solvent molecule
C₀	9.67 × 10²⁶	m⁻³	Number of fixed charge per reference volume
χ_h	0.206	Dimensionless	Flory–Huggins interaction parameter
D₁	4.69 × 10⁻⁸	m²/s	Diffusivity of solvent
$D_{N a^{+}}$	5 × 10⁻¹⁰	m²/s	Diffusivity of sodium ions
$D_{C l^{-}}$	5 × 10⁻¹⁰	m²/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₁, Na⁺ ions $D_{N a^{+}}$ ⁠, and Cl⁻ ions $D_{C l^{-}}$ ⁠. The values are listed in Table 1.

Fig. 6

View large Download slide

Contours of the swelling ratio on the cross section of the cylinder during the shrinking and swelling. (a) At t = 0 h, the gel is in equilibrium with the ionic solution with I = 0.05M. (b) At t =3.97 h, the gel is shrinking after the gel is transferred to the solution with I = 0.2M. The corners shrink faster and generate inhomogeneous deformation. (c) At t = 23.85 h, the gel shrinks everywhere, almost reaching equilibrium state. After this point, the gel is put back to the solution with I = 0.05M. (d) At t = 24.25 h. The gel is swelling with corners and edges swelling faster. (e) At t = 47.7 h, the gel swells everywhere, almost reaching equilibrium state. It is right before the gel is put to I = 0.2M solution again.

Fig. 7

View large Download slide

The cyclic volume change of the gel when it is repetitively transferred between solutions of high and low ionic strengths

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.

Here, we formulate a 1D problem. We start with a cubic dry gel, the length of which is L₀. The gel is put into a container as shown in Fig. 8 and submerged in a NaCl solution of concentration

{\bar{C}}_{N a^{+}} = {\bar{C}}_{C l^{-}} = {\bar{C}}_{ion}

⁠. 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

{\bar{C}}_{ion}

⁠, the initial stretching ratio of the gel in the Z direction λ_z0 needs to be calculated. The deformation gradient of the gel is

F = [\begin{matrix} λ_{f} & 0 & 0 \\ 0 & λ_{f} & 0 \\ 0 & 0 & λ_{z 0} \end{matrix}], J_{0} = λ_{f}^{2} λ_{z 0}

(53)

Fig. 8

View large Download slide

The schematic of the 1D gel. The gel deformation is confined in X and Y directions and is only allowed to move in the Z-direction. The solvent is only allowed to penetrate through the positive Z face that is exposed to external solution.

Substituting Eq. into Eq. , and setting the stress in the Z-direction to be zero, we have

N k T (λ_{z 0}^{2} - 1) - k T (C_{N a^{+}} + C_{C l^{-}}) - J_{0} Π = 0

(54)

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.

The gel together with the container is then put into a new NaCl solution, the ion concentration of which is three times of the initial concentration

{\bar{C}}_{ion}

⁠. Due to the osmotic change, ions and solvent start to flow into or out of the gel through the positive 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

F = [\begin{matrix} λ_{f} & 0 & 0 \\ 0 & λ_{f} & 0 \\ 0 & 0 & \frac{d z}{d Z} \end{matrix}], J = λ_{f}^{2} \frac{d z}{d Z}

(55)

The stress in the Z-direction equals zero:

N k T (\frac{d z}{d Z} - \frac{d Z}{d z}) - k T (C_{N a^{+}} + C_{C l^{-}}) \frac{d Z}{d z} - Π λ_{f}^{2} = 0

(56)

The mass conservation equations become

\frac{d C_{1}}{d t} + \frac{d J_{1}}{d Z} = 0

(57a)

\frac{d C_{m}}{d t} + \frac{d J_{m}}{d Z} = 0, m = {Na}^{+}, {Cl}^{-}

(57b)

where J_m represents the corresponding flux in the Z-direction. The fluxes in the X- and Y- directions are zero. The kinetic relations are

J_{1} = - \frac{D_{1}}{k T} (C_{1} \frac{d μ_{1}}{d Z} + C_{N a^{+}} \frac{d μ_{N a^{+}}}{d Z} + C_{C l^{-}} \frac{d μ_{C l^{-}}}{d Z})

(58a)

J_{m} = \frac{C_{m}}{C_{1}} J_{1} - \frac{D_{m} C_{m}}{k T} \frac{d μ_{m}}{d Z}, m = {Na}^{+}, {Cl}^{-}

(58b)

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. (58a), the direction of solvent flux J₁ depends on three chemical potential gradients. As shown in Fig. 9(a), when ${\bar{C}}_{ion}$ is small, the solvent flux J₁ is dominated by the solvent chemical potential gradient term dμ₁/dZ that drives the solvent to flow out of the gel. As a result, the gel shrinks through the whole kinetic process. Comparatively, if ${\bar{C}}_{ion}$ is large, the solvent flux J₁ becomes dominated by the ion chemical potential gradient term dμ_m/dZ. 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/dZ 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.

Fig. 9

(a) The displacement of the free end of the gel for two different values of initial ion concentrations of external solution. In this case, Dion = 10−2D1. (b) The displacement of the free end of the gel for two different values of Dion/D1. In this case, C¯ion=0.1M. In both plots, the transient behavior is calculated when the ion concentration of external solution is increased to three times of its initial value.

View large Download slide

(a) The displacement of the free end of the gel for two different values of initial ion concentrations of external solution. In this case, D_ion = 10⁻²D₁. (b) The displacement of the free end of the gel for two different values of D_ion/D₁. In this case, ${\bar{C}}_{ion} = 0.1 M$ ⁠. In both plots, the transient behavior is calculated when the ion concentration of external solution is increased to three times of its initial value.

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₁, 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₁, 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.

Fig. 10

The gel behavior for different Dion to D1 ratio and the initial ion concentration. The triangle denotes the case that the gel first swells and then shrinks. The circle denotes the case that the gel only shrinks.

View large Download slide

The gel behavior for different D_ion to D₁ ratio and the initial ion concentration. The triangle denotes the case that the gel first swells and then shrinks. The circle denotes the case that the gel only shrinks.

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.

Fig. 11

View large Download slide

The schematic of the gel bending to the left or to the right and finally reaching equilibrium. The initial chemical potential of solvent is lower in the external solution in the beginning. The external chemical potentials of Na⁺ and Cl⁻ are higher in the beginning. The right surface is permeable to all species while other surfaces are impermeable.

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₂. As shown in Fig. 12(a), the normalized displacement u₂/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₁ = 10⁻³ and dash curve D_ion/D₁ = 10⁻² 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₁ = 10⁻¹ 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)).

Fig. 12

(a) The Y-displacement of the gel with different ion diffusivities when C¯ion=0.2M and (b) the Y-displacement of the gel with different ion concentrations when Dion/D1 = 10−3

View large Download slide

(a) The Y-displacement of the gel with different ion diffusivities when ${\bar{C}}_{ion} = 0.2 M$ and (b) the Y-displacement of the gel with different ion concentrations when D_ion/D₁ = 10⁻³

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

In order to implement the theory into abaqus uel subroutine, all equations are transformed to the current configuration. The mass conservation equations become

{\dot{c}}_{m} + c_{m} \frac{\dot{J}}{J} + div j_{m} = 0, m = 1, 2, \dots, m_{d}

(A1)

The linear momentum equation is

div σ = 0

(A2)

The electroneutrality relation is

\sum_{m = 0}^{m_{d}} e z_{m} c_{m} = 0

(A3)

The kinetics laws become

j_{1} = - \frac{D_{1}}{k T} \sum_{n = 1}^{m_{d}} c_{n} grad μ_{n}

(A4)

j_{m} = - \frac{D_{1} c_{m}}{k T c_{1}} \sum_{n = 1}^{m_{d}} c_{n} grad μ_{n} - \frac{D_{m} c_{m}}{k T} grad μ_{m}, m = 2, 3, \dots, m_{d}

(A5)

In order to avoid divergence due the incompressibility condition, we replace the incompressibility condition J = 1 + ΩC₁ by adding a large bulk modulus to the elastic energy.

W_{el} = \frac{N k T}{2} (F : F - 3 - 2 \ln J) + \frac{κ}{2} (J - 1 - Ω C_{1})^{2}

(A6)

where κ is the bulk modulus of the gel and is much larger than the shear modulus NkT. The constitutive relations are

σ = \frac{N k T}{J} ({F F}^{T} - I) - k T (\sum_{m = 2}^{m_{d}} c_{m}) I + κ (J - 1 - Ω C_{1}) I (A 7 a)

(A7a)

μ_{1} = μ_{1}^{0} + k T (\ln (\frac{Ω C_{1}}{1 + Ω C_{1}}) + \frac{1 + Ω C_{1} + χ_{h}}{(1 + Ω C_{1})^{2}}) - Ω κ (J - 1 - Ω C_{1}) (A 7 b)

(A7b)

μ_{m} = μ_{m}^{0} + k T \ln \frac{c_{m}}{c_{m}^{ref}} + e z_{m} ϕ, m = 2, 3, \dots, m_{d} (A 7 c)

(A7c)

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

Appendix B: Dimensionless Equations

We normalize all stresses by kT/Ω, all lengths by the characteristic length of the system L₀, time by

L_{0}^{2} / D_{1}

⁠, all concentrations by 1/Ω, electric potential by kT/e, all diffusivities by D₁, and all fluxes by D₁/ΩL₀ and define dimensionless chemical potential as

{\hat{μ}}_{1} = \frac{μ_{1}}{k T}

(B1)

{\hat{μ}}_{m} = \exp (\frac{μ_{m}}{k T}), m = 2, 3, \dots, m_{d}

(B2)

All equations in the last section are listed in the dimensionless form. To make the expression consistent with the last section, all quantities in the following text are dimensionless and ^ is ignored.

{\dot{c}}_{m} + c_{m} \frac{\dot{J}}{J} + div j_{m} = 0, m = 1, 2, \dots, m_{d} (B 3 a)

(B3a)

div σ = 0 (B 3 b)

(B3b)

\sum_{m = 0}^{m_{d}} z_{m} c_{m} = 0 (B 3 c)

(B3c)

j_{1} = - (c_{1} grad μ_{1} + c_{2} \frac{grad μ_{2}}{μ_{2}} + \dots + c_{m_{d}} \frac{grad μ_{m_{d}}}{μ_{m_{d}}}) (B 3 d)

(B3d)

\begin{aligned} j_{m} = & - \frac{c_{m}}{c_{1}} (c_{1} grad μ_{1} + c_{2} \frac{grad μ_{2}}{μ_{2}} + \dots + c_{m_{d}} \frac{grad μ_{m_{d}}}{μ_{m_{d}}}) \\ - D_{m} c_{m} \frac{grad μ_{m}}{μ_{m}}, m = 2, 3, \dots, m_{d} \end{aligned} (B 3 e)

(B3e)

σ = \frac{N Ω}{J} ({F F}^{T} - I) - (\sum_{m = 2}^{m_{d}} c_{m}) I + κ (J - 1 - C_{1}) I (B 3 f)

(B3f)

μ_{1} = μ_{1}^{0} + (\ln (\frac{C_{1}}{1 + C_{1}}) + \frac{1 + C_{1} + χ_{h}}{(1 + C)^{2}}) - κ (J - 1 - C_{1}) (B 3 g)

(B3g)

μ_{m} = μ_{m}^{0} \frac{c_{m}}{c_{m}^{ref}} \exp (z_{m} ϕ), m = 2, 3, \dots, m_{d} (B 3 h)

(B3h)

Appendix C: Numerical Solution Procedure

The strong forms of the differential equations on the body

B

and its boundary

\partial B

are

\begin{aligned} {\begin{cases} {\dot{c}}_{m} + c_{m} \frac{\dot{J}}{J} + div j_{m} = 0 & in B \\ μ_{m} = {\bar{μ}}_{m} & on S_{μ} \\ - j_{m} \cdot n = {\bar{j}}_{m} & on S_{j}, m = 1, 2, \dots, m_{d} \end{cases} \end{aligned}

(C1)

\begin{aligned} {\begin{cases} div σ = 0 & in B \\ u = \bar{u} & on S_{u} \\ σ n = \bar{t} & on S_{t} \end{cases} \end{aligned}

(C2)

where S_μ and S_j are the corresponding subsurfaces of the boundary where prescribed chemical potential

{\bar{μ}}_{m}

and prescribed flow

{\bar{j}}_{m}

are applied. S_u and S_t are the corresponding subsurfaces of the boundary where prescribed displacement

\bar{u}

and prescribed traction force

\bar{t}

are applied.

Two test fields w and s are applied to get the weak forms

\int_{B} {\dot{c}}_{m} s + c_{m} \frac{\dot{J}}{J} s - \frac{\partial s}{\partial x_{i}} j_{m i} d v - \int_{S_{j}} s {\bar{j}}_{m} d s = 0, m = 1, 2, \dots, m_{d}

(C3)

\int_{B} σ_{i j} \frac{\partial w_{i}}{\partial x_{j}} d v - \int_{S_{t}} w_{i} {\bar{t}}_{i} d s = 0

(C4)

Applying backward Euler scheme to Eq. to get

\int_{B} \frac{c_{m}^{t + Δ t} - c_{m}^{t}}{Δ t} s + \frac{c_{m}^{t + Δ t}}{J^{t + Δ t}} \frac{J^{t + Δ t} - J^{t}}{Δ t} s - \frac{\partial s}{\partial x_{i}} j_{m i}^{t + Δ t} d v - \int_{S_{j}} s {\bar{j}}_{m}^{t + Δ t} d s = 0

(C5)

Next, the body

B

is discretized into finite elements,

B = ⋃ B^{e}

⁠, and the displacement and chemical potentials are interpolated by

u_{i} = \sum u_{i}^{A} N^{A}, μ_{m} = \sum μ_{m}^{D} M^{D}

(C6)

where the indexes A = 1, 2, … and D = 1, 2, … denote the nodes of the element.

u_{i}^{A}

and

μ_{m}^{D}

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.

The standard Galerkin formulation is applied to interpolate w and s.

w_{i} = \sum w_{i}^{A} N^{A}, s = \sum s^{D} M^{D}

(C7)

Inserting Eqs. and into weak forms Eqs. and to obtain element-level residuals for displacement and chemical potentials.

(R_{u_{i}})^{A} = - \int_{B^{e}} σ_{i j} \frac{\partial N^{A}}{\partial x_{j}} d v + \int_{S_{t}} N^{A} {\bar{t}}_{i} d s (C 8 a)

(C8a)

\begin{aligned} (R_{μ_{1}})^{D} & = \int_{B^{e}} 2 c_{1} M^{D} - c_{1}^{t} M^{D} - \frac{c_{1}}{J} J^{t} M^{D} - \frac{\partial M^{D}}{\partial x_{i}} j_{1 i} Δ t d v \\ - \int_{S_{j}} M^{D} {\bar{j}}_{1} Δ t d s \end{aligned} (C 8 b)

(C8b)

\begin{aligned} (R_{μ_{m}})^{D} & = \int_{B^{e}} 2 c_{m} M^{D} - c_{m}^{t} M^{D} - \frac{c_{m}}{J} J^{t} M^{D} - \frac{\partial M^{D}}{\partial x_{i}} j_{m i} Δ t d v \\ - \int_{S_{j}} M^{D} {\bar{j}}_{m} Δ t d s, m = 2, \dots, m_{d} \end{aligned} (C 8 c)

(C8c)

where all values without superscript t denotes the values at time t + Δt.

In order to solve the problem using Newton's method, the following nine tangents are required in the iterative procedure,

\begin{aligned} {(K_{u_{i} u_{s}})}^{A B} = - \frac{\partial R_{u_{i}}^{A}}{\partial u_{s}^{B}}, {(K_{u_{i} μ_{1}})}^{A E} = - \frac{\partial R_{u_{i}}^{A}}{\partial μ_{1}^{E}}, {(K_{u_{i} μ_{m}})}^{A E} = - \frac{\partial R_{u_{i}}^{A}}{\partial μ_{m}^{E}} \\ {(K_{μ_{1} u_{s}})}^{D B} = - \frac{\partial R_{μ_{1}}^{D}}{\partial u_{s}^{B}}, {(K_{μ_{1} μ_{1}})}^{D E} = - \frac{\partial R_{μ_{1}}^{D}}{\partial {μ_{1}}^{E}}, {(K_{μ_{1} μ_{m}})}^{D E} = - \frac{\partial R_{μ_{1}}^{D}}{\partial {μ_{m}}^{E}} \\ {(K_{μ_{m} u_{s}})}^{D B} = - \frac{\partial R_{μ_{m}}^{D}}{\partial u_{s}^{B}}, {(K_{μ_{m} μ_{1}})}^{D E} = - \frac{\partial R_{μ_{m}}^{D}}{\partial {μ_{1}}^{E}}, {(K_{μ_{m} μ_{n}})}^{D E} = - \frac{\partial R_{μ_{m}}^{D}}{\partial {μ_{n}}^{E}}, \\ m, n = 2, \dots, m_{d} \end{aligned}

(C9)

Then,

K_{u_{i} u_{s}}^{A B} = \int_{B^{e}} \frac{\partial N^{A}}{\partial x_{j}} J^{- 1} F_{j m} F_{l n} \frac{\partial P_{i m}}{\partial F_{s n}} \frac{\partial N^{B}}{\partial x_{l}} d v

(C10)

K_{u_{i} μ_{1}}^{A E} = \int_{B^{e}} \frac{\partial N^{A}}{\partial x_{j}} \frac{\partial σ_{i j}}{\partial μ_{1}} M^{E} d v = - \int_{B^{e}} \frac{\partial N^{A}}{\partial x_{i}} κ \frac{\partial C_{1}}{\partial μ_{1}} M^{E} d v

(C11)

\begin{aligned} K_{u_{i} μ_{m}}^{A E} & = \int_{B^{e}} \frac{\partial N^{A}}{\partial x_{j}} \frac{\partial σ_{i j}}{\partial μ_{m}} M^{E} d v = - \int_{B^{e}} \frac{\partial N^{A}}{\partial x_{i}} (\frac{\partial c_{2}}{\partial μ_{m}} + \dots + \frac{\partial c_{m_{d}}}{\partial μ_{m}}) M^{E} d v, \\ m = 2, \dots, m_{d} \end{aligned}

(C12)

\begin{aligned} K_{μ_{1} u_{s}}^{D B} & = - \int_{B^{e}} M^{D} \frac{\partial N^{B}}{\partial x_{s}} ((2 J - J^{t}) \frac{\partial c_{1}}{\partial J} + \frac{c_{1}}{J} J^{t}) d v - \int_{B^{e}} M^{D} \frac{\partial N^{B}}{\partial x_{s}} (2 c_{1} - c_{1}^{t} - \frac{c_{1}}{J} J^{t}) d v \\ + \int_{B^{e}} (\frac{\partial M^{D}}{\partial x_{s}} \frac{\partial N^{B}}{\partial x_{i}} - \frac{\partial M^{D}}{\partial x_{i}} \frac{\partial N^{B}}{\partial x_{s}}) (c_{1} \frac{\partial μ_{1}}{\partial x_{i}} + \frac{c_{2}}{μ_{2}} \frac{\partial μ_{2}}{\partial x_{i}} + \dots + \frac{c_{m_{d}}}{μ_{m_{d}}} \frac{\partial μ_{m_{d}}}{\partial x_{i}}) Δ t d v \\ + \int_{B^{e}} \frac{\partial M^{D}}{\partial x_{i}} \frac{\partial N^{B}}{\partial x_{i}} (c_{1} \frac{\partial μ_{1}}{\partial x_{s}} + \frac{c_{2}}{μ_{2}} \frac{\partial μ_{2}}{\partial x_{s}} + \dots + \frac{c_{m_{d}}}{μ_{m_{d}}} \frac{\partial μ_{m_{d}}}{\partial x_{s}}) Δ t d v \\ - \int_{B^{e}} \frac{\partial M^{D}}{\partial x_{i}} \frac{\partial N^{B}}{\partial x_{s}} J (\frac{\partial c_{1}}{\partial J} \frac{\partial μ_{1}}{\partial x_{i}} + \frac{\partial c_{2}}{\partial J} \frac{1}{μ_{2}} \frac{\partial μ_{2}}{\partial x_{i}} \dots + \frac{\partial c_{m_{d}}}{\partial J} \frac{1}{μ_{m_{d}}} \frac{\partial μ_{m_{d}}}{\partial x_{i}}) Δ t d v \end{aligned}

(C13)

K_{μ_{1} μ_{1}}^{D E} = \int_{B^{e}} (\frac{J^{t}}{J} - 2) \frac{\partial c_{1}}{\partial μ_{1}} M^{D} M^{E} d v - \int_{B^{e}} c_{1} \frac{\partial M^{D}}{\partial x_{i}} \frac{\partial M^{E}}{\partial x_{i}} Δ t d v + \int_{B^{e}} \frac{\partial M^{D}}{\partial x_{i}} M^{E} \frac{\partial c_{1}}{\partial μ_{1}} \frac{\partial μ_{1}}{\partial x_{i}} Δ t d v

(C14)

\begin{aligned} K_{μ_{1} μ_{m}}^{A B} & = - \int_{B^{e}} (\frac{\partial c_{2}}{\partial μ_{m}} \frac{1}{μ_{2}} \frac{\partial μ_{2}}{\partial x_{i}} + \dots + \frac{\partial c_{m_{d}}}{\partial μ_{m}} \frac{1}{μ_{m_{d}}} \frac{\partial μ_{m_{d}}}{\partial x_{i}} - \frac{c_{m}}{μ_{m}^{2}} \frac{\partial μ_{m}}{\partial x_{i}}) \frac{\partial M^{D}}{\partial x_{i}} M^{E} Δ t d v \\ - \int_{B^{e}} \frac{c_{m}}{μ_{m}} \frac{\partial M^{D}}{\partial x_{i}} \frac{\partial M^{E}}{\partial x_{i}} Δ t d v, m = 2, \dots, m_{d} \end{aligned}

(C15)

\begin{aligned} K_{μ_{m} u_{s}}^{A B} & = - \int_{B^{e}} M^{D} \frac{\partial N^{B}}{\partial x_{s}} ((2 J - J^{t}) \frac{\partial c_{m}}{\partial J} + \frac{c_{m}}{J} J^{t}) d v - \int_{B^{e}} M^{D} \frac{\partial N^{B}}{\partial x_{s}} (2 c_{m} - c_{m}^{t} - \frac{c_{m}}{J} J^{t}) d v \\ + \int_{B^{e}} (\frac{\partial M^{D}}{\partial x_{s}} \frac{\partial M^{B}}{\partial x_{i}} - \frac{\partial M^{D}}{\partial x_{i}} \frac{\partial N^{B}}{\partial x_{s}}) [\frac{c_{m}}{c_{1}} (c_{1} \frac{\partial μ_{1}}{\partial x_{i}} + \frac{c_{2}}{μ_{2}} \frac{\partial μ_{2}}{\partial x_{i}} + \dots + \frac{c_{m_{d}}}{μ_{m_{d}}} \frac{\partial μ_{m_{d}}}{\partial x_{i}}) + \frac{D_{m} c_{m}}{μ_{m}} \frac{\partial μ_{m}}{\partial x_{i}}] Δ t d v \\ - \int_{B^{e}} \frac{\partial M^{D}}{\partial x_{i}} \frac{\partial N^{B}}{\partial x_{s}} J (\frac{\partial c_{m}}{\partial J} \frac{\partial μ_{1}}{\partial x_{i}} + \frac{\partial (c_{m} c_{2} / c_{1})}{\partial J} \frac{1}{μ_{2}} \frac{\partial μ_{2}}{\partial x_{i}} + \dots + \frac{\partial (c_{m} c_{m_{d}} / c_{1})}{\partial J} \frac{1}{μ_{m_{d}}} \frac{\partial μ_{m_{d}}}{\partial x_{i}} + \frac{D_{m}}{μ_{m}} \frac{\partial c_{m}}{\partial J} \frac{\partial μ_{m}}{\partial x_{i}}) Δ t d v \\ + \int_{B^{e}} \frac{\partial M^{D}}{\partial x_{i}} \frac{\partial M^{B}}{\partial x_{i}} (c_{m} \frac{\partial μ_{1}}{\partial x_{s}} + \frac{c_{m} c_{2}}{c_{1} μ_{2}} \frac{\partial μ_{2}}{\partial x_{s}} + \dots + \frac{c_{m} c_{m_{d}}}{c_{1} μ_{m_{d}}} \frac{\partial μ_{m_{d}}}{\partial x_{s}} + \frac{D_{m} c_{m}}{μ_{m}} \frac{\partial μ_{m}}{\partial x_{s}}) Δ t d v, m = 2, \dots, m_{d} \end{aligned}

(C16)

\begin{aligned} K_{μ_{m} μ_{1}}^{D E} & = - \int_{B^{e}} \frac{\partial M^{D}}{\partial x_{i}} c_{m} \frac{\partial M^{E}}{\partial x_{i}} Δ t d v + \int_{B^{e}} \frac{\partial M^{D}}{\partial x_{i}} M^{E} \frac{c_{m}}{c_{1}^{2}} \frac{\partial c_{1}}{\partial μ_{1}} (\frac{c_{2}}{μ_{2}} \frac{\partial μ_{2}}{\partial x_{i}} + \dots + \frac{c_{m_{d}}}{μ_{m_{d}}} \frac{\partial μ_{m_{d}}}{\partial x_{i}}) Δ t d v \\ m = 2, \dots, m_{d} \end{aligned}

(C17)

\begin{aligned} K_{μ_{m} μ_{n}}^{D E} & = - \int_{B^{e}} M^{D} M^{E} (2 - \frac{J^{t}}{J}) \frac{\partial c_{m}}{\partial μ_{n}} d v - \int_{B^{e}} \frac{\partial M^{D}}{\partial x_{i}} M^{E} \frac{c_{m}}{c_{1}} (\frac{\partial c_{2}}{\partial μ_{n}} \frac{1}{μ_{2}} \frac{\partial μ_{2}}{\partial x_{i}} + \dots + \frac{\partial c_{m_{d}}}{\partial μ_{n}} \frac{1}{μ_{m_{d}}} \frac{\partial μ_{m_{d}}}{\partial x_{i}}) Δ t d v \\ - \int_{B^{e}} \frac{\partial M^{D}}{\partial x_{i}} M^{E} \frac{\partial c_{m}}{\partial μ_{n}} (\frac{1}{c_{1}} (c_{1} \frac{\partial μ_{1}}{\partial x_{i}} + \frac{c_{2}}{μ_{2}} \frac{\partial μ_{2}}{\partial x_{i}} + \dots + \frac{c_{m_{d}}}{μ_{m_{d}}} \frac{\partial μ_{m_{d}}}{\partial x_{i}}) + \frac{D_{m}}{μ_{m}} \frac{\partial μ_{m}}{\partial x_{i}}) Δ t d v \\ + \int_{B^{e}} \frac{\partial M^{D}}{\partial x_{i}} M^{E} (\frac{c_{m}}{c_{1}} \frac{c_{n}}{μ_{n}^{2}} + \frac{c_{n} D_{n}}{μ_{n}^{2}} δ_{m n}) \frac{\partial μ_{n}}{\partial x_{i}} Δ t d v - \int_{B^{e}} \frac{\partial M^{D}}{\partial x_{i}} \frac{\partial M^{E}}{\partial x_{i}} (\frac{c_{m} c_{n}}{c_{1} μ_{n}} + \frac{D_{n} c_{n}}{μ_{n}} δ_{m n}) Δ t d v \\ m, n = 2, \dots, m_{d} . \end{aligned}

(C18)

Then all residuals and stiffness matrices are assembled into one residual vector and one stiffness matrix.

R = [\begin{matrix} R_{u_{1}}^{1} & R_{u_{2}}^{1} & R_{u_{3}}^{1} & R_{μ_{1}}^{1} & R_{μ_{2}}^{1} & \dots & R_{μ_{m_{d}}}^{1} & \dots & R_{u_{1}}^{M} & R_{u_{2}}^{M} & R_{u_{3}}^{M} & R_{μ_{1}}^{M} & R_{μ_{2}}^{M} & \dots & R_{μ_{m_{d}}}^{M} \end{matrix}]

(C19)

K = [\begin{matrix} K_{u_{1} u_{1}}^{11} & K_{u_{1} u_{2}}^{11} & K_{u_{1} u_{3}}^{11} & K_{u_{1} μ_{1}}^{11} & K_{u_{1} μ_{2}}^{11} & \dots & K_{u_{1} μ_{m_{d}}}^{11} & \dots & K_{u_{1} u_{1}}^{1 M} & K_{u_{1} u_{2}}^{1 M} & K_{u_{1} u_{3}}^{1 M} & K_{u_{1} μ_{1}}^{1 M} & K_{u_{1} μ_{2}}^{1 M} & \dots & K_{u_{1} μ_{m_{d}}}^{1 M} \\ K_{u_{2} u_{1}}^{11} & K_{u_{2} u_{2}}^{11} & K_{u_{2} u_{3}}^{11} & K_{u_{2} μ_{1}}^{11} & K_{u_{2} μ_{2}}^{11} & \dots & K_{u_{2} μ_{m_{d}}}^{11} & \dots & K_{u_{2} u_{1}}^{1 M} & K_{u_{2} u_{2}}^{1 M} & K_{u_{2} u_{3}}^{1 M} & K_{u_{2} μ_{1}}^{1 M} & K_{u_{2} μ_{2}}^{1 M} & \dots & K_{u_{2} μ_{m_{d}}}^{1 M} \\ K_{u_{3} u_{1}}^{11} & K_{u_{3} u_{2}}^{11} & K_{u_{3} u_{3}}^{11} & K_{u_{3} μ_{1}}^{11} & K_{u_{3} μ_{2}}^{11} & \dots & K_{u_{3} μ_{m_{d}}}^{11} & \dots & K_{u_{3} u_{1}}^{1 M} & K_{u_{3} u_{2}}^{1 M} & K_{u_{3} u_{3}}^{1 M} & K_{u_{3} μ_{1}}^{1 M} & K_{u_{3} μ_{2}}^{1 M} & \dots & K_{u_{3} μ_{m_{d}}}^{1 M} \\ K_{μ_{1} u_{1}}^{11} & K_{μ_{1} u_{2}}^{11} & K_{μ_{1} u_{3}}^{11} & K_{μ_{1} μ_{1}}^{11} & K_{μ_{1} μ_{2}}^{11} & \dots & K_{μ_{1} μ_{m_{d}}}^{11} & \dots & K_{μ_{1} u_{1}}^{1 M} & K_{μ_{1} u_{2}}^{1 M} & K_{μ_{1} u_{3}}^{1 M} & K_{μ_{1} μ_{1}}^{1 M} & K_{μ_{1} μ_{2}}^{1 M} & \dots & K_{μ_{1} μ_{m_{d}}}^{1 M} \\ K_{μ_{2} u_{1}}^{11} & K_{μ_{2} u_{2}}^{11} & K_{μ_{2} u_{3}}^{11} & K_{μ_{2} μ_{1}}^{11} & K_{μ_{2} μ_{2}}^{11} & \dots & K_{μ_{2} μ_{m_{d}}}^{11} & \dots & K_{μ_{2} u_{1}}^{1 M} & K_{μ_{2} u_{2}}^{1 M} & K_{μ_{2} u_{3}}^{1 M} & K_{μ_{2} μ_{1}}^{1 M} & K_{μ_{2} μ_{2}}^{1 M} & \dots & K_{μ_{2} μ_{m_{d}}}^{1 M} \\ \dots & ⋱ & ⋮ & ⋱ & \dots & ⋱ & ⋮ \\ K_{μ_{m_{d}} u_{1}}^{11} & K_{μ_{m_{d}} u_{2}}^{11} & K_{μ_{m_{d}} u_{3}}^{11} & K_{μ_{m_{d}} μ_{1}}^{11} & K_{μ_{m_{d}} μ_{2}}^{11} & \dots & K_{μ_{m_{d}} μ_{m_{d}}}^{11} & \dots & K_{μ_{m_{d}} u_{1}}^{1 M} & K_{μ_{m_{d}} u_{2}}^{1 M} & K_{μ_{m_{d}} u_{3}}^{1 M} & K_{μ_{m_{d}} μ_{1}}^{1 M} & K_{μ_{m_{d}} μ_{2}}^{1 M} & \dots & K_{μ_{m_{d}} μ_{m_{d}}}^{1 M} \\ \dots & ⋱ & \dots \\ K_{u_{1} u_{1}}^{M 1} & K_{u_{1} u_{2}}^{M 1} & K_{u_{1} u_{3}}^{M 1} & K_{u_{1} μ_{1}}^{M 1} & K_{u_{1} μ_{2}}^{M 1} & \dots & K_{u_{1} μ_{m_{d}}}^{M 1} & \dots & K_{u_{1} u_{1}}^{M M} & K_{u_{1} u_{2}}^{M M} & K_{u_{1} u_{3}}^{M M} & K_{u_{1} μ_{1}}^{M M} & K_{u_{1} μ_{2}}^{M M} & \dots & K_{u_{1} μ_{m_{d}}}^{M M} \\ K_{u_{2} u_{1}}^{M 1} & K_{u_{2} u_{2}}^{M 1} & K_{u_{2} u_{3}}^{M 1} & K_{u_{2} μ_{1}}^{M 1} & K_{u_{2} μ_{2}}^{M 1} & \dots & K_{u_{2} μ_{m_{d}}}^{M 1} & \dots & K_{u_{2} u_{1}}^{M M} & K_{u_{2} u_{2}}^{M M} & K_{u_{2} u_{3}}^{M M} & K_{u_{2} μ_{1}}^{M M} & K_{u_{2} μ_{2}}^{M M} & \dots & K_{u_{2} μ_{m_{d}}}^{M M} \\ K_{u_{3} u_{1}}^{M 1} & K_{u_{3} u_{2}}^{M 1} & K_{u_{3} u_{3}}^{M 1} & K_{u_{3} μ_{1}}^{M 1} & K_{u_{3} μ_{2}}^{M 1} & \dots & K_{u_{3} μ_{m_{d}}}^{M 1} & \dots & K_{u_{3} u_{1}}^{M M} & K_{u_{3} u_{2}}^{M M} & K_{u_{3} u_{3}}^{M M} & K_{u_{3} μ_{1}}^{M M} & K_{u_{3} μ_{2}}^{M M} & \dots & K_{u_{3} μ_{m_{d}}}^{M M} \\ K_{μ_{1} u_{1}}^{M 1} & K_{μ_{1} u_{2}}^{M 1} & K_{μ_{1} u_{3}}^{M 1} & K_{μ_{1} μ_{1}}^{M 1} & K_{μ_{1} μ_{2}}^{M 1} & \dots & K_{μ_{1} μ_{m_{d}}}^{M 1} & \dots & K_{μ_{1} u_{1}}^{M M} & K_{μ_{1} u_{2}}^{M M} & K_{μ_{1} u_{3}}^{M M} & K_{μ_{1} μ_{1}}^{M M} & K_{μ_{1} μ_{2}}^{M M} & \dots & K_{μ_{1} μ_{m_{d}}}^{M M} \\ K_{μ_{2} u_{1}}^{M 1} & K_{μ_{2} u_{2}}^{M 1} & K_{μ_{2} u_{3}}^{M 1} & K_{μ_{2} μ_{1}}^{M 1} & K_{μ_{2} μ_{2}}^{M 1} & \dots & K_{μ_{2} μ_{m_{d}}}^{M 1} & \dots & K_{μ_{2} u_{1}}^{M M} & K_{μ_{2} u_{2}}^{M M} & K_{μ_{2} u_{3}}^{M M} & K_{μ_{2} μ_{1}}^{M M} & K_{μ_{2} μ_{2}}^{M M} & \dots & K_{μ_{2} μ_{m_{d}}}^{M M} \\ \dots & ⋱ & ⋮ & ⋱ & \dots & ⋱ & ⋮ \\ K_{μ_{m_{d}} u_{1}}^{M 1} & K_{μ_{m_{d}} u_{2}}^{M 1} & K_{μ_{m_{d}} u_{3}}^{M 1} & K_{μ_{m_{d}} μ_{1}}^{M 1} & K_{μ_{m_{d}} μ_{2}}^{M 1} & \dots & K_{μ_{m_{d}} μ_{m_{d}}}^{M 1} & \dots & K_{μ_{m_{d}} u_{1}}^{M M} & K_{μ_{m_{d}} u_{2}}^{M M} & K_{μ_{m_{d}} u_{3}}^{M M} & K_{μ_{m_{d}} μ_{1}}^{M M} & K_{μ_{m_{d}} μ_{2}}^{M M} & \dots & K_{μ_{m_{d}} μ_{m_{d}}}^{M M} \end{matrix}]

(C20)

References

1.

Brannon-Peppas

,

L.

, and

Peppas

,

N. A.

,

1991

, “

Equilibrium Swelling Behavior of Ph-Sensitive Hydrogels

,”

Chem. Eng. Sci.

,

46

(

3

), pp.

715

–

722

. 10.1016/0009-2509(91)80177-Z

Google Scholar

Crossref

2.

Horta

,

A.

,

Molina

,

M. J.

,

Gómez-Antón

,

M. R.

, and

Piérola

,

I. F.

,

2009

, “

The Ph Inside a Ph-Sensitive Gel Swollen in Aqueous Salt Solutions: Poly (n-Vinylimidazole)

,”

Macromolecules

,

42

(

4

), pp.

1285

–

1292

. 10.1021/ma802204b

Google Scholar

Crossref

3.

Marek

,

S. R.

,

Conn

,

C. A.

, and

Peppas

,

N. A.

,

2010

, “

Cationic Nanogels Based on Diethylaminoethyl Methacrylate

,”

Polymer

,

51

(

6

), pp.

1237

–

1243

. 10.1016/j.polymer.2010.01.060

Google Scholar

Crossref

PubMed

4.

Jabbari

,

E.

,

Tavakoli

,

J.

, and

Sarvestani

,

A. S.

,

2007

, “

Swelling Characteristics of Acrylic Acid Polyelectrolyte Hydrogel in a DC Electric Field

,”

Smart Mater. Struct.

,

16

(

5

), p.

1614

. 10.1088/0964-1726/16/5/015

Google Scholar

Crossref

5.

Lin

,

J.

,

Tang

,

Q.

,

Hu

,

D.

,

Sun

,

X.

,

Li

,

Q.

, and

Wu

,

J.

,

2009

, “

Electric Field Sensitivity of Conducting Hydrogels With Interpenetrating Polymer Network Structure

,”

Colloids. Surf., A

,

346

(

1–3

), pp.

177

–

183

. 10.1016/j.colsurfa.2009.06.011

Google Scholar

Crossref

6.

Sumaru

,

K.

,

Takagi

,

T.

,

Sugiura

,

S.

, and

Kanamori

,

T.

,

2014

,

Soft Actuators

,

K.

Asaka

and

H.

Okuzaki

, eds.,

Springer

,

New York

, pp.

219

–

229

.

7.

Zhang

,

Q. M.

,

Xu

,

W.

, and

Serpe

,

M. J.

,

2014

, “

Optical Devices Constructed From Multiresponsive Microgels

,”

Angew. Chem., Int. Ed.

,

53

(

19

), pp.

4827

–

4831

. 10.1002/anie.201402641

Google Scholar

Crossref

8.

Suzuki

,

A.

, and

Tanaka

,

T.

,

1990

, “

Phase Transition in Polymer Gels Induced by Visible Light

,”

Nature

,

346

(

6282

), p.

345

. 10.1038/346345a0

Google Scholar

Crossref

9.

Drury

,

J. L.

, and

Mooney

,

D. J.

,

2003

, “

Hydrogels for Tissue Engineering: Scaffold Design Variables and Applications

,”

Biomaterials

,

24

(

24

), pp.

4337

–

4351

. 10.1016/S0142-9612(03)00340-5

Google Scholar

Crossref

PubMed

10.

Hu

,

Y.

,

You

,

J.-O.

, and

Aizenberg

,

J.

,

2016

, “

Micropatterned Hydrogel Surface With High-Aspect-Ratio Features for Cell Guidance and Tissue Growth

,”

ACS Appl. Mater. Interfaces

,

8

(

34

), pp.

21939

–

21945

. 10.1021/acsami.5b12268

Google Scholar

Crossref

PubMed

11.

Jen

,

A. C.

,

Wake

,

M. C.

, and

Mikos

,

A. G.

,

1996

, “

Hydrogels for Cell Immobilization

,”

Biotechnol. Bioeng.

,

50

(

4

), pp.

357

–

364

. 10.1002/(SICI)1097-0290(19960520)50:4<357::AID-BIT2>3.0.CO;2-K

Google Scholar

Crossref

PubMed

12.

Hoare

,

T. R.

, and

Kohane

,

D. S.

,

2008

, “

Hydrogels in Drug Delivery: Progress and Challenges

,”

Polymer

,

49

(

8

), pp.

1993

–

2007

. 10.1016/j.polymer.2008.01.027

Google Scholar

Crossref

13.

Li

,

J.

, and

Mooney

,

D. J.

,

2016

, “

Designing Hydrogels for Controlled Drug Delivery

,”

Nat. Rev. Mater.

,

1

(

12

), p.

16071

. 10.1038/natrevmats.2016.71

Google Scholar

Crossref

PubMed

14.

Qiu

,

Y.

, and

Park

,

K.

,

2001

, “

Environment-Sensitive Hydrogels for Drug Delivery

,”

Adv. Drug Delivery Rev.

,

53

(

3

), pp.

321

–

339

. 10.1016/S0169-409X(01)00203-4

Google Scholar

Crossref

15.

Beebe

,

D. J.

,

Moore

,

J. S.

,

Bauer

,

J. M.

,

Yu

,

Q.

,

Liu

,

R. H.

,

Devadoss

,

C.

, and

Jo

,

B.-H.

,

2000

, “

Functional Hydrogel Structures for Autonomous Flow Control Inside Microfluidic Channels

,”

Nature

,

404

(

6778

), p.

588

16.

Dong

,

L.

,

Agarwal

,

A. K.

,

Beebe

,

D. J.

, and

Jiang

,

H.

,

2006

, “

Adaptive Liquid Microlenses Activated by Stimuli-Responsive Hydrogels

,”

Nature

,

442

(

7102

), p.

551

. 10.1038/nature05024

Google Scholar

Crossref

PubMed

17.

Gerlach

,

G.

,

Guenther

,

M.

,

Sorber

,

J.

,

Suchaneck

,

G.

,

Arndt

,

K.-F.

, and

Richter

,

A.

,

2005

, “

Chemical and pH Sensors Based on the Swelling Behavior of Hydrogels

,”

Sens. Actuators, B

,

111

, pp.

555

–

561

. 10.1016/j.snb.2005.03.040

Google Scholar

Crossref

18.

Richter

,

A.

,

Paschew

,

G.

,

Klatt

,

S.

,

Lienig

,

J.

,

Arndt

,

K.-F.

, and

Adler

,

H.-J.

,

2008

, “

Review on Hydrogel-Based Ph Sensors and Microsensors

,”

Sensors

,

8

(

1

), pp.

561

–

581

19.

Maeda

,

S.

,

Hara

,

Y.

,

Sakai

,

T.

,

Yoshida

,

R.

, and

Hashimoto

,

S.

,

2007

, “

Self-Walking Gel

,”

Adv. Mater.

,

19

(

21

), pp.

3480

–

3484

. 10.1002/adma.200700625

Google Scholar

Crossref

20.

Kreuer

,

K.

,

2001

, “

On the Development of Proton Conducting Polymer Membranes for Hydrogen and Methanol Fuel Cells

,”

J. Membr. Sci.

,

185

(

1

), pp.

29

–

39

. 10.1016/S0376-7388(00)00632-3

Google Scholar

Crossref

21.

Keplinger

,

C.

,

Sun

,

J.-Y.

,

Foo

,

C. C.

,

Rothemund

,

P.

,

Whitesides

,

G. M.

, and

Suo

,

Z.

,

2013

, “

Stretchable, Transparent, Ionic Conductors

,”

Science

,

341

(

6149

), pp.

984

–

987

. 10.1126/science.1240228

Google Scholar

Crossref

PubMed

22.

Yang

,

C. H.

,

Chen

,

B.

,

Lu

,

J. J.

,

Yang

,

J. H.

,

Zhou

,

J.

,

Chen

,

Y. M.

, and

Suo

,

Z.

,

2015

, “

Ionic Cable

,”

Extreme Mech. Lett.

,

3

, pp.

59

–

65

. 10.1016/j.eml.2015.03.001

Google Scholar

Crossref

23.

Kim

,

C.-C.

,

Lee

,

H.-H.

,

Oh

,

K. H.

, and

Sun

,

J.-Y.

,

2016

, “

Highly Stretchable, Transparent Ionic Touch Panel

,”

Science

,

353

(

6300

), pp.

682

–

687

. 10.1126/science.aaf8810

Google Scholar

Crossref

PubMed

24.

Xue

,

S.-L.

,

Li

,

B.

,

Feng

,

X.-Q.

, and

Gao

,

H.

,

2017

, “

A Non-Equilibrium Thermodynamic Model for Tumor Extracellular Matrix With Enzymatic Degradation

,”

J. Mech. Phys. Solids

,

104

, pp.

32

–

56

. 10.1016/j.jmps.2017.04.002

Google Scholar

Crossref

25.

Swartz

,

M. A.

, and

Fleury

,

M. E.

,

2007

, “

Interstitial Flow and Its Effects in Soft Tissues

,”

Annu. Rev. Biomed. Eng.

,

9

, pp.

229

–

256

. 10.1146/annurev.bioeng.9.060906.151850

Google Scholar

Crossref

PubMed

26.

Levick

,

J.

,

1987

, “

Flow Through Interstitium and Other Fibrous Matrices

,”

Q. J. Exp. Physiol.

,

72

(

4

), pp.

409

–

437

. 10.1113/expphysiol.1987.sp003085

Google Scholar

Crossref

PubMed

27.

Frank

,

E. H.

, and

Grodzinsky

,

A. J.

,

1987

, “

Cartilage Electromechanics—I. Electrokinetic Transduction and the Effects of Electrolyte Ph and Ionic Strength

,”

J. Biomech.

,

20

(

6

), pp.

615

–

627

. 10.1016/0021-9290(87)90282-X

Google Scholar

Crossref

PubMed

28.

Fraldi

,

M.

, and

Carotenuto

,

A. R.

,

2018

, “

Cells Competition in Tumor Growth Poroelasticity

,”

J. Mech. Phys. Solids

,

112

, pp.

345

–

367

. 10.1016/j.jmps.2017.12.015

Google Scholar

Crossref

29.

Krishna

,

R.

, and

Wesselingh

,

J.

,

1997

, “

The Maxwell-Stefan Approach to Mass Transfer

,”

Chem. Eng. Sci.

,

52

(

6

), pp.

861

–

911

. 10.1016/S0009-2509(96)00458-7

Google Scholar

Crossref

30.

Lai

,

W. M.

,

Hou

,

J.

, and

Mow

,

V. C.

,

1991

, “

A Triphasic Theory for the Swelling and Deformation Behaviors of Articular Cartilage

,”

ASME J. Biomech. Eng.

,

113

(

3

), pp.

245

–

258

. 10.1115/1.2894880

Google Scholar

Crossref

31.

Maxwell

,

J. C.

,

1867

, “

IV. on the Dynamical Theory of Gases

,”

Philos. Trans. R. Soc. London

,

157

, pp.

49

–

88

. 10.1098/rstl.1867.0004

32.

Stefan

,

J.

,

1871

, “

Über Das Gleichgewicht Und Die Bewegung Insbesondere Die Diffusion Von Gasgemengen

,”

Sitzungsber. Akad. Wiss. Wien

,

63

(

2

), pp.

63

–

124

. Cited By 167.

33.

Vanag

,

V. K.

, and

Epstein

,

I. R.

,

2009

, “

Cross-Diffusion and Pattern Formation in Reaction–Diffusion Systems

,”

Phys. Chem. Chem. Phys.

,

11

(

6

), pp.

897

–

912

34.

Grim

,

E.

, and

Sollner

,

K.

,

1957

, “

The Contributions of Normal and Anomalous Osmosis to the Osmotic Effects Arising Across Charged Membranes With Solutions of Electrolytes

,”

J. Gen. Physiol.

,

40

(

6

), pp.

887

–

899

. 10.1085/jgp.40.6.887

Google Scholar

Crossref

PubMed

35.

Toyoshima

,

Y.

,

Kobatake

,

Y.

, and

Fujita

,

H.

,

1967

, “

Studies of Membrane Phenomena. Part 5.—Bulk Flow Through Membrane

,”

Trans. Faraday Soc.

,

63

, pp.

2828

–

2838

. 10.1039/TF9676302828

Google Scholar

Crossref

36.

Sasidhar

,

V.

, and

Ruckenstein

,

E.

,

1982

, “

Anomalous Effects During Electrolyte Osmosis Across Charged Porous Membranes

,”

J. Colloid Interface Sci.

,

85

(

2

), pp.

332

–

362

. 10.1016/0021-9797(82)90003-0

Google Scholar

Crossref

37.

Gu

,

W.

,

Lai

,

W.

, and

Mow

,

V.

,

1997

, “

A Triphasic Analysis of Negative Osmotic Flows Through Charged Hydrated Soft Tissues

,”

J. Biomech.

,

30

(

1

), pp.

71

–

78

. 10.1016/S0021-9290(96)00099-1

Google Scholar

Crossref

PubMed

38.

Gu

,

W.

,

Lai

,

W.

, and

Mow

,

V.

,

1998

, “

A Mixture Theory for Charged-Hydrated Soft Tissues Containing Multi-Electrolytes: Passive Transport and Swelling Behaviors

,”

ASME J. Biomech. Eng.

,

120

(

2

), pp.

169

–

180

. 10.1115/1.2798299

Google Scholar

Crossref

39.

Hong

,

W.

,

Zhao

,

X.

,

Zhou

,

J.

, and

Suo

,

Z.

,

2008

, “

A Theory of Coupled Diffusion and Large Deformation in Polymeric Gels

,”

J. Mech. Phys. Solids

,

56

(

5

), pp.

1779

–

1793

. 10.1016/j.jmps.2007.11.010

Google Scholar

Crossref

40.

Hong

,

W.

,

Zhao

,

X.

, and

Suo

,

Z.

,

2010

, “

Large Deformation and Electrochemistry of Polyelectrolyte Gels

,”

J. Mech. Phys. Solids

,

58

(

4

), pp.

558

–

577

. 10.1016/j.jmps.2010.01.005

Google Scholar

Crossref

41.

Beltran

,

S.

,

Baker

,

J. P.

,

Hooper

,

H. H.

,

Blanch

,

H. W.

, and

Prausnitz

,

J. M.

,

1991

, “

Swelling Equilibria for Weakly Ionizable, Temperature-Sensitive Hydrogels

,”

Macromolecules

,

24

(

2

), pp.

549

–

551

. 10.1021/ma00002a032

Google Scholar

Crossref

42.

Yadav

,

M.

, and

Rhee

,

K. Y.

,

2012

, “

Superabsorbent Nanocomposite (Alginate-g-Pamps/Mmt): Synthesis, Characterization and Swelling Behavior

,”

Carbohydr. Polym.

,

90

(

1

), pp.

165

–

173

. 10.1016/j.carbpol.2012.05.010

Google Scholar

Crossref

PubMed

43.

Tanaka

,

T.

,

Fillmore

,

D.

,

Sun

,

S.-T.

,

Nishio

,

I.

,

Swislow

,

G.

, and

Shah

,

A.

,

1980

, “

Phase Transitions in Ionic Gels

,”

Phys. Rev. Lett.

,

45

(

20

), p.

1636

. 10.1103/PhysRevLett.45.1636

Google Scholar

Crossref

44.

Tanaka

,

T.

, and

Fillmore

,

D. J.

,

1979

, “

Kinetics of Swelling of Gels

,”

J. Chem. Phys.

,

70

(

3

), pp.

1214

–

1218

. 10.1063/1.437602

Google Scholar

Crossref

45.

Dehghany

,

M.

,

Zhang

,

H.

,

Naghdabadi

,

R.

, and

Hu

,

Y.

,

2018

, “

A Thermodynamically-Consistent Large Deformation Theory Coupling Photochemical Reaction and Electrochemistry for Light-Responsive Gels

,”

J. Mech. Phys. Solids

,

116

, pp.

239

–

266

. 10.1016/j.jmps.2018.03.018

Google Scholar

Crossref

46.

Coleman

,

B. D.

, and

Noll

,

W.

,

1963

, “

The Thermodynamics of Elastic Materials With Heat Conduction and Viscosity

,”

Arch. Ration. Mech. Anal.

,

13

(

1

), pp.

167

–

178

. 10.1007/BF01262690

Google Scholar

Crossref

47.

De Groot

,

S. R.

, and

Mazur

,

P.

,

2013

,

Non-Equilibrium Thermodynamics

,

Courier Corporation

,

New York

.

48.

Onsager

,

L.

,

1931

, “

Reciprocal Relations in Irreversible Processes. II.

,”

Phys. Rev.

,

38

(

12

), p.

2265

. 10.1103/PhysRev.38.2265

Google Scholar

Crossref

49.

Flory

,

P. J.

,

1942

, “

Thermodynamics of High Polymer Solutions

,”

J. Chem. Phys.

,

10

(

1

), pp.

51

–

61

. 10.1063/1.1723621

Google Scholar

Crossref

50.

Flory

,

P. J.

, and

Rehner, Jr.

,

J.

,

1943

, “

Statistical Mechanics of Cross-Linked Polymer Networks I. Rubberlike Elasticity

,”

J. Chem. Phys.

,

11

(

11

), pp.

512

–

520

. 10.1063/1.1723791

Google Scholar

Crossref

51.

Flory

,

P.

,

1953

,

Principles of Polymer Chemistry

,

Cornell University Press

,

Ithaca, NY

.

52.

Huggins

,

M. L.

,

1941

, “

Solutions of Long Chain Compounds

,”

J. Chem. Phys.

,

9

(

5

), pp.

440

–

440

. 10.1063/1.1750930

Google Scholar

Crossref

53.

Debanne

,

D.

,

2004

, “

Information Processing in the Axon

,”

Nat. Rev. Neurosci.

,

5

(

4

), p.

304

54.

García-Grajales

,

J. A.

,

Jérusalem

,

A.

, and

Goriely

,

A.

,

2017

, “

Continuum Mechanical Modeling of Axonal Growth

,”

Comput. Methods Appl. Mech. Eng.

,

314

, pp.

147

–

163

. 10.1016/j.cma.2016.07.032

Google Scholar

Crossref

55.

Lodish

,

H.

,

Berk

,

A.

,

Kaiser

,

C. A.

,

Krieger

,

M.

,

Scott

,

M. P.

,

Bretscher

,

A.

,

Ploegh

,

H.

, and

Matsudaira

,

P.

,

2008

,

Molecular Cell Biology

,

Macmillan

,

New York

.

56.

Pinto

,

T. M.

,

Wedemann

,

R. S.

, and

Cortez

,

C. M.

,

2014

, “

Modeling the Electric Potential Across Neuronal Membranes: the Effect of Fixed Charges on Spinal Ganglion Neurons and Neuroblastoma Cells

,”

PLoS One

,

9

(

5

), p.

e96194

. 10.1371/journal.pone.0096194

Google Scholar

Crossref

PubMed

57.

Sun

,

Y.

,

Liu

,

S.

,

Du

,

G.

,

Gao

,

G.

, and

Fu

,

J.

,

2015

, “

Multi-Responsive and Tough Hydrogels Based on Triblock Copolymer Micelles As Multi-Functional Macro-Crosslinkers

,”

Chem. Commun.

,

51

(

40

), pp.

8512

–

8515

. 10.1039/C4CC10094H

Google Scholar

Crossref

58.

Bouklas

,

N.

,

Landis

,

C. M.

, and

Huang

,

R.

,

2015

, “

A Nonlinear, Transient Finite Element Method for Coupled Solvent Diffusion and Large Deformation of Hydrogels

,”

J. Mech. Phys. Solids

,

79

, pp.

21

–

43

. 10.1016/j.jmps.2015.03.004

Google Scholar

Crossref

Kinetics of Polyelectrolyte Gels

Abstract

1 Introduction

2 Thermodynamic Framework

2.1 Kinematics.

2.2 Mass Conservation.

2.3 Linear and Angular Momentum Balance.

2.4 Electrostatic Effects.

2.5 First Law of Thermodynamics.

2.6 Second Law of Thermodynamics.

2.7 Consequences of First and Second Laws.

3 Kinetic Law

4 A Specific Material Model

4.1 Explicit Form of Free Energy Function for Polyelectrolyte Gels.

4.2 Constitutive Relations.

4.3 Electrostatic Double Layer and Electroneutrality.

5 Results and Discussion

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

5.2 Three-Dimensional Swelling and Model Calibration.

5.3 One-Dimensional Swelling and Anomalous Displacement.

5.4 Bending With Abnormal Directions Due to Cross-Diffusion.

6 Conclusion

Acknowledgment

Appendix A: Formulas in Current Configuration

Appendix B: Dimensionless Equations

Appendix C: Numerical Solution Procedure

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

Kinetics of Polyelectrolyte Gels

Abstract

1 Introduction

2 Thermodynamic Framework

2.1 Kinematics.

2.2 Mass Conservation.

2.3 Linear and Angular Momentum Balance.

2.4 Electrostatic Effects.

2.5 First Law of Thermodynamics.

2.6 Second Law of Thermodynamics.

2.7 Consequences of First and Second Laws.

3 Kinetic Law

4 A Specific Material Model

4.1 Explicit Form of Free Energy Function for Polyelectrolyte Gels.

4.2 Constitutive Relations.

4.3 Electrostatic Double Layer and Electroneutrality.

5 Results and Discussion

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

5.2 Three-Dimensional Swelling and Model Calibration.

5.3 One-Dimensional Swelling and Anomalous Displacement.

5.4 Bending With Abnormal Directions Due to Cross-Diffusion.

6 Conclusion

Acknowledgment

Appendix A: Formulas in Current Configuration

Appendix B: Dimensionless Equations

Appendix C: Numerical Solution Procedure

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only