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 <conor.flem...@eng.ox.ac.uk> 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 > fipy@nist.gov > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > >
_______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
