Clock meshes have found increasingly wide applications in today's high-performance IC designs. The inherent routing redundancies associated with clock meshes lead to improved clock skews and reliability. However, the high complexity of clock meshes in modern chip designs has made its verification very challenging. A typical clock distribution network may consist of millions of coupled/interconnected linear elements and hundreds of nonlinear clock drivers attached at different locations on the mesh. Such a large network is often too complex for feasible SPICE-like simulation. In this paper, we present a new simulation methodology which decomposes a clock mesh into linear and nonlinear parts. By exploiting the special matrix property of the linear subsystem resulting from modified nodal analysis (MNA) formulation, the linear subsystem is represented as a matrix-level macro model,which greatly simplifies the overall simulation task. These macromodels can be efficiently computed using Cholesky factorization and significantly speedup the nonlinear Newton-Raphson iterations used inthe transient simulation for the complete clock mesh. Furthermore, a dynamic time step rounding technique is proposed to limit the number of passive macromodels needed in the entire transient simulation which further improves the efficiency of the proposed approach.