In this paper, we design an effective and robust model to solve the 3D single-fluid solar wind plasma magnetohydrodynamics (MHD) problem of low plasma β. This MHD model is formulated on a six-component composite grid system free of polar singularities. The computational domain ranges from the solar surface to the super-Alfvénic region. As common to all MHD codes, this code must handle the physical positivity-preserving property, time-step enlargement, and magnetic field divergence-free maintenance. To maintain physical positivity, we employ a positivity-preserving Harten–Lax–van Leer Riemann solver and take a self-adjusting and positivity-preserving method for variable reconstruction. To loosen the time-step limitation, we resort to the implicit lower–upper symmetric Gauss–Seidel method and keep the sparse Jacobian matrix diagonally dominant to improve the convergence rate. To deal with the constant theme of a magnetic field that is divergence-free, we adopt a globally solenoidality-preserving approach. After establishing the solar wind model, we use its explicit and implicit versions to numerically investigate the steady-state solar wind in Carrington rotations (CRs) 2172 and 2210. Both simulations achieve almost the same results for the two CRs and are basically consistent with solar coronal observations and mapped in situ interplanetary measurements. Furthermore, we use the implicit method to conduct an ad hoc simulation by multiplying the initial magnetic field of CR 2172 with a factor of 6. The simulation shows that the model can robustly and efficiently deal with the problem of a plasma β as low as about 5 × 10−7. Therefore, the established implicit solar wind MHD model is very promising for simulating complex and strong magnetic environments.
