The reason why the force term is neglected is because it is often not needed and it would only complicate the procedure of deriving the macroscopic equations, the Chapman-Enskog expansion. The force term can be easily added afterwards to the incompressible lattice-Boltzmann equation by introducing a clever additional term that will mimic the required force in the relevant equations. I will show you in the last section of my answer how this can be done.
The incompressible lattice-Boltzmann equation can be used to simulate dense fluids as the algorithm can be shown to preserve the incompressible Navier-Stokes equations that take the same form regardless of the fluid (as long as it is Newtonian but there exist even artificial fixes for non-Newtonian fluids), for incompressible gases as well as liquids.
In the following section I will explain you how this works in a continuous setting and how this can be translated to a discrete setting, namely the lattice-Boltzmann methods. I have written my Master's thesis in the field. So feel free to read it as a introduction to the topic, I have put a lot of effort in making it easy to understand. The following section will include several formulas but most of them are just for giving you an idea what this looks like and it is not required to understand them thoroughly in order to follow my comment. I have also linked several of my previous comments on the topic I have given over here. In order to be consistent with my previous comments I will use my own nomenclature.
Boltzmann equation
Around 150 years ago a beardy Austrian guy named Ludwig Boltzmann came up with the Boltzmann equation
$$ \underbrace{ \frac{\partial f}{\partial t} + \vec \xi \boldsymbol{\cdot} \vec \nabla f + \frac{\vec F}{m} \boldsymbol{\cdot} \vec \nabla_{\vec \xi} f }_\text{Propagation} = \underbrace{ \Omega(f) }_\text{Collision} $$
that describes the evolution of a rather unintuitive particle distribution function $f$
$$ f := \frac{d N}{d \vec x \, d \vec \xi}$$
where $N$ is the number of particles. This equation tells you how a distribution function for dilute (mono-atomic) gases evolves over time due to collisions and free motion. Such a distribution function $f$ can be seen as an extended concept of density $\rho$, a density not only in space $\vec x$ but also in velocity space $\vec \xi$. The description seems very complicated with a six-dimensional phase space (in 3D: 3x space and 3x velocity) but the advantage of such a description is that it is more general than continuum-based Navier-Stokes-Fourier equations you are certainly familiar with. Continuum mechanics rests on the assumption that for every point in space and time you can find a suiting macroscopic description given by macroscopic variables such as density, velocity, pressure and viscous stresses. This excludes flow regimes where flow is very dilute, such as in shock-waves around supersonic aircraft or in case of the re-entry of a space shuttle, where the fluid is so dilute atoms can travel a significant distance without collisions. The Boltzmann equation on the other hand still holds in such a case as it considers those collisions. In its general form the collision operator takes following form
$$ \Omega_{Stoß} = df |_{\text{collision}} = \frac{\Delta N_{\text{gain}} - \Delta N_{\text{loss}}}{\Delta V \Delta \vec \xi \Delta t} = \int\limits_{ \vec \xi_1 } \int\limits_{ A_c } |\vec g| (f_1' f' - f f_1 ) d A_c d \vec \xi_1. $$
This basically renders the Boltzmann equation in this form a integro-differential equation as it contains integrals as well as differentials in the variable $f$.
Maxwell-Boltzmann equilibrium distribution
One thing one can derive easily for such a model is what such a distribution $f$ might look like at equilibrium. It takes the form of a Gaussian
$$ f^{(eq)} = n \mathcal{N} (\vec x, \vec \mu, \sigma) = \frac{n}{(2 \pi)^\frac{d}{2} \sigma^d} e^{- \frac{1}{2} \left( \frac{\vec x - \vec \mu}{\sigma} \right)^2}$$
where $n$ is the number density
$$n := \frac{N}{V} = \frac{\rho}{m_P}.$$
Furthermore it can be shown that every distribution would tend to such an equilibrium distribution if external mechanisms such as forces and boundary conditions are removed. As a result one may come up with the idea to replace the traditionally quite involved collision operator with a simple relaxation towards the equilibrium. This is essentially what the BGK collision operator does.
$$ \Omega_{BGK} := \frac{f^{(eq)} - f}{\tau}$$
All the following steps may be performed with the original collision operator as well as with the simplified BGK-operator.
From dilute gas to dense fluids
What the effect of the Boltzmann equation in terms of macroscopic variables precisely is, is not immediately visible: The microscopic distributions are something rather obscure. All we know is that this equation will preserve mass, momentum and energy during collisions but that's about it.
Nonetheless the distribution functions are connected to the macroscopic view through moments. All the particles in a particular point in space must have the same density, momentum and energy as a continuum would have at this position. Thus, we form the corresponding differential contributions and integrate over the velocity space to yield e.g. density and momentum
$$\rho = m_P \int f d \vec \xi, \qquad \rho \vec u_i = m_P \int \xi_i f d \vec \xi.$$
Now developing the unknown distribution $f$ in a perturbation series in a small parameter$\epsilon$ (the later mentioned Knudsen number $\epsilon = Kn$), a procedure similar to a Taylor series, around a Maxwellian equilibrium $f^{(eq)}$
$$f(t_0, t_1, t_2, \ldots) = \underbrace{f^{(0)}(t_0)}_{f^{(eq)}} + \epsilon f^{(1)}(t_1) + \epsilon^2 f^{(2)}(t_2) + \mathcal{O}(\epsilon^3)$$
where
$$ t_n = \epsilon^n t_0$$
and considering the corresponding moments of the Boltzmann equation, one can derive the macroscopic equations that the Boltzmann equation preserves on longer time scales. Surprisingly through this procedure it can be found that in the limit of a continuum (small Knudsen numbers) these are the full Navier-Stokes-Fourier equations for a compressible gas. This basically means the kinetic theory of gases in the limits of dense gases preserves the Navier-Stokes-Fourier equations.
As you are certainly aware flow of liquids and gases in continuum mechanics is generally quite similar, both are governed by a system of partial differential-equations, the Navier-Stokes-Fourier equations
$$\frac{\partial \rho}{\partial t} + \sum\limits_{j \in \mathcal{D}} \frac{\partial (\rho u_j )}{\partial x_j }=0,$$
$$\frac{\partial (\rho u_i )}{\partial t} + \sum\limits_{j \in \mathcal{D}}\frac{\partial (\rho u_i u_j )}{\partial x_j} = \sum\limits_{j \in \mathcal{D}} \frac{\partial \sigma_{ij}}{\partial x_j } + \rho g_i,$$
$$\frac{\partial (\rho e)}{\partial t} + \sum\limits_{j \in \mathcal{D}} \frac{\partial (\rho u_j e)}{\partial x_j} = - \sum\limits_{j \in \mathcal{D}} \frac{\partial q_j}{\partial x_j} + \sum\limits_{i, j \in \mathcal{D}} \frac{\partial (\sigma _{ij} u_i)}{\partial x_j} + \sum\limits_{j \in \mathcal{D}} \rho u_j g_j.$$
Liquids generally can't be compressed though. This is then considered in the equation of state and the Navier-Stokes-Fourier equations may be simplified. While generally gases are not incompressible under certain assumptions even in a gas density fluctuations can be negligible. In this case both, a gas and a liquid are governed by the identical partial differential equations. The only difference is their equation of state (for liquids something like the Tait's equation while for most gases the ideal gas law suffices). This means if choosing similar dimensionless parameters such as the Reynolds number wisely one could even simulate many phenomena with a gas instead of a liquid.
From continuous to discrete: lattice-Boltzmann equation
The previous procedure can be shown for the discrete lattice-Boltzmann equation as well. In this case though the moments will degrade to discrete sums.
Sadly the discretisation you choose (the number of velocities) is connected to artifacts and the preserved differential equations. This means that with a reasonable number of discrete velocities you are unable to preserve the full Navier-Stokes equations. If you choose to follow the derivation through with a single possible velocity (the Boltzmann distribution on the other hand is continuous and has indefinitely many possibilities) a time and space discretisation conveniently chosen as $\Delta x = 1$ and $\Delta t = 1$, you will be left with "handicapped" Navier-Stokes equations with certain artifacts:
$$ \frac{\partial \rho}{\partial t} + \sum\limits_{j \in \mathcal{D}} \frac{\partial (\rho u_j)}{\partial x_j} = 0 $$
$$ \frac{\partial (\rho u_i)}{\partial t} + \sum\limits_{j \in \mathcal{D}} \frac{\partial (\rho u_i u_j)}{\partial x_j} = - \frac{\partial}{\partial x_i} \left( \frac{c^2 \rho}{3} \right) + \sum\limits_{j \in \mathcal{D}} \frac{\partial}{\partial x_j} \left[ \left( \frac{\tau}{\Delta t} - \frac{1}{2} \right) \left( \frac{c^2}{3} \rho \left(\frac{\partial u_j}{\partial x_i} + \frac{\partial u_i}{\partial x_j} \right) - \sum\limits_{k \in \mathcal{D}} \frac{\partial}{\partial x_k} \left( \rho u_i u_j u_k \right) \right) \right] $$
This resembles the familiar continuity and momentum equations if the viscosity is chosen as
$$ \nu = \left( \frac{\tau}{\Delta t} - \frac{1}{2} \right) \frac{c^2}{3} \Delta t = \left( \frac{\tau}{\Delta t} - \frac{1}{2} \right) c_s^2 \Delta t $$
where due to
$$ - \frac{\partial}{\partial x_i} \left( \frac{c^2 \rho}{3} \right) = - \frac{\partial p}{\partial x_i} $$
pressure is linked to density via the equation of state
$$ p = c_s^2 \rho $$
and
$$ - \frac{\partial}{\partial x_k} \left( \rho u_i u_j u_k \right) $$
is an error term emerging from the discretisation that is though of order $\mathcal{O}(Ma^3)$ and hence can be neglected for small local Mach numbers, if $\vec u$ is comparably small to the local speed of sound $c_s$.
As you can see the equation of state is the one of an iso-thermal ideal gas as the speed of sound in that case is given by
$$c_s^T := \sqrt{\left( \frac{\partial p}{\partial \rho} \right)_T } = \sqrt{R_m T} .$$
Furthermore the equations are neither compressible nor truly incompressible as the corresponding dilational term is missing but the density is still inside the temporal and spatial derivative in the continuity and momentum equation. Finally the energy equation would not make sense at all and is thus omitted in the derivation.
The only smart thing you can do with these Frankenstein equations is expose their simple computational nature to simulate transient flows of incompressible fluids by letting the simulation Mach number go to zero (choosing very small simulation velocities). The conversion between the system is done by the law of similarity and setting similar dimensionless numbers in the fictional LBM.
Several authors do not recommend seeing lattice-Boltzmann as a child of the Boltzmann equation as historically it developed from the lattice gas automata (LGA) but I think this is a big mistake. Seeing the equation in the light of its significantly more capable bigger brother might lead to researchers proposing more methods based on kinetic theory that might harvest the enormous theoretical framework more efficiently and to a larger extent. I think actually one should think of lattice-Boltzmann as a fictional system where every the particle has a mass of $1 kg$, a time step is $1 s$ and the length step $1 m$. The velocity might be chosen arbitrarily as long as it is significantly smaller than the speed of sound (low Mach number). There is no direct way to change the temporal resolution directly instead one has to vary the velocity or rescale the system (more cells with a size of $1 m$ each). Thus, one should always list formulas with the space $\Delta x$ and time step $\Delta t$ even though they will be set to 1. These values are important for scaling and give the system its units.
Dimensionless Boltzmann equation
For rendering any equation dimensionless generally you will try to introduce as few variables as possible (so I'd avoid an $F_0$) and you will have to know what a good reference value is. This depends a lot on the flow regime and in this case on the type of force. As a choice for reference values for a force due to gravity (where $g = 9.81 \frac{m}{s^2}$ is the gravitational constant) you would use something like
$$g* = \overbrace{g}^{\frac{F}{m}} \frac{L}{U^2}.$$
In the dimensionless equation then characteristic dimensionless numbers should appear, in the case of the continuous dimensionless Boltzmann equation
$$f^*=f \frac{c_0^3 L^6}{n}, \qquad t^*=t \frac{c_0}{L}, \qquad x_i^*=\frac{x_i}{L}, \qquad \xi_i^*=\frac{\xi_i}{c_0}, \qquad g^*=g \frac{L}{c_0^2}, \qquad \tau^*=\tau \frac{c_0}{\lambda},$$
$$\frac{\partial f^*}{\partial t^*} + \xi_j^* \frac{\partial f^*}{\partial x_j^*} + g_j^* \frac{\partial f^*}{\partial \xi_j^*} = \frac{1}{Kn} \frac{1}{\tau^*} \left( f^{(eq)*} - f^* \right).$$
this is the Knudsen number
$$ Kn := \frac{\lambda}{L} \phantom{spacespace} \frac{\text{mean free path}}{\text{representative physical length scale}}.$$
For the discrete equation you now basically only consider certain directions $\alpha$ meaning. Meaning instead of a continuous variable $f$ we will consider a variable $f_\alpha$ with $\alpha \in \mathcal{L}$ directions. Similarly you would have to decompose your force but instead you might use the following tricks for most relevant forces.
Adding forces to the lattice-Boltzmann equation
Specific forces (such as gravity in the following case) do not pose particular requirements for the conserved momenta and therefore might be modelled through an additionally artificial forcing term on the right hand side
$$ - g_i \frac{(c_{\alpha i} - u_i)}{c_s^2} f_\alpha^{(eq)} $$
that does not contribute to the conservation of mass due to
$$ - \frac{g_i}{c_s^2} \sum_\alpha (c_{\alpha i} - u_i) f_\alpha^{(eq)} = - \frac{g_i}{c_s^2} \Big( \underbrace{\sum_\alpha c_{\alpha i}f_\alpha^{(eq)}}_{\rho u_i} - u_i \underbrace{\sum_\alpha f_\alpha^{(eq)}}_{\rho} \Big) = 0 $$
but leads to a corresponding term in the momentum equation
$$- \frac{g_i}{c_s^2} \sum_\alpha (c_{\alpha i} c_{\alpha j} - u_i c_{\alpha j}) f_\alpha^{(eq)} = - \frac{g_i}{c_s^2} \Big( \underbrace{\sum_\alpha c_{\alpha i} c_{\alpha j} f_\alpha^{(eq)}}_{\Pi_{ij}^{(0)} = c_s^2 \rho \delta_{ij} + \rho u_i u_j} - u_i \underbrace{\sum_\alpha c_{\alpha j} f_\alpha^{(eq)}}_{\rho u_j} \Big) = - \rho g_i \delta_{ i j}.$$
As only these low order moments are physically meaningful the contribution to higher order moments is not of relevance.
