ionized gas in the chemically homogeneous interior of stars are governed by the Navier-Stokes (NS) equations: $$\begin{aligned}
\label{ns1}
\partial\rho/\partial t &=& -\nabla \cdot \rho \vec{v},\\
\label{ns2}
\partial\rho \vec{v}/\partial t &=& -\nabla \cdot \rho \vec{v}\vec{v}
-\nabla p+\nabla\cdot\vec{\Sigma}+\rho\vec{g},\\
\label{ns3}
\partial E/\partial t &=&-\nabla\cdot[(E+p)\vec{v}-\vec{v}\cdot\vec{\Sigma}
+\vec{F}_{\rm d}]+\rho \vec{v}\cdot\vec{g},\end{aligned}$$ where $\rho$ is the density, $\vec{v}$ is the velocity, $p$ is the pressure, $\vec{g}$ is the gravitational acceleration, $E=e+\rho
v^2/2$ is the total energy where $e$ defined in Eq.(\[eos\]) is the internal energy including ionization energy and radiation energy. $$\vec{\Sigma}=2\mu\vec{\sigma}+\lambda(\nabla\cdot\vec{v})\vec{I}$$ is the viscous stress tensor, where $\vec{\sigma}$ is the strain rate tensor, $\vec{I}$ is the identity tensor and $\mu$ represents the sub-grid scale (SGS) eddy viscosity: $$\label{mu} \mu=\rho(c_{\rm
\mu}\Delta)^2(2\vec{\sigma}:\vec{\sigma})^{1/2}.$$ $c_{\rm \mu}$ is an adjustable constant parameter, $\Delta$ is a length scale of the order of the local grid size, and the colon represents the contraction of tensor. The bulk viscosity $\lambda$ is taken to be $-(2/3)\mu$. $$\label{fd} \vec{F}_{\rm d}=\vec{F}_{\rm rad}+\vec{F}_{\rm sgs}$$ is the diffusive flux. $\vec{F}_{\rm rad}$ is the radiative flux which can be accurately computed as diffusion in the deep interior of the star, i.e., $$\vec{F}_{\rm rad}=-\frac{4acT^3}{3\kappa\rho}\nabla T.$$ $a$ is the radiation density constant, $c$ is the speed of light, and $\kappa$ is the Rosseland 