Next Article in Journal
Synthesis, Spectroscopic, and Theoretical Study of Copper and Cobalt Complexes with Dacarbazine
Previous Article in Journal
Zn–0.8Mg–0.2Sr (wt.%) Absorbable Screws—An In-Vivo Biocompatibility and Degradation Pilot Study on a Rabbit Model
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Breakage Strength of Wood Sawdust Pellets: Measurements and Modelling

Institute of Agrophysics, Polish Academy of Sciences, Doświadczalna 4, 20-290 Lublin, Poland
*
Author to whom correspondence should be addressed.
Materials 2021, 14(12), 3273; https://doi.org/10.3390/ma14123273
Submission received: 27 May 2021 / Revised: 10 June 2021 / Accepted: 11 June 2021 / Published: 13 June 2021
(This article belongs to the Section Materials Simulation and Design)

Abstract

:
Wood pellets are an important source of renewable energy. Their mechanical strength is a crucial property. In this study, the tensile strength of pellets made from oak, pine, and birch sawdust with moisture contents of 8% and 20% compacted at 60 and 120 MPa was determined in a diametral compression test. The highest tensile strength was noted for oak and the lowest for birch pellets. For all materials, the tensile strength was the highest for a moisture content of 8% and 120 MPa. All pellets exhibited a ductile breakage mode characterised by a smooth and round stress–deformation relationship without any sudden drops. Discrete element method (DEM) simulations were performed to check for the possibility of numerical reproduction of pelletisation of the sawdust and then of the pellet deformation in the diametral compression test. The pellet breakage process was successfully simulated using the DEM implemented with the bonded particle model. The simulations reproduced the results of laboratory testing well and provided deeper insight into particle–particle bonding mechanisms. Cracks were initiated close to the centre of the pellet and, as the deformation progressed, they further developed in the direction of loading.

1. Introduction

Waste biomaterials, including wood processing residues, are important renewable energy resources [1,2]. Since, in their native loose forms, such materials have low bulk density and their storage and transport are difficult and costly [3], densification is applied to reduce these disadvantages [4]. Biomaterials can be densified into briquettes, pellets, or cubes, among which pelletising is the most popular. The quality of the pellets is related to their physical, chemical, and mechanical properties. Under high-pressure pelletisation, the bulk density of wood sawdust increases from 100–200 to 1000–1200 kg m3, which is the typical density of the pellets [5,6]. The bulk density of pellets is significantly higher than that of natural wood block (400–700 kg m3); however, it is lower than the solid-phase density of sawdust particles (1450 kg m3) [7]. Hydrogen binding at the surfaces of lignocellulosic particles of sawdust provides the main type of binding force used in manufacturing of pellets [1]. With increasing compaction pressure, the contact area between the particles increases and the pellets become denser and more durable. Binding forces are higher when a contact zone between particles covers a larger area; that is, the forces increase with an increase in compaction pressure and a decrease in particle size [8]. To increase binding forces, steam explosion pretreatment is applied; this causes lignin melting and a reduction in particle size [3]. The moisture content (MC) markedly affects the binding forces. It has both antagonistic (with water molecules replacing wood polymer bonds) and protagonistic (decreasing the melting temperature of lignin) actions on pellet durability [2]; therefore, an optimum MC should be maintained during the pelletisation process [9]. Pressure (contact between particles), temperature (plastic deformation of lignin), particle size, and optimum MC are key factors in particle binding [1,3]. Therefore, properly pretreated materials have better physical quality than untreated materials [2,6].
Mechanical durability is a crucial property of pellets [10,11]. However, very simple and convenient standard durability tests provide only a qualitative measure of the mechanical strength of the pellets [12]; therefore, diametral compression can be alternatively applied. This method is also very simple [13] and has been widely used to determine the strength of cylindrical samples by measuring the direct relationship between the load at failure and tensile strength [14]. This test has been used to determine the tensile strength of rocks, concrete, composites, and agglomerates in numerous branches of science and technology, such as civil and chemical engineering [15,16], pharmaceutical science [17,18,19,20], biomass briquettes and pellets engineering [4,6,11,21,22], and food powder engineering [23,24]. Larsson and Samuelsson [11] used diametral compression to show that the tensile strength of pellets was strongly correlated with their density and durability. Using the diametral compression test, Shaw et al. [6] demonstrated that steam treatment leads to an increase of more than tenfold in the tensile strength of poplar wood pellets.
The discrete element method (DEM) proposed by Cundall and Strack [25] has led to new possibilities for modelling interactions and processes between a large number of individual objects. The method has been applied in numerous areas of science and technology, including physics [26], pharmaceutical science [27], civil engineering [28], and agricultural engineering [29]. This method is a very promising tool for modelling biomass post-harvest processing, including handling [30,31,32], compaction [33], and testing biomass bulk properties [34]. The DEM extended with the bonded particle model (BPM) considerably broadened the range of applicability of the DEM to model the internal structure and breakage of agglomerates composed of numerous particles [35,36]. The method has been applied to model the bulk behaviour of deformable biomass fibres [34], packing of flexible fibres [37], compaction of wood chips [38], durability [39], and the strength of pellets [40]. Xia et al. [41] provided a comprehensive discussion of the challenges of modelling the behaviour of the bulk of complex-shaped particles of wood biomass by using the DEM. It seems that DEM modelling of biomass agglomerate breakage is still in an early stage of development, and additional work is needed to demonstrate the validity and limitations of the method.
The aim of this study was to investigate the possibility of numerical reproduction of the pelletisation process of sawdust and the stress–deformation behaviour of pellets in diametral compression using the DEM with a parallel BPM.

2. Experiments

2.1. Materials and Methods

2.1.1. Materials

Sawdust from birch (Betula L.), oak (Quercus L.), and pine (Pinus sylvestris L.) wood obtained from a local sawmill were used for the experiments. The mean particle sizes determined using a set of 0.1 to 2 mm sieves were 0.355 ± 0.028, 0.409 ± 0.024, and 0.375 ± 0.026 mm for birch, oak, and pine sawdust, respectively. The initial MC of all sawdust types was 8% ± 0.7%. To increase the MC to 20 ± 0.9%, sawdust samples were placed in a chamber humidifier at a relative humidity of 70% and a temperature of 21 °C and occasionally shaken for 4 days.

2.1.2. Pellets Preparation

A high-pressure piston-and-mould compaction process was used to densify the sawdust. Samples of sawdust (0.5 g) were compacted in a steel cylindrical mould with a diameter of 10 mm and height of 25 mm up to a pressure of 60 or 120 MPa using a material testing machine (Instron 7782, High Wycombe, UK) with a piston moved at 0.033 mm s−1. To remove the pellet, the base of the mould was removed and the pellet was pressed out [42]. Twelve combinations of pellets were prepared utilising three types of wood, two levels of moisture content, and two levels of compaction pressure. Ten specimens of each variant of pellets were prepared.

2.1.3. Determination of Solid-Phase and Bulk Densities

The solid-phase density of the sawdust particles, ρs, and of the intact wood (40 mm diameter and 70 mm height blocks), ρw, was measured using helium pycnometry (Ultrapyc 1200e, Quantachrome Instruments, Boynton Beach, FL, USA) in five replicates.
The bulk density of the intact wood, BDw, of the pellets prepared at both levels of compaction pressure, BDp, and the initial bulk density of the sawdust (just before densification), BDs, were determined from the mass and volume of the specimens. The volume of the pellets was determined from their diameter and height measured just before the diametral compression test using an electronic calliper with an accuracy of 0.01 mm. The volume of the sawdust was determined from the mould diameter and initial height of the bulk of the sawdust placed in the mould.

2.1.4. Estimation of Porosity and Pore Size Distribution

The total porosity px of the studied materials was estimated from their bulk density BDx as
p x = 1 B D x ρ s
where ρs is the solid density of the sawdust, the index x equals w for the intact wood and p for pellets, and s is the initial bulk of the sawdust. In all cases, the density of the solid phase was assumed to be equal to the density of the sawdust.
Mercury intrusion porosimetry tests were performed to estimate the structural features of the wood and pellets. Prior to the measurements, the specimens were dried overnight at 105 °C. The tests were performed for pressures ranging from ~0.1 to 200 MPa (corresponding to pore radii from ~10.0 µm to 3.8 nm) using an Autopore IV 9500 porosimeter (Micromeritics, Norcross, GA, USA) in two replicates. The intrusion volumes were measured at stepwise increasing pressures, allowing equilibration at each pressure step. The maximum deviations between the replicates were ≤4.5%, and they occurred mainly at low pressures (i.e., for the largest pores). The volume of mercury, V (in m3 kg−1), intruded at a given pressure P (in pascals) gave the pore volume that could be accessed. The intrusion pressure was translated on the equivalent pore radius R (in metres) using the Washburn equation
P = A σ m cos α m R ,
where σm is the mercury surface tension, αm is the mercury–solid contact angle (taken as 141.3° for all studied materials), and A is a shape factor (equal to 2 for the assumed capillary pores).
By knowing the dependence of V on R, a normalised pore size distribution, f(R), was calculated and expressed on a logarithmic scale [43] as:
f ( R ) = 1 V max d V d log ( R ) .

2.1.5. Scanning Electron Microscopy

Scanning electron microscopy (SEM) images of the surfaces of wood blocks, sawdust, and broken pellets were taken in five replicates using a Thermo Scientific Phenom ProX G6 microscope (Thermo Fisher Scientific, Waltham, MA, USA) operating at 10 kV, and the results were analysed using the software provided by the manufacturer. Prior to SEM analysis, all samples were sputtered with gold using a CCU-010 system (Safematic, Zizers, Switzerland).

2.1.6. Diametral Compression Test

The tensile strength of the pellets was determined by diametral compression between two parallel plates using a material testing machine (Lloyd LRX, Advanced Test Equipment Corp., San Diego, CA, USA) in five replicates. The displacement velocity of the plate was 0.033 mm s1. The load was measured with an accuracy of ±0.02 N.
The stress components in the pellet during the diametral compression test are shown in Figure 1. In the centre of the pellets (where the geometrical coordinates are x = 0 and y = 0), the stresses σx and σy are the principal stresses σ1 and σ2, respectively [14]. The tensile strength at failure, σf, is identified by the maximum tensile stress σ1,max in the direction perpendicular to the load in the centre of the object [13]:
σ f = σ 1 , max = P f π R p h p ,
where Pf is the failure load, Rp is the radius of the pellet, and hp is its height.

2.1.7. Statistics

The data were statistically analysed using Statistica v. 10.0 (StatSoft Inc., Kraków, Poland) software and presented as mean values plus or minus the standard deviation (SD). One-way analysis of variance along with Tukey’s test with a significance level of 5% was applied to determine the statistical significance of the differences between the mean values.

2.2. Experimental Results

The densities of wood and sawdust are listed in Table 1. The high solid-phase density of the sawdust was probably the closest to the true solid-phase density of the wood itself. We surmised that the closest packing of the main wood component, cellulose, was similar to that of crystalline sugar (e.g., saccharose) composed of the same chemical subunits. The density of the solid phase of saccharose (1.59 g cm−3 [44]) is close to the densities of all the types of sawdust studied, whereas the density of the wood with its closed pores inside is markedly lower. This may indicate that the vast majority of closed wood pores were damaged during sawdust preparation, with only a small portion eventually remaining. Therefore, the density of the sawdust particles was considered as the true density of the wood solid phase. This was the reason for applying its value for bulk density calculations in Equation (1). The lower solid-phase density of the intact wood indicates that its structure contained closed pores that were not available for gas penetration. The volumetric percentages of these pores, calculated using data from Table 1, were 26.8 ± 0.3 for oak, 2.9 ± 0.2 for pine and 6.5 ± 0.1 for birch.
The bulk density and porosity of the studied intact wood, sawdust, and pellets are presented in Table 2. The porosity of the intact wood ranged from 53.8% (oak) to 65.2% (pine). The initial porosity of the bulk sawdust was ~80% for all the samples. Compaction of the sawdust resulted in a profound reduction of porosity: that of the pelletised samples ranged from 58.4% for pine (MC = 20%, with 60 MPa compaction) to 31.1% for oak (MC = 8%, with 120 MPa compaction). The porosity of the pellets was considerably lower than that of intact wood. The reduction in the porosity of the sawdust during pelletisation can be considered as a combination of two factors: continuous reduction of external pores between particles and closing of the internal pores present inside sawdust particles during the final stages of compaction [38]. At higher MC, the bulk density of the sawdust and the pellets was lower, which was most probably caused by the swelling of wood cells and/or by hydration of cellulose molecules. As a direct consequence of the decrease in the bulk density, the calculated porosity at higher moisture levels was higher. However, the porosity value may be uncertain because the density of the solid phase (sawdust) measured at 8% moisture was used. The water-involving processes described above may also alter the solid-phase density. Other researchers have reported porosities of bulk sawdust, intact wood, and highly compacted pellets of ~85%, ~45–60%, and ~15%–30%, respectively [5,6,7,11], which are consistent with the results of the present study.
The dependencies between the pore volume and pore radius and the pore size distribution functions for the studied wood and pellets are presented in Figure 2 and Figure 3, respectively.
From both figures, one can see that the structures of the different types of wood were quite different. Pine wood had the largest pores, all ~15 µm. Birch wood exhibited a bimodal pore size distribution. A smaller peak was located at ~6 µm and the higher peak was at ~0.8 µm. Oak wood had three groups of pores, developed at 0.08, 0.22, and 6.9 µm. The structures of all pellets appeared to be similar, regardless of the origin of the initial material. The dominant pores in all pellets developed at ~7.9 µm and the second, smaller peak was located at ~1.2 µm. This uniformisation of the pellet structure may have been due to several factors: damage to the wood structure during sawing, the produced sawdust having a similar grain shape regardless of the initial material (as seen in the SEM images of the sawdust of different wood presented below), and the same conditions of the pelletisation process. This similarity of the pellet structure certified the application of the unified model for all pellets together in the further modelling stage.
Figure 4 presents SEM images of the surfaces of the cross-sections of wood blocks, sawdust particles, and broken pellets. Differences in the internal structures of the wood blocks, sawdust particles, and pellets can be easily observed. The surface structure of sawdust particles was more complex compared to that of wood blocks as a result of sawing. Pores between sawdust particles and the majority of their surface pores disappeared during the pelletisation process. On the surfaces of the broken pellets, some bonds between sawdust particles could be identified. Their sizes were markedly smaller than the mean size of the particles. As the size of the bond zone was evaluated to be ~40 μm, the bond radius of 20 μm was used as the baseline input data for calibrations of the DEM modelling parameters.
The compaction pressure–piston displacement relationships (σz versus Δh) for pelletisation are presented in Figure 5a. All dependencies were similar. The noticeable increase in the compaction pressure started at a piston displacement of >10 mm (i.e., >50% of the total displacement). For a higher MC, the maximum compaction pressure was reached at a considerably lower piston displacement (except with oak) and the unloading process lasted longer (Figure 5a inset). This resulted in a lower bulk density of pellets produced from sawdust with a higher MC.
The stress–deformation dependencies for the studied pellets are depicted in Figure 5b. ΔL is the displacement of the loading piston and D is the pellet diameter measured directly before diametral compression. All dependencies of σ1 on ΔL/D were smooth and round without any sudden drops, which is typical for the ductile breakage mode.
The SD of the mean values of the tensile strength σf and the deformation at failure, ΔLf/D, indicate the degree of variability of the experimental relationships resulting from the natural differentiation of pellets (Table 3). Three groups of pellets can be distinguished with respect to the tensile strength: σf > 1 MPa (oak, σz = 120 MPa, and MC = 8%); σf ~ 0.5 MPa (oak, σz = 60 MPa, and MC = 8%; oak, σz = 120 MPa, and MC = 20%; pine, σz = 120 MPa, and MC = 8%); and σf < 0.3 MPa (all others).
An increase in MC and a decrease in σz resulted in a decrease in the tensile strength of the pellets. The highest tensile strength was observed for oak pellets (MC = 8% and σz = 120 MPa) and the lowest for birch pellets (MC = 20% and σz = 60 MPa). An increase in MC and decrease in the compaction pressure resulted in a fourfold decrease in the tensile strength of oak and birch pellets and an eightfold decrease for pine pellets. The range of the determined values of σf is consistent with the findings of other researchers for pellets produced from untreated wood materials [4,6,22]. The range of the deformation at failure, ΔLf/D, found in the present study (0.56–0.82) is very similar to that obtained by Gilvari et al. [40] in the diametral compression tests of biomass pellets (0.05–0.07).

3. Modelling

3.1. DEM Setup

3.1.1. Contact Model

The linear hysteretic spring contact model introduced by Walton and Brown [45] was used for the simulations. The model was extended to account for linear adhesion in the form introduced by Luding [46] to keep particles glued and to minimise changes in the structure of the pellet during unloading and removal from the mould. The contact force–displacement scheme in the normal direction shown in Figure 6 accounts for the plastic contact deformation, elastic unloading and reloading, and attractive adhesion force:
f n = k 1 δ n , k 2 ( δ n δ n , 0 ) k 1 δ n   ( loading ) , k 2 ( δ n δ n , 0 ) , k 1 δ n > k 2 ( δ n δ n , 0 ) > k c δ n   ( loading - unloading ) , k c δ n , k c δ n k 2 ( δ n δ n , 0 )   ( unloading ) ,
where f n is the contact normal force, k1 is the loading (plastic) stiffness, k2 is the unloading (elastic) stiffness, kc is the adhesive stiffness, δn is the overlap in the normal direction, and δn,0 is the residual overlap during unloading when the force f n changes from repulsive to attractive interaction. Adhesive interactions keep the overlap δn among the contacts within the pellet removed from the mould at a level very close to the residual overlap δn,0.
The plastic stiffness k1 is related to the yield strength py of a particle [47],
p y = 2 E π δ n , y r ,
by the following formula [48]:
k 1 = 5 r * min ( p y , i , p y , j ) ,
where E is the Young’s modulus, r is the radius of a particle, δn,y is the yielding overlap, r * = r i r j / ( r i + r j ) is the equivalent radius of the contacting particles, and py,i and py,j are the yield strengths of particles i and j, respectively.
Energy dissipation in the normal direction results from the difference between the loading and unloading stiffnesses. The elastic stiffness k2 for unloading and reloading is related to k1 through the restitution coefficient e as follows [48]:
e = k 1 k 2 .
The particle–particle force in the tangential direction, f t c , was updated incrementally as follows:
f t = ( f t ) 0 + k t Δ δ t + f t d , Δ δ t Δ δ t μ p - p f n , if if f t < μ p p f n , f t μ p - p f n ,
where ( f t ) 0 is the tangential force at the end of the previous timestep, kt and δt are the stiffness and overlap in the tangential direction, respectively, and μp-p is the particle–particle friction coefficient.
The stiffness in the tangential direction, kt, was assumed to be equal to the stiffness in the normal direction. The velocity-dependent dissipative component f t d of the tangential force ft is defined as [48,49]:
f t d = 4 m * k t 1 + π ln e 2 v t ,
where m * = m i m j / ( m i + m j ) is the equivalent mass of the contacting particles, vt is the relative velocity in the tangential direction, and e is the same restitution coefficient used for the energy dissipation in the normal direction.
The dissipative torque Mi associated with rolling friction mr was introduced as
M i = m r f n r i ω i ,
where ωi is the unit angular velocity vector of particle i at the contact point.
The binding mechanism inside pellets is mainly based on the re-solidification of melted lignin during cooling [2]. Therefore, to bring the DEM modelling closer to the experimental conditions, the adhesive interactions in the contacts of the relaxed pellet were replaced by the parallel BPM, as proposed by Potyondy and Cundall [35]. In contrast to the lack of force interactions in the tangent direction assumed in the adhesion model, the BPM provides stiffness in the normal and tangential directions. Additionally, substituting adhesive interactions with the BPM protects the simulations against the appearance of secondary adhesive contacts during deformation, which might falsify the resultant breakage strength of the pellets. The BPM was initiated when the total kinetic energy of the assembly decreased below 10−8 J and the overlap was very close to the residual overlap, δ~δn,0 (Figure 7).
The forces and moments of the bond connections between the two particles were calculated incrementally [48] as follows:
Δ f n b = v n k n b A Δ t , Δ f t b = v t k t b A Δ t , Δ M n b = ω n k t b J Δ t , Δ M t b = ω t k n b J 2 Δ t ,
where vn and vt are the relative velocities in the normal and tangential directions, respectively; k n b and k t b are the stiffnesses in the normal and tangential directions, respectively; A = π r b 2 and J = π r b 4 2 are the area and moment of inertia of the bond cross-section, respectively; rb is the radius of the bond; and Δt is the time increment.
The Young’s modulus of the bond is:
E b = k n b ( r i + r j ) .
The bond is broken when the maximum normal stress σ n , max b exceeds the tension strength σc or the maximum tangent stress τ max b exceeds the shear strength τc:
σ n , max b = f n b A + 2 M t b J r b > σ c , τ max b = f t b A + M n b J r b > τ c .
The critical bond force associated with the tension strength σc is
f c = A σ c ,
as is depicted in Figure 7.

3.1.2. DEM Modelling Parameters

Spherical particles simulating sawdust used in the experiments, with radii normally distributed in the range of 0.14–0.26 mm with the mean = 0.2 mm and SD = 0.02 mm, were generated randomly in a mould of the same size as the mould used in the experiments. The assembly consisted of 40,000 particles. The spherical shape of the particles was used to simplify and fasten the numerical calculations. The DEM simulation parameters are listed in Table 4.
The porosity of random bedding of spherical particles in the DEM simulations was ~37%. The high porosity of the bulk sawdust, ~80%, resulted from the superposition of two components: free spaces between particles and micro- and macro-pores on the surface of sawdust particles and (eventually) inside. Therefore, the porosity of the bedding of spherical particles could represent only the pores located between the sawdust particles. The internal pores inside the particles and pores resulting from the surface roughness of the sawdust particles were not covered by this approach. To model the compaction of deformable wood chips in a cycling loading test using the DEM, Xia et al. [38] set the bulk density of the intact wood as the solid density of particles and reproduced the compaction process well. Similarly, we wanted to set the bulk density of the intact wood as the solid density of the particles. However, we decided to set a particle solid density of 680 kg m−3 as the uniform value for the solid density of all types of sawdust particles. This value was slightly higher than the bulk density of intact oak wood and higher than those of pine and birch. We decided to set the same particle densities for each type of sawdust because we surmised that the structure of the sawdust of all types of wood would be uniform during sawing. Almost identical values of solid-phase density of each type of sawdust (see Table 1), along with almost identical pore size distributions of all pellets (see Figure 4), confirmed this assumption. The assumed high sawdust particle density was a direct consequence of the high solid-phase density. By setting the bulk density of the intact wood as the solid density of particles, in the first stage of compaction, the free spaces between particles were reduced to zero, and, with further compaction, an increase in the density of the pellet above the density of the intact wood reproduced the closing of internal micro- and macro-pores inside sawdust particles.
We set the coefficients of particle–particle friction, μp-p, and particle–wall friction, μp-w, to 0.5 and 0.15, respectively, in accordance with similar studies on biomass pellets [40] and pinewood chips [38]. The effect of the rolling friction coefficient mr on the behaviour of pellets was not considered because, in reality, the rolling of elongated dust particles seems to be absent. Therefore, a very low default value of mr = 0.01 was used in the simulations, as recommended by the EDEM software package [48]. A Poisson ratio of ν = 0.35 and a coefficient of restitution of e = 0.5 were adopted for the DEM simulations, as typical values for wood materials [7,30].
The modulus of elasticity of oak, E = 15.7 GPa [7], was set as the representative value of the modulus of elasticity of sawdust particles. The yield strength py and the resulting plastic stiffness k1 were determined as the values that best fit the σzh/h0) relationship during the compaction process. The bond radius rb evaluated from the SEM images of the cross-sections of the pellets was used as baseline input data for the DEM simulations. The elastic properties of the bonds were treated as adjustable parameters for the calibration process. The bond elasticity modulus Eb, the strength σc, and the final value of the bond radius rb were determined as the values providing the best fit of the σ1L/D) relationship during diametral compression.
To reduce the computational time, the particle density was increased by a factor of 104. To keep the gravitational force unchanged, the gravitational acceleration was reduced by a factor of 104. As also shown in [29], scaling the density by a factor of 106 did not change the shape of the σ1L/D) characteristics and reduced the tensile strength σf by only 1.2% of the strength of the sample without density scaling. Therefore, scaling density by a factor of 104 was considered to introduce an error of <1.2%. Time integration was performed with steps of 2 × 10−6 s, that is, 4% of the Rayleigh timestep [50].

3.1.3. Stages of DEM Simulations

As illustrated in Figure 8, DEM simulations were performed in the following sequence: particle generation and filling of the mould (Figure 8a), compaction (Figure 8b), unloading and mould removal, relaxation (Figure 8c), BPM initiation, and then diametral compression (Figure 8d). Simulated sawdust particles were randomly generated in the space of the mould. After settlement, the particles were compacted with a virtual piston with an axial velocity of 3.3 × 10−5 m s−1 up to the assumed value of the compaction pressure σz of 60 or 120 MPa. The unloading process was performed in two stages: (1) unloading in the z direction by piston removal with the same velocity as during loading and (2) mould removal by increasing the mould radius. Next, after the assembly was relaxed, the BPM was initiated and the adhesive interactions were gradually reduced to zero. The last stage was the diametral compression between parallel plates with the same displacement rate as in the previous stages. The colours in Figure 8 represent the average compressive force exerted on the particles. The force scale ranges illustrate the compressive force differences in the simulated objects: the maximum compressive force in the relaxed pellet was five orders of magnitude lower than that in the compacted pellet (compare the scales “c” and “b” in Figure 8) and the maximum compressive force during diametral compression at σ1 = 0.7 MPa (~0.5σf) was one order of magnitude lower than that in the compacted pellet (scales “d” and “b” in Figure 8).

3.1.4. Fit Quality

The quality of fit of the experimental stress–deformation relationship σ1L/D) obtained from the simulation was evaluated using the relative root mean-square error (RRMSE):
R R M S E = 1 n i = 1 n σ 1 , exp i σ 1 , D E M i σ 1 , exp i 2 ,
where σ 1 , exp i is the tension stress in the centre of a pellet (see Figure 1) obtained from the experiment, σ 1 , D E M i is the tension stress obtained from the DEM simulation, and n is the number of samples in the σ1L/D) relationship.

3.1.5. Calibration of Material Parameters for DEM Modelling

The impact of five material parameters of the DEM simulations was considered: the yield strength of particles py, the coordination number CN, the bond radius rb, the bond elasticity modulus Eb, and the bond tension and shear strength σc = τc. The values of py and CN were calibrated against the experimental data of sawdust compaction and pellet relaxation. The bond radius was evaluated based on SEM images of broken pellets and a preliminary modelling of the diametral compression test fitting to the range of experimental values of tensile strength and deformation at failure. The elastic modulus and strength of the bonds were determined as a set of two independent parameters that provided the best fit to the experimental diametral compression stress–deformation relationship.

3.2. DEM Run

3.2.1. Compaction of Sawdust: Plastic and Elastic Stiffness of Particles

The yield strength of sawdust particles py was determined as the value providing the minimum RRMSE of the DEM approximation of the experimental stress–deformation relationship of compaction process. As presented in Figure 9, the experimental relationships between σz and Δh/h0, where h0 is the initial height of the sawdust in the mould, for compaction of sawdust were very similar for the three wood materials. Therefore, a common value of py was applied for all three materials. The DEM simulation performed for py = 100 MPa fitted all the experimental relationships very well for Δh/h0 < 0.75 (RRMSE < 0.15). The ratio of the yielding overlap to the particle radius, δn,y/r, was 9.63 × 10−3. The plastic stiffness of particles, k1, corresponding to py = 100 MPa, was 1 × 105 N m−1, and the elastic stiffness, k2, was 4 × 105 N m−1; consequently, the assumed value of the restitution coefficient was 0.5. The lower quality of approximation of the experimental relationship for the highest level of compaction (Δh/h0 > 0.75) resulted from the difference in curvature between the experimental and simulated lines. It is possible that a nonlinear contact model could fit the experimental data better.

3.2.2. Coordination Number

Figure 10 illustrates the DEM-simulated relationship of the coordination number CN versus the bulk density ratio of the pellet to the intact wood, BDp/BDw, during compaction. The coordination number CN increased nonlinearly with increasing BDp/BDw, being faster and nonlinear for BDp/BDw < 1 (as free spaces closed between particles) and slower and almost linear for BDp/BDw > 1 (as there remained no free spaces between particles, the surface pores on the sawdust particles and intraparticle pores being closed). In this figure, three values of CN of relaxed pellets (12.9, 12.2, and 10.5) were obtained in the process of unloading and relaxation of the bulk of particles compacted to σz = 120, 60, and 48 MPa, respectively. The values of the BDp/BDw ratio, corresponding to the indicated values of CN, respectively, covered the entire range of variability of the experimental values of the BDp/BDw ratio of 1.62 ± 0.19, 1.32 ± 0.14, and 1.13 ± 0.08, determined as the averaged values for the three wood materials (Table 2) and for three cases: (1) MC = 8%, σz = 120 MPa; (2) MC = 8%, σz = 60 MPa and MC = 20%, σz = 120 MPa; and (3) MC = 20%, σz =120 MPa.
The diametral compression tests were simulated for the three selected values of the coordination number of the pellets (Figure 11a). The shapes of the stress–deformation relationships changed slightly with the change in CN. All shapes were typical for ductile behaviour with progressive loading and gradual collapse. The tensile strength σf increased slightly faster than linearly as the coordination number increased (Figure 11b). However, Gilvari et al. [40] found that the tensile strength of pellets increased linearly with the coordination number.

3.2.3. Bond Radius Estimation

As stated previously, the mean size of bonds of 20 μm was the baseline input data for searching for the bond radius in the DEM simulations. Examples of simulations of the diametral compression performed for different values of the bond radius and constant values of the bond elasticity modulus and bond strength (Eb = 120 MPa and σc = 36 MPa) illustrate the strong dependence of the stress–deformation relationships on the bond radius (Figure 12a). However, simulations performed for the bond parameters fulfilling the conditions E b / r b 2 = constant and E b / σ c = constant had almost identical stress–deformation relationships (Figure 12b). This similarity means that the model can be easily recalibrated to different values of micro-parameters to simulate other stages of compaction (coordination number), other bond radii, or other materials (elastic parameters). This indicates an important need for precise experimental verification of the micro-parameters of bonds to obtain reliable modelling results.
Simulations performed with constant values of Eb and σc and variable bond radii indicated that the breakage strength σf (Figure 13a) and the deformation at failure, ΔLf/D (Figure 13b), increased almost linearly with the bond cross-sectional area A. Experimental values of σf and ΔLf/D for oak pellets (MC = 8% and 20% and σz = 60, 120 MPa) presented in this figure served as reference data to find the best-fitting values of Eb and σc for the selected value of rb of 20 μm and the particular values of MC and σz. By taking into account the fact that the value of the bond radius approximated from SEM images provides a similarity between the simulated and experimental stress–deformation relationships and the breakage profiles, a bond radius of 20 μm was set as the fixed material parameter for all further DEM simulations.

3.2.4. Bond Elasticity and Strength

The bond elasticity modulus Eb and strength σc were determined as values providing the minimum RRMSE of the DEM approximation of the experimental stress–deformation relationship. Details of the procedure for searching for the best-fitting parameters are exemplarily illustrated for oak pellets with MC = 8% compacted to σz = 120 MPa (Figure 14). Generally, 8–10 simulations, divided into two or three groups with sequentially changed Eb (Figure 14a) and σc (Figure 14b) as the constant and variable parameters, were sufficient to find the set of Eb and σc, bringing the RRMSE relatively close to its global minimum. This procedure was applied separately to each variant of the experimental relationship. Eighty simulations were performed to find the best fit for all experimental cases.

4. Simulation Results

Figure 15 presents examples of experimental relationships between the tension stress and the diametral deformation σ1L/D) fitted by the DEM approximations. The DEM simulations approximated very well the individual experimental stress–deformation relationships that were the closest to the mean tensile strength (σf, ΔLf/D).
The best fit with the experimental relationships was obtained for CN corresponding to the degree of compaction of pellets, as given in Table 5.
The best-fitting CN increased with increasing BDp/BDw ratio. This can be clearly seen for the oak pellets. The bond elasticity modulus Eb and its strength σc decreased as MC increased and σz decreased. The range of variability of the bond elasticity modulus applied in the DEM simulations (Table 5) from 9.2 MPa (birch, MC = 20%, and σz = 60 MPa) to 120 MPa (oak, MC = 8%, and σz = 120 MPa) was consistent with a range of values (3–146 MPa) determined experimentally for wood pellets in the puncture test [51,52,53]. The approximately twofold decrease in the bond elasticity modulus with an increase in MC from 8% to 20% applied in our simulations fitted very well to the rate of decrease in pellet elasticity determined experimentally by Gallego et al. [51]. The order of magnitude of the bond elasticity modulus (107–108 Pa) and the corresponding range of the normal stiffness coefficient of the bond (1010–1011 N m−3) applied in our study were similar to the values of the bond elasticity and stiffness applied in the DEM modelling of loading of pinewood chip briquettes by Xia et al. [38], the durability of wood pellets found by Mahajan et al. [39], and the breakage of biomass pellets studied by Gilvari et al. [40]. The maximum bond tensile strength σc of 36 MPa applied in our study for oak (MC = 8% and σz = 120 MPa) corresponded to that applied by Gilvari et al. [40] (σc = 35 MPa) and Mahajan et al. [39] (σc = 40 MPa) to model pellet behaviour. The ~30% decrease in the bond strength with an increase in MC from 8% to 20% was consistent with the results of Whittaker and Shield [2] and Li and Liu [9], who indicated that an excessively high MC reduces the binding forces between particles.
The high value of the coordination number applied in simulations (10.5–12.9) considerably increased the strength of pellets generated by elastic bonds and friction in particle contacts and promoted the rounded shape of the σ1L/D) relationship typical for the ductile breakage mode. As shown by Horabik et al. [54], σc/Eb was <0.1 for the semi-brittle breakage mode and >0.15 for ductile breakage. The σc/Eb values found in this study ranged from 0.25 to 0.5, which is typical for the ductile breakage mode.
Figure 16 presents a comparison of experimental and simulated patterns of cracks of oak pellets.
Cracks were initiated on the flat surfaces of pellets, and progressive deformation developed along the direction of loading toward the interior of the pellet. Two types of breakage patterns were observed in the experiments: a single crack slightly inclined against the direction of loading and diffuse small cracks connected into branches. The first type occurred for pellets of the highest and the lowest strength and the second for the intermediate strengths. The DEM simulations indicated that the breakage profile was strongly influenced by the coordination number and was very slightly influenced by the elastic parameters of the bonds. For CN = 12.9, the DEM simulations reproduced a single-crack breakage well. Simulations performed for CN = 12.2 reproduced the development of diffuse cracks. In general, cracks grew in the central part of the pellet and did not reach the outer surface of the pellet. The earliest appearing cracks were observed in simulations performed for CN = 10.5. In the first stage, breakage occurred as a set of small, separated cracks. As the deformation progressed, these cracks gradually merged into a single zone, which finally separated the pellet into two parts.

5. General Remarks

Bonding mechanisms that join sawdust particles are crucial for the mechanical strength of pellets. The diametral compression test provides a convenient method for determining the tensile strength of pellets. The DEM is a powerful tool for investigating the influence of particle-scale properties on bulk material behaviour. Bonded spheres applied in the DEM modelling proved to be useful for reproducing the breakage behaviour of wood pellets. The main novelty of this research was the good quality of the DEM reproduction of particular phases of pellet formation and their breakage behaviour during diametral compression. Pellet formation was reproduced with acceptable quality, while the response of pellets to diametral compression was reproduced with high accuracy (RRMSE < 0.12). The DEM simulations reproduced the shape of the stress–deformation relationships, tensile strength, and breakage profiles of pellets produced from all three wood materials at two levels of MC and at two levels of compaction pressure. The order of magnitude of the elastic modulus and the strength of bonds were consistent with previously published research results on the application of the BPM for testing the breakage strength of wood pellets [40].
The range of deformation of particles, which was necessary to model the behaviour of sawdust particles during the pelletisation process, was far beyond the typical range of deformation considered in the DEM simulations. To overcome this limitation, a uniform and high density similar to the bulk density of the intact oak wood was set as the bulk density of all sawdust particles (regardless of the type of wood). Similarly, Xia et al. [38] used wood bulk density as the particle density in their DEM simulations. By this means, the large bulk deformation of sawdust particles during pelletisation was considered to result from the reduction of free spaces between particles and compaction of internal pores inside sawdust particles reproduced by a very high level of overlap (δn/r~0.33). The coordination number of relaxed pellets in the range of 10.5–12.9 was consistent with the values obtained by Nordström et al. [27] in their DEM modelling of the formation of microcrystalline cellulose tablets under pressure in the range of 100–300 MPa.
To model the stress–deformation relationship during diametral compression of pellets with a bond radius equal to the radius of contact between particles, the stiffness and strength of the bond should be unrealistically low compared to the experimental findings [50]. It is reasonable to expect that bonding develops only on a fraction of the contact area of the largely deformed sawdust particles. Therefore, supported by the results of the SEM images of the pellets and the results of the calibration procedure, the bond radius applied to the DEM simulations was assumed to be a small fraction (~0.16) of the contact radius. The assumed value of the initial particle density, along with placing the bond radius as a small fraction of the contact radius, allowed the proper modelling of both pelletisation and diametral compression. The simulations yielded good qualitative agreement of the crack profile with the experimental data, reproduction of the stress–deformation relationship during diametral compression, and agreement of the values of the simulated bond elasticity with the experimental data.
Some problems encountered during the DEM modelling of sawdust and their products may have resulted from very limited access to verified experimental data for the micro-properties of sawdust particles [51]. Therefore, more experimental data for micro-properties are necessary to facilitate further development of DEM modelling of these materials.

6. Conclusions

An experimental study and numerical simulations were performed for pressure compaction of sawdust and diametral compression tests of sawdust pellets. The simulation results reproduced the results of laboratory testing well and provided deeper insight into particle–particle bonding mechanisms. The main novelty of the present study was the successful simulation of the diametral compression test of wood pellets and reproduction of the profiles of ductile crack formation using the DEM with the implemented BPM. The following can be concluded from this study:
  • The stress–deformation relationship during pelletisation of wood sawdust comprising a large range of particle deformations was successfully modelled using a DEM equipped with a linear hysteretic contact model. Good agreement between the simulated range of change in the bulk density during pelletisation with experimental data was obtained because of the application density of the intact wood as the input density of particles in the DEM simulations;
  • The highest tensile strength was obtained for oak and the lowest for birch pellets. For all materials, the tensile strength was the highest for MC = 8% of sawdust compacted under a pressure of 120 MPa and the lowest for MC = 20% of sawdust compacted under a pressure of 60 MPa;
  • The breakage processes of pellets of all tested materials were successfully simulated using the DEM with the BPM;
  • All pellets exhibited a ductile breakage mode characterised by a smooth and round stress–deformation relationship without any sudden drops. Cracks were initiated in locations close to the centre of the pellet, and progressive deformation developed in the direction of loading and toward the interior of the pellet;
  • Applying values of the bond elasticity modulus Eb and the tensile strength σc fulfilling the condition σc/Eb > 0.25 in the DEM simulations allowed the stress–deformation relationship and crack formation to be reproduced well for all studied pellets.

Author Contributions

Conceptualisation, J.H.; methodology, J.H., J.W. and P.P.; software, P.P.; validation, R.K.; formal analysis, M.M.; investigation, M.B., G.J., A.A., C.P. and M.S.; resources, J.H.; data curation, J.W.; writing—original draft preparation, J.H.; writing—review and editing, G.J., A.A. and M.M.; visualisation, R.K.; supervision, J.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data are available from the corresponding author for reasonable requests.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Kaliyan, N.; Morey, R.V. Natural binders and solid bridge type binding mechanisms in briquettes and pellets made from corn stover and switchgrass. Bioresour. Technol. 2010, 101, 1082–1090. [Google Scholar] [CrossRef]
  2. Whittaker, C.; Shield, I. Factors affecting wood, energy grass and straw pellet durability—A review. Renew. Sustain. Energy Rev. 2017, 71, 1–11. [Google Scholar] [CrossRef]
  3. Adapa, P.K.; Tabil, L.G.; Schoenau, G.J. Factors affecting the quality of biomass pellet for biofuel and energy analysis of pelleting process. Int. J. Agric. Biol. Eng. 2013, 6, 1–12. [Google Scholar] [CrossRef]
  4. Stasiak, M.; Molenda, M.; Bańda, M.; Wiącek, J.; Parafiniuk, P.; Gondek, E. Mechanical and combustion properties of sawdust—Straw pellets blended in different proportions. Fuel Process. Technol. 2017, 156, 366–375. [Google Scholar] [CrossRef]
  5. Hehar, G.; Fasina, O.; Adhikari, S.; Fulton, J. Ignition and volatilization behavior of dust from loblolly pine wood. Fuel Process. Technol. 2014, 127, 117–123. [Google Scholar] [CrossRef]
  6. Shaw, M.D.; Karunakaran, C.; Tabil, L.G. Physicochemical characteristics of densified untreated and steam exploded poplar wood and wheat straw grinds. Biosyst. Eng. 2009, 103, 198–207. [Google Scholar] [CrossRef]
  7. Kretschmann, D.E. Mechanical Properties of Wood. In Wood Handbook: Wood as Engineering Material, 2nd ed.; U.S. Department of Agriculture: Madison, WI, USA, 2010; pp. 5.1–5.46. [Google Scholar]
  8. Carone, M.T.; Pantaleo, A.; Pellerano, A. Influence of process parameters and biomass characteristics on the durability of pellets from the pruning residues of Olea europaea L. Biomass Bioenergy 2011, 35, 402–410. [Google Scholar] [CrossRef]
  9. Li, Y.; Liu, H. High-pressure densification of wood residues to form an upgraded fuel. Biomass Bioenergy 2000, 19, 177–186. [Google Scholar] [CrossRef]
  10. EN ISO 17831-1:2015. Solid Biofuels-Determination of Mechanical Durability of Pellets and Briquettes—Part 1: Pellets; International Organization for Standardization: Geneva, Switzerland, 2015. [Google Scholar]
  11. Larsson, S.H.; Samuelsson, R. Prediction of ISO 17831-1:2015 mechanical biofuel pellet durability from single pellet characterization. Fuel Process. Technol. 2017, 163, 8–15. [Google Scholar] [CrossRef]
  12. Oveisi, E.; Lau, A.; Sokhansanj, S.; Lim, C.J.; Bi, X.; Larsson, S.H.; Melin, S. Breakage behavior of wood pellets due to free fall. Powder Technol. 2013, 235, 493–499. [Google Scholar] [CrossRef]
  13. Jonsén, P.; Häggblad, H.-Å.; Sommer, K. Tensile strength and fracture energy of pressed metal powder by diametral compression test. Powder Technol. 2007, 176, 148–155. [Google Scholar] [CrossRef]
  14. Timoshenko, S.P.; Goodier, J.N. Theory of Elasticity, 3rd ed.; McGraw-Hill: New York, NY, USA, 1970. [Google Scholar]
  15. Andreev, G. A review of the Brazilian test for rock tensile strength determination. Part I: Calculation formula. Min. Sci. Technol. 1991, 13, 445–456. [Google Scholar] [CrossRef]
  16. Labuz, J.F.; Cattaneo, S.; Chen, L.-H. Acoustic emission at failure in quasi-brittle materials. Constr. Build. Mater. 2001, 15, 225–233. [Google Scholar] [CrossRef]
  17. Alderborn, G. A Novel Approach to Derive a Compression Parameter Indicating Effective Particle Deformability. Pharm. Dev. Technol. 2003, 8, 367–377. [Google Scholar] [CrossRef]
  18. Nordström, J.; Persson, A.-S.; Lazorova, L.; Frenning, G.; Alderborn, G. The degree of compression of spherical granular solids controls the evolution of microstructure and bond probability during compaction. Int. J. Pharm. 2013, 442, 3–12. [Google Scholar] [CrossRef] [Green Version]
  19. Pai, D.A.; Hayes, A.A.; Okos, M.R. Modeling pharmaceutical compacts tensile strength based on viscoelastic properties. Powder Technol. 2013, 239, 441–450. [Google Scholar] [CrossRef]
  20. Han, L.; Elliott, J.; Bentham, A.; Mills, A.; Amidon, G.; Hancock, B. A modified Drucker-Prager Cap model for die compaction simulation of pharmaceutical powders. Int. J. Solids Struct. 2008, 45, 3088–3106. [Google Scholar] [CrossRef] [Green Version]
  21. Okot, D.K.; Bilsborrow, P.E.; Phan, A.N. Effects of operating parameters on maize COB briquette quality. Biomass Bioenergy 2018, 112, 61–72. [Google Scholar] [CrossRef] [Green Version]
  22. Song, X.; Zhang, S.; Wu, Y.; Cao, Z. Investigation on the properties of the bio-briquette fuel prepared from hydrothermal pretreated cotton stalk and wood sawdust. Renew. Energy 2020, 151, 184–191. [Google Scholar] [CrossRef]
  23. Ahmad, M.Z.; Akhter, S.; Anwar, M.; Rahman, M.; Siddiqui, M.A.; Ahmad, F.J. Compatibility and compressibility studies of Assam Bora rice starch. Powder Technol. 2012, 224, 281–286. [Google Scholar] [CrossRef]
  24. Mitchell, W.R.; Forny, L.; Althaus, T.; Dopfer, D.; Niederreiter, G.; Palzer, S. Compaction of food powders: The influence of material properties and process parameters on product structure, strength, and dissolution. Chem. Eng. Sci. 2017, 167, 29–41. [Google Scholar] [CrossRef]
  25. Cundall, P.A.; Strack, O.D.L. A discrete numerical model for granular assemblies. Géotechnique 1979, 29, 47–65. [Google Scholar] [CrossRef]
  26. Azéma, E.; Sánchez, P.; Scheeres, D.J. Scaling behavior of cohesive self-gravitating aggregates. Phys. Rev. E 2018, 98, 030901. [Google Scholar] [CrossRef]
  27. Nordström, J.; Alderborn, G.; Frenning, G. Compressibility and tablet forming ability of bimodal granule mixtures: Experiments and DEM simulations. Int. J. Pharm. 2018, 540, 120–131. [Google Scholar] [CrossRef] [PubMed]
  28. Balevičius, R.; Sielamowicz, I.; Mróz, Z.; Kačianauskas, R. Effect of rolling friction on wall pressure, discharge velocity and outflow of granular material from a flat-bottomed bin. Particuology 2012, 10, 672–682. [Google Scholar] [CrossRef]
  29. Horabik, J.; Wiącek, J.; Parafiniuk, P.; Stasiak, M.; Bańda, M.; Molenda, M. Tensile strength of pressure-agglomerated potato starch determined via diametral compression test: Discrete element method simulations and experiments. Biosyst. Eng. 2019, 183, 95–109. [Google Scholar] [CrossRef]
  30. Rackl, M.; Top, F.; Molhoek, C.P.; Schott, D.L. Feeding system for wood chips: A DEM study to improve equipment performance. Biomass Bioenergy 2017, 98, 43–52. [Google Scholar] [CrossRef]
  31. Xia, Y.; Stickel, J.J.; Jin, W.; Klinger, J. A Review of Computational Models for the Flow of Milled Biomass Part I: Discrete-Particle Models. ACS Sustain. Chem. Eng. 2020, 8, 6142–6156. [Google Scholar] [CrossRef]
  32. Pachón-Morales, J.; Perré, P.; Casalinho, J.; Do, H.; Schott, D.; Puel, F.; Colin, J. Potential of DEM for investigation of non-consolidated flow of cohesive and elongated biomass particles. Adv. Powder Technol. 2020, 31, 1500–1515. [Google Scholar] [CrossRef]
  33. Ilic, D.; Williams, K.; Ellis, D. Assessment of biomass bulk elastic response to consolidation. Chem. Eng. Res. Des. 2018, 135, 185–196. [Google Scholar] [CrossRef]
  34. Guo, Y.; Chen, Q.; Xia, Y.; Westover, T.; Eksioglu, S.; Roni, M. Discrete element modeling of switchgrass particles under compression and rotational shear. Biomass Bioenergy 2020, 141, 105649. [Google Scholar] [CrossRef]
  35. Potyondy, D.O.; Cundall, P.A. A bonded-particle model for rock. Int. J. Rock Mech. Min. Sci. 2004, 41, 1329–1364. [Google Scholar] [CrossRef]
  36. Jiménez-Herrera, N.; Barrios, G.K.P.; Tavares, L.M. Comparison of breakage models in DEM in simulating impact on particle beds. Adv. Powder Technol. 2018, 29, 692–706. [Google Scholar] [CrossRef]
  37. Langston, P.; Kennedy, A.; Constantin, H. Discrete element modelling of flexible fibre packing. Comput. Mater. Sci. 2015, 96, 108–116. [Google Scholar] [CrossRef]
  38. Xia, Y.; Lai, Z.; Westover, T.; Klinger, J.; Huang, H.; Chen, Q. Discrete element modeling of deformable pinewood chips in cyclic loading test. Powder Technol. 2019, 345, 1–14. [Google Scholar] [CrossRef]
  39. Mahajan, A.; Dafnomilis, I.; Hancock, V.; Lodewijks, G.; Schott, D. Assessing the representativeness of durability tests for wood pellets by DEM Simulation—Comparing conditions in a durability test with transfer chutes. EPJ Web Conf. 2017, 140, 15004. [Google Scholar] [CrossRef] [Green Version]
  40. Gilvari, H.; De Jong, W.; Schott, D.L. Breakage behavior of biomass pellets: An experimental and numerical study. Comput. Part. Mech. 2020, 1–14. [Google Scholar] [CrossRef]
  41. Xia, Y.; Chen, F.; Klinger, J.L.; Kane, J.J.; Bhattacharjee, T.; Seifert, R.; Ajayi, O.O.; Chen, Q. Assessment of a tomography-informed polyhedral discrete element modelling approach for complex-shaped granular woody biomass in stress consolidation. Biosyst. Eng. 2021, 205, 187–211. [Google Scholar] [CrossRef]
  42. Molenda, M.; Horabik, J.; Parafiniuk, P.; Oniszczuk, A.; Bańda, M.; Wajs, J.; Gondek, E.; Chutkowski, M.; Lisowski, A.; Wiącek, J.; et al. Mechanical and Combustion Properties of Agglomerates of Wood of Popular Eastern European Species. Materials 2021, 14, 2728. [Google Scholar] [CrossRef] [PubMed]
  43. Sridharan, A.; Rao., G.V. Pore Size Distributions of Soils From Mercury Intrusion Porosimeter Data. Soil Sci. Soc. Am. J. 1972, 36, 980–981. [Google Scholar] [CrossRef]
  44. PubChem—National Center for Biotechnology Information, U.S. National Library of Medicine. Available online: https://pubchem.ncbi.nlm.nih.gov/compound/5988#section=Density (accessed on 10 June 2021).
  45. Walton, O.R.; Braun, R.L. Viscosity, granular-temperature, and stress calculations for shearing assemblies of inelastic, frictional disks. J. Rheol. 1986, 30, 949–980. [Google Scholar] [CrossRef]
  46. Luding, S. Shear flow modeling of cohesive and frictional fine powder. Powder Technol. 2005, 158, 45–50. [Google Scholar] [CrossRef]
  47. Thornton, C.; Ning, Z. A theoretical model for the stick/bounce behaviour of adhesive, elastic-plastic spheres. Powder Technol. 1998, 99, 154–162. [Google Scholar] [CrossRef]
  48. EDEM Solutions. EDEM 2018.2 Documentation; EDEM: Edinburgh, UK, 2018. [Google Scholar]
  49. Tsuji, Y.; Tanaka, T.; Ishida, T. Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe. Powder Technol. 1992, 71, 239–250. [Google Scholar] [CrossRef]
  50. Huang, Y.J.; Nydal, O.J.; Yao, B. Time step criterions for nonlinear dense packed granular materials in time-driven method simulations. Powder Technol. 2014, 253, 80–88. [Google Scholar] [CrossRef]
  51. Gallego, E.; Fuentes, J.M.; Ruiz, Á.; Hernández-Rodrigo, G.; Aguado, P.; Ayuga, F. Determination of mechanical properties for wood pellets used in DEM simulations. Int. Agrophysics 2020, 34, 485–494. [Google Scholar] [CrossRef]
  52. Kocsis, Z.; Csanády, E. Investigation on the Mechanics of Wood Pellet Production from Sawdust and Chips; University of West: Sopron, Hungary, 2017; pp. 1–31. [Google Scholar]
  53. Graham, S.; Eastwick, C.; Snape, C.; Quick, W. Mechanical degradation of biomass wood pellets during long term stockpile storage. Fuel Process. Technol. 2017, 160, 143–151. [Google Scholar] [CrossRef]
  54. Horabik, J.; Wiącek, J.; Parafiniuk, P.; Stasiak, M.; Bańda, M.; Kobyłka, R.; Molenda, M. Discrete Element Method Modelling of the Diametral Compression of Starch Agglomerates. Materials 2020, 13, 932. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Stress components in the pellet during the diametral compression test.
Figure 1. Stress components in the pellet during the diametral compression test.
Materials 14 03273 g001
Figure 2. Dependencies of the pore volume on the pore radius for (a) the studied wood and (b) pellets. Average curves from two replicates are presented.
Figure 2. Dependencies of the pore volume on the pore radius for (a) the studied wood and (b) pellets. Average curves from two replicates are presented.
Materials 14 03273 g002
Figure 3. Pore size distribution functions for (a) the studied wood and (b) pellets. Average curves from two replicates are presented.
Figure 3. Pore size distribution functions for (a) the studied wood and (b) pellets. Average curves from two replicates are presented.
Materials 14 03273 g003
Figure 4. Representative SEM images (at 1000× magnification) of the surfaces of cross-sections of the wood blocks, sawdust particles, and pellets. White ellipses surround examples of potential particle–particle bonds. The 40 µm scale bar shown in the oak pellet image is valid for all other materials.
Figure 4. Representative SEM images (at 1000× magnification) of the surfaces of cross-sections of the wood blocks, sawdust particles, and pellets. White ellipses surround examples of potential particle–particle bonds. The 40 µm scale bar shown in the oak pellet image is valid for all other materials.
Materials 14 03273 g004
Figure 5. Stress–deformation dependencies during (a) the compaction–unloading cycle of sawdust and (b) the diametral compression of pellets.
Figure 5. Stress–deformation dependencies during (a) the compaction–unloading cycle of sawdust and (b) the diametral compression of pellets.
Materials 14 03273 g005
Figure 6. Linear hysteretic contact model with adhesion.
Figure 6. Linear hysteretic contact model with adhesion.
Materials 14 03273 g006
Figure 7. Linear hysteretic contact model with the BPM.
Figure 7. Linear hysteretic contact model with the BPM.
Materials 14 03273 g007
Figure 8. Stages of the DEM simulations: (a) filling, (b) compaction, (c) relaxation, and (d) diametral compression.
Figure 8. Stages of the DEM simulations: (a) filling, (b) compaction, (c) relaxation, and (d) diametral compression.
Materials 14 03273 g008
Figure 9. Fitting the compaction and unloading cycle σzh/h0) of sawdust.
Figure 9. Fitting the compaction and unloading cycle σzh/h0) of sawdust.
Materials 14 03273 g009
Figure 10. Comparison of the coordination number CN versus the bulk density ratio BDp/BDw of sawdust during compaction and for relaxed pellets.
Figure 10. Comparison of the coordination number CN versus the bulk density ratio BDp/BDw of sawdust during compaction and for relaxed pellets.
Materials 14 03273 g010
Figure 11. Impact of the coordination number CN on (a) the shape of the stress–deformation relationship during diametral compression and (b) the tensile strength for rb = 20 μm, Eb = 120 MPa, and σc = 36 MPa.
Figure 11. Impact of the coordination number CN on (a) the shape of the stress–deformation relationship during diametral compression and (b) the tensile strength for rb = 20 μm, Eb = 120 MPa, and σc = 36 MPa.
Materials 14 03273 g011
Figure 12. Stress–deformation relationship during diametral compression: (a) impact of the bond radius for Eb = 120 MPa, σc = 36 MPa, and CN = 12.9; (b) impact of the bond radius compensated by the constant value of the following proportions: Ebrb2 = constant and Eb/σc = constant.
Figure 12. Stress–deformation relationship during diametral compression: (a) impact of the bond radius for Eb = 120 MPa, σc = 36 MPa, and CN = 12.9; (b) impact of the bond radius compensated by the constant value of the following proportions: Ebrb2 = constant and Eb/σc = constant.
Materials 14 03273 g012
Figure 13. DEM simulations fitting the range of experimental values for oak as influenced by the bond cross-sectional area A for CN = 12.9: (a) the tensile strength, σf; (b) deformation at failure, ΔLf/D.
Figure 13. DEM simulations fitting the range of experimental values for oak as influenced by the bond cross-sectional area A for CN = 12.9: (a) the tensile strength, σf; (b) deformation at failure, ΔLf/D.
Materials 14 03273 g013
Figure 14. RRMSE of fitting the experimental stress–deformation relationship during diametral compression of oak pellets (MC = 8% and σz = 120 MPa) performed for CN = 12.9 and rb = 20 μm as influenced by (a) the bond elasticity modulus Eb and (b) the bond strength σc.
Figure 14. RRMSE of fitting the experimental stress–deformation relationship during diametral compression of oak pellets (MC = 8% and σz = 120 MPa) performed for CN = 12.9 and rb = 20 μm as influenced by (a) the bond elasticity modulus Eb and (b) the bond strength σc.
Materials 14 03273 g014
Figure 15. Fitting the σ1(ΔL/D) relationships of the diametral compression for birch, oak, and pine. The bars indicate the SD of the mean value.
Figure 15. Fitting the σ1(ΔL/D) relationships of the diametral compression for birch, oak, and pine. The bars indicate the SD of the mean value.
Materials 14 03273 g015
Figure 16. Experimental and simulated breakage profiles of oak pellets for ΔL/D = 0.1.
Figure 16. Experimental and simulated breakage profiles of oak pellets for ΔL/D = 0.1.
Materials 14 03273 g016
Table 1. Density of the solid phase of the intact wood and of the sawdust measured by helium pycnometry.
Table 1. Density of the solid phase of the intact wood and of the sawdust measured by helium pycnometry.
MaterialWood Solid Phase (Including Closed Pores)Sawdust
Density ρw (kg m−3)Density ρs (kg m−3)
Birch1369.8 ± 0.3 a1465.3 ± 0.8 a
Oak1068.4 ± 1.6 b1459.9 ± 1.1 b
Pine1426.6 ± 0.3 c1468.7 ± 0.6 c
Mean values in a column followed by the same letter are not significantly different at the 5% level.
Table 2. Bulk density and porosity of the studied intact wood, sawdust, and pellets at various values of MC.
Table 2. Bulk density and porosity of the studied intact wood, sawdust, and pellets at various values of MC.
MC
(%)
Intact WoodSawdust (Initial)Pellets
BDw
(kg m−3)
pw
(%)
BDs
(kg m−3)
ps
(%)
Compacted to 60 MPaCompacted to 120 MPa
BDp
(kg m−3)
pp (%)BDp
(kg m−3)
pp (%)
Birch8550 ± 2 a62.5 ± 0.2 a301.7 ± 6.4 a79.4 ± 1.7 a701.2 ± 9.9 a52.1 ± 0.7 a843.9 ± 4.8 a42.4 ± 0.2 a
Birch20 272.1 ± 9.7 b81.4 ± 2.9 b634.6 ± 6.7 b56.7 ± 0.6 b781.9 ± 22 b46.6 ± 1.3 b
Oak8675 ± 3 b53.8 ± 0.2 b312.6 ± 5.4 a78.6 ± 1.4 a824.5 ± 10.4 c43.5 ± 0.5 c1005.8 ± 11.8 c31.1 ± 0.4 c
Oak20 267.8 ± 6.3 b81.6 ± 1.9 b700.3 ± 25.5 a52.0 ± 1.9 a826.4 ± 35.5 a43.3 ± 1.9 a
Pine8510 ± 2 c65.2 ± 0.2 c267.9 ± 6.7 b81.8 ± 2.0 b779.6 ± 12.9 d46.9 ± 0.8 d937.8 ± 17.8 d36.2 ± 0.7 d
Pine20 250.9 ± 3.4 c82.9 ± 1.1 c610.9 ± 44.2 ab58.4 ± 4.2 ab673.1 ± 33.3 b54.2 ± 2.7 e
Mean values in a column followed by the same letter are not significantly different at the 5% level.
Table 3. Tensile strength of pellets.
Table 3. Tensile strength of pellets.
Materialσz (MPa)MC (%)ΔLf/Dσf (MPa)
Birch6080.056 ± 0.005 a0.112 ± 0.012 a
Birch12080.056 ± 0.003 a0.307 ± 0.009 b
Birch60200.057 ± 0.005 a0.081 ± 0.008 a
Birch120200.061 ± 0.008 b0.183 ± 0.018 c
Oak6080.056 ± 0.001 a0.511 ± 0.069 d
Oak12080.061 ± 0.003 b1.346 ± 0.068 e
Oak60200.059 ± 0.004 c0.277 ± 0.039 b
Oak120200.058 ± 0.002 c0.683 ± 0.067 d
Pine6080.064 ± 0.008 b0.262 ± 0.019 b
Pine12080.057 ± 0.003 a0.608 ± 0.034 d
Pine60200.068 ± 0.006 d0.085 ± 0.025 a
Pine120200.082 ± 0.002 e0.182 ± 0.071 c
Mean values in a column followed by the same letter are not significantly different at the 5% level.
Table 4. DEM simulation parameters.
Table 4. DEM simulation parameters.
ParameterSymbolValue
Container:
Diameter (mm)Dc10
Height (mm)Hc25
Solid density (kg m–3)ρ7800
Young’s modulus (MPa)E1.561 × 106
Poisson’s ratioν0.3
Particles:
Number 40,000
Mean particle radius (mm)r0.2
SD of particle radius (mm)rsd0.02
Particle radius range (mm) 0.14–0.26
Particle solid density (kg m–3)ρ680
Young’s modulus (MPa)E1.57 × 104
Poisson’s ratioν0.35
Yield strength (MPa)py100; 150
Mean loading (plastic) stiffness (N m–1)k11 × 105
Mean unloading (elastic) stiffness (N m–1)k24 × 105
Mean adhesion stiffness (N m–1)kc0; 600
Restitution coefficiente0.5
Particle–particle friction coefficientμp-p0.5
Particle–wall friction coefficientμp-w0.15
Rolling friction coefficientmr0.01
Bond radius (μm)rb10; 20; 30; 40
Bond tension strength (MPa)σc3–40
Bond shear strength (MPa)τc3–40
Bond Young’s modulus (MPa)Eb5–120
Table 5. Best-fitting parameters of the DEM simulations.
Table 5. Best-fitting parameters of the DEM simulations.
Materialσz (MPa)MC (%)CNEb (MPa)σc (MPa)RRMSE
Birch60812.212.84.60.027
Birch120812.236.610.10.039
Birch602012.29.240.037
Birch1202012.222.48.050.039
Oak60812.260200.027
Oak120812.9120360.111
Oak602010.552.8130.089
Oak1202012.273.2240.116
Pine60812.227.213.50.052
Pine120812.950.413.50.028
Pine602012.28.64.30.053
Pine1202012.219.79.90.049
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Horabik, J.; Bańda, M.; Józefaciuk, G.; Adamczuk, A.; Polakowski, C.; Stasiak, M.; Parafiniuk, P.; Wiącek, J.; Kobyłka, R.; Molenda, M. Breakage Strength of Wood Sawdust Pellets: Measurements and Modelling. Materials 2021, 14, 3273. https://doi.org/10.3390/ma14123273

AMA Style

Horabik J, Bańda M, Józefaciuk G, Adamczuk A, Polakowski C, Stasiak M, Parafiniuk P, Wiącek J, Kobyłka R, Molenda M. Breakage Strength of Wood Sawdust Pellets: Measurements and Modelling. Materials. 2021; 14(12):3273. https://doi.org/10.3390/ma14123273

Chicago/Turabian Style

Horabik, Józef, Maciej Bańda, Grzegorz Józefaciuk, Agnieszka Adamczuk, Cezary Polakowski, Mateusz Stasiak, Piotr Parafiniuk, Joanna Wiącek, Rafał Kobyłka, and Marek Molenda. 2021. "Breakage Strength of Wood Sawdust Pellets: Measurements and Modelling" Materials 14, no. 12: 3273. https://doi.org/10.3390/ma14123273

APA Style

Horabik, J., Bańda, M., Józefaciuk, G., Adamczuk, A., Polakowski, C., Stasiak, M., Parafiniuk, P., Wiącek, J., Kobyłka, R., & Molenda, M. (2021). Breakage Strength of Wood Sawdust Pellets: Measurements and Modelling. Materials, 14(12), 3273. https://doi.org/10.3390/ma14123273

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop