Solve a set of nonlinear equations for a discrete value in SymPy

Question

I'm trying to solve a system of nonlinear equations for a discrete value in SymPy. The system constists of 7 equations.

Here's my code:

import math
from CoolProp.CoolProp import PropsSI
import sympy as sym

alpha_Fluid, Re_l, v_Fluid, mdot_Fluid, Nu_Fluid, Nu_l_lam, Nu_l_turb = sym.symbols('alpha_Fluid, Re_l, v_Fluid, mdot_Fluid, Nu_Fluid, Nu_l_lam, Nu_l_turb')

################################ Discrete Parameters ################################
fluid = 'Nitrogen';
p_Fluid = 1.1e5; 
T_Fluid = 80; 
T_ZP = 300; 
Pr_Fluid = PropsSI("PRANDTL",'P|gas', p_Fluid, 'T', (T_Fluid+T_ZP)/2, fluid);
eta_Fluid = PropsSI("VISCOSITY",'P|gas', p_Fluid, 'T', (T_Fluid+T_ZP)/2, fluid);
lambda_Fluid = PropsSI("CONDUCTIVITY",'P|gas', p_Fluid, 'T', (T_Fluid+T_ZP)/2, fluid);
rho_Fluid = PropsSI("DMASS",'P|gas', p_Fluid, 'T', (T_Fluid+T_ZP)/2, fluid);
Bi = 0.1;
r_ZP = 7e-3; 
lambda_ZP = 230; 
L_ZP_Biot = r_ZP; 
L_ZP_KONV = (math.pi/2)*2*r_ZP; 
A_Inflow = 4e-3; 
A_Anstr = 2*math.pi*r_ZP*A_Inflow; 


Eq1 = sym.Eq((alpha_Fluid*L_ZP_Biot)/lambda_ZP, Bi)
Eq2 = sym.Eq((alpha_Fluid*L_ZP_KONV)/lambda_Fluid, Nu_Fluid) 
Eq3 = sym.Eq(0.664*(Re_l)**(1/2)*Pr_Fluid**(1/3),Nu_l_lam)
Eq4 = sym.Eq((0.037*Re_l**(0.8)*Pr_Fluid)/(1+2.443*Re_l**(-0.1)*(Pr_Fluid**(2/3)-1)), Nu_l_turb)
Eq5 = sym.Eq(0.3 + (Nu_l_lam**2 + Nu_l_turb**2)**(1/2), Nu_Fluid)
Eq6 = sym.Eq((rho_Fluid*v_Fluid*L_ZP_KONV)/eta_Fluid, Re_l)
Eq7 = sym.Eq((mdot_Fluid*rho_Fluid)/A_Anstr, v_Fluid)

result = sym.solve([Eq1, Eq2, Eq3, Eq4, Eq5, Eq6, Eq7],mdot_Fluid)
print(result)

I expect a discrete value for the variable 'mdot_Fluid'.

But the solution I obtain is the following: {mdot_Fluid: 8.99341199537194e-5*v_Fluid}.

My first guess is, that the whole system is underdetermined, which is weird, because I have 7 equations and 7 unknows.

Answer 1

On the basis of your responses I solved the problem with the scipy.fsolve method. I made the mistake to solve for only the variable 'mdot_Fluid', while I should have solved for the other unknown variables too. Here's my solution:

import math
from scipy.optimize import fsolve

##### Constants ######
p_Fluid = 1.1e5; # 
T_Fluid = 80; 
T_ZP = 300; 
Pr_Fluid = 0.74;
eta_Fluid = 1.23e-5;
lambda_Fluid = 0.0174;
rho_Fluid = 1.96;
r_ZP = 7e-3; 
lambda_ZP = 230; 
L_ZP_Biot = r_ZP; 
L_ZP_KONV = (math.pi/2)*2*r_ZP; 
A_Inflow = 4e-3; 
A_Anstr = (2*math.pi*r_ZP*A_Inflow)/2; 
Bi = 0.01; 

def calc_massflux(g):
    alpha_Fluid, Re_l, v_Fluid, mdot_Fluid, Nu_Fluid, Nu_l_lam, Nu_l_turb = g;
    return [(alpha_Fluid*L_ZP_Biot)/lambda_ZP - Bi, 
            (alpha_Fluid*L_ZP_KONV)/lambda_Fluid - Nu_Fluid, 
            0.664*(Re_l)**(1/2)*Pr_Fluid**(1/3) - Nu_l_lam, 
            (0.037*Re_l**(0.8)*Pr_Fluid)/(1+2.443*Re_l**(-0.1)*(Pr_Fluid**(2/3)-1)) - Nu_l_turb,
            0.3 + (Nu_l_lam**2 + Nu_l_turb**2)**(1/2) - Nu_Fluid,
            (rho_Fluid*v_Fluid*L_ZP_KONV)/eta_Fluid - Re_l,
            mdot_Fluid/(rho_Fluid*A_Anstr) - v_Fluid]
            
alpha_Fluid, Re_l, v_Fluid, mdot_Fluid, Nu_Fluid, Nu_l_lam, Nu_l_turb = fsolve(calc_massflux, (1,1,1,1,1,1,1))
print(mdot_Fluid)

Answer 2

This can be solved as 6 linear equations and a 7th nonlinear. Solve all but equation 4 and then substitute the solution into equation 4 and use nsolve with an initial guess for v_Fluid to get an answer.

import math
import sympy as sym
from sympy import *
math=sym  # don't use math.pi use sympy.pi
alpha_Fluid, Re_l, v_Fluid, mdot_Fluid, Nu_Fluid, Nu_l_lam, Nu_l_turb = sym.symbols('alpha_Fluid, Re_l, v_Fluid, mdot_Fluid, Nu_Fluid, Nu_l_lam, Nu_l_turb')

################################ Discrete Parameters ################################
fluid = 'Nitrogen';
p_Fluid = 1.1e5; 
T_Fluid = 80; 
T_ZP = 300; 
Pr_Fluid = 0.74;
eta_Fluid = 1.23e-5;
lambda_Fluid = 0.0174;
rho_Fluid = 1.96;
Bi = 0.1;
r_ZP = 7e-3; 
lambda_ZP = 230; 
L_ZP_Biot = r_ZP; 
L_ZP_KONV = (math.pi/2)*2*r_ZP; 
A_Inflow = 4e-3; 
A_Anstr = 2*math.pi*r_ZP*A_Inflow; 


Eq1 = sym.Eq((alpha_Fluid*L_ZP_Biot)/lambda_ZP, Bi)
Eq2 = sym.Eq((alpha_Fluid*L_ZP_KONV)/lambda_Fluid, Nu_Fluid) 
Eq3 = sym.Eq(0.664*(Re_l)**(1/2)*Pr_Fluid**(1/3),Nu_l_lam)
Eq4 = sym.Eq((0.037*Re_l**(0.8)*Pr_Fluid)/(1+2.443*Re_l**(-0.1)*(Pr_Fluid**(2/3)-1)), Nu_l_turb)
Eq5 = sym.Eq(0.3 + (Nu_l_lam**2 + Nu_l_turb**2)**(1/2), Nu_Fluid)
Eq6 = sym.Eq((rho_Fluid*v_Fluid*L_ZP_KONV)/eta_Fluid, Re_l)
Eq7 = sym.Eq((mdot_Fluid*rho_Fluid)/A_Anstr, v_Fluid)
req = nsimplify(Tuple(*[Eq1, Eq2, Eq3, Eq5, Eq6, Eq7]))
result = sym.solve(req, dict=True)[0]
nl = Eq4.subs(result)
nsolve(nl.subs(reps), 1)

Alternatively, square both sides of nl, replace v_Fluid**.1 with a positive x and use real_roots on the result

>>> eq=nsimplify(nl.lhs**2-nl.rhs**2).as_numer_denom()[0].subs(root(v_Fluid,10),var('x',positive=True)).n()
>>> [(i**10).n() for i in real_roots(eq)] # x^(1/10) = i so x = i**10
[920.247842938558, 8.55655777550056e-8, 8.55656098540597e-8, 721.920699464287]

You will have to substitute these values into nl.lhs-nl.rhs and result to see which ones (if any) are actually solutions. The scaling of the equations is not what one would desire, however, with these numerical values. I would proceed with caution.