Derivation and implementation of c_v (!= c_p!)
 

I have material properties with different heat capacities, c_p and c_v. My material properties are given as tables. Because some quantities in OpenFOAM depend on the heat capacities, gas constant and other derived quantities I need to calculate and implement the heat capacity. OpenFOAM derives most of these values from the gas constant R=c_p-c_v, so if I know the heat capacities I can implement R from which all further quantities are calculated.

I have tried to derive c_v and R like this: https://www.overleaf.com/read/gmhgqhrxwskk . I think this is correct, but can anyone verify this? Is it safe to implement R in a custom species class?

The specie.H class function R() does not contain the p, T parameters, so it's not possible to implement it there. It might be possible to implement a custom R(p,T) function in the next higher class deriving from it, the equation of state. That class defines a compressibility factor and cpMcv (read: cp minus cv, which is basically the same as R). I hope that it's enough to implement these two properties for a flawless implementation.

I think this applies to ideal gases only. If you want to apply real behavior, you need to use the deviate of the density on pressure. The difference between enthalpy and inner energy relates to this entity. Maybe you can generate the slope of that relationship from your table? Another possibility is to use directly in the program, in particular if you have a mixture.

Other thing, maybe useful too: specific heat is used in differential equations, tabled enthalpies are used for integral considerations. This is such since specific heat should depend only slightly on temperature. Call it a linerization of the enthalpy (evaporation is not considered, of course). So, enthalpies are useful in sofware like Aspen, ChemCAD and the like, specific heat should be used in CFD. If you have an enthalpy table, you can use the same concept as described above by using the slope in T of enthalpy to get the specific heat.

Regards,

Daniel

What exactly do you mean with "this"?
In the tabular material properties library I created I use tables that map enthalpy as a function of pressure and temperature and temperature as a function of pressure and enthalpy (basically just the inverse table, because enthalpy is monotonically increasing in temperature). In my case (plasma simulation) there is a huge dependence of c_p on temperature. The single plasma species are treated as a single fluid and so the heat capacity contains peaks at the ionization energies.

I'm currently implementing the tables based on the formulas above, but I still have some numerical problems. I'm close to some reference data though, so I believe the derivation should be correct.

Your question was how to convert c_p into c_v:

since dh/dT = c_p = de/dT + p (dv/dT) (p= const) = c_V + p (dv/dT) (p= const)

You plot v (= 1/rho) on T at constant pressure (your equation of state should relate p, rho and T). Take that plot and deviate it in T, meaning you take the slope (p is eliminated since it is constant).

If you use ideal gases, you can convert more easily using the gas constant. Note that you will get the molar specific heat, so you will need to convert to mass specific heat by dividing through the molar mass.

Regards,

Daniel

I cannot follow the first step in your derivation. h=e+p/rho as you say, however, c_v is defined as de/dT with constant volume, not pressure, so you cannot equation your de/dT (with constant pressure) with c_v. I do this in my derivation by using the total derivative, which I think you haven't considered.

I just noticed I have a wrong index in line three of the c_v derivation, should be with constant V instead of constant p. Here's the fixed version: https://www.overleaf.com/read/whrpjzqvmnzx

I'm pretty sure my derivation is correct, I used my calculated value to calculate the speed of sound according to http://citeseerx.ist.psu.edu/viewdoc...0.1.1.215.7150 and compared it to the curve they show there and it is nearly equal. I still have some wrong results at room temperature, the gas constant is moving away from the real value, so I get a speed of sound of about 310 m/s instead of 340 m/s. I still need to find the cause for that, but the rest is looking fine at the moment.

You are correct, I missed a term, but this does not change the procedure (since I searched to quickly on the internet).

h = u + p v

dh = du + p dv + v dp

If you the deviation in T you get:

dh/dT = c_p = c_V + p dv/dT (p = const) + v dp/dT (v = const)

For a batch vessel, the volume often does not change, meaning p dv/dT is 0. But since you have stuff going in and out, this is not necessarily the case. So, you have to use your equation of state to convert the 2 partial differentials.

By the way: enthalpy and entropy, chemical potential are all linked using the equation of state in one sense or the other. E.g. in PCSAFT, enthalpy is derived from Helmholz energy and some other terms.

Maybe this is too complicated, but as stated, for sure the conversion of c_p into c_V by just adding the gas constant applies to ideal gas only, since in this case p v = n RT.

Note: I used u instead of your e, but I think this is just a difference of anotation.

Regards,

Daniel