Jonathan,
I have just a few follow up questions. Namely, regarding to working within
cylindrical coordinates. If I switch the mesh to cylindrical, it seems the
temperature profile never changes, unless I go to extremely small time
steps. That strikes me as odd. What do you think?
Also, the shifting profile seems to reverse, with the left boundary
condition dropping to ~300K over time while in Cartesian coordinates, the
temperature rises to 600K over time (a result I would expect). I've
included a copy of the work with the changes you've suggested previously.
If you look at lines 23/24, you can toggle the cylindrical grid on and off,
and in lines 30-34, I've included options for the time steps that work for
each respective coordinate system. Do you have any suggestions?
Finally, in reference to the velocity I would like to incorporate, I
understand the difference between Cell Variables and Face Variables, but
why is the rank of the velocity -1? What does the Face Variable rank
indicate?
Thank you,
Dan
On Thu, Apr 11, 2019 at 8:14 AM Guyer, Jonathan E. Dr. (Fed) via fipy <
fipy@nist.gov> wrote:
>
>
> > On Apr 10, 2019, at 4:42 PM, Daniel DeSantis <desan...@gmail.com> wrote:
> >
> > 1) Set a different Cp, rho, and k in just the PCM domain when the
> temperature in that domain crosses over a certain melt temperature (350K).
> I tried an if statement inside the solving loop and I think that has
> worked, based on some print testing I did in the loop. I just wanted to get
> your opinion and see if that was the right way to go.
>
> Seems reasonable, as long as you're calling `.setValue()` and not
> redefining Cp, rho, and k.
>
> In fact, I wouldn't use an `if`:
>
> >>> for step in range(steps):
> ... T.updateOld()
> ... eq.solve(dt=dt)
> ... rho.setValue(rho_melt, where=PCM & (T > 350))
>
>
> > 2) I'd like to move this to cylindrical coordinates with the lengths
> being radii. I'm not sure if the cylindrical grid is working though. And if
> it isn't, is there a way to convert the diffusion equation into something
> that would solve the problem in Cartesian? More specifically, how would you
> write an appropriate equation for the heat transfer of a radial system in
> Cartesian coordinates? In short, how would you write this:
> > <image.png>
> > in FiPy?
>
> The CylindricalGrids work fine, except for taking the gradient of the
> field, which is an unusual need.
>
> >
> > 3) Ideally, I would move this to a 2D system which is the most confusing
> to me. In a 2D system, heat would be transmitted radially while air would
> be flowing axially. The air would cool as it passes through the tube. The
> diffusion term in the axial direction would be significantly lower than the
> convection term. Can I use the same heat transfer equation you've suggested
> slightly modified? I would guess the equation to be something like this:
> >
> > eq = TransientTerm(coeff = rho*Cp,var = T) +
> PowerLawConvectionTerm(coeff=rho*Cp*v,var = T) == DiffusionTerm(coeff=k,var
> = T)
> >
> > s.t. v = velocity of the air.
>
> Looks OK
>
> > I would also think that I would have to set the value of v to 0 at the
> wall and beyond, which means v would be a CellVariable like rho and k.
>
> Actually, v should be a rank-1 FaceVariable.
>
> Anisotropic diffusion is supported in FiPy (see
> https://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.anisotropy.html),
> but are you really sure you have lower axial diffusivity? It seems to me
> more likely that transport is simply dominated by axial flow.
>
> > I also guess this equation would be subject to any changes that would be
> made in question 2, regarding writing a cylindrical equation in a cartesian
> system. What do you think?
>
> Use a CylindricalGrid and don't change the equation at all.
>
> > Again, thank you so much for your help and I hope I haven't taken up too
> much of your time.
>
> Happy to help.
>
>
> _______________________________________________
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
--
Daniel DeSantis
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 2 15:14:40 2019
@author: ddesantis
"""
from fipy import *
import numpy as np
import matplotlib.pyplot as plt
#%% - Mesh
nx = 100
L = 0.061 #m
Ly = 1 #m
L1 = 0.035 # m, Length of air channel
t_w = 1./1000. # m, wall thickness
L2 = L1+t_w # m, interface between PCM and wall
L3 = L # m, Total Length
mesh = Grid1D(Lx=L,nx=nx)
#mesh = CylindricalGrid1D(Lr = L,nr=nx)
x = mesh.cellCenters[0]
#%% - Steps
# Cartesian coordinates time step
dt = 0.25 #s,time step for sweeps
# Cylindrical coordinates time step
#dt = 0.00005 #s,time step for sweeps
steps = 150
#%% Universal Constants
MW = 2.02 #g H2 mol^-1
R = 8.314e-3 # kj mol^-1 K^-1
#%% - air Properties
T_in = 600. # K
P_in = 1. # bar
v = 0.5 # m s^-2
k_air = 44.41/1000./1000. # kJ s^-1 m^-1 K^-1 (air at ~600K)
Cp_air = 1.051 # kJ kg ^-1 K^-1
rho_air = 0.6 # kg m^-3, At 600K
Re = 1150. # Reynolds Number, dimensionless
Pr = 0.727 # Prandlt Number, dimensionless
Nu = 1.86 # Nusselt Number, dimensionless
#%% - Wall properties
#I'm assuming the wall is steel
rho_w = 8050. # kg m^-3
Cp_w = 0.5 # kJ K^-1 kg^-1
k_w = 16.3/1000. # kJ s^-1 m^-1 K^-1
T_w_initial = (T_in+293)/2
#%% PCM Properties
T_init = 293. #K
T_melt = 300. #K
H_fus = 206. #kJ kg^-1
k_s = 0.18 #W m^-1 K^-1
k_l = 0.19 #W m^-1 K^-1
k_PCM = k_s
Cp_s = 1.8 #kJ kg^-1 K^-1
Cp_l = 2.4 #kJ kg^-1 K^-1
Cp_PCM = Cp_s
# I assume they want solid and liquid presented the same way again?
rho_s = 789. # kg m^-3
rho_l = 750. # kg m^-3
#%% - Domain
air = x<= L1
wall = (x>L1) & (x<=L2)
PCM = x > L2
#%% Equation Variables
T = CellVariable(name = 'Temp', mesh = mesh, hasOld = True)
T.setValue(T_in, where = air)
T.setValue(T_w_initial, where = wall)
T.setValue(T_init, where = PCM)
rho = CellVariable(name = 'rho',mesh=mesh)
rho.setValue(rho_air,where = air)
rho.setValue(rho_w,where = wall)
rho.setValue(rho_s,where = PCM)
Cp = CellVariable(name = 'Cp',mesh=mesh)
Cp.setValue(Cp_air,where = air)
Cp.setValue(Cp_w,where = wall)
Cp.setValue(Cp_s,where = PCM)
k = CellVariable(name = 'k',mesh=mesh)
k.setValue(rho_air,where = air)
k.setValue(rho_w,where = wall)
k.setValue(rho_s,where = PCM)
#%% Boundary Conditions
T.constrain(T_in,mesh.facesLeft)
#%% - Equations
eq = TransientTerm(coeff = rho*Cp,var = T) == DiffusionTerm(coeff=k,var = T)
#%% - Steps
viewer = Viewer(vars=(T),datamin = 50.,datamax = 700.)
viewer.plot()
T_array = np.zeros((steps+1,nx))
T_array[0,:]=T.value
H_array = np.zeros((steps+1,nx))
for step in range(steps):
T.updateOld()
eq.solve(dt=dt)
rho.setValue(rho_l,where = PCM & (T>350) )
k.setValue(k_l,where = PCM & (T>350))
Cp.setValue(Cp_l,where = PCM & (T>350))
T_array[step+1,:] = T.value
H_array[step+1,:] = T.value*Cp.value
print(T.value)
viewer.plot()
percent_complete = float(step)/float(steps)
#print('percent complete: {:0.0f}%'.format(100.*percent_complete))
#print('step {} of {}'.format(step, steps))
#print('elapsed fill time: {:0.2f}s'.format(step*dt))
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
