 Open Access
 Total Downloads : 124
 Authors : Aik Ying Tang, Mohd AlAkhbar Mohd Noor, Airil Yasreen Mohd Yassin, Norsarahaida Saidina Amin
 Paper ID : IJERTV6IS040777
 Volume & Issue : Volume 06, Issue 04 (April 2017)
 Published (First Online): 27042017
 ISSN (Online) : 22780181
 Publisher Name : IJERT
 License: This work is licensed under a Creative Commons Attribution 4.0 International License
OneDimensional Fluid – Structure – Interaction (1DFSI) Steady Flow Stable Formulation
A. Y. Tang1 , N.Amin4
Department of Mathematical Sciences, Faculty of Science,
Universiti Teknologi Malaysia, 81310 Johor Bahru, Johor

A. Mohd Noor2 Faculty of Civil Engineering, Universiti Teknologi Malaysia,
81310 Johor Bahru, Johor

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 onedimensional fluidstructureinteraction (1DFSI) formulation which has been formulated by providing a pressurearea 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 BubnovGalerkin 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 pressurearea constitutive relation. This study can be important for future works in 1DFSI employing pressurearea constitutive relation.
Keywords StreamlineUpwindPetrovGalerkin, Finite Element Method, Biomechanics

INTRODUCTION
Fluidstructureinteraction (FSI) for onedimensional flow can be formulated by providing a pressurearea constitutive relation to complement the mass and linear momentum equations. Such coupling would allow the interaction between volumetric flow rate, Q, crosssectional 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 wellknown stabilization scheme, StreamlineUpwindPetrovGalerkin (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 onedimensional fluidstructureinteraction (1DFSI) flow employing pressure area constitutive relation.

GOVERNING EQUATIONS
1DFSI 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 crosssectional area of the flow. Various detailed forms of () have been proposed in the literature [113] To note, pressurearea constitutive relation as in Eqn. (1) is also termed as tube law elsewhere [56]. 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], BubnovGalerkin 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 pressurearea 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 experimentalfit parameters) is thus considered and termed herein as Model 2 (Eqn. (9)).
Model 1
Pressurearea constitutive relation used in Sherwin et al. [8], Quarteroni and Formaggia [9], and Sochi [1013] is termed as
Model 1 herein and given as
2ln()
2 ])
2 1

SPURIOUS OSCILLATIONS
In numerical analysis of fluid dynamics, spurious oscillations can occur in flows with high Peclet number (for advectiondiffusion problems) and high Reynolds number (for general flows) when solved using either central finite difference method or BubnovGalerkin 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 BubnovGalerkin
() =
0
( 0 ) (6)
formulation detailed in Sochi [12], whilst we obtained similar nonoscillatory 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
Pressurearea 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 BubnovGalerkin formulation to a different pressurearea 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 1DFSI flow as the typical phenomenon of
where is the constant of integration which can be determined from the boundary condition where = at
= 0, thus
BubnovGalerkin 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

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

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 ) + [( )21 1]
1,2
L (H U + B) = 0 (30)
where , are the lefteigenvectors. 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 NewtonRaphson 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)
[11
]
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 asThe above process will be iterated until the satisfaction of some specified convergence criteria is attained.

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.

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

Area

Pressure

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

Area

Pressure

Flow Rate


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


SUMMARY AND CONCLUSIONS
In this study, SUPG formulation has been developed for the onedimensional fluidstructureinteraction (1DFS1) steady flow that employs pressurearea constitutive relation to complement the mass and the momentum equations of NavierStokes. 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 1DFSI employing pressurearea constitutive relation.
REFERENCES

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

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

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

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

R.J. Whittaker, M. Heil, O.E. Jensen, and S.L. Waters, Predicting the onset of highfrequency selfexcited oscillations in elasticwalled tubes. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 2010. 466(2124): p. 36353657.

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

L. Formaggia, D. Lamponi, and A. Quarteroni, One dimensional models for blood flow in arteries. Journal of Engineering Mathematics, 2003. 47(34): p. 251276.

S.J. Sherwin, V. Franke, J. PeirÃ³, and K. Parker, One dimensional modelling of a vascular network in spacetime variables. Journal of Engineering Mathematics, 2003. 47(34): p. 217250.

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

T. Sochi, OneDimensional NavierStokes Finite Element Flow Model, 2013, Imaging Sciences and Biomedical Engineering, King's College London, St Thomas' Hospital, London.

T. Sochi, The flow of Newtonian and power law fluids in elastic tubes. International Journal of NonLinear Mechanics, 2014. 67: p. 245250.

T. Sochi, NavierStokes Flow in Cylindrical Elastic Tubes. Journal of Applied Fluid Mechanics (Accepted), 2015. 8(2): p. 181188.

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

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

A.N. Brooks and T.J.R. Hughes, Streamline upwind/Petrov Galerkin formulations for convection dominated flows with
particular emphasis on the incompressible NavierStokes equations. Computer Methods in Applied Mechanics and Engineering, 1982. 32(13): p. 199259.

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

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

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

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

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

T.E. Tezduyar and Y.J. Park, Discontinuitycapturing finite element formulations for nonlinear convectiondiffusion reaction equations. Computer Methods in Applied Mechanics and Engineering, 1986. 59(3): p. 307325.