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.)

Now that the question has been changed so extensively, the remarks that I made for the old version are no longer of any interest. Here is what I understand the problem to look like now:

First, $\Omega\subset\mathbb{R}^2$ is a convex open set in the plane whose boundary is polygonal, i.e., a union of line segements or rays that meet at a discrete set of 'corners'. In what I am going to discuss, these hypotheses seem a bit restrictive. Most or all of what I have to say would apply if $\Omega$ were simply-connected (and hence contractible) and its boundary $\partial\Omega$ minus a discrete set $C\subset\partial\Omega$ were a union of embedded (smooth) arcs.

Second, we are given a (smooth) function $\lambda>0$ on $\Omega$ that limits to $0$ at every smooth point of $\partial\Omega$. (We don't further specify the behavior of $\lambda$ near the corner points $C\subset\partial\Omega$.) We use $\lambda$ to define a metric $g = \lambda\,(\mathrm{d}{x_1}^2+\mathrm{d}{x_2}^2)$.

Then, we ask the following question: When does there exist a (smooth) foliation $\mathcal{F}$ by curves on an open set in $\mathbb{R}^2$ containing $\Omega\cup (\partial\Omega\setminus C)$ with the following two properties: First, in the interior of $\Omega$, the leaves of $\mathcal{F}$ are $g$-geodesics. Second, for each point $p\in \partial\Omega\setminus C$, the boundary of $\Omega$ and the leaf of $\mathcal{F}$ through $p$ meet orthogonally (at $p$).

Here are a few examples:

1. Let $\Omega$ be the first quadrant, i.e., defined by $x_1>0$ and $x_2>0$. Then $\partial\Omega$ is the union of two rays, the positive $x_i$-axes, which meet at right angles at the origin (which is the unique corner point). Let $\lambda = {x_1}^2{x_2}^2({x_1}^2{+}{x_2}^2)$. Then a solution to the problem is defined by letting $\mathcal{F}$ be the foliation given by the level sets of the function $F = {x_1}^2-{x_2}^2$.

2. Let $\Omega$ be defined by $x_1>0$ and $x_2>{x_1}^2$. Then $\partial\Omega$ is the union of a ray (the positive $x_2$-axis) and the right half of the parabola $x_2-{x_1}^2=0$. These meet at right angles at the origin (which is the unique corner point). Let $\lambda = {x_1}^2(x_2{-}{x_1}^2)^2({x_1}^2{+}(x_2{-}3{x_1}^2)^2)$. Then a solution to the problem is defined by letting $\mathcal{F}$ be the foliation given by the level sets of the function $F = \mathrm{e}^{6x_1}(1-6x_1+18x_2)$.

3. Again, let $\Omega$ be the first quadrant, i.e., defined by $x_1>0$ and $x_2>0$. Then $\partial\Omega$ is the union of two rays, the positive $x_i$-axes, which meet at right angles at the origin (which is the unique corner point). Let $\lambda = {x_1}^2{x_2}^2/({x_1}^2{+}{x_2}^2)^{3}$. Then a solution to the problem is defined by letting $\mathcal{F}$ be the foliation given by the level sets of the function $F = {x_1}^2+{x_2}^2$. (Some interesting features of this example will be explained below, but note, already, that $\lambda$ does not limit to $0$ as one approaches the corner, although it does limit to $0$ as one approaches any smooth boundary point.)

To understand how these examples were constructed and to get an idea of some of the restrictions that the existence of a foliation $\mathcal{F}$ satisfying the two conditions places on the function $\lambda$, suppose that we have a given pair $(\Omega,\lambda)$ satisfying the above conditions and that there exists a foliation $\mathcal{F}$ with the desired properties.

Since $\Omega$ is contractible, we can write $g$ on $\Omega$ in the form $g=\eta_1^2+\eta_2^2$, where $\eta_1$ and $\eta_2$ are $1$-forms and where the leaves of $\mathcal{F}$ are the null curves of $\eta_2$. Then the condition that these leaves be $g$-geodesics is equivalent to the condition that $\mathrm{d}\eta_1 = 0$. Since $\Omega$ is simply-connected, there will exist a function $S:\Omega\to\mathbb{R}$, unique up to an additive constant, such that $\mathrm{d}S = \eta_1$.

Write $\mathrm{d}S = S_1\,\mathrm{d}x_1 + S_2\,\mathrm{d}x_2$. Since $g = \lambda\,({\mathrm{d}x_1}^2+{\mathrm{d}x_2}^2)=\eta_1^2+\eta_2^2$, it follows that $\eta_2 = \pm (S_2\,\mathrm{d}x_1 - S_1\,\mathrm{d}x_2)$, so we have $$ g = \bigl(S_1^2+S_2^2\bigr)\,({\mathrm{d}x_1}^2+{\mathrm{d}x_2}^2) = \lambda\,({\mathrm{d}x_1}^2+{\mathrm{d}x_2}^2). $$ In other words, $\lambda = \bigl(S_1^2+S_2^2\bigr)$. From this, it is not difficult to see that $S$ extends continuously to the smooth part of $\partial\Omega$ and is constant on each connected component of the smooth part of the boundary, while it has no critical points within $\Omega$ itself.

Because of the hypothesis that the foliation $\mathcal{F}$ extends smoothly to an open set $U$ containing $\Omega$ and the smooth part of its boundary, it follows that there is a smooth function $\theta:U\to S^1$ such that the leaves of $\mathcal{F}$ in $U$ are the null curves of the smooth $1$-form $\sin\theta\,\mathrm{d}x_1-\cos\theta\,\mathrm{d}x_2$ and so that $\eta_2 = \sqrt{\lambda}\bigl(\sin\theta\,\mathrm{d}x_1-\cos\theta\,\mathrm{d}x_2\bigr)$ in $\Omega$, and this extends at least continuously to the smooth part of the boundary of $\Omega$ by extending $\lambda$ to be $0$ on the smooth part.

Conversely, the function $S:\Omega\to\mathbb{R}$ gives enough information to recover both the metric $g$ and the foliation $\mathcal{F}$. If, as the OP has asked, one has $\lambda$ already specified, then $S$ must satisfy the first-order scalar PDE $|\nabla S|^2 = \lambda$ on $\Omega$ plus the smooth part of the boundary.

Sometimes, such an $S$ exists (as the examples above show, although they were found by starting with a candidate $S$), and sometimes it does not. I suspect that the OP would like an effectively computable necessary and sufficient condition for a given $\lambda$ to possess an admissable solution $S$, but that is likely to be hard to determine, because it really depends on delicate information about the geodesic flow in $\Omega$ endowed with the metric $g$.

For an example of an definitive condition, suppose that there is a simple closed $g$-geodesic $\gamma$ in $\Omega$. Then $\gamma$ cannot be a leaf of any geodesic foliation $\mathcal{F}$ of $\Omega$, since the function $S$ that corresponds to $\mathcal{F}$ could not then have a critical point on $\gamma$, which it must. Thus, $\gamma$ would have to be transverse to all of the leaves of $\mathcal{F}$, and we would then have a foliation of the disk $D\subset\Omega$ bounded by $\gamma$ that is transverse to the boundary of $D$, an obvious impossibility (by degree theory). (Probably, this argument works even if $\gamma$ is closed but not simple, but I leave that for the interested.)

It is certianly possible to chose an admissable $\lambda$ for which the metric $g$ would have a simple closed geodesic, thus showing that there are $\lambda$ for which there is no $g$-geodesic foliation, and hence no function $S$ with the desired properties. However, for a given $\lambda$, it is not clear how to tell whether there exists a simple closed $g$-geodesic in $\Omega$. There are checkable conditions that guarantee one and checkable conditions that forbid one, but the gap between these is large.

I am dubious that an effectively checkable necessary and sufficient condition for the solvability of the problem of given $\lambda$ exists.

Actually, it's easy to solve this problem for 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.)

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.)