A mixed hybrid finite element framework for the simulation of swelling ionized hydrogels

Ionized hydrogels, as osmoelastic media, swell enormously (1000 times its original volume in unionized water) due to the osmotic pressure difference caused by the presence of the negatively charged ion groups attached to the solid matrix (polymer chains). The coupling between the extremely large deformations (induced by swelling) and fluid permeation is a field of application that regular poroelasticity formulations cannot handle. In this work, we present a mixed hybrid finite element (MHFE) computational framework featuring a three-field (deformation-chemical potential-flux) formulation. This formulation guarantees that mass conservation is preserved both locally and globally. The impact of such a property on the swelling simulations is demonstrated by four numerical examples in 2D. This paper focuses on the implementation aspects of the MHFE model and shows that it stays robust and accurate for a volume increase of more than 3000%.


Introduction
Hydrogels are water-swollen and cross-linked polymeric networks, produced by the simple reaction of one or more monomers. Hydrogels exhibit ability to swell and retain a significant fraction of water within its structure without dissolving. Various environmental stimuli can induce the swelling of a hydrogel. The common stimuli are pH, ion concentration in the solution, electric potential and temperature [1]. Hydrogels have received considerable attention in the past 50 years due to their outstanding application in wide range of fields including biomedicine [2].
Super Absorbent Polymers (SAPs), as a specific type of ionized hydrogels, can absorb and retain extraordinary large amounts of water. They can imbibe deionized water as high as 1000-100,000% of their own mass. In other words, 1 g SAPs can absorb 10-1000 g water. The absorption capacity of common hydrogels, however, is not more than 100% (1 g/g) [3]. SAPs due to its exceptional ability to absorb a large amount of water in a short time are widely used in hygiene products as well as for agricultural applications. Nowadays SAPs are made from partially neutralized, lightly cross-linked polyacrylic acid, which has been proven to give the best performance versus cost ratio [4]. A schematic illustration of the structure of a SAPs gel particle is given in Fig. 1.
There are several mechanisms leading to the swelling of SAPs. First of all, the backbone of the polymer chains (carboxylate acid group) is hydrophilic. As a result, when the polymer is in contact with water, water molecules are attracted to the polymer. Hydration happens and hydrogen bonds are formed as these activities increase the entropy of the system. This process is known as mixing. Secondly, besides mixing, for a partially neutralized gel the positive sodium ions are able to move relatively freely inside the gel as soon as the water molecules weaken the bonding force between them and negatively charged carboxylate groups. Consequently, the gel acts like a semi-permeable membrane and the difference in osmotic pressure arises in and outside of the gel. It has been proven that compared to the mixing  Fig. 1 Schematic structure of a SAPs gel particle part, the difference in osmotic pressure is the main drive contributing to the swelling of a SAPs gel [5].
To study the mechanical behavior of a hydrogel during swelling, a reliable swelling model is needed. Under the framework of mixture theory, Lanir [6] developed an osmoelastic model arguing that the ionic concentrations are always in equilibrium with the external solution. In essence, this is a biphasic (solid and fluid phase) swelling model. Later, Lai et al. [7] developed a triphasic (solid, fluid and ion phases) theory to model the swelling and deformation behavior in articular cartilage. Later, Huyghe and Janssen [8] proposed a finite deformation quadriphasic model where solid, fluid, cation and ion phases are all independently considered for incompressible porous media. Considering a SAPs gel particle facing a gush of urine which can be approximated by the physiological salt solution, i.e. a solution of Na+ and Cl− of concentration 0.154 mol/L, chemical potential responds much faster to the local ions concentration change than solvent permeating into the gel. This claim is justified by the following calculation. We have, on one hand, the diffusion coefficient d ion of ions in free water of the magnitude 10 −9 m 2 /s [9]. On the other hand, the hydraulic pressure diffusion coefficient d p is estimated by the multiplication of hydraulic permeability and Young's modulus. In a SAPs gel particle, its hydraulic permeability is of order 10 −3 mm 4 /(Ns) and Young's modulus is of the order 10 −2 N/mm 2 . Therefore, after some simple calculation, one finds that d ion /d p = 10 2 . For this reason, we assume ions respond infinitely fast to the presence of outer solution compared to the fluid. In other words, Lanir's osmoelastic model was found sufficient to characterize the swelling mechanism.
In order to simulate such a fluid permeation and solid deformation coupled problem, proper numerical schemes need to be implemented. Due to the essential focus on solid deformation, FEM (Finite Element Method) has been a natural choice for such swelling simulations. Lots of effort in numerical implementation are made by various groups over the years. Frijns et al. [10] implemented the quadriphasic model in one dimension to simulate the swelling and shrinking of biological tissues. The extension to the 3D implementation of quadriphasic model in finite deformation is done by Van Loon et al. [11]. Limiting to gel swelling simulations, Hong et al. [12] and Kuang and Huang [13] presented a finite element swelling model for inhomogeneous swelling at equilibrium state. A number of works ( [14][15][16]) take a different approach to model a swelling hydrogel. Namely the swelling front is treated as a sharp transition surface between the solution phase and the gel phase (characterized by different chemical potential). The literature listed above only dealt with spherical gel specimen. The generalization to an arbitrary gel shape is not a trivial job. The application of eXtended Finite Element method in such simulations have received wide attentions [17][18][19]. Extended Finite element is applied to simulate the crack propagation in hydrogel induced by swelling [20].
However, to achieve satisfactory simulation results, it is essential to deploy a discretization method that takes the physics of the problem into account. In swelling simulations, (local) solid deformation is directly related to (local) net in/out flow of the fluid. As a result, it is sensible to value the accurate calculation of the flux field in swelling simulations. In FEM, the flux field is calculated by numerical differentiation of the chemical potential field, which leads to a serious loss of precision [21]. On the other hand, MHFEM (Mixed Hybrid Finite Element Method), which approximates flux as an independent variable using Raviart-Thomas element and resolves the resulting indefinite coefficient matrix by means of hybridization procedure, possesses local mass conservation property and has proven to be effective in solving Darcy type equations [22,23]. Moreover, similar discretization method is applied to solve Biot consolidation problem in geomechanics. Specifically, Jha and Juanes [24] approx-imated displacement, pressure and flux using linear (P 1 ), lowest order Raviart-Thoams (RT 0 ) and element-wise constant (P 0 ) elements respectively and successfully performed reservoir simulations. Ferronato et al. [25] adopted the same discretization method in solving three dimensional Biot consolidation and showed that the mixed formulation alleviates the pressure oscillations at the interface between different permeability.
Extension from linear poroelasticity to finite strain is also investigated. Berger et al. [26] presented a three-field finite strain poroelasticity framework where the pressure and flux are approximated using P 1 /P 0 elements combined with a stabilization technique to avoid spurious pressure oscillation. In the field of biomechanics, where large deformations of biological tissues are often expected, Wall et al. [27] and Levenston et al. [28] presented their own mixed finite element formulations under the theoretical framework of mixture theory.
Motived by the success of MHFEM in solving Darcy's type equations and Biot consolidation problems, we apply MHFEM in swelling simulations in order to achieve more reliable and satisfactory results. Malakpoor et al. [29] applied MHFEM in the simulation of the swelling of cartilaginous tissues. However, his simulations were limited to small deformations. To our knowledge, there is so far no research available that apply MHFEM to the finite swelling of gels. In this study, a nonlinear system of partial differential equations modeling the swelling of SAPs particles is derived and solved using MHFEM. The model is first discretized in time using the first order implicit Euler finite difference scheme and then linearized using the Newton-Raphson strategy. Space discretization is achieved using the lowest order Raviart-Thomas space followed by a hybridization procedure. The numerical model is validated in one dimension by comparing to a semi-analytical solution obtained using MATLAB. By means of numerical examples in two dimensions, we demonstrate that MHFEM is a robust and accurate method for swelling simulations involving large deformations.

Relevant theory
In the framework of mixture theory [30], SAPs are treated as the superposition of two constituents that occupy the same physical domain: fluid phase and solid phase (a ∈ {s, f }). We assume that the fluid phase and the solid phase are nonreactive and incompressible. Body force and inertia are ignored.

Preliminaries
Let Ω s0 ∈ R 2 be the original bounded domains of the solid field. The original domain is transformed to the current (gel) domain Ω by: where x ∈ Ω, X ∈ Ω s0 , χ s is an invertible and continuously differentiable mapping from Ω s0 to Ω. As we are interested in the deformation of the solid field, without explicit indication, the original domain is referred to as Ω 0 = Ω s0 . The velocity of the solid phase is given by: Time derivative D Dt is a material time derivative with respect to the original solid field and is defined as: The deformation gradient F, right Cauchy-Green strain tensor C and volume ratio J , are defined as: We define the apparent density of a−constituent ρ a as the mass of constituent a per unit volume of the mixture, in contrast to the true density of a−constituent γ a , which is defined as the mass per unit volume of the constituent. An important dimensionless quantity in mixture theory of immiscible constituents is volume fraction φ a (x, t) := ρ a /γ a . The physical meaning of φ a (x, t) is the volume of the a−constituent per unit volume of the mixture. In the case of incompressible solid and fluid (γ a is constant), mass conservation of constituent a in the current configuration is derived as: The assumption that the gel is fully saturated yields a constraint: for x ∈ Ω, t ≥ 0. Adding up the mass conservation equation for individual constituents Eq. (7) and making use of the relation Eq. (8), the mass conservation for the mixture at the current configuration is derived: This relationship is the basis of the derivation of mass conservation in Sect. 3 (Field equations).

Helmholtz free energy function
Next, we postulate that Helmholtz free energy per initial mixture volume W consists of two independent parts: ionic part and elastic part. Namely, we have: The form of the elastic part W elastic is taken from Wilson et al. [31]. Explicitly, we have: where K is the bulk modulus and G is the shear modulus. One can recognize that it represents the compressible Neo-Hookean law. Further, we include the following relationship between Poisson's ratio ν and solid volume fraction φ f : ν = 0.5φ s . This relation is justified by the assumption that at the dry state (φ s = 1), the polymer network becomes incompressible (ν = 0.5). Making use of the relationship between the current solid volume fraction and the initial volume fraction φ s = φ s,0 J , the bulk modulus K can be rewritten in terms of G and φ s,0 as: Substituting Eq. (12) into (11), the final form of W elastic is derived: The ionic part W ionic represents the energy related to the affinity between the ions and the water. The form the ionic part is given by: where R is the universal gas constant, T is the absolute temperature, Γ is the osmotic coefficient and c is the molar concentration of the fluid phase and Φ is the fluid volume fraction in initial configuration (Φ = J φ f ). It is related to the Donnan osmotic pressure π via the relation:

Biphasic swelling model
In the biphasic swelling model [6], it is hypothesized that Donnan osmosis equilibrium is achieved instantly. Under the condition that the activity coefficient of the mobile ions is the same in the gel and in the equilibrium solution, the Donnan equilibrium ion concentrations inside the gel are derived [8]: where c f c denotes the (deformation-dependent) fixed charged density and c ex the (constant) ion concentration in the external solution. The osmotic pressure π inside of the gel follows the van 't Hoff law: where R is universal gas constant, Γ is the osmotic coefficient and T is absolute temperature.

Field equations
In Sect. 2, we have presented some preliminaries of mixture theory and the derivation of the osmotic pressure. In this section, we formulate our model problem mathematically using the related theory presented in Sect. 2. The isothermal condition, hyperelasticity and isotropy of the material are assumed. The governing equations contain linear momentum balance, fluid content balance and extended Darcy's law.

Governing equations and constitutive relations
By incorporating general balance principles (mass conservation and linear momentum conservation) and incompressibility constraints into the entropy inequality, we derive the following set of equations Eqs. (19)-(21) that a mixture continuum needs to obey. They are "pulled back" to the reference (initial) configuration [32], since we deal with finite deformation in the current work.
The divergence operator with "0" subscript indicates the divergence is taken with respect to the initial configuration. Equation (19) represents linear momentum balance in the reference configuration. T denotes the first Piola-Kirchhoff stress of the mixture and is calculated by: where p is hydraulic pressure and oe e f f is the effective stress. Assuming proper constitutive relations for the gel, the effective stress is derived from the elastic part of the Helmholtz energy function W elastic (Eq. 13) : Recall that G is shear modulus and φ s,0 is the initial solid volume fraction. Equation (20) represents Darcy's law in an extended sense. Q denotes the flux in the initial configuration and is related to the variables in current configuration by: K is the hydraulic permeability tensor written in the reference configuration and is related to the hydraulic permeability k in the current configuration by: Chemical potential μ inside of the gel consists of two parts: hydraulic pressure p and osmotic pressure π . Specifically we have [33]: where π is given in Eq. (18). Fixed charge is considered to be part of solid matrix. Its density c f c is therefore decreases as the solvent molecules enter the gel. In other words, the fixed charge density is deformation dependent: where c f c 0 denotes the initial fixed charge density and φ f ,0 denotes the initial fluid volume fraction.
At last, Eq. (21) represents the fluid content balance. It is basically Eq. (9) written in the initial configuration. Based on the molecular incompressibility assumption [34], the volume of the solid stays unchanged during swelling. We have thus: which leads to Consequently, we have:

Three-field weak formulation
We adopt mixed formulation to describe fluid permeation: position x, chemical potential μ and flux Q are chosen to be prime variables in the weak formulation. In this section, the corresponding weak form of the governing equations is presented. Firstly, we need to determine proper search/test function spaces for the approximation of the prime variables. For the position field, the continuity requirement demands x at least continuously differentiable. We have thus where Γ μ N denotes the Neumann boundary for μ. To achieve the symmetry of the system, we use the second Piola-Kirchhoff stress S and Green strain tensor E pair to calculate the virtual work. We make use of the following relation: where the bar "-" denotes that the variable is a virtual variable and we define The three-field variational form of the system is given as fol- , such that the following set of equations are satisfied . Note that f ext denotes external surface tension applied to the Neumann part boundary Γ u N . μ ext denotes the chemical potential in the external solution at the Dirichlet boundary part Γ μ D . We have μ ext = −2RT c ex .

Mixed hybrid finite element method
Given the weak form (35)-(37), we present in this section how MHFEM was applied to solve the system. The model is first discretized in time using the first order implicit Euler finite difference scheme and then linearized using the Newton-Raphson strategy. Space discretization is achieved using the lowest order Raviart-Thomas space followed by a hybridization procedure.

Time discretization and linearization
First of all, equations are discretized in time. To achieve unconditional stability, first order implicit Euler time discretization scheme is implemented. The time-discretized system is derived as follows: where subscript n/n−1 denotes variables at time step n/n−1; Δt denotes the time step size. From now on, we will skip subscript n.
The non-linearity stems from both the gel constitutive law (geometric and material non-linearity) and the nonlinear dependency of osmotic pressure on deformation (Eq. 27). Newton-Raphson strategy is applied to solve such a nonlinear system. "δ" in front of a variable indicates the increment of the variable and the tilde above the variables indicates the current (known) value of the variable so far. Linearization is performed as follows: find (δx, δQ, δμ) ∈ where for any (δx, δQ, δμ) where C is the fourth order elasticity tensor.

Spatial discretized form
Given the linearized system of equations (Eqs. 41-43), we further discretize the spatial domain Ω 0 . In this work, we assume two dimensional domain and we use four-node linear quadrilateral elements to discretize the domain. Q h denotes a quadrilateral domain discretization. Finite dimensional approximations of the search function spaces (H 1 D (Ω)) 2 , H N (div, Ω), and L 2 (Ω) are introduced. Specifically, we have: where The finite dimensional function spaces are defined as: Note that P 1 (Ω e ) denotes the space of polynomials of degree less or equal than 1 on one single element domain Ω e . We have: RT 0 (Ω e ) stands for the lowest order Raviart-Thomas element [35] and is defined on one single element Ω e as M 0 (Ω e ) denotes the function space that is constant on element Ω e . Therefore, the spatial discretizations on one single element Ω e are derived: where ϕ i and υ i are basis functions of the function spaces P 1 0 (Q h ) and RT 0 −1 (Q h ) and Ψ is an arbitrary constant. The isoparametric concept has been invoked for the basis functions ϕ i and υ i . The basis functions on the reference domain (x,ŷ) are given as: For Raviart-Thomas element basis function υ i , the transformation between the reference domain (x,ŷ) and the real domain (X , Y ) follows Piola transformation: where M is the affine mapping matrix from the reference domain to the real domain. Such a transformation guarantees that the integration of flux on each edges stays the same as in the reference domain and in the real domain [21]. Finally, based on the discretization details we presented above, the matrix form of the system on the element level is derived as: ⎛ where for i, j = 1, . . . , 4 Remark In poroelasticity, locking (which often manifests itself as spurious pressure oscillation) receives a considerable amount of attention over the years. Locking has been proved to be related to the violation of inf-sup condition for the coupling discrete element spaces [36]. Numerous work has been done to unveil the cause of locking and its numerical remedy [37][38][39][40]. As Haga et al. [40] argued that for (more than two fields) mixed formulation, the coupling between stable element spaces should satisfy individual problems. In our three-field formulation, there are two pairs of coupling in consideration: chemical potential-flux and displacement-chemical potential. For the chemical potentialflux pair, we have chosen the well-known stable pair P 0 /RT 0 for Darcy flow problem, which conveniently also possesses local mass conservation property. As to the displacementchemical potential pair, as heuristically explained by Philips and Wheeler [41], care needs to be taken at the beginning of the deformation of a porous medium when a small time step and low permeability are considered in poroelasticity. Basically, the divergence of the displacement is close to zero due to the incompressibility constraint on the solid matrix. Extension from infinitesimal strain to finite strain is straightforward.
In our formulation, the incompressibility constraint (Eqs. 28) has been implicitly included in the fluid mass conservation equation. As a matter of fact, the infinitesimal volume change of gel δ J is approximated by ∇ · x. Combining with a small time step (δt) and low permeability (k), the fluid mass conservation equation is reduced to a constraint on the position field x at the beginning stage of the swelling simulation: It is well-known that for such a constraint, the current discrete element space combination P 1 /P 0 is unstable [21,36]. Hence, to avoid locking, different element spaces combinations (for example, P 2 /P 0 or P 2 /P 1 ) or stabilization techniques need to be applied. However, since finite swelling simulations stay quite far away from incompressibility constraint except for the very beginning stage, we assess the impact of such an unstable pair on the simulation by means of numerical example (Sect. 6.1).

Hybridization procedure
So far, the search function space we applied in the discretization for the flux field is RT 0 −1 (bigger) instead of RT 0 0 (smaller). For a function that is in RT 0 −1 to be also in RT 0 0 , a necessary and sufficient condition is that the normal flow across the edge between the neighboring elements is continuous [21]. There are several ways to implement such a constraint [42]. Here we adopt the Lagrange multiplier method. By introducing a new variable λ with the physical meaning chemical potential on edges as Lagrange multiplier, the constraint mentioned above is enforced. One of the advantages of such an implementation is that by means of static condensation the total number of unknowns is eventually reduced from (x, Q, μ, λ) to (x, λ), which leads to improvement in computational efficiency. Besides, for such an implementation no extra information like edge orientations is needed and therefore is more desirable when mesh structure is complex. In what follows, we show how the system becomes "hybrid" by introducing a Lagrange multiplier λ.
Firstly, we define the discretized approximation space for λ: where L 0 0 (Q h ) denotes the set of functions that are constant on each edge e of the decomposition Q h . Next, to incorporate boundary conditions, we define which indicates that λ h has to satisfy the Dirichlet boundary condition of μ. The rationale of hybridization can be summarized in the following result: let ) be the solution of the system: for any given function J (X, t) that satisfies D J Dt (X, t) ∈ L 2 (Ω), ∀t. Then, according to [21] Theorem 1.1 p.180, there exists a unique λ h ∈ Λ μ ext ,D satisfying: is the solution of the following "hybrid" system: is also the solution of (82)-(83). Therefore, by introducing the Lagrange multiplier λ, the continuity of the normal flow is guaranteed and we are assured that the solution is found in the correct function space.
In terms of implementation, given the approximation of δλ h on one single element Ω e : where η i 's are constant on e i ∈ ∂Ω e , the system of equations (69) is extended to where After applying static condensation, we end end up solving the following (much smaller) system: where Note that A 1 and A 3 are proven to be symmetric positive definite matrices [29].

Solution verification
So far, we have presented various aspects of the numerical simulation engine. In this section, we focus on the verification of such a simulation engine. Namely, given the partial differential equation system we would like to solve, we calculate the solution in a different way and then compare results with the proposed simulation results. Since there is no analytical solution available for the transient swelling or consolidation simulation with finite deformation even in one dimension. Simulation results are compared with semi-analytical results. These semi-analytical results are solutions calculated by a MATLAB internal partial differential equation solver ("pdepe") in one dimension. On the other hand, the equilibrium state of a circular swollen gel can be calculated analytically (homogeneous swelling) we also compared our simulation results with that.
First of all, the system of nonlinear partial differential Eqs. (19)- (20) are simplified in 1D (y-direction). Basically, deformation tensor F is reduced to volume ratio J , and the simplified equations in terms of J and μ are: where σ e f f and π are nonlinear function of J , given by: (104) Equation (101) expresses the linear momentum balance for the swelling case (no surface compression/traction applied).
In the case of consolidation, Eq. (101) needs to be adapted to: where f ext is the surface tension applied. Equation (102) is the fluid content conservation equation with Darcy's law incorporated. Substituting μ in Eq. (102) using (101), a parabolic (nonlinear) partial differential equation of J is derived: The initial condition (for consolidation as well as swelling) is: Boundary conditions for consolidation simulation are: Their physical meaning are: at y = 1, no flow boundary condition; at y = 0 the sample is in touch with the outer solution whose chemical potential is kept constant and equals −2RT c ex 0 . The value of c ex 0 is chosen so that at the initial state the gel is at "stress-free" state. Therefore, c ex 0 needs to satisfy the following equation: Similarly, boundary conditions for swelling simulation are: Their physical meanings are similar to the one for the consolidation simulation. Parameters given in Table 1 are chosen carefully so that they are within the industrially relevant regime. These parameters are used throughout all simulations presented in this paper unless otherwise indicated. In this section we set c ex to be physiological salt concentration, thus equals 1.54 × 10 −4 mol/ml. The applied (downward) pressure f ext in the consolidation simulation is 0.02 MPa. Mesh size and time steps are chosen to be the same for MHFE and MAT-LAB solutions.
The evolution of swelling ratio and dimensionless chemical potential on the top surface are plotted for both consolidation and swelling simulations (Figs. 2, 3). The characteristic time scale τ in this problem equals l 2 m p /RT k, where l is the characteristic length (dry size), m p is the molar volume of the polymer (taken to be 105 cm 3 /mol). The dimensionless time is found to be t d = t/τ . Similarly, the dimensionless chemical potential is derived as μ d = μ/( RT m p ). Figure 2a shows that the chemical potential at the top surface monotonically decreases over time with some delay at the beginning of the consolidation simulation. This delay can be explained by the fact that it takes time for the fluid inside the gel to drain through the bottom (y = 0). We conclude from the figures that MHFEM solutions match MATLAB solutions very well for both variables of interest in both simulations.
For homogeneous swelling in two dimensions, the following equations of stretch γ must hold: If we let res(γ ) = σ e f f (γ ) − π(γ ) + 2RT c ex , its magnitude stays at 10 −10 for γ calculated by MHFEM under various values of external salt concentrations (c ex ).

Numerical examples
The capability of the MHFEM in the swelling simulations is demonstrated in this section by four two-dimensional numerical examples. We start with the swelling with a low-permeable stripe simulation, in which a heterogeneous permeable domain is introduced. By comparing to the standard FEM, the impact of local mass conservation property is highlighted. Besides, the locking phenomena that are frequently encountered in poroelasticity problems are discussed. Next, we move on to present the simulation results of free swelling of a quarter of a square, which followed by the free swelling of a quarter of a circle using MHFEM. In both cases, the swelling ratios are above 30. At last, we present a complex pattern formation induced by swelling. The results produced by the simulation are compared with the experimental results reported in the literature qualitatively.

Swelling with a low-permeable stripe
Inspired by a numerical example presented in [22], a heterogeneously permeable domain is introduced (Fig. 4).    Simulations are carried out using both MHFEM and standard FEM with different mesh sizes (0.0833 mm, 0.0417 mm, 0.0167 mm, 0.0083 mm). Note that regardless of mesh refinement, the area and location of the low permeable stripe is kept unchanged. Due to the existence of the low perme-ability stripe, the top edge EF is expected to form a curve as the swelling starts and the curve is used to characterize the deformation of the gel (Figs. 5, 6). We conclude that both methods show convergence (to the same curve) as the mesh size decreases. Moreover, with a Both Tables 2 and 3 show that MHFEM produces more accurate results than FEM given the same mesh sizes. As a matter of fact, the average error indicates that MHFEM produced more accurate results than FEM with the mesh 2.5 times coarser when shear modulus G = 0.15 MPa. The accuracy advantages of MHFEM over FEM is more pronounced when the shear modulus is reduced. Table 3 shows that at mesh sizes 0.0833 mm and 0.0417 mm the maximum error of FEM is almost twice as much as the one of MHFEM.
As it is argued that MHFEM possesses local mass conservation property, which should lead to more accurate calculation of flux, we plot the flux vectors and corresponding streamlines originating from the bottom edge for both methods at a given moment Fig. 7) in order to trace back the deformation differences presented above. Note that in standard FEM chemical potential is approximated by linear node-based interpolation functions, as a result of which the flux vector is element-wise constant. As to MHFEM, fluxes are derived edge-wise. Adding up flux on each edge, the net flux of a given element is derived for MHFEM solution and thus also element-wise constant.
A close check on the streamlines near the low permeability area reveals that the streamlines are not completely symmetric in the FEM solution although geometry, permeability distribution and boundary conditions are strictly symmetric (Fig. 7a). On the other hand, the streamlines in MHFEM solution stay strictly symmetric (Fig. 7b). The difference in streamline distribution leads to visible deformation difference. We notice that elements at the corner of the low permeability area in Fig. 7a exhibit non-physical distortions; while the MHFEM solution (Fig. 7b) presents a much more physical deformation. Such non-physical distortions not only affect the accuracy of deformation calculation but also cause converging issues. In fact, FEM simulation failed to converge to the next step in this simulation while MHFEM succeeded. In our simulation practice, it is observed that FEM is more demanding than MHFEM in terms of choosing modeling and discretization parameters to reach convergence.
To investigate the impact of an unstable pair (P 1 /P 0 ) on the swelling simulation, we slightly adapt boundary conditions in Fig. 4. Instead of edge HG in touch with the outer solution, we let EF directly in touch and the rest of edges are with no flow conditions. Also the domain size is changed to 1 mm × 1 mm. This way we generate a sharp gradient in chemical potential especially for the low permeability area. Fig. 8 shows the chemical potential contour of both an unstable pair (P 1 /P 0 ) and a stable pair (P 2 /P 0 ) for displacement-chemical potential coupling. We observe that indeed the unstable pair led to a checkerboard distribution  of pressure while the stable pair not. Also due to the higher degree of approximation for displacement, the deformations generated by the two pairs are slightly different. However, we notice that the oscillation in the pressure field dissipates as the swelling goes on and eventually disappears. Local element refinement also helps with the alleviation of checkerboard distribution. Moreover, the use of higher order element also causes higher computational cost. Hence, unless great impor- tance is attached to the calculation of chemical potential at the initial stage, we keep using P 1 /RT 0 /P 0 pair in our swelling simulations potentially combined with local mesh refinement scheme.

Free swelling of a quarter of a square
In this section, we apply MHFEM to simulate the free swelling of a quarter of a square OABC (Fig. 9). The bound-ary conditions prescribed for free swelling are given as follows: where μ ext denotes the chemical potential in the external solution and equals −2RT c ex ; Γ u D represents OA in ydirection and OC in x-direction; Γ u N contains AB and BC; Γ μ D contains AB and BC; Γ μ N contains OA and OC. The initial conditions are: The outer solution concentration c ex is set to be 4.25×10 −5 mol/ml. The domain is discretized by 20 × 20 elements and we use an adaptive time step size.
The chemical potential contour plots over time are given in Fig. 10. The color indicates the magnitude of the chemical potential. It shows that fluid permeation starts from the two edges (AB and CB) gradually reaching the core. The square shape of OABC is temporarily lost due to the faster swelling area near node B comparing to node A and C. At the equilibrium state (t d → ∞) the square shape is recovered (Fig. 10).
To predict material failure (for example, the initiation of fractures), magnitude, directions and locations of the maximum and minimum stresses must be identified. Between the two parts of stress (effective stress and hydraulic pressure), only the effective stress contributes to the initiation of fractures. Using the proposed numerical model, the magnitude of the maximum and minimum normal stresses and maximum shear stress distribution over the sample can be calculated (Fig. 11). It is observed that the maximum normal stress appears near node A and C along the edges that are in touch with the outer solution. Similar observation can be made for minimum normal stress. The maximum normal stresses are tensile stress (positive values) and the minimum stresses are compressive stresses (negative values). As to the maximum shear stresses, near nodes C and A the sample sustains the highest shear stresses and near origin O and node B the lowest.
The convergence of the simulation as the mesh sizes decreases is demonstrated in  16 Swelling of a gel membrane with periodical holes: chemical potential contour plots at three moments. a t = 100 s. b t = 300 s. c t = 500 s yielded by different meshes are given in the table (e position and e f lux ). As the mesh size decreases, the errors converges towards zero.

Free swelling of a quarter of a circle
Next, we investigate gels with a circular shape. Only a quarter of the circle (OAB) is simulated due to symmetry reasons (Fig. 12). The boundary conditions are the same as for the square except that Γ u D represents OA in x-direction and OB in y-direction; Γ u N and Γ μ D represent curve AB; Γ μ N contains OA and OB. The domain is discretized by 163 quadrilateral elements with time step taken to be same as for the square.
Unlike in the square simulation, the shape of OAB stays unchanged during swelling (Fig. 13). The extremal effective stresses plot (Fig. 14) shows that at the transient state the maximum normal and shear stresses are largest along the outer boundary and gradually decrease towards the core. The opposite holds true for minimum normal stress, where the largest minimum normal stress is at the core and decreases as the radius increases.

Swelling-induced bifurcation
A swelling bifurcation experiment is reported by Zhang et al. [43]. By exposing a hydrogel membrane with periodic circular hole array to a solvent, the holes deform into ellipses and the interaction between them yield a "diamond plate" pattern. In this subsection, we apply the proposed numerical model to simulate the swelling of such a membrane aiming to replicate the deformation of holes qualitatively as reported in the experimental work by Zhang et al. [43].
The domain (Fig. 15) is discretized by 2926 linear quadrilateral elements and the time step is 0.025 s. The outer solution concentration c ex set to be 5.54 × 10 −4 mol/ml. The chemical potential contour plots (Fig. 16) show that the deformation of the originally circular holes. Firstly, the circular holes collapsed into ellipses of alternating directions (Fig. 16a). Next, as the swelling went on, the elliptic holes further elongated themselves in their major axis direction (Fig. 16b). At last, a narrowing "waist" of the elongated slits appeared (Fig. 16c). The simulation using the developed numerical model goes on until contact between the elements happens, since contact was not defined in the numerical model. However, the calculated deformations of the holes by MHFEM model are in good agreement with the experimental results.

Conclusions
In this study, we have developed a MHFE model under the theoretical framework of mixture theory for the simulation of swelling ionized hydrogels in two dimensions. Newton-Raphson strategy was applied to handle the non-linearity arising from hyperelasticity and nonlinear osmotic pressure term. We deployed the lowest order Raviart-Thomas element to approximate flux as an independent variable. Then hybridization procedure was introduced to guarantee the continuity of normal flux across neighboring elements so that the lowest order Raviart-Thomas elements were properly implemented. The calculated solution using proposed model was verified by comparing with semi-analytical solutions calculated by MATLAB in one dimension.
Local mass conservation property of the proposed model guarantees more accurate calculation of flux and deformation, which is crucial for the simulation of extremely large deformation induced by swelling as shown in the lowpermeable stripe simulations. Next, we continue to simulate the free swelling of a square and a circular-shaped gel using MHFEM. In both simulations, the swelling ratio is more than 30. Chemical potential contours and extremal stresses distributions are presented. Chemical potential contours give us a good idea about the fluid permeation progress during swelling and extremal stress distribution plots are useful for failure mechanics. At last, we carried out a swellingbifurcation simulation. The simulation results are verified against experimental results and showed good agreement qualitatively.