Level Set Numerical Surface Tension Calculation
 

Given a phase-function 0<phi(x,y)<1 that varies smoothly from 0 to 1 in a narrow region (approx. 6 grid cells from 0.05<phi<0.95), the gradients are:

dphi/dx
dphi/dy

which are approximated via central difference. This value is located at a grid cell center (i,j) due to central difference formula. The normal vector field is then approximated as

N = gradient/mag(gradient)

This is also located at cell center (i,j). The curvature is the divergence of normal vectors, i.e.

kappa = div(normal)

also located at cell center (i,j). Now...for the surface tension, whose x and y components are located at right edges (i+1/2,j) and top edges (i,j+1/2), respectively...

When I plot gradients, normals, and curvatures, they are all symmetric, i.e., around a circle, all values point perfectly radially outward (normal) to the level set phi=0.5 interface, but when I plot my surface tension, it appears to be asymmetric! I have even taken care to stagger these two vector fields (that are located at different spatial points, remember!) when I plot the field.

Does anybody have any advice? The surface tension (presented in the paper who's work I'm following, Olsson & Kreiss, "A conservative level set method for two phase flow") is a forward difference calculation, whereas normal and gradient calculations are central differences. I'm following their method EXACTLY, I've checked everything...and yet, asymmetric surface tension...help!