Inductively coupled plasmas have been widely used in the etch process due to the high plasma density, simple reactor geometry, etc. Since the plasma characteristics are difficult to understand only via experiments, the numerical study seems to be a valuable and effective tool, which could help us to gain an in-depth insight into the plasma properties and the underlying mechanisms. During the past few years, various models have been employed to investigate inductive discharges, such as global model, fluid model, fluid/Monte Carlo collision hybrid model, biased sheath model, particle-in-cell/Monte Carlo collision hybrid model, etc. Since the plasma parameters are volume averaged in the global model, which effectively reduces the computational burden, it is usually used to study the reactive gas discharges with a complex chemistry set. In order to obtain the spatial distribution, a two-dimensional or three-dimensional fluid model is necessary. However, in the fluid model, the electron energy distribution function is assumed to be Maxwellian, which is invalid under special discharge conditions. For instance, strong electric field and low pressure may result in non-Maxwellian distributions, such as bi-Maxwellian distribution, two-temperature distribution, etc. Therefore, a fluid/Monte Carlo collision hybrid model is adopted to take the electron kinetics into account. Besides, a separate biased sheath model is necessary to study the influence of the sheath on the plasma properties self-consistently. The particle-in-cell/Monte Carlo collision hybrid model is a fully kinetic method based on the first-principles, which could be used to investigate the non-local and non-thermal equilibrium phenomena. In conclusion, the numerical investigation of inductively coupled plasmas has a significant importance for plasma process optimization.