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 ]
