While testing a benchmark for TESPy, we encountered a problem of non-convergence. The error messages from the TESPy are like follows:
ValueError: ('unable to solve 1phase PY flash with Tmin=273.16, Tmax=3000 due to error: HSU_P_flash_singlephase_Brent could not find a solution because Hmolar [0.760517 J/mol] is below the minimum value of 45064.4615885 J/mol : PropsSI("S","P",2,"H",42.21513643,"water")', 'occurred at index <tespy.components.components.pump object at 0x000001CC019BFB38>')
-------------------------------------------------------------------------------------------------------------------After analyzing this error, we realized that during the Newton iteration, a negative pressure value was passed to the CoolProp library, when the enthalpy value for water was needed. In the real physical world, the negative pressure will never happen. This is a numerical issue caused by the Newton iteration, when the primary variable p and T are updated.
By looking into the code and communication with the author:
|
# check properties for consistency |
|
if self.iter < 5: |
|
for c in self.conns.index: |
|
self.solve_check_properties(c) |
|
|
|
for cp in self.comps.index: |
|
cp.convergence_check(self) |
|
|
|
for c in self.conns.index: |
|
self.solve_check_properties(c) |
We learned that the author is aware of this issue, thus decided to add some protection in the code above when solving the process. The author's method is to check primary unknowns and parameters in the first 3 Newton iterations. If either the primary variable or the parameter is beyond the valid range, they are then set to the lower or upper bound set externally by the user. By changing the 3 iterations to 6 iterations, our error message disappeared and the benchmark converged without any further issue. However, I personally believe this is not the best solution to prevent Pressure and Temperature values entering the negative range during Newton iterations.
This non-negative problem also exists in chemical equilibrium calculation, where the primary variables there are the concentrations of each chemical component. Same as the pressure (in Pa) and temperature (in Kelvin), concentrations are always non-negative. During the Newton iteration, the concentrations are never updated directly, but instead with a relaxation factor. Assuming we have the standard newton update procedure as,
$$c_i^{n+1} = c_i^{n} + \Delta c_i$$
The proposed relaxation factor $\delta$ is multiplied before the suggested change
$$c_i^{n+1} = c_i^{n} + \delta \cdot \Delta C_i,$$
and the calculation of relaxation factor $\delta$ follows,
$$ \frac{1}{\delta} = max{ 1, -\frac{\Delta C_i}{0.5 * C_i^{n}} } $$
Assuming we have a concentration value C_i^n equals to 20 mol/kg of water, and the calculated change would be -50 mol/lg of water, then the updated C_i^{n+1} would be -30 if no relaxation is applied. When inserting these values, the relaxation factor can be calculated as 0.2. With this factor multiplied, the updated concentration C_i^{n+1} would be 20 + 0.2 * (-50) = 10 instead of -30 mol/kg water. With the presence of this relaxation factor, the concentration value will be always non-negative through out the Newton iteration.
As the above algorithm has been widely adopted in chemical reaction simulation and has been proven to be very effective, I would suggest to adopt the same treatment in TESPy when updating pressure and temperature values. Just my suggestion for consideration.