Neumann boundary issues


Pneumatik
06-03-2004, 02:13 PM
I'm working on an FMG model for electrical current passing through a thin sheet of metal. The model takes advantage of radial symmetry to only deal with a plane that starts at r=0 and extends out to an arbitrary distance. The current passes through the thickness of the metal sheet. The contact areas on either face are the same size and shape and exactly opposite each other, so I'm only modelling the top half of the metal sheet.

All of these conditions make for an interesting model. I have Dirlecht boundary conditions at the bottom of the model u=0. On either side of the model I have Neumann boundary conditions of du/dr = 0. On the top of the model, roughly half of the boundary has du/dz = 0 and the other half is du/dz = C, where C is a nonzero constant.

As far as I can tell, all the components of my FMG model work correctly. I relax using a symmetric Gauss-Seidel sweep that gives reasonable results. My prolongation is bilinear interpolation and my restriction is full-weighting, both of which are taken right from Numerical Recipies. If I have matrix A, R*(P*A) == A.

I'm having a couple of problems with my model. The one that bothers me the most is that the finer my grid gets, the larger my values at the points along the du/dz = C boundary get. I've discritized that boundary as:

(u(-1) - u(1))/2[delta]z = C --> u(-1) = C*2*[delta]z + u(1)

This particular discretization means that as my mesh gets finer at each level, the C*2*[delta]z term changes (when it's put into the fundamental FDE equation for my problem [similar to 19.5.6 in Numerical Methods], it ends up being divided by ([delta]x)^2, so it ends up growing in magnitude as the mesh shrinks). I've seen this particular discretization used in other places, so I'm assuming that part is correct.

Note: when I say that the values at the points adjacent to the injected electrical current increase, I mean the values that those points asymtotically approach as I constantly relax the mesh increase. If I prolong from one mesh to the next finer one, the values stay very much the same.

I end up having a couple of problems with my FMG model because of the discretization for du/dz = C. The increase in values adjacent to the electric current input means that I can't keep refining my mesh until the values start to converge. Because the values increase each time I refine the mesh, my residual is larger on the /upward/ stroke of each V in my model than on the downward stroke. When I solve exactly (by repeatedly relaxing the coasest mesh) at the bottom of the V, my residuals become very small. When I then prolong to the next finer mesh without increasing my mesh values, the residuals go up. My prolongation operator does not involve any boundary conditions, but the restriction operator does.

Because my residuals don't decrease on the upstroke of my V, I don't know if I should include more relaxation steps at each level. There's a noticable difference in results between if I use a fixed number of relaxation sweeps (say, 4) and if I relax until my relaxation starts to convernge. What makes relaxing until it starts to converge difficult is that as the mesh gets finer, it both takes more sweeps before it starts to converge and each sweep takes longer.

Once I have some idea how many relaxation sweeps I should use, I'll be able to figure out how many partial V's I should use at each level (whether I should have one v, or two - vv, or three - vvv, etc.).

I hope all of that made sense. I'm open to ideas, because I think I'm really stuck here. I'm considering normalizing my results at each level, since that should work well enough for my needs. I've spent a few months working on this, though, so I'd really like it to give the correct answer.

John Spey

hernlund
06-19-2004, 07:59 PM
Well, you don't want all that time to go to waste!

It sounds as if you are using a collocated grid with mixed boundary conditions. I.e. the way you've discretized the problem, the dependent values lie along the boundaries. This works well for Dirichlet problems, but in general it may cause some difficulties when mixed boundary conditions are used. In multi-grid applications, this requires very careful treatments of the prolongations at the boundaries (which is very likely the source of your poor v-cycle performance).

I solve these types of problems quite often, and I've found that using a cell-centered discretization works better with mixed boundary conditions. It is also more forgiving in general for multi-grid problems, see Wesseling's book (free online in pdf) for a discussion of this. For most problems, it also yields a nice "flux-balanced" numerical approximation as well (via the cell measure of divergence). The NR book does not cover this type of approach.

For this discretization, your scalar dependent variables (electrical potential, perhaps) assume values at, say, x=1/2*dx, 3/2*dx, 5/2*dx, etc., and same in other dimensions. I.e. the grid is just offset by one-half grid spacing in all dimensions.

So, say youve done this, then a Dirichlet boundary condition is simply an average (using your notation):
1/2*(u(-1/2)+u(1/2))=boundary value
or
u(-1/2)=2*(boundary value)-u(1/2)

For a specified normal derivative or "flux":
(u(1/2)-u(-1/2))/dz=boundary flux value
or
u(-1/2)=u(1/2)-dz*(boundary flux value)

Both of these are second-order, which preserves the second-order-ness of the overall discretization scheme. This type of thing is especially good for axisymmetric problems, where r=0 is an axis of symmetry.

NOTE: You'll probably have to change the indices for your loops as well as the prolongation/restriction (see Wesseling for good choices).

Good Luck!