Assumptions
As discussed, some assumptions will be made. You mentioned an illiquid stock in your question. I'll get to that, but let's first assume that the stock price $S_u$ satisfies the following SDE under the pricing (risk-neutral) measure:
$$
dS_u = \left[r(u)-q(u)\right]S_u \,du + \sigma_u S_u \left[\rho \,dW_u + \sqrt{1-\rho^2} \,dZ_u\right]
$$
Here $r(u)$ and $q(u)$ are deterministic functions of time, and $W_u$ and $Z_u$ are standard independent Brownian motions.
The instantaneous volatility $\sigma_u$ is assumed be correlated to the stock price and that it satisfies the SDE
$$
d\sigma_u = a(\sigma_u,u)\,du + b(\sigma_u,u)\, dW_u,
$$
where $a(\sigma_u,u)$ and $b(\sigma_u,u)$ deterministic functions of $\sigma_u$ and $u$ only (they do not depend on $S_u$). We won't actually need to further specify these functions $a$ and $b$. The instantaneous volatility can also be driven by fractional noise (eg rough volatility model) without affecting the results below. But for now it suffices to let $\sigma_u$ to be driven by $W_u$ as per above SDE.
For $t\in [0,T]$ the volatility swap price is defined as
$$
v_t := e^{-r(T-t)} E_t \left( \frac{1}{T} \int_0^T \sigma^2_u du \right)^{1/2},
$$
and the variance swap price is given by
$$
V_t := e^{-r(T-t)} E_t \left( \frac{1}{T} \int_0^T \sigma^2_u du \right).
$$
Since both $v_t$ and $V_t$ are discounted expectations of a claim at $T$, both their risk-neutral drifts equals $r$. I am going to make the assumption now that the volatility swap price is lognormal, and that it is driven by the same standard Brownian motion that drives $\sigma_u$. The SDE for $v_u$ is then
$$
dv_u = r(u) v_u \,du + \alpha(u) v_u \, dW_u,
$$
with $\alpha(u)$, which the volatility of the volswap price, a deterministic function of time $u$ and potentially of time to maturity $T$. I don't need to specify the dynamics for $V_u$ as will become clearer shortly.
Now that the main assumptions have been stated, I am going to set $r=q=0$ for simplicity in order to make notation less cluttered. It shouldn't be too difficult for you to slightly modify the results below for nonzero discount rate and dividend yield.
Options on realised volatility
First of all notice that the realised volatility over $[0,T]$ is simply the terminal price of the volatility swap:
$$
v_T = \left( \frac{1}{T} \int_0^T \sigma^2_u du \right)^{1/2}.
$$
Hence, an option on realised volatility is equivalently an option on the terminal price of the volswap. The price of an option on realised volatility with strike $K$, given today's volatility swap price $v_0$, is therefore
$$
F(v_0,K,T) = E_0 \left[ \left (v_T - K \right)_+ \right].
$$
The solution to the SDE for the volswap is
$$
v_T = v_0 e^{ -\frac12 \bar\alpha^2 T + \bar\alpha W_T},
$$
where
$$
\bar\alpha^2 := \frac1T \int_0^T \alpha(u) \,du.
$$
It remains to find $v_0$ and $\bar\alpha$. Once they are known the price of the option on realised volatility is given by the Black-Scholes call formule with spot price $v_0$ and volatility $\bar\alpha$.
Backing out $v_0$ and $\bar\alpha$ from vanilla stock options
The more difficult part is finding $v_0$ from vanilla stock option prices. It is possible however to find $v_0$ in a relatively model-free manner (i.e. only using vanilla stock options). One way is by using Carr and Lee's correlation immunisation strategy. For this you would need a continuous strip of stock options, which is challenging if the stock is illiquid. Another way is to use the `zero vanna approximation' explained here. For that approach you only need one particular implied volatility that is usually close to (but not equal to) the stock ATM implied volatility. Although the difference in approximation accuracy between Carr and Lee and zero vanna has not been systematically studied, numerical experiments I've carried out so far indicates that the differences are really very small.
The next step is to find $\bar\alpha$. Because of the lognormal assumption for the volswap price,
$$
E_0 [V_T] = E_0[v_T^2] = v_0^2 e^{\bar\alpha^2 T}.
$$
Since $E_0 [V_T]$ is just the variance swap strike $V_0$,
$$
\bar\alpha^2 T = \log (V_0/v_0^2).
$$
So if you have $V_0$ then, since you also have an approximation for $v_0$, you have an approximation for $\bar\alpha$.
For liquid stocks, $V_0$ can be backed out using eg the Matytsin formula which uses a continuous strip of stock option implied volatilities with strikes ranging from $0$ to $\infty$. However, for illiquid options, you might want to use the method outlined here which can approximate the volswap price $v_0$ and the varswap price $V_0$ simultaneously using simple matrix inversion and only 3 near the money stock options.
I am not going to explain the papers that I have referenced, otherwise this answer will become a paper. However feel free to ask more and related questions.
Hope this helps.
Last remark: Both Carr and Lee and the zero vanna IV approach for approximating the volswap price tend to slightly underestimate the volswap price relative to its exact model price (approximation error depends on the `exact' model). The upside of this is that the vol of vol is then slightly over-estimated, which if you're a bank selling voloptions is what you want, and if you're a client buying a voloption you know whereabouts the bak's ask-price should be.
