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.
The following graph is a reference which I need to reproduce the same graph with my code.
Here is the graph that I generated.
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:
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")