We present an analytical method for the two-qubit quantum Rabi model. While still operating in the frame of the generalized rotating-wave approximation (GRWA), our method further embraces the idea of introducing variational parameters. The optimal value of the variational parameter is determined by minimizing the energy function of the ground state. Comparing with numerical exact results, we show that our method evidently improves the accuracy of the conventional GRWA in calculating fundamental physical quantities, such as energy spectra, mean photon number, and dynamics. Interestingly, the accuracy of our method allows us to reproduce the asymptotic behavior of mean photon number in large frequency ratio for the ground state and investigate the quasi-periodical structure of the time evolution, which are incapable of being predicted by the GRWA. The applicable parameter ranges cover the ultrastrong coupling regime, which will be helpful to recent experiments.