2022_12_07_8d973a3879ded5997be7g.tex

\documentclass[10pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage[version=4]{mhchem}
\usepackage{stmaryrd}
\usepackage{bbold}
\usepackage{hyperref}
\hypersetup{colorlinks=true, linkcolor=blue, filecolor=magenta, urlcolor=cyan,}
\urlstyle{same}
\usepackage{graphicx}
\usepackage[export]{adjustbox}
\graphicspath{ {./images/} }

\title{A comparative study of mixed least-squares FEMs for the incompressible Navier-Stokes equations }


\author{A. Schwarz ${ }^{1}$, M. Nickaeen ${ }^{2,3}$, S. Serdas ${ }^{1}$, A. Ouazzi $^{2}$, J. Schröder ${ }^{1}$, S. Turek Sc $^{2}$\\
${ }^{1}$ Institute for Mechanics, Faculty of Engineering,\\
University Duisburg-Essen, Universitätsstr. 15, 45141 Essen, Germany\\
alexander.schwarz, serdar.serdas, j.schroeder@uni-due.de\\
${ }^{2}$ Institute for Applied Mathematics, Department of Mathematics,\\
Dortmund University of Technology, Vogelpothsweg 87, 44227 Dortmund, Germany\\
masoud.nickaeen, abderrahim.ouazzi, stefan.turek@math.tu-dortmund.de\\
${ }^{3}$ Department of Mechanical and Aerospace Engineering, Rutgers University\\
Piscataway, New Jersey 08854, USA}
\date{}


\begin{document}
\maketitle


\begin{abstract}
In the present contribution we compare different mixed least-squares finite element formulations (LSFEMs) with respect to computational costs and accuracy. In detail, we consider an approach for Newtonian fluid flow, which is described by the incompressible Navier-Stokes equations. Starting from the residual forms of the equilibrium equation and the continuity condition, various first-order systems are derived. From these systems least-squares functionals are constructed by means of $L^{2}$-norms, which are the basis for the associated minimization problems. The first formulation under consideration is a div-grad first-order system resulting in a threefield formulation with stresses, velocities, and pressure as unknowns. This S-V-P formulation is approximated in $H(\operatorname{div}) \times H^{1} \times L^{2}$ on triangles and for comparison also in $H^{1} \times H^{1} \times L^{2}$ on quadrilaterals. The second formulation is the well-known div-curl-grad first-order velocityvorticity-pressure (V-V-P) formulation. Here all unknowns are approximated in $H^{1}$ on quadrilaterals. Besides some numerical advantages, as e.g. an inherent symmetric structure of the system of equations and a directly available error estimator, it is known that least-squares methods have also a drawback concerning mass conservation, especially when lower-order elements are used. Therefore, the main focus of the work is on performance and accuracy aspects on the one side for finite elements with different interpolation orders and on the other side on the usage of efficient solvers, for instance of Krylov-space or multigrid type. In order to demonstrate the capability of the formulations the results for some well-known benchmark problems are presented and discussed.
\end{abstract}

Keywords: Least-squares FEM, V-V-P formulation, S-V-P formulation, Navier-Stokes, Multigrid

\section{Introduction}
In the last years mixed least-squares finite elements (LSFEMs) were successfully applied to many problems in fluid dynamics and solid mechanics. A reason for the increasing attention of least-squares variational principles is due to some theoretical and computational advantages compared to the Galerkin method as for instance a directly available a posteriori error estimator without additional costs, symmetric positive (semi-)definite system matrices ${ }^{1}$ and the unrestricted choice of the polynomial degree of the finite element spaces (the inf-sup condition is not required), see e.g. Jiang [1998], Bochev and

${ }^{1}$ In general the system matrices of LSFEMs are symmetric positive definite, but due to numerical issues during the Newton iteration steps they may become positive semi-definite only. Gunzburger [2009] and Kayser-Herold and Matthies [2005]. Nevertheless, the Galerkin method is today one of the most employed variational approaches for all kind of boundary value problems in computational mechanics. Despite the mentioned advantages the least-squares method only plays a minor role as variational principle for solving partial differential equations. One reason is the fact that lower-order LSFEMs provide a suboptimal performance. Indeed, there are some contributions with lower-order elements with a satisfying performance in fluid dynamics or solid mechanics, see e.g. Chang and Nelson [1997] and Schwarz et al. [2010], but these approaches are modified least-squares formulations, which may introduce other drawbacks. An illustrative example for the poor solution of lower-order interpolation was given in Pontaza [2003] for the driven cavity problem. A closely related problem is the lack of mass conservation, see e.g. Chang and Nelson [1997] and Deang and Gunzburger [1998]. Different strategies, as e.g. mesh refinement, residual weighting, the introduction of Lagrange multipliers, and improved velocity-pressure coupling can be used in order to overcome the disadvantage, see Bochev et al. [1999], Deang and Gunzburger [1998], Bolton and Thatcher [2005], Chang and Nelson [1997], Pontaza [2006], Heys et al. [2006; 2007], and Heys et al. [2009]. The application of higher-order spectral finite elements demonstrates also an improvement of the mass conservation, see e.g. Pontaza and Reddy [2003; 2004], and Proot and Gerritsma [2006], as well as the recently proposed non-conforming locally conservative LSFEM, see Bochev et al. [2013], with a piecewise divergence-free basis for the velocity.

As it was mentioned, the resulting LSFEM system is symmetric and positive (semi-) definite, see Bochev and Gunzburger [2009]. This permits the design of efficient solvers, which exploit the properties of the least-squares system with respect to both the conjugate gradient (CG) and the multigrid methods. Algebraic multigrid preconditioned CG methods have been widely used to solve the Stokes equation, see Heys et al. [2005], and the NS equations, see Heys et al. [2006; 2007; 2009]. The advantage of this scheme is that the Krylov method, here CG, reduces the error in eigenmodes that are not being effectively reduced by multigrid. A geometric multigrid preconditioned CG solver was used by Ranjan and Reddy [2012] for the Spectral/hp LSFEM solution of the NS equations. They demonstrated superior convergence of the multigrid solver compared to the Jacobi preconditioning. Most recently, Nickaeen et al. [2014] have developed a geometric multigrid solver as a preconditioner for the CG (MPCG) iterations to solve the vorticity formulation of the NS equations with LSFEM. A robust and grid independent behavior is demonstrated for the solution of different flow problems with both bilinear and biquadratic finite elements.

In the present work we study different mixed least-squares formulations, which are based on div-grad and div-curl-grad first-order systems in the Sobolev spaces $H$ (div), $H^{1}$ and $L^{2}$. We show that these formulations for the NS equations may suffer from the well-known mass loss phenomena, especially for low-order elements, but a clear improvement in the mass conservation is obtained for higher order finite elements. Furthermore, the influence of different interpolation orders and the performance of direct and iterative solvers are discussed.

The paper is outlined as follows: In section 2 the least-squares method is derived briefly and the related mixed finite element formulations are proposed. We present and discuss different numerical benchmark examples in section 3 and section 4 closes with a conclusion and an outlook.

\section{Least-squares method and mixed finite element formulations}
In the following we provide the least-squares mixed finite element method for first-order systems, see e.g. Jiang [1998] or Bochev and Gunzburger [2009]. We denote the scalar multiplication of two matrices $\boldsymbol{B}, \boldsymbol{C} \in \mathbb{R}^{d \times d}$ and two vectors $\boldsymbol{b}, \boldsymbol{c} \in \mathbb{R}^{d}$ by the sum of component-wise multiplication, i.e. $(\boldsymbol{B}, \boldsymbol{C})=\boldsymbol{B} \cdot \boldsymbol{C}=\operatorname{tr}\left(\boldsymbol{B} \boldsymbol{C}^{T}\right)$ and $(\boldsymbol{b}, \boldsymbol{c})=\boldsymbol{b} \cdot \boldsymbol{c}=\boldsymbol{b}^{T} \boldsymbol{c}$, respectively. Based on this, we define the $L^{2}$ scalar product and the $L^{2}$-norm of a matrix or a vector by the volume integral of the corresponding scalar multiplication, e.g.

$$
(\boldsymbol{b}, \boldsymbol{c})_{0}=\int_{\mathcal{B}} \boldsymbol{b} \cdot \boldsymbol{c} d v \quad \text { and } \quad\|\boldsymbol{b}\|_{0}=\left\{\int_{\mathcal{B}} \boldsymbol{b} \cdot \boldsymbol{b} d v\right\}^{\frac{1}{2}}
$$

and

$$
(\boldsymbol{B}, \boldsymbol{C})_{0}=\int_{\mathcal{B}} \boldsymbol{B} \cdot \boldsymbol{C} d v \quad \text { and } \quad\|\boldsymbol{B}\|_{0}=\left\{\int_{\mathcal{B}} \boldsymbol{B} \cdot \boldsymbol{B} d v\right\}^{\frac{1}{2}}
$$

The least-squares method is to find the minimizer of a functional $\mathcal{F}(\boldsymbol{b})$ by solving the resulting optimization problem

$$
\boldsymbol{b}_{L S}=\underset{\boldsymbol{b} \in \boldsymbol{X}}{\operatorname{argmin}} \mathcal{F}(\boldsymbol{b}),
$$

where $\boldsymbol{X}$ denotes the minimization space. Such a least-squares functional, which is based on at least one residual $\mathcal{R}(\boldsymbol{b})=\mathbf{0}$, can be constructed by means of the quadratic $L^{2}$-norm as

$$
\mathcal{F}(\boldsymbol{b})=\frac{1}{2}\|\mathcal{R}(\boldsymbol{b})\|_{0}^{2} .
$$

The minimization problem is solved using the calculus of variations to obtain an equivalent variational statement, i.e. with the condition that the first variation of the functional, defined as

$$
\delta_{\boldsymbol{b}} \mathcal{F}(\boldsymbol{b} ; \delta \boldsymbol{b})=\lim _{t \rightarrow 0} \frac{\mathcal{F}(\boldsymbol{b}+t \delta \boldsymbol{b})-\mathcal{F}(\boldsymbol{b})}{t}=(\mathcal{R}(\boldsymbol{b}), \delta \mathcal{R}(\boldsymbol{b}, \delta \boldsymbol{b}))_{0}
$$

equals to zero. All quantities $\boldsymbol{b}$ are defined in $d$ dimensions on a domain $\mathcal{B}$, which is parameterized in $\boldsymbol{x} \in \mathbb{R}^{d}$. In order to describe a boundary value problem for stationary Newtonian fluid flow, we consider the incompressible Navier-Stokes equations consisting of balances of momentum and mass

$$
\begin{array}{r}
-\rho \nabla \boldsymbol{v} \boldsymbol{v}+2 \rho \nu \operatorname{div} \nabla^{s} \boldsymbol{v}-\nabla p=\mathbf{0} \\
\operatorname{div} \boldsymbol{v}=0 .
\end{array}
$$

Here, $\boldsymbol{v}$ denotes the velocities, $p$ the pressure, $\rho$ the density and $\nu$ the kinematic viscosity of a medium flowing through the domain. Furthermore,

$$
\nabla^{s} \boldsymbol{v}=\frac{1}{2}\left(\nabla \boldsymbol{v}+[\nabla \boldsymbol{v}]^{T}\right)
$$

denotes the symmetric velocity gradient.

\subsection{Stress-Velocity-Pressure formulation}
In order to derive the first formulation under consideration we introduce the Cauchy stresses $\boldsymbol{\sigma}=2 \rho \nu \nabla^{s} \boldsymbol{v}-p \mathbf{1}$ as additional variable in system (6). We obtain

$$
\begin{aligned}
\operatorname{div} \boldsymbol{\sigma}-\rho \nabla \boldsymbol{v} \boldsymbol{v} & =\mathbf{0} \\
\boldsymbol{\sigma}-2 \rho \nu \nabla^{s} \boldsymbol{v}+p \mathbf{1} & =\mathbf{0} \\
\operatorname{div} \boldsymbol{v} & =0,
\end{aligned}
$$

the stress-velocity-pressure formulation (S-V-P), see e.g. Cai et al. [2004]. Besides that, similar approaches for stress-based least-squares formulations have been considered and investigated in Bochev and Gunzburger [1995], Chang et al. [1995] and Ding and Tsang [2003]. The resulting least-squares functional $\mathcal{F}_{1}$ is given in terms of the $L^{2}$-norm as

$$
\mathcal{F}_{1}(\boldsymbol{\sigma}, \boldsymbol{v}, p)=\frac{1}{2}\left(\left\|\frac{1}{\sqrt{\rho}}(\operatorname{div} \boldsymbol{\sigma}-\rho \nabla \boldsymbol{v} \boldsymbol{v})\right\|_{0}^{2}+\left\|\frac{1}{\sqrt{\rho \nu}}\left(\boldsymbol{\sigma}-2 \rho \nu \nabla^{s} \boldsymbol{v}+p \mathbf{1}\right)\right\|_{0}^{2}+\|\operatorname{div} \boldsymbol{v}\|_{0}^{2}\right),
$$

where we have used some additional physically motivated weights on the first and second residual, see e.g. Schwarz and Schröder [2011] and Serdas et al. [2013]. The required variation of the functional is explicitly given by

$$
\begin{aligned}
\delta_{\boldsymbol{v}} \mathcal{F}_{1}\left(U_{1} ; \delta \boldsymbol{v}\right)= & \int_{\mathcal{B}} \operatorname{div} \delta \boldsymbol{v} \cdot \operatorname{div} \boldsymbol{v} d V-\frac{1}{\rho \nu} \int_{\mathcal{B}} 2 \rho \nu \nabla^{s} \delta \boldsymbol{v} \cdot\left(\boldsymbol{\sigma}-2 \rho \nu \nabla^{s} \boldsymbol{v}+p \mathbf{1}\right) d V \\
& -\frac{1}{\rho} \int_{\mathcal{B}}(\rho \nabla \delta \boldsymbol{v} \boldsymbol{v}+\rho \nabla \boldsymbol{v} \delta \boldsymbol{v}) \cdot(\operatorname{div} \boldsymbol{\sigma}-\rho \nabla \boldsymbol{v} \boldsymbol{v}) d V=0 \\
\delta_{\boldsymbol{\sigma}} \mathcal{F}_{1}\left(U_{1} ; \delta \boldsymbol{\sigma}\right)= & \frac{1}{\rho} \int_{\mathcal{B}} \operatorname{div} \delta \boldsymbol{\sigma} \cdot(\operatorname{div} \boldsymbol{\sigma}-\rho \nabla \boldsymbol{v} \boldsymbol{v}) d V \\
& +\frac{1}{\rho \nu} \int_{\mathcal{B}} \delta \boldsymbol{\sigma} \cdot\left(\boldsymbol{\sigma}-2 \rho \nu \nabla^{s} \boldsymbol{v}+p \mathbf{1}\right) d V=0 \\
\delta_{p} \mathcal{F}_{1}\left(U_{1} ; \delta p\right)= & \frac{1}{\rho \nu} \int_{\mathcal{B}} \delta p \mathbf{1} \cdot\left(\boldsymbol{\sigma}-2 \rho \nu \nabla^{s} \boldsymbol{v}+p \mathbf{1}\right) d V=0
\end{aligned}
$$

where $U_{1}=(\boldsymbol{\sigma}, \boldsymbol{v}, p)$. In order to seek the minimizer $\mathcal{I}(\boldsymbol{\sigma}, \boldsymbol{v}, p)$ with e.g. the Newton method, we have to linearize (10) with respect to stresses, velocities and pressure. Alternatively, the Newton tangent can also be computed by a standard difference quotient procedure. The interpolation of the unknown variables in $H(\operatorname{div}) \times H^{1} \times L^{2}$ is realized by different conforming finite element spaces, which fulfill directly the boundary conditions on the Dirichlet boundary for $\boldsymbol{v}$ and the Neumann boundary for $\boldsymbol{\sigma}$ and $p$. Finally, the first element formulation results in an triangular element $R T_{m} P_{k} P_{l}$, where $m, k$ and $l$ denote the polynomial order of the interpolation. The stresses are approximated using the so-called Raviart-Thomas space $R T_{m}$, while for the velocities and pressure we use standard Lagrange finite element spaces $P_{k}$ and $P_{l}$. The complete system of equations is obtained by standard finite element assembly operations (similar to Schwarz et al. [2009]) and solved with a standard solver. We further investigate the S-V-P method by using equal-order interpolation functions in $H^{1} \times H^{1} \times L^{2}$. The finite elements in this setting are conforming bilinear and biquadratic nodal elements, which are defined on quadrilateral elements. We solve the resulting system with a multigrid-preconditioned conjugate gradient (MPCG) solver, which is initially developed by Nickaeen et al. [2014] for the solution of vorticity-based NS equations.

\section{$2.2$ Vorticity-Velocity-Pressure formulation}
In order to derive the second investigated formulation we introduce the vorticity $\omega=\nabla \times \boldsymbol{v}$ as additional variable in system (6). We obtain

$$
\begin{aligned}
\rho \nabla \boldsymbol{v} \boldsymbol{v}+\rho \nu \nabla \times \omega+\nabla p & =\mathbf{0} \\
\omega-\nabla \times \boldsymbol{v} & =\mathbf{0} \\
\operatorname{div} \boldsymbol{v} & =0,
\end{aligned}
$$

the vorticity-velocity-pressure formulation (V-V-P), which is the basis for the related least-squares functional $\mathcal{F}_{2}$, see e.g. Jiang [1998]. The resulting least-squares functional $\mathcal{F}_{2}$ is given in terms of the $L^{2}$-norm as

$$
\mathcal{F}_{2}(\omega, \boldsymbol{v}, p)=\frac{1}{2}\left(\left\|\frac{1}{\sqrt{\rho \nu}}(\rho \nabla \boldsymbol{v} \boldsymbol{v}+\rho \nu \nabla \times \omega+\nabla p)\right\|_{0}^{2}+\|(\omega-\nabla \times \boldsymbol{v})\|_{0}^{2}+\|\sqrt{\alpha}(\operatorname{div} \boldsymbol{v})\|_{0}^{2}\right),
$$

where $\alpha$ is a scaling parameter aimed to improve the mass conservation of the LSFEM formulation, see e.g. Deang and Gunzburger [1998] and Bolton and Thatcher [2005; 2006]. The $\mathcal{F}_{2}$ functional is obtained by scaling the momentum balance equations with the inverse viscosity. The required variation of the functional is explicitly given by

$$
\begin{aligned}
\delta_{\boldsymbol{v}} \mathcal{F}_{2}\left(U_{2} ; \delta \boldsymbol{v}\right)= & \alpha \int_{\mathcal{B}} \operatorname{div} \delta \boldsymbol{v} \cdot \operatorname{div} \boldsymbol{v} d V-\int_{\mathcal{B}} \nabla \times \delta \boldsymbol{v} \cdot(\omega-\nabla \times \boldsymbol{v}) d V \\
& +\frac{1}{\rho \nu} \int_{\mathcal{B}}(\rho \nabla \delta \boldsymbol{v} \boldsymbol{v}+\rho \nabla \boldsymbol{v} \delta \boldsymbol{v}) \cdot(\rho \nabla \boldsymbol{v} \boldsymbol{v}+\rho \nu \nabla \times \omega+\nabla p) d V=0 \\
\delta_{\omega} \mathcal{F}_{2}\left(U_{2} ; \delta \omega\right)= & \int_{\mathcal{B}} \delta \omega \cdot(\omega-\nabla \times \boldsymbol{v}) d V \\
& +\int_{\mathcal{B}} \nabla \times \delta \omega \cdot(\rho \nabla \boldsymbol{v} \boldsymbol{v}+\rho \nu \nabla \times \omega+\nabla p) d V=0 \\
\delta_{p} \mathcal{F}_{2}\left(U_{2} ; \delta p\right)= & \frac{1}{\rho \nu} \int_{\mathcal{B}} \nabla \delta p \cdot(\rho \nabla \boldsymbol{v} \boldsymbol{v}+\rho \nu \nabla \times \omega+\nabla p) d V=0,
\end{aligned}
$$

where $U_{2}=(\omega, \boldsymbol{v}, p)$. Similar to the S-V-P formulation, we use the Newton method for the linearization of the convective terms. All variables are in $H^{1}$ and the equal-order interpolation functions are conforming finite element spaces satisfying Dirichlet velocity boundary conditions. We use continuous nodal linear or quadratic finite element spaces for all unknowns and solve the resulting system with the MPCG solver, see Nickaeen et al. [2014].

\section{Numerical examples}
\subsection{Flow around cylinder}
We simulate a laminar flow around a circular cylinder, see Turek and Schäfer [1996] and \href{http://www.featflow.de/en/benchmarks.html}{www.featflow.de/en/benchmarks.html} for further details concerning this benchmark. For the outflow boundary, we impose the zero normal-stress boundary condition defined by

$$
\boldsymbol{\sigma} \cdot \boldsymbol{n}=(\rho \nu \nabla \boldsymbol{v}-p \mathbf{1}) \cdot \boldsymbol{n}=\mathbf{0} \text { on } \Gamma_{\text {out }} .
$$

We incorporate the zero normal-stress boundary condition (14) into the V-V-P functional (12) with an extra $L^{2}$-norm term acting on the outflow, see Nickaeen et al. [2014].

We study the drag/lift coefficient and the pressure drop across the cylinder. For the definition of these flow parameters one should refer to Turek and Schäfer [1996]. Moreover, we study the mass conservation of the LSFEM formulations. We measure the Global Mass Conservation (GMC) in terms of the fractional change of mass flow rate, defined as

$$
\mathrm{GMC}=\frac{\int_{\Gamma_{i}} \rho(\boldsymbol{n} \cdot \boldsymbol{v}) d \Gamma_{i}-\int_{\Gamma_{o}} \rho(\boldsymbol{n} \cdot \boldsymbol{v}) d \Gamma_{o}}{\int_{\Gamma_{i}} \rho(\boldsymbol{n} \cdot \boldsymbol{v}) d \Gamma_{i}} \times 100
$$

where $\Gamma_{i}$ is the inflow boundary of the domain and $\Gamma_{o}$ is any vertical section between the inflow and the outflow boundaries, including the outflow.

Figure 1 shows the computational mesh of the coarsest level. Correspondingly, Table 1 summarizes the mesh information of different levels. We divide every quadrilateral element into two triangles, therefore the number of elements for the case of the quadratic $R T_{0} P_{2} P_{1}$ and cubic $R T_{1} P_{3} P_{1}$ is twice the number of elements presented in Table 1 .

\begin{center}
\includegraphics[max width=\textwidth]{2022_12_07_8d973a3879ded5997be7g-06}
\end{center}

Figure 1: Flow around cylinder, computational grid on level 1.

Table 1: Mesh information for the flow around cylinder problem, the number of elements (NEL) and the number of degrees of freedom.

\begin{center}
\begin{tabular}{rrrrrrrrrr}
\hline
Level & NEL & \multicolumn{4}{c}{D Degrees of freedom} &      \\
 & \multicolumn{2}{c}{$\mathrm{V}-\mathrm{V}-\mathrm{P}$} & \multicolumn{4}{c}{$\mathrm{S}$ S-V-P} \\
\cline { 3 - 4 }\cline { 6 - 10 }
 & $Q_{1}$ & $Q_{2}$ & $Q_{1}$ & $Q_{2}$ & $R T_{0} P_{2} P_{1}$ & $R T_{1} P_{3} P_{1}$ \\
\hline
1 & 346 & 1,564 & 5,896 & 2,346 & 8,844 & 5,505 & 13,989 \\
2 & 1,384 & 5,896 & 22,864 & 8,844 & 34,296 & 21,390 & 54,966 \\
3 & 5,536 & 22,864 & 90,016 & 34,296 & 135,024 & 84,300 & 217,884 \\
4 & 22,144 & 90,016 & 357,184 & 135,024 & 535,776 & 334,680 & 867,576 \\
5 & 88,576 & 357,184 & $1,422,976$ & 535,776 & $2,134,464$ & $1,333,680$     \\
6 & 354,304 & $1,422,976$ &  & $2,134,464$ &      \\
\hline
\end{tabular}
\end{center}

We present the lift and drag coefficients and the pressure drop across the cylinder at Reynolds number $\mathrm{Re}=20$ for the V-V-P and S-V-P formulations in Table 2 and Table 3 , respectively. Using higher order finite elements, both methods show excellent convergence toward the reference solution. Moreover, the two formulations produce similar results with $Q_{1}$ and $Q_{2}$ elements.

Table 2: $V$ - $V$-P formulation: flow parameters in the flow around cylinder at $\operatorname{Re}=20$.

\begin{center}
\begin{tabular}{ccccc}
\hline
Level & Drag coefficient $C_{D}$ & Lift coefficient $C_{L}$ & Pressure drop $\triangle p$ & GMC-value at $x=2.2$ \\
\hline
$Q_{1}, \alpha=1$ &  &      \\
3 & $3.8910464$ & $0.0023588$ & $0.0771009$ & $32.458835$ \\
4 & $4.8914483$ & $0.0043424$ & $0.1009819$ & $12.828753$ \\
5 & $5.3687854$ & $0.0086759$ & $0.1124486$ & $3.871370$ \\
6 & $5.5234496$ & $0.0101034$ & $0.1161682$ & $1.025463$ \\
$Q_{2}, \alpha=1$ &  &      \\
2 & $5.4639017$ & $0.0092029$ & $0.1148151$ & $1.212050$ \\
3 & $5.5668223$ & $0.0104928$ & $0.1172298$ & $0.114906$ \\
4 & $5.5779343$ & $0.0106055$ & $0.1174841$ & $0.010779$ \\
5 & $5.5792792$ & $0.0106169$ & $0.1175144$ & $0.001096$ \\
$Q_{1}, \alpha=100$ &  &      \\
3 & $4.7949823$ & $0.0474048$ & $0.0893839$ & $3.577138$ \\
4 & $5.2716595$ & $0.0245279$ & $0.1044718$ & $1.053974$ \\
5 & $5.4769244$ & $0.0132686$ & $0.1128306$ & $0.283655$ \\
6 & $5.5500584$ & $0.0109745$ & $0.0116144$ & $0.073273$ \\
$Q_{2}, \alpha=100$ &  &      \\
2 & $4.9640652$ & $0.0047990$ & $0.1050202$ & $0.267959$ \\
3 & $5.4744825$ & $0.0094397$ & $0.1157215$ & $0.035024$ \\
4 & $5.5620013$ & $0.0104345$ & $0.1172938$ & $0.004049$ \\
5 & $5.5762393$ & $0.0105925$ & $0.1174909$ & $0.000473$ \\
ref. $: C_{D}=5.57953523384, C_{L}=0.010618948146, \triangle p=0.11752016697$ &  &      \\
\hline
\end{tabular}
\end{center}

We present the GMC values on the outflow (at $x=2.2$ ) in Table 2 and Table 3. The S-V-P method shows much better mass conservation as compared to the V-V-P method. However, by using $\alpha=100$ mass conservation is also improved in the V-V-P method. Moreover, this weighting does not have a pronounced effect on the accuracy of the drag and lift coefficients and the pressure drop especially for the $Q_{2}$ element results. Very good results are obtained by the $R T_{1} P_{3} P_{1} \mathrm{~S}-\mathrm{V}-\mathrm{P}$ formulation. Here, the discretization with 54,966 degrees of freedom produces nearly the same solutions for the drag and lift coefficients and the pressure drop as the $R T_{0} P_{2} P_{1}$ formulation with 1,333,680 degrees of freedom or the $Q_{2}(\mathrm{~S}-\mathrm{V}-\mathrm{P})$ with 135,024 degrees of freedom. Figure 2 shows the horizontal velocities for $R T_{0} P_{2} P_{1}$ level 2 and $R T_{1} P_{3} P_{1}$ level 1 . The mass loss for the $R T_{0} P_{2} P_{1}$ is indicated by the not fully developed outflow, especially if compared to the $R T_{1} P_{3} P_{1}$ result.

Next, we analyze the performance of the MPCG solver for the solution of the V-V-P and S-V-P problems using the conforming nodal finite elements $Q_{1}$ and $Q_{2}$. Table 4 shows the number of nonlinear iterations and the corresponding averaged linear solver (MPCG solver) iterations for different levels. We observe a grid-independent convergence behavior and a constant number of iterations with grid refinement at each $\alpha$ for the V-V-P method. Table 3: $S$ - $V$-P formulation: flow parameters in the flow around cylinder at $\operatorname{Re}=20$.

\begin{center}
\includegraphics[max width=\textwidth]{2022_12_07_8d973a3879ded5997be7g-08}
\end{center}

The optimal number of iterations is obtained at $\alpha=1$. The S-V-P method shows a gridindependent solution behavior as well. However, in this case the solver requires more iterations to reach the same convergence criteria. An important remark about the results of Table 4 is that the MPCG solver remains efficient for both low and high order finite elements.

Table 4: The number of nonlinear iterations and the corresponding averaged number of linear solver iterations for flow around cylinder at $\mathrm{Re}=20$, nonlinear and linear solver relative changes are kept below $1 \mathrm{E}-6$ and $1 \mathrm{E}-3$, respectively.

\begin{center}
\begin{tabular}{ccccccccc}
\hline
  & \multicolumn{5}{c}{V-V-P} & \multicolumn{2}{c}{S-V-P} \\
Level & $Q_{1}$ &  & $Q_{2}$ & $Q_{1}$ & $Q_{2}$     \\
\hline
$\alpha$ & $1.0$ & 10 & 100 & $1.0$ & 10 & 100 & $-$ & $-$     \\
\cline { 2 - 8 }
3 & $8 / 4$ & $8 / 4$ & $8 / 6$ & $6 / 5$ & $6 / 6$ & $6 / 10$ & $7 / 19$ & $6 / 12$     \\
4 & $8 / 4$ & $8 / 4$ & $8 / 6$ & $6 / 5$ & $6 / 6$ & $6 / 9$ & $7 / 17$ & $6 / 12$     \\
5 & $7 / 4$ & $7 / 4$ & $8 / 6$ & $6 / 5$ & $6 / 6$ & $6 / 9$ & $7 / 17$ & $6 / 12$     \\
\hline
\end{tabular}
\end{center}

\begin{center}
\includegraphics[max width=\textwidth]{2022_12_07_8d973a3879ded5997be7g-09}
\end{center}

Figure 2: Horizontal velocities for $R T_{0} P_{2} P_{1}$ level 2 (above) and $R T_{1} P_{3} P_{1}$ level 1 (below).

\section{$3.2$ Lid-driven cavity flow}
We numerically solve the regularized lid-driven cavity flow problem in this section, see Bruneau and Saad [2006] for further details concerning the boundary conditions. We fix the pressure in one point, $p=0$, in the middle of the lower cavity wall.

We investigate two global quantities as defined in Bruneau and Saad [2006]. The first one is the kinetic energy defined as

$$
E=\frac{1}{2}\left\|\boldsymbol{u}_{h}\right\|_{0, \Omega}^{2}
$$

and the other quantity is the enstrophy defined as follows

$$
Z=\frac{1}{2}\left\|w_{h}\right\|_{0, \Omega}^{2} .
$$

For the S-V-P formulation, the enstrophy is obtained by post-processing using the vorticity definition.

Table 5 shows the mesh information of the two LSFEM formulations with different finite element combinations. As in the last example the number of elements for the case of $R T_{m} P_{k} P_{l}$ is twice the number of elements presented in Table 5.

We compare the kinetic energy values of the V-V-P and S-V-P LSFEMs with the mixed finite element method (MFEM) results obtained from FEATFLOW for different Re. We present the V-V-P formulation results for $\alpha=1$ and $\alpha=100$ in Table 6 . The results accurately converge towards the reference solution with mesh refinement at the different $\mathrm{Re}$ Table 5: Mesh information for the regularized driven cavity, the number of elements (NEL) and the number of degrees of freedom.

% \begin{center}
% \begin{tabular}{rrrrrrrrrr}
% \hline
% Level & NEL & \multicolumn{5}{c}{} & \multicolumn{5}{c}{Degrees of Freedom} \\
% \multicolumn{2}{c}{V-V-P} & \multicolumn{5}{c}{S-V-P} &      \\
% \cline { 3 - 4 }\cline { 5 - 10 }
%  & $Q_{1}$ & $Q_{2}$ & $Q_{1}$ & $Q_{2}$ & $R T_{0} P_{2} P_{1}$ & $R T_{1} P_{3} P_{1}$ & $R T_{2} P_{4} P_{1}$ &      \\
% \hline
% 5 & 256 & 1,156 & 4,356 & 1,734 & 6,534 & 4,067 & 10,339 & 19,683 &      \\
% 6 & 1,024 & 4,356 & 16,900 & 6,534 & 25,350 & 15,811 & 40,643 & 77,763 &      \\
% 7 & 4,096 & 16,900 & 66,564 & 25,350 & 99,846 & 62,339 & 161,155 & 309,123 &      \\
% 8 & 16,384 & 66,564 & 264,196 & 99,846 & 396,294 & 247,555 & 641,795 &  &      \\
% 9 & 65,536 & 264,196 & $1,052,676$ & 396,294 & $1,579,014$ & 986,627 &  &      \\
% \hline
% \end{tabular}
% \end{center}

numbers assuring that higher order elements are used. As in the last example good results are obtained by the higher-order S-V-P formulation. In Table 7 , the $R T_{2} P_{4} P_{1}$ discretization with 77,763 degrees of freedom produces a comparable solution as the $R T_{1} P_{3} P_{1}$ formulation with 161,155 degrees of freedom for the kinetic energy at $\operatorname{Re}=1000$. Moreover, the $R T_{0} P_{2} P_{1}$ could not reach the reference solution for $\mathrm{Re}=1000$. Even with 986,627 degrees of freedom the results are not satisfying, which clearly shows the demand for higher-order interpolations in LSFEMs. In addition, we present the nodal element S-V-P formulation results in Table 7 . The accuracy of the results are similar to those obtained with the $\mathrm{V}-\mathrm{V}-\mathrm{P}$ formulation at $\alpha=1$. The poor performance of the $Q_{1}$ finite elements is improved by using higher-order $Q_{2}$ elements.

Table 6: $V$ - $V$ - $P$ formulation: Convergence of the kinetic energy for the regularized cavity problem and comparison with MFEM results by Damanik et al. [2009].

\begin{center}
\includegraphics[max width=\textwidth]{2022_12_07_8d973a3879ded5997be7g-10}
\end{center}

Table 8 shows the mesh convergence for the kinetic energy and the enstrophy at $\operatorname{Re}=1000$. We compare the $Q_{2}$ element results of the V-V-P with two different $\alpha, \alpha=1$ and $\alpha=100$, and the $R T_{2} P_{4} P_{1}$ element results of the S-V-P with the MFEM results by Damanik et al. Table 7: $S$ - $V$-P formulation: Convergence of the kinetic energy for the regularized cavity problem.

\begin{center}
\begin{tabular}{cccccc}
\hline
 & \multicolumn{2}{c}{$\mathrm{LSFEM}$} &      \\
Level & $Q_{1}$ & $Q_{2}$ & $R T_{0} P_{2} P_{1}$ & $R T_{1} P_{3} P_{1}$ & $R T_{2} P_{4} P_{1}$ \\
\hline
$\mathrm{Re}=1$ &  &      \\
6 & $1.624891 \mathrm{E}-02$ & $1.861560 \mathrm{E}-02$ & $1.755107 \mathrm{E}-02$ & $1.862306 \mathrm{E}-02$ & $1.862431 \mathrm{E}-02$ \\
7 & $1.658091 \mathrm{E}-02$ & $1.862344 \mathrm{E}-02$ & $1.827577 \mathrm{E}-02$ & $1.862429 \mathrm{E}-02$ & $1.862438 \mathrm{E}-02$ \\
8 & $1.769969 \mathrm{E}-02$ & $1.862426 \mathrm{E}-02$ & $1.852873 \mathrm{E}-02$ & $1.862438 \mathrm{E}-02$     \\
9 & $1.831796 \mathrm{E}-02$ & $1.862437 \mathrm{E}-02$ & $1.859964 \mathrm{E}-02$ &      \\
$\operatorname{Re}=400$ &  &      \\
6 & $3.242921 \mathrm{E}-02$ & $2.201041 \mathrm{E}-02$ & $1.024579 \mathrm{E}-02$ & $2.164476 \mathrm{E}-02$ & $2.131219 \mathrm{E}-02$ \\
7 & $2.849125 \mathrm{E}-02$ & $2.136759 \mathrm{E}-02$ & $1.400149 \mathrm{E}-02$ & $2.134282 \mathrm{E}-02$ & $2.131461 \mathrm{E}-02$ \\
8 & $2.385371 \mathrm{E}-02$ & $2.131884 \mathrm{E}-02$ & $1.776242 \mathrm{E}-02$ & $2.131732 \mathrm{E}-02$     \\
9 & $2.200378 \mathrm{E}-02$ & $2.131558 \mathrm{E}-02$ & $2.007041 \mathrm{E}-02$ &      \\
$\operatorname{Re}=1000$ &  &      \\
6 & $3.754829 \mathrm{E}-02$ & $3.222452 \mathrm{E}-02$ & $-$ & $2.803115 \mathrm{E}-02$ & $2.325646 \mathrm{E}-02$ \\
7 & $3.853200 \mathrm{E}-02$ & $2.347606 \mathrm{E}-02$ & $7.441780 \mathrm{E}-03$ & $2.325998 \mathrm{E}-02$ & $2.278377 \mathrm{E}-02$ \\
8 & $3.431037 \mathrm{E}-02$ & $2.281841 \mathrm{E}-02$ & $1.037963 \mathrm{E}-02$ & $2.280210 \mathrm{E}-02$     \\
9 & $2.746944 \mathrm{E}-02$ & $2.277026 \mathrm{E}-02$ & $1.551056 \mathrm{E}-02$ &      \\
\hline
\end{tabular}
\end{center}

[2009]. The higher order finite element LSFEM results show very good agreement with the reference solutions. Furthermore, the kinetic energy and the enstrophy results of the V-V-P formulation show better mesh convergence for higher $\alpha$ (Table 6 and Table 8).

Table 8: Convergence of the kinetic energy $E$ and the enstrophy $Z$ for the regularized cavity problem at $\operatorname{Re}=1000$.

\begin{center}
\includegraphics[max width=\textwidth]{2022_12_07_8d973a3879ded5997be7g-11}
\end{center}

\subsection{Solver comparisons for the S-V-P formulations}
In our third example we investigate solver comparisons between iterative and direct solvers, which are available or implemented (Pardiso 10.2.2.025) in the used softwares FEATFLOW and FEAP 8.2. For details about the solvers the reader is referred to Nickaeen et al. [2014], Turek [1999], Schenk and Gartner [2004] and Taylor [2008].

The first benchmark problem is the regularized driven cavity at $\mathrm{Re}=1$ from section 3.2. Here, we compare a conjugate gradient solver (PCG) with preconditioning using the diagonal of the matrix, with two direct solvers for the $R T_{1} P_{3} P_{1}$ formulation with respect to capability and computation time. In the whole section the nonlinear and linear solver relative changes are kept below 1E-4. The results for different mesh levels for the three solvers are shown in Table 9. It could be seen that the direct solvers are much faster than the PCG solver. A reason therefore could be the simple but not very efficient preconditioning. Apart from the computation time, the PCG solver is capable to solve higher mesh levels. Here, the direct solvers could solve only until level 7 and 8, while the PCG solver also completes level 9.

Table 9: Relative residuals and related CPU time for different solvers in FEAP $8.2$ for the regularized driven cavity at $\operatorname{Re}=1$.

\begin{center}
\includegraphics[max width=\textwidth]{2022_12_07_8d973a3879ded5997be7g-12}
\end{center}

In the second case study we compare the before mentioned PCG solver (FEAP 8.2) for the $R T_{1} P_{3} P_{1}$ formulation with the MPCG solver (FEATFLOW), see Nickaeen et al. [2014], used in the last examples for the $Q_{2} \mathrm{~S}-\mathrm{V}-\mathrm{P}$ formulation. The tested benchmarks are the flow around cylinder from section $3.1$ and the regularized driven cavity at $\operatorname{Re}=1$ from section 3.2. In Table 10 the results for regularized driven cavity are depicted. Since for $\operatorname{Re}=1$ the problem exhibits only a weak nonlinearity, there are only two Newton steps necessary to solve the problem. It could be seen that the MPCG solver needs much less linear solver iterations than the PCG solver, which is not surprising due to the above mentioned reasons. Furthermore, the PCG solver shows a mesh level dependent iteration behavior, i.e. the iteration number increases with higher mesh levels. In contrast to that the MPCG solver shows a grid-independent solution behavior. The same findings are obtained for the flow around cylinder with the results shown in Table 11. Due to the higher nonlinearity, more Newton steps are required in this example. Altogether, the MPCG solver clearly outperforms the PCG solver and for future work this technology will be analyzed also for the $H(\operatorname{div})$ based S-V-P formulation.

Table 10: The number of nonlinear Newton iterations and the corresponding averaged number of linear solver iterations for the regularized driven cavity at $\mathrm{Re}=1$ for different S-V-P formulations.

\begin{center}
\includegraphics[max width=\textwidth]{2022_12_07_8d973a3879ded5997be7g-13}
\end{center}

Table 11: The number of nonlinear Newton iterations and the corresponding averaged number of linear solver iterations for flow around cylinder at $\mathrm{Re}=20$ for different $\mathrm{S}-\mathrm{V}-\mathrm{P}$ formulations.

\begin{center}
\begin{tabular}{ccc}
\hline
 & $R T_{1} P_{3} P_{1}$ with PCG & $Q_{2}$ with MPCG \\
Level & Solver (FEAP 8.2) & Solver (FEATFLOW) \\
\hline
2 & $11 / 1729$ & $6 / 12$ \\
3 & $9 / 2914$ & $6 / 12$ \\
4 & $10 / 3541$ & $6 / 12$ \\
\hline
\end{tabular}
\end{center}

\section{Conclusion}
We presented a numerical study regarding the accuracy and the efficiency of two leastsquares finite element formulations of the incompressible Navier-Stokes equations. The first-order system for the first formulation is introduced using the stress, velocity and pressure, known as the S-V-P formulation, and for the second formulation using the vorticity, velocity and pressure, known as the V-V-P formulation. We investigated different continuous nodal finite element spaces of higher and low order for both S-V-P and V-V-P formulations. In addition, the $H$ (div)-conforming Raviart-Thomas finite elements were used for the S-V-P formulation. In combination with the Newton technique to treat the nonlinearity the resulting linear system is solved by a standard solver for the $H$ (div)-conforming Raviart-Thomas approximations, while, for the nodal finite elements an extended multigrid-preconditioned conjugate gradient solver is used.

Two incompressible steady-state laminar flow problems are studied. In the flow around cylinder test case, the flow accuracy and the mass conservation of the LSFEM formulations are investigated. In the lid-driven cavity test case, the results are analyzed with the help of global quantities, namely the kinetic energy and the enstrophy. We summarize the numerical results as follows:

\begin{enumerate}
  \item The mass conservation is investigated on one hand with respect to the different formulations and on the other hand with respect to the order of the interpolation functions. The results show that the S-V-P formulation delivers better mass conservation as compared to the V-V-P formulation. The mass conservation of the V-V-P system is enhanced with the help of an extra weighting parameter, but this comes with a negative side effect on the linear solver performance. We observe that in both formulations using higher order finite elements effectively improves the mass conservation.

  \item The linear MPCG solver performs efficiently for both of the LSFEM formulations with continuous nodal finite elements. We have obtained a grid-independent solver behavior with low and higher order finite elements for the V-V-P as well as the S-V-P formulations. However, the MPCG solver outperforms for the V-V-P system in comparison to the S-V-P system.

\end{enumerate}

We conclude that, highly accurate results are obtained with higher order finite elements. More importantly, we have obtained more accurate results with the higher-order finite elements with less number of degrees of freedom as compared to the lower-order elements. This obviously leads to less computational costs. Thus, the accuracy of the V-V-P and the S-V-P LSFEM formulations depends mainly on the order of the interpolation functions, which will also be investigated in the future for the three-dimensional case and for non-Newtonian fluid flow. Regarding the efficiency aspects, the MPCG solver performs efficiently for both of the LSFEM formulations with continuous nodal finite elements. As mentioned in subsection $3.3$ the next step will also be a multigrid solver not only for $H^{1}$ but also for $H(\operatorname{div})$ finite elements. The availability of such elements combined in one program system seems very promising. It will allow for numerical investigations and comparisons (e.g. with respect to CPU time) between all formulations and is the subject of our ongoing research.

\section{References}
P.B. Bochev and M.D. Gunzburger. Least-squares methods for the velocity-pressurestress formulation of the Stokes equations. Computer Methods in Applied Mechanics and Engineering, 126:267-287, 1995.

P.B. Bochev and M.D. Gunzburger. Least-Squares Finite Element Methods. SpringerVerlag, New York, 1st edition, 2009.

P.B. Bochev, T. Manteuffel, and S. McCormick. Analysis of velocity-flux least squares methods for the Navier-Stokes equations, part II. SIAM Journal on Numerical Analysis, $36: 1125-1144,1999$.

P.B. Bochev, J. Lai, and L. Olson. A non-conforming least-squares finite element method for the velocity-vorticity-pressure stokes equations. International Journal for Numerical Methods in Fluids, 72:375-402, 2013.

P. Bolton and R.W. Thatcher. On mass conservation in least-squares methods. J. Comput. Phys., 203:287-304, 2005.

P. Bolton and R.W. Thatcher. A least-squares finite element method for the Navier-Stokes equations. J. Comput. Phys., 213:174-183, 2006.

C.-H. Bruneau and M. Saad. The $2 \mathrm{D}$ lid-driven cavity problem revisited. Computers and Fluids, 35:326-348, 2006.

Z. Cai, B. Lee, and P. Wang. Least-squares methods for incompressible newtonian fluid flow: linear stationary problems. SIAM Journal on Numerical Analysis, 42:843-859, 2004.

C. Chang, S.Y Yang, and C.H. Hsu. A least-squares finite element method for incompressible flow in stress-velocity-pressure version. Computer Methods in Applied Mechanics and Engineering, 128:1-9, 1995.

C. L. Chang and J. J. Nelson. Least-squares finite element method for the Stokes problem with zero residual of mass conservation. SIAM Journal on Scientific Computing, 34(2): 480-489, 1997.

H. Damanik, J. Hron, A. Ouazzi, and S. Turek. A monolithic FEM-multigrid solver for non-isothermal incompressible flow on general meshes. J. Comput. Phys., 228:3869$3881,2009$.

J. M. Deang and M. D. Gunzburger. Issues related to least-squares finite element methods for the Stokes equations. SIAM Journal on Scientific Computing, 20(3):878-906, 1998.

X. Ding and T.T.H. Tsang. On first-order formulations of the least-squares finite element method for incompressible flows. International Journal of Computational Fluid Dynamics, 17:183-197, 2003.

J. J. Heys, T. A. Manteuffel, S. F. McCormick, and L. N. Olson. Algebraic multigrid for higher-order finite elements. J. Comput. Phys., 204(2):520-532, 2005. J. J. Heys, E. Lee, T. A. Manteuffel, and S. F. McCormick. On mass-conserving leastsquares methods. SIAM Journal on Scientific Computing, 28(3):1675-1693, 2006.

J.J. Heys, E. Lee, T.A. Manteuffel, and S.F. McCormick. An alternative least-squares formulation of the Navier-Stokes equations with improved mass conservation. J. Comput. Phys., 226:994-1006, 2007.

J.J. Heys, E. Lee, T. Manteuffel, S. McCormick, and J. Ruge. Enhanced mass conservation in least-squares methods for navier-stokes equations. SIAM J. Sci. Comp., 31:23032321, 2009.

B.-N. Jiang. The Least-Squares Finite Element Method. Springer-Verlag, Berlin, 1998.

O. Kayser-Herold and H.G. Matthies. Least-squares FEM, literature review. InformatikBericht 2005-05, TU Braunschweig, Institut für Wissenschaftliches Rechnen, 2005.

M. Nickaeen, A. Ouazzi, and S. Turek. Newton multigrid least-squares FEM for the V-V-P formulation of the Navier-Stokes equations. J. Comput. Phys., 256:416-427, 2014

J. P. Pontaza. A least-squares finite element form for unsteady incompressible flows with improved velocity-pressure coupling. Journal of Computational Physics, 217:563-588, 2006.

J. P. Pontaza. Least-squares variational principles and the finite element method: theory, form, and model for solid and fluid mechanics. PhD thesis, Texas A\&M University, 2003.

J. P. Pontaza and J. N. Reddy. Spectral $/ h p$ least-squares finite element formulation for the Navier-Stokes equation. Journal of Computational Physics, 190:523-549, 2003.

J. P. Pontaza and J. N. Reddy. Space-time coupled spectral/hp least-squares finite element formulation for the incompressible Navier-Stokes equation. Journal of Computational Physics, 197:418-459, 2004.

M. M. J. Proot and M. I. Gerritsma. Mass- and momentum conservation of the leastsquares spectral element method for the Stokes problem. Journal of Scientific Computing, 27:389-401, 2006.

R. Ranjan and J.N. Reddy. On multigrid methods for the solution of least-squares finite element models for viscous flows. International Journal of Computational Fluid Dynamics, 26:45-65, 2012.

O. Schenk and K. Gartner. Solving unsymmetric sparse systems of linear equations with pardiso. J. of Future Generation Computer Systems, 20:475-487, 2004.

A. Schwarz and J. Schröder. A mixed least-squares formulation of the Navier-Stokes equations for incompressible Newtonian fluid flow. Proceedings in Applied Mathematics and Mechanics, 11:589-590, 2011.

A. Schwarz, J. Schröder, and G. Starke. Least-squares mixed finite elements for small strain elasto-viscoplasticity. International Journal for Numerical Methods in Engineering, 77:1351-1370, 2009. A. Schwarz, J. Schröder, and G. Starke. A modified least-squares mixed finite element with improved momentum balance. International Journal for Numerical Methods in Engineering, 81:286-306, 2010.

S. Serdas, A. Schwarz, J. Schröder, S. Turek, A. Ouazzi, and M. Nickaeen. Influence of higher interpolation orders in mixed lsfem for the incompressible Navier-Stokes equations. Proceedings in Applied Mathematics and Mechanics, 13(1):301-302, 2013.

R. L. Taylor. Documentation Version 8.2. Robert L. Taylor, Department of Civil and Environmental Engineering, University of California at Berkeley, 2008.

S. Turek. Efficient Solvers for Incompressible Flow Problems: An Algorithmic and Computational Approach. Springer, Berlin, 1999.

S. Turek and M. Schäfer. Benchmark computations of laminar flow around cylinder. In E.H. Hirschel, editor, Flow Simulation with High-Performance Computers II, volume 52 of Notes on Numerical Fluid Mechanics, pages 547-566. Vieweg, 1996.


\end{document}