Re: FiPy Heat Transfer Solution
Thanks Daniel. I think I understand now. What would I do to multiply through by r in writing a FiPy equation? Should I some how create a spatial variable for r? Such as r= Variable? Do I have to somehow give it a range of possible values? (I didn't see an example for this online) Then I could use a standard mesh, correct? mesh = Grid1D(nx=nx,dx=dx) and then the equation becomes? eqT = TransientTerm(r * rho_b * Cp_b*1000/MW_b, var=T) == DiffusionTerm(coeff = L_eff, var=T) + r*S Sorry to have so many questions. On Thu, Aug 30, 2018 at 10:20 AM Daniel Wheeler wrote: > There are two ways to solve on a cylindrical domain in FiPy. You can > either use the standard diffusion equation in Cartesian coordinates (2nd > equation below) and with a mesh that is actually cylindrical in shape or > you can use the diffusion equation formulated on a cylindrical coordinate > system (1st equation below) and use a standard 2D / 1D grid mesh. I think > you're confused about "multiplying through by R". The problem with the > first equation in FiPy (cylindrical) is that FiPy can't have a coefficient > outside the first derivative, which occurs in the diffusion term of the > first equation. Of course that isn't an issue for constant coefficient like > density etc., but r, however, is not constant. We can get around that > problem by multiplying the first equation by r (small r) because the r can > come inside the time derivative in the transient term and just be a > multiplier on the source term, which is fine. Your equation then looks like > this. > > > Hope that helps. > > Cheers, > > Daniel > > On Wed, Aug 29, 2018 at 9:13 AM Daniel DeSantis > wrote: > >> Daniel, >> >> Thanks again for getting back to me. Sorry to keep coming back to the >> same problem. It does seem like there's a problem with the >> CylindricalGrid1D mesh. I noticed that running the most basic 1D diffusion >> example with a cylindrical grid doesn't seem to work correctly. With that >> said, I hope I can impose on you one more time and ask how I should convert >> a cylindrical coordinate PDE into a rectangular coordinate PDE for FiPy, >> essentially turning this >> [image: image.png] >> into this: >> [image: image.png]. >> Previously you said to multiply through by r. Assuming R is the maximum >> of the variable r, am I looking at this equation >> [image: image.png] or should I be using this? >> [image: image.png]If r should still be in the variable format, how would >> I make that work in FiPy? >> >> Thanks, >> Dan DeSantis >> > > > -- > 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
There are two ways to solve on a cylindrical domain in FiPy. You can either use the standard diffusion equation in Cartesian coordinates (2nd equation below) and with a mesh that is actually cylindrical in shape or you can use the diffusion equation formulated on a cylindrical coordinate system (1st equation below) and use a standard 2D / 1D grid mesh. I think you're confused about "multiplying through by R". The problem with the first equation in FiPy (cylindrical) is that FiPy can't have a coefficient outside the first derivative, which occurs in the diffusion term of the first equation. Of course that isn't an issue for constant coefficient like density etc., but r, however, is not constant. We can get around that problem by multiplying the first equation by r (small r) because the r can come inside the time derivative in the transient term and just be a multiplier on the source term, which is fine. Your equation then looks like this. Hope that helps. Cheers, Daniel On Wed, Aug 29, 2018 at 9:13 AM Daniel DeSantis wrote: > Daniel, > > Thanks again for getting back to me. Sorry to keep coming back to the same > problem. It does seem like there's a problem with the CylindricalGrid1D > mesh. I noticed that running the most basic 1D diffusion example with a > cylindrical grid doesn't seem to work correctly. With that said, I hope I > can impose on you one more time and ask how I should convert a cylindrical > coordinate PDE into a rectangular coordinate PDE for FiPy, essentially > turning this > [image: image.png] > into this: > [image: image.png]. > Previously you said to multiply through by r. Assuming R is the maximum of > the variable r, am I looking at this equation > [image: image.png] or should I be using this? > [image: image.png]If r should still be in the variable format, how would > I make that work in FiPy? > > Thanks, > Dan DeSantis > -- 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, Thanks again for getting back to me. Sorry to keep coming back to the same problem. It does seem like there's a problem with the CylindricalGrid1D mesh. I noticed that running the most basic 1D diffusion example with a cylindrical grid doesn't seem to work correctly. With that said, I hope I can impose on you one more time and ask how I should convert a cylindrical coordinate PDE into a rectangular coordinate PDE for FiPy, essentially turning this [image: image.png] into this: [image: image.png]. Previously you said to multiply through by r. Assuming R is the maximum of the variable r, am I looking at this equation [image: image.png] or should I be using this? [image: image.png]If r should still be in the variable format, how would I make that work in FiPy? Thanks, Dan DeSantis On Thu, Aug 16, 2018 at 6:04 PM Daniel Wheeler wrote: > 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 ] > -- 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
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 ] > -- Dani
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 ]
FiPy Heat Transfer Solution
Hello, First off, I'd like to say I really like FiPy and appreciate that you have made such great tool for everyone to use. I am trying to use FiPy to solve a heat transfer problem in 1D and I have a few implementation questions. Forgive me if they seem simplistic. I have looked at a lot of the documentation and examples, but there may be some things that have eluded me. If it seems like I'm writing very simple, obvious things it's just so I can get a full grasp of what is going on. 1) It seems that most, if not all, of the examples are in Cartesian coordinates. Thus, a 1D diffusion equation with a source term would be: [image: image.png] where rho is density, Cb is specific heat, lamba is an effective heat transfer and S is a source term. Of course, T is temperature, t is time, and x is the coordinate. I would normally set this mesh up by using mesh = Grid1D(nx =100., Lx = 1., dx = 1./100.) If I wanted to run this in a cylindrical coordinate grid, the equation would become: [image: image.png] Can I just call the following in FiPy to account for this coordinate system change?: mesh = CylindricalGrid1D(nx = 100.,Lx = 1., dx = 1./100.) 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.? 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)? 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've included a copy of the code below. There are some extraneous details/code in there that I was playing with. Thank you! -- Daniel DeSantis from fipy import * # Constants rho_b_MH = 589. # kg m^-3 Cp_b = 88.25 #kJ (mol*K)^-1, # I'm using L in place of lambda L_eg = 150. # W m^-1 K^-1 L_MH = 1.74 # W m^-1 K^-1 L_eff = 8.4 # W m^-1 K^-1 t_r = 3.7 #min M = 2.02 #g/mol dH = 28.*1000. #kJ/mol h = 150. # W/m^-2 K^-1 T_c = 60 # Coolant Temperature Equation Set up and Boundary Condition R = 1. # There was no radius provided. I just set it to 1 for now nx = 100. # arbitrary number of points dx= R/nx #mesh = Grid1D(nx = nx, dx=dx) # Create a "mesh" for the 1D array mesh = CylindricalGrid1D(nx=nx,dx=dx,Lx=R) T = CellVariable(name ="Temperature",mesh=mesh,value = 60.) #T.constrain(60.,mesh.facesLeft) #T.constrain(120.,mesh.facesRight) T.faceGrad.constrain(0,mesh.facesRight) T.faceGrad.constrain(h*(T.faceValue-T_c),mesh.facesLeft) eq1 = TransientTerm(coeff=rho_b_MH*Cp_b,var=T) == DiffusionTerm(coeff = L_eff, var = T) + rho_b_MH*(dH/M) # Equation Solution vi = Viewer(T,datamin=0,datamax=100) TSD = .0001 #min/step steps = 75 #steps for step in range(steps): eq1.solve(var=T,dt=TSD) vi.plot() ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]