Derivatives applied to thermodynamic properties
When dealing with engineering problems we often face (partial) derivatives. Since most functions can not be differenciated analytically, we need numerical methods. This blog post deals with the numerical differentiation and its application to thermodynamic properties.
Theoretical Basics and Numerical error
First, let us start with a bit of theory. If you are only interested in the application, just continue reading. In case you wanna have a look intothe theory, hit the Read More button.Read More...
Let us have a look at the Taylor series to understand this context. A Taylor series represents a function by an infinite sum of terms. So, lets expand the Taylor row for a function at :
Here, is in the interval (x, x+h). Solving for leads to
Here stands for the remainder of the higher order terms. Hence, the first error term is proportional to . This means, the forward difference method is of first order. Additionally, another first order method is the backward difference method. It is not explained here in detail.
Let us consider a second Taylor series:
By subtracting the first Taylor series from the second we get:
Afterwards, let us solve that for
Since the first error term is proportional to , the central difference method is of second order.
In an exact differentiation, would approach 0. While differentiating numerically, is finite. The truncation error is the difference between a derivative and its finite difference representation.
Furthermore, we need to consider the roundoff error. The roundoff error represents the difference between an exact number and its numerical representation. For example, a variable with the data type double can hold approximately 16 significant figures. Simply put, it arises by chopping off numbers. The smaller the larger the impact of the roundoff error.
Finite difference representation of derivatives
In order to calculate the derivative of a function numerically, we can use its finite difference representation. The forward difference approximation is of first order and can be defined as:
Furthermore, the backward difference approximation is also of first order and can be defined as:
Additionally, the central difference approximation is of second order and can be defined as:
Specific heat capacity
The specific heat capacity at constant pressure is defined as
This means, the enthalpy is derivated with respect to the temperature under constant pressure .
CoolProp is a fantastic open source library to receive thermodynamic properties. I use it extensively and there is a lot of functionality. Since it is even possible to calculate derivatives, we use it for benchmarking reasons.
Applying central difference to specific heat
We can calculate the specific heat capacity by applying the central difference approximation as follows:
Hence, we need to be able to calculate a temperature from the enthalpy and pressure (T = f(p, h)). We use a CoolProp function in order to save us some time.
Furthermore, it is super important to choose an adequate by which the enthalpy is pertubed. Therefore, we use the machine epsilon and define as
Since we use the central difference approximation our order is 2 ().
Writing the Python script
So let us go over the Python script.
- Import the dependencies CoolProp and numpy
- Input parameter: Nitrogen, 100 bar, 300 K
- Get the machine epsilon
- Specify the function specificHeat
- Call the own and the CoolProp function
#License: MIT License (http://opensource.org/licenses/MIT)
import numpy as np
import CoolProp.CoolProp as CP
fluid = "Nitrogen"
eps = np.finfo(float).eps
def specificHeat(h, p):
pert = pow(eps, 0.33333)*h
num = (h + pert)-(h - pert)
denum = CP.PropsSI("T", "P", p, "H", h+pert, fluid)-CP.PropsSI("T", "P", p, "H", h-pert, fluid)
p = 100e5 #Pa
T = 300 #K
h = CP.PropsSI("H", "P", p, "T", T, fluid) #J/kg
cp = specificHeat(h, p)
cp_DerivativeCP = CP.PropsSI('d(H)/d(T)|P','P',p,'T',T,fluid)
First of all, numerical derivatives are quite often needed when dealing with engineering problems. Since derivatives can not be always calculated analytically we need numerical methods.
Furthermore, the CoolProp function returns = 1194.93427991 J/kg, while the own function returns = 1194.93427992 J/kg. This means, the results are identical. Hence, the internal CoolProp function is equal to the own approach. For example, using a forward difference in the own function would lead to a = 1194.93356274 J/kg. This results differs at the third decimal place.
- Be careful which differencing method you choose close to saturated liquid or saturated vapor line.
- In the two-phase region tends to infinity. This is because a pure component is not changing its temperature at constant pressure. A temperature glide can occur for mixtures.
What did you think?
I’d like to hear what you think about this post.
Let me know by leaving a comment below.