# One-Dimensional Fluid – Structure – Interaction (1D-FSI) Steady Flow Stable Formulation

DOI : 10.17577/IJERTV6IS040777

Text Only Version

#### One-Dimensional Fluid – Structure – Interaction (1D-FSI) Steady Flow Stable Formulation

A. Y. Tang1 , N.Amin4

Department of Mathematical Sciences, Faculty of Science,

Universiti Teknologi Malaysia, 81310 Johor Bahru, Johor

1. A. Mohd Noor2 Faculty of Civil Engineering, Universiti Teknologi Malaysia,

81310 Johor Bahru, Johor

1. Y. Mohd Yassin3

Department of Civil and Environmental Engineering, Universiti Teknologi PETRONAS,

32610 Bandar Seri Iskandar, Perak.

Abstract This paper concerns the numerical solution of one-dimensional fluid-structure-interaction (1D-FSI) formulation which has been formulated by providing a pressure-area constitutive relation to complement the mass and linear momentum equations. However, typical spurious oscillations were found for the cases of relatively high pressure difference when Bubnov-Galerkin formulation was employed. In minimizing the oscillation, SUPG stabilization scheme was then formulated and shown as able to stabilize the solutions. For validation purposes, an analytical solution for the limited case of straight vessel has been derived for a specific pressure-area constitutive relation. This study can be important for future works in 1D-FSI employing pressure-area constitutive relation.

Keywords Streamline-Upwind-Petrov-Galerkin, Finite Element Method, Biomechanics

1. INTRODUCTION

Fluid-structure-interaction (FSI) for one-dimensional flow can be formulated by providing a pressure-area constitutive relation to complement the mass and linear momentum equations. Such coupling would allow the interaction between volumetric flow rate, Q, cross-sectional area, A, and pressure, p, of the flow. The constitutive relation can be given in general form as;

0 = () (1)

reported. However, when we repeated the formulation and applied it to high pressure differences that is, in the range higher than reported, spurious oscillations were observed. These oscillations are typical phenomenon of Bubnov- Galerkin formulation hence its shortcoming. In minimizing the oscillation, we then formulated the well-known stabilization scheme, Streamline-Upwind-Petrov-Galerkin (SUPG) for the problem.

Realizing the importance of having stabilized solution in ensuring the attainment of reliable information, it is the interest of this paper, therefore, to report such a formulation for future reference especially in the study of one-dimensional fluid-structure-interaction (1D-FSI) flow employing pressure- area constitutive relation.

2. GOVERNING EQUATIONS

1D-FSI steady flow is governed by the conservation laws of mass and linear momentum as follows [10, 12]:

Mass equation

= 0 (2)

Momentum equation

2

where and 0 are the local and reference pressure respectively, and () highlights the dependency of the

(

) + + = 0 (3)

pressures magnitude and distribution on the cross-sectional area of the flow. Various detailed forms of () have been proposed in the literature [1-13] To note, pressure-area constitutive relation as in Eqn. (1) is also termed as tube law elsewhere [5-6]. Employment of Eqn. (1) thus the formulation has a wide range of applications especially in the field of biomechanics.

Due to the complexity of the governing equations,

where is the fluid density, is the momentum correction

factor and is the axial coordinate. is defined as viscosity friction coefficient as follows

= 2

( 1) (4)

where is the fluid viscosity. With regards to the second term in Eqn. (3), we have

= =

solutions are mostly obtained numerically. In recent works of Sochi [10, 12], Bubnov-Galerkin finite element method has been formulated where good verifications of results were

=

(5)

2.1 Constitutive Relation Equations

Since there are three dependent variables, A, Q and p, a

= ( [1

1ln( )

third equation, that is, the pressure-area constitutive relation, must be provided. Despite the various constitutive relations

1 + 1

(11)

available from literatures, only two relationships are considered in this study since the main purpose is to demonstrate the formulation of a specific numerical technique. Herein, the first constitutive equation is termed as Model 1 and is given by Eqn. (6). Despite its simplistic nature, the equation is chosen because it is the relationship used in Sochi [12] which results we are comparing against. For completeness, a more realistic constitutive equation (as it involves experimental-fit parameters) is thus considered and termed herein as Model 2 (Eqn. (9)).

Model 1

Pressure-area constitutive relation used in Sherwin et al. [8], Quarteroni and Formaggia [9], and Sochi [10-13] is termed as

Model 1 herein and given as

2ln()

2 ])

2 1

3. SPURIOUS OSCILLATIONS

In numerical analysis of fluid dynamics, spurious oscillations can occur in flows with high Peclet number (for advection-diffusion problems) and high Reynolds number (for general flows) when solved using either central finite difference method or Bubnov-Galerkin finite element method. Mathematical wise, both are known to be closely related thus inherit the same numerical difficulty [14, 15, 16].

As mentioned, when we repeated the Bubnov-Galerkin

() =

0

( 0 ) (6)

formulation detailed in Sochi [12], whilst we obtained similar non-oscillatory results for the reported range of pressure

where is known as vessel stiffness, given as

differences (< 1000Pa), we started to observe spurious oscillation for higher pressure differences. These observations are depicted in Figure 1.

=

0 1 2

(7)

and 0 and 0 are the cross sectional area of the flow and vessels wall thickness at reference pressure 0, respectively whilst and are Youngs elastic modulus and Poissons ratio of the vessels wall, respectively. Accordingly, Eq. (5) can be written in expanded form as

3

=

( ) (8)

2

3

Model 2

Pressure-area constitutive relation used in Ku et al. [3] and Downing and Ku [4] is termed as Model 2 herein and given as

Fig. 1 Oscillation due to the employment of Bubnov formulation for Model 1.The tube, fluid and flow

() = ((

1

)

(

2

)

) (9)

parameters used are; = 1060 kgm3, = 0.0035 Pa s, = 1.3333, = 1 m, = 0.1 m and = 5 Ã— 104 Pa m

(same data were used in Sochi [12])

where 1 and 2 are parameters obtained from a fit to experimental data of pressure versus diameter curve for a bovine carotid artery as detailed in Downing and Ku [4]. The vessel stiffness, is defined as

3

In confirming the occurrence of the spurious oscillations, we then employed the same Bubnov-Galerkin formulation to a different pressure-area constitutive relation

i.e. ( Model 2) only to observe similar phenomenon, as depicted in Figure 2.

=

0

12(1 2)3

(10)

where is the mean flow radius. One of the set of values of

1 and 2 proposed in Downing and Ku [4] is used in this study, which are 7 ad 2.5, respectively. Accordingly, for

Model 2, Eqn. (5) can be written in expanded form as

which can be further simplified as

2

1

2

( 2 + (( ) 1 + () 2))

(14)

+ = 0

With some algebraic manipulations, Eq. (14) becomes

2 +

1

2 )

=

(()

1 + () 2

(15)

By integrating Eqn. (15) with respect to , we obtain

ln

2 2 ln

2 1 ln

=

0 + 2 0 1

0

Fig. 2 Oscillation due to the employment of Bubnov formulation for

Model 2. The tube, fluid and flow parameters used are; = 995 kgm3, =

(2 + 2)

+

(2 + 1)

(16)

0.003Pa s, = 1.3333, = 0.1 m, = 0.003 m and = 125 Pa (same data were used in Downing and Ku [4])

Having confirmed the occurrence of the spurious oscillations in this 1D-FSI flow as the typical phenomenon of

where is the constant of integration which can be determined from the boundary condition where = at

= 0, thus

Bubnov-Galerkin formulation and also realized the

ln

2 2 ln

importance for a stabilized solution, we then formulated SUPG stabilization scheme for this particular problem which derivation and results are reported herein.

=

0 2 0

(2 + 2)

2 1 ln

(17)

+ 1 0

4. ANALYTICAL SOLUTIONS

In this work, we verify our SUPG formulation against analytical solutions which are available for limited case of

straight vessel. For Model 1, the analytical solution for

(2 + 1)

Substituting Eqn. (17) into Eqn. (16), we obtain

the volumetric flow rate, has been derived in Sochi [10, 12]

ln

2 2 ln

which is given herein as

=

0 + 2

0

(2 + 2)

5 5

12

1 ln

0

ln

0

+ 22 4 ln( ) 5

(2 2 )

(12)

(2 + 1)

= 2 ln( )

0

2 2 ln

2 0

(2 + 2)

(18)

where is the length of vessel whilst

and

are the

12

+

1 ln

0

flow cross sectional area at the inlet and outlet respectively.

4.1 Derivation of Analytical Solution for Model 2

For Model 2, we have derived the analytical solution which derivation is detailed as follows. By inserting Eq. (11) into Eq. (3) and the fact that Q is spatially constant (refer Eq. (2)), the momentum equation can be restated as

(2 + 1)

Now, to obtain a closed form solution for , we then employ the other boundary condition that is at the outlet which can be given as, = at = . This result in

2

(

+ 1

1 ln

0 2

2 ln

0 )

(1 + 1)

+ = 0

(1 + 2)

(13)

= ( ln + ln )

Eqn. (19) can be rearranged to yield a quadratic polynomial in

0

0

2 2 ln

Q, given as

2

1

+ (

2

(2 + 2)

0

( ln

0

+ ln

0

) +

2

2

ln

0

22 2 ln 0

( (2+2)

2 2 ln 0

+ 2

(2+2)

(20)

+ 2

(19)

(2 + 2)

ln

2 1 0

12 1 ln 0

2 1 ln 0

1 + ) = 0

1

(2+1)

(2+1)

(2 + 1)

12

+

1 ln

0

)

(2 + 1)

We can solve Eqn. (20) for the roots of by applying the quadratic formula, thus obtain

Â± 2 4 ( ln + ln ) ( 2 [2 (( )2 2 ( )2 ) ] + 1 [2 (( 1 2 )1 ) ])

0

0

(2 + 2)

0

0

(2 + 1)

0 )

( 0

=

2 ( ln + ln )

0

0

(21)

If we limit the solution to a specific condition of > , the two roots must be real. Also, in ensuring the flow to be consistent in direction with the pressure gradient, the root with the plus sign should be chosen so that positive flow rate is obtained. This is because, since > it can be shown that the denominator is always positive and the square root is always greater than . So, the flow rate can be given as

+ 2 4 ( ln + ln ) ( 2 [2 (( )2 2 ( )2 ) ] + 1 [2 (( 1 2 )1 ) ] )

0

0

(2 + 2)

0

0

(2 + 1)

0 )

( 0

=

2 ( ln + ln )

0

0

(22)

Eqn. (22) is the equivalent of the Poiseuille equation for rigid tubes but instead of pressure difference, it is expressed in terms of the specified inlet and outlet areas i.e. , which actually represent the specification of pressure at boundaries through the constitutive equation given by Eqn.

(9) (e.g. using Eqn. (9), or is solved for the desired pressure at the boundary) .

compatibility conditions. The non-reflecting boundary conditions are used to project the differential equations in the direction of outgoing characteristic variables at the inlet and outlet in producing the compatibility conditions as proposed in Sochi [10, 12] and Thompson [17] and as detailed next. Based on the method of characteristics and by assuming that > 0, the eigenvalues and left eigenvectors of a matrix defined as

5. FINITE ELEMENT FORMULATION

The weak form of the formulation (after conducting

=

(26)

integration by parts to Eq. (2) and (3)) can be written as

can be obtained as follows. The eigenvalues can be obtained

[ ] =

(23)

by solving

x

( ( x ) + ) x +

|=0 = 0

det( ) = 0 (27)

where is the linear weighting functions whilst and give the vector representation of Eq. (2) and (3), given herein as

where is the eigenvalues, I is identity matrix and is defined in Eq. (26) as the matrix of partial derivative of with respect to . The matrix has two eigenvalues represented as 1,2 which can be obtained by the quadratic

formula. Left eigenvectors are then obtained by solving the

F = [2 +

] (24)

following system

0

B = [ ] (25)

The weak formulation of Eq. (37) is then coupled with the suitable boundary conditions through the introduction of

= (28)

where is the left eigenvectors of and is given as

= [1 0 ] (29)

0 2

Once the left eigenvectors, is obtained, the desired compatibility conditions is then given by

, =

[ Â± 2 (2 ) + [( )2

1 1]

1,2

L (H U + B) = 0 (30)

where , are the left-eigenvectors. The imposition of

2

0

2 + (0)

1]

(35)

boundary condition is accomplished by replacing the continuity equation at the boundary node with Eq. (30). The expanded expressions of the eigenvalues, left eigenvectors and compatibility conditions for each of the constitutive model are given next.

Inserting Eq. (35) into Eq. (30), we then obtain the compatibility conditions for Model which can be as

(

Model 1

2

2

1

The eigenvalues of for Model 1, obtained by

Â± 2 (2 ) +

[(0)

2 + ( )

0

1] )

solving Eq. (27) can be given in expanded form as

2

2

1

= Â± 2 (2 ) + ( ) (31)

+ ( 2 +

[( )

0

2 + ( )

0

1])

1,2

2

20

+ (2 + ) = 0

Inserting Eq. (31) into Eq. (29) and by solving Eq. (28), the left eigenvalues of for Model 1, can be given as

,

5.1 SUPG Formulation

(36)

2

(32)

In employing the SUPG stabilization technique, the

0

= [ Â± 2 (2 ) + (2 ) 1]

Inserting Eq. (32) into Eq. (30), we then obtain the compatibility conditions for Model 1 which can be given as

stabilization term is added to Eq. (23) to give

x

( ( x ) + ) x

+ (() ()) x

[ ] = = 0

+ |=0

(37)

( Â± 2 (2 ) + ( ) )

where () is the operator applied to the test function,

2

20

2

(33)

() is the residual of the governing equation, and is the stabilization parameter, all as detailed in Donea and Huerta

[18] and SoulaÃ¯mani and Fortin [19] and given herein as

+ ( 2 + (20))

w

+ (2 + ) = 0

P(w) = H (38)

R(U) = F + B (39)

x

Model 2

1

= () 2

(40)

The eigenvalues of for Model 2, obtained by solving Eq. (27) can be given in expanded form as

1,2

=

= H (41)

where = () is the actual coordinates whilst refers to normalized local coordinates. To note, the weighting

(34)

function in Eq. (37) is the same as the shape functions,

2

2

1

. The expanded expression of Eq. (37) for each

0

Â± 2 (2 ) + [( )

2 + ( )

0

1]

constitutive model is detailed next for the discretization that uses linear shape functions (i.e. where = 1,2).

Inserting Eq. (34) into Eq. (29) and by solving Eq. (28), the left eigenvalues of for Model 2, can be given as

Model 1

With the descriptions given by Eq. (38) to (41), expanded expression of Eq. (37) can be given for Model 1 as

0

0 0

x { } + [

] x { }

x

x

0

[

3

]

2

0

(

+ ) 0

+

2

2

[]

2

{}

(2

)

(

2 + 2

) (2

)

[

] [

]

2

0

(

+ ) 0 0

2

2

+

[] [0 ] {}

(2

)

[

]

0

=

0

+ ([

])|

{} = { }

where

3

|

0

=0

(42)

1

0 1 0 1 2

= ([

2

] [

2

])

Model 2

2 + 2 2

2 + 2 2

(43)

With the descriptions given by Eq. (38) to (41), expanded expression of Eq. (37) can be given for Model 2 as

0

0 0

1 ln( )

2 ln( )

x {

} + [

] x { }

x

1

0

2

0

0

[

(

1

+ 1

2

1 )

]

2

1

2

0

( [(

)

2 + (

)

1]

)

+

0

0

2

[]

(2

)

[

]

0

2

1

2

{}

[( [(

0 )

2 + (

0 )

1]

2 ) (2 ) ]

2

1

2

0

( [(

)

2 + (

)

1]

) 0 0

0

0

2

+

[] [0 ] {}

(2

)

[

]

0

=

1 ln( )

2 ln( )

0

+

1

0

2

0

|

{} = { }

([

(

1

+ 1

2

1 )

])

0

=0

(44)

where

=

0 1

1

0 1 2

([

2

1

2

] [

2

1

2

])

[(

)

0

2 + (

)

0

1]

2 2

[(

)

0

2 + (

)

0

1]

2 2

(45)

5.2 Nonlinear Solver

This study employs Newton-Raphson as the nonlinear solver. For this, Eqn.(42) or Eqn.(44) can be arranged in matrix form as

= [()]

1

1

1

1

1

1

1

1

1

1

[[11] [12] 0

(46)

1

1

(53)

[21] [22]] {} = {0}

=

which can be further simplified into the general form

[()]{} = 0 (47)

1

1

The residual can be expressed as

{()} [()]{} (48)

[1

1

]

Expanding Eq. (48) by Taylors series about the known

solution gives

The problem is solved by first solving for the change of variable, symbolically given as

= (54)

{()} = 0 = {()} +

{()}

{} {} (49)

Then the variables are updated by

{}+ = {} + {} (55)

where the series has been truncated up to linear terms only.

Rearranging Eq. (49) gives

[()]{} = {() } (50) where [()] is thus the tangent stiffness given as

The above process will be iterated until the satisfaction of some specified convergence criteria is attained.

1. VALIDATION OF FORMULATIONS

In this study, once the formulations are established, the

[() ] =

{()}

{} (51)

corresponding source codes are written in Matlab. Results obtained are then verified against the analytical solutions. The following subsctions detailed such verifications according to the constitutive laws.

For a vessel with n nodes and for residual {} expressed as (from Eq. (24) and (25))

Model 1

Figures 3(a), (b) and (c) show the plotting of results for

{} = {} + {} = [2

]

the cross sectional area, pressure and flow rate distribution,

+

(52)

respectively, along the vessel. The plots are given for various pressure difference, . Based on the figures, it can

+ [ 0 ] = { }

The expanded form of the tangent stiffness given by Eq.

1. can be given as

be seen that, for relatively lower values of pressure differences, (i.e. <1800Pa), no oscillations are observed. On the other hand, for =1800Pa, slight oscillation is observed for Bubnov Galerkin formulation represented by the wiggling-like curve which becomes greater for higher pressure differences (i.e. =2200Pa, =2500Pa). However, these oscillations vanish when SUPG formulation is employed hence the attainment of the stabilized solutions. This observation marks the success of this study.

1. Area

2. Pressure

3. Flow Rate

Fig. 3 Stabilization of solutions with the employment of SUPG formulation for Model 1

Another trend that can be observed is that, despite the wiggling found in the solutions of pressure and area distributions, flow rate solutions seem not being affected. This highlights that pressure and area solutions are more sensitive to instability than the flow rate solution.

Also, it can be observed that somehow both Bubnov and SUPG formulations diverge from the analytical solutions towards the outlet as the pressure difference increases. This is a typical phenomenon (hence shortcoming) of both finite element formulation which occurs due to sharp internal and boundary layers as identified and studied in Hughes et al. [20] and Tezduyar and Park [21].

For future reference, numerical data of pressure taken at = 0.45 related to the employment of constitutive relation of Model 1 are given in Table 1.

Table 1 Numerical data of pressure distributions taken at = 0.45 for

Model 1

Model 2

Figures 4(a), (b) and (c) give the plotting of results for

Model 2. Based on the figures, similar trends are observed where the oscillations are greater for higher pressure difference for Bubnov formulation which is stabilized when SUPG formulation is employed. Again, no oscillation is observed for flow-rate solution thus confirms the insensitivity of the variable to instability problem, at least in the range of pressure differences considered. The typical phenomenon of divergence of the numerical formulations from the analytical solution due to sharp boundary is also observed.

For future reference, numerical data of pressure taken at = 0.045 related to the employment of constitutive relation of Model 2 are given in Table 2.

Table 2 Numerical data of pressure distributions taken at = 0.045 for

 Inlet Pressure, (mmHg) Outlet Pressure, (mmHg) Pressure, P Present (Analytical) (mmHg) Pressure, P Present (Bubnov) (mmHg) Pressure, Present (SUPG) (mmHg) P 70 60 65.7656 65.7697 65.7696 80 60 72.9638 72.9336 72.9593 90 60 81.7944 81.1998 81.6302 80 85.6916 85.6948 85.6949 100 60 91.5684 88.4568 90.7133 80 92.4244 92.4131 92.4241 105 60 96.6077 91.5647 95.2215 80 96.3631 96.3132 96.3560

Model 2

 Inlet Pressure, (Pa) Outlet Pressure, (Pa) Pressure, P -Present (Analytical) (Pa) Pressure, P -Present (Bubnov) (Pa) Pressure, P -Present (SUPG) (Pa) 400 0 219.0113 221.8988 221.9219 900 0 544.7454 550.6878 551.0865 400 679.5921 683.4366 683.4965 800 856.4529 857.0487 857.0489 1400 0 956.2520 960.1980 964.4542 400 1021.7960 1028.1233 1028.7451 800 1147.6963 1152.2351 1152.3028 1200 1311.3055 1312.8957 1312.9524 1800 0 1332.9789 1309.1983 1343.5355 400 1349.2886 1353.3251 1357.4480 800 1418.0922 1424.6723 1425.1408 1200 1540.6901 1545.3126 1545.3213 2200 0 1753.8284 1631.7019 1756.4283 400 1739.3268 1713.7936 1750.2536 800 1755.8363 1760.0212 1763.9044 1200 1821.4187 1828.0692 1828.3641 2500 0 2067.5523 1831.6639 2057.8781 400 2042.1463 1951.5272 2050.8800 800 2033.6771 2022.2246 2041.8587 1200 2061.0896 2067.0662 2069.0130
1. Area

2. Pressure

3. Flow Rate

Fig. 4 Stabilization of solutions with the employment of SUPG formulation for Model 2

2. SUMMARY AND CONCLUSIONS

In this study, SUPG formulation has been developed for the one-dimensional fluid-structure-interaction (1D-FS1) steady flow that employs pressure-area constitutive relation to complement the mass and the momentum equations of Navier-Stokes. For validation purposes, an analytical solution is derived for one of the constitutive relation. From the study, it was found that SUPG able to provide stable solutions to the problem which otherwise would wiggle due to numerical instability. This study is important as it provides the first SUPG formulation and numerical data for future reference for the specific problem of 1D-FSI employing pressure-area constitutive relation.

REFERENCES

1. A.H. Shapiro, Steady Flow in Collapsible Tubes. Journal of Biomechanical Engineering, 197. 99(3): p. 126-147.

2. R.D. Kamm and A.H. Shapiro, Unsteady Flow in a Collapsible Tube Subjected to External Pressure or Body Forces. Journal of Fluid Mechanics, 1979. 95: p. 1-78.

3. D.N. Ku, M.N. Zeigler, and J.M. Downing, One-Dimensional Steady Inviscid Flow Through a Stenotic Collapsible Tube. Journal of Biomechanical Engineering, 1990. 112(4): p. 444- 450.

4. J.M. Downing and D.N. Ku, Effects of Frictional Losses and Pulsatile Flow on the Collapse of Stenotic Arteries. Journal of Biomechanical Engineering, 1997. 119(3): p. 317-324.

5. R.J. Whittaker, M. Heil, O.E. Jensen, and S.L. Waters, Predicting the onset of high-frequency self-excited oscillations in elastic-walled tubes. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 2010. 466(2124): p. 3635-3657.

6. R.J. Whittaker, M. Heil, O.E. Jensen, and S.L. Waters, A Rational Derivation of a Tube Law from Shell Theory. The Quarterly Journal of Mechanics and Applied Mathematics, 2010. 63(4): p. 465-496.

7. L. Formaggia, D. Lamponi, and A. Quarteroni, One- dimensional models for blood flow in arteries. Journal of Engineering Mathematics, 2003. 47(3-4): p. 251-276.

8. S.J. Sherwin, V. Franke, J. PeirÃ³, and K. Parker, One- dimensional modelling of a vascular network in space-time variables. Journal of Engineering Mathematics, 2003. 47(3-4): p. 217-250.

9. A. Quarteroni and L. Formaggia, Mathematical Modelling and Numerical Simulation of the Cardiovascular System, in Handbook of Numerical Analysis. 2004, Elsevier. p. 3-127.

10. T. Sochi, One-Dimensional Navier-Stokes Finite Element Flow Model, 2013, Imaging Sciences and Biomedical Engineering, King's College London, St Thomas' Hospital, London.

11. T. Sochi, The flow of Newtonian and power law fluids in elastic tubes. International Journal of Non-Linear Mechanics, 2014. 67: p. 245-250.

12. T. Sochi, Navier-Stokes Flow in Cylindrical Elastic Tubes. Journal of Applied Fluid Mechanics (Accepted), 2015. 8(2): p. 181-188.

13. T. Sochi, NavierStokes flow in convergingdiverging distensible tubes. Alexandria Engineering Journal, 2015. 54(3): p. 713-723.

14. J.C. Heinrich, P.S. Huyakorn, O.C. Zienkiewicz, and A.R. Mitchell, An upwind finite element scheme for two- dimensional convective transport equation. International Journal for Numerical Methods in Engineering, 1977. 11(1): p. 131-143.

15. 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, 1982. 32(13): p. 199-259.

16. D.W. Kelly, S. Nakazawa, O.C. Zienkiewicz, and J.C. Heinrich, A note on upwinding and anisotropic balancing dissipation in finite element approximations to convective diffusion problems. International Journal for Numerical Methods in Engineering, 1980. 15(11): p. 1705-1711.

17. K.W. Thompson, Time dependent boundary conditions for hyperbolic systems. Journal of Computational Physics, 1987. 68(1): p. 1-24.

18. J. Donea and A. Huerta, Finite Element Methods for Flow Problems. 2003, England: John Wiley & Sons Ltd. 350.

19. A. SoulaÃ¯mani and M. Fortin, Finite element solution of compressible viscous flows using conservative variables. Computer Methods in Applied Mechanics and Engineering, 1994. 118(34): p. 319-350.

20. T.J.R. Hughes, M. Mallet, and M. Akira, A new finite element formulation for computational fluid dynamics: II. Beyond SUPG. Computer Methods in Applied Mechanics and Engineering, 1986. 54(3): p. 341-355.

21. T.E. Tezduyar and Y.J. Park, Discontinuity-capturing finite element formulations for nonlinear convection-diffusion- reaction equations. Computer Methods in Applied Mechanics and Engineering, 1986. 59(3): p. 307-325.