Closed ekremekc closed 1 year ago
andfoy
commented
4 years ago Hi @ekremekc, thanks for opening this issue, however, we cannot do much at this time, as we don't have the script you are trying to run, could you please post it here in order to reproduce your issue?
ekremekc
commented
4 years ago The following script has a problem; For example, script running perfectly in the first place. Then change the T1 from 470 to 600, it will freeze;
from compressibleflow import *
T1 = 470
P1 = 500
v1 = 180
A1 = 0.05
air = Air(T1,P1)
M1 = air.mach_number(v1)
T0 = air.stagnation_temp(M1)
P0 = air.stagnation_pressure(M1)
A_star = air.throat_area(A1, air.mach_number(v1))
m_star = air.m_dot(v1,air.diameter(A1))
area = 0.04
ratio = air.throat_area_ratio(area,A1, M1)
Ma_upward = air.ma_finder("upward",area, ratio,True)
Ma_downward = air.ma_finder("downward",area, ratio, True)
T_upward = air.temperature(Ma_upward,T0)
T_downward = air.temperature(Ma_downward,T0)
P_upward = air.pressure(Ma_upward,P0)
P_downward = air.pressure(Ma_downward,P0)
exit_p = air.exit_pressure(M1)
air.plot(0.5*A1, A1,M1,'Ma','T',1000,method='bisection')
Compressible flow library is shown as below;
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
class Gas():
def __init__(self,T,P,R_u=8.31447):
self.T = T
self.P = P
self.R_u=R_u
def gas_list():
print(" Code\t","Gas","\n",
"----\t","---","\n",
"Air\t","Air" ,"\n",
"CO2\t","Carbon dioxide" ,"\n",
"N2\t\t","Nitrogen" ,"\n"
)
def area(self, diameter):
return (np.pi*(diameter**2))/4
def critical_area(self,massflowrate):
return massflowrate/(self.P*1000*(self.k**(1/2))*(2/(self.k+1))**((self.k+1)/(2*self.k-2))/((self.R*1000*self.T)**(1/2)))
def critical_m_dot(self, diameter=1):
return self.critical_density()*self.area(diameter)*self.critical_speed_of_sound()
def critical_temperature(self):
return self.T*2/(self.k+1)
def critical_pressure(self):
return self.P*(2/(self.k+1))**(self.k/(self.k-1))
def critical_density(self):
return self.rho*(2/(self.k+1))**(1/(self.k-1))
def critical_speed_of_sound(self):
return np.sqrt(self.k*self.R*self.critical_temperature()*1000)
def density(self):
return self.P/(self.R*self.T)
def diameter(self, area):
return np.sqrt(4/np.pi*area)
def exit_temperature(self,Mach):
return self.T/(1+(self.k-1)/2*Mach**2)
def exit_pressure(self,Mach):
return self.P/(1+(self.k-1)/2*Mach**2)**(self.k/(self.k-1))
def exit_density(self, Mach):
return self.rho/(1+(self.k-1)/2*Mach**2)**(1/(self.k-1))
def exit_speed(self, Mach):
return Mach*np.sqrt(self.k*self.R*self.exit_temperature(Mach)*1000)
def exit_area(self, Throat_Area, Mach):
return Throat_Area*(1/Mach)*((2/(self.k+1))*(1+(self.k-1)/2*Mach**2))**((self.k+1)/(2*self.k-2))
def mach_number(self, velocity):
return velocity/self.speed_of_sound()
def m_dot(self, velocity, diameter=1):
return self.density()*self.area(diameter)*velocity
def mfr(self,velocity, diameter):
return self.critical_pressure()*self.area(diameter)*self.mach_number(velocity)*np.sqrt(self.k/(self.R*self.critical_temperature()))
def mass_flowrate(self, velocity, diameter=1):
return (self.area(diameter)*self.mach_number(velocity)*self.stagnation_pressure(velocity)*np.sqrt(self.k/(self.R*self.stagnation_temp(velocity))))\
/((1+(self.k-1)*(self.mach_number(velocity)**2)/2)**((self.k+1)/(2*(self.k-1))))
def ma_finder(self, section, area, area_ratio, show_iterations=False, tolerance=10e-6, method="bisection"):
try:
if section !="upward" and section !="downward":
raise NameError
def finder(Ma):
value = (1/Ma*((1+0.5*(self.k-1)*Ma**2)/(0.5*(self.k+1)))**(0.5*(self.k+1)/(self.k-1)))
if method=='golden' or method=='secant':
target = abs(value - area_ratio)
elif method=='bisection':
target = value - area_ratio
return target
# def check_boundaries(Ma_0, Ma_1):
# if section=="upward":
# if Ma_0>1 or Ma_1>1:
# Ma_0 = 1/Ma_0
# Ma_1 = Ma_0+0.001
# # print("ma kucuk 1 den calisti")
# elif section=="downward":
# if Ma_0<1 or Ma_1<1:
# Ma_0 = 1+Ma_0
# Ma_1 = Ma_0+0.1
# # print("ma buyuk 1 den calisti")
if section=="upward":
if method=='bisection':
Ma=bisection_method( finder,0, 1, tolerance = 10e-6,show_iterations=show_iterations)
elif method=='secant':
Ma=secant_method( finder,0, 1, tolerance = 10e-6,show_iterations=show_iterations)
elif method=='golden':
Ma=golden_section(finder,0, 1, tolerance = 10e-6,show_iterations=show_iterations)
elif section=="downward":
if method=='bisection':
Ma=bisection_method( finder,1, 5, tolerance = 10e-6,show_iterations=show_iterations)
elif method=='secant':
Ma=secant_method( finder,1, 5, tolerance = 10e-6,show_iterations=show_iterations)
elif method=='golden':
Ma=golden_section(finder,1, 5, tolerance = 10e-6,show_iterations=show_iterations)
return Ma
except NameError:
raise NameError("Please specify the flow by using these keywords: \"upward\" or \"downward\"") from None
except ValueError:
raise ValueError("Given area is smaller than throat area. Program has terminated.\n Hint: You could change the division number.") from None
def plot(self,area_start, area_end, Mach_start, y_axis='T', color_bar='Ma', division=250 ,x_axis='A', method="bisection"):
area_upward = np.linspace(area_start, self.throat_area(area_start,Mach_start), division)
area_downward = np.linspace(self.throat_area(area_start,Mach_start), area_end, division)
area_total = np.concatenate((area_upward,area_downward))
ST = self.stagnation_temp(Mach_start)
temp_upward = []
Ma_upward = []
for i in range(division):
ratio = self.throat_area_ratio(area_upward[i], area_start, Mach_start)
Ma=self.ma_finder("upward",area_upward[i],ratio,method=method)
Ma_upward.append(Ma)
temp_upward.append(self.temperature(Ma, ST))
temp_downward = []
Ma_downward = []
for i in range(division):
ratio = self.throat_area_ratio(area_downward[i], area_start, Mach_start)
Ma=self.ma_finder("downward",area_downward[i],ratio,method=method)
Ma_downward.append(Ma)
temp_downward.append(self.temperature(Ma, ST))
temp_total = temp_upward +temp_downward
Ma_total = Ma_upward +Ma_downward
fig = plt.figure(figsize=(10,7.5))
ax = fig.add_subplot(111)
xs = np.linspace(0,1,2*division)
if y_axis == 'T':
y_lbl='Temperature (K)'
if color_bar=='Ma':
color = Ma_total
mp = ax.scatter((xs),(temp_total),c=color,cmap=plt.cm.get_cmap('jet'))
c_lbl = 'Mach Number'
elif color_bar=='T':
mp = ax.scatter((xs),(temp_total),c=temp_total,cmap=plt.cm.get_cmap('jet'))
c_lbl = 'T (K)'
elif y_axis == 'Ma':
y_lbl='Mach Number'
if color_bar=='Ma':
color = Ma_total
mp = ax.scatter((xs),(Ma_total),c=color,cmap=plt.cm.get_cmap('jet'))
c_lbl = 'Mach Number'
elif color_bar=='T':
mp = ax.scatter((xs),(Ma_total),c=temp_total,cmap=plt.cm.get_cmap('jet'))
c_lbl = 'T (K)'
cb = plt.colorbar(mp)
cb.set_label(c_lbl)
ax.set(title=r'Converging- Diverging Nozzle',
xlabel='Area $m^2$', ylabel=y_lbl)
tick_labels=[]
for j in np.linspace(0,(2*division),7):
if j==2*division:
tick_labels.append(round(area_total[-1],4))
else:
tick_labels.append(round(area_total[int(j)],4))
plt.xticks(np.linspace(0,1,7),tick_labels)
plt.show()
def pressure(self, Mach, Stagnation_Pressure):
return Stagnation_Pressure/((1+0.5*(self.k-1)*Mach**2)**(self.k/(self.k-1)))
def speed_of_sound(self):
return np.sqrt(self.k*self.R*self.T*1000)
def stagnation_temp(self,Mach):
return self.T*(1+(self.k-1)/2*Mach**2)
def stagnation_pressure(self,Mach):
return self.P*(1+0.5*(self.k-1)*Mach**2)**(self.k/(self.k-1))
def temperature(self, Mach, Stagnation_Temperature):
return Stagnation_Temperature/(1+(self.k-1)/2*Mach**2)
def throat_area(self,known_area,Mach):
return known_area/((1/Mach)*((2/(self.k+1))*(1+(self.k-1)/2*Mach**2))**((self.k+1)/(2*self.k-2)))
def throat_area_ratio(self,wanted_area, known_area,known_Mach):
return wanted_area/self.throat_area(known_area, known_Mach)
class Air(Gas):
def __init__(self,T,P=101.325):
Gas.__init__(self, T, P)
self.M=28.97
self.k=1.4
self.R=self.R_u/self.M
self.cp=1.9327E-10*self.T**4 - 7.9999E-07*self.T**3 + 1.1407E-03*self.T**2 - 4.4890E-01*self.T + 1.0575E+03
self.rho = self.P/(self.R*self.T)
class CO2(Gas):
def __init__(self,T,P=101.325):
Gas.__init__(self, T, P)
self.M=44.01
self.k=1.289
self.R=self.R_u/self.M
self.cp=0.849
self.rho = self.P/(self.R*self.T)
class N2(Gas):
def __init__(self,T,P=101.325):
Gas.__init__(self, T, P)
self.M=28.01
self.k=1.4
self.R=self.R_u/self.M
self.cp=1.040
self.rho = self.P/(self.R*self.T)
##### NUMERICAL METHODS
def golden_section(func,starting, ending, show_iterations=False, tolerance = 10e-6):
gr=(np.sqrt(5)+1)/2-1
dm=tolerance
a0 = starting+dm
b0 = ending-dm
count=0
while True:
count+=1
# print(finder(Ma_0))
# print(finder(Ma_1))
a1=(b0-a0)*(1-gr)+a0
b1=(b0-a0)*gr+a0
if abs(func((a1+b1)/2))<=tolerance:
if 1>=ending:
print("The Mach number below unity is: ",a1,"\n")
elif starting>=1:
print("The Mach number above unity is: ",a1,"\n")
break
else:
if func(a0)>func(b0):
a0=a1
b1=b1
else:
a0=a0
b0=b1
if show_iterations ==True:
print("Iteration ", count, " :",a1)
return (a1+b1)/2
def secant_method(func, lower_bound, upper_bound, show_iterations=False,tolerance=10e-6):
Ma_0 = (upper_bound+lower_bound)/2
dMa = 0.01
Ma_1 = Ma_0+dMa
count=0
while True:
count+=1
Ma_2 = Ma_1 - func(Ma_1)*(Ma_1-Ma_0)/(func(Ma_1)-func(Ma_0))
if show_iterations ==True:
print("Iteration ", count, " :",Ma_2)
if func(Ma_2)<=tolerance:
if show_iterations ==True:
print("The Mach number below unity is: ",Ma_2,"\n")
break
else:
Ma_0 = Ma_1
Ma_1 = Ma_2
return Ma_2
def bisection_method(func, lower_bound, upper_bound, show_iterations=False,tolerance=10e-6):
if lower_bound==0 :
lower_bound+=tolerance
a=lower_bound
b= upper_bound
count = 0
while True:
count+=1
c = (a+b)/2
if abs(func(c))<=tolerance:
if show_iterations ==True:
print("The Mach number is: ",c,"\n")
break
else:
if func(a)*func(c)>func(b)*func(c):
b=b
a=c
else:
a=a
b=c
if show_iterations ==True:
print("Iteration ", count, " :",c)
return c
I have uploaded the gif.
@ccordoba12
ccordoba12
commented
1 year ago Sorry for not answering in time but I hope the problem is fixed for you in our latest version.
Issue Report Checklist
conda update spyder(orpip, if not using Anaconda)jupyter qtconsole(if console-related)spyder --resetProblem Description
First execution of the script works perfectly. When I change the some parameters in the script, the Spyder IDE freezes till execute all the script. But if I re-run immediately it works fine. Then again, if some parameter changed (e.g. from a=1 to a=2) it freezes about 3 4 seconds again.
What steps reproduce the problem?
Versions
Dependencies
Mandatory:
atomicwrites >=1.2.0 : 1.4.0 (OK) chardet >=2.0.0 : 3.0.4 (OK) cloudpickle >=0.5.0 : 1.4.1 (OK) diff_match_patch >=20181111 : 20181111 (OK) intervaltree : None (OK) IPython >=4.0 : 7.13.0 (OK) jedi =0.15.2 : 0.15.2 (OK) keyring : None (OK) nbconvert >=4.0 : 5.6.1 (OK) numpydoc >=0.6.0 : 0.9.2 (OK) parso =0.5.2 : 0.5.2 (OK) pexpect >=4.4.0 : 4.8.0 (OK) pickleshare >=0.4 : 0.7.5 (OK) psutil >=5.3 : 5.7.0 (OK) pygments >=2.0 : 2.6.1 (OK) pylint >=0.25 : 2.5.0 (OK) pyls >=0.31.9;<0.32.0 : 0.31.10 (OK) qdarkstyle >=2.8 : 2.8.1 (OK) qtawesome >=0.5.7 : 0.7.0 (OK) qtconsole >=4.6.0 : 4.7.4 (OK) qtpy >=1.5.0 : 1.9.0 (OK) rtree >=0.8.3 : 0.9.4 (OK) sphinx >=0.6.6 : 3.0.3 (OK) spyder_kernels >=1.9.1;<1.10.0 : 1.9.1 (OK) watchdog : None (OK) xdg >=0.26 : 0.26 (OK) zmq >=17 : 18.1.1 (OK)
Optional:
cython >=0.21 : 0.29.17 (OK) matplotlib >=2.0.0 : 3.1.3 (OK) numpy >=1.7 : 1.18.1 (OK) pandas >=0.13.1 : 1.0.3 (OK) scipy >=0.17.0 : 1.4.1 (OK) sympy >=0.7.3 : 1.5.1 (OK)