1
$\begingroup$

enter image description here

I am trying to solve this equation for $v_{Sl}$ in the range $[0.001,1]$ to obtain the values of $\delta_l$. I have tried Fsolve with python, which gives results where the values of $v_{Sl}$ are between $[0.001,0.1]$ when I put an initial guess of $\delta_l$ that is close to the solution. But, when I vary the initial guess I don't obtain any solution for the whole range.

I also tried scipy root with all the methods provided and I did not obtain any solution. When I tried with an Evolutionary Algorithm, such as Genetic Algorithm, it seems to give results, but with the limitation of the range of solution provided. For instance, if I choose the range of solution for $\delta_l$ to be from $[0.001,0.065]$, I get the solution I am actually looking for, but if I widen/narrow the range I obtain a different solution. Ideally, I would give a range of solutions that may fall in let's say $[0.0001,0.2]$ and the algorithm would give the same results as when I give a range of $[0.001,0.065]$.

From this I understand that the algorithm is not converging toward the solution of the equation ($0$) and I need a better tool to make it converge. Is there a better algorithm to solve this problem?

Code and Figures

As a matter of fact, the final solution I would be looking to obtain is $v_{sg}$. I will first obtain the $\delta_l$ value from solving the mentioned equation then calculate $v_{sg}$ from the following equation. enter image description here

The following graph is a reference which I need to reproduce the same graph with my code. Reference plot

Here is the graph that I generated.

My Plot

The code I wrote is below. The main problem is the range of solutions for $\delta_l$ that I need to specify to get the desired results, for my case here I put $[0.001,0.065]$, but I should not narrow it as I will not know the solution in another case (different Data). Rather I would put an acceptable range let's say $[0.0001,0.2]$ and would still expect to get similar results, which is not the case.

Plot of $\delta_l$ vs function I am solving:

enter image description here

import pygad
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
import math
###### Constants start here 
D= 0.06  ### Diameter (m)
ROL= 997.9 #### Density liquid  Kg/m3
ROG= 1.2  ##### density gas    kg/m3
UL= 1.1*(10**(-3))### viscosit Pa.S
UG= 0.000018  ### viscosit Pa.S
G= 9.81   # gravitational acceleration  m/s2
sigma= 0.06    #### Interfacial tension  N/m
alpha= 20   #### inclination from horizontal   Degrees
SL=3.14*D
### constants ends

##### list of the values of the variable that we are solving the equation with each time
list_of_liquid_velocities_VSL= np.linspace(0.001,1,80) #va
### ends here

### lists to store the output values in
list_of_derivatives=[]
list_of_interfacial_tension=[]
list_of_interfacial_tension_gas=[]  
list_of_VSG=[]   ###### main values we are looking for
### ends

##################### Main algorithm  ##############

for VSL in list_of_liquid_velocities_VSL: 
    AL=3.14*(D**2)*(0.01-0.01**2)
    
    DL=4*AL/SL
    A=3.14*(0.5*D)**2
    VL=VSL*A/AL
    if (D*ROL*VSL/(UL))<2300:
        n=1
        CL=16
        m=0.2
        CG=0.046    
    else:
        n=0.2
        CL=0.046
        m=0.2
        CG=0.046
        
        
    desired_output = 0
    
    def fitness_func(x,x_idx):
        derivative = ((G*(ROL-ROG)*D*math.sin(alpha*3.14/180)*((1-2*x)**2-2*(x-x**2)))-((CL/16)*ROL*((D*ROL/UL)**(-n))*(VSL**(2-n))*(((x-x**2)+(1-2*x)**2)/(x-x**2)**3)))  ### main equation I am solving

        fitness = (1/np.abs(derivative - desired_output))
        return fitness

    fitness_function = fitness_func

    num_generations = 100
    num_parents_mating = 4

    sol_per_pop = 8
    num_genes = 1
    init_range_low = 0.001
    init_range_high = 0.101

    parent_selection_type = "sss"
    keep_parents = 1

    crossover_type = "single_point"

    mutation_type = "random"
    mutation_percent_genes = 20

    gene_space = [{'low': 0.001, 'high': 0.065}] ### Range I am giving for the algorithm to look for the solution in (try to vary it and see the change in the final plot)
    ga_instance = pygad.GA(num_generations=num_generations,
                       num_parents_mating=num_parents_mating,
                       fitness_func=fitness_function,
                       sol_per_pop=sol_per_pop,
                       num_genes=num_genes,
                       init_range_low=init_range_low,
                       init_range_high=init_range_high,
                       parent_selection_type=parent_selection_type,
                       keep_parents=keep_parents,
                       crossover_type=crossover_type,
                       mutation_type=mutation_type,
                       mutation_percent_genes=mutation_percent_genes,
                       gene_space= gene_space)


    ga_instance.run()
    x, solution_fitness, solution_idx = ga_instance.best_solution()
    
   ### x is the solution of the function above 
    derivative = ((G*(ROL-ROG)*D*np.sin(alpha*3.14/180)*((1-2*x)**2-2*(x-x**2)))-((CL/16)*ROL*((D*ROL/UL)**(-n))*(VSL**(2-n))*(((x-x**2)+(1-2*x)**2)/(x-x**2)**3))) ### recalculation of the equation above to double confirm (which should should be equal to 0)


    list_of_derivatives.append(abs(derivative))
    
    interfacial_tenstion=((G*(ROL-ROG)*D*np.sin(alpha*3.14/180)*(x-x**2)*(1-2*x))+((CL/32)*ROL*(D*ROL/UL)**(-n)*VSL**(2-n)*((1-2*x)/(x-x**2)**2)))  ### should be equal to interfacial_tension_gas

    VSG  =  (interfacial_tenstion*((1-2*x)**4) /((0.5)*(CG)*((D*ROG/UG)**(-m))*(1+300*x)*(ROG)))**(1/(2-m)) ### VSG the main output which will be plotted against VSL in a loglog plot
    
    list_of_interfacial_tension.append(interfacial_tenstion)
    
    interfacial_tension_gas= ((((VSG)**(2-m))*(0.5)*(CG)*((D*ROG/UG)**(-m))*(1+300*x)*(ROG))/((1-2*x)**4))/(G*(ROL-ROG)*D)  ### this should equal to interfacial_tension
    list_of_interfacial_tension_gas.append(interfacial_tension_gas)
    
    list_of_VSG.append(VSG)
    
    
VSLexperiment= [0.02,0.06,0.1]
VSGexperiment=[13.94,14.54,19.09]

plt.plot()
plt.xscale("log")
plt.yscale("log")
plt.xlim(1,100)
plt.ylim(0.001,1)
plt.xlabel("VSG", fontsize=9)
plt.ylabel("VSL",fontsize=9)
plt.grid(True, which="both", ls=":", color='0.002', linewidth=0.4)
plt.scatter(VSGexperiment,VSLexperiment, color="Black", linewidth="0.05", label="VsgExperiment")
plt.plot (list_of_VSG,list_of_liquid_velocities_VSL, label='VBarnea', color='blue')  

I tried with approach of scipy.root with the following codes:

from numpy import sin, linspace
from numpy.polynomial import Polynomial
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
import math

###### Constants start here 
D= 0.06  ### Diameter (m)
ROL= 997.9 #### Density liquid  Kg/m3
ROG= 1.2  ##### density gas    kg/m3
UL= 1.1*(10**(-3))### viscosit Pa.S
UG= 0.000018  ### viscosit Pa.S
G= 9.81   # gravitational acceleration  m/s2
sigma= 0.06    #### Interfacial tension  N/m
alpha= 20   #### inclination from horizontal   Degrees
SL=3.14*D
### constants ends

##### list of the values of the variable that we are solving the equation with each time
list_of_liquid_velocities_VSL= np.linspace(0.001,1,80) #va
### ends here

### lists to store the output values in
list_of_F_Function=[]
list_of_interfacial_tension=[]
list_of_interfacial_tension_gas=[]  
list_of_x_solutions=[]
list_of_VSG=[]   ###### main values we are looking for
rooots=[]
### ends

for VSL in list_of_liquid_velocities_VSL: 
    def fitness_roots(VSL):
        if (D*ROL*VSL/(UL)) < 2300:
            n = 1
            CL = 16
        else:
            n = 0.2
            CL = 0.046

        a1 = G*(ROL-ROG)*D*sin(alpha*3.14/180)
        a2 = -(CL/16)*ROL*(D*ROL/UL)**(-n)*VSL**(2-n)
            #result = a1*((1-2*x)**2-2*(x-x**2)) + a2*((x-x**2)+(1-2*x)**2)/(x-x**2)**3
            #coeffs =[ a2, -3*a2, 3*a2, a1, -9*a1, 27*a1, -37*a1, 24*a1, -6*a1] # from mathematica cause I am lazy
        coeffs=[-6*a1, +24*a1, -37*a1, 27*a1, -9*a1, a1, 3*a2, -3*a2, a2]
        p = Polynomial(coeffs)
        r = p.roots()
            # return r # if you want all roots
        return r [abs(r.imag) < 1e-15].real # if you only want real roots
    
         
    rooots.append(fitness_roots(VSL))
    x=(rooots[0])
    AL=3.14*(D**2)*(x-x**2)
    
    DL=4*AL/SL
    A=3.14*(0.5*D)**2
    VL=VSL*A/AL
    if (D*ROL*VSL/(UL))<2300:
        n=1
        CL=16
        m=0.2
        CG=0.046               
    else:
        n=0.2
        CL=0.046
        m=0.2
        CG=0.046
                
    F = ((G*(ROL-ROG)*D*np.sin(alpha*3.14/180)*((1-2*x)**2-2*(x-x**2)))-((CL/16)*ROL*((D*ROL/UL)**(-n))*(VSL**(2-n))*(((x-x**2)+(1-2*x)**2)/(x-x**2)**3))) ### should be equal to 0
    list_of_F_Function.append((F))
        
    interfacial_tenstion=((G*(ROL-ROG)*D*np.sin(alpha*3.14/180)*(x-x**2)*(1-2*x))+((CL/32)*ROL*(D*ROL/UL)**(-n)*VSL**(2-n)*((1-2*x)/(x-x**2)**2)))  ### should be equal to interfacial_tension_gas
        
    VSG  =  (interfacial_tenstion*((1-2*x)**4) /((0.5)*(CG)*((D*ROG/UG)**(-m))*(1+300*x)*(ROG)))**(1/(2-m))
        
    list_of_VSG.append(VSG)
    list_of_x_solutions.append(x)
    print('VSL: {VSL}\t {fitness_roots(VSL)}')

    # list of the values of the variable that we are solving the equation with each time

for VSL in list_of_liquid_velocities_VSL:
        print(f'VSL: {VSL}\t {fitness_roots(VSL)}')

plt.plot()
plt.xscale("log")
plt.yscale("log")
plt.xlim(0.01,100)
plt.ylim(0.0001,100)
plt.xlabel("VSG", fontsize=9)
plt.ylabel("VSL",fontsize=9)
plt.grid(True, which="both", ls=":", color='0.002', linewidth=0.4)

plt.plot (list_of_VSG,list_of_liquid_velocities_VSL, label='VBarnea', color='blue')    

or this code

```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve, root, minimize
import math

from scipy.interpolate import make_interp_spline, BSpline

###### Constants start here 
D= 0.06  ### Diameter (m)
ROL= 997.9 #### Density liquid  Kg/m3
ROG= 1.2  ##### density gas    kg/m3
UL= 1.1*(10**(-3))### viscosit Pa.S
UG= 0.000018  ### viscosit Pa.S
G= 9.81   # gravitational acceleration  m/s2
sigma= 0.06    #### Interfacial tension  N/m
alpha= 20   #### inclination from horizontal   Degrees
SL=3.14*D
### constants ends

##### list of the values of the variable that we are solving the equation with each time
list_of_liquid_velocities_VSL= np.linspace(0.001,1,80) #va
### ends here

### lists to store the output values in
list_of_F_Function=[]
list_of_interfacial_tension=[]
list_of_interfacial_tension_gas=[]  
list_of_x_solutions=[]
list_of_VSG=[]   ###### main values we are looking for
### ends

VSLexperiment= [0.02,0.06,0.1]
VSGexperiment=[13.94,14.54,19.09]


#################  Barnea Model  ##################


for VSL in list_of_liquid_velocities_VSL:  
    def myF(z):
        x=z[0]
        if (D*ROL*VSL/(UL)) < 2300:
            n = 1
            CL = 16
        else:
            n = 0.2
            CL = 0.046

        a1 = G*(ROL-ROG)*D*np.sin(alpha*3.14/180)       
        a2 = -(CL/16)*ROL*(D*ROL/UL)**(-n)*VSL**(2-n)
            #result = a1*((1-2*x)**2-2*(x-x**2)) + a2*((x-x**2)+(1-2*x)**2)/(x-x**2)**3
            #coeffs =[ a2, -3*a2, 3*a2, a1, -9*a1, 27*a1, -37*a1, 24*a1, -6*a1] # from mathematica cause I am lazy
        F =-6*a1*x**8 +24*a1*x**7 -37*a1*x**6 + 27*a1*x**5 -9*a1*x**7 + a1*x**3+ 3*a2*x**2 -3*a2*x + a2
        return F
    
   
    solution = root(myF, [0.05], method='krylov')  
    x=solution.x   
    AL=3.14*(D**2)*(x-x**2)   
    DL=4*AL/SL
    A=3.14*(0.5*D)**2
    VL=VSL*A/AL
   
    if (D*ROL*VSL/(UL))<2300:
        n=1
        CL=16
        m=0.2
        CG=0.046       
    else:
        n=0.2
        CL=0.046
        m=0.2
        CG=0.046

    
    list_of_x_solutions.append(x)
    
    
    AL=3.14*(D**2)*(x-x**2)
    SL=3.14*D
    DL=4*AL/SL
    #DG=D-a
    #AG=3.14*(0.5*D-a)**2
    A=3.14*(0.5*D)**2
    VL=VSL*A/AL
    #VG=VSG*A/AG
    
    F= ((G*(ROL-ROG)*D*np.sin(alpha*3.14/180)*((1-2*x)**2-2*(x-x**2)))-((CL/16)*ROL*((D*ROL/UL)**(-n))*(VSL**(2-n))*(((x-x**2)+(1-2*x)**2)/(x-x**2)**3)))
    list_of_F_Function.append(F)
    
    if (DL*ROL*VL/(UL))<2300:      
        n=1
        CL=16
        m=0.2
        CG=0.046        
    else:              
        n=0.2
        CL=0.046
        m=0.2
        CG=0.046

for VSL,x in zip(list_of_liquid_velocities_VSL,list_of_x_solutions):
   
    interfacial_tension= ((G*(ROL-ROG)*D*math.sin(alpha*3.14/180)*(x-x**2)*(1-2*x))+((CL/32)*ROL*((D*ROL/UL)**(-n))*(VSL**(2-n))*((1-2*x)/(x-x**2)**2)))

         
    VSG= (interfacial_tension*((1-2*x)**4) /((0.5)*(CG)*((D*ROG/UG)**(-m))*(1+300*x)*(ROG)))**(1/(2-m))
    interfacial_tension_gas= ((((VSG)**(2-m))*(0.5)*(CG)*((D*ROG/UG)**(-m))*(1+300*x)*(ROG))/((1-2*x)**4))
      
    list_of_interfacial_tension_gas.append(interfacial_tension_gas)
    list_of_VSG.append(VSG)  
    list_of_interfacial_tension.append(interfacial_tension)
    
plt.plot()
plt.xscale("log")
plt.yscale("log")
plt.xlim(1,100)
plt.ylim(0.001,1)
plt.xlabel("VSG", fontsize=9)
plt.ylabel("VL",fontsize=9)
plt.grid(True, which="both", ls=":", color='0.002', linewidth=0.4)

plt.plot (list_of_VSG,list_of_liquid_velocities_VSL, label='VBarnea')
plt.scatter(VSGexperiment,VSLexperiment, color="Black", linewidth="0.05", label="VsgExperiment")
$\endgroup$
9
  • $\begingroup$ Can you plot the function on the left for different values of VSL? It should have values where it cross the x-axis (well, the $\delta_l$ axis) and it should be relatively easy to see what the problem may be. $\endgroup$ Nov 5, 2021 at 2:27
  • $\begingroup$ @WolfgangBangerth I edited the post probably now is clear enough. $\endgroup$ Nov 5, 2021 at 3:19
  • $\begingroup$ The edit didn't make it better. I'm now totally confused which equation is the problem. Please make it clear which equation you want to, but can't solve, and what all of the variables in that are. For example, in both of the equations you currently have in the question, I have no idea what variables have fixed values and which ones are free that you try to find values for that solve the equation. $\endgroup$ Nov 5, 2021 at 4:08
  • $\begingroup$ the first equation is the one I am trying to solve by getting the value of delta_l that makes the equation equal to 0. In that equation I have constants such as (g , D, ROL, ROG), each time I will vary VSL and look for delta_l values that make the equation equal to 0 (I will carry VSL 80 times with this np.linspace(0.001,1,80) to obtain 80 values of delta_l that makes the equation equal to 0. The second equation gets that value of delta_l to solve for VSG with the other parameters fixed. the problem is in the first equation which is defined in the code as (def fitness_func...... derivative). $\endgroup$ Nov 5, 2021 at 14:42
  • 1
    $\begingroup$ I tried to reformat your question, but there is a lot more that could be done to clarify your post. For one, you can convert the images of equations into Mathjax, as I have done for several variables you wrote in your question. The code is also very hard read; you could try to change the formatting (reasonable number of spaces between lines, define intermediate quantities rather than massive multiline variables). $\endgroup$
    – Tyberius
    Nov 16, 2021 at 20:35

1 Answer 1

2
$\begingroup$

Your equation essentially reads with coefficients $a_1$ and $a_2(v_{Sl})$ $$a_1 ((1-2x)^2-2(x-x^2))+a_2 \frac{x-x^2+(1-2x)^2}{(x-x^2)^3} =0$$ After multiplying both sides with the denominator $(x-x^2)^3$ you obtain a polynomial of degree 8. Its roots can be easily found by calculating the eigenvalues of its companion matrix, with e.g. numpy.polynomials.roots.

So there is no need for any fancy genetic algorithm.

This is done in the attached code. I didn't get from your question if you have any more criteria to select from multiple roots or what you do if there are no roots.

from numpy import sin, linspace
from numpy.polynomial import Polynomial
# Constants
D = 0.06  # Diameter (m)
ROL = 997.9  # Density liquid  Kg/m3
ROG = 1.2  # density gas    kg/m3
UL = 1.1*(10**(-3))  # viscosit Pa.S
UG = 0.000018  # viscosit Pa.S
G = 9.81   # gravitational acceleration  m/s2
sigma = 0.06  # Interfacial tension  N/m
alpha = 20  # inclination from horizontal   Degrees
SL = 3.14*D



def fitness_roots(VSL):
    if (D*ROL*VSL/(UL)) < 2300:
        n = 1
        CL = 16

    else:
        n = 0.2
        CL = 0.046

    a1 = G*(ROL-ROG)*D*sin(alpha*3.14/180)
    a2 = -(CL/16)*ROL*(D*ROL/UL)**(-n)*VSL**(2-n)
    #result = a1*((1-2*x)**2-2*(x-x**2)) + a2*((x-x**2)+(1-2*x)**2)/(x-x**2)**3
    coeffs =[ a2, -3*a2, 3*a2, a1, -9*a1, 27*a1, -37*a1, 24*a1, -6*a1] # from mathematica cause I am lazy

    p = Polynomial(coeffs)
    r = p.roots()
    # return r # if you want all roots
    return r[abs(r.imag) < 1e-15].real # if you only want real roots


if __name__ == '__main__':

    VSLexperiment = [0.02, 0.06, 0.1]
    VSGexperiment = [13.94, 14.54, 19.09]

    for VSL in VSLexperiment:
        print(f'VSL: {VSL}\t {fitness_roots(VSL)}')

    # list of the values of the variable that we are solving the equation with each time
    list_of_liquid_velocities_VSL = linspace(0.001, 1, 80)  # va

    for VSL in list_of_liquid_velocities_VSL:
        print(f'VSL: {VSL}\t {fitness_roots(VSL)}')
$\endgroup$
2
  • $\begingroup$ I did what you suggested and here is the code. However, when I recalculate the function with the roots it doesn't give me 0 at all. the code was added to the post. $\endgroup$ Nov 6, 2021 at 4:50
  • $\begingroup$ I don't get, why you change the order of coefficients. Of course the now calculated roots are different. $\endgroup$
    – Bort
    Nov 6, 2021 at 8:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.