Analyze the Stirling Cycle with Real-World Materials 

Examine the Stirling cycle for 10 moles of air.

In[1]:=

X mAir = 10 ChemicalData["Air", "MolarMass"]

Out[1]=

Calculate the isothermal expansion and compression curves at 200 °C and 100 °C, respectively.

In[2]:=

X isothermalExpansion = Transpose[{ ThermodynamicData["Air", "SpecificVolume", {"Temperature" -> Quantity[200, "DegreesCelcius"], "Pressure" -> Quantity[Range[100, 200] 1000, "Pascals"]}] mAir, Quantity[Range[100, 200] 1000, "Pascals"]}];

In[3]:=

X isothermalCompression = Transpose[{ ThermodynamicData["Air", "SpecificVolume", {"Temperature" -> Quantity[120, "DegreesCelcius"], "Pressure" -> Quantity[Range[100, 200] 1000, "Pascals"]}] mAir, Quantity[Range[100, 200] 1000, "Pascals"]}];

Create the diagram.

In[4]:=

X ListLinePlot[{isothermalExpansion, isothermalCompression}, Axes -> False, PlotRange -> All, FrameLabel -> Automatic, Frame -> True, GridLines -> {{0.2, 0.3}, {}}]

Out[4]=

Calculate mechanical work through .

In[5]:=

X W = Quantity[ NIntegrate[ Evaluate[ Interpolation[QuantityMagnitude[isothermalExpansion]][V]], {V, 0.2, 0.3}] - NIntegrate[ Evaluate[ Interpolation[QuantityMagnitude[isothermalCompression]][V]], {V, 0.2, 0.3}], "Joules"]

Out[5]=

Find maximal and minimal pressures for expansion and compression.

In[6]:=

X {pMax1, pMin1} = Quantity[Interpolation[QuantityMagnitude[isothermalExpansion]][{0.2, 0.3}], "Pascals"]

Out[6]=

In[7]:=

X {pMax2, pMin2} = Quantity[ Interpolation[QuantityMagnitude[isothermalCompression]][{0.2, 0.3}], "Pascals"]

Out[7]=

Calculate the change in internal energy.

In[8]:=

X U[T_, p_] := ThermodynamicData["Air", "InternalEnergy", {"Temperature" -> T, "Pressure" -> p}]*mAir; \[CapitalDelta]U = U[ Quantity[200, "DegreesCelcius"], pMax1 ] - U[ Quantity[120, "DegreesCelcius"] , pMin2]

Out[8]=

Out[8]=

Using work and the change in internal energy, calculate efficiency.

In[9]:=

X W/\[CapitalDelta]U

Out[9]=

Compare with maximal possible Carnot efficiency.

In[10]:=

X 1. - UnitConvert[Quantity[120, "DegreesCelcius"], "Kelvins"]/ UnitConvert[Quantity[200, "DegreesCelcius"], "Kelvins"]

Out[10]=