An alternative method for field analysis in inhomogeneous domains

1. Introduction

Although final element methods (FEMs) are nowadays the standard tool for electrostatic and magnetostatic calculations, both in homogeneous and inhomogeneous materials, conformal transformations and Schwarz–Christoffel (SC) maps (Driscoll and Trefethen, 2002) in particular still deserve attention because of their capability to cope with convex boundary corners and relevant field singularities and to provide closed domains of computation. Starting from Costamagna and Di Barba’s study (2017) to Costa and Costamagna’s study (2019), in Section 2, we try solving inhomogeneous dielectric problems by matching boundary conditions (BC) at the interfaces between homogeneous subdomains, via two stage SC + BC procedures and refined boundary condition computations using Hilbert transforms. More generally, a question about the accuracy of methods for solving two-dimensional Laplace problems was raised not long ago (Trefethen, 2018), asking to compute the potential ψ at the point c = 0.99 + i0.99 within an L-shaped domain [0,2] × [0,2] minus its upper-right quarter with Dirichlet boundary condition ψ(z = x + iy) = x². For this problem the exact solution is known and is ψ(c) = 1.0267919261073… As reported in Gopal and Trefethen’s study (2019), about 20 experts responded with solutions obtained from FEM codes to 2–4 exact digits; only one reached six digits using 158997 5th-order triangular finite elements near the re-entrant corner. Boundary element method experts did better, up to eight digits, running their own codes exploiting powerful specialized quadrature routines. Corners, even those which are not very critical, still represent a significant computational issue, and the search for alternative methods to compare accuracy is justified indeed. We give our own answer to the question in Section 3.

2. Dielectric interfaces and Hilbert computations

In case of inhomogeneous materials, second stage FEM calculations have been introduced in Costamagna and Di Barba’s study (2017) to provide accuracy and computing speed, exceeding the performances of SC + finite difference (FD) methods already proposed [see for instance (Alfonzetti et al., 2001) and references therein]. It was shown that differences between capacitance values, computed in inhomogeneous dielectric for the geometry in Figure 1 by means of direct FEM, SC + FD and SC + FEM computation, can be reduced to about 0.1%. The geometry, which represents the half cross section of a shielded dielectric supported strip line, was selected as case study due to the singularity at the extreme of the strip line. Relative permittivity was ε_r = 9.5.

Further, the rectangular map of Figure 2, which is similar to the second stage map derived from Figure 1, was introduced as Figure 4 in the cited work, Section 2.3. The reason was to simplify the preparation of input data for FEM codes while maintaining all the significant characteristics of the case study; the dielectric interface at the right of the figure has been replaced by a perfect quarter circle and the other by a very similar vertical straight line segment. Moreover, the lengths of the sides and of the relevant dimensions have been chosen to provide integer ratios among them.

This way, the difference between the SC + FEM and SC + FD results was further reduced in Costamagna and Di Barba’s study (2017) to about 0.03% with a mesh of 548 × 1,189 nodes. This value and the previous one of about 0.1% can therefore be assumed as terms of comparison for the procedures discussed in this work, from which we hope to obtain, however, greater speed of calculation with respect to SC + FDs and greater ease of implementation with respect to FEMs. The map will be used here to exploit both the quarter circle shape of the right interface and reference data derived from it in Costamagna and Di Barba’s study (2017).

The circular Hilbert transform of a real-valued function f(t), 0 ≤ t < 2π:

In Costa and Costamagna’s study (2019), by using voltages obtained from the above inhomogeneous mesh to compute the stream function for the quarter disk subdomain via Hilbert transform, the capacitance agreed with the mesh fluxes at the interface within some 0.6%, and the very uneven behaviour of the local flux evaluated at the mesh nodes was again observed. Further, this error reduces to about 0.06% when comparing the above Hilbert transform value to the total mesh flux computed at the mesh boundary by means of the trapezoidal rule. Of course, formal calculations are possible in homogeneous dielectric, virtually error-free, considering a whole circle geometry obtained first by adding to the quarter circle its symmetrical on the right of a vertical magnetic wall, then the antisymmetric copy of the resulting semicircle, and eventually applying to the obtained full circle copies of the potential function on the original quarter circle with appropriate signs. By comparison, the results obtained with small numbers of discrete potential samples are noteworthy because with only 64 or 32 samples on any quarter circle accuracy was better than 10⁻¹³. This means that, when the geometry is perfect, the sampling of the potentials can be quite coarse.

Consider now the geometry in Figure 3, which still represents the above quarter disk, but now the aim is to simulate a general procedure, abandoning perfect circle geometries and considering instead a polygonal boundary of 128 sides lying on the previous quarter circle, 65 equal sides lying on the magnetic wall, horizontal in the figure, and the closing side of the electric wall, vertical in the figure. This geometry is mapped onto a unit disc via a first SC map onto the upper half plane and then from this onto the unit disk via a Cayley conformal transformation, which can be useful to control the final location on the full circle of all the vertices. Two different SC maps have been used: the first maps to infinity the sector centre, right vertex in Figure 3; the second maps to infinity a point on the left electric wall. Similar polygonal boundaries with different numbers of sides have been used without noticeable differences. Maps on circles are used to avoid vertices at infinity. Periodic potential functions are considered on the boundary to obtain the stream function conjugate to the potential via Hilbert transform computed via fast Fourier transform. The dielectric is homogeneous, and unitary permittivity was imposed to provide reference data by inspection; the potentials at the vertices are supplied by their abscissa, both along the quarter circle and the magnetic wall, and the true stream function is known. Both the small number of vertices and the naive procedure have been chosen precisely to see the possible limits in accuracy. The hypothesis of homogeneous dielectric and unitary permittivity is not limiting the general meaning of the analysis because the conditions are the same when matching BC at any dielectric interface; by imposing potential conditions, fluxes are computed solving the Laplace equation in a homogeneous domain with some well-defined permittivity.

We consider 4,096 sample points on the circle map, and the potentials applied to them have been derived by means of linear interpolation between the potentials at the vertices of the polygonal boundary. They are shown in Figure 4, with a thick dashed line for the first SC map and a thin dashed line for the second.

The derived stream functions, biased to obtain very small values along the magnetic wall, are shown in the same Figure 4 with solid lines. It is worth noting the almost perfect mirror behaviour of the solid line potential and its Hilbert transform, due to an astonishing spectrum with phases of π/4 or −3π/4 rad for any component of important amplitude. From inspection of geometry, the maximum value of the total flux must be equal to the maximum potential, and both stream functions (Figure 4, thick and thin solid lines) show accuracies of the order of 0.5% or 1%. This result is probably due in part to the sparse distribution of the mapped vertices on some region of the circle but mainly to the small oscillations of the Hilbert transforms, which impair accurate evaluations of the difference between their values at the ends of the magnetic wall or of the quarter circle. The oscillations are clearly perceived in Figures 5 and 6, which are enlargements of the left and upper parts of the solid line of Figure 4 for the case of mapping to infinity the right-angle vertex in Figure 3. Considering the behaviour on the magnetic walls (see for instance Figure 5), where a constant value was expected everywhere, a useful measure for both stream functions was found in calculating their average along a short segment, where the smallest oscillations are observed, to obtain a convenient bias and finally replacing the oscillating behaviour with a straight-line segment. This led to accuracies of about 0.5% and 0.35%.

More attention deserved the oscillation near the border between the quarter circle and the electric wall, see for instance Figure 6. A good result was again obtained in this case by means of considering the straight line passing through two points of the average oscillation trend, choosing these points at appropriate distance and replacing the maximum value of the oscillation with the ordinate of the straight line to its abscissa. This reduced the error in the maximum value of the stream function to about 0.014% instead of the previous 0.5%, with an accuracy analogous to that found in Costamagna and Di Barba’s study (2017) when comparing SC + FD and FEM results. Similar results have been obtained in the same way starting from the second SC map, potential and stream functions represented with dashed lines in Figure 4; from the previous 0.35%, the accuracy has been reduced to about 0.03%. The ramp behaviour of the stream function near the maximum value observed in Figures 5 and 6 is due to crowding phenomena of the polygon vertices mapped on the circle. When plotting the samples of the stream functions considering an abscissa running along the perimeter of the polygons of Figure 3, we find trends of the type shown in Figure 7 for the curve in Figure 6: three points with similar value of the function, mapped into the circle in almost coincident points, appear on different position along the perimeter. The right part of the figure shows the trend of the stream function along the electric wall, almost straight as expected.

It should be noted that returning to the potential along the quarter circle from the corrected stream functions leads to an oscillation behaviour analogous to that found in the above stream functions, and the correction in the maximum peak region has no notable influence, being the peak itself an effect of crowding. Therefore, a kind of corrected transform pair is conceivable. Correction of the oscillations in stream functions is of course necessary when they prevent accurate evaluation of partial capacitances along the sides of the original polygon.

A provisional conclusion drawn from the above observations can be that local calculations at the boundary can be more suitable than the integral calculations supplied by the Hilbert stream function if and only if flux values are reliably computed. In any case, the possibility of solving inhomogeneous dielectric problems by means of matching field behaviours at the boundaries of homogeneous subdomains looks very promising, avoiding cumbersome computations but nevertheless ensuring good accuracy, and this is the subject of the next section.

3. Alternative method: introductory considerations

The Conformal Boundary Difference Method (CBDM) is a numerical technique introduced in Costa’s study (2018), conceived for extending to multiply connected and inhomogeneous domains the scope of conformal mapping, insensitive by nature to corner singularities. To illustrate the underpinning idea, we answer the question proposed in the introduction and depicted in Figure 8.

Let f(t) be a conformal transformation mapping the unit disk D in the t-plane onto the physical domain ℙ in the z-plane; let also f(0) = c. For e^iθ ∈ ∂ D the Schwarz integral:

w=[0  2  2+i1+i  1+2i  2i];

spi=@(z)  real(z).  ˆ2;

f= diskmap(polygon (w),scmapopt(′tolerance′,1e-12));

f=centre (f,0.99+0.99i);

fun=@(x)  psi (f(exp(i∗x)));

psi−c=quadl (fun, 0,2∗pi,1e-12)/(2∗pi)

The code runs in one second on a laptop PC, returning ψ(c) = 1.0267919261074 exact to 13 digits.

The idea of determining field values by association of conformal mappings and the Schwarz integral is modified to allow BC other than those of the Dirichlet type, and this leads to the aforementioned CBDM. In short, the method runs as follows; for a potential field problem of the Laplace type:

(where u(z) is the potential to be determined; ℙ is the physical domain, each ℙ_m representing a homogeneous, simply connected subdomain; n⌢ is the outward normal to the boundary ∂ℙ_m; Γ_D and Γ_N ⊆ ∂ℙ_m) we can enforce constraints at all boundaries of the Neumann type Γ_N, namely, borders and interfaces between subdomains, by applying a suitable discretization and FDs there only. We can start with an arbitrary distribution of potential ψ(z) at {z_k}_m ∈ ∂ℙ_m and {z_i}_m ∈ ℙ_m, forming usual 5-points stencils, and converge to a solution via Successive Over-Relaxation (SOR) iterations. We do not really need a grid here to find and update the values of ψ: given a set of SC mappings {F_m} to the unit disk { D_m} such that F_m(ℙ_m) = D_m; at each step, we exploit the images t_i = F_m(z_i) in the unit disk and update the values of the electric potential ψ(z_i) = ψ(t_i) by means of the Schwarz integral. Then go back to the physical domain, update ψ(z_k) via the usual 5-point equation and repeat until some convergence criterion is met, that is, ψ(z) ≈ u(z). Therefore, the CBDM can be thought of as an improved and accelerated FDM; irregular contours pose less problems, and SOR iterations are restricted to boundaries. This is the overall picture, and we refer the reader to Costa’s study (2018) for a detailed description, covering measures to improve accuracy and computing speed.

4. Computations via the CBDM

The physical domain is first split as shown in Figure 1 for the application of the CBDM; then, a set of field solutions is obtained for the homogeneous case (no strip support, i.e. permittivityε = ε₀) with different values of the discretization step h, at the Neumann boundaries and between subdomains (h/8 at L/R interface). The electric potential distributions being determined, the capacitance values C₀ (h) are readily computed from the normal component of their gradient entering the outer conductor and compared with the “exact” value C₀ = 19.22491 pF/m obtained from a direct map to a rectangle (rectmap function of the SC-Toolbox). Notwithstanding the unavoidable loss of precision injected by discretization, the CBDM proves itself fast and accurate. For instance, with a boundary step h = 0.5 units (3rd row in Table 1) and no particular code optimization except for ordinary routine vectorization, the solution [1] takes a few seconds on a laptop PC, including conformal mappings and FD scheme generation, and the error of +0.39% on capacitance is the same order of magnitude of those affecting the best FEM results in (Costamagna and Di Barba, 2017; see Table 1, columns 1– 4 from the left).

Step sizes in Table 1 are by all means intentional; though at first sight they might look somewhat “coarse”, we want first of all guarantee that numerical SC direct and inverse maps be unaffected by the crowding phenomenon, [see (Driscoll and Trefethen, 2002), subsection 2.6 for a detailed description], and meet high accuracy, as reported by the SC-Toolbox diagnostic messages. The CBDM is well suited to the technique of Richardson extrapolation (RE) to zero error, and the capacitance C in particular is conveniently expanded via Taylor series as:

Thus, by taking a logarithm:

Table 2 shows the results obtained from a two-step RE for the inhomogeneous case (dielectric substrate having relative permittivity ε_r = 9.5), where an OoA n = 1 has been used again on the basis of the previous analysis. The final value differs by −0.11% from that obtained by the FEM (32.9654 pF/m) in Costamagna and Di Barba’s study (2017), but the former is likely to be a yet better estimation, given that it does not rely on meshes optimized for a homogeneous case and then transplanted “as is” to an inhomogeneous one. Apart from this small deviation, which does not really invalidate the excellent level of agreement, the point is that CBDM + RE turns the result (capacitance) obtained from conformal analysis in the particular homogeneous case from being a mere numerical reference into a means of determining a critical parameter (the OoA) for the refinement method in general situations. Moreover, this process is quick, neat and stable.

The Hilbert transform (1) represents another route to the calculation of capacitance, as in Laplacian fields the difference in stream function values φ (z₂) − φ (z₁) equals the total dielectric flux crossing the path from point z₁ to z₂, and just as ψ, it is invariant under conformal transformation. The electric potential distribution along the domain boundary corresponds through SC mapping to that along the image circle; here, after a suitable interpolation with evenly distributed points in the range [0, 2π], the Hilbert function available in MATLAB returns the stream values from the electric ones in a single line of code almost instantly. Regretfully, numerical problems arise at the points where the electric potential loses differentiability; refer to Figure 9.

These instabilities vary with the size of h and the number of points for the H-transform in a nearly unpredictable manner, their overall behaviour being dirigible only by the choice of a particular interpolation method. Although these issues are confined in the (wider or narrower) neighbourhood of such well-determined points, they render it almost impossible to pick them, and their closest neighbours, as the beginning/ending of a path for the calculation of charge flux, if accuracy and predictability must be met. Interpolation and regression methods have been tried, but at the cost of introducing sophisticated machinery that is not likely to cope with all possible cases, and they are not worthy of mention here except for the homogeneous Neumann sides where a numerical regression to a zero-order polynomial (straight horizontal line) fits naturally. What really comes to the rescue is the fact that φ is another scalar potential, and as such, it features all the peculiar properties granted by theory to this class of functions. See Figure 10; instead of considering the outer conductor and its uneasy points, the same result is obtained in a smoother and more elegant manner along the path a → b → c → d (or the like), which equally encloses the dielectric flux emanating from the strip. Actually, only the values of φ at the points are relevant; in either case, homogeneous or not, the capacitance C is readily found as:

Tables 3 and 4 show the results obtained from combining the CBDM + H-transform + RE, for the homogeneous and inhomogeneous case, respectively. Compare with Tables 1 and 2; although those results have been derived from the plain, trustworthy technique of integrating the electric potential gradient in regions where it features slow changes (namely, the outer conductor), the latter rely on deep analytical ideas that always work even close to field singularities. The two show high consistency in any aspect indeed; in particular, a slightly lower (−0.12%) value is confirmed for the capacitance in the inhomogeneous case with respect to that obtained from the FEM. In closing, we appreciate that our value is closer (+0.07%) to the best SC + FD result in Costamagna and Di Barba’s study (2017).

5. Discussion

The dielectric supported thin strip line geometry has been considered in the light of a new idea, that is matching conditions at the interfaces between homogeneous subdomains by means of easy calculations, to obtain both fast and undemanding, yet reliable computations. The CBDM has been found well suited to this end; resting on the theory of analytic functions and one-dimensional discretization, its results are easily processed by means of the Hilbert transform for dielectric flux and capacitance calculations, and of the Richardson Extrapolation for accuracy improvements. Problems arising from numerical treatment of non-smooth distributions can be circumvented by resorting to theoretical function properties; it elevates the role of results obtained via conformal mapping in homogeneous cases from side to key. Besides its standing as an independent analysis method, the CBDM can provide a reliable term of comparison for more standard means such as the FEMs.