Analysis of Fully Discrete Finite Element Methods for a System of Differential Equations Modeling Swelling Dynamics of Polymer Gels

# Analysis of Fully Discrete Finite Element Methods for a System of Differential Equations Modeling Swelling Dynamics of Polymer Gels

## Abstract

The primary goal of this paper is to develop and analyze some fully discrete finite element methods for a displacement-pressure model modeling swelling dynamics of polymer gels under mechanical constraints. In the model, the swelling dynamics is governed by the solvent permeation and the elastic interaction; the permeation is described by a pressure equation for the solvent, and the elastic interaction is described by displacement equations for the solid network of the gel. The elasticity is of long range nature and gives effects for the solvent diffusion. It is the fluid-solid interaction in the gel network drives the system and makes the problem interesting and difficult. By introducing an “elastic pressure” (or “volume change function”) we first present a reformulation of the original model, we then propose a time-stepping scheme which decouples the PDE system at each time step into two sub-problems, one of which is a generalized Stokes problem for the displacement vector field (of the solid network of the gel) and another is a diffusion problem for a “pseudo-pressure” field (of the solvent of the gel). To make such a multiphysical approach feasible, it is vital to find admissible constraints to resolve the uniqueness issue for the generalized Stokes problem and to construct a “good” boundary condition for the diffusion equation so that it also becomes uniquely solvable. The key to the first difficulty is to discover certain conservation laws (or conserved quantities) for the PDE solution of the original model, and the solution to the second difficulty is to use the generalized Stokes problem to generate a boundary condition for the diffusion problem. This then lays down the theoretical foundation for one to utilize any convergent Stokes solver (and its code) together with any convergent diffusion equation solver (and its code) to solve the polymer gel model. In the paper, the Taylor-Hood mixed finite element method combined with the continuous linear finite element method are chosen as an example to present the ideas and to demonstrate the viability of the proposed multiphysical approach. It is proved that, under a mesh constraint, both the proposed semi-discrete (in space) and fully discrete methods enjoy some discrete energy laws which mimic the differential energy law satisfied by the PDE solution. Optimal order error estimates in various norms are established for the numerical solutions of both the semi-discrete and fully discrete methods. Numerical experiments are also presented to show the efficiency of the proposed approach and methods.

G

els, soft matters, poroelasticity, Stokes equations, finite element methods, inf-sup condition, fully discrete schemes, error estimates.

{AMS}

65M12, 65M15, 65M60,

## 1 Introduction

A gel is a soft poroelastic material which consists of a solid network and a colloidal solvent. The solid network spans the volume of the solvent medium. The solvent can permeate through the solid network and the permeation can be controlled by external forces. Both by weight and volume, gels are mostly liquid in composition and thus exhibit densities similar to liquids. However, they have the structural coherence of a solid and can be deformed. A gel network can be composed of a wide variety of materials, including particles, polymers and proteins, which then gives different types gels such hydrogels, organogels and xerogels (cf. [9, 13]). Gels have some fascinating properties, in particular, they display thixotropy which means that they become fluid when agitated, but resolidify when resting. In general, gels are apparently solid, jelly-like materials, they exhibit an important state of matter found in a wide variety of biomedical and chemical systems (cf. [9, 10, 19, 20] and the references therein).

This paper develops and analyzes some fully discrete finite element methods for a displacement-pressure model for polymer gels. The model, which was proposed by M. Doi et al in [9, 19, 20], describes swelling dynamics of polymer gels (under mechanical constraints). Let be a bounded domain and denote the initial region occupied by the gel. Let denote the displacement of the gel at the point in the space and at the time , and be the velocity and the pressure of the solvent at . Following [9], the governing equation for the swelling dynamics of polymer gels are given by

 (1) \rm div\,(σ(u)−pI) =0, (2) ξ(vs−ut) =−(1−ϕ)∇p, (3) \rm div\,(ϕut+(1−ϕ)vs) =0.

Here is the friction constant associated with the motion of the polymer relative to the solvent, is the volume fraction of the polymer, denotes the identity matrix, and stands for the stress tensor of the gel network, which is given by a constitutive equation. In this paper, we use the following linearized form of the stress tensor:

 (4) σ(u):=(K−23G)\rm div\,uI+2Gε(u),ε(u):=12(∇u+∇uT),

where and are respectively the bulk and shear modulus of the gel (cf. [9, 7]). We remark that (1) stands for the force balance, (2) states Darcy’s law for the permeation of solvent through the gel network, and (3) describes the incompressibility condition. In addition, if we introduce the total stress , then equation (1) becomes .

Substituting (4) into (1) and (2) into (3) yield the following basic equations for swelling dynamics of polymer gels (see [19])

 (5) α∇\rm div\,u+βΔu =∇p, α:=K+G3,β:=G, (6) \rm div\,ut =κΔp, κ:=(1−ϕ)2ξ,

which hold in the space-time domain for some given .

To close the above system, we need to prescribe boundary and initial conditions. Only one initial condition is required for the system, which is

 (7) u(⋅,0)=u0(⋅)in Ω.

Various sets of boundary conditions are possible and each of them describes a certain type mechanical condition and solvent permeation condition (cf. [19, 20]). In this paper we consider the following set of boundary conditions

 (8) (σ(u)−pI)ν=f,∂p∂ν=0on ΩT:=∂Ω×(0,T),

where denotes the outward normal to . (8) means that the mechanical force is applied on the boundary of the gel. Since , hence, (8) implies that the solvent can not permeate through the gel boundary. We also remark that the force function must satisfy the compatibility condition

 ∫∂ΩfdS=0.

Problem (5)–(8) is interesting and difficult due to its multiphysical nature which describes the complicate fluid and solid interaction inside the gel network. It is numerically tricky to solve because it is difficult to design a good and workable time-stepping scheme. For example, one natural attempt would be at each time step first to solve a Poisson problem for and then to solve a linear elasticity problem for . However, this strategy is difficult to realize because there is no good way to compute the source term for the Poisson equation. In fact, the strategy even has a difficulty to start due to the fact that no initial condition is provided for the pressure .

To overcome the difficulty, in this paper we shall use a reformulation of system (5)–(6), which is now introduced. Define

 (9) q:=\rm div\,u.

Physically, measures the volume change of the solid network of the gel, and often called “elastic pressure” or “volume change function”. Taking divergence on (5) yields

 (α+β)Δq=Δp,

which and (6) imply that satisfies the following diffusion equation

 (10) qt=DΔq,D:=κ(K+43G).

However, the usefulness of the above diffusion equation is hampered by the lack of boundary condition for . We like to note that the above diffusion equation for was first noticed by M. Doi [8], but it was not utilized before exactly because of the lack of boundary condition for . Nevertheless, using the new variable we can rewrite (5)–(6) as

 (11) βΔu =∇˜p, ˜p:=p−αq, (12) \rm div\,u =q, (13) qt =κΔp, p=˜p+αq.

An immediate consequence of the above reformulation is that (11)-(12) implies satisfies the generalized Stokes equations with being the source term at each time , and satisfies a diffusion equation and it interacts with only at the boundary .

This is a key observation because it not only reveals the underlying physical process of swelling dynamics of the gel, but also gives the “right” hint on how the problem should be solved numerically. This indeed motivates the main idea of this paper, that is, at each time step, we first solve the generalized Stokes problem for , which in turn provides (implicitly) an updated boundary condition for , we then use this new boundary condition to solve the diffusion equation for . The process is repeated iteratively until the final time step is reached. However, in order to make this idea work, there is one crucial issue needs to be addressed. That is, for a given the generalized Stokes problem for is only unique up to additive constants. Clearly, how to correctly enforce the uniqueness of the generalized Stokes problem is the bottleneck of this approach. It is easy to understand that one can not use arbitrary constraints to fix because this will lead to bad or even divergent numerical schemes if the exact PDE solution does not satisfy the constraints. Instead, the constraints which can be used to fix should be those satisfied by the exact solution of the PDE system. To the end, we need to discover some invariant (or conserved) quantities for the exact PDE solution. It turns out that the situation is precisely what we anticipated and wanted. We are able to show that the exact PDE solution satisfies the following identities (see Section 2 below for a proof):

 (14) ∫Ωq(x,t)dx ≡Cq:=∫Ωq0(x)dx:=∫Ω\rm div% \,u0(x)dx, (15) ∫Ωp(x,t)dx ≡Cp:=cdCq−1d∫∂Ωf(x,t)⋅xdS, (16) ∫∂Ωu(x,t)⋅νdS ≡Cu:=∫∂Ωu0(x)⋅νdS=Cq,

where denote the dimension of and

 (17) cd:=α+βd={K+5G6if d=2,K+2G3if d=3.

Obviously, the right-hand sides of (14) and (16) are constants. The right-hand side of (15) is also a constant provided that is independent of , otherwise, it is a known function of . In this paper, we shall only consider the case that is independent of . It follows from (14) and (15) that

 (18) ∫Ω˜p(x,t)dx≡C˜p:=Cp−αCq=βCqd−1d∫∂Ωf(x)⋅xdS.

It is clear now that (18) and (16) provide two natural conditions which can be used to uniquely determine the solution to the generalized Stokes problem (11)–(12) for a given source term . This then leads to the following time-discretization for problem (5)–(8):

Algorithm 1:

• Set and .

• For , do the following two steps

Step 1: Solve for such that

 (19) −βΔun+1+∇˜pn+1 =0, in ΩT, (20) \rm div\,un+1 =qn in ΩT, (21) β∂un+1∂ν−˜pn+1ν =f on ∂ΩT, (22) (˜pn+1,1)=C˜p,⟨un+1,ν⟩ =Cu.

Step 2: Solve for such that

 (23) dtqn+1−κΔ(αqn+1+˜pn+1) =0, in ΩT, (24) α∂qn+1∂ν =−∂˜pn+1∂ν on ∂ΩT, (25) (qn+1,1) =Cq,

where .

We note that (23) is the implicit Euler scheme, which is chosen just for the ease of presentation, it can be replaced by other time-stepping schemes. (24) provides a Neumann boundary condition for . Another subtle issue is the role which the initial value plays in the algorithm. Seems is only needed to produce and there is no need to have in order to execute the algorithm. However, to ensure the stability and convergence of the algorithm, it turns out that not only needs to be provided but also must be carefully constructed when the algorithm is discretized (see Sections 3 and 4).

The above algorithm has a couple attractive features. First, it is easy to use. Second, it allows one to make use of any available numerical methods (finite element, finite difference, finite volume, spectral and discontinuous Galerkin) and computer codes for the Stokes problem and the Poisson problem to solve the gel swelling dynamics model (5)–(8). We remark that an almost same model as (5)–(8) also arise from different applications in poroelasticity and soil mechanics and are known as Boit’s consolidation model (cf. [2, 14] and the references therein). [14] proposed and analyzed a standard finite element method which directly approximates under the (restrictive) divergence-free assumption on the initial condition .

This paper consists of four additional sections. In Section 2, we first introduce notation used in this paper. We then present a PDE analysis for the gel swelling dynamics model (5)–(8), which includes deriving a dissipative energy law, establishing existence and uniqueness, and in particular, proving the conservation laws stated in (14)–(16) and (18). In Section 3, we first propose a semi-discrete (in space) finite element discretization for problem (5)–(8) based on the multiphysical reformulation (11)–(13). The well-known Taylor-Hood mixed element and the conforming finite element are used as an example to present the ideas. It is proved that the solution of the semi-discrete method satisfies a discrete energy law which mimics the differential energy law enjoyed by the PDE solution, and the semi-discrete numerical solution also satisfies the conservation laws (14)–(16) and (18). We then derive optimal order error estimates in various norms for the semi-discrete numerical solution. In Section 4, fully discrete finite element methods are constructed by combining the time-stepping scheme of Algorithm and the semi-discrete finite element methods of Section 3. The main results of this section include proving a fully discrete energy law for the numerical solution and establishing optimal order error estimates for the fully discrete finite element methods. Finally, in Section 5, we present some numerical experiments to gauge the efficiency of the proposed approach and methods.

## 2 PDE analysis for problem (5)–(8)

The standard Sobolev space notation is used in this paper, we refer to [4, 6, 18] for their precise definitions. In particular, and denote respectively the standard and inner products. For any Banach space , we let and use to denote its dual space. In particular, we use and to denote the dual products on and , respectively. is a shorthand notation for .

We also introduce the function spaces

 L20(Ω):={q∈L2(Ω);(q,1)=0},X:={v∈H1(Ω);⟨v,ν⟩=0}.

It is well known [18] that the following so-called inf-sup condition holds in the space :

 (26) supv∈X(\rm div\,v,φ)∥∇v∥L2≥α0∥φ∥L2∀φ∈L20(Ω),α0>0.

Throughout the paper, we assume be a bounded polygonal domain such that is an isomorphism; see  [11, 12]. In addition, is used to denote a generic positive constant which is independent of , and the mesh parameters and .

We now give a definition of weak solutions to (5)-(8). {definition} Let , and . Given , a tuple with

 u∈L∞(0,T;H1(Ω)),\rm div% \,u∈H1(0,T;H−1(Ω)),p∈L2(0,T;H1(Ω)),

is called a weak solution to (5)–(8), if there hold for almost every

 (27) ((\rm div\,u)t,φ)dual+κ(∇p,∇φ) =0 ∀φ∈H1(Ω), (28) α(\rm div\,u,\rm div\,v)+β(∇u,∇v)−(p,\rm div\,v) =⟨f,v⟩dual ∀v∈H1(Ω), (29) u(0) =u0.

Similarly, we define weak solutions to problem (11)-(13), (7)-(8).

{definition}

Let , and . Given , a triple with

 u∈L∞(0,T;H1(Ω)), ˜p∈L2(0,T;L2(Ω)), q∈L∞(0,T;L2(Ω))∩H1(0,T;H−1(Ω)), αq+˜p∈L2(0,T;H1(Ω)),

is called a weak solution to (11)–(13), (7)-(8) if there hold for almost every

 (30) β(∇u,∇v)−(˜p,\rm div\,v) =⟨f,v⟩dual ∀v∈H1(Ω), (31) (\rm div\,u,φ) =(q,φ) ∀φ∈L2(Ω), (32) u(0) =u0, (33) (qt,ψ)dual+κ(∇(αq+˜p),∇ψ) =0 ∀ψ∈H1(Ω), (34) q(0) =q0:=\rm div\,u0.
###### Remark 2.1

(a) Clearly, gives back the pressure in the original formulation. What interesting is that both and are only -functions in the spatial variable but their combination is an -function. In other words, the new formulation provides an decomposition for the pressure . It turns out that this decomposition will have a significant numerical impact because it allows one to use low order (hence cheap) finite elements to approximate and but still to be able to approximate the pressure with high accuracy.

(b) (33) implicitly imposes the following boundary condition for :

 (35) α∂q∂ν=−∂˜p∂νon ∂Ω.

Since problem (27)–(29) consists of two linear equations, its solvability should follows easily if we can establish a priori energy estimates for its solutions. The following dissipative energy law just serves that purpose.

{lemma}

Every weak solution of problem (5)–(8) satisfies the following energy law:

 (36) E(t)+κ∫t0∥∇p(s)∥2L2ds=E(0)∀t∈(0,T],

where

 E(t):=12[β∥∇u(t)∥2L2+α∥\rm div\,u(t)∥2L2−2⟨f,u(t)⟩dual].

Moreover,

 (37) ∥(\rm div\,u)t∥L2(H−1)=κ∥Δp∥L2(H−1)≤E(0)12.
{proof}

We first consider the case . Setting in (27) and in (28) yield

 (\rm div\,ut(t),p(t))+κ∥∇p∥2L2 =0, =ddt⟨f,u(t)⟩dual+(p(t),\rm div\,ut(t)).

(36) follows from adding the above two equations and integrating the sum in over the interval for any .

If is only a weak solution, then could not be used as a test function in (28). However, this difficulty can be easily overcome by using as a test function in (28), where denotes a mollification of through a symmetric mollifier (cf. [11, Chapter 7]), and by passing to the limit .

Finally, (37) follows immediately from (27) and (36). The proof is completed.

###### Remark 2.2

Lemma 2 and Theorem 2 can be easily carried over to the reformulated problem (11)–(13), (7)–(8). The only difference is that the energy law (36) now is replaced by the following equivalent energy law:

 (38) J(t)+κ∫t0∥∇[˜p(s)+αq(s)]∥2L2ds=J(0)∀t∈(0,T],

and (37) is replaced by

 (39) ∥qt∥L2(H−1)≤√κ∥∇(˜p+αq)∥L2(L2)≤J(0)12.

Where

 J(t):=12[β∥∇u(t)∥2L2+α∥q(t)∥2L2−2⟨f,u(t)⟩dual].

We also note that the weak solution to problem (11)–(13), (7)–(8) is understood in the sense of Definition 2.

{lemma}

Weak solutions of problem (5)–(8) and problem (11)–(13),(7)–(8) satisfy the conservation laws (14)–(16), and (18).

{proof}

(14) and (16) follows immediately from taking in (33), in both (31) and (27) followed by integrating in and appealing to the divergence theorem.

To prove (15), taking in (28) and using the identities and we get

 α(\rm div\,u,d)+β(∇u,I)−(p,d)=⟨f,x⟩dual

Hence,

 (p,1) =α(\rm div\,u,1)+βd(\rm div\,u,1)−1d⟨f,x⟩dual =(α+βd)(\rm div\,u0,1)−1d⟨f,x⟩dual =cdCq−1d⟨f,x⟩dual,

where and are defined in (17) and (14). So we obtain (15).

Similarly, (18) follows by taking in (30). The proof is complete.

With the help of the above two lemmas, we can show the solvability of problem (5)–(8).

{theorem}

Suppose , and . Then there exists a unique weak solution to (5)–(8) in the sense of Definition 2, and there exists a unique solution to to (11)–(13), (7)–(8).

{proof}

The uniqueness is an immediate consequence of the energy laws (36) and (38) and the conservation laws (14)–(16) and (18).

The existence can be easily proved by using (abstract) Galerkin method and the standard compactness argument (cf. [18]). The energy laws (36) and (38) provide the necessary uniform estimates for the Galerkin approximate solutions. We omit the details because the derivation is standard.

Again, exploiting the linearity of the system, we have the following regularity results for the weak solution.

{theorem}

Let be the weak solution of (30)–(34) Then there holds the following estimates:

 (40) √β∥√t∇ut∥L2(L2)+√α∥√tqt∥L2(L2)+√κ∥√t∇p∥L∞(L2) ≤J(0)12. (41) √β∥t∇ut∥L∞(L2)+√α∥tqt∥L∞(L2)+√κ∥t∇pt∥L2(L2) ≤[2J(0)]12. (42) ∥tqtt∥L2(H−1) ≤[2J(0)]12. (43) κ√α∥√tΔp∥L2(L2)=√α∥√tqt∥L2(L2) ≤J(0)12. (44) κ√α∥tΔp∥L∞(L2)=√α∥tqt∥L∞(L2) ≤[2J(0)]12. (45) ∥∇˜p∥L2(L2)+∥∇˜p∥L∞(L2)+α√κ∥∇q∥L2(L2) +α√κ∥√t∇q∥L∞(L2) ≤C(C˜p,J(0)). (46) β∥Δu∥L∞(L2)=∥∇˜p∥L∞(L2) ≤C(C˜p,J(0)).

Where denotes some positive constant which depends on and .

{proof}

Differentiating each of equations (30), (31) and (33) with respect to yields (note that is assumed to be independent of )

 (47) =0 ∀v∈H1(Ω), (48) (\rm div\,ut,φ) =(qt,φ) ∀φ∈L2(Ω), (49) (qtt,ψ)dual+κ(∇(αq+˜p)t,∇ψ) =0 ∀ψ∈H1(Ω).

Taking in (47), in (48), in (33), and adding the resulted equations we get

 tβ∥∇ut∥2L2+tα∥qt∥2L2+κ2ddt(t∥∇p∥2L2)=κ2∥∇p∥2L2.

Integrating in over gives (40).

Alternatively, setting in (47), in (48) (after differentiating in one more time), in (49), and adding the resulted equations we have

 12ddt[t2β∥∇ut∥2L2+t2α∥qt∥2L2]+t2κ∥∇pt∥2L2=tβ∥∇ut∥2L2+t∥qt∥2L2.

Integrating in over and using (40) yield (41).

(42)–(44) are immediate consequences of (49), (40), and (41).

Choose such that . Then . By an equivalent version of the inf-sup condition (cf. [1]) we have

 ∥∇(˜p−ϕ0)∥L2 ≤α−10supv∈H10(Ω)(∇(˜p−ϕ0),v)∥∇v∥L2 ≤α−10supv∈H10(Ω)−β(∇u,∇v)∥∇v∥L2+α−10√d∥ϕ0∥L2 ≤α−10[β∥∇u∥L2+√d∥ϕ0∥L2].

Hence,

 ∥∇˜p∥L2≤∥∇ϕ0∥L2+α−10[β∥∇u∥L2+√d∥ϕ0∥L2],

which together with (38) and (40) imply that (recall )

 ∥∇˜p∥L2(L2)≤√T∥∇˜p∥L∞(L2)≤√T∥∇ϕ0∥L2+α−10√T[βJ(0)12+√d∥ϕ0∥L2], α√κ∥∇q∥L2(L2)≤J(0)12+√κT∥∇ϕ0∥L2+α−10√κT[βJ(0)12+√d∥ϕ0∥L2], α√κ∥√t∇q∥L∞(L2)≤J(0)12+√κT∥∇ϕ0∥L2+α−10√κT[βJ(0)12+√d∥ϕ0∥L2].

Hence, (45) holds.

Finally, (46) follows immediately from the equation and (45). The proof is complete.

## 3 Semi-discrete finite element methods in space

The goal of this section is to present the ideas and specific semi-discrete finite element methods for discretizing the variational problem (30)–(34) in space based on the multiphysical (deformation and diffusion) approach. That is, we shall approximate using a (stable) Stokes solver and approximate by a convergent diffusion equation solver. The combination of the Taylor-Hood mixed finite element [1] and the conforming finite element is chosen as a specific example to present the ideas.

### 3.1 Formulation of finite element methods

Assume () is a polygonal domain. Let is a quasi-uniform triangulation or rectangular partition of with mesh size such that . Also, let be a stable mixed finite element pair, that is, and satisfy the inf-sup condition

 (50) supvh∈Vh∩X(%divvh,φh)∥∇vh∥L2≥β0∥φh∥L2∀φh∈Mh∩L20(Ω),β0>0.

A well-known example that satisfies (50) is the following so-called Taylor-Hood element (cf. [1]):

 Vh ={vh∈C0(¯¯¯¯Ω);vh|K∈P2(K)  ∀K∈Th}, Mh ={φh∈C0(¯¯¯¯Ω);φh|K∈P1(K)  ∀K∈Th}.

In the sequel, we shall only present the analysis for the Taylor-Hood element, but remark that the analysis can be readily extended to other stable combinations. On the other hand, constant pressure space is not recommended because that would result in no rate of convergence for the approximation of the pressure (see Section 3.3).

Approximation space for variable can be chosen independently, any piecewise polynomial space is acceptable as long as . Especially, can be chosen as a fully discontinuous piecewise polynomial space, although it is more convenient to choose to be a continuous (resp. discontinuous) space if is a continuous (resp. discontinuous) space. The most convenient and efficient choice is , which will be adopted in the remaining of this paper.

We now ready to state our semi-discrete finite element method for problem (30)–(34). Let . We seek and such that for all there hold

 (51) β(∇uh,∇vh)−(˜ph,\rm div\,vh) =⟨f,vh⟩dual ∀vh∈Vh, (52) (\rm div\,uh,φh) =(qh,φh) ∀φh∈Mh, (53) (∇uh(0),∇wh)=(∇u0,∇wh),⟨uh(0),ν⟩ =⟨u0,ν⟩ ∀wh∈Vh, (54) (qht,ψh)dual+κ(∇(˜ph+αqh),∇ψh) =0 ∀ψh∈Wh, (55) (qh(0),χh) =(q0,χh) ∀χh∈Wh.

Where denotes the time derivative of .

###### Remark 3.1

(a) We note that (53) and (55) defines and (see their definitions below). It is easy to see that in general although . But it is not hard to enforce by defining slightly differently (in the case of continuous ) or simply substituting (55) by (in the case of discontinuous ), even such a modification is not necessary for the sake of convergence (see Section 3.3). We also note that