Surface curvature for cylindrical jet for determining surface tension 
Hi,
I have this problem of the stability of a cylindrical jet. The things I do not understand is the expression for the surface curvature for the cylindrical jet, which I need in order to find the difference in pressure across the interface caused by surface tension. I read some text, the curvature is calculated from the divergent of the surface normal vector. Δp = γ ∇ · n. Since it is a cylindrical coordinate problem, so I think I should be using the stardard divergence formula for cylindrical coordinate:
$$
\nabla \cdot \boldsymbol{u}=\frac{1}{r} \frac{\partial\left(r u_{r}\right)}{\partial r}+\frac{1}{r} \frac{\partial u_{\phi}}{\partial \phi}+\frac{\partial u_{z}}{\partial z}
$$
But this is one of the things I would like to clarify is it the correct way, since when it differentiate term involving the coordinate r, r is itself a function of detla and z.
To calculate the unit surface normal, I use the concept of the surface can be represented by a position vector in the r and z components. So my position vector when represented in cylindrical coordinate only has two components r and z, but without $ \theta$ component.
$\vec x =r(z,\theta)\hat e_r +z \hat e_z $
and given in the problem, the surface have an r coordinate:
$r=a+\epsilon(\theta,z)$ where a is "a" constant.
To find the two tangent vectors, I differentiate this position vector with respect to $\theta$ and z. I followed the product ruled that when I differentiate a product of the component and the unit vector, because the unit vector in r and $\theta$ is function of $\theta$.
$d(r\hat e_r)/d\theta=\hat e_r dr/d\theta + rd(\hat e_r)/d\theta=\hat e_r dr/d\theta + r\hat e_\theta$
In general: $d(\hat e_r)/d\theta=\hat e_\theta$  and  $d(\hat e_\theta)/d\theta=-\hat e_r$
Finally I cross product these two tangent vector to find the surface vector, and divide that by its length to make it a unit normal.
$d\vec x/d\theta \times d\vec x/dz=(\hat e_r\partial\epsilon/\partial\theta + r \hat e_\theta )\times(\hat e_r\partial\epsilon/\partial z+  \hat e_z)=\hat e_r r - \hat e_\theta \partial\epsilon/\partial\theta - \hat e_z r \partial\epsilon/\partial z$
I made the assumption that the square of the derivatives in the epsilon with respect to z and $\theta$ is much smaller than one:
$(d\epsilon/d\theta)^2<< 1$  and   $(d\epsilon/dz)^2<< 1$
This make the length of $d\vec x/d\theta \times d\vec x/dz$ just equal to r. Finally I substitute the three components of the unit surface normal into the divergent formula for cylindrical coordinate to find the expression for different pressure caused by surface tension.
$\vec n =\hat e_r - \hat e_\theta (1/r)\partial \epsilon/\partial\theta - \hat e_z \partial \epsilon/\partial z$
The problem is I always ended up without the $-\epsilon/a^2$ that is indicated in the formula shown in the first image. This is what I always get:
$\nabla \bullet \vec n =(1/r)(1+(\frac{1}{r} \frac{\partial \epsilon}{\partial \theta})^2)-(1/r^2)\partial^2 \epsilon/d\theta^2 -\partial^2 \epsilon/dz^2$
The terms of square of the derivative of epsilon with respect to $\theta$ is supposed to be neglected, since it is second order in epsilon.
So I am wondering that maybe something I am doing wrong, or not understanding it right. Apperciate any kind of ideas or advice.
Thank in advance
 A: Ok, this comes one hour after I posted the question. The method used to calculate the curvature is indeed  correct. The difference in the expression is due a linearization process that linearlized the expression to first order in epsilon.
$\nabla \bullet \vec n =(1/r)(1+(\frac{1}{r} \frac{\partial \epsilon}{\partial \theta})^2)-(1/r^2)\partial^2 \epsilon/d\theta^2 -\partial^2 \epsilon/dz^2$
should be simplified to 
$\nabla \bullet \vec n =(1/r)-(1/r^2)\partial^2 \epsilon/d\theta^2 -\partial^2 \epsilon/dz^2$
where $r=a+\epsilon$
The expression have terms involving $1/r$ and $1/r^2$
Writing $1/r$ and $1/r^2$ as r^(-1)=(a+$\epsilon$)^(-1) and r^(-2)=(a+$\epsilon$)^(-2) respectively, and factoring out the $a$, from each expression
r^(-1)=(a+$\epsilon$)^(-1)=a^(-1)(1+$\epsilon$/a)^(-1)
r^(-2)=(a+$\epsilon$)^(-1)=a^(-2)(1+$\epsilon$/a)^(-2)
Using the fact that $a$/$\epsilon$ is a small quantity, expand each  term in the bracket to the first order of $\epsilon$
r^(-1)=(a+$\epsilon$)^(-1)=a^(-1)(1+$\epsilon$/a)^(-1)=a^(-1)(1-$\epsilon$/a)
r^(-2)=(a+$\epsilon$)^(-1)=a^(-2)(1+$\epsilon$/a)^(-2)=a^(-2)(1-2$\epsilon$/a)
Substitute them back into the expression of
$\nabla \bullet \vec n =(1/r)-(1/r^2)\partial^2 \epsilon/d\theta^2 -\partial^2 \epsilon/dz^2$
And neglecting second order terms in $\epsilon$ will give the desired result.
