The reaction rate of heterogeneous reaction is 0

Hello. I am doing research on CFD-DEM. Among them, the model has involved heterogeneous reactions and has defined the usr _ rates _ des file, but the results of monitor and vtk seem to indicate that the reaction did not occur and the reaction rate was 0. This is my relevant documents
3-6-CFD-DEM.zip (30.2 MB), please take the time to help me check what is wrong, thank you very much.

A primitive, but highly effective, method of debugging is to insert WRITE statements into your code.

      if ((K1*x_g(IJK,CO)/mw_g(CO)-x_g(IJK,CO2)/mw_g(CO2)) <= 0) THEN
         write(*,*) "CASE1"
         DES_RATES(REACTION1) = 0
      else if (x1 < 1d-3) THEN
         write(*,*) "CASE2"
         DES_RATES(REACTION1) = 0
      else
         DES_RATES(REACTION1) = x1**3*s*RO_g(IJK)/w*((a3*(a2+b2+b4)+(a2+b2)*b4)*(K1*x_g(IJK,CO)/mw_g(CO)-x_g(IJK,CO2)/mw_g(CO2))-(a3*(b2+b4)+b2*b4)*(K2*x_g(IJK,CO)/mw_g(CO)-x_g(IJK,CO2)/mw_g(CO2))-(a2*b4)*(K3*x_g(IJK,CO)/mw_g(CO)-x_g(IJK,CO2)/mw_g(CO2)))
         write(*,*) "CASE3", DES_RATES(REACTION1)
         if (DES_RATES(REACTION1) < 0) THEN
            write(*,*) "NEGATIVE"
            DES_RATES(REACTION1) = 0
         ENDIF

I’m seeing a lot of negative values when I run this.

This

x1**3*s*RO_g(IJK)/w*((a3*(a2+b2+b4)+(a2+b2)*b4)*(K1*x_g(IJK,CO)/mw_g(CO)-x_g(IJK,CO2)/mw_g(CO2))-(a3*(b2+b4)+b2*b4)*(K2*x_g(IJK,CO)/mw_g(CO)-x_g(IJK,CO2)/mw_g(CO2))-(a2*b4)*(K3*x_g(IJK,CO)/mw_g(CO)-x_g(IJK,CO2)/mw_g(CO2)))

is a very complex expression, and I’m not going to try to analyze it. I would suggest break this down into some intermediate results, which you can print out and check.