Sampling of band-limited white noise The context is communication where we have a front-end that samples a signal and a noise (but here we focus only on the noise). My goal is to determine the noise power that I should use to simulate (in Matlab) the discrete noise after the ADC.
This post is related to this one, and I will refer to it. Below are figures to illustrate the discussion.

Let's consider a (virtual) continuous white noise of PSD $N_0/2$ (see Fig. 1). It passes through an ideal filter of bandwidth $[-B; B]$. So, we obtain a band-limited white noise, of power $N_0 B$ (see Fig. 2, with $B$ = 10 kHz). Here, we are just before the ADC of the front-end. After that, the signal is sampled, and here is how I see the things :

If the signal is sampled at a higher frequency than twice the bandwith (e.g. 30 kHz), the output is no more a white noise, since the PSD is not constant (see Fig. 3, in red is the duplicate of the spectrum around $±f_s$, so the actual spectrum is the grey curve plus the red curve).
Dilip Sarwate illustrates this very clearly using the autocorrelation function of the output noise in this post.

If the signal is sampled at a frequency equal to twice the bandwith (i.e. 20 kHz), the output is a white noise since the PSD is constant and equal to $N_0/2$ (see Fig. 4).
However, this disagrees with the development of Dilip Sarwate still in this post, where with 20 kHz we are not anymore on a zero of the sinc and thus the output samples are correlated...

If the signal is sampled at a frequency lower than the bandwith (e.g. 10 kHz), the output is also a white noise since the PSD is constant; but there is aliasing which increases the amount of noise and the PSD is now $N_0$ (Fig. 5). With lower sampling frequencies, the increase would be higher.

Therefore, for me, if I want to simulate in Matlab the discrete noise after the ADC, I should consider $f_s = 2 B$ to assume the noise as AWGN, and the noise power should be $\sigma^2 = N_0 B = N_0 f_s / 2$, whatever the actual front-end bandwidth is.

Could you tell me where my reasoning is wrong regarding the case of Fig. 4 ?
There are two things to consider here: the noise power after sampling, and the fact whether the sampled noise is uncorrelated or not. Note that the autocorrelation function of the sampled noise equals the sampled autocorrelation function of the continuous-time noise. For a wide-sense stationary (WSS) band-limited white noise signal $X(t)$ with bandwidth $B$ and one-sided noise density $N_0$, the autocorrelation function is given by

$$R_X(\tau)=BN_0\,\text{sinc}(2B\tau)\tag{1}$$

where I use $\text{sinc}(x)=\sin(\pi x)/(\pi x)$.

If you define a discrete time process $Y[n]=X(nT)$, where $f_s=1/T$ is the sampling frequency, then its autocorrelation function is given by

$$R_Y[k]=R_X(kT)=BN_0\,\text{sinc}(2BkT)\tag{2}$$

Since the power of the discrete-time noise $Y[n]$ is given by $R_Y[0]=BN_0$, it is obvious that the sampling frequency has no influence on the noise power. So no matter how fast you sample the original continuous-time noise, the power of the resulting discrete-time noise always equals the power of the noise before sampling.

Now let's have a look at the correlation between the noise samples. The autocorrelation function $R_X(\tau)$ has zeros at

$$\tau_i=\frac{i}{2B},\quad i=1,2,\ldots$$

which for $B=10\,\text{kHz}$ gives $\tau_i=i\cdot 50\mu s$ (this is the small numerical error in Dilip Sarwate's post you referred to in your question). This means that if you choose a sampling interval which is an integer multiple of $50\mu s$, the resulting sampled noise will be uncorrelated. This corresponds to sampling frequencies $f_s=2B/i$, $i=1,2,\ldots$, i.e. $20\,\text{kHz}$, $10\,\text{kHz}$, $6.67\,\text{kHz}$, etc.

So with a sampling frequency of $20\,\text{kHz}$ you do get uncorrelated noise samples.