search

TorrensJoaquin / GERG2008

Code written by Eric Lemmon. Posted on the NIST website. But now in JavaScript, the best EOS for natural gases straight from the browser.
MIT License
11 stars 1 forks source link

Doubt with residual parameters calculation #1

Open fedebenelli opened 3 years ago

fedebenelli commented 3 years ago

I have a doubt on how you calculate the residual Helmholtz energy, specifically the first derivative with the reduced density.

According to the GERG-2008 paper the polinomial term should be: image

In your code, around line 555:

for(let k = 1; k <= kpol[i]; k++){
          ndt = x[i] * delp[doik[i][k]] * taup[i][k];
          ndtd = ndt * doik[i][k];
          ar[0][1] = ar[0][1] + ndtd;

For what I understood, taup[i][k] equals to n_{oi,k} * tau ^ t_{oi,k}. So each term of the derivative is calculated like:

delta^d_{oi,j} * n_{oi,k} * tau ^ t_{oi,k} * d_{oi,j} When the term of the paper is: delta^(d_{oi,j}-1) * n_{oi,k} * tau ^ t_{oi,k} * d_{oi,j}

I'm sure your code work, since I've tested a bit of it and reached good results, could you tell me what I'm reading wrongly? (I'm trying to make my own implementation on Python and Fortran and not getting the expected results)

Thank you for your time!

TorrensJoaquin commented 3 years ago

It's possible that you have the same mistake i had. Let me know if i'm right. Let me use uppercase letters for DELTA = d/dr and lowercase for the partial derivative operator = delta. When you first look at the meaning of the output in the comments. It seems like a[0][1] is delta aij/ delta DELTA. But if you look at it again you will notice that is DELTA * delta aij/ delta DELTA. That cancel out the -1 in the exponent.

//Outputs;
// ar(0,0) - Residual Helmholtz energy (dimensionless, =a/RT)
// ar(0,1) -     delta*partial  (ar)/partial(delta)
// ar(0,2) -   delta^2*partial^2(ar)/partial(delta)^2
// ar(0,3) -   delta^3*partial^3(ar)/partial(delta)^3
// ar(1,0) -       tau*partial  (ar)/partial(tau)
// ar(1,1) - tau*delta*partial^2(ar)/partial(tau)/partial(delta)
// ar(2,0) -     tau^2*partial^2(ar)/partial(tau)^2

PS: I'm aware of two implementation in FORTRAN https://github.com/n-centrix/ISO20765-Part2 https://github.com/usnistgov/AGA8/tree/master/AGA8CODE/FORTRAN (From NIST - Ian Bell and Eric Lemmon) I'm not aware of any implementation on python. Looking forward to see it.

fedebenelli commented 3 years ago

Thanks for the quick response! Ohh, now I understand, thanks for the clarification! Is there any reason you do it this way besides avoiding using the "-1"?

TorrensJoaquin commented 3 years ago

When you run CalculateDensity. This routine is inside the loop of the Newton Raphson method. I understand that Eric tried to avoid expensive computations inside the loop.