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

. Wheat grains can be considered as a natural cemented granular material. They are milled under high forces to produce food products such as ﬂour. The major part of the grain is the so-called starchy endosperm. It contains sti ﬀ 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 e ﬀ ect 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 e ﬀ ects of starch-protein interface adhesion and voids, the granule size distribution has a limited e ﬀ ect on the samples’ yield stress.


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][5][6][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 e-mail: kmheinze@um2.frcan 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.

Sample construction
Dense and homogeneous cemented granular samples were created, consisting of the three bulk phases: starch ρ s in the form of discs, a matrix ρ p , representing the embedding protein network and a void phase ρ v , if matrix porosity is implemented.The condition ρ s +ρ p +ρ v = 1 is fulfilled for all samples.For most samples ρ 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•10 −3 (standard deviation σ A = 0.25r A ) for the A-type discs and a normal distribution around radius r B ≤ r A (σ B = 0.25r B ) for the B-type discs, which was varied.Due to the finite mesh size, a cutoff of A,B = r A,B ± 1.5σ 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 * (πr A 2 ) and length L = √ 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: ng A,B being the respective total number of discs per type.
To define the constant A we assumed a coverage with Btype 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 From equation ( 1) and ( 2) it follows that For R = r B /r A = 0.2 it follows that A = 1.59.This value was kept constant for all samples, while r B was varied.To guarantee homogeneous packing along with a constant 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 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 ≈ 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 μ = 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 (ρ 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 R are shown in Fig. 3.
For R = 0.2 a sub-series of tests was carried out on microstructures containing void spaces (ρ 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%.

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 ρ discX s to node in ρ discX s , all nodes within ρ p to each other), (II) connecting two nodes within the same material phase in different elements (node in ρ discX s to node in ρ discY s ) or (III) connecting nodes of different material phases (node in ρ s to node in ρ 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 quasistatic uniaxial tension in the y-direction, with the lower boundary being fixed, until integrity failure occurred.

Results and Discussion
The increase of R lead to a less dense packing of discs and therefore more free spaces in between the discs, which could be occupied by the matrix material.For the selected variations of R a non-negligible, non-linear increase of the protein matrix volume fraction ρ p from 13% (R min = 0.2) to 19% (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 ρ Y is given as normalized by the yield stress of pure matrix (ρ [m] ): ρ Y =ρ y /ρ [m] ; with ρ y being the absolute yield stress.An increase of R lead to a decrease of the yield stress (Fig. 4).The lowest yield stress was obtained for R = 0.9, which was decreased by about 10% compared to the maximum yield stress obtained for 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 R merges for Series 1 and 2 (see Fig. 4b).The effect of R on yield stress is therefore independent of the specific interface properties.As mentioned before, changes in R are inevitably accompanied by changes of the volume fraction of ρ p , if ρ 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, 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 R=0.2 (ρ v =3.5%, ρ 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 R, but fixed ρ 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 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 ρ Y (Data not shown).

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.

Figure 1 .
Figure 1.SEM image of a sectioned wheat endosperm, showing the dense packing of starch granules.Scale bar: 30μm

Figure 2 .
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 A = 1.59.

Figure 3 .Figure 4
Figure 3. Five examples of square, bi-modal granular samples with increasing 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.

Table 1 .
The model parameters elasticity E in [Pa] and toughness K in [Jm −3 10 −4 ] are summarized for both series.