Hi,
thanks for your reply.
You can find attached a sample problem.
The equation is a heat equation, variable is temperature,other and transient
and diffusion term, there is a spatially varying source term (q) computed by
analytical expression which is and hyperbola.
We should reach about 100°C in 30-35 sec and the shape of the solution should
be similiar to and hyperbola.
As it is now the file works fine (mesh declared with: mesh =
CylindricalGrid1D(dr=dr, nr=(len(DR))) + (r_int,) ) but if you define the mesh
with ( mesh = CylindricalGrid1D(dx=DR) + (r_int,) ) the result you obtain are
wrong.
The 2 declaration of CylindricalGrid1D should be the same, the face and cell
centers are the same, but the result are different.
You can observe the same problem also when solving the same equation but
without transient term.
Thanks
Eduard
> Date: Tue, 11 May 2010 10:45:24 -0400
> From: [email protected]
> To: [email protected]
> Subject: Re: Problem in solving Poisson equation with 1D Cylindrical mesh
>
>
> Can you put together the simplest script possible that demonstrates
> the problem and I'll try and debug it? Thanks.
>
> On Thu, May 6, 2010 at 2:01 PM, Eduard Manley <[email protected]> wrote:
> > Thanks for your reply.
> >
> > I probably found the reason of the problem.
> >
> > As said before I'm trying to solve an eq of this type:
> >
> > A(d phi/d t) = div (D grad phi) + q
> >
> > where A and D are costant coefficient and q is a spatially varying heat
> > source.
> >
> > The problem is in how I create the cylindrical 1D mesh (origin of the mesh
> > is not in 0.).
> >
> > It doesn't matter if the discretization is logarithmic or uniform but how I
> > declare it:
> >
> > (using a uniform spacing:)
> >
> > ** mesh = CylindricalGrid1D(dx=dr, nx=(len(DR))) + (r_int,) **
> > SHOULD BE EQUAL TO:
> >
> > ** mesh = CylindricalGrid1D(dx=DR) + (r_int,) **
> >
> > [dr = 5e-04, nr=58, r_int=0.00125, DR is a list which contains the various
> > dx(58 elements of value dr for uniform grid)]
> >
> > BUT It is NOT.
> >
> > The mesh (cell centers, facecenters) is ok but the results are NOT.
> > The results are right only if I create the mesh using dx=dr and nx=. .
> > And this is why before I thought the problem was the logaritmic
> > discretization (must use DR=[...]).
> >
> > As said before with Cylindrical 2D mesh results are instead correct.
> > Is this a bug?
> >
> >> Date: Thu, 6 May 2010 11:11:10 -0400
> >> From: [email protected]
> >> To: [email protected]
> >> Subject: Re: Problem in solving Poisson equation with 1D Cylindrical mesh
> >>
> >>
> >> Edward, This may be to do with having a very small volume (or area or
> >> line or point) for the inner most element of the domain. It should be
> >> the same whether one is using a 1D or 2D mesh. Since you are getting
> >> differences in the 1D and 2D case, it should be relatively easy to
> >> debug and figure out what's going on It could also be that the
> >> boundary condition on the inner boundary as zero area and this is
> >> causing issues. Try shifting the grid by a small value away from the
> >> zero point and see if things are improved. I have always had this
> >> issue with cylindrical grids and have never really had a satisfactory
> >> solution (other than shifting away from the zero point). If you
> >> discover a better way to handle this, let me know. Cheers.
> >>
> >>
> >> If you can't debug it, then send me the most minimalist scripts that
> >> show the issue and I'll give it a shot.
> >>
> >> On Wed, May 5, 2010 at 9:53 PM, Eduard Manley <celez1...@hotmailcom>
> >> wrote:
> >> > Problem partially solved:
> >> >
> >> > I'm using a logarithmic discretization (first dr= 5e-04 and next
> >> > dr increasing as 1.05)(with internal radius= 0.00125, external
> >> > radius=0.03)
> >> > and,
> >> > for some unknown reason this create problems and wrong result with
> >> > cylindrical 1D mesh. I tried using a uniform discretization (dr=5e-04
> >> > nx=58)
> >> > and now the result is correct.
> >> >
> >> > However I need to use the logarithmic discr. so after some hours of
> >> > sleep
> >> > I'll think about the reason.....
> >> >
> >> > ________________________________
> >> > Hotmail: Trusted email with powerful SPAM protection. Sign up now.
> >>
> >>
> >>
> >> --
> >> Daniel Wheeler
> >>
> >>
> >
> > ________________________________
> > Your E-mail and More On-the-Go. Get Windows Live Hotmail Free. Sign up now.
>
>
>
> --
> Daniel Wheeler
>
>
_________________________________________________________________
Hotmail: Trusted email with Microsoft’s powerful SPAM protection.
https://signup.live.com/signup.aspx?id=60969from fipy import *
dr= 5e-04
DR =[dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr,
dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr,
dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr, dr]
r_int = 0.00125
R_est = len(DR)*dr + r_int
#mesh = CylindricalGrid1D(dx=DR) + (r_int,)
# SHOULD BE EQUAL TO:
mesh = CylindricalGrid1D(dr=dr, nr=(len(DR))) + (r_int,)
x = mesh.getCellCenters()[0]
sigma = (1./3.) #Electrical conductivity
W = 20. #POWER
sp = 0.02 # depth (z direction in cylindrical coord.)
Z = (((1/sigma)*log((R_est/r_int)))/(2*pi*sp)) #electrical impedance
I = sqrt((W/Z)) # electric current
J = I/(2*pi*sp*x) #Current density
q = CellVariable(name = "Resistive Heat", mesh = mesh, value = 0.)
q = (1/(sigma)) * (J**2) #heat source by joule heating
teta = CellVariable(name = "Temperature", mesh = mesh, value = 20.)
totaltime = CellVariable(name = "Time", mesh = mesh, value = 0.)
# BCs
valueLeft = 0.
valueRight = (20.)
facesLeft = mesh.getFacesLeft()
facesRight = mesh.getFacesRight()
BCs = (FixedValue(faces=facesRight, value=valueRight),
FixedFlux(faces=facesLeft, value=valueLeft))
# viewer
if __name__ == '__main__':
viewer = Viewer(vars=teta, datamin=0., datamax=100.)
viewer.plot()
# def time steps and equation
timeStepDuration = 1
steps =40
time = 0.
for step in range(steps):
time += timeStepDuration
totaltime.setValue(time)
OmegaH = 3500000.
GammaH = 0.5
eqH = (TransientTerm(coeff = OmegaH)) == (DiffusionTerm(coeff = GammaH)) +
(q)
#eqH = (DiffusionTerm(coeff = GammaH))+ q == 0
eqH.solve(var=teta, boundaryConditions=BCs, dt=timeStepDuration)
viewer.plot()
