The nano electromechanical system (NEMS) with functionally graded (FG) nanotubes as the core component is usually in the multi field environment leading to multi-source excitation resonance. As a typical and basic excitation form, the nonlinear principal resonance of FG nanotubes under the supercritical pulsatile flow and forced excitation is studied. Geometric nonlinearity originates from the tensile hardening of immovable boundary modeled by von Kármán’s nonlinear geometric relationship. Size-dependencies for both nanofluids and nanosolids are characterized by slip-flow model and nonlocal strain gradient model coupled with surface effects, respectively. Under the framework of Zhang–Fu’s beam theory, the governing equations for coupled resonance based on the initial post-buckling equilibrium state are established. A combination of perturbation-incremental harmonic balance method (IHBM) is established and employed in bifurcation analysis of such a strongly nonlinear coupled resonance system for the first time The softening-type nonlinear behavior and multiple-stable solutions are found and two dominant factors controlling the bifurcation topology are revealed. It is suggested that the size effects of both nanoflow and nanosolid are dual and may change the bifurcation topology through two mechanisms. One is the change of two excitations from strong interaction to weak interaction, and the other is the change of multi-solution threshold of the two-excitation amplitude ratio. The results can guide the nonlinear design and adjustment of the NEMS to obtain available working frequency band maintaining stable motion state.