Earth return

LineCableModels.Engine.EarthReturn

Compute the earth-return impedances and potential coefficients of cables. A formula gives the impedance, the potential coefficients, or both, by the methods it defines.

Formulas

Formula{:unified}

Earth-return impedances and potential coefficients of conductors in two homogeneous half-spaces, the air, medium 1 with $y > 0$, and the earth, medium 2 with $y < 0$, for fields with the dependence $e^{-\Gamma \ell}$ along the line, where $\ell$ is the axial coordinate. Each cable enters as a solid conductor bounded by its exterior circumference, which lies wholly in one medium. The field of each conductor is that of one axial line source, and each conductor takes the circumferential mean of the field on its circumference, the zeroth azimuthal harmonic. The formula contains the solutions of the underground case [3] and of the mixed overhead and underground case [4].

For each medium $t$, the complex conductivity $\widehat\sigma_t$ [S/m], the permeability $\mu_t$ [H/m], the transverse propagation parameter $\kappa_t$ [1/m] and the vertical decay parameter $a_t$ [1/m] are

\[\widehat\sigma_t = \frac{1}{\rho_t} + s \varepsilon_0 \varepsilon_{r,t}, \qquad \mu_t = \mu_0 \mu_{r,t}, \qquad \kappa_t^2 = s \mu_t \widehat\sigma_t - \Gamma^2, \qquad a_t = \sqrt{\lambda^2 + \kappa_t^2},\]

where $s = \mathrm{j} 2 \pi f$ is the complex frequency [1/s], $\rho_t$ is the resistivity [Ω·m], $\varepsilon_{r,t}$ and $\mu_{r,t}$ are the relative permittivity and permeability [dimensionless], $\Gamma$ is the longitudinal wavenumber [1/m], and $\lambda$ is the transverse Fourier variable [1/m]. The constants $\varepsilon_0$ and $\mu_0$ are vacuum_permittivity and vacuum_permeability. The square roots take the decaying branch. The spectral denominators are

\[\Delta_\mu = \mu_2 a_1 + \mu_1 a_2, \qquad \Delta_\sigma = \widehat\sigma_2 a_1 + \widehat\sigma_1 a_2.\]

Conductor $i$ has its center at $(x_i, y_i)$, the distance $h_i = |y_i|$ to the interface, the exterior radius $r_i$ [m] and the medium $t = m(i)$. In a coefficient $ij$, conductor $i$ is the target and conductor $j$ the source, in the medium $s = m(j)$, and $\bar s$ is the other medium. The distance between the centers is $D_{ij}$, the distance from conductor $i$ to the image of conductor $j$ is $D'_{ij} = \sqrt{x_{ij}^2 + (h_i + h_j)^2}$, with $D'_{ii} = 2 h_i$, and $x_{ij} = x_i - x_j$ [m].

For conductors in the same medium, $t = s$, the means of the direct and image fields and the interface corrections are

\[\begin{aligned} \Lambda_{ii} &= K_0(\kappa_s r_i) - I_0(\kappa_s r_i) K_0(2 \kappa_s h_i), \qquad \Lambda_{ij} = I_0(\kappa_s r_i) \left[K_0(\kappa_s D_{ij}) - K_0(\kappa_s D'_{ij})\right], \\ J_{ij} &= \int_0^\infty \frac{\mu_{\bar s}}{\Delta_\mu} e^{-a_s (h_i + h_j)} \cos(\lambda x_{ij}) \, \mathrm{d}\lambda, \qquad M_{ij} = \int_0^\infty \frac{\widehat\sigma_s (\mu_s a_s + \mu_{\bar s} a_{\bar s})}{\Delta_\mu \Delta_\sigma} e^{-a_s (h_i + h_j)} \cos(\lambda x_{ij}) \, \mathrm{d}\lambda, \end{aligned}\]

where $I_\nu$ and $K_\nu$ are the modified Bessel functions of order $\nu$. The axial vector-potential kernel $\widetilde L_{ij}$ [H/m] and the scalar-potential kernel $\widetilde\Phi_{ij}$ [m/F] are

\[\widetilde L_{ij} = \frac{\mu_s}{2 \pi} \left[\Lambda_{ij} + 2 I_0(\kappa_s r_i) J_{ij}\right], \qquad \widetilde\Phi_{ij} = \frac{s}{2 \pi \widehat\sigma_s} \left[\Lambda_{ij} + 2 I_0(\kappa_s r_i) M_{ij}\right]\]

for $t = s$, and, from the transmitted fields,

\[\begin{aligned} \widetilde L_{ij} &= \frac{\mu_t \mu_s}{\pi} I_0(\kappa_t r_i) \int_0^\infty \frac{e^{-a_t h_i - a_s h_j}}{\Delta_\mu} \cos(\lambda x_{ij}) \, \mathrm{d}\lambda, \\ \widetilde\Phi_{ij} &= \frac{s}{\pi} I_0(\kappa_t r_i) \int_0^\infty \frac{\mu_1 a_1 + \mu_2 a_2}{\Delta_\mu \Delta_\sigma} e^{-a_t h_i - a_s h_j} \cos(\lambda x_{ij}) \, \mathrm{d}\lambda \end{aligned}\]

for $t \neq s$. The kernel $\widetilde K_{ij} = s \widetilde L_{ij} - (\Gamma^2 / s) \widetilde\Phi_{ij}$ [Ω/m] is the negative mean axial electric field at conductor $i$ per source amplitude at conductor $j$.

The voltage of a conductor is the circumferential mean of the integral of the vertical electric field from a reference level. A conductor in the earth is referred to infinite depth, and a conductor in the air to the interface below it. The source-normalized voltage coefficients $\widetilde P_{ij}$ [m/F] are

\[\begin{aligned} \widetilde P_{ij} &= \delta_{ts} \frac{s}{2 \pi \widehat\sigma_2} \Lambda_{ij} + \frac{s}{\pi} I_0(\kappa_2 r_i) \int_0^\infty \frac{a_1}{a_2 \Delta_\sigma} e^{-a_2 h_i - a_s h_j} \cos(\lambda x_{ij}) \, \mathrm{d}\lambda, & t &= 2, \\ \widetilde P_{ij} &= \delta_{ts} \frac{s}{2 \pi \widehat\sigma_1} \Lambda_{ij} + \frac{s}{\pi} \int_0^\infty \frac{a_2}{a_1 \Delta_\sigma} \left[I_0(\kappa_1 r_i) e^{-a_1 h_i} - J_0(\lambda r_i)\right] e^{-a_s h_j} \cos(\lambda x_{ij}) \, \mathrm{d}\lambda, & t &= 1, \end{aligned}\]

where $\delta$ is the Kronecker delta and $J_0$ is the Bessel function of order zero, which averages the interface endpoint over the horizontal extent of the circumference.

The enclosed currents follow from the complete field on the circumferences, with the reflected and mutual fields, through the dimensionless current matrix $\mathbf T$, held as L in the state of the formula, and the diagonal matrices $\mathbf W$ [dimensionless] and $\mathbf F$ [m/Ω]. Elimination of the source amplitudes gives the earth-return potential coefficients $\mathbf P_g$ [m/F] and impedances $\mathbf Z_g$ [Ω/m]:

\[\mathbf T = \mathbf W^{-1} - \mathbf F \widetilde{\mathbf K}, \qquad W_{ii} = I_0(\kappa_t r_i), \qquad F_{ii} = \frac{2 \pi \widehat\sigma_t r_i I_1(\kappa_t r_i)}{\kappa_t I_0(\kappa_t r_i)},\]

\[\mathbf P_g = \widetilde{\mathbf P} \mathbf T^{-1}, \qquad \mathbf Z_g = \left(\widetilde{\mathbf K} + \frac{\Gamma^2}{s} \widetilde{\mathbf P}\right) \mathbf T^{-1}.\]

The internal impedances of the conductors are not part of $\mathbf Z_g$: the impedance inside each cable, Zin, gives them.

The integrals run along the complex ray that ray! selects, and integrate evaluates them. The kernel $\widetilde\Phi_{ij}$ enters only for $\Gamma \neq 0$. With $\xi = \kappa_t r_i$, the ratio $I_1(\xi) / (\xi I_0(\xi))$ in $F_{ii}$ is the series $1/2 - \xi^2/16 + \xi^4/96 - 11 \xi^6/6144$ for a nominal $|\xi| < 10^{-3}$. The direct and image terms are evaluated together for a nominal $|\kappa_s D'_{ij}| < 0.5$, and their difference is $\ln(D'_{ij} / D_{ij})$ for $\kappa_s = 0$, with $D_{ii} = r_i$. The implementation scales the column of each source $j$ by $e^{|\operatorname{Re} \kappa_s r_j|}$ with exponentially scaled Bessel functions, and the scaling cancels in $\mathbf P_g$ and $\mathbf Z_g$.

Identifiers of the earth-return formulas.

struct Formula{ID, O<:LineCableModels.Commons.FormulationOptions} <: LineCableModels.Engine.AbstractEarthReturn

Represent an earth-return impedance and potential-coefficient formula. Its options form a record of type FormulationOptions{Formula{ID}}, with ID the formula identifier.

• options::LineCableModels.Commons.FormulationOptions: Options of the formula.

Formula(
    ::Val{ID};
    kw...
) -> LineCableModels.Engine.EarthReturn.Formula

Build the formula named by ID with its declared defaults and supplied options.

Arguments

• Val{ID}: formula identifier from FORMULAS.

Keywords

• kw...: options declared by defaults for Formula{ID}.

Returns

• A Formula with a typed option record.

Errors

• Throws ArgumentError for an unknown identifier or option.
• Throws ArgumentError for invalid integration controls.

defaults(
    _::Type{LineCableModels.Engine.EarthReturn.Formula{:unified}}
) -> @NamedTuple{integral::Pair{typeof(LineCableModels.Engine.integrate), Symbol}}

Return the declared options of the unified earth-return formula.

Arguments

• Formula{:unified}: earth-return formula type.

Returns

• integral: the default integration selection, integrate => :quad.

validate(
    options::LineCableModels.Commons.FormulationOptions{LineCableModels.Engine.EarthReturn.Formula{:unified}}
) -> LineCableModels.Commons.FormulationOptions{LineCableModels.Engine.EarthReturn.Formula{:unified}}

Check the integration controls of the unified earth-return formula.

Arguments

• options: formula options with an integration section.

Returns

• options, unchanged.

Errors

• Throws ArgumentError for a negative or nonfinite tolerance, two zero tolerances, or a nonpositive evaluation limit.

validate(
    problem::LineParametersProblem,
    formula::LineCableModels.Engine.EarthReturn.Formula{:unified}
) -> LineParametersProblem

Check that each outer cable circumference lies strictly within one medium. The circumference-current formula requires an outer radius smaller than the distance from the cable center to the interface.

Arguments

• problem: cable geometry and earth model.
• formula: unified earth-return formula that consumes the problem.

Returns

• problem, unchanged.

Errors

• Throws DomainError for a cable whose outer circumference touches or crosses the interface.

initialize_buffers(
    formula::LineCableModels.Engine.EarthReturn.Formula{:unified},
    _::LineParametersProblem,
    input::LineCableModels.Engine.CoaxialInput{T, L} where L<:(AbstractArray{<:Material{T}, 1}),
    buffers
) -> Union{NamedTuple, DataStructures.SortedMultiDict}

Allocate cable matrices and state vectors under unified, then add the integral arrays for the selected quadrature method.

Arguments

• formula: unified earth-return formula with integration controls.
• LineParametersProblem: problem whose computation needs the arrays.
• input: numerical cable geometry.
• buffers: computation arrays already allocated.

Returns

• buffers extended with the unified and integral blocks. prepare! fills the unified arrays.

prepare!(
    formula::LineCableModels.Engine.EarthReturn.Formula{:unified},
    buffers,
    input::LineCableModels.Engine.CoaxialInput,
    media,
    s,
    k
) -> Any

Compute the earth-return impedances and potential coefficients of the cables at the frequency of index k by the formula Formula{:unified}, whose equations are in the section "Formulas" of EarthReturn. The longitudinal wavenumber is buffers.gamma[k]. The results and the shared state of the two quantities stay in buffers.unified, which impedance! and potential! read.

Arguments

• formula: unified earth-return formula with integration controls.
• buffers: computation arrays with gamma, integral, and unified.
• input: cable centers and outer radii.
• media: air and earth records with resistivity rho in [Ω·m] and relative eps_r, mu_r [dimensionless].
• s: complex frequency, j 2π f, [1/s].
• k: frequency index.

Returns

• buffers, with unified.Ze [Ω/m], unified.Pe [m/F], and their shared state filled.

Errors

• Throws DomainError for a nonfinite weighted earth kernel, a nonfinite integral, or a spectral ray whose geometric weight grows.

impedance!(
    formula::LineCableModels.Engine.EarthReturn.Formula{:unified},
    Zg::AbstractMatrix,
    buffers,
    input,
    media,
    s,
    k
) -> Any

Copy the earth-return impedances computed by prepare! into the destination.

Arguments

• formula: unified earth-return formula.
• Zg: destination for cable impedances in [Ω/m].
• buffers: computation arrays with the prepared unified.Ze matrix.
• input: numerical cable geometry and connections.
• media: air and earth properties at the frequency.
• s: complex frequency, j 2π f, [1/s].
• k: frequency index.

Returns

• Zg, with earth-return impedances [Ω/m].

potential!(
    formula::LineCableModels.Engine.EarthReturn.Formula{:unified},
    Pg::AbstractMatrix,
    buffers,
    input,
    media,
    s,
    k
) -> Any

Copy the earth-return potential coefficients computed by prepare! into the destination.

Arguments

• formula: unified earth-return formula.
• Pg: destination cable potential-coefficient matrix [m/F].
• buffers: computation arrays with the prepared unified.Pe matrix.
• input: numerical cable geometry and connections.
• media: air and earth properties at the frequency.
• s: complex frequency, j 2π f, [1/s].
• k: frequency index.

Returns

• Pg, with earth-return potential coefficients [m/F].

struct SpectralIntegral{F<:Function}

Represent the integral of a complete integrand over $[0, \infty)$. The function includes the weights and Jacobians of its formula.

• integrand::Function: Complete integrand, a function of the nonnegative real coordinate.

integrate(
    integral::LineCableModels.Engine.SpectralIntegral,
    ::Val{:quad},
    section::NamedTuple,
    buffers;
    points
)

Evaluate an integral by adaptive Gauss–Kronrod quadrature with QuadGK. Transform its nonnegative coordinate $\lambda$ into $t$:

\[\int_0^\infty g(\lambda)\,\mathrm{d}\lambda = \int_0^1 g\!\left(\frac{t}{1-t}\right)\frac{\mathrm{d}t}{(1-t)^2}.\]

Here $g$ is the complete integrand. Map each breakpoint $\lambda_b$ to $t_b = \lambda_b / (1 + \lambda_b)$ in $[0, 1)$ and seed one segment between adjacent points. The error norm is the magnitude of the nominal value.

Arguments

• integral: complete integrand and its integral.
• Val(:quad): adaptive quadrature method tag.
• section: controls rtol, atol, and maxevals from defaults.
• buffers: computation arrays with the integral block.

Keywords

• points: finite, nonnegative breakpoints in the coordinate of the integrand.

Returns

• The integral value, in the unit of the integrand times its coordinate.
• The estimated absolute error, in the unit of the integral value.

Notes

An estimated error above max(atol, rtol * abs(nominal(value))), or a nonfinite estimate, produces a warning. The method still returns the value.

Errors

• Throws DomainError for a negative or nonfinite breakpoint, or a nonfinite integral value.

ray!(
    buffers,
    wavenumbers::Tuple,
    poles::Tuple,
    scales::Tuple,
    height,
    reach
) -> Any

Choose a complex spectral ray and write quadrature breakpoints in its dimensionless coordinate $t$. The spectral coordinate is

\[\lambda = R e^{\mathrm{j}\theta}t, \qquad R = \max\!\left(\max_m |k_m|,\frac{1}{h}\right), \qquad h\cos\theta-y\sin\theta > 0.\]

Here $h$ is height, $y$ is reach, and $k_m$ are medium wavenumbers. The angle is

\[\theta = \min\left(\frac{\pi}{6}, \arctan\frac{h}{2 \max(y, \epsilon)}, \min_{\zeta} \frac{\arg \zeta}{2}\right),\]

where $\epsilon$ is the machine epsilon of the coordinate type and $\zeta$ runs over the branch points $\pm\mathrm{j}k_m$ and the poles in the open first quadrant. These choices make the geometric weight $e^{-h\lambda}\cos(y\lambda)$ decay along the ray.

Use nominal values for all decisions. The scales are nonzero wavenumber magnitudes, finite positive formula scales, and the reciprocal geometric decay. Each scale contributes half, once, and twice its value. Each singularity contributes up to five positive points around its projection. Bridge gaps by factors of four, divide points by $R$, include zero and one, then sort and remove duplicates.

Arguments

• buffers: computation arrays with the integral block.
• wavenumbers: tuple of medium wavenumbers [1/m].
• poles: tuple of poles of the integrand [1/m]. Ignore zero entries.
• scales: tuple of formula-specific spectral scales [1/m]. Only finite positive values contribute.
• height: sum of source and receiver distances from the interface [m].
• reach: horizontal separation, with any receiver-radius allowance needed by the integrand [m].

Returns

• The complex ray $R e^{\mathrm{j}\theta}$, in [1/m]. The method also fills buffers.integral.points with dimensionless breakpoints.

Errors

• Throws DomainError when the geometric weight does not decay along the ray.

defaults(
    _::typeof(LineCableModels.Engine.integrate),
    _::Val{:quad}
) -> @NamedTuple{method::Symbol, rtol::Float64, atol::Float64, maxevals::Int64}

Return the declared controls of adaptive Gauss–Kronrod quadrature.

Arguments

• integrate: integral evaluation function.
• Val(:quad): adaptive quadrature method tag.

Returns

• method: :quad.
• rtol: relative tolerance, 1e-8.
• atol: absolute tolerance in the unit of the integral, 0.0.
• maxevals: evaluation limit, 10^7.

validate(
    section::NamedTuple,
    _::typeof(LineCableModels.Engine.integrate),
    _::Val{:quad}
) -> NamedTuple

Check the quadrature controls before an integral evaluation.

Arguments

• section: integration controls.
• integrate: consumer of the controls.
• Val(:quad): adaptive quadrature method tag.

Returns

• section, unchanged.

Errors

• Throws ArgumentError for a negative or nonfinite tolerance, two zero tolerances, or a nonpositive evaluation limit.

select(
    default::Pair{typeof(LineCableModels.Engine.integrate), Symbol},
    given::NamedTuple
) -> NamedTuple

Build an integration section from method defaults and supplied controls, then replace its method identifier with a Val tag.

Arguments

• default: pair of integrate and the default method identifier.
• given: supplied controls, with or without a method key.

Returns

• A named tuple with the method keys and their declared value types.

Errors

• Throws ArgumentError for an unknown method or an unknown control key.

initialize_buffers(
    _::Type{LineCableModels.Engine.SpectralIntegral},
    _::Val{:quad},
    _::Type{T},
    buffers
) -> Union{NamedTuple, DataStructures.SortedMultiDict}

Allocate the integral arrays. Coordinates use the type of the nominal value of T, and integrand values use Complex{T}.

Arguments

• SpectralIntegral: integral type whose arrays the computation needs.
• Val(:quad): adaptive quadrature method tag.
• T: numerical type of the formula.
• buffers: arrays already allocated.

Returns

• buffers extended with the integral block. It holds breakpoint and scale vectors plus mapped points, QuadGK segment arrays, and a zero prototype.

LineCableModels.Engine.EarthProperties

Compute the resistivity, relative permittivity, and relative permeability of earth and air.

Formulas

Formula{:constant}

Properties of a layer of the earth or of the air independent of the frequency $f$:

\[\rho(f) = \rho, \qquad \varepsilon_r(f) = \varepsilon_r, \qquad \mu_r(f) = \mu_r,\]

where $\rho$ is the resistivity of the layer [Ω·m], and $\varepsilon_r$ and $\mu_r$ are its relative permittivity and permeability [dimensionless].

Identifiers of the earth-properties formulas.

struct Formula{ID, O<:LineCableModels.Commons.FormulationOptions} <: LineCableModels.Engine.AbstractEarthProperties

Represent a frequency law for the properties of earth and air. Its options form a record of type FormulationOptions{Formula{ID}}, with ID the formula identifier.

• options::LineCableModels.Commons.FormulationOptions: Options of the formula.

Formula(
    ::Val{ID};
    kw...
) -> LineCableModels.Engine.EarthProperties.Formula

Build the formula named by ID with its declared defaults and supplied options.

Arguments

• Val{ID}: formula identifier from FORMULAS.

Keywords

• kw...: options declared by defaults for Formula{ID}.

Returns

• A Formula with a typed option record.

Errors

• Throws ArgumentError for an unknown identifier or option.

properties(
    _::LineCableModels.Engine.EarthProperties.Formula{:constant},
    layer::EarthLayer{T<:Real},
    f::Real
) -> NamedTuple{(:rho, :eps_r, :mu_r), <:Tuple{Real, Real, Real}}

Return the properties of the earth layer at the frequency f by the formula Formula{:constant}, whose equation is in the section "Formulas" of EarthProperties.

Arguments

• Formula{:constant}: frequency-independent properties law.
• layer: layer of the earth or of the air.
• f: frequency [Hz].

Returns

• rho: resistivity [Ω·m].
• eps_r: relative permittivity [dimensionless].
• mu_r: relative permeability [dimensionless].