Mr Ho

Ph.D Siao Hong

A Mathematical Proof for the Four-Color Problem In honor to the Mathematician August Ferdinand Möbius (1790-1868)

Keywords:

In this manuscript we present a mathematical proof in the traditional manner for the Four-Color Problem, without using any software program. We do so in two parts. In Part One we consider the four-color problem with the following constraint. Each country is a connected region. A proof using mathematical induction on the number of countries of a map is given. In Part Two we prove the Four-Color Problem without this constraint. No attempt is made to generalize the original four-color problem.

I.1. Part One: Preliminary Remarks.

The four-color problem says the following: on any two-dimensional map it is possible to label all countries using just four colors so that any two adjacent countries have dierent colors. In this context, a connected region is in fact simply connected. Two countries are adjacent if they share at least one common border. If they touch just at one point, then they aren't adjacent.

There have been numerous counter examples and proofs since the conjecture dawned in the mind of the mathematician August Ferdinand Möbius in 1840. It was formulated in a letter by another mathematician in 1852.

However, all these proofs and counter examples turned out to contain mistakes.

With the help of digital computer software programs some authors came up with huge number of maps with all kinds of possible country constellations.

Even mathematical proofs using software programs are accepted as correct, we provide a mathematical proof in the traditional manner. 1 2 I.2. Part One: The initial case.

A planar graph only partially describes the four-color problem. Each country is substituted by a vertex, and any two adjacent countries by an edge.

The initial case is a map with four countries using precisely four colors. As a planar graph it has four vertices and six edges. The Euler's formula V -E + R = 2 holds, where

V , E, and R denote the number of vertices, edges, and connected regions respectively.

Four colors are indeed needed in order to ensure that any two adjacent countries have distinct colors. We denote the colors of the four countries by 1, 2, 3, 4.

t 1 t 2 t 3 t 4 @ @ @ @ # " ! 1 3 4 2 I.
3. Part One: The Induction Hypothesis.

Let k be an integer. We assume that we are able to assign four given colors to k countries so that any two adjacent countries have distinct colors using the given four colors 1, 2, 3, 4.

Let l be an integer such that l < (k + 1) and l > (k/2 + 1). Then

(k + 1) -l < k + 1 -(k/2 + 1)≤k/2 < k.
We add one more country to the k existing countries adjacent to the region made up of these k countries. We want to show that we can assign colors to k + 1 countries satisying the coloring conditions.

I.4. Part One: The Induction Conclusion.

By induction hypothesis we are able to assign colors to the l countries satisying coloring conditions, since l≤k. These l countries form a simply connected colored surface area.

Again we apply the induction hypothesis. It is possible to assign colors to the (k + 1) -l countries satisfying the coloring conditions, sincce (k + 1) -l < k. Here we have to start assigning colors to countries adjacent to those countries that already have specic colors. Then we proceed assigning one color to each other remaining country until all (k + 1) -l countries are assigned colors.

Therefore, we are able to assign the given four colors to k + 1 countries, and any two adjacent countries have distinct colors. The Part One is done. q.e.d.

II.1. General Case of the Four-Color Problem.

II.1.1. Denition of a Big Country: A big country consists of more than one nonconnected regions located on a continent. Non-connected regions can be just islands.

They can have adjacent countries located on a dierent continent. The two continents are separated by oceans.

II.1.2. Example: Kaliningrad Oblast is a region of Russia, which is not connected to Russia. Between Kaliningrad Oblast and Russia there are Lithuania, then either Latvia or Belarus.

II.1.3. Remark: We take one region as the starting point of a big country. First, we link those non-connected regions to the big country, using paths that are constructed by linking one country after another like a chain. On each chain we choose countries so that they reach the big country. Note that those countries being used to build a chain don't belong to the big country.

In case a big country has just one non-connected region, we assign the region with the same color as the big country. Then we assign colors to the countries of the chain, to neighboring countries of the chain in order to obtain a simply connected superregion.

When a big country has a region located on a dierent continent, and has adjacent countries, then we assign colors to the latter in order to ensure that coloring conditions are met.

II.1.4. Example: The French oversea territory La Guyane is located in South America. It is assigned the same color as France in Europe. It is a good idea to assign colors to Brasil and Surinam. This precaution ensures that only the given four colors are used and any two adjacent countries have distinct colors.

French oversea territories, which are islands, are easy to assign colors. We simply assign them the same color as France in Europe.

II.1.5. Denition of Superregion: If a big country has more than one non-connected regions, we proceed in the same manner as above until all non-connected regions are linked to the big country and assigned the same color as the original big country.

We assign colors to participating countries in order to eventually construct a simply connected superregion on the same continent.

II.1.6. Remark: We note again that participating countries to form a superrregion don't belong to the big country, and are assigned colors so that only the given four colors are used and any two adjacent countries have distinct colors.

If there are more than one non-connected regions located on a dierent continent with the two continents separated by oceans, then the construction of a superregion on that continent is required to ensure that the coloring conditions are met.

If this situation occurs on several continents and the latter are all separated by oceans, then superregions are constructed on all pertinent continents.

Informally, a simply connected region has no hole. When there is a hole of countries within a superregion, then we cannot be sure that the coloring conditions are met.

II.2. Construction of a Hierarchy of Big Countries.

II.2.1. Remark: For any big country there are more than one possible superregions.

We can add a few not yet colored countries, or take away a few of them. Precisely this exibility enables us to construct superregions in such a manner to suit and meet coloring conditions. In addition, we ensure that any superregion becomes simply connected.

Given a map. We take one region of any big country, and attempt to construct a superregion as described in Part II.1.5. We make a note of this region. This big country corresponds to a vertex of a graph in graph theory.

Using a top-down approach we construct a graph representing a hierarchy of big countries, without yet assigning any colors to them.

If a big country consists of two non-connected regions, then it is always possible to link them using a chain of neighboring countries. None of which belongs to any big country.

But when there are three or more non-connected regions, then some of the participating countries may actually be regions of other big countries. Therefore, we need to establish a hierarchy of related big countries.

We exhaustively search for all participating countries which are actually regions of other big countries. If they exist, then we make a note of them and their corresonding big countries. Next we draw an edge between the vertex at the beginning and all new vertices underneath it.

We continue constructing the graph until all big countries related to each other, by having participating countries belonging to other big countries, are taken into consideration.

Such a graph has a nite number of vertices, because the number of big countries is nite.

What remains are other big countries which have nothing to do with this graph of vertices and their corresponding big countries. They will be processed later on. Now we consolidate the hierarchy of vertices relativ to each other, before we work on assigning colors.

Let us take one simple example to illustrate it and also show how we assign colors using a bottom-up aproach.

II.2.2. Example:

Let us consider a rst big country to begin with. Suppose there are two participating countries A and C which are actually regions of two other big countries. Call them the second and third countries.

For the purpose of coloring we take a region B of the second country. Similarly we take a region D in the third country.

Assume that we have no more big countries as participating countries. We only have two levels of vertices.

If the second country contains a participating country E which is a region of the third big country, then the vertex corresponding to the third country has to be one level lower than the vertex corresponding to the second country. We also need to nd a region F belonging to the third country to work with subsequently. Now instead of two levels of vertices, the graph of big countries has three levels of vertices. This is how we consolidate the hierarchy of vertices.

Next we work on the coloring of regions and countries in this example.

Since the third country is at the lowest level, we start our coloring works with it.

F and E have the same color. We assign colors to participating countries in a chain linking F to E. Likewise, D and C have the same color. We assign colors to a chain of participating countries linking D to C.

In this particular example we don't link A to C because they are just two participating regions of two distinct big countries and therefore have dierent colors. We work on them later.

E is just a participating country in the second country, but is a region of the third country. E and B have dierent colors.

Since F and D belong to the third big country and have the same color, we link them through a chain of participating countries.

Furthermore, when we know that any two non-connected regions belong to this same big country corresponding to the vertex at this lowest level, we connect them through a chain of neighboring countries.

In regard to assigning colors in the construction of a superregion:

1. Starting from one colored non-connected region we connect it to another one, with both belonging to the same big country.

2. Then we must assign colors to adjacent and neighboring countries to form a simply connected surface area. This activity also ensures that the coloring conditions will be met at the end of our works.

We repeat these two steps just mentioned until we are done with all non-connected regions. Then we continue assigning colors in order to construct a superregion of this big country.

Once a superregion is nished, we perform the same tasks to construct a superregion of the big country corresonding to a vertex at the same level.

If there is no more vertex at the same level, we go one level higher in the vertex hierarchy. Then we repeat the identical tasks keeping in mind that any two new superregions form a simply connected colored surface area.

We continue this identical process until all vertices of the given graph are processed.

As a result we obtain a cluster of related big countries, whose superregions together remain simply connected.

II. 4 .

 4 Processing Other Clusters of Big Countries. If there are other clusters of related big countries like in Part II.2 and Part II.3, we work in the same manner until all simply connected superregions are colored satisfying the coloring conditions.

Next we assign the same color to one non-connected region of the third country either from the region F or D. Then we also assign colors of adjacent and neighboring countries to form a simply connected area before we connect to the next non-connected region.

We repeat the same task to all other remaining non-connected regions of the third country and ensure that the growing colored superregion remains simply connected.

The regions F , D, and all other non-connected regions of the third big country have an identical color.

After this superregion is done, we work on coloring the second country which corresponds to a vertex at one higher level.

We link B to A through a chain of participating countries. Starting with the region B we construct a superregion to the second country, ensuring that the two superregions of the third and second big countries together remain simply connected.

Since A and C are just participating regions of the second and third big countries It is fairly simple to construct a superregion to the rst big country, again ensuring that the new colored surface area remains simply connected.

II.3. Coloring a Hierarchy of Big Countries.

From the II.2.2. Example we get an idea how to assign colors in case of related big countries.

After the vertices of the graph are consolidated, we start with vertices at the lowest level in any branch of the graph.

From bottom-up we establish which vertices at higher levels are related to a vertex at a lowest level. Then we link one non-connected region of the big country corresponding to this vertex being considered, through a path of neighboring countries, to a participating region of a big country corresponding to a vertex at a higher level.

We repeat this task until all pertinent vertices at higher levels are taken care of.

When two simply connected colored clusters located on the same continent are going to merge, we have to assign colors so that the new growing merged surface area still remains simply connected.

Islands surrounded by oceans and belonging to any big country are always assigned the same color as the remaining regions of that specic big ccountry.

II.5. The Final Step of the General Proof.

After all individual clusters of big countries are processed, there remain only connected countries that are NOT located inside any simply connected regions.

By Part I we rst assign colors to those countries adjacent to countries that already have a xed color, keeping the coloring conditions in mind. One by one, all these remaining countries are assigned one of the four given colors and any two adjacent countries have distinct colors.

Therefore, the four-color conjecture, rst initiated by the mathematician August Ferdinand Möbius in 1840, is proved. Q. E. D. III. Part Three: Concluding Remarks and Comments.

1. At each step of assigning colors to countries and/or regions of big countries, the growing colored surface must remain simply connected. This is the reason why I construct a hierarchy of related big countries in a top-down manner, then work in a bottom-up approach when it comes to assigning colors.

The notion of simple connectivity is indeed central to this proof of the Four-Color Problem.

2. Since I didn't study any research articles published by any other authors, there is no bibliography at the end in this research manuscript.


Ivan Martinez
email: ivanmartinez@usantotomas.edu.co

Victor Ramos

Choosing a TCP version over static ad hoc wireless networks: wired TCP or wireless TCP?

Keywords: Ad hoc wireless networks, TCP Reno, TCP Westwood, DSDV

The Transmission Control Protocol (TCP) was originally designed to operate on wired networks. However, nowadays the traffic on wireless networks has grown to a point where one must take into account the specific characteristics of such networks when setting up a particular TCP implementation on them. Particularly, the TCP performance has been studied over ad hoc wireless networks leading to several new implementations for TCP. The main issue for TCP on ad hoc wireless networks is to differentiate between losses due to congestion and losses occurred at lower network layers. In this paper, we analyze the TCP performance on two different scenarios of static ad hoc wireless networks over the DSDV routing protocol. Our findings show that choosing the right TCP version for ad hoc wireless networks is a key factor for their performance. We show for the scenarios we study, that TCP Reno outperforms TCP Westwood on the average loss rate as well as on the throughput.

I. INTRODUCTION

Unlike traditional wireless networks that are provided with a centralized infrastructure, wireless ad hoc networks are formed solely of a set of nodes with identical or similar characteristics that perform simultaneously functions of a station as well as of a router [START_REF] Wei | Analysis of the applications and characteristics of ad hoc networks[END_REF]. These characteristics make ad hoc networks very attractive in terms of cost, scalability, and ease of implementation and thus making them well suited for a wide set of applications as in wireless sensor networks (WSN), emergency situation applications, monitoring of specific locations, mobile object tracking, location characterization, characterizing places, and military applications just to mention a few. Figure 1 depicts a simple example of an ad hoc network.

The deployment of wireless protocols such as Bluetooth, IEEE 802.11, and Hyperlan makes possible the implementation of ad hoc networks for commercial purposes. Therefore, several research efforts have been made during the last decade to solve the main problems faced by this type of network architectures. Wireless ad hoc networks may be classified as mobile ad hoc networks (MANETs) and static ad hoc networks (SANETs). On one hand, in MANETs nodes are allowed to move within a given area, and thus mobility introduces additional challenges to overpass for the correct operation of the network. On the other hand, wireless nodes in a SANET are placed within a given area and do not move.

The main problems faced by ad hoc networks include the hidden terminal problem, the exposed terminal problem, path asymmetry, network partition, routing problems, and power constraints. Additionally, in order that applications running over ad hoc networks operate well, a transport layer protocol is often needed. Therefore, several research efforts have been done on the direction of evaluating the performance of the transmission control protocol (TCP) on ad hoc networks [START_REF] Al-Hanbali | A survey of TCP over ad hoc networks[END_REF]- [START_REF] Henna | A throughput analysis of TCP variants in mobile wireless networks[END_REF].

The TCP protocol was originally designed to operate over wired networks. However, the wireless nature of ad hoc networks adds new problems to the operation of TCP. The TCP protocol interprets any packet loss as a sign of congestion; however, wireless networks may induce losses due to interference or collisions, which may be badly interpreted by TCP and thus impacting its performance. This is why several versions of TCP for wireless networks have been proposed in the literature to circumvent the deficiencies of TCP when losses other than those due to congestion take place within a TCP connection.

In this work, we compare the performance of two different versions of TCP on a SANET: TCP Reno and TCP Westwood; the former being originally designed for wired networks and the latter designed to operate on wired as well as on wireless networks. Our aim is to study the performance of the TCP protocol over SANETs arranged in a mesh, which is a classical node placement for several practical applications. The results we obtain during our analysis show that even if TCP Westwood is designed to operate on wired as well as on wireless networks, it is not always the best choice when selecting a transport protocol for the type of wireless networks as those we consider in this paper, because TCP Reno exhibits a better performance for every case we study.

The remainder of this paper is as follows. In Section II, we briefly describe the background and related work of TCP over ad hoc networks. Section III presents the simulation scenarios we consider for our evaluation. In Section IV, we describe the results we obtain by comparing TCP Reno and TCP Westwood regarding the throughput and the loss rate over our scenarios. Finally, Section V sketches our conclusions.

II. BACKGROUND AND RELATED WORK

Several research has been done on TCP over SANETs, the main efforts have been done on how TCP may be well tuned to operate on multi-hop wireless SANETs. In [START_REF] Xu | Performance evaluation of TCP algorithms in multi-hop wireless packet networks[END_REF], the authors compare the performance of five different TCP implementations: Tahoe, Reno, Reno, Sack, and Vegas over a multihop network when using the 802.11 protocol; by adequately tunning the TCP's advertised window the authors find that TCP Vegas outperform the other TCP implementations being besides the fairer TCP version. Moreover, the same author shows in [START_REF] Xu | TCP behavior across multihop wireless networks and the wired Internet[END_REF] that even in scenarios where TCP connections cross wired as well as wireless networks, the TCP unfairness problem persists.

Later, in [START_REF] Chen | On setting TCP's congestion window limit in mobile ad hoc networks[END_REF] the authors analyze the effect of the TCP's congestion window limit on throughput for SANETs. They also focus on a multi-hop scenario for a TCP connection and find that the MAC layer has a big impact on TCP performance. Thus, most of the work on TCP over SANETs is focused on how to well tune TCP to operate well on such kind of networks as well as on facing the problem of how to differentiate between losses due to congestion and losses due to phenomena on lower protocol layers. Therefore, several new versions of TCP for wireless networks have been proposed to date [START_REF] Fu | TCP Veno: TCP enhancement for transmission over wireless access networks[END_REF][START_REF] Wang | Adaptive bandwidth share estimation in TCP westwood[END_REF][START_REF] Akyldiz | TCP-Peach+: enhancement of TCP-Peach for satellite IP networks[END_REF] in order to circumvent such TCP deficiencies.

The main conclusion is that, when implementing an application over an ad hoc network it is often better to choose a TCP version for wireless networks rather than one for wired networks. In this paper, we show that this is not always the best choice. We compare two different TCP implementations: Reno and Westwood. The latter has been proposed to operate on wired as well as on wireless networks, while the former is designed to operate on wired networks. By means of simulation with the NS-2 simulator, we find in the following sections that TCP Reno outperforms better than Westwood for all of the cases we study on two different SANET scenarios. The reader is referred to [START_REF] Henna | A throughput analysis of TCP variants in mobile wireless networks[END_REF][START_REF] Bhanumathi | TCP Variants a comparative analysis for high bandwidth-delay product in mobile ad hoc network[END_REF] for complementary analyses regarding MANETs.

A. Node placement on a SANET

The node placement on a SANET is choice that depends on the application. Most of the applications for SANETs place nodes arranged in a way such that they cover a whole surface of interest. The most widely used node arrangement is a mesh, but different topologies have been proposed which may better adapted to a particular range of applications [START_REF] Biagioni | Wireless sensor placement for reliable and efficient data collection[END_REF].

We show in Figure 2 an example of different topologies for mesh deployment in static ad hoc networks. We can find in the literature several application examples that use different node topologies. For instance, in [START_REF] Kovács | Case study of a simple , low power WSN implementation for forest monitoring 2 overview of WSN systems[END_REF][START_REF] Al-Turjman | Connectivity optimization for wireless sensor networks applied to forest monitoring[END_REF] the authors implement networks with smoke sensors for early fire detection.

The nodes in our analysis are placed in a square mesh with 20 m of separation between nodes. We analyze two different scenarios, one with a TCP connection traversing the whole network in order to stress the routing protocol by having the highest number of hops, and the other one with several simultaneous TCP connections sending data to a receiver located at the center of the mesh.

III. SIMULATION SCENARIOS AND PERFORMANCE

MEASURES

A wide set of applications benefit from placing nodes arranged as a mesh, we thus choose such scenario for evaluation. In [START_REF] Biagioni | Wireless sensor placement for reliable and efficient data collection[END_REF], the authors provide an analytical model allowing to represent such scenario as the one shown in Figure 2(b), their model is represented as follows:

d = 0.71 × |A| N (1) a = |A| N (2) N s = 0.50 × |A| d 2 (3)
N t = |A| a 2 (4)
where n represents the maximum allowed distance between any point from the deployment zone and the nearest sensor, a is the length of a tile's side, N is the number of nodes, and |A| is the total area to cover. Depending of the type of zone to cover and the type of sensors, there is an additional parameter, δ, representing the sampling distance which must accomplish δ ≥ d. In order to guarantee the optimal communication among the whole set of nodes, the transmission range r should fulfill r ≥ a.

We choose the IEEE 802.11 standard for our simulations since it is the most widely used to deploy wireless networks, it is also easily scalable and cost convenient. Additionally, the 802.11 transmission speeds are also according to those used in WSNs, and finally the indoor coverage radio for this standard is 38 m.

We focus our attention on a set of applications benefiting from implementing the scenario chosen. Thus, several WSN implementations where the main application is forest monitoring to warn from fire of from illegal logging, where the δ parameter depends of the coverage radio of each node [START_REF] Lei | Distributed coverage of forest fire border based on WSN[END_REF][START_REF] Papán | WSN for forest monitoring to prevent illegal logging[END_REF]. We thus take these considerations into account and set this value to δ = 20m. Finally, the surface were we deploy the set of nodes is a 160m × 160m square, which results in the set of parameters given in Table I We therefore consider the two following scenarios. The first one consisting of a network with 64 static wireless nodes arranged in a mesh, where nodes are spaced by 20 m. There is a TCP connection between the end nodes in a manner such that we are able to evaluate a scenario where the connection traverses the maximum number of hops, and thus stressing the routing protocol to find out the path with the best conditions.

3 represents the first scenario. We consider the same node mesh arrangement for the second scenario, where we set up 64 simultaneous TCP connections between each node and the one located at the center, which acts as receiver. This scenario is closer to a wireless static sensors network characterizing a real static WSN. We show such zone definition in Figure 4.

We use the NS-2 network simulator for our evaluation. Table II lists the parameters we consider for our scenarios. For each case, we run the simulation 1000 times computing average values for throughput as well as for losses.

IV. RESULTS

We present in this section the results for both scenarios, where we measure the throughput as well as the loss rate for to the first scenario and in Figures 7 and8 the results corresponding to the second scenario.

Regarding our first scenario, we clearly see how in Figures 5 and6 TCP Reno outperforms TCP Westwood regarding the loss rate as well as the throughput. This may be explained by the fact that even if the nodes involved in the TCP connection are opposite end corners in the mesh, the DSDV routing protocol establishes proactively its routing tables for all nodes leading to the case where most of the losses are due to congestion rather than to other factors (e.g., interference, path asymmetry, route changes). TCP Westwood is particularly designed to better handle large bandwidth-delay product paths (large pipes), considering also losses due to transmission.

Finally, regarding our second scenario we see in Figures 7 and8 how, again, TCP Reno outperforms TCP Westwood even for the case where all nodes are sending data to the center node. We note how in this case, the loss rate grows considerably compared to the previous scenario which is clearly explained by the fact that the whole set of nodes is active during transmission. In the case of throughput, we measure the average throughput at the center node, that is why the average throughput is considerably lower than for the first scenario.

V. CONCLUSIONS

In this paper, we evaluated the performance of TCP Reno and TCP Westwood over a SANET arranged in a mesh topology. Two scenarios have been analyzed, the first one allowed us to evaluate the network performance for the nodes located at extreme corners of the mesh. In this case, the routing protocol makes the TCP connection traverse the higher number of hops within this topology.

The second scenario we studied was intended to evaluate the network performance in the presence of several active TCP connections. Our results for both scenarios reveal that choosing a TCP implementation for wireless networks might not always be the better choice for SANETs since during our evaluations TCP Reno behave slightly better than TCP Westwood.

  This work has been done during the Ivan Marino's spring research internship at the Networking and Telecommunications Research Team of UAM-I in Mexico City.

Figure 1 .

 1 Figure 1. Example of an ad hoc network.

2013

  Seventh International Conference on Next Generation Mobile Apps, Services and Technologies 978-0-7695-5090-9/13 $26.00 © 2013 IEEE DOI 10.1109/NGMAST.2013.38

Figure 2 .

 2 Figure 2. Most used topologies for static ad hoc networks: (a) Hexagon, (b) Square, and (c) Triangle.

Figure 3 .

 3 Figure 3. Our first scenario.

Figure 4 .Figure 5 .

 45 Figure 4. Our second scenario.

Figure 6 .

 6 Figure 6. Throughput for the first scenario.

Figure 7 .

 7 Figure 7. Loss rate for the second scenario.

Figure 8 .

 8 Figure 8. Throughput for the second scenario.

Table I PARAMETERS

 I FOR THE NETWORK MODEL.

	A	25600m 2
	a	20m
	d	14.2m
	N	64


Thomas Heuzé
email: thomas.heuze@ec-nantes.fr

Lax-Wendroff schemes for elastic-plastic solids

Keywords: Elastic-plastic solids, Lax-Wendroff, Associated and non-associated plasticity

Two formulations of the Lax-Wendroff scheme are proposed in this paper for its extension to elastic-plastic solids. The unknown vector consists of the strain or stress components in addition to the velocity ones according to the chosen formulation. The Lax-Wendroff scheme is here implemented in its Richtmyer two-step version so that the projection of the elastic trial stresses onto the yield locus is performed twice per time step. It is shown that it allows to reduce the number of cells required for the approximation of discontinuous plastic waves with respect to a classical one-step strategy consisting in an a posteriori projection of the elastic trial stresses onto the yield locus. Elastic-plastic associated and non-associated constitutive models under small strains are considered in examples. In particular, dynamic ratchetting is simulated with the Armstrong-Frederick nonlinear kinematic hardening. Comparisons with respect to finite element numerical solutions show good agreements.

Introduction

The simulation of impacts on solids finds many applications in different domains such as crash or high speed forming processes. They are modeled by means of hyperbolic initial boundary value problems, whose solutions consist of possibly discontinuous waves that are important to be tracked in order to properly understand the physical phenomena occurring in impacted parts. Indeed for solid-type media, waves can propagate both continuous and discontinuous solutions, as well as large strains and irreversible phenomena. Further, wave paths provide information about the loading history undergone by each material particle in time, which in turn allows to compute the residual strains and stresses once the steady state is achieved.

The numerical simulation of impacts on dissipative solids has been and is again mainly performed with the classical finite element method coupled with the centered differences or Newmark finite difference schemes in time [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF][START_REF] Benson | Computational methods in Lagrangian and Eulerian hydrocodes[END_REF], which is implemented in many industrial codes. Indeed, the finite element method enables an easy management of nonlinear partial differential equations, especially for solid-type media for which history-dependent constitutive equations can be implemented by means of appropriate integration algorithms [START_REF] Simo | Computational inelasticity[END_REF], and internal variables are stored at integration points in each element. However, it is first well-known that explicit finite element schemes can exhibit numerical noise near sharp solutions that can be hard to remove with artificial viscosity without loss of accuracy. Second, and especially for elastic-plastic simulations, the stresses and plastic strains converge at a rate that is one order lower than that of the displacement components.

On another side, the finite volume methods initially developed for the simulation of gas dynamics [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF], have gained recently more and more interest for problems involving impacts on solid media (see e.g. [START_REF] Barton | A high-order Eulerian Godunov method for elastic-plastic flow in solids[END_REF][START_REF] Hill | An Eulerian hybrid WENO centered-difference solver for elastic-plastic solids[END_REF][START_REF] Favrie | Mathematical and numerical model for nonlinear viscoplasticity[END_REF][START_REF] Lee | Development of a cell centred upwind finite volume algorithm for a new conservation law formulation in structural dynamics[END_REF][START_REF] Ortega | Numerical simulation of elastic-plastic solid mechanics using an Eulerian stretch tensor approach and HLLD Riemann solver[END_REF][START_REF] Maire | A nominally second-order cell-centered Lagrangian scheme for simulating elastic-plastic flows on two dimensional unstructured grids[END_REF][START_REF] Aguirre | An upwind vertex centred finite volume solver for lagrangian solid dynamics[END_REF][START_REF] Ndanou | Multi-solid and multi-fluid diffuse interface model: Applications to dynamic fracture and fragmentation[END_REF]). This family of methods is constructed based on the characteristic structure of hyperbolic equations, either introduced explicitly by means of the solution of a Riemann problem at cell interfaces for upwind cell-centered schemes, or implicitly for centered schemes through an average process to define the updated quantities. One interesting feature of these schemes is that the same order of convergence is achieved for all components of the unknown vector, and in particular for both the velocity and stress fields in solid mechanics [START_REF] Lee | Development of a cell centred upwind finite volume algorithm for a new conservation law formulation in structural dynamics[END_REF]. Next, some schemes permit to control locally the amount of numerical viscosity introduced as a function of the local regularity of the solution, so as to deliver non-oscillatory numerical solutions while preserving a high order of accuracy in more regular zones [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF].

Since the early work of Wilkins [START_REF] Wilkins | Methods in computational physics[END_REF] and Trangenstein et al. [START_REF] Trangenstein | A higher-order Godunov method for modeling finite deformation in elastic-plastic solids[END_REF], several authors have proposed many ways to simulate impacts on dissipative solid media with this class of methods, both for rate-independent and rate-dependent solid behaviours. Eulerian [START_REF] Miller | A high-order eulerian godunov method for elastic-plastic flow in solids[END_REF][START_REF] Hill | An Eulerian hybrid WENO centered-difference solver for elastic-plastic solids[END_REF][START_REF] Barton | A high-order Eulerian Godunov method for elastic-plastic flow in solids[END_REF][START_REF] Ortega | Numerical simulation of elastic-plastic solid mechanics using an Eulerian stretch tensor approach and HLLD Riemann solver[END_REF][START_REF] Ndanou | Multi-solid and multi-fluid diffuse interface model: Applications to dynamic fracture and fragmentation[END_REF] as well as Lagrangian [START_REF] Maire | A nominally second-order cell-centered Lagrangian scheme for simulating elastic-plastic flows on two dimensional unstructured grids[END_REF][START_REF] Aguirre | An upwind vertex centred finite volume solver for lagrangian solid dynamics[END_REF][START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elasticplastic solids[END_REF] approches have been proposed, each pulling its own profits according to the level of strain reached, mesh entanglement and refined history-dependent constitutive equations. One key point that should be emphasized is that for rate-independent plasticity, no elastic-plastic Riemann solvers exist for two or three-dimensional stress states, even approximate. Only that of Lin & Ballman [START_REF] Lin | A Riemann solver and a second-order Godunov method for elastic-plastic wave propagation in solids[END_REF] was introduced for combined tension-torsion loadings, and is based on the analytical solution provided by Clifton [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF]. Anyway these solvers are generally considered too complex or costly so that upwind schemes are usually coupled to linear elastic/acoustic Riemann solvers to define an elastic trial solution, and then either it is projected onto the yield locus for rate-independent constitutive behaviour [START_REF] Maire | A nominally second-order cell-centered Lagrangian scheme for simulating elastic-plastic flows on two dimensional unstructured grids[END_REF][START_REF] Fridrich | Cell-centered lagangian Lax-Wendroff HLL hybrid method for elastic-plastic flows[END_REF] or a relaxation step associated with a splitting scheme is performed for rate-dependent solids. The latter is often called the Maxwell-type relaxation approach [START_REF] Godunov | Elements of continuum mechanics and conservation laws[END_REF][START_REF] Hill | An Eulerian hybrid WENO centered-difference solver for elastic-plastic solids[END_REF][START_REF] Barton | A high-order Eulerian Godunov method for elastic-plastic flow in solids[END_REF][START_REF] Ortega | Numerical simulation of elastic-plastic solid mechanics using an Eulerian stretch tensor approach and HLLD Riemann solver[END_REF][START_REF] Ndanou | Multi-solid and multi-fluid diffuse interface model: Applications to dynamic fracture and fragmentation[END_REF][START_REF] Favrie | Mathematical and numerical model for nonlinear viscoplasticity[END_REF] which actually refers to an adapted version of Perzyna's elastic-viscoplastic solids [START_REF] Perzyna | Fundamental problems in viscoplasticity[END_REF], and converge to the rate-independent solution when the relaxation parameter tends to zero, provided the solution scheme employed to solve the source term is asymptotic preserving and accurate [START_REF] Heuzé | Simulation of impacts on elasticviscoplastic solids with the splitting finite volume method applied to non-uniform quadrilateral meshes[END_REF][START_REF] Dumbser | Finite volume schemes of very high order of accuracy for stiff hyperbolic balance laws[END_REF]. However, both approaches generally yield a poor solution of discontinuous plastic waves (e.g. plane wave with linear hardening), which are smeared over a lot of cells even though the Courant number is set at one since the elastic precursor travels faster. Conversely, they behave better in more general cases for which plastic waves are simple waves (for nonlinear hardenings or non plane wave propagation).

One way to go around the former approaches while still avoiding elastic-plastic Riemann solvers can be done by looking at the well-known Lax-Wendroff scheme [START_REF] Lax | Systems of conservation laws[END_REF], that has two interesting features for the computation of elastic-plastic solids. First, it belongs to the class of centered schemes [START_REF] Toro | Centred TVD schemes for hyperbolic conservation laws[END_REF] so that no Riemann solver is required. Second, its Richtmyer two-step implementation [START_REF] Richtmyer | A survey of difference methods for non-steady fluid dynamics[END_REF][START_REF] Zwas | On Two Step Lax-Wendroff Methods in Several Dimensions[END_REF] allows to perform the projection of the elastic trial stresses onto the yield locus twice per time step, which can improve the solution of plastic waves in case these are discontinuous. The purpose of this work is to develop this idea. Therefore, two formulations of the Lax-Wendroff scheme are proposed in this paper for its extension to the computation of elastic-plastic solids. The unknown vector consists either of the strain or stress components in addition to the velocity ones according to the chosen formulation, hereafter denoted as 'strain' and 'stress' ones respectively. In the 'strain' formulation, the integration of the elastic-plastic constitutive equations is performed twice per time step, first at nodes, then cellwise. The 'stress' one amounts, for its first step, to solve a small system of equations at nodes to update the stresses, provided some effective elastic-plastic compliance moduli. Then, the stress update is performed cellwise during the Richtmyer second step and still accounts for the yield threshold. At last, the plastic strains and internal variables are updated iteratively by adjusting the elastic convex so that the updated stress state lies on its boundary, consistently with the plastic flow rule and hardening laws. The 'stress' formulation presented here is actually an extension to two space dimensional problems of that proposed in [START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elasticplastic solids[END_REF]. It should be noted that another extension of the Lax-Wendroff scheme to hypoelastic-plastic solids has already been proposed in [START_REF] Fridrich | Cell-centered lagangian Lax-Wendroff HLL hybrid method for elastic-plastic flows[END_REF]. However, only one projection onto the yield locus is performed per time step, and the formulations and strategies presented in this work are completely different.

The paper is organized as follows. First, the elastic-plastic initial boundary value problem is presented in Section 2. Next, the Lax-Wendroff scheme derived in its Richtmyer two-step version is presented in Section 3. The proposed 'strain' and 'stress' formulations are derived in this section. Section 4 is devoted to numerical applications. Several examples conducted in two-dimensional plane strain conditions with associated and non-associated elastic-plastic constitutive models are considered in the linearized geometrical framework, and comparisons with analytical and finite element numerical solutions are presented. In particular, an example of dynamic ratchetting is simulated with the Armstrong-Frederick [START_REF] Frederick | A mathematical representation of the multiaxial Baushinger effect[END_REF][START_REF] Frederick | A mathematical representation of the multiaxial Bauschinger effect[END_REF] nonlinear kinematic hardening to show the possibilities of the proposed approach.

Elastic-plastic initial boundary value problem

Conservation laws

We consider a three-dimensional continuum body Ω ∈ R 3 , of boundary ∂Ω, and current coordinates x = x i e i . Considering the linearized geometrical framework and isothermal setting, the conservation of linear momentum and the geometrical compatibility should be satisfied for the time interval of interest t ∈ τ =]0, T ]:

ρv -∇ • σ = 0 (1) ε - 1 2 (∇v + (∇v) T ) = 0 (2)
which hold ∀(x, t) ∈ Ω×]0, T ], and where v, σ and ε refer to the velocity vector, the Cauchy stress and the linearized strain tensors respectively. Equation (2) can also be written in conservative form:

ε -∇ • v ⊗ 1 + 1 v 2 = 0 ( 3
)
where 1 is the second order identity tensor, and the operator is defined in the cartesian coordinate system as 1 v = δ ik v j e i ⊗ e j ⊗ e k , δ ik being the Kronecker delta symbol. The system formed by Equations ( 1) and (3) appears as a system of conservation laws, written as:

∂U ∂t + ∇ • F = 0 (4)
where U and F denote the arrays of conserved quantities and associated fluxes respectively, defined as follows:

U = v ε ; F =        -σ ρ -v⊗1+1 v 2        (5)
More precisely, System (4) can be derived in full matrix form with cartesian coordinates as:

∂U ∂t + 3 k=1 ∂F k ∂x k = 0 (6)
with vector components such that

U = v {ε • e i } 1≤i≤3 ; F k =        -σ•e k ρ - v p δ ik +v i δ pk 2 e p 1≤i≤3        = F • e k (7)
where e i (1 ≤ i ≤ 3) refer to the cartesian basis vectors, and the redundant equations due to the symmetry of the strain tensor are not considered. The fluxes F depend on both the velocity and stress components, which also defines the following auxiliary vector:

W = v σ (8)
The closure of the above system of conservation laws (4) is performed by means of a set of constitutive equations, introduced thereafter. Besides, Systems (4) or ( 6) must also be supplemented with appropriate initial and boundary conditions.

Elastic-plastic constitutive model

The thermodynamical state of the material is described by a set of state variables [START_REF] Germain | Continuum Thermodynamics[END_REF][START_REF] Lemaitre | Mechanics of Solid Materials[END_REF] that consists of the linearized strain tensor ε and the temperature T , plus some internal state variables Z describing the evolution of the internal microstructure and stored energy due to plastic deformations and other irreversible processes. Under the small strains assumption, the total strain ε is additively decomposed into an elastic strain (ε e ) and a plastic strain (ε p ):

ε = ε e + ε p (9)
The set of internal state variables Z consists here of the plastic strain ε p plus some additional variables α I , 1 ≤ I ≤ N. From a Helmholtz free-energy density potential ψ(ε, T, Z) ≡ ψ(εε p , T, α I, 1≤I≤N ), concave with respect to the temperature and convex with respect to the other variables, the state laws are derived as:

σ = ρ ∂ψ ∂ε e ; s = - ∂ψ ∂T ; A I = ρ ∂ψ ∂α I , 1 ≤ I ≤ N (10)
where s and A I refer to the entropy density and the force conjugated to the variable α I respectively. Thus, the mechanical dissipation reads:

D = σ : εp -A I • αI = Y • Ż ≥ 0 (11) with Y = {σ, A I } ; Z = {ε p , -α I } (12)
where the vector Y denotes the thermodynamic forces conjugated to the internal variables Z, and the dot ( ˙ ) applied to the quantity stands for a time rate. A convex elastic domain

C = {Y| f (Y) ≤ 0} (13)
is defined by means of the yield function f , and plastic flow occurs if the current state lies and remains on its boundary, defined by f = 0. Then, provided a dissipation pseudo-potential Φ(Y), convex with respect to the flux variables Y, positive (Φ(Y) ≥ 0), and containing the origin (Φ(0) = 0), the flow rules read:

Ż = λ ∂Φ ∂Y ; λ ≥ 0 ( 14
)
where λ denotes the plastic multiplier satisfying the complementarity conditions λ ≥ 0, f (Y) ≤ 0, λ f = 0, the identity Φ λ ≥ 0, and is determined by means of the consistency condition ḟ = 0. The convexity of the flow function Φ(Y) ensures that the mechanical dissipation

D = λY • ∂Φ ∂Y ≥ 0 ( 15
)
is non-negative. The definition of both the flow and yield function Φ and f leads to the nonassociative framework for plasticity, which may be useful for instance for the modeling of nonlinear kinematic hardening. The associative case appears as a particular case of the former one, identifying both flow and yield functions.

3. Lax-Wendroff numerical schemes

Mesh notations

The computational domain is divided into N c non-overlapping polygonal cells. The nodes p = 1, . . . , N n are defined at cell vertices. Following the pc notations introduced in [START_REF] Maire | A nominally second-order cell-centered Lagrangian scheme for simulating elastic-plastic flows on two dimensional unstructured grids[END_REF][START_REF] Fridrich | Cell-centered lagangian Lax-Wendroff HLL hybrid method for elastic-plastic flows[END_REF], the set of nodes connected to a cell c is denoted by P(c), that of cells connected to a node p by C(p), that of edges linked to a cell c by E(c), and that of edges linked to a node p by E (p). The outward unit normal to a cell c at an edge e ∈ E(c) is denoted by n (c) e . The dual mesh consists of non-overlapping polygons surrounding a node p, so that the vertices of that polygon consist of the centroids of the cells C(p) connected to a node p, and the midpoints of the primal edges e ∈ E (p) connected to a node p; see Figure 1(a). Thus, the edges of dual cells consist of separator segments connecting the centroids of primal cells to the midpoints of the primal edges E (p). The overlap between the primal and dual cells c and p forms the subcell pc, which is a quadrangle. It is defined by the separators spc-and spc+, of lengths l spc-, l spc+ and outward unit normal n spc-, n spc+ respectively, and by the half of the edges e ∈ {p -p, pp + }. The length of these edges are denoted l e , and their outward unit normals are n pc-, n pc+ respectively. The volumes (actually areas in two-dimensions) of primal and dual cells are hence related to these of subcells by the following relations:

V c = p∈P(c) V pc ; V p = c∈C(p) V pc (16)
Notice that a boundary node may be connected to two or three edges, so that its dual cell consists of one or two subcells respectively (see Figure 1(b)).

• • • • • ∂Ω p p + p - c n spc- n spc+ npc- V p npc+ V p e (a) c • • • • ∂Ω p p - p + V p e n spc- n spc+ n pc+ (b)

Approximation

Consider a polygon of volume V, one defines the following approximation within this domain:

U V = 1 V V U dV (17)
where U V denotes the unknown array, that is the integral average of the array of conserved variables within the polygon V, which can stand for the primal cell c or the dual one p thereafter. Let's consider the time increment [t n , t n+1 ], the time step being defined as ∆t = t n+1t n , for which the solution vector U n V is known at time t n . Then the set of conservation laws (4) written in integral form on the polygon V can be discretized as follows [START_REF] Fridrich | Cell-centered lagangian Lax-Wendroff HLL hybrid method for elastic-plastic flows[END_REF]:

V(U n+α V -U n V ) α∆t = - e∈E(c) l e F n+α-1 2 e • n n+α-1 2 e ( 18
)
where l e is the length of an edge e ∈ E(c), and the time shift index α is either equal to 1/2 for the first step or to 1 for the second step. The Lax-Wendroff scheme defined in its Richtmyer two-step version consists in solving the discrete conservation laws (18) first on the dual mesh at half-time step (α = 1/2), then on the primal mesh at the end of the time step (α = 1).

3.3. Elastic-plastic 'strain' formulation 3.3.1. Lax-Wendroff first step Particularizing Equation ( 18) to a dual cell V p and setting α = 1/2, one gets the updated solution array at time t n + ∆t/2:

V p U n+ 1 2 p = c∈C(p) V pc U n c - ∆t 2 F n c • (l spc-n spc-+ l spc+ n spc+ ) (19)
where the fluxes F n c are evaluated such that F n c = F W n c , W n c refering to the integral average of the array W (8) on the primal cell c at time t n . Notice that for boundary nodes, boundary edges exist in E (p), so that their lengths and normals l e , n pc+ appear in place of separator contributions in Equation [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF]. The fluxes defined on these boundary edges are computed according to prescribed boundary conditions. Ghost cells are used and combined with an elastic Riemann solver to prescribe boundary conditions, and boundary fluxes are computed as the Godunov flux of the Riemann solution.

The unknown vector U p . First, the velocity of a node p is used to compute its displacement at time t n+1 :

u n+1 p = u n p + ∆tv n+ 1 2 p (20)
Second, the updated strain components are used for updating the stresses and internal variables at half-time step. To do so, the strain components, the thermodynamic forces as well as the set of internal variables are approximated in the dual cell V p at time t n by performing a weighted average of their respective values in the primal cells C(p) connected to a node p:

(ε) n p = c∈C(p) V pc ε n c V p ; Y n p = c∈C(p) V pc Y n c V p ; Z n p = c∈C(p) V pc Z n c V p (21)
Then, the set of elastic-plastic constitutive equations can be integrated according to any available integration algorithm [START_REF] Simo | Computational inelasticity[END_REF], so that the update process can be formally written as:

Y n+ 1 2 p , Z n+ 1 2 p = T ((ε) n+ 1 2 p , Y n p , Z n p ) (22)
where T denotes the numerical integrator. In this work, the cutting-plane algorithm introduced in [START_REF] Simo | A unified approach to finite deformation elastoplastic analysis based on the use of hyperelastic constitutive equations[END_REF][START_REF] Ortiz | An analysis of a new class of integration algorithms for elastoplastic constitutive relations[END_REF][START_REF] Simo | Computational inelasticity[END_REF] is used. It consists, when a plastic loading occurs ( f (Y trial n+1/2 ) > 0), in linearizing the yield function f (k) n+1/2 at iteration k about the current values Y (k) n+1/2 . Dropping the index n + 1/2 for the moment, this reads:

f ≈ f (Y (k) ) + ∂ f ∂Y (k) • δY (k) = 0 ( 23
)
where

δY (k) = Y (k+1) -Y (k) .
Then, considering a freezed strain increment at a node p (given ∆ε p ), a plastic relaxation flow is computed consistently with the partition (9), the state laws [START_REF] Ortega | Numerical simulation of elastic-plastic solid mechanics using an Eulerian stretch tensor approach and HLLD Riemann solver[END_REF] and the flow rules [START_REF] Wilkins | Methods in computational physics[END_REF]. The linearization and composition of the latters yield:

δσ (k) = ρ ∂ψ ∂δε e • -δλ (k) ∂Φ ∂σ (k) ( 24
)
δA (k) I = ρ ∂ψ ∂δα I • -δλ (k) ∂Φ ∂A I (k) 1 ≤ I ≤ N (25)
Introducing ( 24)-( 25) into (23) yields, accounting for the linearity of ( 24)-( 25) with respect to (-δλ (k) ), the increment of the plastic multiplier between two iterations:

δλ (k) = f (k) ρ ∂ψ ∂δε e • ∂Φ ∂σ (k) + N I=1 ρ ∂ψ ∂δα I • ∂Φ ∂A I (k) ; λ (k+1) = λ (k) + δλ (k) (26)
Substitution of ( 26) into ( 24)-( 25) gives the updated stress σ (k+1) and the thermodynamic forces

A (k+1)

I

, while its introduction in the flow rules ( 14) yields the updated plastic strain (ε p ) (k+1) and the internal variables α (k+1) I , 1 ≤ I ≤ N. Initial conditions for the above return mapping procedure are provided by the elastic predictor. The iterative process stops when the value of the yield function is lower than a given tolerance.

Lax-Wendroff second step

Let's now write Equation ( 18) on a primal cell c, and set α = 1. The updated solution at time t n+1 is given by the following conservative formula:

U n+1 c = U n c - ∆t V c e∈E(c) (l e F n+1/2 e • n (c) e ) ( 27
)
where l e denotes the length of an edge e ∈ E(c). The fluxes F n+1/2 e defined at these edges are computed by averaging those defined at the end nodes p and p + of an edge e, such that:

F n+1/2 e = F n+1/2 p + F n+1/2 p + 2 (28)
The fluxes F n+1/2 p and F n+1/2 p + are computed from nodal solution arrays defined at time t n + ∆t/2:

F n+1/2 p = F W n+ 1 2 p ; F n+1/2 p + = F W n+ 1 2 p + (29)
where

W n+ 1 2 p
denotes the integral average of the array W (8) on the dual cell V p at half-time step. Next, the set of elastic-plastic constitutive equations is integrated between time t n and time t n+1 in a cell c according to the cutting-plane algorithm [START_REF] Simo | A unified approach to finite deformation elastoplastic analysis based on the use of hyperelastic constitutive equations[END_REF][START_REF] Ortiz | An analysis of a new class of integration algorithms for elastoplastic constitutive relations[END_REF][START_REF] Simo | Computational inelasticity[END_REF]] ( 23)- [START_REF] Toro | Centred TVD schemes for hyperbolic conservation laws[END_REF]. Therefore, the update process can be formally written as:

Y n+1 c , Z n+1 c = T (ε n+1 c , Y n c , Z n c ) (30)
Hence, Equations ( 19) and ( 27) define the Richtmyer two-step algorithm for the elastic-plastic 'strain' formulation, each being supplemented by a projection step of the elastic trial stresses onto the yield locus, ( 22) and ( 30) respectively, to compute the elastic-plastic behaviour.

3.4. Elastic-plastic 'stress' formulation 3.4.1. Quasi-linear form From the system of conservation laws (4), a quasi-linear form is conveniently written as follows [START_REF] Trangenstein | Numerical solution of hyperbolic partial differential equations[END_REF]:

A ∂W ∂t + 3 k=1 B k ∂W ∂x k = 0 ( 31
)
where the auxiliary array W , the matrix A and the k th matrix B k are defined as:

W = v σ ; A = ∂U ∂W = 1 2 0 3 0 3 ∂ε ∂σ ; B k = ∂F k ∂W = - 0 2 1 4 •e k ρ 1 4 • e k 0 4 (32)
where 1 q denotes the identity tensor of order q, F k = F • e k and ∂ε/∂σ ≡ S stands for the elastic-plastic tangent compliance moduli if plastic flow occurs, or simply the elastic compliance tensor for an elastic evolution.

Lax-Wendroff first step

Since the first equation of ( 31) is identical to that of (4), integrating it over the discrete domain V p × t n , t n+ 1 2 yields the same update of the velocity field at a node p at time t n+1/2 as the one given by Equation ( 19):

V p v n+ 1 2 p = c∈C(p) V pc v n c + ∆t 2 l spc ρ σ n c • n spc (33)
where l spc n spc = l spc+ n spc+ + l spc-n spc-. However, the stresses now occur in the second equation of ( 31), thus integrating it over the discrete domain V p × t n , t n+ 1 2 yields:

t n+ 1 2 t n V p S : σ -∇ • v ⊗ 1 + 1 v 2 dVdt = 0 (34)
Considering that the compliance moduli are homogeneous in each cell, and using the divergence theorem for the second term, it comes:

t n+ 1 2 t n c∈C(p) V pc S c : σc dt = ∆t 2 c∈C(p) {spc-∪spc+} v n ⊗ n + n ⊗ v n 2 dl ( 35
)
where have been defined

σ c = 1 V pc V pc σdV ; v n = 1 ∆t 2 t n+ 1 2 t n vdt (36)
Provided the linearized geometrical framework, one gets:

c∈C(p) V pc σ n+ 1 2 p σ n c S c : dσ = ∆t 2 c∈C(p) l spc v n c ⊗ n spc + n spc ⊗ v n c 2 (37)
Equation ( 37) is the two space dimension analog of Equation ( 38) derived for one-dimensional problems in [START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elasticplastic solids[END_REF], which has been shown to be identical to the method of characteristics. If an elastic constitutive response is assumed for the cells C(p), an elastic trial solution can be first computed for the dual cell V p at time t n+1/2 . Assuming a homogeneous medium, the elastic compliance moduli can be factorized in the computation of the left-hand-side of (37), so that the elastic trial stresses read:

V p σ trial p = c∈C(p) V pc σ n c + ∆t 2 l spc C c : v n c ⊗ n spc ( 38
)
where C is the elastic stiffness tensor. Notice that Equation (38) can be deduced from the second equation of [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF], associated with the 'strain' formulation, by means of multiplication by the elastic stiffness tensor.

The yield criterion ( 13) is then computed with the elastic trial stresses (38) in each cell c ∈ C(p). If the yield criterion is satisfied in all cells connected to a node p, then the elastic trial stress state is the actual one, the evolution of the dual cell V p is elastic on the time increment t n , t n+ 1 2 . If the yield criterion is not satisfied in at least one cell of this patch, then an elasticplastic correction is required.

The cells C(p) may have three types of evolution during the time increment t n , t n+ 1 2 . A first subset C E (p) ⊂ C(p) consists of cells having an elastic evolution, for which the elastic trial stresses (38) satisfy the yield criterion. A second subset C EP (p) ⊂ C(p) consists of cells undergoing a transition from elasticity at time t n ((ε p ) n c = 0) to elastoplasticity at time t n+1/2 , for which the yield criterion computed with the elastic trial state is not satisfied at time t n+1/2 while these cells had an elastic evolution at time t n . A third subset C P (p) ⊂ C(p) consists of cells having an elastic-plastic evolution, for which the yield criterion computed with the elastic trial state is not satisfied at time t n+1/2 , while plastic flow occured at time t n . The computation of the left-hand-side of Equation ( 37) is thus performed accordingly, splitting the set of connected cells C(p) into the three subsets C E (p), C EP (p) and C P (p). Further, the integral pertaining to the cells C EP (p) should be split into two parts in order to account for the yield threshold, so that the stress state σ * c leading to the vanishing of the yield function f (σ * c , (A I,1≤I≤N ) n c ) = 0 appears in bounds of integrals. Hence, denoting by S(σ) the elastic-plastic tangent compliance moduli, and by S the elastic compliance moduli, the particularization of Equation (37) to the three above subsets yields:

         c∈C E (p) V pc S c + c∈C EP (p) V pc S c (σ * c ) + c∈C P (p) V pc S c (σ n c )          : σ n+ 1 2 p = c∈C E (p) V pc S c : σ n c + c∈C EP (p) V pc S c (σ * c ) : σ * c - c∈C EP (p)
V pc S c : (σ * c -σ n c )

+ c∈C P (p) V pc S c (σ n c ) : σ n c + ∆t 2 c∈C(p) l spc v n c ⊗ n spc + n spc ⊗ v n c 2 (39)
where the integrals pertaining to plastic flow have been approximated assuming constant elasticplastic tangent compliance moduli over the load increment:

σ n+ 1 2 p σ * c S c (σ) : dσ S c (σ * c ) : σ n+ 1 2 p -σ * c ( 40
)
This approximation is valid if both the plastic strain increment is limited and the normal to the yield function varies little on the time increment t n , t n+ 1 2 . Since the time step ∆t is restricted by the Courant condition, this approximation is reasonable in pratice.

The left-hand-side of Equation (39) appears as some effective compliance moduli S p of the dual cell V p on the time increment, computed by a mixture rule. However, these effective compliance moduli cannot be explicitly inverted. Indeed, the elastic-plastic tangent compliance moduli S(σ) depend explicitly on the unit normal to the yield surface N(σ), whose eigenbasis depends on the plastic flow occurring in each cell connected to a node p, that may be different. Consequently, the computation of the stresses σ n+ 1 2 p has to be performed by solving a small linear system of equations at each node of the mesh, of the form:

S p : σ n+ 1 2 p = B p (41)
where B p stands for the right-hand-side of Equation (39).

3.6. Lax-Wendroff second step 3.6.1. Stress update The integration of the first equation of ( 31) over the discrete domain V c × [t n , t n + 1] gives the update of the velocity field in a cell c at time t n+1 , which can be shown to be identical to the first equation of [START_REF] Richtmyer | A survey of difference methods for non-steady fluid dynamics[END_REF]. Next, the integration of the second equation of ( 31) over the same discrete domain V c × [t n , t n + 1] yields an equation similar to [START_REF] Ortiz | An analysis of a new class of integration algorithms for elastoplastic constitutive relations[END_REF], but with different bounds of integration. Its discretization reads:

σ n+1 c σ n c S c : dσ = ∆t V c e∈E(c) l e         v n+ 1 2 e ⊗ n (c) e + n (c) e ⊗ v n+ 1 2 e 2         (42)
where have been defined

σ c = 1 V c V c σdV ; v n+ 1 2 = 1 ∆t t n+1 t n vdt (43)
and v

n+ 1 2
e is defined at an edge e according to [START_REF] Zwas | On Two Step Lax-Wendroff Methods in Several Dimensions[END_REF], by averaging the velocities defined at the nodes p and p+ at time t n+1/2 :

v n+ 1 2 e = v n+ 1 2 p + v n+ 1 2 p+ 2 (44)
An elastic trial solution can be first computed for a cell c. In that case, Equation (42) reads:

σ trial c = σ n c + ∆t V c C c : e∈E(c) l e v n+ 1 2 e ⊗ n (c) e ( 45
)
and C is the elastic stiffness tensor. Equation (45) can also be deduced from the second equation of formula [START_REF] Richtmyer | A survey of difference methods for non-steady fluid dynamics[END_REF] linked to the Lax-Wendroff second step of the 'strain' formulation, multiplying it by the elastic stiffness tensor. The trial elastic stresses (45) are then tested against the yield criterion (13) in this cell c. If the yield criterion is satisfied, then the elastic trial stresses are the actual ones, the evolution of that cell c is elastic during this step. Conversely, if the yield criterion is not satisfied, an elastic-plastic correction must be performed. The integral appearing in the left-hand-side of Equation ( 42) is thus computed according to the evolution of a cell c on the time increment [t n , t n+1 ]. If it undergoes a transition from elasticity at time t n to elastoplasticity at time t n+1 , the integral is split into two parts in order to account for the yield threshold, so that the stress state σ * c corresponding to the yield surface f (σ * c , (A I,1≤I≤N ) n c ) = 0 appears as an upper and a lower bound of the first and second integrals respectively. Provided that a similar approximation to (40) is done:

σ n+1 c σ * c S(σ) : dσ S(σ * c ) : σ n+1 c -σ * c , (46)
Equation ( 42), premultiplied by the elastic-plastic tangent stiffness tensor C(σ), reads:

σ n+1 c = C c (σ * c ) : S c : σ n c + (1 4 -C c (σ * c ) : S c ) : σ * c + ∆t V c C c (σ * c ) : e∈E(c) l e v n+ 1 2 e ⊗ n (c) e ( 47
)
The above expression can be further simplified in the case of isotropic elasticity and J 2 flow theory. The elastic-plastic tangent stiffness tensor C(σ) can be shown to be of the form [3, Equation (2.3.9)]:

C(σ) = C -βN ⊗ N ( 48
)
where β is a function of the elastic and hardening moduli, as well as hardening variables for certain constitutive models, and N is the unit normal to the yield surface. Introducing (48) into (47) gives the stress update for a cell undergoing an elastic-elastoplastic transition:

σ n+1 c = σ n c + β 2µ N * ⊗ N * : (σ * c -σ n c ) + ∆t V c C c (σ * c ) : e∈E(c) l e v n+ 1 2 e ⊗ n (c) e ( 49
)
where µ is the second Lamé parameter. Conversely, if the stress state in a cell c already lied on the yield surface at time t n , then this cell undergoes an elastic-plastic evolution on the time increment [t n , t n+1 ]. Provided an approximation similar to (46), Equation (42) gives:

σ n+1 c = σ n c + ∆t V c C c (σ n c ) : e∈E(c) l e v n+ 1 2 e ⊗ n (c) e (50)

Update of internal variables

The update of the plastic strains and additional internal variables appears for the 'stress' formulation as a stress-driven process, while the integration algorithms used for elastic-plastic constitutive equations are generally strain-driven [START_REF] Simo | Computational inelasticity[END_REF], provided an increment of strain. The procedure thus amounts to adapt the size of the elastic convex and the location of its center (if both isotropic and kinematic hardenings are accounted for) so that the updated stress state ((49) or (50)) lies on its boundary, consistently with the plastic flow rule and hardening laws [START_REF] Wilkins | Methods in computational physics[END_REF]. Since elastic-plastic equations are nonlinear, the solution process is performed iteratively. One general way to proceed consists in modifying the cutting-plane algorithm [START_REF] Simo | A unified approach to finite deformation elastoplastic analysis based on the use of hyperelastic constitutive equations[END_REF] so that a stress-driven update of internal variables be performed. Thus the yield function f (k) is linearized at iteration k about current values of the sole additional thermodynamic forces A I 1 ≤ I ≤ N, the stresses being given:

f ≈ f (σ n+1 , A (k) I ) + N I=1 ∂ f ∂A I (k) • δA (k) I = 0 ( 51
)
where δA (k)

I = A (k+1) I -A (k) I .
Then introducing Equations ( 25) into (51) and following the same procedure than that of Section 3.3.1 yields the increment of the plastic multiplier between two iterations:

δλ (k) = f (k) N I=1 ρ ∂ψ ∂δα I • ∂Φ ∂A I (k) ; λ (k+1) = λ (k) + δλ (k) (52)
The iterative process stops when the value of the yield function is lower than a given tolerance.

Notice that for the two-dimensional plane strain case, the out-of-plane stress component σ 33 has to be updated together with the internal variables. Therefore, the yield function ( 51) is also linearized with respect to that component. The increment of the plastic multiplier thus reads as [START_REF] Toro | Centred TVD schemes for hyperbolic conservation laws[END_REF] with the sole out-of-plane component of (24).

Applications

Plane waves with Riemann-type initial conditions

Let's consider a rectangular domain Ω : (x, y) ∈]0, L[×]0, H[, made of an elastic-plastic material. Traction free conditions are set at the left and right sides, while zero shear stress and normal velocity components are prescribed on the top and bottom sides, so that combined with the plane strain assumption (zero out-of-plane velocity and strain components), planes waves can be simulated. Riemann-type initial conditions are prescribed: the longitudinal component of the velocity is prescribed to -v in the first half of the medium x ∈ [0, L/2[, and to v in the second half x ∈]L/2, L], while the stresses are considered to be zero everywhere initially. The prescribed velocity is set so that plastic flow occurs:

v = 2 Y H ρc L (53)
where Y H = (λ + 2µ)σ y /2µ denotes the Hugoniot elastic limit, c L is the speed of the elastic pressure waves c L = (λ + 2µ)/ρ, and λ, µ and ρ denote the Lamé parameters and the mass density respectively. Provided linear isotropic or kinematic hardenings, an analytical solution of that problem has been developed in [START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elasticplastic solids[END_REF]. This test case allows to generate a compressive plastic reloading after an initial tensile one, by means of wave interactions and reflexions at free ends. An elastic-plastic material with linear isotropic hardening is considered, associated with the von Mises yield criterion:

f = σ eq -C p -σ y ; σ eq = 3 2 s : s ; C = {σ| f (σ) ≤ 0} ( 54
)
where s is the deviatoric part of the Cauchy stress σ, C is the hardening modulus and p stands for the cumulated plastic strain. Table 1 summarizes the material properties and geometrical parameters used for the simulations. Comparison is here performed between (i) the Lax-Wendroff

Material parameters

Geometrical parameters

E = 200GPa C = 10GPa L = 6m ν = 0.3 σ y = 400MPa H = 0.08m ρ = 7800kg.m -3
Table 1: Parameters solutions obtained with the 'stress' and 'strain' formulations, (ii) a 'projected' solution that consists of a Lax-Wendroff elastic computation (the first step is computed with (38), the second one with (45)) followed by a projection onto the yield locus by means of a radial return [START_REF] Wilkins | Methods in computational physics[END_REF][START_REF] Simo | Computational inelasticity[END_REF], (iii) a finite element solution obtained with the finite element code Cast3M [START_REF]Cast3M, user's manual[END_REF] and computed with an explicit central difference time integrator, Q1 elements and a lumped mass matrix [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF], and (iv) the analytical solution. Computations are performed with a mesh made of 200 × 3 quadrangular grid cells and a Courant number set at 0.9. Figures 2 and3 show for two different times the superposed plots of the longitudinal stress and plastic strain components computed with these different approaches. The solution first consists of two discontinuous elastic predictors travelling leftward and rightward (see Figure 2), captured in few cells by all the methods due to the chosen Courant number. The subsequent discontinuous plastic waves need much more cells to be captured. The numerical solutions obtained with the finite elements and with the Lax-Wendroff 'stress' and 'strain' formulations give results close to each other for the longitudinal stress. They allow to slightly improve the capture of the discontinuous plastic wave enabled by an a posteriori projection of a Lax-Wendroff elastic solution onto the yield locus, although the finite element solution shows oscillations on the stress plateau. However, all the methods show a non-physical peak of plastic strain that appears at the middle of the medium. This non-physical peak can be explained by the combination of the initial discontinuous velocity profile and the fact that neither Lax-Wendroff nor FEM are monotone or Total Variation Diminishing [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] (TVD) schemes. Indeed, the use of a TVD scheme in [START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elasticplastic solids[END_REF] on that example allowed to remove that peak.

The finite element and Lax-Wendroff 'strain' solutions slightly overestimate the plastic strain plateau, while the Lax-Wendroff 'stress' one underestimates it. The 'projected' numerical solution seems closer to that plateau on this example, although its plastic waves are less sharp than these obtained by the other methods. After reflexion of waves at the two ends, their interaction yields a compressive plastic reloading as shown in Figure 3, in which the same trend as depicted above is observed. The Lax-Wendroff 'stress' formulation gives slightly less accurate plastic strains on the compressive plastic plateau than the other methods. Results obtained from computations performed on coarser (100×3 grid cells) and finer meshes (500×3 grid cells) are respectively shown in Figures 4(a) and 4(b). Coarse and underresolved grids do not show any significant difference between a double projection onto the yield locus and a single one. However as the mesh is refined, a clear and net difference between these two procedures is observed in Figure 4(b): fewer cells are required for the Lax-Wendroff 'strain' and 'stress' solutions to capture the plastic discontinuity.

Partial impact on a plane

Let's now consider the bidimensional square domain shown in Figure 5, submitted to an impact on a part of its top face, by means of a step function of a pressure p. A symmetry condition is considered on the left side, and free boundaries are set at the bottom, right, and remaining part of top faces. This problem is treated in the bidimensional plane strain framework. For this test case, a comparison between the Lax-Wendroff numerical solutions and a Q1-finite element one is considered. The finite element numerical solution is obtained with the code Cast3M [START_REF]Cast3M, user's manual[END_REF], computed with an explicit central difference time integrator, and uses a lumped mass matrix. Computations are performed with a mesh that consists of 100 × 100 regular quadrangular grid cells, and a Courant number set at 0.9. An elastic-plastic constitutive model with a linear isotropic hardening is also considered for this example with parameters listed in Table 1, to which are added these 2. Figures 6 and 7 show the comparison between the finite element and the Lax-Wendroff numerical solutions at two different times. These figures consist on the one hand of plots of numerical isovalues, in particular the normal stress σ 22 and the cumulated plastic strain ε eq ≡ p, and on the other hand of superposed plots of these fields along the symmetry line (left side). At time t = 1.88 × 10 -4 seconds (see Figure 6), the pressure and shear waves generated by the step pressure prescribed on a part of the top face propagate downward. An elastic precursor is followed by a discontinuous plastic wave, which propagate plastic strains. As can be seen in Figures 6(a) and 6(c), the Lax-Wendroff 'strain' and 'stress' formulations give numerical results close to these provided by finite elements, though the latter shows much spurious numerical oscillations. In particular, the three former solutions allow to capture well the discontinuous plastic wave, better than the Lax-Wendroff projected solution. This is also seen on the plastic strain front (see Figure 6(c)), though the overall isovalues shown in Figure 6(b) are comparable for all numerical solutions.

Geometry Loading

a = 0.5m p = 2Y H L = 2m
Once the initial waves have done several round trips (see Figure 7), the plastic strains do not evolve anymore in the structure. These have been much spread in the height of the domain, and the four numerical solutions look very close to each other, and should yield almost the same residual state.

Dynamic ratchetting

Let's consider an infinite bidimensional medium, containing a set of holes aligned vertically and equally spaced, as shown in Figure 8. This infinite medium is submitted to given oscillating pressure pulses propagating horizontally. This problem is reduced according to the symmetries of the system so that the computational domain shown in Figure 9 with simplified boundary conditions is considered for the analysis. An elastic-plastic material with a nonlinear kinematic hardening is considered, with the von Mises yield criterion: where γ denotes an additional hardening parameter. The computational domain is submitted on its right side to non-symmetric alternating tension/compression tractions:

f = σ eq (ξ) -σ y ; σ eq = 3 2 ξ : ξ ; ξ = s -X ; C = {(σ, X)| f ≤ 0} (55)
σ ∞ (t) = σ a sin πt 2T + σ 0 ( 57
)
where σ a , σ 0 and T refer respectively to the stress magnitude, the stress offset and a period defined as T = L/c L corresponding to the time required for an elastic pressure wave to travel the length of the domain. Provided these non-symmetric alternating stress pulses, the nonlinear kinematic hardening law is expected to produce plastic strains growing progressively with constant pitch, especially close to the hole due to the geometry of the domain considered. This phenomenon is known as ratchetting [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF], and since it is triggered here by the propagation of waves, we call it dynamic ratchetting. The present test case does not involve any discontinuous solution, but rather aims at showing the ability of Lax-Wendroff schemes to simulate unconventional dynamic plastic flows with 'refined' elastic-plastic constitutive models. Computations are performed in the bidimensional plane strain framework, and comparison is made between the Lax-Wendroff numerical solutions and an explicit Q1-finite element solution obtained with the code Cast3M [START_REF]Cast3M, user's manual[END_REF] as for the previous example. Computations are performed with a mesh that consists of 4200 non-rectangular quadrangular grid cells, and a Courant number set at 0.55. The parameters used for the geometry and the loading are summarized in Table 3, while these pertaining to the material are the same as these listed in Table 1, to which should be added the value of the additional hardening parameter γ = 40 of the Armstrong-Frederick law. Figure 10 shows some comparisons between the finite element and the Lax-Wendroff numerical solutions obtained with the 'strain' and 'stress' formulations at time t = 1.25 × 10 -3 seconds, that is after two cycles of loading/unloading close to the hole area. It can be observed that the isovalues of the longitudinal stress look very close for all numerical solutions, and predict it reaches a peak value at that time close to the top point of the hole. The computed cumulated plastic strain looks also quite close in the three numerical simulations, but shows however small differences close to the top point of the hole. This is better observed in Figure 11 where, at that point, plots of the time evolution of the plastic strain component ε p 11 , and of the stress component σ 11 as a function of this plastic strain component, are made. On the one hand it can be observed that the simulations conducted with the above modeling actually predict the occurence of a ratchet phenomenon as expected. On the other hand, slight differences in the computation of the ratchet appear between the finite element and the Lax-Wendroff numerical solutions. First, the ratchet pitch computed at that point by the finite elements appears slightly greater than these obtained by the two formulations of Lax-Wendroff, which also give results close to each other. Second, stresses computed with Lax-Wendroff are slightly larger than these obtained by finite elements.

Geometry Loading

r = 0.5m σ a = 3Y H L = 1m σ 0 = Y H /2
However, if we look a little bit away from that corner point as shown in Figure 12, the agreement between numerical solutions looks better, the ratchet pitches computed by the finite elements and the Lax-Wendroff 'stress' formulation perfectly match, while that computed by the 'strain' formulation is slightly lower. These differences may be attributed to the different ways of prescribing boundary conditions between the finite element and Lax-Wendroff schemes. However, the results obtained show the ability of the Lax-Wendroff scheme to compute complex plastic flows in dynamic elastic-plastic numerical simulations.

Conclusion

Two formulations were proposed in this work for the extension of the well-known Lax-Wendroff scheme to elastic-plastic solids in the linearized geometrical framework. These formulations are referred to as 'strain' and 'stress' ones since the strain or the stress components are respectively included in the unknown vector in addition to the velocity ones. The Lax-Wendroff scheme was here implemented in its Richtmyer two-step version so that the projection of the 20 elastic trial stresses onto the yield locus is performed twice per time step. It allows to reduce the number of cells required to properly capture discontinuous plastic waves with respect to a classical one-step strategy consisting in an a posteriori projection of the elastic trial stresses onto the yield locus. Besides, it has been shown on bidimensional plane strain examples the ability of these formulations to compute complex plastic flows like dynamic ratchetting in dynamic elastic-plastic numerical simulations.

Figure 1 :

 1 Figure 1: The dual control volumes for interior and boundary nodes are shown in red background, and the edges e ∈ E (p) connected to a node p in red. The normals to the dual cell boundaries are in blue, while the normals to primal edges are in black. The reader is refered to the web version of this article for intepretation of colours.

  at time 0.000437099 sec.

Figure 2 :

 2 Figure 2: Comparison at time 4.37 × 10 -4 seconds of the Lax-Wendroff solutions ('stress' and 'strain' formulations), an elastic Lax-Wendroff solution projected onto the yield locus, a finite element solution and the analytical stress and plastic strain fields.

Figure 3 :

 3 Figure 3: Comparison at time 1.38 × 10 -3 seconds of the Lax-Wendroff solutions ('stress' and 'strain' formulations), an elastic Lax-Wendroff solution projected onto the yield locus, a finite element solution and the analytical stress and plastic strain fields.

Figure 4 :

 4 Figure 4: Comparison at time 4.37 × 10 -4 seconds of the Lax-Wendroff solutions ('stress' and 'strain' formulations), an elastic Lax-Wendroff solution projected onto the yield locus and the analytical for coarser and finer meshes.

2 aFigure 5 :

 25 Figure 5: Partial impact on a plane

Figure 6 :

 6 Figure 6: Comparison at time t = 1.88×10 -4 seconds of the Lax-Wendroff solutions ('stress' and 'strain' formulations), an elastic Lax-Wendroff solution projected onto the yield locus and a finite element solution.

  Plots along the symmetry line

Figure 7 :Figure 8 :Figure 9 :

 789 Figure 7: Comparison at time t = 5 × 10 -3 seconds of the Lax-Wendroff solutions ('stress' and 'strain' formulations), an elastic Lax-Wendroff solution projected onto the yield locus and a finite element solution.

  (a) Isovalues of the stress component σ 11 (b) Isovalues of the cumulated plastic strain p

Figure 10 :

 10 Figure 10: Comparison at time t = 1.25×10 -3 seconds of the Lax-Wendroff solutions ('stress' and 'strain' formulations) and a finite element solution.

Figure 11 :

 11 Figure 11: Plots at the top point of the hole (x 1 = 0, x 2 = r) of the time evolution of the plastic strain component ε p 11 (at the left), and that of the stress component σ 11 as a function of this plastic strain component (at the right).

Figure 12 :

 12 Figure 12: Plots at the location (x 1 = r/5, x 2 = r) of the time evolution of the plastic strain component ε p 11 (at the left), and that of the stress component σ 11 as a function of this plastic strain component (at the right).

Table 2 :

 2 Parameters of the partial impact test case

Table 3 :

 3 Parameters of the dynamic ratchetting test case