Troubleshooting usnistgov/fipy#816

desiccantwheel.py

# Paramaters from the Table 3.1: Transport and thermophysical properties of the adsorbent wheel tested.
T_ci     = 24.0                               # °C                             Exhaust Value in the simulation
T_hi     = 90.0                               # °C                             Inlet Value in the simulation
omega_ci = 0.012                              # kg_w kg_as^-1                  Exhaust Value in the simulation
omega_hi = 0.021                              # kg_w kg_as^-1                  Inlet Value in the simulation
lambda_d = 0.24/1000                          # kW m^-1 K^-1
N        = 0.5                                # rpm
rho_d    = 630                                # kg m^-3
delta    = 0.3/1000                           # m
qst      = 2650                               # kJ kg^-1
Dvs      = 3.5e-07                            # m2 s^-1
e_t      = 0.46                               # -
L        = 0.1                                # m
at_1     = 2.04                               # m^2
nad      = 5200.0                             # -
m_dot_a  = 30.6/3600                          # kg/s
W_max    = 0.35                               # kg_w/kg_d
C_s      = 1.1                                # -

# Properties of the air at atphosferic pressure
from CoolProp.CoolProp import PropsSI
rho_a    = PropsSI('D','P', 101325,'T',(T_ci+273),'Air')           # kg m^-3
c_pa     = PropsSI('C','P', 101325,'T',(T_ci+273),'Air')/1000      # kJ kg^-1 K^-1
lambda_a = PropsSI('L','P', 101325,'T',(T_ci+273),'Air')/1000      # kW m^-1 K^-1
c_pw     = PropsSI('C','P', 101325,'T',(T_ci+273),'Water')/1000    # kJ kg^-1 K^-1
D_va = 2.64e-5                                                     # m^2 s^-1

# The Nusslet and Shewrwood number its extraid from the Table 2.3 : Fully developed Nusselt and Sherwood...
# ...numbers for sinusoidal adsorbent ducts (ducts with adsorbent walls
Nu = 2.98
Sh = 1.39

# From the seccion 3.3 in the thrid paragraph.
# Hydraulic Diameter
A_c = 1.36/(1000**2)                               # m^2
P_f = 3.92/1000                                    # m
D_h = (4.*A_c)/(P_f)                               # m
# Heat transfer coefficient
h = Nu*lambda_a/D_h                                # kW m^-2 K^-1 Eq. (3.3)
# Mass transfer coefficient
k = Sh*D_va/D_h                                    # m s^-1       Eq. (3.4)

A_t = P_f * L* nad                                 # m^2
u_a = m_dot_a/(rho_a*A_c*nad)                      # m s^-1

# In the Table 2.2:  Physical and thermodynamic properties of the honeycomb adsorbent bed.
c_pd = 0.876                                       # kJ kg^-1 K^-1
k_pm = 0.32                                        # s^-1
Kp = W_max/omega_hi                                # Eq. (2.51)

# Coefficients for the equations
k_mu = (60.*k_pm*(omega_hi-omega_ci))/(N*W_max)    # Eq. (3.26)
c_2 = (60.*rho_d*k_pm)/(e_t*rho_a*N)               # Eq. (3.31)
D_vzu = (240.*Dvs)/(N* delta**2*e_t)               # Eq. (3.32)
Bi_m = (k*delta)/(2.*Dvs)                          # Eq. (3.38)
Bi = (h*delta)/(2.*lambda_d)                       # Eq. (3.37)
# Below the Eq (3.40)
c_3   = (N*L)/(60.*u_a)
NTU   = (h*at_1)/(m_dot_a*c_pa)
NTU_m = (rho_a*k*at_1)/(m_dot_a)

# x* value to sample air conditions
air_x = 0.01

# Mesh and equations for the desiccant side
import fipy as fp
n = 100
mesh1 = fp.Grid1D(nx = n, dx = 1/n)

T = fp.CellVariable(mesh=mesh1, name ='T', value=1., hasOld = True)
H = fp.CellVariable(mesh=mesh1, name ='H', value=1., hasOld = True)
water = fp.CellVariable(mesh=mesh1, name ='w', hasOld = True)

# initialize to inlet conditions
Hax = fp.Variable(name='Hax', value=1.)
Tax = fp.Variable(name='Tax', value=1.)
# Transient Parameters  

# Equelibrium Humidity
p = H * (omega_hi - omega_ci) + omega_ci        # Turn dimensional again   Ec. (3.20)
# p = Hax * (omega_hi - omega_ci) + omega_ci      # Turn dimensional again   Ec. (3.20)
RH = p*(fp.numerix.exp(5294/(24+273))/1e6)      # Relative Humidity of the air Ec. (3.16)
w_equ = W_max/(1 - C_s + (C_s/RH))              # Ec. (3.12)
a = w_equ/Kp                                    # Relation betwwen wapter uptake and humidity
omega_equ = (a-omega_ci)/(omega_hi - omega_ci)  # Turn dimensionless

# Coefficients for equation (3.27)
water.value = H*Kp
# c_1 = 43.364*0.32/c_ptot            # Ec. (3.29)
# lambda_zu = 2031.7/c_ptot           # Ec. (3.28)
w = water * 0.35                                                   # Ratio water absorption dimensional again Eq.(3.21)
c_ptot = (c_pd+w*c_pw)                                             # Eq.(3.9)
lambda_zu = (240.*lambda_d)/(N*delta**2*rho_d*c_ptot)              # Eq.(3.28)
c_1 = ((60*k_pm*qst)/(N*c_ptot))*((omega_hi-omega_ci)/(T_hi-T_ci)) # Eq.(3.29)



mask = mesh1.facesRight

# Values of Ha and Ta evaluated at x

# cellDistanceVectors was never implemented for Grid1D
cellDistanceVectors = fp.tools.numerix.zeros((1, mesh1.numberOfFaces), float)
cellDistanceVectors[:] = mesh1.dx
if mesh1.numberOfCells > 0:
    cellDistanceVectors[0, 0] = -mesh1.dx / 2.
    cellDistanceVectors[0, -1] = mesh1.dx / 2.

dPf = fp.FaceVariable(mesh=mesh1,
                      value=mesh1._faceToCellDistanceRatio * cellDistanceVectors)
norm = fp.FaceVariable(mesh=mesh1, value=mesh1.faceNormals, rank=1)

aT = fp.FaceVariable(mesh=mesh1, value=Bi * norm, rank=1)
bT = fp.FaceVariable(mesh=mesh1, value=1., rank=0)
gT = fp.FaceVariable(mesh=mesh1, value=Bi, rank=0) * Tax
Ts_coeff = mask * lambda_zu.faceValue * norm / (-dPf.dot(aT) + bT)

aH = fp.FaceVariable(mesh=mesh1, value=Bi_m * norm, rank=1)
bH = fp.FaceVariable(mesh=mesh1, value=1., rank=0)
gH = fp.FaceVariable(mesh=mesh1, value=Bi_m, rank=0) * Hax
Hs_coeff = mask * D_vzu * norm / (-dPf.dot(aH) + bH)

# Sources Terms
S01 = c_1*omega_equ
S11 = fp.ImplicitSourceTerm(coeff=c_1, var = H) 
S1 = S01 - S11        # Source Term for the Surface Energy equaiton Ec. (3.27)
S02 = c_2 * omega_equ
S22 = fp.ImplicitSourceTerm(coeff=c_2, var = H)
S2 = S02 - S22        # Source Term for the Surface Mass equation Ec. (3.30)
S05 = k_mu * omega_equ
S55 = fp.ImplicitSourceTerm(coeff=k_mu, var = H)
S5 = S05 - S55        # Source Term for the Adsorption equation Ec. (3.25)

# New Diffusion coefficients due to Robin Boundary condition
Gamma_eq1 = lambda_zu.faceValue * ~mask
Gamma_eq2 = D_vzu * ~mask

# Source generated from the BC of Robin Type
BCs1 = ((Ts_coeff * gT).divergence
        - fp.ImplicitSourceTerm(coeff=(Ts_coeff * norm.dot(aT)).divergence, var=T))
BCs2 = ((Hs_coeff * gH).divergence
        - fp.ImplicitSourceTerm(coeff=(Hs_coeff * norm.dot(aH)).divergence, var=H))

eq1 = (fp.TransientTerm(coeff=1, var=T) == fp.DiffusionTerm(coeff=Gamma_eq1, var=T)+ S1 + BCs1)   # Ec (3.27)
eq2 = (fp.TransientTerm(coeff=1, var=H) + S2 == fp.DiffusionTerm(coeff=Gamma_eq2, var=H) + BCs2)  # Ec (3.30)
eq5 = (fp.TransientTerm(coeff=1, var=water) == S5) # Ec (3.25)

eq12 = eq1 & eq2 & eq5

# Mesh and equations for the air side

mesh2 = fp.Grid1D(nx = n, dx = 1/n)

Ha = fp.CellVariable(mesh=mesh2, name ='Ha', value=1., hasOld = True)
Ta = fp.CellVariable(mesh=mesh2, name ='Ta', value=1., hasOld = True)

# Boundary Conditions

# Dirichlet Boundaries
Ha.constrain(1., where=mesh2.facesLeft)     # Ec. (3.6)
Ha.constrain(0., where=mesh2.facesRight)    # Ec. (3.6)
Ta.constrain(0., where=mesh2.facesLeft)     # Ec. (3.7)
Ta.constrain(1., where=mesh2.facesRight)    # Ec. (3.7)

# Values of H and T evaluated at z=z_s (z* = 1)
Hs = fp.Variable(name='Hs')
Ts = fp.Variable(name='Ts')

S03 = NTU * Ts
S33 = fp.ImplicitSourceTerm(coeff = NTU, var = Ta)
S3 = S03 - S33         # Source Term for the air stream energy transfer equation  #  Ec.(3.39)
S04 = NTU_m * Hs
S14 = fp.ImplicitSourceTerm(coeff = NTU_m, var = Ha)    
S4 = S04 - S14         # Source Term for the air strem mass transfer equation  # Ec. (3.40)


# Convective Vector Value
Conv_coeff = fp.FaceVariable(mesh = mesh1, value = [1.])

eq3 = (fp.TransientTerm(coeff = c_3, var = Ta) + fp.UpwindConvectionTerm(coeff = Conv_coeff, var = Ta) == S3)
eq4 = (fp.TransientTerm(coeff = c_3, var = Ha) + fp.UpwindConvectionTerm(coeff = Conv_coeff, var = Ha) == S4)

eq34 = eq3 & eq4

dt = 1/120
steps = 20
sweeps = 20
# solver = fp.LinearLUSolver(tolerance=1e-10, iterations=10)  #fp.LinearCGSSolver(tolerance = 1e-10, iteration = 10)
solver = None

vT = fp.Viewer(vars=(T, Ta))
vH = fp.Viewer(vars=(H, Ha))
input()

for step in range(steps):
    H.updateOld()
    T.updateOld()
    Ha.updateOld()
    Ta.updateOld()

    # Transfer air values to desiccant mesh
    Hax.value = Ha((air_x,), order=1)
    Tax.value = Ta((air_x,), order=1)

#     print("desiccant")
    res = 1e+10
    for sweep in range(sweeps):
        res = eq12.sweep(dt=dt)
#         print(res)

    # Transfer desiccant values to air mesh
    Hs.value = float(H.faceValue[..., mask])
    Ts.value = float(T.faceValue[..., mask])

#     print("air")
    res = 1e+10
    for sweep in range(sweeps):
        res = eq34.sweep(dt=dt)
#         print(res)

    vT.plot()
    vH.plot()

    print(step * dt, Tax * (T_hi-T_ci) + T_ci, Hax * (omega_hi - omega_ci) + omega_ci)
#     input()