Hi, Jonathan.
I didn't know how to submit a pull request on MatForge, so here (attached
and pasted below) is a git diff of mesh1D.py. I added a small section
discussing the import of keeping the original governing equation in mind
when parameters vary with a specific example from heat transfer. I didn't
know where it would fit best in the example; I added it after the
non-uniform diffusivity section. If you have suggestions, I can edit.
Ray
diff --git a/examples/diffusion/mesh1D.py b/examples/diffusion/mesh1D.py
index 71f2da4..3c5209f 100755
--- a/examples/diffusion/mesh1D.py
+++ b/examples/diffusion/mesh1D.py
@@ -450,6 +450,45 @@ And finally, we can plot the result
------
+Note that for problems involving heat transfer and other similar
+conservation equations, it is important to ensure that we begin with
+the correct form of the equation. For example, for heat transfer with
+:math:`\phi` representing the temperature,
+
+.. math::
+ \frac{\partial \rho \hat{C}_p \phi}{\partial t} = \nabla \cdot [ k
\nabla \phi ].
+
+With constant and uniform density :math:`\rho`, heat capacity
:math:`\hat{C}_p`
+and thermal conductivity :math:`k`, this is often written like Eq.
+:eq:`eq:diffusion:mesh1D:constantD`, but replacing :math:`D` with
:math:`\alpha =
+\frac{k}{\rho \hat{C}_p}`. However, when these parameters vary either in
position
+or time, it is important to be careful with the form of the equation used.
For
+example, if :math:`k = 1` and
+
+.. math::
+ \rho \hat{C}_p = \begin{cases}
+ 1& \text{for \( 0 < x < L / 4 \),} \\
+ 10& \text{for \( L / 4 \le x < 3 L / 4 \),} \\
+ 1& \text{for \( 3 L / 4 \le x < L \),}
+ \end{cases},
+
+then we have
+
+.. math::
+ \alpha = \begin{cases}
+ 1& \text{for \( 0 < x < L / 4 \),} \\
+ 0.1& \text{for \( L / 4 \le x < 3 L / 4 \),} \\
+ 1& \text{for \( 3 L / 4 \le x < L \),}
+ \end{cases}.
+
+However, using a ``DiffusionTerm`` with the same coefficient as that in the
+section above is incorrect, as the steady state governing equation reduces
to
+:math:`0 = \nabla^2\phi`, which results in a linear profile in 1D, unlike
that
+for the case above with spatially varying diffusivity. Similar care must be
+taken if there is time dependence in the parameters in transient problems.
+
+------
+
Often, the diffusivity is not only non-uniform, but also depends on
the value of the variable, such that
On Mon, Aug 18, 2014 at 10:54 AM, Guyer, Jonathan E. Dr. <
[email protected]> wrote:
> Conor and Raymond -
>
> Thank you both for posting your findings and interpretation of the issue.
> I agree that this would be a useful issue to make clear in the
> documentation. We would welcome a patch or pull request from either of you
> to illustrate this situation in 'examples.diffusion.mesh1D'.
>
>
> On Aug 18, 2014, at 9:36 AM, Raymond Smith <[email protected]> wrote:
>
> > 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 ]
>
>
> _______________________________________________
> fipy mailing list
> [email protected]
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
diff --git a/examples/diffusion/mesh1D.py b/examples/diffusion/mesh1D.py
index 71f2da4..3c5209f 100755
--- a/examples/diffusion/mesh1D.py
+++ b/examples/diffusion/mesh1D.py
@@ -450,6 +450,45 @@ And finally, we can plot the result
------
+Note that for problems involving heat transfer and other similar
+conservation equations, it is important to ensure that we begin with
+the correct form of the equation. For example, for heat transfer with
+:math:`\phi` representing the temperature,
+
+.. math::
+ \frac{\partial \rho \hat{C}_p \phi}{\partial t} = \nabla \cdot [ k \nabla
\phi ].
+
+With constant and uniform density :math:`\rho`, heat capacity :math:`\hat{C}_p`
+and thermal conductivity :math:`k`, this is often written like Eq.
+:eq:`eq:diffusion:mesh1D:constantD`, but replacing :math:`D` with
:math:`\alpha =
+\frac{k}{\rho \hat{C}_p}`. However, when these parameters vary either in
position
+or time, it is important to be careful with the form of the equation used. For
+example, if :math:`k = 1` and
+
+.. math::
+ \rho \hat{C}_p = \begin{cases}
+ 1& \text{for \( 0 < x < L / 4 \),} \\
+ 10& \text{for \( L / 4 \le x < 3 L / 4 \),} \\
+ 1& \text{for \( 3 L / 4 \le x < L \),}
+ \end{cases},
+
+then we have
+
+.. math::
+ \alpha = \begin{cases}
+ 1& \text{for \( 0 < x < L / 4 \),} \\
+ 0.1& \text{for \( L / 4 \le x < 3 L / 4 \),} \\
+ 1& \text{for \( 3 L / 4 \le x < L \),}
+ \end{cases}.
+
+However, using a ``DiffusionTerm`` with the same coefficient as that in the
+section above is incorrect, as the steady state governing equation reduces to
+:math:`0 = \nabla^2\phi`, which results in a linear profile in 1D, unlike that
+for the case above with spatially varying diffusivity. Similar care must be
+taken if there is time dependence in the parameters in transient problems.
+
+------
+
Often, the diffusivity is not only non-uniform, but also depends on
the value of the variable, such that
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
