03-mod03.Rmd

# Distribution of Residence times

Reactors treated in this course have been modeled under specific flow considerations that lead to simplicity in the mole and energy balances and are valid in many industrial reactors.

Unfortunately, as it is shown in Figure \@ref(fig:fig004000), there are other cases where the deviation of those industrial reactors from ideal behavior may be sufficiently large that the use of an ideal reactor model will give an unacceptably significant error because of the careless treatment that these models do to the fluid flow pattern within the reactors. 

```{r fig004000, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Examples of nonideal flow in a CSTR and packed-bed vessel. [@missen1999]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_25.png")
grid.raster(img)
```

If the flow pattern does not follow the ideal models, it must be necessary to apply a more complex methodology for the reactor performance analysis. Nevertheless, before going deeper into these methodologies, it is essential to emphasize that  the factors that govern the performance of a chemical reactor can be grouped into two broad classifications:

 1. The **amount of time** that a reactant molecule spends in the reactor. The longer that reactant can spend in the reactor, the more chance they can react. The time that molecules spend in the reactor is called **the residence time**, which is a crucial parameter in the flow reactor performance, and the **distribution of residence time** is a key factor in determining the extent of reaction.
 2. The **extent of mixing** within the reactor. When dealing with two or more different molecules that must react, molecular collision must occur as a first step. Therefore, the concentration at the molecular level is essential.

The degree of mixing depends not only on the reactor configuration and operating conditions but also on the fluid and reactants' properties.

Mixing is a complex phenomenon that can be analyzed on both a macroscopic and microscopic scale. The RTD functions provide a measure only of macroscopic mixing. [@missen1999] 

However, non-ideal flow results in irregularities relating to the micromixing of fluid elements of differing ages at the microscopic level, too. Therefore, in this section of the course, we will see how the classical approach of Residence Time Distribution (RTD) provides a method for evaluating the effect of flow patterns within the reactor on chemical conversion levels.

## RTD functions

First of all, the study of non-ideal flow behavior is based on the RTD of non-reacting species in an arbitrary flow system with constant density fluid. From this point, it will be explained its relationship to a chemically reacting species present in the reactor.

Residence time theory deals with the age of particles within a flow system. The particles may be atoms, molecules, or fluid elements conserved as they flow through the system.

:::: {.infobox .info data-latex="information"}
  
The residence time, or age, of a particle in the system is the elapsed time between the time at which the particle enters the system and the time it leaves.
::::

This approach assumes specific essential considerations:

1. If a particle reenters the system after it leaves, the aging of that particle resumes from the value it had as its previous exit.
 2. A particle is not allowed to remain in the system forever.
 3. All particles must have an original entrance and a final departure.
 4. The age of the particle when it finally leaves the system will be its residence time.

Regardless the state of mixing or type of flow, we use the classical definition of mean residence time $\overline{t}$ of an element of fluid volume within a fully contained, continuous flow reactor with a total fixed volumen of $V$ operated at steady state which is:
\begin{equation}
  \overline{t} = \frac{V}{\upsilon}
(\#eq:ec004001)
\end{equation}
where $\upsilon$ is the steady volumetric flow rate.

It is vital to point out that a fully contained vessel has no diffusion or dispersion at the vessel's entrance or exit so that a fluid element once inside the reactor cannot escape back into the input stream or once removed from cannot return to the vessel. 

This assumption is more precisely related to the boundaries of the reactor. It is assumed a plug flow in both the inlet tube and the outlet tube. Such boundary conditions are called "closed-closed." Figure \@ref(fig:fig00400) gives others boundary conditions than closed-closed that they are essential for specific experimental techniques to assess the mixing parameters inside the reactor, but they will not consider in this course. 

```{r fig00400, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Boundary conditions of reactors in relation to the residence time distribution. [@westerterp1984]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_28.png")
grid.raster(img)
```

Now, let us introduce some mathematical concepts. Let $d\,V = \upsilon\, dt$ represent a fluid element which entered the system $t$ time ago. Also, let $\mathbf{P}(t)$ represent the probability that the element remains in the reactor during the period $t$. The total volume of the fluid in the reactor must be

\begin{equation}
  V = \upsilon \int_{0}^\infty {\mathbf{P}(t)\,dt}
  (\#eq:ec004002)
\end{equation}

Now, let us introduce some mathematical concepts. Let $d\,V = \upsilon\, dt$ represent a fluid element which entered the system $t$ time ago. Also, let $\mathbf{P}(t)$ represent the probability that the element remains in the reactor during the period $t$. The total volume of the fluid in the reactor must be

Now, let $\mathbf{F}(t)$ represent the probability of a fluid element leaving the system during the period $t$ so $\mathbf{P}(t) + \mathbf{F}(t) = 1$. The foregoing equation then becomes

\begin{equation}
  \frac{V}{\upsilon} = \int_{0}^\infty {\left(1 - \mathbf{F}(t)\right)\,dt}
(\#eq:ec004003)
\end{equation}
which upon integration by parts gives

\begin{equation}
  \frac{V}{\upsilon} = \left[1 - \mathbf{F}(t)\right] \cdot \infty - \left[1 - \mathbf{F}(t)\right] \cdot 0+\int_{F(0)}^{F(\infty)} {t\,d\mathbf{F}(t)}
(\#eq:ec004004)
\end{equation}
with boundary conditions of $F(t) = 0$ at $t=0$, and 1 at $t=\infty$ (chance of 100% removal in an infinitive time) so that
\begin{equation}
  \frac{V}{\upsilon} = \int_{0}^{1} {t\,d\mathbf{F}(t)} = \overline{t}
(\#eq:ec004005)
\end{equation}

\begin{equation*}
  \mathbf{F}(t) \to 1 \text{ when } t \to \infty
 \end{equation*}
  
$\mathbf{F}(t)$ is called "_Cumulative Distribution Function_''. This function is dimensionless, and **represents the fraction of effluent that has been in the reactor for less time than _t_, and leaves the reactor with age less than _t_**.

The typical shape of $\mathbf{F}(t)$ function is shown on Figure \@ref(fig:fig004000a). From this figure, it is possible to note that 80% of the molecules spend 40 min or less in the reactor, and 20% of the molecules spend longer than 40 min in the reactor.

```{r fig004000a, out.width = "55%", echo=FALSE, fig.align='center', fig.cap='Typical cumulative distribution curve **F**. [@fogler2008]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_16.png")
grid.raster(img)
```

The distribution of residence times is represented by an external residence time distribution or an exit age distribution, $\mathbf{E}(t)$. This function is called **Residence Time Distribution Function** or (RTD-function) which **defines a fraction of exiting fluid elements with residence times lies between time $t$ and $t+dt$.**

Since $\mathbf{E}(t)$ is a fraction, $\mathbf{E}(t)$ has units of time$^{-1}$. In addition, the fraction of fluid that leaves the reactor over all time, this means between $t=0$ and $t=\infty$ must be 1. Therefore
\begin{equation}
 \int_{0}^{\infty} \mathbf{E}(t)\,dt = 1
(\#eq:ec004008)
\end{equation} 

The fraction of fluid **_in the effluent stream_** that was **_in the reactor_** for a time between $t=0$ and $t=t$  si given by
\begin{equation*}
  \int_{0}^{t} \mathbf{E}(t)\,dt 
\end{equation*} 

Another way of saying the same thing is that $\int_{0}^{t} \mathbf{E}(t)\,dt$ is the fraction of fluid in the exit stream with an \textit{exit age} less than $t$. This fractions can be represented graphically as shown in Figure \@ref(fig:fig004019).

```{r fig004019, out.width = "70%", echo=FALSE, fig.align='center', fig.cap='Typical RTD curve **E**. [@roberts2009]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_17.png")
grid.raster(img)
```

The fraction of the fluid that spends a given duration, $t$ inside the reactor, is given by the value of $\mathbf{E}(t)\,dt$. Therefore, $\mathbf{E}(t)$ curve is related to $\mathbf{F}(t)$ because the fraction of the fluid that leaves the reactor with age less than $t$ is equal to the sum over all times less than $t$ of $\mathbf{E}(t) \Delta t$, or expressed continuously. [@roberts2009]

\begin{equation}
  \mathbf{F}(t) = \int_0^{t} \mathbf{E}(t)\,dt
(\#eq:ec004009)
\end{equation}

Equation \@ref(eq:ec004009) could be rewritten as
 \begin{equation}
  \frac{d \mathbf{F}(t)}{dt} = \mathbf{E}(t)
(\#eq:ec004006)
\end{equation} 

Therefore from Eq. \@ref(eq:ec004005)
 \begin{equation}
 \overline{t} = \int_{0}^{\infty} t\,\mathbf{E}(t)\,dt
(\#eq:ec004007)
\end{equation} 

This procedure is called normalizing the distribution, and Figure \@ref(fig:fig004001) shows that.

```{r fig004001, out.width = "70%", echo=FALSE, fig.align='center', fig.cap='The typical exit age distribution curve **E** (RTD) [@levenspiel2004]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_01.png")
grid.raster(img)
```

The RTD function may be used directly to analyze reactors, or related functions may be used depending on the application.

Another useful distribution curve is the Washout function, $\mathbf{W}$, which sometimes is called "Complement function''. This function gives **the fraction of effluent that has been in reactor (with the residence times) for longer than time $t$**. It could be calculated from
\begin{equation}
  \mathbf{W}(t) = \int_{t}^{\infty} \mathbf{E}(t)\,dt \equiv 1 - \mathbf{F}(t) = 1 - \int_0^t{\mathbf{E}(t)\,dt}
(\#eq:ec004010)
\end{equation}
So,
\begin{equation}
  \mathbf{F}(t) + \mathbf{W}(t) = 1
(\#eq:ec004011)
\end{equation}
Furthermore, the relationship between the functions in derivative form is
\begin{equation}
  \frac{d\mathbf{F}(t)}{dt} =- \frac{d\mathbf{W}(t)}{dt} = \mathbf{E}(t)
(\#eq:ec004012)
\end{equation}

A typical $\mathbf{F}(t)$ and the corresponding $\mathbf{W}(t)$ function curve are illustrated in Figure \@ref(fig:fig004002)

```{r fig004002, out.width = "70%", echo=FALSE, fig.align='center', fig.cap='The typical **F(t)** and **W(t)** curves. [@hayes2013]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_02.png")
grid.raster(img)
```

Another distribution function of interest in some application is the "Internal Age Distribution'', $\mathbf{I}$, which is **the fraction of fluid inside the vessel with age between $t$ and $t + dt$**. It has properties similar to $\mathbf{E}(t)$, so
\begin{equation*}
    \int_0^\infty \mathbf{I}(t)\,dt = 1
\end{equation*}
The time $t=0$ refers to an arbitrary initial time and not the start of the flow of fluid into the vessel. In physical terms, this equation states that all fluid has an age between 0 and $\infty$.

As a consequence, the fraction of vessel contents younger that a specified age $t$ is,
\begin{equation*}
     \int_0^{t} \mathbf{I}(t^\prime)\,dt^\prime = 1
\end{equation*}
While the fraction older than $t$ is
\begin{equation*}
     \int_t^{\infty} \mathbf{I}(t^\prime)\,dt^\prime = 1 -  \int_0^{t} \mathbf{I}(t^\prime)\,dt^\prime
\end{equation*}
Similar to $\mathbf{E}(t)$, the mean age fluid elements inside the vessel is,
\begin{equation*}
 \overline{t}_1 =   \int_0^\infty t\,\mathbf{I}(t)\,dt = \text{ mean internal age }
\end{equation*}

As might be expected, there is an interrelationship between $\mathbf{I}$ and $\mathbf{E}$ since the fluid entering a vessel at a given time obviously either leaves it or stays inside. Therefore:
\begin{equation}
  \overline{t}\,\mathbf{I}(t) = 1 - \int_0^t \mathbf{E}(t)\,dt = \mathbf{W}(t)
(\#eq:ec004012a)
\end{equation}
or
\begin{equation}
  \mathbf{I}(t) = \frac{\mathbf{E}(t)}{1 - \mathbf{F}(t)} \Rightarrow \mathbf{E(t)} = - \overline{t}\;\frac{d\mathbf{I}(t)}{dt}
\end{equation}

Finally, the **Intensity function**, $\boldsymbol\Lambda(t)$ is defined as **fraction of fluid in the vessel of age $t$ will leave at time between $t$ and $t + dt$**. The intensity is useful in detecting the existence of dead space and bypassing.

The intensity function can be related to  $\mathbf{E}$ and $\mathbf{I}$ functions from

<center>

![](figures/intensity_E_I.png)

</center>

Therefore
\begin{equation}
  \boldsymbol\Lambda (t) = \frac{1}{\overline{t}} \frac{\mathbf{E}(t)}{\mathbf{I}(t)} 
\end{equation}

Combining the last two equations
\begin{equation}
  \boldsymbol\Lambda (t) = \cfrac{- \overline{t}\;\cfrac{d\mathbf{I}(t)}{dt}
}{\overline{t}\,\mathbf{I}(t)} = - \frac{d\left[\ln\left(\overline{t}\mathbf{I}(t)\right)\right]}{dt}
\end{equation}

## Means and moments

In addition to characterizing the RTD using the three functions $\mathbf{E}(t)$, $\mathbf{F}(t)$ and $\mathbf{W}(t)$, the residence time can also be characterized using the moments of the RTD functions. There are several types of moments. [@hayes2013]

The first type is the moment around the origin, which is defined in terms of density function as
\begin{equation}
  \mu_n = \int_0^{\infty} t^n\,\mathbf{E}(t)\,dt
(\#eq:ec004013)
\end{equation}

The value of **n** determines which moment is calculated. For example, the first moment about the origin is the **mean residence time**, $\overline{t}$
\begin{equation}
  \mu_1 = \overline{t} = \int_0^{\infty} t\,\mathbf{E}(t)\,dt
(\#eq:ec004014)
\end{equation}
As we already know, the mean residence time is the average time that molecules spend in the system. In a constant-density system $\overline{t} = V/\upsilon = \tau$

:::: {.infobox .info data-latex="information"}
Equation \@ref(eq:ec004013) can provide a helpful check on operating data because the volumetric flow rate is relatively easy to set and measure. However, the volume of the reactor filled by a fluid is not always so easy to determine.
::::

Rather than compute the moments about the origin, it is common for the moments of the RTD function to be calculated about the mean residence time. These moments are called central moments, and are defined as
\begin{equation}
  \mu_n^{\prime} = \int_0^{\infty} {\left(t-\overline{t}\right)}^{n}\,\mathbf{E}(t)\,dt
(\#eq:ec004015)
\end{equation}

The variance, $\sigma_t^2$, of the RTD is the second central moment:
\begin{equation}
  \sigma^2 = \mu_2^{\prime} = \int_0^{\infty} {\left(t-\overline{t}\right)}^2\,\mathbf{E}(t)\,dt
(\#eq:ec004016)
\end{equation}

An alternative form may be more convenient to use is obtained by expanding the square in Eq. \@ref(eq:ec004016)
\begin{equation}
  \sigma_t^2 = \left(\int_0^{\infty} {t}^2\,\mathbf{E}(t)\,dt\right) - \overline{t}^{\,2}
(\#eq:ec004016a)
\end{equation}

The variance measures the spread of the distribution about the mean characterized by  the standard deviation, $\sigma_t$. The third central moment, the skewness, measures the symmetry of the distribution about the mean:
\begin{equation}
  \sigma^3 = \mu_3^{\prime} = \int_0^{\infty} {\left(t-\overline{t}\right)}^3\,\mathbf{E}(t)\,dt
(\#eq:ec004017)
\end{equation}

All moments must be calculated meticulously for a complete summary of a distribution.  In practice, these three parameters are typically appropriate for a good RTD characterization.

## Normalized RTD functions
It is common when using residence time theory to work with normalized RTD, defined as the residence time divided by the mean residence time. This dimensionless time is given by [@hayes2013]
\begin{equation*}
  \theta = \frac{t}{\overline{t}}
\end{equation*}
The density functions for the two cases are related by
\begin{equation*}
 \mathbf{E}(\theta)\,d\theta = \mathbf{E}(t)\,dt 
\end{equation*}
Because $\overline{t}\,d\theta = dt$, it follows that
\begin{equation*}
 \mathbf{E}(\theta) = \overline{t}\,\mathbf{E}(t) = \frac{\overline{t}\;C(t)}{\displaystyle \int_0^\infty C(t)\,dt}
\end{equation*}
The other relationships do not change. The normalized cumulative distribution function is
\begin{equation*}
  \mathbf{F}(t) = \int_0^{t} \mathbf{E}(t)\,dt = \int_0^{\theta} \mathbf{E}(\theta)\,d\theta = \mathbf{F}(\theta)
\end{equation*}

The normalized washout function is
\begin{equation*}
  \mathbf{W}(\theta) = \int_{\theta}^{\infty} \mathbf{E}(\theta)\,d\theta
\end{equation*}
Finally, the moments about the origin and the central moments for the normalized RTD functions are
\begin{align*}
  &\nu_{(1,\theta)} = \int_{\theta}^{\infty} \theta \,\mathbf{E}(\theta)\,d\theta = \overline{\theta} = 1 \\
  &\nu_{\theta}^{2} = \int_{\theta}^{\infty} {\left(\theta - \overline{\theta}\right)}^2\,\mathbf{E}(\theta)\,d\theta = \int_{\theta}^{\infty} {\theta}^2\,\mathbf{E}(\theta)\,d\theta - 1 = \frac{\theta^2}{\overline{t}^2}
\end{align*}

For numerical work with tracer data directly, it is common to use the dimensional $\mathbf{E}(t), \overline{t}$, and $\sigma_t^2$. The dimensionless quantities are used with the performance models. Figure \@ref(fig:tabrtd) shows a summary of RTD functions expressed by dimensionless time.


```{r tabrtd, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Relationships between Residence Time Functions', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/tab_rtd.png")
grid.raster(img)
```

## Experimental determination of RTD

The RTD is measured experimentally by injecting a tracer, which is an inert chemical compound into the reactor at some time $t=0$ and then measuring its concentration, $C$, in the effluent stream as a function of time. The tracer must meet the following requirements:

1. it must be a nonreactive species that is easily detectable at low concentrations, e.g., by electrical conductivity, light absorption, or nuclear radiation.
2. it should have physical properties similar to those of the reacting mixture, such as density, and be completely soluble in the mixture to avoid flow changes due to viscosity and molecular diffusion of the tracer in the solution.
3. it should not adsorb on the walls or other surfaces in the reactor.
 4. The tracer should be inexpensive and easy to handle.
 5. The tracer should be stable and conserved so that it can be accounted for by a material balance relating the response to the input; if tracer decays, its half-life should be such that $t_{1/2}\,(\text{tracer}) > 25\overline{t}\,  (\text{fluid})$.

The two most used methods of injection are _pulse input_ and _step input_.

## Determination of **E(t)** from Pulse input

The direct experimental determination of **E(t)** requires a sharp pulse of tracer to be injected into the system. Essentially, an amount of tracer $N_0$ is required to be injected into the vessel's inlet at a single instant of time at $t=0$. Figure \@ref(fig:fig004003) shows typical concentration-time curves at the inlets and outlet of an arbitrary reactor. [@hayes2013]

```{r fig004003, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Typical pulse impulse curves for RTD measurement', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_03.png")
grid.raster(img)
```

Mathematically, the ideal shape of the injected pulse is described by a Dirac delta function. In practice, the tracer should be injected in as short an elapsed time as is experimentally possible to achieve.

The tracer's concentration in the effluent is then measured as a function time, and the density function can be evaluated using the following analysis. The fraction of tracer molecules that leave the system over some time interval $dt$ is

\begin{equation}
  \frac{dN}{N_0} = \mathbf{E}(t)\,dt
(\#eq:ec004018)
\end{equation}
where $dN$ represents the number of tracer molecules that leaves the system over the time interval $dt$, and $N_0$ is the total number of tracer molecules injected.

The number of molecules can be expressed in terms of the concentration and the volumetric flow rate,
\begin{equation}
  dN = \upsilon\,c(t)\,dt
(\#eq:ec004019)
\end{equation}
where $C(t)$ is the concentration of tracer in the effluent streams as a function of time.

Thus, the total number of molecules initially injected can also be calculated from the effluent profile by integration over a long time:

\begin{equation}
  N_0 = \int_0^{\infty} \upsilon\,C(t)\,dt \equiv \frac{N_0}{\upsilon} = \int_0^{\infty} C(t)\,dt
(\#eq:ec004020)
\end{equation}

Combining Eqs \@ref(eq:ec004019) and \@ref(eq:ec004020) into Eq. \@ref(eq:ec004018) for a constant-density system,
\begin{equation}
 \mathbf{E}(t) = \frac{c(t)}{\displaystyle\int_0^\infty{c(t)\,dt}} \approx \frac{c_i(t)}{\displaystyle\sum_i {c_i(t)\,\Delta t_i}}
(\#eq:ec004021)
\end{equation}

Once $\mathbf{E}(t)$ is determined, the mean residence time, $\overline{t}$, and the variance of distribution, $\sigma_t^2$, may be calculated based on Eqs. \@ref(eq:ec004007)  and \@ref(eq:ec004016a)
\begin{gather*}
 \overline{t} = \int_{0}^{\infty} t\,\mathbf{E}(t)\,dt \approx \sum_i t_i E_i(t) \Delta t_i \\
   \sigma_t^2 = \left(\int_0^{\infty} {t}^2\,\mathbf{E}(t)\,dt\right) - \overline{t}^{\,2} \approx \left[ \sum_i t_i^2 E_i(t) \Delta t_i \right] - \overline{t}^2
\end{gather*} 

Advantages of using a pulse input include:

1. requiring only a small amount of tracer.
2. involving usually only a tiny impact on process operation.

Disadvantages include:

1. difficulty in achieving a perfect pulse.
2. difficulty in achieving accurate material balance on traces using equation \@ref(eq:ec004020)

Obtaining $\mathbf{E}(t)$, $\overline{t}$ and $\sigma_t^2$ from experimental tracer data involves determining areas under curves defined continuously or by discrete data. Two simple ways in which discrete data may be treated to obtain the required areas:

1. use of data in histogram form, and
2. use of the trapezoid rule.

In the histogram method, Figure \@ref(fig:fig004003a), the area under the response curve is
\begin{gather*}
  \text{area} = \sum_{i=1}^{n_C} c_i(t)\,\Delta t_i \\
  \Delta t_i = \frac{t_{i+1} - t_{i-1}}{2}
\end{gather*}

```{r fig004003a, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Basis for histogram method of area determination. [@missen1999]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_26.png")
grid.raster(img)
```

In the use of the trapezoid rule, Figure \@ref(fig:fig004003b), the area under the response curve is given by

```{r fig004003b, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Basis for histogram method of area determination. [@missen1999]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_27.png")
grid.raster(img)
```

\begin{equation*}
  \text{area} = \sum_{i=1}^{n_c - 1} c_i(t)_\text{av}\,\Delta t_i = \sum_{i=1}^{n_c - 1} \left(\frac{c_i + c_{i+1}}{2}\right) \left(t_{i+1} - t_i\right)
\end{equation*}

The only way to determine which numerical method provides the best answer is to predict the experimental concentrations using the calculated parameters in an appropriate mixing model.

```{exercise, it004001, name="Calculations of RTD and distribution functions curves from pulse impulse method"}
Consider a vessel through which a fluid is flowing. The RTD function is determined utilizing a pulse injection of the tracer in the inlet stream, and its outlet concentration is monitored as a function of time. The output concentration in mol/L measured as a function of time in seconds is

 t (s) | C (mol/L) | t (s) | C (mol/L) | t (s) | C (mol/L) | t (s) | C (mol/L) 
--|--|--|--|--|--|--|--
0.0 | 0 | 225 | 7.4 | 275 | 8.2 | 375 | 0.5
150 | 0 | 240 | 9.4 | 300 | 5.0 | 400 | 0.2
175 | 1 | 250 | 9.7 | 325 | 2.5 | 450 | 0.0
200 | 3 | 260 | 9.4 | 350 | 1.2 | 500 | 0.0

Determine:
  
1. Plot of tracer concentration vs time.
2. RTD curve, $\mathbf{E}(t)$.
3. Cumulative distribution function plot, $\mathbf{F}(t)$.
4. Washout function plot, $\mathbf{W}(t)$.
5. Mean residence time.
6. Fraction of material leaving the reactor that has spent between 230 and 270 seconds.

```

_Solution_

The solution in a Python script is available at  [link](https://github.com/aliglara/ebook_che5314_ucv/blob/main/codes/lect_021RTDimpulse.ipynb).

(a) By plotting _C_ as a function of time, the curve shown in Figure \@ref(fig:fig004005a) is obtained. From the $C(t)$ curve is possible to calculate the RTD, according to Eq. \@ref(eq:ec004021), for a pulse impulse test.

We just need to divide each point of c(t) by the integral
 \begin{equation*}
  \int_0^\infty c(t)\:dt
\end{equation*}
which is just the area under the $C (t)$ curve.

```{r fig004005a, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Trace concentration curve', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_01_01.png")
grid.raster(img)
```

Figure \@ref(fig:fig004005b) shows the obtained RTD curve

```{r fig004005b, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='RTD curve', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_01_02.png")
grid.raster(img)
```

(b) The cumulative distribution curve, _F(t)_, and washout curve, _W(t)_, were obtained from _E(t)_ curve, and Figure \@ref(fig:fig004006a) shows the results obtained.

```{r fig004006a, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Cumulative curve, F(t)', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_01_03.png")
grid.raster(img)
```

```{r fig004006b, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Washout curve, W(t)', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_01_04.png")
grid.raster(img)
```

(c) The mean residence time is obtained from \@ref(eq:ec004007), thus $\overline{t} = 261.615$ seconds.

(d) To calculate the fraction of material leaving the reactor that has resided in the reactor between 230 and 270 seconds, it is necessary to solve

 \begin{equation*}
   \int_{230}^{270} {\mathbf{E}(t)dt}
 \end{equation*}
 
 However, times 230 and 270 seconds are not presented in the original data. It is necessary to append our data by interpolating _E(t)_ points to correspond times 230 and 270 seconds. Then, the data obtained are plotted in Figure \@ref(fig:fig004007). The shaded area represents the fraction of material leaving the reactor that resided in the reactor between 230 and 270 seconds.
 
```{r fig004007, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Analysis of curve E(t)', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_01_05.png")
grid.raster(img)
```

Then, we learned between 230 and 270 seconds, the 37.53\% of the material left into the reactor.

## Determination of **E(t)** from Step Change

In this case, a step input tracer is a step increase from one steady-state value to another. Usually, as illustrated in Figure \@ref(fig:fig004004), the step increase is from a zero value, and for a $t\to \infty$, the final concentration of the tracer will be equal to its inlet concentration. [@hayes2013]

For a step change, a material-balance criterion is that the steady-steady inlet and outlet tracer concentrations must be equal, both before and after the step change. Then, it may be concluded that the system's response, as shown in Fig. \@ref(fig:fig004004), is linear to the tracer, and there is no loss of tracer because of reaction or adsorption. [@missen1999]


```{r fig004004, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Typical step change curves for RTD measurement. [@fogler2008]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_04.png")
grid.raster(img)
```

In comparison with pulse input, which is related to $\mathbf{E}(t)$, the response to a step change is related to $\mathbf{F}(t)$. The advantages of this experimental method are:

 1. A step change is usually easier to achieve.
 2. A material balance is usually easier to achieve.

While disadvantages include:

 1. Continuous delivery requires a greater amount of tracer.
 2. It may have a significant impact on process operation.
 
Now, considering a constant rate of tracer addition to a feed that is initiated at time $t=0$, the tracer concentration for a system with a constant volumetric flow rate must achieve following conditions:
 \begin{equation*}
   C_0(t) = \left\{ \begin{array}{ll}
   0 & t <0 \\
   {\left(C_0\right)} \,\text{ constant} & t \geq 0 \end{array} \right.
 \end{equation*}
 
Thus, in the case of a step increase from an initial value, $C(0)$, to some final value, $c(\infty)$, the outlet concentration as function of time gives the cumulative distribution function, $\mathbf{F}(t)$, directly from the following equation:
 \begin{equation}
   \mathbf{F}(t) = \frac{c(t) - c(0)}{C(\infty) - c(0)}
   (\#eq:ec004022)
 \end{equation}
 
If the inlet tracer concentration has a step decrease from initial value, $C(0)$, to a final value, $C(\infty)$, then the effluent concentration can be used to calculate the washout function, $\mathbf{W}(t)$:
 \begin{equation}
\mathbf{W}(t) = \frac{c(t) - c(\infty)}{c(0) - c(\infty)}
(\#eq:ec004023)
 \end{equation} 
 
A tracer study may use a step increase followed later by a step decrease; the transient responses in the two cases are the checked consistency. When considered separately, the washout technique has advantages: less tracer is required, and it avoids having to maintain a steady-state value of $C_o$ for a lengthy period. [@missen1999]
 
The normalized response data may be converted to $\mathbf{E}(t)$, since
 \begin{equation*}
   \mathbf{E}(t) = \frac{d\mathbf{F}(t)}{dt}
 \end{equation*}
 Or could be estimated depending on the differencing technique used. From backward differencing:
  \begin{equation*}
   \mathbf{E}(t) = \frac{d\mathbf{F}(t)}{dt} \approx \frac{F_{i+1} - F_{i}}{t_{i+1} - t_{i}}; \quad i = 1, 2, \,\dots \, n-1
 \end{equation*}
 and central differencing:
   \begin{equation*}
   \mathbf{E}(t) = \frac{d\mathbf{F}(t)}{dt} \approx \frac{F_{i+1} - F_{i-1}}{t_{i+1} - t_{i-1}}; \quad i = 2, 3,\, \dots \, n-2
 \end{equation*}
 
Finally, based on the normalized condition of $\mathbf{E}$ curve. The mean residence time could be approximated to:
 \begin{equation*}
   \overline{t} \approx \sum{t_i\,\mathbf{E}_i\,\Delta t}
 \end{equation*}
 and the variance $\sigma^{2}$
 \begin{equation*}
   \sigma_t^{2} \approx \sum{t_i^2\,\mathbf{E}_i\,\Delta t} - \overline{t}^{2}
 \end{equation*}

```{exercise, it004002, name="Calculations of RTD and distribution functions curves from step change method"}
The inlet concentration of helium (tracer A) was step-increased from 1.0 to 2.0 mmol/L to determine the mixing pattern in a fluidized-bed reactor. The response data were as follows:
  
 t (min) | C (mmol/L) | t (min) | C (mmol/L) 
--|--|--|--
0.0  | 1.000 | 30.0  | 1.410 
5.0  | 1.005 | 45.0  | 1.610 
10.0 | 1.020 | 60.0  | 1.770 
15.0 | 1.060 | 90.0  | 1.920
20.0 | 1.200 | 120.0 | 1.960

Determine $F(t)$, $E(t)$, $\sigma^2$ for flow through the vessel, and calculating $\overline{t}$ 

```


_Solution_

The solution in a Python script is available at  [link](https://github.com/aliglara/ebook_che5314_ucv/blob/main/codes/lect_022RTDstep.ipynb).

The figure \@ref(fig:fig004008) shows results obtained for $F(t)$

```{r fig004008, out.width = "80%", echo=FALSE, fig.align='center', fig.cap='Left: Trace concentration curve. Right: F(t) curve', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_02_01.png")
grid.raster(img)
```

 
The RTD was calculated by
 \begin{equation*}
   \mathbf{E}(t) = \frac{d\mathbf{F}(t)}{dt}
 \end{equation*}
 
 Which was estimated using backward differencing for $F(t)$. The figure \@ref(fig:fig004009) shows results obtained for $E(t)$
 
```{r fig004009, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='E(t) curve', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_02_03.png")
grid.raster(img)
```

The mean residence time was estimated by trapezoid rule.
\begin{equation*}
  \overline{t} = \sum\limits_{i=1}^{n_c - 1} {\left[t_i\,E_i(t)\right]}_{av} \Delta t_i
\end{equation*}
From this,

\begin{equation*}
 \overline{t} = 57.537 \text{ min}
\end{equation*}

Finally, the variance was obtained by,
\begin{equation*}
  \sigma^2 = \left\{\sum\limits_{i=1}^{n_c - 1} {\left[t_i^2\,E_i(t)\right]}_{av} \Delta t_i \right\} - \overline{t}^2
\end{equation*}

From this,
\begin{equation*}
  \sigma^2= 530.249 \text{ min}^2
\end{equation*}

# Applications for RTD functions

As it was said before, RTD provides a measure of macroscopic mixing; therefore, it could be helpful as a tool for:

1. detecting and characterizing flow behavior.
2. estimating of values of parameters for nonideal flow models.
3. assessment of the performance of a vessel as a reactor.

## RTD functions for ideal reactors
$\mathbf{F}(t)$ fits the behavior either when $N$ particles are injected individually at different times to a steady-state reactor or if they are entered all at once to the reactor as $\mathbf{F}(t)$ represents the volume fraction of the fluid leaving the reactor, which resided in the system over a time of $t$ or less.

### Plug Flow Reactor (PFR) {-}
In this reactor, there is not back mixing. Hence, every fluid element has the same mean residence time $\overline{t}$ (or $\tau$). This corresponding to a pulse tracer entering the reactor at time $t=0$ during an infinitesimally short time and does not mix with fluid elements running ahead or behind it within the reactor. 
  
The flow pattern in PFR could be represented in terms of the Dirac delta function $\left( \delta \left(t - \overline{t} \right)\right)$.The Dirac delta function is a generalized function, on the real number line that is zero everywhere except at zero, where it is infinite:

\begin{equation*}
\delta \left(t - \overline{t} \right) = \left \{
\begin{aligned}
  & \infty && t = \overline{t} \\
  & 0  && t \neq \overline{t}
\end{aligned}
\right.
\end{equation*}

Therefore, the normalized RTD function is written down in terms of Dirac "delta-function''
\begin{equation}
\mathbf{E}(t) = \delta (t - \overline{t} ) = \left \{
\begin{aligned}
  & \infty \text{ when } t = \overline{t} \\
  & 0 \text{ when } t \neq \overline{t}
\end{aligned}
\right.
(\#eq:ec004025)
\end{equation}

The integral of Dirac $\delta$-function is given as:

\begin{equation*}
\int_a^b \delta {\left(t-\overline{t}\right)} \, f(t) dt = \left \{
\begin{aligned}
& f(\overline{t}) \text{ when } a < \overline{t} < b \\
& 0 \text{ when } a > \overline{t} > b
\end{aligned}
\right.
\end{equation*}
where $a$ and $b$ are any constant.

Therefore,
\[
\mathbf{F}(t) = \left \{
\begin{aligned}
  & 0 \text{ when } 0 < t < \overline{t} \\
  & 1 \text{ when } t \geq \overline{t}
\end{aligned}
\right.
\]

:::: {.infobox .info data-latex="information"}

The moments of the normalized distribution all equal one, while the central moments equal zero.

::::

### Continuous Stirred Tank Reactor (CSTR) {-}

An ideal CSTR implies that as soon as a fluid element (or pulse of tracer) is injected into the reactor, it distributes itself uniformly throughout the reaction zone. In "ideal'' CSTRs, the exit stream is assumed to have the same composition as the reaction mixture within the reaction zone.

It follows that the exit stream will have the highest possible concentration of material from the tagged fluid element (or pulse tracer) at the outset of the experiment (i.e., $t =0$). After that, the tracer concentration in the main reactor would be diluted with the inlet stream, and the concentration of tracer in the exit stream would be expected to decay to zero gradually.

The RTD for a CSTR can be easily analyzed using a pulse injection. Assume, in the first instance, that the inlet stream to the vessel contains a concentration of the tracer equals to $C_{TR,0}$, given by

\begin{equation}
  C_{TR,0} = \frac{N_{TR}}{V_R}
  (\#eq:ec004026)
\end{equation}

The decay of tracer concentration in the reactor is an unsteady state process given by,
\begin{equation}
  -\upsilon\,C_{TR} = V_R\,\frac{dC_{TR}}{dt}
  (\#eq:ec004026a)
\end{equation}

For a constant-density fluid, Eq. \@ref(eq:ec004026a) could be rewritten as

\begin{equation}
  -\frac{\upsilon}{V_R}\,dt=-\frac{dt}{\overline{t}} = \frac{dC_{TR}}{C_{TR}}
  (\#eq:ec004027)
\end{equation}

Integrating Eq. \@ref(eq:ec004027) yields to:
\begin{equation}
C_{TR}(t) = C_{TR,0}\,\exp{\left(-\frac{t}{\overline{t}}\right)}
(\#eq:ec004028)
\end{equation}


Now, substituting Eq. \@ref(eq:ec004028) into Eq. \@ref(eq:ec004023)
 \begin{equation}
\mathbf{W}(t) = \frac{C_{TR}(t) - C_{TR}(\infty)}{C_{TR,0} - C_{TR}(\infty)} = \cfrac{C_{TR,0}\,\exp{\left(-\cfrac{t}{\overline{t}}\right)}
}{C_{TR,0}} = \exp{\left(-\frac{t}{\overline{t}}\right)}
(\#eq:ec004029)
 \end{equation} 
 
And, the RTD function is obtained by differentiation:
 \begin{equation}
\mathbf{E}(t) = -\frac{d\,\mathbf{W}(t)}{dt}  =\frac{1}{\overline{t}} \exp{\left(-\frac{t}{\overline{t}}\right)}\equiv \frac{C(t)}{\displaystyle\int_0^\infty C(t)\, dt}
(\#eq:ec004030)
 \end{equation}
 
 The RTD decays exponentially with time. However, the area under the curve is always equal to unity because of the normalization. The typical residence time distribution for a CSTR is presented in Figure \@ref(fig:fig00406)

```{r fig00406, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='The residence time distribution E(t) in a CSTR and the time interval from t to t+dt', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_07.png")
grid.raster(img)
```

Other tracer experiments can be used to characterize flow patterns for non-ideal reactors. However, the input of a pulse or a step-change is commonly used in practice. Figures \@ref(fig:fig00407a) and  \@ref(fig:fig00407b)  summarize the exit concentration curves and the shapes of the $\mathbf{E}(t)$ curves for a step change and input pulse of tracer, respectively.

```{r fig00407a, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Response of tracer species during a step change experiment', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_08.png")
grid.raster(img)
```


```{r fig00407b, out.width = "55%", echo=FALSE, fig.align='center', fig.cap='Response of tracer species during an input impulse experiment', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_10.png")
grid.raster(img)
```

## RTD as diagnostic tool for characterizing flow behavior

The RTD curve can be used as a diagnostic method to determine the characteristics of reactor flow patterns. The methods of interpretation fall into two main categories:

 1. processes with a relatively small degree of mixing, and
 2. processes with significant deviations from ideal flow patterns, namely, dead space, bypassing, and general nonuniform regions.

Since these maldistributions can cause unpredictable conversion in reactors, they are usually detrimental to reactor operation.

### Relatively Small Degree of Mixing {-}

The age-distribution functions for flow patterns with a relatively small degree of mixing will have the general shapes shown in Figure \@ref(fig:fig004agedistribution)

```{r fig004agedistribution, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Effect of some features of nonideal flow on RTD. [@himmelblau1968]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_29.png")
grid.raster(img)
```

As is shown, the peak of the $\mathbf{E}(\theta)$ curve will almost be at the mean holding time $\theta = \overline{t}/\overline{t} = 1$, and $\mathbf{I}(\theta)$ curve will have a value of about 0.5 at $\theta = 1$. None of the curves will have excessively long tails and, in particular, the $\mathbf{E}(\theta)$ will be almost symmetrical.

These features lead to the use of the type of description based on Gaussian curves, namely, the mean and the variance.
\begin{align*}
  \text{mean} & = \mu_1 = \int_0^\infty \theta\,\mathbf{E}(\theta)\,d\theta = 1 \\
  \text{variance} & = \sigma^2 = \int_0^\infty \left(\theta - \mu_1 \right)^2\,\mathbf{E}(\theta)\,d\theta
\end{align*}

### Detection of Dead Space {-}

A process vessel may exhibit a region where fluid elements are retained for times of an order of magnitudes more significant than the mean residence time of the total fluid, what is called **dead space**. The existence of dead space is most easily verified from the characteristics of the $\mathbf{E}(\theta)$ curve as shown in Figure \@ref(fig:fig004deadspace). The $\mathbf{E}(\theta)$ will have a very long tail indicating that fluid is held in the dead space.

```{r fig004deadspace, out.width = "55%", echo=FALSE, fig.align='center', fig.cap='Identification of the existence of dead space in a process vessel. [@himmelblau1968]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_30.png")
grid.raster(img)
```

The _true_ mean residence time will be equal to the vessel volume divided by the flow rate $\overline{t} = V/\upsilon$. However, the experimental data for times longer than $\theta$ equal to 2 or 3 times mean residence time are seldom of sufficient accuracy to use in calculations. Thus, the calculated apparent mean residence time would be
\begin{equation*}
  \overline{\theta}_a = \int_0^\theta \theta^\prime \, \mathbf{E}(\theta^\prime)\,d\theta^\prime, \qquad \theta \approx 2 \text{ or } 3
\end{equation*}

whereas the true mean would be
\begin{align*}
 \theta & = \int_0^\infty \theta \, \mathbf{E}(\theta)\,d\theta = \int_0^\theta \theta^\prime \, \mathbf{E}(\theta^\prime)\,d\theta^\prime + \int_\theta^\infty \theta^\prime \, \mathbf{E}(\theta^\prime)\,d\theta^\prime = 1 \\
  & = \overline{\theta}_a + \int_\theta^\infty \theta^\prime \, \mathbf{E}(\theta^\prime)\,d\theta^\prime = 1
\end{align*}

In presence of dead space, $\overline{\theta}_a < 1$ or $\overline{t}_a < \overline{t}$. The apparent mean gives a measure of the dead volumen of the vessel:
\begin{equation*}
  \overline{t}_a = \frac{V_\text{active}}{\upsilon_\text{active}} = \frac{V - V_\text{dead}}{\upsilon_\text{active}} = \overline{t}\, \overline{\theta}_a
\end{equation*}

This equation can be rewritten as
\begin{align*}
  \frac{V_\text{dead}}{V} &= \frac{\upsilon_\text{active}}{V} \left(\frac{V}{\upsilon_\text{active}} - \overline{t}\overline{\theta}_a\right) \\
    \frac{V_\text{dead}}{V} &= 1 - \overline{\theta}_a\,\overline{t}\,\frac{\upsilon_\text{active}}{V} \\  
    \frac{V_\text{dead}}{V} &= 1 - \overline{\theta}_a \left(\frac{\upsilon_\text{active}}{\upsilon} \right)\approx 1 - \overline{\theta}_a
\end{align*}

One problem with using these relations is that they assume the true mean holding time, $\overline{t}$, is known. Since it cannot be found from the trace curve (due to the inaccurate tail), it must be known independently and challenging to obtain.

### Detection of Bypassing {-}
If a portion of fluid spends one-tenth of the overall fluid's mean residence time, this fluid can be said to bypass the vessel for all practical purposes. The $\mathbf{E}(\theta)$ curve would then appear as in Figure \@ref(fig:fig004bypass).

```{r fig004bypass, out.width = "55%", echo=FALSE, fig.align='center', fig.cap='Identification of the existence of bypassing in a process vessel. [@himmelblau1968]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_31.png")
grid.raster(img)
```

If there is a large amount of dead space or bypassing, the distinction between the two is not clear-cut but depends on which part of the fluid is considered the "major'' part, a somewhat arbitrary assignment.

However, $\mathbf{E}(\theta)$ can be used to detect dead space or bypassing (as discussed above). If the true mean holding time is not known, this curve cannot be used to determine dead space. Instead, the $\mathbf{I}(\theta)$ curve must be used because it always has the same shape as regards the existence of maxima no matter what the time scale. Therefore, the most significant utility of the intensity function is to determine dead space from experimental data, especially if the true mean holding time is unknown.

Figure \@ref(fig:fig004002a) summarizes the typical shapes of the $\mathbf{E}(\theta)$ curve of (a) stagnation, (b) dispersion and (c) channeling or excessive bypassing.

```{r fig004002a, out.width = "55%", echo=FALSE, fig.align='center', fig.cap='Effect of some features of nonideal flow on RTD. [@missen1999]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_24.png")
grid.raster(img)
```

If large amounts of dead space of bypassing exist, it may not be worthwhile to try to construct a detailed mathematical model because the unit is sure to give poor performance. Thus, some modifications might be called for to eliminate the severe mixing problems.

In order to provide a better picture of the behavior of the distribution functions, a simple mathematical model will be considered that includes the possibility of having well-defined dead space and bypassing. The model consists of two parallel streams, each with a given number of perfectly mixed tanks in series as shown in Figure \@ref(fig:fig004modeldeadzone)

```{r fig004modeldeadzone, out.width = "55%", echo=FALSE, fig.align='center', fig.cap='Scheme of the mathematical representation. [@himmelblau1968]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_32.png")
grid.raster(img)
```

The RTD for this model with $j$ stirred tanks in series is
\begin{equation*}
\mathbf{E}(t) = \frac{j^j}{\left(j-1\right)!} \frac{t^{j-1}}{\overline{t}}\exp\left(-\frac{j\,t}{\overline{t}}\right)
\end{equation*}
where $\overline{t}$ is the mean residence time. If the equation is used for each branch of the model and weighted with respect to the amount of flow in the branch, the RTD for the entire model is:
\begin{equation*}
\mathbf{E}(t) = f\,\left(\frac{n^n}{\left(n-1\right)!} \frac{t^{n-1}}{\overline{t}_1}\exp\left(-\frac{n\,t}{\overline{t}_1}\right)\right) + \left(1-f\right)\,\left(\frac{m^m}{\left(m-1\right)!} \frac{t^{m-1}}{\overline{t}_2}\exp\left(-\frac{m\,t}{\overline{t}_2}\right)\right)
\end{equation*}
where $f$ represents the fraction of flow to branch **1**. The mean residence time is
\begin{equation*}
  \overline{t} = f\, \overline{t}_1 + \left(1 - f\right)\, \overline{t}_2
\end{equation*}

Dimensionless variables based on the total mean residence time, $\overline{t}$ can now be defined:
\begin{equation*}
  \theta = \frac{t}{\overline{t}} \quad \alpha = \frac{\overline{t}_2}{\overline{t}_1} \quad \beta = f + \left(1 - f\right) \alpha
\end{equation*}

Three situations were considered:

1. Fifteen tanks in each branch correspond to a relatively small amount of mixing in the branch.
 2. Two tanks in each branch correspond to a large amount of mixing in the branch.
 3. Four tanks in each branch correspond to a moderate amount of mixing in the branch.  

Scenario 1. is shown in Figure \@ref(fig:fig004RTDdiagnostic). The rest of scenarios is available clicking on the next   [link](https://github.com/aliglara/ebook_che5314_ucv/blob/main/codes/lect_023RTDasDiagnosticTool.ipynb)

```{r fig004RTDdiagnostic, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='RTD curves for parallel stirred-tanks models.', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_07.png")
grid.raster(img)
```

# Models for nonideal flow

The information provided by RTD functions is essential to know if a reactor has dead zones or bypassing. These situations affect the reactor performance. Thus, the main target is answering the question "How much bypassing or dead volume is acceptable?".

For answering, the chemical engineer needs to calculate the performance of the reactor. Furthermore, to do so, two types of information are required: the flow patterns and the kinetic model.

```{r out.width = "65%", echo=FALSE, fig.align='center'}
library(png)
library(grid)
img <- readPNG("figures/004_modelnonidealstructure.png")
grid.raster(img)
```

The information provided by RTD functions is essential to know if a reactor has dead zones or bypassing. These situations affect the reactor performance. Thus, the main target is answering the question "How much bypassing or dead volume is acceptable?''.

For answering, the chemical engineer needs to calculate the performance of the reactor. Furthermore, to do so, two types of information are required: the flow patterns and the kinetic model.

The RTD curves found from tracer studies are only partial information on the details of the flow patterns because all that is known is the length of time-specific fractions of the fluid stay in the vessel (macromixing), but neither what of these fractions do while in the vessel nor in what geometric locations they reside, that is, information on micromixing.

As a result, several mathematical models of reactor performance have been developed to estimate the conversion levels in a nonideal reactor. [@hill1977]

The micromixing behavior has two extremes views:

1. Complete segregation: any fluid element is isolated from all other fluid elements and retains its identity throughout the entire vessel. No micromixing occurs, but macromixing may occur.
2. Complete dispersion: fluid elements interact and mix thoroughly at the microscopic level.

Besides the degree of segregation of fluid elements, another aspect of mixing is the relative time in passage through a vessel at which fluid element mix (''earliness'' of mixing). [@missen1999]

In general, a mathematical model for nonideal flow in a vessel provides a characterization of mixing and flow behavior. Although it may appear to be an independent alternative to the experimental measurement of RTD, the latter may be required to determine the parameter(s) of the model.

For that reason, these mathematical models can be classified in function of the number of adjustable parameters. [@fogler2008]

1. Zero adjustable parameters
	* segregation model
	* maximum mixedness model
2. One parameter model
	* tanks-in-series model
	* dispersion model
3. Two adjustable parameters
	* actual reactor modeled as combinations of ideal reactors

These parameters are correlated as functions of fluid and flow properties, reactor configurations, along with other vital features, and can be used in design calculations.

Generally speaking, the more adjustable parameters used in a model, the better the correlation between predicted and observed RTD would be. These adjustable parameters are often related to RTD moments and are usually expressed in terms of dimensionless moments about the mean.

The usual approach is developing a model flow system with the same RTD as the existing system. Provided that the conservation equations can be written for this model system, the performance of the actual vessel can be predicted.

## Zero adjustable parameters
### Segregational or macrofluid model {-}

This model represents the micromixing condition of complete segregation (no mixing) of fluid elements, which means that fluid flows through a vessel with no mixing between adjacent fluid elements. The feed enters the reactor as little _packets_ of fluid. Each packet retains its identity as they pass through the vessel; there is no mass exchange between individual packets. However, packets can mix each other in the reactor. Thus those that enter at some time, _t_, will not all leave at the same time.
  
Each packet can be treated as a small ideal batch reactor. A reaction (or reactions) takes place as the tiny packets move through the vessel. The composition of each packet changes while it is flowing through the vessel. However, its final conversion will depend on how long the packet has spent in the vessel. [@roberts2009].

After leaving the reactor,  all fluid packets are mixed on a molecular level, and the average composition is measured. This situation is represented in Figure \@ref(fig:fig004010c). [@roberts2009]

```{r fig004010c, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Segregational model scheme. [@roberts2009]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_18.png")
grid.raster(img)
```

When the packets are admitted to a CSTR, we can envisage two scenarios that involve perfect mixing: one with perfect mixing at the macroscale and the other with perfect mixing at the microscale, as is shown in Figure \@ref(fig:fig004010d)

```{r fig004010d, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Different levels of mixing in a CSTR: (a) Segregational model, (b) Ideal CSTR model. [@hayes2013]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_19.png")
grid.raster(img)
```

The conversion in a segregated tank reactor can be calculated using the ideal batch reactor equation design and the RTD curve. Each batch reactor has its conversion, which depends on its time in the tank. While the RTD gives information about the mean residence time.

The initial concentration of reactants in the packet is the same as the inlet concentration of the reactor. While each fluid element  (assuming constant density) behaves as a batch reactor, the total reactor conversion is the average conversion of all the fluid elements at the outlet conditions. That is:

```{r, echo=FALSE, fig.fullwidth=TRUE}
library(png)
library(grid)
img <- readPNG("figures/004_segregationalmodelbalance.png")
grid.raster(img)
```

where the summation is over all fluid elements in the reactor exit stream. This equation can be written analytically as:

\begin{equation}
  \overline{C}_A = \int_0^\infty{C_A(\overline{t})\,\mathbf{E}(t)\,d\overline{t}}
(\#eq:ec004031)
\end{equation}

Where $C_A(t)$ depends on the residence time of the element.

\begin{equation*}
  \frac{d\,C_A(t)}{dt} = - r_A = k\,C_A(t)
\end{equation*}

The equation \@ref(eq:ec004031) is rewritten as

\begin{equation}
  \frac{\overline{C}_A}{C_{A0}} = \int_0^\infty{\left(\frac{C_A(\overline{t})}{C_{A0}}\right)\,\mathbf{E}(t)\,d\,t} \approx \sum_{\substack{\text{\footnotesize{all age}}\\\text{\footnotesize{intervals}}}}\left(\frac{C_A}{C_{A0}}\right)\,\mathbf{E}(t)\,\Delta t
\end{equation}
(\#eq:ec004032)
or in terms of conversions
\begin{equation}
  \overline{X}_A = \int_0^\infty{{\left(X_A(t)\right)}_\text{\footnotesize{element}}\,\mathbf{E}(t)\,dt} 
(\#eq:ec004033)
\end{equation}

Introducing the $\mathbf{E}(t)$ curve for  CSTR, Equation \@ref(eq:ec004032) becomes,
\begin{equation*}
  \frac{\overline{C}_A}{C_{A0}} = \frac{1}{\overline{t}} \int_0^\infty{ \text{e}^{-k\,t}\,\text{e}^{-t/\overline{t}}\,d\,t} = \frac{1}{1 + k\,t} 
 \end{equation*}

Now, introducing the $\mathbf{E}(t)$ curve for a PFR,
\begin{equation*}
  \frac{\overline{C}_A}{C_{A0}} = \frac{1}{\overline{t}} \int_0^\infty{ \text{e}^{-k\,t}\,\delta \left(t - \overline{t}\right)\,d\,t} = \text{e}^{-k\,\overline{t}}
 \end{equation*}
 
The methodology described above can be used without any flow pattern models at all, in that experimental concentration as a function of residence time and the measured RTD can be directly numerically integrated to predict the conversion in the flow reactor.

:::: {.infobox .info data-latex="information"}
The segregational o macrofluid model is essential because it permits reactor behavior to be estimated directly from RTD, $\mathbf{E}(t)$, and the reaction kinetics. No other information is necessary.
::::

However, it is only strictly applicable to isothermal, single-phase systems. Extensions to more complicated behaviors are not straightforward. Therefore, other techniques are required for more general predictive and design purposes.

When we dealt with ideal reactors in series, we learned that:

1. If the order of the reaction is greater than 1 $\left(n > 1\right)$, we must keep the reactant concentration as high as possible for as long as possible to maximize fractional conversion. In order words, we must avoid mixing as long as possible. For example, the optimum arrangement will be a PFR first, then a small CSTR, and a large CSTR.
2. If the order of the reaction is less than 1 $\left(n < 1\right)$, we must drop the reactant concentration as low as possible, which means mix as soon as possible. 
3. Finally, if the order of the reaction is equal to 1, the order of reactors will not affect the final conversion.
    
The macrofluid model represents the latest possible mixing for a given RTD because there is no mixing between fluid elements until the reaction is over. For a reaction with an effective order greater than 1, the macrofluid model represents the _best_ possible situation. It provides an upper bound on conversion. [@roberts2009]

Conversely, for a reaction with a practical order of less than 1, the macrofluid model represents the \textit{worst} possible situation. It provides a lower bound on conversion. In summary,

Reaction order | Segregational model provide
-|-
$>1$ | Upper bound of conversion 
$=1$ | Exact result 
$<1$ | Lower bound of conversion 

```{exercise, it004003, name="Calculation of conversion based on RTD"}

Calculate conversion for a reaction 
\begin{equation*}
A \rightarrow B 
\end{equation*}

for a given RTD and segregation model.

t (min) | c (g/m$^3$) | t (min) | c (g/m$^3$)
- | - | - | -
1 | 1 | 7 | 4
2 | 3 | 8 | 3
3 | 8 | 9 | 2.2
4 | 10 | 10 | 1.5
5 | 8 | 12 | 0.6
6 | 6 | 14 | 0

\begin{equation*}
  r_A = - k\,C_A; \quad k = 0.1 \text{ min}^{-1}
\end{equation*}
  
```

_Solution_

The solution in a Python script is available at  [link](https://github.com/aliglara/ebook_che5314_ucv/blob/main/codes/lect_024RTDconversion.ipynb).

From data problem, the RTD is shown on the left size of Figure \@ref(fig:fig004011)

```{r fig004011, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Graphical solution for the problem. Left: E(t) curve, Right: Fractional conversion curve.', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_03s.png")
grid.raster(img)
```

From kinetic data, the fractional conversion for each batch reactor is
\begin{equation*}
  X_{i,\text{batch}}(t) = 1 - \exp(-k\,t)
\end{equation*}
The profile of conversion is plotted right size of  Figure \@ref(fig:fig004011)

Finally, the mean conversion is estimated by Eq. \@ref(eq:ec004033)
\begin{equation*}
  \overline{X}_i = \int_0^\infty {\left[1 - \exp(-k\,t) \right]\,\mathbf{E}(t)\,dt} = 0.3871 
\end{equation*}

### Maximum Mixedness Model (MMM) {-}

The MMM for a reactor represents the micromixing condition of complete dispersion, where fluid elements are entirely mixed at the molecular level, which is the opposite extreme of the segregational model.

The MMM considers that the fluid packages of different residence times are mixed as early as possible in the flow system consistent with the observed RTD. For this reason, this model is also called **_Early mixing model_**.
 
This behavior can be modeled using a PFR with inlet side streams, as shown in Figure \@ref(fig:fig004010e). The feed flow rate is distributed among the side streams in such a way as to match the required RTD function. It can be seen that in this arrangement, fluid entering the reactor in a given side stream, which has a residence time of zero, will mix radially with the fluid already present in the reactor completely.

```{r fig004010e, out.width = "55%", echo=FALSE, fig.align='center', fig.cap='Scheme for modeling of maximum mixedness model. [@hayes2013]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_20.png")
grid.raster(img)
```

The conversion in a MM model can be calculated by performing a mole balance that considers all side streams. First, it is necessary to define an alternative variable, $\lambda$, defined as the residence time remaining for an entering particle. All particles at a given $\lambda$ have the same time left in the reactor but may have different ages.

```{r fig004010f, out.width = "55%", echo=FALSE, fig.align='center', fig.cap='Control volume used to develop the mole balance for the maximum mixedness. [@hayes2013]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_21.png")
grid.raster(img)
```

The concentration at the reactor outlet can be derived from mole balance equation. Take a small control volume of the reactor, as show in Figure \@ref(fig:fig004010f). The mole balance over the discrete volumen element is
\begin{equation}
q\,C_A\left(\lambda + \Delta \lambda\right) + \Delta q\,C_{A0} + \Delta V \, r_A = \left( q + \Delta q \right)\,C_A (\lambda)
(\#eq:ec004033a)
\end{equation}

The volumetric flow rate at $\lambda$ is the inlet volumetric flow, $q$, plus the inlet side streams, $\Delta q$, thus
\begin{align*}
\Delta q &= \upsilon_0 \, \mathbf{E}(\lambda) \, \Delta \lambda \\
q(\lambda) &= q(\lambda + \Delta \lambda) + \upsilon_0 \, \mathbf{E}(\lambda) \, \Delta \lambda
\end{align*}
where $\upsilon_0$ is the total inlet volumetric flow rate fed to the reactor.

Rewriting
\begin{equation}
\frac{d\,q(\lambda)}{d\lambda} = - \upsilon_0\,\mathbf{E}(\lambda)
(\#eq:ec004033b)
\end{equation}

The inlet volumetric flow rate, $\upsilon_0$ at the inlet of the reactor $(X_A = 0)$ is equal to zero because the fluid is fed only for side streams. Thus, integrating equation \@ref(eq:ec004033b) considering that $q(\lambda) = 0$ at $\lambda = \infty$ and $q(\lambda) = q(\lambda)$ at $\lambda = \lambda$,

\begin{equation}
\boxed{q(\lambda) = \upsilon_0\, \int_\lambda^\infty \mathbf{E}(\lambda)\,d\lambda = \upsilon_0 \left[1 - \mathbf{F}(\lambda)\right]}
(\#eq:ec004033c)
\end{equation}
where $q(\lambda)$ refers to volumetric flow rate inside the reactor at $\lambda$

The volume of fluid with a life expectancy between $\lambda$ and $\lambda + \Delta \lambda$ is 
\begin{equation}
\Delta V (\lambda) = \upsilon_0 \left[ 1 - \mathbf{F}(\lambda)\right] \Delta \lambda
(\#eq:ec004033d)
\end{equation}

The rate of generation of A in this volume is
\begin{equation}
r_A\,\Delta V (\lambda) = r_A\,\upsilon_0 \left[1 - \mathbf{F}(\lambda)\right] \Delta \lambda
(\#eq:ec004033e)
\end{equation}

Now substituting Eqs. \@ref(eq:ec004033c), \@ref(eq:ec004033d) and \@ref(eq:ec004033e) into \@ref(eq:ec004033a)

\begin{align}
\upsilon_0 \left[1 - \mathbf{F}(\lambda)\right]\,C_A\left(\lambda  + \Delta \lambda\right) + \upsilon_0 \, \mathbf{E}(\lambda) \, \Delta \lambda\,C_{A0} & + r_A\,\upsilon_0 \left[ 1 - \mathbf{F}(\lambda)\right] \Delta \lambda = \cdots \\ \cdots &\left(\upsilon_0 \left[1 - \mathbf{F}(\lambda)\right] + \upsilon_0 \, \mathbf{E}(\lambda) \, \Delta \lambda \right)\,C_A (\lambda)
(\#eq:ec004033f)
\end{align}

Dividing by $\upsilon_0 \, \Delta \lambda$,
\begin{equation}
\frac{\mathbf{W}(\lambda)\,C_A\left(\lambda + \Delta \lambda\right)}{\Delta \lambda} + \mathbf{E}(\lambda) \,\left(C_{A0} - C_A\right) + r_A\, \mathbf{W}(\lambda) = \frac{\mathbf{W}(\lambda)\,C_A (\lambda)}{\Delta \lambda}
(\#eq:ec004033g)
\end{equation}

Now, taking the limit as $\Delta \lambda \to 0$ gives
\begin{equation}
-\frac{\mathbf{W}(\lambda)\,d\,C_A}{d\,\lambda} + \mathbf{E}(\lambda) \,\left(C_{A0} - C_A \right) + r_A\, \mathbf{W}(\lambda) = 0
(\#eq:ec004033h)
\end{equation}
or
\begin{equation}
\frac{d\,C_A}{d\,\lambda} = \frac{\mathbf{E}(\lambda)}{\mathbf{W}(\lambda)} \left(C_{A0} - C_A\right) + r_A
(\#eq:ec004033i)
\end{equation}

Equation \@ref(eq:ec004033i) can be written in terms of conversion as
\begin{equation}
-C_{A0}\,\frac{d\,X}{d\,\lambda} = \frac{\mathbf{E}(\lambda)}{\mathbf{W}(\lambda)}\,C_{A0}\,X + r_A
(\#eq:ec004033j)
\end{equation}
or
\begin{equation}
\frac{d\,X}{d\,\lambda} = -\frac{\mathbf{E}(\lambda)}{\mathbf{W}(\lambda)}\,X - \frac{r_A}{C_{A0}}
(\#eq:ec004033k)
\end{equation}

Equation \@ref(eq:ec004033j) is solved to determine the concentration at the reactor outlet. The equation usually requires a numerical solution. An initial condition is also required. In this case, the boundary condition is as $\lambda \to \infty$, then $C_A = C_{A0}$.

Despite this model, it may be impossible to determine the amount of mixing in a reactor and hence obtain an exact prediction of conversion, but it could be helpful to determine the bounds on the expected conversion. This bounding may be represented by the diagram shown in Figure \@ref(fig:fig004010g).

```{r fig004010g, out.width = "55%", echo=FALSE, fig.align='center', fig.cap='The bounding of the micro- and macromixing regions. [@hayes2013]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_22.png")
grid.raster(img)
```

The horizontal axis represents the extent of macromixing as defined by the RTD. The vertical axis represents the amount of micromixing, which may vary from zero to perfect mixing. A PFR is always a segregated flow reactor, while a CSTR can vary between perfect mixing and segregated flow.

## One adjustable parameter
For a first-order reaction or the flow-pattern fluid can be modeled as complete segregation or maximum mixedness, the zero-adjustable parameter models are enough to predict the conversion. However, for non-first-order reactions over fluids with good micromixing, more than just the RTD is needed.

A model of reactor flow patterns is necessary to predict the conversions and product distribution for such systems. We use combinations or modifications of ideal reactors to represent actual reactors to model these patterns.

However, in this course, we are dealing with two of them that are most common, each of which has one adjustable parameter: The first one is _the tanks-in-series model_ and the second is _the dispersion model_ which are widely used. 

In principle, both models can represent the flow pattern in a single vessel between the two extreme ideal models: CSTR and PFR. Additionally, they are roughly equivalent and could be applied to turbulent flow in pipes, laminar flow in very long tubes, flow in a packed bed, etc. However, the dispersion model is mainly used in catalytic reactions carried out within tubular reactors.

### The Tanks-in-Series Model {-}

The tank-in-series (TIS) model for arbitrary flow through a vessel assumes that it can be represented by flow through a series of _N_ equal-sized ideal CSTRs, where _N_ is the adjustable parameter of the model and describes the degree of mixing in the vessel.

The value of $N$ is tailored to match an observed or expected RTD. As $N$ increases from 1 to $\infty$, the flow represented in the vessel varies from back mix flow to plug flow. For practical purposes, a plug flow is virtually achieved for $N>30$. 

The volume of the vessel modeled, $V$, is equal to the sum of the volumes of the $N$ tanks. [@missen1999]
\begin{equation*}
  V = \sum N\,V_i
\end{equation*}

Consider, as an illustration, \@ref(fig:fig004012), three equally sized tanks in series through which a constant-density fluid flows, giving an equal space-time (a mean residence time) for each reactor. 

A pulse tracer is injected into the first tank at time zeroand, instantaneously, becomes uniformly distributed in the first tank, giving an initial concentration of $C_0$. The outlet concentration from the first vessel is given by:
\begin{align}
  V_1 \, \frac{dC_1}{dt} & = - \upsilon \, C_1 \\
  C_1(t) & = C_0 \exp\left(-\frac{t}{\overline{t_1}}\right)
(\#eq:ec004034)
\end{align}
where $\overline{t_i}$ is the mean residence time in a single tank. 

The mole balance equation for the second tank is obtained from the perfectly mixed mole balance, without reaction:
\begin{equation}
  V\frac{dC_2}{dt} = \upsilon\,{\left(C_1-C_2\right)}
(\#eq:ec004035)
\end{equation}

Substituting Eq. \@ref(eq:ec004034) and solving Eq. \@ref(eq:ec004035) using an integration factor $\exp(t/\overline{t}_2)$ along with the initial condition $C_2(0) = 0$:

\begin{equation}
C_2(t) = C_0\,\frac{t}{\overline{t_2}}\,\exp{\left(-\frac{t}{\overline{t_2}}\right)}
(\#eq:ec004036)
\end{equation}

```{r fig004012, out.width = "55%", echo=FALSE, fig.align='center', fig.cap='The tanks-in-series model uses a number of equally sized CSTRs placed in series', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_11.png")
grid.raster(img)
```

Performing the mole balance on the third vessel and solving in a similar way give the outlet concentration with time from the third reactor.
\begin{equation}
C_3(t) = \frac{C_0}{2}\,{\left(\frac{t}{\overline{t}_3}\right)}^2\,\exp{\left(-\frac{t}{\overline{t}_3}\right)}
(\#eq:ec004037)
\end{equation}

Therefore, generalizing this method to a series of $N$ CSTRs gives the RTD for $N$ CSTRs in series:
\begin{equation}
\frac{C(t)}{C(t_0)} = \frac{t^{N-1}}{(N-1)!\,\overline{t}_i^N}\,\exp{\left(-\frac{t}{\overline{t}_i}\right)}
\end{equation}
taking into account that $\overline{t} = N \: \overline{t_i}$ then
\begin{equation}
\frac{C(t)}{C(t_0)} = \mathbf{C}(t) = \frac{t^{N}}{(N-1)!}\,\left(\frac{t}{\overline{t}}\right)^{N-1}\exp{\left(-N\frac{t}{\overline{t}}\right)}
(\#eq:ec004038)
\end{equation}

Equation \@ref(eq:ec004038) could be expressed with normalized time using $\mathbf{C}(t) = \overline{t}\,\mathbf{E}(t)$, leading to the normalized distribution function is:
\begin{equation}
\mathbf{E}_N(\theta) = \frac{N^N}{(N-1)!}\,\theta^{N-1}\,\exp{\left(-N\theta\right)}
(\#eq:ec004041)
\end{equation}
where, again, $\theta = t/\overline{t}$. 

Figure \@ref(fig:fig004013) illustrates the RTDs of various numbers of CSTRs in series. As $N$ increases from 1, the spread of distribution decreases. Furthermore, when $N>2$, the peak would be an infinitive height (behavior of a PFR). 

The mathematical illustration is available on the following [link](https://github.com/aliglara/ebook_che5314_ucv/blob/main/codes/lect_025in-tank-serie.ipynb)

```{r fig004013, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='E(theta) curve using the tanks-in-series model', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_04_01.png")
grid.raster(img)
```

The cumulative RTD function $\mathbf{F}{(\theta)}$ could be calculated from equation \@ref(eq:ec004041), for integral values of $N$
\begin{equation}
\mathbf{F}_N(\theta) = \int_0^\theta \mathbf{E}(\theta)\:d\theta
(\#eq:ec004041a)
\end{equation}
The integral in equation \@ref(eq:ec004041a) can be evaluated analytically:
\begin{equation}
\boxed{\mathbf{F}_N(\theta) = 1 - \exp \left(-N\,\theta\right)\,\sum\limits_{i=0}^{N-1}{\frac{(N\theta)^i}{i!}}}
(\#eq:ec004041b)
\end{equation}

The result for equation \@ref(eq:ec004041a) is shown in Figure \@ref(fig:fig004014) for several values of $N$. For $N=1$, the result is that for a single CSTR, as $N$ increases, the shape becomes more like a step function, corresponding to a plug-flow.

```{r fig004014, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='F(theta) curve using the tanks-in-series model', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_04_02.png")
grid.raster(img)
```

Finally, the washout function, $\mathbf{W}{(\theta)}$, can be obtained from $\mathbf{F}{(\theta)}$
\begin{equation*}
  \mathbf{W}{(\theta)} = 1 - \mathbf{F}{(\theta)}
\end{equation*}

To utilize the tanks-in-series model, we need to determine the number of tanks that will give the same RTD as for the real system. This number may be obtained from the variance of RTD, which is given by the second moment about the mean:
\begin{equation*}
\sigma_{\theta}^2 = \int_0^\infty{{\left(\theta-1\right)}^2 \frac{\theta^{N-1}N^N}{(N-1)!}\,\exp{\left(-N\theta\right)}} = \int_{0}^{\infty} \theta ^2\mathbf{E}(\theta)\,d\theta - 1
\end{equation*}
Solution of the integral equation gives the result
\begin{equation*}
\sigma_{\theta}^2 = \frac{\left(N+2-1\right)!}{N\cdot N!} - 1 = \frac{\left(N+1\right)!}{N\cdot N!} - 1 = \frac{1}{N}
\end{equation*}
The variance of the RTD is seen to be equal to the inverse of the number of tanks used. Therefore, if the RTD curve is known, the appropriated number of tanks in series to use may be calculated, so
\begin{equation*}
  N = \frac{1}{\sigma_{\theta}^2} = \frac{\overline{t}^{ 2}}{\sigma^2}
\end{equation*}

```{exercise, it004004, name="Estimate number of tanks based in Tanks-in-serie model"}
A step increase in the inlet trace A concentration, from 1.0 to 2.0 mmol/L, was used to determine the mixing pattern inside the fluidized-bed reactor. The response data were as follows:

 t (min) | C (mmol/L) | t (min) | C (mmol/L) 
--|--|--|--
0.0  | 1.000 | 30.0  | 1.410 
5.0  | 1.005 | 45.0  | 1.610 
10.0 | 1.020 | 60.0  | 1.770 
15.0 | 1.060 | 90.0  | 1.920
20.0 | 1.200 | 120.0 | 1.960

Estimate values of $N$ based on both central and backward differencing and determine which technique best describes the tracer outflow concentrations.

```

_Solution_

The solution is available on the next [link](https://github.com/aliglara/ebook_che5314_ucv/blob/main/codes/lect_026numberofTanks.ipynb)

For each scenario, the equation \@ref(eq:ec004041b) was applied to calculate the estimated $\mathbf{F}(t)$, then based on Eq. \@ref(eq:ec004022), the output tracer concentration, $C_{out}(t)$, may be calculated from $\mathbf{F}(t)$.

Figure \@ref(fig:fig004015) shows the results of E(t) curve based on several numrical tecniques. On the other hand, Figure \@ref(fig:fig004016) shows the values of $C_{out}(t)$ predicted using forward and central differencing technique, and based on a curve fitting approach. The last one provides a superior representation of the experimental data.

```{r fig004015, out.width = "75%", echo=FALSE, fig.align='center', fig.cap='E(t) curves obtained based on several numerical tecniques applied to the tanks-in-series model', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/lect026_01.png")
grid.raster(img)
```

```{r fig004016, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Selection of the best technique for the tanks-in-series model', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/lect026_02.png")
grid.raster(img)
```

### The Dispersion Model {-}

The dispersion model describes nonideal tubular reactors based on the assumption that exists an axial dispersion of the material, governed by Fick's law diffusion-like, superimposed on the flow. 

```{r fig004017, out.width = "55%", echo=FALSE, fig.align='center', fig.cap='The dispersion model', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_12.png")
grid.raster(img)
```

Figure \@ref(fig:fig004017) shows a sketch of differences between velocity profiles between ideal PFR and actual flow. Fluctuations could be due to different flow and molecular velocities and turbulent diffusion. Although radial dispersion may also be significant, we consider only axial dispersion.

In considering axial dispersion as diffusive flow, it is assumed that Fick's first law applies, with the effective diffusion coefficient replaced by an axial dispersion coefficient, $D_L$. Thus, for unsteady-state behavior concerning a species $T$ (tracer), the molar flux ($N_T$) of $T$ at an axial position $x$ is

\begin{equation*}
  N_T = - D_L \frac{\partial C_T(x,t)}{\partial x}
\end{equation*}
where $C_T(x,t)$ is the concentration of $T$.

Consider a material balance for $T$ around the differential control volume shown in Figure \@ref(fig:fig004018). 

```{r fig004018, out.width = "55%", echo=FALSE, fig.align='center', fig.cap='Control volume for continuity equation for axial dispersion', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_13.png")
grid.raster(img)
```

We assume steady flow overall, but not for tracer ($T$); there is no place. 

```{r, out.width = "75%", echo=FALSE}
library(png)
library(grid)
img <- readPNG("figures/dispersion_balance.png")
grid.raster(img)
```

Under those assumptions, with $C_T \equiv C_T(x,t)$, $A_c$ and $D_L$ constant, for the transient flow of a non-reacting tracer, the mole balance equation is

\begin{equation}
  A_c\,\nu\,C_T - D_L\,A_c\,\frac{\partial C_T}{\partial x} - A_c\left[\nu C_T + \frac{\partial \left(\nu\,C_T\right)}{\partial x}\right] + D_L\,A_c\left[\frac{\partial C_T}{\partial x} + \frac{\partial}{\partial x} \left(\frac{\partial C_T}{\partial x}\right)\,dx\right] = A_c\,\frac{\partial C_T}{\partial t}\,dx
(\#eq:ec004042)
\end{equation}
where $\nu$ is the average axial velocity ($m/s$), $D_L$ is the axial dispersion coefficient, ($m^2/s$).

If system is density-constant and $\nu$ is constant, the equation \@ref(eq:ec004042) becomes to
\begin{equation}
- D_L\,\frac{\partial^2 C_T}{\partial x^2} + \nu\,\frac{\partial C_T}{\partial x} = \frac{\partial C_T}{\partial t}
(\#eq:ec004043)
\end{equation}
If we know the boundary conditions of the problem, solving this equation, we get the profile of tracer concentration vs. time. 

Equation \@ref(eq:ec004043) can be made dimensionless by defining the following scaling factors:
\begin{equation*}
  z = \frac{x}{L}; \qquad \theta = \frac{t}{\overline{t}} = t\,\frac{\nu}{L} \qquad C^\ast = \frac{C_T}{C_0}
\end{equation*}
Here $C_0$ is usually taken as the maximum concentration in the system.

\begin{equation}
\frac{\partial C^\ast}{\partial \theta} = \left(\frac{D_L}{\nu \, L}\right) \frac{\partial ^2 C^\ast}{\partial z^2} -\frac{\partial C^\ast}{\partial z} 
(\#eq:ec004044)
\end{equation}

and $D_L/(\nu\,L)$ is named the dispersion number, when
 \begin{align*}
&\frac{D_L}{\nu\,L} \to 0  && \rightarrow &&\quad \text{large dispersion, hence well-mixed flow} &&\\
&\frac{D_L}{\nu\,L} \to \infty && \rightarrow &&\quad \text{negligible dispersion, hence plug flow}&&
 \end{align*}

We also have a dimensionless group, denoted _Pe_, which is called **Peclet number** based on vessel length. The physical interpretation is that it represents the ratio of convective flux to diffusive flux:
 \begin{equation*}
   \text{Pe} = \frac{\nu\,L}{D_L} \equiv \frac{\nu\,C_T}{D_L(C_T/L)} = \frac{\text{rate of transport by convection}}{\text{rate of transport by diffusion or dispersion}}
 \end{equation*}

 So, we see that
 \begin{equation}
\frac{\partial C^\ast}{\partial \theta} = \frac{1}{\text{Pe}}\frac{\partial ^2 C^\ast}{\partial z^2} -\frac{\partial C^\ast}{\partial z} 
(\#eq:ec004045)
\end{equation}

Now, the results of solutions of Eq. \@ref(eq:ec004045) will depend on the choice of boundary conditions at inlet $(z = 0)$ and outlet $(z = 1)$. There are several possibilities concerning ''closed'' and ''open'' vessels, which do not necessarily satisfy actual positions.

When the flow is undisturbed as it passes the entrance and exit boundaries, we will call this the **open boundary condition**. While, when we have plug flow outside the vessel up to the boundaries, we call this the **closed boundary condition**. There are up to four combinations of boundaries, closed-closed, open-open, or mixed, based on previous conditions. Figure \@ref(fig:fig004020) illustrates the closed and open extremes.

```{r fig004020, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Boundary conditions used with dispersion model', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_14.png")
grid.raster(img)
```

For closed vessels, the flow conditions are illustrated schematically in Figure \@ref(fig:fig004021). $D_L=0$, or $Pe=\infty$ characterizes a plug-flow upstream of the vessel inlet. Since flow inside the vessel is dispersed, and the pulse injection is at inlet $(z=0)$, will have a $C_T$ discontinuity across the inlet. In contrast, a continuity on CT across the outlet $(z=1)$ is usually assumed, even though downstream flow $(z>1)$ is plug-flow.

```{r fig004021, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Flow conditions for closed-closed vessel', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_15.png")
grid.raster(img)
```

In this course, we only deal with a closed-closed case which is illustrated in Figure \@ref(fig:fig004021). Danckwerts formulated the closed configuration, and it means that there is no dispersion before $z=0$ and after $z=1$.

The equation \@ref(eq:ec004045) is a second-order boundary value problem requires two boundary conditions. If the system is taken as a closed system, the dimensionless boundary conditions are:
\begin{align}\label{ec:004_046}
  &C_\text{in}^\ast = C^\ast-\frac{1}{\text{Pe}}\,\frac{\partial C^\ast}{\partial z} && \text{at } z = 0 && \qquad\\
  & \frac{\partial C^\ast}{\partial z} = 0 &&\text{at } z = 1 && \qquad
  (\#eq:ec004046)
\end{align}
Additionally, the system requires an initial condition:
\begin{align}
  &C^\ast = 1  &&\text{at } \theta = 0, \text{ and } z = 0 && \qquad\\
  &C^\ast = 0 &&\text{at } \theta = 0, \text{ and } 0 < z \leq 1 && \qquad
  (\#eq:ec004047)
\end{align}

 An analytical solution to equation \@ref(eq:ec004045) with these boundary conditions is given by Otake and Kunigita (1598), but is rather complex, and it was omitted it here. However, moments and relatively simple an are given by:
 \begin{gather}
   \overline{\theta} = 1 \\
   \sigma_\theta^2 = \frac{2}{\text{Pe}^2}\left(\text{Pe}-1+\exp(-\text{Pe})\right)
   (\#eq:ec004048)
 \end{gather}
 An approximate closed form solutions for Pe > 50:
 \begin{gather*}
   \mathbf{F}(\theta) = \frac{1}{2} \left[1 - \text{erf}\left(\frac{1-\theta}{\sqrt{\theta}}\sqrt{\frac{\text{Pe}}{4}}\right)\right] \\
   \mathbf{E}(\theta) = \frac{\exp\left(-\text{Pe}(1-\theta^2)/4\right)}{\sqrt{4\,\pi/\text{Pe}}}
 \end{gather*}
 where $\text{erf}$ is the error function.
 
As a practical matter, these mathematical complexities should probably not be taken too literally since the precise conditions at the boundaries of actual equipment cannot usually be precisely defined in any event.


**Where to get the dispersion coefficient?**

  (1) **Experimental determination by analyzing pulse response**. The Peclet number can be found experimentally using a pulse injection into the tube inlet and finding the RTD function. The Peclet number depends on the variance of the normalized RTD function:
  
\begin{equation}
  \sigma^2_{\theta} = \frac{2}{\text{Pe}} - \frac{2}{\text{Pe}^2}\left[1 - \exp(-\text{Pe})\right]
  (\#eq:ec004049)
\end{equation}

(2) **Calculation from correlations**.
\begin{align*}
  & \frac{1}{\text{Pe}_R} = \frac{1}{Re \cdot Sc} + \frac{Re \cdot Sc}{192} && \text{ For } 1<Re< 2000 \text{ and } 0.23 < Sc < 1000\\
  & \frac{1}{\text{Pe}_R} = \frac{3\cdot 10^7}{Re^{2.1}} + \frac{1.35}{Re^{1/8}} && \text{ For } Re > 2000 
\end{align*}
where

\begin{align*}
 &\text{Pe}_R = \frac{wd}{D}    	&&:  \text{Peclet number in tube} &&\quad \\
 &Re  = \frac{wd}{\nu}  	&&:\text{Reynolds number} &&\quad \\
 &Sc = \frac{\nu}{D_m}    	&&:  \text{Schmidt number} &&\quad
\end{align*}

```{exercise, it004005, name="Comparison in conversion calculation for a non-ideal flow reactor"}
For a reaction

\begin{equation}
A \rightarrow B; \qquad    r_A = -k\,C_A; \quad k = 0.25 \text{ min}^{-1}
\end{equation}

Compare final conversions based on:

1. axial dispersion model with closed vessel approximations.
2. tanks in series model, and
3. ideal PFR model.

 t (min) | C (g/m<sup>3</sup>) | t (min) | C (g/m<sup>3</sup>) 
--|--|--|--
1.00  | 1.000 | 7.00  | 4.000 
5.00  | 2.000 | 8.00  | 4.000 
3.00  | 8.000 | 9.00  | 2.200 
4.00  | 10.00 | 10.0  | 1.500
5.00  | 8.000 | 12.0  | 0.600
6.00  | 6.000 | 14.0  | 0.000

```

_Solution_

The solution is available on the next [link](https://github.com/aliglara/ebook_che5314_ucv/blob/main/codes/lect_027nonidealmodels.ipynb)


Based on the axial dispersion model. We need to find the Pe number that is following RTD data.

In this case, we calculated from experimental data, the RTD curve, mean time and its variance. Then, the Eq. \@ref(eq:ec004049) was solved:
  \begin{equation*}
  \sigma^2_{\theta} = \frac{\sigma^2}{\overline{t}^2} = \frac{2}{\text{Pe}} - \frac{2}{\text{Pe}^2}\left[1 - \exp(-\text{Pe})\right] \Rightarrow Pe = 8.007
\end{equation*}

When a chemical reaction is included in the dispersion model, a numerical solution is usually required for the resulting differential equation. so, based on the reactions and considering a simple power-law rate expression, the mole balance of _A_ is written as

\begin{equation*}
  D_L\frac{\partial^2 C_A}{\partial x^2} - \nu \, \frac{\partial C_A}{\partial x} - (k\,C_A^{\,n}) = 0
\end{equation*}

Which it could be written in dimensionless variables as:
\begin{equation*}
  \frac{1}{\text{Pe}}\frac{\partial^2 C_A}{\partial z^2} - \nu \, \frac{\partial C_A}{\partial z} - (k\,\overline{t}\,C_A^{\,n}) = 0
\end{equation*}

And, in terms of the fractional conversion
\begin{equation*}
  \frac{1}{\text{Pe}}\frac{\partial^2 X_A}{\partial z^2} - \frac{\partial X_A}{\partial z} - k\,\overline{t}\,C_{A0}^{\,n-1}{\left(1-X_A\right)}^n = 0
\end{equation*}

The fractional conversion is therefore governed by the magnitude of three dimensionless groups, Pe, $n$, and $k\overline{t}\,C_{A0}^{n-1}$. There is an analytical solution to this equation only for a first-order reaction:
\begin{equation*}
  \frac{C_A}{C_{A0}} = 1 - X_A = \frac{4\,a\exp{\left(\text{Pe}/2\right)}}{{\left(1+a\right)}^2 \exp {\left(a\,\text{Pe}/2\right) }- {\left(1-a\right)}^2 \exp{\left(-a\,\text{Pe}/2\right)}}
\end{equation*}
where
\begin{equation*}
  a = {\left(1 + \frac{4\,k\,\overline{t}}{\text{Pe}}\right)}^{0.5}
\end{equation*}
Solving this equation $X_A = 0.6826$

For calculating the fractional conversion based in TIS model, it is necessary to find de number of reactor that is in accordance with given RTD data.
\begin{equation*}
  \text{number of reactors } = \frac{\overline{t}}{\sigma^2} = 4.4
\end{equation*}
Then the conversion is
\begin{equation*}
  X_A = 1 - \frac{C_{AN}}{C_{A0}} =  1 - \frac{1}{\left(1+k\,\overline{t}/N\right)^N} = 0.6736
\end{equation*}

For an ideal PFR. The fractional conversion is:
\begin{equation*}
  X_A = 1 - \exp(-k\,\overline{t}) = 0.7253
\end{equation*}

:::: {.infobox .info data-latex="information"}
  
In summary, **why is it relevant to perform an analysis of RTDs?**

This approach can ease the following tasks:
  1. Diagnostics of problems in reactor operations.
  2. Design of unit operations.
  3. Predictions of conversion and effluent concentrations.
  
::::

We have studied models of several simple, partially mixed reactors. These are shown in Figure \@ref(fig:fig004022) for increasing back mixing.

```{r fig004022, out.width = "65%", echo=FALSE, fig.align='center', fig.cap='Residence time distributions for the nonideal reactors. [@schmidt1998]', fig.lp='fig:'}
library(png)
library(grid)
img <- readPNG("figures/004_draw_23.png")
grid.raster(img)
```