V Ripard

D Goncalves

F Ville

J H O Seabra

J Cavoret

P Charles

Orcid Ids

Grease composition influence on friction & starvation

Keywords: Tribology, grease, lubrication, film thickness, starvation

Nowadays, grease lubrication is frequently used in rolling element applications such as bearings or constant velocity joints. The advantage of grease is to supply lubrication to the application without leaking thanks to its consistency. Nevertheless, starvation can occur leading to damage such as scuffing. In the present study, starvation is analyzed using the Starvation Degree parameter through tests with different oper-ating conditions and different types of greases.

Introduction

Greases are mainly used to lubricate components such as rolling element bearings or drivetrain components. The lubrication mechanism is quite complex. [START_REF] Lugt | A review on grease lubrication in rolling bearings[END_REF] It seems that grease lubrication is not a continuous process but it is characterized by events caused by film breakdown and recovery. [START_REF] Lugt | On the chaotic behavior of grease lubrication in rolling bearings[END_REF] In order to study the lubrication mechanism and possible starvation, film thickness is measured using different techniques such as electrical conductivity [START_REF] Muennich | Elastohydrodynamic lubrication of grease-lubricated rolling bearings[END_REF] or interferometry [START_REF] Gohar R | Optical measurement of Oil film thickness under elasto-hydrodynamic lubrication[END_REF] on a ball on disc tribometer for example. Using these techniques, base oil viscosity appears not to be enough to evaluate film thickness. The thickener must be taken into account especially at low speed. [START_REF] Morales-Espejel | Film thickness in grease lubricated slow rotating rolling bearings[END_REF] The thickener structure in particular plays a key factor on film thickness. [START_REF] Kaneta | Effects of a thickener structure on grease elastohydrodynamic lubrication films[END_REF] Relubrication is a key parameter to extend bearing life at low temperature. In addition, the mechanical design of the component also plays an important role. [START_REF] Cann | Bearing performance limits with grease lubrication: the interaction of bearing design, operating conditions and grease properties[END_REF] If relubrication is poor, starvation could occur. To evaluate starvation risk, a single dimensionless parameter (SD), between the operating parameters and the transition from the fully flooded to starved regime in case of oil lubricating contact was developed. [START_REF] Cann | The transition between fully flooded and starved regimes in EHL[END_REF] Its parameters are surface tension effects σ S , amount of oil present in the vicinity of the track h oil∞ , base oil viscosity η 0 , the track width a and the velocity u. It gives the formula below:

SD = η 0 ua h oil∞ σ S
For greased lubricated contacts, the SD evaluation is quite complicated due to events as mentioned before. [START_REF] Gonçalves | An experimental study on starved grease lubricated contacts[END_REF] Nevertheless, film recovery was often observed even if starvation occurs. [START_REF] Nagata | Track replenishment by lateral vibrations in grease-lubricated EHD contacts[END_REF] It seems to be due to contact oscillations that help film thickness recovery up to fully flooded conditions.

In the present paper, different greases in terms of composition and physical or chemical properties were tested using different tribometers. Starvation and the SD parameter are analysed using different operating conditions such as temperatures and/or speeds.

Material and methods

Grease properties

Four greases were used. Their compositions vary (Table 1). Two base oils are used. The first one (A) is a mix between PAO and mineral oil. The second one (B) is a synthetic one. Two different thickener technologies are used. The first one is a complex soap using Lithium and Calcium. The second one uses polyurea thickener.

However, even if G3 & G4 have same thickener technology, the manufacturing process changes. To finish, two types of additives [START_REF] Dr | The Effect of Friction Modifier additives on CVJ Grease Performance[END_REF] are used : MoDTC (Molybdenum Dithiocarbamates), MoDTP (Molybdenum Dithiophosphate). Concentrations were normalized to 1 for the maximum additive concentration.

Moreover, the oil separation from lubricating grease is determined on greases using IP121 method. [START_REF] Afnor | Produits pétroliers et graisses lubrifiantes -Séparation d'huile au stockage des graisses lubrifiantes -Méthode sous pression -Conditions statiques[END_REF] For 42 h at 40°C, greases are inside cups in order to separate the oil from the grease. A mass is used in order to press the grease into cup. The oil is salvaged in a container.

As expected, G1 and G2 give the same oil separation results because only additives were added between G1 & G2. On the contrary, it is possible to observe a huge gap between G3 and G4. Even if the components are the same, different oil separation results are obtained which will influence their properties. It is due to the completely different manufacturing process of the thickener.

High frequency reciprocating rig

A PCS instrument HFRR is used to study the different available greases. The HFRR allows to measure friction coefficient and separation rate of a lubricant. The test principle is to study an alternative contact with pure sliding of a ball on a disc (Figure 1). It allows also to evaluate the grease ability to supply oil in the contact. It is easy to use thanks to its design. Also balls & discs are inexpensive sample. The applied load can be set from 0.98 N -(0.63 GPa) to 9.81 N (1.37 GPa) which applied a pressure named P h-HFRR .

The upper specimen is a 6 mm diameter ball and the bottom specimen a flat disc. They are made from polished 100Cr6 steel with Ra < 20 nm for the disc and Ra < 50 nm for the ball. The rig can be configured in terms of frequency f HFRR (from 10 Hz to 200 Hz) and stroke d HFRR (from 20 μm to 2 mm). The temperature can be set from ambient temperature up to 150°C thanks to a heater block placed under the lower specimen.

Test conditions, chosen in order to maximize speed on HFRR are the following:

f HFRR = 15Hz ; d HFRR = 1mm ; P h-HFRR = 1.32 GPa
To conduct this study different test durations and operating temperatures are set (Table 2). For this study, is it possible to write SD as a function of SD HFRR which is assumed constant for this study ( a h oil∞ σ S ), η 0 & u which variate between tests:

SD = η 0 ua h oil∞ σ S = η 0 uSD HFRR
This condition gives a mean sliding velocity of 30 mm/ s. A test example is available in Figure 2. For test condition 1 with G4 it is possible to observe repeatability between 2 tests (curves G4-1 & G4-2) and also the film thickness corresponding to first friction test Considering a time t, it is possible to calculate the gap between the friction coefficient G4-1 and G4-2 and divide it by G4-1 value. Taking the mean value, the deviation between these curves is 2%. So, the HFRR test device is a very reliable rig in order to study the friction coefficient with grease.

In addition, to analyse the behaviour of grease during the test, the electrical contact potential is measured to evaluate surface separation.

A second set of tests are considered in Table 3 in order to evaluate starvation. In order to increase the SD parameter, it is possible to increase mean velocity on test rig.

Considering, the minimum SD value as SD G1-3 ,i ti s possible to calculate variation of SD for each condition. It is expose in Table 4 as a ratio between SD SD G1-3 .I ti s important to note that a that a higher SD value is more prone to starvation.

The wear scars induced are measured using a Sensofar PLu neox optical profiler with Nikon objectives. Its vertical resolution is less than 0.1 nm. It uses an interferometry technique in order to determine topography.

Ball-on-disc film thickness

In order to explain the different behaviour of the grease in terms of wear and to complete the rheology study an investigation is carried out using EHD2 from PCS Instruments. It consists of a ball or roller on a disc contact with rolling and/or sliding conditions. This equipment allows to validate first results from HFRR and measure film thickness instead of having film separation. It is also possible to visualize contact and lubricant repartition inside the contact with a camera. The ball and disc are mounted on an independent rotational shaft with speed control S EHD . So, a slide to roll ratio SRR EHD can be introduced inside the contact. It can measure traction coefficient but also the film thickness down to 1 nm. The main goal of that study is to measure the grease film thickness using interferometry. [START_REF] Gohar R | Optical measurement of Oil film thickness under elasto-hydrodynamic lubrication[END_REF] Samples used are:

• A roller made from 100Cr6 with a minor diameter of 19.05 mm and a crowned diameter of 133.35 mm. Its roughness is down to 20 nm. • A disc made from glass (E = 75 GPa v = 0.2) with a roughness of approximately 5 nm. Tests are conducted without grease scoop 9 to go as close as possible to the real contact. To insure repeatability, a 3D printed stencil is used (Figure 3(a)). When the stencil is filled, it allows 7.2 millilitre of grease to be placed on the disc (Figure 3(b)) ensuring that the same volume of grease is used for all the tests.

The tests conditions are available in Table 5. They are set in order to allow maximum load available (so maximum Contact Hertz Pressure P h-EHD ) and a mean velocity equal to conditions 5 & 7.

Compared to HFRR, the speed is almost equivalent to tests 5 & 7. However, EHD2 uses rolling and sliding contact as opposed to HFRR that only uses pure sliding contact. Another difference, is that, even if the load is almost equivalent, one sample is made of glass so the Young's modulus varies a lot. For that reason, hertz pressure drops.

Experimental results

HFRR

2.1.1. Friction at 40°C. The first step is to use HFRR in order to qualify the grease in terms of the friction coefficient and surface separation.

At 40°C, the first observation is the role of additive concentrations. Indeed, the modification of the MODTC concentration and also the addition of MODTP (G1 to G2) allow the friction coefficient to be decreased by 25%, once COF is stabilized (around 16 h). Moreover, the replacement of the thickener and base oil (G2 to G3& G4) allows the friction coefficient (COF) again to be decreased by 30% in order to give a low friction grease. These 2 major modifications give two types of grease: a low friction grease with a COF of 0.06 and a regular grease with a COF of 0.12 after stabilization. Figure 5 gives an indication of the separation during tests measuring the electrical resistance between the 2 surfaces.

Separations for greases G2, G3& G4 vary during tests at the opposite of G1 which is stable from the beginning. The most important is to do not have separation equal to 0 which could mean a dry contact. In this case, separation rate is always above 20%. After the stabilization of the COF, a distinction can be established between the grease with a full additive package (G2, G3, G4) and grease G1 without the same MODTC concentration and MODTP.

In addition, wear scars have been measured and the polyurea grease seems to better protect surfaces.

2.1.2. Friction at 80°C. In Figure 6, the coefficient of friction evolution is represented. For each grease, 2 different durations were set: 3 h and 18 h. The tests were carried out twice and the curves could be superposed. In order to clarify, only 18 h tests were plotted on Figure 6.

The first observation is that the coefficient of friction of greases G1 & G2 converge around 0.08 to 0.09. However, the G1 friction coefficient seems unstable. An explanation can be that the tribofilm is heterogeneous. Using XPS measurements, [START_REF] Ripard | Tribological characterization of greased drive-shaft : Evaluation of constant velocity joint durability[END_REF] it has been proved that Mo intensity can vary all along the scar width depending on the grease composition Figure 7.

At 80°C, the additive concentration does not have a significant influence on the coefficient of friction, when comparing grease G1 and G2. The low friction greases G3 and G4 have a coefficient of friction below 0.08 which is interesting for a drivetrain application. As an example, a low friction grease will lead to reduce generated axial force on a constant velocity joint. [START_REF] Serveto | Joint tripode coulissant de transmission automobile. Effort axial généré : essais et modèles[END_REF] In addition, a friction coefficient increase (G4) is noticed at around 3 h. This increase can be linked to the thickener on G4, the bleed oil from G3 & G4 being the same.

Figure 8 presents the separations of each grease at 80°C which seem to be quite significant. At 80°C the performance gap between greases becomes less. The concentration of additives seems to be less relevant than at 40°C. However, it is important to notice the couple oil + thickener plays a key role on performance at this temperature.

Test operating conditions evolution to provoke starvation.

A first approach is proposed to explain the influence of grease composition on friction. It shows the important role of additives on friction coefficient but also the contribution of the couple thickener + oil. This role is amplified at 40°C as it is possible to observe. Also, the couple oil + thickener is involved in friction coefficient establishment duration.

In order to induce starvation, test conditions available in Table 3 were used on HFRR. These conditions, increasing frequency from 15 to 100 Hz and distance from 1 to 2 mm, induce a mean sliding speed thirteen times higher which can lead to starvation. The first step was to study greases at 40°C using condition number 5. Indeed, 40°C seems to be important to understand the grease behaviour. The friction coefficients are available in Figure 9.

Greases G1, G2 & G4 have exactly the same behaviour as in the previous conditions of Table 2 (Figure 4). Only a variation of friction coefficient between test conditions #1 and #5 is noticed due to speed. In addition, the surface separation is very satisfactory (Figure 10).

However, for grease G3 a major friction coefficient is obtained. Moreover, when this jump occurs, the separation goes to 0%, meaning, there is no lubricant between contacting parts. Starvation seems to occur in that case.

In order to go further, another test at 80 Hz is carried out (test number 6) which reduces the speed and so the SD parameter (Table 4). Also this test is repeated at 100 Hz. (see Figure 11).

This phenomenon appears to be repeatable. A first jump appears, then the lubrication resumes and a second jump occurs. After these failures and due to surface degradation, the grease G3 is now the equivalent of a standard grease as G1. Indeed, it gives a friction coefficient equivalent or superior to G1 due to surface degradation. Also, it is frequency dependent. At 80 Hz, the phenomenon does not occur (Figure 11). Looking at it with the naked eye (figure 12), this major friction seems to involve a lot of particles inside the grease. The grease becomes black between the samples. Now, at 80°C the same tests are performed (Figures 13 and14).

There is no starvation compared to Figures 9 and10. Indeed, here the friction coefficients are almost stable (Figure 13) and the film thickness is always above 50% (Figure 14). Finally, starvation occurs only at 40°C with grease G3.

To explain this wear (Figure 15) a test with condition 5 and G3 is repeated and stopped when the jump in friction coefficient occurs.

Wear scars were measured on the ball and disc after testing with grease G3 as shown in Figures 16 &17 It is possible to observe a material transfer from the disc (Figure 16) to the ball (Figure 17), which shows that, scuffing occurs with grease G3. The disc and ball weld and, with the movement imposed by the actuator, it transfers material onto the ball surface from the disc surface. Now, it could be interesting to study the influence of the grease composition directly on film thickness.

Film thickness study

On an EHL test rig, grease test repeatability is very difficult to obtain due to the chaotic behaviour of the grease 15- 17 even for exactly the same operating conditions and grease volumes. For conditions 7&8, the tests are repeated more than twice to obtain meaningful results. In this section, G1, G3 and G4 are studied as they have meaningful HFRR results.

In order to explain scuffing, the film thickness is measured as a function of time. The condition (8&9) seems relevant compared to previous tests (Table 2) due to velocity which seems to be the key parameter in this study. It can be compared to conditions 5&7 from table 3 of the HFRR section. The results for G1, G3 and G4 are available in Figures 18, 19 First of all, it is interesting to note that for a same test a grease can have some chaotic but at the same time repeatable film thickness behaviour over time. For example, grease G1 at 40°C has 3 distinct behaviours (Figure 18). The first one separates correctly the surface around 250 nm. The second one is the opposite. The film thickness becomes smaller and approaches zero damaging the spaced layer coating of the test track on the disc and another one which oscillates between the 2 curves. This behaviour is very interesting. Indeed, it shows the ability of grease to replenish contact. The same ability is observable for G4 not only at 40°C but also at 80°C.

In order to investigate phenomenon from Figures 18 to 20, interferometric pictures of the contact area have been done. In Figure 21, images from EHD illustrate the contact state before and after replenishment. The rolling direction is indicated by a yellow arrow. When starvation begins in the left of top Figure 21, it could be observed starvation induced vibrations which promoted the replenishment (the sudden rise in film thickness of Figures 18-20).

However, this behaviour is not observed on G3. At 40°C , the film thickness decreases quickly and never recover unlike G4&G1. Finally, after some minutes, it goes to zero and damages the track. It is illustrated in Figure 22.

It is possible to observe the beginning of starvation on the right and left of the contact. However, the grease never replenishes contact and tends to a dry contact. It is similar to the behaviour found by other authors. [START_REF] Huang | Film thickness decay and replenishment in point contact lubricated with different greases: a study into oil bleeding and the evolution of lubricant reservoir[END_REF] On the contrary, at 80°C, the replenishment of this grease G3 greatly improves (temperature drops oil viscosity which improves the oil separation) and the film thickness at 80°C can even be higher than at 40°C, as shown in Figure 19.

Discussions

Grease G1, G2 & G4 were tested by 2 experimental devices. These devices which are complementary allow to study greases having friction coefficient, separation but also film thickness variation information.

It has been shown that G1 & G4 which have completely different composition in terms of additives concentration and couple oil + thickener have completely different friction behaviour. This composition variation which lead to double friction coefficient from G4 to G1, have no influence on scuffing.

Indeed, for G1 & G4, when the contact is at low values of film thickness, the grease can replenish the contact thanks to its properties even against centrifugal effects. For G1 at 80°C, the grease does not show the same ability. Nevertheless, the film thickness is lower than previous values due to the increase of temperature but it is more stable during the test, and the surface separation is improved protecting the surfaces too.

There are sporadic replenishments during tests. It is possible to explain:

• Starvation induces vibrations: these vibrations will allow the grease to move back inside the contact area • Centrifugal effect: during tests, oil is stored on both sides of the contact. Due to centrifugal effect, the storage closer to the rotational shaft can supply the contact with the bleed oil. • Mechanical degradation of the thickener matrix: degraded grease thickener insures a residual film on the rolling track.

All these phenomena improve re-flow until the next starvation event occurs. [START_REF] Cann | Starvation and reflow in a grease-lubricated elastohydrodynamic contact[END_REF][START_REF] Cann | Thin-film grease lubrication[END_REF] On the contrary, although G3 & G4 have close composition, same SD, and performance except at 40°C with condition #5, it has been noticed that G3 at 40°C induces starvation in HFRR. This starvation leads to scuffing and so surface degradation. On EHD2 with G3 at 40°C, the film thickness events do not happen and therefore, the film generated shows a much smaller thickness, operating close to dry contact and hence, damaging the surface which can be the origin of scuffing found in section 2.1.3. However, it is interesting to note that G3 and G4 have the same oil (B) so the same viscosity. Also, in that case, the contact area and velocity are the same. To finish, the film thickness after 500 s at 40°C is almost the same. So, the parameters from starvation degree SD in oil-lubricated contacts is constant (table 4). This starvation criteria is not valid for grease as suggested by authors8. So, it could be interesting to consider oil separation measurements to study greased contact starvation as these are highly different for G3 & G4. Nevertheless, for these last greases the compositions and chemical and physical properties are similar, only the manufacturing process is different.

Conclusions

In The present study, it has been proved that grease composition and in particular additives can decrease friction. Adding MODTC and MOTDP in enough quantity can improve friction coefficient by 25% at 40°C. Also, changing the couple oil-thickener from a lithium-calcium complex with a mineral & PAO base oil to polyurea with Synthetic ester + PAO base oil allows friction to be decreased by 30% at 40°C. However, if the contact frequency increases it is possible to observe starvation and so scuffing on surfaces using greases with the same oil, thickener technology and additives. It is important to note that the manufacturing process of polyurea thickener seems to play a key role on starvation on our case.

Borrowing the SD parameter from oil lubricated contact theory, it is possible to conclude that it is not sufficient to predict starvation on grease lubrication. Indeed, in this study, 2 greases with same oil viscosity, contact parameters leading to equivalent SD parameters give different behaviours due to very different oil separation. Indeed, when one grease insures correct lubrication, the other one induces starvation and so scuffing according to HFRR and EHD2 measurements.

In a future study, it could be interesting in order to predict starvation to modify the SD parameter with a physical notion of oil bleed using the IP121 method for example. To do this, it would be interesting to lead this future study by controlling more finely the parameters between the different greases. Also, it has been shown that temperature and contact frequency play a key role on starvation.

An application of this study can be rolling element lubrication from automotive drivetrains. A constant velocity joint uses grease in order to lubricate the contact between the roller and housing. Drivetrain greases are important for automotive transverse drivetrains. It fills 2 main goals: lubricating the contact and extracting heat. In that case, the friction coefficient can be directly linked to the generated axial force of the tripod [START_REF] Serveto | Joint tripode coulissant de transmission automobile. Effort axial généré : essais et modèles[END_REF] . Friction between rollers and housing creates a force. The strength will depend on the friction coefficient and implies vibration on the car. These functions which are essential, insure the components have a long life 21 avoiding starvation for example.

Figure 1 .

 1 Figure 1. Schematic of the HFRR.

Figure 3 .

 3 Figure 3. (a) EHD stencil, (b) stencil filled.

Figure 4 .

 4 Figure 4. Grease friction coefficient results at 40°C.

Figure 5 .

 5 Figure 5. Greases film thickness at 40°C.

Figure 6 .

 6 Figure 6. Grease friction coefficient results at 80°C.

Figure 7 .

 7 Figure 7. XPS line scan from wear scars with same types of grease.[START_REF] Ripard | Tribological characterization of greased drive-shaft : Evaluation of constant velocity joint durability[END_REF]

Figure 8 .Figure 9 .

 89 Figure 8. Grease film thickness at 80°C

Figure 10 .

 10 Figure 10. Film thickness for test number 5.

Figure 11 .

 11 Figure 11. G3 friction coefficient for test numbers 5 and 6.

  .

Figure 13 .

 13 Figure 13. Friction coefficient for test number 7.

Figure 12 .

 12 Figure 12. Wear G3 on HFRR.

Figure 14 .

 14 Figure 14. Film thickness for test number 7.

Figure 15 .

 15 Figure 15. HFRR test stopped with G3 when jump occurs.

  & 20 respectively.

Figure 17 .

 17 Figure 17. Ball scar.

Figure 16 .

 16 Figure 16. Disc scar.

Figure 18 .

 18 Figure 18. G1 results for conditions 8 & 9.

Figure 20 .

 20 Figure 20.Figure 20. G4 results for conditions 8 & 9.

Figure 20 .

 20 Figure 20.Figure 20. G4 results for conditions 8 & 9.

Figure 19 .

 19 Figure 19. G3 results for conditions 8 & 9.

Figure 21 .

 21 Figure 21. Contact replenishment with G4 at 40°C.

Figure 22 .

 22 Figure 22. Contact starvation with G3 at 40°C.

Table 1 .

 1 composition of greases tested.Four greases were used. Their compositions vary (table1).

		G1	G2	G3	G4
	Oil	A (Mineral + PAO)	A (Mineral + PAO)	B (Synthetic ester +	B (Synthetic ester +
				PAO)	PAO)
	Thickener	Lithium / Calcium	Lithium / Calcium	Polyurea	Polyurea
		Complex	Complex
	Additives	MODTC 0,28	MODTC 1	MODTC 1	MODTC 1
			MODTP 0,4	MODTP 0,4	MODTP 0,4
	Base oil dynamic viscosity	58.2	58.2	64.2	64.2
	μ 0 (cP) 40°C
	Base oil dynamic viscosity	11.1	11.1	13.5	13.5
	μ 0 (cP) 80°C S_Oil (%) IP 121 12	3.6	3.6	0.6	6
	NLGI number 22	2	2	1-2	1-2

Table 2 .

 2 initial HFRR conditions.

	Test number	Duration [h]	Temperature [°C]
	134 0
	21 8	4 0
	338 0
	41 8	8 0

Table 3 .

 3 New HFRR conditions to discriminate grease starvation.

	Test	f HFRR	d HFRR	P h-HFRR	Duration	Temperature
	number	[Hz]	[mm]	[GPa]	[h]	[°C]
	5	100 2	1.32	3	40
	68 0	2	1 .32	3	40
	7	100 2	1.32	3	80

Figure 2. HFRR test example with condition 1.

Table 4 .

 4 SD ratio variation as function of SD G1-3 .

		Test conditions number
		1&2	3&4	5	6	7
	G1	5.2	1	69.9	55.9	13.3
	G2	5.2	1	69.9	55.9	13.3
	G3	5.8	1.2	77.1	61.7	16.2
	G4	5.8	1.2	77.1	61.7	16.2

Table 5 .

 5 EHD test conditions.

	Test	S EHD	SRR EHD	P h-EHD	Duration	Temperature
	number	[mm/s]	[%]	[GPa]	[s]	[°C]
	8	400	10	0.38	2100	40
	9	400	10	0.38	2100	80

Acknowledgements

The authors would like to thank the ANRT (CIFRE N°2017/ 0094), STELLANTIS, the LaMCoS, in particular S. Wegeler, the Universidade do Porto (INEGI and FEUP) and the LTDS in particular C. Minfray and F. Dassenoy for their help and support.

Funding

The author(s) received no financial support for the research, authorship and/or publication of this article.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.


Guilhem Mollon
email: guilhem.mollon@insa-lyon.fr

The soft discrete element method

Keywords: Granular media, Discrete modelling, Soft grains, SDEM

In order to accelerate simulations of assemblies of highly deformable grains, a novel numerical approach, called Soft Discrete Element Method (SDEM) is presented. It consists in extending the classical DEM by introducing the deformability of the grains in a simplified way. Simple kinematics are postulated in order to represent the ovalisation of the grains and their local deformations around interparticle contacts. Adequate equations of motion are derived, and a contact model is proposed in the framework of elastic frictionless grains. Comparisons are made with existing analytical and numerical solutions, and the ability of this method to simulate the compressibility and some micromechanical aspects of this class of materials is validated. Upon implementation in optimized DEM codes, the SDEM is expected to allow large scale simulations of samples of deformable particles.

Introduction

The Discrete Element Method (DEM, [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]) has now become ubiquitous in the field of granular science, because it offers a large number of advantages when compared to a continuous modelling of granular samples: it is conceptually simple, versatile, and particularly informative about local quantities that may be out of reach for experimentalists. Modern refinements of this technique (e.g. use of complex/realistic grain shapes [START_REF] Mollon | 3D generation of realistic granular samples based on random fields theory and Fourier shape descriptors[END_REF][START_REF] Mollon | Discrete modelling of rock avalanches: sensitivity to block and slope geometries[END_REF][START_REF] Mollon | Can friction replace roughness in the numerical simulation of granular materials?[END_REF], enrichment of contact laws [START_REF] Scholtès | On the capillary stress tensor in wet granular materials[END_REF][START_REF] Luding | Cohesive, frictional powders: contact modes for tension[END_REF], coupling with different physics [START_REF] Tran | Numerical modelling of backward front propagation in piping erosion by DEM-LBM coupling[END_REF][START_REF] Renouf | Coupling electrical and mechanical effects in discrete element simulations[END_REF], etc.) allow it to deal with a broad range of scientific and technological applications. A major assumption of DEM, however, is the perfect rigidity of the grains (only slightly degraded by the small interpenetrations allowed at the contacts). This assumption is sensible in most applications, but becomes very questionable when the level of stress applied to the sample is of the same order of magnitude as the stiffness of the material composing the grains. If this material is ductile enough (i.e. if grain breakage is disregarded), such a situation should lead to large deformations of the grains (Fig. 1), which is prevented in the strict framework of DEM.

Several numerical approaches have been proposed to deal with highly deformable grains [START_REF] Gethin | A discrete deformable element approach for the compaction of powder systems[END_REF][START_REF] Güner | Numerical modeling of cold powder compaction using multi particle and continuum media approaches[END_REF][START_REF] Nguyen | Compaction of granular materials composed of deformable particles[END_REF][START_REF] Boromand | Jamming of deformable polygons[END_REF], and all of them rely on some sort of discretization of the grains. This is notably the case of the Multibody Meshfree Approach [START_REF] Mollon | A multibody meshfree strategy for the simulation of highly deformable granular materials[END_REF][START_REF] Mollon | A unified numerical framework for rigid and compliant granular materials[END_REF], which was implemented in the open source code MELODY [START_REF] Mollon | D[END_REF]. In this framework, a large number of field nodes are positioned in each grain and on its contour. Each field node carries two degrees of freedom in displacement, and the displacement field is interpolated between the nodes using meshfree shape functions. This choice was made to improve the robustness of the code when grains are submitted to very large deformations. Grains contours are represented by a piecewise-linear frontier, and contacts are treated by a robust two-pass node-to-segment algorithm. This technique was successfully applied in granular physics [START_REF] Mollon | Mixture of hard and soft grains: micromechanical behavior at large strains[END_REF], tribology [START_REF] Mollon | Solid flow regimes within dry sliding contacts[END_REF][START_REF] Zhang | Significance of third body rheology in friction at a dry sliding interface observed by a multibody meshfree model: influence of cohesion between particles[END_REF], geomechanics [START_REF] Mollon | Can friction replace roughness in the numerical simulation of granular materials?[END_REF], and geophysics [START_REF] Mollon | Simulating melting in 2D seismic fault gouge[END_REF].

A major limitation of this tool, however, is its computational cost. The large number of degrees of freedom and the small time steps typically limit this approach to samples of a few thousands grains at the most, while classical DEM can nowadays deal with several millions of them. To go beyond this limitation, this paper presents a novel approach for the simulation of large samples of soft grains. This approach, called the Soft Discrete Element Method (SDEM), is based on a different philosophy inspired both by the field of Model Order Reduction [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF][START_REF] Amsallem | Nonlinear model order reduction based on local reduced-order bases[END_REF][START_REF] Corigliano | Model Order Reduction and domain decomposition strategies for the solution of the dynamic elastic-plastic structural problem[END_REF] and by Moveable Cellular Automatons [START_REF] Popov | Theoretical principles ofmodeling elastoplastic media by movable cellular automata method. I. Homogeneous media[END_REF][START_REF] Salman | Implementation of MCA in the framework of LIGGGHTS[END_REF]. It consists in simplifying considerably the kinematics of the soft grains by only considering a very limited number of deformation modes. It leads to a dramatic reduction in the computational cost, but keeps the main physics at work in samples of highly deformable grains. Section 2 of this paper presents these simplified kinematics, and Section 3 develops the associated equations of motions. In Section 4, the stress tensor and the contact forces are derived, as well as the techniques associated with the numerical solver. Section 5 proposes a validation of the method by comparing it with existing analytical and numerical solutions.

Kinematics

Degrees of freedom

We first define the main frame O, ��� ⃗ e X , ��� ⃗ e Y of the problem, which remains fixed. For a given grain, an attached rigidbody motion frame is defined as C, � � ⃗ e x , � � ⃗ e y , with C being its center of mass. In the proposed numerical framework a 2D deformable grain is assumed to have only six degrees of freedom. Three of them are related to rigid-body motions: two translations of the grain center of mass x c (t) and y c (t) (which define the position of C with respect to O ) and one rotation (t) (which is the angle of the direction � � ⃗ e x with respect to ��� ⃗ e X ). The remaining three degrees of freedom are related to the deformations of the grain: they are the two principal strains 1 (t) and 2 (t) (with no specific ordering) and the angle (t) of orientation of the principal deforma- tion frame C, � � ⃗ e 1 , � � ⃗ e 2 with respect to the rigid-body frame C, � � ⃗ e x , � � ⃗ e y attached to the grain. Namely, (t) is the angle between the direction � � ⃗ e x and the direction � � ⃗ e 1 of the main stretching 1 (t) . This is summarized in Fig. 2, for a grain with an initial radius R = 0.5 . 1 (t) , 2 (t) and (t) actually define a strain tensor, which is assumed to be constant over the whole domain covered by the grain:

When submitted to such a homogeneous strain field, a circular grain will take the shape of an ellipse with two orthogonal half-axes of respective lengths R ⋅ 1 + 1 (t) and R ⋅ 1 + 2 (t) , oriented along the frame C, � � ⃗ e 1 , � � ⃗ e 2 . These kinematics are thus more complex that those commonly encountered in 2D DEM, as they offer the possibility for each grain to gain a certain amount of ellipticity as a response to the loads it is submitted to. This is, obviously, a strong simplification of the real kinematics of deformable grains which, to be represented in detail, would require a much larger number of degrees of freedom (using for example shape functions attached to field nodes, like in the multibody meshfree approach, [START_REF] Mollon | A multibody meshfree strategy for the simulation of highly deformable granular materials[END_REF][START_REF] Mollon | A unified numerical framework for rigid and compliant granular materials[END_REF]). As we demonstrate in the remainder of this paper, it provides nevertheless a satisfactory description of the main physics of assemblies of deformable grains.

The main question that arises is: how to determine the evolution in time of these six degrees of freedom? Newtonian dynamics can be readily applied to the three components of the rigid body motion (accounting for the possible changes in rotational inertia brought by the deformations), but the case of the three other degrees of freedom is less straightforward. The usual way to predict the deformation of a solid is to employ a constitutive model, i.e. a relation (usually based on thermodynamics and phenomenological observations) between the stress tensor and the strain tensor. Since the strain is assumed to be homogeneous in the grain, so must be the stress. We thus define a homogeneous field of stress tensor within a given grain as follows:

(1) The principal directions of this stress tensor are aligned with those of the strain tensor if the chosen constitutive model is isotropic, which is postulated in the remainder of this paper for the sake of simplicity.

𝜀 = 𝜀 1 0 0 𝜀 2 {� ⃗ e 1 , � ⃗ e 2 } = 𝜀 xx 𝜀 xy 𝜀 xy 𝜀 yy {� ⃗ e x , � ⃗ e y } = 𝜀 XX 𝜀 XY 𝜀 XY 𝜀 YY {� � ⃗ e X , � � ⃗ e Y } (2)
𝜎 = 𝜎 1 0 0 𝜎 2 {� ⃗ e 1 , � ⃗ e 2 } = 𝜎 xx 𝜎 xy 𝜎 xy 𝜎 yy {� ⃗ e x , � ⃗ e y } = 𝜎 XX 𝜎 XY 𝜎 XY 𝜎 YY {� � ⃗ e X , � � ⃗ e Y }
A given grain is submitted to a number of external forces applied on its external contour. These forces can be summarized, as is often done in granular science, by an equivalent stress field called ext . If static equilibrium is reached, the resulting force and torque applied to the grain are null and the stress within the grain is equal to that induced by the external forces ext . However, DEM simulations are dynamic by nature, and so is the proposed method. Hence, static equilibrium is never perfectly enforced. There is always some amount of unbalanced force and, in the present case, of unbalanced stress. Such a system hence relies on equations of motion, which have to be written for each degree of freedom. Any such equation must establish a relation between a generalized force and a generalized mass (both associated to a given degree of freedom) on one hand, and the second time-derivative of the concerned degree of freedom, on the other hand. The six equations of motion are derived in Section 3 for the postulated kinematics.

Deformations around contacts

The general ovalisation of the grain shape captures well some features of the grains deformation shown in Fig. 1, but a major element is still missing. It is, indeed, impossible with such kinematics to represent properly the closure of the intergranular space under a sufficiently large confining stress. The missing ingredient is the local deformation of the grains in the direct neighbourhood of their contacts. We clearly observe in Fig. 1 that, if the grains are squeezed enough, they do not exhibit any more a point contact like in usual rigid granular materials. Instead of this, the matter around the contact deforms to accommodate the load, and the contact point becomes a contact area (or contact line in 2D). As the load increases, these contact lines occupy a larger proportion of the grains contours, until the whole contour of each grain is in contact with other grains and the porosity is closed. We also observe in Fig. 1 that, at least Fig. 2 Illustration of the six degrees of freedom for an initially circular grain. The path to the current configuration can be conceptualized either by a deformation (a) followed by a rigid body motion (b), or by a rigid body motion (c) followed by a deformation (d) in the case where the constitutive model is the same for all grains, the contact line between any two grains is never far from a linear segment.

In classical DEM, a small interpenetration is allowed between contacting grains, in order to compute a repulsive force opposing to contact. In most cases, the normal stiffness used to compute this force is seen as a numerical parameter of regularization rather than as a physical quantity. In some contexts (e.g. small strain behaviours of sands), this interpenetration is seen as a proxy for the local contact-induced deformation. The classical Hertz theory ( [START_REF] Hertz | Uber die Beruhrung fester elastischer Korper[END_REF]) can for example be applied to compute the repulsive normal force between two spheres on physical grounds, and the Hertz-Mindlin model ( [START_REF] Mindlin | Compliance of elastic bodies in contact[END_REF]) can be applied the same way to derive both normal and tangential components of this force. These approaches, however, rely on rather restrictive assumptions, including the necessity for the deformation (and thus for the interpenetration) to remain small with respect to the grain size. For this reason, they cannot be directly applied in the present context.

Precisely describing the large deformation of the matter composing the grains in the neighbourhood of a heavily loaded contact would require a much larger number of degrees of freedom than the six that were defined in the previous section, and would contradict the general spirit of the method. Alternatively, we propose here to account for this phenomenon in a simplified way. We accept that very large interpenetrations take place between the ellipses representing the contacting grains, but we then assume that the real contours of the grains in that case are not completely represented by the ellipses. In the contact area, the contour is replaced by a straight segment linking the two intersection points between the ellipses, as represented in Fig. 3. Hence, in the case of multiple contacts, a given grain contour is a succession of straight segments and of portions of the main ellipse defined by the kinematics of the previous section.

Obviously, this assumption means that a large part of the grains surface areas (i.e. all the area of each ellipse located beyond the contact segment) is ignored, which could be seen as a problem in terms of conservation laws. A simple solution exists to solve this issue: since this loss of surface can be computed, it can be added to the strain tensor as a spherical (i.e. purely volumetric) term, and thus be accounted for when applying the constitutive model to compute the stress tensor. If we assume, for example, that the material composing the grains is incompressible (i.e. any loss of surface is prohibited), then the surface area lost beyond the contact segment will be regained by an increase of 1 and 2 in accordance with the chosen constitutive model. This approach is detailed in Section 4, as well as the calibration of the contact forces associated with such a contact segment.

Equations of motion

Degrees of freedom of the rigid-body part of the motion:

In this paragraph, we consider that a grain with an initial radius R and a unit mass has acquired a certain deformation (defined by the degrees of freedom 1 , 2 , and ), and that it is subjected to a rigid-body motion while keeping a constant shape. The case of the translational degrees of freedom x c and y c is straightforward, as they obey to the classical New- tonian laws of dynamics, written as:

where F x (t) and F y (t) are the resulting forces applied on the grain and:

The case of the rotational rigid motion is just less straightforward, since it requires to account for the current shape of the grain. It writes:

where C (t) is the resulting torque applied on the grain, and M (t) is the current rotational inertia of the deformed elliptic grain (Fig. 4a):

Degrees of freedom " 1 and " 2

In this paragraph and the next one, we consider that the degrees of freedom related to the rigid motion of the grains (i.e. x c , y c , and ) are kept constant and equal to zero. We thus only consider a deformation of the grain defined by two principal strains 1 and 2 , oriented along an angle . Under these assumptions, a point of position (X, Y) in the reference configuration is located at a certain time t at the location (x(t), y(t)) given by:

(3) We first consider that 2 and are constant in time, and we compute the velocity of this point as a function of 1 (t):

⎧ ⎪ ⎨ ⎪ ⎩ ẍc (t) = F x (t) M x ÿc (t) = F y (t) M y (4) M x = M y = R 2 (5) θ(t) = C 𝜃 (t) M 𝜃 (t) (6) M (t) = 1 + 1 (t) 2 ⋅ 1 + 2 (t) 2 ⋅ R 4 2
If we assign a mass dm to every such material point, the total kinetic energy corresponding at a certain time t to the variation of 1 (t) is thus given by: where Ω is, for example, the geometric domain occupied by the body in its reference configuration. Thus, in polar coordinates (r, ) , we have X = r ⋅ cos and Y = r ⋅ sin , giving: This energy can be reorganized under the form:

where M 1 is the generalized mass associated with the degree of freedom 1 , and given by: Which yields: And finally: [START_REF] Tran | Numerical modelling of backward front propagation in piping erosion by DEM-LBM coupling[END_REF] x

(t)= 1 + 1 (t) ⋅ X ⋅ cos 2 (t) + Y ⋅ cos (t) ⋅ sin (t) + 1 + 2 (t) ⋅ X ⋅ sin 2 (t) -Y ⋅ cos (t) ⋅ sin (t) y(t)= 1 + 1 (t) ⋅ X ⋅ cos (t) ⋅ sin (t) + Y ⋅ sin 2 (t) ⋅ sin (t) + 1 + 2 (t) ⋅ -X ⋅ cos (t) ⋅ sin (t) + Y ⋅ cos 2 (t) (8) ̇x(t) = ̇𝜀 1 (t) ⋅ X ⋅ cos 2 𝛼 + Y ⋅ cos 𝛼 ⋅ sin 𝛼 ̇y(t) = ̇𝜀 1 (t) ⋅ X ⋅ cos 𝛼 ⋅ sin 𝛼 + Y ⋅ sin 2 𝛼 (9) E c1 (t) = ∫ Ω ̇x2 (t) + ̇y2 (t)
2 dm (10)

E c1 (t) = 2𝜋 ∫ 0 R ∫ 0 ̇𝜀 2 1 (t) ⋅ r ⋅ cos 𝜃 ⋅ cos 2 𝛼 + r ⋅ sin 𝜃 ⋅ cos 𝛼 ⋅ sin 𝛼 2 2 + ̇𝜀 2 1 (t) ⋅ r ⋅ cos 𝜃 ⋅ cos 𝛼 ⋅ sin 𝛼 + r ⋅ sin 𝜃 ⋅ sin 2 𝛼 2 2 𝜌 ⋅ r ⋅ drd𝜃 (11) E c1 (t) = 1 2 M 1 ⋅ ̇𝜀 2 1 (t) (12) M 1 = ⋅ 2 ∫ 0 R ∫ 0 r ⋅ cos ⋅ cos 2 + r ⋅ sin ⋅ cos ⋅ sin 2 + r ⋅ cos ⋅ cos ⋅ sin + r ⋅ sin ⋅ sin 2 2 ⋅ r ⋅ drd (13) M 1 = ⋅ 2 ∫ 0 R ∫ 0 r 3 ⋅ cos 2 ( -) ⋅ drd
In order to compute the generalized force associated with the same degree of freedom, we postulate that there exists an unbalanced homogeneous stress field s(t) in the body at a certain time t:

This unbalanced stress corresponds to the stress ext (t) induced by external forces applied on the grain, corrected by the stress (t) related to the current strain tensor of the body by the means of a given constitutive model (see section 4). Since this field is homogeneous in the grain, there is an associated resulting force field in each material point of initial coordinates (X, Y):

Hence, the generalized force associated with the degree of freedom 1 is given by the general expression:

where (x(t), y(t)) are the current coordinates of the point of reference coordinates (X, Y) . Still considering 2 and as constant in time, and using Eq. ( 7), we have: These expressions are then injected in F 1 (t) from Eq. ( 17):

(14) M 1 = R 4 ∕4 (15) s(t) = s XX (t) s XY (t) s XY (t) s YY (t) {� � ⃗ e X , � � ⃗ e Y } = s 11 (t) s 12 (t) s 12 (t) s 22 (t) {� ⃗ e 1 , � ⃗ e 2 } (16) F x (X, Y, t) =- s XX (t) X - s XY (t) Y F y (X, Y, t) =- s YY (t) Y - s XY (t) X ( 17
)
F 1 (t) = ∫ Ω x 1 (X, Y, t) ⋅ F x (X, Y, t) + y 1 (X, Y, t) ⋅ F y (X, Y, t) dΩ (18) ⎧ ⎪ ⎨ ⎪ ⎩ x 1 (X, Y, t) = X ⋅ cos 2 + Y ⋅ cos ⋅ sin y 1 (X, Y, t) = X ⋅ cos ⋅ sin + Y ⋅ sin 2 (19)
F 1 (t) =-∫ Ω X ⋅ cos 2 + Y ⋅ cos ⋅ sin ⋅ s XX (t) X + s XY (t) Y + X ⋅ cos ⋅ sin + Y ⋅ sin 2 ⋅ s YY (t) Y + s XY (t) X dΩ
If the grain is considered as isotropic, the generalized force of the degree of freedom 1 (related to the extension/ compression in the direction � � ⃗ e 1 ) should be independent from the value of the angle , and we can thus choose = 0. In that case, the frames O, � � ⃗ e x , � � ⃗ e y and O, � � ⃗ e 1 , � � ⃗ e 2 are identical, and we can rewrite the generalized force:

Integrating this simple expression on the reference domain finally leads to:

This result provides the generalized expression for the dynamics of the degree of freedom 1 : Similar expressions are derived for 2 . We have: with M 2 = M 1 and:

Degree of freedom

We now consider that 1 and 2 are constant in time (Fig. 4b), and we derive the new expressions for the current velocity of a given point of initial coordinates (X, Y) and of current coordinates (x(t), y(t)) . For that purpose, we differentiate in time the Eq. ( 7):

This expression simplifies to:

(20) F 1 (t) =-∫ Ω X ⋅ s 11 (t) X + s 12 (t) Y dΩ (21) F 1 (t) =-R 2 ⋅ s 11 (t) (22) ̈𝜀 1 (t) = F 1 (t) M 1 (23) ̈𝜀 2 (t) = F 2 (t) M 2 (24) F 2 (t) =-R 2 ⋅ s 22 (t) (25
)

̇x(t) =-2 1 + 𝜀 1 X ⋅ sin 𝛼(t) ⋅ cos 𝛼(t) ⋅ ̇𝛼 (t) + 2 1 + 𝜀 2 X ⋅ sin 𝛼(t) ⋅ cos 𝛼(t) ⋅ ̇𝛼 (t) + Y 𝜀 1 -𝜀 2 ⋅ cos 2 𝛼(t) -sin 2 𝛼(t) ⋅ ̇𝛼 (t) ̇y(t) =-2 1 + 𝜀 2 Y ⋅ sin 𝛼(t) ⋅ cos 𝛼(t) ⋅ ̇𝛼 (t) + 2 1 + 𝜀 1 Y ⋅ sin 𝛼(t) ⋅ cos 𝛼(t) ⋅ ̇𝛼 (t) + X 𝜀 1 -𝜀 2 ⋅ cos 2 𝛼(t) -sin 2 𝛼(t) ⋅ ̇𝛼 (t) (26
)

̇x(t) = ̇𝛼 (t) ⋅ 𝜀 1 -𝜀 2 ⋅ -X ⋅ sin 2 𝛼(t) + Y ⋅ cos 2 𝛼(t) ̇y(t) = ̇𝛼 (t) ⋅ 𝜀 1 -𝜀 2 ⋅ Y ⋅ sin 2 𝛼(t) + X ⋅ cos 2 𝛼(t)
The kinetic energy associated with this degree of freedom is thus equal to:

As in the previous paragraph, this energy can be reorganized under the form:

where M (t) is the generalized mass associated with this degree of freedom, and given by: This expression simplifies to: After a simple integration on the reference domain, the generalized mass is:

It is interesting to observe that this mass is not constant in time, since it depends on the values of 1 and 2 (it vanishes if the grain is circular, i.e. if 1 = 2 , since the grain can- not deform along the degree of freedom in that particular case). When implemented numerically, M will thus have to be updated at each time step.

To compute the corresponding generalized force, we use a formula analogous to Eq. ( 17):

where the force field is given in Eq. ( 16). The partial derivatives are computed from Eq. ( 7) and give, at the end of the calculation:

(27) E c𝛼 (t) = ∫ Ω ̇x2 (t) + ̇y2 (t) 2 dm (28) E c𝛼 (t) = 1 2 M 𝛼 (t) ⋅ ̇𝛼 2 (t) (29) M (t) = ∫ Ω 1 -2 2 ⋅ -X ⋅ sin 2 (t) + Y ⋅ cos 2 (t) 2 + 1 -2 2 ⋅ Y ⋅ sin 2 (t) + X ⋅ cos 2 (t) 2 dm (30) M (t) = 1 -2 2 ∫ Ω X 2 + Y 2 dm (31) M (t) = 1 (t) -2 (t) 2 R 4 2 (32) F (t) = ∫ Ω x (X, Y, t) ⋅ F x (X, Y, t) + y (X, Y, t) ⋅ F y (X, Y, t) dΩ
Injecting them and the force field in Eq. ( 32) yields:

For reasons similar to those used in the previous paragraph, this generalized force should be independent on the value of the angle (i.e. it should be frame-independent to remain objective), and we can use an arbitrary value for this parameter, such as = 0 ∶ After integration, we finally obtain: Just like the generalized mass, this force vanishes for a circular shape (i.e. 1 = 2 ) . The general expression of the time-dynamics of the degree of freedom is thus:

Summary

Based on the calculations derived in the previous paragraph, we can summarize the equations of motion of a given grain of initial radius R and of initial unit mass : Using these equations of motion requires to compute both the external forces applied on the grain and its internal stress field. [START_REF] Cantor | Compaction model for highly deformable particle assemblies[END_REF] x (X,

Y, t) = 1 -2 ⋅ -X ⋅ sin 2 (t) + Y cos 2 (t) y (X, Y, t) = 1 -2 ⋅ Y ⋅ sin 2 (t) + X cos 2 (t) (34)
F (t) =-∫ Ω 1 -2 ⋅ -X ⋅ sin 2 (t) + Y cos 2 (t) ⋅ s XX (t) X + s XY (t) Y + 1 -2 ⋅ Y ⋅ sin 2 (t) + X cos 2 (t) ⋅ s YY (t) Y + s XY (t) X dΩ (35) F (t) =-1 -2 ∫ Ω Y ⋅ s 11 (t) X + s 12 (t) Y + X ⋅ s 22 (t) Y + s 12 (t) X dΩ (36) F (t) =-2 R 2 ⋅ 1 (t) -2 (t) ⋅ s 12 (t) (37) ̈𝛼 (t) = F 𝛼 (t) M 𝛼 (t) (38)
⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ẍc (t) = F x (t) 𝜌𝜋R 2 ÿc (t) = F y (t) 𝜌𝜋R 2 θ(t) = 2⋅C 𝜃 (t) (1+𝜀1(t)) 2 ⋅(1+𝜀 2 (t)) 2 ⋅𝜌𝜋R 4
and

⎧ ⎪ ⎨ ⎪ ⎩ ̈𝜀 1 (t) =- 4⋅s 11 (t) 𝜌R 2 ̈𝜀 2 (t) =- 4⋅s 22 (t) 𝜌R 2 ̈𝛼 (t) =- 4⋅s 12 (t) (𝜀1(t)-𝜀2(t))𝜌R 2
4 Forces and stresses

Stress tensors

As exposed in the previous sections, two stress tensors need to be defined. The first one, noted , is related in a univocal way to the strain tensor by the means of the constitutive model, and thus represents the actual stress in the grain.

The second one, noted ext , serves as a compact reduction of the internal efforts induced in the grain by all the contact forces applied on its contour. As a first step into this novel technique, very simple assumptions are to be made. We assume in this paper that a linear relationship exists between the stress and strain tensors. Since we wish to account for the surface loss at the contact segments between grains, we need to distinguish the volumetric and deviatoric parts of the strain tensor. The deviatoric part is:

The volumetric part accounts both for the surface change related to the main deformations 1 and 2 , and for the surface losses S j related to the intersection of the corresponding grain with any other contacting grain j (Fig. 3c):

It should be noted that this evaluation of the volumetric strain does not rely on the trace of , since this approach would only be licit in the limit of small strains. In the more general case of finite strains, which is more likely to be encountered in the present situation, it is necessary to compute directly the area of the deformed grain based on the lengths of its main axes. With this decomposition, we can thus deduce the following stress field:

(39) 1dev (t) = 1 (t) -1 (t) + 2 (t) ∕2 (40) 2dev (t) = 2 (t) -1 (t) + 2 (t) ∕2 (41) V (t) = R 2 ⋅ � 1 + 1 (t) �� 1 + 2 (t) � - ∑ j S j R 2 -1 (42) 1 (t) = (2 + ) ⋅ 1dev (t) + ⋅ 2dev (t) + ⋅ V (t) (43) 2 (t) = ⋅ 1dev (t) + (2 + ) ⋅ 2dev (t) + ⋅ V (t)
where and are the Lamé coefficients and is the bulk modulus, related to Young's modulus and Poisson's coefficient by: Hence, equations of linearized elasticity are applied even though finite strains are expected in the system. Future extensions of the method might include more complex constitutive relations, such as hyper-elasticity or elasto-plasticity.

The external stress tensor associated to the contact forces applied on a given grain is given by the classical Love-Weber formula ( [START_REF] Nicot | On the definition of the stress tensor in granular media[END_REF]):

(44) = E (1 + )(1 -2 ) (45) = E 2(1 + ) (46) = E (3 -6 ) (47) ext ij = 1 S ∑ contacts F i ⋅ r j
where S is the particle surface area, F i are the components of the force ��� ⃗ F c associated to a given contact, and r j are the com- ponents of the corresponding branch vector (i.e. the vector linking the grain center of mass and a representative contact point). When both and ext are determined, the unbalanced stress to be used in the equations of motion of the previous section is simply given by:

Contact algorithm and contact law

Applying Eq. ( 47) requires the computation of contact forces between two grains (and possibly between a single grain and a boundary wall). This force should depend on the local contact conditions (based on quantities such as interpenetration distance, etc.) and on the constitutive behavior of the contacting grains. A contact algorithm is therefore necessary. Several algorithms were proposed in the literature to characterize the intersection between arbitrary ellipses [START_REF] Rothenburg | Numerical simulation of idealized granular assemblies with plane elliptical particles[END_REF][START_REF] Ting | A robust algorithm for ellipse-based discrete element modelling of granular materials[END_REF][START_REF] Ting | An ellipsebased discrete element model for granular materials[END_REF].

they were proposed in a classical DEM framework, they were optimized for the case of small interpenetrations, (48) and relied on the fact that the shape and size of the ellipses were considered as constant.

s(t) = ext (t) -(t)
In order to bring more generality to the method and to its future extensions, the first implementation of SDEM that is presented here is based on a discretization of the contour of each ellipse, with an angular parameterization of the reference (circular and un-rotated) shape of each grain (Fig. 5). For a given grain i of initial radius R i , a number of nodes N n is positioned on the contour of reference, parameterized by the angle k = [0, ,2 , … ,2 -] with an angular step = 2 ∕N n (a value of = 1 • is used in the present work). At a given time in a simulation, if the degrees of freedom of the grains i take the values x ci , y ci , i , 1i , 2i , i , the coor- dinates of these contour nodes become:

With the reference coordinates:

And with:

The transformation defined by Eq. ( 49) allows for a rapid evaluation of the current coordinates of each node on the (49)

x ki y ki = A ⋅ X ki Y ki + B (50) X ki = R i ⋅ cos k Y ki = R i ⋅ sin k (51) B = x ci y ci (52) A = a 11 a 12 a 21 a 22 (53) ⎧ ⎪ ⎨ ⎪ ⎩ a 11 = � 1 cos 2 a + 2 sin 2 a + 1 � cos - � 1 -2 � cos sin sin a 12 = � 1 -2 � cos sin cos - � 1 sin 2 + 2 cos 2 + 1 � sin a 21 = � 1 cos 2 a + 2 sin 2 a + 1 � sin + � 1 -2 � cos sin cos a 22 = � 1 -2 � cos sin sin + � 1 sin 2 + 2 cos 2 + 1 � cos
contour of the grains. A simple intersection algorithm based on a piecewise linear approximation of the current contour of each grain is then applied. For any pair of contacting grains i and j , this allows for a rapid computation of key quanti- ties such as the normal gaps ni and nj (taken as negative is there is an interpenetration of the grains), the interpenetration areas S ij and S ji , and the two intersection points P and P ′ . The contact length L c is then taken as the distance PP ′ , defining a tangential vector ⃗ t = ����� ⃗ PP � ∕L c and its associated normal vector (Fig. 5). A representative point of contact P ij is also defined as the mid-point between P and P ′ , and used as the point of application of the contact force (in particular for the computation of C and ext ). Similar calculations are done for contacts between a grain and an infinite boundary wall.

At that point, it is necessary to establish a relation between the local contact parameters and the contact force. In classical DEM, the standard contact model relies on the concept of contact stiffness, which is usually seen as a numerical regularization of the non-penetration condition. Since it does not bear any physical meaning (except in a limited number of situations, e.g. [START_REF] Gu | DEM simulations of the small strain stiffness of granular soils: effect of stress ratio[END_REF][START_REF] Richefeu | Dissipative contacts and realistic block shapes for modeling rock avalanches[END_REF]), its value is determined in order to ensure limited interpenetrations under the expected average normal load on each contact. In SDEM, however, the non-penetration condition is not enforced, since large overlaps between the deformed grains are expected. Besides, the contact force is seen in this framework as the result of local deformations of the grains, and should therefore be related in some way to their constitutive models. While the Hertzian model is licit at small interpenetrations, its general assumptions are too restrictive to be applied to the present case, especially because of the large deformations of the grains around the contacts. In the general case, it is expected that no closed-form solution can be found for this force, but a simple heuristic model can be proposed under the present assumptions (elastic grains, no friction, no cohesion). We postulate that the contact-related deformation of the grain remains local, and is limited to an area defined by (i) the contact segment and (ii) a curve which is defined as Fig. 6 Sketch of the proposed contact model the symmetrical curve to the initial contour of the grain (i.e. as it would be in the absence of contact) but with a certain scaling factor (Fig. 6). The exact expression of these curves (which are portions of ellipses) does not need to be known. In this area, the deformation is postulated to be constant, and to occur only in the direction normal to the contact. It is thus called n . We also define the normal gap n as the average between the gaps computed for the two grains concerned with the contact: We then assume that the maximum dimension of the deformed area in the normal direction is a linear composition of the two characteristic dimensions of the contact, namely its length L c and its normal gap n . With these assumptions, the normal strain n in the deformed area is related to the ratio between the initial dimensions of the deformed area (in the direction normal to the contact) to its current dimensions.

Considering the point of maximum interpenetration (Fig. 6), we can write:

Considering linear elasticity and integrating the normal stress on the whole contact length, we finally get a repulsive normal force equal to: With a secant normal stiffness k n given by: where a and b are two constants. It is interesting to notice that this expression is independent on the grain radius. To determine a and b , numerical simulations of the isotropic compression of a deformable grain between four walls are performed with the code MELODY [START_REF] Mollon | D[END_REF]. The grain has a diameter equal to 1 and a unit Young modulus, and several values of the Poisson Coefficient (ranging between 0 and 0.5) are tested. This process is summarized in Fig. 7a. As observed in Fig. 7b, the grain shape flattens at the contact, but the stress level remains rather constant on the contact segment, as postulated in the proposed contact model. At each time of each simulation, it is thus possible to measure the value of L c and to compute the value of n (by consid- ering the virtual circle that crosses the square box at the extremities of the contact segments). For given values of a and b , it is thus possible to apply the Eqs. (56-57) and to compare the obtained force with that provided by the simulation. This is done in Fig. 7c with a = 0.58 and b = 2 (i.e. with the same parameters for all values of the Poisson coefficient). The agreement is fairly good, except for low Poisson Coefficients (i.e. 𝜈<0.2 at very large compressive strains, where the model underestimates the normal contact force by up to 15% (Fig. 7d). Additional tests performed with different values of the radius and of the Young modulus (not reported here) confirmed this agreement. This is quite remarkable because the proposed contact model is very

(54) 𝛿 n = 𝛿 ni + 𝛿 nj 2 < 0 (55) n = n a ⋅ L c + b ⋅ n -n (56) F n = E ⋅ n ⋅ L c = k n ⋅ n (57) k n = E ⋅ L c a ⋅ L c + (b -1) ⋅ n
simplistic and relies on assumptions which are obviously wrong. It nevertheless offers a simple way to estimate the repulsive force between deformable grains with large interpenetration, at a very small computational cost.

Solver, stability, damping, and incompressibility

In order to integrate in time the equations of motion summarized in Eq. ( 38), a time step Δt is chosen, and an explicit solver is employed in the same manner as classical DEM. Considering for example the degree of freedom x c , we apply at each time step the following algorithm: Analogous formulations are used for the five other degrees of freedom. Since this solver is explicit, the question of its stability is essential. A critical point is the degree of freedom , because its associated generalized mass M , provided in Eq. ( 29), depends on the current state of deformation of the grain. More importantly, in the circular state of reference of each grain, this generalized mass vanishes to 0. It means that, for very small (in absolute value) levels of deformation 1 -2 , this mass is extremely small, and requires a very small time step. This is clearly a waste of computational resources, since the parameter has only a very small influence on the grain shape when the deviatoric deformation is low. A simple way to bypass this limitation is to ignore the generalized force F is the mass M is too small compared to the generalized mass of the other deformationrelated degrees of freedom. Hence, at each time step and for each grain, the following condition is checked:

where k is a user-defined constant. If this condition is not verified, the force F receives the value 0 before applying Eqs. (58-59), which prevents any instability. Typical values of k of the order of 10 -3 to provide satisfactory results. An alternative approach is to apply a constant mass scaling to M (i.e. to multiply it by a factor 100, for example), but it is only licit in quasi-static situations (i.e. when the dynamic response of this degree of freedom is unimportant), and it does not solve the singular case M = 0 (when 1 = 2 , for example in the initial circular state). In the illustrative simulations described in the next section, a combination of both techniques was implemented.

Another condition to ensure stability is to introduce a certain amount of damping to dissipate kinetic energy. In classical DEM this is introduced in the contact law, and this

(58) ̇xc ← ̇xc +Δt ⋅ F x (t) M x (59) x c ← x c +Δt ⋅ ̇xc (60) M 𝛼 > k 𝛼 ⋅ M 𝜀1
approach is reproduced in the present study. We thus add to each contact force a damping term with the expression:

where is a damping coefficient, k n is the secant contact stiffness provided in Eq. ( 57), M i and M j are the masses of the two grains into contact, and δn is the time-derivative of the normal gap (computed by numerical derivation based on the previous time step). In addition, we must apply a damping on the degrees of freedom involving a deformation of the grains (i.e. 1 , 2 , and ) in order to attenuate their free oscillations in the absence of contact. Hence, we add the following terms to the generalized forces of 1 , 2 , and respectively (Eqs. ( 21) and ( 36)): quasi-static simulations, a value = 1 proved satisfactory.

In many practical and academic situations, soft grains are considered as incompressible (or at least, their compressibility is negligible compared to their deformability). This is problematic, because incompressibility is notoriously difficult to introduce in simulated elastic media (may it be in a discrete or in a continuous framework), due to the divergence of the bulk modulus towards infinity. A common solution for this issue is to introduce a Poisson coefficient close to-but not to-0.5. Typical values could be 0.49, 0.499, etc. This is analogous to a numerical penalization of any shift from perfect incompressibility. This solution is but presents two drawbacks: incompressibility is only enforced in an way, and it strongly reduces the critical time step of explicit solvers (because the bulk modulus, although not infinite, is very high). In the SDEM framework introduced here, this drawback can be circumvented by applying a simple technique. It consists in applying a Poisson coefficient strictly equal to 0.5, but in computing the internal stress field based only on the deviatoric part of the strain tensor. Hence, Eqs. (42-43) are replaced by: (61)

F damp =-𝛾 ⋅ √ √ √ √ k n 1 M i + 1 M j ⋅ δn (62) F 1,damp =-𝛾 ⋅ M 1 ⋅ ̇𝜀 1 (63) F 2,damp =-𝛾 ⋅ M 2 ⋅ ̇𝜀 2 (64) F 𝛼,damp =-𝛾 ⋅ M 𝛼 ⋅ ̇𝛼 (65) 1 (t) = 2 ⋅ 1dev (t)
This avoids any problem related to the singularity of . To enforce incompressibility, though, it is necessary to keep the volume of the grains constant. Hence, after updating the deformation of the grain in Eqs. (58-59) based on the stress tensor of Eqs. (65-66), an additional isotropic deformation should be added to the grain. The added volumetric strain Vadd should verify (for a given grain i in contact with several grains j):

Which yields:

We then simply add Vadd ∕2 to 1dev and 2dev before starting a new time step.

Validation

Isotropic compression

In order to check that the proposed method offers a satisfactory representation of the mechanical response of collections of deformable grains, we first consider the isotropic compaction of two samples of ~ 1000 grains. The sample A follows (66)

2 (t) = 2 ⋅ 2dev (t) (67) R 2 1 + 1dev + Vadd 2 1 + 2dev + Vadd 2 - ∑ j S j = R 2 (68) Vadd =-1dev -2dev -2 + � � � � � 1dev + 2dev + 2 � 2 -4 � � 1 + 1dev �� 1 + 2dev � -1 - ∑ j S j R 2
� Fig. 8 Comparison with the analytical solution of [START_REF] Cantor | Compaction model for highly deformable particle assemblies[END_REF] for incompressible grains 13 a rather narrow Gaussian size distribution, with an average radius of 0.012, a standard deviation of 0.0024 (coefficient of variation of 0.2), and lower and upper size cut-offs at 0.008 and 0.016 respectively. The sample B follows a fractal distribution of fractal dimension 2, with lower and upper cut-offs at 0.004 and 0.038 respectively. Both samples are frictionless, with a Young's modulus equal to 1, and a unit density. They are positioned between four rigid walls forming a unit square box and submitted to a strain-driven isotropic compaction in quasi-static conditions (i.e. at a slow enough pace to avoid any inertial effect).

It was recently shown in [START_REF] Cantor | Compaction model for highly deformable particle assemblies[END_REF] that an analytical solution could be inferred for this kind of physical system in the case of incompressible grains. The mean stress in the sample was indeed reported to be well described by the following expression, E is the Young's modulus of the grains, is the current solid fraction, Z 0 and 0 are the coordination number and the solid fraction at the jamming transition, and max is the maximum solid fraction that can be attained in the (this value is reported to depend mostly on the friction coefficient between the grains, and is usually very close to 1). The pre-factor k and the exponent a are inherited from Zrelations reported in various studies [START_REF] Andreotti | Granular media, between fluid and solid[END_REF][START_REF] Vu | Numerical simulations of the compaction of assemblies of rubberlike particles: a quantitative comparison with experiments[END_REF][START_REF] Nezamabadi | Parallel implicit contact algorithm for soft particle systems[END_REF] and take the values 5.1 and 0.5 respectively. The value of b was fitted in [START_REF] Cantor | Compaction model for highly deformable particle assemblies[END_REF] based on single-grain compaction tests and is equal to 0.12.

Figure 8 shows a comparison between the n -curves obtained with SDEM and with this analytical solution. For this comparison, the values of Z 0 and 0 are first measured at the jamming point for both samples (equal to 4.20 and 0.833 respectively for sample A, 4.22 and 0.858 respectively for sample B), and the value of max is set to 1. The agreement both approaches is very satisfactory. To validate SDEM with different values of the Poisson's coefficient, however, a different approach must be used since the expression of Eq. ( 69) is only valid for = 0.5 . Simulations of isotropic compactions are thus performed on the same samples and in the same conditions with the code MELODY, using the Multibody Meshfree Approach. Qualitative results are provided in Fig. 9 for both samples, in the case = 0.3 . Four loading stages (provided in Fig. 10) are considered. Zoomed views on the same area of the sample are provided, and the grains are represented with a color scale corresponding to the local value of the mean

(69) n =-E ⋅ b 2 ⋅ Z 0 + k -0 a ⋅ ln max - max -0 stress 1 + 2 ∕2
. The agreement between both numerical approaches is acceptable, both in terms of mean stress values (although this quantity is constant in each grain for SDEM while it is allowed to vary spatially in MELODY thanks to a much finer discretization) and in grains motions and deformations. The progressive closure of the void space and the appearance of straight contact segments between grains at larger compaction is also well-rendered, even when the contacting grains have very different radii (Sample B). At loading stage d. (i.e. for a mean stress applied on the sample close to 0.3), the saturation of the sample is almost complete, and the stress level in the grains is slightly over-estimated by SDEM when compared with MELODY. This is probably related to the decrease in accuracy of the simple contact law calibrated in Fig. 7 when the interpenetration becomes too large. It may be corrected in future versions of the method with more accurate contact models. Figure 10 covers the cases = 0.3 and = 0.5 , and pro- vides a quantitative view on this comparison. As shown in the upper part of the figure, SDEM is able to provide a good prediction of the compressibility of all samples, although this compressibility is slightly under-estimated at low pressures (i.e. 𝜎 n < 0.01 ) and slightly over-estimated at large pressures (i.e. 𝜎 n > 0.1 ). The lower part of the figure con- firms that SDEM can also bring fair estimates of microstructural features of the samples, as it provides good first-order evaluations of the evolution of the coordination number Z n during compaction. The prediction of Z is especially good for = 0.3 , but a bit too large for low values of n in the case of incompressible grains. It is interesting to notice that the quality of the SDEM approximation is not degraded by the presence of a very broad size distribution of the grains (Sample B).

Oedometric compression

In order to check the validity of SDEM under a less isotropic loading, samples A and B are submitted to an oedometric compression test. Starting from the isotropic jamming state, the right-hand wall is progressively moved leftwards in quasi-static conditions (i.e. slowly enough to avoid inertial effects in the samples). This operation is performed both with the SDEM and with the Multibody Meshfree Method in MELODY. The whole process is illustrated in Fig. 11 in the case of Sample B. SDEM and MELODY results are qualitatively very similar, and the magnitude of the mean stress inside the grains is well reproduced by SDEM.

During loading, with both approaches, we observe short events of reorganization of the samples, and the grains mostly keep a roughly circular shape (apart from the contact segments). We observe in Fig. 11e, h that some grains seem slightly elongated in the vertical direction, both with SDEM Fig. 9 Zoomed views on isotropic compression tests as performed on samples A and B by MELODY (Multibody Meshfree Approach) and SDEM for ν=0.3 . Loading stages a to d correspond to marks in Fig. 10 ◂ and MELODY, but this effect remains limited. Hence, despite the strongly deviatoric character of the loading path, the stress state in the sample does not seem to deviate much from a spherical stress state. This observation is to be attributed to the absence of friction and to the deformability of the grains. These properties simplify their relative motions and allow them to relax their internal stress field through simple changes in the granular configuration of the sample. It is likely that the presence of intergranular friction will lead the grains to adopt shapes more similar to those represented in Fig. 1b, although this will have to be confirmed when a more advanced version of the method includes this feature.

Figure 12 provides quantitative results provided by SDEM and MELODY during the simulations of samples A and B.

The horizontal xx and yy in the sample are computed based on the forces applied on the rigid walls. This figure shows that the of these stresses during loading and their stabilization at the end of the loading are well-reproduced by SDEM when compared to MELODY, with a relative error which remains below 10%. This figure also shows that the horizontal (in the direction of loading) and vertical (perpendicular to loading) stresses have the same order of magnitude at any loading state, but that the horizontal stress remains slightly larger than the vertical one. It means that both samples can develop a certain amount of deviatoric stress. This is confirmed in Fig. 13, which plots the deviatoric stress xxyy ∕2 against the mean stress xx + yy ∕2 in both samples, as predicted by SDEM and MELODY. Data appear to be scattered, but some trends clearly appear. Both samples are able to develop a deviatoric stress, which increases with the mean stress. The relation between mean and deviatoric stress seems to follow a power law, although this would have to be confirmed by more careful and dedicated simulations. Both SDEM and MELODY predict that sample B presents a resistance to shearing slightly larger than sample A. In general, it appears that SDEM offers a satisfactory prediction of the oedometric response of soft granular samples, at least under the assumptions of the present study.

Conclusion and perspectives

The Soft Discrete Element Method introduced in this work is based on simplified kinematics of soft grains submitted to large contact forces. The main physical assumptions are a constant strain field within each grain (leading initially circular grains to deform into ellipses) and straight contact segments between contacting grains. To close the equations of motion corresponding to these kinematics, a simplistic contact model is proposed in the frictionless case, based on simulations of single-grain compaction. The whole problem is integrated in time by an explicit solver, and solutions are proposed to ensure stability and incompressibility and to introduce damping. The general approach is successfully validated in the cases of isotropic and oedometric compactions, by comparison with an analytical formula and a more accurate numerical approach, for various Poisson's coefficients and grains size distributions.

The conceptual simplicity of this method might make it possible to implement it in existing DEM codes. The main interest of this novel technique when compared to the current state of the art in soft grains modelling lies in the computation time. If we consider a single grain, a reasonable discretization in the Multibody Meshfree Framework requires about 100-150 field nodes, which corresponds to 200-300 degrees of freedom. The critical time step in such a system is driven by the nodes located at the boundaries of the grains, which are characterized by a low mass (because discretization is generally refined in the vicinity of the grains contour) and large stiffness (because the non-interpenetration condition with any contacting body requires a large contact stiffness). It is therefore necessary to adopt very low time steps to ensure stability of the solver. On the other hand, the Soft Discrete Element Method only considers six degrees of freedom per grain, and can be integrated in time with much larger time steps since the contact stiffness is comparatively rather low (Eq. 57). After implementation in an efficient and optimized code, CPU costs are thus expected to drop by 2 to 4 orders of magnitude when using SDEM instead of more finely discretized techniques. Obviously, this gain comes at the expense of a reduced accuracy at the local scale and of a questionable relevance in fully dynamic conditions.

A basic version of the method was presented in this paper, but it is important to stress that this framework is particularly versatile. Future extensions of SDEM may indeed include more complex constitutive behaviours for the grains deformability (such as hyperelasticity, plasticity, etc.) and more complex contact models (friction, cohesion, capillary adhesion, etc.). Kinematics could be either enriched or degraded: if we call the present approach the D2F3-SDEM (2 dimensions, 3 DoF in deformation), standard DEM might then be called D2F0-SDEM. Other versions such as D2F1-SDEM (volumetric deformation only) or D3F6-SDEM could be imagined. A 3D extension seems indeed possible, since the computation of stress and strains tensors will be straightforward. It will however require specific algorithms for rapid computation of large interpenetrations between ellipsoids, and a calibration of contact laws with an accurate 3D code including contacts, in the manner of Fig. 7. A natural and simple extension would also consist in considering ellipses instead of discs, i.e. to modify the elastic state of reference for each body. With such features, SDEM can be expected to have an interesting potential for upscaling simulations of soft grains samples in granular science, complex fluids rheology, soft matter physics, geophysics, and tribology.

Fig. 1

 1 Fig. 1 Close views on numerical experiments performed on soft grains with the Multibody Meshfree Approach; a Frictionless isotropic compaction; b Isochoric deviatoric loading at maximum compaction with friction

Fig. 3 a

 3 Fig. 3 a, b Difference in contact conceptualization between DEM (rigid discs with "point" contact represented by a very limited interpenetration) and SDEM (deformed grains with large deformations

Fig. 4

 4 Fig.[START_REF] Mollon | Can friction replace roughness in the numerical simulation of granular materials?[END_REF] Differences between the degrees of freedom θ (rigidbody rotation) and α (change in the orientation of the principle directions of strain) for ε 1 = 1 and ε 2 = 0 . Note in particular the motion of a given material point when each degree of freedom is continuously increased from 0° to 90° (black arrows)

Fig. 5 a

 5 Fig. 5 a Change of coordinates from the reference to the current configuration =-30 • , = 75 • , 1 = 0.3, 1 = 0.1 ; b Local contact quantities computed based on piecewise-linear contours

Fig. 7

 7 Fig. 7 Calibration of the contact model; a Simulation of single grain isotropic compression; b Fields of mean stress in the grain at different compression stages and for Poisson Coefficients of 0.4, 0.2, and 0 (from left to right); c Normal force in each grain-wall contact as a

Fig. 10

 10 Fig. 10 Comparison with Multibody Meshfree Approach numerical solutions. Marks a to d are loading stages represented in Fig. 9

Fig. 11

 11 Fig. 11 Oedometric compression of sample B with SDEM: a Initial state (isotropic jamming point); b Final state; c-e Zoomed views (black square area in a and b) at t = 5, t = 15, and t = 25, respectively; f-h Same views of the same system simulated with MELODY

Fig. 12 Fig. 13

 1213 Fig. 12 Oedometric compression test results: a Strain-driven loading history; b, c Horizontal and vertical stresses in sample A; d, e Horizontal and vertical stresses in sample B

Acknowledgements

The authors acknowledge that this study contains original material, as a result of a purely academic study without any kind of private funding or conflict of interest. Its publication has been approved tacitly by the responsible authorities at the institutes where the work has been carried out.

Declarations

Conflict of interest The authors declare that they have no conflict of interest.