Modeling of Conductors Catenary in Power Lines: Effects on the Surge Propagation Due to Direct and Indirect Lightning

In this article, the effects of the variable height in overhead power lines are investigated with reference to the surges produced by direct and indirect lightning events. Simulations are carried out both by an implicit finite-difference time-domain scheme and by the commercial tool EMTP-RV. Differences in the waveforms computed by the two approaches are discussed and clarified. Additionally, the validity of the commonly adopted assumption of weakly varying vertical component of the incident electric field is investigated.


I. INTRODUCTION
L IGHTNING strokes represent a severe threat to industrial facilities, infrastructures, or overhead lines [1], [2], [3], [4], being natural hazards that may trigger possible subsequent technological disasters [5].Both direct and indirect lightning strikes pose severe electrical stress to installed power equipment, which should have the proper insulation strength withstanding the lightning impulse voltage to avoid insulation failure and breakdown.
Rachidi [6] proved the equivalence of the different approaches proposed in the literature for the analysis of the interaction between a lightning electromagnetic pulse (LEMP) and a transmission line (TL), namely, those by the authors of [7], [8], and [9].A comprehensive dissertation on the assumptions and derivation of the different field-to-line coupling models is reported in [10], [11], and [12].Other models, e.g., Rusck model and Chowduri model, were revised in [13] and [14].
In the present work, the results of overvoltages produced (direct lightning) and induced (indirect lightning) on multiconductor transmission lines (MTLs) are discussed considering the conductors' sags.In order to show particularly the consequences of nonuniformity of the TL [15], the line is initially considered lossless above a perfectly conducting plane (PEC).Different approximate approaches were proposed to account for the finite conductivity of the soil in the computation of the radial component of the incident electric field (e.g., the wellknown Cooray-Rubinstein correction term [16], [17]); as to the propagation of the vertical component of the incident electric field, the contribution of the soil finite electric conductivity and permittivity is commonly neglected [17].The inclusion of losses in a time-domain algorithm through accurate transient parameters accounting for losses in the conductors and in the ground valid in a wide frequency range (i.e., for wide time intervals) has been presented in [18].The finite-difference time-domain (FDTD) analysis of a lossy nonuniform line has been presented in [19].
The current research is directed towards more efficient implementation of FDTD codes for field-to-line coupling studies [20], improved interface with EMT simulators [21], and assessment of the effects of LEMP [22].The modified telegrapher's equations are usually solved through explicit FDTD algorithms because of their reduced runtime and simplicity; although other approaches, e.g., the one presented in [23], may offer a major accuracy, on the other hand, the computational process is time consuming.In this article, the implicit Crank-Nicolson FDTD scheme with second-order accuracy is implemented: its longer computational time compared with explicit schemes is compensated by the intrinsic stability [24].First, an investigation is carried out with reference to the propagation of overvoltages caused by both direct and indirect lightning events, approximating the line conductors (which have variable height) as a uniform MTL with an average equivalent height.Then, the approximation is removed by considering the height variable along the line.A comparison between the results obtained by the EMTP-RV tool and those obtained by means of the FDTD code is presented.Both implementations are based on the discretization of each span into minor sections.
Differences in the induced overvoltages computed by the FDTD code and the software EMTP-RV-LIOV were observed when modeling the conductors' catenary; the underlying reasons and the accuracy of the results by the software were investigated.An explanation is proposed by showing that a modification of the FDTD code, which alters its standard implementation, allows to reproduce the software results by enforcing a discontinuity in the scattered voltage (which is an intrinsically continuous quantity instead [25]).
The rest of this article is organized as follows.After presenting the modified telegrapher's equations in Section II, the assumptions adopted in the computation of the lightning electromagnetic field are summarized in Section III; Section IV deals with the solution of the telegrapher's equations by the implicit FDTD scheme in combination with different modeling approaches of the conductor's catenary.Numerical results are presented in Sections V and VI.Finally, Section VII is dedicated to conclusions.Appendices A and B include additional numerical results for the interested reader.

II. AGRAWAL ET AL. FIELD-TO-LINE COUPLING MODEL FOR A NONUNIFORM LINE
Agrawal et al. field-to-line coupling model is implemented in the FDTD code to simulate the effects of indirect lightning.It is also adopted by, for example, the LIOV toolbox provided by EMTP-RV [26].
The direction of propagation of voltage and current surges is assumed along the x-axis; the modified telegrapher's equations accounting for the LEMP illumination of a lossless TL read [8] (for the sake of conciseness, reported for a single conductor) The lossless line case has been considered also for comparison purposes since the LIOV code does not include the internal and ground per unit length (p.u.l.) transient impedances [27].In (1), L (x) and C (x) denote the p.u.l.external inductance and capacitance of the single-conductor TL, respectively; E x i is the component in the x-direction of the incident electric field.By means of Agrawal et al. model, propagation is assessed in terms of the scattered voltage v s .A transformation is required to obtain the total voltage at points located at height h(x) or to simulate the interaction between the TL conductors and other devices installed along the line (e.g., grounding of the shield wires (SWs) at the towers, surge arresters, etc.).In case of PEC ground, it reads (2) Herein, subscripts t, s, and i refer to the total, scattered, and incident quantities (fields and voltages), respectively.In (2), v t is the total voltage-to-ground of the conductor and E z i is the vertical component of the incident electric field.
In propagation studies along MTLs, when the points of suspension at the span terminations present the same heights, the geometrical average height h ave = h max − 2 3 S of the conductor is adopted according to the common practice; S is the average sag (which depends on the conductor p.u.l.weight, installation configuration, external temperature, and environmental conditions), and h max (h min ) is the maximum (minimum) height of the conductor within the span, i.e., in correspondence of a tower (at midspan).Nevertheless, the catenary profile of each conductor within the spans (symmetrical with respect to their midpoint) will be considered in the following discussion, approximating its height by a second-order polynomial (i.e., a parabola) [28]: h(x) is the variable height of the conductor, denotes the span length, and x refers to the longitudinal position along the span, being x = 0 and x = the left and right endpoints of the span, respectively.In case, the line includes also inclined spans above the horizontal soil, expression (3) should be replaced by the specific equation in order to account for the different conductor's arrangement and height at the span endpoints [29].

III. FIELD COMPUTATION
The relevant approximation adopted here consists in neglecting the actual tortuosity of the lightning channel in the computation of the fields produced by the return-stroke current [30], [31]; the lightning channel is modeled by a vertical cylinder above a PEC plane [32].The fields are computed by the integral forms given by Rachidi et al. [33].
Many analytical expressions are available in the literature to describe the propagation of the return-stroke current along the lightning channel (comprehensively reviewed in [32] and [34]).The TL model will be considered for comparison purposes since this model is adopted in EMTP-RV-LIOV.In this case, the returnstroke current i rs (z, t), as a function of time t and distance z from the ground, is given by where u(t) is the unit step Heaviside function, ν f is the upwardpropagating front speed (also called return-stroke speed), and ν 0 is the current-wave propagation speed.The value ν 0 = 1.5 • 10 8 m/s is chosen for simulations, being in the range of common values adopted in the literature, along with the hypothesis ν f = ν 0 [32], [35].Since engineering TL models for the propagation of return-stroke currents do not account for reflections at the upper endpoint of the lightning channel, any evaluation of (4) at time t > H/ν 0 may lead to unreliable results (values of the channel height H are expected to be within some kilometers, depending on the specific lightning event and geographical location).

A. Line Model in the Implicit FDTD Code
The Crank-Nicolson implicit scheme, equivalent to the single average alternating direction implicit (ADI) scheme [36], is Fig. 1.Piecewise linear approximation of the line catenary implemented in the implicit FDTD code.The case of a lossless single-conductor TL is displayed.adopted to solve the telegrapher's equations for the lossless line case [24].When the expressions of the FDTD method are modified to account for field-to-line coupling according to the Agrawal et al. model, the following equations are obtained for an MTL with N c conductors: and In ( 5) and ( 6), whose single-conductor equivalent circuit representation is shown in Fig. 1 for three consecutive Δx, [1] denotes the identity matrix of order N c ; L k+ 1 2 and C k are the matrices of p.u.l.external inductances and capacitances of the line; V n s k is the vector of dimension N c × 1 of conductors' scattered voltages at node k and time nΔt; and I n k+ 1   2   is the vector of dimension N c × 1 of conductors' currents flowing from node k toward node (k + 1) at time nΔt.The ith element of the 2 )Δx, y is equal to the average conductors coordinate in the y-direction, z = h i k,k+1 is equal to Expression (7) corresponds to the average height h i k,k+1 of the ith conductor, of variable height (h i ), given by an expression of the type (3) within each span, between nodes located at x = kΔx and x = (k + 1)Δx.
The quantities h i k,k+1 , with i = 1, . . ., N c , are used to derive the conductor matrix of p.u.l.external inductances L k+ 1 2 .Instead, the matrix of p.u.l.capacitances is derived by inversion of the matrix of potential coefficients, which are computed referring to the actual conductors' heights at the defined nodes.The space mesh for the evaluation of the unknown voltages and currents at time t = (n + 1)Δt in ( 5) and ( 6) corresponds to that commonly adopted by the explicit leap-frog scheme since the line nodes for the evaluation of voltages and currents are separated by half a step Δx/2; yet, they are calculated at the same time instants nΔt, disclosing the implicit nature of the method [24].
Zheng et al. proved the unconditional stability of the ADI algorithm [37]; hence, the choice of Δx and Δt for the implicit scheme may be conveniently related to considerations on the numerical burden, on the scheme dispersion error, on the geometry of the configuration under study, and on the known or predicted time constants involved in the specific transient study.
Finally, total voltages are obtained at the nodes of interest by means of the finite-difference counterpart of (2) It should be noted that an approximation is introduced in (8) since the vector of incident voltages at node k is computed by the product of the N c × 1 vector h k (whose elements are the conductors' heights at distance kΔx from the MTL left endpoint) with the quantity E n+1 z k , i.e., the vertical component of the incident electric field at x = kΔx, y (defined previously), and at z = 0.
Piecewise linear (or constant) line modeling approaches neglect any proximity effect that actually exists between the catenary segments.Rigourously, if the segment inclination was considered, the scattered electric field would not be confined in the y-z plane.Nevertheless, similar simplified approaches with two-dimensional approximations are often to be preferred in terms of implementation simplicity and shorter simulation time.

B. Line Model in EMTP-LIOV
The line is simulated in EMTP-RV-LIOV by means of multiple connected LEMP-illuminated lines.The macrostep Δx m = 10 m is chosen to discretize the line catenary.Hence, a number /Δx m of EMTP-LIOV devices connected in series is used to model each span of the line.The simulation time step Δt = 3.33 ns has been set.Since the software adopts a trapezoidal integration method with Courant factor f C = 1, the corresponding space step is equal to Δx = 1 m.The constant Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.height of each conductor within a macrostep is set as the average value of its actual height (accounting for the catenary profile) at x = kΔx m and x = (k + 1)Δx m .
This implementation guarantees the continuity of the total voltages-to-ground of the conductors at the point of connection between LEMP-illuminated line sections with conductors located at different heights (i.e., macrosteps).

C. Modified Line Model in the Implicit FDTD Code
Each span is discretized into segments of length Δx.However, macrosteps Δx m = k m Δx are introduced to model a coarser discretization mesh of the catenary profile, gathering k m small space steps with the same height.The height of the segment within each Δx m is computed as the average height of the catenary nodes limiting each macrostep.The discretization criterion is represented in Fig. 3.
This second approach is implemented as an additional modification to the FDTD solving scheme to compare with results by the LIOV code (as detailed in Section IV-B), which automatically sets the space discretization step to Δx = c 0 Δt (corresponding to f C = 1), c 0 being the velocity of light.The user may either represent the line through a very high number of LIOV lines with different heights (minding possible truncation errors in case the EMTP-computed Δx = c 0 Δt differs from the chosen length of the LEMP-illuminated line segment) or adopt a reduced number of line sections with length equal to Δx m ; the macrostep Δx m is successively discretized by the software into a finer mesh of segments Δx with constant conductors height within each Δx m .
Following the suggestions in [38] and the software documentation, a good approximation of the LIOV results is achieved by means of this modified FDTD approach (as discussed in Section V).The discontinuity introduced by the piecewise constant approximation of the variable conductor height in EMTP-RV is to be enforced into the modified FDTD code through two ideal voltage sources accounting for the different incident voltages at the junction point between adjacent macrosteps at different heights.
The modification to the equivalent circuit at these specific voltage nodes is schematically represented in Fig. 4, where subscripts L and R are employed to simplify the notation and denote the (constant-valued) inductances and capacitances within the

V. PRELIMINARY SIMULATIONS
In this section, the simple configuration sketched in Fig. 5 is considered.The used version of the LIOV toolbox allows for the simulation of lossless MTLs with three or four conductors; hence, three, nonfed, equal conductors (with same h max and S) were simulated for reference, displaying the overvoltages of the central conductor (which lies on the axis y = 0).The line total length is 2 km, consisting of five equal spans with length = 400 m, with symmetric catenary profile with respect to midspan and matched terminations.The average sag is S = 10 m and h max = 27 m.The FDTD code may also be used to simulate more extended TLs; increased computational time is expected due to the larger number of discretization steps (increasing the order of the solving system) and points along the line requiring the computation of the components of the incident electric field at each time iteration.For direct lightning studies (when neglecting the effects of the external field), larger Δx may be adopted to alleviate the total computational burden; nevertheless, for simulations of indirect lightning overvoltages, this choice  should be preceded by the evaluation of the loss of the overall accuracy for larger Δx.In fact, besides a poorer representation of the catenary, the validity of the approximation of a constant E x i along the chosen Δx should be considered (for the presented configuration, the choice of a step equal to 10Δx could also be suitable [40]).EMPT-RV guarantees shorter simulation times, requiring rather cumbersome modeling effort when accounting for the catenary, since a considerable number of line devices with different geometrical characteristics should be individually inserted by the user.
First, overvoltages caused by a lightning current striking the central conductor, as computed by the implicit FDTD code and by EMTP-RV software with Δt = 3.33 ns, are compared when the catenary profile is simulated.In particular, Δx m = 10 m is adopted for simulations in EMTP-RV; as to the implicit FDTD code, the adopted discretization step is Δx = 1 m with Courant factor f C = 1.
Successively, results are compared with reference to an indirect lightning current striking the ground (PEC surface) at x s -y s (the coordinate reference system is shown in Fig. 5).
The same waveform is implemented as the channel base return-stroke current for both direct and indirect lightning simulations, expressed as the sum of two Heidler's functions with parameters in Table I (for a median subsequent stroke current in [39]).

A. Direct Lightning
In the simulated scenario, the lightning current strikes the central conductor at x = 1.2 km.The overvoltages computed by the FDTD algorithm and by EMTP-RV, which are displayed in Fig. 6 at different observation points along the TL and account for the conductor catenary, are practically superimposed, enabling (and validating) results from the developed code.The curves computed by the modified FDTD method are not included since the additional voltage sources in Fig. 4 are commonly neglected for the simulations of direct lightning, and results by the modified and the non-modified methods practically coincide.
The overvoltages computed with the constant height h ave = 20.3m are also shown.In this case, the predicted voltages peak values accounting for the influence of the variable height of the conductor and those computed for the corresponding uniform TL at h ave differ by less than 1%.The effect of the catenary is mainly evident in the low-frequency oscillatory behavior superimposed to the main voltage waveform, which is to be observed more easily at the tail.Its period ≈ 2.67 μs is related to 2 /c 0 , i.e., to the periodicity of the line geometry.Due to the line nonuniformity and discretization with constant Δx, the low-frequency superimposed oscillations may be more or less pronounced depending on the observation point along the line.
As expected, the results support the validity of the average height approximation, which may be adopted to study the propagation of conducted disturbances, e.g., the overvoltages produced by the direct lightning of overhead lines with typical geometric characteristics and values of the conductors sag.

B. Indirect Lightning
Two different channel base locations are simulated as to indirect lightning strokes, namely, x s = 1.0 km, y s = 0.05 km and x s = 1.0 km, y s = 0.10 km.Initially, results from the proposed FDTD code are validated by comparison of the induced overvoltages with results from the LIOV code when accounting for h ave .With reference to the observation points identified in Fig. 7, results are found to agree more than satisfactorily and highlight the natural reduction of the induced surges peak values as the return-stroke channel base is moved farther from the line.
Fig. 8 displays the influence of the conductor catenary on overvoltages computed by the implicit FDTD code (described in Section IV-A), EMTP-RV (in Section IV-B) and the modified FDTD code (in Section IV-C) with x s =1.0, 0.8 km and y s = 0.05 km.Fig. 8(a) refers to the simulated configuration when the lightning channel base is located opposite a span midpoint (i.e., where the conductor has the minimum height).Fig. 8(b) refers to simulations conducted with the lightning channel base located opposite a TL tower (i.e., where the conductor has its maximum height).Overvoltages at x = 0.6 km and x = 1.0 km Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.are superimposed due to the symmetric position with respect to x s ; differences are noted at t 4.7 μs since reflections from the left endpoint reach the closest observation point at x = 0.6 km first. 1 Waveforms obtained with h ave are included for reference.
The results highlight the impact of four relevant aspects.

1) Relative position between the channel base and the line:
By comparing overvoltages in Fig. 8(a) and (b), it may be deduced that harsher overvoltages are expected when the lightning channel is opposite a TL tower.Referring to results computed by the FDTD code, Fig. 9 displays incident and scattered voltages at the observation point closest to the LEMP source and at half-span distance, when the channel base is located at x s = 0.8, 1.0 km and y s = 0.05 km.The scattered voltages obtained with x s = 1.0 km and x s = 0.8 km are less sensitive to the variable conductor's height: in fact, their amplitude is comparable at the observation point opposite the lightning channel (x = 1.0 km and x = 0.8 km, respectively), and the deviation is even reduced at half-span distance (x = 0.8 km and x = 1.0 km, respectively).Hence, the conductor height h(x) contributes significantly to differentiating the total voltages at the two observation points through the incident voltage term, displayed in Fig. 9(a), which holds a predominant role in the determination of the overvoltages pattern (further details in Appendix A). 2) Average height approximation: Depending on the modeling approach, results obtained for the uniform TL with h ave present a deviation, which may be larger than 20%.At the observation point closest to the LEMP source, the average height approach over-(under-)estimates the induced overvoltages when x s -y s are opposite a point where the MTL conductors display the minimum (maximum) height.Similar deviations are not found for simulations of direct lightning that only account for the catenary by means of the height-dependent matrices of the p.u.l.capacitances and inductances of the line (not performing any additional integration of the incident field along the conductor's height).3) Chosen line model: The proposed implementation allows to represent the contribution of the incident voltages more accurately.In fact, it does not resort to any artificial discontinuity in the incident (and, as a consequence, in the scattered) voltages, overcoming this intrinsic limitation of staircase modeling; it considers the line integral of the vertical component of the incident electric field from a point at null potential (z = 0 for the PEC plane case) to the variable conductor height z = h(x).This feature of the FDTD scheme leads to nonnegligible differences against results by the modified FDTD scheme (and EMTP-RV-LIOV), which are enhanced at observation points closer to the LEMP source [e.g., at x = 0.8 km in Fig. 8(b)] and in the proximity of the spans' endpoints [e.g., at x = 0.8 km Fig. 8(a)].With reference to (8), this is justified, respectively, by the larger values of E z i and by the more pronounced height discontinuity between macrosteps (due to the reduced accuracy of the staircase model in representing the actual catenary, where |dh(x)/dx| is larger, i.e., at the span endpoints).These elements accentuate any difference in the amplitude of the additional voltage sources, as introduced in Section IV-C.4) Courant factor: Fig. 10 displays the total voltages already presented in Fig. 8(b), with lightning channel base located at x s = 0.8 km, y s = 0.05 km but with f C = 2, 5, 10, consistent with the time constants typical of lightning studies.Good accuracy of results is retained compared to curves obtained with f C = 1, speeding up the simulation time by a factor 2, 5, and 10, respectively.The differences observed at early times in the inset of Fig. 10 do not substantially affect the voltage front and peak value, relevant to the estimate of the electrical stress caused by lightning overvoltages close to the source (where the most severe stress typically occurs).At observation points located several hundreds of meters from the source, where numerical dispersion is expected to affect more significantly the computed curves, the compromise between accuracy and simulation time suggests to adopt f C ≤ 5 for the analyzed lightning events and configurations [24].

A. Test Line
The same channel base lightning current is used for the computation of induced voltages along the conductors of the lossless MTL, as depicted in Fig. 11, with geometrical characteristics in Table II.The line section has length equal to 2 km.Phase conductors are terminated in their characteristic resistance and the SW is grounded by a small-valued resistance (equal to 1 μΩ, approximating a quasi-ideal grounding condition).Two locations of the lightning channel base are simulated, i.e., x s = 0.8 km and x s = 1.0 km, with y s = 0.05 km.The  II.

TABLE II CHARACTERISTICS OF THE SIMULATED LINE
distance of the channel base was chosen by application of the electrogeometric model [41].
Results from EMTP-RV-LIOV and the implicit FDTD code are compared.In the former case, the line is modeled by the cascade of LEMP-illuminated lines with different heights and length equal to 10 m; the selected time step is Δt = 3.33 ns.As to the original FDTD code, the line is discretized into segments Δx = 5 m, with f C = 0.2.The coarser space mesh adopted in the FDTD code, not affecting the results [40], allows to speed up the computations by reducing the dimension of the solving system.

B. Effect of the Conductors' Catenary
In Fig. 12, the overvoltages are displayed for two different channel base locations, accounting for the conductors catenary (computed by the proposed FDTD method and by the LIOV code) and for the corresponding uniform MTL. 2  Comparing accounting for the catenary or for the average heights, overvoltages are found to be affected predominantly (and nonnegligibly) by the simulated variable conductors' height at observation points closer to the LEMP source.However, regardless the relative position of the lightning channel with respect to the line, peak values of the overvoltages computed along the phase conductor are affected by wider deviations (27%-35%) with respect to those found for the SW (12%-21%).This might be justified by the larger sag characterizing the installation of overhead phase conductors, hence, by the rougher 2 Conversely to results in Figs.7 and 8, overvoltages for the uniform MTL by the FDTD code differ slightly from the ones by EMTP-RV, probably due to the specific derivation of E x along the conductors.In the FDTD scheme, it is evaluated at ȳ at the actual height of the ith conductor; in EMTP-RV, it is calculated at the highest point of the MTL and retrieved along the lower conductors by a linear approximation [27].
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.approximation introduced by the average height approach as to indirect lightning studies.Some different assumptions in the implementation of Agrawal et al. model by the LIOV code (in the first place, the numerical aspects discussed in Section IV-B and IV-C) lead to differences in the range 5%-17% on the predicted peak value of the overvoltages at the observation point nearest to the channel base.Enhanced differences in the tail of the waveforms computed at x = 0.8 km are also noticed for results in Fig. 12(b).This is likely due to a combination of the relative position of the lightning channel base and the wider height gap between subsequent LEMP-illuminated line devices in correspondence of the span endpoints (the adopted piecewise constant approximation for the conductors catenary is less accurate in proximity of the poles, causing harsher discontinuities).

C. Effect of a Different Integration Method for the Incident Voltages
A three-point trapezoidal quadrature rule [42] was included in the code to compute the incident voltages and riser terms at the line terminations (as defined in [43]) in order to assess the validity of (8).The accuracy of the single-point integration method, based on the assumption of weakly varying E z i , may not be satisfactory for larger heights [i.e., longer integration paths in (2)].Hence, the following is implemented: where denotes the Hadamard product; for instance, the total voltage-to-ground V n+1 A,t k of conductor A, with height h A k at node k, may be computed as follows: (10) In ( 9), the quantity E n+1 z k h k denotes a vector with dimensions N c × 1; the ith element corresponds to the vertical component of the incident electric field evaluated at distance kΔx from the line left termination, at z equal to the height of the ith conductor.
Fig. 13(a) displays overvoltages along the SW and conductor B when the MTL is modeled as uniform, with conductors located at their average (constant) height.The simulation is carried out for x s = 1.0 km, y s = 0.05 km for the indirect lightning case; herein, direct lightning is not addressed since only the predominant effect of the injected current is typically considered in this type of simulation.
The waveforms obtained by the single-point (8) or by the three-point (9) quadrature rule converge at observation points farther from the lightning channel; differences of 17% and 8% are observed at x = 1.0 km in the voltage peak values along the SW and conductor B, respectively.Specifically, the improved integration results in the reduced amplitude of the computed overvoltages at points opposite the lightning channel.Already at early times, the single-point integration method overestimates the incident voltage, adopting the value of the vertical component of the electric field attained at z = 0 along the whole integration path and neglecting the time delay necessary for the external field to propagate at 0 < z ≤ h i .
The comparison between results obtained by the two integration methods when the conductors' catenary is simulated and the channel base is located at x s = 1.0 km-y s = 0.05 km or x s = 0.8 km-y s = 0.05 km is displayed in Fig. 13(b) and (c), respectively.
It is confirmed that less accurate integration methods for the calculation of the incident voltages lead to an overestimation of the total voltages-to-ground at points closer to the source.The higher is the local height of the conductors, the more relevant is the (precautionary) overestimation of the conductor voltages by the single-point integration method with respect to the threepoint quadrature; indeed, the deviation of 23.0% of SW peak voltage at x = x s = 0.8 km [in Fig. 13(c)] reduces to 14.3% for x = x s = 1.0 km [in Fig. 13(b)] and is generally wider than that for conductor B, located below the SW.
As expected, the reflected waves produced by the SW grounding at the line endpoints increase the voltage difference between the SW and conductor B, i.e., the stress across the line insulators.

VII. CONCLUSION
This article contributes to assessing the influence of the conductors' catenary as regards overvoltages caused by direct lightning currents striking a line conductor or due to coupling with incident electromagnetic fields.
The relevance of the implemented line model, i.e., uniform line, piecewise constant, or piecewise-linear discretization of the conductors' heights, has been assessed.It has been shown that neglecting the MTL nonuniformity may lead to considerable underestimation of the peak value of the overvoltages induced by external electromagnetic fields.This conclusion holds for lightning channel base locations which are opposite the line towers, where the conductor's height is maximum, due to the contribution of the incident voltages.
Despite the increased computational time, the original FDTD code guarantees several elements of flexibility.Only a few additional code lines allow to arbitrary chose finer discretization meshes in order to model the height variations accurately.Depending on the conductors' height, different numerical integration schemes may be implemented as to the evaluation of the contribution of the incident voltages.
Even though flexibility and accuracy come at the cost of an increased computational burden, measures to speed up the computation can be applied (e.g., reducing the dimension of the solving system by adopting larger discretization step Δx while ensuring results' deviations within a tolerance; adopting larger Courant factors, depending on the time constants involved in the specific study; simulating critical line sections and stroke locations, etc.) Further research may involve analogous studies on power lines with different features (height, sag, etc.); the simulation of the conductors' catenary in dealing with LEMP-induced effects may affect insulation coordination considerations, especially along the line section closest to the lightning channel base.

APPENDIX A
With reference to the configuration analyzed in Section V, overvoltages computed at additional observation points along the line, for the nonuniform and the constant height case, are displayed in Fig. 14(a) with x s = 0.8 km and y s = 0.05 km.Fig. 14(b) and (c) refer to incident and scattered voltages, respectively.
The pattern of the predicted waveforms along the line changes whether accounting or not for the conductors' catenary.For the uniform line, the voltage peak values decrease moving farther from the lightning channel; this trend is not to be found in the computed overvoltages when the catenary is simulated.In fact, in the latter case, the waveform at x 7 = 0.40 km displays larger amplitude with respect to those computed at x 5 = 0.60 km and x 6 = 0.50 km (closer to the LEMP source).This is mainly due to the contribution of incident voltages, being the scattered ones less affected by the catenary; in Fig. 14(c), when the catenary is simulated, larger amplitudes are expected for the incident voltages where the height is larger than h ave .
The percentage deviation between results for the uniform and nonuniform line also depends on x.However, the large deviation of 30% at x = x s is noticeably reduced to less than 10% at half-span distance (x ≤ 0.60 km).At larger distances from the source, a reverse in the polarity of scattered voltages is observed; the amplitude of total voltages is mainly determined by scattered voltages, and the contribution of incident voltages (predominantly affected by the local conductors' height) is limited.
Finally, in Fig. 14(c), the low-frequency oscillations of the scattered voltages, only visible when the catenary is considered, are associated to reflections from subsequent spans, hence, to propagation phenomena along the MTL conductors.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.  in [18] and [24], i.e., by means of a recursive method exploiting suitable approximating functions of the internal and ground p.u.l.impedance obtained by Prony's algorithm.The towers were modeled as lossless TLs with characteristic impedance 175 Ω and propagation velocity 0.85c 0 , with grounding resistance 15.8 Ω.The same lightning current, described by the sum of Heidler's functions with parameters in Table I, was injected at x = 1.2 km, accounting for the lightning channel equivalent resistance 400 Ω [44].The effect of catenary is confirmed to be negligible for direct lightning simulations; the attenuation due to the transient internal and the ground impedance may be already observed at x = 1.0 km, with the multiple reflections due to the periodic grounding of the SW at the towers.
The inception of corona discharge along the line further supports the negligibility of the conductors' catenary in direct lightning studies.It should be noted that accurate simulations require the inception and development characteristics of impulse corona to be dependent on voltage steepness, polarity, weather conditions, etc. [45], [46], [47], [48].In this regard, due to the unavailability of experimental data for the simulated line, the corona model by Mihailescu-Suliciu and Suliciu [49] with symmetrical behavior under different voltage polarities was selected to provide a first estimate of the effect of corona discharge.The parameters given in [ [50], Sec.5.3], comparing the calculated overvoltages against measured data, were used.The Peek formula was implemented to compute the corona inception field (accounting for a position-dependent inception field due to the variable conductors' height, with surface irregularity factor m = 0.75).The corona current in the transverse plane is simulated by means of distributed voltage-controlled current sources.Fig. 15(b) reflects the sizeable impact of corona discharge on traveling voltage waves, practically overshadowing any catenary effect, when plotted against the curves accounting for all the previous simplifying assumptions (lossless MTL, conductors at average height, absence of corona).

Fig. 2 .
Fig.2.Representation of the voltage nodes within a span of a single-conductor TL, with reference to the piecewise linear approximation of the line catenary implemented in the implicit FDTD code.

Fig. 3 .
Fig.3.Piecewise constant approximation of the line catenary implemented in EMTP-RV, and in the modified FDTD code as to reproduce results from the LIOV code; the TL is subdivided into macrosteps with length Δx m (a simple line discretization with k m =Δx m /Δx =3 is shown as an example).

Fig. 4 .
Fig. 4. Modification of the equivalent circuit for discretization steps Δx located at a junction node between macrosteps of the line with different heights, as by the piecewise constant approximation of the line catenary.

Fig. 5 .
Fig. 5. Sketch of the top and lateral view of the central conductor of the simulated three-conductor line; the simulated lightning channel base locations (x s , y s ) and the striking point of the direct stroke are displayed.

Fig. 6 .
Fig. 6.Overvoltages at different observation points along the central conductor of the line in Fig. 5 due to a direct lightning stroke at x = 1.2 km.

Fig. 7 .
Fig. 7. Overvoltages at different observation points along the line central conductor due to indirect lightning strokes with channel base located at x s = 1.0 km and y s = 0.05, 0.10 km; the simulated MTL is uniform, with conductors' height equal to h ave .

Fig. 8 .
Fig. 8. Overvoltages at different observation points along the line central conductor, induced by the indirect lightning current with channel base located at x s -y s .Results account for constant and variable conductor height, computed by the LIOV code, the proposed and the modified FDTD code.Percentage deviations of the peak values are computed with respect to results by the proposed FDTD code when accounting for the catenary.(a) Channel base location point: x s = 1.0 km, y s = 0.05 km.(b) Channel base location point: x s = 0.8 km, y s = 0.05 km.

Fig. 9 .
Fig. 9. Incident and scattered voltages at different observation points along the line central conductor, with channel base location (x s , y s ), accounting for the conductor catenary.(a) Incident voltages.(b) Scattered voltages.

Fig. 10 .
Fig. 10.Overvoltages computed by the FDTD code at different observation points along the line central conductor, induced by the indirect lightning current with channel base located at x s = 0.8 km, y s = 0.05 km.Results account for constant and variable conductor height, and different Courant factors f C .

Fig. 11 .
Fig. 11.Arrangement of the conductors at the towers with geometrical characteristics in TableII.

Fig. 12 .
Fig. 12. Overvoltages at different observation points along the line (for SW and conductor B), induced by the indirect lightning current with channel base located at x s , y s .Results account for constant and variable conductor height, computed by the LIOV code and the proposed FDTD code.Percentage deviations of the peak voltage values are computed with respect to results by the proposed FDTD code when accounting for the catenary.(a) Channel base location point: x s = 1.0 km, y s = 0.05 km.(b) Channel base location point: x s = 0.8 km, y s = 0.05 km.

Fig. 13 .
Fig. 13.Overvoltages at different observation points along the line (for SW and conductor B), induced by the indirect lightning current with channel base located at x s , y s .Results account for constant and variable conductor height, computed by the original FDTD code with a single-or a three-point integration of the vertical component of the incident electric field.Percentage deviations of the peak values are computed with respect to results obtained by the threepoint quadrature.(a) Channel base location point: x s = 1.0 km, y s = 0.05 km.Conductors at constant (average) height.(b) Channel base location point: x s = 1.0 km, y s = 0.05 km.Conductors' catenary is accounted.(c) Channel base location point: x s = 0.8 km, y s = 0.05 km.Conductors' catenary is accounted.

Fig. 14 .
Fig. 14.Total, and scattered voltages at different observation points along the central conductor of the line sketched in Fig. 5, with channel base location x s = 0.8 km, y s = 0.05 km.Results are computed by the FDTD code and account for constant and variable conductor height.(a) Total voltages-toground.(b) Incident voltages.(c) Scattered voltages.

Fig. 15 .
Fig. 15.Overvoltages at different observation points along the SW and conductor B of the line sketched in Fig. 11 due to a direct lightning stroke at x = 1.2 km.Curves in (a) present the effects of losses and catenary, while curves in (b) also account for the corona discharge along the line.