Reputation: 51
Solving non-linear system of equations (issue with sqrt)
I am trying to solve a system of non-linear equations using Python (under-expanded jet with losses - Molkov):
I am trying to use the library Scipy and the module fsolve:
def equations10(p):
u3, T3 = p;
return ((u3 / ABEL_NOBLE_COEF) * (1 - ((np.sqrt(H2_GAMMA * R_H2 * T3)) / u3)) * (np.sqrt((R_H2 * T3) / H2_GAMMA) + ((K + F) / 4) * np.sqrt((2 / K) * (Cp_H2 * TANK_TEMP - Cp_H2 * T3 - ((u3**2) / 2) * ((F / 4) + 1) )) + u3 * ((F / 4) + 1)) - TANK_PRES,
(u3 / ABEL_NOBLE_COEF) * (1 - ((np.sqrt(H2_GAMMA * R_H2 * T3)) / u3)) * (((R_H2 * (TANK_TEMP - (TANK_TEMP - T3 - ((u3**2) / (2 * Cp_H2)) * ((F / 4) + 1)) * ((K + 1) / K))) / (np.sqrt((2 / K) * (Cp_H2 * TANK_TEMP - Cp_H2 * T3 - (u3**2 / 2) * ((F / 4) + 1))) - u3 + np.sqrt(H2_GAMMA * R_H2 * T3))) + ((K / 4) + 1) * np.sqrt((2 / K) * (Cp_H2 * TANK_TEMP - Cp_H2 * T3 - (u3**2 / 2) * ((F / 4) + 1)))) - TANK_PRES)
u3, T3 = fsolve(equations10, (1, 1))
print(equations10((u3, T3)))
print(u3, T3)
res = least_squares(equations10, (10**7, 300), bounds = ((10**4, 200), (10**10, 3000)))
print(res)
I encounter a problem due to the np.sqrt((2 / K) * (Cp_H2 * TANK_TEMP - Cp_H2 * T3 - ((u3**2) / 2) * ((F / 4) + 1) )) part in the last equation. Python is trying to solve the system trying a negative value for T3 (a mathematic error due to negative term with the square-root):
(-40784694.46088256, 48695678.683380365)
1.0 1.0
c:\users\guillaume.thiriet\downloads\test_h2jet.py:65: RuntimeWarning: invalid value encountered in sqrt
return ((u3 / ABEL_NOBLE_COEF) * (1 - ((np.sqrt(H2_GAMMA * R_H2 * T3)) / u3)) * (np.sqrt((R_H2 * T3) / H2_GAMMA) + ((K + F) / 4) * np.sqrt((2 / K) * (Cp_H2 * TANK_TEMP - Cp_H2 * T3 - ((u3**2) / 2) * ((F / 4) + 1) )) + u3 * ((F / 4) + 1)) - TANK_PRES,
c:\users\guillaume.thiriet\downloads\test_h2jet.py:66: RuntimeWarning: invalid value encountered in sqrt
(u3 / ABEL_NOBLE_COEF) * (1 - ((np.sqrt(H2_GAMMA * R_H2 * T3)) / u3)) * (((R_H2 * (TANK_TEMP - (TANK_TEMP - T3 - ((u3**2) / (2 * Cp_H2)) * ((F / 4) + 1)) * ((K + 1) / K))) / (np.sqrt((2 / K) * (Cp_H2 * TANK_TEMP - Cp_H2 * T3 - (u3**2 / 2) * ((F / 4) + 1))) - u3 + np.sqrt(H2_GAMMA * R_H2 * T3))) + ((K / 4) + 1) * np.sqrt((2 / K) * (Cp_H2 * TANK_TEMP - Cp_H2 * T3 - (u3**2 / 2) * ((F / 4) + 1)))) - TANK_PRES)
c:\users\guillaume.thiriet\downloads\test_h2jet.py:69: RuntimeWarning: The iteration is not making good progress, as measured by the
improvement from the last ten iterations.
u3, T3 = fsolve(equations10, (1, 1))
I don't know how to solve the problem. The idea would be to be able to set-up boundaries (min / max values) but it seems fsolve doesn't allow that. The least-squares function of Scipy does not give me an appropriate result. The parameters used in the calculation:
ABEL_NOBLE_COEF = 0.00769 # [m3.kg-1]
H2_GAMMA = 1.4 # [-]
R_GAS = 8.31446261815324 # [J.K-1.mol-1]
M_H2 = 0.002016 # [kg.mol-1]
R_H2 = R_GAS / M_H2 # [J.K-1.kg-1]
Cp_H2 = 14310 # [J.kg-1.K-1]
K = 0.5
D = 0.00075 # [m]
L = 0.015 # [m]
f = 0.049
F = f*L/D
T1 = 300 # [K]
P1 = 400e5 # [Pa]
f is calculated using the following equation (and a Reynolds number of 2000 - arbitrary choice):
The appropriate set of values to be studied should be like that:
u = np.linspace(0, 1500, 200) # Velocity [m.s-1]
t = np.linspace(0, 500, 200) # Temperature [K]
The set of two equations come from a larger set of eight equations / eight unknowns (I am trying to solve these 8 more simple equations but still have some issues):
Reference of the paper I am working with:
- V. Molkov - "Hydrogen jet flames" - International Journal of Hydrogen Energy - 2013
- V. Molkov - "Fundamentals of Hydrogen Safety I"
- V. Molkov - "High pressure hydrogen leak through a narrow channel"
Expected results:
For those who are interested, there is the theory without losses
Upvotes: 1
Views: 248
Answers (3)
Reputation: 6480
OK, for what it's worth I had a go by a completely different method - a bit more physically based. I'm coming up with slightly different results to you - but I'm using the equation of state (and the speed of sound) given in the Molkov paper as Table 1. (Note that the speed of sound is very slightly different in your additional "theory with no losses" text, although it doesn't make much difference here, as b.rho is small.)
From what I see, you start with reservoir conditions P1 and T1. The flow expands with a minor inlet loss (coefficient K) and pipe-friction loss (coefficient F). I have ignored variations in F since your Reynolds number inevitably changes along the pipe and Nikuradze certainly didn't do his pipe-flow measurements in a shock tube.
The key information is that the flow is choked (goes sonic: flow velocity = sound speed sqrt( gamma.p/rho) ) at position 3. There are no further losses to position 4, so I've ignored that.
The fact that the flow is sonic at position 3 allows you to find rho.u^2/p there, and hence find a relationship between the dynamic and static part of the pressure term. Adding the equation of state gives also u^2/T in terms of purely thermodynamic properties here. Then you can do similarly at position 2 using continuity (rho2.u2=rho3.u3).
So, I come up with expressions for P2/P1, P3/P2, T2/T1 and T3/T2 and I simply iterate through those until I reach a steady solution.
Note that if there are no losses then you can ignore state 2 and the connection between 1 and 3, with 3 being sonic, allowing you to find an exact solution without iteration (for the speed of sound as given in the Molkov paper). I've added the analytical solution with no losses to the program.
I'm getting slightly different values to you. Note also that my mass flow rate is in kg/s, not g/s.
Code:
import math
b = 0.007691 # (dimensional) coefficient in Abel-Noble equation of state: p = rho.R.T/(1-b.rho)
R = 8.31446261815324 / 0.002016 # gas constant for hydrogen
gamma = 1.4 # ratio of specific heats (diatomic gas)
cp = R * gamma / ( gamma - 1 ) # specific heat capacity at constant pressure
def solver( P1, T1, F, K ):
# start with uniform pressure and temperature
P2, P3, P4 = P1, P1, P1
T2, T3, T4 = T1, T1, T1
# Iteration
niter = 0
for _ in range( 10000 ):
niter += 1
B2 = b * P2 / ( R * T2 ) # Additional term where Abel-Noble EOS deviates from ideal-gas law
B3 = b * P3 / ( R * T3 )
# Everything fixed by sonic conditions at location 3
# Note that c = sqrt(gamma.R.T/(1-b.rho) = sqrt(gamma.P/rho) in original paper, but c = sqrt(gamma.P/rho/(1-b.rho)) in alternative
ruup3 = gamma # rho.u^2/p at 3, by sonic conditions
# rho3 = P3 / ( R * T3 + b * P3 ) # The speed of sound is (slightly) different in the second paper quoted.
# ruup3 /= ( 1 - b * rho3 ) # This doesn't make a massive difference, but uncomment these 2 lines to try
ruup2 = ruup3 * ( P3 / P2 ) ** 2 * ( T2 / T3 ) * ( 1 + B2 ) / ( 1 + B3 ) # using continuity and equation of state at 2 and 3
uuT3 = ruup3 * R * ( 1 + B3 ) # u^2/T, from the above and the equation of state
uuT2 = ruup2 * R * ( 1 + B2 )
P2new = P1 / ( 1 + ( 1 + K / 4 ) * ruup2 )
T2new = T1 / ( 1 + ( 1 + K ) * uuT2 / ( 2 * cp ) )
P3new = P2new * ( 1 + ( 1 - F / 4 ) * ruup2 ) / ( 1 + ( 1 + F /4 ) * ruup3 )
T3new = T2new * ( 1 + uuT2 / ( 2 * cp ) ) / ( 1 + ( 1 + F / 4 ) * uuT3 / ( 2 * cp ) )
change = abs( P2new - P2 ) + abs( P3new - P3 ) + abs( T2new - T2 ) + abs( T3new - T3 )
P2, P3, T2, T3 = P2new, P3new, T2new, T3new
if change < 1.0e-20: break
print( "Iterations: ", niter )
rho3 = P3 / ( R * T3 + b * P3 )
u3 = math.sqrt( gamma * P3 / rho3 )
return P3, T3, rho3, u3 # sonic (choked) exit conditions
P1 = 40.0e6 # stagnation pressure (Pa)
T1 = 287.65 # stagnation temperature (K)
K = 0.5 # minor loss coefficient at inlet
F = 0.037 * 0.015 / 0.00075 # pipe friction loss coefficient (dubious; not an incompressible or uniform flow)
A = math.pi * 0.00075 ** 2 / 4 # cross-sectional area of pipe (m^2)
P3, T3, rho3, u3 = solver( P1, T1, F, K )
massflow = rho3 * u3 * A
print( "Case with losses:" )
print( "Sonic conditions, P(Pa), T(K), rho(kg/m3) = ", P3, T3, rho3 )
print( "Mass flow rate (kg/s) = ", massflow )
#----------- Special case: no losses. Solution is analytic. -----------
print( "\nCase with no losses:" ) # This is analytical (assumes speed of sound = sqrt(gamma.P/rho) as in original post
P = P1 / ( 1 + gamma )
T = ( T1 - gamma * b * P / ( 2 * cp ) ) / ( 1 + gamma * R / ( 2 * cp ) )
rho = P / ( R * T + b * P )
u = math.sqrt( gamma * P / rho )
massflow = rho * u * A
print( "Sonic conditions, P(Pa), T(K), rho(kg/m3) = ", P, T, rho )
print( "Mass flow rate (kg/s) = ", massflow )
Output:
Iterations: 48
Case with losses:
Sonic conditions, P(Pa), T(K), rho(kg/m3) = 13598666.045112502 219.30168416349343 13.476841457366842
Mass flow rate (kg/s) = 0.00707650210058704
Case with no losses:
Sonic conditions, P(Pa), T(K), rho(kg/m3) = 16666666.666666668 234.52825111058513 15.214676526739288
Mass flow rate (kg/s) = 0.008324001743620746
Upvotes: 2
Reputation: 11002
Solutions
Thank you for publishing a Minimal Complete Verifiable Example.
For your setup, assuming P1=TANK_PRES and T1=TANK_TEMP, and with first version of:
D = 0.001 # [m]
L = 0.001 # [m]
It seems there is at least three solutions for your system:
But only two are accessible by the optimizer:
sol1 = optimize.fsolve(equations10, x0=[1275, 220])
# array([1267.63445497, 221.36788615])
sol2 = optimize.fsolve(equations10, x0=[1225, 210])
# array([1233.28256049, 207.89160204])
The last one lies on the region where square roots make computation fails:
sol3 = optimize.fsolve(equations10, x0=[1325, 240])
# RuntimeWarning: The iteration is not making good progress, as measured by the improvement from the last ten iterations.
If it's this root you need for there might have another solution to catch it numerically. Let me know in the comment.
Modified parameters
For setup:
D = 0.00075 # [m]
L = 0.015 # [m]
There is al least two solutions, one can be found by the solver:
sol1 = optimize.fsolve(equations10, x0=[1275, 220])
#array([1250.20118685, 222.56397776])
I leave the second one as a bonus:
# (1275.0218781410672, 228.98105595734015)
Which has been found by finding the intersection of the isopleths numerically using Shapely.
Your use case without knowing your constants
I'll provide the following set of constants (please update your post with yours, it's really missing for working correctly). I hope they are enough close to reality:
R = 8.31
ABEL_NOBLE_COEF = 0.5
H2_GAMMA = 1.40
R_H2 = R / 2.
K = 1e5
F = 1e1
Cp_H2 = 7/2 * R
TANK_TEMP = 50. + 273.15
TANK_PRES = 3e5
For this set of values, over variable ranges:
u = np.linspace(0, 100, 200)
t = np.linspace(0, 300, 200)
Values of your function equations10 do have zeros separately:
U, T = np.meshgrid(u, t)
X, Y = equations10([U, T])
Lets inspect the shape of the functions (dashed line are zeros):
Now the question is, do they have simultaneous zeros (dashed line is second equation) ?
fig, axe = plt.subplots()
axe.contour(U, T, X, [0], linestyles="-")
axe.contour(U, T, Y, [0], linestyles="--")
axe.grid()
As you can see, for this constant setup, zero isopleths expose quite same behaviours but are a little bit different. The question is then do they intersect?
Zooming on the region where such an intersection seems to be possible, we get not evidence of intersection:
Caution: This is numerical analysis, it inherently has inaccuracies, this does not prove there is no solution.
Anyway it provides a method to check if your equation exhibits an evident solution (zero isopleth intersections). In this setup, it seems the zero isopleth do not intersect, hence there is no solution for such system with the given parameters.
Simpler use case
Lets check it for a simpler system to fix ideas:
def system(x):
return np.array([
x[0] - x[1] + 1.,
x[0] - x[1] + 1.1
])
optimize.fsolve(system, x0=[1, 1])
# Raise warning: The iteration is not making good progress, as measured by the improvement from the last ten iterations.
xlin = np.linspace(-5, 5, 200)
ylin = np.linspace(-5, 5, 200)
X, Y = np.meshgrid(xlin, ylin)
U, V = system([X, Y])
fig, axe = plt.subplots()
axe.contour(X, Y, U, [0], linestyles="-")
axe.contour(X, Y, V, [0], linestyles="--")
The solution is similar to yours, both isopleth are alike but a bit different, so the solver cannot solve the system and warns there is no amelioration in its algorithm.
If we change a bit the system to have a solution:
def system(x):
return np.array([
x[0] - x[1] + 1.,
x[0] - 2 * x[1] + 1.1
])
Optimizer find it easily:
optimize.fsolve(system, x0=[1, 1])
# array([-0.9, 0.1])
And graphical method agrees with this solution:
Conclusion
Without your exact setup, it's simply not possible to tell what is really happening. Update your post accordingly if you wish to continue the analysis.
The method fsolve may return this warning for different reasons, but in this case it seems it is because there is no solution: zero isopleths do not intersects.
After analyzing multiple setup, we can say that your system of equations may none or multiple roots. Indeed for such complex problem fsolve needs a good educated guess to ensure convergence. Here we used numerical geometry (2D space are fine for this) to find such guess.
Upvotes: 4
Reputation: 51
Ok, it seems that the fsolve function may run properly as soon as the appropriate parameters are correctly set up.
In my case, it was difficult to calculate the derivative of the inserted function so I decided to just set up the starting estimate of the roots (as explained in the scipy doc) https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html#scipy.optimize.fsolve
It appears that I had not a correct idea of what my root may be and the solver was not able to run properly.
Then, I encountered very frequently the following error
RuntimeWarning: invalid value encountered in sqrt
This could be easily solved by adding np.abs(X) around the problematic structure like that np.sqrt(np.abs(X)). This method had been proposed by someone above. Thanks to him, it helps a lot.
Finally, I chose to solve the set of 12 equations / 12 variables instead of only the set of two equations / two variables because it seemed easier as it allows to reduce the risk of human errors during the writing of the equation.
The final program is presented below. As far as I know, it works for what I need to do.
import numpy as np
from scipy.optimize import fsolve
""" VARIABLES """
ABEL_NOBLE_COEF = 0.00769 # [m3.kg-1]
H2_GAMMA = 1.4 # [-]
R_GAS = 8.31446261815324 # [J.K-1.mol-1]
M_H2 = 0.002016 # [kg.mol-1]
R_H2 = R_GAS / M_H2 # [J.K-1.kg-1]
Cp_H2 = 14310 # [J.kg-1.K-1]
Diam3 = 0.00075 # [m]
Length = 0.015 # [m]
f = 0.037
F = f*Length/Diam3
K = 0.5
T1 = 287.65 # [K]
P1 = 400e5 # [Pa]
Rho1 = 1 # [kg.m-3]
u1 = 0 # [m.s-1]
# T2 = 1
# P2 = 1
# Rho2 = 1
# u2 = 1
T3 = 1
P3 = 1
Rho3 = 1
u3 = 1
T4 = 1
P4 = 101325
Rho4 = 1
u4 = 1
A4 = 1 # [m2]
""""" WITHOUT LOSSES """""
""" STATE 1 """
rho1 = P1 / (ABEL_NOBLE_COEF * P1 + R_H2 * T1)
""" rho3 """
def equation1(p):
rho3 = p
return ( ((rho3 / (1 - ABEL_NOBLE_COEF * rho3))**(H2_GAMMA - 1)) * (1 + ((H2_GAMMA - 1) / (2 * ((1 - ABEL_NOBLE_COEF * rho3)**2)))) - ((rho1 / (1 - ABEL_NOBLE_COEF * rho1))**(H2_GAMMA - 1)) )
rho3 = fsolve(equation1, 10)[0]
""" T3 """
def equation2(p):
T3 = p
return ( 1 + ((H2_GAMMA - 1) / (2 * ((1 - ABEL_NOBLE_COEF * rho3)**2))) - (T1 / T3) )
T3 = fsolve(equation2, 300)[0]
""" P3 """
def equation3(p):
P3 = p
return ( (P3 / (ABEL_NOBLE_COEF * P3 + R_H2 * T3)) - rho3 )
P3 = fsolve(equation3, 101325)[0]
""" u3 """
u3 = np.sqrt((H2_GAMMA * P3) / (rho3 * (1 - ABEL_NOBLE_COEF * rho3)))
""" T4 """
T4 = ((2 * T3) / (H2_GAMMA + 1)) + ((H2_GAMMA - 1) / (H2_GAMMA + 1)) * (P3 / (rho3 * (1 - ABEL_NOBLE_COEF * rho3) * R_H2))
""" u4 """
u4 = np.sqrt(H2_GAMMA * R_H2 * T4)
""" rho4 """
rho4 = (P4 / (R_H2 * T4))
""" Diam4 / A4 """
Diam4 = Diam3 * np.sqrt((rho3 * u3) / (rho4 * u4))
A4 = np.pi * (Diam4 / 2)**2
Q = rho4 * u4 * A4 * 1000
print("Q NoLosses = " + str(Q))
""""" WITH LOSSES """""
""" REYNOLDS """
H2_Visc_3 = 7.3e-06
#H2_Visc_3 = 8.3969e-05
Re3 = (rho3 * u3 * Diam3) / H2_Visc_3
""" Friction factor - f """
def equation4(p):
f = p
return ( (1 / np.sqrt(f)) - 0.869 * np.log(Re3 * np.sqrt(f)) + 0.8 )
f = fsolve(equation4, 0.01)[0]
F = f*Length/Diam3
print("F coef = " + str(f))
""" EQUATIONS """
def equation5(p):
P2, T2, u2, rho2, P3, T3, u3, rho3, T4, u4, rho4, A4 = p;
return (P2 - P1 + rho2 * (u2**2) * ((K / 4) + 1),
Cp_H2 * T1 - Cp_H2 * T2 - ((K + 1) * (u2**2) / 2),
P2 - ((rho2 * R_H2 * T2) / (1 - ABEL_NOBLE_COEF * rho2)),
P3 - P2 + rho2 * (u2**2) * ((F / 4) - 1) + rho3 * (u3**2) * ((F / 4) + 1),
Cp_H2 * T2 + ((u2**2) / 2) - Cp_H2 * T3 - ((F / 4) + 1) * ((u3**2) / 2),
rho3 - (P3 / (ABEL_NOBLE_COEF * P3 + R_H2 * T3)),
rho2 * u2 - rho3 * u3,
u3 - np.sqrt((H2_GAMMA * R_H2 * np.abs(T3)) / (1 - ABEL_NOBLE_COEF * rho3)),
u4 - np.sqrt((H2_GAMMA * R_H2 * np.abs(T4))),
rho4 - (101325 / (R_H2 * T4)),
Cp_H2 * T3 + ((u3**2)/2) - Cp_H2 * T4 - ((u4**2) / 2),
rho3 * u3 * np.pi * ((Diam3 / 2)**2) - rho4 * u4 * A4)
P2, T2, u2, rho2, P3, T3, u3, rho3, T4, u4, rho4, A4 = fsolve(equation5, (1e7, 200, 1100, 0.1, 1e7, 200, 1100, 0.1, 200, 1100, 0.1, 1e-04))
print(P2, T2, u2, rho2, P3, T3, u3, rho3, T4, u4, rho4, A4)
print("*****")
Q = rho4 * u4 * A4 * 1000
print("Q WithLosses = " + str(Q))
Thanks all for your support and your help.
Regards
Upvotes: 1
Related Questions
- How to fix "Attempted relative import in non-package" even with __init__.py
- Constraining conditional summation in Bounded Optimization in Python
- TypeError when using scipy.optimize.curve_fit with a piecewise function using piecewise functions
- Using openpyxl and scipy to solve non-linear system related to resistor network
- Balancing a System of Non-Linear equations for Chemical Reactions
- Solving single non-linear equation involving sum of numpy array with fsolve
- Solving for a single (non linear) equation in a for loop with substitution?