Heat transfer in fluidized beds with heat production

The heat transfer behavior of particles and gas in an olefin polymerization fluidized bed was numerically analyzed using an in-house developed 3-D computational fluid dynamics discrete element model (CFDDEM). First the implementation of the model was verified by comparing simulation results with analytical results. A constant volumetric heat production rate was implemented in the particle energy equation to mimic the heat production due to the polymerization reaction. It was found that the probability density function (PDF) of the particle temperature becomes more homogeneous with increasing superficial gas velocity. Furthermore, instantaneous snapshots of the thermal driving force (the difference between the single particle temperature and bed-average gas temperature, Tp-) for different heat production rates provide detailed insight in the particle temperature distribution inside the fluidized bed. The timeand bed-averaged particle convective heat transfer coefficient, which was calculated by Gunn’s correlation, was found to be independent of the superficial gas velocity. This is explained by the fact that the relative velocity of gas and particles in the emulsion phase, where most of the particles and gas interact, is hardly influenced when increasing the gas superficial velocity. From the spatial distribution of Nusselt number, it becomes apparent that the high heat transfer regions are found in the wake of rising bubbles, whereas low heat transfer rates are found in the clouds of the bubbles. *This chapter is based on: Li, van Sint Annaland, Kuipers, and Deen (2016a)


Heat transfer in fluidized beds with heat production
A large number of olefin polymerization processes are operated with heterogeneous Ziegler-Natta catalysts in gas-solid fluidized bed reactors. Understanding of the fundamental phenomena encountered in gas-phase polymerization reactors is essential for a proper design and optimal operation. Due to the exothermic nature of this reaction, a large amount of heat needs to be removed from the reactor efficiently to prevent hot-spots formation. In order to reveal the mechanism of formation of hot spots in fluidized beds, the following studies have been described in this thesis: The heat transfer behavior between particles and gas in an olefin polymerization fluidized bed was numerically analyzed using an in-house developed 3-D computational fluid dynamics discrete element model (CFD-DEM). A constant volumetric heat production rate was implemented in the particle energy equation to mimic the heat production due to the polymerization reaction. The effect of various operating conditions on the heat transfer has been systematically studied, viz. gas superficial velocity, operating pressure, particle size distribution, heat source term and the geometry of the bed.
It was found that the probability density function (PDF) of the particle temperature becomes more homogeneous with increasing superficial gas velocity. The time-and bedaveraged particle convective heat transfer coefficient, which was calculated by Gunn s correlation, was found to be independent of the superficial gas velocity. From the spatial distribution of the Nusselt number, it becomes apparent that the high heat transfer regions are found in the wake of rising bubbles, whereas low heat transfer rates are found in the clouds of the bubbles.
The elevated pressure influences the hydrodynamics and significantly enhances the heat transfer performance of a fluidized bed. Operating the bed at an elevated pressure helps to reach a more isothermal condition in the fluidized bed, thereby reducing the opportunity for hot spots to develop. In polymerization reactors the adiabatic temperature rise is relatively high, therefore the effect of the pressure on the average temperature of the particles can be entirely attributed to the increased gas density, whereas changes in the bed hydrodynamics have an effect on the spread of the particle temperature distribution.
Simulations have been carried out for beds with Gaussian particle size distribution with xi xii Summary different standard deviations, but with the same Sauter mean diameter (d 3,2 =1.2 mm). The thermal energy equation of the particles contains a heat source term to mimic the heat production due to the exothermic process and two cases have been considered: a constant volumetric heat production rate (q v , [W/m 3 ]) and a constant heat source per particle (Q, [W]) to represent different reaction systems, viz. "normal" catalytic reactions and polymerization reactions respectively. It was found that small particles with a high heat production rate cause the formation of hot spots in the bed. Furthermore, snapshots of the simulated particle configurations reveals that these small particles have a higher chance to be found near the top of the bed and in the vicinity of the side walls of the reactor. The former is due to minor size segregation in the vertical direction, the latter is caused by preferential particle motion.
To obtain a better understanding of the particle temperature distribution in fluidized beds, a high speed infrared (IR) camera and a visual camera have been coupled to capture the hydrodynamic and thermal behavior of a pseudo-2D fluidized bed. The experimental data were subsequently used to validate an in-house developed computational fluid dynamics and discrete element model (CFD-DEM). In order to mimic the heat effect due to the exothermic polymerization reaction in the pseudo-2D fluidized bed reactor, a model system was used. In this model system, heat is released in zeolite 13X particles (1.8-2.0 mm) due the adsorption of CO 2 . Characteristics of the adsorption kinetics, isotherm and reaction enthalpy have been achieved by performing Thermogravimetric Analysis (TGA) and Simultaneous Thermal Analysis (STA). The combined technique provides insightful information on the particle temperature distribution for different CO 2 concentrations. Furthermore, the comparison of the spatial and temporal distributions of the particle temperature measured in the fluidized bed with the simulation results of CFD-DEM provides qualitative and quantitative validation of the CFD-DEM, in particular concerning the thermal aspects.

Project background
Polyolefins are widely applied in a variety of applications in daily life, such as plastic bottles, grocery bags, home appliances, engineering plastics, automotive parts, medical applications, and prosthetic implants. Polyolefins are made up of monomers composed of only carbon and hydrogen atoms. By simply altering the amount of supplied ethylene, propylene, and higher α-olefins which are bound in the polymer chain, polyolefins with entirely new properties can be produced. Industrial and academic research on olefin polymerization catalysis has been very intense since the original discoveries of Ziegler Natta catalysts, which were produced in the early 1950s by Ziegler and Natta (Cossee, 1964). Polyolefins are made in a variety of reactors, such as autoclave reactors, single-and double-loop reactors, tubular reactors, and fluidized-bed reactors; these processes may be run in solution, slurry or gas phase. This thesis will only focus on the gas-phase olefin polymerization, which is carried out in fluidized beds with heterogeneous catalysts.
Operating heterogeneous catalyzed processes in two-phase or three-phase reactors has an increased level of complexity compared to single phase reactors, due to the fact that these systems are dealing with inter-and intra-particle mass and heat transfer resistances during the polymerization process. The transport resistances may affect the overall polymerization rate and polymer microstructure. If significant, mass and energy transport limitations create nonuniform polymerization conditions within the catalyst particles that lead to non-uniform polymer microstructures. Due to the highly exothermic nature of olefin polymerization a large amount of heat needs to be removed from the particles to prevent the formation of temperature gradients inside the particles. In extreme cases, the accumulation of reaction heat can lead to particle melting, sintering, adherence and agglomeration, which may ultimately cause defluidization. These issues are especially significant in gasphase polymerizations.
In order to prevent hot-spot formation, a full understanding of the heat transfer phenomena of fluidized suspensions is important and hereby essential for optimizing the performance of fluidized bed reactors. The complexity of describing the heat transfer performance in the reactor is mainly associated with the multi-scale phenomena encountered in the process, as shown in Fig. 1.1. The reactions are taking place on the catalyst active sites on the molecular scale, which can be considered as the micro-scale (micro-second). The inter-and intra-particle mass and heat transfer process happen on the scale of particles and are influenced by the scale of reactors, therefore the time scale is described in seconds. From a global scope of the process, it takes 1.5-2.0 hours for a fresh catalyst with a size of 20-25 micron to grow to the size of 1-2 mm spherical polymer product. The mechanisms of hot-spot formation can be related to any of the phenomena at the different time and length scales mentioned above. Therefore, all parameters influencing these process should be investigated.

Motivation of the PhD thesis
The prime motivation of the work described in this thesis is as follows: to investigate the heat transfer phenomena in an olefin polymerization fluidized bed reactor with a high heat production rate, and to explore the mechanism of hot-spot formation.

Fluidized beds
Large amounts of granular material can behave as a liquid when being passed through by sufficiently large amounts of a fluid (either liquid or gas). This phenomena can be explained as follows: when a gas flow is introduced through the bottom of a bed of solid particles, it moves upwards through the bed via the interstitial spaces between the particles. Mainly due to drag force acting on the surface of the particles, they collectively behave as a fluid, when the fluid velocity exceeds a critical value (minimum fluidization velocity). The process is called fluidization and the reactor that is used to contain this type of flow is called a fluidized bed. Intense mixing is one of the most important characteristics of fluidized beds. It is also famous for its excellent heat and mass transfer between the phases. The heterogeneous olefin polymerization is one example of the many applications of fluidized beds. The liquid alike behaviour of catalysts or polymer particles guarantee the continuous operation. Fresh catalyst particles are continuously injected into the bed while the products can be collected at the bottom of the bed and removed. Excess gas can remove a large amount of reaction heat. Up to 90% of the monomer gas (ethylene or propylene) is added in excess and cooled down by an external cooling facility and recycled back into the bed. The monomer gas therefore acts as both reactant and coolant.
Another important fact that needs to be mentioned here is that the commercial olefin polymerization fluidized beds are operated at high pressure, which can be elevated up to above 20 bar. Operating at high pressure can increase the conversion rate of polymerization and meanwhile enhance the cooling capacity of the gas phase. At elevated pressure, several researches have reported that the fluidization pattern will be strongly influenced (Li and Kuipers (2002); Godlieb et al. (2008)). Pressure effects on the heat transfer in fluidized beds is therefore also one of the interesting aspects for further investigation.

Modeling strategy
To improve the fluidization efficiency and heat management of these reactors a better understanding of the gas-particle heat transfer is required. Detailed computational fluid dynamics (CFD) simulations are capable of providing such insight. Several CFD models have already been developed for fluidized bed reactors for olefin polymerization (Khan et al., 2014). However, it is not possible to develop one model that can be used to simulate the whole process covering all the time and length scales aforementioned. A multi-scale modeling strategy is widely used to solve or understand multi-scale multi-physics problems. This idea is far from being new. The constitutive equations that are used to describe the physical phenomena normally have a number of assumptions, and they are solved with discretization methods. Fully resolved simulation methods, such as Lattice-Boltzmann (LB) and direct numerical simulation (DNS), can provide accurate modeling results of transport properties of particulate systems. However, the required computational times becomes impractical for large systems. Therefore typically the number of particles and the size of the simulation domain is too small to capture the macroscopic motion of a large system. These full-resolved models are used to provide closures for larger scale models, such as the drag coefficient, inter-phase heat/mass transfer coefficient, etc. These are required by the popularly applied discrete element model (CFD-DEM) and two fluid model (TFM). The former is an Euler-Lagrange model and the latter is based on the Euler-Euler approach. In the discrete element model the particles are tracked individually and the gas phase is solved as a continuous phase on an Eulerian grid which is much lager than the particle size. In the TFM both the solids and the gas are considered as continuous phases. Hereby the computational time which is spent on sorting and solving particle-particle interactions can be saved in this continuous model. In principle the TFM can be used to simulate bench scale and even large industrial scale reactors and has provided good results that compare well with industrial data. One of the disadvantages of TFM in this application of predicting hot-spot formation is that with this continuous model it is not possible to capture any individual particle motion or particle-scale thermal behaviour. Consequently, this brings difficulties in revealing the mechanism of particle over-heating. Therefore, in this work hot-spot formation is investigated by performing simulations using a discrete model, such as DEM, which normally can simulate up to 1 million particles (at acceptable computational cost).

Experimental investigation
It is mentioned above that the extreme exothermic nature of the olefin polymerization requires excellent heat removal. Therefore, high pressure fluidized bed reactors are commonly used in industry, such as the UNIPOL ® process for producing polyethylene (PE) and polypropylene (PP). Although the gas-phase olefin polymerization fluidized beds have been operated for decades, the complex interplay of hydrodynamics, and heat and mass transfer is not yet fully understood and this still remains a challenge due to the difficulties in the description of these flows. One of the issues related to this process is the local hot-spot formation. Understanding heat transfer in fluidized beds with heat production can help to optimize the process.
Non-intrusive experimental methods such as X-ray tomography, positron emission particle tracking (PEPT) and electrical capacitance tomography (ECT) have been utilized to study the hydrodynamics of fluidized beds for several years. All of the mentioned methods require sophisticated reconstruction of the measured data. Direct visualization of fluidization can be achieved by applying high speed visual cameras to record images obtained from pseudo-2D fluidized beds. When the fluidization domain is purposely designed with a really thin depth (approximately 10-15 particle diameter), it is called a pseudo-2D fluidized bed. It is like a cross section of the 3D fluidized bed and the bubble can be directly imaged by the camera. The recorded images of particles and gas bubbles can later be combined using image processing algorithms such as digital image analysis (DIA) and particle image velocimetry (PIV) to subtract useful information such as solid volume fraction, bubble size, particle velocity and solids flux in the fluidized bed. The temperature field in the a pseudo-2D fluidized bed can be measured non-invasively using infra-red (IR) optical technique. The measurement can only be applied to a pseudo-2D domain as the IR camera can only detect signals from the surface of an object. In some sense the IR and visual camera work in a similar way, in which both cameras detect the same thing: electromagnetic radiation, depending on the range of the spectrum. Infrared radiation is limited in the range of 0.7 µm to 10 3 µm. With proper calibration of the camera, the temperature of the objects can be related to the intensity of the radiation.
Following this basic idea, an integrated IR/DIA technique has been used to investigate the heat transfer in a pseudo-2D fluidized bed. This experiment is used to validate the inhouse developed CFD-DEM. By placing a high speed visual camera and a high speed IR camera in front of a pseudo-2D fluidized bed, the temperature and position of particles are recorded simultaneously.
As aforementioned, the main goal of this work is to investigate the heat transfer in olefin polymerization fluidized beds, and the intensive heat production is one of the crucial causes of hot-spot formation. In literature relatively few experimental studies have been conducted involving the heat production in fluidized beds. It is more common to perform experimental fluid dynamics (EFD) with a "cold flow model" without any reactions taking place using inert phases. It would be very interesting to perform olefin polymerization experiments in a pseudo-2D fluidized bed, and correspondingly to use these for validating the CFD-DEM. However, the olefin polymerization process is operated at elevated pressure, meanwhile the catalysts are extremely sensitive to the ambient conditions. It is difficult to perform the actual polymerization reaction directly in a pseudo-2D fluidized bed. Therefore, an alternative system has been used in this validation, where the heat source term is obtained from the particles. An overview of the main of similarities of these two processes is presented in table 1.1.

Outline of the thesis
The main content of the thesis is summarized in Fig.1.2. In the first two chapters, a constant volumetric heat production rate was implemented in the particle energy equation to mimic the heat production of the polymerization reaction, which was estimated from literature with an averaged value (Choi and Ray, 1985). Verification of the model implementation and a general understanding of the heat transfer process in the fluidized bed are presented in Chapter 2 and 3. In these two chapters, the effects of operating conditions (superficial gas velocity and pressure) are investigated with mono-dispersed particles in a pseudo-2D fluidized bed. Other influential factors such as particle size and the magnitude of heat source terms are also considered in these investigations.
In chapter 4, more complexities of the physical model are added by introducing a particle size distribution and different heat production mechanisms, namely a constant volumetric heat production ([W/m 3 ]) or constant particle heat production ([W]) (mimicking the concept of heat production obtained from a normal catalytic reaction and a polymerization reaction, respectively).
The effect of bed thickness in rectangular fluidized beds is investigated through CFD DEM simulations of the same systems in chapter 5. By keeping the aspect ratio constant, the bed depth is extended stepwise from 10 times the particle diameter to 80 times (i.e. 10, 15, 20, 40, 80d p ). Note that the bed width was kept constant and equal to 80 times the particle 1.6. Outline of the thesis 7 diameter. The influence of the third dimension on both the hydrodynamics and thermal aspects in the fluidized bed is discussed.
Experimental validation of the CFD-DEM model is presented in chapter 6. In this chapter, the particle temperature distribution obtained from experiments and simulation are compared and discussed in detail.
The thesis concludes with an overview of the main conclusion and an outlook for future investigations in chapter 7.
C h a p t e r 2 Effect of superficial gas velocity *

Abstract
The heat transfer behavior of particles and gas in an olefin polymerization fluidized bed was numerically analyzed using an in-house developed 3-D computational fluid dynamics discrete element model (CFD-DEM). First the implementation of the model was verified by comparing simulation results with analytical results. A constant volumetric heat production rate was implemented in the particle energy equation to mimic the heat production due to the polymerization reaction. It was found that the probability density function (PDF) of the particle temperature becomes more homogeneous with increasing superficial gas velocity. Furthermore, instantaneous snapshots of the thermal driving force (the difference between the single particle temperature and bed-average gas temperature, Tp-<Tg>) for different heat production rates provide detailed insight in the particle temperature distribution inside the fluidized bed. The time-and bed-averaged particle convective heat transfer coefficient, which was calculated by Gunn s correlation, was found to be independent of the superficial gas velocity. This is explained by the fact that the relative velocity of gas and particles in the emulsion phase, where most of the particles and gas interact, is hardly influenced when increasing the gas superficial velocity. From the spatial distribution of Nusselt number, it becomes apparent that the high heat transfer regions are found in the wake of rising bubbles, whereas low heat transfer rates are found in the clouds of the bubbles.

Introduction
Polyolefins (PO) are used in a wide variety of applications. A large number of olefin polymerization processes are operated with heterogeneous Ziegler-Natta catalysts in gas-solid fluidized bed reactors. Understanding of the fundamental phenomena encountered in the gas-phase polymerization reactor is essential for a proper design and optimal operation. Although the fluidized bed process for gas-phase olefin polymerization has already been operated for decades, the complex interplay of hydrodynamics, and heat and mass transfer is not yet fully understood and this still remains a challenge due to the difficulties in the description of granular flows (McKenna et al., 2005). Due to the highly exothermic nature of olefin polymerization and especially in gas phase processes a large amount of heat needs to be removed from the particles to prevent the formation of temperature gradients inside the particles and the accumulation of reaction heat that, in extreme cases, can lead to particle melting, sintering, adherence, agglomeration, which may ultimately cause defluidization. The development of catalysts with an increasingly higher activity has made this even more important. To improve the fluidization efficiency and heat management of these reactors a better understanding of the gas-particle heat transfer is required using detailed particle based computational fluid dynamics (CFD) simulations. Several different CFD models have been developed for fluidized bed reactors for olefin polymerization (Khan et al., 2014). These models have been used to investigate the hydrodynamic behavior, such as solids mixing, bubble behavior (van der Hoef et al., 2008), pressure effects Godlieb et al., 2008), and heat and mass transfer characteristics (Behjat et al., 2008). Among all these models, the Discrete Element Model (CFD-DEM) focuses on the scale of individual particles and track individual particle trajectories while accounting for particle-particle interactions, which is key to elucidate the mechanisms governing the complex bed dynamics and heat transfer behavior. For a more comprehensive overview and details of CFD-DEM applications, we refer to recent reviews (Deen et al., 2007;Zhu et al., 2007;van der Hoef et al., 2008).
The first papers that applied CFD-DEM to investigate heat transfer during gas phase polymerization were mainly focused on hop-spot formation (Kaneko et al., 1999), taking both propylene and ethylene as monomer gas. Similar work on fluidized beds with extensive heat production were reported by Behjat et al. (2008); Chang et al. (2013); Eriksson and McKenna (2004); Zhou et al. (2009). According to these studies, hop-spot formation strongly depends on the flow patterns and bed geometry. However, these results were insufficient to arrive at a general conclusion (Behjat et al., 2008;Kaneko et al., 1999). Moreover, most of these works have considered situations where thermal equilibrium in the bed has not yet been reached. Furthermore, although some results on particle and gas temperature distributions have been reported, there is a lack of detailed information on the local distribution of the heat transfer coefficient, which could elucidate the conditions for hop-spot formation and understanding of local heat transfer phenomena. More general methods to analyze hop-spot formation and the relation between hydrodynamics and heat transfer behavior are highly demanded, which are the main purposes of this work. This paper is organized as follows. First, details of the governing equations of CFD-DEM are presented. Subsequently, verification of the numerical implementation of the model was carried out by making a comparison with analytical solutions. A constant heat source for mimicking the propylene polymerization reaction heat has been added to the particle model with an estimated value from literature (Choi and Ray, 1985). In the results and discussion section, the effect of the superficial velocity on the bubble behavior, temperature probability density function (PDF) of the mono-dispersed particles, and instantaneous snapshots of the driving force, Reynolds number (Re p ) and Nusselt number (Nu p ) distribution are presented and discussed.

Governing equations
In the CFD-DEM model, the gas phase is described by the continuity equation (2.1) and the volume-averaged Navier-Stokes momentum equations (2.2): where ε, ρ g , p g and u g are the volume fraction, density, pressure and the velocity vector of the gas phase, respectively. In this work the gas density is calculated from the equation of state for ideal gases. τ g is the gas phase viscous stress tensor, which is assumed to obey the general Newtonian form (2.3): µ g and I are the gas phase viscosity and the unit tensor, respectively. S is the sink term that accounts for the momentum exchange between the gas and the particles: Here r is the position vector to the staggered velocity cell, and r a is the position vector to the particle center of mass. D is the distribution function that distributes the calculated source to the Eulerian grid with grid size V on the basis of volume weighing. β is the inter-phase momentum exchange coefficient, which can be calculated by using a drag force closure. In this work, the Ergun equation (Ergun, 1952) has been used for the dense regime and the Wen & Yu drag equation (Wen and Yu, 1966) for the dilute regime: with: and: The thermal energy equation for the gas phase is given by: where the term Q p accounts for the interface heat exchange and is given by: and T g is the temperature of the gas phase, whereas T p,a is the temperature of particle a.
A p,a is the surface area of particle a and h p is the interface heat transfer coefficient. In equation (2.10). The heat conduction flux q is given by Fourier s law: where k ef f g is the effective thermal conductivity of the gas phase, which can be expressed in terms of the intrinsic fluid thermal conductivity (k g ) as follows:

Discrete particle model
The particles are individually tracked by solving the Newtonian equations of motion accounting for both translation and rotation: The translational motion of each particle is caused by the combined effect of the far field pressure gradient force, the drag force, the gravity force, and contact forces ( ∑ a̸ =b F contact,a ) due to collisions with other particles and the confining walls. A sotft-sphere model is used to simulate the particle-particle interaction (Cundall and Strack, 1979). The rotational motion of particles is also taken into consideration, where ω a is the rotation velocity and T a is the torque due to particle contact. To calculate the heat transfer from the fluid to the particles the gas phase temperature at the particle position (T g ) is calculated by interpolating the gas temperatures from the surrounding Eulerian grid points. The thermal energy equation for every individual particle temperature (T p,a ) is given by: where q v is the volumetric heat production rate due to the exothermic reaction, and h p is the gas-particle interfacial heat transfer coefficient. The latter is calculated from the Gunn s correlation (Gunn (1978)): Nu p is the Nusselt number and Pr is the Prandtl number given by: Finally, the heat transfer coefficient of the individual particle with size d p is obtained from:

Model implementation verification
The gas phase momentum and energy equations are solved numerically with a finite difference method, in which a staggered grid has been utilized. The details of the numerical method and integration method have been reported by Patil et al. (2014) and in earlier literature by Hoomans et al. (1996). Two verification tests have been conducted to verify correct implementation of the model. The physical properties of polypropylene (PP) particles and propylene monomer gas, listed in Table 2.1, have been applied in the CFD-DEM for all the verifications and subsequent simulations. The geometrical characteristics of the bed have been summarized in Table 2.1. 6.25 × 10 −6 m Initial particle temperature T p,0 (Free falling particle test)

K
Initial particle tempertature T p,0 (Fixed bed test) 300 K Normal direction of coefficient of restitution (particleparticle) 0.60 -Tangential direction of coefficient of restitution (particleparticle) A first verification was done with a simple conductive heat transfer test with a single free falling hot PP particle (340 K), with a constant volumetric heat production ofq=1.398 × 10 7 W/m 3 . The particle was placed in a stream of propylene gas (298 K) passing through with a velocity of 0.01 m/s. The CFD-DEM result matches well with the analytical solution, as can be observed in Fig. 2.1. The analytical solution is given in the graph.
The hydrodynamics and convective heat transfer was verified by considering heat transfer in in a fixed bed with a similar geometry as used by van Sint Annaland et al. (2005) and Patil et al. (2014). In the test, initially cold PP particles with a diameter of 0.995 mm were uniformly packed in the bed leading to a uniform gas porosity (voidage) equal to 0.496. At t=0 s, a stream of hot gas was injected from the bottom. Because the particles are uniformly placed in the bed under ideal plug flow conditions, the heat transfer between the two phases can be described by a one-dimensional convection equation that can be solved analytically. The analytical solution can be found in the book of Transport phenomena (Bird et al., 2007). Similar results were reported by Patil et al. (2014), using 3.95 mm copper particles and steam. Hence, the details of the test or results will not be described here. The same conclusion can be drawn here namely, that the CFD-DEM simulation results match the analytical solution of the 1-D model very well. Finally, a case of convective heat transfer with heat production was tested. In that case, the conduction term of the gas and particle phase can be ignored due to the fact that the particles are relatively small and the Péclet number (Pe 70) is relatively high. The model verification was done by comparing CFD-DEM results with the result of 1-D convection equation. The governing equations for the 1-D model are given by: The dimensionless temperature in Fig. 2.2 is defined as difference of the particle and local gas temperature, divided by the difference of particle temperature and the inlet gas temperature((T p − T g )/(T p,0 − T g,in ). The 3-D CFD-DEM heat transfer results and the 1-D approximation for both particle phase and gas phase are in close agreement as shown in Fig. 4.2. From these test cases we draw the conclusion that the implementation of the 3-D CFD-DEM has been verified and can be used for parameter studies of heat transfer in fluidized beds.

Fluidized bed geometry and particle properties
To investigate the effect of the gas superficial velocity on the fluidization behavior, a pseudo-2-D fluidized bed configuration was used. Details of the geometry are shown in Table 2.2. At the top outlet, the pressure is equal to the atmospheric pressure (101,325 Pa), whereas at the bottom a gas inlet with uniformed inlet velocity is prescribed. No-slip boundary and free slip boundary conditions are applied to the side-walls and front and back wall of the bed, respectively. In total 8×10 4 PP particles are initially uniformly packed in the bed. The gas and particle properties are listed in Table 2.1. The time step of solving gas phase equations (∆t) and particle collision (∆t p ) are 5×10 −5 and 5×10 −6 s, respectively.
The simulations were performed by step wise increasing the superficial velocity from slightly larger than the minimum fluidization velocity (u mf = 0.215 m/s) 0.3 m/s to 1.2 m/s, in steps of 0.1 m/s. Initially, the particles are randomly packed in the bed and there is a stream of gas uniformly passing from the bottom of the bed with a temperature of 324 K. These relatively hot or cold particles will be heated or cooled down by convective heat transfer with the gas and eventually reach a state of thermal equilibrium. In order to quickly reach the thermal steady state, a simple CSTR model has been used to estimate the initial temperature

Effect of superficial velocity on bed hydrodynamics
Snapshots showing the voidage pattern of the pseudo-2D fluidized bed at different superficial velocities (u 0 ) can be used to analyze bubble behavior. All the snapshots were taken at the same moment (2 seconds) after the thermal steady state had been attainted in the simulations. The overall solids hold-up and the probability density function (PDF) are also used to quantify the effect of u 0 on the bed hydrodynamics.

Bubble behavior
When the inlet gas velocity is slightly higher than the minimum fluidization velocity, namely 0.3 m/s in the first snapshot in Fig. 2.3, one can observe small spherical bubbles that move towards the bed surface. Increasing the gas velocity u 0 leads to increasing bubble sizes because of increased coalescence rates. The bed exhibits a chaotic movement. For u 0 >0.6 m/s, bubbles in the bed are hard to be distinguished individually. Due to intensive bubble breakup and coalescence particles move in the bed without typical bubble structures. Instead, oddly shaped voids can be identified.

Average voidage distribution
An average voidage (ε) as a function of u 0 was calculated by averaging gas porosity data over 2 s and is shown in Fig. 2 emulsion phase. Notice that in the bubble region (0.85 ≤ ε ≤ 1.00) the peak clearly increases with increasing u 0 , which indicates that the fluidized bed contains more bubbles. The PDF of ε provides quantitative information on the change of bubble behavior in the fluidized bed operated at increasing u 0 and complements the observations from the snapshots of the fluidized bed shown in Fig. 2.3.

Effect on particle temperature
Particle temperature PDF In this section we will investigate the influence of the superficial velocity and volumetric heat production on the particle temperature distribution. To this end the same set of superficial velocities as discussed before were used whereas we used two different volumetric heat production rates to mimic high activity catalysts during polymerization. Fig. 2.6 gives the spatially averaged bed temperature versus the superficial gas velocity. From this figure it can be seen that the average temperatures of both phases decline with increasing superficial velocity and moreover are very close to each other (∼1 K) due to the intense fluid-particle heat transfer. The probability density function (PDF) of the particle temperatures provides the basis for the quantitative comparison. For each case, the temperature of all 8×10 4 particles was made dimensionless by dividing the difference between the particle and inlet gas temperature (T p − T g,0 ) by the difference between the melting point of PP ((T p,m =380 K) and the inlet gas temperature (T p,m − T g,0 ). By increasing u 0 from slightly larger than the minimum fluidization velocity (u mf is equal to 0.215), about 0.3 m/s to 1.0 m/s, a profound change in the particle temperature PDF can be observed, as evident from Fig. 2.7.
From the particle temperature PDF shown in Fig. 2.7a and b, it can be seen that the mean value of the temperature decreases with increasing u 0 . So, increasing u 0 is a method to remove more reaction heat from particles. This conclusion is similar to the results obtained by (Hamzehei et al., 2010;Zhuang et al., 2014). The distribution becomes narrower with increasing u 0 , especially when 0 reaches 1.0 m/s. The symmetry and sharp peak indicates a more homogeneous particle temperature distribution in the bed, and therefore smaller chances of hop-spot formation.
Similar structure and change of the dimensionless particle temperature PDF can be observed whenq is twice as large. However, all dimensionless particle temperatures are approximately doubled in magnitude for the same u 0 . Particles become over-heated, which could lead to melting if u 0 < 0.6 m/s. The peaks in the PDF at corresponding u 0 are less symmetric and wider compared to the case with lowerq. Increased heat production tends to cause a bigger possibility of formation of overheated particles and also a less homogeneous temperature distribution over the bed. The standard deviation (defined in equation 2.21) in the particle temperatures as shown in Fig. 2.8 shows a continuing decrease as a function of u 0 , which indicates that the temperature over the bed is more uniformly distributed at high u 0 . Although a similar trend was found for the higher heat production case, the absolute values of the standard deviation are approximately double the values of the lower heat production case.
It is clearly noticeable that the particle temperature PDF distribution in Fig. 2.7 is less symmetrical when u 0 ≤ 0.6 m/s. In the extreme case where u 0 = 0.3 m/s, there is a long tail reaching to the lower temperature region. This is caused by the fact that at lower u 0 the circulation and mixing of the solids is less intense, leading to the highest cooling rates for particles that are close to the bottom of the bed where the cold gas enters. More information can be obtained by inspecting the snapshots of the particle temperature distribution in the fluidized bed as shown in Fig. 2.9. In order to enable a comparison of the different cases, the particles are colored according to the driving force, which is defined as the difference between the individual particle temperature and the average gas temperature. Large temperature differences can be found in the fluidized bed at lower u 0 , however, the difference becomes less pronounced when u 0 ≥ 0.7 m/s. From the snapshots presented in Fig. 2.9, one can easily trace the route of cold gas passing through the bed, which is closely related to the particle temperatures located in the long tails in the particle temperature PDF in Fig. 2.7.
Doubling the heat production rate in the particle model does not affect the hydrodynamics of the bed, as a similar route of cold gas passing through the bed can be found. However, the temperature differences between the particles are much higher. In the case of a bubbling fluidized bed (u 0 = 0.3 m/s), the temperature difference in the bed was about 25 K for the lower heat production case and about 50 K for the higher heat production case. The fluidization characteristics (visible via the bubble shapes) also affect the range of temperature differences, which is clear from the case where u 0 = 1.0 m/s in which the difference is about 5 K and 10 K for lower and higherq, respectively.

Effect of superficial velocity on heat transfer coefficient
In the previous section, the instantaneous snapshots of the particle temperatures have been used to discuss the effect of u 0 on the thermal behavior of the particles and gas in a fluidized bed. However, the heat transfer rate, which is defined by equations 2.15 and 2.17, has no direct relation to the particle temperature. To discuss the effect of the hydrodynamics on the fluid-particle heat transfer coefficient, the bed average Nu p and Re p numbers and also the instantaneous snapshots of their local distribution were calculated. Quantified Re p and Nu p numbers can provide information both on hydrodynamics and thermal aspects of the fluidized bed. Fig. 2.11 shows the instantaneous distribution of Re p . It is clear that the particles with high Re p (300-500) are in the bubble clouds. When u 0 ≤ 0.6 m/s, the maximum value of Re p in the dense emulsion phase remains approximately constant and in the range between 5 to 100. It implies that an increase of u 0 does not lead to a change in the local gas velocity and local porosity but merely to an increased bubble flow rate following the classical two-phase theory of fluidization. It can also be observed that the gas tends to move through the bed from one bubble to another, leading to regions with somewhat higher Re p in pathways connecting the bubbles.
The effect of the different Re p number on the particle Nusselt numbers can be inferred from Fig. 2.12. In all cases the highest values of Nu p are found in the wake of the bubbles, whereas the lowest values are found in the bubble clouds. Fig. 2.13 shows the time and spatially averaged value of Nu p plotted against u 0 for two different particle diameters.

Effect of superficial velocity on heat transfer coefficient
One main feature to notice is that the overall Nu p is almost independent of the superficial gas velocity, and decreases when the particle size is smaller. With the physical properties of gas and particles and such small particle diameters as we used in the simulations, the first term of the Gunn's empirical equation (2.15) is the dominate term, therefore Nu p is approximately proportional to the particle diameter dp to the power of 0.2. Whereas the heat transfer coefficient is proportional to dp to the power of -0.8. It is expected that smaller particle diameters with smaller Re p lead to smaller values of Nu p and higher heat transfer coefficients. However, it is less trivial that the overall Nu p does not lead to higher convective heat transfer at higher superficial gas velocities. This can be interpreted as the result of no conspicuous particle/gas relative velocity changes after the bed is fluidized (u 0 ≥ u mf ). That is, the excess gas passes through the bed leading to little heat exchange with the particles. Therefore, the heat transfer efficiency reaches a maximum value for superficial gas velocities beyond u mf . We note, however, that from an instantaneous point of view, the local distribution of Nu p is affected by u 0 . More precisely, the distribution of Nu p is affected by the bubble behavior. In the wakes and clouds of bubbles high and low values of Nu p are found, respectively. Our conclusion is confirmed by experimental work by Scott et al. (2004), who heated a phosphor-bronze sphere connected with a thermocouple submerged in a fluidized bed. They proved that the single particle heat transfer coefficient increased with u 0 until u 0 ≤ u mf after it reached a maximum value. Also, these results are consistent with another simulation study by Zhou et al. (2009), although they have explained the observation by assuming that the overall heat transfer coefficient does not depend on u 0 Figure 2.11: Snapshots of distribution of particle Reynolds numbers. Note that the color scales for Re p differ amongst the snapshots. Note that different color scales are used for reasons of clarity.
because of a compromise of increased conductive heat transfer and decreased inter-particle heat conduction. In this work the inter particle heat conduction is ignored due to the fact that the collision time is relatively short (5×10 6 ). Also from their experimental work, Scott et al. (2004) concluded that the heat transfer through conduction to other particles was negligible.

Effect of superficial velocity on particle temperature
Several works explain their observed results of decreased overall average particle temperature via increased superficial velocity by the general idea that the higher superficial velocity has a positive effect on the convective heat transfer coefficient (Behjat et al., 2008;Yiannoulakis et al., 2001). The origin for the lower particle temperature, can be found by considering a steady state approximation of the bed, the derivation of which is given in A. The resulting steady state bed averaged particle temperature is given by: where the terms on the right hand side respectively represent the inlet gas temperature, the temperature difference between the gas and the particles, and the adiabatic temperature rise. The temperature difference between the particles and the gas (second term) is very small if the heat transfer rate is relatively high, at least relative to the heat production rate. In polymerization reactors, the adiabatic temperature rise is much larger, which explains that at higher superficial gas velocities, the bed temperature drops. One would be tempted to say that as a consequence of that, the exact value of the heat transfer coefficient becomes irrelevant. However, since in polymerization processes, the bed is operated close to the melting point of the particles, variations in the heat transfer coefficient, and hence in the individual particle temperature may have a large effect on successful operation. The PDF of the particle temperature indicated that when the bed is operated at high u 0 , its behavior was more isothermal with a lower possibility of hop-spot formation. It turns out that the more isothermal behavior of the bed at increasing u 0 is a result of the enhanced particle mixing, rather than an enhancement of the heat transfer efficiency.

Conclusions
where the terms on the right hand side respectively represent the inlet gas temperature, the temperature difference between the gas and the particles, and the adiabatic temperature rise. The temperature difference between the particles and the gas (second term) is very small if the heat transfer rate is relatively high, at least relative to the heat production rate. In polymerization reactors, the adiabatic temperature rise is much larger, which explains that at higher superficial gas velocities, the bed temperature drops. One would be tempted to say that as a consequence of that, the exact value of the heat transfer coefficient becomes irrelevant. However, since in polymerization processes, the bed is operated close to the melting point of the particles, variations in the heat transfer coefficient, and hence in the individual particle temperature may have a large effect on successful operation. The PDF of the particle temperature indicated that when the bed is operated at high u 0 , its behavior was more isothermal with a lower possibility of hop-spot formation. It turns out that the more isothermal behavior of the bed at increasing u 0 is a result of the enhanced particle mixing, rather than an enhancement of the heat transfer efficiency.

3
Effect of pressure *

Abstract
The effect of elevated pressure on gas-solid heat-transfer behavior in an olefin polymerization fluidized bed was numerically analyzed by using an in-house developed 3-D computational fluid dynamics model coupled with a discrete element model (CFD-DEM). To mimic the heat generation associated with the polymerization reaction, a constant volumetric heat production was incorporated in the particle thermal energy equation. Snapshots of the heat transfer driving force (the difference between single particle temperature and average gas temperature) are presented to provide insight into the particle temperature distribution in the fluidized bed. Furthermore, it was found from the Probability Distribution Function (PDF) of the particle temperature that with increasing operating pressure the average particle temperature drops and the bed becomes more isothermal. Moreover, the average particle-gas heat transfer coefficient, was found to increase with increasing operating pressure. However, it is independent from the superficial gas velocity when the bed was operated above the minimum fluidization condition at the same elevated pressure. Predictive based on a continuous stirred tank reactor (CSTR) approximation of the bed reveals that the average particle temperature at steady state is determined by the inlet gas temperature, the temperature difference of the two phases and the adiabatic temperature rise. In polymerization reactors the adiabatic temperature rise is relatively high, therefore the effect of the pressure on the average temperature of the particles can be entirely attributed to the increased gas density, whereas changes in the bed hydrodynamics have an effect on the spread of the particle temperature distribution.

Introduction
Fluidized bed reactors (FBR) operated at elevated pressure are encountered in a variety of industrial processes such as fluidized bed combustion (Oka and Anthony, 2003), coal gasification (McLendon et al., 2004), Fischer-Tropsch synthesis (Dry, 2002) and heterogeneous olefin polymerization, etc. From these processes the heterogeneous olefin polymerization benefits from operation at elevated pressure (up to 20 bar). Operation at elevated pressure offers several advantages such as improved bed-to-wall heat transfer (Jenkins et al., 1986), less particle segregation, and small reactor size (Chen and Keairns, 1975). To optimize the highly exothermic olefin polymerization process conducted in FBR, fundamental knowledge of the hydrodynamics and heat transfer characteristics as a function of the operating pressure is required. The influence of pressure on the bed hydrodynamics has been investigated by utilizing several experimental methods, i.e. with the use of electrical capacitance tomography (ECT) and X-ray tomography (Sidorenko and Rhodes, 2004;Yates and Simons, 1994). On basis of experimental studies it has been reported that the hydrodynamic characteristics of fluidized beds strongly depend on the operating pressure (Chan et al., 1987;Chitester et al., 1984;Olowson and Almstedt, 1991), i.e. the minimum fluidization velocity (Olowson and Almstedt, 1991), bubble characteristics such as growth, size and pathway in the bed (Sidorenko and Rhodes, 2004;Yates, 1996), etc.
In the past decades, due to the ongoing development of computational technology, the effects of pressure on the hydrodynamics of FBR have been investigated in detailed using CFD approaches. Li and Kuipers (2002) investigated the pressure effect with a CFD-DEM model, where they found that an elevated pressure reduces the incipient fluidization velocity, while the window (i.e. gas velocity range) for uniform fluidization increases. It narrows down the bubbling regime and leads to a quick transition to the turbulent regime. They also observed that an elevated operating pressure suppresses bubble growth and therefore produces more uniform gas-solid flow structures. A number of key features for fluidized beds, such as the bubble behavior, flow patterns, solids mixing and particle granular temperature have been investigated by Godlieb et al. (2008) using state-of-art CFD-DEM simulations. They reported that at elevated pressure, the gas-particle interaction is enhanced and becomes relatively more important as compared to the particle-particle interaction. The bubble size also decreases, whereas bubble breakage becomes more pronounced. Solids mixing is also enhanced at elevated pressure condition. It is noticed that most of these studies describe the effect of elevated pressure on the hydrodynamics of fluidized bed, while very little research has been published on the changed gas-particle heat transfer behavior or the combination of these two aspects. It is known that elevated pressures can enhance the wallto-bed heat transfer (Borodulya et al., 1991;Yew Looi et al., 2002) in fluidized beds, however, there is a lack of detailed information of the pressure effect on gas-particle heat transport. The increase of the heat transfer rate at elevated pressure can directly be related to the increased density of the gas. So, one can directly use the gas properties at elevated pressure to simulate the heat transport in a fluidized bed, like was done in the work of Kaneko et al. (1999). The gas properties (in particular the gas density) at elevated pressure conditions were used to mimic the conditions of elevated pressure to calculate the heat transfer coefficient. However, the effect of the pressure on the hydrodynamics was not considered in these simulations. Therefore, in this chapter we will study the effect of pressure on the heat transport in fluidized beds considering both the effect of the increased density and the effect of the altered hydrodynamics.
This work is organized as follows. First details of the governing equations of the CFD-DEM are given. A constant heat source for mimicking the heat generation due to the strongly exothermic polypropylene polymerization has been incorporated in the particle thermal energy equation using an estimated value from literature (Choi and Ray, 1985). In the results and discussion section, the hydrodynamics is characterized in terms of snapshots of bubble behavior, time-averaged solids circulation patterns and the time-averaged probability distribution function (PDF) of the bed porosity. Subsequently, the particle temperature PDF of mono-dispersed particles, instantaneous snapshots of the thermal driving force, particle Reynolds number (Re p ) and particle Nusselt number (Nu p ) distribution are presented and analyzed to quantitatively assess the pressure effect on the convective gasparticle heat transfer in fluidized beds. The conduction between particles is insignificant compared to the convective heat transfer to the gas phase due to short (almost instantaneous) particle-particle contact time in a vigorously gas-fluidized bed. Within the temperature range studied in this work, the influence of thermal radiation can safely be neglected as well. The Biot number of this system is smaller than one, consequently intra-particle within a unit, hereby the temperature gradients can be neglected. Finally, several simulations with artificial settings of gas properties are discussed to evaluate the effect of solely changing the gas density that is used to describe the gas phase hydrodynamics, while the total heat flux into the bed is kept similar compared to the atmospheric case, making it independent of the gas density. In this way we can identify whether the improved heat transfer is due to the increased heat flux into the bed or due to the change in hydrodynamics.

Gas phase model
In the CFD-DEM model, the gas phase is described by the continuity equation (3.1) and the volume-averaged Navier-Stokes equations (3.2): where ε, ρ g , and u g are the volume fraction, density and the velocity vector of the gas phase, respectively. In this work the gas density is calculated from the equation of state for ideal gases. τ g is the gas phase viscous stress tensor, which is assumed to obey the general Newtonian form (3.3): µ g and I are the gas phase viscosity and the unit tensor, respectively. S is the sink term that accounts for the momentum exchange between the gas and the particles: Here r is the position vector to the staggered velocity cell, and ra is the position vector to the particle center of mass. D is the distribution function that distributes the calculated source to the Eulerian grid with grid size V on the basis of volume weighing. β is the inter-phase momentum exchange coefficient, which can be calculated by using a drag force closure. In this work, the Ergun equation (Ergun and Orning, 1949) has been used for the dense regime and the Wen & Yu drag equation (Wen and Yu, 1966) for the dilute regime: with: and: The thermal energy equation for the gas phase is given by: where the term Q p accounts for the interface heat exchange and is given by: and T g is the temperature of the gas phase, whereas T p,a is the temperature of particle a. In equation (3.10). The heat conduction flux q is given by Fourier s law: where k ef f g is the effective thermal conductivity of the gas phase, which can be expressed in terms of the intrinsic fluid thermal conductivity (k g ) as follows:

Discrete particle model
The particles are individually tracked by solving the Newtonian equations of motion accounting for both translation and rotation: The translational motion of each particle is caused by the combined effect of the far field pressure force, the drag force, the gravity force, and contact forces ( ∑ a̸ =b F contact,a ) due to collisions with other particles and the confining walls. A sotft-sphere model is used to simulate the particle-particle interaction (Cundall and Strack, 1979). The rotational motion of particles is also taken into consideration, where ω a is the rotation velocity and T a is the torque due to particle contact. To calculate the heat transfer from the fluid to the particles the gas phase temperature at the particle position (T g ) is calculated by interpolating the gas temperatures from the surrounding Eulerian grid points. The thermal energy equation for every individual particle temperature (T p,a ) is given by: where q v is the volumetric heat production rate due to the exothermic reaction, and h p is the gas-particle interfacial heat transfer coefficient. The latter is calculated from the Gunn s correlation (Gunn, 1978): Nu p is the Nusselt number and Pr is the Prandtl number given by: Finally, the heat transfer coefficient of the individual particle with size d p is obtained from:

Simulation settings 3.3.1 Fluidized bed geometry and particle properties
In all simulations a pseudo-2D fluidized bed was used to investigate the effect of elevated pressure on the fluidization behavior. Details of the geometry, gas and particle properties are given in Table 3.1. Initially, a total of 80000 polypropylene particles were randomly distributed in the bed. At the bottom gas is introduced with a uniform velocity, whereas at the top of the bed a pressure outlet with a constant pressure of respectively 1 bar, 2 bar, 5 bar, 10 bar or 20 bar is used for different simulation cases, respectively. No-slip boundary conditions are applied to the side-walls, whereas a free-slip boundary conditions are used for the front and back walls of the bed. The inlet gas temperature is equal to 324 K. The ideal gas law is applied to calculate the local gas density. To enable a fair comparison, we compared simulations with a constant excess superficial gas velocity which is defined as the difference between the superficial velocity and the pressure dependent minimum fluidization velocity. The minimum fluidization velocities calculated from the Ergun equation (Ergun and Orning, 1949) (with solid volume fraction equals to 0.62) are 0.215, 0.165, 0.105, 0.085, 0.055 m/s for the different operating pressures. Firstly the simulations with u ex = u 0 u mf = 0.285 m/s are used for the discussion, subsequently more cases that were used to compare the effect of increasing superficial velocities (u 0 ) at constant excess superficial gas velocity (u ex ) are also considered in this work. The initial particle temperatures and simulation conditions are specified in

Porosity distribution
All results discussed in this chapter were analyzed after both hydrodynamic and thermal steady state were reached, and the initial temperatures of the particles and the gas phase are listed in Table B.2. Data was gathered during a period of 2 s using a time interval of 0.0125 s. From the work of Godlieb et al. (2008), it was found that the solids mixing time in a fluidized bed of identical size as used in this work was 1.5-2.0 s. Godlieb et al. (2008) studied the same range of pressures. Based on this, an analysis duration of 2 s is sufficient to characterize the bed behavior. Snapshots of porosity patterns of the pseudo-2D fluidized bed at different operating pressures (p) were used to qualitatively study the fluidization behavior. The time-averaged solids volume fraction distribution and its PDF were used to quantify the effect of the operating pressure on the bed behavior.

Bubble behavior
Fig. 3.1 shows snapshots (taken at simulation time t = 2 s) of the solids distribution in the fluidized bed at different operating pressures. From these snapshots, it is apparent that the bubbles become smaller with increasing pressure. When the bed is operating at relatively low pressure (less than 10 bar) as in the first three snapshots, one can observe several bubbles with clear edges propagating through a dense emulsion phase. The bubble shapes become irregular when the pressure exceeds 10 bar, where after it is difficult to distinguish the bubble phase and the emulsion phase. The bed expansion is higher at elevated pressure because of the dense phase expansion. From animations of the simulations, it was observed that at these pressures due to extensive bubble breakup and coalescence, the gas moves in the bed without typical bubble signature, in contrast to the behavior observed at low pressures. The chaotic movement and increased interaction between particles and gas eliminates the distinct boundary between bubbles and the dense emulsion phase at high pressures. Due to the increased gas density, the interaction between particles and gas is enhanced significantly, and therefore a more homogeneous structure of the two phases at elevated pressure is taken place in the fluidized bed. Consequently a similar behaviour as liquid fluidized bed is observed in Fig. 3.1 e. Smaller bubbles with lower rising speeds and higher break-up rate lead to tortuous bubble trajectories. These typical characteristics of fluidized beds at elevated pressure are in accordance with the findings of others (Li and Kuipers, 2002;Fan et al., 2004;Zhang et al., 2000).

Time-averaged solids motion
The time-averaged particle velocity field has been used to provide information on the solids circulation pattern. This field was obtained by assigning the particles to the Eulerian grid cells in which they reside with subsequent averaging in time and in the depth direction of the bed. The resulting averaged particle velocity fields are shown in Fig. 3 Fig. 3.2 it can clearly be seen how the pressure influences the solids circulation pattern. Particles are dragged up by the gas through the center of the bed and fall down along the side walls, creating a global circulation pattern in the bed. The symmetrical structure of the solids circulation pattern does not change with the pressure although in the bed there is more bubble coalescence and more intense gas-solid interaction, leading to a larger bed expansion as the pressure is increased.

Average voidage
To quantify the effect of pressure on the bubble behavior of the fluidized bed, firstly, the gas voidage was temporally and spatially averaged for the different operating pressures (using a cup-mixing average of the voidage in the Eulerian cells, during a period of 2 s with time intervals of 0.0125 s). The average voidage indicated by ε is reported in Table  B.2. The results of case 2 are visualized in Fig. 3.3, where it can be seen that the average bed voidage increases linearly with increasing operating pressure from 0.58 to 0.75. The elevated pressure leads to an increasing overall bed voidage which is consistent with the observations from the snapshots given in Fig. 3.1. To obtain more quantitative information, the probability distribution function (PDF) of the bed voidage was calculated and the result is shown in Fig. 3.4. For pressures lower than 10 bar, we can observe a peak in the porosity range of 0.40 to 0.55, which is associated with the emulsion phase. However, when the pressure is increased beyond 10 bar, this dominant peak disappears and most of the bed is in a state that is in between those of the dense emulsion phase and the bubble phase. With increasing pressure the bubble phase starts to contain more particles, whereas the emulsion phase contains more gas. Similar observations were reported by Godlieb et al. (2008).

Particle temperature distribution Particle temperature PDF
In this section we will investigate the influence of the elevated pressure on the particle temperature distribution. To this end, the same set of superficial velocities and operating pressures as discussed before were used. Fig. 3.5 shows the driving force for gas-particle heat transfer, which is characterized by the difference between the temperature of the individual particles (T p,a ) and the spatially averaged gas temperature over the domain (<T g >). For the  convenience of discussion, we also show the corresponding gas velocity fields underneath to visualize the main pathways of the gas moving through the bed.
Note that the temperature scales for the driving force shown in Fig. 3.5 are different, i.e. the maximum value of the driving force is around 1 K in all cases, whereas the minimum value changes considerably as a function of the operating pressure. Negative values of the driving force indicate particles with a temperature lower than the average bed temperature. These particles are primarily found near the bottom of the bed where they get in contact with the cold inlet gas. Some particles are dragged along with the gas moving along the main pathway of the gas, leading to the low temperature fingers . These fingers can be used to trace how and where the gas is passing through the bed. With the corresponding instantaneous gas velocity vector graphs given beneath the driving force distribution graphs ( Fig. 3.5 a2-e2), it clearly demonstrated that gas passes through the solids by moving from bubble to bubble. The cold particles (e.g. those with a strong negative driving force) are found between these bubbles. Besides the information obtained from the patterns and temperature differences, it is also noticeable that the differences in driving force are diminishing with increasing pressure, leading to a more isothermal bed. This effect can also be quantified by considering the probability distribution function (PDF) of the normalized particle temperature that is shown in Fig. 3.6. In this figure, the driving force for heat transfer (T p,a − T g,0 ) has been normalized by dividing it by the difference between the melting point of polypropylene and the inlet gas temperature (T p,m -T g,0 ). The melting point of polypropylene is T p,m = 380 K. Fig. 3.6 shows that at low pressures many particles obtain temperatures that are close to the melting point of polypropylene (high values of the temperature ratio is increasing). With increasing pressure, the particles obtain a more constant temperature and are further away from the melting point (the temperature ratio is decreasing). Operating at elevated pressure is not only an essential polymerization condition but also an efficient way to remove excessive polymerization reaction heat, and leads to a more homogeneous temperature distribution in the bed. In the inset in Fig. 3.7 the standard deviation of the particle temperature is also shown, which is defined as equation 2.21 in chapter 2. the standard deviation also reveals that the particle temperature is more uniformly distributed over the bed at increasing pressure. There may be several explanations for this behavior. One explanation is that increasing pressure leads to different bed dynamics, in particular enhanced solids mixing (Godlieb et al., 2008), which leads to quick replacement of the particles in the main pathway of the gas. Another explanation is that by increasing the pressure, the volumetric heat capacity of the gas is considerably increased, which gives the gas a much stronger capacity to cool the particles. We will get back to this issue in one of the next sections.

Particle Reynolds and Nusselt number distribution
In the previous section, the snapshots of the particle temperatures were used to assess the effect of the operating pressure on the hydrodynamics and thermal behavior in a fluidized bed. However, the heat transfer coefficient has no direct relation to the particle temperature. To discuss the effect of the pressure on the heat transfer coefficient, snapshots of Nu p and Re p are shown in Fig. 3.7. Increasing the operating pressure leads to both higher values of Nu p and Re p , which can be observed clearly from the increasing scales alongside each graph of Fig. 3.7. The snapshots of Nu p in graphs Fig. 3.7 (a1-e1) show that most particles with high and low Nu p value are distributed in the wake of bubbles and in the bubble clouds, respectively. These differences cannot be explained from the corresponding distributions of Re p , which are distributed such that low values of Re p are found inside the bubbles, whereas high values of Re p are found in the dense phase. The differing trends in Nu p and Re p indicate that apart from the differences in Re p also differences in the local distribution of the voidage  determines the spatial distribution of the convective heat transfer coefficient.
To assess the effect of the operating pressure on convective heat transfer coefficient between gas and particles, additional simulations at different operation conditions were conducted. Comparing increasing superficial velocities (u 0 ) at constant excess superficial gas velocity (u ex ) for increased elevated pressure in each case. The simulation conditions and the obtained results in terms of Nu p and Re p are listed in the Appendix B in Table. B.2. Case 2 is the one that has been discussed in detail in the previous sections. In Fig. 3.8 we compared the spatial and temporal averaged Re p , emulsion phase voidage (ε e , the threshold for bubble phase equals 0.85) and Nu p as a function of u 0 /u mf at pressures ranging from 1 bar to 20 bar. It is observed that the variation of ε e , Re p and Nu p with respect to u 0 /u mf shows different trends. An overall observation is that the values of εe ( Fig. 3.8 (b)), Re p (Fig. 3.8 (b)) and Nu p ( Fig. 3.8 (c)) all increase as the pressure increases. However, where Re p increases with (u 0 /u mf ) at fixed pressure Nu p weakly depends on (u 0 /u mf ).
This latter observation indicates that the observed increase of Nu p (due to high Re p ) is compensated by the effect of bed expansion which increases with increasing u0. The convective heat transfer coefficient between particle and gas in the fluidized bed reaches Figure 3.7: Snapshots of instantaneous fluid patterns of particle Nu p (a1-e1) and Re p (a2-e2) distribution. Note that different color scales are used for reasons of clarity. . its maximum value above minimum fluidization condition at different elevated pressure. Similar conclusions can be found in a recent simulation study reported by Li et al. (2016a), in which the simulations were performed at atmospheric condition. Apparently the conclusion of the independence of the convective heat transfer coefficient between particles and gas on gas superficial velocity is also valid when the bed is fluidized at elevated pressure.

An artificial case
In the previous sections it was found that the average particle/gas temperature decreased at elevated pressures, which can be used as a means to remove more heat from the particles and to prevent hot spot formation. There are two possible causes of this thermal behavior in the fluidized bed. The first possible explanation is the increased volumetric heat capacity of the gas, which scales directly with the (increased) gas density. The second explanation is the change of the bed hydrodynamics, which as we have discussed in previous sections, involves different bubble behavior, different bed expansion etc. It would be interesting to  investigate these two effects separately. To this end, we decouple the thermal and hydrodynamic bed behavior as much as possible in such a way that at increased pressures the thermal behavior of the bed is only allowed to depend on the different hydrodynamics in the bed. To achieve that, we define an artificial specific gas heat capacity C * p,g to compromise the increasing density by the pressure, keeping an approximately constant convective heat flux (C p,g ρ g,0 u 0,0 T g,0 =3.9×10 5 W/m 2 ) into the simulation domain for all the simulations: Here ρ g,0 u 0,0 and ρ g,1 u 0,1 are the mass flux into the system when the bed is operating at atmospheric and elevated pressure conditions, respectively. The corresponding temperatures are indicated by T g,0 and T g,1 , which are equal in this case. We also need to make sure that the heat transfer coefficient does not depend on the altered mass flux, so we also include the correction factor f in the Gunn's equation as follows: (3.19) note that the product of Re p and Pr remains unaltered: 3.20) and the values of C * p,g and correction factor f are listed in Table3.2 The comparison of the spatially and time-averaged particle temperatures and its standard deviation obtained from the simulations with and without modification of C p,g is given in Fig. 3.9 and Fig. 3.10, respectively.
With the modified specific heat capacity and correction factor, the influence of the gas density vanishes. We can see that the average particle temperature stays more or less unchanged with increasing pressure, despite the differences in bed behavior (e.g. different  Figure 3.9: Spatially and time-averaged particle temperature as a function of pressure for cases with the original heat capacity and with an altered heat capacity that provides the same gas convective heat flux to the bed. voidage distribution and bed expansion). However, the standard deviation of the particles decreases less in the modified cases in comparison with the original ones. To interpret the observation we utilized a continuous stirred tank reactor (CSTR) model to solve the thermal energy equations of both phases. Assuming the solid and gas phases are well-mixed in the bed we can write the governing equation for the entire gas phase: where V B and A B refer to the volume and cross-section of the bed and as is the specific fluid-particle interfacial area, defined as the total surface area of particles per unit of volume in the domain (6(1 − ε)/d p ). The entire solid phase is described by: At steady state we obtain the equations for the temperature of both phases at thermal equilibrium:  Figure 3.10: Standard deviation of particle temperature as a function of pressure for cases with the original heat capacity and with an altered heat capacity that provides the same gas convective heat flux to the bed.
These two equations show that the average particle temperature at steady state can be determined by the sum of inlet gas temperature, the temperature difference between the two phases and the adiabatic temperature rise. In polymerization reactors the adiabatic temperature rise is relatively high, the change of the particle temperature will mainly be effected by the last term of equation (3.24). This explains the reason that the spatially and time-averaged temperature of the particle remain unchanged, when the heat flux is kept constant in the artificial cases. On the other hand, the distributions of the particle temperature characterized by the particle standard deviation do show some response to the hydrodynamics. As mentioned earlier, elevated pressure can lead to faster solids mixing and more effective interaction between bubbles and emulsion phase. All in all these changes of hydrodynamics provide a homogeneous distribution of both phases, which therefore narrows the distribution of the particle temperature preventing local hot spot formation. Consequently, we can conclude that the dominant mechanism for a more favorable average thermal behavior in the fluidized bed can be attributed to the increased gas density and to the associated increased volumetric heat capacity. However, the more pronounced isothermal bed behavior is due to both the change in hydrodynamics and the increased gas density.

Conclusions
In this work, a CFD-DEM model was used to simulate a gas phase polyolefin fluidized bed with a constant heat production in the particle phase. By increasing the operating pressure from 1 bar to 20 bar, with the same excess superficial velocity, it is found that the elevated pressure influences the hydrodynamics and significantly enhances the heat transfer performance of a fluidized bed. The bubble size is reduced and a more flat distribution of bed voidage is achieved. The PDF of the normalized particle temperature shows that the mean value decreases significantly while the standard deviation decreases as the pressure is increased. Operating the bed at an elevated pressure helps to reach a more isothermal condition in the fluidized bed, thereby reducing the opportunity for hot spots to develop.
The convective heat transfer coefficient, which was characterized by the particle Nusselt number increases considerably with increasing operating pressure, which is beneficial for removing more heat from the bed. This can be entirely attributed to the increased density, whereas the changes in the bed hydrodynamics have little effect. Meanwhile the particle Nusselt number is found to be independent from the superficial gas velocity when the bed is operated above minimum fluidization condition at elevated pressure.
The effect of the pressure on the average temperature of the particles can be entirely attributed to the increased density, whereas changes in the bed hydrodynamics has its impact on the variance of the particle temperature distribution.
C h a p t e r 4 Effect of particle size distribution Abstract Many gas-phase olefin polymerization processes are executed in fluidized bed reactors. The particles often have a particle size distribution (PSD) due to a variety of factors, such as residence time distribution, initial catalyst size distribution, different rate of catalyst activity decay, etc. The heat transfer phenomena in poly-disperse fluidized beds with different particle size distributions and different heat sources have been numerically analyzed using an in-house developed 3-D computational fluid dynamics and discrete element model (CFD-DEM). Simulations have been carried out for beds with Gaussian PSD's with different widths but with the same Sauter mean diameter (d32=1.2 mm). The thermal energy equation of the particles contains a heat source to mimic the heat production due to the exothermic process where two cases have been considered: a constant volumetric heat production rate (qv, [W/m 3 ]) and a constant heat source per particle (Q, [W]) to represent different reaction systems, viz. "normal" catalytic reactions and polymerization reactions respectively. The computed probability distribution functions (PDF) of the particle temperature reveal that the temperature distribution in the fluidized bed is strongly affected by the width of the particle size distribution, the magnitude of the heat production rate in the particle phase and the superficial gas velocity. Moreover, the obtained temperature contour maps show the relation between the temperature distribution and the particle size. It was found that small particles with a high heat production rate cause the formation of hot spots in the bed. It was also found that operating the bed with a relatively high superficial velocity cannot suppress the number of particles in the high temperature region. Furthermore, snapshots of the fluidized beds demonstrate that these small particles have a higher chance to be found near the top of the bed and in the vicinity of the side walls of the reactor. The former is due to minor size segregation in the vertical direction, the latter is caused by preferential particle motion.

Introduction
Fluidized beds are often applied for gas-phase olefin polymerization (PO) processes due to their intense solids mixing and excellent heat and mass transfer characteristics. The very large heat of reaction is one of the challenges faced by the polymer industries in terms of optimizing the performance of the fluidized beds. Local hot spots can occur in the bed, and these overheated soft particles can cause sintering, agglomeration and in extreme cases lead to defluidization (Zacca et al., 1996;Mckenna et al., 1995;Zacca and Debling, 2001;Floyd et al., 1987). Although the process for gas-phase olefin polymerization has already been operated for decades, the complex interplay of hydrodynamics, and heat and mass transfer is not yet fully understood and this still remains a challenge due to the multi-scale nature of these phenomena (McKenna et al., 2005). During gas-phase polymerization the fresh catalysts (20-50 µm) are fed into the reactor above the gas distributor, meanwhile the products (1000-1500 µm) are removed from the reactor at the bottom of the bed continuously. A broad particle size distribution can be found in the bed due to several reasons, viz. particle fragmentation due to their rapid growth, residence time distribution, initial catalyst size distribution, different rate of catalyst deactivation, etc. (Zacca et al., 1996;Hatzantonis et al., 1998;Floyd et al., 1986). The motion and thermal behavior of these differently sized particles constitute a key factor in the understanding the formation of overheated particles in the reactor during fluidization.
To improve the fluidization efficiency and heat management in heterogeneous polymerization fluidized bed reactors, a better understanding of the gas-particle heat transfer is required using detailed particle-based computational fluid dynamics (CFD) simulations. Several CFD models have been developed for heterogeneous olefin polymerization processes in fluidized beds; these are reviewed by Khan et al. (2014). Several models have been used to study the hydrodynamics of fluidized beds and their interplay with solids mixing (Godlieb et al., 2008), bubble behavior (Hoomans et al., 1996) and operating pressure (Godlieb et al., 2008;Li and Kuipers, 2002). These models include CFD combined with discrete element model (CFD-DEM) and two-fluid models (TFM). These models have also been extended with heat and mass transfer to study olefin polymerization in fluidized beds (Behjat et al., 2008;Chang et al., 2013;Eriksson and McKenna, 2004). These models can in principle be used to simulate the particle temperature distribution in fluidized beds, however, the reported simulation studies were all limited to systems with a mono-disperse size distribution (Behjat et al., 2008;Kaneko et al., 1999;Zhuang et al., 2014;Li et al., 2016a).
A few publications are related to the experimental investigation of the effect a continuous PSD (Hoffmann and Romp, 1991;Wormsbecker et al., 2005). These studies both deal with the influence of the superficial velocity on size segregation, and are focused on axial and radial segregation respectively. Systems with a binary size distribution behave differently from systems with a continuous PSD (Sun and Grace, 1990). In binary systems, size segregation can be prevented when the bed is operated above the minimum fluidization velocity of the large particles in fluidized beds (Goldschmidt et al., 2003). However, with a continuous PSD, the local concentration of the different sized particles in the bed is affected by the superficial velocity (Hoffmann and Romp, 1991). More information related to the hydrodynamics in a fluidized bed with a continuous PSD can be found from recent modeling studies (Dahl et al., 2003;Dahl and Hrenya, 2005;Chew et al., 2011). Recently, a CFD model combined with population balance equations (PBEs) has been used to predict the PSD evolution for an olefin polymerization process (Fan et al., 2004;Yao et al., 2015;Yiannoulakis et al., 2001;. This model is based on a statistical method and is able to cover around two hours process time of polymerization. However, as with most large scale models, this modeling approach is not able to capture the particle-particle, particle-wall collisions or the temperature evolution of individual particles. Consequently, none of these studies reveal the formation of hot spots in the fluidized bed. Contrary to the aforementioned simulation models, CFD-DEM models do allow us to get a better understanding of the mechanism of hot spots formation in fluidized beds. In this work, the particle temperature distribution in fluidized beds with heat production has been studied by CFD-DEM for fluidized beds with poly-disperse particle size distributions (PSD). Here, we use a Gaussian size distribution with the same Sauter mean diameter (d 32 ), but with different distribution widths. Three standard deviations of the size distribution are used, corresponding to a narrow, a medium and a broad distribution. Two types of constant heat production terms have been considered in the simulations to mimic "normal" catalytic reactions and polymerization reactions, respectively. In the former, the catalytic activity and associated heat production scales with the particle size, whereas in the latter the catalytic activity is the same for all the particles, assuming that each particle contains the same number of catalytic sites. The heat production rate for the reference cases was adopted from Li et al. (2016a) using estimates detailed in Choi and Ray (1985). Details of the model and the simulation parameters are described firstly. In the results and discussion sections, the PDF of the particle temperature is analyzed for all simulations, including a reference case (mono-disperse system), particles with PSD with three different widths, two types of heat source terms with both small and large magnitudes and different superficial gas velocities. Furthermore, the temperature contours are used to relate the particle temperature distribution to the size distribution. The conduction between particles is insignificant compared to the convective heat transfer to the gas phase due to short (almost instantaneous) particle-particle contact time in a vigorously gas-fluidized bed. Within the temperature range studied in this work, the influence of thermal radiation can safely be neglected as well. The Biot number of this system is smaller than one, consequently the internal particle temperature gradients can be neglected. In the last section, the position of overheated particles in the fluidized bed is studied through examination of instantaneous snapshots of particle configurations.

Governing equations 4.2.1 Gas phase model
In the CFD-DEM model, the gas phase is described by the continuity equation 4.1 and the volume-averaged Navier-Stokes equations (equations 4.2): where ε, ρ g , and u g are the volume fraction, density and the velocity vector of the gas phase, respectively. In this work the gas density is calculated from the equation of state for ideal gases. τ g is the gas phase viscous stress tensor, which is assumed to obey the general Newtonian form (4.3): µ g and I are the gas phase viscosity and the unit tensor, respectively. S is the sink term that accounts for the momentum exchange between the gas and the particles: Here r is the position vector to the staggered velocity cell, and r a is the position vector to the particle center of mass. D is the distribution function, which is used to distribute the force exerted by the particles on the gas phase in the Eulerian grid cell with volume V . β is the inter-phase momentum exchange coefficient, which can be calculated by dimensionless drag force closures (F i ). In this work, the Beetstra/van der Hoef drag correlation (Beetstra et al., 2007) has been used.
which depends on the drag force that a mono-disperse system with the local Sauter mean diameter d 32 would have with a correction based on the ratio y i of the diameter and the local Sauter mean diameter: This correlation was derived for static arrays of particles, i.e. all particles have zero velocity. In our case, all particles are moving with an individual velocity. Hence, we calculate the particle Reynolds number on the basis of the local slip velocity: The thermal energy equation for the gas phase is given by: where the term Q p accounts for the interface heat exchange and is given by: and T g is the temperature of the gas phase, whereas T p,a is the temperature of particle a. In equation (4.10), the heat conduction flux q is given by Fourier's law: where k ef f g is the effective thermal conductivity of the gas phase, which can be expressed in terms of the intrinsic fluid thermal conductivity (k g ) as follows:

Discrete particle model
The particles are individually tracked by solving the Newtonian equations of motion accounting for both translation and rotation: The translational motion of each particle is caused by the combined effect of the far field pressure force, the drag force, the gravity force, and contact forces ( ∑ a̸ =b F contact,a ) due to collisions with other particles and the confining walls. A sotft-sphere model is used to simulate the particle-particle interaction (Cundall and Strack, 1979). The rotational motion of particles is also taken into consideration, where ω a is the rotation velocity and T a is the torque due to particle contact. To calculate the heat transfer from the fluid to the particles the gas phase temperature at the particle position (T g ) is calculated by interpolating the gas temperatures from the surrounding Eulerian grid points. The thermal energy equation for every individual particle temperature (T p,a ) is given by: where q v is the volumetric heat production rate due to the exothermic reaction, and h p is the gas-particle interfacial heat transfer coefficient. The latter is calculated from the modified Gunn s (Gunn, 1978) correlation adapted from literature (Tavassoli et al., 2015): Nu ef f p is the average Nusselt number calculated for a poly-disperse system, Re p is the Sauter mean diameter based Reynolds number given in equation (4.9). Pr is the Prandtl number given by: Finally, the heat transfer coefficient of the individual particle with size d p is obtained from:

Simulation settings Simulation geometry and physical properties
Simulations are performed for a pseudo 2D fluidized bed. The main simulation settings are summarized in Table 4.1. At the top outlet, the pressure is set equal to the atmospheric pressure (101,325 Pa), whereas at the bottom an inlet with a uniform inlet gas velocity is prescribed. No-slip boundary conditions for the gas phase are applied at the confining walls. The physical properties of polypropylene particles and propylene gas are used in this study, similar to previous studies reported in literature (Godlieb et al., 2008;Kaneko et al., 1999;Li et al., 2016a). The time step of solving gas phase equations (∆t) and particle collisions (∆t p ) are 5×10 −5 and 5×10 −7 s, respectively. Note that the time step for solving particle collisions is 10 times smaller than it is used in the previous chapters, in order to adjust the situation in which particles collide with extreme size difference. Simulations have been carried out for three different Gaussian particle size distributions (PSD) with the same Sauter mean diameter (d 32 ), but with different widths, where a mono-disperse PSD is used as a reference case (Dahl and Hrenya, 2005). To achieve a fair comparison in terms of the thermal response, the total bed mass for the different cases is kept the same. This was done by adjusting the total number of particles in the system. Detailed information on the different particle systems is given in Table 4.2. The cumulative mass distribution function (CDF) of the particle size for the different systems is illustrated in Fig. 4.1. For practical reasons, the head and the tail of the distribution are clipped beyond 2 times the standard deviation.

Heat source term
To mimic the heat generation associated with the chemical reaction, a constant heat production was incorporated in the particle thermal energy equation. In "normal" catalytic chemical reactions, e.g. CO catalytic oxidation, the catalysts contain a constant number of active sites per unit volume, leading to a constant volumetric heat production rate and henceforward the heat production rate scales with the size of the particle. On the other hand, in heterogeneous polymerization reactions, especially in those where the polymer product grows on or over the catalyst active site (McKenna and Soares, 2001;Weickert et al., 1999), the activity of the polymerization reaction is commonly related to the initial amount of active material inside each particle and does not scale with the size of the particles (Yiannoulakis et al., 2001). Hereby two generic types of the heat source terms have been identified and will be used accordingly in the simulations, i.e.: • for mimicking "normal" catalytic reactions, a constant volumetric heat source has been applied: • for mimicking heterogeneous polymerization reactions, a constant heat source per particle has been applied: Note that in the latter case, the overall heat production rate of all the particles was tuned such that it is the same as the overall heat production in the former case, enabling a fair comparison between the two cases. To summarize, Table 4.3 lists the values of heat sources used in the simulations. A factor of 5 has been used to investigate the situation in which the catalytic reactions proceed faster. Heat source term Reference σ1 σ2 σ3 q v (J/s/m 3 ) 6.7×10 5 6.7×10 5 6.7×10 5 6.7×10

Results and discussion
The average temperature of the particles in fluidized beds strongly depends on the cooling capacity of the supplied gas (ρ g C p,g u 0 ), and it can be well predicted using a simple CSTR model (Li et al., 2016a). The individual particle temperature, on the other hand, will deviate from the average bed temperature due to differences in the local heat transfer rate and differences in driving force. Since in polymerization processes, the bed is operated close to the melting point of the particles, variations in the heat transfer coefficient, and hence in the individual particle temperature may hamper successful operation. The individual particle temperature is influenced as well by the bubble characteristics and varies in the bed due to the spatial distribution of the heat transfer rate and the local temperature differences of the particles and the gas (driving force). Thus it is logical to examine the distribution of the particle temperatures. To achieve that, the PDFs of the particle temperature are used to compare the spread of the temperature in the fluidized bed as a function of the key variables (i.e. heat production rate, PSD, superficial gas velocity). Then, the temperature contours reveal the relation of the size of particles and the temperature distribution. Moreover, instantaneous snapshots of the bed are used to locate overheated particles. To reach a better understanding of the thermal behavior of the overheated particles, one of them has been traced to monitor the temperature evolution along the particle trajectory. Finally, the possible mechanism of local hot-spots formation has been identified.

Particle temperature distribution
The probability distribution function (PDF) of the particle temperatures provides the basis for a quantitative comparison between the cases with different PSDs. For each case, the temperature of all particles was normalized by dividing the difference between the particle and inlet gas temperature (T p − T g,0 ) by the difference between the melting point of PP (T p,m =380 K) and the inlet gas temperature (T p,m − T g,0 ). Fig. 4.2 shows the PDFs for the different cases. The rows indicate the different particle size distributions, whereas the columns show simulation results for different magnitudes of the particle heat sources. The two types of heat source are indicated with different colors: the cases with q v and Q are indicated in red and black, respectively. By comparing the differences between the two columns, it is obvious that the temperature of the bed increases more when the heat production rate increases. The mean value of the dimensionless temperature of the reference case with small heat production rate is located at around 0.5, however in the large heat production rate case it shifts to the right and is close to 0.6. The long tails of the distribution towards low temperatures are caused by the fact that near the inlet the circulation and mixing of the solids is relatively weak, leading to the highest cooling rates for particles that are close to the bottom of the bed where the cold gas enters (Li et al., 2016a;Patil et al., 2015b). It is interesting to notice that for the cases with a constant volumetric heat production rate, both for small and large heat liberation rates, the PDF has the same shape as the corresponding reference case. This indicates that for a normal catalytic chemical reaction (q v ), the temperature distribution in the bed will not be affected much by the particle size distribution. However, for polymerization reactions (Q) the PDF shows a strong dependence on the width of the particle size distribution. The spread in temperature increases with the width of the PSD. Moreover, a tail towards high particle temperatures is observed in those cases. Especially for cases with a large heat production rate, it is clear that more particles have a temperature that exceeds the mean bed temperature. The most extreme case is the one with the widest PSD (case σ3, 5Q), the long tail of the PDF extends to high temperatures, which indicates that some overheated particles have formed and this may trigger hot-spot formation for polymerization reactions. It is well-known that by increasing the superficial gas velocity (u 0 ) the particle temperature in a polymerization fluidized bed reactor can be reduced (Zhuang et al., 2014;Li et al., 2016a) due to the increased cooling capacity of the gas. In Fig. 4.3 the PDF of the dimensionless particle temperature of the bed operating at 3u mf is shown in the same structure as in Fig. 4.2, in which the superficial gas velocity was equal to 2u mf . The reason of only applying a small increase in the superficial gas velocity is to prevent that u 0 exceeds the elutriation velocity of the small particles, which would otherwise be elutriated from the fluidized bed. It is noticeable that the mean temperature of the bed reduces and the particles are colder in all cases compared to the corresponding cases at 2u mf . However, in the cases with broad PSD (σ3) and heat liberation at 5Q, the fraction of hot particles do not reduce compared to the case with superficial velocity at 2u mf , which indicates that operating at higher superficial gas velocity does not reduce the risk of overheated particles.

Temperature contour
In the previous section, we have seen that particles with a broad PSD and/or a high catalyst activity give rise to broad particle temperature distributions with hot particles close to the melting point. To determine which particles get overheated, more detailed and quantitative information is required. Fig. 4.4 shows contour plots of the probability density functions in terms of the particle size and temperature. Fig. 4.4 clearly shows the relationship between these two parameters. Note that all cases are shown with the same color bar to enable a clear comparison. For brevity, only cases with large heat production will be discussed here.
By comparing the two columns in Fig. 4.4, the differences in temperature distribution for the two types of catalytic activity can be clearly observed. With normal catalytic reactions, differences in particle size do not lead to large temperature variations in the bed. However, in polymerization fluidized bed reactors, the small particles containing fresh catalyst with high activity are the hot particles in the bed. In the case with a broad size distribution, only a few particles have the highest temperature, but it is exactly those that can lead to hot-spot formation in the fluidized bed. Fig. 4.5 shows similar results, but in this case for a superficial gas velocity equal to 3u mf . The same trends are observed as for the case with a superficial gas velocity of 2u mf , except that the number of cold particles has increased. From Fig. 4.4 and Fig. 4.5, it is noticed that not all small particles are hot. This suggests that the cause of overheating is also related to the hydrodynamics of the fluidized bed. In the next section, the location of these hot particles will be further investigated.

Spatial distribution of hot particles
As discussed in the previous section, the particle temperature contour plots show the relationship between the particle diameter and temperature. The cases with constant heat production (constant Q and 5Q) mimic polymerization reactions where the small particles have high activities. It is these small particles with a particle diameter of 400 µm that will   From left to right the graphs are shown the particles colored by its radius, the spatial particle temperature distribution and the location of the hot particles by different thresholds (top 1‰). The graphs in the last column are using gas streamlines and bubble surface contours to represent the bed (bubble phase is with a threshold of voidage equal to 0.85). Note that the size of the hot particles has been increased with a factor of 10 to facilitate the visualization. .7: Particle temperature evolution (left) and random single particle trajectory (right) in a fluidized bed (broad size distribution (σ3), heat production term with Q, and u 0 =2u um ).
Due to the fast solids motion in the fluidized bed, the particles are constantly mixed. By only showing particles with a temperature above a particular threshold, the locations of overheated particles can be determined. From the snapshots, it is clear that the particles with the highest temperature are mainly distributed in the freeboard and in the vicinity of the side walls. Since these particles are close to the melting points, both the thermal evolution and the motion of these particles are important aspects in our further discussion.

Particle temperature evolution
In Fig. 4.7, the evolution of the average particle temperature and the standard deviation (red line) together with the lowest and highest instantaneous particle temperatures in the bed (blue and green, respectively) are shown. Included is the thermal behavior of a random particle which had the highest temperature in the bed at some moment in time. As mentioned before, the bed is operated at the steady state, hereby the average particle temperature remains unchanged in time. At the starting time t=0.5 s, a particle which had the maximum particle temperature in the bed has been marked and tracked, and its trajectory in the bed during 5 seconds is shown in the right graph of Fig. 4.7. The position and temperature are marked with capitals in both graphs of Fig. 4.7.
At the starting point (t=0.5 s), the particle with the maximum temperature is found located at the top of the bed (noted with A), then it moves downwards in the bed in the vicinity of the side wall with hardly changing temperature until it reaches the bottom of the bed, where it encounters the cold feed gas. The particle temperature drops due to the heat exchange with the cold feed gas and reaches its minimum temperature at point B. Due to bubble passage it moves towards to the top of the bed while the particle starts to heat up again (point C) and finally reaches the maximum temperature of the bed (point D). From the trajectory and corresponding temperature evolution, the thermal history of a particle in the bed can be clearly observed and the observation confirms the finding of Fig. 4.6 that the hottest particles are situated in the freebroad and in the vicinity of the side walls. Because the density of the gas is quite low compared to the solid phase density, the gas temperature will quickly approach the solids temperature. Only in the bottom region where the particles are getting into contact with the cold stream of feed gas, the hot particles can be sufficiently cooled down.
Possible explanation of the observations can be summarized as follows: i) size segregation in the bed; ii) increased residence time at the side walls; iii) reduced cooling effect near the wall. For the first argument, we found that minor size segregation in the vertical direction of the bed can be expected for particles with a broad PSD, whereas the lateral segregation can be excluded due to the size and geometry of the bed. For the second point, it can be added that this is a direct consequence of the typical solids circulation patterns in fluidized beds with fast upward motion in the core of the bed and slow downward motion in the vicinity of the side walls. These two possible explanations of the formation of hot spots are discussed in more detail in Appendix C. In the next section the analysis of the gas-particle heat transfer coefficient will be discussed.

Heat transfer coefficient distribution
The heat transfer coefficient in this model is calculated from equation 4.19, where Nu ef f p is evaluated from the modified Gunn correlation on the basis of Direct Numerical Simulations (DNS) (Tavassoli et al., 2015). One main feature to notice is that the overall Nu ef f p is almost independent of the superficial gas velocity (Li et al., 2016a), and decreases when the particle size is smaller. With the physical properties of gas and particles and such small particle diameters as we used in the simulations, the first term of the modified Gunn empirical equation is the dominating term. Therefore Nu ef f p is approximately proportional to the particle diameter d p to the power of 0.2. Consequently the heat transfer coefficient is proportional to d p to the power of -0.8. It is expected that smaller particle diameters with smaller Re p lead to smaller values of Nu ef f p but higher heat transfer coefficients. A time-averaged spatial distribution of the particle to gas convective heat transfer coefficient is shown in Fig.  4.8. It can be clearly observed that high heat transfer coefficients are prevailing in the center of the bed, whereas low values are found in the vicinity of the wall with a lateral profile of approximately parabolic shape. In the splash zone of the bed we also see relatively high fluid-particle heat transfer coefficients.
This work is limited to pseudo-2D fluidized beds with adiabatic side walls. The simulation results help to improve the understanding of the hot-spot formation by revealing that the high activity particles experienced less cooling in the vicinity of side walls. However, the particle to wall heat transfer needs to be further investigated. either the direct particle to wall conduction or indirect heat transfer by conduction through the gas layer.

Conclusions
In this work a pseudo-2D CFD-DEM model was extended to study fluidized beds with heat production. Two different relevant cases were studied, i.e. constant volumetric heat production q v and constant heat production for each particle Q, for the purpose of mimicking normal catalytic reaction and polymerization, respectively. Mono-dispersed particles and particles with a continuous Gaussian distribution with three different widths were numerically studied. A fair comparison is made by keeping Sauter mean diameter (d 32 ) constant, and using the same total bed mass for all simulations. It is found that with the same total heat generation and same cooling capacity in the bed, particles with a constant Q show a broader temperature distribution compared to those with a constant volumetric heat production q v . The spread in temperature distribution increases as the heat generation is increased. The largest difference between the highest and lowest particle temperature in the bed occurs in the case with the broadest PSD and constant heat production per particle (i.e. polymerization). Through further inspection of the temperature contours, it is found that the hot particles that are close to the melting point are those small particles with high catalyst activity. Increasing the superficial gas velocity leads to a drop of the mean particle temperature in the bed, however the occurrence of overheated particles is not affected. The hot particles are mostly found in the free-board and near the side walls. Further investigation on particle segregation and particle occupation in the bed revealed that minor size segregation happened in the vertical direction, but not in the axial direction (details can be found in Appendix C). The convective particle-gas heat transfer rate along the side walls is lower than in the center of the bed, so that the over-heated particles undergo less efficient cooling by the gas there, possible leading to hot-spot formation, despite the fact that small particle have higher heat transfer coefficient compared to bigger particles.
This study is limited to a small size fluidized bed and a pseudo-2D geometry. Extrapolation of the findings to large scale industrial situations should be done with great care. Moreover, the simple formulation of the heat source terms used in the simulations may require further refinement for specific reaction systems. Nevertheless, this work can help to improve our understanding of hot-spot formation in fluidized bed polymerization reactors.
C h a p t e r

Abstract
The effect of bed depth of rectangular fluidized beds is investigated with CFD DEM. A heat transfer model with a constant volumetric heat production rate in the particulate phase, which was estimated from an olefin polymerization reaction, is incorporated in this CFD-DEM. The bed depth is increased stepwise from 10 times the particle diameter (a pseudo-2D bed) to 80 times (equal to the bed width). Simulation results for the different bed dimensions have been compared in terms of the hydrodynamics and thermal aspects at adiabatic conditions. The influence of the bed depth on a number of parameters is investigated: i.e. the bubble behavior, bed expansion, bed porosity distribution, solids mixing and motion at different operating conditions (increased superficial gas velocity and elevated pressure). Two-dimensional beds exhibit larger bubbles and more violent bubble breakage than three-dimensional beds of the same height fluidized at the same superficial gas velocity. A more pronounced distinction between bubble and emulsion phases is observed in fluidized beds where the bed depth exceeds 40 times the particle diameter. Moreover, twodimensional (2D) flow and three-dimensional (3D) flow behavior can be clearly observed between beds with depth less than 20 and more than 40 times that particle diameter, respectively. The influence of the third dimension on the solids mixing and associated mixing time is minor compared to the effect of the different operating conditions. Correspondingly, we find that the effect of bed depth does not directly affect the particle temperature or temperature distribution, especially not when the bed is operated in the bubbling fluidization regime at gas velocities above two times the minimum fluidization velocity.

Introduction
Rectangular domains are widely applied when performing Computational Fluid Dynamics (CFD) simulations for the purpose of fundamental studies on fluidization. Especially pseudo-2 dimensional (pseudo-2D) fluidized beds with small bed depth are often applied as they allow for non-intrusive visual techniques to study the fluidization behavior. Several characteristics of fluidization obtained from CFD simulations and experimental approaches performed with 2D or pseudo-2D fluidized beds can be directly compared, such as the particle velocity field (by applying particle image velocimetry (de Jong et al., 2012;Deen et al., 2007;Link et al., 2008), particle mixing and segregation (Dahl and Hrenya, 2005;Fan and Fox, 2008;Feng and Yu, 2010;Olaofe et al., 2011;Goldschmidt et al., 2003), and bubble properties (Busciglio et al., 2009;Deen et al., 2007), etc. Performing fluidization simulations in 2D or pseudo-2D with a relatively small domain and a small number of particles can reduce the computational cost notably (Xie et al., 2008). Moreover, it can also provide valuable information on the model development and fundamental understanding of the fluidization behavior (Feng and Yu, 2010). However, the phenomena observed in 2D or pseudo-2D fluidized beds cannot be generalized to full 3D fluidized beds (Geldart, 1970). The differences are mainly related to the particle and gas flow patterns that are affected by the close proximity of the confining walls (Glicksman and McAndrews, 1985;Kathuria and Saxena, 1987;Li et al., 2010). The bubble shapes, size, moving frequency and velocity in pseudo-2D and 3D gas-solid fluidized beds differ (Geldart and Cranfield, 1972;Geldart, 1970). The bubble coalescence (Lyczkowski et al., 2009), solids velocity field (Glicksman and McAndrews, 1985), particle mixing behavior (Bokkers et al., 2004) are all expected to be affected by the third dimension. The effect of the confining walls on the fluidization behavior can be significant when the bed depth is small. Pronounced 3D flow patterns can be observed with a bed depth larger than 40 times the particle diameter , where the transition zone is found at bed depths of around 20-40 times particle diameter.
Few of the aforementioned works have considered the effect of the bed depth on the heat/mass transfer characteristics in fluidized beds. Patil et al. (2014Patil et al. ( , 2015b recently performed investigations of the single bubble injection in a pseudo-2D and 3D beds operated at minimum fluidization conditions. The heat transfer coefficient for single bubble to the emulsion phase in a 3D domain was found to be twice as large compared to the Davidson and Harrison model. Whereas this coefficient matches well with the theoretical value of Davidson in a pseudo-2D domain. Detailed investigation addressing the effect of the third dimension with respect to the heat transfer in a freely bubbling fluidized bed are still lacking. This is especially relevant as the application of fluidized bed reactors in industries are almost exclusively in the 3D domain. The heterogeneous olefin polymerization process is one of the typical industrial application of fluidized bed reactors that requires careful consideration of the heat management. This process is well-known for its extreme exothermic nature and consequently often hotspots are found during commercial operation, which can cause severe operational issues, i.e. sintering and agglomeration (Khan et al., 2014;McKenna et al., 1998). In the previous chapters, we have summarized multi-scale modeling approaches which have been reported in literature in relation to study gas-particle heat transfer in fluidized bed. Heat transfer and particle temperature distributions in polymerization fluidized beds have been studied in 2D or pseudo-2D fluidized beds using the CFD-DEM approach (Kaneko et al., 1999;Li et al., 2016a,b;Zhou et al., 2009;Zhuang et al., 2014).
In chapter 2 and 3, we studied the effect of the operating conditions, i.e. superficial gas velocity and elevated pressure, on the particle temperature distribution in a pseudo-2D fluidized bed. Moreover, the particle size and size distribution have been also incorporated in the model to study the aspect of hot-spot formation in chapter 4. In the current chapter, we will perform simulations to investigate the effect of increased bed depth on the thermal aspect of these fluidized beds. Based on the conclusions from previous chapters, it is known that the operating conditions can have a significant influence on the particle temperature and temperature distribution. To provide a thorough comparison, the same operating conditions (step-wise increase of the superficial gas velocity (from 0.4 m/s to 0.8 m/s) and elevated pressure (from 1 bar to 20 bar)) are also incorporated in this work. Meanwhile simulations are conducted with a step-wise increase of the bed depth from 10 to 80 times the particle diameter (d p ) (note that the bed width is kept constant and is equal to 80d p ). The results from this simulation matrix can provide insightful information of the influence of the domain geometry on the hydrodynamics, and consequently on the thermal aspects of the fluidized bed.
This chapter is organized as follows: first snapshots of simulated bubble patterns in fluidized beds with increased bed depths operating at different conditions. Consequently the effect of the bed depth on the bed height, gas voidage, and solids mixing and motion is discussed in detail. Finally, the mean and standard deviation of the particle temperature, along with the particle temperature PDF are compared among all cases. In this way we can identify whether the conclusions we have drawn in the previous chapters are restricted to systems with a pseudo-2D geometry.

Gas phase model
In the CFD-DEM model, the gas phase is described by the continuity equation (5.1) and the volume-averaged Navier-Stokes equations (5.2): where ε, ρ g , and u g are the volume fraction, density and the velocity vector of the gas phase, respectively. In this work the gas density is calculated from the equation of state for ideal gases. τ g is the gas phase viscous stress tensor, which is assumed to obey the general Newtonian form (5.3): µ g and I are the gas phase viscosity and the unit tensor, respectively. S is the sink term that accounts for the momentum exchange between the gas and the particles: Here r is the position vector to the staggered velocity cell, and r a is the position vector to the particle center of mass. D is the distribution function that distributes the calculated source to the Eulerian grid with grid size V on the basis of volume weighing. β is the inter-phase momentum exchange coefficient, which can be calculated by using a drag force closure. In this work, the Ergun equation (Ergun and Orning, 1949) has been used for the dense regime and the Wen & Yu drag equation (Wen and Yu, 1966) for the dilute regime: with: and: The thermal energy equation for the gas phase is given by: where the term Q p accounts for the interface heat exchange and is given by: T g is the temperature of the gas phase, whereas T p,a is the temperature of particle a. In equation (5.10). The heat conduction flux q is given by Fourier s law: where k ef f g is the effective thermal conductivity of the gas phase, which can be expressed in terms of the intrinsic fluid thermal conductivity (k g ) as follows:

Discrete particle model
The particles are individually tracked by solving the Newtonian equations of motion accounting for both translation and rotation: The translational motion of each particle is caused by the combined effect of the far field pressure force, the drag force, the gravity force, and contact forces ( ∑ a̸ =b F contact,a ) due to collisions with other particles and the confining walls. A sotft-sphere model is used to simulate the particle-particle interaction (Cundall and Strack, 1979). The rotational motion of particles is also taken into consideration, where ω a is the rotational velocity and T a is the torque due to particle contact. To calculate the heat transfer from the fluid to the particles the gas phase temperature at the particle position (T g ) is calculated by interpolating the gas temperatures from the surrounding Eulerian grid points. The thermal energy equation for every individual particle temperature (T p,a ) is given by: where q v is the volumetric heat production rate due to the exothermic reaction, and h p is the gas-particle heat transfer coefficient. The latter is calculated from Gunn s (Gunn, 1978) correlation: Nu p = (7−10ε g +5ε 2 )(1+0.7 Re 0.2 p Pr 0.33 )+(1.33−2.40ε+1.20ε 2 ) Re 0.7 p Pr 0.33 (5.15) Nu p is the Nusselt number and Pr is the Prandtl number given by: Finally, the heat transfer coefficient of every individual particle with size d p is obtained from:

Physical model and simulation geometry
The particle and gas phase physical properties applied in the simulations are listed in table 2.1, and are the same as used in the previous chapters. Simulation cases at the same operating conditions discussed in Chapter 2 and 3 are reported here with increased bed depth from 10 (in chapter 2 and 3) to 15, 20, 40 and 80 d p . Note that the bed width is constant and kept equal to 80 d p . The aspect ratio based on the static bed height is kept the same for all cases, therefore the total particle numbers used in simulations increase linearly as a function of the original bed depth, i.e. these are 80,000, 120,000, 160,000, 320,000, and 640,000, respectively. Details of the bed geometry are shown schematically in Fig Initially, the particles are randomly packed in the bed and the gas uniformly enters from the bottom of the bed at a temperature of 324 K. The initial temperature of both gas and particles is obtained from the equilibrium temperature obtained from the reference cases (D=10d p ) discussed in Chapters 2 and 3. Data are obtained for 10 seconds of simulation time after the thermal steady state is being achieved. No-slip boundary conditions are applied and all simulations are conducted using adiabatic walls.  Similar instantaneous snapshots of the bubble patterns are shown in Fig. 5.3, however this time at an elevated pressure of 20 bar. The bubble shapes become irregular compared to those presented in Fig. 5.2 at atmospheric conditions. Consequently, the bubble phase and the emulsion phase cannot be clearly distinguished. From animations of the simulations, it is observed that at high pressures, the gas phase moves in the bed without a typical bubble signature (Li et al., 2016b;Godlieb et al., 2008;Li and Kuipers, 2002) due to extensive  bubble breakup and coalescence. The chaotic movement and increased interaction between particles and gas eliminates the distinct boundary between bubbles and the dense emulsion phase at high pressures. Smaller bubbles with lower rising velocity and higher break-up rate lead to tortuous bubble trajectories. As it was discussed in chapter 3, the gas density increases at elevated pressure. A more homogeneous structure of the gas and particle phases is found in the fluidized bed due to the intensive particle-gas interaction. Hereby, the bubble features are strongly affected by the high pressure, whereas they are less influenced by the bed geometry.

Bed height
In Fig. 5.4, the bed height of the fluidized bed is monitored for different bed depths. Inspection of the four graphs in Fig. 5.4 indicates that the bed expansion reduces whereas fluctuations become less intense as the bed depth is increased. Especially when the gas velocity exceeds 3 times u mf (u 0 ≥ 0.6 m/s), vigorous fluctuation of the bed height can be clearly noticed for bed depths of 10-20d p , whereas the fluctuation is less intense in the 3D beds (D= 40 d p and 80d p ). When comparing the simulation results at different superficial gas velocities, one can conclude that increasing u 0 leads to larger bubbles and higher a coalescence rate, consequently the bed exhibits a more chaotic behavior. The increased third dimension in the 3D beds provides more space to form small bubbles, and thus is less prone to bubble coalescence. The breakage of large bubbles leads to vigorous splash zone dynamics in pseudo-2D beds and explains the highly fluctuating bed height in pseudo-2D beds. However, at elevated pressure the fluctuations in the bed expansion show no significant dependence on the bed depth, as shown in Fig. 5.5. As previously discussed, the bubbles are smaller at high pressure, and the bubble shapes become irregular when the pressure exceeds 10 bar, after which it is difficult to distinguish the bubble phase and the emulsion phase (Li et al., 2016b). The chaotic movement and increased interaction between particles and gas eliminates the distinct boundary between bubbles and the dense emulsion phase at high pressures.
The time-averaged bed height (over 8 s) of all simulation cases are shown in Fig. 5.6. The left and right graph are the average bed height as a function of u 0 and pressure, respectively. Overall speaking, it is observed that the average bed height is reduced in 3D beds in comparison to pseudo-2D beds. There are 10-20% differences between the cases with depth of 10d p and 80d p , however only very minor differences can be observed between the cases with a bed depth of 40d p and 80d p . 3D flow is achieved when the bed depth is larger than 40d p , but a further increase of the bed depth hardly affects the bubble behavior. It is

Bed porosity PDF
To quantify the effect of bed depth on the bubble behavior in the fluidized bed, the probability density function (PDF) of the bed voidage was calculated and presented in Fig. 5.7. The data were collected in 50 bins ranging from 0.25 to 1.00 utilizing 8 s of simulation time. In Fig. 5.7, the peak at ε=0.45 in the dense emulsion phase (0.4 ≤ ε ≤ 0.55) increases (slightly) with increasing bed depth. Notice that in the bubble phase (0.8 ≤ ε ≤ 1.0), the peak is slightly higher with larger bed depth. In other words, a slightly better distinction between the bubble and emulsion phase is found in 3D fluidized beds. The PDF of bed porosity ε in the fluidized bed operated at high pressure (20 bar) with increased bed depth is shown in Fig. 5.7 (b). In pseudo-2D beds (D < 20d p ), ε is distributed in the range of dense emulsion phase and the bubble phase. This finding is consistent with the snapshots presented in Fig. 5.3. In 3D beds (D ≤ 40d p ), a dominant peak is observed in the emulsion phase (0.4 ≤ ε ≤ 0.55), similar as observed at atmospheric conditions. With increasing bed depth the bubble phase starts to contain more particles, meanwhile the emulsion phase contains less gas.

Solids motion Solids circulation pattern
The time-averaged particle velocity field (data collected during 8 s simulation time with a time step of 0.1 s) has been computed to provide information on the solids circulation pattern (the velocity field only considers x and z direction). This field was obtained by assigning the particles to the Eulerian grid cells in which they reside with subsequent averaging in time and in the depth direction of the bed. Fig. 5.8 presents the influence of the increased bed depth on the solids circulation pattern. Generally speaking, in a fluidized bed particles move upwards through the center of the bed and then fall down along the side walls, creating a global circulation pattern in the bed. The symmetrical structure of the solids circulation pattern can be found when the bed depth is relatively thin (D=10 and 20d p ), as shown in Fig. 5.8 (a). On the contrary, an asymmetrical structure appears when the bed depth is increased. Especially when the bed depth is equal to 80d p as shown in Fig. 5.8 (d), the circulation of the particles split into top and bottom parts. The increased third dimension provides particles freedom to move in the depth direction. Consequently, small bubbles coalesce to provide multiple larger ones instead of forming one extremely large bubble (as often observed in pseudo-2D beds). More intense gas-solid interaction prevails in the central region of bed because the particles move freely in horizontal directions in 3D bed. Thus instead of being pushed up by one big bubble, particles move along with multiple smaller bubbles leading to the asymmetrical structure of the solids circulation pattern observed in 3D beds.
In Fig.5.9 the time-averaged lateral profiles of the solids phase velocity for the different bed depths are presented at four different heights in the bed. The solid lines and symbols indicate the results of pseudo-2D and 3D beds with a bed depth of 15d p and 80d p , respectively. The averaged axial velocity profiles are compared at the same heights (correspondingly 0.025, 0.05, 0.075 and 0.1 m) above the gas distributor. It can be observed that the lateral solids velocity profiles differ strongly between the cases with different bed depths. A symmetric velocity profile with up flow in the center and down flow along the walls are clearly shown in the case of pseudo-2D beds, whereas an asymmetric profile is found in

Mixing in three directions
The effect of the depth on solids mixing is briefly discussed in this section. Investigating the mixing behavior of particles helps to understand the thermal homogeneity of fluidized beds. The average bed height method, which is based on the average height of a group of coloured particles, reported by Hoomans et al. (1996) and Godlieb et al. (2009) is adapted to study the mixing index of the fluidized beds here. The average bed height method is used to calculate the mixing index in the x, y and z direction. Initially, half of the particles are indexed by a colour, subsequently the average position of all particles is monitored. In the first step of the algorithm the vertical positions of all particles are sorted to determine the median height. The lower half of the particles is colored white, while the upper half is colored black. For each time step the average height of the white particles is monitored and normalized with the average height of all particles: where whitez white is the height normalized by the average vertical position of the white particles. The initial height of the white particles is 0.5 (half of the bed height). The average bed height is calculated based on the z position of all the particles: where j indicates the color of particles (white, black or all). Hereby the mixing index based on the average height method is defined as below: Initially the mixing index M is equal to 0 as the white and black particles are completely segregated. Eventually the two colored particles will be fully mixed with mixing index M equal to 1. However, the average height method shows periodic overshoots, which is caused by the global circulation pattern of particles, as shown in Fig. 5.10. Therefore the mixing index curve is fitted with a damped harmonic function (equation shown in the graph) as proposed by Godlieb et al. (2009). By only considering the dampened part and excluding the oscillatory part, a smooth monotonously increasing curve is obtained. In Fig. 5.10, we observe M with an overshoot in the z direction (vertical direction), whereas it evolves in the other two directions (x and y) rather smoothing. The fitting equations are applied in all three directions and the resulting fitted mixing indexes as a function of time for bed depths of 10d p and 80d p are illustrated in Fig. 5.11 for three different superficial velocities.
In Fig. 5.11, we present the mixing index evolution monitored in the x, y and z direction in the bed. Data were obtained for bed depths of 10d p and 80d p to represent pseudo-2D and 3D situations. Both cases are distinguished by using solid lines and symbols respectively in the graphs. The effect of the superficial gas velocity on the mixing index is included in the analysis. For the sake of simplicity, simulation results at u 0 equal to 0.4, 0.5 and 0.8 m/s are presented by different colors, i.e, respectively black, red and green. A general observation is that the mixing index is strongly affected by the superficial gas velocity, which indicates that faster global solids mixing can be achieved when the bed is operated at higher u 0 . Meanwhile the mixing index in both x and y directions (width and depth) increases much faster in 3D beds compared to pseudo-2D beds. However, the mixing index in the z direction shows little dependence on the bed depth. The vertical solids mixing shows a similar trend in both beds with 10d p and 80d p . Apparently the macroscopic movement of the solids in the vertical direction is not affected by the proximity of the front and back walls.

Mixing in Eulerian cells
From the previous section, it can be concluded that the solids movement in the vertical direction of the bed is not affected by the other two dimensions, however, faster mixing in the horizontal directions (x and y) has been observed when increasing the bed depth. Correspondingly it is interesting to have a closer look into the effect of the bed depth on solids "micro-mixing". The microscopic mixing of the solids in a fluidized bed can be analyzed by monitoring the Lacey mixing index (Lacey, 1954), which is based on a statistical analysis and was developed in order to estimate the chemical component proportion in mixtures. The compositions are correspondingly defined as colors of particles. Initially the particles are colored based on their vertical position (top white and bottom black, for instance). To obtain the Lacey mixing index, first the standard deviation of the proportion of white particles is calculated as equation 5.21: (5.21) where N is the number of samples (total number of cells in the domain), and ϕ i is the proportion of white particle in one cell, whereas ϕ m is its mean. The initial variance of white particles is denoted as S 0 and a random mixing status of mixture noted as S R , which are respectively defined as: where n is the total number of particles. By normalizing the standard deviation of the white particles, the Lacey index is monitored by comparing the variance of proportion evolution during the mixing process: Due to the use of grid cells the Lacey index is grid-dependent and the applied cell size has a direct influence on the results of mixing. One can consider two extremes. If the cell size is equal to the domain size, the particles are always well-mixed, and if the cell size is equal to (or smaller than) the particle size, particles are never mixed by the definition of Lacey mixing index. In this work, the cells used to calculate the Lacey index are the same as the Eulerian cells used to solve the constitutive equations of the gas phase. To compare the different simulation cases in a quantitative way, the mixing time is recorded at the moment that the mixing index reaches a value of 0.80. To prevent noise from influencing the results, a dampened exponential function has been adapted to fit the mixing index curve (Deen et al., 2010, equation 35). The results of the mixing time at different operating conditions in fluidized beds with different bed depth are presented in Fig. 5.12. Significant influence of the operating conditions can be observed as with increased u 0 and operation pressure, the mixing time is decreased. However, there appears to be no systematic effect of the bed depth on the solids mixing time.

Particle temperature
In this section the influence of the bed depth on the thermal aspects of the fluidized beds is investigated. As described in chapter 2 and 3, a CSTR model has been applied to estimate the thermal steady state of the gas and particle phase (more details are presented in Appendix A). The initial temperatures used in the pseudo-2D fluidized bed (D=10d p ) (Li et al., 2016a,b) are directly applied in the other simulations. The spatially averaged bed temperature at a thermal steady state only varies with bed height, but not with bed depth as expected (see Fig. 5.13). For brevity, the results of other bed depths studied in this work are not presented here. As it was previously concluded from simulation results obtained from pseudo-2D beds, the averaged particle temperature decreases as u 0 increases. Correspondingly, this conclusion can also be extended to 3D fluidized beds.

Deviation of particle temperature
The standard deviations (defined in equation 2.21) of the particle temperatures in fluidized beds are shown in Fig. 5.14 for different bed depths. The different colors are used to represent cases with different superficial gas velocity and the line styles are used to distinguish the different bed depths.
The variance of particle temperature fluctuates around a mean value, meanwhile it decreases as u 0 is increased, correspondingly the fluctuation becomes less intense. Moreover, the bed depth does not appear to have an influence on the variance, which is evident from the overlapping curves. Due to the vigorous movement of bubbles, particles are continuously mixed, hence a better thermal homogeneity is achieved as the mixing time reduces. Increasing the superficial gas velocity and operation pressure can both lead to improved solids mixing in the bed, at least to a certain extent. To obtain a better comparison among all cases, the time-averaged variance of the particle temperature is plotted as a function of bed depth and shown in Fig. 5.15. It can be clearly noticed that the bed thermal homogeneity is significantly affected by the operating conditions. However, there is no dependence on the bed depth. More detailed information is obtained by inspecting the probability density function (PDF) of the particle temperature.
The probability density function (PDF) of the particle temperatures provides the basis for a detailed quantitative comparison, as shown in Fig. 5.16. Examples are taken to compare the effect of bed depth on the particle temperature PDF at different superficial gas velocity (Fig. 5.16a-c) and elevated pressures (Fig. 5.16d). The particle temperature PDF obtained from the cases with bed depth D = 80d p are shown with red lines whereas the cases with D = 10d p are shown in black, respectively. It is observed that there is only a marked difference between both bed depths at low superficial gas velocity (u 0 =0.4 m/s). The distribution of the particle temperature shows a symmetrical bell shape in the 3D bed. However, with a further increased superficial velocity, the PDFs of the particle temperature obtained with different bed depths are in close agreement. In both cases long tails reaching to the lower temperatures can be observed. These are related to the particles with the highest cooling rates, which can be found close to the bottom of the bed where the cold gas enters. So far, it can be concluded that the bed depth has a minor influence on either the spatially averaged particle temperature and the temperature distribution. Especially when the fluidized bed is  Figure 5.14: Particle temperature variance as a function of time in the fluidized bed with different bed depth. Different line style and color of lines are used to distinguish cases with varied bed depth and superficial gas velocity in m/s, respectively.
operated at higher superficial gas velocity, hardly any difference can be observed between the thermal aspects in a pseudo-2D or 3D bed. This can be further confirmed by examining the axial profile of the gas temperature. The latter is averaged over the cross section as a function of bed height and shown in Fig. 5.17. Colors are used to differentiate between different operating conditions, whereas line styles are used to differentiate between the different bed depths. From the graph, a rapid rise of the gas temperature versus the bed height can be observed in all cases. The sharp gradient of gas temperature can be seen up to a height of around 15d p (0.015 m) above the gas distributor. As aforementioned, due to the low thermal capacity of the monomer gas, the temperature quickly responds to the surrounding (hot) particles. Although the cold gas stream passes through the bed as bubbles, particles that are actually carried by the bubbles are only few in number compared to those in the dense emulsion phase. Therefore, the efficient cooling zone of the fluidized bed is confined to the bottom of the bed. Consequently, the vertical movement of the particles becomes crucial to establish a certain temperature distribution.
A further supportive evidence of this conclusion is provided by the simulations performed at elevated pressure. The gas axial temperature profiles versus bed height are shown in Fig. 5.17b, where it can be clearly observed that it takes relatively long (z ≈ 80d p ) for the gas to reach a steady temperature compared to the simulations at atmospheric condition. Due to the elevated pressure, the gas phase heat capacity increases linearly with respect 10 dp 15 dp 20 dp 40 dp 80 dp 10 dp 15 dp 20 dp 40 dp 80 dp to the density, and therefore a longer time is needed for the gas phase to be heated up by the particles, thus a narrower particle temperature distribution can be achieved at elevated pressure.
To summarize, all the observations show that the operating conditions have a large effect on the various thermal aspects of the bed. However, the influence of the bed depth is very minor, which implies that the results obtained in chapters 2 and 3 for pseudo-2D beds can be extended to 3D beds as well.

Conclusions
The effect of the bed depth of rectangular fluidized beds has been investigated with the aid of CFD DEM simulations. A heat transfer model with constant volumetric particle heat production rate has been incorporated in the CFD-DEM to mimic an exothermic polymerization reaction. Simulation results of various bed depthes have been compared in terms of the hydrodynamics and thermal aspects in fluidized beds.
Instantaneous snapshots of the bubble phase in fluidized beds with stepwise increased bed depth were used to qualitatively show the influence of the third dimension on the bubble behavior. Increasing the bed depth leads to a decrease in the mean bed height and its fluctuations. A more clear distinction between the bubble and emulsion phases is observed in fluidized beds when the bed depth is more than 40 times the particle diameter. Increasing the third dimension can accelerate the solids mixing in the lateral directions to a certain extent, however it only has a minor influence on the mixing behavior in the vertical direction. From the mixing times derived from the Lacey index no clear evidence was found that the bed depth has any influence on the mixing time. However, it is found that the solids .16: Comparison of particle temperature PDF in pseudo-2D and 3D geometry simulations when changing the operating condition. Note that the minimum fluidization velocity for all cases can be found in Table B.2 in Appendix B.
mixing behavior is strongly related to the operating conditions. Two-dimensional (2D) flow behavior is observed in the bed having a depth up to 20 times the particle diameter, whereas three-dimensional (3D) flow is found in beds with depths greater than 40 times the particle diameter.
Regarding the thermal aspects, it is found that when the fluidized bed is operated at a low superficial gas velocity (less than twice the minimum fluidization velocity), the influence of the third dimension can be significant. The Probability Density Function (PDF) of the particle temperature is more symmetric in a 3D than in a pseudo-2D bed at the same operating conditions. However, when the fluidized bed is operated at a superficial gas velocity larger than two times the minimum fluidization velocity, a good agreement is achieved between the particle temperature distributions obtained in pseudo-2D and 3D beds. In the fluidized bed the cold gas entering the bed is quickly heated up by the particle phase due to the relatively low thermal capacity of the gas phase, leading to the fact that the zone with  noticeable temperature difference between the phases of the fluidized bed is confined to the bottom of the bed. However, there is no effect of the bed depth on the rate of gas phase cooling. Henceforward, we conclude that to study the thermal behavior of fluidized beds, it is possible to conduct simulations using a pseudo-2D domain. This has the benefit of reduced computational costs and the opportunity of using non-intrusive imaging techniques for validation purposes.
C h a p t e r

Experimental validation
Abstract As a result of highly exothermic reactions during gas-phase olefin polymerization in fluidized bed reactors, difficulties with respect to the heat management play an important role in the optimization of these reactors. To obtain a better understanding of the particle temperature distribution in fluidized beds, a high speed infrared (IR) camera and a visual camera have been coupled to capture the hydrodynamic and thermal behavior of a pseudo-2D fluidized bed. The experimental data were subsequently used to validate an in-house developed computational fluid dynamics and discrete element model (CFD-DEM). In order to mimic the heat effect due to the exothermic polymerization reaction, a model system was used. In this model system, heat is released in zeolite 13X particles Geldart D type) due to the adsorption of CO 2 . All key aspects of the adsorption process (kinetics, equilibrium and heat effect) were studied separately using Thermogravimetric Analysis (TGA) and Simultaneous Thermal Analysis (STA), and subsequently fluidized bed experiments were conducted, by feeding gas mixtures of CO 2 and N 2 with different CO 2 concentrations to the bed, where the total heat of liberation could be controlled. The combined infrared/visual camera technique generated detailed information on the thermal behavior of the bed. Furthermore, the comparison of the spatial and temporal distributions of the particle temperature measured in the fluidized bed with the simulation results of CFD-DEM provides qualitative and quantitative validation of the CFD-DEM, in particular concerning the thermal aspects.

Introduction
Polyolefins (PO) are used in a wide variety of applications. A large number of olefin polymerization processes are operated with heterogeneous Ziegler-Natta catalysts in gas-solid fluidized bed reactors. Due to the highly exothermic nature of olefin polymerization and especially in gas phase processes a high rate of heat removal from the particles needs to be attained to prevent the formation of temperature gradients inside the particles and the accumulation of reaction heat that, in extreme cases, can lead to particle melting, sintering, adherence and agglomeration, which may ultimately cause severe operational problems Zacca et al. (1996);McKenna et al. (2005). The development of catalysts with increasingly higher activity has motivated the research on this topic. To improve the fluidization efficiency and heat management, a better understanding of the gas-particle heat transfer is required which can be obtained from detailed particle-based computational fluid dynamics (CFD) simulations.
Several different CFD models have been developed for fluidized bed reactors for olefin polymerization Khan et al. (2014). Early work that applied CFD-DEM to study heat transfer in a gas phase polymerization process was mainly focused on hot-spot formation Kaneko et al. (1999), taking both propylene and ethylene as monomer gas. Similar work on fluidized beds with extensive heat production were reported by Behjat et al. (2008); Chang et al. (2013); Eriksson and McKenna (2004). Recently the effect of operation conditions, such as superficial gas velocity and elevated pressure have been investigated in detail with CFD-DEM models Li et al. (2016a,b). In these studies a constant heat production rate estimated from polymerization data has been applied in the energy balance of the particle phase. The CFD-DEM model can reveal all the details of the spatial distribution of the particle temperature. It can thereby provide information to investigate the mechanism of hot-spot formation.
In literature relatively few studies exist addressing the experimental validation of CFD-DEM models for heat transfer. In this work a high resolution infra-red optical technique combined with Particle Image Velocimetry and Digital Image Analysis (PIV/DIA) has been applied to study fluid-particle heat transfer in gas-fluidized beds. Experimental methods to study the hydrodynamics aspects of fluidized beds have been developed for years, such as Electrical Capacitance Tomography (ECT) Banaei et al. (2015) and X-ray tomography Kantzas and Kalogerakis (1996); Yates and Simons (1994) and magnetic resonance particle tracking or positron emission particle tracking Pore et al. (2015) and particle image velocimetry (PIV) van Buijtenen et al. (2011); de Jong et al. (2012). Infrared tomography is also a mature and commonly used measurement technique Fischer Verlag et al. (1998), however, it is limited by the fact that infrared measurements can only be applied to measure the surface temperature of objects. Patil et al. (2015a) studied hydrodynamics and heat transfer in a pseudo 2D fluidized bed using a combined CFD-DEM modeling approach and non-intrusive monitoring using an integrated PIV/DIA/IR technique. It was found that the main heat removal mechanism from the particle phase was due to the gas phase convective transport and heat loss to the environment through the confining walls. The heat loss by direct particle and wall conduction is minor compared to these two effects. However, the validation was conducted by studying a cooling process of glass beads without heat production. In this work, we extend the work of Patil et al. (2015a) to include heat production in the particles. In order to mimic the heat production in the exothermic polymerization, a model system was used. In this model system, heat is released in zeolite 13X particles (1.8-2.0 mm, Geldart D type) due to the adsorption of CO 2 Siriwardane et al. (2005). This is the first time (as far as the authors' knowledge) that the validation of a CFD-DEM with heat transfer by experimental technique has been done by involving heat production directly from the particle phase. The kinetic model of adsorption and the corresponding heat production rate has been implemented in our in-house developed CFD-DEM.
Step-by-step characterization of the heat source term, including understanding of the adsorption kinetics and associated heat production was performed.
This paper consists of two main parts, which are the CFD-DEM modeling and the experimental validation. First, details of the governing equations of CFD-DEM are presented. Subsequently, the used optical apparatus and image processing technique are described. In the results and discussion section, we first characterize the heat loss to the reactor walls, employing a simple CSTR model. The simple model was mainly used to fit the gas-to-wall heat transfer coefficient required in our CFD-DEM to account for heat isolateds to the environment. We will also show the instantaneous snapshots of the particle temperature and report the overall particle temperature evolution measured by the IR camera and CFD-DEM model. The monitoring time is up to 80 seconds, covering the entire duration of the experiments. In order to keep the main structure of the paper simple, the equations of the CSTR model and the characterization of the heat source term of adsorption by zeolite 13X particles are presented in the appendix D.

Modeling method
In the CFD-DEM model, the gas phase is described by the continuity equation (6.1) and the volume-averaged Navier-Stokes equations (6.2): where ε, ρ g , and u g are the volume fraction, density and the velocity vector of the gas phase, respectively. In this work the gas density is calculated from the equation of state for ideal gases. τ g is the gas phase viscous stress tensor, which is assumed to obey the general Newtonian form (equation (6.3)): µ g and I are the gas phase viscosity and the unit tensor, respectively. S is the sink term that accounts for the momentum exchange between the gas and the particles: Here r is the position vector of the staggered velocity cell, and r a is the position vector of the particle center of mass. D is the distribution function, which is used to distribute the force exerted by the particles on the gas phase to the Eulerian grid with cell volume V . β is the inter-phase momentum exchange coefficient, which can be calculated by dimensionless drag force closures (F i ). In this work, the Ergun equation (Ergun and Orning, 1949) has been used for the dense regime and the Wen & Yu drag equation (Wen and Yu, 1966) for the dilute regime: with: and: The thermal energy equation for the gas phase is given by: where the term Q p accounts for the interface heat exchange and is given by: and T g is the temperature of the gas phase, whereas T p,a is the temperature of particle a.
In equation (6.10), the heat conduction flux q is given by Fourier's law: where k ef f g is the effective thermal conductivity of the gas phase, which can be expressed in terms of the intrinsic fluid thermal conductivity (k g ) as follows:

Discrete particle model
The particles are individually tracked by solving the Newtonian equations of motion accounting for both translation and rotation: The translational motion of each particle is caused by the combined effect of the far field pressure force, the drag force, the gravity force, and contact forces ( ∑ a̸ =b F contact,a ) due to collisions with other particles and the confining walls. A sotft-sphere model is used to simulate the particle-particle interaction (Cundall and Strack, 1979). The rotational motion of particles is also taken into consideration, where ω a is the rotation velocity and T a is the torque due to particle contact. To calculate the heat transfer from the fluid to the particles the gas phase temperature at the particle position (T g ) is calculated by interpolating the gas temperatures from the surrounding Eulerian grid points. The thermal energy equation for every individual particle temperature (T p,a ) is given by: where q v is the volumetric heat production rate due to the exothermic reaction, and h p is the gas-particle interfacial heat transfer coefficient. The latter is calculated from the Gunn s (Gunn, 1978) correlation: Nu p = (7−10ε g +5ε 2 )(1+0.7 Re 0.2 p Pr 0.33 )+(1.33−2.40ε+1.20ε 2 ) Re 0.7 p Pr 0.33 (6.15) Nu p is the Nusselt number and Pr is the Prandtl number given by: Finally, the heat transfer coefficient of the individual particle with size d p is obtained from:

Experimental approach 6.3.1 Experimental set-ups
The experiments have been carried out in a pseudo-2D fluidized bed. A schematic overview of the position of the fluidized bed and the cameras is shown in Fig. 6.1. The fluidized bed is 8 cm in width, 1.5 cm in depth and 25 cm in height. The cameras are positioned at a distance of 1.0 m from the bed to ensure that a full image of the bed can be recorded by both cameras while reflections are avoided. The cameras are fixed on two tripods with an 15 • angle to the vertical position of the front window, and in the same plane 0.8 m above the ground. The front window of the fluidized bed column is made of sapphire glass with 3 mm thickness. Sapphire covers the IR spectrum measured by the camera (middle range, 1.5-5.1 µm) and allows up to 0.82-0.85 transmittance in the measurement temperature range (below 100 • C). The side and back walls consist of PMMA of 20 mm thickness. Its low thermal conductivity (λ=0.17-0.19 [W/(m · K)]) helps to limit the heat loss through conduction to the environment. PMMA can also reduce the interference relate to an increase of the temperature of the back plate during the experiments (C p =1466 [J/(kg · K)]) compared to the aluminum plate (λ=205 [W/(m · K)] and C p =900 [J/(kg · K)]) used in the earlier work of Patil et al. (2015a). Particles are heated-up due to the heat of adsorption, while the whole column also increases its temperature from room temperature. As long as the temperature of the particles and the background remains at certain difference, the contrast can be measured by the IR camera (up to 0.1 • C sensitivity). It takes a relatively long time to cool down the whole bed to its original state (room temperature, 18 • C), meanwhile only minor differences of the particle phase and the back wall can be observed, thus from the total of 120 s measurement, only the first 80 s data will be used for detailed analysis.

Experimental procedure
The described pseudo-2D fluidized bed was filled with spherical zeolite 13X beads of 1.8-2.0 mm diameter (supplier: SIGMA-ALDRICH ® , 8-12 mesh). These particles are fluidized by a gas entering uniformly from the bottom of the bed by a gas distributer fabricated from a porous metal plate. Initial fluidization was achieved by feeding pure nitrogen gas (N 2 ) at room temperature (18-19 • C) controlled by a mass flow controller (up to 120 l/min). Subsequently a mixture of nitrogen and carbon dioxide (N 2 /CO 2 ) with a known composition was fed to the bed, shown as Fig. 6.2. Due to the CO 2 adsorption, the particles heat up by the adsorption heat, where the bed temperature rises up to 40-80 • C above room temperature depending on the CO 2 concentration. Once the particles are fully saturated the heat production is terminated and consequently the bed will (eventually) cool down to room temperature. We have used a CO 2 concentration ranging from 20% to 100%. Feeding with pure CO 2 gas causes initial defluidization due to the rapid adsorption, whereas feeding with 20% of CO 2 does not lead to a sufficiently high bed temperature. After the particles became saturated, they were taken out from the column and were replaced by a fresh batch. All experiments with different CO 2 concentrations have been repeated three times.

Optical apparatus
Two cameras are used for recording the fluidization images. The visual and infrared cameras were connected to the computer system via a trigger box. The cameras are triggered to take images simultaneously according to a preselected signal, of which details can be found in Patil et al. (2015a) and Sutkar et al. (2015). The two images taken by IR and visual cameras can be combined to use the position of particles to correct the thermal image such that only IR radiation from the particle phase is included. This mapping technique will be discussed in the next section.

IR camera and calibration
The temperature of the particles in the fluidized bed was determined from infrared recordings obtained from a FLIR ® X8400sc high speed IR camera with a maximum resolution of 1280×1024 pixels. Using an MW 50 mm 2.0 HD lens without camera filter, pictures were taken at a frame rate of 40.0 Hz. The integration time of 541 µs has been set for all measurements, to enable accurate detection in the range of 5-300 • C with a maximum error of 1 • C. As output of the camera, the radiation intensity (noted as digital level (DL)) of each pixel is transferred to a 14 bit digital value. The infrared signal received by the detector of the camera has propagated through the window of the fluidized bed, the air and the lens of the camera. The influences of window and air on the DL should be reduced as much as possible. The absorption by the air is reduced by placing the camera relatively close to the set-up, so that the transmittance of the air can be assumed to be close to 1. As described before, we used a sapphire glass front window for its high transmittance Pishchik et al. (2009). The reflection of environment, on the other hand, is minimized by performing the experiment in a closed unit without interference of any warm objects. White LED lamps have been used to illuminate the set-up and enhance the signal for the visual camera. It does not interfere with the infrared part of the electromagnetic spectrum, as their peaks are in the blue and yellow part of the visual spectrum. However the lamps can get warm after a long measurement and provide an IR source of the reflection by the front window.
The IR camera measures the radiation from the surface of an object. The radiation is either absorbed, reflected or transmitted. The transmittance of the front window is close to 1, meaning that the major part of the detected radiation is thus emitted by particles IR counts [-] 4000 5000 6000 7000 8000 9000 10000 11000 12000 13000 14000  behind it. Although the infrared camera was calibrated using perfect black-body radiation and contains a correction for gray-body objects, developed by the manufacturer (FLIR ® ), using values for the emissivity of the objects and for the room temperature, these results did not reflect the actual temperature of the object present at the inside of the pseudo-2D fluidized bed. Both emissivity and shape of the object can affect the radiation emitted. We adopted a similar way of camera calibration as proposed by Patil et al. (2015a). A bed of zeolite 13X particles of 1.8-2.0 mm diameter was heated up to 200 • C in an external oven and transferred into the bed, where one zeolite 13X particle attached to a thermocouple was placed close to the front observation window. During the process of cooling down to room temperature, by conduction only, the temperature of the particle was measured with a thermometer. Simultaneously, the DL detected by the infrared camera was recorded. Multiple repeated calibration experiments according to the described procedure have led to the calibration curve presented in Fig. 6.3 in which the detected number of infrared counts is shown as a function of the correlated temperature [K]. The black line indicates the third order correlation function to estimate the temperature from the detected number of counts with a correlation factor R 2 =99.7%. T = 7.117 · 10 −11 · DL 3 − 1.997 · 10 −6 · DL 2 + 2.575 · 10 −2 · DL + 219.7 (6.18) in which T is the temperature [K] and DL stands for the digital level of the received infrared signal. Figure 6.4: Schematic of the general processing procedure used to determine the particle phase in the fluidized bed reactor. The used filter technique may differ and can include the usage of a post-processed visual image.

Visual camera
Digital images of the fluidized bed were recorded using a high-speed visual camera of LaVision ® (Imager pro HS) with a resolution of 560×1280 pixels. LED lamps illuminating the set-up were used to minimize reflections and to provide sufficient light to use short shutter times. An integration time of 300 µs was used, after which a 100 µs delay is used before a second recording was taken with the same integration time. These two consecutive recordings are required for the particle image velocimetry (PIV) analysis (this part of analysis is not included in the paper), as explained before, in which the time difference between the two frames (i.e. 400 µs) is used.

Images mapping and data processing Detect particles phase from the images
Difference of the resolution of the IR and visual cameras brings difficult to map two images taken instantaneously. A schematic overview of the image processing is improved, in which is elaborated the steps to subtract information from the recording. Because the visual and infrared cameras are both positioned at an angle to avoid reflections of either one of both. Images taken by both cameras contain larger area than the desired part, the fluidized bed. Transformation of the images is necessary. After the images are cropped, the infrared DL in the leftover pixels -the particle phase -are then converted to their respective temperatures, conform the calibration method as explained in the previous section.

Mapping of the images
A dummy image of 320×800 pixels, representing the front window of 8×20 cm, is used as a blueprint for both images. By selecting the corner points of the images of the front window from any frame recorded by both cameras, using the built in interface in MATLAB ® 2015b image processing toolbox, the image of interest can be stretched, rotated and cropped to fit the dummy matrix which contains all digit pixels. The raw IR and visual images will not be presented here for the sake of brevity. The processed visual and IR images and the Figure 6.5: A schematic illustration of processing and overlapping visual and IR images, by following the steps described in Fig. 6.4.
illustration of overlapping these two are shown in Fig. 6.5, which is also a supplementary explanation of Fig. 6.4. As can be seen in Fig. 6.5 step (3) and (4), some isolated particles are not overlapping well, because the images were recorded from different angles. The mentioned mask C in Fig. 6.5 is used to extract particle phase only from IR images which will be introduced in the next section.

Different filtering techniques
As a base case, we will first discuss the results when no filter is used to detect the particle phase. By simply taking the raw data of the counts, exporting this data to MATLAB ® , subtracting the background image and converting the detected number of counts to the respective temperatures, the background is also taken into account, leading to a bias in the obtained temperature. To efficiently correct for the bias three measures can be taken of which the outcome is presented in Fig. 6.6 in this section. Firstly, the region of interest (ROI) is defined by restricting the image to a region coinciding with the initial packed bed (before introduction of the gas mixture to the fluidized bed), which covers most of the dense phase. Hereby the relative numbers of pixels containing the background (not the gas phase) can be reduced.
Secondly, a threshold value for the infrared image can be used, determining whether a pixel belongs to the particle phase or the background, denoted as "IR". A threshold value of 0.32 (converted image intensity from 0-1, where 0 and 1 indicate the minimum and maximum, respectively) can be used to differentiate between the particles and the background. As reported by Patil et al. (2015a), this value can be found by trial and error. However, by using this method for each particle in the bed, a small halo of pixels with a higher temperature than the surroundings is observed. Once the cooling phase in the experiment is approached, the background will have a profound influence on the measured particle temperatures because of the reduced temperature difference.
As reported by Patil et al. (2015a), another method is to use a visual mask, obtained by using a threshold value from the visual image, denoted as "VIS-IR" (same concept as mask C mentioned in Fig. 6.5). However, due to the angle between the two cameras, overlaying of these images causes distortions in the cases of few isolated particles (above the free-board). This is related to 3D-effects and cannot be avoided unless the two cameras are placed at the same position, as shown Fig. 6.5. Fortunately the images capture a large part of the dense phase and loss of a few particles can be tolerated and therefore this mapping approach can still provide sufficiently accurate results.
Therefore, the pixels that are considered to belong to the particle phase passing both the visual and the infrared filter are then selected to be the pixels of interest. Although not all particles will be included using this technique, the results are an accurate representation of the particle phase, excluding halo effects in the infrared image.
From Fig. 6.6, it can be observed that direct application of a threshold on the IR images to subtract particle data will suffer from interference with the background. Both ROI and IR-VIS filter provide almost overlapping curves. However, during the heating stage of the experiments, vigorous fluctuations in the temperature curve obtained by applying ROI can be clearly noticed. It is a direct result of gas bubbles passing through the bed. The back plate becomes partly "visible" to the camera and hence its temperature will be recorded. During the final phase of the adsorption process, the particles cool down and the fluctuation in the temperature curve is reduced, due to smaller differences of the temperature of particles and background. Both ROI and IR-VIS can capture the average bed temperature evolution, however to obtain sufficient accuracy in the data for the instantaneous spatial distribution of the particle temperatures requires usage of the IR-VIS filter. More details of the results of image processing after applying this filtering method to determine the spatial temperature distribution can be found in Fig. 6.7.

Results and discussion
The experiments have been performed with different CO 2 concentrations in a CO 2 /N 2 gas mixture. The evolution of the averaged bed temperature is shown in Fig. 6.8 for different CO 2 concentrations in the feed, viz. 60%, 40% and 20% of CO 2 in the gas mixture. The  Figure 6.6: Average particle phase temperature < T p >, determined from the number of infrared counts after application of the mentioned filtering technique. Experiment using 60% CO 2 in the gas feed. evolution of the average particle temperature can be divided roughly in three stages. In the first stage, T p quickly increases due to fast adsorption; during the particle saturation phase, adsorption slows down and a short thermal equilibrium can be achieved when the amount of heat produced is equal the amount being removed. In the final (third) stage, the bed cools down by the gas. We have performed experiments with pure CO 2 as well, however a sudden defluidization occurs when the gas is switched from pure N 2 to CO 2 , which is caused by the very fast initial adsorption. Our CFD-DEM is not accounting for any CO 2 supply limitation and therefore it cannot capture the defluidization happening in the experiments with pure CO 2 . When feeding a CO 2 /N 2 mixture we still assume that there is no CO 2 supply limitation due to fast mixing in the relatively small bed and constant CO 2 partial pressure on the particles surface. However when the concentration of CO 2 is lower than 40%, for instance at 20% as shown in the graph, the largest temperature rise in the bed is only around 10 • C which causes larger experimental errors due to the increased disturbance of the background. For simplicity in the following discussion we will focus mostly on the cases with 40% and 60% CO 2 . The simulation conditions and parameters used in the adsorption model based on A are listed in the

IR-VIS filter
(f) t=30 s, using an infrared-visual mask threshold as filter (VIS-IR) Figure 6.7: Temperature distribution of the particle phase after using different filtering techniques, 30 seconds after addition of CO 2 in a gas feed of 40% CO 2 /60% N 2 . Figure 6.8: Spatially averaged particle temperature as a function of time measured with the infrared camera. The experiments were performed with a CO 2 concentration in the feed equal to 60%, 40%, and 20%.

Evolution of the spatially averaged particle temperature Heat loss to surroundings
In our experiments, the heat is produced by CO 2 adsorption in the zeolite 13X particles. The adsorption heat is removed from the particle phase by convective gas-particle heat transfer. Due to the relatively small volumetric heat capacity of the gas phase (ρ g C p,g ) and the high gas-particle heat transfer rate in the fluidized bed, the temperature difference between the gas and particles is small. The produced heat due to adsorption is transferred to the gas, and removed via the outlet or lost to the environment via the surrounding walls. The later contribution can be significant and should be taken into consideration in the CFD-DEM. Before detailed CFD-DEM simulations were performed, a CSTR model was developed to provide a simple tool (learning model) to study the process. The energy balance for both the gas and particle phase as well as the walls are all considered. The heat loss term has been incorporated in the energy equation of the gas phase. More details of the model can be found in C. We have fitted the kinetic model in non-isothermal condition from the TGA data, then correlated the contribution of the temperature change during adsorption with an Arrhenius equation (details are shown in A). Combined with the fitted gas-to-wall heat transfer coefficient, the CSTR reaches a good agreement with the experimental data regarding to the temperature evolution, as shown in Fig. 6.9. In the graph the black lines indicate the IR measured averaged particle temperature, and blue lines with different dot- dash styles represent CSTR fitting results. From the CSTR model, it is estimated that the gas-to-wall heat transfer coefficient is in the range of 200-300 [W/(m 2 · K)]. This range is in accordance with the general model proposed by Kunii and Levenspiel (1991). It shows that the maximum heat transfer coefficient of a bubbling fluidized bed at a heat-exchange surface is affected by the size of the particles. With the given particle size of our system, a suggested h max of 200-250 [W/(m 2 · K)] can be estimated from Figure 8 of Kunii and Levenspiel (1991). The heat liberation rate due to CO 2 adsorption as well as the increase in particle weight have been accounted for in our  Figure 6.9: Results of the average particle temperature evolution comparing experimental and the CSTR simulation data. For the CSTR model data results for different values for the heat transfer coefficient from the wall to the surroundings are shown. CFD-DEM model. During the whole process, particles are heated-up by adsorption, where the heat is transferred to the gas phase by convection. Heat transfer by particle-particle conduction is neglected due to the rapid particle-particle collision under fluidization conditions. The Biot number of the zeolite 13X particle used in our system is smaller than 0.1, consequently there is no prominent internal temperature gradient. The bed temperature at the highest point is approximately equal to 20-25 • C above ambient temperature, therefore the heat transfer from particle to surroundings by radiation can be also ignored. In summary, in our CFD-DEM we applied the same thermal model as in the CSTR model, excepted the heat loss to the environment.
This heat loss is accounted for in the CFD-DEM by a boundary condition imposed on the gas phase thermal energy equations. Within CFD-DEM the grid size typically exceeds the thickness of the thermal boundary layer and the unresolved heat loss should be fitted through an effective gas-to-wall heat transfer coefficient h w which can be related to the effective heat conductivity of the gas and the near-wall gas film thickness. The thickness of this thin gas film was suggested by Patil et al. (2006) to amount one or two particle diameters. The boundary condition used in the CFD-DEM is: with: Applying linear interpolation, this proposed boundary condition can be discretized and leads to the formula to calculate gas temperature at the boundary cell (T g,0 ): So far, the heat loss to the wall is handled in largely the same manner as in the work of Patil et al. (2015a). However, we found that the effective gas phase conductivity, which is calculated from equation 6.11 (≈ 0.01 W/(m · K)), in combination with a thermal boundary layer thickness δ h ≈ 10 −3 leads to a much too low value for h w (10 W/(m 2 · K)) in comparison with the fitted value on basis of the experiments. The reason is that only the gas phase contribution is accounted for in equation 6.11.
An estimate of the effective conductivity of the bed can be obtained from the Zehner-Schlünder model for stagnant thermal conductivity, which has been used in several studies (Kuipers et al., 1992;Armstrong et al., 2010). The rate of k eff g /(1−(1−ε) 0.5 )k g for our system (CO 2 and N 2 gas mixture and zeolite 13X particle, with k g =0.02 [W/(m · K)] and ε=0.43) is around 20 which leads to h w =200 [W/(m 2 · K)] in agreement with the experimentally fitted value. Moreover, when the bed is fluidized, the contacting surface of the gas and particles are continuously renewed by random particle fluctuations. This kinetic contribution can contribute to the effective thermal conductivity of the bed. We can therefore anticipate that the ratio of k eff g /(1 − (1 − ε) 0.5 )k g is larger than 20 for fluidization conditions. Hence, we can conclude that this approximation is reasonable. However, more accurate correlations for the effective thermal conductivity for fluidization conditions should be developed by fully resolved computations using DNS. .10: Results of the spatial averaged particle temperature evolution feeding with 60% of CO 2 comparing with experimental and simulation results. (The CSTR results are shown with blue lines at adiabatic condition as well as with heat loss to surroundings applying a heat transfer coefficient of 300 [W/(m 2 · K)]. The CFD-DEM simulation results presented in the graph are with anticipated effective thermal conductivity (k ef f g ), which is shown as a result of 5, 50 and 100 times of the value calculated from equation 6.11.)

Comparison of experimental data and CFD-DEM
In Fig. 6.10 we show the average bed temperature obtained from the experiments and CFD-DEM simulation results together with the predictions obtained from the CSTR model. In this figure we have included the results obtained for adiabatic operation. The experimental data are indicated by black solid lines, the CSTR model and CFD-DEM results are shown with blue lines and red symbols, respectively.
As evident from Fig. 6.10, a good agreement between CSTR and CFD-DEM simulation results is obtained in this case. In this figure we also show the CFD-DEM simulation results using the boundary condition discussed in the previous section in combination with the corrected effective thermal conductivity.
The simulated temperature curves closely follow the trend observed in the experiments. Initially (first 0-8 s), the temperature of the particle phase increases rapidly. During this period the energy which is used to heat-up the particles is produced by fast adsorption. At a certain moment however the adsorption rate reduces and the effect of cooling of the gas becomes more pronounced. A short equilibrium of heat production and heat removal can be observed after around 8-11 s (roughly speaking). Finally, the cooling effect by the gas phase plays a dominant role and the bed temperature declines gradually. Once the same wall temperature and gas-to-wall heat transfer coefficient applied in the CFD-DEM as well  Figure 6.11: Results of the spatial averaged particle temperature evolution (feeding with 40% of CO 2 ) comparing with experimental and simulation results. (The CSTR results are shown with blue lines at adiabatic condition as well as with heat loss to surroundings applying heat transfer coefficients of 250 and 300 [W/(m 2 · K)]. The CFD-DEM simulation results presented in the graph are with anticipated effective thermal conductivity (k ef f g ), which is shown as a result of 50 times of the value calculated from equation 6.11.) as the CSTR model (T w = 300 K and h w = 300 [W/(m 2 · K)], blue solid line in the graph), the CFD-DEM simulation results show good agreement with the experimental data, provided that the effective conductivity in equation 6.11 is given a value of 50-100 times the original value. Apparently, according to the CSTR model, the heat transfer resistance should not reside in the bulk of the gas phase but in the gas-wall boundary. With increasing effective conductivity the resistance in the gas phase reduces.
These two parameters used in CFD-DEM (h w and k eff g ) are directly linked to the thickness of the thermal boundary layer. So far, there are no empirical correlations that can be directly applied to CFD-DEM to estimate these parameters. The effective conductivity of the gas phase in fluidization condition requires further investigation, e.g. by more detailed models using direct numerical simulation (DNS).
The CSTR and CFD-DEM model also reach a good agreement with experimental data when the bed was fluidized by feeding a gas mixture containing 40% of CO 2 . For brevity, the data from the CSTR model with adiabatic condition and with gas-to-wall heat transfer coefficient equal to 250 and 300 [W/(m 2 · K)] are shown in Fig. 6.11. Again, the simulation results of CFD-DEM matches with experimental data provided that the corrected effective thermal conductivity of the gas has been applied.

Spatial particle temperature distribution
We have used the CSTR model as a learning model to facilitate interpretation of the experimental data. However, the limitation is that only the global data can be obtained from the CSTR model. More details, especially the spatial distribution of the properties in the fluidized bed, such as the particle temperature distribution are of interest as well. Building on the successful simulation of the particle temperature evolution, the spatial temperature distribution will be discussed in this section. First the snapshots obtained from IR images and CFD-DEM simulations are compared.

Snapshots
The IR camera can generate high resolution images of the particles in the pseudo-2D fluidized bed. On average, each particle occupies 9-10 pixels, therefore the spherical shape of the particles can be well detected from the IR images. This feature is convenient to permit a direct comparison with the CFD-DEM thermal data of the particle phase. In Fig. 6.12, the IR images and snapshots of CFD-DEM simulated particle temperature field are shown at the same moment, top and bottom, respectively. Note that the relative time started after the fluidization gas was switched from pure N 2 to CO 2 and N 2 mixture. Prior to CO 2 introduction into the bed, there was no temperature difference between the particles and the column, therefore the particles can hardly be observed. The shown snapshots are at 1.0, 1.5, 4.0, 9.0, 11.0 and 21.0 s after CO 2 introduction into the bed. The first three snapshots represent the phase that the bed heats up. At 9.0-10.0 s, particle temperature reaches its maximum, and in the last snapshot, (at 21 s) the saturated zeolite particles cool down by the gas phase.
Generally speaking, the snapshots of the CFD-DEM simulations exhibit a similar pattern as the temperature field obtained from experiments. Both images reveal the same cooling area that is found at the bottom of the bed, where the fresh cold gas enters the bed. However the "fingering" of cold particles observed in the CFD-DEM simulation can hardly be seen in the snapshots of IR images. It is speculated that when the bed temperature increases by adsorption, the front window is heated along with the bulk of the bed. The heated front window transfers its heat to the cold stream of the particles through the gas phase, hereby the "fingering-effect" is not observed in the experimental IR images. We have computed the front window temperature resulting from the CSTR model (see in D) and from Fig. D.7 it can be seen that the temperature of the sapphire window increases at a rate of 0.5 K/s. The recorded IR radiation by the camera from the bed certainly contains the contribution from the front window during the measurement, however it should be minor because of the low emittance of sapphire in the middle range of the IR spectrum. Sapphire has a high transmittance up to 82% and its emittance is reported in the range of 0-10% at 300 K Pishchik et al. (2009). Figure 6.12: Snapshots of the images taken by the IR camera (top) and the CFD-DEM simulation results (bottom) at different moments after introducing CO 2 (60% in the gas feed) into the system. From left to right the corresponding snapshots are taken at 1, 1.5, 4.0, 9.0, 11.0 and 21.0 s, covering the initial heating-up, equilibrium and cooling down three stages of the measurement.

Probability density distribution
Following the qualitative observation of the spatial distribution of the particle temperature from the snapshots, a quantitative analysis will be presented by examining the PDF of the particle temperature at different moments. In Fig. 6.13, the instantaneous particle temperature PDF are shown together with CFD-DEM simulation results. Experimental and simulation results are shown for 40% CO 2 (left) and 60% CO 2 (right) in the figures. The distributions of particle temperature with 60% CO 2 correspond to Fig. 6.12 at time 4, 11 and 21 s. The experimental PDF presented here is post-processed after application of a mask obtained from the visual image (see section 3.4.3). In comparison to the CFD-DEM results, the experimental data include a slightly broader temperature range, both at low and high temperatures. The peak in the PDF corresponds to the bed average temperature. Overall the agreement between the experimental and simulated PDF's is very good. The tails observed in the experimental PDF is most likely due to noise and ineluctable reflection from the environment. Figure 6.13: Particle temperature PDF obtained from CFD-DEM (black line) and IR images after applying IR-VIS filter mentioned in previous section. The left and right column of graphs are corresponding with CO 2 concentration of 40% and 60%. From top to bottom rows the graphs are shown the PDF obtained at 4 s, 11 s and 21 s, respectively.
To estimate the effect of the boundaries in the depth (y-) direction on the PDF of the particles near the front wall (i.e. those particles present in the front layer of the cells) and the PDF of particles in the model (i.e. those particles presented in the third layer of cells) are compared to the original PDF (using all particles) in Fig. 6.14. From the figure it can be seen that the particles presented in the front layer have a relatively low mean temperature compared to those present in the middle. The PDF peak of the front layer locates slightly to the left compared to the peak of middle layer. The differences between the PDFs are not large during the initial stage of the absorption process (t=4 s). Clear differences can be found when the bed reaches its highest temperature (approximately, at t=11 s). Moreover, Figure 6.14: Particle temperature PDF obtained from different layers of the simulation domain in y-axis direction at different moments. Left and right are at 4 s and 11 s, respectively. the long tail of the PDF of all particles in the low temperature region partly overlaps with the one of the middle layer, meanwhile the PDF of particles present in the front layer is more symmetrical. The long tails of the PDF signify the "fingering-effect" observed in Fig.  6.12 in CFD-DEM simulations. Due to the low thermal capacity of the gas phase (ρ g C p,g ), the cooling zone of the fluidized bed is confined to the bottom of the bed, where the fresh gas is entering. Particles can cool down in this zone and are subsequently carried by the bubbles towards to the top of bed.

Summary and conclusions
In this work, a CFD-DEM was used to simulate heat transfer in a pseudo-2D fluidized bed with heat production in the particle phase. In order to mimic the heat production of an exothermic reaction, a model system is used, in which heat is released by adsorption of CO 2 in zeolite 13X particles (1.8-2.0 mm, Galdart D type). This system was studied experimentally using a non-intrusive optical technique combining an IR and a visual camera. The temperature field is captured through IR thermography, simultaneously the particle position is recorded by the visual camera. Data of interest processed from the IR measurements are accurately extracted by applying filters obtained from visual images. A simple CSTR model is used as a learning model for the purpose of fitting the gas-towall heat transfer coefficient (h w ) in order to match the experimental data. In the CSTR model, the adsorption kinetic model and heat loss (through the gas phase) are introduced. Incorporated with the estimated h w from the CSTR by applying it on the boundary condition, the CFD-DEM simulations are subsequently conducted and compared with respect to the experimental data as well. The comparison shows good agreement when the effective gas thermal conductivity is much larger than its original value, which was estimated by only considering the contribution from the gas phase. However, the accurate value of the effective thermal conductivity of the gas at fluidization conditions is unknown. Therefore, for a further application of the CFD-DEM on investigating the heat transfer in a fluidized bed at non-adiabatic boundary condition, an accurate value or correlation for the effective thermal conductivity of the gas is critical to lead to successful simulations, which can be obtained, for example by fully resolved computations using DNS.
Based on the matching results of average particle temperature and temperature evolution profiles of experiments and simulations, quantitative comparison of the particle temperature distribution is presented, by examining the PDF obtained from both approaches. Overall, a good agreement between the simulations and experiments is achieved with respect to the spatial temperature distribution. However, the "fingering-effect" found in the CFD-DEM simulations was not clearly captured by the IR measurements. This can be attributed to the noise in the measurement data and interference of the front window during the heat-up stage. Furthermore, the effect of the boundary condition in the depth (y-) direction on the particle temperature was also discussed, and it was found that particles present in the front layer appear to have a (slightly) lower mean temperature compared to those present in the middle of the bed, and their temperature distribution is more symmetrical compared to all particles. Thus the differences between the PDF of all particles and those present in the font layer shows that the IR measurement can capture the temperature field of the whole bed rather than only the front layer.

Conclusions
In this thesis, the heat transfer phenomena in a fluidized bed with heat production have been studied using an in-house developed CFD-DEM and an experimental approach. The project is part of a large research topic which aims to establish a complete model in order to obtain an accurate description of gas-phase polymerization processes under actual operation conductions, including all the important features encountered in these processes, such as the hydrodynamics of the multi-phase flow, heat/mass transfer, heat production, liquid injection and so on. The olefin polymerization is an extremely exothermic process, and therefore a large amount of heat needs to be removed from the reactor efficiently to prevent hot-spot formation, consequently the heat management of the reactor is crucial for a successful operation. Understanding and modeling of gas-phase polymerization is particularly challenging because phenomena prevailing at different length and time scales need to be properly accounted for. This thesis is focused on a relatively small lab system and investigated the particle temperature distribution in a fluidized bed with heat production and the influence of operating conditions on this distribution. In addition, the mechanism of hot-spot formation was studied.
To mimic the polymerization heat liberation, a constant volumetric heat production rate [W/m 3 ] was implemented in the particle energy equation. The effect of various operating conditions on the heat transfer has been systematically studied: viz. superficial gas velocity (Chapter 2), elevated pressure (Chapter 3), bed dimension (Chapter 5), particle size distribution and heat source term (Chapter 4). In addition, the model was validated with experiments (Chapter 6). Simulation and experimental results show that the CFD-DEM is an effective approach to study the temporal evolution and spatial distribution of the particle temperature in a small lab scale fluidized bed.
The simulations were first performed in a pseudo-2D fluidized bed in order to conduct the parameter study and to obtain a general understanding of the hydrodynamics and thermal aspects of the fluidized bed associated with this particular physical model. Later the simulations were performed in 3D rectangular domains by step-wise increasing the bed depth. From the pseudo/2D bed simulations we can conclude that the time-and spatialaveraged particle convective heat transfer coefficient, which was calculated by Gunn s correlation, is independent of the superficial gas velocity (Chapter 2) but increases with increasing pressure (Chapter 3). The elevated pressure influences the hydrodynamics and significantly enhances the heat transfer performance of a fluidized bed. Operating the bed at elevated pressure helps achieving a more isothermal condition in the fluidized bed, thereby reducing the opportunity for hot spots to develop. The conclusions made from the pseudo-2D bed simulations were confirmed by the 3D bed simulations as it was found that the particle temperature and distribution are hardly affected by the increased bed depth. Thus, to study the heat transfer behavior in a fluidized bed at adiabatic conditions, it is effective to conduct the simulations using a pseudo-2D domain reducing the computational cost and compare with the experiments using non-intrusive analysis techniques.
Additional complexity is introduced to the CFD-DEM model by implementing a particle size distribution. The simulations have been carried out in beds with Gaussian PSD's with different distribution widths but with the same Sauter mean diameter (d 32 =1.2 mm). Different heat production mechanisms, namely a constant volumetric heat production ([W/m 3 ]) or constant particle heat production ([W]) are used to mimic the heat production obtained from a normal catalytic reaction and polymerization reaction, respectively. Compared to simulations performed with a mono-dispersed particle diameter, much wider temperature distributions were found in the simulations with a wide Gaussian distribution. It was found that the hot particles that are close to the melting point are small particles with high catalyst activity, and they are mostly found in the free-board and near the side walls. However, this study is limited to a small size fluidized bed and a pseudo-2D geometry. Extrapolation of the findings to large scale industrial situations should be done with great care. Nevertheless, this work can help to improve our understanding of hot-spot formation in fluidized bed polymerization reactors. The next step was to validate the CFD-DEM model extended with heat production in the particle phase with experiments. Due to the obvious technical limitations we did not directly perform the actual olefin polymerization with a pseudo-2D fluidized bed. In order to mimic the heat effect due to the exothermic polymerization reaction in the pseudo-2D fluidized bed reactor, a model system was used. In this model system, heat is released in zeolite 13X particles (1.8-2.0 mm) due the adsorption of CO 2 . Characteristics of the adsorption kinetics, isotherm and reaction enthalpy have been achieved by performing Thermogravimetric Analysis (TGA) and Simultaneous Thermal Analysis (STA). Model validation was carried out with non-intrusive experiments in a pseudo-2D set-up combining a high speed infrared (IR) camera and a visual camera in order to record the particle positions and particle temperature simultaneously. The combined technique provides insightful information on the particle temperature distribution for different CO 2 concentrations. The comparison shows good agreement when the effective gas thermal conductivity is properly estimated. However, an accurate value of the effective thermal conductivity of the gas under fluidization conditions is unknown. Therefore, for a further application of the CFD-DEM on investigating the heat transfer in a fluidized bed with non-adiabatic boundary conditions, an accurate value or correlation for the effective thermal conductivity of the gas is critical to lead to successful simulations, which can be obtained, for example, by fully resolved computations using DNS.

Discussion and outlook
Simulations have been performed with an in-house developed CFD-DEM in order to study the heat transfer in fluidized bed with extensive heat liberation. The investigation covers several aspects, including the effect of operating conditions, bed geometry, particle properties and heat source terms. The model validation has been successfully conducted in a pseudo-2D bed with an optical technique. However, due to limitations associated with computational cost, this discrete model can only be used effectively to simulate a relatively small system (up to 1 million particles). Hence here in this section, we are going to discuss how to improve or extend the existing model in order to facilitate future investigations on this specific topic. Based on the existing model, the following aspects can be considered for further studies: • We have extended the particle model by considering a Gaussian distribution with different widths. However, in practice, the particle size distribution obtained from commercial reactors are different for different catalyst types, and moreover, it is related to the particle residence time. Gas phase olefin polymerization is a typical coordination polymerization process, which is composed of (multiple) active sites of Ziegler-Natta catalyst. The complexity of the coordination polymerization is due to the multi-step reaction chains: i.e. the formation of active sites, initiation of active sites, propagation, chain transfer to monomer, transfer to hydrogen, transfer to co-catalysts, and deactivation reactions . It still remains a challenge to develop a complete kinetic model and incorporate this in CFD models (Zhang and Luo, 2014;Yao et al., 2015;Schneiderbauer et al., 2015;Hamzehei, 2011). The simple formulation of the heat source terms used in our simulations may require further refinement for specific reaction systems. Thus, the future success on conducting more realistic simulations through CFD-DEM is of interest to incorporate the PSD obtained from industrial catalyst particle growth models, as well as to implement a kinetic model to derive the corresponding heat production.
• The soft sphere approach has been utilized to compute the particle-particle and particlewall collisions in the CFD-DEM. However, the interaction model does not account for particle breakage or agglomeration after collisions. In real system, particles can agglomerate and swarm into a big chunk when the particles are overheated or charged, while on the other hand they can also break into small fragments due to collisions or other external forces. From our study, we have found that the small but high activity particles have the possibility to develop to hot spots. They can be found above the free-board and in the vicinity of the side walls. The former is caused by minor vertical segregation and the latter was found to be attributed to the solids movement. Moreover, due to the static force, particles are expected to be charged during the contact with the reactor wall (Sowinski et al., 2012;Moughrabiah et al., 2012). The static force will further cause segregation that may lead to the development of hot spots at the reactor walls, where the heat transfer coefficient was found to have a smaller value than in the center of the bed (Chapter 4). Therefore it is interesting in the future models to incorporate the static force and also consider particle agglomeration and breakage into the particle interaction model of the CFD-DEM.
• When the olefin polymerization is operated in condensed mode, a major part of reaction heat is removed by evaporating injected liquid. The liquids are injected into the reactor, evaporate in the interior of the fluidized bed, and the generated gases undergo the respective catalytic reactions (Bruhns and Werther, 2005). The existence of the liquid phase affects the solids movement significantly due to the solid-liquid contact. Therefore, it is of importance to monitor the solids flux when there is a liquid phase in the fluidized bed. However, the complexity of simulating the liquid phase is greatly enhanced due to change of liquid configuration, which can be briefly split into three regions, viz. jetting region, fragmentation and formation of ligaments, and the third, breakdown of ligaments to drops (Bruhns and Werther, 2005). Moreover, besides the exiting of "free liquid" in the bed, the liquid can also be attached to the particles and form a "bridge" between particles or form a coating layer covering the surface before it is evaporated (Heinrich et al., 2003). Therefore, it is also of importance to extend the existing CFD-DEM with liquid injection, and to monitor solids flux movement and associated thermal effects.
• As discussed before, the CFD-DEM is limited by the fact that it can only be applied for the simulation of a relative a small lab scale system with particle numbers up to 1 million. To simulate a bench scale or even a larger industrial scale fluidized bed, coarser models which make use of more assumptions are required. The two fluid model (TFM) in which both the gas and particle phase are computed as continuous phases is able to simulate much larger systems. Therefore, it is possible for such a type of model to make a direct comparison to a bench scale and even an industrial scale fluidized bed (Chang et al., 2013;Yao et al., 2015). The CFD-DEM, which was applied as a learning model, has provided insightful guidance for future investigations on heat transfer and hot-spot formation in gas-phase olefin polymerization reactors. It was found that the bed with a wide PSD can easier develop hot spots. Moreover, the small and high activity particles are found above the free board and in the vicinity of the walls. Correspondingly, it is hereby of importance to introduce a more realistic PSD to the continuous type of model to lead further investigations on this issue. The TFM that use statistic models such as the population balance model and polymerization kinetics have been used to predict the PSD and molecular weight of polymer product, and the results have been compared with experimental data from industry directly (Chang et al., 2013;Yao et al., 2015). A further extended TFM with a heat/mass transfer model, polymerization kinetics as well as a more realistic PSD is recommended to further investigate hot-spot formation.

Particle size segregation
Particles with size distribution or density differences may exhibit segregation in fluidized beds. Size segregation in a fluidized bed with a continuous PSD is not as pronounced as in a bed with a binary size distribution. It is found that after fluidization, the local PSD in a fluidized bed with continous size distribution preserves a simular distribution pattern (Dahl et al., 2003). In this appendix the discussion will be focused on the possible occurrence of size segregation in the fluidized bed. The spatial and time-averaged ratio of the local particle radius and the averaged radius of the all particles can provide detailed information on this size segregation. Moreover, the occupation of the different sized particles in the fluidized bed are also discussed.
The simulation results of particle size related information with narrow (σ1) and broad particle PSD (σ3) are shown in Fig. C.1 and Fig. C.2, respectively. To evaluate the local size distribution, the bed domain is divided into equally sized cells (the cell size is 4 times larger than the Eulerian grids size, which allows a maximum of 10 3 particles in a cell), and subsequently the local average radius and the standard deviation of the radius of particles are calculated in each cell and time-averaged from 1 s to 5 s simulation time using a time step of 0.05 s.
The time-averaged local ratio of particle radius and the mean radius of all particles (R p,i /< R p >) is shown in a1 and a2 in Fig. C.1, indicates clearly a uniform distribution of the particle size in the bulk of the fluidized beds. Values of this ratio below 1 (slightly) can be found at the top of the bed. The local standard deviation of the particle size distribution Figure C.1: Spatial distribution of dimensionless particle size in fluidized beds with narrow size distribution (σ1) at different gas superficial velocity. Graphs noted with a are the time-averaged local ratio of particle radius and mean radius of all particles (R p,i /< R p >), whereas the graphs b are the ratio of local particle size deviation and size deviation of the whole bed (σ i /< σ >). Graphs c and d show the occupancy of small and large particles in the domain, respectively. Note the arab numbers 1 and 2 behind a, b, c, and d are indicating different superficial velocities, namely, 2u mf and 3u mf .
from graphs b1 and b2 in Fig. C.1 suggests that the largest size difference can be found in the free-board. Meanwhile, larger bed expansion can be observed when the bed is operated at 3u mf , however the minor segregation of small particles on the top of the bed cannot be prevented by operating with higher u 0 . Similar conclusions can be drawn from graphs a1, a2 and b1, b2 in Fig. C.2 for particles with a broader size distribution (σ3). Combined with the previous discussion that most of Figure C.2: Spatial distribution of dimensionless particle size in fluidized beds with broad size distribution (σ3) at different gas superficial velocity. Graphs noted with a are the time-averaged local ratio of particle radius and mean radius of all particles (R p,i /< R p >), whereas the graphs b are the ratio of local particle size deviation and size deviation of the whole bed (σ i /< σ >). Graphs c and d show the occupancy of small and large particles in the domain, respectively. Note the Arab numbers 1 and 2 behind a, b, c, and d are indicating different superficial velocities, namely, 2u mf and 3u mf .
the particles with highest temperature are found in the free-board and in the vicinity of the side walls, the formation of these overheated particles residing near the free-board can be attributed to the minor size segregation in the vertical direction of the bed. However, there is no evidence for lateral segregation as evident from all the graphs (a). To provide more information on analyzing the distribution of the hot particles, it is also interesting to investigate the occupancy of differently sized particles in the fluidized beds. The occupancy of the particles in a cell is recorded when they are found to reside in that particular cell of the domain during the simulation time. The definition of the occupancy is shown in equation (C.1).
Before detailing the occupancy of particles, it is necessary to define "small" and "large" particles in a continuous PSD. Two limits can be indicated in the cumulative PSD curves shown in Fig. 4.1, one is at (bottom) 10% of the mass fraction, and the other is at (top) 90%. All particles with a size below the lower limit are defined as small particles, whereas all particles with a size exceeding the upper limit are defined as large particles. From the c and d graphs of Fig. C.1 and Fig. C.2, the occupancy of both small and large particles reveals that the particles are residing by preference near the two side walls of the bed. This behavior is logical because of the fact that the gas phase passes through the center of the bed as bubbles. When increasing the superficial velocity, i.e at 3u mf , larger bubbles are expected to be formed in the bed compared to the bed at 2u mf (Zhuang et al., 2014), hereby the bed is expanded and correspondingly particles can occupy a larger area in the bed, but still mainly reside at the left and right sides of the bed.

Supplementary of chapter 6
Characterization of the heat source Zeolites are natural or synthetic crystalline aluminosilicates with a structure of repeating pore networks, which can capture molecules (only molecules of certain sizes and shapes that can pass through). The 13X type, among all the natural types of zeolites, is reported to have the largest physical adsorption capacity of CO 2 at room temperature and atmospheric conditions (Yu, 2012). In our work, the heat of adsorption of CO 2 on zeolite 13X particles is utilized as a source term of heat production from the particle phase during fluidization. The general adsorption heat production can be described by the equation shown below: Q ads = ∆H ads k ads (q E − q) n (D. 1) This equation indicates all the critical variables and unknown parameters required for a quantitative description of the heat production. The adsorption kinetics (k ads (q E − q) n ) is established by fitting data from the Thermogravimetric Analysis (TGA), whereas the adsorption enthalpy (∆H ads ) is obtained by performing a Simultaneous Thermal Analysis (STA).

Adsorption kinetics and equilibrium
This process of CO 2 adsorption on zeolite 13X has been extensively studied in the literature (Siriwardane et al., 2005;Lee et al., 2002;Mérel et al., 2006;Cavenati et al., 2004). Nevertheless, the resulting descriptions of the adsorption kinetics vary over a wide range and Figure D.1: Weight change of zeolite 13X particles due to adsorption (increase of mass) and desorption (decrease of mass) of CO 2 /N 2 (60%/40%) over time.
therefore are applied on different approaches. Significantly differing values for the equilibrium adsorption of CO 2 on zeolite 13X are found in the literature (Siriwardane et al., 2005;Hefti et al., 2015). Therefore, a detailed analysis of the adsorption process is crucial to lead a success of the modeling work.

TGA measurement
The adsorption experiments performed in the TGA have been carried out with a CO 2 /N 2 mixture with different CO 2 concentrations. The data obtained from the TGA are with CO 2 concentrations of 20%, 40%, 60%, 80% and 100%. A fixed gas inflow (480 ml/min) is introduced to the sample basket which has been placed in a closed chamber. Particles were pretreated by placing them in an external oven overnight at a temperature of 400 • C. The entire procedure is as follows: a CO 2 /N 2 gas mixture inflow at room temperature (18-19 • C) is introduced to the sample basket, and subsequently the adsorption starts. Rapidly the particles are saturated, then the inflow gas is switched to pure N 2 at a temperature of 450 • C in order to lead a fast desorption of CO 2 , until the weight of the sample reaches a stable value. Subsequently the TGA temperature is decreased to the room temperature. In Fig.  D.1 a typical profile (raw data) for two cycles of TGA experiments are shown using a gas mixture with 60% CO 2 mixed with N 2 . The weight of the sample experiences a change up to 15% as a result of adsorption (increase) and desorption (decrease).

Pseudo-n th order fitting
Among all the different fitting methods to establish an adsorption kinetic model, which are summarized in the literature (Azizian and Fallah, 2010), the general n th order kinetics model k ads (q E − q) n is chosen to fit the TGA data in our work. Fig. D.2 presents the data of one cycle of measurement, whereas Fig. D.3 shows a fitting results by applying the pseudo n th order (PnO) fitting method proposed by Azizian and Fallah (2010). A fairly good agreement between the fitting curve and TGA data is achieved. Nevertheless, the temperature influence on the overall kinetic model was not identified during non-isothermal TGA measurements. The estimated rate constant (k ads = 1.342 [L n−1 mol 1−n s −1 ]) is an averaged value over time. Influence of the temperature on the rate constant cannot be ignored, since the temperature increases 15-20 • C rapidly during fluidization experiments, as shown in Fig. 6.8. Therefore, an Arrhenius equation (D.2) is proposed in order to account for the rapid temperature change during adsorption: Here k 0 is the infinite adsorption constant and E A is the isosteric heat of adsorption at zero loading which is equal to -29.3 ×10 3 J/kg. The values of these variables are obtained from literature (Sumida et al., 2012). Different k 0 value have been used in the Arrhenius equation to estimate k ads as a function of particle temperature during adsorption. The results shown in Fig. D.4 indicate that the best approximation value of k 0 is equal to 1.5 ×10 −5 . The mean of k is consistent with k ads fitted from TGA data.
The adsorption equilibrium have also been measured by TGA for different CO 2 concentrations, ranging from 20% to 100%. In Fig. D.5, the mass of CO 2 adsorbed by sample particles as a function of time is shown with data obtained from two repeated measurements. The concentration of 80% CO 2 in the gas inlet were left out for the sake of clarity. Increased adsorption equilibrium at increased partial pressure of CO 2 can be clearly noticed.
It can also be noticed that there are differences between the repeated and the first batch of measurements. The variation is attributed to the different exposure time of the sample particles to the environment during tests. The feature of containing an extremely large surface area makes zeolites response sensitively to the ambient conditions, i.e. temperature,  humidity, etc. The value of the adsorption equilibrium q E used in the model are averaged over the two repeated measurements, which are equal to 0.13 [g/g] and 0.15 [g/g], in a CO 2 gas concentration of 40% and 60%, respectively.

Adsorption enthalpy measured by STA
Zeolite 13X particles of 1.8-2.0 mm diameter are placed in a crucible to perform a simultaneous thermal analysis in the Netzsch STA 449F1 analyzer. The sample particles are first pre-heated up to 200 • C to indicate the time, and subsequently they are cooled down to the selected temperature. A CO 2 gas flow is introduced into the system when thermal stability has been achieved. Adsorption starts at the moment when CO 2 gas reaches the crucible, and after the adsorption finishes the sample is flushed with N 2 gas at 200 • C. Flushing the sample with high temperature N 2 gas leads to a fast desorption and therefore the sample particles can be regenerated and prepared for a new experiment of CO 2 adsorption. The heat of adsorption is measured at three different temperatures: 30 • C, 50 • C and 100 • C. The same procedure has been applied on a blank test with an empty crucible to correct for the interference due to switching of the inlet gases. The mass and power deviations due to adsorption of CO 2 on zeolite 13X after correction with the blank experiment have been measured four times at the selected temperatures of 30 • C, 50 • C and 100 • C. For the sake of brevity, the STA data obtained by following the description introduced above is shown in shows the mass change of the sample as a function of time, meanwhile the bottom one shows the temperature and the amount of energy required to remove from the system to keep the system isothermal. Moreover, the integrated area of the energy generated during adsorption is marked with gray colour blocks. The resulting value for the energy release is normalized by the adsorbed mass of CO 2 . An overview of the obtained values for the adsorption enthalpy at different temperatures is shown in Table D

Energy balances
In the CSTR model, both gas and particle phases are assumed to be well-mixed.The thermal energy balances reduce for the CSTR model to: Particle phase: Gas phase: Walls: by which the particle phase temperature T p , the gas phase temperature T g and the wall temperature T w can be calculated, respectively. ρ p stands for the particle density and the specific heat capacities is represented by C p,p for the particle phase and C p,g for the gas phase. V b is the bed volume and ∆H ads the adsorption enthalpy. The heat transfer coefficient h w and h c are related to the gas-to-wall and wall-to-environment heat transfer. More details on this model can be found in the literature (Li et al., 2016a). The physical properties of the gas mixtures are calculated as a linear combination of the properties of the respective mixed gases, which are CO 2 and N 2 . For example, the gas density is calculated as follow: N2 (D.6) in which f i represents the volume fraction of the respective compound i. The other physical properties such as specific heat capacity (C p,g ) and conductivity of gas (k g ) are calculated similarly. The results of the particle and gas phase temperature obtained from the CSTR model have been presented in the aforementioned sections. Correspondingly, the wall temperature T w evolution, which was applied in the CFD-DEM heat loss boundary conditions as a constant value, is shown in Fig. D.7. It is interesting to point out that the temperature of the front sapphire window is influenced by the heat production during the experiment. It gets heated up by the adsorption heat released by the particles and cooled down after the particles became saturated. However, the back PMMA wall has not been significantly affected by the adsorption heat. As mentioned before, the PMMA back wall helps avoiding extra IR radiation emission because of its large heat capacity. Although the front window can be heated up during the experiment, the interference which acted on the recorded IR image is minor due to the low emissivity of sapphire. Hence the assumption of constant window temperature applied in the CSTR and CFD-DEM is reasonable. already understood where I was going with my question or problem and was never short of a new idea. By then I thought his mentoring was like the carrot in front of a donkey. No matter how hard I worked, I always left our meeting with a follow-up question, a remark on something I wrote or a suggestion for further research. Although this frustrated me at times, now I realize it was not carrot but light at the end of the tunnel. I could not have achieved what I have achieved now without him! Aside from the two professors aforementioned, I would also like to extent my gratitude to my daily supervisor Niels Deen. Niels made me feel like a PhD is more than just work. He always asked me about my life and the things I had been up to lately. Because of his approach, I felt like I could ask him any question. In the beginning we always had long discussions about how the project should go and which new fields of research we should explore. He guided me when I needed it the most. Later he helped me with writing my papers while at the same time challenged me to go into even more detail. I was also impressed by how efficiently he could work and it is something I want to incorporate into my own working style. During my last year as a PhD student, Niels got promoted to being a full professor. Despite his full professorship and all the new tasks and responsibilities that came with it, his attention on my project never fainted.
This project would not have been possible without the continuous financial support of the Dutch Polymer Institute (DPI). I would therefore like to thank them for the opportunity to do this project. I always enjoyed the DPI meetings, which I have attended frequently. The project sponsors were often engaged in the success of this project and it made me want to perform even better. I would like to specially express my gratitude to Jan Stamhuis, professor Vincenzo Busico for reviewing the project and organizing the cluster meetings, where I enjoyed a lot with networking and nice DPI evenings. I also would like to thank my industrial partners: Erik Delsman (SABIC), who is also the committee member of my defense ceremony; Mpho Setlhaku (SABIC); Frederik Gemoets (Dow); Jeffrey Brinen (ExxonMobil) and John Walzer (ExxonMobil), for all the discussions we had and insightful input to the project.
In the 3 rd year of my PhD I worked in the lab on chapter 6 of experimental validation. At that time I received lots of help from our professional technicians. They are Joris Garenfeld, Joost Kors and Tijs van Moll. Joris is a very kind helper, he patiently listened to what I needed to modify the set-up and then came up with effective solutions. Joost and Tijs helped me without any hesitation whenever I needed them. Without their efficient and effective assistance in the lab, that part of work would have never been completed in such a good shape.
There were two other contributors who are also specially thanked for their input into this project -Marvin Arts and Tom Janssen. Together we worked on one of the most interesting and also the most difficult part of this project-the experimental validation. There were hard times when nothing seemed to work and no good results could be used to present in your project meetings, and yet you were both still hard working and persistent. I am really happy that in the end you both received a high score on your bachelor/master graduation project. I think it was well-deserved and above all, I really enjoyed the time I spent with either of you! Thank you Tom for all the nice "office-cookies" and all the "best wishes" (Sinterklaas, Christmas, New Year, my birthday, Carnaval, etcetera).
I would like to mention some persons who were not involved in this project but gave constructive suggestions and help on it. They are Frank Peters, Tom Kolkman and Kai Coenen. Frank, many thanks for all the discussions related to the code and especially in the last few months, when I was struggling on finding the right way, your office door was always open when I came to ask for advice. Tom, thank you for all the nice discussions we had about my progress and yours. You never kept any useful information back when you were helping me on solving my problem and these advices from you inspired me to innovate more. Kai, your kind assistance on conducting the TGA experiments and interpreting the results are highly appreciated.
Since I spent lots of time on doing modeling, many of my days have been spend in the office. I started out in the office together with Arash Helmi Siasi Farimani, Ildefonso Campos Velarde and Mariët Slagter. You made my starting days in a foreign country not so lonely and helpless: you organized office outings, dinners and borrels, which were always superfun and full of laughter. Later me and Mariët moved to the Matrix for one year, and Vinayak Sutkar and Mohammad Banaei joined us. Vinayak was left with a private office in the Matrix, while the rest of us moved back to the Helix again. Mariët, me and Mohammad shared an office conveniently close to entrance. From then on we were known as the "DPI team". Later Carin Dietz also joined the office meanwhile Mariët has the other challenge from her new job. Now Mariët and Mohammad are both my paranimphs. I cannot begin to express my gratitude towards all of you. Special thanks to Mohammad, who is not only my officemate but also my project partner, I enjoy a lot to work with you. I could not imagine a better partner to work together, your kindness and knowledge on our field impresses me the most. I am so proud to have shared my working hours with all of my officemates aforementioned. It made those difficult days in which nothing seemed to go right, more enjoyable. Instead of only being colleagues, you have become close friends and this friendship will continue after my PhD.
Aside from my office mates, there are more people that have made my days more colorful. They might started in the group earlier than me, and now they already work somewhere else, however I will always remember the nice time we had together during or out of the office time. They are Yukman Lau, Nhi Dang, Deepak Jain, Paul Hamers, Amit Patil, Yupeng Xu, Dang Chen, NingNing Lu, Lianghui Tan and so on. The others who started together or later than me, like Krushna, Kay, Alvaro, Vincenzo, Luuk, Maaike, Saurish, Shauvik, Marian, Mechelle, Jose, Rommon, Salomon, Alexandro, Giulia, Maxim and so on, thank you for the nice time we had during my days in the group.
We are gathered here in the Netherlands from different cities, my nice fellow country folks, who are still in the group: Yali, Jiangtao, Lei, Lijing, Gaopan. These years we had quite a lot of fun, together we made dumpings, moon cakes and various dishes from different regions of China. Yali, Lei, Lijing, we had a fantastic time together in Canada! These are all my precious memories. Special thanks to Jiangtao, who once got up 6:30am and went to the office just for helping me to print out my draft thesis in order to hand in on time.
I also would like to specially thank two other persons, Judith Wachters-van Gemert and Ada Rijnberg. Both of you always helped me to book the last minute meeting/talk with supervisors. In the last period of completing the draft thesis, thank you Ada for taking care of it when I was in China.
There are many other friends from outside of the office, who also make my life here in the Netherlands much more fun than I could have ever imagined. Tianwei and Lennart, you are "the perfect couple" of ours. We spent lots of nice time together during these years, we went camping, skiing, visiting each other's hometown, and being there for each other's important moments! Huanghuang and Dinesh, girls from Blue, all my co-workers in Hangtang Chinese school, because of you I could take my mind away from the struggling, frustrating times when things went not as good as they supposed to, and also because of you, I could share my joy from the achievement.
I am so glad that my cousin Kai and his family, who are also living here in Eindhoven, can always support me whenever I need them. A great coincidence that we both started our studies here in the Netherlands, and later your wife also became a doctor candidate in TU/e. More joy is among us since your baby girl Mia came to the world and she is the most adorable baby in the world.
No words can express my gratitude to my beloved parents.
Finally, this is for you, Sander, my life partner. Thank you for always being there for me, especially those bad days, you were there for me and help me to deal with my negative emotions. Your family, especially Lena, Alex and Lineke, showed great care to me and your positive attitude towards life not only encouraged me to go through the difficulties but also motivated me to explore more. Without you I could not achieve what I have achieved today. As I mentioned, the 4 years journey of PhD was not an easy trip, with every step I have went through, I was full of joy, because I know, you are also walking with me.