CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   Fluent UDF and Scheme Programming (https://www.cfd-online.com/Forums/fluent-udf/)
-   -   Enthalpy Calculation with UDF (https://www.cfd-online.com/Forums/fluent-udf/245825-enthalpy-calculation-udf.html)

erginbayrak October 28, 2022 03:49
Enthalpy Calculation with UDF
 
Hi,

I have struggled to solve specific heat and enthalpy values of supercritical CO2 in the pipe using the following code (only for 80 bar). Although my code is running in Fluent, the solution does not converge. Also, there is a warning in Fluent as C4715: not all control paths return a value". I want to learn your recommendations about my problem. Regards,

Code:
#include "udf.h"

DEFINE_SPECIFIC_HEAT(supercritical_cp, T, Tref, h, yi)
{
  real cp;

 
    if (300. >= T)
      {
              cp = 2.256999118844760E-02 * pow(T, 4.) - 2.606618818487370E+01 * pow(T, 3.) + 1.129129898467780E+04 * pow(T, 2.) - 2.174235523392760E+06 * T + 1.570278261279480E+08;
                }
    if (306. >= T > 300.)
      {
                cp = 1.147018354719330E+01 * pow(T, 4.) - 1.385519872072390E+04 * pow(T, 3.) + 6.276090044658930E+06 * pow(T, 2.) - 1.263530241490370E+09 * T + 9.539285116131690E+10;
      }
    if (307.8 >= T > 306.)
      {
                cp = -1.611612838918710E+04 * pow(T, 4.) + 1.978085636985280E+07 * pow(T, 3.) - 9.104578256081260E+09 * pow(T, 2.) + 1.862478067464320E+12 * T - 1.428738938943020E+14;       
                }
    if (310. >= T > 307.8)
                {
                cp = -1.576135476848110E+03 * pow(T, 4.) + 1.945866130859760E+06 * pow(T, 3.) - 9.008641337104630E+08 * pow(T, 2.) + 1.853611584329810E+11 * T - 1.430227164258690E+13;
                }
    if (320. >= T > 310.)
                {
                cp = 1.770529718817220E+00 * pow(T, 4.) - 2.242833202341920E+03 * pow(T, 3.) + 1.065432793165650E+06 * pow(T, 2.) - 2.249454519629040E+08 * T + 1.781003821884030E+10;
                }
      if (T > 320)
                {
                cp = 2.086961941877520E-03 * pow(T, 4.) - 2.851496372954040E+00 * pow(T, 3.) + 1.461548345109060E+03 * pow(T, 2.) - 3.330816449664890E+05 * T + 2.848088978782360E+07;
                }
  return cp;
}

#include<stdio.h>
float y(float T)
        {
    if (300. >= T)
       
    return 2.256999118844760E-02 * pow(T, 4.) - 2.606618818487370E+01 * pow(T, 3.) + 1.129129898467780E+04 * pow(T, 2.) - 2.174235523392760E+06 * T + 1.570278261279480E+08;
       
    if (306. >= T > 300.)
       
    return 1.147018354719330E+01 * pow(T, 4.) - 1.385519872072390E+04 * pow(T, 3.) + 6.276090044658930E+06 * pow(T, 2.) - 1.263530241490370E+09 * T + 9.539285116131690E+10;
       

    if (307.8 >= T > 306.)
       
    return -1.611612838918710E+04 * pow(T, 4.) + 1.978085636985280E+07 * pow(T, 3.) - 9.104578256081260E+09 * pow(T, 2.) + 1.862478067464320E+12 * T - 1.428738938943020E+14;
     
    if (310. >= T > 307.8)
       
    return -1.576135476848110E+03 * pow(T, 4.) + 1.945866130859760E+06 * pow(T, 3.) - 9.008641337104630E+08 * pow(T, 2.) + 1.853611584329810E+11 * T - 1.430227164258690E+13;
       

      if (320. >= T > 310.)
       
    return 1.770529718817220E+00 * pow(T, 4.) - 2.242833202341920E+03 * pow(T, 3.) + 1.065432793165650E+06 * pow(T, 2.) - 2.249454519629040E+08 * T + 1.781003821884030E+10;
     
      if (T > 320)
       
    return 2.086961941877520E-03 * pow(T, 4.) - 2.851496372954040E+00 * pow(T, 3.) + 1.461548345109060E+03 * pow(T, 2.) - 3.330816449664890E+05 * T + 2.848088978782360E+07;
        }
// Function to evaluate the value of integral
        float trapezoidal(float Tref, float T, float n, float k, float *h)
        {

        k = (T-Tref)/n;
        *h = y(Tref)+y(T);
        for (int i = 1; i < n; i++)
        *h += 2*y(Tref+i*k);
        return (k/2)*(*h);
        }

        // Driver program to test above function
int main()
{
    // Range of definite integral (initial values)
    float T0 = 295;
    float Tn = 300;
    int n = 2000;
    return 0;
}

AlexanderZ November 2, 2022 23:00
first of all this is your code
Code:
#include "udf.h"

DEFINE_SPECIFIC_HEAT(supercritical_cp, T, Tref, h, yi)
{
  real cp;

 
    if (300. >= T)
      {
              cp = 2.256999118844760E-02 * pow(T, 4.) - 2.606618818487370E+01 * pow(T, 3.) + 1.129129898467780E+04 * pow(T, 2.) - 2.174235523392760E+06 * T + 1.570278261279480E+08;
                }
    if (306. >= T > 300.)
      {
                cp = 1.147018354719330E+01 * pow(T, 4.) - 1.385519872072390E+04 * pow(T, 3.) + 6.276090044658930E+06 * pow(T, 2.) - 1.263530241490370E+09 * T + 9.539285116131690E+10;
      }
    if (307.8 >= T > 306.)
      {
                cp = -1.611612838918710E+04 * pow(T, 4.) + 1.978085636985280E+07 * pow(T, 3.) - 9.104578256081260E+09 * pow(T, 2.) + 1.862478067464320E+12 * T - 1.428738938943020E+14;       
                }
    if (310. >= T > 307.8)
                {
                cp = -1.576135476848110E+03 * pow(T, 4.) + 1.945866130859760E+06 * pow(T, 3.) - 9.008641337104630E+08 * pow(T, 2.) + 1.853611584329810E+11 * T - 1.430227164258690E+13;
                }
    if (320. >= T > 310.)
                {
                cp = 1.770529718817220E+00 * pow(T, 4.) - 2.242833202341920E+03 * pow(T, 3.) + 1.065432793165650E+06 * pow(T, 2.) - 2.249454519629040E+08 * T + 1.781003821884030E+10;
                }
      if (T > 320)
                {
                cp = 2.086961941877520E-03 * pow(T, 4.) - 2.851496372954040E+00 * pow(T, 3.) + 1.461548345109060E+03 * pow(T, 2.) - 3.330816449664890E+05 * T + 2.848088978782360E+07;
                }
  return cp;
}
remove everything below ( starting from #include<stdio.h> )

second point here is your condition sequence:
(300. >= T) will always be true, so only first equation would be executed, also Cp is not specified at temperature below 300
-> go from top bottom

erginbayrak November 3, 2022 03:17
Thank you for your reply. I couldn't understand the logic of the sequence. Is it important for initialization? I have just tried it by obeying your advice, but I take some invalid cp errors (I checked the invalid cp errors at related temperatures by hand, and there is no problem with my equation). My other question is that the calculation of the enthalpy value. According to the manual, we must describe the equation of enthalpy. In my case, it should be integration from Tref to T, shouldn't it?

So, my new code:

Code:
#include "udf.h"

DEFINE_SPECIFIC_HEAT(supercritical_cp, T, Tref, h, yi)
{
  real cp;

   
 if (T > 320.)
                {
                cp = 2.086961941877520E-03 * pow(T, 4.) - 2.851496372954040E+00 * pow(T, 3.) + 1.461548345109060E+03 * pow(T, 2.) - 3.330816449664890E+05 * T + 2.848088978782360E+07;
                *h = (2.086961941877520E-03 / 5.) * (pow(T, 5.) - pow(Tref, 5.)) - (2.851496372954040E+00 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) + (1.461548345109060E+03 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) - (3.330816449664890E+05 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) + 2.848088978782360E+07 * (T-Tref);
                }
 if (320. >= T > 310.)
                {
                cp = 1.770529718817220E+00 * pow(T, 4.) - 2.242833202341920E+03 * pow(T, 3.) + 1.065432793165650E+06 * pow(T, 2.) - 2.249454519629040E+08 * T + 1.781003821884030E+10;
              *h = (1.770529718817220E+00 / 5.) * (pow(T, 5.) - pow(Tref, 5.)) - (2.242833202341920E+03 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) + (1.065432793165650E+06 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) - (2.249454519629040E+08 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) + 1.781003821884030E+10 * (T-Tref);
                }
 if (310. >= T > 307.8)
                {
                cp = -1.576135476848110E+03 * pow(T, 4.) + 1.945866130859760E+06 * pow(T, 3.) - 9.008641337104630E+08 * pow(T, 2.) + 1.853611584329810E+11 * T - 1.430227164258690E+13;
                *h = (-1.576135476848110E+03 / 5.) * (pow(T, 5.) - pow(Tref, 5.)) + (1.945866130859760E+06 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) -  (9.008641337104630E+08 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) + (1.853611584329810E+11 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) - 1.430227164258690E+13 * (T-Tref);
                }

 if (307.8 >= T > 306.)
                {
                cp = -1.611612838918710E+04 * pow(T, 4.) + 1.978085636985280E+07 * pow(T, 3.) - 9.104578256081260E+09 * pow(T, 2.) + 1.862478067464320E+12 * T - 1.428738938943020E+14;       
                *h = (-1.611612838918710E+04 / 5.) * (pow(T, 5.) - pow(Tref, 5.)) + (1.978085636985280E+07 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) - (9.104578256081260E+09 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) + (1.862478067464320E+12 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) - 1.428738938943020E+14 * (T-Tref);
                }

 if (306. >= T > 300.)
                {
                cp = 1.147018354719330E+01 * pow(T, 4.) - 1.385519872072390E+04 * pow(T, 3.) + 6.276090044658930E+06 * pow(T, 2.) - 1.263530241490370E+09 * T + 9.539285116131690E+10;
                *h = (1.147018354719330E+01 / 5.) * (pow(T, 5.) - pow(Tref, 5.)) - (1.385519872072390E+04 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) + (6.276090044658930E+06 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) - (1.263530241490370E+09 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) + 9.539285116131690E+10 * (T-Tref);
                }
 if (300. >= T)
                {
              cp = 2.256999118844760E-02 * pow(T, 4.) - 2.606618818487370E+01 * pow(T, 3.) + 1.129129898467780E+04 * pow(T, 2.) - 2.174235523392760E+06 * T + 1.570278261279480E+08;
              *h = (2.256999118844760E-02 / 5.) * (pow(T, 5.) - pow(Tref, 5.))- (2.606618818487370E+01 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) + (1.129129898467780E+04 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) - (2.174235523392760E+06 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) + 1.570278261279480E+08 * (T-Tref);
                }
  return cp;
}

erginbayrak November 3, 2022 03:19
Quote:
Originally Posted by AlexanderZ (Post 838703)
first of all this is your code
Code:
#include "udf.h"

DEFINE_SPECIFIC_HEAT(supercritical_cp, T, Tref, h, yi)
{
  real cp;

 
    if (300. >= T)
      {
              cp = 2.256999118844760E-02 * pow(T, 4.) - 2.606618818487370E+01 * pow(T, 3.) + 1.129129898467780E+04 * pow(T, 2.) - 2.174235523392760E+06 * T + 1.570278261279480E+08;
                }
    if (306. >= T > 300.)
      {
                cp = 1.147018354719330E+01 * pow(T, 4.) - 1.385519872072390E+04 * pow(T, 3.) + 6.276090044658930E+06 * pow(T, 2.) - 1.263530241490370E+09 * T + 9.539285116131690E+10;
      }
    if (307.8 >= T > 306.)
      {
                cp = -1.611612838918710E+04 * pow(T, 4.) + 1.978085636985280E+07 * pow(T, 3.) - 9.104578256081260E+09 * pow(T, 2.) + 1.862478067464320E+12 * T - 1.428738938943020E+14;       
                }
    if (310. >= T > 307.8)
                {
                cp = -1.576135476848110E+03 * pow(T, 4.) + 1.945866130859760E+06 * pow(T, 3.) - 9.008641337104630E+08 * pow(T, 2.) + 1.853611584329810E+11 * T - 1.430227164258690E+13;
                }
    if (320. >= T > 310.)
                {
                cp = 1.770529718817220E+00 * pow(T, 4.) - 2.242833202341920E+03 * pow(T, 3.) + 1.065432793165650E+06 * pow(T, 2.) - 2.249454519629040E+08 * T + 1.781003821884030E+10;
                }
      if (T > 320)
                {
                cp = 2.086961941877520E-03 * pow(T, 4.) - 2.851496372954040E+00 * pow(T, 3.) + 1.461548345109060E+03 * pow(T, 2.) - 3.330816449664890E+05 * T + 2.848088978782360E+07;
                }
  return cp;
}
remove everything below ( starting from #include<stdio.h> )

second point here is your condition sequence:
(300. >= T) will always be true, so only first equation would be executed, also Cp is not specified at temperature below 300
-> go from top bottom
Thank you for your reply. I couldn't understand the logic of the sequence. Is it important for initialization? I have just tried it by obeying your advice, but I take some invalid cp errors (I checked the invalid cp errors at related temperatures by hand, and there is no problem with my equation). My other question is that the calculation of the enthalpy value. According to the manual, we must describe the equation of enthalpy. In my case, it should be integration from Tref to T, shouldn't it?

So, my new code:

Code:
#include "udf.h"

DEFINE_SPECIFIC_HEAT(supercritical_cp, T, Tref, h, yi)
{
  real cp;

   
 if (T > 320.)
                {
                cp = 2.086961941877520E-03 * pow(T, 4.) - 2.851496372954040E+00 * pow(T, 3.) + 1.461548345109060E+03 * pow(T, 2.) - 3.330816449664890E+05 * T + 2.848088978782360E+07;
                *h = (2.086961941877520E-03 / 5.) * (pow(T, 5.) - pow(Tref, 5.)) - (2.851496372954040E+00 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) + (1.461548345109060E+03 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) - (3.330816449664890E+05 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) + 2.848088978782360E+07 * (T-Tref);
                }
 if (320. >= T > 310.)
                {
                cp = 1.770529718817220E+00 * pow(T, 4.) - 2.242833202341920E+03 * pow(T, 3.) + 1.065432793165650E+06 * pow(T, 2.) - 2.249454519629040E+08 * T + 1.781003821884030E+10;
              *h = (1.770529718817220E+00 / 5.) * (pow(T, 5.) - pow(Tref, 5.)) - (2.242833202341920E+03 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) + (1.065432793165650E+06 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) - (2.249454519629040E+08 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) + 1.781003821884030E+10 * (T-Tref);
                }
 if (310. >= T > 307.8)
                {
                cp = -1.576135476848110E+03 * pow(T, 4.) + 1.945866130859760E+06 * pow(T, 3.) - 9.008641337104630E+08 * pow(T, 2.) + 1.853611584329810E+11 * T - 1.430227164258690E+13;
                *h = (-1.576135476848110E+03 / 5.) * (pow(T, 5.) - pow(Tref, 5.)) + (1.945866130859760E+06 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) -  (9.008641337104630E+08 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) + (1.853611584329810E+11 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) - 1.430227164258690E+13 * (T-Tref);
                }

 if (307.8 >= T > 306.)
                {
                cp = -1.611612838918710E+04 * pow(T, 4.) + 1.978085636985280E+07 * pow(T, 3.) - 9.104578256081260E+09 * pow(T, 2.) + 1.862478067464320E+12 * T - 1.428738938943020E+14;       
                *h = (-1.611612838918710E+04 / 5.) * (pow(T, 5.) - pow(Tref, 5.)) + (1.978085636985280E+07 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) - (9.104578256081260E+09 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) + (1.862478067464320E+12 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) - 1.428738938943020E+14 * (T-Tref);
                }

 if (306. >= T > 300.)
                {
                cp = 1.147018354719330E+01 * pow(T, 4.) - 1.385519872072390E+04 * pow(T, 3.) + 6.276090044658930E+06 * pow(T, 2.) - 1.263530241490370E+09 * T + 9.539285116131690E+10;
                *h = (1.147018354719330E+01 / 5.) * (pow(T, 5.) - pow(Tref, 5.)) - (1.385519872072390E+04 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) + (6.276090044658930E+06 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) - (1.263530241490370E+09 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) + 9.539285116131690E+10 * (T-Tref);
                }
 if (300. >= T)
                {
              cp = 2.256999118844760E-02 * pow(T, 4.) - 2.606618818487370E+01 * pow(T, 3.) + 1.129129898467780E+04 * pow(T, 2.) - 2.174235523392760E+06 * T + 1.570278261279480E+08;
              *h = (2.256999118844760E-02 / 5.) * (pow(T, 5.) - pow(Tref, 5.))- (2.606618818487370E+01 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) + (1.129129898467780E+04 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) - (2.174235523392760E+06 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) + 1.570278261279480E+08 * (T-Tref);
                }
  return cp;
}

erginbayrak November 4, 2022 03:56
My following code is compiling very well, but I have some convergence problems. There is divergence detected and floating point exception warnings. If I decrease the URF of the energy equation, iteration starts, but I take some errors after a certain iteration again. This situation didn't change although used a finer mesh structure. I am continuing to study on it. Any comments or suggestions are welcome. Regards. EB

My revised code is as follows:

Code:

#include "udf.h"

DEFINE_SPECIFIC_HEAT(supercritical_cp, T, Tref, h, yi)
{
  real cp;

   
 if (T > 320.)
  {
               
                cp = 2.086961941877520E-03 * pow(T, 4.) - 2.851496372954040E+00 * pow(T, 3.) + 1.461548345109060E+03 * pow(T, 2.) - 3.330816449664890E+05 * T + 2.848088978782360E+07;
                *h = (2.086961941877520E-03 / 5.) * (pow(T, 5.) - pow(Tref, 5.)) - (2.851496372954040E+00 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) + (1.461548345109060E+03 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) - (3.330816449664890E+05 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) + 2.848088978782360E+07 * (T-Tref);
               
                }
 else if (320. >= T > 310.)
  {
               
                cp = 1.770529718817220E+00 * pow(T, 4.) - 2.242833202341920E+03 * pow(T, 3.) + 1.065432793165650E+06 * pow(T, 2.) - 2.249454519629040E+08 * T + 1.781003821884030E+10;
              *h = (1.770529718817220E+00 / 5.) * (pow(T, 5.) - pow(Tref, 5.)) - (2.242833202341920E+03 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) + (1.065432793165650E+06 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) - (2.249454519629040E+08 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) + 1.781003821884030E+10 * (T-Tref);
               
                }
 else if (310. >= T > 307.8)
  {
               
                cp = -1.576135476848110E+03 * pow(T, 4.) + 1.945866130859760E+06 * pow(T, 3.) - 9.008641337104630E+08 * pow(T, 2.) + 1.853611584329810E+11 * T - 1.430227164258690E+13;
                *h = (-1.576135476848110E+03 / 5.) * (pow(T, 5.) - pow(Tref, 5.)) + (1.945866130859760E+06 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) -  (9.008641337104630E+08 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) + (1.853611584329810E+11 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) - 1.430227164258690E+13 * (T-Tref);
               
}
 else if (307.8 >= T > 306.)
  {
               
                cp = -1.611612838918710E+04 * pow(T, 4.) + 1.978085636985280E+07 * pow(T, 3.) - 9.104578256081260E+09 * pow(T, 2.) + 1.862478067464320E+12 * T - 1.428738938943020E+14;       
                *h = (-1.611612838918710E+04 / 5.) * (pow(T, 5.) - pow(Tref, 5.)) + (1.978085636985280E+07 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) - (9.104578256081260E+09 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) + (1.862478067464320E+12 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) - 1.428738938943020E+14 * (T-Tref);
               
}
 else if (306. >= T > 300.)
  {
               
                cp = 1.147018354719330E+01 * pow(T, 4.) - 1.385519872072390E+04 * pow(T, 3.) + 6.276090044658930E+06 * pow(T, 2.) - 1.263530241490370E+09 * T + 9.539285116131690E+10;
                *h = (1.147018354719330E+01 / 5.) * (pow(T, 5.) - pow(Tref, 5.)) - (1.385519872072390E+04 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) + (6.276090044658930E+06 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) - (1.263530241490370E+09 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) + 9.539285116131690E+10 * (T-Tref);
             

            } 
 else
  {
             
              cp = 2.256999118844760E-02 * pow(T, 4.) - 2.606618818487370E+01 * pow(T, 3.) + 1.129129898467780E+04 * pow(T, 2.) - 2.174235523392760E+06 * T + 1.570278261279480E+08;
              *h = (2.256999118844760E-02 / 5.) * (pow(T, 5.) - pow(Tref, 5.))- (2.606618818487370E+01 / 4.) * (pow(T, 4.) - pow(Tref, 4.)) + (1.129129898467780E+04 / 3.) * (pow(T, 3.) - pow(Tref, 3.)) - (2.174235523392760E+06 / 2.) * (pow(T, 2.) - pow(Tref, 2.)) + 1.570278261279480E+08 * (T-Tref);
            } 
  return cp;
}


All times are GMT -4. The time now is 07:01.