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

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.

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)