In the dem chemical reaction simulation, an error is reported when the two solid phases particles encounter

[lixiang0426]

Hello everyone！
I encountered a problems in the dem chemical reaction simulation.
I set up separate reactions for the two solid phase particles. But when the particles of solid phase 2 meet the particles of solid phase 1, an error will be reported as shown below：

It seems that there is a problem with the heat conduction. I try to reduce the inlet wind speed (for example, 100m/s), it will not have this error, but the wind speed does not seem to be the root cause.
I want to know where the source of this error is and how to solve it.
I hope you can answer my doubts, I am very grateful again.
My case is attached below.
2_2021-11-03T235134.746095.zip (1.0 MB)

[cgw]

Hello Li -

This is a floating-point exception, that is, a math error.

Here’s the offending line (496) in des_thermo_cond_mod.f

496	         TERM3 = TWO*PI*(ONE-OLAP)*log((ONE-OLAP-sqrt(ONE-Rout**2))/Kn)

I added a write(*,*) statement on the previous line to report the values of OLAP, ROUT, and KN. Right before failing it prints:

OLAP=   3.1569290598721614E-002 
ROUT=   1.0000000000000000      
KN=  -1.0024507434496799E-006

Note that Kn is negative, so the logarithm is being attempted on a negative value:

>>> print( (ONE-OLAP-sqrt(ONE-Rout**2))/Kn)
-966063.1364975298

Looking at L465 in des_thermo_cond_mod.f:

     DOUBLE PRECISION FUNCTION EVAL_H_PFP(RLENS_dim,S,OLAP_dim,RP)
      USE CONSTANT
      USE PARAM1

      IMPLICIT NONE
      ! Note: Function inputs dimensional quantities
      DOUBLE PRECISION, intent(in) :: RLENS_dim, S, OLAP_dim, RP
      ! BELOW VARIABLES ARE NONDIMENSIONALIZED BY RP
      DOUBLE PRECISION :: RLENS, OLAP, KN
      DOUBLE PRECISION :: TERM1,TERM2,TERM3
      DOUBLE PRECISION :: Rout,Rkn
      DOUBLE PRECISION, PARAMETER :: TWO = 2.0D0

      RLENS = RLENS_dim/RP
      KN = S/RP
      OLAP = OLAP_dim/RP

A debugger session with gdb (or more write statements) shows us that S is getting passed in with a negative value:

#1  0x00007fe08357376e in des_conduction::eval_h_pfp (
    rlens_dim=1.0001875441682284, s=-2.0048488274173704e-08, 
    olap_dim=0.00063136922839164208, rp=0.019999474692573838)

|
One frame up the stack:

#2  0x00007fe083574486 in des_thermo_cond::des_conduction (i=394, j=420, 
    center_dist=0.018987945901915475, im=<optimized out>, 
    iijk=<optimized out>)
    at /home/cgw/mfix/model/des/des_thermo_cond_mod.f:184

184            H_2 = EVAL_H_PFP(RLENS_2, S_2, OLAP_2,MAX_RAD)*MAX_RAD*k_gas

and a few lines further up in the same file:

163        RATIO = (OLAP_actual + DES_MIN_COND_DIST) / MAX_RAD
164        S_1 = (MAX_RAD*MAX_RAD*(RATIO-0.5D0*RATIO*RATIO)/&
165            &(CENTER_DIST_CORR-DES_MIN_COND_DIST)) - OLAP_1
166        S_2 =  DES_MIN_COND_DIST - S_1

This is where the negative value is coming in - S_1 is greater than DES_MIN_COND_DIST, so S_2 is negative. I suppose MFiX should have a check for this, we will take a look at the solver code and see if more checks are needed.

But DES_MIN_COND_DIST is an MFiX keyword so you could try setting it to a higher value, to keep S_2 from going negative. However if you set it too high, S_1 could become negative.

It looks like the default value for DES_MIN_COND_DIST is 1E-6. This is not well-documented, see check_data/check_solids_dem.f lines 165

I suggest you try adjusting DES_MIN_COND_DIST and let me know if it helps.

[cgw]

[jeff.dietiker]

I agree with Charles we should have a data check to avoid the fatal FPE. However, this kind of error could be indicative of excessive overlap, and it may crop up at some other stage in the calculation (you could see particles shoot out of the domain). Since you are running with fairly large velocities, the particles may be too soft. You can try increasing the spring constant see if this helps.

There is another issue you need to solve first though: the mesh is too coarse to resolve the inlet tube so you get a bunch of particles that are not injected in the domain.

[lixiang0426]

Thank you very much for your answers, it helped me a lot.

1. I try to increase DES_MIN_COND_DIST, but s_1 will increase along with it, always being a positive value, s_2 is still assigned, and the gap between s_1 and s_2 increases as DES_MIN_COND_DIST increases.like this:
   
   Then I study the expression of s_1:
   
   
   
   Substituting the expression of RATIO into the expression of S_1 and simplifying it can be obtained:
   
   
   
   Since a, b, and c are not simple numbers, the monotonic relationship between S_1 and DES_MIN_COND_DIST is not well studied. The result of the attempt is that within a certain range, s_1 will increase with the increase of DES_MIN_COND_DIST.So adjusting DES_MIN_COND_DIST is not very effective.
2. I got some other hints from Jeff’s answer. The root cause of the problem may be too much overlap, so I followed the suggestion to increase the spring constant. This time it worked.
   
   
   
   Thank you very much again for your help！