View Issue Details

IDProjectCategoryView StatusLast Update
0003368OpenFOAMContributionpublic2019-10-23 14:28
Reportervijay Assigned Towill  
PrioritynormalSeverityminorReproducibilityalways
Status closedResolutionno change required 
PlatformGNU/LinuxOSUbuntuOS Version15.04
Summary0003368: Correct enthalpy transport even when specie specific heats are a function of temperature
DescriptionThis is a continuation of bug report "0003363: Solved unphysical temperature fields in multi-species gas tranport solver". Through the commit 54f379f668485bd2af0d05518c6a9eb5e20ca5a9, the unphysical temperature fields issue in reactingFoam solver was overcome. This is Lewis Number = 1 formulation. In the bug report mentioned above, I had said that this rectification will work only as long as specie specific heats are not a function of temperature. "Will" didn't agree with my argument and had asked for a case to justify my argument. I was not able to come up with a OpenFOAM case, but here I have attached the necessary files to test the effect mathematically.

In conclusion, I was wrong and Will was right. The current implementation in reactingFoam solver gives correct enthalpy transport even when specie specific heats are a function of temperature.

Math Expressions.pdf
- Discusses the expressions in Exact form and the form implemented in OpenFOAM

Excercise:
- A one dimensional domain of length = 1 m. Two gases Hydrogen and Air. The mass fraction of Hydrogen varies linearly from 0 at x=0 to 1 at x=1. For Air it varies from 1 at x=0 to 0 at x=1.
- Temperature varies linearly from 275 K at x=0 to 2000K at x=1.

Case1:
- Cp_Hydrogen and Cp_Air constant with temperature
- Result: CpiConstWithTemperature.png
Case2:
- Cp_Hydrogen and Cp_Air varies with temperature (polynomial function).
- Result: CpiVaryWithTemperature.png
Steps To ReproduceNone
Additional InformationNone
TagsNo tags attached.

Activities

vijay

2019-10-11 14:18

reporter  

Math Expressions.pdf (518,029 bytes)
CpiConstWithTemp.py (1,454 bytes)   
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

# Domain size 
x_start, x_end = 0,1

# number of nodes
nx  = 101

# values of x
x = np.linspace(x_start, x_end, nx)

# Temperature
T = np.linspace(275,2000,nx)

# Specific heats 
# Hydrogen
Cp1 = 14000
# Air
Cp2 = 1005

# Mass fractions
Y1 = np.linspace(0,1,nx)
Y2 = 1-Y1

# Thermal conductivity
kappa1 = 0.4
kappa2 = 5e-2

kappa = np.zeros(nx)
kappa = kappa1*Y1+kappa2*Y2

Cp = np.zeros(nx)
Cp = Cp1*Y1+Cp2*Y2

h = np.zeros(nx)
h = Cp*T

alpha = np.zeros(nx)
alpha = kappa/Cp

# Values at cell centres
Exact       = np.zeros(nx-1)
OpenFoam    = np.zeros(nx-1)
x_c         = np.zeros(nx-1)

for i in range(0,nx-1):
    kappa_c = 0.5*(kappa[i]+kappa[i+1])
    alpha_c = 0.5*(alpha[i]+alpha[i+1])
    Cp1_c    = Cp1
    Cp2_c    = Cp2
    T_c     = 0.5*(T[i]+T[i+1])
    
    gradT   = (T[i+1]-T[i])/(x[i+1]-x[i])
    gradY1   = (Y1[i+1]-Y1[i])/(x[i+1]-x[i])
    gradY2   = (Y2[i+1]-Y2[i])/(x[i+1]-x[i])
    gradH   = (h[i+1]-h[i])/(x[i+1]-x[i])

    x_c[i]      = x[i]+0.5*(x[i+1]-x[i])
    Exact[i]    = kappa_c*gradT + alpha_c* (Cp1_c*T_c*gradY1 + Cp2_c*T_c*gradY2)
    OpenFoam[i] = alpha_c*gradH 


plt.plot(x_c,Exact,'b-',label="Exact")
plt.plot(x_c,OpenFoam,'r--',label="OpenFOAM")
plt.legend()
plt.xlabel('X [m]')
plt.ylabel('Enthalpy from Heat conduction and Specie diffusion')
plt.title('Cpi constant with temperature')
plt.grid()
plt.savefig('CpiConstWithTemperature.png')
CpiConstWithTemp.py (1,454 bytes)   
CpiVarWithTemp.py (2,733 bytes)   
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

# Domain size 
x_start, x_end = 0,1

# number of nodes
nx  = 101

# values of x
x = np.linspace(x_start, x_end, nx)

# Temperature
T = np.linspace(275,2000,nx)

# Specific heats 
# Hydrogen
Cp1 = np.zeros(nx)
Cp1 = 1000*(-1.64044e-18*pow(T,6) +1.22802e-14*pow(T,5) -3.70152e-11*pow(T,4) +5.6503e-8*pow(T,3) -4.4297e-5*pow(T,2) +0.0172878*T +11.8523)
# Air
Cp2 = np.zeros(nx)
Cp2 = 1000*(7.36435e-22*pow(T,6) -8.89851e-17*pow(T,5) +5.73158e-13*pow(T,4) -1.44058e-9*pow(T,3) +1.64424e-6*pow(T,2) -0.000626846*T +1.07935)

# Mass fractions
Y1 = np.linspace(0,1,nx)
Y2 = 1-Y1

# Thermal conductivity
kappa1 = 0.4
kappa2 = 5e-2

kappa = np.zeros(nx)
kappa = kappa1*Y1+kappa2*Y2

Cp = np.zeros(nx)
Cp = Cp1*Y1+Cp2*Y2

h1 = np.zeros(nx)
for i in range(0,nx):
    Te      = T[i]
    upLim   = 1000*(-1.64044e-18*pow(Te,7)/7 +1.22802e-14*pow(Te,6)/6 -3.70152e-11*pow(Te,5)/5 +5.6503e-8*pow(Te,4)/4 -4.4297e-5*pow(Te,3)/3 +0.0172878*pow(Te,2)/2 +11.8523*Te)
    Te      = T[0]
    lowLim  = 1000*(-1.64044e-18*pow(Te,7)/7 +1.22802e-14*pow(Te,6)/6 -3.70152e-11*pow(Te,5)/5 +5.6503e-8*pow(Te,4)/4 -4.4297e-5*pow(Te,3)/3 +0.0172878*pow(Te,2)/2 +11.8523*Te)
    h1[i]   = Cp1[0]*T[0]+upLim-lowLim

h2 = np.zeros(nx)
for i in range(0,nx):
    Te      = T[i]
    upLim   = 1000*(7.36435e-22*pow(Te,7)/7 -8.89851e-17*pow(Te,6)/6 +5.73158e-13*pow(Te,5)/5 -1.44058e-9*pow(Te,4)/4 +1.64424e-6*pow(Te,3)/3 -0.000626846*pow(Te,2)/2 +1.07935*Te)
    Te      = T[0]
    lowLim  = 1000*(7.36435e-22*pow(Te,7)/7 -8.89851e-17*pow(Te,6)/6 +5.73158e-13*pow(Te,5)/5 -1.44058e-9*pow(Te,4)/4 +1.64424e-6*pow(Te,3)/3 -0.000626846*pow(Te,2)/2 +1.07935*Te)
    h2[i]   = Cp2[0]*T[0]+upLim-lowLim

h = np.zeros(nx)
h = h1*Y1+h2*Y2

alpha = np.zeros(nx)
alpha = kappa/Cp

# Values at cell centres
Exact       = np.zeros(nx-1)
OpenFoam    = np.zeros(nx-1)
x_c         = np.zeros(nx-1)

for i in range(0,nx-1):
    kappa_c = 0.5*(kappa[i]+kappa[i+1])
    alpha_c = 0.5*(alpha[i]+alpha[i+1])
    Cp_c    = 0.5*(Cp[i]+Cp[i+1])
    h1_c    = 0.5*(h1[i]+h1[i+1])
    h2_c    = 0.5*(h2[i]+h2[i+1])
    
    gradT   = (T[i+1]-T[i])/(x[i+1]-x[i])
    gradY1   = (Y1[i+1]-Y1[i])/(x[i+1]-x[i])
    gradY2   = (Y2[i+1]-Y2[i])/(x[i+1]-x[i])
    gradH   = (h[i+1]-h[i])/(x[i+1]-x[i])

    x_c[i]      = x[i]+0.5*(x[i+1]-x[i])
    Exact[i]    = kappa_c*gradT + alpha_c* (h1_c*gradY1 + h2_c*gradY2)
    OpenFoam[i] = alpha_c*gradH 


plt.plot(x_c,Exact,'b-',label="Exact")
plt.plot(x_c,OpenFoam,'r--',label="OpenFOAM")
plt.legend()
plt.xlabel('X [m]')
plt.ylabel('Enthalpy from Heat conduction and Specie diffusion')
plt.title('Cpi variation with temperature')
plt.grid()
plt.savefig('CpiVaryWithTemperature.png')
CpiVarWithTemp.py (2,733 bytes)   
CpiConstWithTemperature.png (34,726 bytes)   
CpiConstWithTemperature.png (34,726 bytes)   
CpiVaryWithTemperature.png (34,350 bytes)   
CpiVaryWithTemperature.png (34,350 bytes)   

will

2019-10-23 14:28

manager   ~0010847

Great. Thanks for confirming that it works with variable specific heats.

Issue History

Date Modified Username Field Change
2019-10-11 14:18 vijay New Issue
2019-10-11 14:18 vijay File Added: Math Expressions.pdf
2019-10-11 14:18 vijay File Added: CpiConstWithTemp.py
2019-10-11 14:18 vijay File Added: CpiVarWithTemp.py
2019-10-11 14:18 vijay File Added: CpiConstWithTemperature.png
2019-10-11 14:18 vijay File Added: CpiVaryWithTemperature.png
2019-10-23 14:28 will Assigned To => will
2019-10-23 14:28 will Status new => closed
2019-10-23 14:28 will Resolution open => no change required
2019-10-23 14:28 will Note Added: 0010847