### ppt - The Institute for Computational Engineering and Sciences (ICES)

Structure-Preserving B-spline Methods for the
Incompressible Navier-Stokes Equations
John Andrew Evans
Institute for Computational Engineering and Sciences, UT Austin
Stabilized and Multiscale Methods in CFD
Spring 2013
Motivation
So why do we need another flow solver?
Incompressibility endows the Navier-Stokes equations with
important physical structure:
• Mass balance
• Momentum balance
• Energy balance
• Enstrophy balance
• Helicity balance
However, most methods only satisfy the incompressibility
constraint in an approximate sense.
Such methods do not preserve structure and lack robustness.
Motivation
Consider a two-dimensional Taylor-Green vortex. The vortex is a
smooth steady state solution of the Euler equations.
1
0.9
0.8
¶u
+ Ñ × ( u Ä u ) + Ñp = f
¶t
Ñ×u = 0
0.7
y
0.6
0.5
0.4
with Symmetry BCs
0.3
0.2
0.1
0
0
0.2
0.4
x
0.6
0.8
1
Motivation
Conservative methods which only weakly satisfy incompressibility
do not preserve this steady state and blow up in the absence of
artificial numerical dissipation.
2-D Conservative Taylor-Hood Element
Q2/Q1 and h = 1/8
Motivation
Methods which exactly satisfy incompressibility are stable in the
Euler limit and robust with respect to Reynold’s number.
Re = ¥
Re = 10
Increasing
Reynold’s
number
2-D Conservative Structure-preserving B-splines
k’ = 1 and h = 1/8
Motivation
Due to the preceding discussion, we seek new discretizations that:
• Satisfy the divergence-free constraint exactly.
• Harbor local stability and approximation properties.
• Possess spectral-like stability and approximation properties.
• Extend to geometrically complex domains.
Structure-Preserving B-splines seem to fit the bill.
The Stokes Complex
The classical L2 de Rham complex is as follows:
¾¾
® H ¾¾¾
® H(curl) ¾¾® H(div) ¾¾® L2 ¾¾
®0
1
curl
div
From the above complex, we can derive the following smoothed
complex with the same cohomology structure:

   
     V   Q  
0
 : H
V : H
1
1
curl
{
div
F := v ÎL2 : curl v ÎH1
Q : L
2
}
The Stokes Complex
The smoothed complex corresponds to viscous flow, so we
henceforth refer to it as the Stokes complex.

   
     V   Q  
0
Scalar
Potentials
curl
Vector
Potentials
 

div
Flow
Pressures
Flow
Velocities
 
curl
 
div
The Stokes Complex
For simply-connected domains with connected boundary, the Stokes
complex is exact.
• grad operator maps onto space of curl-free functions
grad      : curl   
• curl operator maps onto space of div-free functions
curl  v  V : divv   
• div operator maps onto entire space of flow pressures
divV  Q
The Stokes Complex
 grad ( )  d s   ( (1))   ( (0 ))

 curl( )  d a     d s
Curl Theorem:
S
S
 div(v ) dV
Divergence Theorem:

V
V
 

 
curl
 v da
 
div
The Discrete Stokes Subcomplex
Now, suppose we have a discrete Stokes subcomplex.

  h  
  h    Vh   Q h  
0
Discrete
Scalar
Potentials
Discrete
Vector
Potentials
curl
div
Discrete Flow
Velocities
Discrete Flow
Pressures
Then a Galerkin discretization utilizing the subcomplex:
• Does not have spurious pressure modes, and
• Returns a divergence-free velocity field.
The Discrete Stokes Subcomplex
P roposition. If v h  V h satisfies
div
v h , q h L2  0,  q h  Q h
then div v h  0.
Proof. Let q h  div v h  Q h .
 d iv v h
2
L
2
 0
The Discrete Stokes Subcomplex
Proposition. For every qh Î Xh3 there exists a v h Î Xh2
such that
di v v h = qh
and
v h 1 £ b h-1 qh
0
where b h only depends on the norm of P2h .
The Discrete Stokes Subcomplex
Proof. For qh ÎXh3 , let v ÎH1 be such that
div v = qh
and
v 1 £ b -1 qh 0 .
Define v h = P2h v. Then
div v h = div P2h v = Ph3 div v = qh
and
v h 1 £ P2h v 1 £ b -1 P2h qh
0
º b h-1 qh 0 .
Structure-Preserving B-splines
Review of Univariate (1-D) B-splines:
Knot vector on (0,1) and k-degree
Ξ = {0, 0, 0, 0.2, 0.4, 0.6, 0.8, 0.8, 1, 1, 1}, k = 2
Knots w/
multiplicity
B-spline basis on (0,1) by recursion:
ì1 if xi £ x < xi +1 ,
N i,0 (x ) = í
î0 otherwise
x - xi
N i,k (x ) =
N (x ) +
xi+k - xi i,k-1
x i+k+1 - x
N i+1,k-1 (x )
xi+k+1 - xi+1
Start w/ piecewise
constants
Bootstrap
recursively
to k
Structure-Preserving B-splines
Review of Univariate (1-D) B-splines:
• Open knot vectors:
• Multiplicity of first and last knots is k+1
• Basis is interpolatory at these locations
• Non-uniform knot spacing allowed
• Continuity at interior knot a function of knot repetition
k=2
Structure-Preserving B-splines
Review of Univariate (1-D) B-splines:
• Derivatives of B-splines are B-splines
1
1
d/dx
onto
0
0
1/3
2/3
k=4
1
0
0
1/3
2/3
k=3
1
Structure-Preserving B-splines
Review of Univariate (1-D) B-splines:
• We form curves in physical space by taking weighted sums.
4
5
3
2
0
“control mesh”
1
- control points
- knots
Structure-Preserving B-splines
Review of Multivariate B-splines:
• Multivariate B-splines are built through tensor-products
• Multivariate B-splines inherit all of the aforementioned
properties of univariate B-splines
In what follows, we denote the space of n-dimensional tensorproduct B-splines as
polynomial degree in direction i
k1 ,k2 ,...,ki ,...,kn
a1 ,a 2 ,...,a i ,...,a n
S
i th continuity vector
Structure-Preserving B-splines
Review of Multivariate B-splines:
• We form surfaces and volumes using weighted sums of
multivariate B-splines (or rational B-splines) as before.
Control mesh
Mesh
Structure-Preserving B-splines
Two-dimensional Structure-Preserving B-splines
Define for the unit square:
k ,k
Sˆ h : S  11 ,22
k ,k 1
k 1,k
Rˆ h : S  11 ,22 1  S  11 1,22
k 1,k 1
Wˆ h : S  11 1,22 1
In the context of fluid flow:
Sˆ h  Sp ace o f Stream Fu n ctio n s
Rˆ h  Sp ace o f Flo w V elo cities
Wˆ h  Sp ace o f Pressu res
and it is easily shown that:
curl
ˆ  div

 Sˆ h    R
 Wˆ h  
0
h
Structure-Preserving B-splines
Two-dimensional Structure-Preserving B-splines:
Mapped Domains
On mapped domains, the Piola transform is utilized to map
flow velocities. Pressures are mapped using an integral
preserving transform.
Structure-Preserving B-splines
Two-dimensional Structure-Preserving B-splines
We associate the degrees of freedom of structure-preserving
B-splines with the control mesh. Notably, we associate:
Sˆ h : co n tro l p o in ts
Rˆ h : co n tro l faces
Wˆ h : co n tro l cells
Structure-Preserving B-splines
Two-dimensional Structure-Preserving B-splines,
k1 = k2 = 2:
Structure-Preserving B-splines
Two-dimensional Structure-Preserving B-splines,
k1 = k2 = 2:
Structure-Preserving B-splines
Two-dimensional Structure-Preserving B-splines,
k1 = k2 = 2:
Structure-Preserving B-splines
Two-dimensional Structure-Preserving B-splines,
k1 = k2 = 2:
Structure-Preserving B-splines
Two-dimensional Structure-Preserving B-splines,
k1 = k2 = 2:
Structure-Preserving B-splines
Two-dimensional Structure-Preserving B-splines,
k1 = k2 = 2:
Structure-Preserving B-splines
Two-dimensional Structure-Preserving B-splines,
k1 = k2 = 2:
Structure-Preserving B-splines
Two-dimensional Structure-Preserving B-splines,
k1 = k2 = 2:
Structure-Preserving B-splines
Three-dimensional Structure-Preserving B-splines
Define for the unit cube:
k ,k ,k
Sˆ h : S  11 ,22 ,3 3
k 1,k ,k
k ,k 1,k
k ,k ,k 1
Nˆ h : S  11 1,22 ,3 3  S  11 ,22 1,3 3  S  11 ,22 ,3 3 1
k ,k 1,k 1
k 1,k ,k 1
k 1,k 1,k
Rˆ h : S  11 ,22 1,3 3 1  S  11 1,22 ,3 3 1  S  11 1,22 1,3 3
k 1,k 1,k 1
Wˆ h : S  11 1,22 1,3 3 1
It is easily shown that:
curl
div

 Sˆ h  
 Nˆ h    Rˆ h   Wˆ h  
0
Structure-Preserving B-splines
Three-dimensional Structure-Preserving B-splines
Flow velocities: map w/ divergence-preserving transformation
Flow pressures: map w/ integral-preserving transformation
Vector potentials: map w/ curl-conserving transformation

h
h
S h : w : w oF  Sˆ h

: u
N h : v :  D F 
Rh
W h :

h
h
T
v
h

oF   Nˆ h
1
: det  D F  D F 

h
u
 o F   Rˆ h
h
h
p : det( D F ) p oF  Wˆ h


Structure-Preserving B-splines
Three-dimensional Structure-Preserving B-splines
We associate as before the degrees of freedom with the
control mesh.
Sˆ h
:
control points
Nˆ h
:
con trol ed ges
Rˆ h
:
con tro l faces
Wˆ h
:
control cells
Structure-Preserving B-splines
Three-dimensional Structure-Preserving B-splines
Control points:
Scalar potential DOF
Control edges:
Vector potential DOF
Structure-Preserving B-splines
Weak Enforcement of No-Slip BCs
Nitsche’s method is utilized to weakly enforce the no-slip
condition in our discretizations. Our motivation is threefold:
• Nitsche’s method is consistent and higher-order.
• Nitsche’s method preserves symmetry and ellipticity.
• Nitsche’s method is a consistent stabilization procedure.
Furthermore, with weak no-slip boundary conditions, a
conforming discretization of the Euler equations is obtained
in the limit of vanishing viscosity.
Structure-Preserving B-splines
Weak Enforcement of Tangential Continuity Between Patches
On multi-patch geometries, tangential continuity is enforced
weakly between patches using a combination of the
symmetric interior penalty method and upwinding.
!
1
F1
!
F2
"2
!
3
2
1
x2
F3
0
0
1
"1
x1
Summary of Theoretical Results
(
(
)
)
Ñ × u Ä u + Ñp = Ñ × 2nÑsu + f in W
Ñ ×u = 0
u=0
in W
on ¶W
• Well-posedness for small data
• Optimal velocity error estimates and suboptimal, by one order,
pressure error estimates
• Conforming discretization of Euler flow obtained in limit of
vanishing viscosity (via weak BCs)
• Robustness with respect to viscosity for small data
Summary of Theoretical Results
(
)
(
)
¶u/ ¶t + Ñ × u Ä u + Ñp = Ñ × 2nÑsu + f in W ´ (0,T )
Ñ ×u = 0 in W ´ (0,T )
u = 0 on ¶W ´ (0,T )
u = u0 on W
• Existence and uniqueness (well-posedness)
• Optimal velocity error estimates in terms of the L2 norm for
domains satisfying an elliptic regularity condition (local-in-time)
• Convergence to suitable weak solutions for periodic domains
• Balance laws for momentum, energy, enstrophy, and helicity
• Balance law for angular momentum on cylindrical domains
Spectrum Analysis
Spectrum Analysis
Consider the two-dimensional periodic Stokes
eigenproblem:


   2  u   p   u in 
s
 u  0
in 
w / P eriodic B C s
We compare the discrete spectrum for a specified
discretization with the exact spectrum. This analysis sheds
light on a given discretization’s resolution properties.
Spectrum Analysis
Spectrum Analysis:
Structure-Preserving B-splines
Spectrum Analysis
Spectrum Analysis:
Taylor-Hood Elements
Spectrum Analysis
Spectrum Analysis:
MAC Scheme
Selected Numerical Results
Numerical Confirmation of Convergence Rates
é
ù
2e x (-1+ x)2 x 2 (y 2 - y)(-1+ 2y)
u(x,y) = ê
x
2 2 ú
ë(-e (-1+ x)x(-2 + x(3 + x))(-1+ y) y ) û
-424 + 156e + (y 2 - y)(-456 + e x (456 +
p(x,y) = x 2 (228 - 5(y 2 - y)) + 2x(-228 + (y 2 - y)) +
2x3 (-36 + (y 2 - y)) + x 4 (12 + (y 2 - y))))
2-D Manufactured Vortex Solution
Selected Numerical Results
Numerical Confirmation of Convergence Rates
2-D Manufactured Vortex Solution
Selected Numerical Results
Numerical Confirmation of Convergence Rates
2-D Manufactured Vortex Solution
Selected Numerical Results
Numerical Confirmation of Convergence Rates
Re
0
1
10
100
1000
10000
Energy
1.40e-2
1.40e-2
1.40e-2
1.40e-2
1.40e-2
1.40e-2
H1 error u
1.40e-2
1.40e-2
1.40e-2
1.40e-2
1.40e-2
1.40e-2
L2 error - u
2.28e-4
2.28e-4
2.28e-4
2.28e-4
2.28e-4
2.28e-4
L2 error - p
3.49e-4
3.49e-4
1.98e-4
1.96e-4
1.96e-4
1.96e-4
Robustness with respect to Reynolds number
k’ = 1 and h = 1/16
2-D Manufactured Vortex Solution
Selected Numerical Results
Numerical Confirmation of Convergence Rates
Re
0
1
10
100
1000
10000
H1 error u
6.77e-4
6.77e-4
7.11e-4
2.26e-3
2.16e-2
X
L2 error u
6.54e-4
6.54e-6
6.79e-6
1.97e-5
1.86e-4
X
L2 error p
1.96e-4
1.96e-4
1.96e-4
1.96e-4
1.96e-4
X
Instability of 2-D Taylor-Hood with respect to
Reynolds number
Q2/Q1 and h = 1/16
2-D Manufactured Vortex Solution
Selected Numerical Results
Lid-Driven Cavity Flow
U
u ÎW1,r , p ÎW 0,r , 1< r < 2
H
H
Selected Numerical Results
Lid-Driven Cavity Flow at Re = 1000
Selected Numerical Results
Lid-Driven Cavity Flow at Re = 1000
Method
umin
vmin
vmax
k’ = 1, h = 1/32
-0.40140
-0.39132
0.54261
k’ = 1, h = 1/64
-0.39399
-0.38229
0.53353
k’ = 1, h = 1/128
-0.39021
-0.37856
0.52884
k’ = 2, h = 1/64
-0.38874
-0.37715
0.52726
k’ = 3, h = 1/64
-0.38857
-0.37698
0.52696
Converged
-0.38857
-0.37694
0.52707
Ghia, h = 1/156
-0.38289
-0.37095
0.51550
Selected Numerical Results
Flow Over a Cylinder at Re = 100
L
H
D
U
Lout
Selected Numerical Results
Flow Over a Cylinder at Re = 100
Patch 1
Patch 4
Patch 2
Patch 3
Patch 5
Selected Numerical Results
Flow Over a Cylinder at Re = 100
Simulation Details:
- Full Space-Time Discretization
- Method of Subgrid Vortices
- Linears in Space and Time
- D = 2, H = 32D, L = 64D
- Time Step Size: 0.25
- Trilinos Implementation
- GMRES w/ ILU Preconditioning
- 3000 time-steps (shedding
initiated after 1000 time-steps)
- Approximately 15,000 DOF/time step
Selected Numerical Results
Flow Over a Cylinder at Re = 100
Selected Numerical Results
Flow Over a Cylinder at Re = 100
Quantity of Interest:
Strouhal Number: St = fD/U
Computed Strouhal Number: 0.163
Accepted Strouhal Number: 0.164
Selected Numerical Results
Three-Dimensional Taylor-Green Vortex Flow
3-D Periodic Flow
Simplest Model
of Vortex Stretching
No External Forcing
Selected Numerical Results
Three-Dimensional Taylor-Green Vortex Flow
Time Evolution of
Dissipation Rate
Reproduced with
permission from
[Brachet et al. 1983]
Selected Numerical Results
Three-Dimensional Taylor-Green Vortex Flow at Re = 200
Enstrophy isosurface at
time corresponding to
maximum dissipation
Selected Numerical Results