That did it. Thanks! ~Dan
On Thu, Oct 29, 2015 at 10:18 AM, Dan Brunner <[email protected]> wrote: > Thanks for the help Jon. I did not realize that the coeff in the > TransientTerm was assumed to only have a time dependence. I agree that > derivations of the heat equation are often not careful about non-linear > coefficients. I'll try your reformulation of it later this week. We have > some other tools in-house to benchmark it with. > > ~Dan > > On Wed, Oct 28, 2015 at 5:08 PM, Guyer, Jonathan E. Dr. < > [email protected]> wrote: > >> Dan - >> >> Interesting problem. What I actually find is that the residual decays or >> at least holds fairly steady until some point at which the residual rises >> rapidly, regardless of the number of sweeps. >> >> I was ultimately able to make some headway by looking at the form of your >> equations. By writing >> >> TransientTerm(coeff=rho*cp) >> >> you are saying that cp (and rho) have time dependence, but not spatial >> dependence. Clearly, in your case, the do have spatial dependence. In >> theory, sweeping should resolve some of this, but is apparently inadequate. >> >> Some derivations lump everything together as >> >> alpha = kappa/(rho*cp) >> >> and solve >> >> TransientTerm(coeff=1.) == DiffusionTerm(coeff=alpha) >> >> but I don't think this is right. >> >> The derivations of the heat equation I've seen that are careful about >> non-uniform coefficients (which is not most derivations of the heat >> equation that I've seen) end up with >> >> \rho c_p \frac{\partial T}{\partial t} = \nabla\cdot(\kappa\nabla T) >> >> This can be rearranged as >> >> \frac{\partial T}{\partial t} = \nabla\cdot(\kappa/(\rho c_p)\nabla T) >> - \nabla(1/(\rho >> c_p))\cdot\kappa\nabla T >> >> which FiPy can express as >> >> eq = (fp.TransientTerm(coeff=1.) >> == fp.DiffusionTerm(coeff=kappa/(rho*cp)) >> - (1/(rho*cp)).grad.dot(kappa*temperature.grad)) >> >> >> I don't know if it's right, but it solves. >> >> You also should be able to increase the time step considerably. >> >> - Jon >> >> >> >> On Oct 27, 2015, at 9:11 AM, Dan Brunner <[email protected]> wrote: >> >> > I'm working on a nonlinear diffusion problem (specifically finding how >> long a spatially-dependent heat flux can be held on a surface before it >> reaches the melting point, 'simple' example at the end) and I am seeking >> assistance in setting up the sweeping. I've looked through all of the FiPy >> examples to learn what I can. If I set it up to sweep until a certain >> residual threshold is reached, each iteration on the time loop takes longer >> and longer because the initial residual keeps increasing and the residual >> decrease between sweeps keeps decreasing. This usually goes on until the >> residual bounces up and down in a sweep loop. If I try some logic to >> decrease the sweep dt, the sweep residuals blow up. >> > >> > Does anyone have guidance on what to do here and/or can refer me to a >> resource where I may learn how to do this better? >> > >> > Thanks >> > ~Dan >> > >> > from fipy import * >> > import numpy as np >> > >> > q0=5e8 >> > qLambda=1e-2 >> > Lx=qLambda*4.0 >> > dx=qLambda/20.0 >> > Ly=0.02 >> > dy=Ly/2.0 >> > T0=300.0 >> > peakTemp=3695.0 >> > resMin=1e-2 >> > >> > mesh=Gmsh2D(''' >> > Lx=%(Lx)g; >> > Ly=%(Ly)g; >> > dx=%(dx)g; >> > dy=%(dy)g; >> > Point(1) = { 0, 0, 0, dx}; >> > Point(2) = {Lx, 0, 0, dx}; >> > Point(3) = {Lx, Ly, 0, dy}; >> > Point(4) = { 0, Ly, 0, dy}; >> > Line(1) = {1, 2}; >> > Line(2) = {2, 3}; >> > Line(3) = {3, 4}; >> > Line(4) = {4, 1}; >> > Line Loop(5) = {1, 2, 3, 4}; >> > Plane Surface(6) = {5}; >> > ''' % locals()) >> > X, Y = mesh.faceCenters >> > >> > time = Variable() >> > >> > temperature = CellVariable(name = "temperature", >> > mesh = mesh, >> > value = T0, >> > hasOld=1) >> > # thermal properties of tungsten >> > rho = 19300.0 >> > cp = (5.33e+01*temperature**0.0+ >> > 4.50e-01*temperature**1.0- >> > 8.30e-04*temperature**2.0+ >> > 7.07e-07*temperature**3.0- >> > 2.73e-10*temperature**4.0+ >> > 3.91e-14*temperature**5.0) >> > kappa = (2.32e+02*temperature**0.0- >> > 2.67e-01*temperature**1.0+ >> > 2.42e-04*temperature**2.0- >> > 1.13e-07*temperature**3.0+ >> > 2.50e-11*temperature**4.0- >> > 1.99e-15*temperature**5.0) >> > >> > eq = TransientTerm(coeff=rho*cp) == DiffusionTerm(coeff=kappa) >> > >> > dTime = 0.9*dx**2.0/(2.0*5e-5) >> > >> > if __name__ == '__main__': >> > viewer = Viewer(vars=temperature, datamin=0., datamax=4000.) >> > viewer.plot() >> > >> > while temperature().max() < peakTemp: >> > time.setValue(time()+dTime) >> > temperature.updateOld() >> > >> > facesDecayBottom = (mesh.facesBottom & (X>=(qLambda))) >> > valueDecayBottom = >> -q0*np.exp(-(X-qLambda)/qLambda)/kappa.getFaceValue() >> > >> > temperature.faceGrad.constrain(0.0, mesh.facesBottom) >> > temperature.faceGrad.constrain(valueDecayBottom, facesDecayBottom) >> > temperature.faceGrad.constrain(0.0, mesh.facesTop) >> > >> > res=1e9 >> > while res>resMin: >> > res = eq.sweep(var=temperature,dt=dTime) >> > print res >> > >> > print time, temperature().max() >> > if __name__ == '__main__': >> > viewer.plot() >> > _______________________________________________ >> > 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 ]
