Modeling the devolatilization and fragmentation of biomass pellets with the bonded particle method for fluidized bed applications

By extending the bonded particle method, the major structural changes during the devolatilization of a wood pellet in a fluidized bed and the resulting mechanical behavior have been successfully reproduced. A comparison with experiments from the literature showed that the implemented particle-based pyrolysis model enables predicting the entire pellet’s kinetics with a high agreement. The developed shrinkage model for the particles and bonds further allowed to emulate the reported formation of a large-scale pore network inside the pellet. The simulation of a radial compression test with the predicted structure showed good agreement with experimental data and could confirm the importance of the pores for the mechanical behavior. The results demonstrated that the large pores cause the fragmentation of agglomerates already at low mechanical loads which could promote attrition. In general, the results have shown that the developed extension of the bonded particle method allows studying and predicting the behavior of a single pellet during conversion inside a fluidized bed gasification reactor in more detail.


Introduction
Fluidized bed gasification of biomass is a promising approach to produce synthesis gas sustainably. Synthesis gas, or short syngas, consists mainly of carbon monoxide, hydrogen, methane, and carbon dioxide. It can be used as an energy source or for producing chemicals, biofuels, and hydrogen. Syngas can be produced from biomass by gasification. Gasification is a thermochemical conversion process typically run at elevated temperatures of above 700 • C. By partially oxidizing the biomass, the process can run autothermally, allowing a climate-neutral conversion. To improve fuel properties, reduce transportation costs and simplify material feeding, the biomass is usually pelletized. Often fluidized bed reactors are used for gasification to achieve high efficiency and enable greater flexibility in the materials and capacity. One investigated gasification technology is chemical looping gasification (CLG). CLG allows reaching a high concentration of the produced syngas and in-bed tar removal by providing lattice oxygen from a circulating oxygen carrier. The oxygen carrier is undergoing oxidation-reduction cycles as it circulates between two fluidized bed reactors. In a so-called air reactor, the carrier is oxidized and then transported to the fuel reactor, where it provides heat for the conversion and is reduced. The separation into two reactors allows supplying oxygen from the air without diluting the product gas.
However, a common problem of all approaches using fluidized bed reactors is the increased attrition and breakage of the material inside the reactor [35,42]. In the case of the biomass pellets, their attrition and breakage can result in fines leaving the reactor unreacted, carbon slipping into the air reactor, or a change in their mixing behavior. [20] have shown that after devolatilization in a fluidized bed reactor, the pellet shape typically remains, but large internal cavities and macro-scale pores are formed. As a result, the pellets' compressive strength is massively reduced, Fig. 1 Schematic representation reprinted from [20] of microstructural modification occurring in the wood pellet during devolatilization based on a micro-computed tomography image analysis resulting in their fragmentation under minimal load. In this work, the described changes of wood pellets during devolatilization are modeled to deepen the knowledge about the pellets' behavior. The goal is to reproduce their structural changes and mechanical behavior. A representative model for these changes allows predicting the pellets' behavior under different conditions and facilitates a more accurate model of the entire reactor.
For modeling the pellet, the bonded particle method (BPM) was chosen. The BPM has been developed as an extension of the discrete element method (DEM) for modeling granular materials [30]. Therefore, the BPM is also referred to as bonded DEM or DEM-BPM. It represents an object as a set of primary particles connected via solid or liquid bonds. Each bond connects two primary particles and transmits forces and moments depending on the relative motion of the connected particles. All forces and moments acting on particles are integrated using Newton's law of motion leading to the particles' translational and rotational motion. In each time step, breakage criteria are evaluated for the bonds. If one of the criteria is met, the bond is considered broken and is removed from the simulation domain. The BPM has successfully been used for modeling wood pellets [14], and the discrete nature of the approach allows the simple handling of internal discontinuities like pores and reproducing complex failure behavior as seen for devolatilized wood pellets.
The BPM is foremost a mechanical modeling approach. However, recent studies have shown the capability of the DEM and the BPM to study phenomena including heat transfer [28,39,40]. The underlying DEM is often used for simulating heat transfer phenomena including also thermochemical conversion where it is typically combined with computational fluid dynamics (CFD) [27,33]. Furthermore, in such CFD-DEM investigations, the bonded DEM has been successfully implemented to model non-spherical shapes when simulating wooden particles conversion in fluidized bed reactors [22]. In the work at hand, we are now extending the BPM to enable modeling the thermochemical conversion of wood pellets, mimic the caused structural changes and predict the resulting mechanical behavior. As solver, the open-source framework MUSEN [10] was used due to its customizability.

Methods
The motivation for developing the model is to improve the capability to predict and increase the understanding of the biomass pellets' behavior inside a fluidized bed reactor. Therefore, the goal is to develop a comprehensive microscale model which will consider various phenomena occurring during devolatilization to reproduce the microstructural changes and to predict the resulting mechanical behavior.
As shown by [20], during devolatilization, the pellet is losing roughly 80% of its mass causing massive enlargement of internal pores and the shrinkage of the structure. The main changes found are shown in Fig. 1. Mainly, largescale pores (dimension: 150-500 µm) arise in areas where pores were visible in the raw state and build up into internal cavities (dimension: 0.1-2 µm). These pores in the raw state result from manufacturing and are located between different wooden chips. After devolatilization, the wood chips are covered by small pores (dimension > 150 mm), which are connected to the global pore network. Images showing the structural changes can be found in [20] and in Sect. 3.3.
Despite the large changes, the cylindrical shape of the pellet remains after devolatilization, but the strength is heavily reduced. Under compressive load, a complex failure behavior with high fluctuation in the measured force-displacement curve was found. The behavior was observed even for low strains and is likely caused by pores collapse and resulting local dislocations.
To reproduce all these phenomena, a new formulation of bonded particle method is proposed and applied for the biomass pellet. The entire model will be described in the following sections.

Bonded particle model
The BPM is a mathematical model which is applied to represent real objects and to mimic their behavior. In the following, it is briefly introduced. Further information about the method can be found in [11,12,29].
In general, a bonded particle model can be subdivided into three parts: -structural model: positions and radii of primary particles, connection between primary particles, bonds radii, etc.; -functional model: dependencies describing forces and moments acting on primary particles caused by the particle-particle or particle-wall interactions as well as due to the existence of bonds; -model parameters: specific parameters of the functional models to describe plasticity, softening, or other effects.
The structural model can either represent the real physical microstructure consisting of spherical particles and cylindrical bonds [11,37] or it can be a mechanical approximation [12,14,21]. If the real structure cannot be reproduced exactly, a coarsening approach is typically applied and some kind of homogeneity is assumed. In the work at hand, a hybrid approach is chosen which will be explained in the next section (Sect. 2.2). For the functional model, mainly three different interactions need to be modeled: particle-wall interactions, particleparticle interaction, and interaction via bonds. Depending on the desired behavior, bonds can superimpose particleparticle interactions or replace them. For modeling devolatilized pellets in this work, the superposition of bonds and particle-particle interactions is used to ensure continuity in case of bond breakage.
All functional models typically have several model parameters, which need to be set. Especially in the case of a coarsening approach, they are determined by calibrating the bonded particle model with experimental data. An in-depth explanation of how the calibration can be efficiently realized is shown by [19]. The proposed linearization methods for a faster calibration were also applied in the work at hand.

Basic functional model
For all simulations in this work, particle-particle and particle-wall interactions are modeled with the well-known Hertz-Mindlin contact model based on the work by [25] and [41]. The bond forces are modeled by a linear elastic model proposed by [30].
The movement of particles resulting from all forces is described by Newton's or rather Euler's laws of motion. Time integration of these equations is achieved via the leapfrog algorithm, which is an explicit integration scheme [10].
Based on the explicit time integration, increments in the bond forces in normal ΔF n and tangential ΔF t directions as well as torsional ΔM n and bending ΔM t moments are calculated in each time step as: where L b and A b are the bond's length in relaxed state and its cross-sectional area. E b and ν b are the Young's modulus and Poisson ratio of the bond. I and J are the moments of inertia. v r and w r are the translational and rotational relative velocities, Δt is the simulation time step. The breakage criteria for all bonds are: where σ n and τ n are the bond's normal and tangential stresses. σ b,max and τ b,max are the corresponding strengths. F n is the signed absolute value of the normal force, where the sign is defined to be negative for compression. R b is the bonds radius. If the condition (5) is fulfilled, the corresponding bond is considered to be broken and removed from the simulation domain.

Structural model
The goal of the structural model is to reproduce the major changes in the structure of the pellet during devolatilization. Thus, large-scale pores and their formation must be resolved by the structural model, i.e., by the particles and bonds.
However, the small matter pores in the devolatilized state are too small and widespread to be truly resolved by primary particles. Instead, the matter pores are considered mainly in the functional models. Therefore, the solid matter is assumed to be homogeneous, since pores in the biomass matter are smaller than the chosen scale of the model. Furthermore, the wood chips are not represented individually but the solid areas are represented by randomly distributed homogeneous primary particles. This reduces the computational complexity heavily and since the behavior of individual chips is not analyzed, the homogenization should give sufficiently accurate results. However, the pores between chips in the raw state have to be considered, since they are suspected to act as origin for the formation of larger cavities. Pores in the structural model are represented by the absence of particles and bonds in the specific area. The actual pore formation and enlargement during devolatilization will be emulated by bond breakage due to shrinkage. The bond breakage is used as an approximation to represent the loss of connectivity due to the large mass reduction. In the model, only primary particles are considered to have a mass and the conversion causes them to shrink. Bonds shrink in accordance to the particles they connect. The bond shrinkage causes local stress resulting in the failure of some bonds. Therefore, failures will mainly happen in areas with a low amount of bonds and fast shrinkage. By this means, an enlargement of the raw pores as well as the generation of new pores is achieved. The approach will allow us to mimic real-world behavior without having to model the conversion on the micro-scale.
For generating the structural model at raw conditions, we will use the micro-computed tomography (microCT) data from [20]. Moreover, the images at devolatilized state will be used to evaluate if the model is capable of mimicking the true structural changes.

Conversion of BPM structural model to voxel representation
To allow the comparison of the structural model with the microCT images, an equivalent voxel representation of the structural model is developed. For generating the voxel representation, the same voxel size of 12 µm as in the microCT was used. The entire simulation domain is represented by voxels which can either be empty or occupied by bonds or particles. The resulting three-dimensional image has roughly 550 · 550 · 1100 voxels, i.e., pixels. The voxelization of a BPM structure is schematically shown for two dimensions in Fig. 2.
To compute the voxel representation, simple conditions are evaluated for every particle and bond. They are shortly explained in the following. For each primary particle, the closest voxel to its center is used as the center in the voxel representation named x part . Each voxel x voxel fulfilling the condition is assumed to be part of the particle, where r part is the corresponding particle radius. For a bond connecting the primary particles x part [, 1] and x part [, 2], all voxels fulfilling the following conditions are assumed to be part of the bond: x voxel − x part [, 1] > r part,1 (10) The first two conditions ensure that the voxel is between the centers of the two particles and the third condition ensures that it is inside the cylindrical area of the bond. The last two conditions exclude the voxels occupied by the primary particles themselves. The voxel representation allows a straightforward identification of areas with no bonds and particles. Through these areas, no forces are transferred and no heat is conducted. Therefore, these void areas are used to represent pores in the structural model. The voxelized structural model will also be used to approximate the external surface of the structural model which is described later in Sect. 2.5.2.

Generating the structural model
The generation of the bonded particle structure consists of several steps which are shown Fig. 3 and described below: Fig. 3 Schematic visualization of the procedure to generate the bonded particle structure of a wood pellet in raw conditions randomly placed. To reduce the overlap between particles, a force-biased algorithm is applied (see e.g., [12]). All primary particles are generated with the same radius of 130 µm, and a porosity for filling the pellet with spheres of 0.36 was set. 2. Delete particles overlapping with pores in voxel representation: This step's purpose is to reproduce the holes found in the microCT data of raw pellets. A random approach would also be possible, but we are using the available data here as a template. First, a voxel representation is generated for all created primary particles. This voxel representation is compared to the microCT image of a wood pellet. The microCT image is binarized as described in [20]. If more than 90% of a primary particle corresponding voxels are occupied by pores, the particle is deleted. 3. Generate bonds between neighboring primary particles: Between every pair of primary particles where the particles are less than 130µm apart from each other, a bond is generated. This results in a bond to particle ratio of about 10 : 1 ensuring a homogeneous reaction under load in all directions. 4. Reduce bonds diameter to avoid overlap with pores: All bonds and particles are transferred again into a voxel representation and compared to the binarized microCT image. The goal of this step is to generate the defects in structure from the initial pore structure. Every bond's diameter is reduced until less than 5% of its voxels are occupied by a pore. If the final bond's diameter is less or equal to two voxels, it is deleted. 5. Adjust bond diameter for a consistent cross-sectional area: To ensure that the cross-sectional area of the entire pellet is approximately equal to the internal area used for the bond models, it is necessary to adjust bonds diameters. Otherwise, the number of bonds would for example influence the thermal conductivity of the pellet.
The cause for this problem is that material properties defined for a continuum are used in a discrete model. For more detailed information about the problem of heat fluxes in the DEM, see e.g., [39] or [28]. In the case of the BPM, quasi-randomly distributed cylinders are used to describe fluxes. For a perfectly uniform flux, each point in space should be crossed by exactly three orthogonal bonds. This is of course not possible in the model, as the bonds are cylindrical objects connecting particles. Therefore, the number of bonds and their orientation crossing each point vary (see Fig. 2). Thus, the goal is to achieve uniformity on average inside the pellet. For evaluating uniformity, the voxel representation is used as an approximation. For each bond, the voxels occupied by that bond are computed. Due to their orientation, the bonds do not transfer fluxes equally in all directions. Therefore, three different values for the different coordinate axes are computed. Each value represents the percentage that the bond would transfer in that respective direction. By summing the values over all bonds for each voxel, global values are obtained. The average of these global values for all voxels inside the pellet which are not particles should be approximately one for each direction. Values below or above one would mean that the model pellet transfers less or more heat than the thermal conductivity implies. To achieve nearly uniform transport, the bond diameter is changed until the average values in each direction are approximately one. The most important effect of this step is to achieve consistency between different models if the amount of bonds changes. Otherwise, a change of the bonds would require to re-calibrate the thermal material properties.
A resulting structure will be shown and described later in Sect. 3.3.

Extension of the BPM with a pyrolysis model
The general idea for modeling the thermochemical conversion follows the concept of the bonded particle method: By applying a simple model for the individual primary particles and bonds, complex behavior is achieved due to the interaction of the network of particles and bonds. Following this idea, each primary particle has its unique, uniform temperature. Although the individual primary particles are considered isothermal, the combined structure allows for a nonuniform temperature of the pellet. Each primary particle consists of three species: moisture (M), biomass (B), and char (C). The mass of a primary particle m part is defined by summing the mass of each individual species m i The mass fraction of a species is, thus, defined by The char content X C is initially zero and increases during devolatilization until X C = 1. The three reactions considered are: Biomass B → Primary gases G & Tars T Such a conversion model is called a single component model in literature because biomass is considered as one component. Models considering the biomass composition of lignin, cellulose, and hemicellulose are typically more accurate but the computational costs are also much higher [9]. In our model, only bonds are transferring heat. Bonds do not have a temperature and their thermal conductivity is derived from averaging the thermal conductivities of the connected primary particles. Furthermore, no bound water diffusion or other mass diffusion is considered.

Heat transfer
External heat transferQ E xt to a primary particle is modeled by a mean convective heat transfer coefficient (HTC),h conv , and an approximate radiative contribution,h rad , witḣ T part is the particles temperature and T env is the temperature of the environment, in this case, that is the bed temperature. A E xt is the primary particles external surface area. Only for particles at the boundary of the structure, the value A E xt is bigger than zero. The calculation of this area is explained later together with the definitions of the heat transfer coefficients in Sect. 2.5. The heat flux between two primary particles, having temperatures T part,1 and T part,2 and connected by a bond with cross-sectional area A b , is calculated by a discretized Fourier's law written as: x part,1 and x part,2 are the coordinates of the particles center and k mean is the mean thermal conductivity defined by the individual thermal conductivities with 0.5 k part,1 + k part,2 .

Drying and devolatilization
The mass of a primary particles changes due to two reactions: drying and devolatilization. Similar to the work by [44] or [38], drying, i.e., moisture vaporization, is assumed to be only limited by the availability of heat. Therefore, the mass flow of generated vaporṁ M→W is calculated bẏ T vap is the temperature at which moisture vaporization starts, assumed to be 373K. c M is the specific heat capacity of moisture at constant pressure and h M the heat of moisture vaporization. Devolatilization is modeled by a typical Arrhenius-based correlation (see e.g., [9]). The conversion from biomass B to char C is calculated witḣ (20) and the conversion from biomass to primary gases G and tars T is calculated witḣ The correlation factors, i.e., the frequency factors A X →X and activation energies E X →X , and all other used parameters are given in the next section (Sect. 2.4). R is the specific gas constant, which is rounded for all calculations to 8.314 J K −1 mol −1 .

Conservation of mass
As we are focusing on the pyrolysis and its impact on the solid phase, char conversion is not considered and tars, primary gases, as well as vapor are not monitored after leaving the primary particle. Therefore, the conservation of mass can be written as:

Conservation of energy
Each primary particle is defined as an open system, where the primary particles enthalpy is equal to the sum of the enthalpy of the moisture, biomass and char it contains. The total enthalpy changes due the external and internal heat flow and the mass flow of vapor, tars, and gases leaving the particle. Therefore, the equation for energy conservation is: where h M , h B and h C are the specific enthalpy for each species.
By applying the product rule for derivatives, the equation can be reformulated to: The right-hand side can be simplified by, firstly, introducing the heat of vaporization as Secondly, the change in enthalpy due to devolatilization is described by the common term of the heat of devolatilization: Assuming a homogeneous temperature inside the primary particle and applying the definition of the specific heat capacity, the energy equation becomes By introducing as the particles specific heat capacity, the equation can be written as:

Material properties and models
The thermochemical conversion model requires additional material properties and model parameters to close the set of equations. All used properties and models are summarized in Table 1. Furthermore, each parameter choice is shortly explained in the following section. As described above, the primary particle specific heat capacity is calculated by interpolating the individual specific heat capacities by the mass fraction analogously to the approach shown by [15]. For the individual specific heat capacities, the values derived by [15] are used.
The kinetics of devolatilization are described by the model of [7]. Although the used type of wood and the particle size is not given in their work, the model has proven to give reasonable agreement for different conditions and is widely used as a single component model [3,9,17,38]. For modeling, the heat of devolatilization different approaches exist. For example, [17] are integrating the specific heat capacities. Often the heat of devolatilization is also neglected assuming that the endothermic generation of volatiles and the exothermic generation of char balance out (see e.g., [5]). However, depending on the conditions, the heat of devolatilized can be significant for accurately predicting the kinetics. A more detailed model for the heat of reaction in biomass pyrolysis was derived by [2]. They are using a detailed reaction scheme to get a more accurate prediction of the required heat. To get an estimate of the heat, there are also existing correlations based on the char yield for a limited range of conditions [9,32]. They require determining the char yield experimentally or by other models. [38] have shown that a linear fitting based on the temperature slightly improves the accuracy when modeling thermally thick particles but a fixed value already gives accurate results. Since we are working with a simple devolatilization model with no adaption to the type of wood, it appears reasonable to use also a very simple model for the heat required. Furthermore, no accurate data for all different conditions and pellets used in this work are available. Therefore, we are using for the required heat a proposed fixed value by [38] to get a rough approximation.
The thermal conductivity of the primary particles is calculated by weighting the conductivity of wood pellets and char based on the current mass fraction. The approach is equal to the approach for calculating the specific heat capacity and is often used in literature to model a mixture of materials (e.g., [15,24]).
The influence of the moisture on the thermal conductivity is neglected due to the low moisture content in pellets and for this application comparatively small influence of the moisture which was measured by [36]. The value for the wood pellet's thermal conductivity is based on measurements by [23]. They measured the conductivity of different homogenized wood pellets. As the density of the wood pellets they used are similar to the reference pellets in this work (i.e., the pellets from [4,20,34]), the values should give a good approximation for all cases. For the char thermal conductiv-ity, the softwood char particle conductivity measured by [16] is used.

Heat transfer
The heat transfer has a significant influence on the conversion time in a fluidized bed reactor. Predicting the heat transfer coefficient for the pellet is, therefore, crucial to accurately predict the conversion rate of the pellet. [4] compared different models for their agreement with experiments and their performance in simulations. They found that the best agreement was achieved with the model by [1], which will also be used in this work. Agarwal's model is a mechanistic model for a particle inside a hot fluidized bed. It accounts for particle convective components as well as gas convective components in the emulsion and bubble phase.
In Section 3, three different experimental conditions will be analyzed: the experiments by [4], the experiments by [34] and the experiments by [20]. For all conditions, a mean heat transfer coefficient was calculated with the Agarwal model. Since the change of the heat transfer coefficient with the pellet's temperature is comparatively small, a fixed value was assumed during simulation. Simulations for comparing with the results from [4], [34] and [20] are based on a convective heat transfer coefficienth conv = 450 W m −2 K −1 ,h conv = 550 W m −2 K −1 andh conv = 550 W m −2 K −1 , respectively. The values are also in good agreement with the correlation derived by [26] and appear to be reasonable compared to the given experimental data by [4].

Thermal radiation
Thermal radiation is not considered by the model of [1]. Radiative heat transfer has mostly a limited influence for fluidized bed devolatilization due to the high convective heat transfer [4,8,18,26]. Nevertheless, since in this work bed temperatures of up to 1173 K are considered, a small influence of thermal radiation is expected. Therefore, a simple approximation of the radiative heat flux is included. The radiative heat transfer coefficient can be calculated by: where b is the bulk emissivity, s is the pellets emissivity, T b is the bulk temperature, i.e., the bed temperature, and σ is the Stefan-Boltzmann constant equal to 5.670 × 10 −8 W m −2 K −4 [26]. According to [13], the bulk emissivity can be approximated by where B is the back-scatter fraction and e p the emissivity of the bed material particles. For circulating fluidized beds B can be approximated with 0.667 [6]. The emissivity for the bed material in the experiments, i.e., quartz sand, is 0.9 [43].
By approximating the emissivity of the biomass pellets with the emissivity of carbon, s = 0.8 [31], the bulk emissivity can be calculated to be approximate b = 0.9.
During calculation, it is assumed that thermal radiation occurs at all time. This appears a sufficiently accurate approximation since the influence of radiation is small. Furthermore, the Agarwal model [1] predicts a probability of above 90% for a comparable sphere to be in the emulsion phase.

External surface area
Convective and radiative heat transfer are both assumed to only impact the primary particles with contact to the surrounding. However, the surfaces of the spheres are not representative for the surface of the solid area they are representing. Therefore, the surface in contact with the surrounding is calculated based on the voxel representation.
For determining the object's boundary, three steps are performed on the voxelized structural model: 1. Spherical dilation of occupied voxels, 2. Flood filling starting from the domain boundary, 3. Spherical erosion of occupied voxels, 4. Identification of boundary voxels.
Due to the spherical dilation followed by later erosion, only voxels that can be touched by a sphere with the specified diameter are potentially identified as boundary voxels. We are using a sphere diameter corresponding to the minimal diameter of the bed material used by [20] or [34] (i.e., 100 µm or 8 voxels). Therefore, all areas reachable by the bed material are considered as external boundaries. Thus, pore surfaces that can be reached by the bed material are considered to directly receive heat from the fluidized bed. The assumption is that if the pore is big enough that bed material can freely enter the pore, the heat transfer inside the pore is of a similar dimension as at the outer surface due to particle convective heat transfer. To reduce the complexity, all areas are assumed to have the same heat transfer coefficient. Thus, the heat flux inside the considered pores is likely overestimated. However, since at the same time the heat flux in smaller pores is neglected, the approach is considered as a reasonable approximation. Further investigations of the internal heat transfer mechanisms and flow behavior will be needed, to facilitate more accurate modeling.
For calculating the surfaces, all boundary voxels are assigned to the nearest primary particle. Each boundary voxel is assumed to represent the same external surface equal to the squared voxel size (12 2 µm 2 ).

Shrinkage model
The shrinkage of primary particles and bonds is crucial for coupling the thermochemical conversion model to the structural and mechanical model. In principle, the mass loss can be simply coupled to a volumetric shrinkage of the primary particles. However, since the smallest pores are not resolved by the structural model, we need to assume an increase in porosity of the primary particles. Therefore, the time-dependent porosity of each primary particle φ(t) is defined by with φ 0 as initial porosity, α as porosity factor and i=M,B X i,0 − X i (t) as the relative mass loss. In this work, φ 0 is considered to be zero. The primary particle volumes and radius are calculated by For the bonds, we want that the ratio of the bond length to the radii of connected primary particles remains constant. A constant relative distance ensures the similarity of the mechanical behavior because it keeps the influence of particle-particle contacts in the mechanical model at a similar level. Furthermore, the bond diameter should shrink comparably to the particles radii as well. This ensures a consistent shrinkage of the cross-sectional area and volume. Therefore, the bond length l bond and the bond diameter d bond are constantly updated by the following equations: The reduction in the bond length will induce forces, leading to the movement of primary particles. Globally, this will cause the entire pellet to shrink. Hence, the shrinkage factor α has to be adjusted to match the real-world shrinkage of the pellet. The local bond forces induced by shrinkage can exceed the bond's strength, especially due to the hightemperature gradients and the fast reaction. This mechanism of the resulting bond failure is used to model the formation of macro-scale pores.

Results and discussion
In the following section, the capabilities and limitations of the developed method are presented by comparing simulations to experimental results. In the first step, the calibration of unknown model parameters is explained. Secondly, the thermochemical conversion is validated against different experimental data by [34] and [4]. Thirdly, the structural changes and the three-dimensional transient behavior are analyzed based on the results shown by [20]. In the last step, compression tests were simulated with the predicted structures and also compared to experimental data from [20].
All simulations were carried out using the before described method which is based on the implementation of MUSEN in Version v1.71.2 [10]. The main part of the computations were performed on graphic cards using CUDA (versions 11.2-11.5).

Calibration
Since the BPM cannot exactly represent the material microstructure an additional calibration of the model parameter is needed. For the proposed model three different adjustment steps during simulation need to be considered.
Firstly, there are the parameters that define the kinetics of thermochemical conversion. These parameters do not require any calibration as they can be obtained directly from experimental data (Sect. 1). For the thermal conductivity, this is only true due to the developed approach for generating the structural model (Sect. 2.2).
Secondly and thirdly, there are the mechanical model parameters that play an important role during the conversion and during the compression tests. Ideally, these parameters would be identical and change continually based on a single model for the transient behavior of the pellet's properties during conversion at elevated temperatures. However, no experimental data for such a complex model exists to our knowledge. Instead, two independent sets of constant parameters are used during the conversion and the subsequent compression testing. Since the goal of this work is to evaluate the general suitability of such a modeling approach, this approximation was considered sufficient. If the model proves to be suitable and more detailed data about the conversion becomes available, a more sophisticated approach could be implemented.
For achieving an agreement between the structural model and the experimentally obtained pore structure, the strength of the bonds and the stiffness were calibrated. The used values for the Young's modulus of the bonds and particles as well as the breakage strength are E = 160 MPa and σ b,max = τ b,max = 30 MPa, respectively. Reducing the breakage strength results in higher porosity and can potentially lead to the fragmentation of the pellet. A higher stiffness slightly reduces the structure's shrinkage and increases the porosity because of the higher stresses (assuming a constant strength). The bonds were assumed to be soft during conversion; therefore, the relatively low stiffness was chosen for mimicking the structural changes. In the shrinkage model, a porosity factor of α = 0.375 was used, which was calibrated against the real-world pellets shrinkage. The physical meaning of these parameters is, however, limited, as they would describe a kind of average behavior of the pellet during the conversion from wood into char.
After determining the parameters for the mechanics during conversion, the parameters had to be adjusted to represent the pellet in the devolatilized state. Due to the large scattering of the properties of devolatilized pellets, the goal was to match the general behavior and not the reaction of individual pellets. Since the raw material properties of wood pellets just like the pore structure after devolatilization fluctuate between different pellets, we cannot expect a single structural model and functional model parameter set to represent the entirety of experiments. Therefore, multiple sets of parameters were used for simulating the compression testing. They are discussed later in Sect. 3.4.
The accuracy of the model predictions based on the described calibration is evaluated in the following sections.

Thermochemical conversion
Predicting the thermochemical conversion is crucial for all subsequent steps. Therefore, multiple simulations were carried out under different conditions to compare the results to the experimental conditions from [34] and [4]. For the comparison, the predicted mass change and temperature are used. A more detailed insight into the internal spatial distribution of the simulation results is given in Sect. 3.3. In the experiments by [4] and [34], multiple pellets of different lengths were used for determining the mass loss to account for fluctuations. [4] found a limited influence of pellet length on the results. Therefore, a single pellet with a length of 15 mm was used in all simulations. Based on the microCT images from [20], multiple generic pellets were generated. Afterward, the thermochemical conversion of these pellets was simulated. Figure 4 shows the comparison between experimental data and the model predictions of the normalized remaining mass and temperature over time inside fluidized bed reactor at 673 K. The comparison shows a good agreement of the model with the experimental data. The expected slower mass loss during initial heat up, the slow-down of the temperature increase during drying, the subsequent increase in the devolatilization rate and the decay at the end are all reproduced. However, the heat up appears to be slightly slower in the simulation. This could be caused by the homogeneous Fig. 4 Comparison of the model predictions with the measured relative remaining mass and the pellet core temperature from different experiments by Reschmeier [34] and Berg et al. [4] in fluidized bed reactors at 673 K. Shaded areas represent the standard deviation as given by the authors Fig. 5 Comparison of the model predictions with the measured relative remaining mass from different experiments by Reschmeier [34] and Berg et al. [4] in fluidized bed reactors at different bed temperatures. Shaded areas represent the standard deviation as given by the authors initialization of the model, while in the experiments, the initial time point is likely determined by offsetting. Figure 5 shows a comparison between experimental data and the model predictions of the relative remaining mass of the pellet for different bed temperatures. The figures show that the scaling of the model with temperatures in general agrees well with the experimental data. However, the remaining mass at the end of the devolatilization, i.e., the char content, is not correctly predicted. The model predicts rather similar values of 0.19, 0.18, and 0.17 for bed temperatures of 873 K, 923 K, and 1123 K. In contrast, in the experiments, the values are approximately 0.18, 0.13, and 0.11. This behavior is expected since an accurate prediction of the char content requires an adjustment of the parameters of the single component conversion model or a more sophisticated model.
The computation time for only evaluating the thermochemical conversion within the BPM framework is in the Fig. 6 Image comparison of the generated structural model and its voxelized representation with the visualized linear attenuation coefficient measured via microCT (measurements from [20]) order of minutes for simulating 60 seconds with 0.3 million particles and 3 million bonds (on an NVIDIA GeForce GTX 1080 Ti). However, when coupled to the mechanical model, the computation time is roughly 100 times higher since lower time steps are needed for the mechanical model, and the computational cost of the mechanical model is much higher.
In summary, we can conclude that for the goals of this work, the conversion model predictions are more than sufficiently accurate. For mimicking the structural changes, a basic approximation of the processes is needed to model the internal shrinkage appropriately. In fact, the very good agreement of the predicted values with the experimental data shows that the proposed mesh-free BPM might even be a viable option when more accurate modeling of the conversion is needed. In principle, it appears that the BPM allows using conversion models developed for the conversion of small particles for predicting the behavior of pellets manufactured from those small particles.

Structural changes and transient behavior
Based on the thermochemical conversion model, the structural changes are emulated. The goal is to reproduce the main phenomena found by [20]. First, a pellet structure is generated based on the microCT data. Note that for the experiments, both ends of the imaged pellet were sanded, and one end was sanded obliquely to allow identifying the orientation of the pellet after devolatilization. The corresponding experimental conditions of a fluidized bed reactor at 1173 K were used for the simulations. The simulations were carried out with a time step of Δt = 6 × 10 −7 s. Simulating one second takes around 1.5 h of computation on an NVIDIA Quadro GV100 with the used 0.13 million particles and 1.3 million bonds.

Initial structure and conversion
In Fig. 6, the initial voxelized structural model is compared to the microCT image data available from [20]. When comparing the microCT image and the derived structural model, the generated larger pores between different chips are clearly visible. The chosen bond to particle ratio of 10:1 and the chosen bond diameter result in mainly filled areas beside the voids. However, due to the random placement of the homogeneous spheres, some smaller randomly generated voids also exist. When regenerating the model, the random placement of spheres will lead to small deviations in the displayed structure.
During simulation, the temperature of the outer particles rises rapidly causing the expected internal temperature gradients as visualized in Fig. 7. Connected particles can have temperature differences above 800 K resulting in local temperature gradients above 4 × 6 K m −1 . The temperature difference between core and surface of the pellet increases the first few seconds until about 750 K and then starts to decay.
Due to the fast heat transfer and temperature gradients, also the process of the conversion varies inside the pellet. Figure 7 shows that the gradient of the char content is very steep and mainly two areas exist: a shrinking unconverted core and an extending converted shell. In the transition zone between those two, most of the bond breakages occur. Around 85% of the bond breakages happen within 0.5 s after a significant change in their char content (ΔX C > 0.1)). The bond breakages are caused by local stresses due to shrinkage and result in the targeted formation of pores. With increasing time and conversion rate, more and more pores arise. In the end, around 40% of the bonds are broken for the shown simulation. The increasing porosity has also a slight impact on the conversion since larger pores are considered open to the fluidized bed. Therefore, the heat flux is increased close to these pores (e.g., visible in the top left for X C ≈ 0.2).

Final structure
In Fig. 8, the voxelized structural model is compared to the microCT image after full devolatilization. The voxelized structural model shows that pores have successfully formed during simulation. Since the model only emulates the realworld behavior, the pores do not match one to one to the real-word pellet. A full match is also not likely, because the predicted pores depend on the initial structural model, where the exact particle position is quasi-random. Therefore, also the occurrence of breakage and, thus, the pore formation slightly vary between different structures. With corresponding experimental data, this fact could even be used to reproduce the measured stochastic of the fragmentation inside the fluidized bed reactor. Although the pores are not exactly reproduced, the general trends and characteristics of their formation are: Large pores are formed and build an entire network including cavities in the inner part. Micropores are, by design, not reproduced and were modeled by an increase in porosity of the primary particles.
Although the formation of a large-scale pore network with macro-scale pores is visible, the formation of cavities is less pronounced. This results in an about 30% lower overall porosity (including the porosity of the particles) compared to the porosity obtained from the microCT image. Moreover, the predicted pore inlets are bigger than in the measured pellet. The microCT images in Fig. 8 show that the inlets are rather small and open up toward the core. It was suspected that the pore shape could be influenced by an over-pressure inside the pores [20]. Since the model is currently not resolving the gaseous phase, the pressure from the gaseous release is not taken into account. However, an influence of the cooling stage by the thermal contraction on the imaged pore shape can also not be ruled out because the pellet had to be cooled down before the microCT imaging. Nevertheless, the general trends of the pore formation are reproduced.

Summary of structural changes
In summary, the developed approach for reproducing the structural changes with the BPM shows a good agreement with the structural changes revealed by microCT. The formation of a large pore network is clearly visible in the model and the similarity of the structure appears high enough to allow accurate mechanical modeling. However, cavity formation is underestimated resulting in a smaller porosity. Not all physical phenomena are considered and more experiments are needed to better understand the pore formation and improve the accuracy of the model.

Mechanical behavior under compressive load
In the final step, we compare the mechanical behavior of the predicted structural model compared to experimental data. Compression tests in the radial direction of the devolatilized model were simulated and compared to the measured forcedisplacement curves by [20].
For the simulations, one structural model of a 1 mm long pellet was generated, and, subsequently, the devolatilization was simulated. Afterward, compression was simulated with different parameters. All simulations were carried out with a time step below 2 × 10 −7 s. Simulating one second takes around 3 h of computation on an NVIDIA Quadro GV100 for the 0.14 million particles and 1.1 million bonds used in this section. The speed of the compression stamp was set equal to the experimental tests (i.e., 0.1 mm s −1 ). Figure 9 shows the measured force on the stamp over the displacement both in the experiments and simulations. Three different parameter sets are presented since no single set of parameters can reproduce the diverse behavior of the pellets under compression. Furthermore, the different sets visualize the influence of different parameters on the pellet's behavior.
In Fig. 9, one can clearly see that the high-frequency fluctuations in the measured force are reproduced by the model. In the model, they are caused by the failure of bonds followed by local dislocations inside the structure. Breakages of bonds are occurring from the beginning on and are dominating the behavior of the model. The bond breakage depends of course on the stiffness and the strength. For the shown Fig. 9 Comparison between experimental and simulation data of force-displacement curves during radial compression. The experimental data originates from [20] and five representative curves of identically treated and tested wood pellets are shown. The simulation data shows the results from one raw structure with three different material parameters Fig. 10 Monochrome slices through structural model during compression in relation to the force-displacement curve and the number of broken bonds. In the top, overlapping radial slices before (red) and after (black) a sudden force drop are shown. At the bottom, the evolution of the structure is shown by axial slices. Case 1 of Fig. 9 is used as example parameter sets, the typical order of magnitude of the total bonds breaking per one percent strain (displacement divided by raw diameter) is 5%. The constant failure of bonds is also limiting the slope of the force-displacement curve. Therefore, the bond breakage strength has a very high influence on the apparent stiffness as shown by the different parameter sets in Fig. 9b).
As Fig. 9a) also illustrates, there does not exist a clear breakage point for devolatilized pellets (see also [20]). With increasing strain, higher force drops might occur in some cases. In other cases, the slope of the force-displacement reduces strongly and the fluctuations become much stronger. These tipping points can be used to define a breakage point. Despite the diversity, the failure of the simulation model appears to be quite similar to the measured behavior. The same structural model can have, depending on the model parameter, a slow plastic failure, a sudden larger force drop, or a kink in the perceived stiffness followed by multiple force drops.
To get a better understanding of the model's reaction and failure under load, Fig. 10 visualizes the change in the voxelized structural model during compression for the first case in Fig. 9. The figure reveals how the breakage of bonds leads to the fragmentation of the model. The comparison of the structure before and after a rapid reaction force drop in the top shows how an agglomerate breaks apart from the structure causing the rapid change in the measured force. In both shown cases only around 0.2% of the total bonds are failing, but their failure is causing a fragment to break apart from the pellet.
In general, the number of broken bonds raises steadily as Fig. 10 visualizes. With increasing strain, more bonds fail and agglomerates, start to move horizontally. This leads to large void zones and the pellet loses its cylindrical shape and structural integrity. Although the pellet structure eventually fails, as also the rapidly increasing number of broken bonds indicates, there is still a reaction force due to contact forces between agglomerates.
The behavior after breakage is highly influenced by particle-particle interactions. In the beginning, only around 5% of particles have significant particle-particle interaction forces increasing to around 15% at the end. Particles with a significant interaction force have on average triple the number of broken bonds. Hence, the particle-particle forces are slightly compensating for the increasing amount of breaking bonds. This results in the reaction force after pore collapses and the final failure of the structure, which was also observed in the experimental data.
The mechanical behavior and the breakage shown in Fig. 9 and Fig. 10 matches the described behavior by [20]. Thus, the failure of bonds and the resulting local movement of agglomerates inside the structure are possible explanations for the reported behavior of devolatilized pellets. In a real-world pellet, this could correspond to the movement of devolatilized chips.
In general, the simulation results of the compression test demonstrate the significance of the pore structure for describing and understanding the mechanical behavior of devolatilized wood pellets. The performed simulations further could show how the pore network promotes attrition under load due to fragmentation. The large pores allow agglomerates to break apart from the structure when the pellet hits a wall.

Conclusions
Fluidized bed gasification of biomass pellets might be one possible approach to generate synthesis gas sustainably and mitigate climate change. By modeling the devolatilization of a single pellet, we want to deepen the understanding of the processes to facilitate a better reactor design and enable more sophisticated simulations for an improved upscaling of the technology.
In this work, we explored the capabilities of a bonded particle method to simulate a single pellet's behavior during devolatilization. For this purpose, the bonded particle method was extended with models to predict the thermochemical conversion of a pellet. The comparison with experimental data from literature showed that this approach is capable of predicting the conversion characteristics, like mass loss and core temperature, very accurately. Simple conversion models based on the particles' behavior were employed, and by the interaction of many interconnected particles, the complex behavior of the pellet was reproduced. Furthermore, by implementing a shrinkage model for the particles and bonds, the structural changes of a single pellet during devolatilization were successfully emulated with the BPM. The good agreement of the results gives an indication that the decomposition of real-world bonds could be a driving force for the pore generation during devolatilization. However, overpressure in the pores and the diffusion limitation through the pore network are phenomena that were neglected and would need to be included for an accurate analysis.
The predicted structural changes of the developed model were shown to be accurate enough to reproduce the pellet's mechanical behavior with standard functional models (elastic bond and Hertz-Mindlin contact model). The complex reaction force and failure behavior of the devolatilized pellets under radial compression were successfully reproduced. The simulation showed that fragmentation and the following movement of agglomerates can explain the reported behavior under load. Therefore, the simulations could confirm the assumption that the mechanical behavior of wood pellets devolatilized in fluidized bed is dominated by the forming pore structure.
The derived multi-physics model showed that the BPM is capable of reproducing the mechanical behavior of devolatilized wood pellets based on a rough approximation of the raw structure of the pellets. Therefore, it enables studying in future different load scenarios on a pellet in a reactor and, to a small extent, the effects of different reaction conditions on the mechanical behavior.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.