* Replacing the T.faceGrad constraint with a flux divergence source helps: eq = DiffusionTerm(coeff=D) + (alph * (Tinf-T.faceValue) * mesh.faceNormals * mesh.facesTop).divergence
Max,Min 500.0 450.943343505 49.0566564949 450.0 expected surface value 50.0 150.0 0.333333333333 * Increasing the mesh resolution to dx = dy = 0.000025 helps further: Max,Min 500.0 450.093747663 49.9062523366 450.0 expected surface value 50.0 150.0 0.333333333333 Surprisingly, replacing T.faceValue.value with T.faceValue doesn't make any difference, but you should do it anyway. Wheeler may be able to speculate why the T.faceGrad constraint doesn't do a very good job in this case. On Jul 25, 2014, at 8:03 AM, Burak Atakan <[email protected]> wrote: > Dear all, > again I have a problem with a simple convection boundary condition in a heat > transfer problem. The problem will get more complicated later, but for > verification I set up a 2 D Mesh set to a fixed T at the bottom and having > convection boundary conditions at the top. The equation is just the steady > state conduction equation. > 0= k (d^2T/dx^2+d^2T/dy^2) > the convection boundary condition should lead to: > dT/dy=h/k(Tinf-T) > > (Tinf is a constant fluid Temperature, h the convection coeff. and k the > thermal conductivity) > > The calculated surface temperature in my script depends on the value I > assign initially to T (T_bot -53 in the script). The simple analytical > result is calculated at the end, but I never get the correct result. > > What am I doing wrong? (I first tried solve() before trying sweep()...) > > Thanks for your help! > Burak > -------------- > from fipy import * > # > L=0.02 > H=0.005 > dx=0.00025 > 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=100. #convection coefficient > D=1.5 #thermal conductivity > > > # > X, Y = mesh.faceCenters > #facesTopLeft = ( mesh.facesTop )#& (X < L *.1)) > facesBottomLeft = ( mesh.facesBottom)# & (X <=2*dx)) > T = CellVariable(name = "solution variable", > mesh = mesh,value=T_bot-53, hasOld=1) > # > #One heated position at the bottom (later, at the moment all is heated) > T.constrain(T_bot, facesBottomLeft) > #cooling on the top due to convection > > T.faceGrad.constrain((alph/D*(Tinf-T.faceValue.value)* > mesh.faceNormals),where=mesh.facesTop) > eq = DiffusionTerm(coeff=D) > # > viewer=Viewer(T) > # > res=1e10 > co=0 > while co<1e1: > res=eq.sweep(var=T,dt=1e5) > co+=1 > print res,co#eq.solve(var=T) > # > 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 > > _______________________________________________ > 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 ]
