The effect of particle shape on porosity of swelling granular materials: Discrete element method and the multi-sphere approximation

Theporosityvalueofswellinggranularmaterialsdependsonthedegreeofswellingandthestate-of-stressinside the packing. This research focusses on the swelling of a bed of Super Absorbent Polymer (SAP) particles, which are absorbent particles that can absorb 30 to 40 times their initial weight of saline ﬂ uids. Due to the complexity of swelling materials, measurements of porosity, and also hydraulic parameters, are complex and often not feasible, unless a quasi-static approach is employed. An alternative is not only to conduct measurements, but in ad- dition to generate arti ﬁ cial packings using a grain-scale model from which porosity values can be derived. The combination ofexperimentsand simulationsreducetheexperimentalefforts whilst allowing for anunderstand- ing of particle-scale phenomena. A suitable grain-scale model is the Discrete Element Method (DEM). In DEM, particles are often represented as spheres, while SAP particles have very irregular shapes, which strongly affect their hydraulic properties. Therefore, the shapes of particles should be included in DEM simulations, which is achieved by using sets of overlap- ping spheres (often referred to as clumps). In this research, an algorithm was employed to make realistic clumps that are representative of measured SAP particles, using micro-CT scans of individual SAP particles. In DEM, particles were randomly generated based on 20 clumps and the particle size distribution of real SAP parti-cles.Theparticleswerethencompactedtoobtainapacking,fromwhichtheporosityvalueswereobtained.These porosity values were close to experimental values and thus we concluded that our simulations of particle shape in DEM does improve the quality of simulations in case of granular materials. Finally, the relationship between porosity, state-of-stress and the degree of swelling is developed based on our grain-scale model.


Introduction
The pore geometry inside a granular material depends on the arrangement of particles, particle sizes, and particle shapes. As a result, both the properties of particles and the particle packing determine the bulk mechanical properties, porosity values and hydraulic parameters. For example, the permeability of a particle packing decreases with decreasing particle size and decreasing porosity value, which indeed is the basis of the Kozeny-Carman equation that relates particle size and porosity to permeability (e.g. [2]). On top of that the permeability is lower when particles are irregularly shaped than when particles are spherical, for particle packings that have the same porosity value and same particle size distribution [11,36]. Moreover, the entry pressure of air to invade a saturated granular material increases with decreasing porosity values or decreasing particle sizes (e.g. [28,35]).
Regarding non-deforming granular materials, the pore geometry and the corresponding hydraulic parameters can be obtained experimentally. Column experiments can be conducted to measure various hydraulic properties while the pore space can be visualized using for example X-ray tomography or micro-CT. However, measuring those properties is difficult or impossible for granular materials that are being deformed over time, unless a quasi-static approach is employed, where deformation is applied in small discrete steps while measuring various properties.
An example of a fast deforming granular material is the swelling of Super Absorbent Polymer (SAP) particles, which are absorbent particles that are frequently used in hygienic products (e.g. Buchholz and Graham [4]). SAP particles can absorb up to 1000 times their initial weight within a timeframe of b10 min [4,8,17]. During this process, SAP particles soften and increase in size [18]. The characterization of a bed of swelling SAP particles is required to enable continuum-scale modelling, as described in Diersch et al. [6]. For example, the following three parameters all change with swelling and thus with time: porosity, permeability, and the capillary pressure-saturation curve. Therefore, experimental characterization may be complex. An alternative is to employ a grain-scale model that can obtain dependencies of hydraulic parameters on the swelling degree, but also on the state-ofstress.
A suitable grain-scale model is the Discrete Element Method (DEM), which can simulate the movement of individual particles in a packing during deformation [34]. DEM can be used to construct packings of spheres that represent specific granular materials (e.g. [3,5,37]), by matching the particle size distribution and either the porosity or the state-of-stress of that material. The hydraulic properties can be measured independently from DEM simulations, by employing a variety of techniques. For example, the capillary pressure-saturation relationship can be constructed using the Pore Unit Assembly (PUA) method [12,24,25,35,38] or by using the pore morphology method [14,39]. The permeability can be found by solving for water pressure using the pore finite volume method [5] or using direct simulations by solving Stokes equation (e.g. [5,11]); the permeability can be computed afterwards using Darcy's law.
For swelling SAP particles, Sweijen et al. [35] have constructed capillary pressure-saturation curves for various porosity values and degrees of swelling, using packings of spheres generated in DEM. The largest porosity value obtained was approximately 0.43, which corresponds to a loose packing of spheres. However, Mirnyy et al. [26] found a maximum porosity value for dry SAP particles of 0.90 and Diersch et al. [7] reported a porosity value of 0.70. These large porosity values are most likely a result of the irregularity in particle shapes.
Therefore, in this research, we study whether the inclusion of particle shape in DEM would improve the estimation of porosity of generated packings. Particle shapes are included in DEM using sets of overlapping spheres that act as one particle, which are referred to as clumps (see e.g. [9,23]). Other methods exist in DEM to account for particle shape, in context of particle mechanics. For example, Höhner et al. [15,16] employed polyhedra as particles to account for the angularity of particles, in case of a rotating drum. Pasha et al. [40] studied the effect of particle shape on particle flow in a rotary batch seed coater, they found that using spherical particles that have a rolling friction yields similar results to as using clumps. However, the shape of particles is important for determining hydraulic parameters [36]. Therefore, we generated various packings of clumps that represent packings of real SAP particles. Ideally, the shape of each SAP particle in a packing should be considered, which is not feasible for packings that typically have N1000 particles of usually slightly different shapes. Therefore, we constructed a set of 20 clumps in a random fashion, based on shapes of real SAP particles. From 20 SAP particles, a 3-dimensional reconstruction was made using micro-CT, to which clumps were fitted using a newly developed algorithm based on that by Price et al. [30]. The set of clumps was used to construct various packings in DEM for different friction angles and confining stresses and for each packing the porosity was determined. In future work, the hydraulic properties of the packings can be determined, such as permeability, relative permeability and capillary pressure-saturation curves.
The aims of this research are as follows: i) to construct packings of SAP particles in DEM, using clumps; ii) to obtain the relation of porosity as function of the degree of swelling, iii) to study the effect of confining stress and how this can be parameterized in terms of fitting equations and iv) to compare porosity values of simulations with those of experiments.

Realization of clumps
Each particle has its own unique shape and surface irregularities, which is a result of the production process of these SAP particles. Typically, SAP is produced in large blocks which are then ground to form particles having an irregular shape [13,29]. Including all unique particle shapes of a particle packing having N1000's of particles is not feasible during simulations, because it would require an excessive amount of time for particle characterization. Therefore, a library of 20 particles is generated, in a random fashion, whose shapes are representative for most of the particles. Thus each shape represents a class of particles, whilst small heterogeneities of the shapes are ignored.
For each particle, a set of overlapping spheres is constructed, to which we refer to as a clump. A clump represents the main shape and features of a class of particles (see Fig. 1); e.g. cubical, elongated, or plate-like particles. Note that the aim is not to fully reproduce a measured particle shape, but that an averaged approximation of the shape is deemed sufficient (i.e. micro-roughness, small tight corners are ignored). When constructing a packing of SAP particles, shapes are chosen randomly from the library (following a uniform distribution) and are given a size following the particle size distributions (usually known from experiments, this typical entails a Gaussian distribution).
In this work, we follow the procedure and algorithms as in Sweijen et al. [33], which is based on the work by Price [30]. More comprehensive algorithms exist in the literature to generate clumps (see for example: [10,21]), however these algorithms typically result in very detailed reconstructions of particles, which would result in too many spheres per clump. In what follows, the methodology of generating a clump library is explained.

Imaging of individual SAP particles
Individual SAP particles were scanned using a Micro-CT setup. First, SAP particles were put into a small petri-disk in one single layer such that no SAP particle touched another particle. Then, a stack of 2dimensional images was taken with a 6 μm resolution. The images were post-processed using image analysing software Aviso (FEI, 2016). The images were converted to a 3-dimensional image of the layer of SAP particles from which 20 individual SAP particles were randomly chosen (based on a uniform distribution) and their shapes were extracted into surface plot files. The surface plots were subsequently used in the algorithm to construct multi-sphere approximations, as described in the following section. Previous studies on particle shapes have employed similar workflows, for example, Latham et al. [20] used a 3-dimensional laser setup to acquire 3-dimensional reconstructions of pebbles, with a resolution of 100 to 500 μm.

Creating a library of clump shapes
A fitting algorithm was developed based on the work by Price et al. [30]. The aim was to find a set of overlapping spheres whose outer surface would lie close on the surface of the particle, but where the number of spheres to represent the shape was minimal to minimize computational time in DEM. For example, an elongated particle should be represented by an elongated set of spheres.
The first step of the fitting algorithm was to import the surface mesh of a SAP particle, which was formed by a web of triangles. Next, four vertices on the surface of that particle were selected randomly (based on a uniform distribution) and a sphere was fitted to pass through those four points. Then, a selection procedure was employed to check whether the fitted sphere was admissible or not.
For the sake of simplicity, we explain the selection procedure in two dimensions, where the shape of a particle is shown as a closed curve (see Fig. 2), which is a string of points instead of a web of vertices (in three dimensions), and we use circles (instead of spheres) to represent a clump. Thus, circles will be chosen to fill out the particle shape such that their outer edges coincide with the particle boundary as much as possible. From the string of points, we randomly select three points (or 4 points in 3-dimensions) based on a uniform distribution, fit a circle to them, and check whether the circle is admissible.
If the selected points are on an almost straight line, the resulting circle will be much larger than the particle (circle I in Fig. 2a). To avoid this scenario, we require the radius of a fitted circle to be smaller than the smallest principal axis of the particle. A fitted circle can stick too far out of the particle (circle II in Fig. 2a); this is quantified in terms of bulginess. The bulginess is defined as the distance between a point on the perimeter of a particle, which falls inside a fitted circle, and the nearest point on the boundary of that circle. The maximum allowable bulginess for a fitted circle is defined to be a fraction (α) of the radius of that circle. Obviously, the bulginess can be bigger for larger spheres, which results in a tendency to generate large spheres rather than many small spheres. A circle can also be dismissed if its centre falls outside the particle (circle III in Fig. 2a). To identify whether a centre of a circle k with coordinates X k is located inside a particle, the following statement was employed which has to be true for all points on the boundary of a particle: withn being the outward pointing normal vector of a surface element and X i being a point on the boundary. Note that Eq. (1) disqualifies spheres at convex curvatures, which is beneficial to avoid small spheres representing small spike-like features of the particle shape. In addition, to reduce the number of circles (in fact spheres) of a clump for irregular particles, a small fraction (φ) of them are allowed to have their centres outside of the particle, provided they comply with the maximum allowable bulginess. Once an admissible circle had been found (similar to circle IV in Fig.  2a), the points on the surface of the particle that fall within that circle were blocked for further usage. Then, other circles were fitted to the remaining points. Obviously, it is not possible to have all points on the surface of a particle to fall within admissible spheres; many points may still fall outside the spheres (see e.g. Fig. 2b). We required a pre-specified fraction ϑ of the total number of points to fall within admissible spheres. This completed the generation of a clump of spheres. However, it is still possible to have the clump to fit the particle more closely. This was done by adjusting the radii of all spheres to have their outer surface as close as possible to the surface of the particle. Thus, we minimized the distance between the surface of the particle and the outer surface of the set of overlapping spheres. The centres of spheres were considered fixed, to avoid areas with a large vertex density (number of vertices per surface area) to have more weight during fitting of a sphere than areas with a low density, for example, to avoid that highly irregular areas have more weight than smooth areas. Vertices i on the surface of the particle that fall within a sphere k were assigned to that sphere. Vertices i on the surface of the particle that did not fall within any sphere, were assigned to the nearest sphere k (see Fig. 2b). The total number of vertices that were assigned to sphere k was denoted as N ik . The shortest distance between a vertex and its corresponding sphere was defined as d ik . The following sum of residuals (Θ k ) was calculated for each sphere: The optimization algorithm was as follows: 1) for each vertex, the closest sphere was identified; 2) for each sphere, the radius was incrementally increased from 0 to a maximum value while calculating Θ k ; 3) for each sphere, the lowest value of Θ k gave the new radius of the sphere; 4) after updating all radii, step 1, 2, and 3 were repeated until X k Θ k converged to a minimum value. Finally, the clump was checked for internal porosity, which did not occur for the clumps generated in this work. One set of values for φ, α, and ϑ was used to generate all clumps. We used the values: φ = 0.025, α = 0.05, and ϑ = 0.85 as those gave the best fit of the volume-to-area ratio of clumps to that of particles, which is evaluated in Appendices A and B. Fig. 1. Illustration of: a) one SAP particle, b) one SAP particle and its corresponding clump realization, and c) the clump realization. Fig. 2. Two-dimensional schematic overview of the fitting procedure: a) possible scenarios that may occur during the fitting of a sphere in a particle; I) a circle that is too large, II) a circle with too much bulginess, III) a circle that has its centre outside of the particle, and IV) a circle that is considered a good fit and b) overview of the optimization procedure.

Multi-sphere approximation in the discrete element method
Here, we use open-source software Yade-DEM [32] which is a 3dimensional model based on the discrete element method (DEM). To model contact mechanics, Hertz-Mindlin contact mechanics are employed, which has previously been implemented in Yade-DEM by Modenese et al. [27].
To include the shape of real SAP particles, we employed the multisphere approximation. The advantage of using the multi-sphere approximation is that the sphere-to-sphere contact detection is the same as for spherical particles, which is faster than that of irregularly shaped particles [9,19]. In the multi-sphere approximation, the sub spheres of a clump can have contacts with other clumps giving rise to contact forces acting on the sub spheres. The contact forces and torques are then projected to the centre of mass of the clump. For the centre of mass, the force balance is summed and integrated to find its acceleration. The acceleration is subsequently projected back to the sub spheres of the clump in order to have the relative positions of the sub spheres constant during the movements of the clump. The reader is referred to Šmilauer et al. [32] for more information on the implementation of clumps in Yade-DEM.

Swelling of clumps
Swelling was included by increasing the particle sizes and by softening the particles. In this study, we assume that the shape of particles remains constant while they grow due to swelling. This assumption is a simplification of the physical process of swelling, which is a complex dynamic process. In experiments, the initial stage of swelling is characterized by particles having a cauliflower-like shape because the edges will swell faster than centre of the particle. However during swelling the initial shape is obtained again (due to isotropic swelling of the material). In addition, we assume that the stresses at the particle contact points are small compared to the stresses involved with swelling of a particle; thus the stresses at the contact points do not interfere with the process of swelling. Of course, when a particle bed is constrained, the stresses at particle contact points will become large and will alter the shape of particles as they grow into the pore-space rather than pushing away other particles.
To define swelling, first the degree of swelling is defined in terms of the absorption ratio, Q i abs [4]: in which M i w is the mass of water inside particle i and M i s is the mass of dry SAP in particle i. Based on Eq. (3), we define a volumetric grow factor f v as: in which V i w and V i s are the water and dry SAP volume, respectively, and ρ w and ρ s are the water and dry SAP density, respectively. The initial location of spheres in a clump, denoted as a position vector Y k 0 , was changed according to f v such that the location of a sphere after swelling (Y k ) is given by: in which Y 0 is the position of the centre of mass of the clump. The radii of all sub spheres were updated by simply applying: is the initial radius of a sub sphere. Finally, to include the softening of the SAP particle during water absorption, the following equation was employed [4,18,34]: in which G i is the shear modulus of particle i and β is a constant. In DEM, changing the size and shear modulus of the particles will effectively result in only changing a stiffness that is used in the Hertz-Mindlin contact mechanics. Thus by using DEM, we study the effect of stiffness on porosity, where the stiffness is determined by the shear modulus and particle size. The shear modulus and particle size are related to experimental work as they are both physical properties.

Clump insertion into the modelling domain
To generate a dense packing of clumps, first a loose packing was generated by randomly inserting clumps into a large cubic modelling domain of 5 × 5 × 5 cm 3 . Clumps were inserted into the model domain as follows: i) a clump was randomly chosen from the shape library, ii) it was resized to match a randomly chosen size from the particle size distribution, where the smallest transect of a clump was used as particle size, iii) the clump was resized and its shear modulus adjusted according to a predefined absorption ratio as described in Section 3.2, and iv) the clump was randomly rotated and randomly placed in the modelling domain, such that it did not touch any boundary or other clump. The porosity of the cloud of clumps was typically larger than 0.90. For small particle sizes, defined as particle sizes smaller than half of the D50 of the particle size distribution, a sphere was inserted rather than a clump. This was done because small irregularities of small particles were assumed to have no effect on the mechanics of the packings.
To compact the cloud of particles, we applied a confining pressure of 0.3 Psi (2068 Pa) to all boundaries of the modelling domain. All boundaries were periodic. The model was run until force equilibrium was achieved inside the packing. Then the porosity value was determined by converting the packing of clumps to a 3-dimensional voxel image, from which the porosity was determined.

Parameterization
The mechanical parameters were chosen to resemble those of SAP particles. The shear modulus was set to 7100 Pa at Q i abs = 30g/g; i.e. β was 22.1 kPa, following Knaebel et al. [18] and in Sweijen et al. [34]. The Poisson ratio was chosen to be that of incompressible water, namely 0.5. The friction angle was set to that of acrylic acid, namely, 5.5°following Lorenz et al. [22], but the friction angle was also varied from 0°to that of Labenne sand, namely 35° [3]. The particle size distribution was set to that of SAP particles, which were the same batch as measured using the micro-CT setup in Section 2.1; D50 was set to 460 μm and a standard deviation to 69 μm.

Experiments on a bed of particles
Experiments were performed on a bed of SAP particles to obtain the porosity as function of the absorption ratio, for a constant confining stress. Two pressures were used, namely: 0.3Psi (2068 Pa) and 0.015Psi (104 Pa). First, a predefined mass of dry SAP particles was put in a column having a diameter of 6 cm. At the bottom of the column, a filter was located that kept the particles in place while allowing water to pass through. At the top of the column, an impermeable weight was placed that corresponds to the predefined confining stress. The weight was able to move up or down freely. Then, the column was put in a container with water, such that water could flow into the particle bed. Then, after 4 h, the column was taken out of the water and the total mass of the swollen particles and the interstitial water was measured. Then, the excess water was removed gently using a paper towel. This removed all interstitial water and only left swollen SAP particles. The mass was then measured again, which allowed for the determination of the porosity.

Minimum number of particles and shapes
Prior to any other simulations, the influence should be determined of the number of particles and shapes on the porosity. In other words, the porosity value should be independent of both the number of clumps and shapes and therefore a minimum number of particles and shapes is required. Therefore, we generated packings of 4000 dry SAP particles (Q abs = 1 g/g) for various numbers of shapes, ranging from 1 to 25 (note that 5 extra random clump realizations were made based on 5 particles from the measured 20). Once a packing was generated, various porosity values were measured for varying subdomains in a packing. The subdomain was gradually increased from 10% of the dimensions of the modelling domain to 90% of the dimension. The porosity value of the subdomain was plotted as function of the number of particles in that subdomain (see Fig. 3). Fig. 3 shows that increasing the number of particles resulted in the porosity value converging to a constant value, which was typically obtained after 1500 particles. To ensure that the reported porosity values are independent of the number of particles in a simulation, we used approximately 3000 particles to determine the porosity value. By increasing the number of shapes per simulation, the porosity value also converged to a single value for N10 shapes. Therefore, 20 shapes were used, making the simulations independent of the number of shapes. Note that the minimum number of 10 shapes is specific to this study, for the shape of SAP particles and particle size distribution. Other particles, having different particle size distribution and particle shapes could have a different minimum number of particles. Thus, in all subsequent simulations, 4000 particles and 20 shapes were used.

The effect of absorption ratio on the porosity value
Using a confining stress of 0.3Psi (2068 Pa), different packings were constructed for various absorption ratios [1-40 g/g]. The conceptual experiment for these simulations was as follows: 1) SAP particles were swollen to a given absorption ratio; 2) the SAP particles were put into a triaxial cell; 3) the particles were compacted using an isotropic compaction; 4) the porosity of the packing was calculated. Using a triaxial compaction on already swollen particles, a simple packing was generated that had no particular loading history (i.e. the degree of anisotropy is minimalized). But the packing was considered to be dense as triaxial compression enhances rotation and interlocking of the clumps. Fig. 4 shows that the porosity value monotonically decreases with increasing absorption ratio, which is caused by a decrease in shear modulus with increasing absorption ratio. This phenomenon was also reported by Mirnyy et al. [26] and Diersch et al. [7], who studied the swelling under a constant confining pressure. Fig. 4 also shows the effect of the friction angle on the porosity -absorption ratio curve. The porosity value decreases with decreasing friction angle, which is a wellknown phenomenon in the DEM literature (e.g. [31]). Friction angles smaller than or equal to 5.5°gave results that correspond well with experimental data. In particular, a friction angle of 1°gave a good match. However, porosity values are overestimated when using the friction angle of sand (35°). We therefore conclude that the friction angle for SAP particles should be relatively low (≤5.5°), which is in the order of magnitude of the friction angle of acrylic acid, 5.5° [22].

Mechanical considerations
Swelling and deformation affects the structure of the particle packing, as shown in Fig. 5. Dry SAP particles, under a loading stress of 0.3Psi, form a relatively loose packing of particles with a porosity of 0.39. In Fig. 5a, it is clearly visible that high porosity zones exist and that the pores form an open and easily accessible network for water to  flow. However, high particle density zones exist for swollen particles (see Fig. 5b). Pores in Fig. 5b show less connectivity than the open pore geometry of dry SAP particles, which is a result of more particleparticle connections and larger contact areas in packings of swollen particles. This is required to support the applied load for relatively soft swollen particles.
SAP particles are relatively soft compared to the confining pressure; i.e. the shear moduli of SAP vary between 22 kPa to 6.5 kPa during swelling, while the confining stress is 2068 Pa. In other words, the shear moduli are 10.6 times and 3.3 times higher than the confining stress, respectively, thus the particles are relatively soft. Therefore, contact mechanics may become non-linear during swelling. Therefore, we compute a threshold value for the absorption ratio (Q linear ) above which the contact mechanics is expected to become non-linear. Following Hertzian contact mechanics, the stress-strain (σ − ε) relationship at a contact point is given by (Popov [42]; Johnson [41]): where σ is the confining stress, ν the Poisson ratio and E is the Young's modulus, which is related to the shear modulus G by E = 2G(1 + ν). Using Eq. (6), we can rewrite Eq. (7) to: In this work, we assume linear elasticity to be valid for strains up to 1% (i.e. ε b 0.01). Therefore, we evaluate Eq. (8) at ε = 0.01, β = 22.1 kPa, σ = 2068 Pa, which results in Q linear = 20.2 g/g. For Q abs lower than Q linear the deformation at the contact points is linear (see dashed line in Fig. 4), where we assume that the stress at contact points does not exceed the mean stress in a packing too much. Thus, for values of Q abs below 20.2 g/g the porosity predictions are reliable whereas for larger values of Q abs the porosity prediction may not be reliable.

Gravity effect in case of small confining pressures
At relatively low confining pressures, the stress inside a particle packing is also significantly affected by the weight of the SAP particles, which is not included in the DEM simulations as the modelling domain is only a small fraction of the experimental volume. Therefore, the effect of overlaying particles should be considered in the confining stress at the boundaries of the DEM simulations.
In Fig. 6 the porosity is shown as function of absorption ratio for a confining stress σ 0 of 104 Pa (0.015 Psi). In that Figure two types of simulations are presented, namely one where we only use 0.015 Psi (σ 0 ) as confining pressure and one where we use σ 0 and a correction term for the average overburden stress inside a particle bed (σ c ), which is the weight of the overlying particles (in the middle of the particle bed). The overburden stress is evaluated for a point in the middle of the SAP particle bed such that the confining stress (σ) at the boundary conditions becomes:   where ρ is the density of swollen SAP particles, which is that of water, and h is the height of the particle bed in experiments which is given by h(Q i abs ) = 7.1 × 10 −4 Q abs [L], which is valid for 3 b Q abs b 28 [g·g −1 ]. Fig. 6 shows that without using the correction term, the model overestimates porosity values from experiments; after all, the overburden is not included in DEM simulations whilst it is present in experiments. However, including the overburden stress results in a closer match in porosity values between simulations and experiments.
Albeit that the experimental data point at Q i abs = 5 g/g shows a different behaviour than expected from modelling results and experimental results for 0.3 Psi (Fig. 4). This phenomenon is attributed to the fact that the packing history in the model and experiments are not the same, which is true for all data points reported in this work. But at Q i abs = 5 g/g, the particles are stiffer and smaller than at higher swelling degrees, which make the packing prone to structural changes due to small vibrations during handling in experiments (which also explains the large error bar for this point). At higher swelling degrees or at larger loads, the sensitivy of particle packings to vibrations becomes less pronounced and therefore the measurments match better with simulations. Fig. 8. Porosity values of SAP particle beds in a dry-state, for three different shear moduli, using a confining stress of 2 kPa. The shear moduli of 22.1 kPa is in accordance with Eq. (6) while the shear moduli 1 and 100 MPa are estimated values for real SAP particles in drystate. Fig. 9. Porosity as function of absorption ratio for different confining stresses. Symbols represent data from clump simulations, the solid lines represent the fitting of Eq. (10) and the colours indicate the confining stress, using a) the shear modulus given by Eq. (6) and b) the shear modulus given by Eq. (6) but the initial porosity is set to 0.68. Note that for confining pressures 3000 and 4000 Pa, the porosity values should be considered as a indicative value as contact mechanics are showing non-linear deformation (following Eq. (8)) as well as absorption ratio's higher than 20 g/g for confining stress 2068 Pa.

Spheres versus clumps
To test the effect of shape, we have simulated SAP particle packings for varying absorption ratios (1 to 40 g/g) for both clumps and spherical particles, using a friction angle of 5.5°and a confining stress of 0.3 Psi. The particle size distribution and the solid volume was identical for both packings (spherical particles and clumps). Fig. 7 shows the difference in porosity values between spheres and clumps for a confining pressure of 0.3 Psi. The porosity value of clumps, ranging from 0.10 to 0.34, is lower than that of spheres, ranging from 0.21 to 0.34.
Thus, usage of clumps can result in lower porosity values than for spheres, in case of swelling particles. This phenomenon is a result of interlocking of irregularly shaped particles. However, Fig. 3 shows that the porosity increases with the number of particle shapes, indicating that interlocking does not occur for dry SAP particles. The explanation of dry SAP particles behaving differently than swollen SAP particles is that swollen SAP particles are much softer than dry SAP particles and as a result the confining stress is more significant for swelling particles than for dry particles. The confining stress forces swollen particles to compact, while the confining stress is very small compared to the stiffness of dry SAP particles and therefore has little effect on dry packings. Therefore, we observed an increase of porosity with the number of particle shapes (interlocking is not trigged for dry SAP particles), but for swelling particles the confining stress forces the particles to interlock.

Initial porosity
SAP particles in a dry-state, where Q i abs equals unity, have a much larger shear modulus than particles that are slightly swollen [1]. Therefore, the shear modulus predicted by Eq. (6), which is valid for swollen particles, is an underestimation of the shear modulus of dry particles. To investigate the porosity value of particles in a dry-state, we have constructed packings with a shear modulus of 1 MPa and 100 MPa for a variety of friction angles. Fig. 8 shows the porosity value as function of friction angle, for a 2 kPa confining pressure. A shear modulus of 22.1 kPa, which is predicted by Eq. (6) and using β of 22.1 kPa, results in lower porosity values than when using a shear modulus of 1 MPa or higher. A shear modulus of 100 MPa results in a porosity value of 0.68 that is almost independent of the friction angle. This is close to the porosity value found in experiments, namely 0.69, which is for dry SAP particles without using a confining stress. The insensitivity of the porosity value on friction angle is caused by the extremely stiff particles, compared to a confining pressure of 2 kPa, as well as the presence of particle shape. Thus, these packings have already the loosest arrangement resulting in a maximum value of porosity. For shear moduli 1 MPa and 22.1 kPa, the porosity value increases with increasing friction angle.

Relationship of porosity, absorption ratio, and confining stress
To identify the dependency of porosity on both the absorption ratio and confining stress, we have conducted a series of simulations with varying values of Q abs [1-40 g/g] and σ[10-4000 Pa], using a friction angle of 5.5°and Eq. (6) to describe the shear modulus in combination with a value for β of 22.1 kPa. Fig. 9 shows that the porosity decreases as function of absorption ratio, but this decrease becomes more pronounced with increasing confining stress. In the work by Diersch et al. [6], the following fitting equation was employed to describe the porosity (ϕ) as function of absorption ratio for a constant confining stress: in which ϕ scale and ϕ exp are fitting parameters and ϕ max is the maximum porosity value, which is that of dry SAP particles. Our model results are fitted by Eq. (10) using the following fitting equations: The fits are shown in Fig. 9a as solid lines. The quality of fitting was described by a R 2 value of 0.98. However, as discussed in Section 3.4, SAP have a different modulus when in a dry-state than predicted by Eq. (6). Let us assume that the transition from dry-state to slightly swollen occurs between Q abs of 1 g/g and 5 g/g. If we replace the initial porosity in Fig. 9a by that of dry SAP particles, namely 0.68, we can obtain fitting parameters for more realistic porosity and absorption ratio curves (see Fig. 9b): Using Eq. (12) to fit Eq. (10) to simulated data resulted in a reasonable fit although the fitted curves deviate significantly from their corresponding data. Fitting Eqs. (11) and (12) can be employed for continuum-scale simulations where the confining stress can vary during simulations. Note, that the parameterization of Eqs. (11) and (12) may change for various types of SAP particles.

Conclusions
In this research, the effect of particle shape is investigated on the porosity of packings of Super Absorbent Polymer (SAP). The discrete element method (DEM) is employed as grain-scale model in which particle shape is represented by sets of overlapping spheres, so called clumps or multi sphere approximations. Clumps were constructed based on 3-dimensional reconstructions of real SAP particles that were obtain from micro-CT imaging. Using clumps, particle packings were confined using 0.3Psi (2068 Pa) confining pressure. A variety of packings were generated for varying degree of swelling, friction angle, and confining stresses. The minimum number of particles required for porosity values to be independent on the number of particles was found to be 3000 particles. The minimum number of shapes (i.e. clump realizations) was found to be 10. This indicates that the variation in particle shapes is a statistical problem, which allows to simplify simulations because only a small number of particles shapes are needed. The porosity versus absorption ratio plots were generated for SAP particles for which we have experimental data. We found a good agreement with experiments if a relatively low friction angle was used (b5.5°). A 3-dimensional surface of porosity, absorption ratio, and confining pressure was generated that can be employed as a closure equation in continuum-scale simulations.
Appendix A. Determination of φ, α, and ϑ In this Appendix, the values of parameters φ, α, and ϑ are determined (see Section 2.2 for their definition). Clump should represent SAP particles as good as possible, while not using too many spheres. To identify appropriate values of φ, α, and ϑ for SAP particles, we have chosen one SAP particle for which different clump realizations were generated by varying: i) φ values from 0.025 to 0.075, which defines the fraction of spheres having their centres outside the particle; ii) α values from 0.05 to 0.25, which defines the bulginess of the spheres, and iii) ϑ from 0.75 to 0.95, which defines the fraction of the surface of a particle that is overlaid by spheres. After all realizations were constructed, the quality of the clump was evaluated in terms of sphericity. Sphericity is computed for both the SAP particles and the clump realizations; it is defined as: in which V is the volume and A is the surface area of the particle or clump. Fig. A.1a shows the sphericity of all clump realizations of one SAP particle as function of the number of spheres. The sphericity decreases with increasing number of spheres per clump, which is related to the increase in surface area when using more spheres in a clump. The sphericity tends to go to unity when the number of spheres goes to unity, which is indeed the sphericity of a sphere. A broad range of clump realizations give a reasonable approximation of the sphericity of the SAP particle. In particular, using more than 9 particles per clump yields good results. The closest sphericity value was obtained using φ = 0.025, α = 0.05, and ϑ = 0.85, which gave a sphericity of 0.70 which is close to that of the SAP particle, which 0.72.
By matching the sphericity of a clump to that of a SAP particle, the ratio of volume to surface area is maintained. However, the magnitude of both the surface area and the volume of a clump can deviate from that of a particle. To investigate this phenomenon, we plot the normalized volume and normalized surface area in Fig. 1b, in respect to their corresponding values of the SAP particles. Thus, in that Figure, unity represents a perfect match (see dashed lines). Results show that the normalized volume increases roughly linearly with the normalized surface area ratio. The best realization in terms of sphericity (deviation of sphericity within 5% deviation) yields a deviation in normalized volume and normalized surface area of 1.17 and 1.15, respectively. Even though, better approximations of both normalized values are found using different parameters (φ = 0.025, α = 0.05, and ϑ = 0.75), we use the value of φ, α, and ϑ that correspond to the best sphericity value. This is to preserve the volume-tosurface area ratio of a particle, which is deemed more important than to mimic the absolute values of surface area and volume. Appendix B. Clump library Fig. B.1a shows the statistics of clump realizations for all SAP particles, using φ = 0.025, α = 0.05, and ϑ = 0.85. Clump realizations overestimate both the volume and surface area of the SAP particles for round-like SAP particles, by a factor up to 1.6. For plate-like particles, the clump realization overestimates the volume by a factor up to 2.8. This is an artefact of using spheres to represent a particle, because the likelihood of generating a sphere that sticks out of a thin particle is large and on top of that the fact that a sphere sticks out has a more pronounced effect for thin particles than for round-like particles. Subsequently, the volume gets overestimated. As a consequence, the sphericity of clumps that represent plate-like particles is larger than the particle, see Fig. B.1b.