|Multiphysics Model Setup||Solve the PDE Model|
|Physics Domains Coupling||Nomenclature|
Thermal decomposition is a chemical reaction where heat is a reactant. Since heat is a reactant, these reactions are endothermic meaning that the reaction requires thermal energy to break the chemical bonds in the molecule. A classic example is the calcination process:
Calcium carbonate, , will decompose into carbon dioxide, , and calcium oxide, , when heated above at a pressure of atmosphere. For endothermic reactions in this case the heat of reaction, , is negative. This reaction is known as the common way to make quicklime, , that is an industrially important product.
This study is to model the decomposition process of calcium carbonate, and simulate its concentration field within an expanding reactor. The calcium carbonate is transported into the reactor by suspending it in water that enters the domain from the bottom, flows across a heated tube, and exits the reactor through the top. During this process the calcium carbonate particles decompose and absorb heat from the system.
The symbols and corresponding units used for this model are summarized in the Nomenclature section.
Since this problem considers more than one physics domain, a multiphysics model is to be constructed. As a first step we look at the contributing physics domains independently. A fluid dynamics model, a heat transfer model and a mass transport model are set up to solve for the flow field , temperature field and concentration field within the reactor.
The Navier-Stokes equation (1) is used to solve for the steady state flow field within the reactor. Note that we use a 2D domain to simulate the 3D reactor, the reason for this model order reduction will be explain in the later section about the modeling domain setup.
Since the thermal decomposition is an endothermic process, the heat absorbed by the reaction is modeled by the heat sink term , which depends on both the chemical reaction rate, , and the heat of reaction, :
The steady state mass transport equation (3) is used to solve for the concentration field within the reactor. The mass transported by the flow field is modeled by the convective term .
As described in the previous section, the temperature field within the reactor is computed with the heat equation (6) and depends on both the fluid flow velocity computed with the Navier-Stokes equation (7) and the concentration through reaction rate of the mass transport equation (8). Furthermore, the reaction rate itself also depends on the temperature from the heat equation (9), resulting in nonlinear, two-way relation between each physical mode.
To solve a two-way coupled nonlinear PDE system like this, a possible strategy is to first solve a simpler model with a weaker nonlinearity and use the solution of the simpler model as an initial guess for the full nonlinear model. Ignoring the heat generated by the reaction gives a good simpler model to start with. This approach is detailed in the following sections.
By ignoring the heat of reaction (i.e. , ) from the decomposition, the heat equation (10) is simplified to: . This leads to an one-way coupled multiphysics model. For clarity we have highlighted the dependent variable for each physical mode:
This one-way coupled multiphysics model can either be solved sequentially or in one go. Solving it sequentially would mean that each of the physics domains is solved with a separate call to NDSolve. This approach of splitting the coupled equations in chunks and solve them independently and sequentially is called a segregated solution process. In case the single physics fields can not be split easily an initial guess for the coupling variable of the first PDE is made. Then the solution of the first PDE is take to drive the next PDE and so forth. Once all PDEs have been solved, the process is started over again. This will improve the initial guess for the initial coupling variable. This procedure is executed until the solution converges.
Note that the temperature field and the concentration field would affect with each other, and should therefore be solved as a system of PDEs in NDSolve. This solving process will be shown in the later section: "Solve the PDE Model: Fully Coupled Model Solution".
Since the depth ( direction) of the reactor is considerably more than its width ( direction) and its height ( direction) it is reasonable to neglect the variation of in the direction, so a two dimensional model is sufficient to represent the 3D reactor.
As a further step to simplify the geometry note the symmetry about the vertical () axis. It is effective to only make use of the right half of the reactor as the simulation domain . Here and denote the flow inlet and flow outlet. The wall boundary and the symmetric axis are symbolized as and , respectively.
In order to get a good result, a finer than the default grid is used for the mesh generation. Here the maximum grid size is set to , which means that there will be about a hundred elements in the direction (height).
When a simulation involves fluid flow, it is advisable to calculate the Reynolds number to determine the type of the flow. Empirically, the flow is considered as laminar when and as turbulent when . For a Reynolds number between roughly an unsteady separation of flow may form around blunt bodies , resulting in an oscillating vortex shedding in the flow field. Whenever this happens, the steady state fluid dynamics model will failed to converge and a steady state solution is not achievable.
The parameters for the reactor are chosen so that the Reynolds number falls into the laminar flow regime, which means the chemical reactor in this case is a laminar flow reactor (LFR). One feature of the LFR is the higher residence time  (the time interval during which the chemicals stay in the reactor), and this is desirable for the calcination to achieve higher production rate.
The wall boundaries are assumed to be perfectly thermally insulated and are modeled by a thermally insulated boundary condition. An outflow boundary condition and a symmetric boundary condition are applied on the flow outlet and the symmetric axis , respectively.
However, since the outflow boundary condition, symmetric boundary condition and thermally insulated boundary condition are all Neumann zero conditions, they are implicitly applied without further setup.
On the walls of the reactor the mass flux of the species is assumed to be zero, and is modeled by an impermeable boundary condition. An outflow boundary condition and a symmetric boundary condition are applied on the flow outlet and the symmetric axis , respectively.
A stable solution can be found if the velocities are interpolated with a higher order than the pressure. NDSolve allows an interpolation order for each dependent variable to be specified. Here the velocities and are set to be interpolated with second order and the pressure with first order.
The flow enters the reactor at a temperature of and is heated as it passes by the heating tube . Note that since we ignore reaction heat generated in this case, the temperature variance is caused by the tube heating only.
As described in the previous section, particles within the reactor are decomposing at a reaction rate .
The concentration of is held at at the inlet, and then deceases as it flows over the reactor due to decomposition. Note that the concentration drops significantly in the lower part of the reactor because of the higher reaction rate. We can verify this by visualizing the reaction rate field: .
Note that a significant decomposition occurs in the lower part of the reactor around the flow inlet. This explains why the previous plot showed a significant drop of concentration in this region.
Making use of the one-way coupling model, we've solved for the flow velocity , temperature and the concentration in a sequential manner. Now we are ready to solve for a fully coupled model and investigate the difference.
In the previous section we've shown the heat source term depends the reaction rate by:
To efficiently solve the coupled PDE model, the result from the previous step is used as an initial seed for the new solution. In fact not doing so will result in a model where the nonlinear solver does not converge. It is a general strategy to first solve a less nonlinear model and use the solution there of as an initial seed for more nonlinear model.
To measure the effectivity of the reactor, the net of reduction is calculated between the flow inlet and outlet . For this the mean concentration at the boundary is computed with a boundary integration:
1. Tansley, Claire E. and Marshall, David P. Flow past a Cylinder on a Plane, with Application to Gulf Stream Separation and the Antarctic Circumpolar Current, Journal of Physical Oceanography. 31 (11): 3274–3283. (2001).