# Numerical modeling of the tensile strength of a biological granular aggregate: Effect of the particle size distribution

Karsta Heinze, Xavier Frank, Valerie Lullien-Pellerin, Matthieu George, Farhang Radjai, Jean-Yves Delenne

► To cite this version:

Karsta Heinze, Xavier Frank, Valerie Lullien-Pellerin, Matthieu George, Farhang Radjai, et al.. Nu-merical modeling of the tensile strength of a biological granular aggregate: Effect of the particle size distribution. International workshop on Powders and Grains, Jul 2017, Montpellier, France. 10.1051/epjconf/201714008013 . hal-01539472

HAL Id: hal-01539472

https://hal.science/hal-01539472

Submitted on 2 Jun 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical modeling of the tensile strength of a biological granular aggregate: Effect of the particle size distribution

Karsta Heinze<fnref n="1" />,<fnref n="2" />,*, Xavier Frank<fnref n="2" />, Valérie Lullien-Pellerin<fnref n="2" />, Matthieu George<fnref n="1" />, Farhang Radjai<fnref n="3" />,<fnref n="4" />, and Jean-Yves Delenne<fnref n="2" />

**Abstract.** Wheat grains can be considered as a natural cemented granular material. They are milled under high forces to produce food products such as flour. The major part of the grain is the so-called starchy endosperm. It contains stiff starch granules, which show a multi-modal size distribution, and a softer protein matrix that surrounds the granules. Experimental milling studies and numerical simulations are going hand in hand to better understand the fragmentation behavior of this biological material and to improve milling performance. We present a numerical study of the effect of granule size distribution on the strength of such a cemented granular material. Samples of bi-modal starch granule size distribution were created and submitted to uniaxial tension, using a peridynamics method. We show that, when compared to the effects of starch-protein interface adhesion and voids, the granule size distribution has a limited effect on the samples' yield stress.

## 1 Introduction

Wheat is one of the three most important cereal crops worldwide. The major part of the grain is the starchy endosperm, which consists of a multi-modal distribution of stiff starch granules [1, 2], that are embedded in a softer protein matrix (see Fig.1). Therefore, the endosperm can be considered as a cemented granular material.

The milling of the grains is a crucial step in the processing of wheat. It separates the different grain parts and results in a transformation of the grains into flour. The milling performance is determined by endosperm texture, which depends on genetics and environment. The genetic component influences grain hardness, which was suggested to result from the starch-protein adhesion strength. The environmental conditions impact the porosity of the protein matrix, commonly measured as endosperm vitreousness [3]. Both components have been studied by experimental and numerical approaches [4–7], while other possible sources for differences in milling behavior, such as the starch granule size distribution (SGSD), are still less well understood. Experimental studies on the connection of SGSD and endosperm texture are not numerous and even not convergent [8, 9]. Studies on the milling performance as a function of SGSD are contradicting, too [2, 10]. However, one point of general agreement is that SGSD is genetically controlled [8]. Therefore, the study of SGSD using numerical modeling appears essential to understand its specific effect on the fragmentation behavior and they also

can help the selection of wheat varieties for breeding.

The fragmentation behavior of wheat endosperm has already been the focus of DEM or LEM studies [6, 7, 11], but with a focus on starch-protein adhesion and the volume fractions of the different material phases. We use a peridynamic approach [12, 13] to study how the variation of the starch granules polydispersity affects the yield stress and failure properties. Here we present the methodology in detail along with some first results.

## 2 Methods

### 2.1 Sample construction

Dense and homogeneous cemented granular samples were created, consisting of the three bulk phases: starch $\rho_s$ in the form of discs, a matrix $\rho_p$, representing the embedding protein network and a void phase $\rho_v$, if matrix porosity is implemented. The condition $\rho_s+\rho_p+\rho_v = 1$ is fulfilled for all samples. For most samples $\rho_v$ was chosen to be 0 and therefore all non-starch spaces were filled with the protein phase.

Square volumes were filled with discs of a bi-modal distribution of larger A-type discs and smaller B-type discs, embedded in a matrix phase. Two distributions of disc radii were implied: a normal distribution around a constant radius $r_A = 5 \cdot 10^{-3}$ (standard deviation $\sigma_A = 0.25r_A$) for the A-type discs and a normal distribution around radius $r_B \le r_A$ ($\sigma_B = 0.25r_B$) for the B-type discs, which was varied. Due to the finite mesh size, a cutoff of

<footnote n="1">¹Laboratoire Charles Coulomb, UMR 5221, CNRS-Université de Montpellier, 163 Rue Auguste Broussonet, Montpellier, France</footnote>

<footnote n="2">²UMR IATE, CIRAD, INRA, Montpellier SupAgro, Université de Montpellier, Montpellier, France</footnote>

<footnote n="3">³LMGC, UMR 5508, CNRS-Université de Montpellier, 163 Rue Auguste Broussonet, Montpellier, France</footnote>

<footnote n="4">⁴MSE2, UMI 3466, CNRS-MIT, DCEE, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, USA</footnote>

<footnote>*e-mail: kmheinze@um2.fr</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** SEM image of a sectioned wheat endosperm, showing the dense packing of starch granules. Scale bar: 30µm</figcaption></figure>

$c_{A,B}^{+,-} = r_{A,B} \pm 1.5\sigma_{A,B}$ had to be applied to avoid the occurrence of very small discs, which could not be meshed sufficiently. Consequently, the smallest allowed grain would be meshed by about 10 nodes along the diameter.

In the first step, a square box of the area $A_{box}=900*(\pi r_A^2)$ and length $L = \sqrt{A_{box}}$ was filled with A-type discs using a Fast Poisson Disc Sampling (FPDS) method [14]. The first disc with a random radius within the normal distribution around $r_A$ was placed at a random position within the box. Then a high number of preliminary discs with random radii within the distribution around $r_A$ were randomly placed in the periphery of the already existing disc.

All preliminary discs were tested consecutively for overlap with the existing discs within a limited distance. If no such overlap exists, the preliminary disc was permanently placed into the sample volume. Additionally, a minimum distance $d_{min}$ between the discs was required. The steps of randomly selecting a disc from the sample box, creating preliminary discs and checking for overlap were repeated. If after a finite number of iterations, no more discs could be permanently placed, the box was regarded as being filled. In the second step, smaller discs with a random radius following the distribution around radius $r_B$ were placed into the free spaces remaining in-between the bigger discs. The same FPDS approach as for the A-type discs was used, only that no minimum distance between discs was required and which therefore enabled a dense packing of discs.

Samples with fixed radius $r_A$ and different radii $r_B$ were created to investigate the effect of varying B-type disc radius on the mechanical properties. The radius $r_B$ was varied from $r_B = 0.2r_A$ to $r_B = r_A$ in steps of $0.1r_A$. The resulting nine different radius ratios $R = r_B/r_A$ were therefore [0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]. For every radius ratio $R$ five independent samples were created.

An increase in $r_B$ inevitably resulted in a different disc packing, including a significant change in the number of A- and B-type discs and the total number of discs. Therefore, the sole parameter that was chosen to be kept constant between samples was the ratio of total area occupied by A-type discs to total area occupied by B-type discs:

$$ \mathcal{A} = \sum_{ng_A} (\pi r_A^2) / \sum_{ng_B} (\pi r_B^2) \quad (1) $$

<figure><img src="image_3.png" /><figcaption>**Figure 2.** The schematic of an A-type disc of radius $r_A$ surrounded by B-type discs of radius $r_B$ illustrates the calculation of the area ratio to be $\mathcal{A} = 1.59$.</figcaption></figure>

$ng_{A,B}$ being the respective total number of discs per type. To define the constant $\mathcal{A}$ we assumed a coverage with B-type discs of 83% of a perimeter with the radius $r_{equ} = r_A + r_B$ around one A-type disc for $R = 0.2$ in two dimensions. An illustrating scheme is shown in Fig.2. This coverage allowed some "sharing" of B-type discs between different A-type discs. The number of B-type granules $n_B$ that cover 83% of the perimeter was estimated by

$$ n_B = 0.83 \frac{2\pi r_{equ}}{2r_B} = 0.83\pi \left( 1 + \frac{r_A}{r_B} \right) \quad (2) $$

From equation (1) and (2) it follows that

$$ \mathcal{A} = \frac{\pi r_A^2}{n_B \pi r_B^2} = \frac{1}{0.83\pi (r_A r_B + r_B^2)} \quad (3) $$

For $R = r_B/r_A = 0.2$ it follows that $\mathcal{A} = 1.59$. This value was kept constant for all samples, while $r_B$ was varied. To guarantee homogeneous packing along with a constant $\mathcal{A}$ a sufficient number of B-type discs was needed to fit in the spaces between the A-type discs. Therefore, the minimum distance $d_{min}$ between A-type discs was increased along with an increasing $r_B$. Iterations of a step-wise increase of $d_{min}$ were implied in the disc packing computations, until the required area ratio $\mathcal{A} = 1.59$ was met for every $R$.

With a limited computation power and time, FPDS was able to achieve dense, but not the densest possible packing of discs. Therefore, a classical Discrete Element Methods (DEM) approach [11] was used to apply a slight compression over 100.000 time steps of $10^{-4}$ s. The contact law of a classical linear spring-dash pot model was used with the relation of normal stiffness of the granules $k_n$ to the confining pressure $P$ being $k_n/P \approx 10^4$. As a consequence the overlap of the rigid granules during the compression was smaller than the size of the applied mesh [15]. Coulombs law of friction was applied with the sliding friction coefficient $\mu = 0.3$ to avoid segmentation of the small granules towards the lower sample border. The confining pressure was set to achieve maximum dense packing and the number of time steps was chosen sufficiently high to result in a mechanical equilibrium. Finally, all free spaces that were not occupied by discs of type A or B, were considered to be filled with matrix ($\rho_v = 0$). A square box of the length
$l_{x,y} = 24 r_A = 0.12$ was cropped out of the center of the sample box to avoid boundary effects. Examples for five different samples with increasing $\mathcal{R}$ are shown in Fig.3.

For $\mathcal{R} = 0.2$ a sub-series of tests was carried out on microstructures containing void spaces ($\rho_v > 0$). The matrix filling was computed following the procedure described by Topin et al. [6] and resulted in samples with an average void volume of 3.5%.

## 2.2 Peridynamic mechanical simulations

The peridynamic approach is an alternative to continuum mechanics. It is based on integral equations and is more accurate than classical LEM [6, 7]. Peridynamics allows taking into account heterogeneous distributions of mechanical properties using bond-based schemes with brittle linear springs. The stress and strain tensors are computed on a rectilinear grid by integrating long-range interactions over a limited neighborhood. As the algorithm is based on non-local interactions, failure occurs in the form of damaged zones with a characteristic thickness.

The square sample volumes were meshed onto a rectilinear grid of $1024^2$ nodes. Starch discs were considered as distinct elements, whereas the protein matrix was considered as one single element. Nodes were connected via brittle linear springs within a certain horizon, which was set to be three nodes in each direction. The bonds were classified to three cases: (I) connecting two nodes within the same material phase in the same element (node in $\rho_s^{discX}$ to node in $\rho_s^{discX}$, all nodes within $\rho_p$ to each other), (II) connecting two nodes within the same material phase in different elements (node in $\rho_s^{discX}$ to node in $\rho_s^{discY}$) or (III) connecting nodes of different material phases (node in $\rho_s$ to node in $\rho_p$). Depending on these cases, the bonds were characterized by an elastic modulus $E_{s,p,i}$ and toughness $K_{s,p,i}$. For simplification in this first study, granules of type A and B were assumed to have the same elastic and toughness properties. All bonds of case (II) and (III) were characterized by interface properties. Therefore a small interface layer was considered to be present around all starch discs, even if two discs are in contact with each other.

Two distinct series with the same sample micro-structures but different interface properties were investigated. The parameter values for elasticity and toughness of the different phases and series are summarized in Tab.1. Series 1 represents the case of the strongest possible adhesion between starch and protein, where protein matrix and interface have the same toughness ($K_p = K_i$). In Series 2 the toughness of the interface $K_i$ was decreased by 50%. The parameters were chosen similar as was reported in [7], taking into account nano-mechanical measurements of starch and protein inside wheat endosperm by atomic force microscopy[16]. Samples were submitted to quasi-static uniaxial tension in the y-direction, with the lower boundary being fixed, until integrity failure occurred.

## 3 Results and Discussion

The increase of $\mathcal{R}$ lead to a less dense packing of discs and therefore more free spaces in between the discs, which

<table><caption>**Table 1.** The model parameters elasticity $E$ in [Pa] and toughness $K$ in [$Jm^{-3}10^{-4}$] are summarized for both series.</caption><thead><tr><th>Series</th><th>Protein</th><th>Starch</th><th>Interface</th></tr></thead><tbody><tr><td>S1</td><td>$E_p=1$<br>$K_p=1 \cdot 10^{-3}$</td><td>$E_s=4E_p$<br>$K_s=3.5K_p$</td><td>$E_i=E_p$<br>$\mathbf{K_i}=\mathbf{K_p}$<br>$E_i=E_p$<br>$\mathbf{K_i}=\mathbf{0.5K_p}$</td></tr><tr><td>S2</td><td></td><td></td><td></td></tr></tbody></table>

could be occupied by the matrix material. For the selected variations of $\mathcal{R}$ a non-negligible, non-linear increase of the protein matrix volume fraction $\rho_p$ from 13% ($\mathcal{R}^{min} = 0.2$) to 19% ($\mathcal{R}^{max} = 1.0$) was observed, the starch volume fraction decreasing accordingly. In a cemented granular material, such as the endosperm, it is not only the SGSD that determines the amount of protein fraction, but also other factors such as packing density and presence of voids. Experimental studies on wheat are inconclusive whether or not SGSD is correlated with grain protein content [8, 17]. Though our first numerical samples could mimic the endosperm of vitreous wheat grains, which have low porosity, a direct comparison with the experimental results is unfeasible at this point. Future numerical simulations, which could investigate the relation of SGSD and protein content further, would need to take into account more parameters, such as void spaces and packing density.

All tested sample volumes showed a brittle behavior. Yield stress $\rho^Y$ is given as normalized by the yield stress of pure matrix ($\rho^{[m]}$): $\rho^Y = \rho_y/\rho^{[m]}$; with $\rho_y$ being the absolute yield stress. An increase of $\mathcal{R}$ lead to a decrease of the yield stress (Fig.4). The lowest yield stress was obtained for $\mathcal{R} = 0.9$, which was decreased by about 10% compared to the maximum yield stress obtained for $\mathcal{R} = 0.2$, for both interface toughness values respectively.

When comparing samples with the same micro-structures, but with a differing interface toughness $K_i$, a more drastic decrease of the yield stress by 50% is observed. In the range of the observed variation, which cover the range of biological micro-structures, the effect of the SGSD is of secondary order compared to the effect of the interface. The high significance of the interface properties has already been established by numerical simulations[6, 7] and our results are consistent with this. When yield stress is additionally normalized by the interface toughness $K_i$, the evolution of yield stress with $\mathcal{R}$ merges for Series 1 and 2 (see Fig.4b). The effect of $\mathcal{R}$ on yield stress is therefore independent of the specific interface properties.

As mentioned before, changes in $\mathcal{R}$ are inevitably accompanied by changes of the volume fraction of $\rho_p$, if $\rho_v=0$. Yield stress was found to decrease with increasing protein volume (data not shown). Because the increase in matrix volume fraction is the consequence of significant changes in the micro-structures, $\mathcal{R}$ and matrix volume fraction have a joint effect on the yield stress, which cannot easily be separated for the presented sample selection.

The introduction of voids into the samples with $\mathcal{R}=0.2$ ($\rho_v=3.5\%$, $\rho_p=8.6\%$, data not shown) had a significant effect on the yield stress in the same range as was observed for the 50% decrease of $K_i$. Samples with varying
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Five examples of square, bi-modal granular samples with increasing $\mathcal{R}$ are shown after integrity failure due to uniaxial tension. Starch discs are shown in blue. All voids are filled with protein matrix (yellow). The interface around the starch granules is highlighted in red. The fracture pathway, signaled by broken bonds, is shown in light blue.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** *a:* The evolution of yield stress $\rho_Y$ with increasing $\mathcal{R}$ is shown for two series with differing interface toughness: $K_i = 1.0$ (circles, S1) and $K_i = 0.5$ (triangles, S2). Yield stress was normalized by the yield stress of pure matrix. *b:* Both series merge, when yield stress is normalized by $K_i$.</figcaption></figure>

$\mathcal{R}$, but fixed $\rho_p$ are difficult to construct and it can be hypothesized that the necessary introduction of voids would strongly dominate over the effect of a changing $\mathcal{R}$.
It is visible in Fig. 3 that the fracture does not propagate through starch granules. This is a direct result of the value chosen for $K_s$, which affects the resulting amount of damage to the starch granules, but not $\rho_Y$ (Data not shown).

## 4 Conclusion

We found that, when compared to the effects of interface adhesion and void inclusions, the SGSD presents a limited effect on the yield stress in tension. The material is strongest for the smallest size ratio (small/large) between the two granule populations. The effect of SGSD is coupled to changes in the volume fractions of the starch and protein phase, but is fully independent of the interface properties. An ongoing study investigates the effect of interface adhesion and void occurrence on the tensile properties for different volume ratios of A- and B-type granules. Additionally, further studies will investigate fracture under dynamical conditions. This is possible with the presented peridynamic assay, which enables simulations under conditions closer to the milling process than quasi-static tests.

## Acknowledgment

K. Heinze received a French grant from INRA (CEPIA dept.) and Labex Numev as a PhD student. We thank the bioinformatics platform Toulouse Midi-Pyrenees (Bioinfo Genotoul) for providing computing resources.

## References

[1] D.B. Bechtel, I. Zayas, L. Kaleikau, Y. Pomeranz, Cereal Chem **67**, 59 (1990)

[2] M. Edwards, B. Osborne, R. Henry, J Cereal Sci **48**, 180 (2008)

[3] F.X. Oury, P. Lasme, C. Michelet, M. Rousset, J. Abecassis, V. Lullien-Pellerin, Theoretical and Applied Genetics **128**, 913 (2015)

[4] V. Greffeuille, J. Abecassis, M. Rousset, F.X. Oury, A. Faye, C.B. L'Helgouac'h, V. Lullien-Pellerin, Theor Appl Genet **114**, 1 (2006)

[5] Y. Haddad, F. Mabille, A. Mermet, J. Abecassis, J. Benet, Powder Technol **105**, 89 (1999)

[6] V. Topin, J.Y. Delenne, F. Radjai, L. Brendel, F. Mabille, Eur Phys J E **23**, 413 (2007)

[7] E. Chichti, V. Lullien-Pellerin, M. George, F. Radjai, R. Affès, J.Y. Delenne, J of Food Eng **190**, 40 (2016)

[8] G. Igrejas, B. Faucher, D. Bertrand, D. Guibert, P. Leroy, G. Branlard, J Cereal Sci **35**, 103 (2002)

[9] D.B. Bechtel, I. Zayas, R. Dempster, J.D. Wilson, Cereal Chem **70**, 238 (1993)

[10] M. Edwards, B. Osborne, R. Henry, J Cereal Sci **52**, 314 (2010)

[11] J. Delenne, Y. Haddad, J. Bénet, J. Abecassis, J Cereal Sci **47**, 438 (2008)

[12] S. Silling, J Mech Phys Solids **48**, 175 (2000)

[13] P. Perré, G. Almeida, M. Ayouz, X. Frank, Ann For Sci **73**, 147 (2016)

[14] R. Bridson, *Fast Poisson Disk Sampling in Arbitrary Dimensions*, in ACM SIGGRAPH 2007 Sketches (ACM, New York, NY, USA, 2007), SIGGRAPH '07

[15] J.N. Roux, G. Combe, C R Phys **3**, 131 (2002)

[16] E. Chichti, M. George, J.Y. Delenne, F. Radjai, V. Lullien-Pellerin, Eur Polym J **49**, 3788 (2013)

[17] S.H. Park, J.D. Wilson, B.W. Seabourn, J Cereal Sci **49**, 98 (2009)