Convection in superimposed layers of fluids heated from below is commonly observed in many industrial and natural occurrences, such as crystal growth, co-extrusion processes and atmospheric flow. The stability analysis of this system reveals a complex dynamic behavior, with multiplicity of stationary states and temporal oscillations. In the present study, a linear stability analysis is applied, allowing the determination of the convection onset as a function of imposed boundary conditions, different geometric configurations and perturbations with specific wavenumbers. However, non-linear effects are neglected in this method of analysis and it can induce to deviations in comparison with experimental observations. To investigate the differences from a fully non-linear model, a direct numerical simulation approach using CFD techniques is used, reveling excellent agreement near the convection onset and an increase in the deviation as the Rayleigh number increases beyond the critical value.