\documentclass[reqno]{amsart}
\AtBeginDocument{{\noindent\small
{\em Electronic Journal of Differential Equations},
Vol. 2001(2001), No.~07, pp. 1--12.\newline
ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.swt.edu \quad ejde.math.unt.edu (login: ftp)}
\thanks{\copyright 2001 Southwest Texas State University.}
\vspace{1cm}}
\begin{document}
\title[\hfilneg EJDE--2001/07\hfil A penalty method for approximations]
{ A penalty method for approximations of the
stationary power-law Stokes problem }
\author[ Lew Lefton \& Dongming Wei \hfil EJDE--2001/07\hfilneg]
{Lew Lefton \& Dongming Wei }
\address{Lew Lefton \hfill\break\indent
School of Mathematics\\
Georgia Institute of Technology \hfill\break\indent
Atlanta, Georgia 30332, USA}
\email{llefton@math.gatech.edu}
\address{Dongming Wei \hfill\break\indent
Department of Mathematics\\
University of New Orleans \hfill\break\indent
New Orleans, Louisiana 70148, USA}
\email{dwei@math.uno.edu}
\date{}
\thanks{Submitted October 31, 2000. Published January 10, 2001.}
\subjclass{65N30, 65N12, 65N15, 35J70 }
\keywords{Power-law flow, penalty method, stationary Stokes problem,
\hfill\break\indent
non-Newtonian flows, convergence and error estimates, LBB condition}
\begin{abstract}
We study approximations of the
steady state Stokes problem governed by the power-law model for
viscous incompressible non-Newtonian flow using the penalty
formulation. We establish convergence and find error estimates.
\end{abstract}
\maketitle
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\theoremstyle{remark}\newtheorem{remark}[theorem]{Remark}
\numberwithin{equation}{section}
\newcommand{\abs}[1]{\lvert#1\rvert}
\newcommand{\bu}{\ensuremath{\mathbf{u}}}
\newcommand{\bv}{\ensuremath{\mathbf{v}}}
\newcommand{\bx}{\ensuremath{\mathbf{x}}}
\newcommand{\by}{\ensuremath{\mathbf{y}}}
\newcommand{\bX}{\ensuremath{\mathbf{X}}}
\newcommand{\bD}{\ensuremath{\mathbf{D}}}
\newcommand{\bsigma}{\mbox{\boldmath\ensuremath{\sigma}}}
\newcommand{\bphi}{\mbox{\boldmath\ensuremath{\phi}}}
\newcommand{\bW}{\ensuremath{\mathbf{W}}}
\newcommand{\sn}[3]{\ensuremath{\|#3\|_{#1,#2}}} % Sobolev norm
\section{Introduction}
Let $\Omega $ be a convex bounded domain in
$\mathbb{R}^d$, $d\ge 2$. We consider the
steady state flow of a fluid in $\Omega$, where
$\bu(\bx)=(u_1(\bx),\dots ,u_d(\bx))$ denotes the velocity of a fluid
particle at $\bx=(x_1,\dots ,x_d)\in \Omega$. Let $\bsigma\in
\mathbb{R}^d\times \mathbb{R}^d$ denote the stress tensor for the fluid.
The momentum equations for an isothermal steady state flow are
\begin{equation}\label{eq:int1}
\rho(\bu\cdot \nabla ) \bu =
\nabla \cdot \mbox {\bsigma} + \mathbf f \quad \mathrm{in}~\Omega,
\end{equation}
where $\rho$ is the density of the fluid, $\mathbf f =(f_1,\dots ,f_d)$ the
body force, and $\nabla = (\frac{\partial} {\partial x_1},\ldots,
\frac{\partial} {\partial x_d})$.
The $j^{th}$ component of $(\bu \cdot \nabla) \bu$ is $
\sum_{i=1}^d {u_i}\frac{\partial u_j} {\partial x_i}$ and
$\nabla \cdot \mbox {\bsigma}$ is obtained by applying the divergence
operator, defined by $\nabla \cdot \bu = \sum_{i=1}^d
\frac{\partial u_i}{\partial x_i}$,
to each row of $\mbox {\bsigma}$. We further assume that the
fluid is incompressible so it satisfies the continuity equation
\begin{equation}\label{eq:int2}
\nabla \cdot \bu = 0 \quad \mathrm{in}~ \Omega.
\end{equation}
The rate of deformation tensor $\bD(\bu)\in \mathbb{R}^d\times\mathbb{R}^d$ is the symmetric gradient of
$\bu$ with components $D_{ij}(\bu) = \frac 1 2 \left(\frac{\partial
u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$. For
incompressible fluids, the second invariant of $\bD(\bu)$ denoted
$\Pi_{\bD}(\bu)$ satisfies \[ -2 \Pi_{\bD}(\bu)= \bD(\bu):\bD(\bu) =
\sum_{i,j=1}^d {D_{ij}(\bu)^2}= \abs{\bD(\bu)}^2, \]
where $|\cdot|$ denotes the Euclidean matrix norm; that is for ${\bf K}$, a
$ d \times d$ real matrix, $|{\bf K}|=[\sum_{i,j=1}^d (k_{ij})^2]^{\frac {1}{2}}.$
In the power-law model for non-Newtonian fluid flows, it is assumed
that viscosity $\eta$ varies as a power of the shear-strain rate, giving a
stress tensor $\mbox {\bsigma}$ of the form
\begin{equation}\label{eq:int3}
\bsigma =-p {\mathbf I} +
\eta(\Pi_{\bD}(\bu)) \bD(\bu), \end{equation}
where $p$ is the scalar pressure,
$\mathbf I\in \mathbb{R}^d\times\mathbb{R}^d$ is the $d$ dimensional identity matrix, and
\[
\eta(z)=\eta_0 |z|^{(r-2)/2}, z \in \mathbb{R}.
\]
Here we assume $1< r < +\infty$, and
$\eta_0 > 0$. Substituting (\ref{eq:int3}) into (\ref{eq:int1}), we
obtain the steady state power-law Navier-Stokes equation
\begin{equation}\label{eq:int4}
- k \nabla \cdot (\abs{\bD(\bu)}^{r-2}
\bD(\bu) )+\rho(\bu\cdot \nabla) \bu+ \nabla p = \mathbf f, \end{equation}
where $k=\eta_0/2^{\frac{r-2}{2}}$.
The \textit{power-law Stokes equation} is obtained by neglecting the
inertial term $(\bu\cdot \nabla )\bu$ in (\ref{eq:int4}). This model
of non-Newtonian flow is very popular in chemical engineering \cite{BCH} as well
as in geophysics \cite{SE}. It has also been used
in applications for the design of extrusion dies \cite{BBM},
\cite{LWT}, and for the study of the lithosphere \cite{EJ}, \cite{EM1},
\cite{EM2}. To make the Stokes problem well posed, we assume that the solution satisfies the continuity
equation (\ref{eq:int2}) and, for simplicity, a homogeneous boundary condition
of Dirichlet type. The resulting problem is
\begin{equation}\label{eq:int5}
\gathered
- k \nabla \cdot(\abs{\bD(\bu)}^{r-2} \bD(\bu) ) +\nabla p = \mathbf f
\quad \mbox{in}~ \Omega,\\
\nabla\cdot \bu = 0 \quad \mbox{in}~ \Omega,\\
\bu = \mathbf 0 \quad \mbox{on}~ \partial \Omega,
\endgathered
\end{equation}
where $\partial \Omega$ denotes the boundary of $\Omega$.
\begin{remark}
If $d \ge 2$ and the flow is unidirectional then $\bu(\bx)=
(u(\bx),0,\dots ,0)$, where $u(\bx)$ is a scalar valued function and ${\mathbf f}=(f_1(\bx), 0, \ldots,0)$. For
$\bu$ to satisfy the continuity equation (\ref{eq:int2}), we have
$\frac{\partial u}{\partial x_1}=0$ which implies that $u(\bx)=
u(x_2,\dots ,x_d) $. Substituting $\bu(\bx)$ into equation (\ref{eq:int5}) and writing it as a system shows that $(\frac{\partial p}{\partial x_2}, \ldots, \frac{\partial p}{\partial x_n}) ={\mathbf 0}$ and $\tilde f \equiv f_1(\bx) - \frac{\partial p}{\partial x_1}$ is independent of $x_1$.
Thus, we are left with
a scalar quasilinear elliptic Dirichlet problem in the $\mathbb{R}^{d-1}$ domain $\tilde \Omega = \Omega \cap \{x_1=0\}$
\[
-{\tilde k} \Delta_r u = {\tilde f} \quad
\mbox{in}~ {\tilde\Omega},\quad \quad u = 0\quad \mbox{on}~ \partial {\tilde\Omega}
\]
where $\Delta_r u = \nabla \cdot(\abs{\nabla u}^{r-2}\nabla u)$ is
the quasilinear generalization of the Laplacian known as the $r$-Laplacian, and
$\nabla u = (\frac{\partial u}{\partial x_2},\dots , \frac{\partial u}{\partial x_d})$.
(We note that $\Delta_p$ is frequently called the $p$-Laplacian in the
literature, but we are using $p$ to denote fluid pressure here.) There has been
a great deal of analytical (e.g.,~\cite{Bh}, \cite{D}, \cite{MNZ}) and
numerical work (e.g.,~\cite{BL3}, \cite {GM}, \cite{LW2}, \cite{W}) devoted to problems
involving $\Delta_r u$.
\end{remark}
\begin{remark} In our notation, the power-law index is the value
$n=r-1$. We recall that a fluid is considered to be Newtonian if it
has power-law index $n = 1$, and non-Newtonian if $n \ne 1$. The
power-law equation (\ref{eq:int3}) is also called the Ostwald-deWaele
equation. When $ 0 < n < 1$, which corresponds to our parameter $ 1 <
r < 2$, power-law fluids exhibit a decrease in viscosity with
increasing shear stress and they are known as pseudoplastic or
shear-thinning fluids. When $ 1 < n < \infty $, corresponding to $ 2 <
r < \infty$, power-law fluids exhibit an increase in viscosity with
increasing shear stress and they are known as dilatant or
shear-thickening fluids \cite{Wi}. For a specific example, we cite \cite{J} where the value of $n$ for a certain tomato paste is given as $n=0.257$.
\end{remark}
One way of studying a stationary power-law Stokes flow is to
consider the velocity field $\bu$ as the minimizer of an appropriate
energy functional. In order to enforce the constraint of a divergence
free flow, we seek a minimizer $\bu$ in the space of divergence free
vector fields, that is, we solve the problem:
\begin{equation}\label{eq:int6}
\gathered
\text{Find the minimum of}~J(\bu)
= \frac{k}{r}\int_\Omega |\bD(\bu )|^ {r}\; d\bx- \int_\Omega
\mathbf f \cdot \bu \;d\bx, \\
\text{where }~\bu\in\bX= \left\{ \bu \in \bW_0^{1,r}(\Omega):
\nabla \cdot \bu =0\right\}.
\endgathered\end{equation}
This is a
minimization problem with a constraint and we refer to it as the
\textit{variational formulation} of power-law flow. The Euler equation
corresponding to
this minimization problem describes a solution for the velocity field.
Another common way of studying power-law flow
is to simultaneously find a pair $(\bu,p)$ that satisfies (\ref{eq:int5}). This is called the \textit{mixed weak
formulation} and it is written down precisely in (\ref{eq:mwf1}) and (\ref{eq:mwf2}). The connection bewteen these two formulations comes from a technical
inf-sup condition (see Theorem \ref{thm:lbb1}-\ref{thm:lbb2} below) which is frequently called the \textit{LBB condition} named after Ladyzhenskaya, Brezzi, and Bab\v uska. This condition can be stated as the following: $\exists \, \beta > 0$ such that
\begin{equation}\label{eq:int7}
\inf_{ q \in
L^{r'}_0(\Omega)}\sup_{ \bv \in \bW_0^{1,r}(\Omega)} \frac{\langle \nabla \cdot
\bv , q\rangle}{\|q\|_{0,r'}\|\bv \|_{0,r}}\ge \beta.
\end{equation}
When this condition holds, the variational formulation and the mixed weak
formulation are equivalent in the sense that $\bu$ is a solution of
the variational formulation if and only if $(\bu ,p)$ is a solution of
the mixed weak formulation with $p$ being solved in terms of $\bf u$
using the inverse of the gradient operator. The pressure can then be computed after the velocity is known provided that the \textit{LBB} condition holds.
The \textit{LBB} condition is well known to hold for the linear problem ($r=2$) on
Lipschitz domains in any dimension \cite{Ba}, \cite{Br},
\cite{La}. Since any bounded convex
domain has a Lipschitz continuous boundary \cite[Corollary
1.2.-2.3]{G} we conclude that the \textit{LBB} condition holds for our $\Omega$ when $r=2$. For the nonlinear problem $1 < r \le 2$, the \textit{LBB} has been
shown (\cite{AG1}, \cite{BN}) for smooth domains in dimension $d=2$.
This is generalized in \cite{AG2} where it is shown that the \textit{LBB} condition holds for all
dimensions
$d>1$ and for the full range $ 1< r <\infty$ in Lipschitz domains. Thus, the mixed weak formulation makes sense and it is
equivalent to the variational formulation in our setting in this work.
We note that, in the variational formulation, $\bu$ is defined as the
minimum of a convex functional on a separable Banach space, thus (\ref{eq:int6}) always has a unique
solution $\bu$ for any $ 1 < r < \infty$, whether or not the \textit{LBB}
condition holds. However, the pressure $p$ is not necessarily well
defined, so that results involving the pressure function $p$ typically
require the additional assumption of an \textit{LBB} condition and it is not known if this condition holds in a nonconvex domain.
The two formulations discussed above for the power-law Stokes problem
are both useful for the numerical analysis of the problem. A finite
element analysis of the power-law Stokes problem using the mixed weak
formulation has been studied by several authors, for example, in \cite{BN}, \cite{BGN}, \cite{BL1}, \cite{BL2}.
Finite element analysis using the direct variational formulation
(\ref{eq:int6}) requires one to solve a constrained minimization problem
and construct finite element spaces with approximately divergence free interpolation
functions. Thus, the variational formulation and its associated
constrained minimization problem is more difficult for both analysis
and numerical approximation. A natural way to overcome this
difficulty is to introduce a penalty functional that eliminates the
constraint. The first use of the penalty function method in
conjuction with the finite element method is due to Bab\v uska
\cite{Ba}. The method was quickly adopted as a standard tool for the
finite element analysis of viscous, incompressible fluid flows
\cite {Zi}. Extensive studies of the penalty method applied to fluid flow
problems, both experimentally and mathematically, have appeared from
the late seventies to the present day. Here we cite only several
important articles among them \cite{F1}, \cite{F2}, \cite {FKg},
\cite{H}, \cite{HLB}, \cite{OKS}, \cite{Od1}, \cite{Od2}, \cite{R},
\cite{R2}, and \cite{Zi}. A very general mathematical analysis of the penalty method applied to nonlinear
problems including a class of non-Newtonian fluid flow problems was
presented by Oden \cite{Od2}. His work provides some important convergence
results. It appears that when the penalty method is applied to the Newtonian Stokesian flow problems, the resultant matrices are ill-conditioned. However, this
deficiency can be overcomed by the use of reduced integration techniques.
For a given penalty parameter $\epsilon>0$, the penalty formulation
requires the unconstrained minimization of the nonlinear convex
functional
\[ J_{\epsilon}(\bu)=J(\bu)+\frac{1}{r\epsilon}\int_{\Omega}|\nabla\cdot \bu|^r d\bx \]
over the Sobolev space
${\bW}_0^{1,r}(\Omega)$.
The corresponding pressure $p^\epsilon$ is defined in terms of the minimizer $\bu^\epsilon$. We prove that the penalty
approximation $\bu^\epsilon$ of the unconstrained minimization problem
min$\{J_\epsilon(\bu): \bu \in \bW_0^{1,r}(\Omega)\}$ converges to the
true solution $\bu$ of min$\{J(\bu): \bu \in X \}$ as $\epsilon
\to 0$ for any $1 < r < \infty$ without assuming that the domain $\Omega$ is convex. This convergence result is only
for the velocity field since the pressure may be undefined. However,
because of the more general variational setting, this result
establishes the validity of a penalty approximation even when the \textit{LBB} condition fails to hold. This is a
convergence result, not an error estimate, but it doesn't require the
\textit{LBB} condition (\ref{eq:int7}).
Here we are writing $(\bu, p)$ for the unique solution of the mixed
formulation (\ref{eq:mwf2}), and $\bu^\epsilon$ the penalty solution.
When the \textit{LBB} condition holds, we obtain error
estimates for the velocity field $\sn{1}{r}{\bu - \bu^\epsilon} = O(\epsilon^{g_1(r)})$, where
$g_1(r) = \frac{1}{(r-1)(3-r)}$ for $1 < r \le 2$, and $g_1(r) =
\frac{1}{(r-1)^2}$ for $2 \le r < \infty $. Let $\phi(z) = |z|^{r-2}z, z \in \mathbb{R}$ and let $p^{\epsilon}
= c-\phi(\nabla\cdot\bu^{\epsilon})$, where $ c =
\int_{\Omega}\phi(\nabla\cdot\bu^{\epsilon})d\bx$. We also show
error estimates for the pressure $ \|p - p^{\epsilon}\|_{0,r'} =
O(\epsilon^{g_2(r)})$ where $g_2(r) = \frac{1}{(3-r)}$ for $ 1 < r \le
2$ and $g_2(r) = \frac{1}{(r-1)^2}$ for $2 \le r < \infty $. These
rates of convergence reduce to known results $\sn{1}{2}{\bu -
\bu^\epsilon}+\|p - p^{\epsilon}\|_{0,2} = O(\epsilon)$ for the Newtonian
case $ r = 2 $ as discussed in \cite{GR}, \cite{OKS}, \cite{R}, and
\cite{T}.
This is an important feature of the unconstrained
penalty minimization formulation which makes it convenient for error
analysis and numerical implementation. In contrast, the mixed weak
formulation requires the solution of a system of nonlinear equations and a
discrete \textit{LBB} condition.
Thus, we provide a
mathematical analysis of the penalty
method applied to the power-law Stokes problem in the variational
formulation. To our knowledge, this is a generalization of the analysis available in the literature for Newtonian flows.
Numerical experiments have been performed on power-law flow
problems using the penalty method in the engineering literature \cite{LLH}, \cite{LWT}, \cite{R}, \cite{RP}, and \cite{Zi}.
Since pressure must then be calculated from the computed
velocity field, the accuracy of the pressure is lower than
that of the velocity as shown in our error estimates. It is interesting to note that the penalty term $\frac{1}{2\epsilon}\int_{\Omega}|\nabla\cdot \bu|^2 d\bx $ was used in \cite{RP} to approximate power-law flows instead of
$\frac{1}{r\epsilon}\int_{\Omega}|\nabla\cdot \bu|^r d\bx$, which is used in this work and reasonable numerical results were obtained without a mathematical analysis.
\section{Preliminaries}
We begin by establishing some notation. Let $ L^r(\Omega)$ for $ 1 <
r < \infty$ be the space of real scalar functions defined on $\Omega$ whose
$r^{\hbox{th}}$ power is absolutely integrable with respect to
Lebesgue measure $d{\bx}=dx_1 \ldots dx_d$. This is a Banach space with
the norm $\sn{0}{r}{u} = (\int_\Omega |u({\bx})|^r\, d{\bx})^{1/ r}$. The
Sobolev space $W^{k,r}(\Omega)$ is the space of functions in $
L^r(\Omega)$ with distributional derivatives up to order $k$ also in $
L^r(\Omega)$. The norm for this space is $ \sn{k}{r}{u} =
(\int_\Omega \sum_{|j|\le k}|D^j u({\bx})|^r \, d{\bx})^{1/ r}$, where we use
the standard multi-index notation. That is, for $j = (j_1 \ldots
j_d)\in {\mathbb N}^d$, where ${\mathbb N}$ is the set of natural
numbers, define $|j|=\sum\limits_{i=1}^{d}j_i$ and write the partial
derivative $D^{j}
u({\bx})=\frac{\partial^{|j|}u}{\partial^{j_1}{x_{1}}\ldots
\partial^{j_d}{x_d}}$.
The closure of $ C^{\infty}_0(\Omega)$ in $W^{k,r}(\Omega)$ is denoted
by $W_0^{k,r}(\Omega)$. For systems of equations, we need the product
spaces defined by ${\mathbf L}^r(\Omega)= [L^r(\Omega)]^d$, ${\mathbf
W}^{k,r}(\Omega)=[W^{k,r}(\Omega)]^d$, and ${\mathbf
W}_0^{k,r}(\Omega)=[W_0^{k,r}(\Omega)]^d$. The norm for $\bv = (v_1,
\ldots, v_d)\in {\mathbf L} ^{r}(\Omega)$ is $\sn{0}{r}{\bv} = (
\int_{\Omega} \sum_{i=1}^d |v_i|^r\, d\bx )^{1/r}$. For
$\bv\in \bW^{k,r}(\Omega)$ we have norm $\sn{k}{r}{\bv} = (
\int_{\Omega} \sum_{i=1}^d \sum_{|j|\le k}|D^j v_i|^r\,
d\bx)^{1/r}$. It is well known \cite{A}, that the seminorm $|\bv|_{k,r}
=( \int_{\Omega}(\sum_{i=1}^d \sum_{|j| = k}|D^j v_i|^r\,
d\bx)^{1/r}$ is equivalent to $\sn{k}{r}{\bv}$ for $\bv\in
\bW_0^{k,r}(\Omega)$. In addition, by Korn's inequality, see \cite {A},
the norm $\| \bD(\bu)\|_{0,r}=( \int_\Omega \sum_{i,j=1}^d
|D_{ij}(\bu)|^r\, d\bx)^{1/r}$ is equivalent to $\|\bu\|_{1,r}$ in
$\bW_0^{1,r}(\Omega)$. For $10$ is independent of $\bx$ and $\by$.
\begin{equation}
\gathered
|\bx - \by|^2 \leq C (\bphi(\bx) - \bphi(\by))\cdot(\bx - \by)
(|\bx| + |\by|)^{2-r}, \\
|\bphi(\bx) - \bphi(\by)| \leq C
|\bx-\by|^{r -1},\, \text{for~} 1 < r < 2;
\endgathered \label{eq:pre1}
\end{equation}
\begin{equation}\gathered
|\bx -\by|^{r} \leq C \ (\bphi(\bx) - \bphi(\by))\cdot(\bx - \by), \\
|\bphi(\bx) - \bphi(\by)| \leq C |\bx-\by|(|\bx| + |\by|)^{r-2},\,
\text{for~} 2 \leq r < \infty.
\endgathered\label{eq:pre2}
\end{equation}
They were proved for the case $d=2$ in Glowinski and Marroco \cite{GM} and were generalized by Barrett and Liu \cite{BL2}, see also \cite{Ch}. A simple proof for
general $d$ is shown in DiBenedetto \cite{Di}.
Finally, for completeness, we quote some important
results from convex analysis and functional analysis which will apply to our problem. See \cite{C} or \cite{FK} for further details.
Let $X$ be a reflexive Banach space with dual space $X^*$. Suppose the operator $T: X\to X^*$. Let $\langle u^*, u\rangle$ denote the duality pairing between $u\in X$ and $u^*\in X^*$. We say $T$ is \textit{bounded} if $\exists\, C>0$ such that $\Vert T u\Vert_{X^*} \leq C \Vert u\Vert_{X}$ for all $u\in X$. The operator
$T$ is \textit{monotone} if $\langle T u - T v, u-v\rangle \ge 0$ for all $u, v\in X$. $T$ is \textit{strictly monotone} if the inequality is strict for all $u, v\in X$ with $u\ne v$. A \textit{coercive} operator $T$ satisfies $\lim_{\Vert u \Vert \to \infty} \frac{\langle T u, u\rangle}{\Vert u\Vert} = \infty$.
Finally, we say
$T$ is \textit{potential} if $\exists$ a functional $J:X\to \mathbb{R}$ such that $J'(u) = T u$ for all $u\in X$ (i.e. $\langle J'(u), v\rangle = \langle T(u), v\rangle$
for all $u, v\in X$).
\begin{theorem}\label{thm:pre1}
Let $T: X\to X^*$ be a bounded, monotone, coercive, potential operator. Then $TX= X^*$. Thus, $Tu=f$ has a solution for every $f\in X^*$. Moreover, if $T$ is strictly monotone, then $Tu=f$ has a unique solution.
\end{theorem}
\begin{theorem}\label{thm:pre2}
Let $J$ be a functional defined on $X$ such that $\lim\limits_{\Vert u\Vert \to \infty} J(u) = \infty$. If $J$ is either (i) continuous and convex on $X$ or (ii) weakly lower semicontinuous on $X$ then $\inf_{u\in X} J(u) > -\infty$ and there exists at least one $u_0\in X$ such that $J(u_0) = \inf_{u\in X} J(u)$. Moreover, if $J$ is continuous and strictly convex on $X$ then there is precisely one such $u_0$.
\end{theorem}
\begin{theorem}\label{thm:pre3}
Let $J:X\to \mathbb{R}$ be a functional with a local extremum at $u_0\in X$.
If $\langle J'(u_0), v\rangle$ exists for some $v\in X$ then
$\langle J'(u_0), v\rangle =0$.
\end{theorem}
\section{The Variational Formulation (VF) of the Stokes Problem}
In order to apply the penalty method, we first consider the
variational formulation of (\ref{eq:int5}) given in (\ref{eq:int6}).
The Euler-Lagrange equation associated to problem
(\ref{eq:int6}) is
\begin{equation}\label{eq:vf1}
\langle A(\bu),\bv
\rangle = \langle{\mathbf f},\bv \rangle\, \quad \forall~\bv \in \bX
\end{equation}
where
$A:\bW_0^{1,r}(\Omega)\to \bW^{-1,r'}(\Omega)$ is defined by
\[
\langle A(\bu),\bv \rangle = k \int_{\Omega}
\abs{\bD(\bu)}^{r-2}\bD(\bu):\bD(\bv )d\bx \quad \forall \bv\in
\bW_0^{1,r}(\Omega).
\] We use $\langle \cdot, \cdot\rangle$ to denote the
duality pairing between $\bW_0^{1,r}(\Omega)$ and
$\bW^{-1,r'}(\Omega)$ as well as between ${\bf L}^{r}(\Omega)$ and ${\bf L}^{r'
}(\Omega)$. In particular, we only need to assume
${\mathbf f}\in \bW^{-1,r'}(\Omega)$ for this formulation.
By using (\ref{eq:pre1}) and (\ref{eq:pre2}), we obtain the following where
$C>0$
denotes a generic constant independent of $\bu$ and $\bv$.
\begin{equation} \gathered
\|\bu - \bv\|_{1,r}^2 \leq C\langle A(\bu) - A(\bv ),{\bu-\bv}
\rangle (\|\bu\|_{1,r} + \|\bv\|_{1,r})^{2-r}, \\
\| A(\bu) -
A(\bv )\|_{-1,r'} \le C \|\bu - \bv\|_{1,r}^{r-1},\text{ for~}
1 < r \leq 2; \endgathered \label{eq:vf2} \end{equation}
\begin{equation}\gathered
\sn{1}{r}{\bu - \bv}^r \le C\langle A(\bu) - A(\bv ),
{\bu-\bv}\rangle,\\
\sn{-1}{r'}{A(\bu) - A(\bv )} \le C
\sn{1}{r}{\bu - \bv} (\sn{1}{r}{\bu} + \sn{1}{r}{\bv})^{r-2},
\text{ for~} 2 \leq r < \infty.
\endgathered\label{eq:vf3}
\end{equation}
From (\ref{eq:vf2}) and (\ref{eq:vf3}), $A$ can be shown to be a bounded, monotone, coercive, potential operator on $\bX=\{\bu\in\bW^{1,r}_0(\Omega) : \nabla\cdot\bu = 0\}$. We conclude using Theorem \ref{thm:pre1} that (\ref{eq:vf1}) has a unique solution~$\bu$. It follows that $J:\bX \to \mathbb{R}$
is a continuous, strictly convex functional on $\bX$ and
$$\lim\limits_{\|\bu\|_{1,r} \to \infty} J(\bu) = \infty.$$ Thus, by
Theorem \ref{thm:pre2} problem (\ref{eq:int6}) has one and
only one solution $\bu$ and hence (\ref{eq:vf1}) and (\ref{eq:int6}) are equivalent.
Note that in this formulation, the pressure
function $p$ does not appear.
\section{The Mixed Weak Formulation (MWF) of the Stokes Problem}
For the mixed weak formulation of the Stokes problem (\ref{eq:int5}),
we suppose ${\mathbf f}\in \bW^{-1,r'}(\Omega)$. The problem is
then to simultaneously find $\bu \in \bW_0^{1,r}(\Omega)$ and $ p \in
L^{r'}_0(\Omega)$, such that
\begin{equation}\label{eq:mwf1}
\gathered
k \int_\Omega \abs{\bD(\bu)}^{r-2}\bD(\bu):\bD(\bv)\,d\bx -
\int_\Omega p\nabla\cdot \bv \,d\bx
= \int_\Omega \mathbf f\cdot \bv \,d\bx\,\quad
\forall \bv \in \bW_0^{1,r}(\Omega) \\
\int_\Omega q \nabla \cdot \bu \,d\bx
=0\,\quad\quad \forall q \in L^{r'}_0(\Omega). \endgathered
\end{equation}
If we let $b(p,\bv)$ be the bilinear form defined on
$L^{r'}_0(\Omega)\times \bW_0^{1,r}(\Omega)$ by $ b(p,\bv )
=\int_\Omega p\nabla \cdot\bv \,d\bx$, then the weak formulation
(\ref{eq:mwf1}) can be rewritten as the problem of finding $(\bu,p) \in
\bW_0^{1,r}(\Omega) \times L^{r'}_0(\Omega)$ such that
\begin{equation}\label{eq:mwf2}
\begin{split} \langle A(\bu),\bv
\rangle - b(p,\bv ) &= \langle{\mathbf f},\bv \rangle\, \quad
\forall~\bv \in \bW_0^{1,r}(\Omega),\\ b(q,\bu) &= 0\, \quad \forall~q
\in L^{r'}_0(\Omega).\\ \end{split} \end{equation}
The existence and uniqueness of solutions of (\ref{eq:mwf1}) and
(\ref{eq:mwf2}) was studied by J.~Ba\-ran\-ger and Najib \cite{BN} and
J.~W. Barrett and W. B. Liu \cite{BL2}.
This mixed weak formulation requires the pressure $p\in
L^{r'}_0(\Omega)$ and the
\textit{LBB} condition is a sufficient condition to guarantee $p\in L^{r'}_0(\Omega)$.
It is possible that the velocity field is well defined within $
\bW_{0}^{1,r}(\Omega)$, but the mixed weak formulation is not well-posed when, e.g., when
the domain is nonconvex and the \textit{LBB} condition fails to hold.
In this case, good approximations of the pressure from the velocity field are not expected from the penalty method.
\section{The \textit{LBB} Condition and Equivalence of (VF) and (MWF)}
Baranger and Najib prove in \cite{BN} that (\ref{eq:mwf2}) is equivalent to
(\ref{eq:vf1}) (and hence (\ref{eq:int6})) for any $1 0$ such that
\begin{equation}\label{eq:lbb1}
0< \beta \le \inf_{ q \in
L^{r'}_0(\Omega)}\sup_{ \bv \in \bW_0^{1,r}(\Omega)} \frac{\langle B
\bv , q\rangle}{\|q\|_{0,r'}\|\bv \|_{1,r}}.
\end{equation}
\end{theorem}
This allows us to conclude the equivalence of (\ref{eq:vf1}) and
(\ref{eq:mwf2}) in our more general setting.
\section{An \textit{A Priori} Bound}
Using Theorems \ref{thm:lbb1} and \ref{thm:lbb2}, we conclude, (\ref{eq:mwf2}) has an unique
solution $(\bu,p)$ in which $\bu$ is the unique solution of (\ref{eq:int6}).
See, e.g, \cite{BN} and \cite {GR}.
\begin{lemma} \label{thm:apb1} Let $\bu$ be the solution of
(\ref{eq:int6}), then $ \|\bu \|_{1,r} \le C$. Suppose further that
(\ref{eq:lbb1}) holds, and let $(\bu,p) $ be the solution of
(\ref{eq:mwf2}). Then $ \|\bu \|_{1,r} \le C$ and $\|p\|_{0,r'} \le C
$. These constants $C>0$ depend only on $r$, $\Omega$ and $\mathbf f$.
\end{lemma}
\begin{proof} In (\ref{eq:vf1}), let $\bv = \bu$. Then
$\langle A(\bu),\bu\rangle = \langle {\mathbf f}, \bu\rangle$. We have,
by (\ref{eq:vf2}), $\sn{1}{r}{\bu}^r \le C
\langle A(\bu),\bu\rangle = C \langle {\mathbf f},\bu\rangle \le C
\|{\mathbf f}\|_{-1,r'}\sn{1}{r}{\bu}$ which implies
\begin{equation}\label{eq:apb1}
\|\bu\|_{1,r} \le C\|{\mathbf f}\|_{-1,r'}^\frac{1}{(r-1)},
\end{equation}
for $2\le r\le\infty$.
Similarly, by (\ref{eq:vf3}), we have $\sn{1}{r}{\bu}^2 \le C \langle
A(\bu),\bu\rangle \sn{1}{r}{\bu}^{2-r}= C \langle {\mathbf
f},\bu\rangle \sn{1}{r}{\bu}^{2-r}$ which gives (\ref{eq:apb1}) for
$1< r\le2$. By (\ref{eq:lbb1}) and (\ref{eq:mwf2})
\begin{align}\label{eq:apb2}
\|p\|_{0,r'} &\le C \sup_{\bv \in
\bW^{1,r}_{0}(\Omega)}\frac{\langle \nabla\cdot
\bv,p\rangle}{\|\bv\|_{1,r}} \notag\\
&= C\sup_{\bv \in \bW^{1,r}_0(\Omega)}
\frac{\langle A(\bu), \bv \rangle - \langle \mathbf{f},
\bv\rangle }{ \| \bv \|_{1,r}} \\
&\le C( \|A(\bu)\|_{-1,r'} + \|{\mathbf f}\|_{-1,r'}). \notag
\end{align}
Upon applying (\ref{eq:vf2}) and
(\ref{eq:vf3}) to the right hand side of (\ref{eq:apb2}) and using
(\ref{eq:apb1}) we conclude that $\|p\|_{0,r'} \le C $ for $10$ be given and suppose that $\bu$ is the solution of (\ref{eq:int6}) and $ {\bu}^{\epsilon}$ is the solution of (\ref{eq:pf2}). Then $ {\bu}^{\epsilon}$ converges strongly to $\bu$ in $\bW_0^{1,r}(\Omega)$
as $\epsilon \to 0$. Furthermore, $\exists\;C$ independent of $\epsilon$ such that $\| {\bu}^{\epsilon}\|_{1,r}\le C$ and
$\|\nabla \cdot {\bu}^{\epsilon}\|_{0,r}\le C \epsilon^{\frac{1}{r}}$. Therefore, it follows that $\nabla \cdot {\bu}^{\epsilon} \to 0$ in $L^{r}(\Omega)$
as $\epsilon \to 0$.
\end{theorem}
\begin{proof} The proof that $ {\bu}^{\epsilon}$ converges strongly to $\bu$ in
$\bW_0^{1,r}(\Omega)$
as $\epsilon \to 0$ follows along the lines of \cite[Theorem 46.D. and Corollary 46.7]{Z}. We only need to
check that the functional $G(\bv)=\|\nabla \cdot \bv\|_{0,r}$ is weakly sequentially continuous in
$\bW_0^{1,r}(\Omega)$. To this end, let $ \bv_{n} \rightharpoonup \bv$ in $\bW^{1,r}_0(\Omega)$.
Then $\left<\nabla \cdot \bv_{n}, \eta\right> \to \left<\nabla \cdot \bv, \eta\right>$ as $ n \to
\infty$ for any $ \eta \in L^{r'}(\Omega)$. Indeed, $C^\infty(\Omega)$ is dense in $L^{r'}(\Omega)$ and
\[
\langle \nabla \cdot \bv_{n}, \eta\rangle = -
\langle \bv_{n},\nabla \eta\rangle \to -\langle\bv, \nabla \eta\rangle = \langle\nabla \cdot \bv, \eta\rangle, \forall \eta\in C^\infty(\Omega).
\]
Therefore, $\nabla \cdot \bv_{n} \rightharpoonup
\nabla
\cdot \bv$ in $L^{r}(\Omega)$. It is well known that the norm $\|\cdot \|_{0,r}$ is
weakly
sequentially continuous in $L^{r}(\Omega)$. We conclude that
$\lim_{n\to
\infty}
\|\nabla \cdot \bv_{n}\|_{0,r} = \|\nabla \cdot \bv\|_{0,r}$ and hence
$\bu^\epsilon\to \bu$ in $\bW_0^{1,r}(\Omega)$ because of \cite[Theorem 46.D. and Corollary 46.7]{Z}. To complete the proof,
let $\bu=\bv = \bu^\epsilon$ in (\ref{eq:pf1}). We have
\[
\sn{1}{r}{\bu^\epsilon}^r + \frac{1}{\epsilon} \sn{0}{r}{\nabla\cdot\bu^\epsilon}^r \le
C\sn{-1}{r'}{\mathbf f} \sn{1}{r}{\bu^\epsilon}
\]
which implies $\sn{1}{r}{\bu^\epsilon}^{r-1} \le
C\sn{-1}{r'}{\mathbf f} $ and $\sn{0}{r}{\nabla\cdot\bu^\epsilon}^r \le
{\epsilon}C\sn{-1}{r'}{\mathbf f}^{r'}$.
Therefore $\| {\bu}^{\epsilon}\|_{1,r}\le C$ and $\nabla \cdot
{\bu}^{\epsilon} \to 0$ in $L^r(\Omega)$
as $\epsilon \to 0$.
\end{proof}
\noindent Note that when $r=2$, Theorem \ref{thm:pf1} is a well-known result \cite{T}. A generalized version of it for convergence in higher order derivative norms can be found in \cite{AG2}.
\begin{lemma}\label{lemma:pf2}
For each $ \epsilon >0$, let $\bu^{\epsilon}$ denote the unique minimizer of (\ref{eq:pf2}) and let
$p^{\epsilon} = c -\frac{1}{\epsilon}\bphi(\nabla \cdot \bu^\epsilon)$, where $ c=\frac{1}
{\epsilon|\Omega|} \int_{\Omega}\bphi(\nabla \cdot \bu^\epsilon)d\bx$.
If (\ref{eq:lbb1}) holds (which is the
\textit{LBB} condition), then there exists $C>0$ which
depends only on $r$, $\Omega$ and
$\mathbf f$ such that $\|\nabla \cdot {\bu}^{\epsilon}\|_{0,r}\le C \epsilon^{\frac{1}{r-1}}$ and $\|p^{\epsilon}\|_{0,r'}\le C.$
\end{lemma}
\begin{proof}
The pair $(\bu^{\epsilon},p^{\epsilon})$ satisfies
\begin{equation}\label{eq:pf3}
\begin{split}
\langle A({\bu^{\epsilon}}),\bv \rangle - b( p^{\epsilon},\bv ) &=
\langle {\mathbf f},\bv \rangle,\quad \forall \bv \in
\bW_0^{1,r}(\Omega),
\end{split}
\end{equation}
since $\bu^{\epsilon}$ satisfies (\ref{eq:pf1}) and
$\int_\Omega c\nabla \cdot \bv d\bx =c \int_{\partial\Omega}
\bv\cdot {\bf n} ds =0$ by Gauss' Theorem.
Moreover, $p^{\epsilon}\in L^{r'}_0(\Omega)$, and
by (\ref{eq:lbb1}) and (\ref{eq:pf3})
\begin{align*} \sn{0}{r'}{p^\epsilon}
&\le C\sup_{\bv \in \bW^{1,r}_0(\Omega)} \frac{\langle \nabla\cdot\bv,
p^\epsilon\rangle }{ \| \bv \|_{1,r}}=C \sup_{\bv \in \bW^{1,r}_0(\Omega)} \frac{
b(p^\epsilon,\bv) }{ \| \bv \|_{1,r}} \\
&= C\sup_{\bv \in \bW^{1,r}_0(\Omega)}
\frac{\langle A(\bu^\epsilon), \bv\rangle
-\langle{\mathbf f}, \bv\rangle}{ \| \bv \|_{1,r}} \\
&\le C(\sn{-1}{r'}{A(\bu^\epsilon)} + \sn{-1}{r'}{\mathbf f}).
\end{align*}
Using (\ref{eq:vf2}) and (\ref{eq:vf3}) we conclude $\|p^{\epsilon}\|_{0,r'}\le C
(\sn{1}{r}{\bu^\epsilon}^{r-1} + \sn{-1}{r'}{\mathbf f})$, therefore
by Theorem \ref{thm:pf1}
\[ \|p^{\epsilon}\|_{0,r'}\le C \sn{-1}{r'}{\mathbf f}. \]
Let $\bv = \bu - \bu^{\epsilon}$ in (\ref{eq:mwf2}) and
(\ref{eq:pf3}) and then subtract (\ref{eq:pf3}) from
(\ref{eq:mwf2}) to get
\begin{equation}\label{eq:pf4}
\langle A(\bu)-A(\bu^{\epsilon}),\bu-\bu^{\epsilon}\rangle = b(p -
p^{\epsilon},\bu-\bu^{\epsilon}).
\end{equation}
Since $\nabla\cdot\bu =0$ and $\int_\Omega \nabla\cdot\bu^\epsilon\, d\bx =0$, the above equation gives
\begin{equation*}
\langle A(\bu)-A(\bu^{\epsilon}),\bu-\bu^{\epsilon}\rangle + \frac{1}{\epsilon} \langle \bphi(\nabla \cdot \bu^{\epsilon}), \nabla \cdot \bu^{\epsilon}
\rangle =- b(p,\bu^{\epsilon}).
\end{equation*}
Since $ \langle A(\bu)-A(\bu^{\epsilon}),\bu-\bu^{\epsilon}\rangle\ge 0$ and $ |b(p,\bu^{\epsilon})|\le \|p\|_{0,r'} \|\nabla \cdot \bu^{\epsilon}\|_{0,r}$
we have $$\frac{1}{\epsilon}\|\nabla \cdot \bu^{\epsilon}\|_{0, r}^r =\frac{1}{\epsilon} \langle \bphi(\nabla \cdot \bu^{\epsilon}), \nabla \cdot \bu^{\epsilon}\rangle \le \|p\|_{0,r'} \|\nabla \cdot \bu^{\epsilon}\|_{0,r},$$
which gives $\|\nabla \cdot \bu^{\epsilon}\|_{0,r}\le C \epsilon^{\frac{1}{r-1}}$ by using Lemma \ref{thm:apb1} to bound $\|p\|_{0,r'}$.
\end{proof}
\begin{theorem}\label{thm:pf3} Suppose (\ref{eq:lbb1})
(the \textit{LBB} condition) holds.
Let $(\bu,p)$ be the unique solution of
(\ref{eq:mwf2}) and $\bu^{\epsilon}$ be a solution of
(\ref{eq:pf2}). Let
$p^{\epsilon} = c -\frac{1}{\epsilon}\bphi(\nabla \cdot \bu^\epsilon)$, where $ c=\frac{1}
{\epsilon|\Omega|} \int_{\Omega}\bphi(\nabla \cdot \bu^\epsilon)d\bx$.
Then there exists $C>0$ which depends only on $r$,
$\Omega$, and $\mathbf f$ such that
\begin{align*} \|\bu -
\bu^{\epsilon}\|_{1,r}&\le C \epsilon^\frac{1}{ (r-1)(3-r)}\quad
\text{for~} 1 < r \le 2\\ \text{and}\quad \|\bu -
\bu^{\epsilon}\|_{1,r}&\le C \epsilon^\frac{1}{(r-1)^2} \quad
\text{for~} 2 \le r <\infty.
\end{align*}
Furthermore,
\begin{align*} \|p - p^{\epsilon}\|_{0,r'}&\le C
\epsilon^\frac{1}{(3-r)}\quad \text{for~} 1 < r \le 2\\
\text{and}\quad \|p - p^{\epsilon}\|_{0,r'}&\le C
\epsilon^\frac{1}{(r-1)^2}\quad \text{for~} 2 \le r < \infty.
\end{align*}
\end{theorem}
\begin{proof}
By (\ref{eq:lbb1}), we have
\begin{equation}\label{eq:pf6}
\sn{0}{r'}{p - p^\epsilon} \le C \sup_{\bv \in \bW^{1,r}_0(\Omega)}
\frac{b( p- p^{\epsilon}, \bv) }{\sn{1}{r}{\bv}}. \end{equation}
Subtracting (\ref{eq:pf3}) from (\ref{eq:mwf2}) gives
\begin{equation}\label{eq:pf7} b(p-
p^{\epsilon}, \bv) = \langle A(\bu),\bv \rangle -\langle
A(\bu^{\epsilon}),\bv \rangle.
\end{equation}
Using (\ref{eq:pf6})
and (\ref{eq:pf7}) we get
\begin{equation}\label{eq:pf8}
\sn{0}{r'}{p
- p^\epsilon} \le C \|A(\bu) - A(\bu^{\epsilon})\|_{-1,r'}.
\end{equation}
By (\ref{eq:pf8}), (\ref{eq:vf2}), Lemma \ref{thm:apb1} and Theorem \ref{thm:pf1}, we have
\begin{equation}\label{eq:pf9}
\sn{0}{r'}{p -
p^\epsilon} \le C \sn{1}{r}{\bu - \bu^\epsilon} \left(
\sn{1}{r}{\bu} + \sn{1}{r}{\bu^\epsilon} \right)^{r-2}\le
C\sn{1}{r}{\bu - \bu^\epsilon}
\end{equation}
for $ 2 \le r <\infty$.
Similarly, using (\ref{eq:pf8}) and (\ref{eq:vf3}) we get for $ 1 < r \le 2$
\begin{equation}\label{eq:pf10}
\sn{0}{r'}{p - p^\epsilon} \le
C\sn{1}{r}{\bu - \bu^\epsilon}^{r-1}.
\end{equation}
Since $(\bu, p)$ solves (\ref{eq:mwf2})
we have $\nabla\cdot \bu = 0$ and hence by (\ref{eq:pf4})
\begin {equation*}
\langle A(\bu) -A(\bu^\epsilon), \bu - \bu^\epsilon\rangle = -b(p-p^\epsilon,
\bu^\epsilon).
\end{equation*}
Therefore, for $1