Bauer & Schlunder thermal conductivity model

MFiX
#1

Dear MFiX forum,
I have been checking the equations of Bauer & Schlunder model following the theoretical review: Theoretical Review of the MFIX Fluid and Two-Fluid Models (osti.gov)

However, I am unable to replicate the results of Figure 3-4 using the same values of conductivities stated in the review. I have in fact copy-pasted (and adapted) the code in MFiX and I am still not replicating the values in the graph.

I am doing something wrong or there is a typo in the equations? Also, I was wondering about the validity of that model for high thermal conductivity particles (values around 200-300 W/(m·K).

Thank you for your help!

#2

Can you please attach the input files for your simulation?

#3

Hello,

there is no associated simulation. I am just trying to replicate the curves in Figure 3-4 by manually implementing the equations in the Theory Guide in Matlab. However, the bed conductivities I am obtaining are not coherent with the curves in Figure 3-4. I have tripled-checked everything and I am still getting some different curves, which is why I wonder if I am doing something wrong or if there is something weird with the equations. I can attach the Matlab files if needed.

Thank you!

#4

Please attach your Matlab files, but I cannot guarantee that we have the resources to follow up.

#5

Hello Charles,
please see attached the Matlab file. I checked it again, it is probably something silly, but I don’t see where the difference can come from.
Thanks!
Bauer_cond_MFiX.zip (690 Bytes)

#7

I was able to run your MATLAB code and for reference this is what the plot looks like:

Figure_1
Figure_11269×1019 18.7 KB

and here’s figure 3-4 from the cited paper:

plot
plot949×602 102 KB

I wrote a small Python script to do the calculations from the paper, and here’s my result:

plot2

This looks a lot closer to your result than figure 3-4 in the cited paper.

Here’s the script that produced above plot:

#!/usr/bin/env python
#  Attempt to reproduce figure 3-4 from reference 1
# References:
# 1: https://www.osti.gov/biblio/1604993
# 2: https://mfix.netl.doe.gov/forum/t/bauer-schlunder-thermal-conductivity-model/4069

from math import *
ln = log

# table 2-3
𝑇_𝑟𝑒𝑓 = 300
𝜅_𝑟𝑒𝑓 = 0.0252

# eq 3.87
𝜌_𝑘 = sqrt(3.5E-4)

# eq 3.86
𝜑 = 23 * 𝜌_𝑘**2 / (1 + 22*𝜌_𝑘**(4/3))

for 𝑇_𝑔 in range(300, 2101, 300):
    output = open('T%d.dat' % 𝑇_𝑔, 'w')
    𝜀_𝑔 = 0.15
    while 𝜀_𝑔 <= 1.0:
        # eq 3-88
        𝐵 = 1.25 * pow((1-𝜀_𝑔)/𝜀_𝑔, 10/9)

        # eq 2-28
        𝜅_𝑔 = 𝜅_𝑟𝑒𝑓 * sqrt(𝑇_𝑔 / 𝑇_𝑟𝑒𝑓)

        # eq 3-89
        𝑅 = 1 / 𝜅_𝑔

        # eq 3-86
        𝛤 = (-2/(1-𝐵/𝑅)) * ( 𝐵*(𝑅-1) / (𝑅*(1-(𝐵/𝑅))**2) * ln(𝐵/𝑅) +
                             (𝐵-1) / (1-(𝐵/𝑅)) + (𝐵+1)/2 )

        # eq 3-85
        𝜅_𝑏 = 𝜅_𝑔 * sqrt(1 - 𝜀_𝑔) * (𝜑*𝑅 + (1-𝜑)*𝛤)

        output.write("%g %g\n" % (𝜀_𝑔, 𝜅_𝑏))
        𝜀_𝑔 += 0.001
    output.close()

and gnuplot was used for producing the plot

bs.py (1.1 KB)

More investigation is called for.

#8

Hello Charles,

thank you for your answer and for the code. I will in parallel also dig a bit on why this difference may come from. I will get in touch in case I find some useful information.

Kind regards,
Eduardo

#9

The most complex equation is 3-86 for 𝛤 and I suspected this was most likely to contain a typo. However lines 120-121 of calc_k_s.f read

               BoR  = BB/R_km
               L_rm = -(2.d0/(ONE-BoR)) * &
                      ( ((R_km-ONE)/(ONE-BoR)**2)*BoR*LOGo(BoR) + &
                        (BB-ONE)/(ONE-BoR) + (BB+ONE)/2.d0 )

L_rm is the variable referred to as 𝛤 in the paper.

I was able to copy this into the Python and run it, with only trivial modifications (2.d0 → 2)

        𝛤 = (-2/(1-𝐵/𝑅)) * ( 𝐵*(𝑅-1)* ln(𝐵/𝑅) / (𝑅*(1-(𝐵/𝑅))**2)  +
                             (𝐵-1) / (1-(𝐵/𝑅)) + (𝐵+1)/2 )

        BoR = B/R
        ONE = 1
        BB = B
        R_km = R
        LOG = log
        L_rm = ( -(2/(ONE-BoR)) *
                 ( ((R_km-ONE)/(ONE-BoR)**2)*BoR*LOG(BoR) +
                   (BB-ONE)/(ONE-BoR) + (BB+ONE)/2. ))


        if (abs(𝛤-L_rm)) > 1e-12:
            print(𝛤, L_rm, 𝛤-L_rm)

I’m seeing differences on the order of 1e-15, due to different order-of-operations. But this shows that our 𝛤 at least is correct.
This

#10

My guess is what is plotted on Fig 3.4 is actually Ks (Eq. 3.84) not Kb (Eq. 3.85). Can you please try that and see if it makes more sense?

#11

No, I’m afraid that’s not it.

Eq 3.48 reads 𝜅_𝑠 = 𝜅_𝑏 / (1 - 𝜀_𝑔)

So as 𝜀_𝑔 goes to 1, this increases rather than decreasing … totally different shape, as in these plots:
plot3

#13

I communicated with one of the authors of the paper and we believe the published Fig 3-4 is in error.

The published equations are correct and the above plots above produced by @edcanop and myself are also correct. Fig 3-4 is incorrect due to an error in the computation of gamma,

img2

instead of

img1