COMSOL 6.4 - Thermal Decomposition

In this tutorial, the heat and mass transport equations are coupled to laminar flow in order to model exothermic reactions in a parallel plate reactor. It exemplifies how you can use COMSOL Multiphysics to systematically set up and solve increasingly sophisticated models using predefined physics interfaces.

Model Definition

In this model you investigate the unimolecular decomposition of a chemical passing through a parallel plate reactor. A heat–sensitive compound is present in a water solution. After entering the reactor, the liquid first experiences expansion – due to a step in the bottom plate. Before exiting, the fluid also passes a heated cylinder.

The full 3D representation of the reactor geometry is given in Figure 1.

Figure 1: 3D geometry of a parallel plate reactor. The reacting fluid is heated as it passes the cylinder.

The short inlet section of the reactor is considerably wider than it is high. With such a geometry, it is reasonable to assume that the laminar flow develops a parabolic velocity profile between the top and bottom plate. At the same time, the velocity between the side walls is expected to be close to constant (Ref. 1). As a consequence, you can reduce the modeling domain to 2D without dramatically reducing the validity of the simulation (see Figure 2).

Figure 2: Neglecting edge effects, the modeling geometry can be reduced to 2D.

Chemistry

A heat–sensitive chemical (A) undergoes thermal decomposition into fragments (F) according to the following unimolecular reaction in water:

The reaction rate (SI unit: mol/(m3·s)) is given by:

where rate constant k (SI unit: s−1) is temperature–dependent according to the Arrhenius equation:

(1)

In Equation 1, A is the frequency factor (1 × 1010 1/s), E the activation energy (72 × 103 J/mol), Rg the gas constant (8.314 J/(mol·K)), and T the temperature (SI unit: K).

In addition, the decomposition reaction is exothermic, and the rate of energy expelled is given by:

where H is the heat of reaction (−100 kJ/mol).

The reaction kinetics are set up with the Chemistry interface. The fragments’ concentration is for simplicity set as constant in the model, thus only the decomposing species is modeled.

The conversion of species A in the reactor is a function of the residence time; that is, it depends on the detailed fluid flow. Furthermore, the decomposition is influenced by the temperature distribution.

Momentum Transport

The Navier–Stokes equations, which are solved by the single-phase flow interface, consist of the continuity equation:

(2)

and the momentum equations:

(3)

Here, μ denotes the dynamic viscosity (SI unit: Ns/m2), u the velocity (SI unit: m/s), ρ the density of the fluid (SI unit: kg/m3), p the pressure (SI unit: Pa), and F a body force term (SI unit: N/m3). This particular model contains the solution to a steady-state problem, so the first term in each of the equations above disappears.

Apart from the domain equations you also need to select proper boundary conditions. At the inlet you specify a velocity vector normal to the boundary:

(4)

At the outlet boundary you specify a pressure

. Finally, at the surfaces of the reactor plates and the heating cylinder you set the velocity to zero, that is, a no slip boundary condition:

(5)

By selecting the Laminar Flow interface you can easily associate the momentum balance (Equation 2) and boundary conditions (Equation 4, Equation 5, and Equation 5) with your modeling geometry.

Energy Transport

The energy balance equation applied to the reactor domain considers heat transfer through convection and conduction:

(6)

In Equation 6, Cp denotes the specific heat capacity (SI unit: J/(kg·K)), k is the thermal conductivity (SI unit: W/(m·K)), and Q is a sink or source term (SI unit: W/m3).

At the inlet and at the surface of the heating cylinder you set a Temperature boundary condition:

(7)

(8)

At the outlet you set an Outflow boundary condition. This prescribes that all energy passing through this boundary does so by means of convective transport. Equivalently, this means that the heat flux due to conduction across the boundary is zero:

(9)

so that the resulting equation for the total heat flux becomes:

(10)

This is a useful boundary condition, particularly in convection-dominated energy balances where the outlet temperature is unknown.

Finally, assume that no energy is transported across the reactor plates, that is, apply a Thermal Insulation boundary condition:

(11)

Using the Heat Transfer in Fluids interface, you can associate the energy balance (Equation 6) and boundary conditions (Equation 7 to Equation 11) with the modeling geometry.

Mass Transport

The mass transfer in the reactor domain is given by the stationary convection and diffusion equation:

(12)

where Di denotes its diffusion coefficient (SI unit: m2/s), and Ri denotes the reaction term (SI unit: mol/(m3·s)).

Equation 12 assumes that the species i is diluted in a solvent. The mass transport of the fragments are neglected, due to the constant concentration setting of these in the Chemistry interface.

For the boundary conditions, specify the concentration of compound A at the inlet:

(13)

At the outlet, specify that the mass flow through the boundary is dominated by convection. This assumes that any mass flux due to diffusion across this boundary is zero

(14)

and that

(15)

Finally, at the surfaces of the reactor plates and the heating cylinder, assume that no mass is transported across the boundaries — that is, an insulation boundary condition:

(16)

By selecting the Transport in Diluted Species interface you can easily associate the mass balance (Equation 12) and boundary conditions (Equation 13 to Equation 16) with the modeling geometry.

Preparing for Modeling

Before you can start modeling you need to gather the physical data that characterize your reacting flow. For instance, flow modeling requires you to supply the fluid density and viscosity. Mass transport requires knowledge of diffusivities and the reaction kinetics.

Another part of the preparations involves selecting the appropriate physics interfaces and investigating the couplings between different transport equations.

Transport Properties

The term transport properties refers to the physical properties occurring in the transport equations (see the previous section). The momentum and heat transfer equations (Equation 3 and Equation 6) require the following fluid-specific transport properties:

•

Viscosity (μ)

•

Density (ρ)

•

Thermal conductivity (k)

•

Heat capacity (Cp)

The mass transport equation (Equation 12) requires the following species-specific property:

•

Diffusivities (Di)

You need to supply appropriate values of the transport properties to the physics interfaces in order to ensure accurate simulation results. In the present example, water with the dissolved compound A enters the reactor at 300 K. Because water is the solvent, you can assume that its physical properties are representative for the entire fluid. The warmest part of the reactor is held at 325 K. Table 1 lists the transport properties of water as well as the diffusivity of A in water at 300 K and 325 K.

Table 1: Physical properties of liquid water.
Property	at 300 K	At 325 K
Density (kg/m3)	997	987
Viscosity (Ns/m2)	8.5·10-4	5.3·10-4
Thermal conductivity (W/(m·K))	0.62	0.66
Heat capacity (J/(kg·K))	4180	4182
Diffusivity (m2/s)	2.0·10-9	2.0·10-9

When you build this model you make use of the built-in material databases of COMSOL Multiphysics, which automatically provides temperature-dependent properties.

The Flow Regime

The Reynolds number indicates whether a flow is in the laminar or turbulent regime:

As a rule of thumb, a Reynolds number between of 2000 and 2500 marks the transition from stable streamlines to stable turbulent flow. It is always good practice to evaluate the Reynolds number related to the specific flow conditions of the model, because its magnitude guides you to choose the appropriate flow model and corresponding physics interface.

In the present example, you can evaluate the Reynolds number using values from Table 1 and setting the velocity to 5 × 10-4 m/s and the characteristic length to 0.007 m:

Calculating the Reynolds number at 325 K produces a nearly identical result.

The Reynolds numbers are well within the limits of the laminar flow regime.

Dilute or Concentrated Mixtures

When modeling mass transport, it is advisable to discriminate between dilute and concentrated mixtures. For dilute mixtures, Fick’s Law is adequate to describe the diffusional transport. Furthermore, you can assume that the transport properties of the fluid are those of the solvent. For concentrated mixtures, on the other hand, other diffusion models, such as the Maxwell-Stefan model, may be required. The transport properties of the fluid then depend on the mixture composition.

In COMSOL Multiphysics, the Transport of Diluted Species interface is appropriate for dilute mixtures, while the Transport of Concentrated Species interface is recommended for concentrated mixtures.

As a rule of thumb, you can consider concentrations of up to 10 mol% of a solute in a solvent as a dilute mixture.

In the example at hand, the compound A is dissolved in water at a concentration of 1000 mol/m3. As the concentration of pure water is 55,500 mol/m3, the molar fraction of A is approximately 2%. Because the mixture is dilute, it is appropriate to select the Transport of Diluted Species interface for mass transport and to select the transport properties of water as representative values for the mixture.

Solving Coupled Models

As noted previously, the chemistry occurring in the reactor depends both on the fluid flow and the temperature distribution in the reactor. More explicitly, the mass transport equation

(17)

depends on the velocity vector, u, which is solved for in the momentum transfer equation (Equation 3).

Furthermore, the source term Ri in Equation 17 is a function of the temperature, which in turn is the dependent variable of the energy transport equation

When attempting to solve a coupled system of equations such as the one illustrated above, it is often a good idea to analyze the couplings involved and then approach the solution in a stepwise fashion.

In the current model, you first neglect the heat of reaction, Q. This leads to a loose bidirectional coupling between the transport equations:

•

The momentum transport is weakly dependent on the energy and mass transport through the material properties.

•

The energy transport depends only on the momentum transport.

•

The mass transport depends on both the momentum transport and the energy transport.

This structure suggests that it is possible to solve the problem sequentially in the following order: First solve the momentum transport and energy transport problem. Then add the mass transport and investigate the difference.

The last step leads to a fully coupled problem:

•

The momentum transport depends on the energy transport.

•

The energy transport depends on both momentum and mass transport (added heat of reaction).

•

The mass transport depends on both the momentum transport and the energy transport.

In this case you must solve the equations describing all transport phenomena simultaneously.

Results and Discussion

Figure 3 shows the velocity field in the reactor domain along with arrows indicating the velocity magnitude.

Figure 3: Velocity field (m/s) in the reactor.

The cross-sectional area of the fluid increases at the step and decreases at the cylinder, leading to a corresponding local reduction and then increase in the fluid velocity. Recirculation zones appear after the step and the cylinder.

The water solution enters the reactor at a temperature of 300 K and is heated as it passes the cylinder (325 K). Figure 4 shows the temperature distribution in the reactor domain at steady state.

Figure 4: A water solution enters the reactor at 300 K and is heated by a cylinder kept at 325 K.

At the reactor inlet, the concentration of A is 1000 mol/m3. Figure 5 shows the concentration of A as the compound undergoes decomposition.

Figure 5: Concentration of the heat sensitive chemical (A) (mol/m3) as function of position in the reactor.

These plots make it possible to identify some general trends. It is clear that decomposition occurs mainly after the liquid has been heated by the cylinder. In the first half of the reactor, where the temperature is relatively low, decomposition is still fairly advanced near the wall and after the step. This is due to the longer residence times in these areas. In the second part of the reactor, where heating takes place, regions with relatively high concentrations of compound A are visible. This also makes physical sense because the water velocity is relatively high.

The temperature distribution in the entire reactor is affected by the heat of reaction. As shown in Figure 6, the maximum fluid temperature now exceeds the temperature of the heating cylinder. Furthermore, the water temperature is higher than 300 K in the region between the inlet and the cylinder.

Figure 6: Reactor temperature (K) when the heat of reaction is taken into account.

Figure 7 plots the rate of reaction as a function of the position in the reactor. Clearly, significant reaction now occurs in the first part of the reactor, before the heating cylinder.

Figure 7: Decomposition rate of compound A. Significant decomposition occurs in the first half of the reactor.

Reference

1. H. Schlichting, Boundary Layer Theory, 4th ed., McGraw Hill, p. 168, 1960.

Chemical_Reaction_Engineering_Module/Reactors_with_Mass_and_Heat_Transfer/thermal_decomposition

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Import some model parameters from a text file.

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file thermal_decomposition_parameters.txt.

Geometry 1

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type W1.

In the text field, type H1.

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type W2.

In the text field, type H2.

Circle 1 (c1)

In the toolbar, click

In the window for , locate the section.

In the text field, type R1.

Locate the section. In the text field, type xpos.

In the text field, type ypos.

Add a line segment extending from the back of the cylinder. This will be used to refine the mesh in the wake region.

Line Segment 1 (ls1)

In the toolbar, click

and choose .

On the object , select Point 3 only.

In the window for , locate the section.

From the list, choose .

In the text field, type xpos+10*R1.

In the text field, type ypos.

Click

Difference 1 (dif1)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

Click to select the

toggle button for .

Select the objects and only.

Click

Mesh Control Edges 1 (mce1)

In the toolbar, click

and choose .

On the object , select Boundary 6 only.

In the toolbar, click

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Materials

Water (mat1)

By default, the first material you add applies on all domains.

Definitions

In preparation for defining boundary conditions, it is practical to define some named selections.

Inlet

In the toolbar, click

In the window for , type Inlet in the text field.

Locate the section. From the list, choose .

Select Boundary 1 only.

Outlet

In the toolbar, click

In the window for , type Outlet in the text field.

Locate the section. From the list, choose .

Select Boundary 6 only.

Heater

In the toolbar, click

In the window for , type Heater in the text field.

Locate the section. From the list, choose .

Select Boundaries 7–9 only.

Laminar Flow (spf)

Follow the instructions below to set up the interface. The fluid properties are automatically taken from the material assigned to the reactor domain, so all you need to do is to define inlet and outlet boundary conditions.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Inlet 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Locate the section. In the Uav text field, type v_inlet.

Outlet 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. Select the checkbox.

For a more accurate simulation, edit the default mesh by refining it in the region around and behind the heater.

Mesh 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Size

In the window, under > click .

In the window for , locate the section.

Click the button.

Locate the section. In the text field, type 1.1.

Size 1

The node provides mesh refinement on wall boundaries. Include the internal boundary behind the cylinder as well. When done all boundaries except the inlet and outlet should be selected.

In the window, click .

Select Boundaries 2–5 and 7–11 only.

Boundary Layers 1

Due to the heating, the reaction will be promoted in vicinity of the cylinder. Improve the resolution by increasing the number of boundary layer elements on the cylinder surface.

Boundary Layer Properties 2

In the window, expand the > > node.