Sliding and multi-fluid velocities in Staggered Mesh Eulerian and MMALE codes

1. Sliding and multi-fluid velocities in Staggered Mesh Eulerian and MMALE codes Gabi Luttwak1 1Rafael, P.O. Box 2250, Haifa 31021, Israel

2. The Velocity in Eulerian and MMALE Codes • Most Eulerian and MMALE codes assume a common velocity in multimaterial cells • Is this right and what are the consequences? • Let explore the price and the benefits of defining separate velocities for each fluid • MMALE (Luttwak & Rabie 1985,Luttwak1989) , had multifluid velocities. • We review its scheme. • We introduce the multifluid velocities into The SMG scheme (Luttwak & Falcovitz 2005) • Show some results • Discuss the advantages price & limitations and look ahead for further improvements

3. Multi-fluid velocities • What happens at an interface? • In 1D both pressure (or the normal stress) and velocity are continuous • In 2D (3D) slip lines (planes) form at density discontinuities: • Even in a single fluid • And almost always at material interfaces • In Lagrangian codes • Sliding-Impacting techniques are used to handle such situations

4. Sliding and Impact in Lagrange and SMALE • In Lagrangian and almost-Lagrangian ALE or SMALE codes a sliding/impacting contact algorithm enables the simulation of sliding, impact and separation • Without it: • impact/separation cannot be described • Due to sliding, a row of cells near the interfaces undergoes “non-physical” shear. This leads to a mesh size dependent and un-physical friction. If the relative sliding is large, this mesh deformation stops the calculation.

5. In Eulerian and MMALE codes: • The fluids flow through the cells • The mesh does not deform due to sliding/impact • The calculation can proceed • But, with a common interface velocity: • Sliding introduces an “effective” friction and viscosity • Therefore the strain field will be wrong near the interface • in elastic-plastic flow this affects failure • The rate of growth of instabilities (like R-T,R-M) might be slower.

6. MM-Euler and MM-ALE codes using multi-fluid velocities • Most MM-Euler and MM-ALE codes use a common vertex velocities, for all species. Remarkable exceptions are: • Luttwak & Rabie-1985 One of the first MMALE codes … • Walker & Anderson 1994 noticed the importance of this issue for elastic plastic flow and introduced cell centered multifluid velocities into the MM-Eulerian code CTH. • Vitali & Benson 2006 – use Belytschko X-FEM method to locally add multifluid velocity using an extended FEM method for staggered mesh MMALE formulations. • Our purpose is to incorporate multi-fluid velocities into the SMG scheme for ALE hydrodynamics (Luttwak & Falcovitz 2005)

7. The “old” 1985 2D MMALE code • Like the SALE code, it used a staggered mesh in space and time: • pressure, stress, density, internal energy and Q are cell centered • The velocity is defined at the vertices at Tn+1/2 • These variables (except the Q) were also defined separately for each fluid. Unlike to SALE and to most other codes, each fluid had a separate vertex velocity and vertex mass. • Both the cell centered and vertex material variables used a LIFO stack logic so that memory was used only at cells were a fluid is present • Each fluid was accelerated separately: • A common normal to interface velocity was enforced by momentum conservation at the vertex.

8. 2 1 2 1 3 1 The (1985) 2D MMALE code • The normal to the interface and the wet volume in cell uniquely defines the wet polygon for 2 materials. • For more than 2 materials the result depend on the order: • The wet polygon of the first fluid is found as before. • The wet volume for the next fluid defined by its normal and wet volume inside the polygon not filled by the previous fluids.

9. 4 5 3 1 2 MMALE (1984) – The Advection Phase • The advection phase was based on a remap: • The vertices of the wet polygon in cell c filled with a fluid were moved with the fluid velocity there . • The mesh was advanced with the grid velocity. • The overlap volume was found and it is the flux of fluid i from cell c to a nearby cell n. (Luttwak 1994,Luttwak&Cowler 1999) • The densities in the flux volume were obtained by linear interpolation between donor and acceptor cell values with the flux values limited to prevent the generation of new extreme values. n c

10. Vertex Variables Advection in the “old” MMALE-2D code • Vertex variable update used SALE’s method A cell centered mass and momentum was first advected • For single material cells, the influx of mass and momentum to a cell were equally distributed to its vertices • For multi-material cells, the fluxes were split according to the relative wet volume fraction of each fluid in the corner zone. • The method was excessively diffusive and the distribution of the mass to the vertices for multi-material vertices caused problems

11. The Staggered Mesh Godunov (SMG/Q) Scheme • SMG schemes for 3D ALE hydrodynamics were presented at Oxford 2005,APS 2003,2005. • They use Riemann problem (RP) solutions instead of pseudo-viscosity to capture shocks. • The vertex defined velocities may have jumps mid-way, on the faces separating the corner zones, near the cell center. • These are simpler, “collision” RP. (with constant there). • They are solved in the normal to shock direction, i.e. along the velocity difference .

12. The Staggered Mesh Godunov (SMG/Q) Scheme • Let p* be the RP solution pressure. We look at as a uni-axial tensor pseudo-viscosity acting on the corner zone internal face along the shock direction • Its impulse is exchanged between the two neighboring vertices. • Its work, which must be dissipative and is added to the zone internal energy. • The above procedure is similar to Caramana et al. compatible hydro-scheme. The SMG scheme has some similarities with Christensen-1990 split-Q and Caramana 1997 et al. edge viscosity.

13. Limited slope Original slope Velocity jump for RP L C R The Staggered Mesh Godunov (SMG/Q) Scheme • The velocities on both sides of the face, serving as data for the RP are evaluated from the vertex values using a, cell centered velocity gradient limited, to preserve a smooth velocity distribution. This serves as a shock detector. • The resulting scheme captures sharp and smooth shocks, with an additional strong and natural mesh stabilizing effect.

14. Multi-material velocities • The fluids have distinct vertex defined velocities and masses • On a material interface a common normal to interface velocity component is enforced. • When appropriate (fluid flowing away from the interface, negative pressure or normal stress, large void fraction) these normal velocities are left unchanged enabling interface separation • The fluid vertex velocities are only stored for vertices sharing at least one cell containing that fluid.

15. The Set-up of the SMG-MMALE code • Geometry Update (not necessary for pure Euler) • compute face areas, volumes and volume changes. • Lagrangian Phase • The SMG Scheme • Definition of pressure and stress in a MM zone • Advection Phase (not present for pure Lagrange) • Second order advection for cell-centered variables • Staggered Mesh Momentum Advection • HIS (Benson) • Interface tracking: VOF (Luttwak 1984,2002) similar to Youngs scheme • Grid Motion Algorithm • Almost all parts need some change….

16. Lagrangian Phase • Strain, velocity gradient and volume changes are specific for each fluid in the cell: • To find the velocity is integrated over the full cell. • If a fluid has a very small vertex mass, an “extrapolated” velocity is used obtained from the fluid velocity average in the cell.This prevents unphysical large gradients when the wet volume is small. • A common vertex acceleration is computed integrating the pressure (without Q) over the (dual) cell using Wilkins scheme. The specific fluid accelerations are found using:

17. SMG • The “collision” RP is solved separately for each fluid • The RP data uses the fluid velocities and densities • The forces act only on the particular fluid momentum • The associated work is added to that fluid internal energy

18. Interface Tracking • For each fluid we find thenormal to the interface: • The normal is found by integrating the relative partial volumes over the cell faces. • The normal and the wet volume Videtermine the interface position and the wet part of the cell. • Unlike to the 1985 MMALE 2D , we are ignoring possible overlap when several materials present in a cell. • The normal at a vertex is found by integrating over the faces passing through the vertex (like for acceleration)

19. R L The Advection Phase • The volume flux between two face neighbor cells is defined by the “wet” part of the volume enclosed by the common face old and new positions. • The new face is obtained by moving its vertices with the fluid velocity relative to the grid. This is different for each fluid. • The “wet” part is defined by cutting it with the interface plane found by the tracking algorithm advanced to Tn+1. • The remap algorithm used in the older 2D code is superior as it took into account diagonal advection but is more expensive

20. 8 7 5 6 4 3 1 2 Vertex Variable Advection • Similar to Benson HIS • Treating I,J,K pair of faces separately, let be the mass influx from left and right into a cell, is advected from the left face vertices (1,2,6,5 in the picture) to the corresponding vertices on the right. Each vertex takes an equal share. (A quarter for the hexa cell shown). Momentum advection is done separately for each pair of vertices, the velocity of the fluxed mass is evaluated using the donor vertex velocity and the cell limited velocity gradient. • We keep the same scheme for the fluid fluxes. Thus each vertex shares an equal part of the cell fluid mass. However, until the vertex is wet (it is in the cell wet part) its velocity will be modified only by advection.

21. LIFO stack memory set-up • A structure (type in Fortran 90 ) is defined with all the fluid specific vertex variables. • An array of this structure is allocated. The unused elements are connected in a linked list. When a fluid leaves all the cells surrounding the vertex its element is freed and put on the top of the stack to be first used. • Each vertex points to an element of the array storing the data for the first fluid. The element itself points to the element containing the next fluid data or to a Null element (All empty vertices also point to the Null element) . • When there are no more free elements, another array is allocated with the elements linked in a stack. • A similar logic is used to store the vertex variable fluxes. This is done in a separate structure as it must exist whenever there is a fluid in any vertex around.

22. Preliminary Tests • Some results showing the capabilities of the method are presented. • 2D problems run in the 3D code with one layer of cells and plane symmetry. • The first is a trivial frictionless sliding of a copper block on steel. It exhibits some difficulties, which were mainly overcome as shown in the following picture • The second are two test cases of normal plate penetration: • 2m/ms impact of copper on steel • 2m/ms impact of aluminum on steel

23. T=0 T=0 T=9.6 2m/ms Cu Cu 2mm/ Cu Cu Steel Steel Steel Steel A (not so) trivial test

24. A (not so) trivial test of sliding • The interface tracking produces some corner rounding. Taking a common normal velocity component there, starts an un-physical interaction. Checking if the moment of the impact is within the time step helps here, but makes the normal impact less smooth. The local interface location is found from the tracking algorithm. When this check is restricted to cells including more than 2 materials (including void) both problems run well. • A better interface normal definition for corner cells sharing more than 2 materials will help. Perhaps the interface direction should be determined by the neighboring cells which share only 2 materials

25. 2m/ms Normal Impact of Copper or Aluminum on Steel • The impact was run both with the standard common vertex velocity and the new multifluid vertex velocity versions. • A Steinberg EOS was used with the material parameters for Copper (Cu-OFHC), Steel (SS-21-6-9) aluminum(AL-2024-T4), taken from the AUTODYN library. • The 3D code was again run in 2D mode using one layer of cells

26. 2m/ms Normal Impact of Copper on Steel T=0 T=5

27. 2m/ms Normal Impact of Copper on Steel T=15

28. Comparing T=5, 2mm/ms impact Cu on steel a) Multifluid velocities b) Common vertex velocities

29. Comparing T=7.5, 2mm/ms impact Cu on steel a) Multifluid velocities b) Common vertex velocities

30. Comparing T=10, 2mm/ms impact Cu on steel a) Multifluid velocities b) Common vertex velocities

31. Comparing T=15, 2mm/ms impact Cu on steel a) Multifluid velocities b) Common vertex velocities

32. Comparing T=5, 2mm/ms impact Al on steel

33. Comparing T=7.5, 2mm/ms impact Al on steel a) Multifluid velocities b) Common vertex velocities

34. Comparing T=10, 2mm/ms impact Al on steel a) Multifluid velocities b) Common vertex velocities

35. Comparing T=12.5, 2mm/ms impact Al on steel a) Multifluid velocities b) Common vertex velocities

36. The results of the normal impact tests • The early time results are similar. • At later stages the multi-fluid velocities allows separation of the interface which is a feature we look for. • However, the multifluid-velocity version seems to produce too much fragmentation; with the common vertex velocity we get smoother results.

37. Concluding Remarks • We have set-up a multi-fluid velocities version of the SMG scheme. • The present multi-fluid velocity version of the code is not yet stable. Further improvements and more debugging is needed • A better normal definition for cells sharing more than 2 materials • The re-use of some of the features of the “old” MMALE 2D code should help: • Remap type advection near the interfaces • Multi-fluid velocities better describe the physics of the interface dynamics: • at low speed, where strength and failure are important • at higher speed where the growth of instabilities might be relevant

38. References • G. Luttwak, R.L. Rabie, ”The Multimaterial Arbitrary Lagrangian Eulerian code MMALE and its application to some problems of penetration and impact”,LA-UR-85-2311,(1985) • Luttwak G., "Numerical Simulation of Jet Formation Using an Adaptive Grid in the Multi-Material ALE Code MMALE", p45-52 in Proc. of the 3d HDP (High dynamic Pressure) Symp., La Grande Motte, France, (1989) • J.D. Walker, C.E. Anderson,” Multimaterial Velocities for Mixed Cells”, p 1773-1776,High Pressure Science and Technology-1993,ed. Schmidt et al.,AIP,(1994) • G. Luttwak, J. Falcovitz, ”Staggered Mesh Godunov (SMG) Schemes for ALE Hydrodynamics”, presented at the “Numerical Methods for Multi-material Flows”, Oxford, Sept. (2005) • E. Vitali, D. Benson, Int. J. Numer. Meth. Engng. 67,p142-1444,(2006)

39. More 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). • Wilkins M. L., "Calculation of Elastic-Plastic Flow”, Meth. Comp. Phys., 3, p211, B. Alder et al. eds., Academic Press (1964). • 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. • Luttwak, G. ”Staggered Mesh Godunov (SMG) Schemes for Lagrangian Hydrodynamics”, presented at the APS Topical Conference on Shock Compression of Condensed matter,Baltimore,MD,July31-Aug5,2005

40. More References • Luttwak G. "Second Order Discrete Rezoning", in High-Pressure Science and Technology -1993, Schmidt S.C. et al (editors) , pp1777-1780 AIP-Press (1994) • Luttwak G., Cowler M. S., "Remapping Techniques for Three Dimensional Meshes", presented at the 1999 Int. Workshop for New Models and Hydro-codes for Shock Wave Processes, held at the Univ. of Maryland, College Park, MD,(1999) • Luttwak G., "Interface Tracking in Eulerian and MMALE Calculations",p283-286, Shock Compression of Condensed Matter-2001,Furnish M. D., et al. eds., AIP, CP620, (2002) • Luttwak G., Cowler M. S., Birnbaum N. , "Virtual Memory Techniques in Eulerian Calculations", p363-366, Shock Compression of Condensed Matter-1999, Furnish M. D., et al. eds.,AIP,CP505,(2000)