Unsteady Flow - Shock Capturing Methods
Download
Report
Transcript Unsteady Flow - Shock Capturing Methods
Solution of the St Venant
Equations / Shallow-Water
equations of open channel flow
Dr Andrew Sleigh
School of Civil Engineering
University of Leeds, UK
www.efm.leeds.ac.uk/CIVE/UChile
Shock Capturing Methods
Ability to examine extreme flows
Changes between sub / super critical
Other techniques have trouble with trans-critical
Steep wave front
Front speed
Complex Wave interactions
Alternative – shock fitting
Good, but not as flexible
More recent
Developed from work on Euler equations in
the aero-space where shock capturing is very
important (and funding available)
1990s onwards
Euler equations / Numerical schemes:
Roe, Osher, van Leer, LeVeque, Harten, Toro
Shallow water equations
Toro, Garcia-Navarro, Alcrudo, Zhao
Books
E.F. Toro. Riemann Solvers and Numerical Methods
for Fluid Dynamics. Springer Verlag (2nd Ed.) 1999.
E.F. Toro. Shock-Capturing Methods for FreeSurface Flows. Wiley (2001)
E.F. Toro. Riemann Problems and the WAF Method
for Solving Two-Dimensional Shallow Water
Equations. Philosophical Transactions of the Royal
Society of London. A338:43-68 1992.
Dam break problem
The dam break problem can be solved
It is in fact integral to the solution technique
Conservative Equations
As earlier, but use U for vector
U t F U
A
U
Q
Q
F U Q 2
gI 1
A
x
S U
0
S U
gI gA S S
o
f
2
I1 and I2
Trapezoidal channel
Base width B, Side slope SL= Y/Z
h dS L
1 dB
I2 h
2
dx
3
dx
S
B
I1 h h L
3
2
2
2
Rectangular, SL = 0
2
I1
h B
2
Source term
A
2
I2
2B
Sf
QQn
2
A R
2
(4 / 3)
A
2
2B
dB
2
dx
Rectangular Prismatic
Easily extendible to 2-d
h
U
hu
0
S U
gh S S
o
f
hu
F U 2 1
2
hu gh
2
Integrate over control volume V
U dx F U dt S U dU
V
V
2-dimensions
In 2-d have extra term:
U t F U
h
U hu
hv
hu
2 1
2
F U hu gh
2
huv
friction
S fx
n
h
x G U y
hu
G U
huv
2 1
2
hv
gh
2
2
(1 / 3 )
S U
u u v
2
2
0
S U gh S o x S f x
gh S o S f
y
y
Normal Form
Consider the control volume V, border Ω
n
V
x
U
t
V
dV
H U d S U dU
V
H U cos F U sin G U n1 F U n 2 G U
Rotation matrix
H(U) can be expressed
H U n1 F U n 2 G U T F T U
1
1
T 0
0
0
cos
sin
sin
cos
0
T
1
1
0
0
0
cos
sin
sin
cos
0
h
T U h uˆ
h vˆ
Finite volume formulation
Consider the homogeneous form
i.e. without source terms
U dx F U dt 0
V
And the rectangular control volume in x-t space
t
tn + 1
F i+ 1 /2
F i-1 /2
x i-1 /2
i
x i+ 1 /2
x
Finite Volume Formulation
The volume is bounded by
xi+1/2 and xi-1/2
tn+1 and tn
The integral becomes
t n 1
xi 1 / 2 U x , t n 1 dx xi 1 / 2 U x , t n dx tn F U x i 1 / 2 , t dt
xi 1 / 2
xi 1 / 2
t n 1
tn
F U x i 1 / 2 , t dt
Finite Volume Formulation
We define the integral averages
U
n
i
F i 1 / 2
1
x
x i 1 / 2
1
t
t n 1
tn
U x , t n dx
xi 1 / 2
n 1
Ui
Fi 1 / 2
F U x i 1 / 2 , t dt
1
x
1
t
U x , t n 1 dx
xi 1 / 2
t n 1
x i 1 / 2
tn
F U x i 1 / 2 , t dt
And the finite volume formulation becomes
U
n 1
i
U
n
i
t
x
F i 1 / 2 F i 1 / 2
Finite Volume Formulation
Up to now there has been no approximation
The solution now depends on how we
interpret the integral averages
In particular the inter-cell fluxes
Fi+1/2 and F1-1/2
Finite Volume in 2-D
The 2-d integral equation is
U
t
V
dV
H U d S U dU
V
H(u) is a function normal to the volume
Fn s
A s 1
As
H U d
n F U n G U d
A s 1
As
1
2
Finite Volume in 2-D
Using the integral average
Ui
1
V
U dV
V
Where |V| is the volume (area) of the volume
then
Fn
U
n 1
i
U
n
i
t
V
N
Fn
s 1
s
Finite Volume in 2-D
If the nodes and sides are labelled as :
A2
L1
L2
A1
V
A3
L4
L3
U
n 1
i
U
n
i
t
V
A4
Fn s1 L1 Fn s 2 L 2
Fn s 3 L 3 Fn s 4 L 4
Where Fns1 is normal flux for side 1 etc.
FV 2-D Rectangular Grid
y
For this grid
G j+ 1 /2
i+ 1 /2 ,j+1 /2
i-1 /2 ,j+ 1/2
F i-1 /2
F i+ 1 /2
i, j
i+ 1 /2 ,j-1 /2
i-1 /2 ,j-1 /2
G j-1 /2
x
Solution reduces to
n 1
Ui
Ui
n
t
x
F
i 1 / 2 , j
Fi 1 / 2 , j
t
y
G
i , j 1 / 2
G i , j 1 / 2
Flux Calculation
We need now to define the flux
Many flux definitions could be used to that
satisfy the FV formulation
We will use Godunov flux (or Upwind flux)
Uses information from the wave structure of
equations.
Godunov method
Assume piecewise linear data states
n+ 1
F i-1 /2
F i+ 1 /2
C ells
n
i-1
i-1
i
D ata state s
U (0 ) i-1 /2
U (0 ) i+ 1 /2
n+ 1
n
Means that the flux calc is solution of local
Riemann problem
Riemann Problem
The Riemann problem is a initial value
problem defined by
U t F U
x
0
U in if x x i 1 / 2
U x, tn n
U i 1 if x x i 1 / 2
Solve this to get the flux (at xi+1/2)
FV solution
We have now defined the integral averages
of the FV formulation
U
n 1
i
U
n
i
t
x
F i 1 / 2 F i 1 / 2
The solution is fully defined
First order in space and time
Dam Break Problem
The Riemann problem we have defined is a
generalisation of the Dam Break Problem
D a m w all
D eep w a ter at rest
S hallo w w a ter at rest
Dam Break Solution
Evolution of solution
h
W ater le v els at tim e t= t*
x
v
V elo city at tim e t= t*
x
R are factio n
t
S ho ck
t*
x
Wave structure
Exact Solution
Toro (1992) demonstrated an exact solution
Considering all possible wave structures a
single non-linear algebraic equation gives
solution.
Exact Solution
Consider the local Riemann problem
U t F U x 0
h
U hu
hv
U L if x 0
U x ,0
U R if x 0
hu
2 1
2
F U hu gh
2
huv
Wave structure
S tar regio n, h *, u *
L e ft w a ve ( u -c)
S hea r w a v e
t (u)
vL
R ight w a ve, ( u+ c)
vR
h L, u L, v L
hR, uR, vR
0
x
Possible
Wave structures
Across left and
right wave h, u
change v is
constant
S hea r w a v e
L e ft R are fa ctio n
R ight S ho c k
t
x
L e ft S ho c k
S hea r w a v e
R ight R are fa ctio n
t
x
L e ft R are fa ctio n
S hea r w a v e
R ight R are fa catio n
t
Across shear
wave v changes,
h, u constant
x
L e ft S ho c k
R ight S ho c k
S hea r w a v e
t
x
Determine which wave
Which wave is present is determined by the
change in data states thus:
h* > hL
h* ≤ hL
left wave is a shock
left wave is a rarefaction
h* > hR
h* ≤ hR
right wave is a shock
right wave is a rarefaction
Solution Procedure
Construct this equation
f h f L h , hL f L h , hL u
And solve iteratively for h (=h*).
The functions may change each iteration
f(h)
The function f(h) is defined
f h f L h , hL f L h , hL u
u uR uL
2 gh gh L
fL
1 h hL
h
h
g
L
2 hh L
And u*
u
*
if
2 gh gh R
fR
1 h hR
h
h
g
R
2 hh R
1
2
u L
uR
h h L ( rarefactio n )
if
h h L ( shock )
h h R ( rarefactio n )
if
if
1
2
h h R ( shock )
f h
R
*
, hR f L h , hL
*
Iterative solution
The function is well behaved and solution by
Newton-Raphson is fast
(2 or 3 iterations)
One problem – if negative depth calculated!
This is a dry-bed problem.
Check with depth positivity condition:
u u R u L 2 c L c R
Dry–Bed solution
Dry bed on one side
of the Riemann
problem
Dry bed evolves
Wave structure is
different.
t
D ry b ed
W et b ed
x
t
D ry b ed
W et b ed
x
D ry b ed
t
W et b ed
W et b ed
x
Dry-Bed Solution
Solutions are explicit
Need to identify which applies – (simple to do)
Dry bed to right
u*
3
u L
2cL
c*
1
c*
1
3
Dry bed to left
u*
1
1
3
u R
2cR
Dry bed evolves
u L
3
2cL
u R
h* c / g
2cR
h* = 0 and u* = 0
Fails depth positivity test
2
*
Shear wave
The solution for the shear wave is straight
forward.
If vL >
Else
0
v* = vL
v* = vR
Can now calculate inter-cell flux from h*, u*
and v*
For any initial conditions
Complete Solution
The h*, u* and v* are sufficient for the Flux
But can use solution further to develop exact
solution at any time.
i.e. Can provide a set of benchmark solution
Useful for testing numerical solutions.
Choose some difficult problems and test your
numerical code again exact solution
Complete Solution
h
W ater le v els at tim e t= t*
x
v
V elo city at tim e t= t*
x
R are factio n
t
S ho ck
t*
x
Difficult Test Problems
Toro suggested 5 tests
Test No.
hL (m)
uL (m/s)
hR (m)
uR (m/s)
1
1.0
2.5
0.1
0.0
2
1.0
-5.0
1.0
5.0
3
1.0
0.0
0.0
0.0
4
0.0
0.0
1.0
0.0
5
0.1
-3.0
0.1
3.0
Test 1: Left critical rarefaction and right shock
Test 2: Two rarefactions and nearly dry bed
Test 3: Right dry bed problem
Test 4: Left dry bed problem
Test 5: Evolution of a dry bed
Exact Solution
Consider the local Riemann problem
U t F U x 0
h
U hu
hv
U L if x 0
U x ,0
U R if x 0
hu
2 1
2
F U hu gh
2
huv
Wave structure
S tar regio n, h *, u *
L e ft w a ve ( u -c)
S hea r w a v e
t (u)
vL
R ight w a ve, ( u+ c)
vR
h L, u L, v L
hR, uR, vR
0
x
Returning to the Exact Solution
We will see some other Riemann solvers that
use the wave speeds necessary for the exact
solution.
Return to this to see where there come from
Possible
Wave structures
Across left and
right wave h, u
change v is
constant
S hea r w a v e
L e ft R are fa ctio n
R ight S ho c k
t
x
L e ft S ho c k
S hea r w a v e
R ight R are fa ctio n
t
x
L e ft R are fa ctio n
S hea r w a v e
R ight R are fa catio n
t
Across shear
wave v changes,
h, u constant
x
L e ft S ho c k
R ight S ho c k
S hea r w a v e
t
x
Conditions across each wave
t
Left Rarefaction wave
t
Right bounding characteristic
Left bounding characteristic
dx / dt u * c *
dx / dt u L c L
hL, uL
h*, u*
x
Smooth change as move in x-direction
Bounded by two (backward) characteristics
Discontinuity at edges
Crossing the rarefaction
We cross on a forward characteristic
u 2 c constant
States are linked by:
u L 2 c L u * 2 c*
or
u * u L 2 c L c *
Solution inside the left rarefaction
dx
The backward characteristic equation is dt u c
For any line in the direction of the rarefaction
Crossing this the following applies: u L 2 c L u 2 c
1
dx
c u L 2cL
3
dt
Solving gives
On the t axis dx/dt = 0
c
1
3
u L
2cL
u
u
1
3
u L
1
dx
u
2
c
L
L
3
dt
2cL
Right rarefaction
dx
Bounded by forward characteristics
Cross it on a backward characteristic
u R 2 c R u * 2 c*
In rarefaction
If Rarefaction crosses axis
c
1
3
c
u R
2cR
dt
u * u R 2 c * c R
1
dx
u
2
c
R
R
3
dt
uc
u
1
3
u
u R
1
dx
u R 2cR 2
3
dt
2cR
Shock waves
Two constant data states are separated by a
discontinuity or jump
Shock moving at speed Si
Using Conservative flux for left shock
UL
hL
h
u
L L
h*
U*
h
u
* *
Conditions across shock
Rankine-Hugoniot condition
F U * F U L S i U * U L
Entropy condition
i U L S i i U *
λ1,2 are equivalent to characteristics.
They tend towards being parallel at shock
Shock analysis
Change frame of reference, add Si
uˆ * u * S L
uˆ L u L S L
hL
Uˆ L
ˆ
h
u
L L
h*
ˆ
U*
ˆ
h
u
* *
Rankine-Hugoniot gives
h* uˆ * h L uˆ L
h* uˆ
2
*
1
2
gh
2
*
h L uˆ
2
L
1
2
2
gh L
Shock analysis
Mass flux conserved
From eqn 1
M
L
Using
uˆ * M
L
/ h*
uˆ L M
L
L
M
L
also
uˆ * uˆ L u * u L
M
h*2 h L2
g
2 uˆ * uˆ L
1
M
h* uˆ * h L uˆ L
/ hL
1
2
h*2 h L2
g
2 u * u L
1
L
g h* h L h* h L
Left Shock Equation
Equating gives
u * u L f L h* , h L
Also
S L uL aLqL
f L h* , h L h* h L
qL
1 h* h L h*
2
2
hL
h hL
g *
2 h* h L
1
Right Shock Equation
Similar analysis gives
u * u R f R h* , h R
Also
S R uR aRqR
f R h* , h R h* h R
qR
1 h* h R h*
2
2
hR
h hR
g *
2 h* h R
1
Complete equation
Equating the left and right equations for u*
u * u L f L h* , h L
u * u R f R h* , h R
u R u L f L h , hL f R h , hL 0
Which is the iterative of the function of Toro
f h* f L h* , h L f R h* , h L u 0
Approximate Riemann Solvers
No need to use exact solution
Expensive
Iterative
When other equations, exact may not exist
Many solvers
Some more popular than others
Toro Two Rarefaction Solver
Assume two rarefactions
Take the left and right equations
u * u L 2 c L c *
u * u R 2 c * c R
Solving gives
c*
uL uR
4
cL cR
2
For critical rarefaction use solution earlier
Toro Primitive Variable Solver
Writing the equations in primitive variables
Non conservative
h
W
u
W t A W W x 0
h
u
Approximate A(W) by a constant matrix
W
1
2
u
A W
g
W R
WL
Gives solution
h*
h L
hR
2
u R
u L h L h R
4 c L c R
u*
u L
uR
2
h R
h L c L c R
4 h L h R
Toro Two Shock Solver
Assuming the two waves are shocks
h*
u*
qL
1
2
q L hL q R hR u L u R
qL qR
1
u L u R h* h R q R h* h L q L
2
g ho h L
2 ho h L
qR
g ho h R
2 ho h R
Use two rarefaction solver to give h0
Roe’s Solver
Originally developed for Euler equations
Approximate governing equations with:
U t F U
Where
uL
u~
~
A
hL u R
hL
hL
x
U t AU
x
~
U t AU x
is obtained by Roe averaging
hR
~
h
hL hR
c~
1
2
c
2
L
cR
2
Roe’s Solver
Properties of matrix
Eigen values
Right eigen vectors
Wave strengths
Flux given by
~
1 u~ c~
~
1 u~ c~
~ (2) 1
R
~ ~
u c
~ (1 ) 1
R ~ ~
u c
~
1
h
~1 h ~ u
2
c
Fi 1 / 2
1
2
F
n
i
F
~
1
h
~ 2 h ~ u
2
c
n
i 1
2
~
2
1
j 1
~ ~ j
j
j R
h hR hL
HLL Solver
Harten, Lax, van Leer
Assume wave speed
Construct volume
Integrate round
SL
t
xL
SR
t
xR
t
U*
FL
FR
F*
t
UL
UR
x
xL
xR
x LU L x RU R x L x R U * tU R tU L 0
U * S R U R S LU L F L F R / S R S L
Alternative gives flux
F* S R F L S L F R S L S R U R U L / S R S L
HLL Solver
What wave speeds to use?
Free to choose
One option:
S L min u L c L , u TR c TR
For dry bed (right)
Simple, but robust
S R min u R c R , u TR c TR
S L u L cL
S R u R 2cR
Higher Order in Space
Construct Riemann problem using cells
further away
Danger of oscillations
Piecewise
reconstruction
L
R
i-1
i
Limiters
i+ 1
i+ 2
i+ 1
i+ 2
L
R
i-1
i
Limiters
Obtain a gradient for variable in cell i, Δi
UL Ui
2
x i
U R U i 1
1
2
x i 1
Gradient obtained from Limiter functions
Provide gradients at cell face
i 1 / 2 a
1
u i 1 u i
xi 1 / 2 xi
Limiter Δi =G(a,b)
i 1 / 2 b
u i u i 1
xi xi 1
Limiters
A general limiter
max 0 , min a , b , min a , b
G a , b
min 0 , max a , b , max a , b
for a 0
for a 0
β=1 give MINMOD
β=2 give SUPERBEE
van Leer
G a , b
ab ab
a b
Other SUPERBEE expression
G a , b s max 0 , min 2 b , s a , min b , 2 s a
s = sign
Higher order in time
Needs to advance half time step
MUSCL-Hancock
Primitive variable
W A W W
Limit variable
Evolve the cell face values 1/2t:
t
Wi
L,R
Wi
L,R
x
0
1 t
2 x
A Wi
n
W
Update as normal solving the Riemann
problem using evolved WL, WR
n 1
Ui
Ui
n
t
x
F i 1 / 2 F i 1 / 2
L
i
Wi
R
Alternative time advance
Alternatively, evolve the cell values using
n 1 / 2
Ui
Ui
1 t
2 x
F W F W
R
i
L
i
Solve as normal
n 1
Ui
n
Ui
n
t
x
F i 1 / 2 F i 1 / 2
Procedure is 2nd order in space and time
Boundary Conditions
Set flux on boundary
Directly
Ghost cell
Wall u, v = 0.
Transmissive
Ghost cell un+1=-un
Ghost cell hn+1 = hn
un+1 = un
Wet / Dry Fonts
Wet / Dry fronts are difficult
Source of error
Source of instability
Common
near tidal boundaries
Flooding - inundation
Dry front speed
We examined earlier dry bed problem
S *L u c
SL
t
u L 2cL u 2c
W et b ed
D ry b ed
x
c0
Front is fast
S *L u L
Faster than characteristic.
Can cause problem with time-step / Courant
2cL
Solutions
The most popular way is to artificially wet bed
Give a small depth, zero velocity
Loose a bit of mass and/or momentum
Can drastically affect front speed
E.g. a 1.0 dambeak with 1cm gives 38% error
1mm gives 25% error – try it!
Conservation Errors
The conserved variable are h and hu
Often require u
Need to be very careful about divide by zero
Artificial dry-bed depths could cause this
Source Terms
“Lumped” in to one term and integrated
Attempts at “upwinding source”
Current time-step
Could use the half step value
E.g.
0
S z u u
2
gh
x K
2 B
2 B x
K CR
K
R
1/ 2
2/3
n
Main Problem is Slope Term
Flat still water over uneven bed starts to
move.
z
Problem with discretisation of
gh
i-1
i
x
i+ 1
w ater surfa ce
hi 1
hl
hi
hr
hi 1
bed level
gh i
z, h
zi 1
zl
zi
zr
z i 1
x
datum le vel
zr
zl
x
Discretisation
Discretised momentum eqn
hu i
n 1
hu i
2
2
t
gh l
gh r
2
2
hu l hu r
x
2
2
For flat, still water
hu i
n 1
t
x gh i z r z l
Require
2
t hl
g
hi z l
x 2
2
hl
2
h r2
2 hi z r
2
hi z l
hr
2
hi z r
0
A solution
Assume a “datum” depth, measure down
2
g hi
z
x
g hi
2
hi
x
1 hi
2
g
2 x
and
g
1 h i z i z r
h i z i z l
2
x
2
2
Momentum eqn:
hl hi z i z l
Flat surface
hu i
'
g t
2 x
h
2
l
hr hi z i z r
h r hi z i z r hi z i z l
2
2
2
Some example solutions
Weighted mesh gives more detail for same number of cells
Dam Break - CADAM
Channel with 90° bend
Secondary Shocks