You might model PLL in the way for a digial PD:
PD+ChargePump Gain = I/2PI
It's a state average modeling method, which involves use the pulse width to esimate the average signal amplitude as practised by power electronics. Assume high frequency components will be filled out---as in this case by the load cap of the charge pump.
Quantitively, you can do this as long as your PLL BW is below 1/5 of carrier freq according to F. Garner.
For high BW above 1/5, you need go into discrete time modeling to do simulation. The best paper of Z domain model I found is one published over 10 years ago by one IEEE life fellow from UC-Irvine. And for recently results, check out the papers by Michael Perrott from MIT.