Hi,
I am using FiPy to determine the depth of heat penetration into concrete
structures due to fire over a certain period of time. I am solving the unsteady
heat equation on a 1D grid, and modelling various scenarios, e.g.
time-dependent temperature boundary condition, temperature-dependent diffusion
coefficient. For these cases, the model compares very well to results from
other solvers (for example the commercial finite element solver TEMP/W).
However I am having trouble modelling problems with a spatially-varying
diffusion coefficient, in particular, a two layer model representing a concrete
wall covered with a layer of insulating material.
I have attached a number of images to clarify the issue. The first,
'FiPy_TEMPW_insulation_only.png' shows the temperature distribution in a single
material (the insulation) with constant diffusivity, D_ins, when a constant
temperature of 1200 C is applied at one boundary for 3 hours. The FiPy result
agrees very well with the analytical solution
1 - erf(x/2(sqrt(D_ins t))*1200 +20,
taken from the examples.diffusion.mesh1D example (scaled and shifted
appropriately) and with a numerical solution calculated using TEMP/W (using the
same spatial and temporal discretisation, material properties and boundary
conditions). The three results agree well, showing that the FiPy model is
performing as expected.
The second figure, 'FiPy_TEMPW_insulated_concrete.png', presents the
temperature distribution through an insulated concrete wall (where for
simplicity the concrete is also modelled with constant diffusivity, D_con) for
the same surface temperature and period. There is now a considerable difference
between the FiPy and TEMP/W predictions. The FiPy model predicts a lower
temperature gradient in the insulation layer, which leads to higher
temperatures throughout the domain.
I am confident that the TEMP/W result is accurate, as it agrees perfectly with
a simple explicit finite difference solution coded in FORTRAN. I have tried to
identify any coding errors I have made in my FiPy script. I am aware that
diffusion terms are solved at the grid faces, so when the diffusion coefficient
is a function of a CellVariable, an appropriate interpolation scheme must be
used to obtain sensible face values. However, in my case the diffusion
coefficients, D_ins and D_conc, are created as FaceVariables and assigned
constant values. I have also examined the effects of space and time
discretisation, implicit and explicit DiffusionTerm, multiple sweeps per
timestep, but these have made no significant difference.
I would be very interested to hear anyone's opinion on what I might be doing
wrong here. Also, does anyone think it is possible for FiPy to deliver an
accurate result, and for the finite difference and finite volume solvers to be
wrong? Below this email I have written out a minimal working example, which
reproduces the 'FiPy' curve from figure 'FiPy_TEMPW_insulated_concrete.png'.
Many thanks,
Conor
--
Conor Fleming
Research student, Civil Engineering Group
Dept. of Engineering Science, University of Oxford
###################################################################
# FiPy script to solve 1D heat equation for two-layer material
#
import fipy as fp
nx = 45
dx = 0.009 # grid spacing, m
dt = 20. # timestep, s
mesh = fp.Grid1D(nx=nx, dx=dx) # define 1D grid
# define temperature variable, phi
phi = fp.CellVariable(name="Fipy", mesh=mesh, value=20.)
# insulation thermal properties
thick_ins = 0.027 # insulation thickness
k_ins = 0.212
rho_ins = 900.
cp_ins = 1000.
D_ins = k_ins / (rho_ins * cp_ins) # insulation diffusivity
# concrete thermal properties
k_con = 1.5
rho_con = 2300.
cp_con = 1100.
D_con = k_con / (rho_con * cp_con) # concrete diffusivity
valueLeft = 1200. # set temperature at edge of domain
phi.constrain(valueLeft, mesh.facesLeft) # apply boundary condition
# create diffusion coefficient as a FaceVariable
D = fp.FaceVariable(mesh=mesh, value=D_ins)
X = mesh.faceCenters.value[0]
D.setValue(D_con, where=X>thick_ins) # change diffusivity in concrete region
# unsteady heat equation
eqn = fp.TransientTerm() == fp.DiffusionTerm(coeff=D)
# specify viewer
viewer = fp.Viewer(vars=phi, datamin=0., datamax=1200.)
# solve equation
t = 0.
while t < 10800: # simulate for 3 hours
t += dt # advance time
eqn.solve(var=phi, dt=dt) # solve equation
viewer.plot() # plot result_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
