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

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 the L2-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.


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 firstorder 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 related Galerkin finite element method GFEM , both leading to a nonsymmetric, nonpositive definite matrices as Jiang has illustrated in his text 1 .
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 and p-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.

The least-square finite element method
Consider the transient advection equation given as 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: In the least-squares approach, the L 2 -norm of the differential equation is minimized with respect to unknown coefficients over the solution domain Ω.Applying the L 2 -norm to 2.2 and minimizing the functional with respect to U n 1 leads to the weak statement where the row vector {N} contains the basis functions N j used to approximate the solution over the domain as U x, y, z j N j x U j {N} T {U}.

D i fferential Equations and Nonlinear Mechanics
The weak statement can be expanded and written in matrix form where the individual matrix contributions are given by

2.5
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.

The least-square finite element formulations
For the sake of simplicity, let us consider 1D scalar advection equation ∂U ∂t a ∂U ∂x 0.

3.1
The three least-square finite element formulations tried are as follows.

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 Using forward difference for time derivative term and θ-method for approximating U in convective term gives Let the unknown U be defined as where U j is the solution at the jth node and N j is the interpolation function.Taking the derivative with respect to U n 1 , 3.3 leads to the Crank-Nicolson LSFE formulation

Space-time LSFEM
In space-time formulation, both time and space derivatives are discretized the finite element way and the unknown U becomes function of both spatial and temporal variables, that is,

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 by P α,β p .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 P α,β p as 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

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

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 where x 0 2/15, l 7 √ 2/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

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 The solution was convected by t 0.6 using a mesh of uniform size h 1/50.

Convection of a concentration spike
A concentration spike, given by U x, y, 0 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.

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 where X x − x 0 /σ and Y y − y 0 /σ, and the boundary condition is U 0 on Γ in .The initial positions of the center and the radius of the cosine hill are x 0 , y 0 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.

Calculation of flow parameters
Important flow parameter, Courant number, is given as C 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 where u ui vj is the velocity vector and Courant number was given as 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 ωr peak , where ω is the angular velocity.

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 of 1.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.

Convection of Gaussian hill
Results of the Gaussian hill problem are presented in Figure 2 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.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.
For C 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.For C 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.

Propagation of a steep front
Discontinuity was propagated by t 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 U max and U min causing the overshoot and undershoot representatives of the dispersive error.All the comparative results were based on the results from Crank-Nicolson LSFEM.Space-time LSFEM is more dissipative than CNLSFEM for all the three Courant numbers as can be seen in Figure 3.However, it is more dispersive than Crank-Nicolson LSFEM for C 0.75 and 1.0 and less dispersive for C 2.
At C 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.
At C 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 .
At C 2.0, BELSFEM fails to better the slope of Crank-Nicolson LSFEM, although it is less dispersive.

Convection of a concentration spike
The concentration spike was convected linearly by t ≈ 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.
At C 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 U max and 13.4% decrease in U min compared to Crank-Nicolson LSFEM.
Results further improved for C 1.0 as 42.3% increase in U max and 20% decrease in U min accrued see Figures 4, 5, and Table 3 .And finally for C 1.47, about 22% increase in U max and 10.4% decrease in U min were recorded.

Rotating cosine hill problem
Results for rotating cosine hill problem are shown in Figures 6, 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.
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.

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 of U max and U min 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 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.Figure 9 shows the effect of number of bubbles on the performance of BELSFEM.Typical variation of U max and U min 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.
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.

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 spacetime elements was investigated.Orthogonal Jocobi polynomials were chosen as the 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 in x-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 of C 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.

Figure 2 :
Figure 2: Propagation of Gaussian hill by time t 0.6 for Courant numbers, C 0.5 a , C 1.0 b and C 1.5 c for continuous LSFEM.

Figure 3 :
Figure 3: Propagation of a steep front by time t 0.6 for Courant numbers, C 0.75 a , C 1.0 b , and C 2.0 c for continuous LSFEM.

Figure 4 :Figure 5 :
Figure 4: Variation of maximum and minimum concentrations with time for advection of concentration spike : comparison of results over the time of advection.

Figure 6 :Figure 7 :
Figure 6: Variation of maximum and minimum concentrations with time for advection of cosine hill in rotation : comparison of results over one rotation.

1 :
Convection of Gaussian hill by t 0.6.

2 :
Propagation of discontinuity by t 0.6.

Table 3 :
Advection of concentration spike by t 1.3.

Table 4 :
Advection of cosine hill in rotation.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.
‡Courant number based on velocity of the peak.*