  (Abridged) We describe the results of three-dimensional (3D) numerical
simulations designed to study turbulent convection in the stellar interiors,
and compare them to stellar mixing-length theory (MLT). Simulations in 2D are
significantly different from 3D, both in terms of flow morphology and velocity
amplitude. Convective mixing regions are better predicted using a [dynamic
boundary condition] based on the bulk Richardson number than by purely local,
static criteria like Schwarzschild or Ledoux. MLT gives a good description of
the velocity scale and temperature gradient for a mixing length of $\sim 1.1
H_p$ for shell convection, however there are other important effects that it
does not capture near boundaries. Convective "overshooting" is best described
as an elastic response by the convective boundary, rather than ballistic
penetration of the stable layers by turbulent eddies. We find that the rate at
which material entrainment proceeds at the boundaries is consistent with
analogous laboratory experiments as well as simulation and observation of
terrestrial atmospheric mixing. In particular, the normalized entrainment rate
E=$u_E/\sigma_H$, is well described by a power law dependence on the bulk
Richardson number $Ri_B = \Delta b L/\sigma_H^2$ for the conditions studied,
$20\lesssim Ri_B \lesssim 420$. We find $E = A Ri_B^{-n}$, with best fit
values, $\log A = 0.027 \pm 0.38$, and $n = 1.05 \pm 0.21$. We discuss the
applicability of these results to stellar evolution calculations.
