Consider a finite sum of a function $f(x)$ over discrete values of $x$. $$S=\sum_{x=a}^{b} f(x)\tag{1}$$
Now suppose that, instead of having only certain values of $i$, this variable can vary continuously in the interval $[a,b]$, i.e. $x \in [a,b] \subset \mathbb{R}$.
In many occasion, studying physics mainly, I read on textbooks that such sum can be in that case "approximated" with an integral. Nevertheless I'm quite sure that I cannot write
$$S \approx \int_{a}^b f(x) dx\tag{2}$$
Since this would not work dimensionally. I could imagine to divide by a constant factor for dimensions but I don't think it would work either, because the integral of a function itself is proportional to its integral average on the interval $[a,b]$, which I do not think can be proportional to the "sum" $S$.
So I can imagine to divide $[a,b]$ in sub-intervals and creating a density function $\frac{dn}{dx}$ which sholud give the "relative number of allowed values of $x$ in the interval $dx$" and then finally
$$S \approx \int_{a}^b \frac{dn}{dx} f(x) dx\tag{3}$$
But again if there are no reason for $x$ not to be spread uniformly I should set $\frac{dn}{dx}=constant$ and I would get the same of above.
So what integral should I think of when I read on a physics book that a discrete sum like $(1)$ can be subsituted by an integral if the variable is allowed to vary continously?
If you insist on everything proven rigorously, you could raise a number of objections to what follows below (and some have already done so in comments to your question), but since you referred to applications in physics, here's how physicists justify that. In a nutshell, the answer is: the correct version is (3).
Split up the interval $[a,b]$ into $k$ equal pieces of length $\Delta x = \frac{b-a}{k}$. Let's agree that when we write $\sum_{x = a}^b$ we mean that $x$ takes the discrete allowed values in $(a, b]$, to avoid double inclusion of any that fall exactly on a boundary. So $$ \sum_{x=a}^b f(x) = \sum_{j=0}^{k-1}\left[ \sum_{x = a+j\Delta x}^{a+(j+1)\Delta x} f(x)\right]. $$ Now, if $f(x)$ doesn't change much as $x$ takes the different values in $(a+j\Delta x,a+(j+1)\Delta x]$ we can approximate the sum for each sub-interval by the number of terms in it, $n_j$, times an average, $f(x_j)$ for some $x_j \in (a+j\Delta x,a+(j+1)\Delta x]$. (Actually, for a continuous function, this is not an approximation at all, it is always possible to pick an $x_j$ so that it is exact, and that doesn't depend on how much $f$ varies over the interval). Introducing a cumulative function for the number of allowed values $x \leq s$, $N(s)$, so that $n_j = N(a+(j+1)\Delta x)-N(a+j\Delta x)$, we have $$ \sum_{x=a}^b f(x) \approx \sum_{j=0}^{k-1} \left[\frac{N(a+(j+1)\Delta x)-N(a+j\Delta x)}{\Delta x} f(x_j) \right]\Delta x. $$ This will be a Riemann sum if we can also give the big fraction the meaning of a density, $dN(x)/dx$ as $\Delta x$ becomes small. The trouble is that $N$ isn't really continuous (in other words, we don't really want $\Delta x$ approaching zero because of the discreteness), but supposing this can be done (at least as an approximation), we can then write $$ \sum_{x=a}^b f(x) \approx \int_a^b \frac{dN(x)}{dx}f(x)dx. $$ This all depends on being able to choose sub-intervals small enough so that the integral approximates the Riemann sum well, but large enough so that the "density" can be reasonably approximated by some continuous function $dN/dx$. If the discrete values of $x$ are equally spaced, the latter is not a problem; we can pick $\Delta x$ to be the spacing and then $dN/dx = 1/\Delta x$; this will generally work so long as $\Delta x \ll b-a$ (from the first condition).
One common example of this is an approximation for $\ln(n!)$ for large integers often found in statistical physics books: $$ \ln(n!) = \sum_{x=1}^n \ln x \approx \int_1^n \ln x dx = n\ln n -n + 1 $$ where $x$ takes integer values, so we use $\Delta x = 1$ and $dN/dx = 1$ (see the Wikipedia article Stirling's approximation for more technical details on the validity of this). Now let's try an example where the $x$'s aren't evenly spaced. Let $x$ take values from the set $\{\sqrt{n}|n\in\mathbb{N}\}$ and let $f(x) = e^{-x^2}$ (yes, this is rather contrived; we might as well do this one with uniformly spaced values and $f(x) = e^{-x}$ instead, but this is in the interest of being able to evaluate the sum algebraically). Suppose we want to sum from $x = 10$ to $20$.
To illustrate the whole process, let's split the interval into 10 parts, so that $\Delta x = 1$. What is $n_j$? The last value that is no greater than an integer $s$ is $\sqrt{s^2}$, thus there are $s^2$ values of $x \leq s$, so $N(s) = s^2$. Then $n_j = N(10 + (j+1)) - N(10+j) = (11+j)^2 - (10+j)^2 = 21+2j$. On the other hand, extending the expression to $N(x) = x^2$ for a continuous variable $x$ produces $dN/dx = 2x$. For any $x_j \in (10+j, 10+j+1]$ the difference $|n_j - (dN/dx)_{x=x_j}| \leq 1 \ll n_j$. So this step of the approximation works ok.
Let's calculate then: $$\int_{10}^{20} \frac{dN}{dx}f(x)dx = \int_{10}^{20} 2xe^{-x^2}dx = \int_{100}^{400}e^{-u}du \approx e^{-100},$$ where I have omitted the negligible in comparison $e^{-400}$. The sum, on the other hand, is a geometric series: $$\sum_{x=10}^{20} e^{-x^2} = \sum_{n=100}^{400} e^{-n} \approx e^{-100}\frac{1}{1-1/e},$$ again neglecting a very small term in the numerator. These are significantly different results, and the reason here is that with this function and this partition, the Riemann sum is not a good approximation for the integral. But we can't make $\Delta x$ much smaller, because then the approximation from the first step, $n_j \approx dN/dx$ will be poor. The discrete values simply aren't dense enough in this example. If we take instead $x$'s from the set $\{\sqrt{n/10}|n\in\mathbb{N}\}$, we will get much better agreement (note $N(x) = 10x^2$ now): $$ \int_{10}^{20} 20x e^{-x^2} dx \approx 10e^{-100},\\ \sum_{x = 10}^{20} e^{-x^2} = \sum_{n=1000}^{4000} e^{-n/10} \approx \frac{e^{-100}}{1-e^{-0.1}} \approx 10.51e^{-100}. $$