Re: FiPy Heat Transfer Solution
Apologies for the slow reply. I hadn't noticed your response until now. See below. On Wed, Aug 8, 2018 at 1:09 PM Daniel DeSantis wrote: > > Daniel, > > Thank you very much for your help.The code does run much better. > > I was curious about these lines of code. I'm hoping to understand what this > means so that I can use it if needed, next time. > > Thank you again, > Dan DeSantis > > constraint_value = FaceVariable(mesh=mesh) > T.faceGrad.constrain(constraint_value,mesh.facesLeft) > > for sweep in range(sweeps): > constraint_value[:] = h * ((60. + 273.15) - T.faceValue) I only used the above code as it happened to work. It should work without having to do the update in the loop, but something isn't right with the way FiPy is updating constraints in the case when the constraint is not a bare FaceVariable. In other words, we make a face variable and then say that the constraint depends on that face variable. We have to explicitly update the face variable in the loop. The expression "h * ((60. + 273.15) - T.faceValue" should also be a face variable and so we should be able to use "T.faceGrad.constrain(h * ((60. + 273.15) - T.faceValue, mesh.facesLeft)" and not update in the loop. That doesn't work, however. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: FiPy Heat Transfer Solution
Daniel, Thank you very much for your help.The code does run much better. I was curious about these lines of code. I'm hoping to understand what this means so that I can use it if needed, next time. Thank you again, Dan DeSantis constraint_value = FaceVariable(mesh=mesh) T.faceGrad.constrain(constraint_value,mesh.facesLeft) for sweep in range(sweeps): constraint_value[:] = h * ((60. + 273.15) - T.faceValue) On Wed, Aug 8, 2018 at 12:05 PM Daniel Wheeler wrote: > Hi Daniel, > > I updated your code a bit and it seems to do better. See > https://pastebin.com/51xfqWH3 for the full version. > > I changed the "eqX" equation to > > param = k_c * (P - Peq) > maxX = 0.9 > largeValue = 1e9 > constraint = largeValue * (X > maxX) > big_value = 1e9 > eqX = TransientTerm(var = X) == param - ImplicitSourceTerm(param, > var=X) + maxX * constraint - ImplicitSourceTerm(constraint, var=X) > > to have a constraint on the maximum value of X. I'm not sure whether > this is physical, but it does constrain the max value of X. > > I changed the "eqY" equation to be, > > eqT = TransientTerm(rho_b_MH * Cp_b, var=T) == DiffusionTerm(coeff > = L_eff, var=T) + rho_b_MH*(dH/M)*rc > > to have the transient coefficient inside the term. Again, not sure if > that matters, but did it anyway. > > I changed the constraint to be > > constraint_value = FaceVariable(mesh=mesh) > T.faceGrad.constrain(constraint_value,mesh.facesLeft) > > and then updated the constraint in the loop with > > for sweep in range(sweeps): > constraint_value[:] = h * ((60. + 273.15) - T.faceValue) > > If FiPy worked correctly that wouldn't be necessary, but evidently it > doesn't. > > Anyway that all seems to work at least as far as I can tell. > > Cheers, > > Daniel > > -- > Daniel Wheeler > ___ > fipy mailing list > fipy@nist.gov > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > -- Daniel DeSantis ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: FiPy Heat Transfer Solution
Hi Daniel, I updated your code a bit and it seems to do better. See https://pastebin.com/51xfqWH3 for the full version. I changed the "eqX" equation to param = k_c * (P - Peq) maxX = 0.9 largeValue = 1e9 constraint = largeValue * (X > maxX) big_value = 1e9 eqX = TransientTerm(var = X) == param - ImplicitSourceTerm(param, var=X) + maxX * constraint - ImplicitSourceTerm(constraint, var=X) to have a constraint on the maximum value of X. I'm not sure whether this is physical, but it does constrain the max value of X. I changed the "eqY" equation to be, eqT = TransientTerm(rho_b_MH * Cp_b, var=T) == DiffusionTerm(coeff = L_eff, var=T) + rho_b_MH*(dH/M)*rc to have the transient coefficient inside the term. Again, not sure if that matters, but did it anyway. I changed the constraint to be constraint_value = FaceVariable(mesh=mesh) T.faceGrad.constrain(constraint_value,mesh.facesLeft) and then updated the constraint in the loop with for sweep in range(sweeps): constraint_value[:] = h * ((60. + 273.15) - T.faceValue) If FiPy worked correctly that wouldn't be necessary, but evidently it doesn't. Anyway that all seems to work at least as far as I can tell. Cheers, Daniel -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: FiPy Heat Transfer Solution
Hello,
First off, I appreciate your previous responses. I do want to follow up on
this question because I was running into an issue with the same/updated
code I sent before. I'm not sure it's working the way it should. I am
simultaneously solving two variables. X is maxing out around 1, which I
would expect, though I would like to constrain it to 0.9 but I can't seem
to figure out how. T, however, doesn't seem to be following constraints
very well. I am trying to use a robin constraint on both the left and right
faces but for testing purposes I've also tried to constrain the faces to
set points and that doesn't seem to work either. I would expect T to vary
along the radial direction, which is what should be showing in the viewer
plot. Do you have any advice?
Thank you,
#%% - Imports
from numpy import *
import matplotlib.pylab as plt
from fipy import *
import pandas as pd
#%% - Constants
rho_b_MH = 589. # kg m^-3
Cp_b = 181.170 #J (mol*K)^-1, Thermal conductivity and specific heat
measurements of metal hydrides
# I'm using L in place of lambda
L_eff = 8.4 # W m^-1 K^-1
E_c = 45.*1000 # J/mol
Beta = 13.5 #%
P = 100. #atm
M = 2.02 #g/mol
dH = 28.*1000. #J/mol
h = 150. # W/m^-2 K^-1
T_c = 60.+273.15 #K, Coolant Temperature
To = 60+273.15 # K, Bed Temperature
R = 8.314 # J/mol/K
Peq = ((E_c/dH)/(1+E_c/dH))*P #atm
A = 1.1319e25
#%% - Mesh
L = 1.
nx = 10.
dx = L/nx
# Create the mesh
mesh = Grid1D(nx=nx,dx=dx)
#mesh = CylindricalGrid1D(nx=nx,dx=dx)
# Sorbtion Setup #
X = CellVariable(name = "Sorbtion",mesh = mesh, value=0.1, hasOld=True)
CC = mesh.cellCenters[0]
PDE Variable Setup #
#T = CellVariable(name = "Temperature", mesh = mesh,value = 60.+273.15)
T = CellVariable(name = "Temperature", mesh = mesh, value = 60.+273.15,
hasOld=True)
# I want to use the robin boundary condition but in either case, the
constraint doesn't seem to be working
#T.constrain(60.+273.15,mesh.facesLeft)
T.faceGrad.constrain(h*((60.+273.15)-T.faceValue),mesh.facesLeft)
# I want to use the robin boundary condition
#T.constrain(110.+273,mesh.facesRight)
T.faceGrad.constrain(0.,mesh.facesRight)
k_c = A*T**-9.73625
rc = k_c*(1-X)*(P-Peq)
eqX = TransientTerm(var = X) == k_c*(1.-X)*(P-Peq)
eqT = rho_b_MH * Cp_b * TransientTerm(var = T) == DiffusionTerm(var =
T,coeff = L_eff) + rho_b_MH*(dH/M)*rc
eq = eqT & eqX
#%% - Equation Solution
TSD = .005 # This seems to be a stable time step for this function
sweeps = 150 # identify the number of sweeps to run, running 150 for now
ts = np.linspace(0,TSD*sweeps,sweeps+1) # ts is the number of time slots
plotted for the ODE. create an array to plot for the ode
solver = LinearLUSolver() # Trying a different solver
#Prepare an array for storing calculated values of X
Xsolution = np.zeros((sweeps+1,int(nx))) #Create a matrix a row for initial
value + 1 row for each sweep. There will be nx columns per row.
#Store the initial value of X at all points in the first row of sweeps
Xsolution[0,:] = X.value
#Prepare an array for storing calculated values of T
Tsolution = np.zeros((sweeps+1,int(nx))) #Create a matrix a row for initial
value + 1 row for each sweep. There will be nx columns per row.
#Store the initial value of T at all points in the first row of sweeps
# This matrix has the form with the rows are T|t=0...T|t=t_max, and columns
T|r=0T|r=R
Tsolution[0,:] = T.value
viewer = Viewer(vars = T, datamin = 330., datamax = 450., title =
"Temperature")
viewer2 = Viewer(vars = X, datamin = 0., datamax = 1.05, title = "X")
for sweep in range(sweeps):
X.updateOld()
T.updateOld()
eq.sweep(dt=TSD,solver=solver)
viewer.plot()
viewer2.plot()
Xsolution[sweep+1] = X.value #Store the recently computed X value in
the current sweep
Tsolution[sweep+1] = T.value #Store the recently computed T value in
the current sweep
#print("X values",X.value[-1])
#print("T values",T.value[-1])
#print(rc[-1])
df=pd.DataFrame([X.value[-1],T.value[-1],k_c[-1],rc[-1]],index =
["X","T","k_c","rc"],columns=[""])
print(df.T)
#print(Xsolution[:,0])
#plt.plot(CC.value,X.value,'ro') # plot the points where calculations are
made in the fipy solver.
fig,ax = plt.subplots()
ax.plot(ts,Xsolution[:,0],'bo')
ax.set_xlabel("time")
ax.set_ylabel(X.name)
ax.set_title("ODE Solution")
fig,ax = plt.subplots()
ax.plot(ts,Tsolution[:,0],"bo")
ax.set_xlabel("time")
ax.set_ylabel(T.name)
ax.set_title("Time vs. Temp")
"""
fig,ax = plt.subplots()
ax.plot(ts,Tsolution[-1,:],"bo")
ax.set_xlabel("time")
ax.set_ylabel(T.name)
ax.set_title("Time vs. Temp")
"""
On Tue, Jul 24, 2018 at 2:37 PM Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:
> > On Jul 24, 2018, at 9:12 AM, Daniel Wheeler
> wrote:
> >
> > Jon, is that correct?
>
> I honestly don't know
>
>
> ___
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
--
Re: FiPy Heat Transfer Solution
> On Jul 24, 2018, at 9:12 AM, Daniel Wheeler wrote: > > Jon, is that correct? I honestly don't know ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: FiPy Heat Transfer Solution
Thank you very much! On Tue, Jul 24, 2018 at 10:13 AM Daniel Wheeler wrote: > On Wed, Jul 18, 2018 at 4:36 PM, Daniel DeSantis > wrote: > > > > > > If I wanted to run this in a cylindrical coordinate grid, the equation > would become: > > > > > > Can I just call the following in FiPy to account for this coordinate > system change?: > > > > mesh = CylindricalGrid1D(nx = 100.,Lx = 1., dx = 1./100.) > > Yes, that should work. There is a bug (see > https://github.com/usnistgov/fipy/issues/547) when calculating > gradients, but I don't think that applies to diffusion terms. > > Jon, is that correct? > > > If not, how do I account for the dependent variable when it is not in > the differential? Essentially, how do can I put lambda/r in the coefficient > etc.? > > You can multiply through by r and then split the transient term into > an implicit TransientTerm and a source term. It gives an extra source > term that is dependent on the time step. > > > 2) Is the correct way to add a source term simply to write something > like: > > > > eq = TransientTerm() == DiffusionTerm() + S > > > > or > > > > eq = TransientTerm() == DiffusionTerm() + ImplicitSourceTerm(S)? > > Either one. ImplicitSourceTerm implies that the source term is "var * > S", not "S" only, but the "var" is the variable that you're solving > for. To use ImplicitSourceTerm the S should be negative. In your case, > the source does not depend on temperature so you can't use an > ImplicitSourceTerm. > > > 3) Regarding boundary conditions, I am trying to use a boundary > condition that has my dependent variable (T) in the boundary condition. The > condition is dT/dr = h(T - Tc), where Tc is a constant. Would this be the > correct way to do that? I'm mostly inquiring about writing the dependent > variable as T.faceValue... > > > > T.faceGrad.constrain(h*(T.faceValue - Tc), mesh.facesLeft) > > I think that should work. Test it carefully. It's never obvious to me > what the sign should be. Don't trust it though. > > -- > Daniel Wheeler > > ___ > fipy mailing list > fipy@nist.gov > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > -- Daniel DeSantis ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: FiPy Heat Transfer Solution
On Wed, Jul 18, 2018 at 4:36 PM, Daniel DeSantis wrote: > > > If I wanted to run this in a cylindrical coordinate grid, the equation would > become: > > > Can I just call the following in FiPy to account for this coordinate system > change?: > > mesh = CylindricalGrid1D(nx = 100.,Lx = 1., dx = 1./100.) Yes, that should work. There is a bug (see https://github.com/usnistgov/fipy/issues/547) when calculating gradients, but I don't think that applies to diffusion terms. Jon, is that correct? > If not, how do I account for the dependent variable when it is not in the > differential? Essentially, how do can I put lambda/r in the coefficient etc.? You can multiply through by r and then split the transient term into an implicit TransientTerm and a source term. It gives an extra source term that is dependent on the time step. > 2) Is the correct way to add a source term simply to write something like: > > eq = TransientTerm() == DiffusionTerm() + S > > or > > eq = TransientTerm() == DiffusionTerm() + ImplicitSourceTerm(S)? Either one. ImplicitSourceTerm implies that the source term is "var * S", not "S" only, but the "var" is the variable that you're solving for. To use ImplicitSourceTerm the S should be negative. In your case, the source does not depend on temperature so you can't use an ImplicitSourceTerm. > 3) Regarding boundary conditions, I am trying to use a boundary condition > that has my dependent variable (T) in the boundary condition. The condition > is dT/dr = h(T - Tc), where Tc is a constant. Would this be the correct way > to do that? I'm mostly inquiring about writing the dependent variable as > T.faceValue... > > T.faceGrad.constrain(h*(T.faceValue - Tc), mesh.facesLeft) I think that should work. Test it carefully. It's never obvious to me what the sign should be. Don't trust it though. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
