Ciao a tutti,
Innanzitutto chiarisco che utilizzo Anaconda ed in particolare l'interfaccia Spyder .


Per un progetto universitario avevo intenzione di modellare la camera di pompaggio di un compressore, descrivibile matematicamente come un sistema di ODE nel dominio del tempo. Una volta scritte in Python le equazioni pero mi accorgo che qualcosa non mi funziona nello script; in particolare noto che provando a stampare durante l'integrazione il tempo questo non prosegue e rimane fisso ad un dato valore. In altri termini sembra che non prosegua con gli intervalli temporali successivi e che non segua il "linespace" che ho inserito. Trascurando possibili errori nelle formule (che posso ricontrollare autonomamente) se provaste a dare un occhio allo script (e sulla metodologia di integrazione con odeint) dandomi una vostra opinione sul malfunzionamento mi sarebbe davvero d'aiuto . A seguire lo script:






import numpy as np




# y sarà la variabile lista formata da [Temperatura(t),Massa(t)]




def Compressor(y,t): #INPUT:t=tempo s ,

#T=temperatura al tempo t-h con h=passo d'integrazione K ,

#m=massa al tempo t-h con h=passo d'integrazione kg




T , m = y #Assegno i simboli T e m rispettivamente al primo e secondo valore di y



D=55. #Alesaggio [mm^2]

N=1500. #Velocità angolare rpm

r=84.845 #lunguezza biella mm

d=0.2817 #Rapporto lunghezza manovella / Lunghezza biella dimensionless

Ra=287.05*10**(6) #Costante aria [mm^2/(s^2 * K)]

Vc=5678.23 #Volume morto cilindro [mm^3]

T_amb=15 + 273.15 #Temperatura ambiente K

ro_amb=1.225*10**(-9) #Densità aria a 15°C e 1atm [kg/mm^3]

alpha=0.9 #Flow coefficient dimensionless

c=0.9 #Coeff. di valvola dimensionless

Cg=0.45 #Coeff. di superficie utile valvola dimensionless

kg=0.028*10**(3) #Coefficiente di conducibilità termica dell’aria [(kg*mm)/(s^3 *K)]

Tw=78. + 273.15 #Temperatura cilindro a regime ipotizzata costante K

a=0.7 #Coeff. per coefficiente convettivo

b=0.7 #Coeff. per coefficiente convettivo

qlimit=2.5 #Distanza massima imposta dal limitatore mm

kmol=1.3*10**(-3) #Coeff. elastico N/mm

Av=40*10 #Area valvola [mm^2]

ps=1.01325*10**(-1) #Pressione suction chamberMPa

pd=10. *10**(-1) #Pressione discharge chamberMPa

Ts=T_amb #Temperatura suction chamber K

Td=570.5 #Temperatura discharge chamberK

Cvs=0.7187529739145476 #Calore specifico volume cost. suction chamber

Cvd=0.7580872067793204 #Calore specifico volume cost. discharge chamber



ro_s=ro_amb



def x(t): #Spostamento del pistone

teta=((2*np.pi*N)/60)*t

return r*(1-np.cos(teta)+(d)**(-1)*(1-np.sqrt(1-(d*np.sin(teta))**2)))




def xdot(t): #Velocità del pistone

w=((2*np.pi*N)/60)

teta=w*t

return(r*(w*np.sin(teta)+((d*w*np.sin(teta)*np.cos(teta))/

(np.sqrt(1-(d*np.sin(teta))**2)))))



def V(t): #Volume

return(Vc+((np.pi*D**2)/4)*x(t))



def Vdot(t): #Derivata temporale del volume

return((np.pi*D**2)/4)*xdot(t)



def ro(m,t): #Densità nella camera

return(m/V(t))



def p(T,m,t): #Calcolo la pressione con la legge dei gas perfetti

return(ro(m,t)*Ra*T)



def Cp(T): #calore speciico a pressione costante con relazione polinomiale

a0=0.1034090*10

a1=-0.2848870*10**(-3)

a2=0.7816818*10**(-6)

a3=-0.4970786*10**(-9)

a4=0.1077024*10**(-12)



return(a0 + a1*T + a2*T**2 + a3*T**3 + a4*T**4)



def Cv(T): #calore specifico a valume costante con relazione di Mayer

return(Cp(T)-Ra)



def k(T):

return(Cp(T)/Cv(T))



def dmvdt(pu,pds,ro_u): #pds=pressione downstream, pu=pressione upstream

#portata massica attraversante la generica valvola

q=(1/kmol)*(Cg*Av*(pu-pds))



if q>qlimit:

q=qlimit

A0=2*45*q

eps=1-(c/1.41)*((pu-pds)/pu)

return(alpha*eps*A0*(2*ro_u*(pu-pds))**(1/2))



elif q<0:

return(0.)



else:

A0=2*45*q

eps=1-(c/1.41)*((pu-pds)/pu)

return(alpha*eps*A0*(2*ro_u*(pu-pds))**(1/2))



def dmdt(T,m,t):

#derivata temporale della massa nel volume di controllo

dmsdt=dmvdt(ps,p(T,m,t),ro_s)

dmddt=dmvdt(p(T,m,t),pd,ro(m,t))

return(dmsdt-dmddt)



def dQdt(T,t):

#Derivata temporale del calore scambiato con il cilindro

A=((np.pi*D**2)/4)+(4*Vc/D)+(np.pi*D*x(t))

up=np.abs(dmdt(T,m,t)*(4/(ro(m,t)*np.pi*D**2)))

h=a*(kg/D)*((ro(m,t)*up*D*Cp(T))/kg)**b

return(h*A*(Tw-T))




def dTdt(T,m,t):

return((1/(m*Cv(T))) * (k(T)*Cvs*Ts*dmvdt(ps,p(T,m,t),ro_s) -

k(T)*Cvd*Td*dmvdt(p(T,m,t),pd,ro(m,t)) +

dQdt(T,t) - T*(ro(m,t)*Ra*Vdot(t) + Cv(T)*dmdt(T,m,t))))




# p.append(p0)



dydt=[dTdt(T,m,t),dmdt(T,m,t)]

return dydt




y0 = [570.5,3.4673667*10**(-] #[Temperatura al tempo 0, Densità al tempo 0]

t = np.linspace(0, 40*10**(-3), 1000)




from scipy.integrate import odeint

sol = odeint(Compressor, y0, t)







import matplotlib.pyplot as plt

plt.plot(t, sol[:, 0], 'b', label='T(t)')

plt.legend(loc='best')

plt.xlabel('t')

plt.grid()

plt.show()











Grazie in anticipo per l'aiuto! Se serve direttamente il file ditemi che lo allego

Dato che ho visto che il copia e incolla ha avuto scarsi risultati posto direttamente il file xD