We present results of a lattice QCD calculation of $B\to D^*$ and $B_s\to D_s^*$ axial vector matrix elements with both states at rest. These zero recoil matrix elements provide the normalization necessary to infer a value for the CKM matrix element $|V_{cb}|$ from experimental measurements of $\bar{B}^0\to D^{*+}\ell^-\bar{\nu}$ and $\bar{B}^0_s\to D_s^{*+}\ell^-\bar{\nu}$ decay. Results are derived from correlation functions computed with highly improved staggered quarks (HISQ) for light, strange, and charm quark propagators, and nonrelativistic QCD for the bottom quark propagator. The calculation of correlation functions employs MILC Collaboration ensembles over a range of three lattice spacings. These gauge field configurations include sea quark effects of charm, strange, and equal-mass up and down quarks. We use ensembles with physically light up and down quarks, as well as heavier values. Our main results are $\mathcal{F}^{B\to D^*}(1)= 0.895\pm 0.010_{\mathrm{stat}}\pm{{0.024}_{\mathrm{sys}}}$ and $\mathcal{F}^{B_s\to D_s^*}(1)= 0.883\pm 0.010_{\mathrm{stat}}\pm{0.028_{\mathrm{sys}}}$. We discuss the consequences for $|V_{cb}|$ in light of recent investigations into the extrapolation of experimental data to zero recoil.