I am doing biomass pyrolysis study with Debiagi scheme in MFIX-Exa-PIC,
the case looks like below:
Please check if anything went wrong.
@YupengXu just to clarify there are two sepearate issues here:
is that correct?
Yes, you are right! Since they are both about the biomass pyrolysis work, so I put them together.
unit-test-PIC.tar (100 KB)
@wfullmer Here is a unit test for the kinetic model, and the compare.py just sums up the total mass of gas and solids from monitors. Please try it out and let me know if you have any questions.
@YupengXu I ran your inputs with your reactions. When I use your python script to post, I get a similar plot:
When I add ep_g to your density monitor (and see it’s constant at roughly 0.931), and calculate total mass as lagrangian mass + ep_g * ro_g * V_tot(=1 m^3), I see the total mass decrease. I see this in grace w/ ep_g as a function of time:
and in excel with ep_g $\approx$ 0.931:
Not that losing mass is any better than gaining mass, but I need to figure out if we’re actually looking at the same thing first. Can you tell me how you’re calculating total mass? I tried to look at your python script but I just 
hate python so much.
I have tow monitors: one for the VolumeIntegral of ro_g to get gas mass, and another is the summation of Lagrangian mass, then the python script just add these 2 together.
I did a quick try, as I multiply the ep_g, I got the same figure as yours.
OK. So issue 1, you were not including volume fraction in your total mass calculation. Issue 2 is that the integral of ep_g * the integral of rho_g is not the same as the integral of their product. To do that, change your eulerian density monitor to type Eulerian::VolumeIntegral::MassWeightedIntegral and set the variable to ones (because all you want is ep_g * rho_g, see docs here). That gets you closer. When you sum that with your lagrangian mass, you should get something like this:
As you can see, it depends somewhat on your precision. If you bump setprecision down to 4, you see a lot of noise that’s just truncation; if you bump it up to 8, the truncation error is gone but there’s still a little behavior there that’s not hard zero. Two things come to mind: 1) summing the eulerian and lagrangian masses separately like this is still an approximation, could be some error in that or 2) probably more likely, we’re down in the 6th sig fig here and I see you have hacked into the code m_mass_balance_tolerance(1.e-5) (instead of the default 1.e-12), so you’re probably just running into your own mass imbalance issues at this point.
Which, BTW, @jmusser points out you can do from the inputs with chemistry.mass_balance_tolerance = 1e-5
So based on this unit test, the kinetic part and the stiff solver is OK.
How about the implementation to calculate the solids outlet flow, did it print the correct particle information?
amr.max_grid_size_x = 6
amr.max_grid_size_y = 6
amr.max_grid_size_z = 4
So many idle nodes available, and the simulation can finish 120s within a day (very helpful for me to check different settings).
Ok, you’re right. It’s pretty slow on 43 cores. I don’t know what Oil, Gas and Char, those are not listed in your species list. I adjusted the total fluid mass monitor as above, added a monitor for lagrangian mass flow at the outlet and moved the eulerian monitor for outlet mass flow down to match, just a little bit at z = 0.42. (It’s probably totally fine to monitor at the exit. Just that old CFD habits die hard. I never want to monitor adjacent to BCs.) Looking at the change in total mass inside the system vs. what I expect to see by summing and integrating the inlet/outlet mass flow rates, I do see a difference.
If we really need to dig into this, I should probably start with a much simpler case and add in complexity to see when it breaks. But I’m not entirely sure if this 
is the issue you are asking about. Actually, I’m not sure exactly what you’re asking at all.
To underscore @wfullmer’s comments – one of the chemical reactions is unbalanced.
By default, MFIX-Exa throws an error when the mass of products does not equal the mass of reactants. You got around this by commenting out the check in chemical reactions initialize routine, but similarly, adding the following setting to your inputs also skirts the issue.
chemistry.mass_balance_tolerance = 1e-5
Running the simplified 16^3 grid closed box case shows a (quasi) steady-state mass error of O(1e-4) kg. If the offending reaction formula is modified so that it balances, the error drops to O(1e-11) which coincides with the sig-figs being written into the monitor files.
This is to say, there are checks within the code to provide guard rails, and by removing them, you risk introducing error into your simulation.
