Hi Conor,
Just to add to your observations, I'm guessing FiPy is "correct" here in
both situations, and as you noticed, the original input with a thermal
diffusivity is simply not the correct representation of your physical
situation. The actual conservation equation for thermal energy is
(neglecting convection and source terms)
\frac{\partial rho * c_p * T}{\partial t} = \nabla\cdot(k\nabla T)
So if you have spatially varying rho and c_p, you cannot sweep them into
the coefficient for the flux term, else they will appear *inside *FiPy's
divergence operator. They are often combined with k to form alpha for the
convenience of writing it like the mass conservation (diffusion) equation
(when chemical diffusivity is constant), but it's always important to start
from the correct conserved quantities. I would guess that if you run it
with any rho, c_p that are the same in the concrete and insulator, you will
get the correct result with the "alpha form" because, being both constant
and uniform, they can be swept into and out of either type of partial
derivative.
Cheers,
Ray
On Mon, Aug 18, 2014 at 8:51 AM, Conor Fleming <[email protected]>
wrote:
> 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]>
> 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]
>> 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 ]
>
>
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
