Simulation of the mechanical behaviour of coated particles using DEM‐BPM

In previous years, research into the use of granular media as crash absorbers in a ship double hull has led to an increased interest in further optimisation to improve their energy absorption. One way is to add a coating layer on the outside surface of these particles. This work provides an introduction to the numerical modelling of coated particles to simulate their crushing behaviour. Depending upon the type of coating material, the compressive behaviour differs. Therefore, a sound methodology needs to be developed. For that purpose, the discrete element method in combination with the bonded particle method was utilised using the open‐source code MUSEN. The bonded particle method helps to simulate the crushing behaviour of coated particles. However, for this high‐fidelity model, the number of parameters needs to be reduced and afterwards optimised. A sensitivity analysis is performed for that objective followed by optimisation using particle swarm algorithm. The results illustrate a good approximation of the experimental force‐displacement curve by the optimised model on a uniaxial single particle compression test.


INTRODUCTION
The collision of ships is still a major contributor to shipping accidents [1].To prevent the loss of cargo and damage to the environment, for example in case of oil spillage, the crashworthiness of ship design therefore plays an important role.For that purpose, previous studies have shown that granular media can be introduced inside the ship double hull [2].Expanded glass, Poraver R , in particular was found to be a suitable candidate [3].The breakage of particles helps to absorb the kinetic energy from the collision.Additionally, the particles enable the transfer of load from the outer to the inner hull.Consequently, the whole structure works like a sandwich plate.To further optimise the energy absorption and particle breakage under dynamic loads, a coating layer made of environmentally friendly material properties was added and the mechanical behaviour of coated and uncoated particles was compared [4].The results show differences in breakage behaviour and energy absorbed depending on the type of coating.To simulate this behaviour independent of the coating material, a robust methodology needs to be developed.For this objective, discrete element method (DEM) combined with the bonded particle method (BPM) is used [5].However, the increased number of parameters due to the presence of both coating material and Poraver R particles adds to the computational cost and time.To overcome this problem, the number of parameters needs to be reduced and optimised.The methods used to perform these simulations are described in Section 2 followed by the corresponding results in Section 3.

MATERIALS AND METHODS
Each Poraver R particle has its unique shape and microstructure which adds to the complexity in determining its modelling parameters.Compression of each particle results in a different stress-strain relationship.Nevertheless, they can be divided into different diameter fractions to illustrate a better distribution of mechanical properties [6].However, the scattering of data remains.The addition of an environmentally friendly coating on Poraver R increases complexity.Although it does improve the energy absorption capabilities of the particles to an extent, the scattering of results still remains [4].Figure 1 highlights the microstructure of a wax coated Poraver R particle.The thin wax coating layer is approximately 0.25-0.5 % the diameter of Poraver R .Furthermore, it is a continuous medium in comparison to the porous microstructure of Poraver R .The type of coating material affects the mechanical behaviour as well [4].This is illustrated in Figure 3 where the wax coating is brittle and breaks away from the Poraver R particle.In contrast, the silicone coating is elastic and keeps the particles held together.
In order to simulate the mechanical behaviour of the coated particles under compression, an open-source code MUSEN is used [7].It allows the macroscopic modelling of particles using DEM where each particle is a discrete quasi rigid body F I G U R E 3 Different mechanical behaviour under uniaxial compression of coated particles.Wax coated Poraver R particle a and silicone coated Poraver R particle b. [8].Once the neighbouring particles have been determined, the Hertz-Mindlin contact model is used to generate contact forces.These are fed into Newton-Euler equations to calculate the velocity and displacement of each particle using an explicit time integration scheme.Additionally, MUSEN allows the extension of DEM with BPM.This enables the particles to be connected to each other via cylindrical bonds forming an agglomerate.To allow for conservation of mass, these bonds are weightless and are removed from the simulation domain once they are broken.In order to simulate the bond breakage, the stress-strain relationship in each bond is calculated using a material model.This allows the user flexibility in modifying well known existing material models, such as elasticity or elastoplasticity, to the required application.To simulate the micro-cracking behaviour of coated particles, a micro-cracking elastic-plastic material model is used.This particular model is well suited to simulate the micro-cracking in Poraver R particles and the energy dissipation caused by them [9].
Figure 2 depicts the cross-section of the numerical model of wax-coated Poraver R particle.The wax particles on the outside are coloured with red and the larger Poraver R particles on the inside are in blue.As mentioned above, the particles are connected to each other with bonds.Wax particles are bonded to themselves and Poraver R particles in light red whereas the Poraver R particles are connected in light blue.In comparison to the actual microstructure in Figure 1, the numerical model is approximated on the mesoscale.It does not take into account each and every pore which might have been possible using a CT scan.This simplification helps in reducing the computational cost and time later on.Nevertheless, it still provides a good approximation of the stress-strain relationship of each particle, as seen in the later sections.
A total of 43 different parameters (geometric and material) need to be defined in MUSEN to simulate the numerical model shown in Figure 2.For simplification, the model for Poraver R particles described in previous works is adapted [9][10][11].In addition to this simplification, the parameters related to the geometry of the wax coating layer are also fixed and adapted in a similar way.The thickness of the layer is derived from the observation of microscopic scans [4].Furthermore, the bond model between wax-wax particles and wax-Poraver R particles is assumed to be similar.Therefore, these assumptions lead to a reduction of parameters from 43 to 17.A sensitivity analysis is then performed on the remaining parameters to determine their influence on the breakage force and breakage displacement when changed by a factor 10×.This results in seven sensitive parameters, as shown in Figure 5.It includes the breakage strain for wax -Poraver R bonds (BS w -P b), Young's modulus for wax -Poraver R bonds and particles (YM w -b/P) and sliding friction for wax -Poraver R particles (SF w -P), wax -steel plate (SF w -s) and for wax -wax particles (SF w -w).
For the purpose of optimisation, the particle swarm algorithm as part of the MATLAB optimisation toolbox is used.The difference between the experiments and the simulation was used to define an objective function where  exp and  sim are vectors containing the function values from experiments and simulation at similar displacement values, respectively.The optimisation is done on experimental results from uniaxial single particle compression test found in [4].Due to computational costs, the optimisation is done on a single particle.Figure 4 highlights the difference in computational times for coated and uncoated particles on a GPU based on the results from single particle compression F I G U R E 4 Computational times for coated and uncoated particles for uniaxial compression of single particle (SPT) and multi particle (MPT).
F I G U R E 5 Variation of sensitive parameters by a factor of 10×.test (SPT) on GeForce RTX 2080TI, 11 GB GDDR6 memory, 4352 CUDA cores and multi particle compression test (MPT) on GeForce RTX 3090, 24 GB GDDR6X memory, 10496 CUDA cores.Firstly, the simulation times are higher for coated particles for both single and multi particle compression by factor 2.8× and 6.7×, respectively.Secondly, for each particle type, the simulation times are significantly higher for the multi particle compression.This leads to the realisation that optimisation is to be done on simple particle compression based on the high-fidelity numerical model depicted before in Figure 2.

RESULTS AND DISCUSSION
The results for uniaxial compression for the particle diameter between 2.0 to 2.15 mm are shown in Figure 6.The experimental results depict a significant variation for each test run since every particle has its unique microstructure.Consequently, two experimental results can never be repeated.However, an average curve, highlighted in colour red, can still be derived based on the mean breakage force and mean crushing displacement.This experimental curve is fed to the objective function and compared with the simulation curve in blue, generated by the parameters optimised by the optimisation algorithm mentioned in Section 2. The optimised set of parameters is listed in Table 1.The result of this Difference between experimental and optimised simulation force-displacement curve for wax coated Poraver R particles of diameter 2.0-2.5 mm.

TA B L E 1
Optimised set of parameters for a wax coated Poraver R particle with diameter 2.28 mm.set shows a good agreement between the experimental and simulation curves.Similar to the experiments, it does depict microcracking which is explained by the small dips in force before the main particle breakage characterised by the drop in force to less than 1 N.However, there are differences.Up to 0.05 mm displacement, the slope of the curve remains constant.Afterwards, a drop in stiffness is observed.This is explained by the fact that the numerical model assumes the coating layer to be a set of bonded particles whereas in reality, it is a continuous layer.As the coating layer in the numerical model breaks, it results in a loss of stiffness till the steel plate further moves downwards and comes into direct contact with the Poraver R particles.This means that majority of resistance to compression is provided by the the Poraver R particles.

Sliding friction coefficient -
Nonetheless, the simulation model provides an overall reasonable approximation to the experimental averaged curve.Similar methodology can be repeated for other diameter sizes of coated particles.One example is depicted in Figure 7 with its optimised set of parameters listed in Table 2.In this case, the larger particle size leads to a higher mean breakage force and mean crushing displacement.Consequently, the averaged experimental curve differs from the previous one.Nonetheless, it can still be fed into the optimisation algorithm in a similar way to generate the simulation forcedisplacement curve.Similar to the previous result, the initial stiffness up to 0.1 mm displacement matches the average curve.Subsequently, there is a drop in the value of slope of the curve as the bonds are broken between the coated and Poraver R particles.The slope starts to increase again once more Poraver R particles come into contact with the steel plate up to a point where a significant portion of bonds are broken.
The set of optimised parameters for both results are also comparable.The sliding friction coefficient is more than 0.5 for all material interactions for both outcomes.Moreover, the differences between the two results are negligible.Therefore, F I G U R E 7 Difference between experimental and optimised simulation force-displacement curve for wax coated Poraver R particles of diameter 3.125 -4.0 mm.

TA B L E 2
Optimised set of parameters for a wax coated Poraver R particle with diameter 3.48 mm. the sliding friction coefficient parameter does not play a significant role in the compression of a single particle.However, there are differences in the values for Young's modulus and breakage strain.In contrast to the smaller particle, the larger particle has approximately twice as high values.This is reasonable since the higher values help to compensate for the higher breakage force and crushing displacement.

Sliding friction coefficient -
Once the optimisation is done on a single particle, the next step is to test the optimised numerical model in a multi particle compression test setting to determine its suitability in simulating the actual goal of simulating the compression of a ship double hull.The setup for the test is similar to experiments [4].The results of compression of wax-coated Poraver R particles are depicted in Figure 8.In contrast to the single particle compression, the variation between experimental results is less.This means if a certain number of agglomerates are present inside the cylinder, the results remain consistent.The total average number of coated particles in the experiments was 2306.To replicate this, a same number of optimised agglomerates were randomly generated inside the cylinder with similar geometric and material parameters.However, the simulation crashed as soon as it began.This is due to the lack of GPU memory available on the GeForce RTX 3090 graphics card.However, in a second attempt only 1300 agglomerates were generated inside the cylinder which is more than 50% of the average.The results of this simulation are shown in blue colour with the stress and strain normalised to take into account the lack of particles.During the initial part of compression till 15% of strain, the slope of the simulation curve matches that of the experiment.However, it starts to increase exponentially beyond this point.This is explained by the lack of unbroken bonds present inside the model to allow for energy absorption through bond breakage.Consequently, the stress increases as the particles inside are densely packed.This shows that the high-fidelity model of a F I G U R E 8 Uniaxial multi particle compression test of wax coated Poraver R particles.
single particle is able to capture the initial deformation of coated particles in a multi particle compression test.However, to model large deformation with strains greater than 15%, more number of agglomerates need to be present inside the cylinder.Unfortunately, that approach does not work due to the lack of computational power available.

CONCLUSION
Based on the results described above, the optimised model based on DEM in combination with BPM provides a reasonable approximation of the compressive behaviour of a coated particle.Although the physical microstructure is different from the numerical mesoscopic model, it is still able to predict particle breakage, especially the small microcracking behaviour due to the presence of bonds.Additionally, the methodology can be repeated for other diameter sizes of coated particles.However, the real suitability of this model is seen during a uniaxial multi particle compression test.The current highfidelity model is limited due to lack of computational power available.Therefore, a high fidelity model in combination with a simplified model needs to be developed.This will decrease the computational cost as well as provide the ability to simulate the numerical model on an actual ship double hull simulation.

A C K N O W L E D G M E N T S
The

F I G U R E 1 F I G U R E 2
Microstructure of a wax coated Poraver R particle with 500× magnification.Agglomerate consisting of primary particles of wax (red), wax-wax-Poraver R bonds (light red), Poraver R particles (blue), and Poraver R -Poraver R bonds (light blue).
presented research is funded by the German Research Foundation (Deutsche Forschungsgemeinschaft, DFG) in the framework of the research training group GRK 2462 "Processes in natural and technical Particle-Fluid-Systems" (PintPFS), which is gratefully acknowledged.Open access funding enabled and organized by Projekt DEAL.Open access funding enabled and organized by Projekt DEAL.
R E F E R E N C E S