Generalized multiscale finite element method for a strain-limiting nonlinear elasticity model

TL;DR

This paper develops a multiscale finite element method for nonlinear elasticity with strain-limiting behavior, efficiently capturing multiscale features and nonlinearities through offline and online basis functions.

Contribution

It introduces a GMsFEM framework tailored for strain-limiting nonlinear elasticity, incorporating adaptive basis function updates based on material property changes.

Findings

Accurate solutions with few basis functions per coarse region.

Adaptive basis updates improve solution accuracy.

Method effectively handles multiscale nonlinear elasticity problems.

Abstract

In this paper, we consider multiscale methods for nonlinear elasticity. In particular, we investigate the Generalized Multiscale Finite Element Method (GMsFEM) for a strain-limiting elasticity problem. Being a special case of the naturally implicit constitutive theory of nonlinear elasticity, strain-limiting relation has presented an interesting class of material bodies, for which strains remain bounded (even infinitesimal) while stresses can become arbitrarily large. The nonlinearity and material heterogeneities can create multiscale features in the solution, and multiscale methods are therefore necessary. To handle the resulting nonlinear monotone quasilinear elliptic equation, we use linearization based on the Picard iteration. We consider two types of basis functions, offline and online basis functions, following the general framework of GMsFEM. The offline basis functions depend nonlinearly on the solution. Thus, we design an indicator function and we will recompute the offline basis functions when the indicator function predicts that the material property has significant change during the iterations. On the other hand, we will use the residual based online basis functions to reduce the error substantially when updating basis functions is necessary. Our numerical results show that the above combination of offline and online basis functions is able to give accurate solutions with only a few basis functions per each coarse region and updating basis functions in selected iterations.

Equations

$$\mathbf{E}=\frac{\mathbf{T}}{1+\beta (\mathbf{x})|\mathbf{T}|}.$$

$$\mathbf{T}=\frac{\mathbf{E}}{1-\beta (\mathbf{x})|\mathbf{E}|}.$$

$$\mathbf{E}:=\frac{1}{2}(\nabla \mathbf{u}+\nabla {\mathbf{u}}^{\text{T}}).$$

$$\mathbf{E}=\mathbf{D}(\mathbf{u})=\mathbf{Du}={\nabla }_s\mathbf{u}.$$

$$\mathbf{T}=\frac{\mathbf{D}(\mathbf{u})}{1-\beta (\mathbf{x})|\mathbf{D}(\mathbf{u})|}.$$

$$|\mathbf{E}|=\frac{|\mathbf{T}|}{1+\beta (\mathbf{x})|\mathbf{T}|}<\frac{1}{\beta (\mathbf{x})}.$$

$$\Vert \mathbf{v}{\Vert }_{{\mathbf{W}}_0^{1,2}(\Omega )}=(\Vert \mathbf{v}{\Vert }_{{\mathbf{L}}^2(\Omega )}^2+\Vert \nabla \mathbf{v}{\Vert }_{{\mathbb{L}}^2(\Omega )}^2{)}^{\frac{1}{2}}.$$

$$\Vert \mathbf{b}{\Vert }_{{\mathbf{H}}^{-1}(\Omega )}=\underset{\mathbf{v}\in {\mathbf{H}}_0^1(\Omega )}{\sup }\frac{|⟨\mathbf{b},\mathbf{v}⟩|}{\Vert \mathbf{v}{\Vert }_{{\mathbf{H}}_0^1(\Omega )}}.$$

$$\boldsymbol{f}\in\boldsymbol{H}^{1}_{*}(\Omega)=\left\{\boldsymbol{g}\in\boldsymbol{H}^{1}(\Omega)\biggr{|}\int_{\Omega}\boldsymbol{g}\,dx=\boldsymbol{0}\right\}\subset\boldsymbol{L}^{2}(\Omega)\subsetneq\boldsymbol{H}^{-1}(\Omega)\,.$$

$$\begin{array}{rl}-\text{div}(\mathbf{T})& =\mathbf{f}\text{in }\Omega ,\\ \mathbf{Du}& =\frac{\mathbf{T}}{1+\beta (\mathbf{x})|\mathbf{T}|}\text{in }\Omega ,\\ \mathbf{u}& =0\text{on }\partial {\Omega }_D,\\ \mathbf{Tn}& =\mathbf{G}\text{on }\partial {\Omega }_T,\end{array}$$

$$\mathbf{S}\text{ as }\mathbf{T}\text{and}\mathbf{D}(\mathbf{u})=\mathbf{Du}\text{ as }\mathbf{E}=\mathbf{E}(\mathbf{u}).$$

$$-\text{div}\left(\frac{\mathbf{D}(\mathbf{u})}{1-\beta (\mathbf{x})|\mathbf{D}(\mathbf{u})|}\right)$$

$$\mathbf{u}$$

$$\kappa (\mathbf{x},|\mathbf{D}(\mathbf{u})|)=\frac{1}{1-\beta (\mathbf{x})|\mathbf{D}(\mathbf{u})|},\mathbf{a}(\mathbf{x},\mathbf{D}(\mathbf{u}))=\kappa (\mathbf{x},|\mathbf{D}(\mathbf{u})|)\mathbf{D}(\mathbf{u}),$$

$$0\le |\mathbf{D}(\mathbf{u})|<\frac{1}{M}\le \frac{1}{\beta (\mathbf{x})}\le \frac{1}{m}.$$

$$\mathbf{F}(\mathbf{\xi })=\frac{\mathbf{\xi }}{1-\beta (\mathbf{x})|\mathbf{\xi }|},$$

$$0\le |\mathbf{\xi }|<\frac{1}{M}\le \frac{1}{\beta (\mathbf{x})}.$$

$$\mathcal{Z}:=\left\{\boldsymbol{\zeta}\in\mathbb{R}^{2\times 2}\;\biggr{|}\;0\leq|\boldsymbol{\zeta}|<\dfrac{1}{M}\right\}\,.$$

$$\mathbf{\xi }\in \mathcal{Z}\mapsto \mathbf{F}(\mathbf{\xi }):=\frac{\mathbf{\xi }}{1-\beta (\mathbf{x})|\mathbf{\xi }|}\in {\mathbb{R}}^{2\times 2}.$$

$$|\mathbf{F}({\mathbf{\xi }}_1)-\mathbf{F}({\mathbf{\xi }}_2)|$$

$$(\mathbf{F}({\mathbf{\xi }}_1)-\mathbf{F}({\mathbf{\xi }}_2))\cdot ({\mathbf{\xi }}_1-{\mathbf{\xi }}_2)$$

$$\mathbf{F}({\mathbf{\xi }}_1)-\mathbf{F}({\mathbf{\xi }}_2)$$

$$={\int }_0^1\frac{d}{dt}\left(\frac{t{\mathbf{\xi }}_1+(1-t){\mathbf{\xi }}_2}{1-\beta (\mathbf{x})|t{\mathbf{\xi }}_1+(1-t){\mathbf{\xi }}_2|}\right)dt$$

$$={\int }_0^1\frac{{\mathbf{\xi }}_1-{\mathbf{\xi }}_2}{(1-\beta (\mathbf{x})|t{\mathbf{\xi }}_1+(1-t){\mathbf{\xi }}_2|{)}^2}dt.$$

$$|t{\mathbf{\xi }}_1+(1-t){\mathbf{\xi }}_2|\le \text{max}\{|{\mathbf{\xi }}_1|,|{\mathbf{\xi }}_2|\}\le |{\mathbf{\xi }}_1|+|{\mathbf{\xi }}_2|.$$

$$|\mathbf{F}({\mathbf{\xi }}_1)-\mathbf{F}({\mathbf{\xi }}_2)|$$

$$(\mathbf{F}({\mathbf{\xi }}_1)-\mathbf{F}({\mathbf{\xi }}_2))\cdot ({\mathbf{\xi }}_1-{\mathbf{\xi }}_2)={\int }_0^1\frac{|{\mathbf{\xi }}_1-{\mathbf{\xi }}_2{|}^2}{(1-\beta (\mathbf{x})|t{\mathbf{\xi }}_1+(1-t){\mathbf{\xi }}_2|{)}^2}dt\ge |{\mathbf{\xi }}_1-{\mathbf{\xi }}_2{|}^2.$$

$$\mathcal{U}=\{\mathbf{w}\in {\mathbf{H}}^1(\Omega )|\mathbf{Dw}\in {\mathcal{Z}}^{\prime }\},$$

$${\int }_{\Omega }\mathbf{a}(\mathbf{x},\mathbf{D}\mathbf{u})\cdot \mathbf{D}\mathbf{v}dx={\int }_{\Omega }\mathbf{f}\cdot \mathbf{v}dx,\forall \mathbf{v}\in \mathbf{V}.$$

$$\text{ Find }(\mathbf{u},\mathbf{Du})\in \mathbf{V}\times {\mathbb{L}}^{\infty }(\Omega ),\text{that is, find }\mathbf{u}\in \mathbf{V}\text{ such that }(\text{ref}w8.2)\text{ holds for each }\mathbf{v}\in \mathbf{V}.$$

$${\int }_{\Omega }\mathbf{T}\cdot \mathbf{D}(\mathbf{w})dx$$

$$\mathbf{D}(\mathbf{u})$$

$$\mathbf{u}$$

Full text

Keywords. Generalized multiscale finite element method; Strain-limiting nonlinear elasticity; Adaptivity; Residual based online multiscale basis functions

Mathematics Subject Classification. 65N30, 65N99

Shubin Fu $\cdot$ Eric Chung

Department of Mathematics, The Chinese University of Hong Kong, Shatin, Hong Kong

E-mail: shubinfu89@gmail.com (Shubin Fu); tschung@math.cuhk.edu.hk (Eric Chung)

Tina Mai

Institute of Research and Development, Duy Tan University, Da Nang 550000, Vietnam

E-mail, corresponding author: maitina@duytan.edu.vn; tinagdi@gmail.com

1. Introduction

Even though multiscale methods for linear equations are positively growing up, their applications to nontrivial nonlinear problems are still hard. In addition, many nonlinear elastic materials contain multiple scales and high contrast. To overcome the challenge from nonlinearity, the main idea is to linearize it in Picard iteration, as employed in [16] and references therein. To deal with the difficulties from multiple scales and high contrast, we apply the GMsFEM ([15]) for the linear equation at the current iteration.

The motivation for our nonlinear elasticity problem is a remarkable trend of studying nonlinear responses of materials based on the recently developed implicit constitutive theory (see [22, 23, 25, 24]). As Rajagopal notes, the theory offers a framework for establishing nonlinear, infinitesimal strain theories for elastic-like (non-dissipative) material behavior. This setting is different from classical Cauchy and Green approaches for modeling elasticity which, under the condition of infinitesimal strains, lead to traditional linear models. Moreover, it is significant that the implicit constitutive theory provides a firm theoretical foundation for modeling fluid and solid mechanics variously, in engineering, physics, and chemistry.

We consider herein a special sub-class of the implicit constitutive theory in solids, namely, the strain-limiting theory, for which the linearized strain remains bounded, even when the stress is very large. (In the traditional linear model, the stress blows up when the strain blows up and vice versa.) Therefore, the strain-limiting theory is helpful for describing the behavior of fracture, brittle materials near crack tips or notches (e.g. crystals), or intensive loads inside the material body or on its boundary. Either case leads to stress concentration, regardless of the small gradient of the displacement (and thus infinitesimal strain). Our considering materials are science-non-fiction and physically meaningful. These materials can sustain infinite stresses, and do not break (because of the boundedness of strains).

To solve the multiple scales, rather than direct numerical simulations on fine scales, model reduction methods are proposed in literature to reduce the computational expensiveness. Generally, model reduction techniques include upscaling and multiscale approaches. On coarse grid (which is much larger than fine grid), upscaling methods involve upscaling the media properties based on homogenization to capture macroscopic behavior, whereas multiscale methods additionally need precomputed multiscale basis functions.

Within the framework of multiscale methods, the multiscale finite element method (MsFEM), as in [18, 14], has effectively proved notable success in a variety of practical applications, but it requires scale separation nevertheless. To overcome this requirement, we use the generalized multiscale finite element method (GMsFEM), as in [15], to systematically construct multiple multiscale basis. More specifically, in the GMsFEM, the computation is divided into two stages: offline and online. One constructs, in the offline stage, a small dimension space, which can be used effectively in the online stage to construct multiscale basis functions, to solve the problem on coarse grid. The construction of offline and online spaces is based on the selection of local spectral problems as well as the selection of the snapshot space. In [8], the GMsFEM was applied to handle the linear elasticity problem. In [6], a model reduction method was introduced to solve nonlinear monotone elliptic equations. In [2], thanks to [16] (for handling nonlinearities), the GMsFEM was used to solve nonlinear poroelasticity problems.

Here, our paper will combine the ideas of Picard iteration and the GMsFEM in [9] to solve a strain-limiting nonlinear elasticity problem [21, 20]. At each Picard iteration, we will either apply the offline GMsFEM or the residual based online adaptive GMsFEM. In the latter approach, we study the proposed online basis construction in conjunction with adaptivity ([7, 9, 12]), which means that online basis functions are added in some selected regions. Adaptivity is an important step to obtain an effective local multiscale model reduction as it is crucial to reduce the cost of online multiscale basis computations. More specifically, adaptive algorithm allows one to add more basis functions in neighborhoods with more complexity without using a priori information. Given a sufficient number of initial basis functions in the offline space, our numerical results show that the adaptive addition of online basis functions substantially decreases the error, accelerates the convergence, and reduces computational cost of the GMsFEM.

Our strategy is that after some Picard iterations, when the relative change of the permeability coefficient is larger than a given fixed tolerance, we need to update (either offline or online, context-dependently) basis functions. This updating procedure ends when we obtain desired error, and these new basis functions is kept the same in the next Picard iterations until we need to compute them again.

The next Section contains the mathematical background of our considering strain-limiting nonlinear elasticity problem. Section 3 is for some preliminaries about the GMsFEM, including fine-scale discretization and Picard iteration for linearization. Section 4 is devoted to the general idea of GMsFEM, including some existing results regarding offline GMsFEM, for the current nonlinear elasticity problem. Section 5 is about the existing method of residual based online adaptive GMsFEM for computing online multiscale basis functions, in our context. In Section 6, some numerical examples will be shown. The last Section 7 is for conclusions.

2. Formulation of the problem

2.1. Input problem and classical formulation

Let us consider, in dimension two, a nonlinear elastic composite material $\Omega=\Omega^{1}\times\Omega^{2}\in\mathbb{R}\times\mathbb{R}$.

We assume that the material is being at a static state after the action of body forces $\boldsymbol{f}:\Omega\to\mathbb{R}^{2}$ and traction forces $\boldsymbol{G}:\partial\Omega_{T}\to\mathbb{R}^{2}$. The boundary of the set $\Omega$ is denoted by $\partial\Omega$, which is Lipschitz continuous, consisting of two parts $\partial\Omega_{T}$ and $\partial\Omega_{D}$, where the displacement $\boldsymbol{u}:\Omega\to\mathbb{R}^{2}$ is given on $\partial\Omega_{D}$. We are considering the strain-limiting model of the form (as in [21])

[TABLE]

Equivalently,

[TABLE]

In the equations (2.1) and (2.2), $\boldsymbol{T}$ denotes the Cauchy stress $\boldsymbol{T}:\Omega\to\mathbb{R}^{2\times 2}$; whereas, $\boldsymbol{E}$ denotes the classical linearized strain tensor,

[TABLE]

We can write

[TABLE]

Then, by (2.2), it follows that

[TABLE]

The strain-limiting parameter function is denoted by $\beta(\boldsymbol{x})$, which depends on the position variable $\boldsymbol{x}=(x^{1},x^{2})$. We notice from (2.1) that

[TABLE]

This means that $\displaystyle\frac{1}{\beta(\boldsymbol{x})}$ is the upper-bound on $|\boldsymbol{E}|$, and taking sufficiently large $\beta(\boldsymbol{x})$ raises the limiting-strain small upper-bound, as desired. However, we avoid $\beta(\boldsymbol{x})\to\infty$. (If $\beta(\boldsymbol{x})\to\infty$, then $|\boldsymbol{E}|<\displaystyle\frac{1}{\beta(\boldsymbol{x})}\to 0$, a contradiction.) For the analysis of the problem, $\beta(\boldsymbol{x})$ is assumed to be smooth and have compact range $0<m\leq\beta(\boldsymbol{x})\leq M$. Here, $\beta(\boldsymbol{x})$ is chosen so that the strong ellipticity condition holds (see [21]), that is, $\beta(\boldsymbol{x})$ is large enough, to prevent bifurcations arising in numerical simulations.

2.2. Function space

The preliminaries are the same as in [13]. Latin indices vary in the set $\{1,2\}$. The space of functions, vector fields in $\mathbb{R}^{2}$, and $2\times 2$ matrix fields defined over $\Omega$ are respectively denoted by italic capitals (e.g. $L^{2}(\Omega)$), boldface Roman capitals (e.g. $\boldsymbol{V}$), and special Roman capitals (e.g. $\mathbb{S}$). The space of symmetric matrices of order 2 is denoted by $\mathbb{S}^{2}$. The subscript $s$ appended to a special Roman capital denotes a space of symmetric matrix fields.

Our considering space is $\boldsymbol{V}:=\boldsymbol{H}_{0}^{1}(\Omega)=\boldsymbol{W}_{0}^{1,2}(\Omega)$. However, the techniques here can be used in more general space $\boldsymbol{W}_{0}^{1,p}(\Omega)$, where $2\leq p<\infty$. The reason we consider the space $\boldsymbol{W}_{0}^{1,2}(\Omega)$ is that we can characterize displacements that vanish on the boundary $\partial\Omega$ of $\Omega$. The dual space (also called the adjoint space), which consists of continuous linear functionals on $\boldsymbol{H}_{0}^{1}(\Omega)$, is denoted by $\boldsymbol{H}^{-1}(\Omega)$, and the value of a functional $\boldsymbol{b}\in\boldsymbol{H}^{-1}(\Omega)$ at a point $\boldsymbol{v}\in\boldsymbol{H}_{0}^{1}(\Omega)$ is denoted by $\langle\boldsymbol{b},\boldsymbol{v}\rangle$.

The Sobolev norm $\|\cdot\|_{\boldsymbol{W}_{0}^{1,2}(\Omega)}$ is of the form

[TABLE]

Here, $\|\boldsymbol{v}\|_{\boldsymbol{L}^{2}(\Omega)}:=\||\boldsymbol{v}|\|_{\boldsymbol{L}^{2}(\Omega)}\,,$ where $|\boldsymbol{v}|$ denotes the Euclidean norm of the 2-component vector-valued function $\boldsymbol{v}$; and $\|\nabla\boldsymbol{v}\|_{\mathbb{L}^{2}(\Omega)}:=\||\nabla\boldsymbol{v}|\|_{\mathbb{L}^{2}(\Omega)}\,,$ where $|\nabla\boldsymbol{v}|$ denotes the Frobenius norm of the $2\times 2$ matrix $\nabla\boldsymbol{v}$. We recall that the Frobenius norm on $\mathbb{L}^{2}(\Omega)$ is defined by $|\boldsymbol{X}|^{2}:=\boldsymbol{X}\cdot\boldsymbol{X}=\text{tr}(\boldsymbol{X}^{\text{T}}\boldsymbol{X})\,.$

The dual norm to $\|\cdot\|_{\boldsymbol{H}_{0}^{1}(\Omega)}$ is $\|\cdot\|_{\boldsymbol{H}^{-1}(\Omega)}$, i.e.,

[TABLE]

Let $\Omega$ be a bounded, simply connected, open, Lipschitz, convex domain of $\mathbb{R}^{2}$. Let

[TABLE]

be bounded in $\boldsymbol{L}^{2}(\Omega)$. We consider the following problem: Find $\boldsymbol{u}\in\boldsymbol{H}^{1}(\Omega)$ and $\boldsymbol{T}\in\mathbb{L}^{1}(\Omega)$ such that

[TABLE]

where $\boldsymbol{n}$ stands for the outer unit normal vector to the boundary of $\Omega$.

Benefiting from the notations in [3], we will write

[TABLE]

The considering model (2.1) is compatible with the laws of thermodynamics [26, 27], which means that the class of materials are non-dissipative and are elastic.

We assume $\partial\Omega_{T}=\emptyset$. Using (2.7), we write the considering formulation in the form of displacement problem: Find $\boldsymbol{u}\in\boldsymbol{H}_{0}^{1}(\Omega)$ such that

[TABLE]

We denote

[TABLE]

in which $\boldsymbol{u}(\boldsymbol{x})\in\boldsymbol{W}_{0}^{1,2}(\Omega)$. In this setting, $\boldsymbol{a(\boldsymbol{x},\boldsymbol{\xi})}\in\mathbb{L}^{1}(\Omega)$, $\boldsymbol{\xi}\in\mathbb{L}^{\infty}(\Omega)$, as in [1].

2.3. Existence and uniqueness

Notice from (2.5) that $\displaystyle\frac{1}{\beta(\boldsymbol{x})}$ is the upper-bound on $|\boldsymbol{D}(\boldsymbol{u})|$, and $0<m\leq\beta(\boldsymbol{x})\leq M$ such that

[TABLE]

We define

[TABLE]

where $\boldsymbol{\xi}\in\mathbb{L}^{\infty}(\Omega)$, and

[TABLE]

With these facts and thanks to [3], we derive the following results, which were also stated in [4] (p. 19).

Lemma 2.1**.**

Let

[TABLE]

For any $\boldsymbol{\xi}\in\mathcal{Z}$ such that $0\leq|\boldsymbol{\xi}|<\dfrac{1}{M}$, consider the mapping

[TABLE]

Then, for each $\boldsymbol{\xi}_{1},\boldsymbol{\xi}_{2}\in\mathcal{Z}$, we have

[TABLE]

Proof.

We present here the proof in details, thanks to [3]. Notice first that

[TABLE]

For the proof of (2.15), we observe that

[TABLE]

Therefore,

[TABLE]

Then (2.15) follows.

For (2.16), we notice that $(1-\beta(\boldsymbol{x})|t\boldsymbol{\xi}_{1}+(1-t)\boldsymbol{\xi}_{2}|)^{2}\leq 1$. Thus,

[TABLE]

∎

Remark 2.2**.**

The condition (2.16) also implies that $\boldsymbol{F}(\boldsymbol{\xi})$ is a monotone operator in $\boldsymbol{\xi}$.

Remark 2.3**.**

For the rest of the paper, without confusion, we will use the condition $\boldsymbol{\xi}\in\mathbb{L}^{\infty}(\Omega)$ with the meaning that \boldsymbol{\xi}\in\mathcal{Z}^{\prime}=\left\{\boldsymbol{\zeta}\in\mathbb{L}^{\infty}(\Omega)\;\biggr{|}\;0\leq|\boldsymbol{\zeta}|<\dfrac{1}{M}\right\}.

2.3.1. Weak formulation

Let

[TABLE]

with the given $\mathcal{Z}^{\prime}$ in Remark 2.3.

Remark 2.4**.**

For the rest of the paper, without confusion, we will use the condition $\boldsymbol{u},\boldsymbol{v}\in\boldsymbol{H}_{0}^{1}(\Omega)\text{ or }\boldsymbol{H}^{1}(\Omega)$ (context-dependently) with the meaning that $\boldsymbol{u},\boldsymbol{v}\in\mathcal{U}$.

Now, for $\boldsymbol{u}\in\boldsymbol{V}=\boldsymbol{H}_{0}^{1}(\Omega)$, we multiply equation (2.8) by $\boldsymbol{v}\in\boldsymbol{V}$ and integrate the equation with respect to $\boldsymbol{x}$ over $\Omega$. Integrating the first term by parts and using the condition $\boldsymbol{v}=\boldsymbol{0}$ on $\partial\Omega$, we obtain

[TABLE]

By the weak (often called generalized) formulation of the boundary value problem (2.8)-(2.9), we interpret the problem as follows:

[TABLE]

2.3.2. Existence and uniqueness

In [3], the existence and uniqueness of weak solution $(\boldsymbol{u},\boldsymbol{T})$ for (2.7) have been proved. Also, see [28], for further reference.

Similarly, in our paper, we consider the problem: Find $(\boldsymbol{u},\boldsymbol{T})\in\boldsymbol{H}_{0}^{1}(\Omega)\times\mathbb{L}^{1}(\Omega)$ such that

[TABLE]

This problem is equivalent to problem (2.19).

As noticed in [4] (Section 4.3), the identity (2.1) can be equivalently rewritten as

[TABLE]

where $\boldsymbol{a}$ is a uniformly monotone operator (2.16) with at most linear growth at infinity (2.15). Hence, the existence and uniqueness of the solution $\boldsymbol{T}\in\mathbb{L}^{1}(\Omega)$ and $\boldsymbol{u}\in\boldsymbol{H}_{0}^{1}(\Omega)$ (or, in higher regularity, $\boldsymbol{T}\in\mathbb{L}^{2}(\Omega)$ and $\boldsymbol{u}\in\boldsymbol{W}_{0}^{1,2}(\Omega)$) to (2.20) - (2.22), or $\boldsymbol{u}\in\boldsymbol{W}_{0}^{1,2}(\Omega)$ to (2.18), is guaranteed by [1].

In the case with $\beta(\boldsymbol{x})$ in (2.7), these results are still valid, arriving from similar argument as in [3], thanks to Lemma 2.1.

3. Fine-scale discretization and Picard iteration for linearization

The solution $\boldsymbol{u}\in\boldsymbol{V}$ to (2.8) satisfies

[TABLE]

where

[TABLE]

Starting with an initial guess $\boldsymbol{u}^{0}=\boldsymbol{0}$, to solve equation (2.8), we will linearize it by Picard iteration, that is, we solve

[TABLE]

where superscripts involving $n\text{ }(\geq 0)$ denote respective iteration levels.

To discretize (3.3), we next introduce the notion of fine and coarse grids. Let $\mathcal{T}^{H}$ be a conforming partition of the domain $\Omega$. We call $H$ the coarse mesh size and $\mathcal{T}^{H}$ the coarse grid. Each element of $\mathcal{T}^{H}$ is called a coarse grid block (patch). We denote by $N_{v}$ the total number of interior vertices of $\mathcal{T}^{H}$ and $N$ the total number of coarse blocks. Let $\{\boldsymbol{x}_{i}\}^{N_{v}}_{i=1}$ be the set of vertices in $\mathcal{T}^{H}$ and $w_{i}=\cup\{K_{j}\in\mathcal{T}^{H}\;|\;\boldsymbol{x}_{i}\in\bar{K_{j}}\}$ be the neighborhood of the node $\boldsymbol{x}_{i}$. The conforming refinement of the triangulation $\mathcal{T}^{H}$ is denoted by $\mathcal{T}_{h}$, which is called the fine grid, where $h>0$ is the fine mesh size. We assume that $h$ is very small so that the fine-scale solution $\boldsymbol{u}_{h}$ (to be founded in the next paragraph) is sufficiently close to the exact solution. The main goal of this paper is to find a multiscale solution $\boldsymbol{u}_{\textup{ms}}$ which is a good approximation of the fine-scale solution $\boldsymbol{u}_{h}$. This is the reason why the GMsFEM is used to obtain the multiscale solution $\boldsymbol{u}_{\textup{ms}}$.

On the fine grid $\mathcal{T}_{h}$, we will approximate the solution of (3.1), denoted by $\boldsymbol{u}_{h}$ (or $\boldsymbol{u}$ for brevity). To fix the notation, we will use Picard iteration and the first-order (linear) finite elements for the computation of the fine-scale solution $\boldsymbol{u}_{h}$. In particular, we let $\boldsymbol{V}_{h}\text{ }(\subset\boldsymbol{V}=\boldsymbol{H}_{0}^{1}(\Omega))$ be the first-order Galerkin finite element basis space with respect to the fine grid $\mathcal{T}_{h}$. Toward presenting the details of the Picard iteration algorithm, we define the bilinear form $a(\cdot,\cdot;\cdot)$

[TABLE]

and the functional $J(\cdot)$

[TABLE]

Given $\boldsymbol{u}_{h}^{n}$, the next approximation $\boldsymbol{u}_{h}^{n+1}$ is the solution of the linear elliptic equation

[TABLE]

This is an approximation of the linear equation

[TABLE]

We reformulate the iteration (3.7) in a matrix form. That is, we define $\boldsymbol{A}^{n}$ by

[TABLE]

and define vector $\boldsymbol{b}$ by

[TABLE]

Then, in $\boldsymbol{V}_{h}$, the equation (3.7) can be rewritten in the following matrix form:

[TABLE]

Furthermore, at the $(n+1)$-th Picard iteration, we can solve (3.11) for the multiscale solution $\boldsymbol{u}^{n+1}_{\textup{ms}}\in\boldsymbol{V}^{n}_{\textup{ms}}$ by using the GMsFEM (to be discussed in the next Sections 4 and 5), with multiscale basis functions in $\boldsymbol{V}^{n}_{\textup{ms}}$ computed for $|\boldsymbol{D}\boldsymbol{u}_{\textup{ms}}^{n}|$ in each coarse region $w_{i},i=1,\cdots,N$.

Each of $\boldsymbol{u}_{h}$ and $\boldsymbol{u}_{\textup{ms}}$ is computed in a separate Picard iteration procedure, whose termination criterion is that the relative $\boldsymbol{L}^{2}$ difference is less than $\delta_{0}$, which can be found in Subsection 5.3 and Section 6 ($\delta_{0}=10^{-7}$).

4. GMsFEM for nonlinear elasticity problem

4.1. Overview

We will construct the offline and online spaces. Being motivated by [2], we will concentrate on the effects of the nonlinearities. From the linearized formulation (3.8), we can define offline and online basis functions following the general framework of GMsFEM.

Given $\boldsymbol{u}^{n}$ (which can stand for either $\boldsymbol{u}_{h}^{n}$ or $\boldsymbol{u}^{n}_{\textup{ms}}$, context-dependently), at the current $(n+1)$-th Picard iteration, we will obtain the fine-scale solution $\boldsymbol{u}_{h}^{n+1}\in\boldsymbol{V}_{h}$ by solving the variational problem

[TABLE]

where

[TABLE]

At the $n$-th Picard iteration, we equip the space $\boldsymbol{V}_{h}$ with the energy norm $\|\boldsymbol{v}\|^{2}_{\boldsymbol{V}_{h}}=a_{n}(\boldsymbol{v},\boldsymbol{v})$.

4.2. General idea of GMsFEM

For details of GMsFEM, we refer the readers to [16, 12, 9, 5]. In this paper, at the current $n$-th Picard iteration, we will consider the continuous Galerkin (CG) formulation, having a similar form to the fine-scale problem (4.1).

First, we start with constructing snapshot functions. Then, by solving a class of specific spectral problems in that snapshot space, for each coarse node $\boldsymbol{x}_{i}$, we will obtain a set of multiscale basis functions $\{\boldsymbol{\psi}_{k}^{i}\,|\,k=1,2,\cdots,l_{i}\}$, such that each $\boldsymbol{\psi}_{k}^{i}=\boldsymbol{\psi}_{k}^{w_{i}}$ is supported on the coarse neighborhood $w_{i}$. Furthermore, the basis functions satisfy a partition of unity property, that is, there exist coefficients $\alpha^{i}_{k}$ such that $\sum_{i=1}^{N_{v}}\sum_{k=1}^{l_{i}}\alpha^{i}_{k}\boldsymbol{\psi}^{i}_{k}=1$. With the constructed basis functions, their linear span (over $i=1,\cdots,N_{v},k=1,\cdots l_{i}$) defines the approximate space $\boldsymbol{V}^{m}_{\textup{ms}}$ (at the $m$-th inner iteration, to be specified in Section 5). The GMsFEM solution $\boldsymbol{u}^{m}_{\textup{ms}}\in\boldsymbol{V}^{m}_{\textup{ms}}$ can then be obtained via CG global coupling, which is given through the variational form

[TABLE]

where $a_{n}$ is defined by (4.2).

In summary, one observes that the key component of the GMsFEM is the construction of local basis functions. First, we will use only the so called offline basis functions, which can be computed in the offline stage. Second, to improve the accuracy of the multiscale approximation, we will construct additional online basis functions that are problem-dependent and computed locally and adaptively, based on the offline basis functions and some local residuals. As in [9], our results show that the combination of both offline and online basis functions will give a rapid convergence of the multiscale solution $\boldsymbol{u}^{m}_{\textup{ms}}$ to the fine-scale solution $\boldsymbol{u}_{h}$.

4.3. Construction of offline multiscale basis functions

At the current $n$-th Picard iteration, we will present the construction of the offline basis functions ([9]). We start with constructing, for each coarse subdomain $w_{i}$, a snapshot space $\boldsymbol{V}^{i}_{\textup{snap}}$. For simplicity, the index $i$ can be omitted when there is no confusion. The snapshot space $\boldsymbol{V}^{i}_{\textup{snap}}$ is a set of functions defined on $w_{i}$ and contains all or most necessary components of the fine-scale solution restricted to $w_{i}$. A spectral problem is then solved in the snapshot space to extract the dominant modes, which are the offline basis functions and the resulting reduced space is called the offline space.

4.3.1. Snapshot space

The first choice of $\boldsymbol{V}^{i}_{\textup{snap}}$ is the restriction of the conforming space $\boldsymbol{V}_{h}$ in $w_{i}$, and the resulting basis functions are called spectral basis functions. Note that $\boldsymbol{V}^{i}_{\textup{snap}}$ contains all possible fine-scale functions defined on $w_{i}$.

The second choice of $\boldsymbol{V}^{i}_{\textup{snap}}$ is the set of all $\kappa$-harmonic extensions, and the resulting basis functions are called harmonic basis functions. More specifically, we denote the fine-grid function $\delta_{j}^{h}(\boldsymbol{x}_{k}):=\delta_{jk}$ for $\boldsymbol{x}_{k}\in J_{h}(w_{i})$, where $J_{h}(w_{i})$ denotes the set of all nodes of the fine mesh $\mathcal{T}_{h}$ belonging to $\partial w_{i}$. The cardinality of $J_{h}(w_{i})$ is denoted by $J_{i}$. At the $n$-th Picard iteration, for each $j=1,\cdots,J_{i}$, the snapshot function $\boldsymbol{\psi}^{i}_{j}$ is defined to be the solution to the following system

[TABLE]

For each coarse region $w_{i}$, the corresponding local snapshot space $\boldsymbol{V}^{i}_{\textup{snap}}$ is defined as $\boldsymbol{V}^{i}_{\textup{snap}}:=\textup{span}\{\boldsymbol{\psi}^{i}_{j}:j=1,\cdots,J_{i}\}$. Then, one may define the global snapshot space $\boldsymbol{V}_{\textup{snap}}$ as $\boldsymbol{V}_{\textup{snap}}:=\oplus_{i=1}^{N_{v}}\boldsymbol{V}^{i}_{\textup{snap}}$.

For simplicity, in this paper, we will use the first choice of $\boldsymbol{V}^{i}_{\textup{snap}}$ consisting of the spectral basis functions. We also use this choice in our numerical simulations, and still use $J_{i}$ to denote the number of basis functions of $\boldsymbol{V}^{i}_{\textup{snap}}$.

4.3.2. Offline multiscale basis construction

To obtain the offline basis functions, we need to perform a space reduction by a spectral problem. The analysis in [9] motivates the following construction. The spectral problem that is needed for the purpose of space reduction is as follows: find $(\boldsymbol{\psi}^{i}_{j},\lambda^{i}_{j})\in\boldsymbol{V}^{i}_{\textup{snap}}\times\mathbb{R},j=1,2,\cdots,J_{i}$ such that

[TABLE]

where the weighted function $\tilde{\kappa}(\boldsymbol{x},|\boldsymbol{Du}^{n}_{\textup{ms}}|)$ is defined by (see [9])

[TABLE]

and $\{\chi_{i}\}$ is a set of standard multiscale finite element basis functions, which is a partition of unity, for the coarse node $\boldsymbol{x}_{i}$ (that is, with linear boundary conditions for cell problems). Specifically, $\forall K\in w_{i}\,,$ $\chi_{i}$ is defined via

[TABLE]

where $\Phi_{i}$ is linear and continuous on $\partial K$. That is, the multiscale partition of unity is defined as $\chi_{i}=\tilde{\Phi}_{i}=(\boldsymbol{\zeta}_{i})_{1}$.

After arranging the eigenvalues $\lambda_{j}^{i}$, $j=1,2,\cdots,J_{i}$ from (4.4) in ascending order, we choose the first $l_{i}$ eigenfunctions from (4.4), and denote them by $\boldsymbol{\Psi}_{1}^{\textup{off}},\cdots,\boldsymbol{\Psi}_{l_{i}}^{\textup{off}}$. Using these eigenfunctions, we can establish the corresponding eigenvectors in the space of snapshots via the formulation

[TABLE]

for $k=1,\cdots,l_{i}$, where $(\Psi_{k}^{i,\textup{off}})_{j}$ denotes the $j$-th component of the vector $\boldsymbol{\Psi}_{k}^{i,\textup{off}}$. At the final step, the offline basis functions for the coarse neighborhood $w_{i}$ is defined by $\boldsymbol{\psi}_{k}^{i,\textup{off}}=\chi_{i}\boldsymbol{\phi}_{k}^{i,\textup{off}}$, where $\{\chi_{i}\}$ is a set of standard multiscale finite element basis functions, which is a partition of unity, for the coarse neighborhood $w_{i}$. We now define the local auxiliary offline multiscale space $\boldsymbol{V}_{\textup{off}}^{i}$ as the linear span of all $\boldsymbol{\psi}_{k}^{i,\textup{off}},k=1,2,\cdots,l_{i}$.

Using the notation in (4.3), one can take $\boldsymbol{V}^{m}_{\textup{ms}}$ as $\boldsymbol{V}^{m}_{\textup{off}}:=\textup{span}\{\boldsymbol{\psi}_{k}^{i,\textup{off}}\,|\,1\leq i\leq N_{v}\,,1\leq k\leq l_{i}\}$. We refer to [17] for the convergence of the method within the current Picard iteration.

5. Residual based online adaptive GMsFEM

As we mentioned in the previous Sections, some online basis functions are required to obtain a coarse representation of the fine-scale solution and give a fast convergence of the corresponding adaptive enrichment algorithm. In [9], such online adaptivity is proposed and mathematically analyzed. More specifically, at the current $n$-th Picard iteration, when the local residual related to some coarse neighborhood $w_{i}$ is large (see Subsection 5.2), one may construct a new basis function $\boldsymbol{\phi}_{i}\in\boldsymbol{V}_{i}=\boldsymbol{H}^{1}_{0}(w_{i})\cap\boldsymbol{V}_{h}$ (with the equipped norm $\|\boldsymbol{v}\|^{2}_{\boldsymbol{V}_{i}}=\int_{w_{i}}\kappa(\boldsymbol{x},|\boldsymbol{Du}^{n}_{\textup{ms}}|)|\boldsymbol{Dv}|^{2}\,dx$), and add it to the multiscale basis functions space. It is further shown that if the offline space contains sufficient information in the form of offline basis functions, then the online basis construction results in an efficient approximation of the fine-scale solution $\boldsymbol{u}_{h}$.

At the considering $n$-th Picard iteration, we use the index $m\text{ }(\geq 1)$ to stand for the adaptive enrichment level. Thus, $\boldsymbol{V}^{m}_{\textup{ms}}$ denotes the corresponding GMsFEM space, and $\boldsymbol{u}^{m}_{\textup{ms}}$ represents the corresponding solution obtained in (4.3). The sequence of functions $\{\boldsymbol{u}^{m}_{\textup{ms}}\}_{m\geq 1}$ will converge to the fine-scale solution $\boldsymbol{u}_{h}$. In this Section, we remark that our space $\boldsymbol{V}^{m}_{\textup{ms}}$ can consist of both offline and online basis functions. We will establish an approach for obtaining the space $\boldsymbol{V}^{m+1}_{\textup{ms}}$ from $\boldsymbol{V}^{m}_{\textup{ms}}$.

In the following paragraphs, based on [9], we present a framework for the construction of online basis functions. By online basis functions, we mean basis functions that are computed during the adaptively iterative process, which contrasts with offline basis functions that are computed before the iterative process. The online basis functions are computed based on some local residuals for the current multiscale solution, that is, the function $\boldsymbol{u}^{m}_{\textup{ms}}$. Hence, we realize that some offline basis functions are crucial for the computations of online basis functions. We will also obtain the sufficiently large number of offline basis functions, which are required in order to get a rapidly converging sequence of solutions.

At the current $n$-th Picard iteration, we are given a coarse neighborhood $w_{i}$ and an inner adaptive iteration $m$-th with the approximation space $\boldsymbol{V}^{m}_{\textup{ms}}$. Recall that the GMsFEM solution $\boldsymbol{u}^{m}_{\textup{ms}}\in\boldsymbol{V}^{m}_{\textup{ms}}$ can be obtained by solving (4.3):

[TABLE]

Suppose that we need to add a basis function $\boldsymbol{\phi}_{i}\in\boldsymbol{V}_{i}$ on the $i$-th coarse neighborhood $w_{i}$. Initially, one can set $\boldsymbol{V}^{0}_{\textup{ms}}=\boldsymbol{V}_{\textup{off}}$. Let $\boldsymbol{V}^{m+1}_{\textup{ms}}=\boldsymbol{V}^{m}_{\textup{ms}}+\textup{span}\{\boldsymbol{\phi}_{i}\}$ be the new approximation space, with $\boldsymbol{u}^{m+1}_{\textup{ms}}\in\boldsymbol{V}^{m+1}_{\textup{ms}}$ being the corresponding GMsFEM solution from (4.3). Let

[TABLE]

The argument from [9] deduces that the new online basis function $\boldsymbol{\phi}_{i}\in\boldsymbol{V}_{i}$ is the solution of

[TABLE]

and $\|\boldsymbol{\phi}_{i}\|_{\boldsymbol{V}_{i}}=\|R_{i}\|_{\boldsymbol{V}^{*}_{i}}\,.$ This means the residual norm $\|R_{i}\|_{\boldsymbol{V}^{*}_{i}}$ (using $\boldsymbol{H}^{-1}(w_{i})$ norm) gives a measure on the quantity of reduction in energy error. Also, it holds that

[TABLE]

for $\alpha=a_{n}(\boldsymbol{u}_{h}-\boldsymbol{u}^{m}_{\textup{ms}},\boldsymbol{\phi}_{i})$. This algorithm is called the online adaptive GMsFEM because only online basis functions are used.

The convergence of this algorithm is discussed in [9], within the current Picard iteration.

5.1. Error estimation in a Picard iteration

At the current $n$-th Picard iteration, we show a sufficient condition for reduction in the error. Let $I_{p}\subset\{1,2,\cdots,N_{v}\}$ be the index set over coarse neighborhoods $w_{i}$ $(i\in I_{p})$, which are non-overlapping. For each $i\in I_{p}$, we define the online basis functions $\boldsymbol{\phi}_{i}\in\boldsymbol{V}_{i}$ by the solution to the equation (5.1). Set $\boldsymbol{V}^{m+1}_{\textup{ms}}=\boldsymbol{V}^{m}_{\textup{ms}}\oplus\textup{span}\{\boldsymbol{\phi}_{i}:i\in I_{p}\}\,.$ Let $r_{i}=\|R_{i}\|_{\boldsymbol{V}^{*}_{i}}$. Let $\Lambda_{p}=\displaystyle\min_{i\in I_{p}}\lambda^{i}_{l_{i}+1}$, where $\lambda^{i}_{l_{i}+1}$ is the $(l_{i}+1)$-th eigenvalue ($k=l_{i}+1$) from (4.4) in the coarse region $w_{i}$.

From ([9], equation (15)), we obtain the following estimate for the energy norm error reduction

[TABLE]

where $C$ is a uniform constant, independent of the contrast in $\kappa(\boldsymbol{x},|\boldsymbol{Du}^{n}_{\textup{ms}}|)$. Toward error decreasing, for a large enough $\Lambda_{p}$, there is a need to pick sufficiently many offline basis functions, which means the online error reduction property ([9]) as below.

Definition 5.1**.**

We say that $\boldsymbol{V}_{\textup{off}}$ satisfies Online Error Reduction Property (ONERP) if

[TABLE]

for some $\theta_{0}>\gamma>0$, where $\gamma$ is independent of physical parameters such as contrast.

Theoretically, the ONERP is required in order to obtain fast and robust convergence, which is independent of the contrast in the permeability, for general quantities of interest.

5.2. Online adaptive algorithm

Set $m=0$. Pick a parameter $\theta\in(0,1]$ and denote $\boldsymbol{V}^{m}_{\textup{ms}}=\boldsymbol{V}^{0}_{\textup{ms}}=\boldsymbol{V}_{\textup{off}}$. Choose a small tolerance $\boldsymbol{tol}\in\mathbb{R}_{+}$. For each $m\in\mathbb{N}$, assume that $\boldsymbol{V}^{m}_{\textup{ms}}$ is given. Go to the following Step 1.

Step 1: Solve for $\boldsymbol{u}^{m}_{\textup{ms}}\in\boldsymbol{V}^{m}_{\textup{ms}}$ from the equation (4.3).

Step 2: For each $i=1,\cdots,N_{v}$, compute the residual $r_{i}$ for the coarse region $w_{i}$. Assume that we obtain

[TABLE]

Step 3: Pick the smallest integer $k_{p}$ such that

[TABLE]

Now, for $i=1,\cdots,k_{p}$, add basis functions $\boldsymbol{\phi}_{i}$ (by solving (5.1)) to the space $\boldsymbol{V}^{m}_{\textup{ms}}$. The new multiscale basis functions space is defined as $\boldsymbol{V}^{m+1}_{\textup{ms}}$. That is

[TABLE]

Step 4: If $\displaystyle\sum_{i=1}^{N_{v}}r_{i}^{2}\leq\boldsymbol{tol}$ or the dimension of $\boldsymbol{V}^{m+1}_{\textup{ms}}$ is sufficiently large, then stop. Otherwise, set $m\leftarrow m+1$ and go back to Step 1.

Remark 5.2**.**

If $\theta=1$, then the adaptive enrichment is said to be uniform.

In practice, we do not update the basis function space $\boldsymbol{V}^{n+1}_{\textup{ms}}$ ($n\geq 0$) at every Picard iteration step $(n+1)$-th. We consider a simple adaptive strategy to update the basis space. More specifically, at each Picard iteration $(n+1)$-th, after updating the multiscale solution $\boldsymbol{u}^{n+1}_{\textup{ms}}$ in $\boldsymbol{V}^{n}_{\textup{ms}}$, we compute the coefficient $\kappa(\boldsymbol{x})=\dfrac{1}{1-\beta(\boldsymbol{x})|(\boldsymbol{Du}^{n+1}_{\textup{ms}})|}$, and the relative $\boldsymbol{L}^{2}$ change of this updated coefficient and the coefficient corresponding to the last step that the basis was updated. If the change is larger than a predefined tolerance, then we recompute the offline and online basis functions. In particular, $\delta=0$ implies that we update the basis functions in every Picard iteration, while $\delta=\infty$ implies that we do not update the basis functions.

5.3. GMsFEM for nonlinear elasticity

We summarize the major steps of using the GMsFEM to solve problem (2.8-2.9): pick a basis update tolerance value $\delta\in\mathbb{R}_{+}$ and Picard iteration termination tolerance value $\delta_{0}\in\mathbb{R}_{+}$ (where $\delta_{0}$ and $\delta$ to be specified in Section 6). We also take an initial guess of $\boldsymbol{u}^{\text{old}}_{\textup{ms}}$, and compute $\kappa^{\text{old}}(\boldsymbol{x})=\dfrac{1}{1-\beta(\boldsymbol{x})|(\boldsymbol{Du}^{\text{old}}_{\textup{ms}})|}$ and the multiscale space $\boldsymbol{V}^{\text{old}}_{\textup{ms}}$, then we repeat following steps:

Step 1: Solve for $\boldsymbol{u}^{\text{new}}_{\textup{ms}}\in\boldsymbol{V}^{\text{old}}_{\textup{ms}}$ from the following equation (as (4.3)):

[TABLE]

If $\dfrac{\|\boldsymbol{u}^{\text{new}}_{\textup{ms}}-\boldsymbol{u}^{\text{old}}_{\textup{ms}}\|_{\boldsymbol{V}_{h}}}{\|\boldsymbol{u}^{\text{old}}_{\textup{ms}}\|_{\boldsymbol{V}_{h}}}>\delta_{0}$, let $\boldsymbol{u}^{\text{old}}_{\textup{ms}}=\boldsymbol{u}^{\text{new}}_{\textup{ms}}$ and go to Step 2.

Step 2: Calculate $\kappa^{\text{new}}(\boldsymbol{x})=\dfrac{1}{1-\beta(\boldsymbol{x})|(\boldsymbol{Du}^{\text{new}}_{\textup{ms}})|}$. If $\dfrac{||\kappa^{\text{old}}(\boldsymbol{x})-\kappa^{\text{new}}(\boldsymbol{x})||_{\boldsymbol{L}^{2}(\Omega)}}{||\kappa^{\text{old}}(\boldsymbol{x})||_{\boldsymbol{L}^{2}(\Omega)}}>\delta$, compute the new basis functions space $\boldsymbol{V}^{\text{new}}_{\textup{ms}}$, let $\boldsymbol{V}^{\text{old}}_{\textup{ms}}$=$\boldsymbol{V}^{\text{new}}_{\textup{ms}}$ and $\kappa^{\text{old}}(\boldsymbol{x})=\kappa^{\text{new}}(\boldsymbol{x})$. Then go to Step 1.

6. Numerical examples

In this section, we will present several test cases to show the performance of our GMsFEM. In our simulations, we consider two choices of $\beta(\boldsymbol{x})$, which are shown in Figure 1. In these two test cases, the blue region represents $\beta(\boldsymbol{x})=1$ and the red regions (channels) represent $\beta(\boldsymbol{x})=10^{-4}$. In addition, the coefficient $\beta(\boldsymbol{x})$ is defined on a $200\times 200$ fine grid. For the coarse grid size, we choose $H=1/20$. We take the source term $f=(\sqrt{x^{2}+y^{2}+1},\sqrt{x^{2}+y^{2}+1})$, $\delta_{0}=10^{-7}$, and $\delta$ as in the following Tables 1-8.

At the $(n+1)$-th Picard iteration, to quantify the quality of our multiscale solutions, we use the following relative $\boldsymbol{L}^{2}$ norm error and weighted $\boldsymbol{H}^{1}$ norm error:

[TABLE]

where the reference solution $\boldsymbol{u_{h}}$ is computed on the fine grid, and the bilinear form $a_{n}$ is defined in (4.2).

We will first present the results using only offline basis functions. We tested the performance on the use of various number of basis functions, which is denoted by $Nb$. The error history with different number of basis functions and update tolerance $\delta$ is shown in Table 1 for the Model 1 and in Table 3 for the Model 2. First of all, we observe that by updating the basis functions more frequently, one can obtain better approximate solutions. On the other hand, we observe that the errors stay around a fixed level. This motivates us to use the online basis functions in order to improve the approximation quality.

In our next test, we consider the addition of online basis functions. In this case, we will construct both offline and online basis functions. More precisely, we will first find the new offline basis functions by the updated coefficient. After that, we solve the PDE using the new offline basis functions. Then, using the residual, we construct online basis functions. The error history is shown in Table 2 for the Model 1 and in Table 4 for the Model 2. In these tables, we use $1+2$ in the $Nb$ column to represent the use of $1$ offline basis function and $2$ online basis functions. From these tables, we observe that updating the basis functions will produce more significant improvement in the approximate solutions. Also, we observe that the use of online basis functions is able to produce more accurate solutions.

In our second set of simulations, we repeat the above tests but with $\beta(\boldsymbol{x})=10^{4}$ in the red regions (channels). We take $f=10^{-4}(\sqrt{x^{2}+y^{2}+1},\sqrt{x^{2}+y^{2}+1})$ to ensure the convergence of the Picard iteration procedure. From the results in Tables 5-8, we observe similar results to the first set of simulations. These results indicate that our method is robust with respect to the contrast in the coefficient, and is able to give accurate approximate solution with a few local basis functions per each coarse neighborhood.

7. Conclusions

In this paper, we propose a GMsFEM framework for a strain-limiting nonlinear elasticity model. The main idea here is the combination of Picard iteration procedure (for linearization) and the two types (offline and online) of basis functions within GMsFEM (for handling the multiple scales and high contrast of materials). This means that at each Picard iteration, the problem is linear; and in each coarse neighborhood, we use offline multiscale basis functions (whose number is determined by a local error indicator) or combine them with residual based online adaptive basis functions (which are only added in regions with large errors, and can capture global features of the solution).

Our numerical results show that the combination of offline and online basis functions is able to give accurate solutions, accelerate the convergence, and reduce computational cost with only a small number of Picard iterations as well as basis functions per each coarse region. In a future contribution, we will address the development of this GMsFEM using the constraint energy minimization approach [10, 11, 19], for nonlinear problems.

Acknowledgments. EC’s work is partially supported by Hong Kong RGC General Research Fund (Projects: 14317516 and 14304217) and CUHK Direct Grant for Research 2017-18.