Differential Equations and Nonlinear Mechanics
Volume 2008 (2008), Article ID 267454, 21 pages
doi:10.1155/2008/267454

Copyright © 2008 Rajeev Kumar and Brian H. Dennis. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

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 the -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 -version of STLSFEM, where they have used -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 -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 -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)where is the property being convected at a velocity with , and as its components in , and 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)In the least-squares approach, the -norm of the differential equation is minimized with respect to unknown coefficients over the solution domain Ω. Applying the -norm to (2.2) and minimizing the functional with respect to leads to the weak statement (2.3) where the row vector contains the basis functions used to approximate the solution over the domain as .

The weak statement can be expanded and written in matrix form (2.4)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.

3. The Least-Square Finite Element Formulations

For the sake of simplicity, let us consider 1D scalar advection equation (3.1)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, , given by , where is the approximate solution. For sake of simplicity, we will use in place of . The LSFEM formulation based on minimization of square of residual leads to (3.2)Using forward difference for time derivative term and θ-method for approximating U in convective term gives (3.3)Let the unknown be defined as (3.4)where is the solution at the jth node and is the interpolation function. Taking the derivative with respect to , (3.3) leads to the Crank-Nicolson LSFE formulation (3.5) For , it becomes Crank-Nicolson LSFEM formulation as (3.6)

3.2. 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, (3.7)where is bilinear interpolation function for 1D and 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)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)Similarly, Trilinear shape functions for trilinear elements are given by (3.9b) where and are the linear shape functions and and 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 by . 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 with respect to the function . Bubble modes were generated from as (3.10)where p is the order of the Jacobi polynomial. Jacobi polynomials with 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.

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)where , and the boundary condition as and convection velocity . The solution was convected by over a uniform mesh of size . The exact solution is given by (4.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 of the domain ]0,1[.

The discontinuity is given as (4.3) The solution was convected by using a mesh of uniform size .

4.2. Two-Dimensional Problems

4.2.1. Convection of a Concentration Spike

A concentration spike, given by (4.4) was convected by with a velocity given by and at an angle of to the -axis. A mesh in 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) where and , and the boundary condition is on . The initial positions of the center and the radius of the cosine hill are and , respectively. The angular velocity is given by . A uniform mesh of four-node elements over the unit square was used in the computations.

5. Calculation of Flow Parameters

Important flow parameter, Courant number, is given as , where u is the convection velocity, is the time step, and is the characteristic length in the direction of the convection. In one-dimensional problems, h is simply taken as and . In the first problem of convection of Gaussian hill, and in the second problem of propagation of discontinuity 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)where is the velocity vector and Courant number was given as (5.2)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 , 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 of 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 Figure 2 and Table 1. The initial profile shown in dotted line was propagated till , 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, , BELSFEM with one bubble in and 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 , 8 bubbles in 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 , BELSFEM with 8 and 10 bubbles in and , 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 by 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, , 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 (), the exact solution had a slope, . Also considered were the values of and 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 and 1.0 and less dispersive for .

At , 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 , the 8/10 bubble combination shows a significant improvement in the results as slope reaches very close to the exact value of (see Table 2) and the dispersion error completely disappears and the solution looks almost like the exact solution (see Figure 3).

At , 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 by at a unit velocity given by and and making an angle of 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 (the same 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 and 13.4% decrease in compared to Crank-Nicolson LSFEM.

Results further improved for as 42.3% increase in and 20% decrease in accrued (see Figures 4, 5, and Table 3). And finally for , about 22% increase in and 10.4% decrease in were recorded.

6.2.2. 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 , and 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 , 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 , 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 (which corresponds to , 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 of and 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 . 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 and 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.

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 in - 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 , 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 .

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

Acknowledgment

The authors would like to acknowledge the partial support for this research by the Texas Space Grant Consortium through New Investigations Grant UTA-06-685.

References

1. B.-N. Jiang, The Least-Squares Finite Element Method: Theory and Applications in Computational Fluid Dynamics and Electromagnetics, Scientific Computation, Springer, Berlin, Germany, 1998.
2. I. Christie, D. F. Griffiths, A. R. Mitchell, and O. C. Zienkiewicz, “Finite element methods for second order differential equations with significant first derivatives,” International Journal for Numerical Methods in Engineering, vol. 10, no. 6, pp. 1389–1396, 1976.
3. J. J. Westerink and D. Shea, “Consistent higher degree Petrov-Galerkin methods for the solution of the transient convection-diffusion equation,” International Journal for Numerical Methods in Engineering, vol. 28, no. 5, pp. 1077–1101, 1989.
4. A. N. Brooks and T. J. R. Hughes, “Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations,” Computer Methods in Applied Mechanics and Engineering, vol. 32, no. 1–3, pp. 199–259, 1982.
5. J. Donea, “A Taylor-Galerkin method for convective transport problems,” International Journal for Numerical Methods in Engineering, vol. 20, no. 1, pp. 101–119, 1984.
6. B.-N. Jiang and L. A. Povinelli, “Least-squares finite element method for fluid dynamics,” Computer Methods in Applied Mechanics and Engineering, vol. 81, no. 1, pp. 13–37, 1990.
7. B.-N. Jiang, T. L. Lin, and L. A. Povinelli, “Large-scale computation of incompressible viscous flow by least-squares finite element method,” Computer Methods in Applied Mechanics and Engineering, vol. 114, no. 3-4, pp. 213–231, 1994.
8. J. Donea and L. Quartapelle, “An introduction to finite element methods for transient advection problems,” Computer Methods in Applied Mechanics and Engineering, vol. 95, no. 2, pp. 169–203, 1992.
9. G. F. Carey and B.-N. Jiang, “Least-squares finite elements for first-order hyperbolic systems,” International Journal for Numerical Methods in Engineering, vol. 26, no. 1, pp. 81–93, 1988.
10. C. W. Li, “Least-squares characteristics and finite elements for advection-dispersion simulation,” International Journal for Numerical Methods in Engineering, vol. 29, no. 6, pp. 1343–1358, 1990.
11. N.-S. Park and J. A. Liggett, “Taylor-least-squares finite element for two-dimensional advection-dominated unsteady advection-diffusion problems,” International Journal for Numerical Methods in Fluids, vol. 11, no. 1, pp. 21–38, 1990.
12. N.-S. Park and J. A. Liggett, “Application of Taylor-least squares finite element to three-dimensional advection-diffusion equation,” International Journal for Numerical Methods in Fluids, vol. 13, no. 6, pp. 759–773, 1991.
13. H. Nguyen and J. Reynen, “A space-time least-square finite element scheme for advection-diffusion equations,” Computer Methods in Applied Mechanics and Engineering, vol. 42, no. 3, pp. 331–342, 1984.
14. T. J. R. Hughes, L. P. Franca, and G. M. Hulbert, “A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations,” Computer Methods in Applied Mechanics and Engineering, vol. 73, no. 2, pp. 173–189, 1989.
15. K. S. Surana and J. S. Sandhu, “Investigation of diffusion in p-version ‘LSFE’ and ‘STLSFE’ formulations,” Computational Mechanics, vol. 16, no. 3, pp. 151–169, 1995.
16. C. Baiocchi, F. Brezzi, and L. P. Franca, “Virtual bubbles and Galerkin-least-squares type methods (Ga.L.S.),” Computer Methods in Applied Mechanics and Engineering, vol. 105, no. 1, pp. 125–141, 1993.
17. F. Brezzi, L. P. Franca, and A. Russo, “Further considerations on residual-free bubbles for advective-diffusive equations,” Computer Methods in Applied Mechanics and Engineering, vol. 166, no. 1-2, pp. 25–33, 1998.
18. L. P. Franca, A. Nesliturk, and M. Stynes, “On the stability of residual-free bubbles for convection-diffusion problems and their approximation by a two-level finite element method,” Computer Methods in Applied Mechanics and Engineering, vol. 166, no. 1-2, pp. 35–49, 1998.
19. G. E. Karniadakis and S. J. Sherwin, Spectral/hp Element Methods for CFD, Numerical Mathematics and Scientific Computation, Oxford University Press, New York, NY, USA, 1999.
20. J. Donea and A. Huerta, Finite Element Methods for Flow Problems, John Wiley & Sons, Chichester, UK, 2003.
21. C.-C. Yu and J. C. Heinrich, “Petrov-Galerkin method for multidimensional, time-dependent, convective-diffusion equations,” International Journal for Numerical Methods in Engineering, vol. 24, no. 11, pp. 2201–2215, 1987.