Hi!
Daniel, thanks for your reply!
I had already seen CLAWPACK and also considered it as an alternative if I don't
manage to solve this problem with FiPy. Thanks for pointing it out, however I
don't want to try too many approaches at once.
I did in fact START to implement something which looks as follows
--- (snip) ---
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from fipy import Grid1D, CellVariable, FixedValue, TransientTerm,
ConvectionTerm, Viewer
# constants
# specific heat of carrier gas
c_P_G = 1.
# specific heat of vapor
c_P_V = 1.86
# specific heat of adsorbed water
c_P_W_ads = 2.2
# specific heat of dry sorbent
c_P_S = 1.
# specific surface of adsorbent
a_S = 903.
# define a mesh
length = 1.
nx = 100
dx = length/nx
mesh = Grid1D(dx=dx, nx=nx)
# define solution variables
u1 = CellVariable(mesh = mesh, hasOld = True, name = 'X_G')
u2 = CellVariable(mesh = mesh, hasOld = True, name = 'X_S')
u3 = CellVariable(mesh = mesh, hasOld = True, name = 'T_G')
u4 = CellVariable(mesh = mesh, hasOld = True, name = 'T_S')
# initial conditions: start with a dry sorbent at ambient temperature
X_G_init = 6.2e-06 # eq. water content air above sorbent with ...
# following parameters
X_S_init = 0.09 # initial water content of sorbent
T_G_init = 20. + 273.15 # initial temperature of the carrier gas in K
T_S_init = T_G_init # initial temperature of sorbent
u1.setValue(X_G_init)
u2.setValue(X_S_init)
u3.setValue(T_G_init)
u4.setValue(T_S_init)
# boundary conditions: humid carrier gas enters so that sorbent adsorbs water
X_G_in = 0.0147 # water content inc. air, corr. to dewpoint of 20°C
X_S_in = 0.3484 # eq. water content of sorbent with given carrier gas ...
# parameters
T_G_in = 20. + 273.15 # temp. of inc. carrier gas
T_S_in = T_G_in # temp. of sorbent at entrance equal to carrier gas temp.
# (use ramp when values T_S_init and T_G_in differ a lot)
u1.constrain(X_G_in, mesh.facesLeft)
u2.constrain(X_S_in, mesh.facesLeft)
u3.constrain(T_G_in, mesh.facesLeft)
u4.constrain(T_S_in, mesh.facesLeft)
# coefficients in the pde
# filling ratio of carrier gas/sorbent
eps_G = 0.4
eps_S = 1. - eps_G
# density of carrier gas/sorbent
rho_G = 1.1
rho_S = 1150.
# effective specific heat of carrier gas/sorbent
c_P_G_ast = c_P_G + X_G_in * c_P_V
c_P_S_ast = c_P_S + X_S_in * c_P_W_ads
# mass flow carrier gas per area
m_G_dot = 0.1
# water mass flow from carrier gas to sorbent per volume
def m_GS_dot(X_G):
# mass transfer coefficient
beta_GS = 0.0197
X_G_ast = X_G/1.5
return beta_GS * a_S * rho_G * (X_G - X_G_ast)
# heat flow from carrier gas to sorbent per volume
def q_GS_dot(T_G,T_S):
#heat transfer coefficient
alpha_GS = 0.0739
return alpha_GS * a_S * (T_S - T_G)
# heat of adsorption
def h_ads(X_S):
L = 2501.
return L + X_S * (-1.e4) + 6000
# partial differential equations
eqn1 = TransientTerm(coeff=1., var=u1) + \
ConvectionTerm(coeff=(m_G_dot/(eps_G*rho_G),), var=u1) == \
-m_GS_dot(u1)/(eps_G*rho_G)
eqn2 = TransientTerm(coeff=1., var=u2) == m_GS_dot(u1)/(eps_S*rho_S)
eqn3 = TransientTerm(coeff=1., var=u3) + \
ConvectionTerm(coeff=(m_G_dot/(eps_G*rho_G),), var=u3) == \
q_GS_dot(u3,u4)/(eps_G*rho_G*c_P_G_ast)
eqn4 = TransientTerm(coeff=1., var=u4) == \
-q_GS_dot(u3,u4)/(eps_S*rho_S*c_P_S_ast) + \
m_GS_dot(u1) * h_ads(u2) / (eps_S*rho_S*c_P_S_ast)
# convection term coefficient has units: kg/(s m²) / (kg/m³) = m/s so indeed is
a velocity
# the CFL number is defined as C = u * dt/dx and typically C_max = 1, so u *
dt/dx < C_max yields
# dt < dx/u
dt_CFL = dx/(m_G_dot/(eps_G*rho_G))
# solving by sweeping to 'convergence'
# for t in range(10):
# # only move forward in time after a number of sweeps is completed
# u1.updateOld()
# u2.updateOld()
# u3.updateOld()
# u4.updateOld()
# # sweep several times per time step
# res = 1e+10
# #while res > 1e-5:
# for sweep in range(3):
# res1 = eqn1.sweep(var=u1, dt=dt_CFL)
# res2 = eqn2.sweep(var=u2, dt=dt_CFL)
# res3 = eqn3.sweep(var=u3, dt=dt_CFL)
# res4 = eqn3.sweep(var=u4, dt=dt_CFL)
# res = max(res1, res2, res3, res4)
# print res
# solving 'directly'
for t in range(100):
u1.updateOld()
u2.updateOld()
u3.updateOld()
u4.updateOld()
eqn1.solve(dt=dt_CFL)
eqn2.solve(dt=dt_CFL)
eqn3.solve(dt=dt_CFL)
eqn4.solve(dt=dt_CFL)
print u2
if __name__ == '__main__':
viewer = Viewer(vars = (u2),datamin=0., datamax=0.4)
viewer.plot()
--- (snap) ---
You realize that I tried two different ways of finding a solution. I assume
that using the sweep method is what you recommended in your answer to my first
post. However, when I run the script with the call to sweep, the iPython
notebook kernel dies immediately. Also from the command line version of iPython
I have no luck but get the following error message
/home/markus/.venv/fipy/local/lib/python2.7/site-packages/fipy/solvers/pysparse/linearLUSolver.py:84:
RuntimeWarning: divide by zero encountered in double_scalars
L = L * (1 / maxdiag)
/home/markus/.venv/fipy/local/lib/python2.7/site-packages/fipy/solvers/pysparse/linearLUSolver.py:85:
RuntimeWarning: divide by zero encountered in double_scalars
b = b * (1 / maxdiag)
/home/markus/.venv/fipy/local/lib/python2.7/site-packages/fipy/solvers/pysparse/linearLUSolver.py:85:
RuntimeWarning: invalid value encountered in multiply
b = b * (1 / maxdiag)
Speicherzugriffsfehler (Speicherabzug geschrieben)
(i.e.: memory access violation or segmentation faults)
If the lines using the direct call to solve are active, the script runs.
However, comparing the solution to what I get with the MATLAB version set to
the same input parameters there are no similarities whatsoever. This might be
ok, since you already suggested that the solutions have to be swept to
convergence but what puzzles me is that looking at e.g. u2, the cell value at
the left boundary is not equal to 0.3484 as I would expect it from the
constraint. Granted, the correct solution isn't found, yet, but why is the
constraint not reflected at the left boundary of the grid? What am I not
getting?
Thanks for your time!
Markus
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
