1 Introduction

Hydrogels swell by the permeation of solvent (water) molecules into polymer networks. Due to the existence of cross-links in the polymer networks, water permeation leads to the expansion and rearrangement of polymer chains as a whole. There are various external stimuli that induce the swelling of a hydrogel: temperature, pH value and ionic strength in the outer solution [39]. The responsiveness and bio-compatibility of hydrogels make them popular materials for diverse applications including carriers for drug delivery, tissue engineering matrices [21], actuator and sensors [2].

Over the years, many works are dedicated to the development of theoretical frameworks that enable the coupling between fluid permeation and solid deformation. Following the linear and finite poroelasticity approach established by Biot [4, 5], Hong et al. [23] proposed a theoretical framework that enables finite deformations with combination of Flory theory [19] to model the swelling of a polymeric gel. Many more works have been done to extend this framework to include features like non-Gaussian constitutive relations for solid [12], viscoelasticity [11, 43] and temperature-sensitivity [13, 16]. Besides the thermodynamic approach, which treats the gel as one single phase, mixture theory [8], which allows the study of multiple phases (fluid, solid, ions and other substances) at the same time, gains popularity in the biomechanical studies. A biphasic swelling model is proposed by Lanir [31] assuming the existence of fixed charge groups and the instantaneous equilibrium of the ion phase. Later, Lai et al. [29] developed a triphasic theory in which fluid, solid and ion phases were considered separately. A further separation of cation and anion phases leads to the quadriphasic theory [25].

Development in the theoretical frameworks fostered many numerical investigations into the swelling process. Hong et al. [22] and Kuang and Huang [27] simulated the inhomogeneous swelling of a gel at equilibrium state. Later on, transient swelling simulations are the topic of many studies [7, 14, 32, 34, 40, 46]. We note that all the studies listed above adopted standard FEM in which deformation and chemical potential were chosen to be prime variables implemented in commercial softwares like ABAQUS or COMSOL. Böger et al. [6] proposed an alternative variational formulation which treated deformation and flux as prime variables. On the other hand, studies [3, 23] on the solvent molecules migration kinetics suggest that the mobility tensor (the coefficient tensor in the Darcy’s equation) is not homogeneous during transient swelling. Basically, the part of a gel that gets in touch with the external solvent earlier increases dramatically its hydraulic permeability. This inhomogeneity in the mobility tensor may require special attention in terms of numerical treatment. As pointed out by Kaasschieter and Huijben [26] and Mosé et al. [36], in a Darcy flow problem with heterogeneous domain, the standard FEM suffers poor accuracy in computing flux field due to the violation of local mass conservation. In contrast, a mixed formulation that treats flux and pressure as separate independent variables greatly improves the accuracy of the flux field as mass conservation is preserved locally and globally. As a matter of fact, using proper numerical methods to solve Darcy flow problem is a well-established topic in porous media studies [15, 33]. Inspired by the success of mixed formulations in Darcy flow problems, a mixed hybrid finite element method (MHFEM) was implemented to simulate the transient swelling of a hydrogel in one of our previous works [45].

In this work, we aim to evaluate the effectiveness of MHFEM comparing with standard FEM in terms of robustness, accuracy and efficiency the solution method when applied to simulate the finite transient swelling of hydrogels. As the theoretical basis of the MHFEM is elaborated in [45], in this work, the focus lies on comparing and furthermore demonstrating the advantages of MHFEM over standard FEM in swelling simulations involving large deformations. Results and discussions of numerical examples form the key part of this work. Depending on the hydraulic permeability laws we apply, our numerical examples fall into three categories: homogeneous permeability (the mobility tensor is constant), heterogeneous permeability (the mobility tensor is a strain-dependent, thus non-uniform but change continuous during transient swelling) and heterogeneous with interface of discontinuous permeability (mobility tensor is thus discontinuous). The effectiveness of MHFEM is evaluated for all three categories. The paper is organized as follows. In Sect. 2, we recap the swelling mechanism and the resulting governing equations. In Sect. 3 and 4, the weak form and discretization details of standard FEM and MHFEM are described respectively. Section 5 presents the comparison results and discussions of four selected numerical examples. At last, conclusions are drawn in Sect. 6.

Fig. 1
figure 1

SAP particle structure illustration

2 The swelling mechanism and governing equations

We limit ourselves to one specific type of hydrogels: superabsorbent polymers (SAP) or in other words, partially neutralized sodium polyacrylate hydrogels. A schematic illustration of its structure is given in Fig. 1. Once the particle gets in touch with the outer solution, the water interacts (1) with the hydrophillic polychains and (2) with the solium ions electrostatically linked to those polychains. The latter interaction is represented by the ionic part of the free energy while the former is represented by the mixing part of the free energy. Moreover, it has been suggested that in the case of monovalent solvent, the contribution from the mixing energy can be ignored [24].

We make the following assumptions in the swelling model. The gel is placed in isothermal condition and responsive only to the ion concentration change in the outer solution due to Donnan osmosis. Both the external and the solid skeleton are assumed to be incompressible and the gel as a whole is assumed to be hyperelastic. The ion equilibrium in and outside of the gel is established instantly, as a consequent of which the biphasic swelling model [31] applies.

The swelling model consists of two governing equations (linear momentum balance and (fluid) mass conservation) plus constitutive relation (a Darcy type equation which relates flux and chemical potential gradient). Since deformation of the solid field is of our interest and is in the finite deformation regime, the equations are given in Lagrangian form. First of all, the linear momentum balance in Lagrangian form is given as:

$$\begin{aligned} \nabla _X\cdot {{\varvec{T}}}={\varvec{0}}, \end{aligned}$$
(1)

where \({{\varvec{T}}}\) is the first Piola-Kirchoff stress tensor. It is related to Cauchy stress tensor \({\varvec{\sigma }}\) by the relation: \({{\varvec{T}}}=J{\varvec{\sigma }}{\varvec{F}}^{-T}\), where J is the volume ratio \(J=\det ({\varvec{F}})\) and \({\varvec{F}}\) is the deformation tensor. Cauchy stress in a hydrogel following the convention in poroelasticity is given by:

$$\begin{aligned} {\varvec{\sigma }}={\varvec{\sigma }}_{eff}-p{\varvec{I}}, \end{aligned}$$
(2)

where \({\varvec{\sigma }}_{eff}\) denotes the effective stress and p stands for the hydraulic pressure in the gel. We further specify the concrete form of \({\varvec{\sigma }}_{eff}\) following the work of Wilson et al. [44]. Although the solid skeleton and fluid phase are both assumed to be incompressible, the gel experiences large volume change by imbibing external fluid. Therefore, the material law of the gel is assumed to be compressible Neo-Hookean. Besides, we make the assumption that Poisson ratio \(\nu \) is proportional to the solid volume fraction \(\phi _s\), \(\nu =0.5\phi _{s}\). This assumption ensures that when there is no fluid (\(\phi _s=1\)), the incompressibility of the solid skeleton is recovered. Following the derivation presented in [44], the final form of \({\varvec{\sigma }}_{eff}\) is given by:

$$\begin{aligned}&{\varvec{\sigma }}_{eff}=-\frac{G}{6}\frac{\ln (J)}{J}\left[ -1 +\frac{3(J+\phi _{s,0})}{(-J+\phi _{s,0})}+\frac{3\ln (J)J\phi _{s,0}}{(-J+\phi _{s,0})^2}\right] {\varvec{I}}\\ \nonumber&+\frac{G}{J}({\varvec{F}}\cdot {\varvec{F}}^T-J^{2/3}{\varvec{I}}), \end{aligned}$$
(3)

where G is shear modulus, which indicates the softness of the gel and \(\phi _{s,0}\) denotes the initial solid volume fraction and related to \(\phi _s\) by \(\phi _s=\phi _{s,0}/J\). The other part of Cauchy stress, p is derived by the relation [3]:

$$\begin{aligned} p=\mu +\varPi , \end{aligned}$$
(4)

where \(\mu \) denotes the chemical potential and \(\varPi \) is the osmotic pressure due to Donnan osmosis. Assuming biphasic swelling model [31], the osmotic pressure is related to fixed charge density \(c^{fc}\) as:

$$\begin{aligned} \varPi =RT\sqrt{(c^{fc})^2+4{\bar{c}}}, \end{aligned}$$
(5)

where R and T are universal gas constant and absolute temperature respectively; \({\bar{c}}\) is the external ion concentration. The current fixed charge density \(c^{fc}\) is related to the initial fixed charge density \(c^{fc}_0\) by the relation: \(c^{fc}=c^{fc}_0\phi _{f,0}/(J-\phi _{s,0})\) [44].

The second equation represents mass balance of the fluid phase:

$$\begin{aligned} \frac{D\varPhi _f}{D t}+\nabla _X\cdot {\varvec{Q}}=0, \end{aligned}$$
(6)

where D/Dt denotes the material time derivative of a given material field variable. \(\varPhi _f\) denotes the nominal porosity and related to the porosity at the current configuration \(\phi _f\) by \(\varPhi _f=J\phi _f\). \({\varvec{Q}}\) denotes the flux vector in Lagrangian form. Due to the incompressibility conditions, the conservation of the solid phase can be written as:

$$\begin{aligned} \phi _{s,0}=1-\phi _{f,0}=J(1-\phi _f)=J-\varPhi _f. \end{aligned}$$
(7)

Therefore, Eq. (6) can be rewritten as:

$$\begin{aligned} \frac{D J}{D t}+\nabla _X\cdot {\varvec{Q}}=0. \end{aligned}$$
(8)

At last, we relate the flux Q with the chemical potential \(\mu \) with a Darcy’s type equation:

$$\begin{aligned} {\varvec{Q}}=-{\varvec{K}}\nabla _X\mu , \end{aligned}$$
(9)

where \({\varvec{K}}\) is the Lagrangian form of the intrinsic permeability k. They are related to each other by \(\varvec{K}=J{\varvec{F}}^{-1}{\varvec{F}}^{-T}k/\eta \) [35], where \(\eta \) is the viscosity of the fluid. The intrinsic permeability k is assumed to be constant, strain-dependent or discontinuous in this work. The specific values of k and strain-dependent forms of k are discussed in Sect. 5. Unlike the original Darcy’s law, Eq. (9) features the chemical potential instead of pressure. It has two contributions: hydraulic pressure p and osmotic pressure \(\varPi \), as given in Eq. (4).

3 Two-field standard finite element formulation

Standard finite element formulation is by far the most widely used formulation in swelling simulations. Commercial software such as ABAQUS or COMSOL is used in combination with user-defined subroutines. Regardless of the swelling mechanisms and constitutive relations of the gel, flux field \({\varvec{Q}}\) in the mass conservation equation (8) is substituted using Darcy type equation (9) in this formulation, yielding a system of equations of \(({\varvec{x}},\mu )\):

$$\begin{aligned} \nabla _X\cdot {{\varvec{T}}}&=0, \end{aligned}$$
(10)
$$\begin{aligned} \frac{D J}{D t}-\nabla _X\cdot ({\varvec{K}}\nabla _X \mu )&=0. \end{aligned}$$
(11)

Due to geometric and material nonlinearity in the material law for the gel, Newton–Raphson procedure is carried out. Therefore, instead of solving \(({\varvec{x}}, \mu )\) directly, its linear increment \((\delta {\varvec{x}},\delta \mu )\) is solved from the linearized version of Eqs. (10) and (11). Let \(\varGamma _D^x\) and \(\varGamma _D^\mu \) denote the Dirichlet boundary for \({\varvec{x}}\) and \(\mu \). We define the function space \(H_D^1(\varOmega )\):

$$\begin{aligned} H_D^1(\varOmega )=\{\nu \in L^2(\varOmega ):\frac{\partial \nu }{\partial x}\in L^2(\varOmega ),\nu \vert _{\varGamma _D^{x/\mu }=0}\}. \end{aligned}$$
(12)

The linearized weak form of Eqs. (10) and (11) is thus derived: find \((\delta {\varvec{x}},\delta \mu )\in ((H_D^1(\varOmega ))^2,H_D^1(\varOmega ))\), such that for any \((\bar{{\varvec{x}}}.{\bar{\mu }})\in ((H_D^1(\varOmega ))^2,H_D^1(\varOmega ))\),

$$\begin{aligned} a(\delta {\varvec{x}},\bar{{\varvec{x}}})+b(\bar{{\varvec{x}}},\delta \mu )&=R_1({\varvec{x}},\mu ), \end{aligned}$$
(13)
$$\begin{aligned} \frac{D}{D t}b(\delta \varvec{x},{\bar{\mu }})+c(\delta \mu ,{\bar{\mu }})&=R_2({\varvec{x}},\mu ), \end{aligned}$$
(14)

with

$$\begin{aligned} a(\delta {\varvec{x}},\bar{{\varvec{x}}})&=\int _\varOmega \nabla \bar{{\varvec{x}}}{\varvec{S}}:\nabla \delta {\varvec{x}}+\bar{{\varvec{E}}}:{\mathbb {C}}:\delta {\varvec{E}} d\varOmega , \end{aligned}$$
(15)
$$\begin{aligned} b({\bar{x}},\delta \mu )&=-\int _{\varOmega }J{\varvec{F}}^{-1}:\nabla \bar{{\varvec{x}}}\delta \mu d\varOmega , \end{aligned}$$
(16)
$$\begin{aligned} c(\delta \mu ,{\bar{\mu }})&=-\int _\varOmega ({\varvec{K}}\nabla \delta \mu )\cdot \nabla {\bar{\mu }} d\varOmega , \end{aligned}$$
(17)
$$\begin{aligned} R_1({\varvec{x}},\mu )&=-\int _\varOmega {\varvec{S}}:\bar{{\varvec{E}}}d\varOmega ,\end{aligned}$$
(18)
$$\begin{aligned} R_2({\varvec{x}},\mu )&=\frac{D}{D t}\int _\varOmega J{\bar{\mu }}d\varOmega +\int _\varOmega \varvec{K}\nabla \mu \cdot \nabla {\bar{\mu }}d\varOmega , \end{aligned}$$
(19)

where \({\varvec{S}}\) is the second Piola-Kirchoff stress tensor; \({\mathbb {C}}\) is the fourth-order elasticity tensor; \(\bar{{\varvec{E}}}\) and \(\delta {\varvec{E}}\) are the Green strain tensor based on \(\bar{{\varvec{x}}}\) and \(\delta {\varvec{x}}\). Specifically, we have:

$$\begin{aligned} \bar{{\varvec{E}}}&=\frac{1}{2}(\nabla \bar{{\varvec{x}}}^T{\varvec{F}}+{\varvec{F}}^T\nabla \bar{{\varvec{x}}}), \end{aligned}$$
(20)
$$\begin{aligned} \delta {\varvec{E}}&=\frac{1}{2}(\nabla \delta {\varvec{x}}^T{\varvec{F}}+{\varvec{F}}^T\nabla \delta {{\varvec{x}}}). \end{aligned}$$
(21)

The finite element method is introduced by a partition \({\mathcal {T}}_h=\{\varOmega ^e\}\) of the initial domain \(\varOmega \). A natural choice for the discretized sub-spaces for \(H_0^1(\varOmega )\) is

$$\begin{aligned} P_0^1({\mathcal {T}}_h)=\{\varphi \in L^2(\varOmega ):\varphi \vert _{\varOmega ^e}\in P^1(\varOmega ^e)\}\cap H_0^1(\varOmega ), \end{aligned}$$
(22)

where \(P^1(\varOmega ^e)\) denotes constant or linear polynomial on \(\varOmega ^e\). Let \(\varphi _i\) be the basis interpolation functions in \(P_0^1({\mathcal {T}}_h)\), we have:

$$\begin{aligned} \delta {\varvec{x}}&=\sum \varphi _i {\varvec{x}}_i, \end{aligned}$$
(23)
$$\begin{aligned} \delta \mu&=\sum \varphi _i\mu _i, \end{aligned}$$
(24)

where i runs over all the nodes introduced by the spatial discretization and \({\varvec{x}}_i\) and \(\mu _i\) denotes the value of the position and chemical potential at node i.

To avoid any instability due to time discretization, we choose backward Euler scheme for time discretization. Let subscript \(n/n-1\) denote the variable at time step \(n/n-1\) and \(\varDelta t\) the time step size, the fully discretized version of Eqs. (13) and (14) is:

$$\begin{aligned} a_n(\delta {\varvec{x}},\bar{{\varvec{x}}})+b_n(\bar{{\varvec{x}}},\delta \mu )&=r^1_n({\varvec{x}},\mu ), \end{aligned}$$
(25)
$$\begin{aligned} b_n(\delta {\varvec{x}},{\bar{\mu }})+\varDelta t\;c_n(\delta \mu ,{\bar{\mu }})&=\varDelta t\;r^2_n(\varvec{x},\mu )+b_{n-1}(\delta {\varvec{x}},{\bar{\mu }}). \end{aligned}$$
(26)

For each time step, flux field \({\varvec{Q}}\) is computed by numerical differentiation of the chemical potential field. As a result, the fluid mass balance is only satisfied in a weak sense. At the element level, this leads to the violation of mass conservation.

4 Three-field mixed hybrid finite element formulation

To circumvent the violation of local mass conservation, a mixed formulation is proposed. In the mixed formulation, flux field and the chemical potential are solved simultaneously. Raviart–Thomas elements [37] guarantees the continuity of flux in the normal direction between neighboring elements; therefore, mass conservation is preserved locally and globally. Further, to implement Raviart–Thomas elements, a Lagrange multiplier \(\lambda \) is introduced [1]. Finally, a static condensation technique is applied which reduces the number of degree of freedom from four to two. To begin with, the system of equations of \((\varvec{x},\mu ,{\varvec{Q}})\) is given by:

$$\begin{aligned} \nabla _X\cdot {{\varvec{T}}}&=0, \end{aligned}$$
(27)
$$\begin{aligned} {\varvec{Q}}&=-{\varvec{K}}\nabla _X \mu , \end{aligned}$$
(28)
$$\begin{aligned} \frac{D J}{D t}+\nabla _X\cdot {\varvec{Q}}&=0. \end{aligned}$$
(29)

We introduce the function space \(H_0(\text {div},\varOmega )\):

$$\begin{aligned} H_0(\text {div},\varOmega ){=}\{{\varvec{q}}\in (L^2(\varOmega ))^2, \text {div}{\varvec{q}}\in L^2(\varOmega ), \varvec{q}\cdot {\varvec{n}}\vert _{\varGamma _N^\mu }{=}0 \}, \end{aligned}$$
(30)

where \(\varGamma _N^\mu \) denotes the Neumann boundary for the chemical potential. Like in the standard FEM, the linearized weak form of Eqs. (27), (28) and (29) is derived: find \((\delta {\varvec{x}},\delta \mu , \delta \varvec{Q})\in ((H_D^1(\varOmega ))^2,L^2(\varOmega ),H_0(\text {div},\varOmega ))\),

$$\begin{aligned} a(\delta {\varvec{x}},\bar{{\varvec{x}}})+b(\bar{{\varvec{x}}},\delta \mu )&=r^1({\varvec{x}},\mu ), \end{aligned}$$
(31)
$$\begin{aligned} d(\delta {\varvec{Q}},\bar{{\varvec{Q}}})+e(\delta \mu ,\bar{{\varvec{Q}}})&=r^3({\varvec{Q}},\mu ), \end{aligned}$$
(32)
$$\begin{aligned} \frac{D}{D t}b(\delta {\varvec{x}},{\bar{\mu }})+e(\delta \varvec{Q},{\bar{\mu }})&=r^4({\varvec{x}},{\varvec{Q}},\mu ), \end{aligned}$$
(33)

for any \((\bar{{\varvec{x}}},{\bar{\mu }},\bar{{\varvec{Q}}})\in ((H_D^1(\varOmega ))^2,L^2(\varOmega ),H_0(\text {div},\varOmega ))\). The form of ab and \(r^1\) has been given in Eqs. (15), (16) and (18). Here we give the specific form of \(d,e,r^3\) and \(r^4\):

$$\begin{aligned} d(\delta {\varvec{Q}},\bar{{\varvec{Q}}})&=\int _\varOmega {\varvec{K}}^{-1}\delta {\varvec{Q}}\cdot \bar{{\varvec{Q}}}d\varOmega , \end{aligned}$$
(34)
$$\begin{aligned} e(\delta \mu ,\bar{{\varvec{Q}}})&=\int _\varOmega \nabla \cdot \bar{{\varvec{Q}}}\delta \mu d\varOmega , \end{aligned}$$
(35)
$$\begin{aligned} r^3({\varvec{Q}},\mu )&=-\int _\varOmega {\varvec{K}}^{-1} {\varvec{Q}}\cdot \bar{{\varvec{Q}}}d\varOmega \nonumber \\&\quad +\int _\varOmega \mu \nabla \cdot \bar{{\varvec{Q}}}d\varOmega +{\tilde{\mu }}\int _{\varGamma _N^\mu }\bar{{\varvec{Q}}}\cdot {\varvec{n}}{\varvec{d}}\varGamma , \end{aligned}$$
(36)
$$\begin{aligned} r^4({\varvec{x}},{\varvec{Q}},\mu )&=\frac{D}{D t}\int _\varOmega {\bar{\mu }}Jd\varOmega \nonumber \\&\quad +\int _\varOmega {\bar{\mu }}\nabla \cdot \varvec{Q}d\varOmega , \end{aligned}$$
(37)

where \({\tilde{\mu }}\) is the prescribed value for chemical potential at the Dirichlet boundary condition, \({\varvec{n}}\) is the outward normal vector at \(\varGamma _N^\mu \). The spatial discretization spaces are chosen to be \((\delta {\varvec{x}}_h,\delta \mu _h,\delta {\varvec{Q}}_h)\in (P_0^1({\mathcal {T}}_h), M_{-1}^0({\mathcal {T}}_h), RT_0^0({\mathcal {T}}_h))\), where

$$\begin{aligned} M_{-1}^0({\mathcal {T}}_h)&=\{l\in L^2(\varOmega ):l\vert _{\varOmega ^e}=c,\forall c\in {\mathbb {R}}\}, \end{aligned}$$
(38)
$$\begin{aligned} RT_0^0({\mathcal {T}}_h)&=\{{\varvec{v}}\in H_0(\text {div},\varOmega ):{\varvec{v}}\vert _{\varOmega ^e}\in RT^0({\mathcal {T}}_h)\}. \end{aligned}$$
(39)

Note that the combination of approximation spaces \(M_{-1}^0({\mathcal {T}}_h)\) and \(RT_0^0({\mathcal {T}}_h)\) are proven to satisfy the LBB condition and are thus a stable combination [9]. However, to facilitate the implementation, we first consider a bigger function space \(RT_{-1}^0({\mathcal {T}}_h)=\{{\varvec{v}}\in L^2(\varOmega ):\varvec{v}\vert _{\varOmega ^e}\in RT^0({\mathcal {T}}_h)\}\) and later limit it to the desired space \(RT_{0}^0({\mathcal {T}}_h)\) by introducing a Lagrange multiplier. This way, we do not need to design certain procedure to make sure the orientation of the edges are compatible. The interpolation of the independent variables is given by:

$$\begin{aligned} \delta {\varvec{x}}&=\sum \varphi _i{\varvec{x}}_i, \end{aligned}$$
(40)
$$\begin{aligned} \delta \mu&=\sum l_k\mu _k, \end{aligned}$$
(41)
$$\begin{aligned} \delta {\varvec{Q}}&=\sum {\varvec{v}}_j Q_j, \end{aligned}$$
(42)

where ik and j loop over each node, element and edge respectively. Again, we use backward Euler as time discretization scheme for stability consideration. The procedure of introducing the Lagrange multiplier and static condensation is known as hybridization technique. Its relevant details can be found in [45]. The result of hybridization is the simple and efficient implementation of lowest-order Raviart–Thomas elements with a reduction of degrees of freedom. In the end, the system of linear equations we solve is:

$$\begin{aligned} \left[ \begin{array}{cc} {\varvec{A}}_1 &{} {\varvec{A}}_2\\ {\varvec{A}}_2^T&{}{\varvec{A}}_3 \end{array}\right] \left[ \begin{array}{c} \delta {\varvec{x}}\\ \delta \lambda \end{array}\right] =\left[ \begin{array}{c} {\varvec{F}}_1\\ {\varvec{F}}_2 \end{array}\right] . \end{aligned}$$
(43)

The chemical potential and flux field can be reconstructed from \({\varvec{x}}\) and \(\lambda \) without loss of accuracy.

Table 1 Model parameters
Fig. 2
figure 2

Initial state geometry illustration

Table 2 FEM solution convergence performance summary
Table 3 MHFEM solution convergence performance summary

5 Numerical examples

In this section, four numerical examples solved by standard FEM as well as MHFEM are presented. Depending on different permeability distributions over the domain, the four examples can be categorized into three cases. In the first example, we consider constant permeability k over the whole domain. The second example treats permeability k to be strain-dependent and continuously varying its value over the domain as soon as the swelling starts. Example 3 and 4 deal with permeability with discontinuous jumps. Specifically, the swelling is confined to one direction in the third example and a low permeability stripe is added to the geometry. Example 4 allows swelling in two dimensions and represents a soft high-permeable core surrounded partly by a stiff low-permeable shell structure.

For all four examples, comprehensive comparisons are carried out between the two numerical methods in terms of and robustness, accuracy and computational cost. First, we look into the solution convergence behavior of both methods. In other words, the robustness of both methods for a given range of material and simulation parameters is investigated. Next, we examine the accuracy of the deformation field by means of comparing mesh convergence rates (provided that solution convergence is achieved). In addition, we present the contour plots of the chemical potential with streamlines on top to gain good insights on the accuracy of the chemical potential and flux field calculation at transient states. At last, computational cost is summarized in tables with focus on the computing time and the number of degree of freedoms. All examples are two dimensional with the standard set of parameters given in Table 1 (note the very low shear modulus of 15 kPa [17]). The simulation time span covers the whole swelling process (0 s to \(10^{6}\) s) and is divided into 48 adaptive steps with smaller time step size at the initial stage and much larger step size near equilibrium. All the simulations are done using MATLAB 2016b on a windows machine with Intel Core 17-6700 CPU and 64GB RAM.

Fig. 3
figure 3

Chemical potential contour plot after one Newton iteration with mesh \(20\times 20\), \(\hbox {G}=0.015\) MPa

Fig. 4
figure 4

Situations where the wavy edges disappears for FEM solutions

Fig. 5
figure 5

Chemical potential contour plot after one time step (0.1 s) with mesh \(10\times 10\), \(\hbox {G}=0.15\) MPa

Table 4 \(\hbox {FEM}(8\times 4\mu )\) solution convergence performance summary

5.1 Swelling of a square-shaped gel with constant permeability

In the first example, the geometry is a quarter square as given in Fig. 2. To avoid rigid body motions and to account for symmetry constraints, node O is fixed in both x and y directions and edge OA and OC are fixed in y and x direction respectively (\(\varGamma _D^u\)). Edge BC and BA are traction free (\(\varGamma _N^u\)) and in touch with outer solution (\(\varGamma _D^\mu \)). Fluid enters the sample as the arrows indicate. Edge OA and OC are equipped with no flow boundary conditions (\(\varGamma _N^\mu \)) because of symmetry constraints. To summarize, the initial and boundary conditions for free swelling simulation are given as:

$$\begin{aligned} {\varvec{x}}({\varvec{X}},0)&={\varvec{X}}, \end{aligned}$$
(44)
$$\begin{aligned} \mu ({\varvec{X}},0)&= -\pi _0, \end{aligned}$$
(45)
$$\begin{aligned} {\varvec{x}}({\varvec{X}},t)&={\varvec{0}},\;\;\; {\varvec{X}}\in \varGamma _D^u, \end{aligned}$$
(46)
$$\begin{aligned} {{\varvec{T}}}({\varvec{X}},t) \cdot {\varvec{n}}_R&={\varvec{0}},\;\;\; {\varvec{X}}\in \varGamma _N^u, \end{aligned}$$
(47)
$$\begin{aligned} \mu ({\varvec{X}},t)&=\mu _D,\;\;\; {\varvec{X}}\in \varGamma _D^\mu , \end{aligned}$$
(48)
$$\begin{aligned} {\varvec{Q}}({\varvec{X}},t)\cdot {\varvec{n}}_R&=0,\;\;\; {\varvec{X}}\in \varGamma _N^\mu , \end{aligned}$$
(49)

where \(\pi _0=RT\sqrt{(c^{fc}_0)^2+4{\bar{c}}^2}\) and \(\mu _D=-2RT{\bar{c}}\).

Fig. 6
figure 6

Domain area calculation error of FEM and MHFEM solution with different mesh sizes at two different moments of swelling

First of all, we investigate the solution convergence behavior of each methods. Considering three (key) simulation parameters: mesh size, time step size and shear modulus, the solution convergence results of FEM and MHFEM are summarized in Tables 2 and 3 respectively. Note that \(\varDelta t\) indicates the standard time step we use throughout this study with the initial time step 0.1 s. From Table 2, we see that as the shear modulus decreases (larger degree of swelling expected), FEM shows severe convergence issues. As a matter of fact, five mesh sizes and time steps combinations fail to converge for shear modulus 0.055 MPa and none of the mesh sizes and time steps lead to convergence when shear modulus is further decreased to the desirable value 0.015 MPa.

In contrary to the convergence issues suffered by FEM, MHFEM are much robust in reaching solution convergence. Table 3 shows that the only case that MHFEM fails to converge is at shear modulus 0.015 with the mesh \(65\times 65\) and time step \(\varDelta t\). A further refinement of time steps (dividing \(\varDelta t\) by 10) leads to the convergence again.

Fig. 7
figure 7

Contour plots of chemical potential and streamlines at time \(t=5.32\) s

Fig. 8
figure 8

Contour plots of chemical potential and streamlines at time \(t=2740.1\) s

Table 5 Standard FEM computing performance summary
Table 6 MHFEM computing performance summary

In order to better ideas about the performance difference between the two methods, we examine the chemical potential contour plot after the very first Newton iteration of FEM and MHFEM. Using a relatively coarse mesh (\(30\times 30\)) with shear modulus 0.015 MPa, FEM (Fig. 3b) solution shows certain wave-like instability on the edge. Such instability later may stop propagate (for example in some cases with shear modulus 0.055 MPa) and eventually lead to convergence. However, with shear modulus 0.015 MPa the instability persists and within several Newton iterations excessive distortions in elements appear. The wavy edges tend to disappear when the mesh is refined or the shear modulus is increased as shown in the Fig. 4. In the case of large shear modulus (G = 0.15 MPa), the eventual convergence is indeed achieved. In contrast, MHFEM generates no wavy edges at all after the first Newton iteration and indeed yields convergence for a broader combinations parameters.

Bouklas et al. [7] argue that at the beginning of the swelling, due to the incomprehensibility of solid, Taylor–Hood elements should be applied to satisfy the LBB condition between the displacement and chemical potential field. It was shown that less oscillation appears in chemical potential field when the displacement field is interpolated with one order higher function spaces than that of the chemical potential field. In order to preclude the influence of such oscillations on the convergence results, we carried out FEM simulations using quadratic elements for displacement and linear elements for the chemical potential (\(8\times 4\mu \)). As shown in Fig. 5, the oscillation in chemical potential shown in standard FEM (\(4\times 4\mu \)) is indeed disappeared when the higher order interpolation functions for displacement is used. However, the less oscillation does not improve the general solution convergence performance. As a matter of fact, the standard FEM (\(8\times 4\mu \)) delivers the worst convergence robustness results among the three as indicated in Table 4.

Fig. 9
figure 9

Fitted curves between total time and DOF for each method

To summarize, the appearance of wavy curves in FEM acts like a source of instability and causes the failure to converge. Refinement in mesh size and increase in shear modulus help to diminish the wavy curves. Satisfying LBB condition for displacement filed and chemical potential field in FEM simulations reduces the oscillation of the chemical potential field but does not improve its solution convergence performance. In swelling simulations that involves extremely large deformations, MHFEM shows more robust performance than standard FEM. This might be related to the local mass conservation property which is possessed by MHFEM only.

For accuracy analysis, the shear modulus value is taken to be 0.08 MPa to guarantee the solution convergence for both methods. To compare quantitatively the deformations calculated by the two simulations, the total area of the domain (sum of all elements areas) is calculated with different mesh sizes. Mesh sizes h are taken to be 0.1 mm, 0.06 mm, 0.047 mm, 0.033 mm, 0.025 mm and 0.02 mm. Since there is no analytical solution available in two dimensions for a nonlinear transient swelling problem, we take the averaged solution between MHFEM and FEM obtained with the finest mesh (\(50\times 50\) linear quadrilateral elements of 0.02 mm each) as the canonical solution. In Fig. 6, we plot the dimensionless error e (absolute value of area canonical minus area FEM/MHFEM) and mesh size h in a double log figure at two transient swelling moments. We observe that both methods converge to the same results as the mesh is refined. At the early stage (\(16.16\%\) swollen) of welling, FEM (blue continuous line) gives smaller errors than MHFEM when the mesh is coarse. The difference between the error gets smaller at the later stage (\(78.15\%\) swollen). The convergence rate for MHFEM and standard FEM are comparable.

In order to gain insight into the chemical potential distribution and net flux in each element calculated by each method, we plot the chemical potential contour with streamlines (stemming from the edge of the sample) on top for both MHFEM and standard FEM solution at two selected moments (Figs. 7, 8). The chemical potential and streamlines distribution derived by the two methods are almost identical for both moments. Due to the symmetry in boundary conditions and geometry, both chemical potential and streamlines exhibit symmetry for both methods.

The computing cost of both methods is summarized in Tables 5 and 6. For different mesh sizes, total computing time, total numbers of iterations, time per iteration to build the tangent stiffness matrix and residual vector and degree of freedoms (DOF) are recorded. In addition, we fitted the curve between the total time and degree of freedoms using second order polynomials for each method (Fig. 9). We conclude that for the same mesh FEM requires less degree of freedoms but more iterations to reach equilibrium. MHFEM is more efficient in terms of computing time per degree of freedom.

Table 7 FEM solution convergence performance summary
Table 8 MHFEM solution convergence performance summary

5.2 Swelling of a square-shaped gel with strain-dependent permeability

During transient swelling, permeability k is typically not constant over the domain. The part of the gel that is in touch with the outer solution earlier has higher hydraulic permeability due to the enlarged pore. This observation has been reported by several studies [3, 12, 23]. Depending on swelling materials, different permeability-strain dependencies are observed and corresponding mathematical relationships are proposed [10, 20, 28, 30, 41, 42]. In this example, we adopt the law from [18] which is developed for polyacrylamide gels:

$$\begin{aligned} k(\phi _f)=k_0\frac{(1-\phi {f,0})^\beta }{\phi _{f,0}}\frac{\phi _f}{(1-\phi _f)^\beta }, \end{aligned}$$
(50)

where \(\beta =1.5\). Note that the current porosity \(\phi _f\) is related to strain by the relation

$$\begin{aligned} \phi _f=1-\frac{\phi _{s,0}}{J}, \end{aligned}$$
(51)

on element level.

Like in the first example, we carry out solution convergence investigations. The results is summarized in Tables 7 and 8. We notice that the convergence results are the same as in the previous example. Therefore, similar conclusions can be drawn and will not be repeated here.

Next, we investigate the accuracy of the each method in terms of mesh convergence rate. Again, we take shear modulus to be 0.08 MPa to guarantee solution convergence for both methods. In Fig. 10, errors and the mesh sizes at two different moments are plotted in a double-log axis. Standard FEM shows accuracy advantage at both moments with the advantage at an earlier moment (\(24.73\%\) swollen) is slightly larger. Both methods converge to the same solution at a comparable rate. The chemical potential and streamlines distributions over the sample at two given moments are plotted in Figs. 11 and 12. Similar conclusions can be drawn as in the previous example.

Fig. 10
figure 10

Domain area calculation error of FEM and MHFEM solution with different mesh sizes at two different moments of swelling

Fig. 11
figure 11

Contour plots of chemical potential and streamlines at time \(t=5.32\) s

Fig. 12
figure 12

Contour plots of chemical potential and streamlines at time \(t=2740.1\) s

At last, the computing cost summary is given in Tables 9 and 10. Comparing to the constant permeability case, the total computing time increases because of higher numbers of iterations. FEM once again demands less computing time compared to MHFEM due to the smaller DOF; however, MHFEM needs less iterations in total to reach equilibrium.

Table 9 Standard FEM computing performance summary
Table 10 MHFEM computing performance summary

5.3 Swelling of a gel with a low-permeability stripe

In this example, we deal with heterogeneous permeability with discontinuous jumps. The initial and typical transient swelling geometry is given in Fig. 13. The yellow zone represents high permeability area with permeability \(p_1\) and the blue zone low permeability area with permeability \(p_2\). The sample is confined in x direction thus swelling is only allowed in y direction. Only edge AB is in touch with the outer solution. Translated into initial and boundary conditions applied in this simulation, they are:

$$\begin{aligned} {\varvec{x}}({\varvec{X}},0)&={\varvec{X}}, \end{aligned}$$
(52)
$$\begin{aligned} \mu ({\varvec{X}},0)&= -\pi _0, \end{aligned}$$
(53)
$$\begin{aligned} {\varvec{x}}({\varvec{X}},t)&={\varvec{0}},\;\;\; {\varvec{X}}\in \varGamma _D^u, \end{aligned}$$
(54)
$$\begin{aligned} {{\varvec{T}}}({\varvec{X}},t) \cdot {\varvec{n}}_R&={\varvec{0}},\;\;\; {\varvec{X}}\in \varGamma _N^u, \end{aligned}$$
(55)
$$\begin{aligned} \mu ({\varvec{X}},t)&=\mu _D,\;\;\; {\varvec{X}}\in \varGamma _D^\mu , \end{aligned}$$
(56)
$$\begin{aligned} {\varvec{Q}}({\varvec{X}},t)\cdot {\varvec{n}}_R&=0,\;\;\; {\varvec{X}}\in \varGamma _N^\mu , \end{aligned}$$
(57)

where \(\varGamma _D^u\) contains edge AB in both x and y direction plus edge CA and edge BD in x direction only; \(\varGamma _N^u\) contains only edge CD; \( \varGamma _D^\mu \) contains only edge AB; and \(\varGamma _N^\mu \) contains edges AC, CD and BD. Note that \(\pi _0=RT\sqrt{(c^{fc}_0)^2+4{\bar{c}}^2}\) and \(\mu _D=-2RT{\bar{c}}\). This is the same as in the previous simulations.

Due to the presence of the low permeability stripe, top surface CD forms typically a curve with the middle point lowest in y direction during transient swelling. In this case, the deformations of the gel can be completely characterized by curve CD.

Fig. 13
figure 13

Initial and typical transient swelling geometry of the sample. Yellow(light): high permeability area, Blue(dark): low permeability area

Fig. 14
figure 14

Convergence rate of FEM and MHFEM solutions with different shear modulus at \(\hbox {t}=120~\hbox {s}\)

Fig. 15
figure 15

Convergence rate of FEM and MHFEM solutions with different permeability contrast at \(\hbox {t}=120~\hbox {s}\)

Following the mesh convergence analysis in the previous two examples, we investigate the changes of errors of each method as the mesh sizes decrease. In this example, mesh sizes h are chosen to be 1/12, 1/24, 1/36/, 1/48, 1/60, 1/120 mm. The averaged solution between FEM and MHFEM at \(h=1/120\) mm is used as the canonical solution. In this example, the time step size is fixed to be 10 s and we run 12 steps (thus \(t=120\) s). The standard value for \(p_1\) and \(p_2\) are \(10^{-2}\) mm\(^{4}\)/(Ns) and \(10^{-6}\) mm\(^{4}\)/(Ns).

Fig. 16
figure 16

Contour plots of chemical potential at time \(t=120\) s with mesh size 0.0417 mm

Fig. 17
figure 17

Contour plots of chemical potential at time \(t=120\)  s with mesh size 0.0417 mm shear modulus 0.15 MPa

Fig. 18
figure 18

Contour plots of chemical potential at time \(t=120\) s with mesh size 0.0417 mm small permeability contrast

Fig. 19
figure 19

Contour plots of chemical potential at time \(t=120\) s with fine mesh (0.0083 mm)

Table 11 standard FEM computing performance summary
Table 12 MHFEM computing performance summary

In Fig. 14, errors of both methods with different shear moduli (thus different swelling degrees) are plotted against mesh sizes. We see that a lower shear modulus (implies larger swelling degree) generates larger errors than higher shear modulus for both methods. MHFEM shows smaller errors than standard FEM with shear moduli 0.03 MPa and 0.08 MPa for all mesh sizes. As to shear modulus 0.15 MPa, MHFEM shows smaller error when the mesh is coarse; the error performance is very comparable as the mesh is refined.

Figure 15 shows convergence results under different high/low permeability contrast ratios (100, 1000 and 10,000). A larger contrast ratio implies larger errors for both methods. We observe that MHFEM generates smaller errors than standard FEM in all three cases. Moreover, the error differences between standard FEM and MHEFEM are not affected by the contrast ratio. In other words, a smaller permeability contrast ratio does not necessarily suggest a more comparable accuracy performance between FEM and MHFEM.

The advantages of MHFEM in accuracy over standard FEM in deformation field are extended to the chemical potential and net flux. Figure 16 compares MHFEM and FEM in terms of chemical potential contour with streamlines (stemming from the bottom side) with a relatively coarse mesh (mesh size about 0.0417 mm) at time \(t=120\) s. The top surface of MHFEM solution appears smooth, whereas nonphysical distortion of elements near the low permeability stripe in FEM solution is observed. A closer look at the streamlines near the low permeable area, one can find that symmetry of streamlines are lost in the FEM solution despite of symmetric geometry and boundary condition; on the other hand, the MHFEM solution produces perfectly symmetric streamlines.

After increasing the shear modulus (from 0.015 to 0.15 MPa), we observe that the degree of swelling is decreased (Fig. 17). At the same time, the distortion of the local elements as well as the asymmetry of the streamlines of the FEM solutions are alleviated.

Next we decrease the permeability contrasts between high and low permeable areas. Specifically, the low-permeable area is increased from \(10^{-6}\) to \(10^{-5}\) mm\(^4\)/(Ns) and the high-permeable area is decreased from \(10^{-2}\) to \(10^{-3}\) mm\(^4\)/(Ns). Figure 18 shows that the degree of swelling is further reduced because of the decrease of permeability in the high- permeable area. In this case, top surface deformation, the chemical potential and streamlines distribution of FEM solution resemble very much with the ones of MHFEM solution.

At last, we investigate the influence of mesh size on the solution. By refining the mesh, the chemical potential and streamline plot as well as the top surface deformation in the FEM solution are much close to the MHFEM solution (Fig. 19) compared to Fig. 16. Although the complete symmetry in streamlines in the FEM solution is still missing, mesh refinement appears to have improved the situation greatly.

The solution convergence behavior in this case is consistent with the first two examples. Besides the conclusions we have drawn in Example 1, we notice that for coarse mesh in the FEM solution excessive element distortion tends to appear near the discontinuity interface (e.g. Fig. 16b) and that is another source contributing to convergence failure for FEM.

Finally, the computing cost summary is given in Tables 11 and 12. Standard FEM shows advantage (especially) in total computing time for a dense mesh due to lower degrees of freedom. The number of iterations, unlike in the previous two examples, is almost identical for both FEM and MHFEM and does not decrease as the mesh is refined. In addition, the fitted curves (Fig. 20) are very close to each other. We believe that the similar numbers of iterations between FEM and MHFEM is due to the fact that the swelling simulation in this example is of one-dimensional nature (only one edge is in touch of outer solution, namely edge AB). Therefore, to reach equilibrium smaller volume change is involved and thus the advantage of MHFEM in solving the equations more efficiently is less evident in this example.

Fig. 20
figure 20

Fitted curves between total time and DOF for each method

Fig. 21
figure 21

Schematic illustration of the core-shell structure

Fig. 22
figure 22

MHFEM solution of crack opening angle \(\theta \) with different mesh sizes

Fig. 23
figure 23

FEM solution of crack opening angle \(\theta \) with different mesh sizes

Fig. 24
figure 24

MHFEM solution of crack opening angle \(\theta \) zoomed in at the initial stage

Fig. 25
figure 25

FEM solution of crack opening angle \(\theta \) zoomed in at the initial stage

Fig. 26
figure 26

MHFEM solution of shell depth d with different mesh sizes

Fig. 27
figure 27

FEM solution of shell depth d with different mesh sizes

5.4 Swelling of gels with core-shell structure

The last numerical example features an industrially-relevant gel structure (see Fig. 21). To prevent the swelling SAPs from blocking the passage of inter-particles fluid, nowadays the outer layer (shell) of a SAP is more heavily crosslinked than the inner core to achieve maximum absorbing efficiency [38]. As a result, the outer layer is stiffer and with lower permeability than the inner core. As the swelling continues, the stiffer shell fractures and develops into isolated islands. This example aims to simulate the scenario after the isolated islands are formed and the swelling dynamics from then on.

In the simulations, the shear modulus of the shell is assumed to be 0.15 MPa (thus 10 times larger than that of the core). The permeability of the shell is \(10^{-7}\) mm\(^{4}\)/(Ns) (4 orders of magnitude smaller than of the core). There are crack openings on the core between the isolated shell islands present. The crack does not propagate in the simulation. Further, to avoid rigid body motions, we fix O in both x and y directions and node A in y direction. Curve AB (including the outer contour of the shell) is in touch with the other solution. Edge OA and OB are with no flow boundary conditions. The initial and boundary conditions can be summarized as:

$$\begin{aligned} {\varvec{x}}({\varvec{X}},0)&={\varvec{X}}, \end{aligned}$$
(58)
$$\begin{aligned} \mu ({\varvec{X}},0)&= -\pi _0, \end{aligned}$$
(59)
$$\begin{aligned} {\varvec{x}}({\varvec{X}},t)&={\varvec{0}},\;\;\; {\varvec{X}}\in \varGamma _D^u, \end{aligned}$$
(60)
$$\begin{aligned} {{\varvec{T}}}({\varvec{X}},t)\cdot {\varvec{n}}_R&={\varvec{0}},\;\;\; {\varvec{X}}\in \varGamma _N^u, \end{aligned}$$
(61)
$$\begin{aligned} \mu ({\varvec{X}},t)&=\mu _D,\;\;\; {\varvec{X}}\in \varGamma _D^\mu , \end{aligned}$$
(62)
$$\begin{aligned} {\varvec{Q}}({\varvec{X}},t)\cdot {\varvec{n}}_R&=0,\;\;\; {\varvec{X}}\in \varGamma _N^\mu , \end{aligned}$$
(63)
Fig. 28
figure 28

trajectory of node B in x direction with different mesh sizes

Fig. 29
figure 29

trajectory of node B in y direction with different mesh sizes

Fig. 30
figure 30

Contour plots of chemical potential at time \(t=0.1\) s with mesh size 0.013 mm, also see movies crack_MHFEM_0.013 and crack_FEM_0.013

Fig. 31
figure 31

Contour plots of chemical potential at time \(t=76.8\) s with mesh size 0.013 mm, also see movies crack_MHFEM_0.013 and crack_FEM_0.013

where \(\varGamma _D^u\) contains edge OA in y direction plus origin O in both x and y directions; \(\varGamma _N^u\) contains only edge OB and the irregular surface between A and B (including the shell surface plus the crack opening surface); \( \varGamma _D^\mu \) contains the irregular surface between A and B (including the shell surface plus the crack opening surface); and \(\varGamma _N^\mu \) contains edges OA, and OB. Note that \(\pi _0=RT\sqrt{(c^{fc}_0)^2+4{\bar{c}}^2}\) and \(\mu _D=-2RT{\bar{c}}\). This is the same as in the previous simulations.

Fig. 32
figure 32

Contour plots of chemical potential at time \(t=0.1\) s with mesh size 0.0052 mm, also see movies crack_MHFEM_0.013 and crack_FEM_0.013

Fig. 33
figure 33

Contour plots of chemical potential at time \(t=76.8\) s with mesh size 0.0052 mm, also see movies crack_MHFEM_0.013 and crack_FEM_0.013

The displacement boundary in this numerical example is not symmetric with respect to edge OA and OB. This was chosen on purpose by the authors considering the fact that by relaxing edge OB’s movement in x-direction, the difference in swelling dynamics between FEM simulation and MHFEM simulation can be better highlighted. This is especially true for the crack opening angles. By allowing the x-displacement of edge OB, the domain loses its symmetry during transient swelling and the crack angle (\(\theta \)) will be opened up as the results of this x-displacement. In this example, mesh sizes are chosen from small to large: 0.013 mm, 0.0087 mm, 0.0065 mm and 0.0052 mm. Note that the simulation results are also available in video form. They are included in the electronic supplement materials.

To characterize the swelling dynamics of such a structure, we examine the time evolution of crack opening angle \(\theta \) and shell depth d, trajectory of node B (in both x and y directions), the chemical potential and flux chart over the sample with different meshes sizes.

Following previous examples, the canonical solution is derived by averaging FEM and MHFEM solution with the smallest mesh size (0.0052 mm). In Figs. 22 and 23, we plot the MHFEM and standard FEM solution of crack opening angle \(\theta \) against time with different mesh sizes respectively. We observe that MHFEM is able to generate reasonably accurate solution even with the most coarse mesh (0.013 mm); while FEM shows large deviation from the canonical solution with the most coarse mesh (0.013 mm). In Figs. 24 and 25, we zoomed into the early stage of swelling to highlight the difference between the FEM and MHFEM simulation. Using a denser mesh alleviates the problem but the disagreement with the canonical solution stays relatively large compared to the MHFEM solutions. As a matter of fact, only with the most dense mesh (0.0052 mm) FEM solution is close enough to the canonical solution (less than 0.5 degree deviation).

Next, we plot the MHFEM and FEM solution of time evolution of shell depth d in Figs. 26 and 27. An interesting phenomenon showed by both simulations is that the depth of the shell first decrease and then increase, which suggests that the shell first shrinks and then swells. This can be explained by the large permeability difference between the shell and the core. Since the core is 10,000 times more permeable than the shell, any change of chemical potential will lead to a much more significant amount flow. As a result, as soon as the swelling starts, the slight change in chemical potential in the shell leads to the flow into the core elements underneath the shell. Even if at the same time there is external fluid flowing into the shell itself, since the core is much more permeable than the shell, the flow that goes out into the core beneath the shell is much larger than the flow goes into the shell. Therefore, only after the core is nearly fully swollen the shell starts to swell too. Similar to the observations made in the crack opening investigation, FEM solution with the most coarse mesh (0.013 mm) fails to capture the this dynamic of the depth of the shell. Further refinement of the mesh leads to convergence to the canonical solution but the deviation from the canonical solution of FEM is larger than that of MHFEM solution.

Besides crack opening angle \(\theta \) and shell depth d, the trajectory of node B (which is free to move in both x and y directions) as an indication of the soft core geometry evolution is studied. Figure 28 plots the MHFEM and FEM solution of the trajectory of node B in x direction and Fig. 29 in y direction. They show that the disagreement with the canonical solution happens mostly in x direction with FEM with the most coarse mesh (0.013 mm) yielding the largest deviation. The calculation in y direction seems more unified compared with x direction regardless of mesh sizes and numerical methods. Overall, the accuracy advantages of MHFEM over FEM are apparent in the calculation of trajectory B.

The difference in the deformations calculation is caused by the calculation of the net flux in/out of each element. In what follows four chemical potential contour plots featuring comparison between MHFEM and FEM solutions with two mesh sizes (dense and course) at two selected moments (right after swelling starts and a short while after) are presented. The black arrows in each element indicates the direction of the net flux in the given element and its length indicates the magnitude of the flux.

For a coarse mesh (Fig. 30), the FEM solution generates a visibly deformed (swollen) shell with the highest chemical potential at the edge of the shell; while there is no visible deformation of the shell observed in MHFEM solution and the chemical potential stays relatively uniform except from the part of the soft core directly in touch with the outer solution. Given a closer look at the net flux over the sample, one can find that the net flux given by MHFEM solution is more physically plausible, since the largest influx comes from the edge of the soft core that is in direct contact with the outer solution and barely no flux exchange at the core-shell interface due to the low permeability of the shell. On the other hand, due to the absence of local mass conservation, FEM shows a visible influx from the core to the shell right at the beginning of swelling which is clearly not physical.

The earlier discrepancy in chemical potential and flux calculation becomes more prominent as the swelling goes on. At time \(t=76.8\) s, the chemical potential at the shell stays the lowest over the sample for MHFEM solution as the consequence of little influx till this moment; however, the shell in the FEM solution is fully swollen and with the highest chemical potential over the sample (Fig. 31). The flux distribution is also different between the two solutions especially in the element right beneath the shell. The difference in flux distribution clearly has contributed to the aforementioned difference in crack opening angle, depth of the shell and node B trajectories.

By deploying a denser mesh, FEM simulations results are greatly improved (Figs. 32, 33) and resemble MHFEM solutions at both moments. However, FEM still suffers slightly from local nonphysical distortion of elements at the edge of shell as shown in Fig. 32b.

6 Concluding remarks

In this work, we first reviewed the Donnan osmosis based swelling mechanism. Then two computation formulations was introduced in detail. The standard FE formulation uses deformation \({\varvec{x}}\) and chemical potential \(\mu \) as prime variables. The flux \({\varvec{Q}}\) is calculated by numerical differentiation. On the other hand, MHFEM adopts mixed formulation in which \({\varvec{Q}}\) and \(\mu \) are considered as independent variables. Lowest order Raviart–Thoams element was used in the mixed formulation, which guarantees the conservation of mass locally. Such a property has proven to be important in the calculation of flux field solving Darcy’s type equations with heterogeneous permeability [26, 36]. By examining four numerical examples, we evaluate the importance of using MHFEM in swelling simulations involving extremely large deformations.

The four numerical examples falls into three categories (depending on the permeability distributions): constant permeability, strain-dependent permeability and permeability with discontinuous jumps. For each example, aspects such as robustness in solution convergence, accuracy (mesh convergence rate) and computational cost are investigated.

The first two examples involve the swelling of a quarter of a square-shaped gel with constant and strain-dependent permeability respectively. FEM demands much more stringent condition on mesh size and time step size than MHFEM to reach convergence especially when the swelling degree increases. Providing solution convergence, FEM shows certain advantage in the accuracy of deformation calculation than MHFEM; on the other hand, chemical potential and flux field calculations are quite similar between the two. The third and fourth example are about permeability with discontinuous jumps. Example 3 deals with swelling confined in y direction with a low-permeable stripe. By measuring the deformation of the top surface, we conclude that MHFEM generates more accurate results than FEM especially when the mesh is coarse. The chemical potential and flux calculation also supports a similar conclusion since the symmetry of flux streamlines was lost in FEM but not in MHFEM. Moreover, FEM generates nonphysical distortion of elements near the discontinuous jump with a coarse mesh which contributes to a more unfavorable position when the robustness of solution convergence is considered. Example 4 studies the two dimensional swelling of a core-shell structure that is widely used in disposable diaper industry. The soft high-permeable core is surrounded by hard low-permeable shell islands. Simulations of such a structure have shown that FEM using a coarse mesh can lead to a severely wrong prediction of the swelling dynamics; while MHFEM is able to generate reasonably accurate solution with the same (coarse) mesh. Further investigations on chemical potential and flux also favor MHFEM in terms of accuracy. The computational cost favors FEM slightly due to less degrees of freedom. However, if one takes the accuracy argument into account, MHFEM is not necessarily in disadvantage when the computational costs are calculated to achieve certain accuracy.

To summarize, the local mass conservation property of MHFEM has been proven to improve the accuracy of the solution in swelling simulations when there is discontinuities in permeability present. Moreover, MHFEM largely outperforms FEM in terms of solution convergence robustness.