Dynamic simulation of reciprocating compressor-expansor.

© 2026 Universidad de Salamanca, Department of Applied Physics, Thermodynamics, Statistical Physics and Energy Optimization Research Group. 
Author: David Pérez Gallego. https://gredos.usal.es/browse?authority=487d6919-66c1-4547-b165-abc0ce69d7d6&type=author
ORCID: https://orcid.org/0000-0003-2519-6014
Google Scholar : https://scholar.google.com/citations?user=MuzFipMAAAAJ&hl=es
ResearchGate: https://www.researchgate.net/profile/David-Perez-Gallego

This work is licensed under the Creative Commons Attribution–NonCommercial–ShareAlike 4.0 International License (CC BY - NC - SA 4.0) . You are free to use, reproduce, and adapt this work for non - commercial purposes, provided that proper attribution is given to the author and the same license is maintained for any derivative works .

Start and results

Clear[CompReciprocante];CompReciprocante:=Module[{Tx,px,mx,Vx,mdotin,mdotout,n,ϕ,ValveAdm,ValveEsc},
ϕ=ω*(t-1)*Δt;

Clear[Tx,px,mx,Vx,mdotin,mdotout,ValveAdm,ValveEsc];
(*padm>px*)If[padm>p[[t-1]]*1.25,ValveAdm=Cd*((AValve*padm)/(Sqrt[γ R Tadm]))γ((((px/padm)))^(1/γ))((((2/(γ-1))(1-((((px/padm)))^((γ-1)/γ)))))^0.5);,ValveAdm=0.;];
(*px>pesc*)If[p[[t-1]]>pesc*1.6,ValveEsc=-0.02*Cd*((AValve px)/(Sqrt[γ R Tx]))γ((((pesc/px)))^(1/γ))((((2/(γ-1))(1-((((pesc/px)))^((γ-1)/γ)))))^0.5);,ValveEsc=0.;];


{T[[t]],p[[t]],m[[t]],V[[t]],mdotin,mdotout}={Tx,px,mx,Vx,mdotin,mdotout}/.NSolve[
{
((Vx-V[[t-1]])/(Δt))==(V0/2)(r-1)(Sin[ϕ]+((Sin[ϕ]Cos[ϕ])/(Sqrt[(1/(f^2))-((Sin[ϕ])^2)])))ω,
((Tx-T[[t-1]])/(Δt))==(1/(mx *Cpf[T[[t-1]],p[[t-1]]/101325.]))((Qdot[ϕ]*(Tp-Tx))(*Qdot*)+mdotin*hadm+mdotout*hesc-((mx-m[[t-1]])/(Δt))*hTP[T[[t-1]],p[[t-1]]/101325.]+((px-p[[t-1]])/(Δt))*Vx),
mdotin==ValveAdm,
mdotout==ValveEsc,
((mx-m[[t-1]])/(Δt))==mdotin+mdotout,


(*((px-p[[t-1]])/(Δt))==γ px((1/mx)*((mx-m[[t-1]])/(Δt))+(1/(mx Cp Tx))(Qdot+mdotin*hin+mdotout*hout-((mx-m[[t-1]])/(Δt))*h)-(1/Vx)((Vx-V[[t-1]])/(Δt)))*)
px*Vx==(mx/Mm) R Tx

}

,{Tx,px,mx,Vx,mdotin,mdotout},Reals][[1]];




];

Off[NSolve::illcnd];
SetDirectory[NotebookDirectory[]];
fluid="Air";
(*Rango: 1<P<15;  300<T<1700*)
EntropiaTP="Air_properties/sTP_"<>fluid<>".txt";
EntalpiaTP="Air_properties/hTP_"<>fluid<>".txt";
EntalpiaPS="Air_properties/hPS_"<>fluid<>".txt";
TemperaturaPh="Air_properties/TPh_"<>fluid<>".txt";
TemperaturaPS="Air_properties/TPS_"<>fluid<>".txt";
Clear[sTP,hTP,hPS,TPh,TPS,T,P,p,S,h];
sTP[T_,P_]:=Evaluate[ToExpression[Import[EntropiaTP,"Text"]]];
hTP[T_,P_]:=Evaluate[ToExpression[Import[EntalpiaTP,"Text"]]];
hPS[P_,S_]:=Evaluate[ToExpression[Import[EntalpiaPS,"Text"]]];
TPh[P_,h_]:=Evaluate[ToExpression[Import[TemperaturaPh,"Text"]]];
TPS[P_,S_]:=Evaluate[ToExpression[Import[TemperaturaPS,"Text"]]];
ConductividadFluido="Air_properties/Conductividad_"<>fluid<>".txt";
DensidadFluido="Air_properties/Densidad_"<>fluid<>".txt";
ViscosidadFluido="Air_properties/Viscosidad_dinamica_"<>fluid<>".txt";
CapacidadFluido="Air_properties/Cpf_"<>fluid<>".txt";

Clear[Kfaire];Kfaire[T_,P_]:=Evaluate[ToExpression[Import[ConductividadFluido,"Text"]]];
Clear[ρfaire];ρfaire[T_,P_]:=Evaluate[ToExpression[Import[DensidadFluido,"Text"]]];
Clear[μfaire];μfaire[T_,P_]:=Evaluate[ToExpression[Import[ViscosidadFluido,"Text"]]];
Clear[Cpf];Cpf[T_,P_]:=Evaluate[ToExpression[Import[CapacidadFluido,"Text"]]];










rpm=2950.;(*rpm*)
tf=1/(rpm/60.)*6;(*tf*)
Δt=0.0001;
Nt=IntegerPart[tf/Δt]+1;

Clear[V,m,p,T];
V=m=p=T=ConstantArray[0,Nt];
T0=63.6+273.;(*K*)
p0=62600.;(*Pa*)

padm=59157.36;(*Pa*)
Tadm=32+273.;
hadm=hTP[Tadm,padm/101325.];

pesc=620000.;
Tesc=116.85+273.;
hesc=hTP[Tesc,pesc/101325.];


Cd=0.6;
l=33/1000.;(*m  longitud biela*)
a=15.4/1000/2.;(*m   Radio cigüeñal*)
f=a/l;
V0=85./(1000^3);(*(m^3) espacio nocivo*)
ω=2π*rpm/60;

Diam=21.1/1000;(*m diámetro del pistón*)
At=((π (Diam^2))/4);

DValve=21.1/4000;(*m  diámetro valve escape y admisión*)
AValve=((π (DValve^2))/4);(*m  área de valvula de admisión/escape*)

Vdt=At*2a;
r=((Vdt+V0)/V0);(*relación de presión*)


γ=1.4;
R=8.314;


(*Aire*)
Mm=28.97/1000.;(*kg/mol*)






Tp=82.8+273;(*T pared*)

Clear[Agas];Agas[ϕ_]:=π*(Diam^2)/2+2a*Sin[ϕ]*π*Diam;(*área de intercambio de calor gas-paredes*)

k=0.024;(*aire*)

L=2a;(*longitud característica*)
ρ=531;u=1.;μ=1.7*(10^(-5));
Reyd=((ρ u L)/μ);
Nu=0.035*(Reyd^0.7);(*Woschni*)
hconv=Nu*k/L;(*h convectivo*)

Clear[Qdot];Qdot[ϕ_]:=hconv*Agas[ϕ](*(Tp-Tx) *);





V[[1]]=V0;
p[[1]]=p0;
T[[1]]=T0;
m[[1]]=((p[[1]]*V[[1]])/(R*T[[1]]))*Mm;

For[t=2,t<=Nt,t++,
CompReciprocante
];

ListPlot[Table[{(t-1)*Δt,V[[t]]},{t,1,Length[V]}],Frame->True,FrameStyle->Directive[12,Black],FrameLabel->{"t [s]","V [(m^3)]"},Joined->True,PlotRange->All]

ListPlot[Table[{(t-1)*Δt,T[[t]]-273.},{t,1,Length[T]}],Frame->True,FrameStyle->Directive[12,Black],FrameLabel->{"t [s]","T [ºC]"},Joined->True,PlotRange->All]

ListPlot[Table[{(t-1)*Δt,m[[t]]},{t,1,Length[m]}],Frame->True,FrameStyle->Directive[12,Black],FrameLabel->{"t [s]","m [kg]"},Joined->True,PlotRange->All]


ListPlot[Table[{(t-1)*Δt,p[[t]]},{t,1,Length[p]}],Frame->True,FrameStyle->Directive[12,Black],FrameLabel->{"t [s]","p [Pa]"},Joined->True,PlotRange->All]

ListPlot[Table[{V[[t]],p[[t]]},{t,1,Length[p]}],Frame->True,FrameStyle->Directive[12,Black],FrameLabel->{"V [(m^3)]","p [Pa]"},Joined->True,PlotRange->All]