Most perturbations of overdetermined, involutive systems, such as $\mathrm{d}F = \delta\,F = 0$, have no nontrivial solutions.
For example, consider the first nontrivial case: A $2$-dimensional (oriented) Riemannian manifold $(M^2,g)$ with a nonvanishing vector field $X$. Locally can choose an oriented orthonormal coframe field $(\omega_1,\omega_2)$ so that $g = {\omega_1}^2 + {\omega_2}^2$ and $\omega_2(X)=0$ while $\omega_1(X) = x$ for some function $x>0$. Write $F = f\,\omega_1\wedge\omega_2$. Of course, $\mathrm{d}F=0$ is an identity, but one has
$$
(\delta - \iota_X)F = f_2\,\omega_1 - (f_1{+}xf)\,\omega_2
$$
where $\mathrm{d}f = f_1\,\omega_1 + f_2\,\omega_2$, so this is two equations
$$
f_1 + xf = f_2 = 0.
$$
However, since $\delta^2 = 0$, we also get $\delta(\iota_XF) = -(f_2x + f x_2) = 0$ (where $\mathrm{d}x = x_1\,\omega_1 + x_2\,\omega_2$), which, together with the above equations becomes $x_2\,f = 0$. Thus, if $x_2$ is nonvanishing, we must have $f = 0$, i.e., $F = 0$ as the only solution.
Now, in the special case that $x_2$ vanishes identically, the above equations become
$$
\mathrm{d}f = -fx\,\omega_1\qquad\text{and}\qquad \mathrm{d}x = x_1\,\omega_1\,.
$$
If there is any nonzero solution $f$ to the first equation, then we must have $\mathrm{d}(x\omega_1) = 0$, and, coupling this with the second equation, we find that $0 = \mathrm{d}(x\omega_1) = x\,\mathrm{d}\omega_1$, so, since $x$ is nonvanishing, we must have $\mathrm{d}\omega_1 = 0$. In particular, there must exist local coordinates $(u,v)$ so that $\omega_1 = \mathrm{d}u$ and $\omega_2 = h(u,v)\,\mathrm{d}v$. Now, the local solutions can all be written in the form
$$
F = c\,\mathrm{e}^{a(u)}\,h(u,v)\,\mathrm{d}u\wedge\mathrm{d}v,\qquad
X = a'(u)\,\frac{\partial}{\partial u},\qquad
g = \mathrm{d}u^2 + h(u,v)^2\,\mathrm{d}v^2,
$$
where $a'(u) > 0$ and $c$ is a constant.
A similar analysis can be done in any dimension, and you'll get nontrivial conditions on the pair $(g,X)$ in order for there to be nontrivial solutions. In general, though, the 'extra' first order equation $\delta(\iota_XF)=0$ on the unknown $2$-form $F$ is not compatible with the (already overdetermined) system $\mathrm{d}F = (\delta - \iota_X)F = 0$.