Transcript The Application of a General Higher

LA-UR 11-03531
The Implementation of a General Higher-Order Remap
Algorithm
Vincent P. Chiravalle
MultiMat 2011 Conference
Arcachon, France
Slide 1
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
Outline
•
Cercion code structure
–
–
–
–
•
Linked list data structures
Lagrangian solution
Calculation of stress tensor
Remap algorithm
Test Problems
– Riemann shock tube
– Aluminum flyer plate
– Imploding cylindrical shell
•
Conclusions
Slide 2
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
The hydrodynamics is solved on a block structured mesh
•
Block structured mesh
–
–
–
Four-sided cells
Fixed connectivity among mesh blocks
Velocities at the cell vertices
A simple mesh with three structured
blocks (2 zones each)
Slide 3
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
Cercion makes extensive use of cell oriented data structures
•
Cell oriented data structures
– C implementation
– Fortran-type array indexing
•
Cell_t data structure
–
–
–
–
–
Storage for vertex quantities such as position and velocity components
Cell-centered quantities such as density, pressure and volume
Pointer to a linked list of materials in the cell (mat_t)
Pointer to a linked list of material fluxes for the top boundary (flux_t)
Pointer to a linked list of material fluxes for the right boundary (flux_t)
Slide 4
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
Each cell has a linked list to store information about all the materials in
the cell
Material Linked List for a Cell Containing 1 Material
cell_t
header
mat_t
bas_t
•
•
•
tail
str_t
A single mat_t element for each material in the cell
Ordering of the elements according to onion skin method
Two pointers for each element
–
–
Basic material properties in bas_t data structure
• volume fraction, mass and energy
• density, energy density and associated derivatives
• parameters for interface reconstruction
Strength properties in str_t data structure
• stress deviator components and derivatives
• equivalent plastic strain
Slide 5
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
Each cell has two additional linked lists to store material fluxes through
the top and right boundaries
Material Flux Linked List for a Cell Boundary
with 1 Material Crossing the Boundary
cell_t
header
flux_t
tail
sflux_t
•
A flux_t element for each material crossing the boundary
– volume, mass and energy fluxes
– pointer to sflux_t data structure for stress energy and strain fluxes
Slide 6
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
A three phase approach is used to solve the hydrodynamics equations
•
Lagrangian Phase
–
–
–
–
•
Grid Relaxation Phase
–
•
Uses the algorithm from HEMP (Wilkins 1963)
Spatially staggered grid with vertex velocities
Material strength included
• four stress deviator components
• yield stress correction
• Margolin method for calculating strain rates and flow divergence
Margolin anti-hourglass treatment
Simple finite difference mesh relaxer
• nine point stencil
Remap Phase
–
–
–
donor cell procedure from the SALE code (Amsden et al. 1980) for cell-centered quantities
second order correction to the remapped fluxes from the donor cell method
• derivatives using Barth-Jespersen slope-limited method
Mass fluxes from the density remap for the momentum remap
Slide 7
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
The Lagrangian equations of motion are solved using the HEMP approach
Vertex (i,j) and surrounding
cells A,B,C, and D
Vertex quantities for the equations of motion
1
AA  AB  AC  AD 
4
   xy A    xy A    xy A    xy A  
 
 
 
 
 
2  M  A  M  B  M C  M  D 
   yy      yy      yy      yy    
 
A   
A   
A   
A  
2 
M
M
M
M
A 
B 
C 
 D 
i , j 
 i, j
i, j
Axial Velocity Equation
t
A
 y4  y1    Bxx  y1  y2   Cxx  y2  y3    Dxx  y3  y4 
uin, j 1  uin, j 
 xx
2i , j





t
A
x4  x1    Bxy x1  x2   Cxy x2  x3    Dxy x3  x4   t i , j
 xy
2i , j
Radial Velocity Equation
vin, j 1  vin, j 



t
 Ayy x4  x1    Byy x1  x2   Cyy x2  x3    Dyy x3  x4 
2i , j


t
A
 y4  y1    Bxy  y1  y2   Cxy  y2  y3    Dxy  y3  y4   ti, j
 xy
2i , j
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
Slide 8
Strain rates are calculated using the Margolin method
Cell (i,j) and its four
neighboring vertices
Flow divergence
 U   xx   yy  
Equations for the strain rate components
1
u3r234  u1r412  y2  y4   u4 r341  u2 r123 y3  y1 
xx 
6V
1
v1r412  v3r234 x2  x4   v2 r123  v4 r341x3  x1   v1 A412  v2 A123  v3 A234  v4 A341
 yy 
6V
1
v1 A412  v2 A123  v3 A234  v4 A341
 
2V
1
v3r234  v1r412  y2  y4   v4 r341  v2 r123 y3  y1 
xy 
6V
1
u1r412  u3r234 x2  x4   u2 r123  u4 r341x3  x1   u1 A412  u2 A123  u3 A234  u4 A341

6V
V 
1
r123A123  r234A234  r341A341  r412A412
6


 

1
yk  y j xi  x j  yi  y j xk  x j
2
rijk  yi  y j  yk
Aijk 

Slide 9
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
Material fluxes are used to redistribute materials among cells during
the remap phase
Donor Cell Method
k
MFzik, j  upwind
Fzik, j
•
•
•
Mesh relaxation moves (xp,yp) to (x,y)
Total volume fluxes (Fr, Fz) for donor cell method
Mass fluxes (MFr,MFz) for each material leaving the cell
based on material volume fluxes and upwind densities
Update of material k mass from Lagrangian value, M, to remapped value M*

M ik, *j  M ik, j  MFzik, j  MFzik1, j  MFrik, j  MFrik, j 1

Slide 10
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
A higher order correction to the donor cell method is used for cells with
a single material
MFri , j
 


 upwindFri , j  Fri , j dx
 dy


x

y
upwind

upwind 

donor cell method
correction
Slide 11
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
Riemann shock tube problem was used to test the code
•
•
An initial discontinuity between two ideal gas regions (g=1.4)
Region 1
– density of 1 kg/m3
– pressure of 105 N/m2
•
Region 2
– density of 0.01 kg/m3
– pressure of 103 N/m2
•
•
•
A uniform box mesh with 500 (axial) by 2 (radial) zones
Eulerian calculation
The solution is obtained at 0.01s
Slide 12
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
The locations of the shock and contact discontinuity are captured
(a) Density
(b) Pressure
Slide 13
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
The velocity and sound speed spatial profiles are reasonably well
represented
(a) Velocity
(b) Sound Speed
Slide 14
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
Cercion calculates the density as well as other numerical
methods for the Riemann problem
•
•
•
•
•
•
An initial discontinuity
between two ideal gas
regions (g=1.4)
Region 1
– density of 1 kg/m3
– pressure of 1 N/m2
Region 2
– density of 0.125
kg/m3
– pressure of 0.1 N/m2
A uniform box mesh with
100 axial zones
Eulerian calculation
The solution is obtained at
0.14s
Slide 15
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
Cercion compares well with other numerical methods in
capturing the velocity profile
Slide 16
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
Flyer plate problem tests the material strength routines in the code
•
•
Cylindrical geometry with 1 radial zone
Aluminum target
–
•
Aluminum projectile
–
•
r0=2.707
C0=0.5386
S1=1.339
g0=1.97
b=0.48
Material strength model
–
–
•
•
0.2 cm thick with 40 axial zones
Gruneisen EOS
–
–
–
–
–
•
1 cm thick with 200 axial zones
yield strength of 0.0004 Mbar
shear modulus of 0.271 Mbar
Lagrangian calculation
Companion FLAG calculation
Slide 17
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
The temporal velocity profile at the target-vacuum interface from the
code and FLAG agree fairly well
Slide 18
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
An imploding steel shell is a good test of energy conservation in
converging cylindrical geometry
•
Two steel shells each with 10 radial
zones
–
–
•
•
•
5.0cm
8.0cm
r0=7.9
C0=0.457
S1=1.49
g0=1.93
b=0.5
P=0.588 Mbar
=1.84 g/cc
2.0cm
Material strength model
–
–
•
Outer shell is 0.25 cm thick
Inner shell is 0.5 cm thick
Gruneisen EOS for steel
–
–
–
–
–
•
Initial Geometry
yield strength of 0.05 Mbar
shear modulus of 0.895 Mbar
Ideal gas EOS for the pressure driver
material (80 zones)
Companion FLAG calculation with
same mesh
5 ms run time
steel
P=0.0 Mbar
=0.001 g/cc
r
z
Slide 19
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
Cercion and FLAG give similar kinetic energies for the inner shell
when material strength is not included
Kinetic Energy
Internal Energy
Slide 20
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
When material strength is included the agreement between Cercion
and FLAG does not change
Kinetic Energy
Internal Energy
The Cercion calculation is insensitive
to the use of ALE
Slide 21
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
FLAG and Cercion predict virtually the same total energy for the inner
steel shell
Without Strength
With Strength
Slide 22
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
Conclusions
•
•
•
•
•
•
•
Cell-centered data structures simplify the programming of Cercion and allow for
efficient memory allocation for multi-material problems
Cercion uses proven numerical methods for the solution of the Langrangian
equations of motion and the calculation of material strength properties
A second-order accurate remap method is implemented for density, energy and
momentum, enabling ALE calculations to be performed
Cercion shows excellent agreement with the analytic solution for the Riemann shock
tube problem
In pure Langrangian mode both Cercion and FLAG give similar velocity profiles at the
target-vacuum interface for the flyer plate test problem
The cylindrical implosion test problem illustrates that the Cercion solution is relatively
insensitive to the use of ALE in the calculation
Cercion calculations with and without strength for the cylindrical implosion problem
generally agree with the corresponding FLAG results
Slide 23
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
The Sedov blast wave problem tests the symmetry of the code
solution in two dimensional cylindrical geometry
•
•
An ideal gas (g=1.4)
A uniform box mesh (6 cm by 3 cm)
– 300 axial zones
– 150 radial zones
•
Energy source at the origin
– 85 kJ
– 4 zones (2 axial by 2 radial)
•
•
•
Eulerian calculation
The solution is obtained at 10.0 ms
Comparison RAGE calculation
Slide 24
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA
The code and RAGE give different pressures behind the blast wave
and both have a lower density at the front than the self-similar solution
Pressure
Density
Slide 25
Operated by Los Alamos National Security, LLC for the U.S. Department of Energy’s NNSA