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 ]
