FlUENT并行UDF实现：自定义NOx速率模型，炉膛温度及其相对标准差，空燃比等

#include "udf.h"


real average_species_fraction(Domain *d,int id,char *species)
{
    
    real vol_tot=0.,speciesavg=0.;
    #if !RP_HOST
   
    Material *mixture;
    Thread *t;
    int species_index;
    cell_t c;
    mixture=mixture_material(d);
    t=Lookup_Thread(d,id);
  
    species_index=mixture_specie_index(mixture,species);
    
    begin_c_loop(c,t)
    {
        vol_tot +=C_VOLUME(c,t);
        speciesavg +=C_YI(c,t,species_index)*C_VOLUME(c,t);
         //speciesavg +=Pdf_Yi(c,t,species_index)*C_VOLUME(c,t);
             }
    end_c_loop(c,t)
    
    #if RP_NODE
      speciesavg=PRF_GRSUM1(speciesavg);
      vol_tot=PRF_GRSUM1(vol_tot); 
    #endif
    #endif


    node_to_host_real_2(speciesavg,vol_tot);


     
    speciesavg /= vol_tot;
        
    return speciesavg;   
        
}




real averagetemperature(Domain *d,int id)
{
real vol_tot=0.,tavg=0.; 
 
 #if !RP_HOST
 Thread *t;  
 cell_t c;
 t=Lookup_Thread(d,id);
 
     begin_c_loop(c,t)
       {
         vol_tot += C_VOLUME(c,t);
         tavg += C_T(c,t)*C_VOLUME(c,t);
       }    
     end_c_loop(c,t)
  #if RP_NODE
     tavg=PRF_GRSUM1(tavg);
     vol_tot=PRF_GRSUM1(vol_tot); 
  #endif
  #endif 
    
  node_to_host_real_2(tavg,vol_tot); 
   
  return tavg /= vol_tot;


}

real standarderror(Domain *d,int id,real tavg)
{
    
    int totalelement=0;
    real totaltemp,stderror;
    
    #if !RP_HOST
    Thread *t;
    cell_t c;
    t=Lookup_Thread(d,id);
    begin_c_loop(c,t)
       {
        totalelement++;
        totaltemp+=SQR(C_T(c,t)-tavg);
      }
    end_c_loop(c,t)  
    


    #if RP_NODE
    totalelement=PRF_GRSUM1(totalelement);
    totaltemp=PRF_GRSUM1(totaltemp); 
    #endif
    #endif 
    node_to_host_real_1(totaltemp);
    node_to_host_int_1(totalelement);
    stderror=pow(totaltemp/totalelement-1,0.5);
    
    return stderror;
}


real relativestandarderror(real tavg,real stderror)
{
    return stderror*100/tavg;
}


DEFINE_NOX_RATE(user_nox, c, t, Pollut, Pollut_Par, NOx)
{
    #if !RP_HOST
    Pollut->fluct.fwdrate = 0.0;
    Pollut->fluct.revrate = 0.0;


    switch (Pollut_Par->pollut_io_pdf) {
    case IN_PDF:
    /* Source terms other than those from char must be included here*/


       if (POLLUT_EQN(Pollut_Par) == EQ_NO) {


           /* Thermal NOx */
           if (NOx->thermal_nox && NOx->thermal_udf_replace) {
          
          }
       }


    case OUT_PDF:
    /* Char Contributions, that do not go into pdf loop must be 
       included here */


       break;


    default:
       ;
    }
    #endif
}




DEFINE_ADJUST(my_adjust,d)
{
 
 int frontid=7,backid=6;
 int crossid=5,circulationid=4; 
 
 real fronttavg,backtavg;
 real crosstavg,circulationtavg;
 
 real fronterr,backerr;
 real frontrsd,backrsd;
 
 real frontfo2,frontfco2,backfo2,backfco2,frontfco,backfco;
 real k=3.469;
 real frontratio,backratio;
 
 fronttavg=averagetemperature(d,frontid);
 backtavg=averagetemperature(d,backid);
 fronterr=standarderror(d,frontid,fronttavg);
 backerr=standarderror(d,backid,backtavg);
 frontrsd=relativestandarderror(fronttavg,fronterr);
 backrsd=relativestandarderror(backtavg,backerr); 
 
 crosstavg=averagetemperature(d,crossid);
 circulationtavg=averagetemperature(d,circulationid);
 
 
 frontfo2=average_species_fraction(d,frontid,"o2");/*//air excess N2+O2 mass fraction*/
 frontfco2=average_species_fraction(d,frontid,"co2");
 frontfco=average_species_fraction(d,frontid,"co");
 frontratio=(frontfo2+k*(frontfco2))/(k*(frontfco2));
 ////frontratio=(frontfo2-0.5*frontfco+k*(frontfco2+frontfco))/(k*(frontfco2+frontfco));


 backfo2=average_species_fraction(d,backid,"o2");/*//air excess N2+O2 mass fraction*/
 backfco2=average_species_fraction(d,backid,"co2");
 backfco=average_species_fraction(d,backid,"co");
 backratio=(backfo2+k*(backfco2))/(k*(backfco2));
 ////backratio=(backfo2-0.5*backfco+k*(backfco2+backfco))/(k*(backfco2+backfco));
 
  
  #if !RP_NODE     
    Message("\n\nup quirk temperature %.5f,down quirk temperature %.5f.\n",fronttavg,backtavg);//tavg /= vol_tot);
    Message("up quirk relative standard error %.5f.\n",frontrsd);
    Message("down quirk relative standard error %.5f.\n",backrsd);
    Message("up quirk air-fuel ratio %.5f,down quirk air-fuel ratio %.5f.\n",frontratio,backratio);
    Message("cross hole temperature %.5f, circulation ports temperature %.5f.\n\n",crosstavg,circulationtavg);
  #endif
}