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}$$
New Comments referring to progress made on this problem
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.