Always finding non-negative solution for P and T

[HBShaoUFZ]

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: https://github.com/oemof/tespy/blob/27dd934fa8035b6fa6ba8a658ef584d3ff8443cc/tespy/networks.py#L1175-L1184 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.

[fwitte]

Hey, thanks for your hint, I will have a look into that in the next week or two! Have a nice day!

[fwitte]

Hi,

I have started looking into it and tried to implement the relaxation factors (I pushed the changes to the features/improve_convergence_stability branch: 9fa04b304b7661a7634ecfa5e5abd94167972506).

Unfortunately I was not always able to find feasible results without the manipulation of the pressure values as implemented at the moment. Negative values for pressure are prevented, but in some cases the values exceed upper (a lot higher than reasonable) or fall below lower bounds of the allowed range in the fluid property database (CoolProp).

Also, I implemented the relaxation factor for the fluid mass fractions. It works fine for the lower bounds, but analogously the upper values are not limited (in this case limit should be 1 = 100%).

Maybe we can talk about these issues in detail, if you have some spare time.

Thanks a lot for your suggestion!

[fwitte]

I merged an improved convergence check into dev amongst other new features. Now the convergence check makes sure that the fluid properties of pure fluids always remain within the CoolProp boundaries. For mixtures the user defined ranges are still applied. For the moment, this seems to work very well.

I applied your consideration for pressure, which is working very good. For mass flow and enthalpy it is possible to have negative values, thus I did not change the newton iteration for those properties. I want to use this for fluid composition, too, but have not looked into maximum allowed values at the moment, as the mass fraction of a fluid in a mixture should always be less or equal than 1.

The corresponding lines are https://github.com/oemof/tespy/blob/79a1177d5db130465676907d9a51b680005bdc7a/tespy/networks.py#L1270-L1279, https://github.com/oemof/tespy/blob/79a1177d5db130465676907d9a51b680005bdc7a/tespy/networks.py#L1296 and https://github.com/oemof/tespy/blob/79a1177d5db130465676907d9a51b680005bdc7a/tespy/networks.py#L1335-L1380

If you have an idea, how to limit the maximum value as well as the minimum value for the fluid composition, you are very welcome to post it here!

Have a nice week!