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 ]
