Possible error in the polar contribution

[DiegoTMelfi]

Dear,

I believe there might exist an error in the polar section of the pcsaft_eletrolyte.cpp code.

When calculating the terms dA2_det and dA3_det in lines 308-309 and 763-764) the derivative of the density wrt the reduced density is not being accounted for.

    A2 = -PI*den*A2;
    A3 = -4/3.*PI*PI*den*den*A3;
    dA2_det = -PI*den*dA2_det;
    dA3_det = -4/3.*PI*PI*den*den*dA3_det;

I believe the product rule should be applied there. So these lines would be substituted by the following in the compressibility and fugacity functions:

    dden_det = den/eta;
    dA2_det = -PI*(dden_det*A2 + den*dA2_det); //Here A2 is summation part of Eq. (8) from Gross and Vrabec (2006)
    dA3_det = -4/3.*PI*PI*(2*den*dden_det*A3 + den*den*dA3_det); //A3 is the summation part only           
    A2 = -PI*den*A2;
    A3 = -4/3.*PI*PI*den*den*A3;

Also, it seems to me that when calculating the terms J2, J3, and their derivatives we should use e_ij[i times ncomp+j] instead of e_ij[j times ncomp+j].

For example in lines 269 and 270

                J2 += (adip[l] + bdip[l]*e_ij[j*ncomp+j]/t)*pow(eta, l);        // j*ncomp+j needs to be used for e_ij because it is formatted as a 1D vector
                dJ2_det += (adip[l] + bdip[l]*e_ij[j*ncomp+j]/t)*l*pow(eta, l-1);

But the comment makes me wonder if I am missing something.

Thank you for your attention, Diego

[zmeri]

I finally found some time to take a look at this. Yes, it looks like you are right that the dA2_det derivative was not quite correct. Thank you for pointing that out. Indeed, eta in the density term was not taken into account.

The chain rule actually isn't needed because eta can be factored out of the density term and combined with J2, like so: detJ2_det += (adip[l] + bdip[l]*e_ij[i*ncomp+j]/t)*(l+1)*pow(eta, l); Then the final dA2_det value is calculated by multiplying that sum by den/eta, as such: dA2_det = -PI*den/eta*dA2_det; detJ3_det can be calculated in a similar manner. I update the variable names to reflect that now it is the derivative of eta*J2 (detJ2_det).

And yes, I think you are right that in lines 269 and 270 (and later in the fugcoef function) the index for e_ij should be i*ncomp + j. That makes more sense to me as well. All the tests pass with either j*ncomp + j or i*ncomp + j, which is probably why I never noticed this before. I think problems would occur once there are more than one polar compound in the mixture, but I don't have any data handy for such a mixture and Gross and Vrabec (2006) didn't give any parameters or data for such a system either. I'll update the index to i*ncomp + j.

It will take me a little while to get an update with these changes. I'm also trying to make some changes to the flash algorithm before I commit an update.

Thanks again for mentioning this!

[DiegoTMelfi]

I just sent you an e-mail about this.

Thank you! Diego

[zmeri]

I released an updated version (1.4.1) with corrections to those polar derivatives. Flash improvements will have to come later.