I am attempting to generate breakthrough curves under different feeding concentrations and loading flow rates. I modified the codes from the tutorial (05-simple chromatographic processes, 01-chromatography_demonstration_solution, Example 2). This example is modified to one component and one section (only loading), with the parameters fixed. The breakthrough results are pretty weird and the curves are attached here. The codes are also attached.
The x-axis ‘time’ is transformed to feeding volume (time* flowrate/column volume) and y-axis ‘outlet concentration’ is transformed to breakthrough (outlet conc./feeding conc.*100%). The theoretical breakthrough curves should be the same shape as the following figure (Shi, C, et al, 2020, JCA, 1619, 460936).
The breakthrough when c_feed is 0.01 or 0.05 exceeds 100%, while is less than 100% when c_feed is 0.1. What’s more, the order of curves under different concentrations is also not correct, as breakthrough happens first under higher concentration. The curves do not change when flowrate changes.
Thanks for the prompt reply. You are right and I did overwrite the c_feed. I got the same breakthrough curves as you under different feeding concentrations. Sorry for missing that. However, the breakthrough curves under different loading flowrates are still overlapping, which is obviously not correct. Any advice/suggestions on resolving this would be much appreciated!
You’re normalizing time with respect to column volume. Hence, it is expected that the breakthrough comes at a similar point on your plot. If you plot the breakthrough over time, you can see the effect of different flow rates.
Thanks for your explanation. First sorry about that. Yes, the difference can be seen if we directly plot the outlet concentration (breakthrough) versus the time. But if we take look at the breakthrough curves under different flowrates from the literature, it can be found that the breakthrough curve at a higher flowrate becomes more broadening. In this case, even if time is normalized to column volume, the difference in the shape of breakthrough curves still can be observed. But in my simulation, the shape of breakthrough curves under different flowrates seems to be the same. I am very confused about this. Many thanks for your help!
In reality, the dispersion depends on the interstitial velocity / volumetric flow rate.
Does your model take this into account (e.g., via some correlation between axial dispersion coefficient and flow rate)?
Have you managed to replicate the breakthrough curves in the papar (Shi, C, et al, 2020, JCA, 1619, 460936) using CADET? Would you please share the code if possible? I have tried a lot of times and still have no luck to generate the same breakthrough curves.
Thanks for checking. My idea is that the Langmuir model form used in this paper is a little different from the one provided in CADET, so I guess we need to perform parameter estimation by ourselves. In addition, the breakthrough curves are from the twin-column process data, instead of the batch column data. I only modified the parameters from this paper and tried to simulate the batch breakthrough curve. I didn’t build the twin column model and perform parameter estimation by myself. Technically speaking, I am not sure where the problem is as I didn’t fully replicate this work.
Sorry for the late response, here I have attached the code and let’s find out how to use CADET to replicate the paper’s result.
First things first, here’s the link to the paper that we are discussing:
@Chaoying Hi Chaoying, thank you for your opinions. I think the problem is not related with the twin-column system as the author indicates that the Fig.5 is for the breakthrough curves of the second column. Fig4, which is the one you showed for this topic is indeed the breakthrough curves for the first column. Therefore, CADET should be able to replicate the curves by treating the first column as the batch column.
Some thoughts I have in my mind and maybe useful for determining the deviations:
Unit conversions. But I don’t think it’s a big issue as long as the units are consistent for all the parameters, for example I use g/L for c0, qmax, kd. Should I worry about converting the units?
The author uses eqn.10 to approximate De using the values of Dp, Ds, ep, qmax, kd, c0. I think this approximation is reasonable and should not cause large deviations from the CADET results.
For the GRM model, the governing equation from paper eqn(7) is slightly different from the CADET eqn(4). Maybe the value of Dp and Ds should be adjusted for the CADET by Dp_cadet = Dp_paper/ep and Ds_cadet = Ds_paper/(1-ep) ?
(By doing so, the CADET result is much closer to the paper’s result, but still not the same, check: Replicate attempt-V2.ipynb (391.1 KB))
For the Langmuir binding, the governing equation from paper eqn(3) is under the assumption of quasi-stationary (dq/dt=0) binding mode. By converting the Langmuir model from Cadet, I get the result as q = (qmax*c)/(kd/ka + c). The author uses kd instead of kd/ka, so I assume value for ka as 1 in the simulation. Maybe we can look into that as well.
In the Chinese version the author uses 42.5e-6m for rp while in this paper author uses 85e-6m for rp. I am not sure which one is used for the corresponding graph, so try both values to see if it helps.
Here’s the code for my simulation and validation, I am really looking forward to find out the reason why I failed to replicate the result. You can use the paper result file to validate the breakthrough curve. As you can see, for each breakthrough curve, the CADET return approximately the same time to start breakthrough, but the time of reaching 100% is quite different.
But I feel it’s not possible to fix by simply rescaling the Diffusion coefficients since the entire adsorption term is affected which is by nature nonlinear.
To proceed further, I’d recommend fitting CADET to the experimental data. It might also be worth contacting the author and ask if we’re missing something here.
Shi et al. (2020) use different units for the pore phase concentration than CADET. Let’s denote them by \hat{q} (Shi et al., 2020) and q (CADET). Both are concentrations but \hat{q} refers to the whole bead volume while q only to the pore volume within the bead. Hence, \hat{q}=(1-\varepsilon_p)q. Consequently, also \hat{q}_{max}=(1-\varepsilon_p)q_{max}. With this parameter transform the units should be the same. However, the results still don’t match. There is a slight flaw (or typo) in the Danckwerts boundary condition (eq. 6a) of Shi et al. (2020). The \varepsilon is not required here. However, this should not make a big difference. In rapid adsprption equilibrium, the relation in eq. 8 of Shi et al. (2020) should yield exact results within numerical accuracy. I have contacted the corresponding author. Hopefully, he can contribute to rectifying the observed differences.
The first author of the paper discussed, Ce Shi, couldn’t reply to this post because of some technical issues, so he reached out to me to help post his responses. The original note from him is provided below.
Dear Yunhao and Chaoying,
This is from Ce Shi, the author of the above paper discussed. It is a great pleasure to see that the paper is getting noticed in the continuous chromatography modeling community and that you are attempting to replicate my results using CADET. My supervisor Dong-Qiang Lin saw the email from CADET team and asked me to address your questions and perhaps to provide some suggestions.
Answers to the above discussions:
Governing equation difference: Eqn 7 in the paper is not consistent with the equation used in CADET as our definition towards porosity is different, just like what Eric Von Lieres has identified. So you may focus on how to transfer the parameter based on such a difference in definition, and I think this would be the primary cause of issues.
Boundary condition difference: The difference of Danckwert boundary condition is present because u is in fact a superficial velocity, and you can refer to https://doi.org/10.1002/aic.690470224 for detail. This may not be critical just as Eric suggested, but it would be good to check if the right velocity value is used.
The radius information is a typo, and you can look up the resin information on the internet.
Some suggestions
If you can provide your results with jpg, png version instead of CADET code, maybe I can understand your questions better.
A paper has successfully reproduced my work using CADET, and they obtained very good fitting results. You can refer to https://doi.org/10.1002/jctb.6922 for more information.
If you can understand Chinese, perhaps reading my PhD thesis would also be helpful for getting a better perspective. If you wish, you can contact me through WeChat: 13777887731, and we can figure the problem out together.
Best Regards,
Ce Shi
The mentioned paper of Deulgaonkar et al. (2021) presents CADET simulations. I assume the authors have converted the parameters as discussed above. Maybe they can share their configuration files.
We have used the cadet equation framework for simulating the breakthrough data in Deulgaonkar et al. (2021). Shi’s paper has a slightly different equation framework than a cadet’s. I have not reproduced Shi’s result with the cadet.
Now, I have tried to reproduce the results of Shi’s paper to answer on a forum at 2 min RT and 0.5, 1 mg/mL feed titer.
The attached file contains a unit conversion excel file, an m file for code, and obtained results in jpg format at 0.5 and 1 mg/mL.
If we use superficial velocity instead of interstitial velocity then we can get same lifting point in the BT curve as in Shi’s paper but I think the broadening of the curve is different due to variations in the equations. I have assumed the 150 kDa molecular weight of the molecule for unit conversion.
OK, so it seems there are two major differences between the equations used by Shi et al. (2020) and CADET.
When \hat{u} is the superficial velovity, it needs to be divided by \varepsilon to determine the interstitial velocity u. So if equation 6a is correct there must be a typo in eq. 1 of Shi et al. (2020).
When \hat{q} refers to the whole particle volume, q_{max} also needs to be adjusted as described above.
