\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2017 (2017), No. 300, pp. 1--19.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2017 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2017/300\hfil Growing sandpile problem] {Growing sandpile problem with Dirichlet and Fourier boundary conditions} \author[E. Nassouri, S. Ouaro, U. Traor\'e \hfil EJDE-2017/300\hfilneg] {Estelle Nassouri, Stanislas Ouaro, Urbain Traor\'e} \address{Estelle Nassouri \newline Laboratoire de Math\'ematiques et Informatique (LAMI), UFR, Sciences Exactes et Appliqu\'ees, Universit\'e Ouaga I Pr Joseph Ki-Zerbo, 03 BP 7021 Ouaga 03, Ouagadougou, Burkina Faso} \email{nassouristella@yahoo.fr} \address{Stanislas Ouaro \newline Laboratoire de Math\'ematiques et Informatique (LAMI), UFR, Sciences Exactes et Appliqu\'ees, Universit\'e Ouaga I Pr Joseph Ki-Zerbo, 03 BP 7021 Ouaga 03, Ouagadougou, Burkina Faso} \email{ouaro@yahoo.fr, souaro@univ-ouaga.bf} \address{Urbain Traor\'e \newline Laboratoire de Math\'ematiques et Informatique (LAMI), UFR, Sciences Exactes et Appliqu\'ees, Universit\'e Ouaga I Pr Joseph Ki-Zerbo, 03 BP 7021 Ouaga 03, Ouagadougou, Burkina Faso} \email{urbain.traore@yahoo.fr} \dedicatory{Communicated by Vicentiu D. Radulescu} \thanks{Submitted September 9, 2017. Published December 6, 2017.} \subjclass[2010]{35D30, 35K86, 35R37, 35B40} \keywords{Growing sandpile; Fourier boundary condition; nonlinear semi-group; \hfill\break\indent Dirichlet boundary condition; Euler discretization in time} \begin{abstract} In this work, we study the Prigozhin model for growing sandpile with mixed boundary conditions and an arbitrary time dependent angle of repose. On one part of the boundary the homogeneous Dirichlet boundary condition is provided, on the other one the Robin condition is used. Using the implicit Euler discretization in time, we prove the existence and uniqueness of variational solution of the model and for the numerical analysis we use a duality approach. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{remark}[theorem]{Remark} \newtheorem{definition}[theorem]{Definition} \allowdisplaybreaks \section{Introduction} In this work, we study the Prighozin type growing sandpile problem, with mixed boundary conditions. On one part of the boundary the homogeneous Dirichlet boundary condition is applied. This condition means that, on this boundary, the sand can fall down. On the other part of the boundary we use the Robin boundary condition. It is well known that the sand has a limit angle, the so-called angle of repose. It corresponds to the steepest angle which the surface of a mass of particles in bulk make with the ground. In the paper \cite{PV,Ig2, Ig4}, the authors have worked with an angle of stability equal to $\pi/4$. Furthermore, if we assume that the moisture of the material is changing in time, we can assume that the angle of repose is a given time dependent function $ c:t\in [0,T)\to c(t)\in \mathbb{R}^{*}_{+}$. Hence, the model can be described by the following PDE. \begin{equation} \label{eP} \begin{gathered} u_{t}-\nabla(m\nabla u) = f \quad\text{in } [0,T]\times\Omega\\ | \nabla u | \leq c(t), \quad m \geq 0, \quad m(| \nabla u | - c(t))=0 \quad\text{in } [0,T]\times\Omega\\ u = 0 \quad\text{on } [0,T]\times\Gamma_{D}\\ m\frac{\partial u}{\partial \nu} +\lambda u = g \quad\text{on } [0,T]\times\Gamma_{N}\\ u(0)=u_{0}, \end{gathered} \end{equation} where $ \Omega\subset \mathbb{R}^s$, is a bounded open domain. The solution $ u$ is the height of the surface and $ f$ is the source. Let us give a brief description of the model \eqref{eP} under the following assumptions. \begin{itemize} \item The flow of material is confined in a fine layer on the surface of the sand pile, \item the density of the material is constant, \item the dynamic effects are assumed to be negligible. \end{itemize} The conservation law allows us to write $$ \frac{\partial u}{\partial t}+\nabla . q=f \quad \text{in } \Omega; $$ $f$ being the source and $q$ the orthogonal projection of the material flow. Moreover, since the flow of the material flows along the largest slope, then $$ q=-m\nabla u , $$ with $m$ an unknown scalar function verifying with $u$, the following properties. \begin{itemize} \item The angle formed by the sandpile surface and the horizontal is never greater than the angle of stability; i.e. $$ | \nabla u(t,x) |\leq c(t); $$ \item since there is not a flow of material when the slope is not very inclined, we have $$ | \nabla u(t,x) | < c(t) \; \Rightarrow\; m(x,t)=0. $$ \end{itemize} The initial free surface is given by $$ u(x,0)= u_{0}(x). $$ Concerning the boundaries conditions, we have split the boundary of the domain denoted $\Gamma$ into two part $\Gamma_{N}$ et $\Gamma_{D}$ such that $\Gamma_{N} \cap\Gamma_{D}=\emptyset$; \begin{itemize} \item in $\Gamma_{D}$, we apply homogeneous Dirichlet boundary condition. In other words, we assume that at this boundary, the sand falls down, that is $$ u=0 \quad\text{on } [0, T]\times \Gamma_{D}; $$ \item in $\Gamma_{N}$, we consider Fourier boundary condition. In other words, we assume that there exists a wall that prevents the flowing of sand, that is $$ m\frac{\partial u}{\partial \nu} +\lambda u = g \quad\text{on } [0,T]\times\Gamma_{N}, $$ with $\lambda$ being a positive real number. \end{itemize} The PDEs with a mixed boundary appear in the modeling of many phenomena like fluid flows in domains with a boundary, which includes several parts, these parts differ by their physical properties \cite{Mik,Hao}. Let us recall that in the literature, many models of sandpile, with critical slope model (cf. \cite{PV, Ig1,Ig2,Ig3,Ig4}), two layers models (cf. \cite{BCRE94}) were studied. In particular, in \cite{Ig2}, the authors used nonlinear semi-group theory (see \cite{Bre}) to prove the existence and uniqueness of solution to a sandpile problem like \eqref{eP}. And for numerical study of the problem they show that the Euler implicit time discretization problem associated with the model is equivalent to some kind of minimization problem. So, using the dual problem associated with this minimization problem the authors have showed how to compute the solution of the model knowing the solution to the dual problem. Our aim in this work is to generalize the ideas introduced in \cite{Ig2}. This article is organized as follows. In the next Section, we use the nonlinear semi-group theory (see \cite{Bre}) to get the existence and uniqueness of a variational solution of \eqref{eP} and the convergence of the approximate Euler discretization in time solutions to problem \eqref{eP}. In Section 3, we show how to compute the solution of Euler implicit time discretization of \eqref{eP} using duality argument and in Section 4, some results of numerical simulations of \eqref{eP} are given. \section{Theoretical study of problem \eqref{eP}} Let us start by introducing the Euler discretization in time of problem \eqref{eP}. Given $ \varepsilon >0$, we say that $ (t_i,f_i,g_i)_{i=1,\dots,n}$ is an $ \varepsilon$-discretization for problem \eqref{eP} if $ 0= t_00. \end{equation} \begin{definition} \label{def1} \rm For any $\varepsilon >0$, $u_{\varepsilon}$ is an $\varepsilon$-approximate solution of \eqref{eP}, if there exists $(t_i,f_i,g_i)_{i=1,\dots ,n}$ an $\varepsilon-$ discretization of problem \eqref{eP} such that \begin{equation} \label{Ueq2} u_{\varepsilon}= \begin{cases} u_0 & \text{for } t \in (0, t_1] \\ u_i & \text{for } t \in (t_{i-1}, t_i ], \; i=1,\dots ,n \end{cases} \end{equation} and $u_i$ solves the following Euler implicit time discretization of \eqref{eP}: \begin{equation}\label{Ueq3} \begin{gathered} u_i-\varepsilon\nabla(m_i\nabla u_i) = \varepsilon f_i + u_{i-1} \quad \text{in } \Omega\\ | \nabla u_i | \leq c(t_i), \quad m_i \geq 0, \quad m_i (| \nabla u_i | - c(t_i))=0 \quad\text{in } \Omega\\ u_i = 0 \quad\text{on } \Gamma_{D}\\ m\frac{\partial u_i}{\partial \nu} + \lambda u_i = g_i \quad \text{on } \Gamma_{N}. \end{gathered} \end{equation} \end{definition} To simplify the analysis, we introduce the following generic problem associated with \eqref{Ueq3} \begin{equation} \label{eS} \begin{gathered} v-\nabla(m\nabla v) = \widetilde{f} \quad\text{in } \Omega\\ | \nabla v | \leq d , \quad m \geq 0 , \quad m(| \nabla v | - d)=0 \quad \text{in } \Omega\\ v = 0 \quad\text{on } \Gamma_{D}\\ m\frac{\partial v}{\partial \nu} +\lambda v = \widetilde{g} \quad \text{on } \Gamma_{N}, \end{gathered} \end{equation} with $ \widetilde{f}\in L^2(\Omega)$, $\widetilde{g}\in L^2(\Gamma_{N})$, $d=c(t_i)$ and $ \lambda$ a strictly positive real number. For convenience, we note the $ L^2(\Omega)$ scalar product by $ (\cdot,\cdot)$ and the euclidean norm on $ \mathbb{R}^s$ by $ \|\cdot \|$. We introduce the convex set $ K(d)$ given by $$ K(d)=\{z\in W^{1,\infty} (\Omega) \cap H^{1}_{D}(\Omega): | \nabla z (x) |\leq d \text{ a.e } x\in\Omega\} $$ with $$ H^{1}_{D}(\Omega)=\{z\in H^{1}(\Omega); z=0 \ on \ \Gamma_{D}\}. $$ We also introduce the following function defined on $ H^{1}_{D}$ by $$ F_{\widetilde{g}}(z) = \begin{cases} -\int_{\Gamma_{N}} (\widetilde{g}-\lambda Z)zds & \text{if } z\in K(d),\\ +\infty & \text{otherwise}, \end{cases} $$ with $Z$ a variational solution of the following Laplace problem. \begin{equation} \label{eL} \begin{gathered} -\Delta Z=0 \quad\text{in } \Omega\\ Z=0 \quad\text{on } \Gamma_{D}\\ \frac{\partial Z}{\partial\nu} +\lambda Z=\widetilde{g} \quad \text{on } \Gamma_{N}. \end{gathered} \end{equation} The sub-differential of $ F_{\widetilde{g}}$ is defined in $ L^2(\Omega)$ by $$ \partial F_{\widetilde{g}} (v)=\{ w \in K(d): \forall z\in K(d), \, F_{\widetilde{g}}(z)\geq F_{\widetilde{g}}(v) + (w, z-v)\}. $$ We are now in a position to define our notion of solution to problems \eqref{eS} and \eqref{eP}. \begin{definition} \label{def2} \rm For given $ \widetilde{f}\in L^2(\Omega)$, $ \widetilde{g}\in L^2(\Gamma_{N})$, we say that $ v$ is a variational solution of \eqref{eS} if $ v\in K(d)$ and \begin{equation} \label{Ueq4} \int_{\Omega}(\widetilde{f}-v)(z-v)dx + \int_{\Gamma_{N}}(\widetilde{g}-\lambda Z)(z-v)ds \leq 0, \quad\text{for all } z\in K(d). \end{equation} \end{definition} \begin{definition} \label{def3} \rm Given $f\in L^{\infty} (0,T; L^2(\Omega))$, $g\in L^2_{loc}(\Gamma_N)$ and $ u_0 \in K(c(0))$, a variational solution (resp. $ \varepsilon$-approximate solution) of \eqref{eP} is a function $ u\in W^{1,1} (0,T; L^2(\Omega))$ satisfying for any $ t\in (0,T)$, $u(t)\in K(c(t))$ and $$ \int_{\Omega} (f-u_t (t))(z-u(t))dx + \int_{\Gamma_{N}} (g-\lambda Z)(z-v) ds\leq 0, \quad \text{for any } z\in K(c(t)). $$ (resp. $u_{\varepsilon}$ given by \eqref{Ueq2} with $ u_i$ a variational solution of \eqref{Ueq3}). \end{definition} Using the same method given in \cite{Ig2}, the following result can be proved. \begin{lemma} \label{lem1} \begin{itemize} \item[(i)] $ \partial F_{\widetilde{g}}$ is a maximal monotone graph in $ L^2(\Omega)$. \item[(ii)] $ v$ is a solution of \eqref{eS} if and only if $ v$ is a solution of $ v+\partial F_{\widetilde{g}}(v)\ni \widetilde{f}$. \end{itemize} \end{lemma} Since $ \partial F_{\widetilde{g}}$ is a maximal monotone graph in $ L^2(\Omega)$, then, thanks to \cite{Bre} for any $ \widetilde{f}\in L^2(\Omega)$ there exists a unique solution $ v$ of \begin{equation} \label{Ueq5} v+\partial F_{\widetilde{g}}(v)\ni \widetilde{f}, \end{equation} which is equivalent to saying that \eqref{eS} admits a unique variational solution. Moreover, if $ v_i$ is the solution corresponding to $ \widetilde{f}_i$ for $ i = 1, 2$, then \begin{equation} \|v_1-v_2\|_2\leq \|\widetilde{f}_1-\widetilde{f}_2\|_2. \end{equation} To prove the existence and uniqueness of the solution to problem \eqref{eP}, we need the following lemma. \begin{lemma} \label{lem2} Let $ f\in L^{\infty}(0,T; L^2(\Omega))$, $ g\in L^2_{loc}(\Gamma_{N})$, $c:t\in [0,T)\to c(t)\in \mathbb{R}^{*}_{+}$ and $u_{0}\in K(c(0))$. Then, $u$ is a solution of \eqref{eP} if and only if $ v(t,.)=\frac{u(t,.)}{c(t)}$ is a solution of \begin{equation} \label{eP1} \begin{gathered} v_{t}+\partial G_{g/c} (v)\ni h\\ v(0)=\frac{ u(0)}{ c(0)}\,, \end{gathered} \end{equation} with $h=\frac{f}{ c(t)}-\frac{ c^{'}(t)}{ c(t)}v(t)$, for any $t\in [0, T)$ and $$ G_{g/c}(z):=\begin{cases} -\int_{\Gamma_{N}}\frac{ (g-\lambda Z)}{ c(t)} z & \text{if } z\in K(1)\\ 0 & \text{otherwise}. \end{cases} $$ \end{lemma} \begin{proof} If $ u$ is a solution of \eqref{eP}, then \begin{equation} \int_{\Omega}(f-u_{t})(z-u)dx +\int_{\Gamma_{N}}(g-\lambda Z)(z-u)ds\leq 0. \end{equation} For any $t\in [0,T)$, we denote $$ A:=\int_{\Omega}(f-u_{t})(z-u)dx +\int_{\Gamma_{N}}(g-\lambda Z)(z-u)ds. $$ We take $t\in [0,T)$, $v(t,\cdot)=\frac{ u(t,\cdot)}{ c(t)}$ to obtain \begin{align*} A &= \int_{\Omega}[f-(v_{t}c(t)+v(t)c^{'}(t))](z-v(t)c(t))dx +\int_{\Gamma_{N}}(g-\lambda Z)(z-v(t)c(t))ds \\ &= \int_{\Omega}[f-v_{t}c(t)-v(t)c^{'}(t)](z-v(t)c(t))dx +\int_{\Gamma_{N}}(g-\lambda Z)(z-v(t)c(t))ds \\ &= \int_{\Omega}(c(t))^2(\frac{f}{c(t)}-v_{t} -\frac{v(t)c^{'}(t)}{c(t)} )(\frac{z}{c(t)}-v(t))dx \\ &\quad +\int_{\Gamma_{N}}(c(t))^2 \frac{g-\lambda Z}{c(t)}(\frac{z}{c(t)}-v(t))ds. \end{align*} Since $A\leq 0$ for any $t\in [0,T)$, we have \[ \int_{\Omega}(\frac{f}{c(t)}-v_{t}-\frac{v(t)c^{'}(t)}{c(t)} ) (\frac{z}{c(t)}-v(t))dx +\int_{\Gamma_{N}} \frac{g-\lambda Z}{c(t)}(\frac{z}{c(t)}-v(t))ds\leq 0, \] for any $t\in [0,T)$. Since $z\in K(c(t))$, then $w=\frac{ z}{ c(t)}\in K(1)$ for any $t\in [0,T)$ and, it follows that $$ \int_{\Omega}(h-v_{t})(w-v(t))dx +\int_{\Gamma_{N}} \frac{g-\lambda Z}{c(t)}(w-v(t))ds\leq 0; $$ which is equivalent to saying that $ v$ is a solution of \eqref{eP1}. Now, we suppose that $v$ is a solution of \eqref{eP1}, then $$ \int_{\Omega}(h-v_{t})(w-v(t))dx +\int_{\Gamma_{N}} \frac{g-\lambda Z}{c(t)}(w-v(t))ds\leq 0,\quad \text{for any } z\in K(1). $$ Using $v(\cdot,t)=\frac{ u(\cdot,t)}{ c(t)}$ and denoting $$ B=\int_{\Omega}(h-v_{t})(w-v(t))dx +\int_{\Gamma_{N}} \frac{g-\lambda Z}{c(t)}(w-v(t))ds, $$ we get \begin{equation} \begin{aligned} B &= \int_{\Omega}(h-v_{t})(w-v(t))dx +\int_{\Gamma_{N}} \frac{g-\lambda Z}{c(t)}(w-v(t))ds \\ &= \int_{\Omega}(\frac{f}{c(t)}-\frac{c^{'}(t)}{c(t)^2}u(t) -\frac{u_{t}c(t)-u(t)c^{'}(t)}{c(t)^2})(w-\frac{ u(t)}{ c(t)})dx \\ &\quad +\int_{\Gamma_{N}} \frac{g-\lambda Z}{c(t)}(w-\frac{ u(t)}{ c(t)})ds \\ &= \int_{\Omega}(\frac{f}{c(t)}-\frac{u_{t}}{c(t)})(w-\frac{ u(t)}{ c(t)})dx +\int_{\Gamma_{N}} \frac{g-\lambda Z}{c(t)}(w-\frac{ u(t)}{ c(t)})ds. \end{aligned} \end{equation} Since $B\leq 0 $ and $(c(t))^2B\leq 0$, we get $$ \int_{\Omega}( f-u_{t})(c(t)w- u(t) )dx +\int_{\Gamma_{N}} (g-\lambda Z)(c(t)w-u(t))ds\leq 0, \quad \text{for any } w\in K(1). $$ Since $w\in K(1)$, then $c(t)w\in K(c(t))$; therefore $$ \int_{\Omega}( f-u_{t})(z- u(t) )dx +\int_{\Gamma_{N}} (g-\lambda Z)(z-u(t))ds\leq 0, \quad \text{for any } z\in K(c(t)). $$ So, $u$ is a solution of \eqref{eS} \end{proof} \begin{theorem} \label{thm2.6} Let $ f\in L^{\infty}(0,T; L^2(\Omega))$ and $ g\in L^2_{loc}(\Gamma_{N})$, then \eqref{eP} has a unique variational solution $u\in W^{1,1}(0,T;L^2(\Omega))$; and for any subsequence $ \varepsilon\rightarrow 0$, if $ u_{\varepsilon}$ is a $\varepsilon$-approximate solution of \eqref{eP}, then $$ u_{\varepsilon}\to u\in \mathcal{C}([0,T); L^2(\Omega)). $$ Moreover, if for $i=1, 2$, $ u_i$ is a solution corresponding to $ f_i$, then $$ \frac{d}{dt}\int_{\Omega}(u_1-u_2)^{+} \leq \int_{\Omega} (f_1-f_2)^{+}\quad \text{in }\mathcal{D}^{\prime}(0,T). $$ In particular, if $ f\geq 0$ and $ g\geq 0$, then $ u\geq 0$ a.e. in $\Omega$. \end{theorem} \begin{proof} Using \cite[Proposition 3.13]{Bre}, we deduce that for any $v_{0}\in K(1)$, \eqref{eP1} has a unique solution $ v\in W^{1,1}(0,T;L^2(\Omega))$, then from Lemma \ref{lem2}, we deduce that \eqref{eP} has a unique variational solution. To prove the convergence of the $\varepsilon$-approximate solution, let us consider, for $ i= 1,\dots,n$, \begin{equation} \label{Ueq7} u_i + \partial F_{g, c(t_i)}(u_i ) \ni \varepsilon f_i + u_{i-1}, \end{equation} where $$ F_{g,c(t_i)}(z)= \begin{cases} -\int_{\Gamma_{N}} (g-\lambda Z)zds & \text{if } z\in K(c(t_i))\\ +\infty &\text{otherwise.} \end{cases} $$ For $i=0,1,\dots,n$, setting $$ z_i = \frac{u_i}{c(t_i)} \quad \text{with } c(t_i)> 0. $$ We have \[ \int_{\Omega}( \frac{\varepsilon f_i }{c(t_i)} +\frac{c(t_{i-1})}{c(t_i)}z_{i-1}-z_i)( w -z_i )dx + \int_{\Gamma_{N}}\frac{g_i-\lambda Z_i}{c(t_i)} (w-z_i )dx \leq 0 , \] for $w \in K(1)$, which is equivalent to saying \begin{equation} \label{Ueq8} z_i + \partial G_{\frac{g}{c(t_i)}} (z_i )\ni \frac{\varepsilon f_i}{c(t_i)} + \frac{c(t_{i-1})}{c(t_i)}z_{{i-1}}. \end{equation} Now, we introduce the Euler implicit discretization in time associated with \eqref{eP1}, \begin{equation}\label{Ueq6} v_i + \partial G_{\frac{g}{c(t_i)}} (v_i ) \ni \varepsilon h_i + v_{i-1}, \end{equation} for $i=0,1,\dots,n$ with $$ G_{\frac{g}{c(t_i)}}(z):= \begin{cases} -\int_{\Gamma_{N}} \frac{g-\lambda Z}{c(t_i)} zds & \text{if } z\in K(1)\\ +\infty &\text{otherwise}. \end{cases} $$ Let $ v_{\varepsilon}$, defining by $$ v_{\varepsilon}:= \begin{cases} v_0 \quad \text{for } t \in (0, t_1] \\ v_i \quad \text{for } t \in (t_{i-1}, t_i ], \; i=1,\dots,n. \end{cases} $$ Thanks to the nonlinear semigroup theory (see for instance \cite[Theorem 4.6]{Bre}), it follows that $$ v_{\varepsilon} \to v \in \mathcal{C}([0,T); L^2(\Omega)) \quad \text{as } \varepsilon \to 0, $$ where $v$ is a solution of \eqref{eP1}. Problem \eqref{Ueq6} is equivalent to saying \begin{equation} \label{Ueq9} v_i + \partial G_{\frac{g}{c(t_i)}} (v_i ) \ni \frac{\varepsilon f_i}{c(t_i)}-\frac{c(t_i)-c(t_{i-1})}{c(t_i)}(v(t_{i-1})-v_{i-1}) +\frac{c(t_{i-1})}{c(t_i)}v_{{i-1}} , \end{equation} for $i=1,\dots ,n$. We take $z_i$ as a test function in \eqref{Ueq9} to obtain \begin{equation} \label{Ueq10} \begin{aligned} &\int_{\Omega}\Big( \frac{\varepsilon f_i}{c(t_i)} -\frac{c(t_i)-c(t_{i-1})}{c(t_i)}(v (t_{i-1})-v_{i-1}) + \frac{c(t_{i-1})}{c(t_i)}v_{i-1}-v_i\Big)( z_i -v_i )dx \\ &+ \int_{\Gamma_{N}}\frac{g_i-\lambda Z_i}{c(t_i)} (z_i-v_i )dx \leq 0. \end{aligned} \end{equation} Now, we take $ v_i$ as a test function in \eqref{Ueq8} to get \begin{equation} \label{Ueq11} \int_{\Omega}( \frac{\varepsilon f_i}{c(t_i)} +\frac{c(t_{i-1})}{c(t_i)}z_{i-1}-z_i)( v_i -z_i )dx + \int_{\Gamma_{N}}\frac{g_i-\lambda Z_i}{c(t_i)} (v_i-z_i )dx \leq 0 . \end{equation} Combining \eqref{Ueq10} and \eqref{Ueq11}, we obtain \begin{align*} &\int_{\Omega}\Big(-\frac{c(t_i)-c(t_{i-1})}{c(t_i)}(v(t_{i-1})-v_{i-1}) +\frac{c(t_{i-1})}{c(t_i)}(v_{i-1}-z_{i-1}) \\ &+ (z_i -v_i)\Big)( z_i -v_i )dx\leq 0 , \end{align*} which implies \begin{equation} \label{Ueq12} \begin{aligned} &\int_{\Omega}\Big(-\frac{c(t_i)-c(t_{i-1})}{c(t_i)}(v(t_{i-1})-v_{i-1}) +\frac{c(t_{i-1})}{c(t_i)}(v_{i-1}-z_{i-1})\Big)(z_i -v_i)dx \\ &\quad + \int_{\Omega} ( z_i -v_i)^2dx\leq 0. \end{aligned} \end{equation} So, using the H\"{o}lder inequality in \eqref{Ueq12}, we deduce that \begin{align*} &\| z_i -v_i \|^2_{L^2(\Omega)} \\ & \leq \int_{\Omega}(\frac{c(t_i)-c(t_{i-1})}{c(t_i)}(v(t_{i-1})-v_{i-1}) -\frac{c(t_{i-1})}{c(t_i)}(v_{i-1}-z_{i-1}))( z_i -v_i )dx\\ &\leq \big\|\frac{c(t_i)-c(t_{i-1})}{c(t_i)}(v(t_{i-1})-v_{i-1}) -\frac{c(t_{i-1})}{c(t_i)}(v_{i-1}-z_{i-1})\big\|_{L^2(\Omega)} \| z_i -v_i \|_{L^2(\Omega)}. \end{align*} It follows that \begin{equation} \label{Ueq13} \begin{aligned} \| z_i -v_i \|_{L^2(\Omega)} &\leq \big| \frac{c(t_{i-1})}{c(t_i)} \big|\, \| v_{i-1}-z_{i-1} \|_{L^2(\Omega)} +\big| \frac{c(t_i)-c(t_{i-1})}{c(t_i)}\big|\\ &\quad\times \| v(t_{i-1})-v_{i-1} \|_{L^2(\Omega)}, \end{aligned} \end{equation} for $ i=1,\dots,n$. Since $ v_{0}=z_{0}$, then iterating \eqref{Ueq13} for $i=k,\dots,1$, we obtain \begin{align*} \| v_{k}-z_{k} \|_{L^2(\Omega)} &\leq \sum_{i=1}^{k} \big| \frac{c(t_{k-i+1})-c(t_{k-i})}{c(t_{k})}\big| \, \| v(t_{k-i})-v_{k-i} \|_{L^2(\Omega)}\\ &\leq \frac{1}{\delta} \sum_{i=1}^{k} | c(t_{k-i+1})-c(t_{k-i})|\, \| v(t_{k-i})-v_{k-i} \|_{L^2(\Omega)}\\ &\leq \frac{\| c' \|_{\infty} }{\delta} \sum_{i=1}^{k-1}(t_{k-i+1}-t_{k-i}) \| v(t_{k-i})-v_{k-i} \|_{L^2(\Omega)}. \end{align*} Setting $\widetilde{v}_{\varepsilon}(t)=v(t_i)$ if $t\in [t_i, t_{i+1} )$ for $ i=0,1,\dots,n-1$, we get \begin{align*} \| v_{k}-z_{k} \|_{L^2(\Omega)} &\leq \frac{\| c' \|_{\infty} }{\delta} \sum_{i=0}^{k-1}\int_{t_i}^{t_{i+1}} \| \widetilde{v}_{\varepsilon}(t) -v_{\varepsilon}(t) \|_{L^2(\Omega)} dt \\ &\leq \frac{\| c' \|_{\infty} }{\delta}\int_{0}^{T} \| \widetilde{v}_{\varepsilon}(t) -v_{\varepsilon}(t) \|_{L^2(\Omega)}dt \end{align*} and $$ \| v_{\varepsilon}(t) -z_{\varepsilon}(t) \|_{L^2(\Omega)} \leq \frac{\| c' \|_{\infty} }{\delta}\int_{0}^{T} \| \widetilde{v}_{\varepsilon}(t) -v_{\varepsilon}(t) \|_{L^2(\Omega)}dt, \quad \text{for any } t\in [0,T). $$ Since $\varepsilon \to 0$ and $\widetilde{v}_{\varepsilon} \to v$ in $\mathcal{C}([0,T); L^2(\Omega))$, we have \begin{equation} \label{Ueq14} \lim_{\varepsilon\to 0}\sup_{t\in [0,T)} \| v_{\varepsilon}(t) -z_{\varepsilon}(t) \|_{L^2(\Omega)}=0. \end{equation} We define $ c_{\varepsilon}$ by $ c_{\varepsilon}(t)=c(t_i )$ for $t\in [t_i, t_{i+1} )$ and $ i=0,1,\dots,n,$. We have \begin{align*} &\| u(t) -u_{\varepsilon}(t) \|_{L^2(\Omega)} \\ &\leq c_{\varepsilon}(t) \big\| \frac{u(t)}{c_{\varepsilon}(t)} -z_{\varepsilon}(t) \big\|_{L^2(\Omega)} \\ &\leq \| c \|_{\infty} ( \big\| \frac{u(t)}{c_{\varepsilon}(t)} - v(t) \|_{L^2(\Omega)} + \| v(t) -v_{\varepsilon}(t) \|_{L^2(\Omega)} + \| v_{\varepsilon}(t) - z_{\varepsilon}(t)\big\|_{L^2(\Omega)} ). \end{align*} Combining \eqref{Ueq14}, the fact that $v_{\varepsilon} \to v$ in $\mathcal{C}([0,T); L^2(\Omega))$ and $c_{\varepsilon} \to c$ in $\mathcal{C}([0,T))$, and as $\varepsilon\to 0$, we deduce that \begin{equation} \lim_{\varepsilon\to 0}\sup_{t\in [0,T)} \| u(t) -u_{\varepsilon}(t) \|_{L^2(\Omega)}=0. \end{equation} So, $u_{\varepsilon} \to u$ in $\mathcal{C}([0,T); L^2(\Omega))$ as $\varepsilon\to 0$. At last, for the contraction property, if for $ i = 1, 2$, $u_i$ is the solution corresponding to $ f_i$, then we have \begin{align*} &\frac{d}{dt}\int_{\Omega}(u_1(t)-u_2(t))^{+}dx \\ &=\frac{d}{dt}\int_{[u_1(t)>u_2(t)]}(u_1(t)-u_2(t))dx\\ &=\int_{[u_1(t)>u_2(t)]}(f_1(t)-f_2(t))dx\\ &\quad - \int_{[u_1(t)>u_2(t)]}\operatorname{div}(m_1(t)\nabla u_1(t)) -\operatorname{div}(m_2(t)\nabla u_2(t)))dx\\ &\leq \int_{\Omega}(f_1-f_2)^{+}dx \\ &\quad - \int_{\Omega}\text{sign}_{0}(u_1-u_2)\operatorname{div} (m_1(t)\nabla u_1(t))-\operatorname{div}(m_2(t)\nabla u_2(t)))dx\\ &\leq \int_{\Omega}(f_1-f_2)^{+}dx\\ &\quad -\int_{\Gamma_{N}}\text{sign}_{0}(u_1-u_2)(m_1(t)\nabla u_1(t)-m_2(t) \nabla u_2(t)).\nu d\sigma\\ &\leq \int_{\Omega}(f_1-f_2)^{+}dx \end{align*} \end{proof} \section{Dual formulation and numerical approximation} Following the ideas developed in \cite{Ig1, Ig2, Ig3}, we show in this section how to approximate the solution of problem \eqref{eP}. In \cite{Ig3}, the authors introduced a numerical method based on Fenchel duality theory (see, \cite{Ek}) for numerical approximation of problem \eqref{eS}, in the special case where the homogeneous Dirichlet boundary condition is applied. In fact, in the homogeneous Dirichlet boundary condition case, it is not difficult to see that \begin{equation} \label{Ueq15a} v= \mathbb{P}_{K(d)}(\tilde{f}), \end{equation} where $ \mathbb{P}_{K(d)}$ denotes the projection with respect to the $ L^2$ norm onto the convex $ K(d)$ and, thanks to \cite{Ek}, the associated dual problem is given by \begin{equation} \label{Ueq15} \sup\{-G(q): q\in (\mathcal{C}(\Omega))^{N}\}, \end{equation} where \begin{equation} \label{Ueq16} G(q)=\frac{1}{2}\int_{\Omega} (\operatorname{div}(q))^2dx + \int_{\Omega} \widetilde{f} \operatorname{div}(q) dx + d\int_{\Omega} |q| dx. \end{equation} In \cite{Ig3}, the authors have given the connection between \eqref{Ueq15a} and \eqref{Ueq15} by studying $ G$ in the space $ H_{\operatorname{div}}(\Omega)$, defined by $$ H_{\rm div} (\Omega):= \{w\in (L^2(\Omega))^{N}; \operatorname{div}(w)\in L^2(\Omega) \}, $$ which has allowed them to compute the projection $ v$ on $ K(d)$ of $ \tilde{f}$ by computing the solution of the dual problem \eqref{Ueq15}. Since in the present work, we consider a mixed boundary condition, we cannot see problem \eqref{eS} as a projection problem on a convex set. Furthermore, we have the following result. \begin{lemma}\label{lem3} Let $ v\in K(d)$. $ v$ is a solution of \eqref{eS} if and only if $ v$ is also a solution of the following minimization problem \begin{equation}\label{Ueq17} \min_{z\in K(d)}\big\{\frac{1}{2} \int_{\Omega} | z-\widetilde{f}|^2dx -\int_{\Gamma_{N}}(\widetilde{g}-\lambda Z)zds \big\}. \end{equation} \end{lemma} \begin{proof} Let $ v\in K(d)$ be a solution of \eqref{eS}. Then, we have \begin{equation}\label{Ueq18} \int_{\Omega} (\widetilde{f}-v)(z-v)dx + \int_{\Gamma_{N}} (\widetilde{g}-\lambda Z)(z-v) ds \leq 0, \quad \forall z\in K(d). \end{equation} Since $$ \| v-\widetilde{f} \|^2_{L^2(\Omega)} - \| z-\widetilde{f} \|^2_{L^2(\Omega)} =2 (\widetilde{f}-v,z-v) - \| v-z \|^2_{L^2(\Omega)}, $$ we have $$ \frac{1}{2} \int_{\Omega} | v-\widetilde{f} |^2dx - \frac{1}{2} \int_{\Omega} | z-\widetilde{f}|^2dx = \int_{\Omega} (\widetilde{f}-v,z-v) dx - \frac{1}{2} \int_{\Omega} | v-z |^2dx. $$ Thus, \begin{align*} &\frac{1}{2} \int_{\Omega} | v-\widetilde{f} |^2dx - \frac{1}{2} \int_{\Omega} | z-\widetilde{f} |^2dx + \int_{\Gamma_{N}} (\widetilde{g}-\lambda Z)(z-v) ds\\ &=\int_{\Omega} (\widetilde{f}-v,z-v) dx + \int_{\Gamma_{N}} (\widetilde{g}-\lambda Z)(z-v) ds - \frac{1}{2} \int_{\Omega} | v-z |^2dx. \end{align*} Therefore, by using \eqref{Ueq18}, we deduce that \begin{equation} \frac{1}{2} \int_{\Omega} | v-\widetilde{f} |^2dx - \frac{1}{2} \int_{\Omega} | z-\widetilde{f} |^2dx - \int_{\Gamma_{N}} (\widetilde{g}-\lambda Z)(z-v) ds \leq 0; \end{equation} which is equivalent to saying $$ \frac{1}{2} \int_{\Omega} | v-\widetilde{f} |^2dx -\int_{\Gamma_{N}} (\widetilde{g}-\lambda Z)v ds \leq \frac{1}{2} \int_{\Omega} | z-\widetilde{f} |^2dx - \int_{\Gamma_{N}} (\widetilde{g}-\lambda Z)z ds, $$ for any $z \in K(d)$. So, $ v\in K(d)$ is a solution of \eqref{Ueq17}. Now, we suppose that $v\in K(d)$ is a solution to the minimization problem \eqref{Ueq15a}. Let $z_0 \in K(d)$. Since $K(d)$ is convex, then for any $t\in [0,1)$, $z=(1-t)v + tz_0 \in K(d)$. We have \begin{align*} &\frac{1}{2} \int_{\Omega} | v-\widetilde{f} |^2dx -\int_{\Gamma_{N}} (\widetilde{g}-\lambda Z)v ds \\ &\leq \frac{1}{2} \int_{\Omega} | \widetilde{f}- ((1-t)v + tz_0) |^2dx - \int_{\Gamma_{N}} (\widetilde{g}-\lambda Z)((1-t)v + tz_0) ds \\ &\leq \frac{1}{2} \int_{\Omega} | (\widetilde{f}-v)-t(z_0 -v) |^2dx - \int_{\Gamma_{N}} (\widetilde{g}-\lambda Z)v ds - t\int_{\Gamma_{N}}(\widetilde{g}-\lambda Z)(z_0 - v) ds \\ &\leq \frac{1}{2} \| (\widetilde{f}-v)-t(z_0 -v) \|^2_{L^2(\Omega)} - \int_{\Gamma_{N}} (\widetilde{g}-\lambda Z)v ds - t\int_{\Gamma_{N}}(\widetilde{g}-\lambda Z)(z_0 - v) ds. \end{align*} It follows that \begin{align*} &\frac{1}{2} \int_{\Omega} | v-\widetilde{f} |^2dx -\int_{\Gamma_{N}} (\widetilde{g}-\lambda Z)v ds \\ &\leq \frac{1}{2} \int_{\Omega}| v-\widetilde{f} |^2dx - t(\widetilde{f}-v,z_{0}-v) + \frac{t^2}{2}\int_{\Omega} | z_0 -v |^2dx\\ &\quad - \int_{\Gamma_{N}} (\widetilde{g}-\lambda Z)v ds - t\int_{\Gamma_{N}}(\widetilde{g}-\lambda Z)(z_0 - v) ds. \end{align*} So, \begin{equation} \label{Ueq19} (\widetilde{f}-v,z_{0}-v)+ \int_{\Gamma_{N}}(\widetilde{g}-\lambda Z)(z_0 - v)ds \leq \frac{t}{2}\int_{\Omega} | z_0 -v |^2dx. \end{equation} Hence, by letting $t\to 0$ in \eqref{Ueq19}, we obtain $$ \int_{\Omega}(\widetilde{f}-v)(z_{0}-v)dx + \int_{\Gamma_{N}}(\widetilde{g}-\lambda Z)(z_0 - v)ds \leq 0, $$ for any $z_0 \in \ K(d)$, which implies that \begin{equation} \label{Ueq20} \int_{\Omega}(\widetilde{f}-v)z_{0}dx + \int_{\Gamma_{N}}(\widetilde{g}-\lambda Z)z_0 \, ds \leq \int_{\Omega}(\widetilde{f}-v)vdx + \int_{\Gamma_{N}}(\widetilde{g}-\lambda Z)v \, ds, \end{equation} for any $z_0 \in K(d)$. From \eqref{Ueq20}, we deduce that $v$ is a solution of \eqref{eS}. \end{proof} Inspired by the works of Igbida et al.\ (see \cite{Ig2}), we consider as a dual problem associated with problem \eqref{Ueq17}, the following optimization problem \begin{equation}\label{Ueq22} \sup\{-G(w): w\in H_{{\rm div},\tilde{g}}(\Omega)\}, \end{equation} where \begin{equation}\label{Ueq23} G(w)=\frac{1}{2}\int_{\Omega} (\operatorname{div}(w))^2dx + \int_{\Omega} \widetilde{f} \operatorname{div}(w) dx + d\int_{\Omega} | w | dx \end{equation} and \begin{align*} H_{{\rm div},\widetilde{g}}(\Omega) =\Big\{&w\in H_{\rm div} (\Omega):\int_{\Omega} [-\operatorname{div}(w)]\xi dx = \int_{\Omega}w \nabla\xi dx -\int_{\Gamma_{N}}(\widetilde{g}-\lambda Z)\xi ds, \\ &\forall\xi\in H^{1}_{D}(\Omega) \Big\}. \end{align*} It is not clear that the dual problem associated with problem \eqref{Ueq17} is problem \eqref{Ueq22}. In the following, we prove the connection between problems \eqref{Ueq22} and \eqref{Ueq17}. We also present a method for numerical approximation of the solution of problem \eqref{eS} by computing $ \sup\{-G(w): w\in H_{{\rm div},\widetilde{g}}(\Omega)\}$. We first show the following result. \begin{lemma} \label{lem4} For any $ \widetilde{f}\in L^2(\Omega)$, $ \widetilde{g}\in L^2(\Gamma_N)$, $\omega\in H_{{\rm div},\widetilde{g}}(\Omega)$ and $z\in K(d)$, we have $$ -G(w)\leq J(z), $$ where $J(z)=\frac{ 1}{ 2}\int_{\Omega} | z-\widetilde{f} |^2dx - \int_{\Gamma_N}(\widetilde{g}-\lambda Z)z ds$. \end{lemma} \begin{proof} Let $w\in H_{{\rm div},\widetilde{g}}(\Omega)$ and $z\in K(d)$ be fixed. Since \begin{equation} \frac{1}{2}(\operatorname{div}(w)-(z-\widetilde{f}))^2\geq 0, \end{equation} we deduce that $$ -\frac{1}{2}\int_{\Omega} (\operatorname{div}(w))^2 dx - \int_{\Omega} \operatorname{div}(w)\widetilde{f} dx + \int_{\Omega} \operatorname{div}(w)z dx \leq \frac{1}{2}\int_{\Omega} (z-\widetilde{f})^2 dx, $$ which implies $$ -\frac{1}{2}\int_{\Omega} (\operatorname{div}(w))^2 dx - \int_{\Omega} \operatorname{div}(w)\widetilde{f} dx \leq \frac{1}{2}\int_{\Omega} (z-\widetilde{f})^2 dx - \int_{\Omega} \operatorname{div}(w)z dx. $$ Using the fact that $w\in H_{{\rm div},\widetilde{g}} (\Omega)$, we get $$ -\frac{1}{2}\int_{\Omega} (\operatorname{div}(w))^2 dx - \int_{\Omega} \operatorname{div}(w)\widetilde{f} dx \leq \frac{1}{2}\int_{\Omega} (z-\widetilde{f})^2 dx + \int_{\Omega} w\cdot\nabla z dx - \int_{\Gamma_N} (\widetilde{g}-\lambda Z)z\,ds. $$ Therefore, \begin{equation} -G(w)\leq \frac{1}{2}\int_{\Omega} (z-\widetilde{f})^2 dx - \int_{\Gamma_N} (\widetilde{g}-\lambda Z)zds + \int_{\Omega} w\cdot \nabla z dx - d\int_{\Omega} | w | dx. \end{equation} Since $z\in K(d)$, we have \begin{equation} \int_{\Omega} w\cdot \nabla z dx - d\int_{\Omega} | w | dx \leq 0. \end{equation} Thus, $$ -G(w)\leq \frac{1}{2}\int_{\Omega} (z-\widetilde{f})^2 dx - \int_{\Gamma_N} (\widetilde{g}-\lambda Z)zds $$ \end{proof} Our main result of this section reads as follows. \begin{theorem} \label{theo1} Let $ \widetilde{f}\in L^2(\Omega)$, $ \widetilde{g}\in L^2(\Gamma_N)$ and $v$ the variational solution of \eqref{Ueq17}. Then, there exists a sequence $ (w_{\varepsilon})_{\varepsilon >0}$ in $ H_{{\rm div},\widetilde{g}}(\Omega)$ such that, as $\varepsilon\to 0$, \begin{gather} \label{Ueq24} d\int_{\Omega} | \omega_{\varepsilon}| dx \to \int_{\Omega} v(\widetilde{f}-v)dx + \int_{\Gamma_N} (\widetilde{g}-\lambda Z)v\,ds, \\ \label{Ueq25} \operatorname{div}(\omega_{\varepsilon}) \to v-\widetilde{f} \in \ L^2(\Omega) \end{gather} and \begin{equation} \label{Ueq26} \begin{aligned} \lim_{\varepsilon\to 0} G(\omega_{\varepsilon}) &= \inf_{\omega\in H_{{\rm div},g}(\Omega)} G(\omega) \\ &= -\min_{z\in K(d)} J(z)\\ &= -\Big[\frac{1}{2}\int_{\Omega} | \widetilde{f}-v |^2 dx - \int_{\Gamma_N}(\widetilde{g}-\lambda Z)v\,ds\Big] . \end{aligned} \end{equation} \end{theorem} To prove the above result, we consider the following elliptic problem. \begin{equation} \label{eSe} \begin{gathered} v_{\varepsilon} - \nabla w_{\varepsilon} = \widetilde{f} \quad\text{in } \Omega\\ w_{\varepsilon} = \phi_{\varepsilon} (\nabla v_{\varepsilon}) \quad\text{in } \Omega\\ v_{\varepsilon} = 0 \quad\text{on } \Gamma_D \\ w_{\varepsilon} \cdot \eta = \widetilde{g}-\lambda Z \quad\text{on } \Gamma_N, \end{gathered} \end{equation} where for any $ \varepsilon >0$, $\phi_{\varepsilon}:\mathbb{R}^s \to \mathbb{R}^s$ is given by $$ \phi_{\varepsilon}(r)=\frac{1}{\varepsilon} (| r | -d)^{+}\frac{r}{| r |} \quad \text{for all } r \in\mathbb{R}^{s} $$ and satisfies the following properties. \begin{itemize} \item[(i)] for any $r_1$, $r_2 \in\mathbb{R}^s $, $(\phi_{\varepsilon}(r_1 ) - \phi_{\varepsilon}(r_2)).(r_1 - r_2)\geq 0$; \item[(ii)] there exists $\varepsilon_{0}>0$ and $A>d$ such that $\phi_{\varepsilon}(r).r\geq | r |^2$ for any $| r | > A$ and $\varepsilon <\varepsilon_{0}$; \item[(iii)] for any $\varepsilon >0$ and $r\in \mathbb{R}^s$, $d| \phi_{\varepsilon}(r)| \leq \phi_{\varepsilon}(r)\cdot r$. \end{itemize} Then, from the same method as in the proof of \cite[Lemma 3.6]{Ig2}, we can easily prove the following result. \begin{lemma} \label{lem5} There exists a unique weak solution $v_{\varepsilon}$ to problem \eqref{eSe} with $0<\varepsilon<\varepsilon_0$ in the sense that $v_{\varepsilon}\in H^{1}_{D}(\Omega)$, $w_{\varepsilon}=\phi_{\varepsilon} (\nabla v_{\varepsilon})\in (L^2_{D}(\Omega))^{s}$ and for all $z\in H^{1}_{D}(\Omega)$, \begin{equation} \label{Ueq27} \int_{\Omega} v_{\varepsilon} z dx + \int_{\Omega} \phi_{\varepsilon}(\nabla v_{\varepsilon}) \nabla z \,dx =\int_{\Omega} \widetilde{f}zdx + \int_{\Gamma_N} (\widetilde{g}-\lambda Z)z\,ds. \end{equation} Moreover, $(v_{\varepsilon})_{0<\varepsilon<\varepsilon_{0}}$ is bounded in $H^{1}_{D}(\Omega)$ and for any Borel set $B\subset\Omega$, we have \begin{equation} \label{Ueq28} \lim\inf_{\varepsilon\to 0}\int_{B} | \nabla v_{\varepsilon} | dx \leq d| B |. \end{equation} \end{lemma} \begin{proof}[Proof of the Theorem \ref{theo1}] Thanks to Lemma \ref{lem5}, one has that the sequence $(v_{\varepsilon})_{\varepsilon>0}$ is bounded in $ H^{1}_{D}(\Omega)$. Therefore, we can extract a subsequence (still denoted by $ (v_{\varepsilon})_{\varepsilon>0}$) such that $$ v_{\varepsilon}\to\widetilde{v} \quad\text{in } H^{1}_{D}(\Omega)\text{-weak and in } L^2(\Omega). $$ Hence, by problem \eqref{eSe}, we deduce that $$ \operatorname{div}(w_{\varepsilon})\to \widetilde{v}-\widetilde{f} \quad \text{in } L^2(\Omega). $$ We define the set $A_{\delta}$ by $A_{\delta}=[| \nabla \widetilde{v} |\geq d+\delta]$, with $\delta>0$. Using the fact that $\nabla v_{\varepsilon} \to \nabla \widetilde{v}$ in $(L^{1}(\Omega))^s$-weak as $\varepsilon\to 0$, it follows that \begin{equation} (d+\delta)| A_{\delta}| \leq \int_{A_{\delta}} | \nabla\widetilde{v} | dx \leq \lim\inf_{\varepsilon \to 0}\int_{A_{\delta}} | \nabla v_{\varepsilon} | dx. \end{equation} From Lemma \ref{lem4}, one has $ B= A_{\delta}$ and then \begin{equation} (d+\delta)| A_{\delta}|\leq | A_{\delta}|. \end{equation} Therefore $| A_{\delta}|=0$ since $\delta>0$. So, $| \nabla\widetilde{v} | \leq d$ a.e. in $\Omega$ and $\widetilde{v}\in K(d)$. Now we must show that the solution $\widetilde{v}$ is also a solution to problem \eqref{Ueq17}. For any $ z\in K(d)$, one has \begin{equation}\label{Ueq37} \begin{aligned} &\int_{\Omega} (\widetilde{f}-\widetilde{v})(z-\widetilde{v})dx + \int_{\Gamma_N} (\widetilde{g}-\lambda Z)(z-\widetilde{v})ds \\ & = \lim_{\varepsilon\to 0} \int_{\Omega}-\nabla\phi_{\varepsilon}(\nabla v_{\varepsilon})(z-\widetilde{v})dx + \int_{\Gamma_N} (\widetilde{g}-\lambda Z)(z-\widetilde{v})ds \\ &= \lim_{\varepsilon\to 0} \int_{\Omega}\phi_{\varepsilon}(\nabla v_{\varepsilon}) \nabla(z-\widetilde{v})dx-\int_{\Gamma_N} (\widetilde{g}-\lambda Z) (z-\widetilde{v})ds \\ &\quad + \int_{\Gamma_N} (\widetilde{g}-\lambda Z)(z-\widetilde{v})ds \\ &= \lim_{\varepsilon\to 0} \int_{\Omega}\phi_{\varepsilon}(\nabla v_{\varepsilon}) \nabla(z-\widetilde{v})dx \\ &= \lim_{\varepsilon\to 0}\int_{\Omega}(\phi_{\varepsilon}(\nabla v_{\varepsilon}) - \phi_{\varepsilon} (\nabla z))\nabla(z-\widetilde{v})dx \leq 0, \end{aligned} \end{equation} which proves that $ \widetilde{v}$ is a also an variational solution of \eqref{eS}. Then, it follows that $ v= \widetilde{v}$ and $ \widetilde{v}$ also satisfies \eqref{Ueq17}. Now, we must show that $\omega_{\varepsilon}$ satisfies \eqref{Ueq24}. From the property (iii) of $\phi_{\varepsilon}$, it follows that \begin{equation}\label{Ueq38} \begin{aligned} \limsup_{\varepsilon\to 0}d\int_{\Omega} | \omega_{\varepsilon}| dx &= \lim\sup_{\varepsilon\to 0}\int_{\Omega}d| \phi_{\varepsilon}(\nabla v_{\varepsilon}) | dx \\ &\leq \limsup_{\varepsilon\to 0}\int_{\Omega} \phi_{\varepsilon}(\nabla v_{\varepsilon}) \nabla v_{\varepsilon}\,dx \\ &\leq \limsup_{\varepsilon\to 0} \Big( \int_{\Gamma_N}\phi_{\varepsilon} (\nabla v_{\varepsilon}).\eta v_{\varepsilon}dx - \int_{\Omega}\nabla\phi_{\varepsilon}(\nabla v_{\varepsilon})v_{\varepsilon} dx \Big) \\ &\leq \lim\sup_{\varepsilon\to 0}\Big(\int_{\Omega} (\widetilde{f}-v_{\varepsilon}) v_{\varepsilon} dx + \int_{\Gamma_N}(\widetilde{g}-\lambda Z)v_{\varepsilon}\, ds \Big) \\ &\leq \int_{\Omega} (\widetilde{f}-\widetilde{v})\widetilde{v} dx + \int_{\Gamma_N} (\widetilde{g}-\lambda Z)\widetilde{v} ds. \end{aligned} \end{equation} and since $ v$ is a variational solution to problem \eqref{eS}, \eqref{Ueq38} becomes \begin{equation} \label{Ueq39} d\limsup_{\varepsilon\to 0}\int_{\Omega} | w_{\varepsilon}| dx \leq \int_{\Omega} (\widetilde{f}-v)v dx +\int_{\Gamma_N}(\widetilde{g}-\lambda Z )v ds. \end{equation} Using the fact that $\widetilde{v}$ is a variational solution to problem \eqref{eS}, it follows that \begin{equation} \label{Ueq40} \begin{aligned} \int_{\Omega} (\widetilde{f}-v)v dx +\int_{\Gamma_N}(\widetilde{g}-\lambda Z)v ds &\leq \int_{\Omega} (\widetilde{f}-\widetilde{v})\widetilde{v} dx + \int_{\Gamma_N} (\widetilde{g}-\lambda Z)\widetilde{v} ds \\ &= \lim_{\varepsilon\to 0} \int_{\Omega} (\widetilde{f}-v_{\varepsilon}) \widetilde{v} dx + \int_{\Gamma_N} ( \widetilde{g}-\lambda Z)\widetilde{v} ds \\ &=\lim_{\varepsilon\to 0}\int_{\Omega} -\nabla w_{\varepsilon}.\widetilde{v} dx + \int_{\Gamma_N}(\widetilde{g}-\lambda Z)\widetilde{v} ds \\ &=\lim_{\varepsilon\to 0}\int_{\Omega} w_{\varepsilon}\cdot\nabla\widetilde{v} dx \\ &\leq d\lim_{\varepsilon\to 0}\int_{\Omega} | w_{\varepsilon} | dx. \end{aligned} \end{equation} From \eqref{Ueq38} and \eqref{Ueq40}, it follows that \begin{equation} \label{Ueq41} d\lim_{\varepsilon\to 0 } \int_{\Omega} | w_{\varepsilon} | dx =\int_{\Omega} (\widetilde{f}-v)v dx +\int_{\Gamma_N}(\widetilde{g}-\lambda Z)v ds. \end{equation} To complete the proof, it remains to prove \eqref{Ueq26}. For that, we use \eqref{Ueq24}, \eqref{Ueq26} and \eqref{Ueq41} to get \begin{align*} \lim_{\varepsilon\to 0} (-G(\omega_{\varepsilon})) &= \lim_{\varepsilon\to 0} \Big( -\frac{1}{2} \int_{\Omega} [\operatorname{div}(w_{\varepsilon})]^2 dx -\int_{\Omega} \operatorname{div}(\omega_{\varepsilon})\widetilde{f}dx -d\int_{\Omega} | w_{\varepsilon} | dx \Big)\\ &= -\frac{1}{2}\int_{\Omega} (\widetilde{v}-\widetilde{f})^2 dx -\int_{\Omega} (\widetilde{v}- \widetilde{f})\widetilde{f} dx - \int_{\Omega} (\widetilde{f}-\widetilde{v})\widetilde{v} dx\\ & \quad - \int_{\Gamma_N}(\widetilde{g}-\lambda Z)\widetilde{v} ds\\ &= -\frac{1}{2} \int_{\Omega} (\widetilde{v}-\widetilde{f})^2 dx +\int_{\Omega} (\widetilde{v}-\widetilde{f})^2 dx - \int_{\Gamma_N} (\widetilde{g}-\lambda Z)\widetilde{v} ds\\ &= \frac{1}{2} \int_{\Omega} (\widetilde{v}-\widetilde{f})^2 dx - \int_{\Gamma_N}(\widetilde{g}-\lambda Z)\widetilde{v} ds. \end{align*} Therefore, \begin{equation}\label{Ueq42} \lim_{\varepsilon\to 0}(-G(w_{\varepsilon})) = J(\widetilde{v})= J(v). \end{equation} Lemma \ref{lem4} allows us to write \begin{equation}\label{Ueq43} \sup_{\omega\in H_{{\rm div},\widetilde{g}}(\Omega)} (-G(w))\leq J(v) =\lim_{\varepsilon\to 0}(-G(w_{\varepsilon})). \end{equation} Since \begin{equation} \label{Ueq44} J(v)=\lim_{\varepsilon\to 0}(-G(w_{\varepsilon})) \leq \sup_{\omega\in H_{{\rm div},\widetilde{g}}(\Omega)} (-G(w)), \end{equation} one concludes that \begin{equation}\label{Ueq45} \lim_{\varepsilon\to 0}(-G(w_{\varepsilon})) =\sup_{\omega\in H_{{\rm div},\widetilde{g}}(\Omega)} (-G(w))=J(v) \end{equation} \end{proof} \begin{remark} \label{rem1} \rm As consequences of Theorem \ref{theo1}, we have the following characterization of the solution $v$ to problem \eqref{eS}: \begin{equation} \label{Ueq46} \begin{gathered} v-\operatorname{div}(\omega)=\widetilde{f} \quad\text{in } (H^{1}_{D}(\Omega))^{*} \\ d| \omega |(\Omega) = \int_{\Omega}v(\widetilde{f}-v)dx + \int_{\Gamma_N }(\widetilde{g}-\lambda Z)v\,ds, \end{gathered} \end{equation} where $ w\in (\mathcal{M}_{b}(\Omega))^s$ is the weak limit of $ (w_{\varepsilon})_{\varepsilon>0}$, in $(\mathcal{M}_{b}(\Omega))^s$. \end{remark} For numerical tests, we need to construct a subset of $ H_{{\rm div},\widetilde{g}}(\Omega)$. For that, we introduce the following result. \begin{lemma} \label{lem6} $\sup\{-G(y): y\in H_{{\rm div},\widetilde{g}}\} = \sup\{\mathcal{G}(\varphi):\varphi\in\mathcal{H}\}$, where $\mathcal{G}(\varphi)=G(\varphi+\nabla Z)$ with $Z$ the variational solution to the Laplace problem \eqref{eL} and $$ \mathcal{H}=\Big\{\varphi\in H_{\rm div}(\Omega): -\int_{\Omega} \operatorname{div}(\varphi)\xi =\int_{\Omega} \varphi\cdot\nabla \xi, \ \forall \xi\in H^{1}_{D}(\Omega) \Big\}. $$ \end{lemma} This lemma can be proved, by the same method given in \cite{Ig2}. Consequently, the dual problem \eqref{Ueq22} is equivalent to \begin{equation}\label{Ueq47} \inf \{G(\phi + \nabla Z): \phi \in\mathcal{H}\}. \end{equation} To approximate the numerical minimization problem \eqref{Ueq47}, we use the finite element method and we make the following assumptions. \begin{itemize} \item $ \Omega$ is an open and bounded subset of $ \mathbb{R}^2$. \item $ T_h$ is a regular partition (quadrangulation) of $ \overline{\Omega}$ by $N$ disjoint open simplex $ \tau$ of diameter no greater than a given real $ h$, with $ \overline{\Omega}=\cup_{\tau\in T_h} \overline{\tau}$. \end{itemize} Let $ V_h$ be a finite dimensional subspace of $ RT_{0}(T_h)\cap\mathcal{H}$ with dimensional equal to $N=N(h)$, with $ RT_{0}(T_h)$ the space of lowest-order Raviart-Thomas finite elements (see \cite{JE}) defined by \begin{align*} RT_{0}(T_h)=\Big\{&q_{h} \in (L^2(\Omega))^2 : q_{\tau}^{h}=a_{\tau} + b_{\tau}x, \; a\in\mathbb{R}^2 , b\in\mathbb{R}, \; \forall \tau \in T_h,\\ &\text{and } (q_{\tau}^{h} -q_{\tau '}^{h}).\nu_{\partial_{\tau}} = 0 \text{ on } \partial_{\tau} \cap \partial_{\tau '} \Big\}, \end{align*} 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 $ V_h$ given in \cite[Theorem 6.1]{JE}. Then, thanks to \cite{JE}, for all $ w\in H_{\rm div}(\Omega)$, we have \begin{equation}\label{Ueq48} r_{h} (w)\to w \text{ in } (L^2(\Omega))^2 \text{ and } \operatorname{div}(r_{h} (w)) \to \operatorname{div}(\omega) \text{ in } L^2(\Omega), \text{ as } h\to 0. \end{equation} We have the following convergence properties as $h$ tends to $0$. \begin{theorem} \label{theo2} Let $ \widetilde{f}\in L^2(\Omega)$, $ \widetilde{g}\in L^2(\Gamma_N)$, $ v$ a solution to the minimization problem \eqref{Ueq17} and $ w_{h}$ a solution of the optimization problem \begin{equation} \label{Ueq49} \sup\{- \mathcal{G} ( q_h ) : q_h \in V_h \}. \end{equation} Then, as $h\to 0$, \begin{equation}\label{Ueq50} \operatorname{div}(\omega_{h}) \to v-\widetilde{f} \ in \ L^2(\Omega) \end{equation} and \begin{equation}\label{Ueq51} - \mathcal{G} (\omega_{h}) \to \min_{z\in K} J(z) =\frac{1}{2} \int_{\Omega} | v-\widetilde{f}|^2 dx - \int_{\Gamma_N} (\widetilde{g}-\lambda Z)v\,ds. \end{equation} \end{theorem} The proof of this result is similar to the proof of \cite[Theorem 3.8]{Ig2}; therefore, we omit it. \section{Numerical simulations} We start by approximating the term $\int_{\tau}|q_{h} + \nabla Z|dx$ on each element $ \tau$ of the partitioning $ T_{h}$. We take $$ \int_{\tau}|q_{h} + \nabla Z\,|dx \sim |\tau||q_{h} + \nabla Z|(P_{\tau}) , $$ where $ |\tau| $ represents the area of simplex $ \tau$ and $ (P_{\tau})$ 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 \eqref{Ueq49} is a minimizer of the non-differentiable functional $\mathcal{G}:\mathbb{R}^{s}\to\mathbb{R}$, \[ \mathcal{G}(w_{h}):=\frac{1}{2}(Aw_{h},w_{h}) + (dtf_{h}^{n}+u_{n-1}, \operatorname{div}(w_{h}))+ c(ndt)\sum_{\tau \in T_{h}}|\tau|\, | w_{h} + \nabla Z |(P_{\tau}). \] The minimization of the functional $\mathcal{G}$ is done according to the Gauss-Seidel type algorithm which is described in the following way. \begin{itemize} \item We start the algorithm by the initial vector $q_{0}\in\mathbb{R}^{s}$, for any $k\geq 0$ until convergence; \item we choose a canonical direction $e_{j}\in\mathbb{R}^{s}$ and we find $\rho_{jk}$ minimizing the functional $\phi_{jk}:\mathbb{R}\to\mathbb{R}$, \[ G_{h}(q_{k}+\rho e_{j}), \] \item we take $q_{k+1}=q_{k}+\omega\rho_{jk}e_{j}$, with $\omega$ over relaxation parameter, \item algorithm is performed until $\| q_{k+1}-q_{k} \|_{L^2(\Omega)} \leq \varepsilon$, $\varepsilon$ is the convergence criterion. Afterwards, take $ w_{h}=q_{k}$. \end{itemize} Knowing a minimizer $ w_{h}$, the solution $ u_{n}$ of the Euler implicit time discretization of \eqref{eP} is computed by using extremality relation \eqref{Ueq46} in a weak sense with piecewise finite elements $ P_{0}$. Here, we study the growing sandpile problem with time evolution that is $c(t)$ which is connected with time. So, we take $$ c(t_{k})=\frac{t_{k}}{T}. $$ For numerical tests, we consider $\Omega=[-1,1]\times [-1,1]$, $\omega=0.5$, $N=60$, $dt=10^{-3}$ and $\varepsilon=10^{-3}$ the convergence criterion. We consider firstly $\Gamma_{N}=\{y\in[-1,1], \ x=-1\}$ that is we suppose that there exists a wall at the boundary $x=-1$. \begin{figure}[ht] \centering \includegraphics[width=0.45\textwidth]{fig1a} % f=0g=1L=10^-60 \includegraphics[width=0.45\textwidth]{fig1b} \\ % f=0g=1L=10^55} $f=0$, $g=1$, $\lambda=10^{-60}$\hfil $f=0$, $g=1$, $\lambda=10^{55}$ \caption{} \label{fig1} \end{figure} In Figure \ref{fig1}(left), $g$ has the same role as the source. In Figure \ref{fig1}(right), we take $\lambda=10^{55}$ and we see that for this value the boundary value $g-\lambda Z$ is near zero, that is the reason why the height of the sandpile is very small. \begin{figure}[ht] \centering \includegraphics[width=0.45\textwidth]{fig2a} % f=5g=0L=10^-70 \includegraphics[width=0.45\textwidth]{fig2b} \\ % f=5g=0L=10^55 $f=5$, $g=0$, $\lambda=10^{-70}$. \hfil $f=5$, $g=0$, $\lambda=10^{55}$ \caption{} \label{fig2} \end{figure} We remark that the two graphs in Figure \ref{fig2} have the same configuration. This is so because when $g=0$, the solution $Z$ of problem \eqref{eL} is zero regardless of the choice of $\lambda$. \begin{figure}[ht] \centering \includegraphics[width=0.45\textwidth]{fig3a} % f=3g=2L=10^-70 \includegraphics[width=0.45\textwidth]{fig3b} \\ % f=3g=2L=10^30 $f=3$, $g=2$, $\lambda=10^{-70}$\hfil $f=3$, $g=2$, $\lambda=10^{30}$ \caption{} \label{fig3} \end{figure} It can be seen that two graphs in Figure \ref{fig2}, and the one in Figure \ref{fig3}(right) have the same configuration because for $\lambda$ is very large; the boundary value is similar to the case when $g=0$. For the rest, we consider $\Gamma_{N}=\{y\in[-1,1], \, x=-1, \text{ and } x=1\}$, that is we suppose that there exist two opposites sides on $x=-1$ and $x=1$. \begin{figure}[ht] \centering \includegraphics[width=0.45\textwidth]{fig4a} % 2Bf=0g=40L=10^-20 \includegraphics[width=0.45\textwidth]{fig4b} \\ % 2Bf=0g=40L=10^70 $f=0$, $g=40$, $\lambda=10^{-20}$ \hfil $f=0$, $g=40$, $\lambda=10^{70}$ \caption{} \label{fig4} \end{figure} From Figure \ref{fig4}(left) one sees that $g$ has the same role as the source. One also sees that for $\lambda=10^{70}$, the height of sandpile is very small because for this value, the solution $Z$ of \eqref{eL} goes to zero. \begin{figure}[ht] \includegraphics[width=0.45\textwidth]{fig5a} % 2Bf=5g=0L=10^-60 \includegraphics[width=0.45\textwidth]{fig5a} \\ % 2Bf=5g=0L=10^50 $f=5$, $g=0$, $\lambda=10^{-60}$ hfil $f=5$, $g=0$, $\lambda=10^{50}$ \caption{} \label{fig5} \end{figure} The two graphs in Figure \ref{fig5} have the same profile because when $g=0$, the solution $Z$ to problem \eqref{eL} is near to zero. Hence, the formation of sandpile depends only on the source $f$. \begin{figure}[ht] \includegraphics[width=0.45\textwidth]{fig6a} % 2Bf=5g=3L=10^-45 \includegraphics[width=0.45\textwidth]{fig6b} \\ % 2Bf=5g=3L=10^60 $f=5$, $g=3$, $\lambda=10^{-45}$ \hfil $f=5$, $g=3$, $\lambda=10^{60}$ \caption{} \label{fig6} \end{figure} The two graphs in Figure \ref{fig5}, and the one in Figure \ref{fig6}(right) are similar because for $f\neq0$, $g=0$ and $f\neq0$, $g\neq 0$ and $\lambda=10^{60}$ ($\lambda$ very large), boundary value $g-\lambda Z$ is near to zero. \begin{remark} \label{rmk4.1} \rm (1) In the simulations, when the source $f$ is equal to zero and $\lambda$ is very large, it sees that the evolution of the sandpile is very low; i.e., the height of sandpile is near to zero. (2) When the source $f\neq 0$, $g\neq 0$ and $\lambda$ is very large ($\lambda\to \infty$) we get same sandpile configuration as when $g=0$, regardless of the choice of $\lambda$. \end{remark} \begin{thebibliography}{00} \bibitem{Mik} M. A. Artemov, E. S. Bararanovskii; \emph{Mixed boundary-value problems for motion equations of viscoelastic medium.} Electronic Journal of Differential Equations, Vol. 2015 (2015), No. 252, pp. 1-9. \bibitem{Ek} I. Ekeland, R. Temam; \emph{Convex analysis and variational problems.} In Classics in Applied Mathematics, 28. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA ,1999. \bibitem{Ev} L. C. Evans; \emph{Partial Differential Equations.} Grad. Stud. Math. 19, Amer. Math. Soc., 1998. \bibitem{Bre} H. Brezis; \emph{Op\'erateurs maximaux monotones et s\'emi-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, 1973. \bibitem{Hao} R. Haoua, A. Medeghri; \emph{Robin boundary value problems for elliptic operational differential equations with variable coefficients.} Electronic Journal of Differential Equations, Vol. 2015 (2015), No. 87, pp. 1-9. \bibitem{JE} J. E. Robert, J. M. Thomas; \emph{Mixed and Hybrid Methods.} In: P. G. Ciarlet and J. L. Lions, (editors), Handbook of numerical Analysis, Finite Elements Methods (Part 1), Vol. II, North-Holland, Amsterdam (1991). \bibitem{BCRE94} J. P. Bouchaud, M. E. Cater, J. Ravi Prakash, S. E. Edwards; \emph{A model for the dynamics of sandpile surfaces.} J. Phys. FR.1, \textbf{4}, 1383-1410 (1994). \bibitem{LB} L. Prigozhin, J. W. Barrett; \emph{Dual formulation in critical state problems.} Interface Free Bound. \textbf{8}, 349-370 (2006). \bibitem{PV} L. Prigozhin; \emph{Variational model of sandpile growth.} Euro. J. Appl. Math. \textbf{7}, 225-236 (1996). \bibitem{Ig1} N. Igbida; \emph{A generalized collapsing sandpile model.} Arch. Math. (Basel), \textbf{94}, pp. 193-200 (2010). \bibitem{Ig2} N. Igbida, F. Karami, S. Ouaro and U. Traor\'e; \emph{On a dual formulation for growing sandpile problem with mixed boundary conditions.} Appl. Math. Optim. DOI 10.1007/s00245-017-9408-2, 2017. \bibitem{Be} P. B\'enilan, M. G. Crandall, A. Pazy; \emph{Evolution Equations Governed by Accretive Operators}, to appear. \bibitem{Ig3} S. Dumont, N. Igbida; \emph{On the collapsing sandpile problem.} Commun. Pure Appl. Anal., \textbf{10} (2011) no. 2, pp. 625-638. \bibitem{Ig4} S. Dumont, N. Igbida; \emph{On a dual formulation for growing sandpile problem.} European J. Appl. Math. \textbf{20} (2009), no 2, 169-185. \bibitem{SH89} S. B. Savage, K. Hutter; \emph{The motion of a finite mass of granular material down a rough incline.} J. Fluid Mech., \textbf{199} (1989), 177-215. \end{thebibliography} \end{document}