In this work, we present an adaptive multifidelity approach for the efficient uncertainty quantification of probabilistic reliability analysis, with applications in aeroelastic flutter prediction across a flight envelope. Our approach is based on a multifidelity Gaussian process model, which may integrate multiple high- and low-fidelity models, with several key extensions to existing methodology. First, we propose a new batch-sequential adaptive sampling method for parallel computing environments. The method balances estimated information gain against computational cost for each model while accounting for heterogeneous batch sizes resulting from the parallel execution of lower-fidelity models on multi-core computing systems, and selects the optimal model and sampling points. Second, we introduce a multifidelity kernel with relaxation to improve hyperparameter tuning consistency and reduce model selection bias. This addresses challenges that arise in the presence of sparse high-fidelity samples and many evaluations of poorly correlated lower-fidelity models. Third, we extend an entropy-based acquisition function to the selection of model fidelity and heterogeneously sized sampling points. We apply the framework to probabilistic flutter boundary prediction in both subsonic and transonic flow regimes. We incorporate a hierarchy of aeroelastic models that differ in both physical modeling and discretization resolution: an Euler-based finite-volume CFD solver with various time step sizes, a low-fidelity model based on the doublet-lattice method, and a low-fidelity model based on Theodorsen's unsteady airfoil theory. We demonstrate that the multifidelity approach reduces the computational cost required to achieve a given accuracy in both subsonic and transonic regimes.