Turbulent Rayleigh–Bénard convection (RBC), triggered in a fluid layer uniformly heated from below and cooled from above, is considered a paradigm for buoyancy-driven turbulence, ranging from stellar convection through atmospheric and oceanic turbulence to cooling blankets in nuclear engineering and the storage of renewable energy in liquid metal batteries. The flow dynamics and heat transfer in regimes with extremely large or small Prandtl numbers (Pr), corresponding to geophysical and astrophysical flows and to liquid metals, respectively, have not been fully verified, and much remains to be learned. While exploring such regimes poses a challenge for both direct numerical simulations and laboratory experiments, it opens the door to large-eddy simulation (LES) due to its ability to balance accuracy and computational cost. LES resolves turbulent structures larger than a specified filter scale while employing subgrid-scale (SGS) models to represent smaller, unresolved scales \cite{pope20}. This work utilizes the minimum-dissipation LES model to model both the SGS momentum flux and heat flux. Combined with symmetry-preserving discretization, the current simulations are stable (without explicitly introducing numerical dissipation) on highly stretched computational grids. This non-error-based stability enables improved understanding of the Prandtl-number dependence across a wide range of Rayleigh numbers (Ra), particularly in the ultimate regime. The proposed LES scheme has been applied to a differentially heated cavity at Rayleigh numbers up to $10^{15}$ \cite{sun26}. The results show that the LES accurately captures flow features, heat transfer, and transition points across a wide range of Rayleigh numbers, even on extremely coarse meshes. The predicted power law, $Nu_{LES}=0.157Ra^{0.281}$, is very close to the DNS result, $Nu_{DNS} =0.182Ra^{0.275}$, and the relative error remains below $1.6\%$. The eddy viscosity also $\overline{\nu_e}$ exhibits a power-law dependence on the Rayleigh number. At $\mathrm{Ra} = 2 \times 10^9$, the results show weak dependence on grid resolution, demonstrating the robustness of the minimum-dissipation model and symmetry-preserving discretization for simulating thermally driven turbulent flows on relatively coarse meshes. This LES framework is currently being applied to Rayleigh–Bénard convection at extreme Prandtl numbers and high Rayleigh numbers in the ultimate regime. The results are expected to be presented at the