Hi Ray, Thanks for your additional explanation of this situation. I think your suggested edit is ideal.
Kind regards, Conor On 18 Aug 2014, at 19:00, Raymond Smith <[email protected]<mailto:[email protected]>> wrote: 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]<mailto:[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]<mailto:[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]<mailto:[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]<mailto:[email protected]> > [[email protected]<mailto:[email protected]>] on behalf of Conor > Fleming [[email protected]<mailto:[email protected]>] > Sent: 08 August 2014 17:24 > To: [email protected]<mailto:[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]<mailto:[email protected]> > [[email protected]<mailto:[email protected]>] on behalf of Kris > Kuhlman [[email protected]<mailto:[email protected]>] > Sent: 08 August 2014 17:04 > To: [email protected]<mailto:[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]<mailto:[email protected]> > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > > _______________________________________________ > 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]<mailto:[email protected]> http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] <mesh1D.diff>_______________________________________________ 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 ]
