Guyer, Jonathan E. Dr. <jonathan.guyer@...> writes:
>
>
> On Jul 28, 2014, at 3:18 AM, Burak Atakan <herzige@...> wrote:
>
> > After playing around I also found that putting the T.faceGrad constraint
> > into the sweep also resoved the problem.
>
> That sounds like a bug. The constraint should update automatically. Can
you provide the code that you ended
> up using?
>
> _______________________________________________
> fipy mailing list
> fipy@...
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
Here we are:
#---------------------
from fipy import *
#
part=0 # heating the whole bottom or only part of it (select: 0/1)
L=0.025
H=0.005
hb=0.002 #Size of the Heating element (relevant when part>0)
dx=0.000125/2.
dy = dx
nx = int(L/dx)
ny = int(H/dy)
L = dx * nx
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
#
T_bot = 500.
Tinf=300.
alph=60. # convection coefficient
D=1.6e1 # thermal conductivity
T0=zeros((nx,ny))
for x in range(nx): #initial guess T-Profile
T0[x,0:ny]=linspace(T_bot-1,Tinf+1,ny)
#
X, Y = mesh.faceCenters
T = CellVariable(name = "solution variable",
mesh = mesh, hasOld=1)
T.setValue(reshape(T0,nx*ny))
#
eq = DiffusionTerm(coeff=D)
#
viewer=Viewer(T)
#
res=1e10
co=0
if part>0:
facesBottomLeft = ( mesh.facesBottom & (X <=(L+hb)/2.)&(X >=(L-hb)/2.))
else:
facesBottomLeft = ( mesh.facesBottom)
T.constrain(T_bot, facesBottomLeft)
while res>1e-4:
#faktor=.01
T.faceGrad.constrain((alph/D*(Tinf-T.faceValue)*
mesh.faceNormals),where=mesh.facesTop)
#T.faceGrad.constrain((alph/D*(Tinf-T.faceValue.value)*
mesh.faceNormals),where=mesh.facesTop)
# the latter is faster for more complicated problems...
res=eq.sweep(var=T)
print co,res
co+=1
#
ma,mi=max(T.faceValue),min(T.faceValue)
print "Max,Min",ma,mi,ma-mi
viewer.plot()
Bi=alph*H/D
Ts=(T_bot+Tinf*Bi)/(1+Bi)
print Ts,"=expected surface value; ",T_bot-Ts,Ts-Tinf,Bi
woX=where(X==min(X)) #links
woY=where(Y==max(Y)) #oben
woY0=where(Y==min(Y))
figure(17)
plot(Y[woX],T.faceValue.value[woX])
plot(X[woY],T.faceValue.value[woY])
plot(X[woY0],T.faceValue.value[woY0])
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
