The accurate prediction of transition is relevant for aerodynamic analysis and design applications. Extending the laminar flow region over airframes is a potential way to reduce the skin friction drag, which in turn reduces fuel burn and greenhouse gas emissions. This paper introduces a numerical framework that includes the modeling of transition effects for high Reynolds number flows in a high-fidelity, Reynolds–averaged Navier–Stokes (RANS) aerodynamic design framework. The CFD solver uses a discontinuous Galerkin (DG) finite element approach and includes goal-oriented adaptation. The Spalart–Allmaras (SA) turbulence model is used for the closure of the governing equations. In the flow stability analysis, the nonlocal, nonparallel effects that characterize boundary layers are accounted for by using the parabolized stability equations (PSE). Transition onset is obtained through an eN method based on the PSE computations, while a smooth intermittency function includes the transition region length. Numerical results for the NLF(1)-0416 airfoil present good agreement with experimental data, improving the computations when compared to fully-turbulent ones.