Part I (An Example): It's easy to construct examples of such pairs $(g,X)$ on an unbounded polygonal domain whose boundary is connected:
Let $\Omega$ be such a domain, and let $w = u+iv:\Omega\to\mathbb{C}$ be a Schwarz-Christoffel transformation that maps $\Omega$ one-to-one onto the upper half plane and sends the point at infinity to the point at infinity. Then $v$ vanishes on the boundary of $\Omega$ and is positive on the interior of $\Omega$. Let $g = \mathrm{d}u^2 + \mathrm{d}v^2$, which is a positive definite form on the interior of $\Omega$ (in fact, conformal to the standard metric on $\Omega$).
Let $X=-\partial/\partial v$. This solves the problem, since $X$ is perpendicular to the boundary of $\Omega$ and it clearly has divergence zero and length $1$ with respect to $g$.
Of course, $X$ is not continuous at the corner points. However, the metric $g$ is smooth and conformal (in fact, flat) in the interior of $\Omega$ and it extends smoothly to an open set $\Omega^+$ that contains the smooth points of the boundary of $\Omega$ (but not the corners, in general). In fact, you don't even need that the boundary be polygonal; just being a finite union of smooth, embedded curves that meet at 'corners' is sufficient to construct an example.
Part II (The PDE system):
Any (smooth) metric $g$ admits unit divergence-free vector fields locally: If $p$ is a point in $\Omega$, let $c:[-\epsilon,\epsilon]\to\Omega$ be an embedded smooth curve with $c(0)=p$ and let $C:[-\epsilon,\epsilon]\times[-\delta,\delta]\to\Omega$ be defined so that, for each $t\in[-\epsilon,\epsilon]$, the curve $c_t:[-\delta,\delta]\to\Omega$ defined by $c_t(s) = C(t,s)$ is the unit speed geodesic satisfying $c_t(0)=c(t)$ and with initial velocity $c'_t(0)$ perpendicular to $c'(t)$. Then, for $\delta>0$ sufficiently small, $C$ is a smooth embedding of a rectangle into $\Omega$, and the $1$-form $\mathrm{d}s$ is a closed unit $1$-form on the image of $C$. As such, the metric $g$ can be written as $g = \mathrm{d}s^2 + \alpha^2$ for some $1$-form $\alpha$. Now, let $X$ be the vector field that satisfies $\mathrm{d}s(X) = 0$ and $\alpha(X)=1$. Then $X$ is a unit vector field and, since the area form is $dA = \mathrm{d}s\wedge\alpha$, we can compute that $$ \mathrm{div}_g(X)\,dA = \mathscr{L}_X\,dA = \mathrm{d}(\iota_X\,dA) = \mathrm{d}\bigl(\iota_X\,(\mathrm{d}s\wedge\alpha)\bigr) = \mathrm{d}(-\mathrm{d}s) = 0. $$ Thus, $X$ is divergence-free and unit size with respect to $g$. (This is the standard argument in dimension $2$.)
Conversely, if $X$ is a unit, divergence-free vector field on $\Omega$, then $\sigma = -\iota_X\,dA$ is a closed unit $1$-form on $\Omega$, and, since $\Omega$ is simply-connected, $\sigma = \mathrm{d}s$ for some function $s$. Moreover, the leaves of the foliation by the curves orthogonal to $X$ are geodesics that are calibrated by $\mathrm{d}s$.
In particular, if $g$ is smooth and positive definite along the 'edges' (except possibly at the corners) and $X$ extends smoothly to the 'edges' (but not necessarily to the corners), then the edges have to be geodesics in the $g$-metric. (If $g$ goes to zero on the edges, then $X$ has to blow up as you approach the boundary, so I'm not sure what it would mean for $X$ to be 'orthogonal' to the boundary.)