# Bubble-Enriched Least-Squares Finite Element Method for Transient Advective Transport

Authors: Rajeev Kumar; Brian H. Dennis
Journal: Differential Equations and Nonlinear Mechanics
 (2008)

---

## Abstract

The least-squares finite element method (LSFEM) has received increasing attention in recent years due to advantages over the Galerkin finite element method (GFEM). The method leads to a minimization problem in theL2-norm and thus results in a symmetric and positive definite matrix, even for first-order differential equations. In addition, the method contains an implicit streamline upwinding mechanism that prevents the appearance of oscillations that are characteristic of the Galerkin method. Thus, the least-squares approach does not require explicit stabilization and the associated stabilization parameters required by the Galerkin method. A new approach, the bubble enriched least-squares finite element method (BELSFEM), is presented and compared with the classical LSFEM. The BELSFEM requires a space-time element formulation and employs bubble functions in space and time to increase the accuracy of the finite element solution without degrading computational performance. We apply the BELSFEM and classical least-squares finite element methods to benchmark problems for 1D and 2D linear transport. The accuracy and performance are compared.


---

## Body

## 1. Introduction

In an age of  increasing atmospheric pollutions, air
pollution  modeling is getting increasingly important. Air pollution models are
generally based on atmospheric advection-diffusion equation. Major part of
uncertainty in the model predictions is due to the presence of first-order
advective transport term which causes serious numerical difficulties. However,
the nature of difficulties seems to be substantially different in steady and
unsteady advection.In steady state advection problems, the difficulty in
the form of oscillations or wiggles is a consequence of negative (numerical)
diffusion that is inherent in use of centered type discretization for the
convective terms. This applies to central finite difference method as well as
the closely relate	d Galerkin finite element method (GFEM), both leading to a
nonsymmetric, nonpositive definite matrices as Jiang has illustrated in his
text [1]. These asymmetric matrices give rise to odd even decoupling, which
causes node-to-node oscillations in the solution. This can be tackled by severe
refinement of the mesh that greatly undermines the utility of the scheme.Numerical difficulties of different types are
encountered in the time-dependent advection problems. Transient convection
problems are governed by hyperbolic differential equations. The characteristic
lines now assume great importance. The discretization in space now influences
discretization in time and vice versa as they are now interlinked through the
characteristics. One can circumvent the issue by resorting to a Lagrangian
(moving coordinates) formulation in which the convective term vanishes. However,
the formulation is difficult and thus not very popular. The popular Eulerian
formulation, therefore, must properly accommodate the flow physics of
information propagation along the characteristic line, while discretizing in
space and time.Over the years, the Galerkin method in form of its
variants has been used extensively to solve convection problems. Classical GFEM
is very dispersive in nature due to inherent generation of the negative
diffusion. Its popular variant Petrov-Galerkin provides stabilized solutions by
generating numerical diffusion. Petrov-Galerkin method using higher degree
polynomial as weighting function (Christie et al. [2];
Westerink and  Shea [3])
and the streamline upwind Petrov-Galerkin method (SUPG) by Brooks and Hughes
[4] both have at least one free parameter or an intrinsic time function that
has to be tuned in order to control the amount of artificial diffusion. This is
the disadvantage of Petrov-Galerkin methods.Donea [5] proposed Taylor-Galerkin (TG) method, where
Taylor
series
for time discretisation is used before applying space discretisation. The
resulting Taylor-Galerkin methods do not introduce any free parameter but they
require the use of higher-order derivatives.LSFEM which is based  on minimizing theL2-norm
of the residuals is naturally suited for a first order system of differential
equations. Unlike GFEM, LSFEM formulation leads to symmetric positive definite
(SPD) matrices that can be effectively solved using matrix-free iterative
methods like preconditioned conjugate gradient method.Jiang and  Povinelli [6] pointed out the advantages of
LSFEM by demonstrating and validating the method for a variety of compressible
and incompressible flow problems. Jiang et al. [7] also developed a matrix-free
LSFEM for three-dimensional, steady state lid-driven cavity flow.Donea and Quartapelle [8] classified the following
four different least square finite element approaches: the LSFEM proposed by
Carey and  Jiang [9] based on Crank-Nicolson approximation across the time step;
characteristic LSFEM by Li [10];
Taylor-LSFEM by Park and Liggett [11, 12]; and
space-time finite element method, STLSFEM by Nguyen and Reynen [13]. The first
three approaches rely on a quadratic functional associated with time
discretized version of governing
equation, whereas the last one extends the least square formulation and
its finite element representation to space-time domain. Donea and Quartapelle
pointed out that the LSFEM proposed by Carey and Jiang [9] was the most
interesting least square method for advective transport problems presumably
because of simplicity of its formulation and accuracy, and its close
relationship with the SUPG, Galerkin least square (Hughes et al. [14]), and
Taylor Galerkin method. They also found the space-time LSFEM very inaccurate
and diffusive; therefore, not worth recommending for advective transport
problems.The numerical difficulties faced in the form of “wiggles”
can be tackled by resorting to severe mesh refinement which forces the use of
very small time steps, thereby undermining the utility of GFEM. In a study,
Surana and  Sandhu [15] have demonstrated that these oscillations can be
completely
eliminated by using p-version of STLSFEM, where they have used p-values as high
as 7 in space and 11 in time to completely recover the exact solution even
after convecting the Gaussian distribution profile to some distance in the
domain. But the p-version, especially in 2- and 3-dimensional problems, becomes
computationally very expensive and difficult to program.In the present work, we have used space-time LSFEM
with linear elements enriched with bubble modes to get reasonably accurate
solutions to advective transport equation without resorting to severe mesh
refinement andp-version of LSFEM. We term this approach the bubble-enriched least-squares
finite element method (BELSFEM). The Space-time LSFEM as described by
Donea and Quartapelle [8] is second-order accurate and unconditionally stable. Results from
STLSFEM applied to pure advection problems are less accurate and more
dissipative compared to the one obtained from LSFEM using Crank-Nicolson time
discretization. Notwithstanding that STLSFEM has  been chosen as it has finite element
discretization both in space and time domains essential
for applying bubble modes. Results were also generated using Crank-Nicolson
LSFEM proposed by Carey and Jiang, deemed most interesting by  Donea and Quartapelle in
their 1992 article, in order to be used as baseline for comparison.


## 2. The Least-Square Finite Element Method

Consider the transient advection equation given as(2.1)∂U∂t+(V→•∇)U=0,where U is the
property being convected at a velocity V→ with u,v, and w as its
components in x,y, and z directions,
respectively. To illustrate the main benefits of LSFEM, consider the
application of a simple least-squares finite element method to the transient
advection equation. Before application of the finite element method in space,
the time derivative of (2.1) is discretized with a simple backward-Euler method: (2.2)Un+1−UnΔt+v→⋅∇Un+1=0.In the least-squares approach, the L2-norm
of the differential equation is minimized with respect to unknown coefficients
over the solution domain Ω.
Applying the L2-norm to (2.2) and minimizing the functional with
respect to Un+1 leads to the weak statement (2.3)∫Ω({N}Δt+(v→⋅∇){N})({N}Δt+(v→⋅∇){N})TdΩ{Un+1}=∫Ω({N}Δt+(v→⋅∇){N}){N}TΔt{Un}dΩ, where
the row vector {N} contains the
basis functions Nj
used to approximate the solution over the domain as U(x,y,z)=∑jNj(x)Uj={N}T{U}.The
weak statement can be expanded and written in matrix form(2.4)([M]Δt+([C]+[C]T)+Δtv→Tv→[K]){Un+1}=([M]Δt){Un},where the individual matrix contributions are given
by (2.5)[M]=∫Ω{N}{N}TdΩ,[C]=∫Ω{N}{(v→⋅∇)N}TdΩ,[K]=∫Ω[∇N][∇N]TdΩ.
Equation
(2.4) clearly shows
that the resulting system of equations is symmetric, a quality that is not
achievable for Galerkin finite element methods or even finite difference or
finite volume methods. In addition, one can notice an upwind diffusion term
that is implicit to the least-squares approach. The upwind diffusion is often
useful for smoothing nonmonotone solutions that occur before and after any
sharp gradients that appear in the flow direction. We also wish to emphasize
that there are no tunable parameters in the LSFEM approach, such parameters
often appear in stabilized Galerkin methods and are difficult to determine in
general.


## 3. The Least-Square Finite Element Formulations

For the sake of simplicity, let us consider 1D scalar
advection equation(3.1)∂U∂t+a∂U∂x=0.The three least-square finite element formulations
tried are as follows.


### 3.1. Crank-Nicolson LSFEM

In least-square finite element formulation, we
minimize the square of the residual,R, given by R=∂U˜/∂t+a(∂U˜/∂x), where U˜ is the
approximate solution. For sake of simplicity, we will use U in place of U˜.
The LSFEM formulation based on minimization of
square of residual leads to (3.2)∂∂Un+1∫Ω(∂U∂t+a∂U∂x)2dxdt≈0.Using forward difference for time derivative term and
θ-method for approximating U in
convective term gives (3.3)∂∂Un+1∫Ω(Un+1−UnΔt+ad(θUn+1+(1−θ)Un)dx)2=0.Let the
unknown U be defined as (3.4)U(x)=∑jNj(x)Uj,where Uj
is the solution at the jth node and Nj is the interpolation function. Taking the derivative with
respect to Un+1, (3.3) leads to the Crank-Nicolson LSFE formulation (3.5)∑i∫Ω{Ni(x)+aΔtθdNi(x)dx}{Ni(x)+aΔtθdNi(x)dx}TUin+1dx=∑i∫Ω{Ni(x)+aΔtθdNi(x)dx}{Ni(x)−aΔt(1−θ)dNi(x)dx}TUindx.  For θ=1/2, it becomes
Crank-Nicolson LSFEM formulation as (3.6)∑i∫Ω{Ni(x)+aΔt2dNi(x)dx}{Ni(x)+aΔt2dNi(x)dx}TUin+1dx=∑i∫Ω{Ni(x)+aΔt2dNi(x)dx}{Ni(x)−aΔt2dNi(x)dx}TUindx.


### 3.2. Space-Time LSFEM

In space-time formulation, both time and space
derivatives are discretized the finite element way and the unknownU becomes function of both
spatial and temporal
variables, that is, (3.7)U(x,t)=∑jNj(x,t)UjorU(x,y,t)=∑jNj(x,y,t)Uj,where Nj(x,t) is bilinear interpolation function for 1D and
Nj(x,y,t) is the
trilinear interpolation function for 2D formulation. Equations (3.2) and
(3.7)
lead to simple space-time least square finite element formulation (3.8)∑i∫Ω{∂Ni(x,t)∂t+a∂Ni(x,t)∂x}{∂Ni(x,t)∂t+a∂Ni(x,t)∂x}TUin+1dxdt=0.Linear elements of 1D domain transform to 2D bilinear
elements and 2D quadrilateral element transform to trilinear elements in the
space-time formulation. For bilinear elements, the bilinear shape functions are
given in terms of natural  coordinates by
(3.9a)N(ξ,τ)={L1(ξ)L1(τ),L2(ξ)L1(τ),L2(ξ)L2(τ),L1(ξ)L2(τ)}T.Similarly, Trilinear shape functions for trilinear
elements are given by (3.9b)N(ξ,η,τ)={L1(ξ)L1(η)L1(τ)L2(ξ)L1(η)L1(τ)L2(ξ)L2(η)L1(τ)L1(ξ)L2(η)L1(τ)L1(ξ)L1(η)L2(τ)L2(ξ)L1(η)L2(τ)L2(ξ)L2(η)L2(τ)L1(ξ)L2(η)L2(τ)},
where L1(ξ)=(1/2)(1−ξ),L2(ξ)=(1/2)(1+ξ),L1(η)=(1/2)(1−η),L2(η)=(1/2)(1+η),L1(τ)=(1/2)(1−τ), and L2(τ)=(1/2)(1+τ) are the linear shape functions and ξ=x/Δx,η=y/Δy and τ=t/Δt the natural coordinates.


### 3.3. Bubble-Enriched LSFEM

Since space-time formulation has finite element
discretization for both time and space derivative it has been selected for
application of bubble modes in this work. In this approach, bubble functions
are used to enrich the function space of the finite element. We refer this new
approach as the bubble-enriched least-squares finite element method (BELSFEM).
Bubbles are the functions defined in the interiors of the finite elements that
vanish on the element boundaries. Baiocchi et al.
[16] were the first to point out that
the enrichment of the finite element space by summation of polynomial bubble
functions results in stabilized procedures for convection-diffusion problems
formally similar to SUPG and GLS. Brezzi et al. [17]
and Franca et al. [18] introduced more general framework for the
discretization of problem involving multiscale phenomena.In bubble enrichment method, we add bubble functions
to the set of nodal shape functions of the linear elements in space and time
direction and their tensor product gives the set of bilinear shape functions.
We include only the modes falling inside the bilinear element (excluding the
modes falling on the edges). Bubble functions take zero value on the element
boundaries. This property of bubble functions allows the use of classical
static condensation procedure to condense the bubble modes out and include
their effect in the basic element matrix.Bubble functions were taken from orthogonal set of
Jacobi polynomials denoted byPpα,β.
Jacobi polynomials are a family of polynomial solutions to the singular
Sturm-Liouville problem. A significant feature of these polynomials is that
they are orthogonal in the interval [−1,1] with respect to the function (1−x)α(1+x)β(α,β>−1).
Bubble modes were generated from Ppα,β as (3.10)ψp(x)=(1−ξ2)(1+ξ2)Pp−11,1(ξ),0<p,where p is the order of the Jacobi polynomial. Jacobi polynomials with α=β=1
were chosen as they produce symmetric and diagonally strong matrices for second-order
differential equations (Karniadakis and  Sherwin  [19]). First few of the Jacobi polynomials
used are shown in Figure 1. A pseudo code  outlining the whole process is shown in
Algorithm 1.Algorithm 1
Pseudo Code:(1)
Formulate and initialize STLSFEM(2)
Generate bubble fns using Jacobi
polynomials(3)
Introduce bubble fns into original set of
nodal shape function using tensor product,
element stiffness matrix
size goes up from originalm to p=m+bn. Wheren is number of dimensions.(4)
While(p≥m){/∗to get the original size of element
stiffness matrix back/∗Apply Static
Condensation()p=p−1;}(6)
Set the time limit and convect the solution
using linear solver(7) endFigure 1
First few bubble modes generated using Jacobi polynomials(Pp−11,1(x),0<p).


## 3.1. Crank-Nicolson LSFEM

In least-square finite element formulation, we
minimize the square of the residual,R, given by R=∂U˜/∂t+a(∂U˜/∂x), where U˜ is the
approximate solution. For sake of simplicity, we will use U in place of U˜.
The LSFEM formulation based on minimization of
square of residual leads to (3.2)∂∂Un+1∫Ω(∂U∂t+a∂U∂x)2dxdt≈0.Using forward difference for time derivative term and
θ-method for approximating U in
convective term gives (3.3)∂∂Un+1∫Ω(Un+1−UnΔt+ad(θUn+1+(1−θ)Un)dx)2=0.Let the
unknown U be defined as (3.4)U(x)=∑jNj(x)Uj,where Uj
is the solution at the jth node and Nj is the interpolation function. Taking the derivative with
respect to Un+1, (3.3) leads to the Crank-Nicolson LSFE formulation (3.5)∑i∫Ω{Ni(x)+aΔtθdNi(x)dx}{Ni(x)+aΔtθdNi(x)dx}TUin+1dx=∑i∫Ω{Ni(x)+aΔtθdNi(x)dx}{Ni(x)−aΔt(1−θ)dNi(x)dx}TUindx.  For θ=1/2, it becomes
Crank-Nicolson LSFEM formulation as (3.6)∑i∫Ω{Ni(x)+aΔt2dNi(x)dx}{Ni(x)+aΔt2dNi(x)dx}TUin+1dx=∑i∫Ω{Ni(x)+aΔt2dNi(x)dx}{Ni(x)−aΔt2dNi(x)dx}TUindx.


## 3.2. Space-Time LSFEM

In space-time formulation, both time and space
derivatives are discretized the finite element way and the unknownU becomes function of both
spatial and temporal
variables, that is, (3.7)U(x,t)=∑jNj(x,t)UjorU(x,y,t)=∑jNj(x,y,t)Uj,where Nj(x,t) is bilinear interpolation function for 1D and
Nj(x,y,t) is the
trilinear interpolation function for 2D formulation. Equations (3.2) and
(3.7)
lead to simple space-time least square finite element formulation (3.8)∑i∫Ω{∂Ni(x,t)∂t+a∂Ni(x,t)∂x}{∂Ni(x,t)∂t+a∂Ni(x,t)∂x}TUin+1dxdt=0.Linear elements of 1D domain transform to 2D bilinear
elements and 2D quadrilateral element transform to trilinear elements in the
space-time formulation. For bilinear elements, the bilinear shape functions are
given in terms of natural  coordinates by
(3.9a)N(ξ,τ)={L1(ξ)L1(τ),L2(ξ)L1(τ),L2(ξ)L2(τ),L1(ξ)L2(τ)}T.Similarly, Trilinear shape functions for trilinear
elements are given by (3.9b)N(ξ,η,τ)={L1(ξ)L1(η)L1(τ)L2(ξ)L1(η)L1(τ)L2(ξ)L2(η)L1(τ)L1(ξ)L2(η)L1(τ)L1(ξ)L1(η)L2(τ)L2(ξ)L1(η)L2(τ)L2(ξ)L2(η)L2(τ)L1(ξ)L2(η)L2(τ)},
where L1(ξ)=(1/2)(1−ξ),L2(ξ)=(1/2)(1+ξ),L1(η)=(1/2)(1−η),L2(η)=(1/2)(1+η),L1(τ)=(1/2)(1−τ), and L2(τ)=(1/2)(1+τ) are the linear shape functions and ξ=x/Δx,η=y/Δy and τ=t/Δt the natural coordinates.


## 3.3. Bubble-Enriched LSFEM

Since space-time formulation has finite element
discretization for both time and space derivative it has been selected for
application of bubble modes in this work. In this approach, bubble functions
are used to enrich the function space of the finite element. We refer this new
approach as the bubble-enriched least-squares finite element method (BELSFEM).
Bubbles are the functions defined in the interiors of the finite elements that
vanish on the element boundaries. Baiocchi et al.
[16] were the first to point out that
the enrichment of the finite element space by summation of polynomial bubble
functions results in stabilized procedures for convection-diffusion problems
formally similar to SUPG and GLS. Brezzi et al. [17]
and Franca et al. [18] introduced more general framework for the
discretization of problem involving multiscale phenomena.In bubble enrichment method, we add bubble functions
to the set of nodal shape functions of the linear elements in space and time
direction and their tensor product gives the set of bilinear shape functions.
We include only the modes falling inside the bilinear element (excluding the
modes falling on the edges). Bubble functions take zero value on the element
boundaries. This property of bubble functions allows the use of classical
static condensation procedure to condense the bubble modes out and include
their effect in the basic element matrix.Bubble functions were taken from orthogonal set of
Jacobi polynomials denoted byPpα,β.
Jacobi polynomials are a family of polynomial solutions to the singular
Sturm-Liouville problem. A significant feature of these polynomials is that
they are orthogonal in the interval [−1,1] with respect to the function (1−x)α(1+x)β(α,β>−1).
Bubble modes were generated from Ppα,β as (3.10)ψp(x)=(1−ξ2)(1+ξ2)Pp−11,1(ξ),0<p,where p is the order of the Jacobi polynomial. Jacobi polynomials with α=β=1
were chosen as they produce symmetric and diagonally strong matrices for second-order
differential equations (Karniadakis and  Sherwin  [19]). First few of the Jacobi polynomials
used are shown in Figure 1. A pseudo code  outlining the whole process is shown in
Algorithm 1.Algorithm 1
Pseudo Code:(1)
Formulate and initialize STLSFEM(2)
Generate bubble fns using Jacobi
polynomials(3)
Introduce bubble fns into original set of
nodal shape function using tensor product,
element stiffness matrix
size goes up from originalm to p=m+bn. Wheren is number of dimensions.(4)
While(p≥m){/∗to get the original size of element
stiffness matrix back/∗Apply Static
Condensation()p=p−1;}(6)
Set the time limit and convect the solution
using linear solver(7) endFigure 1
First few bubble modes generated using Jacobi polynomials(Pp−11,1(x),0<p).


## 4. Test Problems

Standard test problems taken in one and two dimensions
are as follows.


### 4.1. One-Dimensional Problems

#### 4.1.1. Convection of Gaussian Hill

This one-dimensional problem was taken from
Donea and  Huerta  [20]. A
Gaussian distribution profile was convected over 1D domain ]0,1[ with the
initial condition (4.1)U(x,0)=57exp{−(x−x0l)2},where x0=2/15,l=72/300,  and the boundary condition as U(0,t)=U(1,t)=0 and convection velocity a=1. The
solution was convected by t=0.6 over a uniform mesh of size h=1/150. The exact solution is given by (4.2)U(x,t)=57exp{−(x−x0−atl)2}.


#### 4.1.2. Propagation of a Steep Front

This 1D problem also taken from Donea and  Huerta [20] considers the
convection at unit speed of a discontinuous initial data. The discontinuity
occurs over one element and is initially located at position x=0.2
of the domain ]0,1[.The discontinuity is given as(4.3)U(x,0)={1ifx<0.2,0ifx≥0.2. The solution was convected by t=0.6
using a mesh of uniform size h=1/50.


### 4.2. Two-Dimensional Problems

#### 4.2.1. Convection of a Concentration Spike

A concentration spike, given  by(4.4)U(x,y,0)={exp{−[(x−0.175)2+(y−0.175)2](0.00125)}0ifU(x,0)≤10−10, was convected by t=1.3 with a velocity given by u=0.25 and v=0.1166 at an angle of 25° to the x-axis. A 40×20 mesh in
0≤x≤1,0≤y≤0.5 was used and this problem was
picked from Yu and  Heinrich [21]. Profile was convected for Courant numbers of 0.73 (same as in
Yu and  Heinrich [21]), 1.0, and 1.47.


#### 4.2.2. Rotating Cosine Hill Problem

This classical test problem for 2D convection schemes taken
from Donea and  Huerta [20] considers the convection of a product cosine hill in a pure
rotational velocity field. The initial data is given by (4.5)U(x,y,0)={14(1+cosπX)(1+cosπY)ifX2+Y2≤1,0otherwise, where X=(x−x0)/σ  and  Y=(y−y0)/σ,
and the boundary condition is U=0
on Γin. The initial positions of the center
and the radius of the cosine hill are (x0,y0)=(1/6,1/6) and σ=0.2,
respectively. The angular velocity is given by ω(x)=(−y,x).
A uniform mesh of 30×30
four-node elements over the unit square [−0.5,0.5]×[−0.5,0.5]  was used in the computations.


## 4.1. One-Dimensional Problems

### 4.1.1. Convection of Gaussian Hill

This one-dimensional problem was taken from
Donea and  Huerta  [20]. A
Gaussian distribution profile was convected over 1D domain ]0,1[ with the
initial condition (4.1)U(x,0)=57exp{−(x−x0l)2},where x0=2/15,l=72/300,  and the boundary condition as U(0,t)=U(1,t)=0 and convection velocity a=1. The
solution was convected by t=0.6 over a uniform mesh of size h=1/150. The exact solution is given by (4.2)U(x,t)=57exp{−(x−x0−atl)2}.


### 4.1.2. Propagation of a Steep Front

This 1D problem also taken from Donea and  Huerta [20] considers the
convection at unit speed of a discontinuous initial data. The discontinuity
occurs over one element and is initially located at position x=0.2
of the domain ]0,1[.The discontinuity is given as(4.3)U(x,0)={1ifx<0.2,0ifx≥0.2. The solution was convected by t=0.6
using a mesh of uniform size h=1/50.


## 4.1.1. Convection of Gaussian Hill

This one-dimensional problem was taken from
Donea and  Huerta  [20]. A
Gaussian distribution profile was convected over 1D domain ]0,1[ with the
initial condition (4.1)U(x,0)=57exp{−(x−x0l)2},where x0=2/15,l=72/300,  and the boundary condition as U(0,t)=U(1,t)=0 and convection velocity a=1. The
solution was convected by t=0.6 over a uniform mesh of size h=1/150. The exact solution is given by (4.2)U(x,t)=57exp{−(x−x0−atl)2}.


## 4.1.2. Propagation of a Steep Front

This 1D problem also taken from Donea and  Huerta [20] considers the
convection at unit speed of a discontinuous initial data. The discontinuity
occurs over one element and is initially located at position x=0.2
of the domain ]0,1[.The discontinuity is given as(4.3)U(x,0)={1ifx<0.2,0ifx≥0.2. The solution was convected by t=0.6
using a mesh of uniform size h=1/50.


## 4.2. Two-Dimensional Problems

### 4.2.1. Convection of a Concentration Spike

A concentration spike, given  by(4.4)U(x,y,0)={exp{−[(x−0.175)2+(y−0.175)2](0.00125)}0ifU(x,0)≤10−10, was convected by t=1.3 with a velocity given by u=0.25 and v=0.1166 at an angle of 25° to the x-axis. A 40×20 mesh in
0≤x≤1,0≤y≤0.5 was used and this problem was
picked from Yu and  Heinrich [21]. Profile was convected for Courant numbers of 0.73 (same as in
Yu and  Heinrich [21]), 1.0, and 1.47.


### 4.2.2. Rotating Cosine Hill Problem

This classical test problem for 2D convection schemes taken
from Donea and  Huerta [20] considers the convection of a product cosine hill in a pure
rotational velocity field. The initial data is given by (4.5)U(x,y,0)={14(1+cosπX)(1+cosπY)ifX2+Y2≤1,0otherwise, where X=(x−x0)/σ  and  Y=(y−y0)/σ,
and the boundary condition is U=0
on Γin. The initial positions of the center
and the radius of the cosine hill are (x0,y0)=(1/6,1/6) and σ=0.2,
respectively. The angular velocity is given by ω(x)=(−y,x).
A uniform mesh of 30×30
four-node elements over the unit square [−0.5,0.5]×[−0.5,0.5]  was used in the computations.


## 4.2.1. Convection of a Concentration Spike

A concentration spike, given  by(4.4)U(x,y,0)={exp{−[(x−0.175)2+(y−0.175)2](0.00125)}0ifU(x,0)≤10−10, was convected by t=1.3 with a velocity given by u=0.25 and v=0.1166 at an angle of 25° to the x-axis. A 40×20 mesh in
0≤x≤1,0≤y≤0.5 was used and this problem was
picked from Yu and  Heinrich [21]. Profile was convected for Courant numbers of 0.73 (same as in
Yu and  Heinrich [21]), 1.0, and 1.47.


## 4.2.2. Rotating Cosine Hill Problem

This classical test problem for 2D convection schemes taken
from Donea and  Huerta [20] considers the convection of a product cosine hill in a pure
rotational velocity field. The initial data is given by (4.5)U(x,y,0)={14(1+cosπX)(1+cosπY)ifX2+Y2≤1,0otherwise, where X=(x−x0)/σ  and  Y=(y−y0)/σ,
and the boundary condition is U=0
on Γin. The initial positions of the center
and the radius of the cosine hill are (x0,y0)=(1/6,1/6) and σ=0.2,
respectively. The angular velocity is given by ω(x)=(−y,x).
A uniform mesh of 30×30
four-node elements over the unit square [−0.5,0.5]×[−0.5,0.5]  was used in the computations.


## 5. Calculation of Flow Parameters

Important flow parameter, Courant number, is given asC=∥u∥(Δt/h),
where u is the convection
velocity, Δt is the
time step, and h is the characteristic length in the direction of the
convection. In one-dimensional problems, h is simply taken as h=Δx and ∥u∥=a. In the first problem of convection
of Gaussian hill, Δx=1/150 and
in the second problem of propagation of discontinuity Δx=1/50 was taken. Different values of Courant number were obtained by varying
Δt values.For the 2D test problems, the flow parameters were
calculated as done in the source papers. For the concentration spike test
problem,h was calculated as (5.1)h=1∥u∥(|u|Δx+|v|Δy),where u=ui+vj is the velocity vector and Courant number was given
as (5.2)C=(|u|Δx+|v|Δy)Δt.For the second test problem, since the flow field is
rotational, the velocity is changing throughout the cone; therefore, the
Courant number based on the velocity at the peak of the cone is given by ωrpeak, where
ω is the angular velocity.


## 6. Results and Discussion

The least-squares  methods previously described were
implemented in C++ on uniform quadrilateral and hexahedral meshes. Integration
was performed using Gaussian quadrature. A sparse matrix data structure was
used to conserve memory. Linear systems of equations were solved efficiently
using a Jacobi preconditioned conjugate gradient (PCG) method. An  absolute
tolerance of1.0E−6 was used for all PCG iterations. Inaccurate results of
STLSFEM were considerably improved by introduction of bubble functions. Results
improved gradually with increase in number of bubble functions until a number
beyond which the effect seems to saturate. Results for the number of bubble
functions giving best performance have been discussed.


### 6.1. One-Dimensional Problems

#### 6.1.1. Convection of Gaussian Hill

Results of the Gaussian hill problem are presented in
Figure2 and Table 1. The initial profile shown in dotted line was propagated
till t=0.6, for three Courant numbers of 0.5, 1.0, and
1.5. All the results
have been compared with results from Crank-Nicolson LSFEM as baseline. Results
of the space-time LSFEM are far more dissipative and dispersive compared to the
Crank-Nicolson LSFEM for all the three Courant numbers. However, results show
significant improvement with BELSFEM.Table 1
Convection of Gaussian hill byt=0.6.
CN-LSFEMST-LSFEMBE-LSFEMCourant no.UminUmaxUminUmaxUminUmax%redn.
inUmin%gain
inUmax0.5−0.00550.6861−0.01860.6784−
0.0013†0.6967†76.91.551.0−0.04900.6606−0.10040.6196≈00.7140≈1008.081.5−0.11960.6210−0.15360.5532−0.10490.640112.23.1†with one bubble in
both x and t.Propagation
of Gaussian hill by timet=0.6
for Courant numbers, C=0.5 (a),C=1.0
(b) and C=1.5
(c)
for continuous LSFEM.
(a)(b)(c)For Courant  number,C=0.5, BELSFEM with one bubble
in x and
t direction gives
1.5% increase in maximum
value and 77% decrease in dispersion error compared to Crank-Nicolson
LSFEM. More than one bubble in fact degraded the results.ForC=1.0, 8 bubbles in
x and 10 in time completely remove
the dispersion
error and increase
the peak by around 8% leading to complete recovery of the exact solution.ForC=1.5, BELSFEM with 8 and 10 bubbles in x and t, respectively, causes 12.2% reduction in dispersion error and about 3% increase
in the peak value.


#### 6.1.2. Propagation of a Steep Front

Discontinuity was propagated byt=0.6
and the results  presented in Figure 3 and
Table 2 were computed for Courant
numbers of 0.75, 1.0, and 2.0. Few parameters were considered for comparative
quantification of the results. Slope, m, of the solution at the discontinuity
which indicates the amount of dissipation in the solution was measured across
the two nodes that capture the discontinuity in exact solution. Since the
discontinuity spanned one element (h=1/50), the exact solution had a slope, m=−50.
Also considered were the values of Umax
and Umin
causing the overshoot and undershoot representatives
of the dispersive error. All the comparative results were based on the results
from Crank-Nicolson LSFEM.Table 2
Propagation
of discontinuity byt=0.6.
CN-LSFEMST-LSFEMBE-LSFEMCourant no.slopemUminUmaxSlopemUminUmaxSlopemUminUmax%gain inm%redn. inUmin%redn.
inUmax0.75−12.66−0.00051.1341−9.78901.1740−14.64−0.1791.000115.6−356.611.81.0−10.33none1.1684−7.9650.00011.193−48.3101.0109367.7—13.52.0−5.947none1.2934−4.9070.00541.2232−5.6110.00251.245−5.65—3.8Propagation of a steep front by timet=0.6 for Courant numbers, C=0.75 (a),
C=1.0
(b), and C=2.0
(c) for continuous
LSFEM.
(a)(b)(c)Space-time LSFEM is more dissipative than CNLSFEM for
all the three Courant numbers as can be seen in Figure3. However, it is more
dispersive than Crank-Nicolson LSFEM for C=0.75 and 1.0 and less dispersive
for C=2.AtC=0.75, BELSFEM with 8 bubbles in
x and 10
bubbles in time causes 15.6% increase in the slope (meaning reduced dissipative
error) but a large increase in dispersive error in the form of a deep
undershoot. Although results are much better with one bubble each in x and t directions with 40% increase in the slope and much smaller undershoot, as can
be seen in Figure 3.AtC=1.0, the 8/10
bubble combination shows a
significant improvement in the results as slope m reaches very close to the
exact value of −50 (see Table 2) and the dispersion error completely disappears
and the solution looks almost like the exact solution (see Figure 3).AtC=2.0, BELSFEM fails to better the slope of
Crank-Nicolson LSFEM, although it is less dispersive.


### 6.2. Two-Dimensional Problems

#### 6.2.1. Convection of a Concentration Spike

The concentration spike was convected linearly byt≈1.3 at a
unit velocity given by u=0.25 and v=0.1177 and making an angle of 25° with the
x-direction for Courant numbers of 0.73, 1.0, and
1.47. Results are presented in Figures 4, 5, and
Table 3.
Figure 4 presents the
variation of maximum and minimum concentrations with time and
Figure 5 shows
typical plot of concentration profile before and after being convected. For all
the Courant numbers, tested Space-time LSFEM is far more dissipative and
dispersive compared to Crank-Nicolson LSFEM (see
Figures 4, 5, and
Table 3).
However, there is a marked improvement in the results with bubbles. In
addition, the maximum number of PCG iterations per time step required to
achieve tolerance remained consistent as the number of bubble functions was
increased as shown in Table 3. This clearly indicates the ability of the
BELSFEM to increase accuracy without dramatically increasing computational
effort.Table 3
Advection
of concentration spike byt=1.3.
CN-LSFEMST-LSFEMBE-LSFEMImprovementsRange of
PCG iterations*Courant no.UminUmaxUminUmaxUminUmax%redn.
inUmin%gain
inUmax0.73−0.06950.6119−0.08650.5692−0.05310.694123.613.46-71.0−0.08430.5663−0.10790.5143−0.04860.678042.319.735-61.47−0.13910.5164−0.12500.4462−0.10840.570422.110.46-7Variation of maximum and minimum
concentrations with time for advection of concentration spike : comparison of results over the time of advection.
(a)(b)(c)(d)(e)(f)Convective transport of the concentration spike (initial condition shown by the left cone) with flow at25° to
x-axis for C=1.0.
(a)(b)(c)AtC=0.73 (the same C used by Yu and  Heinrich [21]
in convecting the same profile with Petrov-Galerkin formulation), BELSFEM with
6 bubbles each in spatial and time directions results in 23.6% increase in
Umax
and
13.4% decrease in Umin
compared to Crank-Nicolson LSFEM.Results further improved forC=1.0 as 42.3%
increase in Umax
and 20%
decrease in Umin
accrued (see Figures 4, 5, and
Table 3). And finally
for C=1.47, about
22% increase in Umax
and
10.4% decrease in Umin
were recorded.


#### 6.2.2. Rotating Cosine Hill Problem

Results for rotating cosine hill problem are shown in
Figures6, 7
and Table 4. The variation of maximum and minimum values of
concentration over one rotation for ▵t=2π/120,2π/60, and 2π/30 is shown in
Figure 6. A typical profile after
one rotation is shown for the three formulations in Figure 7. Again,
Crank-Nicolson LSFEM serves as the baseline for comparison.Table 4
Advection
of cosine hill in rotation.
CN-LSFEMST-LSFEMBE-LSFEMImprovementsRange of
PCG iterations*ΔtCourant no.‡UminUmaxUminUmaxUminUmax%redn.
inUmin%gain
inUmax2π/1200.2618−0.02650.9691−0.04050.958−0.01890.976928.80.817-82π/600.5236−0.06150.9165−0.11380.8872−0.02700.971356.15.988-82π/301.047−0.20090.8369−0.21920.7398−0.21290.8418−5.970.615-16‡Courant number
based on velocity of the peak.
*Range of PCG iterations/time step: number of iterations with one bubble—the number with six bubbles.
% Reduction and %
gain calculated on Crank-Nicolson LSFEM results as baseline.Variation of maximum and minimum concentrations with time for advection of cosine hill in rotation : comparison of results over one rotation.
(a)(b)(c)(d)(e)(f)Convection of a cosine hill in a pure rotational velocity
field with▵t=2π/60 : comparison of results after a complete revolution.
(a)(b)(c)For▵t=2π/120, BELSFEM with 6 bubbles each in spatial
and time directions shows about 29% reduction in dispersive error and about 1%
increase in the peak value. This improvement in the peak value is significant
considering the fact that the baseline value from Crank-Nicolson LSFEM itself
was high at 0.9691 (see Table 4).For▵t=2π/60, there is more improvement in the results
as the dispersion error declines by 56% and the peak value goes up by around
6%. Typical profiles after one rotation for this case are shown in Figure 7.For▵t=2π/30 (which corresponds to C≈1, based on velocity at the peak of the
profile), however, there is only 0.6% improvement in peak value and the
dispersive error is worse than CNLSFEM, as can be seen in Figure 6.


### 6.3. Effect of Mesh Size and Number of Bubbles

Two-dimensional benchmarks problems were run on
different sizes of mesh and also on basic meshes with different number of
bubble functions in order to investigate the effect of mesh size and number of
bubble functions on the performance of BELSFEM. Mesh size parameter,h,
was varied from 0.01 to 0.1 (where h = side-length/number of elements per side). In all the cases,
h in x and y directions was the same.Typical  comparative plots ofUmax
and Umin
from  the three least-square methods for different mesh sizes are shown in
Figure 8. For this part of study, four bubbles each in space and time were used. The
maximum and minimum values for the cosine hill are recorded after one full
rotation and those for linear convection of concentration spike have been taken
after being convected by t=1.3. The bubbles seem to be most effective for
moderately coarse meshes as can be observed from the figure where large gain over
both CNLSFEM and STLSFEM can be seen in this region. However, for very coarse
and very fine meshes the benefits of bubbles seem to diminish.Effect of mesh sizeh on the performance of BELFEM and LSFEM.
(a)(b)Figure9 shows the effect of number of bubbles on the
performance of BELSFEM. Typical
variation of Umax
and Umin
for the two problems
is displayed. Results improve sharply with the number of bubble functions
initially but the improvements diminish with further increase in the number and
beyond 3-4 bubbles the effect saturates. It, therefore, can be stated that
generally good improvements in the results can be achieved with 4–6 bubbles.Effect of number of bubbles on the performance of BELFEM. (*pure LS method-results (not functions of n)).
(a)(b)These results show the clear benefit of bubble
functions for linear transport problems, which are purely hyperbolic in nature.
Extensions of this work to mixed problems, such as Navier-Stokes equations, are
of great practical interest and a topic of further research. In addition, there
likely exist optimal bubble functions that will achieve highly accurately
solutions with a small number of functions. The form of these functions is also
a topic of further research.


## 6.1. One-Dimensional Problems

### 6.1.1. Convection of Gaussian Hill

Results of the Gaussian hill problem are presented in
Figure2 and Table 1. The initial profile shown in dotted line was propagated
till t=0.6, for three Courant numbers of 0.5, 1.0, and
1.5. All the results
have been compared with results from Crank-Nicolson LSFEM as baseline. Results
of the space-time LSFEM are far more dissipative and dispersive compared to the
Crank-Nicolson LSFEM for all the three Courant numbers. However, results show
significant improvement with BELSFEM.Table 1
Convection of Gaussian hill byt=0.6.
CN-LSFEMST-LSFEMBE-LSFEMCourant no.UminUmaxUminUmaxUminUmax%redn.
inUmin%gain
inUmax0.5−0.00550.6861−0.01860.6784−
0.0013†0.6967†76.91.551.0−0.04900.6606−0.10040.6196≈00.7140≈1008.081.5−0.11960.6210−0.15360.5532−0.10490.640112.23.1†with one bubble in
both x and t.Propagation
of Gaussian hill by timet=0.6
for Courant numbers, C=0.5 (a),C=1.0
(b) and C=1.5
(c)
for continuous LSFEM.
(a)(b)(c)For Courant  number,C=0.5, BELSFEM with one bubble
in x and
t direction gives
1.5% increase in maximum
value and 77% decrease in dispersion error compared to Crank-Nicolson
LSFEM. More than one bubble in fact degraded the results.ForC=1.0, 8 bubbles in
x and 10 in time completely remove
the dispersion
error and increase
the peak by around 8% leading to complete recovery of the exact solution.ForC=1.5, BELSFEM with 8 and 10 bubbles in x and t, respectively, causes 12.2% reduction in dispersion error and about 3% increase
in the peak value.


### 6.1.2. Propagation of a Steep Front

Discontinuity was propagated byt=0.6
and the results  presented in Figure 3 and
Table 2 were computed for Courant
numbers of 0.75, 1.0, and 2.0. Few parameters were considered for comparative
quantification of the results. Slope, m, of the solution at the discontinuity
which indicates the amount of dissipation in the solution was measured across
the two nodes that capture the discontinuity in exact solution. Since the
discontinuity spanned one element (h=1/50), the exact solution had a slope, m=−50.
Also considered were the values of Umax
and Umin
causing the overshoot and undershoot representatives
of the dispersive error. All the comparative results were based on the results
from Crank-Nicolson LSFEM.Table 2
Propagation
of discontinuity byt=0.6.
CN-LSFEMST-LSFEMBE-LSFEMCourant no.slopemUminUmaxSlopemUminUmaxSlopemUminUmax%gain inm%redn. inUmin%redn.
inUmax0.75−12.66−0.00051.1341−9.78901.1740−14.64−0.1791.000115.6−356.611.81.0−10.33none1.1684−7.9650.00011.193−48.3101.0109367.7—13.52.0−5.947none1.2934−4.9070.00541.2232−5.6110.00251.245−5.65—3.8Propagation of a steep front by timet=0.6 for Courant numbers, C=0.75 (a),
C=1.0
(b), and C=2.0
(c) for continuous
LSFEM.
(a)(b)(c)Space-time LSFEM is more dissipative than CNLSFEM for
all the three Courant numbers as can be seen in Figure3. However, it is more
dispersive than Crank-Nicolson LSFEM for C=0.75 and 1.0 and less dispersive
for C=2.AtC=0.75, BELSFEM with 8 bubbles in
x and 10
bubbles in time causes 15.6% increase in the slope (meaning reduced dissipative
error) but a large increase in dispersive error in the form of a deep
undershoot. Although results are much better with one bubble each in x and t directions with 40% increase in the slope and much smaller undershoot, as can
be seen in Figure 3.AtC=1.0, the 8/10
bubble combination shows a
significant improvement in the results as slope m reaches very close to the
exact value of −50 (see Table 2) and the dispersion error completely disappears
and the solution looks almost like the exact solution (see Figure 3).AtC=2.0, BELSFEM fails to better the slope of
Crank-Nicolson LSFEM, although it is less dispersive.


## 6.1.1. Convection of Gaussian Hill

Results of the Gaussian hill problem are presented in
Figure2 and Table 1. The initial profile shown in dotted line was propagated
till t=0.6, for three Courant numbers of 0.5, 1.0, and
1.5. All the results
have been compared with results from Crank-Nicolson LSFEM as baseline. Results
of the space-time LSFEM are far more dissipative and dispersive compared to the
Crank-Nicolson LSFEM for all the three Courant numbers. However, results show
significant improvement with BELSFEM.Table 1
Convection of Gaussian hill byt=0.6.
CN-LSFEMST-LSFEMBE-LSFEMCourant no.UminUmaxUminUmaxUminUmax%redn.
inUmin%gain
inUmax0.5−0.00550.6861−0.01860.6784−
0.0013†0.6967†76.91.551.0−0.04900.6606−0.10040.6196≈00.7140≈1008.081.5−0.11960.6210−0.15360.5532−0.10490.640112.23.1†with one bubble in
both x and t.Propagation
of Gaussian hill by timet=0.6
for Courant numbers, C=0.5 (a),C=1.0
(b) and C=1.5
(c)
for continuous LSFEM.
(a)(b)(c)For Courant  number,C=0.5, BELSFEM with one bubble
in x and
t direction gives
1.5% increase in maximum
value and 77% decrease in dispersion error compared to Crank-Nicolson
LSFEM. More than one bubble in fact degraded the results.ForC=1.0, 8 bubbles in
x and 10 in time completely remove
the dispersion
error and increase
the peak by around 8% leading to complete recovery of the exact solution.ForC=1.5, BELSFEM with 8 and 10 bubbles in x and t, respectively, causes 12.2% reduction in dispersion error and about 3% increase
in the peak value.


## 6.1.2. Propagation of a Steep Front

Discontinuity was propagated byt=0.6
and the results  presented in Figure 3 and
Table 2 were computed for Courant
numbers of 0.75, 1.0, and 2.0. Few parameters were considered for comparative
quantification of the results. Slope, m, of the solution at the discontinuity
which indicates the amount of dissipation in the solution was measured across
the two nodes that capture the discontinuity in exact solution. Since the
discontinuity spanned one element (h=1/50), the exact solution had a slope, m=−50.
Also considered were the values of Umax
and Umin
causing the overshoot and undershoot representatives
of the dispersive error. All the comparative results were based on the results
from Crank-Nicolson LSFEM.Table 2
Propagation
of discontinuity byt=0.6.
CN-LSFEMST-LSFEMBE-LSFEMCourant no.slopemUminUmaxSlopemUminUmaxSlopemUminUmax%gain inm%redn. inUmin%redn.
inUmax0.75−12.66−0.00051.1341−9.78901.1740−14.64−0.1791.000115.6−356.611.81.0−10.33none1.1684−7.9650.00011.193−48.3101.0109367.7—13.52.0−5.947none1.2934−4.9070.00541.2232−5.6110.00251.245−5.65—3.8Propagation of a steep front by timet=0.6 for Courant numbers, C=0.75 (a),
C=1.0
(b), and C=2.0
(c) for continuous
LSFEM.
(a)(b)(c)Space-time LSFEM is more dissipative than CNLSFEM for
all the three Courant numbers as can be seen in Figure3. However, it is more
dispersive than Crank-Nicolson LSFEM for C=0.75 and 1.0 and less dispersive
for C=2.AtC=0.75, BELSFEM with 8 bubbles in
x and 10
bubbles in time causes 15.6% increase in the slope (meaning reduced dissipative
error) but a large increase in dispersive error in the form of a deep
undershoot. Although results are much better with one bubble each in x and t directions with 40% increase in the slope and much smaller undershoot, as can
be seen in Figure 3.AtC=1.0, the 8/10
bubble combination shows a
significant improvement in the results as slope m reaches very close to the
exact value of −50 (see Table 2) and the dispersion error completely disappears
and the solution looks almost like the exact solution (see Figure 3).AtC=2.0, BELSFEM fails to better the slope of
Crank-Nicolson LSFEM, although it is less dispersive.


## 6.2. Two-Dimensional Problems

### 6.2.1. Convection of a Concentration Spike

The concentration spike was convected linearly byt≈1.3 at a
unit velocity given by u=0.25 and v=0.1177 and making an angle of 25° with the
x-direction for Courant numbers of 0.73, 1.0, and
1.47. Results are presented in Figures 4, 5, and
Table 3.
Figure 4 presents the
variation of maximum and minimum concentrations with time and
Figure 5 shows
typical plot of concentration profile before and after being convected. For all
the Courant numbers, tested Space-time LSFEM is far more dissipative and
dispersive compared to Crank-Nicolson LSFEM (see
Figures 4, 5, and
Table 3).
However, there is a marked improvement in the results with bubbles. In
addition, the maximum number of PCG iterations per time step required to
achieve tolerance remained consistent as the number of bubble functions was
increased as shown in Table 3. This clearly indicates the ability of the
BELSFEM to increase accuracy without dramatically increasing computational
effort.Table 3
Advection
of concentration spike byt=1.3.
CN-LSFEMST-LSFEMBE-LSFEMImprovementsRange of
PCG iterations*Courant no.UminUmaxUminUmaxUminUmax%redn.
inUmin%gain
inUmax0.73−0.06950.6119−0.08650.5692−0.05310.694123.613.46-71.0−0.08430.5663−0.10790.5143−0.04860.678042.319.735-61.47−0.13910.5164−0.12500.4462−0.10840.570422.110.46-7Variation of maximum and minimum
concentrations with time for advection of concentration spike : comparison of results over the time of advection.
(a)(b)(c)(d)(e)(f)Convective transport of the concentration spike (initial condition shown by the left cone) with flow at25° to
x-axis for C=1.0.
(a)(b)(c)AtC=0.73 (the same C used by Yu and  Heinrich [21]
in convecting the same profile with Petrov-Galerkin formulation), BELSFEM with
6 bubbles each in spatial and time directions results in 23.6% increase in
Umax
and
13.4% decrease in Umin
compared to Crank-Nicolson LSFEM.Results further improved forC=1.0 as 42.3%
increase in Umax
and 20%
decrease in Umin
accrued (see Figures 4, 5, and
Table 3). And finally
for C=1.47, about
22% increase in Umax
and
10.4% decrease in Umin
were recorded.


### 6.2.2. Rotating Cosine Hill Problem

Results for rotating cosine hill problem are shown in
Figures6, 7
and Table 4. The variation of maximum and minimum values of
concentration over one rotation for ▵t=2π/120,2π/60, and 2π/30 is shown in
Figure 6. A typical profile after
one rotation is shown for the three formulations in Figure 7. Again,
Crank-Nicolson LSFEM serves as the baseline for comparison.Table 4
Advection
of cosine hill in rotation.
CN-LSFEMST-LSFEMBE-LSFEMImprovementsRange of
PCG iterations*ΔtCourant no.‡UminUmaxUminUmaxUminUmax%redn.
inUmin%gain
inUmax2π/1200.2618−0.02650.9691−0.04050.958−0.01890.976928.80.817-82π/600.5236−0.06150.9165−0.11380.8872−0.02700.971356.15.988-82π/301.047−0.20090.8369−0.21920.7398−0.21290.8418−5.970.615-16‡Courant number
based on velocity of the peak.
*Range of PCG iterations/time step: number of iterations with one bubble—the number with six bubbles.
% Reduction and %
gain calculated on Crank-Nicolson LSFEM results as baseline.Variation of maximum and minimum concentrations with time for advection of cosine hill in rotation : comparison of results over one rotation.
(a)(b)(c)(d)(e)(f)Convection of a cosine hill in a pure rotational velocity
field with▵t=2π/60 : comparison of results after a complete revolution.
(a)(b)(c)For▵t=2π/120, BELSFEM with 6 bubbles each in spatial
and time directions shows about 29% reduction in dispersive error and about 1%
increase in the peak value. This improvement in the peak value is significant
considering the fact that the baseline value from Crank-Nicolson LSFEM itself
was high at 0.9691 (see Table 4).For▵t=2π/60, there is more improvement in the results
as the dispersion error declines by 56% and the peak value goes up by around
6%. Typical profiles after one rotation for this case are shown in Figure 7.For▵t=2π/30 (which corresponds to C≈1, based on velocity at the peak of the
profile), however, there is only 0.6% improvement in peak value and the
dispersive error is worse than CNLSFEM, as can be seen in Figure 6.


## 6.2.1. Convection of a Concentration Spike

The concentration spike was convected linearly byt≈1.3 at a
unit velocity given by u=0.25 and v=0.1177 and making an angle of 25° with the
x-direction for Courant numbers of 0.73, 1.0, and
1.47. Results are presented in Figures 4, 5, and
Table 3.
Figure 4 presents the
variation of maximum and minimum concentrations with time and
Figure 5 shows
typical plot of concentration profile before and after being convected. For all
the Courant numbers, tested Space-time LSFEM is far more dissipative and
dispersive compared to Crank-Nicolson LSFEM (see
Figures 4, 5, and
Table 3).
However, there is a marked improvement in the results with bubbles. In
addition, the maximum number of PCG iterations per time step required to
achieve tolerance remained consistent as the number of bubble functions was
increased as shown in Table 3. This clearly indicates the ability of the
BELSFEM to increase accuracy without dramatically increasing computational
effort.Table 3
Advection
of concentration spike byt=1.3.
CN-LSFEMST-LSFEMBE-LSFEMImprovementsRange of
PCG iterations*Courant no.UminUmaxUminUmaxUminUmax%redn.
inUmin%gain
inUmax0.73−0.06950.6119−0.08650.5692−0.05310.694123.613.46-71.0−0.08430.5663−0.10790.5143−0.04860.678042.319.735-61.47−0.13910.5164−0.12500.4462−0.10840.570422.110.46-7Variation of maximum and minimum
concentrations with time for advection of concentration spike : comparison of results over the time of advection.
(a)(b)(c)(d)(e)(f)Convective transport of the concentration spike (initial condition shown by the left cone) with flow at25° to
x-axis for C=1.0.
(a)(b)(c)AtC=0.73 (the same C used by Yu and  Heinrich [21]
in convecting the same profile with Petrov-Galerkin formulation), BELSFEM with
6 bubbles each in spatial and time directions results in 23.6% increase in
Umax
and
13.4% decrease in Umin
compared to Crank-Nicolson LSFEM.Results further improved forC=1.0 as 42.3%
increase in Umax
and 20%
decrease in Umin
accrued (see Figures 4, 5, and
Table 3). And finally
for C=1.47, about
22% increase in Umax
and
10.4% decrease in Umin
were recorded.


## 6.2.2. Rotating Cosine Hill Problem

Results for rotating cosine hill problem are shown in
Figures6, 7
and Table 4. The variation of maximum and minimum values of
concentration over one rotation for ▵t=2π/120,2π/60, and 2π/30 is shown in
Figure 6. A typical profile after
one rotation is shown for the three formulations in Figure 7. Again,
Crank-Nicolson LSFEM serves as the baseline for comparison.Table 4
Advection
of cosine hill in rotation.
CN-LSFEMST-LSFEMBE-LSFEMImprovementsRange of
PCG iterations*ΔtCourant no.‡UminUmaxUminUmaxUminUmax%redn.
inUmin%gain
inUmax2π/1200.2618−0.02650.9691−0.04050.958−0.01890.976928.80.817-82π/600.5236−0.06150.9165−0.11380.8872−0.02700.971356.15.988-82π/301.047−0.20090.8369−0.21920.7398−0.21290.8418−5.970.615-16‡Courant number
based on velocity of the peak.
*Range of PCG iterations/time step: number of iterations with one bubble—the number with six bubbles.
% Reduction and %
gain calculated on Crank-Nicolson LSFEM results as baseline.Variation of maximum and minimum concentrations with time for advection of cosine hill in rotation : comparison of results over one rotation.
(a)(b)(c)(d)(e)(f)Convection of a cosine hill in a pure rotational velocity
field with▵t=2π/60 : comparison of results after a complete revolution.
(a)(b)(c)For▵t=2π/120, BELSFEM with 6 bubbles each in spatial
and time directions shows about 29% reduction in dispersive error and about 1%
increase in the peak value. This improvement in the peak value is significant
considering the fact that the baseline value from Crank-Nicolson LSFEM itself
was high at 0.9691 (see Table 4).For▵t=2π/60, there is more improvement in the results
as the dispersion error declines by 56% and the peak value goes up by around
6%. Typical profiles after one rotation for this case are shown in Figure 7.For▵t=2π/30 (which corresponds to C≈1, based on velocity at the peak of the
profile), however, there is only 0.6% improvement in peak value and the
dispersive error is worse than CNLSFEM, as can be seen in Figure 6.


## 6.3. Effect of Mesh Size and Number of Bubbles

Two-dimensional benchmarks problems were run on
different sizes of mesh and also on basic meshes with different number of
bubble functions in order to investigate the effect of mesh size and number of
bubble functions on the performance of BELSFEM. Mesh size parameter,h,
was varied from 0.01 to 0.1 (where h = side-length/number of elements per side). In all the cases,
h in x and y directions was the same.Typical  comparative plots ofUmax
and Umin
from  the three least-square methods for different mesh sizes are shown in
Figure 8. For this part of study, four bubbles each in space and time were used. The
maximum and minimum values for the cosine hill are recorded after one full
rotation and those for linear convection of concentration spike have been taken
after being convected by t=1.3. The bubbles seem to be most effective for
moderately coarse meshes as can be observed from the figure where large gain over
both CNLSFEM and STLSFEM can be seen in this region. However, for very coarse
and very fine meshes the benefits of bubbles seem to diminish.Effect of mesh sizeh on the performance of BELFEM and LSFEM.
(a)(b)Figure9 shows the effect of number of bubbles on the
performance of BELSFEM. Typical
variation of Umax
and Umin
for the two problems
is displayed. Results improve sharply with the number of bubble functions
initially but the improvements diminish with further increase in the number and
beyond 3-4 bubbles the effect saturates. It, therefore, can be stated that
generally good improvements in the results can be achieved with 4–6 bubbles.Effect of number of bubbles on the performance of BELFEM. (*pure LS method-results (not functions of n)).
(a)(b)These results show the clear benefit of bubble
functions for linear transport problems, which are purely hyperbolic in nature.
Extensions of this work to mixed problems, such as Navier-Stokes equations, are
of great practical interest and a topic of further research. In addition, there
likely exist optimal bubble functions that will achieve highly accurately
solutions with a small number of functions. The form of these functions is also
a topic of further research.


## 7. Conclusions

A study of Crank-Nicolson least square finite element
method, space-time least square finite element method, was done and the effect
of the bubble modes applied to linear space-time elements was investigated.
Orthogonal Jocobi polynomials were chosen as the bubble functions. Convection
of a Gaussian hill and propagation of a discontinuity in one-dimension and
linear convection of a concentration spike and convection of a cosine hill in
rotation inx-y plane were the standard test problems considered.Emphasis of the current study was to prove the
effectiveness of bubble modes towards generating improved solution for the
linear convection equation without resorting to expensive higher order elements
and severe mesh refinement which undermines the utility of a scheme. Additional
computational work was required on element level due to introduction of bubble
modes and keeping more or less same amount of computation on global level
overall. This was to great extent achieved due to the fact that bubble modes are
easily condensed out using the classic static condensation procedure.It was observed that bubbles greatly improve the
accuracy of the least-squares method compared to the otherwise dissipative and
dispersive space-time least square finite element formulation. The results thus
achieved were compared with the results from Crank-Nicolson least square
formulation. It was observed that the addition of bubble modes increasingly
improves the performance of STLSFEM till about 8 bubble modes when the effect
seems to saturate. It was recorded that for convection of Gaussian hill the
peak value of the profile improves in the range of 1.5%–8% for the CFL numbers of 0.5, 1.0, and 1.5.
Decline of the order of 12%–100% in the dispersion error was seen. In case ofC=1.0, the dissipation and dispersion errors were almost completely removed.
Similar trends were observed in the problem of propagation of discontinuity,
where considerable steepening of profile was observed along with decrease in
the dispersive error almost for all the cases. Here too exact solution was
almost completely recovered for C=1.More interesting results were obtained in two
dimensional test cases. In case of linear convection of concentration spike, an
increase in peak profile value in the range of 10%–20% and a decrease in
dispersive error in the range 22%–43% were recorded for the three Courant
numbers tested. In the second test problem of rotation of cosine hill, also an increase
in peak value of the order of 1%–6% and a decrease in dispersion error in the
range 20%–56% were recorded although in case of▵t=2π/60; a 5% increase in dispersive error
occurred.Overall, the bubble enriched least-squares finite
element method (BELSFEM) seems to be very promising though further work is
required to determine the optimal form of the bubble functions.



---

*Source: 267454-2008-11-26.xml*