© A. Skarlatos et al., Published by EDP Sciences, 2025

1 Introduction

The calculation of the transient field response in eddy-current problems involving arbitrary excitation is conventionally carried out by direct integration of the governing diffusion equation via a time-stepping scheme. Such approaches comprise single or multistep background differentiation (BDF) schemes or the Runge-Kutta method. The time-stepping approaches are well-established and have been successfully applied in combination with grid-based numerical methods based on volume meshes such as the finite elements method (FEM) for the solution of various partial differential equations. Nevertheless, their application is not straightforward when combined with other type of methods, including integral equation methods like the boundary element method (BEM) [1–3] or modal approaches [4,5]. As far as the latter is concerned, the conventional surface development bases must be complemented by an appropriate volumetric basis in order to account for the solution of the previous time-steps, a generalization which complicates the solution, in particular when the geometry departs from the simple multilayer pieces. In addition, time-stepping approaches are in principle sequential schemes, hence their parallelization is normally restricted to the spatial part through domain decomposition methods¹.

An alternative to direct time-integration are frequency domain (FD) methods. The calculation is carried out for each frequency of the excitation spectrum by solving the corresponding Helmholtz equation, which makes these approaches directly applicable to methods like BEM or mode-matching, and which were originally designed for the FD. Given that the monochromatic solutions are independent to each other, this class of methods is straightforwardly parallelizable both in space and time, thus offering an optimal use of the modern computer architectures.

The most well-known FD approach involves the Fourier Transform (FT) and its fast variant, the so-called fast Fourier Transform (FFT), where the problem is solved for a number of real frequencies. Nonetheless, the determination of the optimal spectrum discretization turns to be a challenging task for problems with high dissipation as the eddy-current problem [6,7]. Approaches based on a sampling in the Laplace or Z plane, on the contrary, prove to be more robust for such problems. Common point of these approaches is that the sample frequencies are in general located in the complex plane. The main challenge with these methods is how to invert the solution spectrum back into time domain in an efficient way requiring the least number of frequency samples (recall each frequency is associated to the solution of a field problem).

In this contribution, two alternative approaches based on the inversion of the Laplace and Z transform, respectively, will be compared for the solution of the eddy-current problem under step or pulse excitation. This type of excitations are important since it corresponds to the majority of excitation signals used in pulsed eddy-current applications, and because, once the response under step excitation is known, one can easily compute the corresponding signal under arbitrary excitation using Duhamel's integral [8]. The two approaches will be compared against the results obtained with a first order BDF scheme (BDF1), which will be used as reference.

2 The governing equation: weak formulation

The underlying physical problem is described by the diffusion equation for the magnetic vector potential A: $$\nabla \times \nu \nabla \times \text{A}+\sigma \frac{d\text{A}}{\mathrm{d}t}={\text{J}}_e$$(1)

where ν, σ is the magnetic reluctivity and the electrical conductivity of the medium, respectively, and J_e the excitation current. The magnetic vector potential is defined in terms of the magnetic flux density via the relation $$B=\nabla \times A.$$(2)

Since A is defined up to a gradient, it must be gauged. The usual choice is the Coulomb gauge, which imposes zero gradient of the magnetic potential.

Introducing a finite basis w_m(x) for the development of A(x, t) $$\text{A}\left(\text{x};t\right)\approx {\displaystyle {\displaystyle \sum }_{m=0}^{N_e}}{a}_m{\text{w}}_m\left(\text{x}\right)$$(3)

and projecting (1) onto the same basis according to the Galerkin approach, we end up with a set of ordinary differential equations which constitute the weak form of (1) $$\left(\text{K}+{\text{M}}_{\sigma }\frac{d}{\mathrm{d}t}\right)\text{a}(t)={\text{j}}_e(t)$$(4)

where a = [a₁, ..., a_Ne]^T, j_e= [j_e1,...,j_eNe]^T are column vectors compiling the development coefficients of A and J, respectively, and K, M_σ stand for the corresponding stiffness and mass matrices. The explicit form of the approximation basis w_m(x) as well as the details of M_σ can be found in the literature and will not be repeated here. In the context of the of studies considered in this article we choose to use the modal approach reported in [8,9]. The same mathematical analysis applies, however, equally well in the context of the FEM or the BEM method [1–3].

3 Solution in the Laplace plane

We assume that the initial value of the magnetic potential is zero throughout the solution domain. With application of the Laplace transform $${\displaystyle \tilde{\text{a}}}\left(s\right)={\displaystyle \underset{0}{\overset{\infty }{{\displaystyle \int }}}}\text{a}\left(t\right)e^{-st}dt$$(5)

(4) becomes $$\left(\text{K}+{\text{s M}}_{\sigma }\right){\displaystyle \tilde{\text{a}}}\left(s\right)={{\displaystyle \tilde{\text{j}}}}_e\left(s\right)\text{.}$$(6)

All variables are assumed have been calculated in the complex (Laplace) domain s. The problem is equivalent to the harmonic one with complex radial frequencies $$\omega =s/i\text{.}$$(7)

The transient response can then be obtained by the harmonic solution by applying the inverse Laplace transform, which is defined via the contour integral $$\text{a}\left(t\right)=\frac{1}{2\pi i}{\displaystyle \underset{c-i\infty }{\overset{c-i\infty }{{\displaystyle \int }}}}{\displaystyle \tilde{\text{a}}}\left(s\right)e^{st}ds\text{.}$$(8)

The integration path in (8) is known as the Bromwich path, where c is chosen such that c > maxRe[p_ℓ], with p_ℓ standing for the F(s) poles, in order the integral to converge.

With an appropriate sampling of the integration path, the weak formulation of (6) can be readily solved for each complex frequency s_n, n = 1, ..., N_s. Once the frequency domain response has been calculated one can get the transient response by simply calculating the inverse Laplace transform given in (8).

The direct computation of the integral can be however computationally inefficient since a large number of points may be needed for good accuracy. Since each frequency is translated to the solution of a field problem, it is easily understood that the computational time can rise quickly. Furthermore, the estimation of the path truncation is not apriori known.

To avoid the complexities associated with the numerical evaluation of the contour integral and given that in many practical situations, where the field transients are required, we are primarily interested to the short-time response, other methods may be more efficient. This is particularly the case when the impulse or step response of a particular variable (such as the coil EMF) is required. The response to an arbitrary excitation can then be easily obtained via convolution with the input signal.

A family of methods that works particularly well for the computation of short-time signals, is the Zakian class of methods, which was originally proposed in the 60s by Zakian and has been subjected to several variations thereafter [10]. The main idea behind these methods consists in computing approximations of the delta function via exponentials and calculating the convolution with the sought time function, which reduces to the computation of a weighted sum of function samples in the Laplace plane. It can be shown that this method is a special case of the Gaussian quadrature method.

Starting from the identity $$\text{a}\left(t_n\right)\equiv {\text{a}}_n={\displaystyle \underset{0}{\overset{\infty }{{\displaystyle \int }}}}\delta \left(t-t_n\right)\text{a}\left(t\right)dt$$(9)

where t_n is an arbitrary sample of the function domain and assuming that the delta function can be approximated by the following series of exponentials $$\delta \left(t-t_n\right)\approx t_n^{-1}{\displaystyle {\displaystyle \sum }_{j=1}^N}{K}_j{e}^{-{\text{a}}_{j}t/t_n}$$(10)

we obtain upon substitution $${\text{a}}_n\approx t_n^{-1}{\displaystyle {\displaystyle \sum }_{j=1}^N}{K}_j{\displaystyle \tilde{a}}\left(a_j/t_n\right)$$(11)

where ã (a_j/t_n) is the Laplace transform of the solution evaluated at the sampling frequencies a_j/t_n. The approximation can attain with arbitrary accuracy since A_mn → A_m(t) as N → ∞. The series coefficients K_j and the exponents a_j in (10) admit different values depending on the way of choosing the quadrature points resulting in different variants of the method.

In the original Zakian formulation, K_j and a_j are determined via Padé approximations for the exponential function [11]. The resulting values for the first 5 terms are given in Table 1. It turns out that 5 quadrature points provide a satisfactory accuracy for most of the practical cases.

A popular variant of the original Zakian method, which is commonly used for the calculation of impulse or step response is the so-called Gaver-Stehfest method, which proposes the evaluation sum [10,12]: $${\text{a}}_n\left(x\right)\approx \frac{\ln 2}{t_n}{\displaystyle {\displaystyle \sum }_{j=1}^N}{K}_j{\displaystyle \tilde{\text{a}}}\frac{\left(j\text{ln}2\right)}{t_n}$$(12)

where K_j are calculated by the formula $$K_j={\left(-1\right)}^{\frac{N}{2}+1}{\displaystyle {\displaystyle \sum }_{k=\frac{j+1}{2}}^{min\left(j,N/2\right)}}\frac{k^{N/2}\left(2k\right)!}{\left(N/2-k\right)!k!\left(k-1\right)!\left(j-k\right)!\left(2k-j\right)!}.$$(13)

The number of points N used for the approximation varies between 10 and 20 depending on the precision one wishes to achieve. For the most cases N = 10 provides a sufficient accuracy, and this choice is also adopted in this work. Although this number of points is higher than the one of the original Zakian formulation, and hence it needs (at least) the double number of field calculations, both variants may have different performances in terms of precision and long-time behavior depending on the signal.

An important point to highlight here is that the choice of the time samples is completely arbitrary. This is particularly convenient when calculating time signals such as the Dirac or step response where the main information is stored in the short time part of the signal. In such cases, one needs to proceed to inhomogeneous sampling of the time axis (e.g. logarithmic) for optimal performance. An illustration of the frequency point distribution for logarithmic sampling of the time axis using the two approaches is given in Figure 1. This point will become more clear through some examples in the results section.

4 Solution using the Z-transform

Instead of transforming directly the diffusion equation (4) into the Laplace domain, we can introduce an intermediate step by discretizing the temporal derivative of (4) using a homogeneous time grid t_j= jΔt with j = 0,1,...,N_t. Thereupon, we will be able to work in the Z-plane instead of the Laplace plane.

Following common practice in case of parabolic problems, the time derivative is discretised using a nth order BDF scheme yielding the following recursive formula $${\text{Ma}}_n+\frac{\sigma }{\mathrm{\Delta }t}{\displaystyle {\displaystyle \sum }_{j=0}^n}{\gamma }_{n-j}{\text{a}}_j=\frac{\sigma }{\mathrm{\Delta }t}{\text{j}}_{e_n}$$(14)

where γ_j are the coefficients of the multistep scheme. In the particular case of a first order scheme (implicit Euler), we have γ₀= 1, γ₁= 0, and γ_j= 0 for j < 1.

Introducing the unilateral Z-transform $${\displaystyle \tilde{\text{a}}}\left(z\right)={\displaystyle {\displaystyle \sum }_{n=0}^{\infty }}{\text{a}}_n{z}^{-n}$$(15)

and applying (15) in (14), we yield $$(\text{M}+\frac{\sigma }{\mathrm{\Delta }t}\sum _{j=0}^n{\gamma }_{n-j}{z}^{-j})\stackrel{}{}\text{a}~(z)={\stackrel{}{}\text{j}~}_e\text{(}z)$$(16)

which is the discrete analogue of the harmonic problem with complex radial frequencies $$\omega =\frac{\sigma }{i\mathrm{\Delta }t}{\displaystyle \sum _{j=0}^n}{\gamma }_{n-j}{z}^{-j}\text{.}$$(17)

Having calculated A(x; z) the corresponding time sequence A_n(x) can be retrieved by applying the inverse Z-transform, which in its general form reads $${\text{a}}_n=\frac{1}{2\pi i}\underset{C}{\int }\stackrel{}{}\text{a}~(z^{n-1})dz$$(18)

with C being a closed contour in the complex z plane, which encloses all function poles. For the open diffusion problem (conductive pieces in air) we are interested here and for finite length excitation signals it can be shown that all poles are included in the unit circle when the modal approach is applied for the inversion of the topological operators.

Under these circumstances, (18) can be evaluated in a straightforward manner by choosing the unit circle |z| = 1 as integration contour C and carrying out the integration using the trapezoidal rule. Note that by making use of the z-plane symmetry, i.e. ${\displaystyle \overline{F}}\left(z\right)=F\left({\displaystyle \overline{z}}\right)$, where the bar symbol denotes complex conjugate, only the points at the upper half-circle (Re[z] > 0) need to be calculated, the remaining ones at the lower half-circle being obtained by simply applying the aforementioned symmetry property. Choosing a discretization of the integration contour $z_k=e^{2\pi i\frac{k}{N_f}}$ with k_f= 1, ..., N_f and applying the trapezoidal rule for the numerical computation of the integral, we obtain [13] $${\text{a}}_n(\text{x})=\frac{2}{N_f}\sum _{k=1}^{N_f}\tilde{\text{a}}(z_k)z_k^n.$$(19)

The sampling points for the evaluation of (19) are shown in Figure 2.

It is worth mentioning that (19) is the equivalent of the discrete Fourier transform on the unit circle, therefore the FFT algorithm can be applied for accelerating the calculation.

Fig. 1 Location of the complex frequencies for the numerical evaluation of the inverse Laplace transform with the Zakian and the Stehefest method for a logarithmic sampling of the time axis.

Fig. 2 Sampling points for the numerical evaluation of the inverse Z-transform. Only the samples at the upper half-circle need to be numerically evaluated (filled markers), which reduces the number of problems solutions by a factor of two.

5 Numerical results

The above-discussed approaches have been tested for the calculation of the transient eddy-current field in an infinite conducting plate. The plate is electromagnetically excited by a cylindrical coil positioned above its upper interface with its axis being normal to the plate interface as shown in Figure 3. The plate electrical conductivity σ = 5 MS/m and its relative magnetic permeability μ_r is allowed to vary between 1 (non-magnetic) and 150, which corresponds to the average permeability of a carbon structural steel. The plate thickness is taken equal to 1 mm. The coil characteristics are summarized in Table 2. The coil lift-off is taken equal to e = 1 mm.

The coil inner and outer radius is r₁= 1 mm and r₂= mm respectively, its length is l = 2 mm and is wound with 336 turns. The field is observed at a sample point underneath the plate as shown in the figure. The point coordinates are r_o= 4.65 mm and z_o= −3.0 mm. The specific point has been selected far from the coil and the plate, where the field is weaker, in order to provide a good precision test.

In the first example, we calculate the eddy-current response upon step current excitation. Although the step function does not present a physically feasible current waveform, it makes an interesting stress-test for the tested algorithms due to the steep transient it produces at short times. It is also a case of practical importance since the step response can be used for the calculation of the system response upon arbitrary excitation by evaluation of the Duhamel integral, as already mentioned [8]. Note that in order to get a very short raising time, the plate is considered being non-magnetic for this particular case.

The magnetic field is calculated at a sample point underneath the plate at 0.25 mm from the coil axis and 2 mm below its lower interface. The transient response has been calculated using the three inversion approaches presented above: contour integration in the Z plane, which will be conventionally referred to as IZT henceforth, and calculation of the inverse Laplace transform using the Zakian and the Stehfest methods. As reference solution, we use a calculation based on the finite integration technique (FIT) implementation with a BDF1 scheme for numerical integration [14]. The comparison of the Laplace/Z domain results with the reference solution is given in Figure 4.

In order to test the impact of the discretisation step on the IZT solution several evaluations have been carried out using a different number of sampling points, namely 30, 100 and 300 samples. In order to achieve an adequate precision, it is recommended to proceed to an oversampling in the Z-domain, with the number of sample frequencies being higher than the time samples by a factor of at least 2. For the computations of this article, an oversample factor of 5 was used, which is translated to 150, 500 and 1500 frequency points repsectively. The reference BDF1 solution is obtained using the finest discretisation with 300 time-steps. For the Zakian and the Stehfest approaches a logarithmic sampling with 10 time samples is applied. The number of frequency samples for the two methods was 5 and 10 scomplex frequencies per time sample respectively.

From the comparison of the results, it becomes clear that the steep decay of the field signal at short times poses restrictive constrains to the time discretisation for the IZT approach. Such kind of situations are however well adapted for the Zakian/Stehfest schemes, since the evolution of the time signal is concentrated in a narrow zone in the close to zero of, which makes the use of inhomogeneous time discretisation very efficient.

The second example involves a current waveform with finite time constant, which represents a typical case of a feed current from a real circuit. The pulse duration is 1 ms, its duty-cycle is 50 % and its time constant is 0.1 ms. The plate in this case is ferromagnetic with a relative permeability equal to 150. The plate conductivity is the same with the previous example. The comparison of the field transients, calculated with the three approaches, against the reference solution is shown in Figure 5.

The situation is inverted with respect to the previous case we can see the Zakian methods fail to reproduce the slope change at the negative pulse edge. The IZT signal on the contrary lies very close to the reference using only 30 time samples. The important performance difference of the IZT approach between the first and the latter are attributed to much smoother slopes of the second example. This inefficiency is of the same nature with the difficulties encountered in time-stepping schemes when treating such problems, which is not surprising given that the method resides on a BDF scheme as explained above.

Fig. 3 Problem configuration. Infinite conducting plate excited by a cylindrical coil with normal to the interface axis. The field is calculated at a point at the opposite side of the plate.

Fig. 4 (Up) radial Bρ and (down) axial Bz magnetic flux density at the observation point as function of time. The transient response calculated with the Stehfest and the IZT approaches are compared against the BDF1 solution.

Fig. 5 (Up) radial Bρ and (down) axial Bz magnetic flux density at the observation point as function of time. The transient response calculated with the Stehfest and the IZT approaches are compared against the BDF1 solution.

6 Conclusions

Two different methods for the solution of the diffusion equation by passing from the complex (Laplace/Z) plane are compared. The approach is applicable to FD solvers without modifications and allows a seamless parallelization. The transient signals are reconstructed by applying an inverse transformation, which for the case of the Laplace plane is based on an approximate development of the dirac function, whereas the inverse Z-transformation is performed numerically via integration on the unit circle. The two approaches can be considered as complementary since the first is proved to be better suited for calculating the short-time response whereas the latter is capable of treating arbitrary signals in intermediate and long times. Overall the Z-based approach is more robust. Hybridization of the two methods can be of interest in cases where short-time features are important without compromising the long-time robustness.

Funding

All authors report no external funding for this research.

Conflicts of interest

All authors declare that they have no conflicts of interest.

Data availability statement

No data are associated with this article.

Author contribution statement

Anastassios Skarlatos: Conceptualization, Methodology, Investigation, Formal analysis, Validation, Visualization, Writing original draft, Writing review & editing. Mark Bakry: Conceptualization, Methodology, Investigation, Formal analysis, Writing original draft, Writing review & editing.

References

All Tables

All Figures

	Fig. 1 Location of the complex frequencies for the numerical evaluation of the inverse Laplace transform with the Zakian and the Stehefest method for a logarithmic sampling of the time axis.
In the text

	Fig. 2 Sampling points for the numerical evaluation of the inverse Z-transform. Only the samples at the upper half-circle need to be numerically evaluated (filled markers), which reduces the number of problems solutions by a factor of two.
In the text

	Fig. 3 Problem configuration. Infinite conducting plate excited by a cylindrical coil with normal to the interface axis. The field is calculated at a point at the opposite side of the plate.
In the text

	Fig. 4 (Up) radial Bρ and (down) axial Bz magnetic flux density at the observation point as function of time. The transient response calculated with the Stehfest and the IZT approaches are compared against the BDF1 solution.
In the text

	Fig. 5 (Up) radial Bρ and (down) axial Bz magnetic flux density at the observation point as function of time. The transient response calculated with the Stehfest and the IZT approaches are compared against the BDF1 solution.
In the text