Short-Circuit Analytical Model for Modular Multilevel Converters Considering DC Cable Capacitance

Developing analytical short-circuit models for Modular Multilevel Converters (MMC) is not straightforward due to their switching and blocking characteristics. Short-circuit models for MMCs have been developed previously in the literature. However, there is a lack of understanding regarding the dynamics in the short-circuit model when the DC cable capacitance is taken into account. Therefore, this work proposes an analytical pole-to-pole short-circuit model for half-bridge MMCs that considers the cable capacitance and terminal capacitors and accounts their contribution to fault dynamics. An approximated analytical model has been derived separating the system solutions in different natural frequencies. The proposed model provides an excellent approximation for a vast range of realistic system parameters. The analytical model reproduced the behaviour of the variables in the time domain and provided a clear basis for interpreting the dynamics of the voltages and currents involved.


I. INTRODUCTION
Power transmission based on Voltage Source Converters and High-Voltage Direct Current (VSC-HVDC) have drawn much attention from the scientific community and industry over the last two decades. This is mainly due to the VSCs controllability and easiness to form DC grids [1].
Among the VSCs, the Modular Multilevel Converter (MMC) has been widely studied due to its advantages over two-and three-level converters. These benefits include quality of the generated waveform, full or partial blocking capability and modularity [2]. However, in order to design reliable DC grids, the short-circuit phenomenon must be better understood, which is crucial to design appropriate protection systems [3]- [6]. Moreover, the system's components can be sized using parametric EMT simulations, as in [7]- [9].
The associate editor coordinating the review of this manuscript and approving it for publication was Chandan Kumar .
However, though precise and essential, EMT simulations are time-consuming and do not reveal how each system parameter participates in the fault behaviour as analytical models do.
Developing short-circuit models for MMCs is not a straightforward task. Due to its switching characteristic, control-dependency and blocking capability, the short-circuit model must be separated in different fault periods where the equations defining the fault dynamics can be simplified [10], [11]. Moreover, the model complexity increases when other DC capacitances are considered, such as the stray capacitance of DC cables.
Concerning short-circuit models, [12] presents a poleto-pole fault calculation method for DC grids. However, in this study the DC conductors were simplified to an RL model, considering only the converter capacitance. The same simplification was adopted by [11], which analysed the influence of limiting reactors, AC grid strength and fault location on the short-circuit currents. In [13], the influence of wind power plant control methods in the DC short-circuit behaviour was analysed. The study was performed considering detailed models, but focusing on the stage after the converter blocking. The influence of the control on DC grid fault currents was analysed by [14], using the RL simplification for the DC conductors. Other studies such as [10], [15]- [17] also increased the understanding of the DC fault phenomena, but they only considered the converter capacitance in the equations.
Although these previous studies advanced in specific areas of short-circuit models, they neglected the presence of additional capacitances in the system. However, multiterminal DC systems are also intended to interconnect renewable generation, such as offshore wind farms, and this connection is generally made via cables. The submarine and underground cables capacitance, different from overhead lines, cannot be neglected [18]. Moreover, thought not essential, other DC capacitors may exist in the system, such as DC filter and DC link capacitors [19]- [23]. In this context, there is a lack of research that considers the presence of additional capacitors in the short-circuit model. Therefore, the present work proposes a short-circuit model for half-bridge MMCs (HB-MMC) that considers the presence of multiple capacitors in the DC side and accounts their contribution to pole-to-pole faults.
In order to address the aforementioned problem, precise models using ordinary differential equations (ODEs) or frequency-based models can be used [24]- [26]. However, depending on the system complexity, the models can have complex analytical solution or lose interpretability, failing to provide a clear basis for the understanding of the phenomena being analysed. Thus, to avoid both numerical and extensive analytical solutions, this research provides approximated analytical models which can interpret the short-circuit dynamics with remarkable accuracy. The proposed approximated analytical models preserve the frequency content and waveform of the variables and provide a clear basis for the understanding of the first stage of the fault. We prove that separating the solution in different time constants represents accurately the dynamic of the fault and allows the complete model to be solved, which otherwise would have extensive analytical solution.
A precise yet simple model is key to size DC circuit breakers, size fault current limiters [27] and design reliable protection algorithms, specially those based on transient and frequency components [28]- [30].
The paper is organized as follows. Section II introduces the HB-MMC operating principle and conventional fault analysis. Section III describes the proposed analytical model considering the model with cable capacitance. Section IV presents a comparative analysis between the full numerical solution and the proposed model. The proposed model is extended including a DC terminal capacitor in Section V and its assessment is presented in Section VI. Finally, the discussions and conclusions of this research are presented in Section VII and Section VIII, respectively.

II. CONVENTIONAL MMC FAULT ANALYSIS
The basic structure of a three-phase MMC is shown in Fig. 1. The three-phase MMC comprises six arms, and each arm consists of N series-connected submodules (SMs) and one arm inductor. The SMs can be configured as half-bridge or full-bridge. This work analyses the half-bridge configuration. The MMC fault analysis can be divided into three different stages [10]: the capacitive discharge stage; the AC transient infeed stage and the AC steady-state infeed stage. The capacitive discharge stage is the period between the occurrence of the fault and the instant when the converter blocks. In this stage, the energy stored in the SM capacitors is discharged in the short-circuit until the arm currents reach the limit value of the converter's internal protection. In this stage, if conventional control methods are used, the total number of capacitors in the legs is kept the same [31], as the control has no substantial time to change references [32]. However, if special control strategies as in [33], [34] are used, their influence on the number of inserted capacitors is not neglectable. In this work, the general case is considered, with constant number of inserted capacitors across converter legs. The second stage, AC transient infeed stage, is characterized by the sum of the currents flowing through the freewheeling diodes and the currents resulting from the release of the stored energy in the arm inductors and DC inductors. After the inductors' discharge, the fault is fed solely by the AC grid in the AC steady-state infeed stage. This work analyses the first stage.
Regarding the capacitive discharge stage, Fig. 2 shows the conventional equivalent circuit of the MMC during a poleto-pole fault, without considering the DC cable capacitance, where L a and R a are the arm inductance and resistance, respectively, L dc and R l are the inductance and resistance of the fault path, respectively, R f is the fault resistance and C leg is the equivalent leg capacitance in this fault stage. L dc can be considered a sum of the DC reactor (L fcl ) and the line/cable series inductance. C leg = C SM /N [35], where C SM is the capacitance of the SM, as seen in Fig. 1. For the sake VOLUME 8, 2020 FIGURE 2. MMC circuit at the capacitive discharge stage [15].
of simplicity, the DC fault path resistance can be joined to R dc = R l + R f /2. Although pole-to-pole faults might also involve earth when the cable's sheath is grounded, the models used in this paper do not consider the ground as the system is assumed to be in a balanced steady-state condition before the fault.
Applying Kirchhoff's voltage law in the equivalent circuit of Fig. 2 and assuming Moreover, Considering L eq = 2L a +6L dc , R eq = 2R a +6R dc , the system of differential first order equations is written as Let V dc and I dc be the converter voltage and total current values at t = 0. The solution for (3) is where is the circuit damped natural frequency, ω c0 = 1/ L eq C leg is the undamped natural frequency and τ = 2L eq /R eq is the circuit decaying time constant. If I dc = 0, the expression for i f (t) reduces to the one presented in [15]. It is worth highlighting that the sum of cos ω c t and sin ω c t/(ω c τ ) in (4) and (5) does not represent two oscillations, but a single oscillation at the frequency ω c with an equivalent angle and amplitude. This can be verified using appropriate trigonometric transformation.

III. MMC FAULT ANALYSIS CONSIDERING THE CABLE CAPACITANCE
Although the circuit presented in Fig. 2 is a good representation of the sole MMC during the first stage of a fault, it considers only the converter capacitors, neglecting the DC cable capacitance (Fig. 3) and other capacitive elements, such as DC pole and DC filter capacitors [19]- [23]. The DC cable capacitance can be modelled using an approximated lumped capacitor (C p ), as shown in Fig. 4, where v p and i p represent the capacitor voltage and current, respectively.
The cable T-model was used, as this model requires one capacitor less than the Pi-model. The transformation from the Pi-model to the T-model produces a negative conductance, which is numerically correct but physically meaningless [36], therefore this conductance was neglected. The addition of the cable capacitance splits the previous cable series inductance and resistance into L l1 + L l2 and R l1 + R l2 , respectively.
It is convenient to define equivalent parameters for the circuit in Fig. 4 to reduce the size of the system's Ordinary Differential Equations (ODEs): Using the equivalent parameters, the system of ODEs becomes: As (7) is linear, the system has an unique solution. However, obtaining an exact closed-form solution involves calculating the roots a forth-order polynomial, which produces lengthy expressions. To circumvent this, (7) can be transformed to the frequency domain using the Laplace transform and then simplified to obtain an approximated analytical solution. For the Laplace transform, the following initial conditions are used: i c (0) = I dc /3, v c (0) = V dc , i p (0) = 0 and v p (0) = V dc . These initial conditions represent the DC side circuit in steady state (C leg and C p fully charged and constant DC voltage equal V dc ) for a given power flow determined by I dc (see nomenclature in Fig. 4).
The solution in the frequency domain is presented in (8) and (9), in a compressed form called low-entropy format. The low-entropy format is useful for giving meaningful information about how the parameters influence the dynamics of the variable in a compact format [37]. Thus, the variables are expressed in terms of frequencies and time constants, rather than expressing them in terms of basic RLC parameters.
where s is the Laplace variable and: and ω p0 = 2L eq 3C p L dc1 L dc2 (11) where ω p0 is a new undamped natural frequency resulting from the addition of C p to the circuit, representing the cable capacitive discharge. The voltages v c (t) and v p (t) are obtained using (7b) and (7d) or its equivalents in the frequency domain. The addition of a new capacitor introduces a second oscillation defined by the frequency ω p0 and splits the previous circuit time constant τ in the time constants: The equivalent inductances and resistance are now: Inverting (8) and (9) back to the time-domain requires the roots of P(s) to be known, so that P(s) can be factored and the partial fraction expansion can be applied. As P(s) is a polynomial of forth-degree, the same issue with lengthy expressions appears when obtaining an exact solution. In order to obtain an analytical solution, the common denominator P(s) can be approximated seeing that C p C leg , which means ω c ω p . Now, ω c is still obtained by (6) but using the new L eq and R eq defined in (12). Using the low-Q approximation, which is a technique to find approximate analytical solutions for the roots of n th -order polynomials (see [38]), P(s) can be rewritten in the form: (13) to (10) produces the system of equations: Using the assumption that the natural frequency related to C leg is much smaller than the natural frequency related to C p (ω c ω p ), the approximation for the time constants and frequencies are: and This approximation means that the total response of the circuit was separated into two different oscillations, one high-frequency oscillation related to the cable, of frequency ω p , and one low-frequency oscillation related to the MMC, of frequency ω c . The approximation also means that the decaying time constant of the oscillation governed by C p (τ p ) is influenced by C leg which possibly means that the voltage in C leg sustains the voltage in C p , reducing the speed of voltage drop. The damped frequency ω p is calculated as ω p = ω 2 p0 − 1/τ 2 p . With the approximation for P(s), the voltages and currents in (8) and (9) in the frequency domain can now be inverted back to the time domain, but the amount of terms produced is enormous. Therefore, to present expressions that can show good interpretations, the smallest terms, those involving multiplication by small numbers, were neglected. The solutions are presented in (15), as shown at the bottom of the next page, where: Analysing (15) it can be observed that the total response of each variable is a sum of a low-frequency term, governed by ω c , and a high-frequency term, governed by ω p . The contribution of C p to the low-frequency part of i c (t) and i p (t) is minor as C p plays a major role in the dynamics of the high-frequency part. If C p is neglected, then L eq = 2L s and the low-frequency part of i c (t) will become the same as in the conventional model, given by (5).
As ω c ω p , the addition of a new small capacitor has a minor impact on v c (t), which means that the energy stored in the cable capacitance is not sufficient to refrain the voltage drop in the MMC. The contribution of i p (t) to the low-frequency content of the short-circuit current is minor as τ rp τ and C p C leg . The dynamics of v p (t) is largely influenced by the ratio between L dc1 /L eq and R dc1 /R eq , expressed in the parameter k τ and in the terms that multiply V dc and I dc . This is a consequence of the system itself, as C p discharges partially on the MMC and partially on the fault.
The superposition of low-and high-frequency terms will only be precise when ω c and ω p were clearly separable in the frequency spectra. If both frequencies are close in value, the dynamics will not be separable, and the equivalent frequency will be defined by a contribution of both capacitors. We observed that if the frequencies are separated at least by one octave, the proposed solution presents well-approximated results.

IV. ASSESSMENT OF THE MODEL WITH CABLE CAPACITANCE
In order to assess the validity and precision of the proposed model, the approximated and numerical solution of (7) were compared in time domain. Fig. 5 shows both responses using the parameters summarised in Table 1. We should highlight that although the converters might be blocked before 5 ms, we have shown this time span only to prove that the model is precise for a long time interval.
As can be observed in Fig. 5, there was no noticeable difference between the full numerical solution of the system of ODEs and the approximated analytical solution. It also confirms what was revealed by the equations, that the voltages and currents are formed by a superposition on a low and a high-frequency oscillation with defined frequencies and decay constants.
To assess the validity of the approximated model not only for the parameters presented but for a wide range of feasible values, a parameter-sensibility test was conducted.
The system basic parameters L a , R a , L l1 + L l2 , R l1 + R l2 , C p , C leg and R f were varied in a range from 1/10 to 10 times the value shown in Table 1. Then, the average difference between the full and the approximated model was calculated. To avoid divisions by zero, the following relative error was calculated: difference between the analytical model and the numerical model divided by I dc for i c (t) and i p (t) and divided by V dc for v c (t) and v p (t). The average relative errors are shown in Fig. 6. In Fig. 6, it can be observed that for the tested range of parameters the proposed analytical model approximated the full solution with low error. From 1/10 to 10, there was a difference of 100 times the parameters values, which highlights how the analysed range was wide. In most cases, the errors were less than 1%. The cases of greater error were the ones of higher ratio C p /C leg . This happened because when the difference between C p and C leg reduces, ω c and ω p get closer and the assumption ω p ω c starts to lose validity. Based on the results shown, the proposed model was considered valid for the whole range of parameters, which confirms its precision while maintaining the desired interpretability.

V. MMC FAULT ANALYSIS CONSIDERING THE CABLE CAPACITANCE AND TERMINAL CAPACITOR
When additional DC capacitors are connected to the system, such as the DC pole or DC filter capacitor [19]- [23], they interact with the system during a pole-to-pole fault adding another oscillatory component to the fault current. Their influence and interaction with the MMC and with the DC cable capacitance can be modelled by a lumped capacitor C b connected in parallel between the MMC and the cable, as shown in Fig. 7.
The inclusion of C b to the circuit separates the previous equivalents L dc1 and R dc1 from L a and R a and adds two new variables to the problem. To distinguish these parameters from the parameters used with the cable model only, the new equivalent parameters are renamed as: Hence, the system of ODEs becomes: To find the exact solution of (16), the roots of a sixth-order polynomial must be obtained algebraically, which is not possible for polynomials of order greater than four VOLUME 8, 2020 (see Abel-Ruffini Theorem). However, the same approach used in (7) can be applied for (16): solve the system in the frequency domain, approximate the dynamics considering the frequency of oscillation imposed by C b (ω b ) separated from ω c and ω p and then approximate the solution using the low-Q approximation. For the Laplace transform, the following initial conditions are considered for the new variables: As ω p is imposed by the cable/line capacitance and this capacitance is a function of cable length and properties, ω p can be either bigger or smaller than ω b . Here, ω p was considered bigger than ω b .
The solution in the frequency domain for I b (s) is presented in (17). I c (s) and V c (s) were considered unaffected by C p because C b actuates as a filter for the high-frequency oscillation associated with C p , as can be verified from the simulations depicted in Fig. 9. I p (s) and V p (s) were not presented because the focus of the analysis was on the components at the converter terminals. The components associated with ω c were not presented because they were considered unaffected, hence they are the same as in (15a) and (15b). V b (s) is obtained through (16d).
where (18) and the new values for R eq and L eq are: and the undamped natural frequencies become As C b splits the system series parameters the new time constants are: (18), differently from P(s) in (10), has two frequency components influencing the terms s 4 and s 3 . This suggests that the dynamics of the components governed by ω p and ω b will have higher cross coupling than the previous case with only ω c and ω p . Summed to this, here we can no longer assume ω p ω b as we assumed ω p ω c , because as ω p is related to the cable/line capacitance, this capacitance can be even bigger than C b depending on the cable/line length and configuration. Therefore, the low-Q approximation, in this case, is expected to be less precise than the case with just ω c and ω p , and the quality of the approximation will be directly dependent on the difference between ω p and ω b .
The approximated Q(s) is: Using the approximated Q(s) and inverting (17) back to time yields the solution in (22a), as shown at the bottom of the next page. Then, using (16d) yields (22b), as shown at the bottom of the next page.
As can be observed in (22), the components in i b (t) and v b (t) related to ω b and to ω p are mutual and cross linked, showing that their dynamics were not totally separated. This happened because we considered ω p > ω b instead of ω p ω b . It should be noted that even in systems where there is no terminal capacitor and the DC side capacitance is given only by the cable capacitance, the model with two capacitors can also be used by splitting the total cable capacitance in two equal parts between C b and C p . This increases the model precision as well as the model complexity.

VI. ASSESSMENT OF THE MODEL WITH CABLE CAPACITANCE AND TERMINAL CAPACITOR
To assess the precision of the short-circuit model that considers both the DC terminal capacitor and cable capacitance, the approximated solutions for v c (t), i c (t), v b (t) and i b (t) were compared to simulations of a two-terminal MMC-HVDC system modeled in PSCAD/EMTDC. The system was adapted from [39]. The DC link was modelled using the frequency-dependent model. The cable parameters were based on [40]. Pole-to-pole faults were applied along the system indicated in Fig. 8 with fault resistance of 5 . The converters were intentionally not blocked during the faults to show that the proposed model is precise even for long time intervals.  The cable equivalent series resistance and inductance were calculated and used separately for each frequency (ω c , ω b and ω p ) because these parameters depend greatly on frequency. The travelling wave delay was shifted when plotting the graphs to have the t = 0 when the wave reaches the converter. The initial current and voltage were I dc = −1.24 kA and V dc = 624 kV. The curves referred to as conventional in Fig. 9 were calculated using the model presented in Section II.
In Fig. 9, it can be observed that even not considering the travelling wave effect, the approximated analytical solution had a similar response to the simulated waveforms. The differences can be justified by some non-considered effects. Besides the travelling wave, the approximated model does not consider the system beyond the fault and other possible connected converters. As other converters also feed the fault, the voltage at the fault point will be partially raised by their fault current contribution. Moreover, even considering two capacitors in the DC side, the lumped model is still an approximation of the cable distributed capacitance. A third non-considered effect is the modulation action. It can be observed in the v c (t) curve that the converter modulation and submodule balancing algorithm still operate during the fault, inserting capacitors into the circuit that were not previously being discharged. This instantly raises v c (t) above the predicted value and appears in the v b (t) curve as a difference in the low-frequency component. For i c (t), the analytical solution was close to the simulated value, and for i b (t) the frequency was correct, but a deviation was seen in the exponential decay. The component associated to ω p vanished around 3 ms in i b (t). In the presented fault case, the damped natural frequencies were: ω c = 211 rad/s, ω b = 3973 rad/s and ω p = 8326 rad/s. The equivalent lumped capacitances were C p = 4.85 µF and C b = 2.0 µF. In Fig. 10 it can be confirmed that the proposed model has greater precision compared to the conventional model described by (4)-(6).

A. MAXIMUM FAULT CURRENT ESTIMATION
The proposed model is useful to better understand the fault phenomena and to discriminate the participation of the system's parameters on the fault currents. In addition, it can be used to calculate fault currents in the time domain using (15) and (22) and also to estimate the maximum fault currents given a set of parameters. What is commonly intended from a maximum fault current equation is a simple approximated expression to be used as a starting point in a design phase, before time-consuming EMT simulations. In addition, the maximum fault current should be a conservative estimate, that provides a sufficiently precise value but with a safety margin. Therefore, as the fault current through a DC circuit breaker is calculated as i f (t) = 3i c (t) + i b (t), and the low-frequency part of i c (t) has monotonic increasing behaviour during the first milliseconds, a good estimation can be performed considering the maximum values of the high-frequency parts of i c (t) and i b (t). To achieve this, the exponential terms and sine terms of the high-frequency parts are considered equal to 1 and the cosine terms are neglected as they are much smaller than the sine terms. Even considering these approximations, as the low-frequency part of i c (t) has monotonically increasing behaviour, the maximum current value will be defined after the time instant is specified. Thus, the maximum fault current i f ,max is still a function of time: where, As an example, Fig. 11 compares the maximum fault current calculated by the proposed model and simulated in PSCAD, for a 5 fault 30 km from MMC 2. As can be observed in Fig. 11, the proposed model was very precise in calculating the maximum fault current with a sufficient safety margin for all time instants, confirming that it can be used to estimate maximum fault currents.

VII. DISCUSSION
Compared to previous research, in this study it could be observed that when other capacitances in the HVDC system are considered, such as the DC cable capacitance, the complexity of the problem increases and some assumptions are required to solve it.
The assumption used either for deriving the expressions considering the DC cable capacitance or both the DC terminal capacitor and DC cable capacitance was that the natural frequencies are clearly separable in the frequency spectrum. This assumption can be better visualised if the capacitors are substituted by voltage sources at the given natural frequencies and with the given damping. If this substitution is made, the circuit can be solved by superposition of the excitations of each separable source. When the frequencies are not well separated, this approximated superposition loses its validity because the charge in each capacitor influences the charge of the others. Thus, it can be shown that each additional capacitor in the circuit that introduces a new node in the circuit will add a new natural oscillation, and this new natural oscillation will have components depending on V dc and I dc , except for cases where they are neglected.
Short-circuit models are useful to define circuit breaker requirements. As DC faults occur in the millisecond scale and C leg has higher stored energy than C b and C p , the components related to ω c are seen as low-frequency trends. Therefore, according to the proposed model, the maximum current a circuit breaker must withstand can be approximated as the value of the components related to ω c at the given instant plus the maximum of the components related to ω b and ω p . This is another benefit of the analytical model over numerical models and simulations, that the influence of each parameter on one variable, e.g. the maximum fault current, can be better analysed separately.
For this reason, the proposed approximated analytical model is useful, providing meaningful insight into the components involved in the DC fault when the capacitive contributions are considered, with good accuracy as verified in the comparison with the numerical simulation.

VIII. CONCLUSION
Analytical models of MMCs were previously developed considering only the converter contribution to the fault. However, there was a lack of research that considered the presence of the DC cable capacitance in the analytical short-circuit model. Therefore, this paper proposed an analytical poleto-pole short-circuit model for HB-MMCs that considers the cable capacitance and terminal capacitor and accounts their influence on the fault current waveform.
As the complete analytical solution involves lengthy expressions in time and the numerical solution does not comprehend the variables governing the dynamics of the problem, two approximated analytical models were derived. The models were derived separating the problems in independent natural oscillations, which have shorter analytical solutions, and then summing the separate oscillations obtaining an approximated analytical model. The proposed model including cable capacitance provided excellent approximation (less than 1% of error) for a vast range of realistic parameters, compared with the ODEs numerical solution. In addition to its precision, the model also provided a clear basis for interpreting the dynamics of the voltages and currents involved, showing that the variables are formed by a superposition of oscillations with defined frequencies and decay constants.
The expanded model that considered both DC cable capacitance and DC terminal capacitor also proved to be useful for increasing the understanding of the influence of additional capacitive elements on the fault. The comparison between the proposed model to simulations confirmed the precision of the solution.
The presented derivations extend those found in the literature, providing a more in-depth understanding of the transient short-circuit period of MMC-HVDC systems with respect to the converter's voltages and currents.