Chi square approximation of the likelihood test ratio I wasn't able to find any satisfying answer about that topic. I hope someone who understand correctly the subject could enlighten this shadow.

This is not very important, just for the sake of notations: In the
  following we work with exponential families:  $$ \exp \{ \frac {y
 \theta - b( \theta ) } { a(\phi)} + c(y,\phi) \} $$ the deviance is
  defined as the difference in likelihood between the model and the saturated model:  $$D = 2 \phi \left ( l(y, \phi , y ) - l ( \hat{\mu},
 \phi, y ) \right ) $$ the scaled deviance as : $$ D^* = \frac D {\phi} $$

When one compares two nested models, in particular we are trying to check the goodness of fit, and we are trying to test the null hypothesis (smaller model is doing a better approximation), we can compute $D_0^* - D_1^* $ where $D_0^*$ is the scaled deviance of the smaller model. This, should be approximately, be following a $ \chi^2 $ distribution, of the difference in number of parameters.
My question is the following, how to test whether that approximation is good ? My professor is not very precise on that topic, he mentions doing "simulations". What does it mean please?
 A: That is a large topic, and the quality of the approximation must be studied, generally, on a case-by-case basis. That's why simulation is a useful approach. So the best way to answer your question is by an example, which you can adapt for your data and models.
So I will simulate, in R, some data for a logistic regression, and I will simulate binomial responses from the null model, that is, all regression parameters (except the intercept) assumed to be zero. Then I fit the logistic model, calculate the deviance, repeat many times, and plot the histogram, with the chi-square approximation overlaid.  For my case the result is (based on 10000 replications)

which looks rather good. A qqplot could also be useful:

This looks like a very good fit. Maybe some clumpiness (discreteness), but no systematic deviations. But also, the far tail is poorly represented, to check quality of fit there, maybe a larger simulation.
You can modify the code below for your experiments:
set.seed(7*11*13)  # My public seed
X <- data.frame(x1=rnorm(200, 20, sd=3),
                x2=factor(rep(1:3, c(100, 50, 50))),
                y=rbinom(200, 1, 0.25))

N <- 10000 # Simulating the null distribution

sim_dev <- function(N) {
    res <- numeric(N)
    for (i in seq_along(res)) {
        X$y <- rbinom(200, 1, 0.25)
        mod <- glm(y  ~ x1+x2, data=X, family=binomial)
        res[i] <- with(mod,  null.deviance-deviance)
    }
    res
}

deviance <- sim_dev(N)

mean(deviance); var(deviance)
 3.016633
 6.090653  # Rather close to theoretical values 3, 6

hist(deviance, prob=TRUE, breaks="FD")
plot( function(x) dchisq(x, df=3), col="red", add=TRUE, 
                         from=0, to=22) 
qqplot(qchisq(ppoints(N), df=3), deviance)
qqline(deviance, distribution=function(x)qchisq(x, df=3))
title("qqplot against chisquare distribution (df=3)") 

