Solving the First Law of Thermodynamics using Python

The first law of thermodynamics is related to energy conservation, i.e., if energy in form disappears then the energy will in appear in some other form. In thermodynamics majorly we are concerned about heat, work and internal energy. Heat and work are the form of energy which are also called as "energy in transit", hence they are path functions. They cannot be stored by a system, whereas internal energy is the property of the system which can be stored by system.

For a closed system, the first law is written as −

$$\mathrm{\Sigma Q=\Sigma W}$$

But this is only valid for a closed system which is undergoing a cycle. For a process it can be written as −

$$\mathrm{Q_{1-2}=\triangle U+W_{1-2}}$$

Before we can start modelling these one should know how to model work and heat for different processes. So let us do that first.

If the only impact on the system's external environment can be reduced to weight increase, then the system is said to have performed thermodynamic work. In actual circumstances, the weight may not be lifted, but the effect may still be that of raising weights. The work might be carried out by or on the system depending on the direction. The former is regarded as favourable and the latter as negative.

Work is a route function since it depends on the path taken between two states as well as the initial and ending states. We shall concentrate on non-dissipative/displacement work in this piece. It refers to the region under the p-v graph for a closed system.

For a general polytropic process represented by process $pv^{n}$=Constant, the work done depends on the index 𝑛. Table shown below will tells about the different processes based on the index 𝑛.

n	Process	𝒑−𝒗 relation
0	Isobaric process (p=constant)	p = c
1	Isothermal process for an ideal gas (T=constant)	pv =c
1.4	Adiabatic process (no heat interaction)	$\mathrm{pv^{y}\:=c}$
∞	Isochoric process (constant volume)	v = c

To derive the work interaction, as given in Eq. (3), the typical technique is to integrate the process from the beginning state to the final state.

$$\mathrm{W=\int ^{2}_{1}pdv}$$

As a result, the Table shown below contains a list of the work interaction for different process.

Process	p-v relation
Isobaric process (p = c)	$\mathrm{p(v_{2}-v_{1})}$
Isothermal process for an ideal gas (pv=c)	$\mathrm{p_{1}v_{1}\mathbb{I}n\frac{v_{2}}{v_{1}}=p_{1}v_{1}\mathbb{I}n\frac{p_{1}}{p_{2}}}$
Polytropic process ($pv^{n}=c)$	$\mathrm{\frac{\mathrm{p_{1}v_{1}-p_{2}v_{2}}}{n-1}}$
Adiabatic process ($pv^{y}=c$)	$\mathrm{\frac{\mathrm{p_{1}v_{1}-p_{2}v_{2}}}{\gamma -1}}$
Isochoric process (v=c)	0

Heat transfer is defined as the energy transition between a system and its surroundings which comes into play due to the difference of temperature. Conduction, convection, and radiation are its three main classifications.

The heat that enters the system is regarded as positive, whereas the heat that exits the system is regarded as negative. When a system doesn't exchange heat with its surroundings, it is said to be undergoing an adiabatic process. The quantity of energy in the form of heat needed to increase the temperature of a unit mass of a substance by a unit degree is known as specific heat.

The expression for specific heat is shown in equation (4). The heat capacity of a substance is represented by the symbol 𝐶 (𝐶=𝑚𝑐), which is the product of mass and specific heat. For solids and liquids, the specific heat is unaffected by the process; however, for gases, the specific heat must be qualified by the method. Thus, it differs for operations involving constant pressure and volume for a gas.

$$\mathrm{Q=mc\Delta t}$$

Latent heat (L) is the quantity of energy in the form of heat needed to alter a substance's phase measured per unit mass. Latent heat transformation does not cause a change in temperature or pressure.

The total heat interaction for a unit mass of a substance during the phase shift is given by −

$$Q=mL$$

Evaluating Heat and Work using Python Programming

The table shown below shows the different Python functions to find work and heat interaction. Moreover, they will also be used to plot the process and area under it.

from pylab import *

# Font and size
font = {'family' :'Times New Roman','size' : 16}
rc('font', **font)

# Constant Volume Process
def w_cv(p1,p2,v):

   """
   Function for constant volume process
   Input: Pressure limits p1 & p2 (in kPa) and volume (m^3)
   Output: Work output and p-v plot
   """
   
   figure(1,figsize=(7.20,4.80))
   
   p=linspace(min(p1,p2),max(p1,p2),30)
   v=zeros(30)+v
   
   plot(v,p,'k-',linewidth=4)
   xlabel('Volume (m$^3$)')
   ylabel('Pressure (kPa)')
   
   show()
   return 0
   
# Constant Pressure Process
def w_cp(p,v1,v2):

   """
   Function for constant pressure process
   Input : Volume limits v1 & v2 (m^3) and pressure (kPa)
   Output: Work output and p-v plot
   """
   
   work = p*(v2-v1)
   figure(2,figsize=(7.20,4.80))
   v=linspace(min(v1,v2),max(v1,v2),30)
   p=zeros(30)+p
   
   plot(v,p,'k-',linewidth=4)
   fill_between(v,p,color='red',alpha=0.3)
   xlabel('Volume (m$^3$)')
   ylabel('Pressure (kPa)')
   
   show()
   return round(work,3)
# Constant Temperature Process
def w_ct(p1,v1,v2):
   """
   Function for constant temperature process
   Input : Initial pressure (kPa) and volume (m^3) and final volume (m^3)
   Output: Work output and p-v plot
   """
   
   work = p1*v1*log(v2/v1)
   figure(2,figsize=(7.20,4.80))
   
   v=linspace(min(v1,v2),max(v1,v2),30)
   c=p1*v1
   p=c/v
   
   plot(v,p,'k-',linewidth=4)
   fill_between(v,p,color='red',alpha=0.3)
   xlabel('Volume (m$^3$)')
   ylabel('Pressure (kPa)')
   
   show()
   return round(work,3)
   
# Polytropic Process
def w_pol(p1,p2,v1,v2,n):
   """
   Function for polytropic process
   Input : p1 & p2 (kPa), v1 & v2 (m^3), and n (m^3)
   If the process is adiabatic write n=1.4
   Output: Work output and p-v plot
   """
   
   work = (p1*v1-p2*v2)/(n-1)
   figure(2,figsize=(7.20,4.80))
   
   v=linspace(min(v1,v2),max(v1,v2),30)
   c=p1*v1**n
   p=c/v**n
   
   plot(v,p,'k-',linewidth=4)
   fill_between(v,p,color='red',alpha=0.3)
   xlabel('Volume (m$^3$)')
   ylabel('Pressure (kPa)')
   
   show()
   return round(work,3)
# Sensible Heat Transfer
def q_sh(m,c,T1,T2):
   """
   
   Function for the evaluation of sensible heat transfer
   Input : mass (m), sp. heat (c), initial temp (T1), final temp (T2)
   Output: heat transfer
   """
   
   return m*c*(T2-T1)
   
# Heat transfer during phase change
def q_mel_sol(m,Ti,Tf,T_pc,c_bpc,c_apc,L):
   """
   Function for the evaluation of heat transfer during phase change
   Input: mass (m), initial temp (Ti), final temp (Tf),
   phase change temp (T_pc), sp. heat below phase change (c_bpc),
   sp. heat above phase change (c_apc), latent heat (L)
   Output: heat interaction
   """
   
   if Ti>Tf:
      print('Process is either freezing or condensation')
      return m*(c_bpc*(T_pc-Ti)-L+c_apc*(Tf-T_pc))
   else:
      print('Process is either melting or vaporization')
      return m*(c_bpc*(T_pc-Ti)+L+c_apc*(Tf-T_pc))

Let us see how to implement these functions with the help of two examples.

Example 1

Air having a mass of 1.5 kg is compressed from 0.2 MPa to 1 MPa according to the isothermal process. The air density at the initial point is 1.21 kg/m3. Find the work done and plot it in p-v diagram.

Solution

The program will be as follows −

# Input Data
m=1.50 # kg
p1=0.2E3 # kPa
p2=1.0E3 # kPa
ρ1=1.21 # kg/m^3

# Volume at the beginning
v1=m/ρ1

# Volume at end of process
v2=(p1*v1/p2)

# Given Process is isothermal
# The function will be: w_ct
Work = w_ct(p1,v1,v2)

print('Work output = ', Work)

On execution, you will get the following output −

Work output = -399.034

Example 2

650 kg of Tuna at 5°C are to be frozen and stored at -12°C. The specific heat above and below the freezing points (-2°C) are 3.182 and 1.717 kJ/kg-K, respectively. For a latent heat of fusion of 234.5 kJ/kg, the net heat removed from Tuna is to be evaluated.

Solution

The Python program to solve above problem will be as follows −

# Input data
m=650 # kg
Ti=5 # deg C
Tf=-12 # deg C
T_pc=-2 # deg C
c_bpc=3.182 # kJ/kg
c_apc=1.717 # kJ/kg
L=234.5

# Function calling
Q_net=q_mel_sol(m,Ti,Tf,T_pc,c_bpc,c_apc,L)
print("Q = ",Q_net)

Output

The program output will be −

Process is either freezing or condensation
Q = -178063.6

Now the question comes that how to model Eq. 1 and 2 in the code. As you have already evaluated the functions for heat and work so you can very easily model the first law. Lets take two example problems to solve it.

Example 3

In Example 1, the system has lost the heat 5 kJ of heat then evaluate the change in internal energy of the system.

Solution

First the code for work will remain the same (as in Example 1) the only change will be that of implementation of Eqution-1, so the full code will be −

# Input Data
m=1.50 # kg
p1=0.2E3 # kPa
p2=1.0E3 # kPa
ρ1=1.21 # kg/m^3

# Volume calculation at point 1
v1=m/ρ1

# Calculation of volume at point 2
v2=(p1*v1/p2)

# Evaluation of W_12
# The function used is w_ct
w = w_ct(p1,v1,v2)

# Q_12
Q=-5

# Equation 1 says:
Δu=Q-w

print(f"Δu = {Δu} kJ")

Output

Program output will be −

Δu = 394.034 kJ

Therefore, the internal energy of the system rises by 394.034 kJ

Example 4

A closed system undergoes three heat and two work exchanges as follows −

Q (kJ)	+10	-30	-40
W (kJ)	+100	-50	?

Find the missing energy interaction.

Solution

from numpy import *

nq=int(input("Enter number of heat interactions known: "))
nw=int(input("Enter number of work interactions known: "))
Q=empty(nq)
W=empty(nw)

for i in range(nq):
   Q[i]=float(input(f' Enter {i+1} heat interaction: '))
   
ΣQ=sum(Q)

for i in range(nw):
   W[i]=float(input(f' Enter {i+1} work interaction: '))
ΣW=sum(W)

# Missing Work interaction
W=ΣQ-ΣW

print(f' Missing work interaction : {W} kJ')

Output

The program output will be −

Enter number of heat interactions known: 3
Enter number of work interactions known: 2
Enter 1 heat interaction: 10
Enter 2 heat interaction: -30
Enter 3 heat interaction: -40
Enter 1 work interaction: 100
Enter 2 work interaction: -50
Missing work interaction: -110.0 kJ

Conclusion

In this tutorial, we explained how to model and plot work interaction and heat interaction. Also, we have seen via two problems that how to implement Equation 1 and 2 of the first law for open system.