  In recent years, coupled with traditional turbulence models, the second-order
gas-kinetic scheme (GKS) has been used in the turbulent flow simulations. At
the same time, high-order GKS has been developed, such as the two-stage
fourth-order scheme (S2O4) GKS, and used for laminar flow calculations. In this
paper, targeting on the high-Reynolds number engineering turbulent flows, an
implicit high-order GKS with Lower-Upper Symmetric Gauss-Seidel (LU-SGS)
technique is developed under the S2O4 framework. Based on Vreman-type LES model
and $k - \omega$ SST model, a turbulent relaxation time is obtained and used
for an enlarged particle collision time in the implicit high-order GKS for the
high-Reynolds number turbulent flows. Numerical experiments include
incompressible decaying homogeneous isotropic turbulence, incompressible
high-Reynolds number flat plate turbulent flow, incompressible turbulence
around NACA0012 airfoil, transonic turbulence around RAE2822 airfoil, and
transonic high-Reynolds number ARA M100 wing-body turbulence. Comparisons among
the numerical solutions from current implicit high-order GKS, the explicit
high-order GKS, the implicit second-order GKS, and experimental measurements
have been conducted. Through these examples, it is concluded that the
high-order GKS has high accuracy in space and time, especially for smooth
flows, obtaining more accurate turbulent flow fields on coarse grids compared
with second-order GKS. In addition, significant acceleration on computational
efficiency, as well as super robustness in simulating complex flows are
confirmed for current implicit high-order GKS. This study also indicates that
turbulence modeling plays a dominant role in the capturing physical solution,
such as in the transonic three-dimensional complex RANS simulation, in
comparison with numerical discretization errors.
