Hi Kris,
I have identified the cause of this difference - I had not specified the
governing equation correctly in FiPy. Originally, I had prescribed the heat
equation as follows:
eqn = TransientTerm() == DiffusionTerm(coeff=alpha)
where the thermal diffusivity alpha is
alpha = k / (rho * c_p),
and FiPY gives an erroneous answer. Additionally, solving the steady equation
'DiffusionTerm(coeff=alpha)==0' also yields an incorrect result where (rho*c_p)
has a value other than unity. I have realised that the equation should be
specified as
eqn = TransientTerm(coeff=C) == DiffusionTerm(coeff=k)
where the volumetric heat capacity is C = rho * c_p and k is the thermal
conductivity. For a steady case, this equation reduces to
'DiffusionTerm(coeff=k)==0' and gives a correct result.
I have attached a figure with the updated comparison of FiPy and TEMP/W for an
insulated concrete slab, showing perfect agreement.
I hope this result is useful to other FiPy users. It may be helpful to add a
note to the documentation page 'examples.diffusion.mesh1D', explaining that for
some applications, e.g. the heat equation, it is appropriate to separate the
(thermal) diffusivity into two portions, which act on the TransientTerm and
DiffusionTerm respectively.
Many thanks,
Conor
________________________________
From: [email protected] [[email protected]] on behalf of Conor Fleming
[[email protected]]
Sent: 08 August 2014 17:24
To: [email protected]
Subject: RE: Unexpected result, possibly wrong, result solving 1D unsteady heat
equation with spatially-varying diffusion coefficient
Hi Kris,
Thank you for the prompt response. You are right - altering the insulation
conductivity in the FiPy model to k_ins=0.1 improves the agreement. I have just
checked TEMP/W again, and can confirm that k_ins=0.212 in that model. Specific
heat capacity and density match as well. I will read up on FE and FD
implementations generally to see if I can spot any issues on that side.
Many thanks,
Conor
________________________________
From: [email protected] [[email protected]] on behalf of Kris Kuhlman
[[email protected]]
Sent: 08 August 2014 17:04
To: [email protected]
Subject: Re: Unexpected result, possibly wrong, result solving 1D unsteady heat
equation with spatially-varying diffusion coefficient
Conor,
if you reduce the thermal conductivity in the insulation to about 0.1, the fipy
solution looks about like the other model (the knee in T is about at 400
degrees C). Is there an issue with how your compute or specify this in fipy or
the other model?
Kris
On Fri, Aug 8, 2014 at 9:35 AM, Conor Fleming
<[email protected]<mailto:[email protected]>> wrote:
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]<mailto:[email protected]>
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
