We improve our previous variational method based nonrelativistic quark model by introducing a complete set of three-dimensional harmonic oscillator bases as the spatial part of the total wave function. To assess the validity of our approach, we compare the binding energy thus calculated with the exact value for the hydrogen model. After fitting to the masses of the ground state hadrons, we apply our new method to analyzing the doubly heavy tetraquark states $q{q}^{\prime }\overline{Q}\overline{Q^{\prime }}$ and compare the results for the binding energies to results in other works. We also calculate the ground state masses of $T_{sc}\left(ud\overline{s}\overline{c}\right)$ and $T_{sb}\left(ud\overline{s}\overline{b}\right)$ with $(I,S)=(0,1),(0,2)$. We find that $T_{bb}\left(ud\overline{b}\overline{b}\right)$ and $us\overline{b}\overline{b}$, both with $(I,S)=(0,1)$, are stable against the two lowest threshold meson states with binding energies $-145$ and $-42\mathrm{MeV}$, respectively. We further find that $T_{cb}\left(ud\overline{c}\overline{b}\right)$ is near the lowest threshold. The spatial sizes for the tetraquarks are also discussed.