1

I want to solve the following problem numerically and first defined it as a function:

def F(x, na=1, nb=2, L=6, A=0.09, V=0.54, dh=0.03, k=10**-2, dB=0.005, fB=38, ur=0.39, rhoW=998, rhoG=1.2, g=9.81, pu=10**5, myW = 10**-3):
    
    #variables
    uWa = x[0]
    uWb = x[1]
    ta = x[2]
    tb = x[3]
    ea = x[4]
    eb = x[5]
    rhoa = x[6]
    rhob = x[7]
    ua = x[8]
    ub = x[9]
    pa0 = x[10]
    pb0 = x[11]
    Rea = x[12]
    Reb = x[13]
    lambdaa = x[14]
    lambdab = x[15]

    #nonlinear equations in form f=0:
    f = np.zeros(16)
    f[0] = ta - L/(ur+uWa)
    f[1] = tb - L/(ur-uWb)
    f[2] = ea - na*ta*fB*np.pi*dB**3/(6*V)
    f[3] = eb - nb*tb*fB*np.pi*dB**3/(6*V)
    f[4] = rhoa - rhoW + ea*(rhoW-rhoG)
    f[5] = rhob - rhoW + eb*(rhoW-rhoG)
    f[6] = rhoa*ua - ea*rhoG*(ur+uWa) - (1-ea)*rhoW*uWa
    f[7] = rhob*ub + eb*rhoG*(ur-uWb) - (1-eb)*rhoW*uWb
    f[8] = (-rhoa*ua + rhob*ub + rhoG*(na+nb)*fB*np.pi*dB**3/(6*A**2))*A
    f[9] = pb0 - pa0 - rhoa*rhob*(ua**2-ub**2)/(rhoa+rhob)
    f[10] = pa0 - pu - rhoa*((ua**2*lambdaa*dh*np.pi)/(8*A)+g)*L    
    f[11] = pb0 - pu + rhob*((ub**2*lambdab*dh*np.pi)/(8*A)+g)*L
    f[12] = Rea - rhoa*ua*dh/myW
    f[13] = Reb - rhob*ub*dh/myW  
    f[14] = lambdaa - fixpunktiteration(colebrook, 0.05, args=(Rea,k,dh)) 
    f[15] = lambdab - fixpunktiteration(colebrook, 0.05, args=(Reb,k,dh))

    return f

Following functions I defined for the equations f[14] and f[15]:

def colebrook(lam,args):
    Re, k, d = args
    return ( 1/(4*np.log10( 2.51/(Re*lam**0.5) + (k/d)/3.71 )**2) )

def fixpunktiteration(f,xalt,args,kmax=500,tol=1.e-7):
    for k in range(1,kmax):
        xneu = f(xalt,args)
        xdiff = xneu - xalt
        if abs(xdiff/xneu) < tol:
            break
        xalt = xneu
    else:
        xneu = None
    return xneu

Then I defined all the parameters and initial values of the function F and tried to solve with fsolve (after importing it from scipy):

xinit = np.array([uWa_0, uWb_0, ta_0, tb_0, ea_0, eb_0, rhoa_0, rhob_0, ua_0, ub_0, pa0_0, pb0_0, Rea_0, Reb_0, lambdaa_0, lambdab_0]) 
loesung = fsolve(func=F,x0=xinit)

It works without error but it returns all initial values as solution which I don't understand. If I substitute f[14]=64/Rea the program gives a solution which is different from the initial values but physical nonsense. So I think the problem must have something to do with equations f[14] and f[15] or colebrook but I don't find it.

mkrieger1
  • 19,194
  • 5
  • 54
  • 65
Nikolas
  • 11
  • 1
  • - what is the numerical value of each entry in `xinit`? - have you tried the `full_outputbool`option (from docs)[https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html] to see what additional output might say? – neonwatty Sep 21 '22 at 22:32
  • 1
    Add the argument `full_output=True` to the call of [`fsolve`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html), and inspect the return value `ier`. If `ier` is not 1, it means `fsolve` failed. In that case, print the value of `mesg`. – Warren Weckesser Sep 21 '22 at 22:42
  • 1
    Why do you employ a fixed-point iteration inside an iterative numerical solver? Why not just give the fixed-point equations in f[14] and f[15]? – Lutz Lehmann Sep 22 '22 at 06:02
  • For ier I get "5" an for mesg "The iteration is not making good progress, as measured by the improvement from the last ten iterations." – Nikolas Sep 23 '22 at 19:07
  • @neonwatty this are all numerical values i defined in the main programm Rea_0 = 3000 Reb_0 = 3000 na = 2 nb = 1 s = 0.03 #m b = 0.03 #m L = 6 #m A = s*b #m**2 querschnittsflaeche V = s*b*L #m**3 volumen dh = 2*s*b/(b+s) #m hydraulischer durchmesser k = 5*10**-2 dB = 0.005 #m fB = 38 #Hz ur = 0.39 #m/s uWa_0 = 0.1 uWb_0 = 0.1 ta_0 = 10.5 tb_0 = 20.5 – Nikolas Sep 23 '22 at 19:12
  • ea_0 = 0.01 eb_0 = 0.01 rhoa_0 = 900 rhob_0 = 900 ua_0 = 0.15 ub_0 = 0.15 pa0_0 = rhoa_0*9.81*L pb0_0 = rhob_0*9.81*L lambdaa_0 = 0.05 lambdab_0 = 0.05 – Nikolas Sep 23 '22 at 19:12
  • @Lutz Lehmann I think fixed point iteration is the fastest way to solve equations of f[14] and f[15], first I tried as you proposed but I got for solution the inital values and "RuntimeWarning: invalid value encountered in double_scalars f[14] = 1/lambdaa**0.5 + 2*np.log10(2.51/(Rea*lambdaa**0.5) + (k/dh)/3.71)", same for f[15] – Nikolas Sep 23 '22 at 19:16
  • Yes, that, domain violation, is one reason to stay with an explicit solution. Sometimes one can find a parameter change or a transformation of the equation that removes such a domain restriction. – Lutz Lehmann Sep 23 '22 at 19:51
  • 1
    @Nikolas, the Colebrook equation can be solved in terms of [Lambert's W function](https://en.wikipedia.org/wiki/Lambert_W_function). SciPy has `scipy.special.lambertw`, and I have an implementation of the Colebrook function in [colebrook.py](https://github.com/WarrenWeckesser/experiments/blob/master/python/scipy/colebrook/colebrook.py). – Warren Weckesser Sep 28 '22 at 11:35

0 Answers0