Dev elias 2026 04 29

Compressor_model_CP.py

@@ -0,0 +1,390 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jan 31 17:46:25 2019
+
+@author: roskosch
+"""
+
+from pylab import *
+import Fluid_CP as fl
+
+
+Rm=8.3145 #Gaskonstante J/mol/K
+Tu=25. #Umgebungstemperatur
+Ver0=[34e-3,34e-3,2.,.04]# Fit-Verdichter: D, H, Zylinder, Außenfläche
+
+
+###########################################
+#  Verdichterspezifische Parameter
+# pV:    0 - D, Kolbendurchmesser, m
+#        1 - H, Hub, m
+#        2   Verhältnis aus Kurbelwellenradius und Pleillänge
+#        3   A_a, Oberfläche Verdichter, m²
+#        4 - c1, Totvolumen relativ zum Hubvolumen
+#        5 - preib, Reibdruck, kPa
+#        6 - f, elektrische Verdichterfrequenz Hz
+#        7 - Verdichterdrehazhal Hz
+#        8 - Zylinderanzahl 
+
+
+
+# pZ: 
+# Index S: Saugseite
+# Index D: Druckseite
+#       0 - TS, °C
+#       1 - pS, kPa
+#       2 - vS, m³/kg
+#       3 - uS, kJ/kg
+#       4 - hS, kJ/kg
+#       5 - sS, kJ/kg/K
+#       6 - pD, kPa
+
+#pZyk: Größen, die während eines Zyklusses konstant bleiben bzw. über die
+#      Zyklen iteriert werden
+#       0 - Aeff_e, effektiver Ströumgsquerschnitt Eintritt, m²
+#       1 - Aeff_a, effektiver Strömungsquerschnitt Austritt, m²
+
+
+# z_it: Nimmt für jeden Iterationsschritt folgende Werte auf
+#       0 - tet, Kurbelwinkel
+#       1 - x, Kolbenposition
+#       2 - V, Zylindervolumen
+#       3 - Fläche Wärmeübertragung Zylinderwand 
+#       4 - Zyklusschritt, 0=Verdichten, 1=Ausschieben, 2=Expandieren, 3=Ansaugen
+#       5 - T, Temperatur im Zylinder °C
+#       6 - p, Druck im Zylinder kPa
+#       7 - v, spezifisches Volumen im Zylinder m³/kg
+#       8 - u, innere Energie im Zylinder kJ/kg
+#       9 - h, Enthalpie im Zylinder kJ/kg
+#      10 - s, Entropie im Zylinder kJ/kg/K
+#      11 - m, Masse im Zylinder kg
+#      12 - T_th, Temperatur der thermischen Masse, °C
+#      13 - alpha, Wärmeübergangskoeffizient innen
+#      14 - dm, Eingeströmte bzw. ausgeströmte Masse im Iterationsschritt kg
+#      15 - Übertragene Wärme
+
+
+
+
+
+
+IS=5000 # Anzahl der differentiellen Schritte für einen Zyklus
+pV=zeros(8,float)
+pZ=zeros(7,float)
+pZyk=zeros(2,float)
+z_it=zeros([IS,16])
+fluid=[]
+
+############################################################################################################################################################
+def init(p_e_,fluid_in):
+    global v_fit, u_fit, T_fit, p_fit, h_fit, initialisiert
+    er1=fl.get_fluid_info(fluid_in)
+    Tc=er1["T_crit"]
+    Tmin=er1["T_min"]
+    Tmax=er1["T_max"]
+
+    T_e=fl.state(['p','x'],[p_e_,1],fluid_in,Eh="CKPa")["T"]-5.
+
+    if T_e+200.>=Tmax:
+        T_a=Tmax-5.
+    else:
+        T_a=T_e+200.
+
+    if (T_e-20.)>=Tmin:
+        T_q0_min=T_e-20.
+    else:
+        T_q0_min=Tmin
+
+    T_q0=linspace(T_q0_min,Tc-10.,100)
+    v_fit=zeros(100,float)
+    u_fit=zeros([20,100],float)
+    T_fit=zeros([20,100],float)
+    p_fit=zeros([20,100],float)
+    h_fit=zeros([20,100],float)
+    for i1 in range(100):
+        v_fit[i1]=fl.state(['T','x'],[T_q0[i1],1],fluid_in,Eh="CKPa")["v"]
+        TT=linspace(T_q0[i1],T_a,20)
+        for i2 in range(20):
+            T_fit[i2,i1]=TT[i2]
+
+            erg=fl.state(['T','v'],[TT[i2],v_fit[i1]],fluid_in,Eh="CKPa")
+            u_fit[i2,i1]=erg["u"]
+            p_fit[i2,i1]=erg["p"]
+            h_fit[i2,i1]=erg["h"]
+
+    return
+
+
+def Stoffdaten(u_su,v_su):
+
+    #print u_su, v_su
+    pos=where(v_fit>=v_su)[0][-1]
+    u_pos1=where(u_su>=u_fit[:,pos])[0][-1]
+    u_pos2=where(u_su>=u_fit[:,pos+1])[0][-1]
+
+    #### Temp
+
+    T1=(T_fit[u_pos1+1,pos]-T_fit[u_pos1,pos])/(u_fit[u_pos1+1,pos]-u_fit[u_pos1,pos])*u_su+T_fit[u_pos1+1,pos]-(T_fit[u_pos1+1,pos]-T_fit[u_pos1,pos])/(u_fit[u_pos1+1,pos]-u_fit[u_pos1,pos])*u_fit[u_pos1+1,pos]
+    T2=(T_fit[u_pos2+1,pos+1]-T_fit[u_pos2,pos+1])/(u_fit[u_pos2+1,pos+1]-u_fit[u_pos2,pos+1])*u_su+T_fit[u_pos2+1,pos+1]-(T_fit[u_pos2+1,pos+1]-T_fit[u_pos2,pos+1])/(u_fit[u_pos2+1,pos+1]-u_fit[u_pos2,pos+1])*u_fit[u_pos2+1,pos+1]    
+    T_pos=(T2-T1)/(v_fit[pos+1]-v_fit[pos])*v_su+T2-(T2-T1)/(v_fit[pos+1]-v_fit[pos])*v_fit[pos+1]
+
+
+    #### Druck
+    p1=(p_fit[u_pos1+1,pos]-p_fit[u_pos1,pos])/(u_fit[u_pos1+1,pos]-u_fit[u_pos1,pos])*u_su+p_fit[u_pos1+1,pos]-(p_fit[u_pos1+1,pos]-p_fit[u_pos1,pos])/(u_fit[u_pos1+1,pos]-u_fit[u_pos1,pos])*u_fit[u_pos1+1,pos]    
+    p2=(p_fit[u_pos2+1,pos+1]-p_fit[u_pos2,pos+1])/(u_fit[u_pos2+1,pos+1]-u_fit[u_pos2,pos+1])*u_su+p_fit[u_pos2+1,pos+1]-(p_fit[u_pos2+1,pos+1]-p_fit[u_pos2,pos+1])/(u_fit[u_pos2+1,pos+1]-u_fit[u_pos2,pos+1])*u_fit[u_pos2+1,pos+1]    
+    p_pos=(p2-p1)/(v_fit[pos+1]-v_fit[pos])*v_su+p2-(p2-p1)/(v_fit[pos+1]-v_fit[pos])*v_fit[pos+1]
+
+
+
+
+    #Enthalpie
+    h1=(h_fit[u_pos1+1,pos]-h_fit[u_pos1,pos])/(u_fit[u_pos1+1,pos]-u_fit[u_pos1,pos])*u_su+h_fit[u_pos1+1,pos]-(h_fit[u_pos1+1,pos]-h_fit[u_pos1,pos])/(u_fit[u_pos1+1,pos]-u_fit[u_pos1,pos])*u_fit[u_pos1+1,pos]    
+    h2=(h_fit[u_pos2+1,pos+1]-h_fit[u_pos2,pos+1])/(u_fit[u_pos2+1,pos+1]-u_fit[u_pos2,pos+1])*u_su+h_fit[u_pos2+1,pos+1]-(h_fit[u_pos2+1,pos+1]-h_fit[u_pos2,pos+1])/(u_fit[u_pos2+1,pos+1]-u_fit[u_pos2,pos+1])*u_fit[u_pos2+1,pos+1]    
+    h_pos=(h2-h1)/(v_fit[pos+1]-v_fit[pos])*v_su+h2-(h2-h1)/(v_fit[pos+1]-v_fit[pos])*v_fit[pos+1]
+
+
+    return [T_pos,p_pos,v_su,u_su,h_pos,0.]
+
+
+
+
+#########################################################################################################################################
+
+def getalp():
+    global z_it
+    '''
+    Berechnet den Wärmeübergangskoeffizient Gas/Zylinderwand
+    Woschni Korrelation
+    '''
+    if z_it[i,4]==0 or z_it[i,4]==2: #Ventile geschlossen
+        k=2.28
+    else: #Ventile offen, Ansaugen oder Ausstoßen
+        k=5.18
+    v_p=abs(z_it[i,1]-z_it[i-1,1])/((z_it[i,0]-z_it[i-1,0])/(2.*pi*pV[7])) # dX/dt  
+    alp=127.93*pV[0]**(-.2)*(z_it[i-1,6]/100.)**.8*(z_it[i-1,5]+273.15)**(-.55)*(k*v_p)**.8
+    z_it[i,13]=alp
+
+    return alp
+
+def Zustand_th_Masse(Q):
+    global z_it
+    '''Berechnet die Temperaturänderung der thermischen Masse als Funktion des
+    Wärmeübergang Innen (Q) und des Wärmeübergangs zur Umgebung (Q_u)
+    '''
+    ### Masse und cv der thermischen Masse sind im stationären Betrieb nicht 
+    ### entscheidend, die Parameter sind so gewählt, dass das Modell schnell konvergiert
+    ### aber keine Schwingungen auftreten
+    m=.0001 #kg
+    cv=.502 #kJ/kg/K
+    alp_a=6. #Wärmeübergangskoeffizient zur Umgebung
+    A=Ver0[3]*pV[8]/Ver0[2]*pV[0]/Ver0[0]*pV[1]/Ver0[1] # Außenfläche Zylinder abgeschätzt über Geometrie bezogen auf Fitting-Verdichter
+    Q_u=alp_a*A*(Tu-z_it[i-1,12])*((z_it[i,0]-z_it[i-1,0])/(2.*pi*pV[7]))*1e-3# kJ
+    z_it[i,12]=(Q+Q_u)/cv/m+z_it[i-1,12]
+    return
+
+
+def Kompression():
+    global z_it
+    Schritt=0
+    W=-z_it[i-1,6]*(z_it[i,2]-z_it[i-1,2])#Kompressionsarbeit, kJ
+    Wr=-pV[5]*(z_it[i,2]-z_it[i-1,2])#Reibarbeit, kJ
+    getalp()
+    Q=z_it[i,13]*z_it[i,3]*(z_it[i-1,12]-z_it[i-1,5])*((z_it[i,0]-z_it[i-1,0])/(2.*pi*pV[7]))*1e-3# kJ
+    Zustand_th_Masse(-Q)
+    dm=0. #Keine Ein- oder Ausströmende Masse
+    mi=z_it[i-1,11] #Masse im Zylinder, kg
+    ui=(Q+W+Wr)/mi+z_it[i-1,8]# kJ/kg
+    vi=z_it[i,2]/mi #Spezifisches Volumen im Zylinder, m³/kg
+    # Abfrage Stoffdaten Zustand i mit u und v, zurück: T, p, v, u, h, s
+
+    try:
+        zi=Stoffdaten(ui,vi)
+    except:
+        print("Fehler_1")
+        zi=fl.state(['u','v'],[ui,vi],fluid,Eh="CKPa")[0:6]
+
+    # Setzen der Daten in z_it
+    z_it[i,4]=Schritt
+    z_it[i,5:11]=zi
+    z_it[i,11]=mi
+    z_it[i,14:16]=dm,Q
+    return
+
+def Ausschieben():
+    global z_it
+    Schritt=1
+    W=-z_it[i-1,6]*(z_it[i,2]-z_it[i-1,2])#Kompressionsarbeit, kJ
+    Wr=-pV[5]*(z_it[i,2]-z_it[i-1,2])#Reibarbeit, kJ
+    getalp()
+    Q=z_it[i,13]*z_it[i,3]*(z_it[i-1,12]-z_it[i-1,5])*((z_it[i,0]-z_it[i-1,0])/(2.*pi*pV[7]))*1e-3# kJ
+    Zustand_th_Masse(-Q)
+    m_dot=pZyk[1]/z_it[i-1,7]*sqrt(2.*(z_it[i-1,6]-pZ[6])*1000.*z_it[i-1,7])#Massenstrom der den Zylinder verlässt, kg/s
+    dm=m_dot*((z_it[i,0]-z_it[i-1,0])/(2.*pi*pV[7])) #Ausgeschobene Masse, kg
+    mi=z_it[i-1,11]-dm #Resultierende Masse im Zylinder, kg
+    vi=z_it[i,2]/mi #Resultierendes spezifisches Volumen im Zylinder, m³/kg
+    ui=(Q+W+Wr-dm*z_it[i-1,9]+z_it[i-1,11]*z_it[i-1,8])/mi #Energiebilanz, resultierende Innere Energie im Zylinder, kj/kg
+    # Abfrage Stoffdaten Zustand i mit u und v
+
+    try:
+        zi=Stoffdaten(ui,vi)
+    except:
+        print("Fehler_2")
+        zi=fl.state(['u','v'],[ui,vi],fluid,Eh="CKPa")[0:6]
+    # Setzen der Daten in z_it
+    z_it[i,4]=Schritt
+    z_it[i,5:11]=zi
+    z_it[i,11]=mi
+    z_it[i,14:16]=dm,Q
+    return
+
+
+def Expansion():
+    global z_it
+    Schritt =2            
+    W=-z_it[i-1,6]*(z_it[i,2]-z_it[i-1,2])#Kompressionsarbeit, kJ
+    Wr=pV[5]*(z_it[i,2]-z_it[i-1,2])#Reibarbeit, kJ
+    getalp()
+    Q=z_it[i,13]*z_it[i,3]*(z_it[i-1,12]-z_it[i-1,5])*((z_it[i,0]-z_it[i-1,0])/(2.*pi*pV[7]))*1e-3# kJ
+
+    Zustand_th_Masse(-Q)
+    dm=0. #Keine Ein- oder Ausströmende Masse
+    mi=z_it[i-1,11] #Masse im Zylinder, kg
+    ui=(Q+W+Wr)/z_it[i-1,11]+z_it[i-1,8]#Energiebilanz
+    vi=z_it[i,2]/mi #Resultierendes spezifisches Volumen im Zylinder, m³/kg
+    # Abfrage Stoffdaten Zustand i mit u und v
+
+    try:
+        zi=Stoffdaten(ui,vi)
+    except:
+        print("Fehler_3")
+        zi=fl.state(['u','v'],[ui,vi],fluid,Eh="CKPa")[0:6]
+    # Setzen der Daten in z_it
+    z_it[i,4]=Schritt
+    z_it[i,5:11]=zi
+    z_it[i,11]=mi
+    z_it[i,14:16]=dm,Q
+    return
+
+
+def Ansaugen():
+    global z_it
+    Schritt=3
+    W=-z_it[i-1,6]*(z_it[i,2]-z_it[i-1,2])#Kompressionsarbeit, kJ
+    Wr=pV[5]*(z_it[i,2]-z_it[i-1,2])#Reibarbeit, kJ
+    getalp()
+    Q=z_it[i,13]*z_it[i,3]*(z_it[i-1,12]-z_it[i-1,5])*((z_it[i,0]-z_it[i-1,0])/(2.*pi*pV[7]))*1e-3# kJ
+    Zustand_th_Masse(-Q)
+
+    m_dot=pZyk[0]/pZ[2]*sqrt(2.*(pZ[1]-z_it[i-1,6])*1000*pZ[2]) #Massenstrom der den Zylinder verlässt, kg
+    dm=m_dot*((z_it[i,0]-z_it[i-1,0])/(2.*pi*pV[7])) #Ausgeschobene Masse, kg
+    mi=z_it[i-1,11]+dm #Resultierende Masse im Zylinder, kg
+    vi=z_it[i,2]/mi #Resultierendes spezifisches Volumen im Zylinder, m³/kg
+    ui=(Q+W+Wr+dm*pZ[4]+z_it[i-1,11]*z_it[i-1,8])/mi #isenthalpe Drosselung auf Zylinderdruck
+    # Abfrage Stoffdaten Zustand i mit u und v
+
+    try:
+        zi=Stoffdaten(ui,vi)
+    except:
+        print("Fehler_4")
+        zi=fl.state(['u','v'],[ui,vi],fluid,Eh="CKPa")[0:6]
+
+    # Setzen der Daten in z_it
+    z_it[i,4]=Schritt
+    z_it[i,5:11]=zi
+    z_it[i,11]=mi
+    z_it[i,14:16]=dm,Q
+    return
+
+
+
+def ProzessIter():
+    global pZyk, z_it, i
+
+    #Setzen von Aeff_e, explizite Funktion
+    M=fl.get_fluid_info(fluid,"CKPa")["molar_mass"] #Molmasse kg/mol
+    pZyk[0]=2.0415e-3*(Rm/M)**(-.9826)*pV[0]**2./Ver0[0]**2. # Effektiver Strömungsquerschnitt Eintritt, m²
+
+    #Setzen von Aeff_a, implizite Funktion bzgl der mittleren Massenstromdichte über Ventil
+    #Bei der 1. Iteration ist Massenstromdichte unbekannt, typischer Wert wird gesetzt
+    pZyk[1]=1.5e-5*pV[0]**2./Ver0[0]**2.
+
+
+    while 1:
+        for i in range(1,IS):
+            if z_it[i,0]<=pi:
+                if z_it[i-1,6]<=pZ[6]:
+                    Kompression()
+                else:
+                    Ausschieben()
+            else:
+                if z_it[i-1,6]>=pZ[1]:
+                    Expansion()
+                else:
+                    Ansaugen()
+
+        #Fehlerquadrahtsumme T, p, T_th_mittel          
+        error=sqrt((z_it[-1,5]-z_it[0,5])**2.)+sqrt((z_it[-1,6]-z_it[0,6])**2.)+\
+                sqrt((average(z_it[-1,12])-average(z_it[0,12]))**2.)
+        if error<.01:
+            break
+        else:
+            #print("New loop")
+            Zellen_ausschieben=where(z_it[:,4]==1)[0]
+            m_aus=sum(z_it[Zellen_ausschieben,14]) #Insgesamt ausgeschobene Masse
+            t_aus=(z_it[Zellen_ausschieben[-1],0]-z_it[Zellen_ausschieben[0],0])/(2.*pi*pV[7])#Zeit des Ausschiebens
+            m_dichte=m_aus/t_aus/pZyk[1] #Massenstromdichte kg/s/m²
+            pZyk[1]=5.1109e-4*(m_dichte)**(-.486)*pV[0]**2./Ver0[0]**2. #Aeff_a neu
+            z_it[0,5:14]=z_it[-1,5:14] #Endwerte letzter Zyklus = Anfangswerte nächster Zyklus
+
+
+    # Auswertung Wirkungsgrade    
+    Zellen_ausschieben=where(z_it[:,4]==1)
+    m_aus=sum(z_it[Zellen_ausschieben,14]) #Insgesamt ausgeschobene Masse
+    m0=pi*pV[0]**2.*pV[1]/pZ[2]/4. #Angesaugte Masse idealer Verdichter
+    Liefer=m_aus/m0 #Liefergrad
+
+    h_aus=sum(z_it[Zellen_ausschieben,9]*z_it[Zellen_ausschieben,14])/m_aus #Mittlere Enthalpie Ausschieben
+    h_aus_s=fl.state(['p','s'],[pZ[6],pZ[5]],fluid,Eh="CKPa")["h"]#Austrittsenthalpie isentrop
+    Isentr=(h_aus_s-pZ[4])/(h_aus-pZ[4])#Isentroper Wirkungsgrad
+
+    return Isentr, Liefer
+
+
+def getETA(T_e,p_e,p_a,fluid_in):
+    global pV, pz, z_it, fluid
+    fluid=fluid_in
+
+###############################     Parametersatz spezifisch für Verdichter   ###################################
+
+    pV=[34e-3,34e-3,3.5,.04,.1,80.23,50.,50./2.,2.] # Parameter siehe oben
+
+#################################################################################################################
+    pZ[0:6]=fl.state(['T','p'],[T_e,p_e],fluid,Eh="CKPa")[0:6] #Zusatnd Saugleitung
+    pZ[6]=p_a #Druck in Druckleitung
+
+    ############### Setze Geometrie ##################################
+    z_it[:,0]=linspace(0.,2*pi,IS)
+    z_it[:,1]=-(pV[1]/2.*(1.-cos(z_it[:,0])+pV[2]*(1.-sqrt(1.-(1./pV[2]*sin(z_it[:,0]))**2.))))+pV[4]*pV[1]+pV[1]#Kolbenpositionen, x=0 bei UT
+    A_kopf=pi/4.*pV[0]**2.
+    z_it[:,2]=A_kopf*z_it[:,1]# Zylindervolumina
+    z_it[:,3]=pi*pV[0]*z_it[:,1]+2.*A_kopf #Wärmeübertragungsflächen
+
+    ########## Setze Startbedingungen im Zylinder ##################
+    z_it[0,5:11]=pZ[0:6]    
+    z_it[0,11]=z_it[0,2]/z_it[0,7] #V/v, Zylinder vollständig gefüllt mit Sauggas
+    ##### Starttemp thermische Masse, pro Zyklusdurchlauf gibt es stehts nur eine Temperatur
+    ##### Startwert frei wählbar, beeinflusst maßgeblich Iterationszeit.
+    z_it[:,12]=42.
+    Isentr,Liefer=ProzessIter()
+    return Isentr,Liefer
+
+
+
+
+
+
+
+
+
+
+