Neumann boundary condition on axis in cylindrical coordinates

335
views
3
10 months ago by
I'm trying to solve the eigenvalue problem from a linear stability analysis about a steady Navier-Stokes solution:

$s\hat{\mathbf{u}} + (\mathbf{u}_b . \nabla) \hat{\mathbf{u}} + (\hat{\mathbf{u}}.\nabla)\mathbf{u}_b + \nabla\hat{p} - \nu \nabla^2 \hat{\mathbf{u}} = 0$

In particular I'm looking at solutions in cylindrical coordinates where $\hat{\mathbf{u}}(r,z,\theta) = \tilde{\mathbf{u}}(r,z)e^{mj\theta}$ with $m=\pm1$.  An asymptotic analysis shows that the boundary conditions for the $m\pm-1$ modes as $r\rightarrow 0^+$  are $\frac{\partial \tilde{u}_\theta}{\partial r}=\frac{\partial \tilde{u}_r}{\partial r}=0$ .

I'm using axisymmetric finite elements with $\int_\Omega d\Omega \rightarrow \int_{0}^L \int_{0}^R \int_{0}^{2\pi} r d\theta dr dz$ and a test function of the form $\mathbf{v}(r,z,\theta) = \tilde{\mathbf{v}}(r,z)e^{-mj\theta}$. Following recommendation from 'Global mode interaction and pattern selection in the wake of a disk: a weakly nonlinear expansion, Meliga et al. JFM 633 [159-189]' I multiply the governing equations by $r$ to ensure that there are no $\frac{1}{r}$ terms left in the variational form.

Unfortunately I can't work out how I should set the radial derivative boundary conditions on the axis. All the terms involving $\frac{\partial \tilde{u}_\theta}{\partial r}$ in boundary integrals on the axis look like:
$\int_{\Gamma_{\text{axis}}} r^2 \frac{\partial \tilde{u}_r}{\partial r} 2\pi dz$

Normally I would set Neumann boundary conditions by just setting this boundary integral to zero. Unfortunately, because on the axis $r=0$,  this term is automatically zero regardless of what the radial derivative is.

Is there something I'm missing? Are there any other common methods of enforcing a Neumann boundary condition? Is there a better way to set up the axisymmetric formulation?

Thanks for any help.
Community: FEniCS Project

Your observation is correct. Since  $r=0$r=0 on the axis, the boundary term goes away naturally. You shouldn't need to set a boundary condition there and it's unphysical to try to set a boundary condition on the axis since you are really solving a problem in a cylinder.