Consortium for Computational Physics and Chemistry

A research collaboration of national laboratories for the U.S. DOE Bioenergy Technologies Office

Liquid phase catalysis (PNNL)

Computational modeling of the upgrading of fast pyrolysis oil

Barriers addressed:

  1. Yields of hydrocarbons
  2. Improve utilization of hydrogen

The CPC has targeted the modeling of two aspects of liquid phase upgrading (LPU) of bio-oil: the engineering of the reactors and the engineering of the reactions. The performance (yield) of a chemical reactor depends on three factors:

  • intrinsic kinetics of the reactions
  • residence time distribution of the reactants
  • micromixing of the contents, for reactions that are not 0 or 1st order,
    Yield = kinetics x RTD x micromixing

Intrinsic reaction kinetics can be measured in reactors that afford rates of mass transport and residence times sufficient to span the relevant ranges of rates and conversions. The residence time distribution can, at least in principle, be measured using a bolus of an inert but detectable tracer. While micromixing can be difficult to characterize experimentally, modern simulations of the fluid dynamics can help define it.

Kinetics. Because measurements of detailed kinetics of the reactions involved in LPU are still pending, for the sake of exercising the differences occasioned by nonidealized residence time distributions, we set the kinetics represent two parallel reactions, a bimolecular deactivation and a first order hydrogenation, both as macroscopic expressions that included Langmuir-type adsorption of the reactant, A (a highly lumped representation of the bio-oil) and the product, B (the precursor of “gunk”, which fouls the catalyst) and the product, C (the desired upgraded product):


The rate constants, ki and equilibrium constants, Ki, (Table 1; note that all concentrations were normalized by the concentration of the feed, A) were chosen to reflect that, at the typical operating temperature of the stabilization reactor of the LPU process, the deactivation reaction is about 10-times faster than the hydrogenation.

Table 1. Parameters employed to model the nominal reaction kinetics.
Parameter Value
k1, k2 1 h-1
k3 0.1 h-1
K 1, K3 1
K2 10
τ 1.0

The resulting concentration profiles (Figure 1) were calculated using the Matlab differential equation solver, ODE45.

Figure 1. Calculated Concentration profiles for the assumed kinetics.

Residence time distributions. Figure 3 compares the normalized Residence Time Distributions for the two scales of reactor with those of three idealized RTDs. The characteristics of the two scaled reactors were chosen to approximate those employed in a bench scale reactor (2.5 cm diameter) with the demonstration scale reactor (10.5 cm diameter) in operation at PNNL. The fluid dynamics were modeled the continuum multi-fluid flow and transport simulator STOMP (White and Oostrom, 2006).

Figure 2. Comparison of Residence Time Distributions from the CFD of two packed bed reactors (top) with those of idealized reactors (botttom). The time scale has been normalized such that the mean residence time for each reactor = 1.0.
Table 2. Parameters employed to model the pilot (bench) and demonstration scale reactors.

In lieu of detailing the geometry of the beds, we employed an empirical correlation devised by Cohen and Metzner (1981) that describes the radial changes in void fraction in the bed (Figure 4), given the diameter of the packing and the diameter of the cylindrical reactor.

Figure 4. Radial variation in the density of random packings in the pilot and demonstration scale reactor.

The bimodal residence time distributions for the packed bed reactors arise from channeling at the wall, which is evident from the flow fields calculated by the fluid dynamics code (Figure 5).

Figure 5. While both scales of reactor exhibit gas bypass along the walls, channeling is severe in the pilot scale reactor, as would be expected from the radial dependence of packing density (vide infra).

Results. The residence time distributions were then used to estimate the conversion of bio-oil to gunk (product B), according to the simplified kinetics scheme described above. The model results show marked differences between the large and small reactors (Table 2). Happily, it appears that the larger reactor, in which less of the volume is channeled, could perform better. Yet, all the reactors appear to be susceptible to the formation of gunk because they all will maintain the polymerizable reaction mixture at temperatures sufficient to promote gunking for times long enough to complete the desired reactions. Evidently, all the results need to be validated to assure that we have employed faithful representations of the reactor geometry and, most importantly, of the reaction kinetics.

Table 3. Results of convolving the Residence Time Distributions of the reactors with the time dependent conversions for macromixing.
Reactor Estimate yield of "gunk"
Bench scale 54%
Idealized plug flow 34%
Demonstration scale 33%
Idealized laminar flow 30%
Idealized continuous stirred tank 25%

Next Steps: We are currently approaching the reaction engineering of the catalysts by calculating the effects of the solvophilicity of the support on the approach of the reactants to the vicinity of the surfaces of supported, small metal particles that catalyze the hydrodeoxygenation reactions. (See

Open source code and tools

Computational models and functions developed by consortium members.

Surface Phase Explorer
Create interactive and downloadable surface phase diagrams from ab initio data.



The CCPC is an enabling project in the ChemCatBio consortium


ChemCatBio is part of DOE’s Energy Materials Network

U.S. DOE Bioenergy Technologies Office

Billion Ton Report
2016 Billion-Ton Report: Advancing Domestic Resources for a Thriving Bioeconomy

NREL Thermochemical Users Facility
Home to thermochemical reactors and pilot plants that CCPC models

PNNL Bioproduct, Sciences, and Engineering Laboratory
Home to upgrading reactors and pilot plants that CCPC models