Suppose we want to solve the steady low-Mach-number Navier-Stokes equations coupled with a passive scalar $\xi$, which read:
\begin{align} \nabla \cdot \rho \mathbf{v} & = 0,\\ \rho \mathbf{v} \cdot \nabla \mathbf{v} + \nabla p & = \nabla \cdot \tau + \mathbf{F}(\rho) \\ \rho \mathbf{v} \cdot \nabla\xi & = \nabla \cdot \alpha \nabla \xi, \end{align} where $\rho$ is density defined as a function of $\xi$ through some state equation $\rho = \rho(\xi)$, $\alpha$ is a diffusivity coefficient, $\mathbf{F}$ is some function accounting for buoyancy, $p$ is pressure and $\tau$ is the usual stress tensor. A vaporizing boundary condition must be imposed at $x = 0$, $ r \leq 1$ in an axisymmetric domain:
\begin{align} v_r = 0, \quad \xi = 1, \quad \rho v_x = - \alpha \partial_x \xi, \end{align}
where the last one accounts for the rate of vaporization of the fuel at that surface. Note that we have a Dirichlet BC for $\xi$ and a boundary condition that relates the normal part of the velocity at the inlet, $v_x$, and the normal gradient of $\xi$.
I am interested in solving this problem using a finite element solver (for instance, FreeFem++) for which we need to introduce these equations in weak form. The problem comes when I try to figure out the corresponding terms for the boundary condition written above. If $q$, $\mathbf{w}$ and $\varphi$ are test functions for each equation, I believe the problem is the following:
\begin{align} \int_\Omega q \nabla \cdot \rho \mathbf{v} \, \mathrm{d}\Omega & = 0, \\ \int_\Omega \left[ \rho (\mathbf{v} \cdot \nabla \mathbf{v})\cdot \mathbf{w} - p \nabla \cdot \mathbf{w} + \tau : \nabla \mathbf{w} - \mathbf{F} \cdot \mathbf{w} \right] \, \mathrm{d} \Omega & = 0, \\ \int_\Omega \left[ \rho (\mathbf{v} \cdot \nabla \xi) \, \varphi + \alpha \nabla \xi \cdot \nabla \varphi \right] \, \mathrm{d} \Omega & = 0, \end{align} where one must take into account the corresponding boundary terms (we assume stress-free boundary conditions and no-slip).
Note that $\xi = 1$ and $v_r = 0$ are Dirichlet boundary conditions, for which no boundary term arises at the inlet. But, I wonder what kind of boundary condition is $v_x = -\partial_x \xi$ and how to impose it.
I have tried writing continuity equation in weak form and then substitute the boundary term $\int_{\partial \Omega_{in}} \rho \, \mathbf{v} \cdot \mathbf{n} \, q \, \mathrm{d}\Gamma$ by the boundary condition, but this has been proved unsuccessful. I have to rely on some iterative procedure to satisfy the prescribed balance, but this is not very elegant...
Any hints or thoughts are greatly appreciated.