As far as I know, Dynare only supports stochastic simulation based on normally distributed shocks, since we only specify the standard deviations in “shocks” block. Is there any way to simulate the model with non-normal shocks?

I know that some simple non-normal distribution can be expressed in terms of normal distribution, e.g. log-normal distribution. However, in my model, I have some shocks following beta distribution, which is a more general distribution which cannot be expressed in terms of normal distribution. Though, normal distribution is a special case of beta distribution with both two parameters = inf.

The problem is not really that Dynare does only support normally distributed shocks but that perturbation is used to solve the models and up to first order, only the mean and the variance of the shock distribution play a role. Thus, even a beta distribution is effectively approximated with a normal distribution for computing the model solution, because the distribution is fully characterized by its first two moments. A limitation of Dynare currently is that it does not allow for skewed distributions, i.e. even at order 3 the skewness of the shock distribution is assumed to be 0.

I am not sure what you are exactly trying to do. If you are OK with having a solution based on the first two moments of a beta distribution and just want to generate simulations with shock drawn form the beta distribution, you can hack Dynare’s corresponding routines to use shocks drawn from a different distribution.
See