Erroneus results in no-linear coupled equations.

[AitorCendoya]

I am trying to model a diseccant wheel. The equations system, it's strong coupled and no-linear. Is as following are attached the picture of equations. The results are erroneous compared to the literature example..

• I used different typs of solvers, steps time and nodes.
• When i aply the Dirichlet boundanry conditions and print the mesh.faces, the value it's 100 times greater.
• Is it correct to place Robin's boundary condition in a function? (Because the parameters changes with a variable to be solved).
  • When i change the numbers of iteration or increase a steptime, most of the time print de "RuntimeError : Factor is exactly singular". The script

[AitorCendoya]

Try to organize more the solution of the problem, attach a PDF with the development of the problem. GitHub.pdf

I still have the same questions.

• Are the equations well formulated? since by means of different solver and sweeps and the results are not similar to those in the literature.
• When i decresed the steptime the matrix tourn singular.
• I make an arrangement for the diffusive term corresponding to $ \lambda_{zu} $.
• The limits for Dirichelt's conditions are not being respected in the resolution of the problems.

I apologize for the length of the PDF, but it was printing a lot results for analysy what happend .

[guyer]

I apologize for the length of the PDF, but it was printing a lot results for analysy what happend .

You weren't kidding about length: "14.23 × 949.1 inches"!

I'll take a look and get back to you shortly.

[guyer]

I've skimmed through the chapter you cited and looked at Fig. 3.4 and I'm not sure I understand the geometry. I think the air is modeled in 1D along the x-axis and the desiccant is modeled perpendicular to that, in 1D along the z-axis?

You've written your equations as though and are equivalent (likewise and ). This can't be true. If I understand Fig 3.4 correctly, then air and desiccant are only in contact at ().

Letting is clear enough for Eq. (3.39) (similarly, for Eq. (3.40)). That takes care of coupling the desiccant solution to the air equations.

The coupling in the other direction, from the air solution to the desiccant equations happens via the Robin conditions (3.35) and (3.36), but I don't understand where and are supposed to be evaluated. , but . I would think the Robin condition at the inlet would be different from the outlet. Ahah, that is precisely the point of Fig. 3.9. To get a solution at the inlet, then use in (3.35); to get a solution at the outlet, then use instead.

FiPy was never designed to solve different problems on different meshes. It can be done, but it certainly doesn't support anything like your equations that have variables from two different meshes on them. The result from one mesh would need to be extracted and used as the boundary condition for the other mesh. Updates probably have to be manual (FiPy's lazy evaluation isn't likely to work).

Note: w isn't coupled to anything else. It looks from the reference that Eq. (3.25) has already been substituted into (3.27) and (3.30). If you want to know about water, you could independently integrate (3.25) as an ODE, and FiPy might even do a reasonable job of it, but I would leave it out for now.

[AitorCendoya]

I've skimmed through the chapter you cited and looked at Fig. 3.4 and I'm not sure I understand the geometry. I think the air is modeled in 1D along the x-axis and the desiccant is modeled perpendicular to that, in 1D along the z-axis?

This is correct.

You've written your equations as though image and image are equivalent (likewise image and image). This can't be true. If I understand Fig 3.4 correctly, then air and desiccant are only in contact at image (image).

You're absolutely right, I didn't appreciate it this way.

FiPy was never designed to solve different problems on different meshes. It can be done, but it certainly doesn't support anything like your equations that have variables from two different meshes on them.

What if I put all the variables in a 2D cell grid, making the dimensionless coefficients a vector, corresponding to $ x^{}$ and $ z^{}$?

The result from one mesh would need to be extracted and used as the boundary condition for the other mesh. Updates probably have to be manual (FiPy's lazy evaluation isn't likely to work).

I will try to model with a for between the sweep equation.

Note: w isn't coupled to anything else. It looks from the reference that Eq. (3.25) has already been substituted into (3.27) and (3.30). If you want to know about water, you could independently integrate (3.25) as an ODE, and FiPy might even do a reasonable job of it, but I would leave it out for now.

w is required to calcualted the $ c_{tot}$ (3.9) which is substituted in the dimensionless parameters (3.28 and 3.29)

Jhonatan, thanks for yours answer. I will try to fix the approch of the equations and posted here

[guyer]

What if I put all the variables in a 2D cell grid, making the dimensionless coefficients a vector, corresponding to $ x^{}$ and $ z^{}$?

You could certainly do that, but it would be much more expensive to solve. If, ultimately, you're interested in solutions at any point in (x,z), then that's definitely what to do, but if inlet and outlet are sufficient, it's much more efficient to solve the 1D problems and transfer values at the boundaries.

[AitorCendoya]

You could certainly do that, but it would be much more expensive to solve. If, ultimately, you're interested in solutions at any point in (x,z), then that's definitely what to do, but if inlet and outlet are sufficient, it's much more efficient to solve the 1D problems and transfer values at the boundaries.

I am trying to do it this way, but the results seem strange. Attached the pdf of the jupyter nootebok, the behavior is more reasonable, but the values it takes are excessively high. GitHub_2.pdf

The scipt is

It s god the approch of equation for knows of the air state and inlet and outlet?

[guyer]

I've posted a revised script that get's closer, but I suspect is still not right.

The major changes I made are:

• Reimplemented Robin condition as in our documentation.
• Applied Dirichlet conditions to the air solutions, as per your reference.
• Changed eq3 and eq4 to use UpwindConvectionTerm, since they are purely convective (note that FiPy is not great at purely hyperbolic equations).
• Revised the sweep/solve loop to appropriately communicate between the desiccant and air solutions.

I suggest you check over all of your parameters. It seems like some of the boundary conditions are getting swamped and it could be because coefficients in the equations are the wrong order of magnitude. For example, Table 3.1 says that N = 0.5 rpm, but I can only get D_vzu = 4058 if N = 0.5 s**-1. I don't know where \rho_a is documented; if I use \rho_a = 1.225, I get c_2 to be close to, but not the same as, your script. In general, I recommend you calculate your coefficients in the script and don't just put in random numbers; it'll be easier to troubleshoot any errors.

[guyer]

What did you base your implementation of the Robin conditions on? It looks like something we plausibly could have written, but it's not how we document to write it now. I want to make sure that we're not putting out any misguided code.

[AitorCendoya]

Table 3.1 says that N = 0.5 rpm, but I can only get D_vzu = 4058 if N = 0.5 s**-1.

In the references. While explained the equations, comment what is the dimensional unit for each variable. Almost all coefficients, depend on N, also in dimensionless time \tau = \frac{tN}{60}, which is supposed to be dimensionless, below Eq (2.22) this N dimension.

suggest you check over all of your parameters. It seems like some of the boundary conditions are getting swamped and it could be because coefficients in the equations are the wrong order of magnitude.

Thank you very much, the inlet velocity and hydrodynamic diameter was wrong.

• The revised script, model the Adsorption and Regeneration zone at the sime time? why not appear a inflexion point in the graphic where the slope changes ( at tau = 0.5).
• In the Dirichilet boundaries condition, we limit the limit of the variable, but in the graphs this does not happen, what is the reason?.

  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)

I put all variables in the script and made the calculus corresponding for some parameters. In this chapter, appear the Table for Nusslet , Sherwood Number and internal mass transfer coefficient,

• I apologize for the mistak interpretation for p = H * (0.0090000000000000) + 0.012 # Turn dimensional again Eq.(3.20) RH = p*(fp.numerix.exp(5294/(T_ci+273))/1e6) # Relative Humidity of the air Eq.(3.16) The correct variable is Ha, i put Hax and the solutions turn divergent, because i put Ha it is not possible because the variable are in different mesh) Here is the script

[AitorCendoya]

What did you base your implementation of the Robin conditions on? It looks like something we plausibly could have written, but it's not how we document to write it now. I want to make sure that we're not putting out any misguided code.

I based in the matrixial from of the faceGrad, it is not in the FIPY document, just try to get a matrix that will represent the diffusion coefficient, except where the mask is. I did not use set.Value, because lambda_zu, already a matrix, depended on a Variable and if I applied it when solving, it gave me a dimension error in the shape of associated DiffusionTerm lambda_zu.(In the first-second version of scipt)

[guyer]

Table 3.1 says that N = 0.5 rpm, but I can only get D_vzu = 4058 if N = 0.5 s**-1.

In the references. While explained the equations, comment what is the dimensional unit for each variable. Almost all coefficients, depend on N, also in dimensionless time \tau = \frac{tN}{60}, which is supposed to be dimensionless, below Eq (2.22) this N dimension.

Ah, I see. Putting on my NIST hat, that's a misguided way for them to have written it. N = 0.5 rpm is identical to N = (1/120) s**-1. You shouldn't change equations when you change units. I understand it's what you have to work with.

[guyer]

I based in the matrixial from of the faceGrad, it is not in the FIPY document, just try to get a matrix that will represent the diffusion coefficient, except where the mask is.

Oh, OK. It was a good effort.

[guyer]

• In the Dirichilet boundaries condition, we limit the limit of the variable, but in the graphs this does not happen, what is the reason?.

Because Eqs. (3.39) and (3.40) are purely convective, there is nothing in these equations to "feel" the constraint at the outlet boundary.

If you remove the constraints at x*=1, you get a strange jump in Ta and Ha at the boundary, which I believe is pile-up to satisfy the default no-flux boundary condition. I'm not sure about this, though; there's no diffusion in the equation, so I don't know why counter-diffusion would be introduced at the boundary.

Ultimately, I think it's physical that these boundary conditions aren't reflected in the solution. The outlet temperature and humidity will be determined by the exchange with the desiccant along the flow channel, not by how dry or how warm you want the air to be when it comes out.

[AitorCendoya]

I try to approach a basic model that has more hypotheses, based on. Also, to simplify further, consider that the storage terms for air are negligible.( $ \frac{\partial Y}{\partial t}$ and $ \frac{\partial T}{\partial t}$).

In the first case I try to solve the equations, a 2 PDE and 2 EDO, but I get the "RuntimeError: Factor is exactly singular", based on the following forum. I try to put all variable and SourceTerms in the Implicit Fom, but still giveme the error.

# Parameteres listed on the Table 1.
p_atm = 101325                                       # Pa,  air pressure
T_reg = 100+273                                      # K, Regeneration temperature
T_ad = 34.3+273                                      # K, Adsorption temperature
RH = 56.2                                            # %, relative humidity of the adsorption air
omega_ad = 0.019                                     # kgw/kgair, humidity in adsorption inlet
V_ad = V_reg = 2.                                    # m s^-1, velocity of the adsorption air
t0 = 256                                             # K
p0 = 0.98e5                                          # Pa
D = 2.302e-5*(p0/p_atm)*(T_ad/t0)**(1.81)            # m s^-2, Mass difussion of the vapor in the air
#D_va = 2.64e-5 
c_pd = 921                                           # J kg^-1 K^-1, Specific heat of silica gel 
c_pm = 900                                           # J kg^-1 K^-1, Specific heat of matrix
# Property of the air
lambda_a = 0.0321                                    # W m^-1 K^-1, Thermal conductivity
c_pa = 1009                                          # J kg^-1 K^-1, Specific heat of the air
c_pv = 2028.1                                        # J kg^.1 K^-1, Specific heat of the vapor 
c_pl = 4179                                          # J kg^-1 K^-1, Specific heat of the water
rho_a = 1.2754                                       # kg m^-3
h_v = 2358*1000                                      # J kg^-1, Laten Heat of evaporation
# Geometry
a = 0.0015                                           # m, Hight of the duct
b = 0.0015                                           # m, Width of the duct
gamma = 2                                            # aspect ratio
Nu = 2.45                                            # Nusslet number
Sh = 2.45                                            # Sherwood number
L = 0.2                                              # m, Thickness of desiccant
r_0 = 0.2                                            # m, Radius of the desiccant wheel
omega = 0.0123                                       # rad s^-1, Rotational speed
f_d = 0.005                                          # kg m^-1, Specific desiccant mass 
f_m = 0.003                                          # kg m^-1, Specific matrix mass
#########################################################################################################
from fipy import (numerix,Grid1D, Viewer, LinearLUSolver, CellVariable, FaceVariable, ImplicitSourceTerm, TransientTerm, UpwindConvectionTerm)

n = 100
mesh = Grid1D(nx = n, Lx = L)

# Call of Variables
Y = CellVariable(mesh = mesh, name = 'Ha', value = omega_ad, hasOld = True)
W = CellVariable(mesh = mesh, name = 'Wd', value = 0.0, hasOld = True)
T = CellVariable(mesh = mesh, name = 'Ta', value = T_ad, hasOld = True)
T_as = CellVariable(mesh = mesh, name = 'Td', value = 0.0, hasOld = True)

A = 2*a*b                                                                                            # m^2, Cross sectional area of the duct
P = 2*b+2*numerix.sqrt(b**2+(a*numerix.pi)**2)*(3+(2*b/a*numerix.pi)**2)/(4+(2*b/a*numerix.pi)**2)   # Perimeter of the duct

# Transient Parameters     
Q = h_v*(1.0+0.284*numerix.exp(-10.28*W))                                        # J kg^-1, Adsorption Heat Eq.(8)
phi_w = 0.0078-0.05759*W+24.16554*W**2-124.78*W**3+204.226*W**4                  # Relative Humidity in silica gel Eq.(9)
p_ws = numerix.exp(23.196-(3816.44/(T_as-46.13)))                                # Pa, pressure of satued water vapor Eq.(10)
Y_w = (0.62188*phi_w)/(p_atm/p_ws-phi_w)                                         # Humidity ratio of air in equilibrium desiccant Eq.(11)
alpha = (Nu*lambda_a*P)/(4*A)                                                    # W m^2 K^-1, Heat transfer coefficient Eq.(12)
Ky = rho_a*Sh*D*P/(4*A)                                                          # kg m^-2 s^-1, Mass transfer coefficient Eq.(12)

# Boundary Conditions
Y.constrain(omega_ad, where = mesh.facesLeft)
T.constrain(T_ad, where = mesh.facesLeft)

# Coefficients for the equations
aux_1 = A*rho_a*(c_pa + c_pv*Y)
aux_2 = f_d * (c_pd + c_pl*W) +f_m * c_pm
omega_1 = f_d/(2*A*rho_a)                  # -
omega_2 = 2*Ky*P/f_d                       # s^-1
omega_3 = (aux_2)/(2*aux_1)                # -
omega_4 = (Ky*P*Q)/(aux_1)                 # m K
omega_5 = (2*alpha*P)/(aux_2)              # s^-1
omega_6 = (2*Ky*P*Q)/(aux_2)               # K s^-1
omega_7 = (2*Ky*P*c_pv)/(aux_2 )           # s^-1

# Sources Terms
S02 = omega_2*Y_w
S22 = ImplicitSourceTerm(coeff = omega_2, var = Y)
S2 = S02-S22
S03 = omega_4*Y_w
S33 = ImplicitSourceTerm(coeff = omega_4, var = Y)
S3 = S33-S03
S04 = omega_6*Y_w 
S44 = (ImplicitSourceTerm(coeff = omega_5, var = T_as) - ImplicitSourceTerm(coeff = omega_5, var = T)
       - ImplicitSourceTerm(coeff = omega_6, var = Y) + ImplicitSourceTerm(coeff = omega_7*Y_w, var = T)
       - ImplicitSourceTerm(coeff = omega_7*Y_w, var = T_as) - ImplicitSourceTerm(coeff = omega_7*Y, var = T)
       + ImplicitSourceTerm(coeff = omega_7*Y, var = T_as))
S4 = S04+S44

# Convective coefficient
V_ad1 = FaceVariable(mesh = mesh, value = [V_ad])

# Equations 

eq1 = TransientTerm(coeff = omega_1, var = W) + UpwindConvectionTerm(coeff = V_ad1, var = Y) == 0
eq2 = TransientTerm(coeff = 1, var = W) + S2 == 0
eq3 = TransientTerm(coeff = omega_3, var = T_as) + UpwindConvectionTerm(coeff = V_ad1, var = T) == S3
eq4 = TransientTerm(coeff = 1, var = T_as) + S4 == 0

eq1234 = eq1 & eq2 & eq3 & eq4

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

vT = Viewer(vars=(T, T_as))
vH = Viewer(vars=(Y, W))

for step in range(steps):
    Y.updateOld()
    W.updateOld()
    T.updateOld()
    T_as.updateOld()

    print("Humidity")
    res = 1e+10
    for sweep in range(sweeps):
        res = eq1234.sweep(dt=dt)
        print(res)

    vT.plot()
    vH.plot()

[guyer]

I think I've gotten your previous script to work. The major changes are:

• Made sure that Gamma_eq1 updates when lambda_zu changes. Previously it was held at the initial value.
• Introduced the water evolution equation.

Solutions are stable and make some kind of sense, although they don't agree with Fig. 3.9, etc. I believe that's because the initial conditions aren't the same.

[AitorCendoya]

Solutions are stable and make some kind of sense, although they don't agree with Fig. 3.9, etc. I believe that's because the initial conditions aren't the same.

I have changed the input conditions of the variables. I think that in the variables T and H, they must be less than 0, because the dimensionless of $\omega^{*} values are between 0 and 1 ($\omega{ci} and \omega{hi}$).

I place the input conditions and Dirichlet boundary conditions for the adsorption zone.

T = fp.CellVariable(mesh=mesh1, name ='T', value= (-T_ci/(T_hi-T_ci)), hasOld = True)  
H = fp.CellVariable(mesh=mesh1, name ='H', value=(-omega_ci/(omega_hi-omega_ci)), hasOld = True)  
water = fp.CellVariable(mesh=mesh1, name ='w', hasOld = True)

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

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

The behavior of the curve is accordingly, but the values are lower than normal.

However, when I set the input value for the regeneration zone, the behavior of the curve is inappropriate. The Ha should increase, exceeding the input condition(Ha =1), and it does not happen.

T = fp.CellVariable(mesh=mesh1, name ='T', value= ((T_ci+T_hi)/2-T_ci/(T_hi-T_ci)), hasOld = True)   
H = fp.CellVariable(mesh=mesh1, name ='H', value=((W_max/Kp)-omega_ci/(omega_hi-omega_ci)), hasOld = True)  
water = fp.CellVariable(mesh=mesh1, name ='w',value= 1., hasOld = True)

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

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

The values of water uptake in the solid are low, this term must be the one that is causing the other values to be low. I will check again, all the parameters and the interaction between them and comment to you. Because there may be discrepancies between the values of $\omega{min} and \omega{max}$.