A novel VIP (Vector Image Polygon) slope limiter for scalar

Document Sample

```					     “Numerical Methods for multi-Material Fluid Flows”, Arcachon, France,5-9 Sept 2011

A novel VIP (Vector Image Polygon)
slope limiter for scalar variables in
multidimensional grids
Gabi Luttwak1
and
Joseph Falcovitz2

1Rafael,   P.O. Box 2250, Haifa 31021, Israel
2Institute   of Mathematics, The Hebrew University of Jerusalem, Israel
VIP-Convex Hull based limiters
 What is VIP?
– the natural extension of the monotonicity from scalar to vector variables
 VIP => for the Lagrangian phase of the SMG/Q scheme (Pavia
2009)
 VIP => for the advection phase of ALE/Euler SMG
calculations (Paris 2010)
 VIP => for advection of scalar variables ?
– In the advection phase of ALE/Euler SMG/Q calculations
 Test cases:
– Noh 2D, Sedov 2D/3D, Standing wave in Saltzmann mesh
– New “Noh Corner” Test
In Finite Volume Schemes
 The conservation laws solved over a control
volume :
– (e.g. the computational zone - for zone centered
d                   
 dV    (u  u
variables)
g   )  ds
dtV          V

 Fluxes evaluated on the control volume faces
 Face-centered values of variables are required
Slope Limiters
 Second order Godunov and some other high resolution
schemes use the gradients of the variables to get face-value:

 

 f  c   c  rf  rc 
 

 In hyperbolic PDE discontinuities (shock or slip lines) are
present, or can be generated, during the flow.
 At discontinuities gradients are unbounded.
 Using them produces unphysical fluctuations in the solution
 Thus either gradients or the related fluxes must be limited to
preserve a monotonic flow field
Slope Limiters in 1D
 Limiters were first formulated in 1D:
– Either the slope-extrapolated values of a variable at the zone faces
must lie in the range spanned by the values of that variable in
neighboring cells:
1                             
min( n )  0   lim  r  r0   max(n )
 Or the limited gradient defined as some kind of monotonic
averaging of the gradients in the zone c and its neighbors .
e.g. van Leer’s double minmod:                       
a   ; b  
n

     ab
min(     ,2a,2b) if ab  0
2    ave(a, b)        2
0
               if ab  0
Limiters for scalars in 2D/3D
 In 1D the two approaches are equivalent
 Equation (1) is the usual way to extend the
limiters to 2D/3D
– See Dukowicz, Barth&Jesperson
 Equation (2) is van Leer’s double minmod
– The gradient of scalars are vectors
– What would be monotonic averaging of vectors
in 2D/3D ?
How is (1) extended to 2D/3D
 Look at zone centered scalar variables
 Zone face in 1D=>either face/edge center or zone
vertex in 2D/3D
 Extrapolate    using (1) from the zone center to each of its
vertices
 Find a limited gradient            
          
 lim   
– Such that all vertex extrapolated values of the density
remain within the density range spanned by values in
the zones around that vertex
Extending Equation (1) to 2D/3D
(a) Scalar Limiting    lim   


 For each vertex v find the largest  v such
that:
 
3   min(n )  0  v   r  r0   max(n )
 And take:
4   min v 
v

 Straightforward but dissipative !!
 This was our default scheme for scalar variables
to
Extending Equation (1)  2D/3D
(a) Non scalar α       lim   


Find  lim obeying Eq. (1) at each

           
vertex ʋ “closest” in some sense to the

unlimited   .
 Limiting should be preferentially normal to
shock direction, so that the gradient direction
may change:  lim   
          

 Berger et.al. - Linear programming
 Ridder&Kothe – Least square reconstruction
Let us extend Eq.(2) to 2D/3D !
 We apply the VIP convex hull based monotonicity
criterion to the slopes:
ab
 In 1D:     min(     ,2a,2b) if ab  0
2
 In 2D/3D:
– If a slope is outside the VIP generated by twice
the slopes of its neighbors=>move into the hull
• Factor 2 like in the double minmod
• e.g. we move it to the nearest point on the VIP
What is the analog of ab  0 ?
 In 1D   ab  0 corresponds to an extremum

– (a) In 2D/3D extremum is when       0 .
• In the vector space of the slope-vectors,
extremum exists if the origin (0,0,0) is inside
the VIP => set slope to zero
• Difficulty: “Ridge configuration” a small
slope component in some direction may
prevent detecting a discontinuity normal to a
VIP slope limiter for scalars

(b) At a surface discontinuity   0
n
ˆ
n is the normal to the surface

                      ˆ
•       changes sign along n
n
• Take n    in the direction of
ˆ
VIP slope limiter for scalars
    
 
 Let ai   i ; i  0,1..k the unlimited slopes
in cell c and its k neighbors.
                
 Let n  a0 a0 and si  ai  ni ; i  1..k
ˆ                          ˆ

 with smn  min si ; smx  max si
1i  k   1i  k

 If ( smn smx  0) => zero the slope: a0   0  0
 lim  lim
                               k
 Else if a 0 lies outside the VIP 2ai i 1
(convex hull ) move it to the closest point
on hull
A Convex Hull (CH)
3

1
4
5           2

6           7
9       12

11
8
10
The VIP convex hull limiter
V0    “Move to
hull”

V1

Vc   V5

V2
V0     V6

V4

V3         • if V0 inside VIP => no limiting:
• if V0 outside, VIP limiter is
applied and V0 => Vc
The Staggered Mesh Godunov (SMG/Q) Scheme

 SMG schemes for 3D Lagrangian and ALE hydrodynamics were presented at
Oxford 2005, Prague 2007, Pavia 2009,APS 2003,2005.
 Use Riemann problem (RP) solutions to capture shocks.
 Vertex-centered velocities have jumps on the in-cell corner zone faces.
 These are simplified, “impact” RP or (IRP) with continuous p,     .

   The IRP are solved in the normal to shock direction, i.e. along the

velocity difference u L  u R  . This defines a uniaxial tensor pseudo-
viscosity.
The Staggered Mesh Godunov (SMG/Q) Scheme

 The velocities on either side of the face,           Limited slope

serving as data for the RP, are evaluated
from the vertex velocities using a cell    
u
centered velocity gradient, limited to                              Original
preserve a monotonic velocity                   uL                    slope

distribution.
 A limited velocity gradient is also used                              Velocity
uR        jump for
to compute the momentum fluxes during                                    RP
L     C          R
 The resulting scheme captures shocks

with sharp monotonic profiles. It also has           u
a strong inherent mesh stabilizing effect.   u                
FQ

 The rotation invariant VIP limiter is used                   FQ
to preserve the symmetries present           u        
u
Test Problems
 Saltzman, Sedov, Noh problems test the propagation of 1D
shocks over non-aligned or distorted meshes
 Our corner Noh (a 2D/3D collision Riemann problem)
tests 2D/3D sector-wise converging flows
 VIP for scalars is used only in the advection phase
– ALE could smooth the Saltzman mesh making it regular
– Sedov test will lack resolution at the shock in an Eulerian or
 In the examples VIP for scalars is compared with using the
convex hull limiter only for vectors and “scalar limiting”
(Eqs.1,3,4) for scalar variables
The Standing Shock Saltzman Test
 p, , uL  [133.33,4,10]; p, , uR  [.00667,1,0];  5 / 3 ;US  13.33

 In the original test, a piston moves with a velocity
of u=10 into an almost cold ideal gas. This 1D
problem is solved on a 2D 100x10 mesh as shown.
 Instead, we look at the same mesh fixed in a
frame of reference in which the shock is stationary
with the gas flowing in from the right and out to
the left
 p, , uL  [133.33,4,3.333]; p, , uR  [.00667,1,13.33];  5 / 3 ;US  0
Saltzman with standing wave T=0.4
isodensity plot
(a) VIP vectors & scalars (b) VIP vectors
Saltzman with standing wave T=0.4
Profiles of density and pressure
(a) VIP vectors & scalars (b) VIP vectors
2D Noh Test
 A disk of cold ideal gas   5 / 3 , e  0
collapses to its axis with [ p,  , ur ]  [0,1,1]
 A 100x100x2 disk is meshed with 50x50x1 zones
 The analytic solution is a shock expanding
with U S  1 / 3
16         ; 0  r  t /3
 and density profile:  (r )  1  t / r ; t / 3  r  t

 ALE motion is limited to preserve resolution at
and behind the shock
Noh 2D T60 iso-density plot
(a) VIP vectors & scalars   (b) VIP vectors
Noh 2D T60 iso-density plot inside
(a) VIP vectors & scalars   (b) VIP vectors
“Noh Corner” Test

 A disk of cold   5 / 3 , e  0 ideal gas moves
toward the origin [ p,  , un ]  [0,1,1] with the
velocity directed in each quadrant along the
diagonal as shown
 Due to the symmetry only a quadrant has to be
considered. We carry out this simulation in a
frame in which the gas is initially stationary and 

the walls are moving                          u    u
 In a “hot gas” corner flow e  0.3 and
[ p,  , un ]  [0.2,1,1]
     
u     u
“Cold” Noh Corner Test Iso-densities at T=75.0
100x100x1 zones
(a) VIP vectors and scalars      (b) VIP vectors only

a                                                       b
“Cold” Noh Corner Test Iso-densities at T=75.0
VIP vectors and scalars
(a) 200x200x1 zones     (b) 100x100x1 zones
“Hot” Noh Corner Test Iso-densities at T=70.0
100x100x1 zones
(a) VIP vectors and scalars              (b) VIP vectors only
Sedov-Taylor blast wave problem

 The initial data is:  p, , u0  [0,1,0];   5 / 3, e  0
cold ideal gas with a “point source” e  5027.7 in a
single zone at the origin
 The physical space is a 253 box divided into a
uniform 903 mesh with 3D computation conducted
in the [ x  0, y  0, z  0] octant with the symmetry
planes at x  0, y  0, z  0
 The exact solution has a post-shock density   4
and the front radius is R  1 at T  1.
 To keep the mesh fine near the shock, the ALE
motion was restricted to start at densities   0.25
Sedov-Taylor 3D Test T=1
VIP vectors & scalars isodensity plot
2D Sedov-Taylor test

 The initial data is:  p, , u0  [0,1,0];   5 / 3, e  0
 To remove any effect of boundary conditions we
ran the simulations in all 4 quadrants:
 25  x, y  25
 We set e  5027.7 at the 4 central zones
 The stronger cylindrical shock reaches R  1
at T  0.24
 The ALE motion was again restricted to              0.25
Sedov-Taylor 2D Test T=0.24
isodensity plot
(a) VIP vectors and scalars   (b) VIP vectors only
Sedov-Taylor 2D Test T=0.24
isodensity plot inside
(a) VIP vectors and scalars   (b) VIP vectors only
Conclusions
 We have devised a new VIP based limiter for the gradient
of scalar zone centered variables
 This limiter is applied if either the unlimited gradient lies
outside the VIP spanned by (twice) the gradients in the

neighbor zones or if  n changes sign
 Near a shock the second criterion dominates, so that we do
not have to compute the convex hull and we zero the
limited gradient. Thus for this case the results are identical
to the scalar limiter (Eq. (1),(3),(4) )
More Conclusions

 In principle the new limiter should be less
dissipative than the “scalar limiter”
 While we see some differences in the finer details
pointing to this, in the test cases shown the results
using both limiters are quite similar, virtually
identical.
 On one hand the new limiter performs well in all
the test cases, on the other hand we will consider
more test cases where it may have a clear
Future Perspectives
 This limiter is applied once to the zone
gradient is used for fluxes at all zone faces.
 Other variants, e.g. a directional limiter,
should also be considered
References
 Luttwak G. , Falcovitz J., ”Slope Limiting for vectors: A Novel Vector
Limiting Algorithm”, Int.J.Num.Meth.Fluids,65,p1365-1375,(2011), presented
at the Conf. Num. Meth. for Multi-Material Flows , Pavia, Italy, (2009)
   Luttwak G. , Falcovitz J., “Vector Image Polygon (VIP) limiters in ALE
Hydrodynamics”, EPJ Web of Conferences,10,00020,(2010), presented at the
8th Int. Conf. on New Models and Hydro-codes, Paris, France, (2010)
   M.Berger, M.J.Aftosmis, ”Analysis of slope limiters on Irregular Grids”, NAS
Technical Report NAS-05-007,[2005]
   W. J. Rider, D. B. Kothe, “Constrained Minimization for Monotonic
Reconstruction,” AIAA-97–2036, 1997.
   Dukowicz J. K. and Kodis J. W. “Accurate conservative remapping (rezoning)
for ALE computations, SIAM J. Scient. Stat. Comp., (8),p305-321,1987
   Luttwak G. , Falcovitz J., ”Applying the SMG Scheme to Reactive Flow”, J.
Energetic Materials, 28(1), p279 - 302,(2010), presented at the 7th Int. Conf. on
New Models and Hydro-codes, Lisbon,Portugal,2008
References
 Luttwak, G. ”Staggered Mesh Godunov (SMG) Schemes for
Lagrangian Hydrodynamics”, p339-342, Shock Compression of
Condensed Matter-2005, Furnish M. D. et al Eds, (2006), AIP, CP
 Luttwak G., “ Sliding and Multifluid Velocities in Staggered Mesh
(MMALE) Codes ", Conference/Workshop On Numerical Methods
For Multi-Material Flows, Prague , Sept. 2007, www-
troja.fjfi.cvut.cz/~multimat07
 Luttwak G. , Falcovitz J. , "Staggered Mesh Godunov (SMG)
Schemes for ALE Hydrodynamics", Workshop On Numerical
Methods For Multi-Material Flows , Oxford, UK, Sept. 2005 ,
www.extra.rdg.ac.uk/ifcd/Multimaterial-Workshop.htm
 Luttwak, G., "Comparing Lagrangian Godunov and Pseudo-Viscosity
Schemes for Multi-Dimensional Impact Simulations”, p255-258,
Shock Compression of Condensed Matter-2001, Furnish M. D. et al
Eds, (2002), AIP, CP620.
References
 Christensen R.B., "Godunov Methods on a
Staggered Mesh. An Improved Artificial
Viscosity", L.L.N.L report UCRL-JC-105269,
(1990).
 Caramana E.J., Shaskov M.J., Whalen P.P., J.
Comp. Phys. 144, p70, (1998).
 D.J.Benson,”Computational Methods in
Lagrangian and Eulerian Hydrocodes”,
Comp.Meth.Appl.Mech. Engn.,99 (1992)
,p235-391
How does a slope limiter work?
 Familiar algorithms are formulated for scalar variables.
 The slope extrapolated values of a variable must lie in
the range defined by the values of that variable in
neighboring cells ν:
          
min( )  0   lim  r  r0   max( )

 For vectors, limiting is usually applied separately to
each component.
 Such procedure is frame-dependent:
– Rotating the coordinate axes produces different results
– Component limiters do not preserve problem
symmetry
– They do not transform like vector entities should
The SMG VIP Vector Limiter
 Both for the Lagrangian and Advection phases of the SMG scheme
(in 2D/3D).
 The extension to a general connectivity (FEM) mesh is
straightforward
 The limiting is done separately along each edge of a cell (e.g. 1-2 in
the figure) :
– The outward extrapolated velocities must lie in the VIPs of edge nodes
1,2
– The limited gradient assumed to lie along the velocity difference
4                                                            
6                                     v1, 2  v2  v1 ; v1, 2  v1, 2 / v1, 2
ˆ

7               1
2          3

8
5
SMG VIP vector limiter
 The velocities of 1,2 are extrapolated outward along
the edge.
 e.g. velocity of 2 extrapolated toward middle of 23:
         
r2  0.5(r23  r12 )r12
ˆ ˆ

 The cell centered velocity gradient is used, taking its

component along the velocity difference v1, 2 :
4
                         
v2,e  v2   2 v i  r2   v12 v12
6
ˆ ˆ
7       1
2   3

8
5
SMG VIP limiter

 The value of 0   2  1 chosen to keep v2 ,e inside the VIP of
the edge neighbors (1,5,3,4)      
 The same way  1 chosen to keep v1,e inside the VIP of the edge
neighbors (2,6,7,8)
 The VIP limited slope along 12 generates a velocity jump
between the corner zones :
v   v

12

 v12
12
 red
                         
where: v red  min  , v  (v  r ) v
12        1 2 ˆ12          12
ˆ12
   
 serves as the data for the IRP in the
 The above jump v
Lagrange phase 12
4
6

7             1
2     3

8
5

 Momentum integrated around corner zones
around vertex O (12345678 in the figure)
d                  

 udV  V u (u  ug )  ds
dt Vst         st

2
3                     1

O
4                 8

5
6   7
 Momentum fluxes at the mid-faces between the
corner zones found using the limited velocity
profile
 The vertex mass flux between these corner zones
found from the mass flux through the zone faces
(similar to Benson HIS)
 The upstream weighted velocity at the middle of
the fluxed volume is:         mda   red
       
v f  vd  0.51     vda
    md 
here d,a are the donor and acceptor sides (1or 2)

 And the momentum flux => mda v f
The Staggered Mesh Godunov (SMG/Q) Scheme

 Let p* be the RP solution pressure. We take q  p *  pcell
as a uni-axial tensor pseudo-viscosity acting along the
shock direction.
 Its impulse is imparted to the two neighboring vertices.
 Its work, which must be dissipative, is added to the zone
internal energy.
– Some of these ideas are related to Christensen’s split-Q, and to the
edge viscosity and compatible hydro-scheme by Caramana et al.
When is a vector “monotonic” relative to its
neighbors?
 We seek a monotonicity preservation criterion for a
vector, analogous to the scalar criterion:
min      max (  )
               

 We have defined VIP - Vector Image Polygon in 2D
(Vector Image Polyhedron in 3D) as the convex hull of
the vector-space points corresponding to the neighbor
vectors.
 If a slope-extrapolated vector lies inside the VIP, the
slope is monotonicity preserving. Otherwise, slope
limiting is required.
 If the original vector lies entirely outside the VIP, we zero
its slopes, in analogy to the scalar extremum case
The VIP Monotonicity Criterion
 A point inside a convex hull can be nexpressed as:

n

v    i vi ;0   i  1;  i  1
i 1                 i 1
 In 1D (n=2) this reduces to:
                
v  v1  (1   )v2 ;0    1
which is equivalent to the scalar criterion (1)
 A vector point lying on the line connecting two
neighbor points should be monotonicity complying.
Since VIP is convex any segment connecting two
vertices lies inside the VIP (or on boundary).                 1
  The VIP monotonicity criterion is a natural
extension of (scalar) monotonicity to vectors.
A
2
3
Falcovitz: Cold Noh Corner Test
1000x1000 Eulerian GRP calculation

```
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
 views: 18 posted: 2/13/2012 language: Latin pages: 49
How are you planning on using Docstoc?