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 ]
