# Dynamic mass balance in plug flow reactor

So I would like to ask you about a problem I encountered. The case study goes like this: I want to study species distribution inside a packed bed reactor, at first at steady state but eventually what interests me is the dynamic behavior of the reactor. The reactants are gases. For the smimle steady state case the mass balance of each species if defined by:
$$u_s\frac{\partial C_j}{\partial z}=\rho_b \sum_{k=1}^{N_{reac}} \nu_{j,k} R_{k}$$ where $u_s$ is the superficial velocity (m/s), $C_j$ the $j$ species concentration (mol/m^3), $z$ the axial direction, $\rho_b$ the catalyst density in the reactor ($\text{kg}_\text{cat}/\text{m}^3$), $\nu_{j,k}$ the $j$ species stoichiometric index in $k$ reaction and $R_k$ the $k$th reaction rate ($\text{mol}/(\text{kg}_\text{cat}*\text{s})$). Reaction rates are known functions of concentrations, temperature and pressure and superficial velocity is given by: $$u_s=\frac{Q}{A}=\frac{\sum_{j=1}^{N_c}F_j \frac{Rg T}{P}}{A}$$ where $Q$ is the volumetric flowrate ($\text{m}^3/\text{s}$), $A$ the cross-section area, $F_j$ the $j$th molar flowrate ($\text{mol}/\text{s}$), $Rg$ universal gas constant, $T$ is the temperature, $P$ is the pressure and $N_c$ is the number of components. While, when dealing with gases it is not correct to assume that velocity does not change along the reactor because of mole changes, in this case study I found (in terms of mass balance between reactors inlet and outlet) that the error produced is insignificant.
Now the easiest way to solve this is by discretizing the differential equation via finite difference scheme and you get a set of algebraic equations like:
$$u_s\frac{(C_{i,j}-C_{i-1,j})}{\delta z}=\rho_b \sum_{k=1}^{N_{reac}} \nu_{j,k} R_{i,k}$$ where $i$ is the spatial node number (not sure that's the correct terminology).
For the case where I want to study the dynamic behavior of the system the equation would be like:
$$\frac{\partial C_{i,j}}{\partial t}=-u_s\frac{(C_{i,j}-C_{i-1,j})}{\delta z}+\rho_b \sum_{k=1}^{N_{reac}} \nu_{j,k} R_{i,k}$$ or discretized: $$\frac{(C_{t,i,j}-C_{t-1,i,j})}{\delta t}=-u_s\frac{(C_{t,i,j}-C_{t,i-1,j})}{\delta z}+\rho_b \sum_{k=1}^{N_{reac}} \nu_{j,k} R_{t,i,k}$$ Please correct me if I am wrong

NOW, here is where my problem starts. The assumption of constant velocity (calculated at the inlet of the reactor) is not safe anymore. Explanation:
1. Velocity contant both in space and in time (calculated by $F$ at inlet of reactor at $t=0$): Not good. Gives large deviations in mass balance when inlet flowrate changes (checked it when the system reached steady state).
2. Velocity constant in space but not time (calculates by $F$ at inlet of reactor at each $t$): Good in terms of mass balance but not good in terms of dynamic behavior because at each time interval the velocity might be correct BUT changes for the whole length of the reactor thus does not take into account the residence time.

What are my choices now? At first I thought to calculate the velocity in space like: $$u_{s,i}=\frac{Q_i}{A}=\frac{\sum_{j=1}^{N_c}F_{i,j} \frac{Rg T}{P}}{A}$$ but $F_j=u_s C_j A\Rightarrow u_s=f(u_s)$.
And also the mass balance equation that i showed earlier does not apply. It must be like: $$\frac{\partial (u_sC_j)}{\partial z}=\rho_b \sum_{k=1}^{N_{reac}} \nu_{j,k} R_{k}$$ and discretized $$\frac{((u_sC_j)_i-(u_sC_j)_{i-1})}{\delta z}=\rho_b \sum_{k=1}^{N_{reac}} \nu_{j,k} R_{i,k}$$ or $$u_{s,i}\frac{(C_{i,j}-C_{i-1,j})}{\delta z} +C_{i,j}\frac{(u_{s,i}-u_{s,i-1})}{\delta z} =\rho_b \sum_{k=1}^{N_{reac}} \nu_{j,k} R_{i,k}$$ In the first case i can solve it but what happens when it gets dynamic? I have more variables than equations as $(u_sC_j)_{t,i}$ is one variable and $(C_j)_{t,i}$ is another. Is there a way to overcome this problem?
In the second case I am almost sure that this system of equations does not have one unique solution. Thus my solver keeps giving me wrong (in terms of realistic) solutions (without reaction both velocity and concentration changes in an arithmetically correct way but no reaction should give no change).
Is this statement true or false? Should this last equation along with inlet boundary conditions give me one and only one solution?
Also should I solve along with the mass balance another equation like continuity equation? (I guess along with an equation of state, like $PV=nRT$). If so, I will also face the problem with unique solution because my equations will be: $$u_{s,i}\frac{(C_{i,j}-C_{i-1,j})}{\delta z} +C_{i,j}\frac{(u_{s,i}-u_{s,i-1})}{\delta z} =\rho_b \sum_{k=1}^{N_{reac}} \nu_{j,k} R_{i,k}$$ and $$u_{s,i}\frac{(\rho_i-\rho_{i-1})}{\delta z} +\rho_i\frac{(u_{s,i}-u_{s,i-1})}{\delta z}=0$$ and $$\rho_i=MW_{mean}\frac{P}{Rg T}$$

So, today I tried to solve the dynamic mass balance equation (DMBE) thinking that each quantity inside a differential should be treated like one. Explanation: in the DMBE we have actually 2 unknown variables, concentration $C$ and velocity $u$. One way to solve it is, as mentioned by nluigi bellow, to solve it along with an equation to calculate velocity (will test it tomorrow and get back to you).
The way I tried (and I think I made it) is to think my unknown variable as concentration $C$ and an other quantity like $G=C*u$ (molar flux). In this case the discretized equation do have one unique solution, given initial (in time) and boundary (in space) conditions. But still I need one more equation. The one that I thought is shown below along with the mass balance: $$\frac{\partial C_i}{\partial t}=-\frac{\partial G_i}{\partial z}+\rho_b \sum_{k=1}^{N_{reac}} \nu_{j,k} R_{i,k}$$ $$C_i=\frac{G_i}{\sum_{j=1}^{N_{comp}} G_i}\frac{P}{RgT}$$ Now, the solution of these equations seems to give perfect results both in terms of mass balance as well as in terms of dynamic behavior, but I would like someone to make a comment regarding the validity of this set of equations.

UPDATE
Dear all,
After a lot of tests on solving the last system of equations (last two in this post) I am able to say that the solution does not give good results in terms of dynamic behavior. The truth is that my hopes where really high for this one but they didn't pay off. I may be doing an error during modelling and simulation which may result to this problematic behavior.. but I dont know. I also tested what nluigi mentioned in his answer ($(u \rho A)_{in}=(u \rho A)_{out}$ for each volume) but still the dynamic behavior is not good enough. You see, when I make a step change on one of the inlet component's flowrate, velocity changes at the same timestep all over the reactor's length resulting to an error on outlet molar flowrate calculation $(F_i=C_iuA)$. Here is a figure that might help you understand what I am saying, look closer at $t=5s$ which if the time when the step change is made on one of the inlet flowrates.

The steady state $(t>>time of step change)$ values are correct but the dynamic of the system is not what one would expect. Actually the problem I am experiencing is (in my opinion) attributed to the fact that a term like $\frac{\partial u}{\partial t}$ does NOT appear in the equation where velocity is calculated.
For this reason I decided to solve the mass balance equation $$\frac{\partial C_{i,j}}{\partial t}=-u_i\frac{(C_{i,j}-C_{i-1,j})}{\Delta z}+\rho_b \sum_{k=1}^{N_{reac}} \nu_{j,k} R_{i,k}$$ along with momentum balance: $$\rho _i\frac{\partial u_i}{\partial t}+ \rho _iu_i\frac{u_i - u_{i-1}}{\Delta z}+=0$$ and an equation to calculate density: $$\rho_i=\sum_{j=1}^{N_{comp}}C_{i,j}mw_j$$ Until now, I thing that the results produced this way are looking good but.. let's say that I hope I will not need to update this question again... The next figure shows results of inlet (blue line, at reactor's inlet) and outlet (red line, at reactor's outlet) velocity versus time while tree step changes are imposed at one of the inlet flowrates. As you can see residence time is taken into account and the calculated (this one you can't see but I can assure you that it does) outlet flowrates are smooth and hopefully correct. (Do not mind about the different scale. the last results where for different flowrates)

Now, what is the question (one might ask): I really need a comment on the correctness and validity of this set of equations.

• Your question is very detailed, while that is generally a good thing, I am having difficulties understanding what specifically you are asking. Can you summarize at the end of your post, your key issues? – nluigi Jan 27 '16 at 10:33
• good catch! indeed if the velocity is a function of the axial direction, i.e. $u(z)$ then the convection term needs to be $\partial_z c_i u_z$ rather than $u_z\partial_z c_i$. – nluigi Jan 29 '16 at 14:42
• Can you model it as a series of mixed vessel reactors? – mart Feb 12 '16 at 6:34
• @mart, You mean like a large number of CSTRs in series? I guess I could test this one also. You think that the problem with velocity will not still bother me this way? – ASK22 Feb 12 '16 at 7:56
• Yes, that's what I meant. Velocity will become mass flow from one reactor to the next. Concentration of species at z will become C in a given reactor. – mart Feb 12 '16 at 8:10

A couple of points:

1. Have you considered the stability of your discretization? Specifically are you sure that: $$C=\frac{u\Delta t}{\Delta x}<1$$ This may be difficult in the case $u$ is variable but perhaps you can make an estimate of its average value?
2. Can you calculate the new velocity based on the mass balance? i.e. the total mass at any point in the reactor must be equal: $$\rho A u = \rho_{in} A_{in} u_{in} = \rho_{out} A_{out} u_{out}$$ You could calculate at each point in the reactor: $$u= \frac{\rho_{in} A_{in} u_{in}}{\rho A}$$ where $\rho$ is the total density at that point in the reactor determined from $$\rho = \sum_i^N c_i m_i$$ where $c_i$ is the species concentration and $m_i$ is the molecular mass.
• I understand that this is too long for a comment, but it really isn't an answer. – hazzey Jan 27 '16 at 13:57
• @hazzey, his question is 'what are my choices now'? I give him two options to look into the stability and a method to calculate the velocity. His question is so long, it is difficult to pinpoint exactly what he wants, but i invite you to provide a better answer – nluigi Jan 27 '16 at 15:01
• @nluigi, I am really glad that you spend some of your time to post an answer to my question and, in a way, I think you helped me. My main problem is that (referring to the dynamic form of mass balance, MB) I need a second equation, knowing that the unknown variables are 2*NumberOfzDirectionGridPoint*NumberOfTimeGridPoint. I have one (MB) and I need another one to calculate the velocity. I haven't thought this equation that you just mentioned and i will probably test it tomorrow and get back to you. This actually describes the change of velocity due to density? – ASK22 Jan 27 '16 at 20:18
• @nluigi, I do not have any stability-caused problem. Already looked into this. – ASK22 Jan 27 '16 at 20:21
• @ASK22 - yes locally the density of the gas will change because of the change in number of moles due to the chemical reaction. This way it gives you a method to calculate the velocity based on the density change and makes sure the massflux is everywhere the same in the reactor – nluigi Jan 27 '16 at 20:35

I think u are wrong assuming that velocity is variable along the reactor. The velocity of gas is only the effect of pressure difference according to Darcy-Weissbach putted into Stoke's momentum balance.

For mass balance to be correct, I'd assume that density is changed over the reactor with constant velocity. Namely it depends on the pressure loss, but it is mainly dependent on the ratio of L/d (reactor length to r. diameter), so if it's long enough, the velocity can be changed somewhere.

I haven't considered deeply your issue. Could u attend?