Acta Polytechnica https://doi.org/10.14311/AP.2021.61.0068 Acta Polytechnica 61(SI):68–76, 2021 © 2021 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague CELL-CENTERED LAGRANGIAN LAX-WENDROFF HLL HYBRID SCHEME ON UNSTRUCTURED MESHES David Fridricha, Richard Liskaa,∗, Ivan Taranta, Pavel Váchala, Burton Wendroffb a Czech Technical University in Prague, Faculty of Nuclear Sciences and Physical Engineering, Břehová 7, 115 19 Prague, Czech Republic b Retired Fellow, Los Alamos National Laboratory, Los Alamos, NM, USA ∗ corresponding author: liska@siduri.fjfi.cvut.cz Abstract. We have recently introduced a new cell-centered Lax-Wendroff HLL hybrid scheme for Lagrangian hydrodynamics [Fridrich et al. J. Comp. Phys. 326 (2016) 878-892] with results presented only on logical rectangular quadrilateral meshes. In this study we present an improved version on unstructured meshes, including uniform triangular and hexagonal meshes and non-uniform triangular and polygonal meshes. The performance of the scheme is verified on Noh and Sedov problems and its second-order convergence is verified on a smooth expansion test. Finally the choice of the scalar parameter controlling the amount of added artificial dissipation is studied. Keywords: Lagrangian hydrodynamics, Lax-Wendroff, HLL. 1. Introduction Lagrangian hydrodynamical methods are crucial for simulations of high speed compressible fluid flows as e.g. in astrophysics or inertial confinement fusion (ICF). Lagrangian hydrodynamics solves Euler equa- tions for compressible fluid flow on a computational mesh moving with the fluid under the assumption of no mass flux between the cells. Staggered Lagrangian methods, e.g. [1–3], using an artificial viscosity (orig- inally introduced by von Neumann and Richtmyer [4]) as a dissipative mechanism, have been standard for many years. Staggered methods define thermody- namical quantities (density, pressure, internal energy) in the cells, while velocity is defined at mesh nodes. A class of high order staggered finite element methods on meshes with curvilinear cells has been developed in [5]. Another form of the dissipative mechanism, which comes from an approximate Riemann solver, has been developed in the cell-centered methods [6–12], which assign both velocity and thermodynamical quantities to the cells. Our recently introduced Lax-Wendroff HLL hybrid scheme [13] is also a cell-centered method. The artifi- cial dissipation, in the form of dissipative part of the HLL approximate Riemann solver flux [14], is added to the momentum and energy equations. The original method is improved by applying another weighting in the predictor. The original paper [13] has pre- sented numerical results only on logical rectangular, quadrilateral meshes. Here we present results on un- structured meshes, including uniform triangular and hexagonal meshes and non-uniform triangular and polygonal meshes. 2. Lagrangian Finite Volume In Lagrangian coordinates we define the finite volume method on a moving finite volume V with a surface S. The mass mV , momentum Mv and total energy EV of the moving volume V are given by mV = ∫ V ρ dV, MV = ∫ V ρu dV, EV = ∫ V ρE dV, (1) where ρ is density, u velocity and E total specific energy (i.e., total energy per unit mass). From the Reynolds transport theorem one can derive the finite volume formulation of Lagrangian compressible gas dynamics equations for conserved quantities WV = (V,MV , EV ) dWV d t = ∫ S F · n dS, (2) where n is the outward normal to the surface S and the flux F is given by F = (u,−pI,−pu). (3) Here I is the unit matrix and p is pressure, for which we employ the ideal gas equation of state p = (γ − 1)ρ(E−u2/2). The system (2) is solved for all cells V of the moving unstructured computational mesh. For all nodes of the moving mesh one solves the ordinary differential equations dX d t = u = (u, v), (4) where X = (x, y) are coordinates of a node of the mesh and u is the velocity of the node. 68 https://doi.org/10.14311/AP.2021.61.0068 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en vol. 61 Special Issue/2021 Lagrangian Lax-Wendroff HLL Scheme The specific quantities wV = WV /mV (i.e. quanti- ties per unit mass) ηV = V mV = 1 ρV , uV = MV mV , EV = EV mV (5) are the specific volume ηV , the velocity uV and the specific total energy EV . The specific quantities wV = (ηV ,uV , EV ) will be our variables in the Lagrangian schemes. 3. Meshes For the computational mesh we follow the p-c notation [15]. The unstructured computational mesh contains cells c, points p and edges e. Each cell is a polygon. A cell with n vertices can be split into n quadrilateral subcells Spc by separators, which are the segments connecting the cell center to the edge midpoints. The edge e(ca) separates cells c and a. The outer normal to the edge e(ca) (w.r.t. the cell c) is denoted by ne(ca) and its length is equal to the length of the edge e(ca). Around each node p, we construct a dual cell as the union of the subcells Spc of all cells attached to vertex p. The edges of the dual cell p which connect the cell center and the edge centers are the separators spc±. The outer normals nspc± to the separators spc± have the same length as the separators spc±, see Fig. 1. The corner vector npc = ne(ca) + ne(cb) 2 = −(nspc+ + nspc−) (6) is the average of two neighboring edge normals ne(ca) and ne(cb) (see Fig. 2) and is equal to minus the sum of separator normals nspc+ and nspc− of the corresponding subcell Spc. We define the cell center Xc = (xc, yc) as the arith- metic average of the cell nodes positions Xp. Each subcell Spc is a quadrilateral with volume (area) Vpc = 1 2 4∑ i=1 (yi+1 − yi)(xi+1 + xi), where the subscript i denotes four nodes of subcell Spc and X5 ≡ X1. The nodal volume Vp and the cell volume Vc are given by sums of subcells having node p or cell center c as their vertex, respectively: Vp = ∑ c(p) Vpc, Vc = ∑ p(c) Vpc. The subzonal mass mpc is defined by mpc = ρcVpc. The mass of primary cell mc and the mass of dual cell mp are defined similarly as volumes: mc = ∑ p(c) mpc, mp = ∑ c(p) mpc. The Lagrangian approach assumes no mass flux be- tween subcells. Thus the subzonal masses, as well as cell and nodal masses, remain constant. 4. Schemes We employ Richtmyer’s finite volume formulation [16] of the Lax-Wendroff method. Each step of the method consists of a predictor and a corrector. 4.1. Predictor The predictor computes the nodal estimates of con- served quantities at time level n + 1 2 from cellular quantities at time level n wn+ 1 2 p = ∑ c(p) wn c /Vpc∑ c(p) 1/Vpc + ∆t 2mp ∑ c(p) Fnc ·(nnspc++nnspc−), (7) where the summation goes over c(p), which is the set of all cells sharing the vertex p. In the original paper [13] the conservative cell quantities in the first term on the right hand side of (7) were weighted by the Cartesian subzonal masses mpc, giving the standard Lax-Friedrichs (LF) approximation. Here, being inspired by Wendroff and White [17], we weight them by the inverses of the subzonal volumes 1/Vpc. The standard LF weighting on a 1D non-uniform Eulerian mesh looks like (∆xjwj + ∆xj+1wj+1)/(∆xj + ∆xj+1), while Wendroff and White [17] proposed to use (∆xj+1wj + ∆xjwj+1)/(∆xj + ∆xj+1), which is the same as the linear interpolation from cell centers to the nodes. We call this weighting Wendroff- White (WW) weighting. On a 2D orthogonal mesh the WW weighting can be generalized to the weighting by the inverses of subzonal volumes, which is the same as the bi-linear interpolation from cell centers to nodes. In (7) we apply the same weighting to general Lagrangian meshes. 4.2. Corrector The nodal estimates from the predictor (7) are used to compute the final LW cell values at time level n+ 1 mcwn+1 c = mcwn c + ∆t ∑ e(ca) Fn+ 1 2 e(ca) · n n+ 1 2 e(ca), (8) where the summation goes over e(c) which is the set of all edges of cell c, the edge e(ca) is the edge between cells c and a, and where Fn+ 1 2 e(ca) is the average of the two nodal fluxes at the endpoints of the edge e(ca) Fn+ 1 2 e(ca) = Fn+ 1 2 pe(ca)+ + Fn+ 1 2 pe(ca)− 2 . The corrector (8) can be rewritten as mcwn+1 c = mcwn c + ∆t ∑ p(c) Fn+ 1 2 p · nn+ 1 2 pc (9) using the summation over p(c), which is the set of all points, i.e. vertices, of cell c and the corner vectors nn+ 1 2 pc . 69 D. Fridrich, R. Liska, I. Tarant et al. Acta Polytechnica Figure 1. The primary (here triangular) mesh drawn by solid segments with an edge normal to the primary edge eca and one dual cell delineated by the separators (dotted segments), shown along with their normals. Figure 2. Averaging of the edge normals ne(ca) and ne(cb) of the cell c to define the corner vector npc. 70 vol. 61 Special Issue/2021 Lagrangian Lax-Wendroff HLL Scheme 4.3. Mesh positions update The node positions are updated using the nodal ve- locity estimates from the predictor: Xn+1 p = Xn p + ∆tun+ 1 2 p , Xn+ 1 2 p = (Xn+1 p + Xn p )/2 4.4. Artificial dissipation The dissipative terms Dc are added to the corrector (8) mc wn+1 c −wn c ∆t = ∑ e(ca) Fn+ 1 2 e(ca) · n n+ 1 2 e(ca) + Dc, so that the scheme can compute shock waves. We use the dissipative parts of the HLL approximate Riemann solver fluxes [14] as the artificial dissipation terms Dc = Dτ · ∑ e(c) σcaσac σca + σac |nne(ca)|(wn a −wn c ), (10) where the modified signal velocity is defined by σca = ρnc ( cs,nc + |(unc − una) · nne(ca)| |nne(ca)| ) , and Dτ =  τη 0 0 0 0 τu 0 0 0 0 τu 0 0 0 0 τE  (11) is the diagonal matrix of dimensionless coefficients, which control how much artificial dissipation is added. We employ the method with τη = 0, which we call LW+2 and which satisfies the geometric conservation law [13]. In all these tests we use τη = 0 and τu = τE = τ . In [13] we have used also the scheme LW+3 with τη 6= 0. 5. Numerical results In the original paper [13], all 2D numerical results were shown on a structured quadrilateral mesh, i.e., a logically rectangular (Cartesian) grid, although the actual method was already formulated for general polygonal meshes, using the p-c notation. Here we present the results on a variety of general meshes. In particular, the following topologies are used: First, we consider a honeycomb mesh, that is, a tessellation of the computational domain by regular hexagons. The boundary cells are trimmed accordingly to preserve the straight domain boundaries, which makes them non-hexagonal. Below, we call this a regular hexag- onal mesh. Splitting each hexagon into six triangles meeting at the hexagon’s center, and similarly split- ting the boundary cells, we obtain a structure which we will refer to as a regular triangular mesh. Note that the regular hexagonal (honeycomb) mesh can be seen as a Voronoi tessellation with regularly dis- tributed generators. Perturbing the positions of these generators and generating a new Voronoi tessellation, we obtain a mesh with cells being generally convex polygons. At this point, the mesh may need cleaning from too short edges and similar degeneracies. Fi- nally, splitting each of these polygons of N vertices into N triangles meeting at the polygon’s center, we obtain an irregular triangular grid. Note that while in the regular triangular grid each of the interior nodes was connected to six neighbors, now the connectivity varies and thus we have an unstructured triangular mesh. In the process of generating these test grids, we utilized the Voronoi tessellator ShaPo [18] and its advanced tools. 5.1. Time step control The time step ∆t is adaptively controlled by the CFL condition and limited by the multiple of the previous time step, so that it does not grow too fast. At time level n we set ∆tn = min ( CCFL min c √ V nc cs,nc , Cg∆tn−1 ) , (12) where cs,nc is the sound speed and we use Cg = 1.05 to limit the time step growth. In the numerical tests we set CCFL = 0.15. The Noh problem on the triangular mesh with CCFL = 0.5 generated small oscillations on the plateau of the solution. The time step is further constrained by the requirement that the cell volumes do not change much during one time step. If |V n+1 c /V nc −1| > Cv for some cell c, then the time step is reduced to a half ∆tn = ∆tn/2 and recomputed. In the tests we set Cv = 0.1. 5.2. Noh problem The Noh implosion test problem [19] is a spherical infinite strength shock propagating out from the origin. Its cylindrical version is modeled in the Cartesian (x, y) plane as follows. Initially, the density is set to unity and the pressure to zero everywhere, the velocity has a magnitude of 1 and is directed toward the origin, that is, u = −X/|X|, and the ideal polytropic gas equation of a state with adiabatic index γ = 5/3 is used. The exact solution at time t inside the shocked region (|X| ≤ t/3) is density ρ = 16, pressure p = 16/3, zero velocity, while outside the shock (|X| > t/3) we should have the exact density ρ = 1 + t/|X|, zero pressure and the velocity should remain u = −X/|X|. In our simulations, we approximate the zero initial pressure by p = 10−6. The results are shown at the usual final time tfinal = 0.6. Figures 3 and 4 show the density color maps result- ing from calculations, where the initial computational domain [−1, 1]× [−1, 1] was covered by a regular tri- angular mesh (described above) consisting of 10,108 cells, resp. by a regular hexagonal mesh of 10,083 cells. These meshes have a comparable numbers of cells and both will be further referred to as “coarse”. 71 D. Fridrich, R. Liska, I. Tarant et al. Acta Polytechnica 0.0 0.2 0.4 0.0 0.2 0.4 0 2 4 6 8 10 12 14 16 18 Figure 3. Density and mesh for Noh problem on coarse (10108-cell) regular triangular mesh with τ = 2. 0.0 0.2 0.4 0.0 0.2 0.4 0 2 4 6 8 10 12 14 16 18 Figure 4. Density and mesh for Noh problem on coarse (10083-cell) regular hexagonal mesh with τ = 1.25. Note that for each mesh we are showing the results with a different dissipation parameter τ . Indeed, as we can see in Figs. 5 and 6, the optimal dissipation coefficient, which prevents the oscillation near the shock but still does not smooth it too much, seems to depend on the mesh topology. For the regular trian- gular mesh a reasonable choice seems to be around τ = 2, while for hexagons it is about τ = 1.25. Let us remark, that for a regular quadrilateral grid, which is not shown here, the optimum is about τ = 1.5. This is in good agreement with the results shown in [13], which however used the default LF mass weighting in the predictor instead of (7). This mesh dependence of the optimal dissipation parameter deserves further study. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 2 4 6 8 10 12 14 16 Figure 5. Scatter plot of density for Noh problem on coarse regular triangular mesh with different values of dissipation parameter τ . 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 2 4 6 8 10 12 14 16 Figure 6. Scatter plot of density for Noh problem on coarse regular hexagonal mesh with different values of dissipation parameter τ . As for the symmetry and value of density on the plateau (inside the shocked region), the contour plots suggest that the regular hexagons in Fig. 8 do better than regular triangles in Fig. 7. The comparison of results of our LW+2 method and the standard stag- gered method with tensor artificial viscosity on logical 72 vol. 61 Special Issue/2021 Lagrangian Lax-Wendroff HLL Scheme 0.00 0.05 0.10 0.15 0.20 0.00 0.05 0.10 0.15 0.20 12.0 12.4 12.7 13.1 13.5 13.8 14.2 14.5 14.9 15.3 15.6 16.0 Figure 7. Contour plot of density for Noh problem on coarse regular triangular mesh with τ = 2. 0.00 0.05 0.10 0.15 0.20 0.00 0.05 0.10 0.15 0.20 12.0 12.4 12.7 13.1 13.5 13.8 14.2 14.5 14.9 15.3 15.6 16.0 Figure 8. Contour plot of density for Noh problem on coarse regular hexagonal mesh with τ = 1.25. rectangular quadrilateral mesh in [13] shows that our method preserves polar symmetry much better. Finally, to demonstrate the numerical convergence, we show the sequence of density scatter plots on meshes with the same topology and different reso- lutions, referred to as “coarse”, “medium” and “fine”, namely in Fig. 9 for the regular grid of 10,108, 40,280 and 160,816 triangles and in Fig. 10 for the regular hexagonal grid of 10,083, 40,037 and 159,561 cells. The mesh resolutions in both cases are chosen so that the cell numbers approximately correspond to those of regular quadrilateral meshes of 50×50, 100×100, resp. 200×200 initally square cells per quadrant, shown in the original paper [13]. (But again, keep in mind that the results shown there use the standard LF weighting.) A more thorough test of convergence will be pre- sented in Sec. 5.4. Note that we obtain decent results also on unstruc- tured meshes of 9,914 triangles in Fig. 11 and of 10,083 Figure 9. Convergence of scatter plots of density for Noh problem on regular triangular mesh with τ = 2. Figure 10. Convergence of scatter plots of density for Noh problem on regular hexagonal mesh with τ = 1.25. general polygons (ranging from triangles to nonagons) in Fig. 12. 5.3. Sedov problem The Sedov blast wave test problem [20] describes the evolution of a blast wave in a point-symmetric ex- plosion. Again, we consider its cylindrical version in Cartesian geometry. The ideal polytropic gas, this time with the adiabatic index γ = 7/5, is initially at 73 D. Fridrich, R. Liska, I. Tarant et al. Acta Polytechnica 0.0 0.2 0.4 0.0 0.2 0.4 0 2 4 6 8 10 12 14 16 18 Figure 11. Density and mesh for Noh problem on irregular triangular mesh of 9914 cells with τ = 2. 0.0 0.2 0.4 0.0 0.2 0.4 0 2 4 6 8 10 12 14 16 18 Figure 12. Density and mesh for Noh problem on ir- regular polygonal mesh of 10083 cells with τ = 1.25. rest, its density is set to 1 everywhere and the total energy E = 0.979264 is stored at the origin in the form of internal energy. At the final time tfinal = 1.0 the exact solution is a point-symmetric diverging shock whose front is at radius |X| = 1 and has a density peak ρ = 6. The exact solution shown in the scatter plots as a reference was obtained by the code Cococubed [21]. Note that while we are only showing one quadrant of the computational domain, the calculation is per- formed on the entire domain. This allows us to use a general polygonal mesh which is, unlike the regular meshes, non-symmetric. In the literature, one often encounters this test set up on a quadrilateral mesh with the calculation performed only on one quadrant, so that the total energy E = 0.244816, one quarter of ours, is stored in a single corner cell. In our im- 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0 1 2 3 4 5 6 Figure 13. Density and mesh for Sedov problem on regular triangular mesh of 10108 cells with τ = 2. plementation, all of the initial energy is stored is one regular hexahedral cell centered at (0, 0) in the case of the honeycomb and general polygonal meshes. Six triangular cells are created by its subdivision for the structured and unstructured triangular meshes. The initial computational domain [−1.2, 1.2] × [−1.2, 1.2] was covered by the same meshes as ear- lier for the Noh test, just now scaled to this bigger domain. In Fig. 13 we see the density color map on the coarse regular triangular mesh, in Fig. 14 on the coarse regular hexagonal grid, and in Fig. 17 on the same number of general polygons. The scatter plots in Fig. 15 show that with increas- ing resolution on regular triangles (coarse mesh at left and fine at right) we approach to the correct diverging shock with the front at radius |X| = 1 and a density peak ρ = 6 and to correct low density around the origin. Similarly in Fig. 16 we again compare the solutions of Sedov problem on the lowest and high- est resolution meshes of regular hexagons introduced earlier for the Noh test. 5.4. Smooth expansion problem To check the numerical order of accuracy of the LW scheme we define a smooth, radially symmetric, ex- pansion problem for a gas with γ = 5/3 on the cir- cular domain with the initial radius 3, centered at the origin. The initial density has Gaussian profile ρ0 = exp(−|X|2), the initial specific internal energy is constant ε0 = E0 = 3/4 and the gas is initially at rest. The final time was set to tfinal = 1. This polar symmetric problem has been also solved in 1D on a high resolution mesh with 5,000 cells. The obtained high resolution solution is used as a reference solution and instead of errors we are looking at the deviations of the 2D numerical solutions from this reference. 74 vol. 61 Special Issue/2021 Lagrangian Lax-Wendroff HLL Scheme 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0 1 2 3 4 5 6 Figure 14. Density and mesh for Sedov problem on regular hexagonal mesh of 10083 cells with τ = 1.25. 0 0.5 1 0 1 2 3 4 5 6 7 (a). (b). Figure 15. Density scatter plot for Sedov problem on regular triangular meshes with τ = 2 on coarse (a) and fine mesh (b). The 2D solutions were calculated on honeycomb meshes, which were constructed by the Voronoi tessel- lation of the initial circular domain. Since the solution is smooth, the calculations were performed by the pure LW method without any artificial dissipation. At the final time, the deviations from the linearly interpo- lated 1D reference solution were compared at the cell centers; i.e., for each cell we evaluate the radius of its center and compare the 2D value at the center with the value of the 1D reference solution linearly inter- polated to this radius. Table 1 shows the maximum deviations of density for consecutively refined meshes. Each finer mesh has about four times as many cells as 0 0.5 1 0 1 2 3 4 5 6 7 (a). 0 0.5 1 0 1 2 3 4 5 6 7 (b). Figure 16. Density scatter plot for Sedov problem on regular hexagonal meshes with τ = 1.25 on coarse (a) and fine mesh (b). 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0 1 2 3 4 5 6 Figure 17. Density and mesh for Sedov problem on irregular polygonal mesh of 10083 cells with τ = 1.25. than the previous one, which amounts to doubling the resolution in each direction. The ratio of the devia- tions is about 4, which confirms that the LW method is 2nd order accurate. 6. Conclusions We have presented the cell-centered hybrid Lax- Wendroff HLL Lagrangian scheme for compressible hy- drodynamics on unstructured computational meshes. The results of Noh and Sedov problems on uni- form honeycomb hexagonal, uniform triangular, non- 75 D. Fridrich, R. Liska, I. Tarant et al. Acta Polytechnica No. of cells Max. Ratio w. r. t. in the mesh deviation prev. resol. 75 0.036 309 0.0088 4.1 1 251 0.0021 4.2 5 012 0.00051 4.1 20 094 0.00013 4.0 Table 1. Maximum density deviations of 2D solutions from the interpolated 1D reference solution and the ratios of the deviations. uniform polygonal and non-uniform triangular meshes confirm the quality of the approach. The method preserves polar symmetry of the solution very well. Acknowledgements The first four authors have been supported by the Czech Science Foundation project 18-20962S, by the project CZ.02.1.01/0.0/0.0/16_019/0000778 from European Re- gional Development Fund, by the Czech Technical Univer- sity in Prague project SGS19/191/OHK4/3T/14 and by the Czech Ministry of Education project RVO 68407700. References [1] E. J. Caramana, D. E. Burton, M. J. Shashkov, P. P. Whalen. The construction of compatible hydrodynamics algorithms utilizing conservation of total energy. Journal of Computational Physics 146(1):227–262, 1998. doi:10.1006/jcph.1998.6029. [2] E. J. Caramana, M. J. Shashkov, P. P. Whalen. Formulations of artificial viscosity for multi-dimensional shock wave computations. Journal of Computational Physics 144(1):70–97, 1998. doi:10.1006/jcph.1998.5989. [3] E. J. Caramana, M. J. Shashkov. Elimination of artificial grid distortion and hourglass-type motions by means of Lagrangian subzonal masses and pressures. Journal of Computational Physics 142(2):521–561, 1998. doi:10.1006/jcph.1998.5952. [4] J. vonNeumann, R. D. Richtmyer. A method for the numerical calculation of hydrodynamic shocks. Journal of Applied Physics 21(3):232–237, 1950. doi:10.1063/1.1699639. [5] V. A. Dobrev, Tz. V. Kolev, R. N. Rieben. High-order curvilinear finite element methods for Lagrangian hydrodynamics. SIAM Journal on Scientific Computing 34(5):B606–B641, 2012. doi:10.1137/120864672. [6] B. Després, C. Mazeran. Lagrangian gas dynamics in two dimensions and Lagrangian systems. Archive for Rational Mechanics and Analysis 178(3):327–372, 2005. doi:10.1007/s00205-005-0375-4. [7] P.-H. Maire, R. Abgrall, J. Breil, J. Ovadia. A cell-centered Lagrangian scheme for two-dimensional compressible flow problems. SIAM Journal on Scientific Computing 29(4):1781–1824, 2007. doi:10.1137/050633019. [8] P.-H. Maire. A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes. Journal of Computational Physics 228(7):2391–2425, 2009. doi:10.1016/j.jcp.2008.12.007. [9] P.-H. Maire, B. Nkonga. Multi-scale Godunov-type method for cell-centered discrete Lagrangian hydrodynamics. Journal of Computational Physics 228(3):799–821, 2009. doi:10.1016/j.jcp.2008.10.012. [10] D. E. Burton, N. R. Morgan, T. C. Carney, M. A. Kenamond. Reduction of dissipation in Lagrange cell- centered hydrodynamics (CCH) through corner gradient reconstruction (CGR). Journal of Computational Physics 299:229–280, 2015. doi:10.1016/j.jcp.2015.06.041. [11] J. Cheng, C.-W. Shu. Positivity-preserving Lagrangian scheme for multi-material compressible flow. Journal of Computational Physics 257(A):143–168, 2014. doi:10.1016/j.jcp.2013.09.047. [12] W. Boscheri, M. Dumbser, R. Loubere, P.-H. Maire. A second-order cell-centered Lagrangian ADER-MOOD finite volume scheme on multidimensional unstructured meshes for hydrodynamics. Journal of Computational Physics 358:103–129, 2018. doi:10.1016/j.jcp.2017.12.040. [13] D. Fridrich, R. Liska, B. Wendroff. Some cell-centered Lagrangian Lax-Wendroff HLL hybrid schemes. Journal of Computational Physics 326:878–892, 2016. doi:10.1016/j.jcp.2016.09.022. [14] A. Harten, P. D. Lax, B. van Leer. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Review 25(1):35–61, 1983. doi:10.1137/1025002. [15] R. Loubère, M. Shashkov. A subcell remapping method on staggered polygonal grids for arbitrary-Lagrangian-Eulerian methods. Journal of Computational Physics 209(1):105–138, 2005. doi:10.1016/j.jcp.2005.03.019. [16] R. D. Richtmyer. A survey of difference methods for non-steady fluid dynamics. NCAR Technical Notes 63-2, National Center for Atmospheric Research, Boulder, Colorado, 1962. doi:10.5065/D67P8WCQ. [17] B. Wendroff, A. B. White. A supraconvergent scheme for nonlinear hyperbolic systems. Computers & Mathematics with Applications 18(8):761–767, 1989. doi:10.1016/0898-1221(89)90232-0. [18] J. Pouderoux, F. Chevassu, M. Kenamond, M. Shashkov. ShaPo - Recent advances on the Voronoi mesh generation toolkit. Presented at the 9-th International Conference on Numerical Methods for Multi-Material Fluid Flow (MULTIMAT 2019), Trento, Italy, September 9-13, 2019. [19] W. F. Noh. Errors for calculations of strong shocks using an artificial viscosity and artificial heat flux. Journal of Computational Physics 72(1):78–120, 1987. doi:10.1016/0021-9991(87)90074-X. [20] L. I. Sedov. Similarity and Dimensional Methods in Mechanics. Academic Press, New York, 1959. [21] J. R. Kamm, F. X. Timmes. Cococubed. Sedov verification code. LA-CC-07-020, 2007. http: //cococubed.asu.edu/research_pages/sedov.shtml. 76 http://dx.doi.org/10.1006/jcph.1998.6029 http://dx.doi.org/10.1006/jcph.1998.5989 http://dx.doi.org/10.1006/jcph.1998.5952 http://dx.doi.org/10.1063/1.1699639 http://dx.doi.org/10.1137/120864672 http://dx.doi.org/10.1007/s00205-005-0375-4 http://dx.doi.org/10.1137/050633019 http://dx.doi.org/10.1016/j.jcp.2008.12.007 http://dx.doi.org/10.1016/j.jcp.2008.10.012 http://dx.doi.org/10.1016/j.jcp.2015.06.041 http://dx.doi.org/10.1016/j.jcp.2013.09.047 http://dx.doi.org/{10.1016/j.jcp.2017.12.040} http://dx.doi.org/10.1016/j.jcp.2016.09.022 http://dx.doi.org/10.1137/1025002 http://dx.doi.org/10.1016/j.jcp.2005.03.019 http://dx.doi.org/10.5065/D67P8WCQ http://dx.doi.org/10.1016/0898-1221(89)90232-0 http://dx.doi.org/10.1016/0021-9991(87)90074-X http://cococubed.asu.edu/research_pages/sedov.shtml http://cococubed.asu.edu/research_pages/sedov.shtml Acta Polytechnica 61(SI):68–76, 2021 1 Introduction 2 Lagrangian Finite Volume 3 Meshes 4 Schemes 4.1 Predictor 4.2 Corrector 4.3 Mesh positions update 4.4 Artificial dissipation 5 Numerical results 5.1 Time step control 5.2 Noh problem 5.3 Sedov problem 5.4 Smooth expansion problem 6 Conclusions Acknowledgements References