The invention relates to a power grid harmonic power flow phasor matrix calculation method supporting rapid development and model expansion. The method comprises the steps of S1, setting harmonic times, statistical model initial values of various harmonic sources and an iterative error limit; S2, generating a fundamental wave and harmonic wave admittance matrix according to the input network parameters; S3, calculatng fundamental wave power flow; S4, detecting a node with a harmonic source, generating an index matrix, predicting harmonic data by using a harmonic source statistical model according to a fundamental wave power flow result, and performing harmonic power flow calculation; S5, updating the fundamental wave load power by utilizing the harmonic network loss power; and S6, judging whether the front and back node injection power errors are smaller than an iterative error limit, if so, calculating convergence, outputting a harmonic power flow calculation result, and otherwise, returning to the step S3. According to the power grid harmonic power flow phasor matrix calculation method based on support of rapid development and model expansion and program architecture based on thinking of the phasor matrix, rapid development and efficient calculation of large-scale power grid harmonic power flow simulation are realized.