Growing Sandpile Problem with Local and Non-Local Boundaries Conditions

Urbain Traoré

Laboratoire de Mathématiques et d’Informatique (LAMI), UFR, Sciences Exactes et Appliquées, Université Ouaga I Pr Joseph.

DOI: 10.4236/jamp.2021.910153                

Abstract

In this paper, we study a class of Prigozhin equation for growing sandpile problem subject to local and a non-local boundary condition. The problem is a generalized model for a growing sandpile problem with Neumann boundary condition (see [1]). By the semi-group theory, we prove the existence and uniqueness of the solution for the model and thanks to a duality method we do the numerical analysis of the problem. We finish our work by doing numerical simulations to validate our theoretical results.

Keywords

Sandpile, Collapsing, Avalanche, Nonlinear Semi-Group, Nonlocal Boundary Conditions, Euler Discretization in Time, Time Dependent Gradient Constraints

1. Introduction

The mathematical modeling of the dynamics of granular materials is a fascinating subject as evidenced by the numerous studies which treat this subject in the literature. In this work, we are interested in the formation of a pile of sand under the effect of a vertical source. Regarding this phenomenon, there are mainly two models based on partial differential equations. The first model is that of Hadeler and Kuttler (see [2] [3] [4] ) where the authors model the pile of sand using two layers: a fixed layer above which a mobile layer moves. The second approach to which this work relates is the variational model of Prigozhin (see [5] [6] ) where the surface flow of sand is supposed to exist only at critical slope, that is the maximal admissible slope $\alpha$ for any stationary configuration of sand (here for simplicity it is assumed $\alpha =1$ ).

In the particular case of the Dirichlet boundary condition; for that model there are a sufficiently developed theory and efficient numerical schemes which use duality arguments (see [7] [8] [9] ). However, very few results exist in the literature on the model of Prighozin in the case where one considers other types of boundary conditions other than the Dirichlet boundary condition, since it raises many difficult questions: how to define the notion of the solution for the model and how to approach the problem numerically. We have recently partially provided solutions to these questions where we study the Prigozhin model with Neumann and Robin boundaries conditions (see [10] ). These types of conditions are interpreted by the presence of an obstacle for example a wall at the boundary and which prevents the sand escaping. In this work, we look at the model with a non-local boundary condition which generalizes the Neumann and Robin boundary condition:

$\begin{array}{l}{u}_t-\nabla \left(m\nabla u\right)=f\text{in}\left(\mathrm{0,}T\right)\times \Omega \\ \left|\nabla u\right|\le \mathrm{1,}m\ge \mathrm{0,}m\left(\left|\nabla u\right|-1\right)=0\text{in}\left(\mathrm{0,}T\right)\times \Omega \\ u=0\text{on}\left[\mathrm{0,}T\right]\times {\Gamma }_0\end{array}$

$\begin{array}{l}m\frac{\partial u}{\partial \nu }=g\text{on}\left(\mathrm{0,}T\right)\times {\Gamma }_1\\ {\displaystyle {\int }_{{\Gamma }_2}}m\frac{\partial u}{\partial \nu }=a\left(t\right)\left(\text{given}\text{function}\text{of}t\right)\text{on}\left[\mathrm{0,}T\right]\times {\Gamma }_2\\ u=c\left(t\right)\left(\text{unknown}\text{constant}\text{of}t\right)\text{on}\left(\mathrm{0,}T\right)\times {\Gamma }_2\\ u\left(0\right)=u_0\mathrm{,}\end{array}$ (1)

with $\Omega \subset {ℝ}^N$ a bounded open domain with a boundary $\Gamma ={\Gamma }_0\cup {\Gamma }_1\cup {\Gamma }_2$ such that ${\Gamma }_0\cap {\Gamma }_1\cap {\Gamma }_2=\varnothing$. The solution u is the height of the surface andf is the source, $g\in L^2\left({\Gamma }_1\right)$ and a is a function defined on $ℝ$.

In this work, we consider a mixed boundary condition: Dirichlet, Neumann and a non-local boundaries condition:

· on ${\Gamma }_0$, we apply homogeneous Dirichlet boundary condition. In other words, we assume that at this boundary, the sand falls down, that is

$u=0\text{on}\left(\mathrm{0,}T\right)\times {\Gamma }_0\mathrm{<mo>\; ;\; </mo>}$

· on ${\Gamma }_1$, there is a wall that prevents the sand escaping:

$m\frac{\partial u}{\partial \nu }=g\text{on}\left(\mathrm{0,}T\right)\times {\Gamma }_1\mathrm{<mo>\; ;\; </mo>}$

· we control the outgoing flow through the boundary ${\Gamma }_2$, which results in the following non-local condition:

${\displaystyle {\int }_{{\Gamma }_2}}m\frac{\partial u}{\partial \nu }\text{d}s=A\left(t\right)\text{on}\left[\mathrm{0,}T\right]\times {\Gamma }_2\mathrm{<mo>\; ;\; </mo>}$

Beside the mathematical interest of non-local conditions, it seems that this type of boundary condition is also encountered in other physical applications. For example, in petroleum engineering model for well modelling in a 3D stratified petroleum reservoir with arbitrary geometry; this kind of boundary condition also arises in petroleum engineering, in the simulation of wells performance, since a nonlinear relation exists between the performance pressure tangential gradient and the fluid velocity along the well (see [11] [12] for details). Another application of this type of the boundary condition is in the study of the heat conduction within linear thermoelasticity (see [13] [14] [15] ), and for the reaction-diffusion equation (see [16] [17] ). One could find other applications in other fields of physics in the papers [18] [19] [20].

The present work is organized as follows. In the next section, we use the nonlinear semi-group theory (see [21] ) to get the existence and uniqueness of a variational solution of (1) and the convergence of the approximate Euler discretization in time solutions to problem (1). In Section 3, we show how to compute the solution of Euler implicit time discretization of (1) using duality argument and in the last section, some results of numerical results are given.

2. Global Existence

Given $\epsilon >0$, we say that ${\left(t_i,f_i\right)}_{i=1,\cdots ,n}$ is an $\epsilon$ -discretization for problem (1) if $0=t_0<t_1<\cdots <t_{n-1}<T=t_n$ with $t_i-t_{i-1}\le \epsilon$, $f_1\mathrm{,}\cdots \mathrm{,}{f}_n\in L^2\left(\Omega \right)$, $g_1\mathrm{,}\cdots \mathrm{,}{g}_n\in L^2\left({\Gamma }_1\right)$ such that

${\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle {\int }_{t_{i-1}}^{t_i}}{\Vert f\left(t\right)-f_i\Vert }_{L^2\left(\Omega \right)}\le \epsilon ;{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle {\int }_{t_{i-1}}^{t_i}}{\Vert g\left(t\right)-g_i\Vert }_{L^2\left({\Gamma }_1\right)}\le \epsilon$

Definition 1. For any $\epsilon >0$, $\left(u_{\epsilon }\right)$ is an $\epsilon$ -approximate solution of (1), if there exists ${\left(t_i\mathrm{,}{f}_i\mathrm{,}{g}_i\right)}_{i=1,\cdots ,n}$ an $\epsilon$ -discretization of problem (1) such that

$u_{\epsilon }=(\begin{array}{ll}{u}_0\hfill & \text{for}t\in \left(\mathrm{0,}{t}_1\right]\hfill \\ u_i\hfill & \text{for}t\in \left(t_{i-1}\mathrm{,}{t}_i\right]\mathrm{,}i=\mathrm{1,}\cdots \mathrm{,}n\hfill \end{array}$ (2)

and $u_i$ solves the following Euler implicit time discretization of (1):

$\begin{array}{l}{u}_i-\epsilon \nabla \left(m_i\nabla u_i\right)=\epsilon f_i+u_{i-1}\text{in}\Omega \\ \left|\nabla u_i\right|\le \mathrm{1,}{m}_i\ge \mathrm{0,}{m}_i\left(\left|\nabla u_i\right|-1\right)=0\text{in}\Omega \\ u_i=0\text{on}{\Gamma }_0\\ m_i\frac{\partial u_i}{\partial \nu }=g_i\text{on}{\Gamma }_1\\ {\displaystyle {\int }_{{\Gamma }_2}}m\frac{\partial u_i}{\partial \nu }=a_i\text{on}{\Gamma }_2\\ u_i=c_i\left(\text{unknown}\text{constant}\right)\text{on}{\Gamma }_2\end{array}$ (3)

The generic problem associated with (3) is given by

$\begin{array}{l}v-\nabla \left(m\nabla v\right)=F\text{in}\Omega \\ \left|\nabla v\right|\le \mathrm{1,}m\ge \mathrm{0,}m\left(\left|\nabla v\right|-1\right)=0\text{in}\Omega \\ v=0\text{on}{\Gamma }_0\\ m\frac{\partial u}{\partial \nu }=G\\ {\displaystyle {\int }_{{\Gamma }_2}}m\frac{\partial v}{\partial \nu }=A\left(\text{given}\text{constant}\right)\text{on}{\Gamma }_2\\ v=C\left(\text{unknown}\text{constant}\right)\text{on}{\Gamma }_2\end{array}$ (4)

with $F\in L^2\left(\Omega \right)$, $G\in L^2\left({\Gamma }_1\right)$.

For convenience, we note the $L^2\left(\Omega \right)$ scalar product by $\left(\cdot \mathrm{,}\cdot \right)$ and the euclidean norm on ${ℝ}^N$ by $\Vert \cdot \Vert$. We introduce the convex set $\mathcal{K}$ :

$\mathcal{K}=\left\{z\in W^{\mathrm{1,}\infty }\left(\Omega \right)\cap H_D^1\left(\Omega \right)\mathrm{<mo>\; :\; </mo>}\left|\nabla z\left(x\right)\right|\le 1\text{a}\text{.e}x\in \Omega \right\}$

where

$H_D^1\left(\Omega \right)=\left\{z\in H^1\left(\Omega \right)\mathrm{<mo>\; :\; </mo>}z=0\text{on}{\Gamma }_{0}z\equiv \text{constant}\text{on}{\Gamma }_2\right\}\mathrm{.}$

We also introduce the function ${\Pi }_{\mathcal{K}}$ defined on $L^2(\; \Omega \; )$

${\Pi }_{\mathcal{K}}\left(z\right)=(\begin{array}{ll}-{\displaystyle {\int }_{{\Gamma }_1}}Gz\text{d}s-{Az|}_{{\Gamma }_2}\hfill & \text{if}z\in \mathcal{K}\mathrm{,}\hfill \\ +\infty \hfill & \text{otherwise}\mathrm{,}\hfill \end{array}$

The sub-differential of ${\Pi }_{\mathcal{K}}$ is defined in $L^2\left(\Omega \right)$ by

$\partial {\Pi }_K\left(v\right)=\left\{w\in \mathcal{K}\mathrm{<mo>\; :\; </mo>}\forall z\in K\mathrm{,}{\Pi }_{\mathcal{K}}\left(z\right)\ge {\Pi }_{\mathcal{K}}\left(v\right)+\left(w\mathrm{,}z-v\right)\right\}\mathrm{.}$

Remark. In order to define our notions of variational solutions for the above problems, we make the following observations:

(1) Assuming that $v\in \mathcal{K}$ is a solution of (4) in the following sense: there exists a measurable function m satisfying the properties $m\nabla v\in {\left(L^1\left(\Omega \right)\right)}^d$, $m\left(\left|\nabla v\right|-1\right)=0$ a.e. in $\Omega$ and

${\displaystyle {\int }_{\Omega }}vz\text{d}x+{\displaystyle {\int }_{\Omega }}m\nabla v\cdot \nabla z\text{d}x={\displaystyle {\int }_{\Omega }}Fz\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}Gz\text{d}s+{Az|}_{{\Gamma }_2}\forall z\in H_D^1\left(\Omega \right)\mathrm{<mo>\; ;\; </mo>}$ (5)

then, v is also a solution of the following optimization problem

$\underset{z\in \mathcal{K}}{\max }\left\{{\displaystyle {\int }_{\Omega }}\left(f-v\right)z\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}Gz\text{d}s+{Az|}_{{\Gamma }_2}\right\}\mathrm{.}$ (6)

Indeed, taking $z\in \mathcal{K}$ as test function in (5), we get

${\displaystyle {\int }_{\Omega }}\left(F-v\right)z\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}gz\text{d}s+{Az|}_{{\Gamma }_2}={\displaystyle {\int }_{\Omega }}m\nabla v\cdot \nabla z\text{d}x\mathrm{.}$ (7)

Thus, using the fact that $\left|\nabla z\right|\le 1$ a.e. $x\in \Omega$, we get

${\displaystyle {\int }_{\Omega }}\left(F-v\right)z\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}Gz\text{d}s+{Az|}_{{\Gamma }_2}\le {\displaystyle {\int }_{\Omega }}m\left|\nabla v\right|\text{d}x$ (8)

and taking $z=v$ in (7), we obtain

$-{\displaystyle {\int }_{\Omega }}\left(F-v\right)v\text{d}x-{\displaystyle {\int }_{{\Gamma }_1}}Gz\text{d}s-{Az|}_{{\Gamma }_2}=-{\displaystyle {\int }_{\Omega }}m{\left|\nabla v\right|}^2\text{d}x\mathrm{.}$ (9)

So, adding (8) and (9) and using the fact that $m\left(\left|\nabla v\right|-1\right)=0$ a.e. in $\Omega$, we obtain

${\displaystyle {\int }_{\Omega }}\left(F-v\right)z\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}Gz\text{d}s+{Az|}_{{\Gamma }_2}-\left({\displaystyle {\int }_{\Omega }}\left(F-v\right)v\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}Gv\text{d}s+{Av|}_{{\Gamma }_2}\right)\le 0.$ (10)

(2) It is not clear that there exists $\left(m\mathrm{,}v\right)$ such that (5) is fulfilled, however we can find $v\in \mathcal{K}$, solution of the problem (6), which is a consequence of the following result.

Lemma 2. $\partial {\Pi }_{\mathcal{K}}$ is a maximal monotone graph in $L^2\left(\Omega \right)$.

Proof. The proof is the same as in our previous paper. $\square$

We are now in a position to define our notion of solution to problems (4) and (1).

Definition 3 (see [21] ). For given $F\in L^2\left(\Omega \right)$, $G\in L^2\left({\Gamma }_1\right)$, we say that v is a variational solution of (4) if $v\in \mathcal{K}$ and

${\displaystyle {\int }_{\Omega }}\left(F-v\right)\left(z-v\right)\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}g\left(z-v\right)\text{d}s+A\left({z|}_{{\Gamma }_2}-{v|}_{{\Gamma }_2}\right)\le \mathrm{0,}\text{for all }z\in \mathcal{K}\mathrm{.}$ (11)

Definition 4 (see [21] ). Given $f\in L^{\infty }\left(\mathrm{0,}T\mathrm{<mo>\; ;\; </mo>}{L}^2\left(\Omega \right)\right)$, and $u_0\in \mathcal{K}$, a variational solution (resp. $\epsilon$ -approximate solution) of (1) is a function $u\in W^{\mathrm{1,1}}\left(\mathrm{0,}T\mathrm{<mo>\; ;\; </mo>}{L}^2\left(\Omega \right)\right)$ satisfying for any $t\in \left(\mathrm{0,}T\right)$, $u\left(t\right)\in \mathcal{K}$ and

${\displaystyle {\int }_{\Omega }}\left(f-u_t\left(t\right)\right)\left(z-u\left(t\right)\right)\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}g\left(z-u\left(t\right)\right)\text{d}s+A\left({z|}_{{\Gamma }_2}-{u\left(t\right)|}_{{\Gamma }_2}\right)\le \mathrm{0,}$

for any $z\in \mathcal{K}$ (resp. $u_{\epsilon }$ given by (2) with $u_i$ a variational solution of (3)).

Since $\partial {\Pi }_{\mathcal{K}}$ is a maximal monotone graph in $L^2\left(\Omega \right)$, then, thanks to [21] for any $F\in L^2\left(\Omega \right)$ there exists a unique solution v of

$v+\partial {\Pi }_{\mathcal{K}}\left(v\right)\ni F\mathrm{,}$ (12)

which is equivalent to saying that (4) admits a unique variational solution. Moreover, if $v_i$ is the solution corresponding to $g_i$ for $i=1,2$, then

${\Vert v_1-v_2\Vert }_2\le {\Vert F_1-F_2\Vert }_2\mathrm{.}$ (13)

As in the works of Igbida and al, we get the following result

Theorem 5. Let $u_0\in \mathcal{K}$ , $T>0$ , $f\in L^2\left(\mathrm{0,}T\mathrm{<mo>\; ;\; </mo>}{L}^2\left(\Omega \right)\right)$ . Then,

1. for any $\epsilon >0$ and any $\epsilon$ -discretization of (1), there exists a unique $\epsilon$ -approximate variational solution of (1).

2. There exists $u\in \mathcal{C}\left(\left[\mathrm{0,}T\right)\mathrm{<mo>\; ;\; </mo>}{L}^2\left(\Omega \right)\right)$ such that $u\left(0\right)=u_0$ and as $\epsilon \to 0$,

$u_{\epsilon }\to u\text{in}\mathcal{C}\left(\left[\mathrm{0,}T\right)\mathrm{<mo>\; ;\; </mo>}{L}^2\left(\Omega \right)\right)\mathrm{.}$

3. The function u given by (2) is the unique variational solution of (1). Moreover, if for $i=\mathrm{1,2}$, $u_i$ is a solution corresponding to $f_i$, then

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_{\Omega }}{\left(u_1-u_2\right)}^{+}\le {\displaystyle {\int }_{\Omega }}{\left(f_1-f_2\right)}^{+}\text{in}{\mathcal{D}}^{\prime }\left(\mathrm{0,}T\right)\mathrm{.}$

In particular, if $f\ge 0$, then $u\ge 0$ a.e. in $\Omega$.

3. Numerical Analysis

Following the ideas developed in our previous work we show in this section how to approximate the solution of problem (1). In our previous works, we are introduced a numerical method based on Fenchel duality theory (see [22] ) for numerical approximation of problem (4), in the special case where the homogeneous Dirichlet boundary condition is applied. In order to use the same approach, we introduce the following result.

Lemma 6. The problem (4) is equivalent to: find $v\in \mathcal{K}$ solution of the following of

$\arg \min {}_{z\in \mathcal{K}}\left\{\frac{1}{2}{\displaystyle {\int }_{\Omega }}{\Vert z-g\Vert }_{L^2\left(\Omega \right)}^2+{\Pi }_{\mathcal{K}}\left(z\right)\right\}.$ (14)

Proof. Let $v\in K$ be a solution of (4). Then, we have

${\displaystyle {\int }_{\Omega }}\left(F-v\right)\left(z-v\right)\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}G\left(z-v\right)\text{d}s+A\left({z|}_{{\Gamma }_2}-{v|}_{{\Gamma }_2}\right)\le \mathrm{0,}\forall z\in K\mathrm{.}$ (15)

We use the well know inequality

${\Vert v-F\Vert }_{L^2\left(\Omega \right)}^2-{\Vert z-F\Vert }_{L^2\left(\Omega \right)}^2=2\left(F-v\mathrm{,}z-v\right)-{\Vert v-z\Vert }_{L^2\left(\Omega \right)}^2\mathrm{,}$

to obtain

$\begin{array}{l}\frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left|v-F\right|}^2\text{d}x-\frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left|z-F\right|}^2\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}G\left(z-v\right)\text{d}s+A\left({z|}_{{\Gamma }_2}-{v|}_{{\Gamma }_2}\right)\\ ={\displaystyle {\int }_{\Omega }}\left(F-v\right)\left(z-v\right)\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}G\left(z-v\right)\text{d}s+A\left({z|}_{{\Gamma }_2}-{v|}_{{\Gamma }_2}\right)-\frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left|v-z\right|}^2\text{d}x.\end{array}$

From (15), we deduce that

$\frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left|v-F\right|}^2-\frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left|z-F\right|}^2+{\displaystyle {\int }_{{\Gamma }_1}}G\left(z-v\right)\text{d}s+A\left({z|}_{{\Gamma }_2}-{v|}_{{\Gamma }_2}\right)\le \mathrm{0\; <mo>\; ;\; </mo>}$ (16)

i.e.

$\frac{1}{2}{\Vert v-F\Vert }_{L^2\left(\Omega \right)}^2+{\Pi }_{\mathcal{K}}\left(v\right)\le \frac{1}{2}{\Vert z-F\Vert }_{L^2\left(\Omega \right)}^2+{\Pi }_{\mathcal{K}}\left(z\right)\text{for}\text{any}z\in \mathcal{K}\mathrm{.}$

Which end the first part of the proof. Now assume v is the solution of the minimization problem (14). Let $z_0\in \mathcal{K}$, since $\mathcal{K}$ is convex, then for any $t\in \left[\mathrm{0,1}\right)$, $z=\left(1-t\right)v+t{z}_0\in \mathcal{K}$. Then, we have

$\begin{array}{l}\frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left|v-F\right|}^2\text{d}x-{\displaystyle {\int }_{{\Gamma }_1}}Gv\text{d}s-{Av|}_{{\Gamma }_2}\\ \le \frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left|F-\left(\left(1-t\right)v+t{z}_0\right)\right|}^2\text{d}x-{\displaystyle {\int }_{{\Gamma }_1}}G\left(\left(1-t\right)v+t{z}_0\right)\text{d}s-A{\left(\left(1-t\right)v+t{z}_0\right)|}_{{\Gamma }_2}\\ \le \frac{1}{2}{\Vert F-v-t\left(z_0-v\right)\Vert }_{L^2\left(\Omega \right)}^2-{\displaystyle {\int }_{{\Gamma }_1}}Gv\text{d}s-{Av|}_{{\Gamma }_2}\\ -t{\displaystyle {\int }_{{\Gamma }_1}}G\left(z_0-v\right)\text{d}s-tA\left({z_0|}_{{\Gamma }_2}-{v|}_{{\Gamma }_2}\right)\end{array}$

Which implies that

$\begin{array}{l}\frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left|v-F\right|}^2\text{d}x-{\displaystyle {\int }_{{\Gamma }_1}}Gv\text{d}s-{Av|}_{{\Gamma }_2}\\ \le \frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left|v-F\right|}^2\text{d}x-t\left(F-v,z_0-v\right)+\frac{t^2}{2}{\displaystyle {\int }_{\Omega }}{\left|z_0-v\right|}^2\text{d}x\\ -{\displaystyle {\int }_{{\Gamma }_1}}Gv\text{d}s-{Av|}_{{\Gamma }_2}-t{\displaystyle {\int }_{{\Gamma }_1}}G\left(z_0-v\right)\text{d}s-tA\left({z_0|}_{{\Gamma }_2}-{v|}_{{\Gamma }_2}\right).\end{array}$

Therefore,

$\left(F-v\mathrm{,}{z}_0-v\right)+{\displaystyle {\int }_{{\Gamma }_1}}G\left(z_0-v\right)\text{d}s+A\left({z_0|}_{{\Gamma }_2}-{v|}_{{\Gamma }_2}\right)\le \frac{t}{2}{\displaystyle {\int }_{\Omega }}{\left|z_0-v\right|}^2\text{d}x\mathrm{.}$ (17)

Hence, by letting $t\to 0$ in (17), we obtain

${\displaystyle {\int }_{\Omega }}\left(F-v\right)\left(z_0-v\right)\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}G\left(z_0-v\right)\text{d}s+A\left({z_0|}_{{\Gamma }_2}-{v|}_{{\Gamma }_2}\right)\le \mathrm{0,}$

for any $z_0\in \mathcal{K}$, so v is a solution of (4). $\square$

Inspired by the works of Igbida and al, we consider as a dual problem associated with problem (14), the following optimization problem

$\sup \left\{-G\left(w\right)\mathrm{<mo>\; :\; </mo>}w\in H_{\text{d}iv\mathrm{,}\mathcal{K}}\left(\Omega \right)\right\}\mathrm{,}$ (18)

where

$G\left(w\right)=\frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left(\text{div}\left(w\right)\right)}^2\text{d}x+{\displaystyle {\int }_{\Omega }}g\text{div}\left(w\right)\text{d}x+{\displaystyle {\int }_{\Omega }}\left|w\right|\text{d}x$ (19)

and

$\begin{array}{l}{H}_{\text{div},\mathcal{K}}(\Omega )=\{w\in H_{\text{div}}\left(\Omega \right):{\displaystyle {\int }_{\Omega }}\left[-\text{div}\left(w\right)\right]\xi \text{d}x={\displaystyle {\int }_{\Omega }}w\nabla \xi \text{d}x-{\displaystyle {\int }_{{\Gamma }_1}}G\xi \text{d}s-{A\xi |}_{{\Gamma }_2},\\ \forall \xi \in H_D^1\left(\Omega \right)\},\end{array}$

$H_{\text{div}}\left(\Omega \right)\mathrm{<mo>\; :\; </mo>}=\left\{w\in {\left(L^2\left(\Omega \right)\right)}^N\mathrm{<mo>\; ;\; </mo>}\text{div}\left(w\right)\in L^2\left(\Omega \right)\right\}\mathrm{.}$

In the following, we prove the connection between problems (18) and (14). We also present a method for numerical approximation of the solution of problem (4) by computing $\sup \left\{-G\left(w\right)\mathrm{<mo>\; :\; </mo>}w\in H_{\text{div}\mathrm{,}\mathcal{K}}\left(\Omega \right)\right\}$. We first show the following result.

Lemma 7. Let $F\in L^2\left(\Omega \right)$ , $G\in L^2\left({\Gamma }_1\right)$ , $w\in H_{\text{div}\mathrm{,}\mathcal{K}}\left(\Omega \right)$ and $z\in \mathcal{K}$ , we have

$-G\left(w\right)\le J\left(z\right)\mathrm{,}$

where $J\left(z\right)=\frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left|z-F\right|}^2\text{d}x+{\Pi }_{\mathcal{K}}\left(z\right)$.

Proof. Taking $w\in H_{\text{div}\mathrm{,}\mathcal{K}}\left(\Omega \right)$, $z\in \mathcal{K}$ and using the fact that $\frac{1}{2}{\left(\text{div}\left(w\right)-\left(z-g\right)\right)}^2\ge 0$, we get

$-\frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left(\text{div}\left(w\right)\right)}^2\text{d}x-{\displaystyle {\int }_{\Omega }}\text{div}\left(w\right)g\text{d}x\le \frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left(z-g\right)}^2\text{d}x-{\displaystyle {\int }_{\Omega }}\text{div}\left(w\right)z\text{d}x\mathrm{.}$

Since $w\in H_{\text{div}\mathrm{,}K}\left(\Omega \right)$, we have

$\begin{array}{l}-\frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left(\text{div}\left(w\right)\right)}^2\text{d}x-{\displaystyle {\int }_{\Omega }}\text{div}\left(w\right)g\text{d}x\\ \le \frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left(z-g\right)}^2\text{d}x+{\displaystyle {\int }_{\Omega }}w\cdot \nabla z\text{d}x-{\displaystyle {\int }_{{\Gamma }_1}}Gz\text{d}s-{Az|}_{{\Gamma }_2}\mathrm{.}\end{array}$

Thus,

$-G\left(w\right)\le \frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left(z-g\right)}^2\text{d}x-{\displaystyle {\int }_{{\Gamma }_1}}Gz\text{d}s-{Az|}_{{\Gamma }_2}+{\displaystyle {\int }_{\Omega }}w\cdot \nabla z\text{d}x-{\displaystyle {\int }_{\Omega }}\left|w\right|\text{d}x\mathrm{.}$ (20)

As z belongs to $\mathcal{K}$, it follows that

${\displaystyle {\int }_{\Omega }}w\cdot \nabla z\text{d}x-{\displaystyle {\int }_{\Omega }}\left|w\right|\text{d}x\le 0.$ (21)

Then,

$-G\left(w\right)\le \frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left(z-g\right)}^2\text{d}x-{\displaystyle {\int }_{{\Gamma }_1}}Gz\text{d}s-{Az|}_{{\Gamma }_2}$

The connection between the problem (14) and (18) is given by the following results.

Theorem 8. Let $F\in L^2\left(\Omega \right)$ , $G\in L^2\left({\Gamma }_1\right)$ and v the variational solution of (14). Then, there exists a sequence ${\left(w_{\epsilon }\right)}_{\epsilon >0}$ in $H_{\text{div}\mathrm{,}\mathcal{K}}\left(\Omega \right)$ , $\omega \in {\left(\mathcal{M}\left(\Omega \right)\right)}^s$ such that, as $\epsilon \to 0$ ,

${\displaystyle {\int }_{\Omega }}\left|{\omega }_{\epsilon }\right|\text{d}x\to {\displaystyle {\int }_{\Omega }}v\left(g-v\right)\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}Gv\text{d}s+{Av|}_{{\Gamma }_2}\mathrm{,}$ (22)

$\text{div}\left({\omega }_{\epsilon }\right)\to v-F\in L^2\left(\Omega \right)$ (23)

${\omega }_{\epsilon }\to \omega \in {\left(\mathcal{M}\left(\Omega \right)\right)}^s-{\text{weak}}^{\ast }\mathrm{.}$ (24)

Moreover, we have

$\begin{array}{c}\underset{\epsilon \to 0}{\lim }G\left({\omega }_{\epsilon }\right)=\underset{\omega \in H_{\text{div}\mathrm{,}\mathcal{K}}\left(\Omega \right)}{\inf }G\left(\omega \right)=-\underset{z\in K}{\min }J\left(z\right)\\ =-\left[\frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left|v-F\right|}^2\text{d}x-{\displaystyle {\int }_{{\Gamma }_1}}Gv\text{d}s-{Av|}_{{\Gamma }_2}\right]\mathrm{.}\end{array}$ (25)

Let recall that, $\mathcal{M}\left(\Omega \right)$ denote the space of all bounded Radon measures in $\Omega$.

For the proof of this result, we analyze the following problem

$\begin{array}{l}{v}_{\epsilon }-\nabla w_{\epsilon }=F\text{in}\Omega \\ w_{\epsilon }={\varphi }_{\epsilon }\left(\nabla v_{\epsilon }\right)\text{in}\Omega \\ v_{\epsilon }=0\text{on}{\Gamma }_0\\ \frac{\partial w_{\epsilon }}{\partial \eta }=G\text{on}{\Gamma }_1\mathrm{,}\\ {\displaystyle {\int }_{{\Gamma }_2}}\frac{\partial w_{\epsilon }}{\partial \eta }\text{d}s=A\text{on}{\Gamma }_2\mathrm{,}\\ v_{\epsilon }\equiv \text{constanton}{\Gamma }_2\end{array}$ (26)

where for any $\epsilon >0$, ${\varphi }_{\epsilon }\mathrm{<mo>\; :\; </mo>}{ℝ}^N\to {ℝ}^s$ is given by

${\varphi }_{\epsilon }\left(r\right)=\frac{1}{\epsilon }{\left(\left|r\right|-1\right)}^{+}\frac{r}{\left|r\right|}\text{forall}r\in {ℝ}^N$

and satisfies the following properties.

(i) for any $r_1,r_2\in {ℝ}^N$, $\left({\varphi }_{\epsilon }\left(r_1\right)-{\varphi }_{\epsilon }\left(r_2\right)\right)\cdot \left(r_1-r_2\right)\ge 0$ ;

(ii) there exists ${\epsilon }_0>0$ and $\rho >d$ such that ${\varphi }_{\epsilon }\left(r\right)\cdot r\ge {\left|r\right|}^2$ for any $\left|r\right|>\rho$ and $\epsilon <{\epsilon }_0$ ;

(iii) for any $\epsilon >0$ and $r\in {ℝ}^N$, $\left|{\varphi }_{\epsilon }\left(r\right)\right|\le {\varphi }_{\epsilon }\left(r\right)\cdot r$.

Lemma 9. There exists a unique weak solution ${\left(v_{\epsilon }\right)}_{0<\epsilon <{\epsilon }_0}$ to problem (26) in the sense that $v_{\epsilon }\in H_D^1\left(\Omega \right)$ , $w_{\epsilon }={\varphi }_{\epsilon }\left(\nabla v_{\epsilon }\right)\in {\left(L^2\left(\Omega \right)\right)}^N$ and for all $z\in H_D^1\left(\Omega \right)$ ,

${\displaystyle {\int }_{\Omega }}{v}_{\epsilon }z\text{d}x+{\displaystyle {\int }_{\Omega }}{\varphi }_{\epsilon }\left(\nabla v_{\epsilon }\right)\nabla z\text{d}x={\displaystyle {\int }_{\Omega }}Fz\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}Gz\text{d}s+{Az|}_{{\Gamma }_2}\mathrm{.}$ (27)

Moreover,

(1) ${\left(v_{\epsilon }\right)}_{0<\epsilon <{\epsilon }_0}$ is bounded in $H_D^1\left(\Omega \right)$.

(2) ${\left(w_{\epsilon }\right)}_{0<\epsilon <{\epsilon }_0}$ is bounded in ${\left(L^1\left(\Omega \right)\right)}^N$.

(3) For any Borel set $B\subset \Omega$,

$\lim \underset{\epsilon \to 0}{\inf }{\displaystyle {\int }_B}\left|\nabla v_{\epsilon }\right|\text{d}x\le d\left|B\right|\mathrm{.}$ (28)

Proof. (1) We define the operator $A_{\epsilon }\mathrm{<mo>\; :\; </mo>}{H}_D^1\left(\Omega \right)\to {\left(H_D^1\left(\Omega \right)\right)}^{\ast }$ by

$\langle A_{\epsilon }v\mathrm{,}\xi \rangle ={\displaystyle {\int }_{\Omega }}v\xi \text{d}x+{\displaystyle {\int }_{\Omega }}{\varphi }_{\epsilon }\left(\nabla v\right)\cdot \nabla \xi \text{d}x\mathrm{.}$

It’s not difficult to see that the operator $A_{\epsilon }$ is monotone, coercive, hemicontinous and bounded.

Defining $C\mathrm{<mo>\; :\; </mo>}{H}_D^1\left(\Omega \right)\to ℝ$ defined by

$C\left(z\right)={\displaystyle {\int }_{\Omega }}G\xi \text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}G\xi \text{d}s+{A\xi |}_{{\Gamma }_2}\mathrm{.}$

then thanks to [23], there exists $v_{\epsilon }\in H_D^1\left(\Omega \right)$ such that

$\langle A{v}_{\epsilon }\mathrm{,}\xi \rangle =\langle C\mathrm{,}\xi \rangle \mathrm{,}\text{for}\text{all}\xi \in H_D^1(\; \Omega \; )$

i.e.

${\displaystyle {\int }_{\Omega }}{v}_{\epsilon }\xi \text{d}x+{\displaystyle {\int }_{\Omega }}{\varphi }_{\epsilon }\left(\nabla v_{\epsilon }\right)\cdot \nabla \xi \text{d}x={\displaystyle {\int }_{\Omega }}F\xi \text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}G\xi \text{d}s+{A\xi |}_{{\Gamma }_2}\mathrm{.}$

Which end the proof of the existence. Now, if (26) admits two solutions u and v, then we subtract the two equations obtained by replacing respectively $v_{\epsilon }$ by u and v in (27) to get

${\displaystyle {\int }_{\Omega }}\left(u-v\right)\xi \text{d}x+{\displaystyle {\int }_{\Omega }}\left({\varphi }_{\epsilon }\left(\nabla u\right)-{\varphi }_{\epsilon }\left(\nabla v\right)\right)\nabla \xi \text{d}x=\mathrm{0,}\forall \xi \in H_D^1\left(\Omega \right)\mathrm{.}$ (29)

We take $u-v$ as a test function in (29) to obtain

${\displaystyle {\int }_{\Omega }}{\left(u-v\right)}^2\text{d}x+{\displaystyle {\int }_{\Omega }}\left({\varphi }_{\epsilon }\left(\nabla u\right)-{\varphi }_{\epsilon }\left(\nabla v\right)\right)\cdot \left(\nabla u-\nabla v\right)\text{d}x=0.$ (30)

Thus, using the property (i) of ${\varphi }_{\epsilon }$, we deduce that $u=v$. a.e. $\Omega$

(1) Taking $v_{\epsilon }$ as test function in (27), we have

${\displaystyle {\int }_{\Omega }}{v}_{\epsilon }^2\text{d}x+\frac{1}{\epsilon }{\displaystyle {\int }_{\Omega }}{\left(\left|\nabla v_{\epsilon }\right|-1\right)}^{+}\left|\nabla v_{\epsilon }\right|\text{d}x={\displaystyle {\int }_{\Omega }}F{v}_{\epsilon }\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}G{v}_{\epsilon }\text{d}s+{A{v}_{\epsilon }|}_{{\Gamma }_2}\mathrm{.}$

Since $v_{\epsilon }$ is constant on ${\Gamma }_1$, there exists $C\in ℝ$ such that

$\begin{array}{c}{\displaystyle {\int }_{\Omega }}{v}_{\epsilon }^2\text{d}x\le {\displaystyle {\int }_{\Omega }}F{v}_{\epsilon }\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}G{v}_{\epsilon }\text{d}s+{A{v}_{\epsilon }|}_{{\Gamma }_2}\\ \le {\Vert F\Vert }_{L^2\left(\Omega \right)}{\Vert v_{\epsilon }\Vert }_{L^2\left(\Omega \right)}+{\Vert G\Vert }_{L^2\left({\Gamma }_N\right)}{\Vert v_{\epsilon }\Vert }_{L^2\left({\Gamma }_N\right)}+\left|AC\right|\\ \le {\Vert F\Vert }_{L^2\left(\Omega \right)}{\Vert v_{\epsilon }\Vert }_{H^1\left(\Omega \right)}+{\Vert G\Vert }_{L^2\left({\Gamma }_N\right)}{\Vert v_{\epsilon }\Vert }_{L^2\left({\Gamma }_N\right)}+\left|AC\right|\mathrm{.}\end{array}$

Using the following trace inequality (see [24] )

${\displaystyle {\int }_{\Gamma }}{\left|v_{\epsilon }\right|}^2\text{d}s\le C_1\left({\displaystyle {\int }_{\Omega }}{\left|v_{\epsilon }\right|}^2\text{d}x+{\displaystyle {\int }_{\Omega }}{\left|\nabla v_{\epsilon }\right|}^2\text{d}x\right)\mathrm{,}$

where $C_1$ is a positive constant. One obtains

${\Vert v_{\epsilon }\Vert }_{L^2\left(\Omega \right)}^2\le \left({\Vert F\Vert }_{L^2\left(\Omega \right)}+C_1{\Vert G\Vert }_{L^2\left({\Gamma }_N\right)}\right){\Vert v_{\epsilon }\Vert }_{H^1\left(\Omega \right)}+\left|AC\right|\mathrm{.}$ (31)

On the other hand

$\frac{1}{\epsilon }{\displaystyle {\int }_{\Omega }}{\left(\left|\nabla v_{\epsilon }\right|-1\right)}^{+}\left|\nabla v_{\epsilon }\right|\text{d}x\le \left({\Vert F\Vert }_{L^2\left(\Omega \right)}+C_1{\Vert G\Vert }_{L^2\left({\Gamma }_N\right)}\right){\Vert v_{\epsilon }\Vert }_{H^1\left(\Omega \right)}+\left|AC\right|\mathrm{.}$ (32)

Combining (32) and property (ii) of ${\varphi }_{\epsilon }$, for any $0<\epsilon <{\epsilon }_0$, we get

$\begin{array}{c}{\displaystyle {\int }_{\Omega }}{\left|\nabla v_{\epsilon }\right|}^2\text{d}x\le {\displaystyle {\int }_{\left[\left|\nabla v_{\epsilon }\right|\le \rho \right]}}{\left|\nabla v_{\epsilon }\right|}^2\text{d}x+{\displaystyle {\int }_{\left[\left|\nabla v_{\epsilon }\right|>\rho \right]}}{\left|\nabla v_{\epsilon }\right|}^2\text{d}x\\ \le {\displaystyle {\int }_{\left[\left|\nabla v_{\epsilon }\right|\le \rho \right]}}{\left|\nabla v_{\epsilon }\right|}^2\text{d}x+\frac{1}{\epsilon }{\displaystyle {\int }_{\Omega }}{\left(\left|\nabla v_{\epsilon }\right|-1\right)}^{+}\left|\nabla v_{\epsilon }\right|\text{d}x\\ \le {\left|\rho \right|}^2\left|\Omega \right|+\left({\Vert F\Vert }_{L^2\left(\Omega \right)}+C_1{\Vert G\Vert }_{L^2\left({\Gamma }_N\right)}\right){\Vert v_{\epsilon }\Vert }_{H^1\left(\Omega \right)}+\left|AC\right|.\end{array}$ (33)

Thus, adding (31) and (33), we obtain

${\Vert v_{\epsilon }\Vert }_{H^1\left(\Omega \right)}^2\le {\left|\rho \right|}^2\left|\Omega \right|+2\left({\Vert F\Vert }_{L^2\left(\Omega \right)}+C_1{\Vert G\Vert }_{L^2\left({\Gamma }_N\right)}\right){\Vert v_{\epsilon }\Vert }_{H^1\left(\Omega \right)}+2\left|AC\right|\mathrm{.}$

Using the Young’s inequality, we get

${\Vert v_{\epsilon }\Vert }_{H^1\left(\Omega \right)}^2\le {\left|\rho \right|}^2\left|\Omega \right|+2\left|AC\right|+\frac{1}{2}{\left[2\left({\Vert F\Vert }_{L^2\left(\Omega \right)}+C_1{\Vert G\Vert }_{L^2\left({\Gamma }_N\right)}\right)\right]}^2+\frac{1}{2}{\Vert v_{\epsilon }\Vert }_{H^1\left(\Omega \right)}^2\mathrm{,}$

Consequently

${\Vert v_{\epsilon }\Vert }_{H^1\left(\Omega \right)}^2\le 2{\left|\alpha \right|}^2\left|\Omega \right|+4\left|A\ast C\right|+4{\left({\Vert F\Vert }_{L^2\left(\Omega \right)}+C_1{\Vert G\Vert }_{L^2\left({\Gamma }_N\right)}\right)}^2\mathrm{.}$ (34)

Then $\left(v_{\epsilon }\right)$ is bounded in $H_D^1\left(\Omega \right)$.

(2) Using (32) and the property (iii) of $\left({\varphi }_{\epsilon }\right)$, we have

$\begin{array}{c}{\displaystyle {\int }_{\Omega }}\left|w_{\epsilon }\right|\text{d}x={\displaystyle {\int }_{\Omega }}\left|{\varphi }_{\epsilon }\left(\nabla v_{\epsilon }\right)\right|\text{d}x\\ \le {\displaystyle {\int }_{\Omega }}{\varphi }_{\epsilon }\left(\nabla v_{\epsilon }\right)\cdot \nabla v_{\epsilon }\text{d}x\\ \le \frac{1}{\epsilon }{\displaystyle {\int }_{\Omega }}{\left(\left|\nabla v_{\epsilon }\right|-1\right)}^{+}\left|\nabla v_{\epsilon }\right|\text{d}x\\ \le \left({\Vert F\Vert }_{L^2\left(\Omega \right)}+C_1{\Vert G\Vert }_{L^2\left({\Gamma }_N\right)}\right){\Vert v_{\epsilon }\Vert }_{H^1\left(\Omega \right)}\end{array}$ (35)

and since $\left(v_{\epsilon }\right)$ is bounded in $H_D^1\left(\Omega \right)$, it follows that $\left(w_{\epsilon }\right)$ is bounded in ${\left(L^1\left(\Omega \right)\right)}^N$.

(3) Now, let $B\subseteq \Omega$ be a fixed Borel set. We have

$\begin{array}{c}{\Vert \nabla v_{\epsilon }\Vert }_{L^1\left(B\right)}\le {\Vert {\left(\left|\nabla v_{\epsilon }\right|-1\right)}^{+}+1\Vert }_{L^1\left(B\right)}\\ \le {\Vert {\left(\left|\nabla v_{\epsilon }\right|-1\right)}^{+}\Vert }_{L^1\left(B\right)}+{\Vert 1\Vert }_{L^1\left(B\right)}\\ \le {\displaystyle {\int }_B}\left|{\left(\left|\nabla v_{\epsilon }\right|-1\right)}^{+}\right|\text{d}x+\left|B\right|\\ \le {\displaystyle {\int }_B}{\left(\left|\nabla v_{\epsilon }\right|-1\right)}^{+}\left|\nabla v_{\epsilon }\right|\text{d}x+\left|B\right|\\ \le \epsilon \left({\Vert F\Vert }_{L^2\left(\Omega \right)}+C_1{\Vert G\Vert }_{L^2\left({\Gamma }_N\right)}\right){\Vert v_{\epsilon }\Vert }_{H^1\left(\Omega \right)}+\left|B\right|.\end{array}$

Letting $\epsilon \to 0$ and using the fact that $v_{\epsilon }$ is bounded in $H^1\left(\Omega \right)$, we obtain

$\underset{\epsilon \to 0}{\lim \inf }{\displaystyle {\int }_{\Omega }}\left|\nabla v_{\epsilon }\right|\text{d}x\le \left|B\right|.$

$\square$

Proof of the Theorem 8. Thanks to Lemma 9, there exists $\tilde{v}$ in $H_D^1\left(\Omega \right)$, $\omega \in {\left({\mathcal{M}}_b\left(\Omega \right)\right)}^N$ and a subsequence that we denote again by $\epsilon$, such that

$w_{\epsilon }\to \omega \text{in}{\left({\mathcal{M}}_b\left(\Omega \right)\right)}^N-{\text{weak}}^{\ast }\mathrm{,}$

$v_{\epsilon }\to \tilde{v}\text{in}{H}_D^1\left(\Omega \right)-\text{weakandin}{L}^2\left(\Omega \right)\mathrm{.}$

Hence, by problem (26), we deduce that

$\text{div}\left(w_{\epsilon }\right)\to \tilde{v}-F\text{in}{L}^2\left(\Omega \right)\mathrm{.}$

We define the set $A_{\delta }$ by $A_{\delta }=\left[\left|\nabla \tilde{v}\right|\ge 1+\delta \right]$, with $\delta >0$. Using the fact that $\nabla v_{\epsilon }\to \nabla \tilde{v}$ in ${\left(L^1\left(\Omega \right)\right)}^s$ -weak as $\epsilon \to 0$, it follows that

$\left(1+\delta \right)\left|A_{\delta }\right|\le {\displaystyle {\int }_{A_{\delta }}}\left|\nabla \tilde{v}\right|\text{d}x\le \lim \underset{\epsilon \to 0}{\inf }{\displaystyle {\int }_{A_{\delta }}}\left|\nabla v_{\epsilon }\right|\text{d}x\mathrm{.}$ (36)

Taking $B=A_{\delta }$ in Lemma 7, we obtain

$\left(1+\delta \right)\left|A_{\delta }\right|\le \left|A_{\delta }\right|$ (37)

which implies that $\left|A_{\delta }\right|=0$. Since $\delta >0$ is arbitrary, then, it follows that $\left|\nabla \tilde{v}\right|\le 1$ a.e. in $\Omega$. Let us now show that $\tilde{v}$ is also a solution of (4). For all $z\in \mathcal{K}$, we have

$\begin{array}{l}{\displaystyle {\int }_{\Omega }}\left(F-\tilde{v}\right)\left(z-\tilde{v}\right)\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}G\left(z-\tilde{v}\right)\text{d}s+A\left({z|}_{{\Gamma }_2}-{\tilde{v}|}_{{\Gamma }_2}\right)\\ =\underset{\epsilon \to 0}{\lim }\left[{\displaystyle {\int }_{\Omega }}{\varphi }_{\epsilon }\cdot \left(\nabla v_{\epsilon }\right)\left(z-\tilde{v}\right)\text{d}x-{\displaystyle {\int }_{{\Gamma }_1}}G\left(z-v_{\epsilon }\right)\text{d}s+A\left({z|}_{{\Gamma }_2}-{v_{\epsilon }|}_{{\Gamma }_2}\right)\right]\\ +{\displaystyle {\int }_{{\Gamma }_1}}G\left(z-\tilde{v}\right)\text{d}s+A\left({z|}_{{\Gamma }_2}-{\tilde{v}|}_{{\Gamma }_2}\right)\end{array}$

$\begin{array}{l}=\underset{\epsilon \to 0}{\lim }{\displaystyle {\int }_{\Omega }}{\varphi }_{\epsilon }\left(\nabla v_{\epsilon }\right)\nabla \left(z-\tilde{v}\right)\text{d}x-{\displaystyle {\int }_{{\Gamma }_1}}G\left(z-\tilde{v}\right)ds+A\left({z|}_{{\Gamma }_2}-{\tilde{v}|}_{{\Gamma }_2}\right)\\ +{\displaystyle {\int }_{{\Gamma }_1}}G\left(z-\tilde{v}\right)\text{d}s+A\left({z|}_{{\Gamma }_2}-{\tilde{v}|}_{{\Gamma }_2}\right)\\ =\underset{\epsilon \to 0}{\lim }{\displaystyle {\int }_{\Omega }}{\varphi }_{\epsilon }\left(\nabla v_{\epsilon }\right)\nabla \left(z-\tilde{v}\right)\text{d}x\\ =\underset{\epsilon \to 0}{\lim }{\displaystyle {\int }_{\Omega }}\left({\varphi }_{\epsilon }\left(\nabla v_{\epsilon }\right)-{\varphi }_{\epsilon }\left(\nabla z\right)\right)\nabla \left(z-\tilde{v}\right)\text{d}x\le \mathrm{0,}\end{array}$ (38)

i.e. $\tilde{v}$ is also a solution of (4). Since (4) admits a unique solution then $v=\tilde{v}$. Now, we must show that $w_{\epsilon }$ satisfies (22). By (iii), we have

$\begin{array}{c}\underset{\epsilon \to 0}{\lim \sup }{\displaystyle {\int }_{\Omega }}\left|{\omega }_{\epsilon }\right|\text{d}x=\lim \underset{\epsilon \to 0}{\sup }{\displaystyle {\int }_{\Omega }}\left|{\varphi }_{\epsilon }\left(\nabla v_{\epsilon }\right)\right|\text{d}x\\ \le \underset{\epsilon \to 0}{\lim \sup }{\displaystyle {\int }_{\Omega }}{\varphi }_{\epsilon }\left(\nabla v_{\epsilon }\right)\cdot \nabla v_{\epsilon }\text{d}x\\ \le \lim \underset{\epsilon \to 0}{\sup }\left({\displaystyle {\int }_{\Omega }}\left(F-v_{\epsilon }\right)v_{\epsilon }\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}G{v}_{\epsilon }\text{d}s+{A{v}_{\epsilon }|}_{{\Gamma }_2}\right)\\ \le {\displaystyle {\int }_{\Omega }}\left(F-v\right)v\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}Gv\text{d}s+{Av|}_{{\Gamma }_2}.\end{array}$ (39)

Moreover taking v as test function in (27) and having in mind that $v=\tilde{v}$, we get

${\displaystyle {\int }_{\Omega }}\left(F-v\right)v\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}Gv\text{d}s+{Av|}_{{\Gamma }_2}$ (40)

$=\underset{\epsilon \to 0}{\lim }{\displaystyle {\int }_{\Omega }}\left(F-v_{\epsilon }\right)v\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}Gv\text{d}s+{Av|}_{{\Gamma }_2}$

$=\underset{\epsilon \to 0}{\lim }{\displaystyle {\int }_{\Omega }}{w}_{\epsilon }\cdot \nabla v\text{d}x$

$=\underset{\epsilon \to 0}{\lim }{\displaystyle {\int }_{\Omega }}{w}_{\epsilon }\cdot \nabla v\text{d}x$

$\le \underset{\epsilon \to 0}{\lim }{\displaystyle {\int }_{\Omega }}\left|w_{\epsilon }\right|\text{d}x.$ (41)

Combining (39) and (40), we get

$\underset{\epsilon \to 0}{\lim }{\displaystyle {\int }_{\Omega }}\left|w_{\epsilon }\right|\text{d}x={\displaystyle {\int }_{\Omega }}\left(F-v\right)v\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}Gv\text{d}s+{Av|}_{{\Gamma }_2}\mathrm{.}$ (42)

Thanks to the Lemma 7 we have $-G\left(w_{\epsilon }\right)\le {\min }_{z\in \mathcal{K}}J\left(z\right)$, hence from (22) and (23), we deduce that, as $\epsilon$, $G\left(w_{\epsilon }\right)\to -{\min }_{z\in \mathcal{K}}J\left(z\right)$. $\square$

Remark. As consequences of Theorem 8, we have another characterization of the solution v of the problem (4):

$\begin{array}{l}v-\text{div}\left(\omega \right)=g\text{in}{\left(H_D^1\left(\Omega \right)\right)}^{*}\\ \left|\omega \right|\left(\Omega \right)={\displaystyle {\int }_{\Omega }}\left|{\omega }_{\epsilon }\right|\text{d}x\to {\displaystyle {\int }_{\Omega }}\left(F-v\right)v\text{d}x+{\displaystyle {\int }_{{\Gamma }_1}}Gv\text{d}s+{Av|}_{{\Gamma }_2}.\end{array}$ (43)

Moreover, we have another characterization of $H_{\text{div}\mathrm{,}\mathcal{K}}\left(\Omega \right)$ which will be useful for the numerical part:

$H_{\text{div}\mathrm{,}K}\left(\Omega \right)=\left\{\varphi +\nabla L\mathrm{,}\varphi \in \mathcal{H}\right\}$

where L is the variational solution to the Laplace problem

$\begin{array}{l}-\Delta L=0\text{in}\Omega \\ L=0\text{on}{\Gamma }_0\\ \frac{\partial L}{\partial \nu }=G\text{on}{\Gamma }_1\mathrm{,}\\ L\equiv \text{constanton}{\Gamma }_2\\ \frac{\partial L}{\partial \nu }=\frac{A}{\left|{\Gamma }_2\right|}\text{on}{\Gamma }_2\mathrm{.}\end{array}$ (44)

and

$\mathcal{H}=\left\{\phi \in H_{\text{div}}\left(\Omega \right)\mathrm{<mo>\; :\; </mo>}-{\displaystyle {\int }_{\Omega }}\text{div}\left(\phi \right)\xi ={\displaystyle {\int }_{\Omega }}\phi \cdot \nabla \xi \mathrm{,}\forall \xi \in H_D^1\left(\Omega \right)\right\}\mathrm{.}$

Consequently, the dual problem (18) is equivalent to

$\inf \left\{\mathcal{G}\left(\varphi \right)=G\left(\varphi +\nabla L\right)\mathrm{<mo>\; :\; </mo>}\varphi \in \mathcal{H}\right\}\mathrm{.}$ (45)

For the numerical approximation of the problem (45), we use the finite element method, then we assume that:

· $\Omega$ is an open and bounded subset of ${ℝ}^2$ ;

· $T_h$ will be a regular partitioning of $\bar{\Omega }$ by n disjoint open simplex $\tau$, with $\bar{\Omega }={\displaystyle {\cup }_{\tau \in T_h}\bar{\tau }}$.

Let ${\mathcal{H}}_h$ be a finite dimensional subspace of $R{T}_0\left(T_h\right)\cap \mathcal{H}$, with a finite dimension equal to $N=N\left(h\right)$, where $R{T}_0\left(T_h\right)$ is the space of lowest-order Raviart-Thomas finite elements (see [25] ) defined by

$\begin{array}{l}{\mathcal{H}}_h=\{q_h\in {\left(L^2\left(\Omega \right)\right)}^2\mathrm{<mo>\; :\; </mo>}{q}_{\tau }^h=a_{\tau }+b_{\tau }x\mathrm{,}a\in {ℝ}^2\mathrm{,}b\in ℝ\mathrm{,}\forall \tau \in T_h\mathrm{,}\\ \text{and}\left(q_{\tau }^h-q_{{\tau }^{\prime }}^h\right)\cdot {\nu }_{{\partial }_{\tau }}=0\text{on}{\partial }_{\tau }\cap {\partial }_{{\tau }^{\prime }}\}\mathrm{,}\end{array}$

where ${\nu }_{{\partial }_{\tau }}$ represents the outward unit normal to ${\partial }_{\tau }$, the boundary of $\tau$.

By $r_h$ we denote the interpolation operator onto ${\mathcal{H}}_h$ given in ( [25], Theorem 6.1). Then, thanks to [25], for all $w\in H_{\text{div}}\left(\Omega \right)$, we have

$r_h\left(w\right)\to w\text{in}{\left(L^2\left(\Omega \right)\right)}^2\text{and}\text{div}\left(r_h\left(w\right)\right)\to \text{div}\left(\omega \right)\text{in}{L}^2\left(\Omega \right)\mathrm{,}\text{as}h\to 0.$ (46)

We have the following convergence properties as h tends to 0.

Theorem 10. Let $F\in L^2\left(\Omega \right)$ , $G\in L^2\left({\Gamma }_1\right)$ , v a solution to the minimization problem (14) and $w_h$ a solution of the optimization problem

$\sup \left\{-\mathcal{G}\left(q_h\right)\mathrm{<mo>\; :\; </mo>}{q}_h\in {\mathcal{H}}_h\right\}\mathrm{.}$ (47)

Then, as $h\to 0$,

$\text{div}\left({\omega }_h\right)\to v-\tilde{f}\text{in}{L}^2\left(\Omega \right)$ (48)

and

$-\mathcal{G}\left({\omega }_h\right)\to \underset{z\in K}{\min }J\left(z\right)=\frac{1}{2}{\displaystyle {\int }_{\Omega }}{\left|v-g\right|}^2\text{d}x-{\displaystyle {\int }_{{\Gamma }_N}}\left({\displaystyle {\int }_{\Omega }}K\left(s\mathrm{,}y\right)L\left(y\right)\text{d}y-L\left(s\right)\right)v\left(s\right)\text{d}s\mathrm{.}$ (49)

The proof of this result is similar to the proof of Theorem 3.8 in the paper of Igbida et al.; therefore, we omit it.

4. Numerical Simulations

We start by approximating the term ${\displaystyle {\int }_{\tau }}\left|q_h+\nabla Z\right|\text{d}x$ on each element $\tau$ of the partitioning $T_h$. We take

${\displaystyle {\int }_{\tau }}\left|q_h+\nabla L\right|\text{d}x\sim \left|\tau \right|\left|q_h+\nabla L\right|\left(P_{\tau }\right)\mathrm{,}$

where $\left|\tau \right|$ represents the area of simplex $\tau$ and $\left(P_{\tau }\right)$ is one of the vertices of $\tau$. Using this approximation, at each time $n\times dt$, $n\in \mathbb{N}$ and $dt$ the time step, the solution of (47) is a minimizer of the non-differentiable functional $\mathcal{G}\mathrm{<mo>\; :\; </mo>}{ℝ}^s\to ℝ$,

$\mathcal{G}\left(w_h\right)\mathrm{<mo>\; :\; </mo>}=\frac{1}{2}\left(A{w}_h\mathrm{,}{w}_h\right)+\left(dt{f}_h^n+u_{n-1}\mathrm{,}\text{div}\left(w_h\right)\right)+{\displaystyle \underset{\tau \in T_h}{\sum }}\left|\tau \right|\left|w_h+\nabla L\right|\left(P_{\tau }\right)\mathrm{.}$

The minimization of the functional $\mathcal{G}$ is done according to the Gauss-Seidel type algorithm.

In all the tests the discretization in time $dt=0.001$ and convergence criterion $\epsilon ={10}^{-6}$, $\Omega ={\left(-\mathrm{1,1}\right)}^2$. For the discretization of the problem, we use the Raviart-Thomas elements of the lowest order [25] on a regular square grid.

After having obtained the minimizer $w_h$ of the functional $\mathcal{G}$, the solution $u_n$ of the Euler implicit time discretization of (1) is computed by using extremality relation (43) in a weak sense with piecewise finite elements $P_0$.

In the first test, the source $f\left(x\mathrm{,}t\right)$ is a constant equal to 1, for all time and for all $x\in \Omega$. The Neumann boundary conditions $q\cdot \nu =m\frac{\partial \nabla u}{\partial \nu }=0$ are considered on the boundary $\left\{\left(-\mathrm{1,}y\right)|y\in \left]-\mathrm{1,1}\right[\right\}$ which simulated the existence of a wall on the boundary. The outgoing flow across the border ${\Gamma }_1=\left\{\left(\mathrm{1,}y\right)\mathrm{<mo>\; |\; </mo>}y\in \left]-\mathrm{1,1}\right[\right\}$ at any time t is is ${\displaystyle {\int }_{{\Gamma }_1}}m\frac{\partial \nabla u}{\partial \nu }\text{d}s=1$.

Figure 1 shows the configuration of the sandpile at time $t=1.520$ when the solution u becomes stationary. On sees also the effect of the non-local boundary condition of the configuration of the sandpile.

In the second test (see Figure 2), one analyzes the impact of the condition of Neumann on the configuration of the sandpile in the absence of the source f. We can also see that the flux is concentrated around ${\Gamma }_0=\left\{\left(-\mathrm{1,}y\right)\mathrm{<mo>\; |\; </mo>}y\in \left]-\mathrm{1,1}\right[\right\}$. This example allows us to conclude that the function g behaves likes a source of sand.

In the last simulation, the three data f, g and A act at the same time. Figure 3 and Figure 4 show respectively an intermediate and the final configuration

Figure 1. Final configuration of the sandpile surface at $t=1.520$ for $f\equiv 1$, $A\equiv 1$ and $g=0$.

Figure 2. Sandpile surface and the flux at $t=15$ for $f\equiv 0$, $A\equiv 0$ and $g\equiv 20$.

Figure 3. Sandpile surface and the flux at $t=0.787$ for $f\equiv 1$, $A\equiv -1$ and $g\equiv 0.1$.

Figure 4. Sandpile surface and the flux at $t=15$ for $f\equiv 1$, $A\equiv -1$ and $g\equiv 0.1$.

of the sandpile and the behavior of the flow. We can see here that the final configuration of the sandpile is slightly different from that of the test 1 (see Figure 1). This is due to the sign of the total outflow through ${\Gamma }_1$. We also remark that the flux has the same direction as the gradient of the surface and almost zero on the axis $y=0$.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Igbida, N., Karami, F., Ouaro, S. and Traoré, U. (2017) On a Dual Formulation for Growing Sandpile Problem with Mixed Boundary Conditions. Applied Mathematics & Optimization, 78, 329-359. https://doi.org/10.1007/s00245-017-9408-2
[2]	Bouchaud, J.P., Cates, M.E., Ravi Prakash, J. and Edwards, S.F. (1994) A Model for the Dynamics of Sandpile Surfaces. Journal de Physique, 4, 1383-1410. https://doi.org/10.1051/jp1:1994195
[3]	Boutreux, T. and de Gennes, P.G. (1996) Surface Flows of Granular Mixture, I. General Principles and Minimal Model. Journal de Physique, 6, 1295-1304. https://doi.org/10.1051/jp1:1996138
[4]	Savage, S.B. and Hutter, K. (1989) The Motion of a Finite Mass of Granular Material Down a Rough Incline. Journal of Fluid Mechanics, 199, 177-215. https://doi.org/10.1017/S0022112089000340
[5]	Prigozhin, L. and Barrett, J.W. (2006) Dual Formulation in Critical State Problems. Interface Free Bound, 8, 349-370. https://doi.org/10.4171/IFB/147
[6]	Prigozhin, L. (1996) Variational Model of Sandpile Growth. European Journal of Applied Mathematics, 7, 225-236. https://doi.org/10.1017/S0956792500002321
[7]	Dumont, S. and Igbida, N. (2011) On the Collapsing Sandpile Problem. Communications on Pure & Applied Analysis, 10, 625-638. https://doi.org/10.3934/cpaa.2011.10.625
[8]	Dumont, S. and Igbida, N. (2009) On a Dual Formulation for Growing Sandpile Problem. European Journal of Applied Mathematics, 20, 169-185. https://doi.org/10.1017/S0956792508007754
[9]	Igbida, N. (2010) A Generalized Collapsing Sandpile Model. Archiv der Mathematik, 94, 193-200. https://doi.org/10.1007/s00013-009-0038-z
[10]	Nassouri, E., Ouaro, S. and Traoré, U. (2017) Growing Sandpile with Dirichlet and Fourier Boundary Conditions. Electronic Journal of Differential Equations, 2017, 1-19.
[11]	Giroire, J., Ha-Duong, T. and Moumas, V. (2005) A Non-Linear and Non-Local Boundary Condition for a Diffusion Equation in Petroleum Engineering. Mathematical Methods in the Applied Sciences, 28, 1527-1552. https://doi.org/10.1002/mma.622
[12]	Ta-Tsien, L. (1989) A Class of Non-Local Boundary Value Problems for Partial Differntial Equations and It’s Application in Numerical Analysis. Journal of Computational and Applied Mathematics, 28, 49-62. https://doi.org/10.1016/0377-0427(89)90320-8
[13]	Day, W.A. (1982) A Decreasing Property of Solution of a Parabolic Equation with Applications to Thermoelasticity. Quarterly of Applied Mathematics, 40, 319-330. https://doi.org/10.1090/qam/678203
[14]	Day, W.A. (1985) Heat Conduction within Linear Thermoelasticity. Springer-Verlag, New York. https://doi.org/10.1007/978-1-4613-9555-3
[15]	Friedman, A. (1986) Monotonic Decay of Solution of Parabolic Equations with Nonlocal Boundary Conditions. Quarterly of Applied Mathematics, 44, 401-407. https://doi.org/10.1090/qam/860893
[16]	Pao, C.V. (1995) Reaction Diffusion Equations with Nonlocal Boundary and Nonlocal Initial Conditions. Journal of Mathematical Analysis and Applications, 195, 702-718. https://doi.org/10.1006/jmaa.1995.1384
[17]	Pao, C.V. (1995) Dynamics of Reaction-Diffusion Equations with Nonlocal Boundary Conditions. Quarterly of Applied Mathematics, 53, 173-186. https://doi.org/10.1090/qam/1315454
[18]	Burghardt, R. (2015) Collapsing Schwarzschild Interior Solution. Journal of Modern Physics, 6, 1895-1907. https://doi.org/10.4236/jmp.2015.613195
[19]	Shan, M., Zhu, Y., Yao, C., Han, Q. and Zhu, C. (2017) Modeling for Collapsing Cavitation Bubble Near Rough Solid Wall by Mulit-Relaxation-Time Pseudopotential Lattice Boltzmann Model. Journal of Applied Mathematics and Physics, 5, 1243-1256. https://doi.org/10.4236/jamp.2017.56106
[20]	Yang, J., Shen, Z., Zheng, X. and Li, L. (2015) Simulation on Cavitation Bubble Collapsing with Lattice Boltzmann Method. Journal of Applied Mathematics and Physics, 3, 947-955. https://doi.org/10.4236/jamp.2015.38116
[21]	Brezis, H. (1973) Opérateurs maximaux monotones et sémi-groupes de contractions dans les espaces de Hilbert. (French). North-Holland Mathematics Studies, No. 5. Notas de Matermatica (50). North-Holland Publishing Co., Amsterdam-London, American Elsevier Publishing Co., Inc., New York.
[22]	Ekeland, I. and Temam, R. (1999) Convex Analysis and Variational Problems. Classics in Applied Mathematics, 28. Society for Industrial and Applied Mathematics (SIAM), Philadelphia. https://doi.org/10.1137/1.9781611971088
[23]	Lions, J.L. (1996) Quelques methodes de rsolution des problmes aux limites non linaires. Dunod, Gauthier-Villars, Paris.
[24]	Evans, L.C. (1998) Partial Differential Equations. Graduate Studies in Mathematics, Vol. 19, American Mathematical Society, Providence.
[25]	Robert, J.E. and Thomas, J.M. (1991) Mixed and Hybrid Methods. In: Ciarlet, P.G. and Lions, J.L., Eds., Handbook of Numerical Analysis, Finite Elements Methods (Part 1), Vol. 2, North-Holland, Amsterdam, 523-639. https://doi.org/10.1016/S1570-8659(05)80041-9