FACTA UNIVERSITATIS Series: Electronics and Energetics Vol. 34, No 1, March 2021, pp. 141-156 https://doi.org/10.2298/FUEE2101141D © 2021 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper HIGH-ACCURACY QUASISTATIC NUMERICAL MODEL FOR BODIES OF REVOLUTION TAILORED FOR RF MEASUREMENTS OF DIELECTRIC PARAMETERS Antonije Djordjević1,2, Dragan Olćan1, Jovana Petrović1, Nina Obradović3, and Suzana Filipović3 1University of Belgrade – School of Electrical Engineering, Belgrade, Serbia 2Serbian Academy of Sciences and Arts, Belgrade, Serbia 3Institute of Technical Sciences of the Serbian Academy of Sciences and Arts, Belgrade, Serbia Abstract. We have developed rotationally symmetrical coaxial chambers for measurements of dielectric parameters of disk-shaped samples, in the frequency range from 1 MHz to several hundred MHz. The reflection coefficient of the chamber is measured and the dielectric parameters are hence extracted utilizing a high-accuracy quasistatic numerical model of the chamber and the sample. We present this model, which is based on the method- of-moments solution of a set of integral equations for composite metallic and dielectric bodies. The equations are tailored to bodies of revolution. The model is efficient and accurate so that the major contribution of the measurement uncertainty comes from the measurement hardware. Key words: dielectric measurements, electromagnetic modeling, method of moments, bodies of revolution 1. INTRODUCTION The key parameter for characterization of a linear, isotropic dielectric material is the relative complex permittivity and its dependence on frequency. There exist many methods for measuring the permittivity [1], [2]. It is beyond the scope of this paper to present and compare these techniques, so that we give only a brief overview. For measurements at frequencies up to several hundred megahertz, the most commonly used technique is based on the measurement of the capacitance of a parallel-plate capacitor, where the sample is inserted between the capacitor electrodes. This method assumes that the electromagnetic field within the measured sample is quasistatic, which imposes a high- Received September 2, 2020; received in revised form October 16, 2020 Corresponding author: Dragan Olcan School of Electrical Engineering, Kralja Aleksandra Blvd. 73, 11000 Belgrade, Serbia E-mail: olcan@etf.rs 142 A. DJORDJEVIĆ, D. OLĆAN, J. PETROVIĆ, N. OBRADOVIĆ, S. FILIPOVIĆ frequency limit on the method. At high frequencies, this technique has a drawback due to the strong electromagnetic coupling with the environment. Hence, shielding is required. Another potential drawback is that commercially available meters [1] require large- diameter samples (15 mm or more). For broadband measurements at microwave frequencies (above around 1 GHz), open coaxial lines [3] or waveguides [4] can be used. Parameters of sheet materials can be estimated by measuring the transfer between two antennas [5]. All these techniques require relatively large samples. In yet another set of techniques, a material sample is inserted into a coaxial line or a waveguide [6]. The sample is relatively large and has to be machined according to the shape of the coaxial line, viz. the waveguide. For measurements of dielectric substrates, other techniques can be used (e.g., [7]), which also require a special shape of the dielectric or a particular metallization pattern on it. Narrowband measurements are performed in resonators. They are convenient for low- loss materials and can be used for measurements of anisotropic dielectric materials [8], but they provide data only for discrete frequencies. In our research, we primarily deal with ceramic materials. We utilize disk-shaped samples, which are relatively small due to the restricted available mass of starting components used for sintering: the diameter (d) is in the range from 4 mm to 12 mm, whereas the height (h) is between 1 mm and 4 mm. The samples are too small for the standard test equipment based on the parallel-plate capacitor. Further, their shape and size do not fit into the available test equipment for standard measurements at microwave frequencies. Hence, for measurements in a wide frequency range (1 MHz–5 GHz), we have developed several coaxial chambers. The first prototype is described in [9], whereas an improved design of the chamber is shown in Fig. 1. The dielectric sample is pressed between a plate and a plunger, both made of brass. Using a vector network analyzer (VNA), we measure the reflection coefficient at the SMA (SubMiniature version A) connector and, hence, evaluate the input admittance of the chamber. On the other hand, we utilize a numerical electromagnetic model of the chamber with the sample. In the model, we optimize the dielectric parameters of the sample in order to match the measured admittance. SMA connector BrassTeflon Air Dielectric sample Plunger Plate d h VNA calibration plane Shifted reference plane Fig. 1 Cross-section of coaxial chamber High-accuracy Quasistatic Model for Bodies of Revolution 143 We have selected a coaxial structure because it is electromagnetically closed, and thus well shielded from the environment. Note that the chamber, with the inserted sample, is practically a rotationally symmetrical structure, i.e., a body of revolution (BoR). In practice, the measurement structure does not have a perfect rotational symmetry. This is not critical for lower frequencies because the chamber input admittance, as a function of the sample position, has a stationary point when the sample is in the middle of the chamber. However, at higher frequencies, when resonances of the chamber and the sample occur, the positioning is critical because modes with asymmetrical field distributions can be excited. In order to facilitate positioning of the sample, we use three thin screws that protrude through the chamber wall. After inserting the dielectric sample, the screws hold it in the required position. The screws are removed once the plunger presses the sample so that they do not influence the measurements. At lower frequencies, up to around 500 MHz, which we consider in this paper, the dimensions of the coaxial structure are relatively small compared to the wavelength. (The frequency limit depends on the dimensions and the relative permittivity of the sample.) Hence, the quasistatic approximation can be used for the analysis. Assuming a time-harmonic electromagnetic field [10], the equations involved in the analysis are formally the same as for the electrostatic fields. The differences from the electrostatics are: (a) phasors are involved, for the field sources (charges), electric scalar-potential, and the electric-field vector, and (b) the complex permittivity is used to characterize dielectrics. Such an approach enables analysis of lossy dielectrics. Note that losses in the metallic parts of the chamber have a negligible influence on the overall results of our measurements, which we have verified experimentally and computationally. The structure shown in Fig. 1 belongs to the class of structures that consist of metallic (conductive) regions and piecewise-homogeneous dielectric regions [11]. The electrostatic (quasistatic) analysis of the chamber cannot be performed analytically, but only numerically. To that purpose, various methods can be used, like the method of moments (MoM) [12], the finite-element method (FEM) [13], the method of fictive charges [14], the method of equivalent electrode [15], etc. Based on the MoM and the FEM, methods have been developed for the analysis of arbitrary two-dimensional (2-D) and three-dimensional (3-D) structures. Also, several commercial electrostatic solvers for arbitrary 2-D and 3-D structures are available, e.g., [16]–[18]. In the implementation of such general 3-D solvers, the analyzed structure is segmented without taking into account the rotational symmetry. Consequently, the required computer resources (memory and CPU time) are substantially larger than if a BoR solver were used (where the rotational symmetry of the sources and fields is utilized), resulting in non-optimal running time and even jeopardizing the accuracy due to oversized systems of equations. Unfortunately, there is no commercial simulator for the electrostatic analysis of BoRs. Also, in the open literature we could not find papers devoted to the electrostatic analysis of arbitrary BoRs (which consist of metallic and dielectric regions). Only very few older papers partly deal with this topic, e.g., [19] and [20], but their scope is limited because they are related to the analysis of slender conductors, viz. oblate dielectric bodies. In both papers, uniform asymptotic expansions are used. Tailoring the analysis method to BoRs is important because it can be substantially faster compared to the conventional analysis of 3-D structures. The speed is important for our applications, because many analysis cycles are involved in the optimization. The accuracy of 144 A. DJORDJEVIĆ, D. OLĆAN, J. PETROVIĆ, N. OBRADOVIĆ, S. FILIPOVIĆ the method is even more important because we want to achieve negligible influence of numerical errors on the overall measurement uncertainty. Hence, we have been motivated to develop a new method for precise and efficient quasistatic (electrostatic) analysis of arbitrary BoRs, which is described here. The paper is organized as follows. In Section 2, the numerical method is described. In Section 3, some benchmark numerical results are presented. Section 4 illustrates the implementation of the proposed method on actual measurements. The paper is concluded with Section 5. 2. NUMERICAL METHOD We consider a BoR structure (Fig. 2) that consists of charged conducting (metallic) bodies (electrodes) and electrically-neutral dielectric bodies (which collectively constitute a piecewise-homogeneous, isotropic dielectric medium). The axis of symmetry (revolution) is z. The generatrix of the BoR is in the right-hand part of the Oxz plane (x  0). Hence, x is the distance from the axis of symmetry. The operating frequency is f. Based on our past experience, as the preferred technique for the numerical analysis, we have selected the integral-equation approach, along with the MoM. We follow a similar path as in [16], [18]. x0 Conductor 2 Conductor 1 Vacuum (air) Axis of revolution (z) Dielectric 1 er1 Dielectric 2 er2 sb sb s sb s 'r 'dl r Source element Field point sb 1V 2V xu zu Fig. 2 Example of BoR consisting of conductors and dielectrics 2.1. Integral Equations First, we replace the conducting bodies by free surface charges (whose density is s) and the dielectric bodies by bound surface charges (whose density is sb). All these charges are assumed to be in a vacuum. The electric scalar-potential and the electric-field High-accuracy Quasistatic Model for Bodies of Revolution 145 vector of these charges are the same as in the original (analyzed) system. The reason for introducing these surface charges is to homogenize the medium, so that the potential and the electric field can be evaluated using the standard integral relations for a vacuum. We collectively refer to these surface charges as the total charges, whose (phasor) surface density is st. At an interface (boundary) between a conducting body and a vacuum, the total charges comprise only the free charges, i.e., st = s. At an interface between a conducting body and a dielectric body, we have st = s + sb. At an interface between a dielectric body and a vacuum, there are only bound charges, so that st = sb. Finally, at an interface between two dielectric bodies, there are also only bound charges. In this case, we write st = sb and assume that sb is the sum of the densities of the bound charges of these two dielectrics. Assuming rotationally-symmetric charge distributions, their (phasor) electric scalar- potential at a field point defined by the position-vector r = xux + zuz (where ux and uz are the unit vectors of the Cartesian coordinate system) is given by [21]  e = ' BoRst 0 'd)',()'( 1 )( C lgV rrrr , (1) where 12 0 F 8.8541878128 10 m e −   is the permittivity of a vacuum and )( ' )',(BoR mK q x g  =rr (2) is BoR Green’s function. Further, r' = x'ux + z'uz defines the location of the source (i.e., the element dl'), K(m) is the complete elliptic integral of the first kind, q = (x + x')2 + (z − z')2, and m = 4xx' / q. The BoR generatrix line C ' defines boundaries of all conducting and dielectric bodies and st is an unknown function of the position along C ', i.e., a function of a local coordinate l'. The reference point for the potential is at infinity. Note that the kernel of (1) becomes singular when r = r'. The singularity is logarithmic and integrable. The corresponding (phasor) electric field is Vgrad−=E . We formulate a set of integral equations for st based on the boundary conditions. The first part of this set is based on the boundary condition for the potential at the surfaces of electrodes. Each conducting body is equipotential. We denote the number of the conducting bodies by Nc and assume to know their potentials, Vi, i = 1,...,Nc. Consequently, when the field point is on the surface of a conducting body whose potential is Vi, we have an integral equation of the form V(r) = Vi, i.e., i C VlmK q x =   e  ' st 0 'd)( ' )'( 1 r , c,...,1 Ni = . (3) The second part of this set of integral equations is based on the boundary condition for the normal component of the electric field at the dielectric-to-dielectric interfaces (Fig. 3). We include here interfaces between any two dielectric bodies, as well as between a dielectric body and the surrounding vacuum. This boundary condition yields 146 A. DJORDJEVIĆ, D. OLĆAN, J. PETROVIĆ, N. OBRADOVIĆ, S. FILIPOVIĆ 0 st 2r 1r 211 1 e  =         e e −uE , (4) where E1 is the electric field in the first dielectric just at the boundary, u21 is the unit vector perpendicular to the boundary surface (directed from the second dielectric toward the first dielectric), st = sb is the charge density at the interface, and er1 and er2 are the relative complex permittivities of the two dielectrics. Note that E1  u21 = En1 is the normal component of the electric field in the first dielectric, so that 21n11n uE E= . Dielectric 1 er1 Dielectric 2 er2 stsb = 21u1 E A B D n Fig. 3 Boundary surface between two dielectrics The Cartesian components of the electric field in Fig. 2 are [21] ( )  − +−−−  e = ' 2/3st 0 'd )1( )'('2)()1()('2 ')'( 1 )( C x l mmq xxmxmEmmKx xE rr , (5)  − −  e = ' 2/3st 0 'd )1( )()'( ')'( 1 )( C z l mq mEyy xE rr , (6) where E(m) is the complete elliptic integral of the second kind. These expressions contain harder singularities compared to (1), because the kernels in (5) and (6) come from the derivative of Green’s function. Nevertheless, the technique for the evaluation of integrals described in Subsection 2.3 handles well even these integrals. An alternative approach is to compute the electric field by numerical differentiation. We evaluate V(r) at two points (A and B) on u21, which are close to the boundary surface (Fig. 3) and calculate the normal component of the electric field as En1  (VA − VB) / Dn. The distance Dn has to be carefully chosen in order to maximize the accuracy of computations. If Dn is too small, the error of subtracting two similar numbers (VA and VB) dominates. If nD is too large, the error of replacing the differentiation by differencing becomes pronounced. Our primary goal is to numerically evaluate the matrix of electrostatic-induction coefficients [B] [22]. We consider here a system that consists of two conductors (such as the one shown in Fig. 2). The conductor free charges (Q1 and Q2) and potentials (V1 and V2) are related as Q1 = b11V1 + b12V2, Q2 = b21V1 + b22V2, where bij, i, j =1, 2, are the electrostatic-induction coefficients, so that       = 2221 1211 ][ bb bb B . High-accuracy Quasistatic Model for Bodies of Revolution 147 (Due to reciprocity, bij = bji.) A generalization to a system with an arbitrary number of conductors is straightforward. In order to compute the elements of the matrix [B], we assume that one conductor is at a certain non-zero potential (e.g., 1 V) and all other conductors are at a zero potential. We numerically evaluate the free charges of the conductors. Hence, the elements of one column of [B] are easily calculated. This procedure is repeated for all conductors. For our measurements, we also need the matrix of partial capacitances [C] of the analyzed system. For a two-conductor system, the matrix [C] is defined in terms of the elements of the matrix [B] as       +− −+ =       212221 121211 2221 1211 ][ bbb bbb cc cc C . 2.2. Method-of-Moments Solution The complete set of integral equations is solved numerically using the MoM. As the basis (expansion) functions, we implement one of the simplest approximations for the distribution of the total surface charges: the piecewise-constant (staircase, pulse) approximation. To that purpose, we divide the contour C ' into a number of straight-line segments. (In the general case, each segment corresponds to a right conical frustum, which may degenerate into a right cylindrical frustum or a flat circular ring.) We assume that st is constant along a segment, though yet unknown. In order to provide a high accuracy and at the same time minimize the number of unknowns, we take the segments to be denser in regions where we expect faster variations of st, e.g., near edges of conductor and dielectric bodies. The distribution of the segments is defined in a way similar to [16]. For testing, we implement the Galerkin procedure: we integrate the left-hand side and the right-hand side of each integral equation over the surface of one-by-one frustum. As the result, the elements of the part of the MoM matrix that corresponds to the boundary condition (3) have the form              e = i jC ii C j j ij lxlmK q x Z d2d)( 1 0 , sc,...,1 Ni = , s,...,1 Nj = , (7) where the index i corresponds to the field segment (i.e., the segment where the boundary condition is implemented) and j corresponds to the source segment. Further, Ci denotes the field segment and Cj denotes the source segment. Nsc is the total number of segments for conductors and Ns = Nsc + Nsd is the total number of segments (unknowns) for the whole structure, where Nsd is the total number of segments for dielectric-to-dielectric interfaces. Finally, in (7), q = (xi + xj) 2 + (zi − zj) 2 and m = 4xixj / q. The elements of the remaining part of the MoM matrix, which corresponds to the boundary condition (4), i.e., Zij, i = Nsc + 1,..., Ns, j = 1,..., Ns, have a similar form, except that, in their derivation, the integrals in (5) and (6) are used instead of the integral in (3). We have found that high-contrast dielectrics (e.g., if the relative permittivity of one dielectric is 1000 and the other dielectric is a vacuum) tend to destabilize the system. In 148 A. DJORDJEVIĆ, D. OLĆAN, J. PETROVIĆ, N. OBRADOVIĆ, S. FILIPOVIĆ order to solve this problem, we add an equation that requires that the total bound charge of the system is zero [23]. We solve the resulting system of linear equations by the LU decomposition and back substitution. Thus we obtain the total charge densities on the segments. All BoR conductors are assumed to have finite thicknesses. Hence, we evaluate the free-charge density of a segment simply as s = erst, where er is the relative complex permittivity of the adjacent dielectric. Knowing the free-charge densities, we evaluate the free charges of the conductors and, hence, calculate the matrices [B] and [C]. 2.3. Evaluation of Integrals We have devoted particular care to the evaluation of integrals, in order to soften the influence of singularities, yet obtain a good accuracy and high computational speed. We use double precision arithmetic. We evaluate the elliptic integrals using library functions [24]. The inner integration in (7), along the source segment Cj, is performed numerically in the following way. Let us consider the source segment shown in Fig. 4. In the coordinate system Oxz, the endpoints of the segment are P1(x1, z1) and P2(x2, z2). A local coordinate system is attached to the segment, so that its origin (Ouv) is in the middle of the segment, the u-axis is along the segment, and the v-axis is perpendicular to it. Let us assume that the global coordinates of the field point are P(x, z). The local (u, v) coordinates of the field point are evaluated and the point P is projected onto the u-axis to obtain 'P . The minimal distance between P and the segment is calculated. Two distinct cases are considered: first, when P ' lies on the segment, and second, when it is out of the segment (either towards negative u-coordinates or towards positive u-coordinates). In the first case, the distance is equal to 'PP and the segment is divided into two integration intervals, bounded by P '. The integration is further carried on these two parts separately. In the second case, the minimal distance is the distance between P and the closer end of the segment (P1 or P2), and the integration is carried out on the whole segment as one integration interval. xO z u v x1 x2 z1 z2 Ouv P P' P2 P1 Fig. 4 Local coordinate system for evaluation of integrals Based on our experience, if the minimal distance is greater than one half of the length of the integration path ( 1 2 P P ), the integration is performed on the whole path as a unique integration interval, using a Gauss-Legendre integration formula. Otherwise, the integration interval is divided into nonuniform subintervals (at most 30), whose lengths progressively increase away from P '. Each increase is by the factor of 2. The same integration formula is used for all subintervals, both for the potential and for the field components. High-accuracy Quasistatic Model for Bodies of Revolution 149 The outer integration in (7), along the field segment Ci, is also performed numerically using Gauss-Legendre integration. If the electric field is evaluated using differentiation, numerical experiments have shown that the optimal choice for the evaluation of the electric field is to take ABn =D (Fig. 3), where  = 10−6. As the result, all the integrals (and their derivatives) are calculated to at least 5 significant digits. 3. BENCHMARK RESULTS The analysis method was tested on various examples where analytical solutions exist (Fig. 5). 3.1. Conducting Sphere Shown in Fig. 5a is a conducting sphere located in a vacuum. The radius of the sphere is a = 10 mm. The theoretical capacitance of the sphere is Cth = 4e0a = 1.112650 pF. The cross section of the sphere is a circle, which is approximated in our computations by a regular polygon with np sides. Hence, the generatrix of the sphere is a semi-circle, which is approximated by a polygonal line with ns = np / 2 uniform segments. In the numerical model, the actual sphere is approximated by a set of right conical frustums. In order to reduce the error of the geometrical modeling, we use the same strategy as in [25]: the radius of the given sphere (a) is the mean value of the radius of the circle inscribed into the polygon (rin) and the radius of the circle circumscribed around the polygon (rout). Hence, rout = 2a / (1 + cos (/np)) and the generatrix is easily constructed. 10 30 er 20 4 30 32 10 10 62 60 50 16 14 4 er (a) (b) (c) (d) (e) Fig. 5 Longitudinal cross sections of benchmark structures: (a) conducting sphere, (b) dielectric-covered conducting sphere, (c) conducting prolate ellipsoid, (d) spherical capacitor, and (e) coaxial-line section; all dimensions are in millimeters 150 A. DJORDJEVIĆ, D. OLĆAN, J. PETROVIĆ, N. OBRADOVIĆ, S. FILIPOVIĆ This strategy is in accordance with the theorem, due to Maxwell, that the capacitance of a conducting body is larger than the capacitance of an inscribed body and smaller than the capacitance of a circumscribed body [26]. The numerical result for the capacitance obtained with 20 pulses is Cnum = 1.112098 pF, which corresponds to a relative error with respect to Cth of around 0.0005. The relative error is reduced below 10−6 when the number of pulses is increased to 150. 3.2. Dielectric-Covered Conducting Sphere Fig. 5b shows a conducting sphere, whose radius is a = 10 mm, covered by a concentric dielectric layer. The outer radius of the dielectric is b = 30 mm and the relative permittivity is er = 10 4. The remaining space is a vacuum. The theoretical capacitance of the sphere is Cth = 4e0 / ((b − a) / erab + 1/b) = 3.336459 pF. The computed value, obtained with 20 segments per spherical surface, is Cnum = 3.336185 pF, so that the relative error is around 0.0005. The same low relative error is obtained for any other er ranging from 1.000000 to 1018. Similar results are obtained for a sphere with several concentric dielectric layers. 3.3. Conducting Prolate Ellipsoid Fig. 5c shows a conducting prolate ellipsoid, located in a vacuum. The longer semi-axis of the spheroid (which is the axis of rotational symmetry) is a = 10 mm and the shorter semi-axis is b = 2 mm. The theoretical capacitance is Cth = 8e0c / ln((a + c) / (a − c)) = 0.4755518 pF, where 2 2 c a b= − . In order to keep the relative error below 0.001, at least 30 segments are needed. 3.4. Spherical Capacitor A spherical capacitor, which consists of two concentric conducting spherical shells, is shown in Fig. 5d. The radius of the inner conductor is a = 10 mm, the inner radius of the outer conductor is b = 30 mm, and the outer radius of the outer conductor is c = 32 mm. The medium is a vacuum. The theoretical matrix of electrostatic induction coefficients for this system is           e+ − e − e − − e − − e = c ab ab ab ab ab ab ab ab 0 00 00 t h 4 44 44 ][B and the corresponding capacitance matrix is pF 3.5604801.668975 1.6689750 4 4 4 0 ][ 0 0 0 t h       =           e − e − e = c ab ab ab ab C . Using 20 segments per spherical surface (i.e., a total of 60 unknowns), the computed capacitance matrix is High-accuracy Quasistatic Model for Bodies of Revolution 151 pF 3.5587161.668181 1.66818110048.2 ][ 6 num         − = − C . If 50 segments are used, then pF 3.5601921.668842 1.66884210684.2 ][ 9 num         − = − C . Theoretically, c11 = 0 because the inner conductor is completely shielded by the outer conductor. The numerical result for c11 is very small, indicating a high accuracy of computations. 3.5. Coaxial Line The last example considered here, for which an analytical solution exists, is a section of a coaxial line (Fig. 5e), whose dielectric is Teflon, of relative permittivity er = 2.1. The radius of the inner conductor is a = 2 mm, the inner radius of the outer conductor is b = 7 mm, and the outer radius of the outer conductor is c = 8 mm. The coaxial line is open-circuited at both ends and the width of both gaps between the conductors is 5 mm. The length of the inner conductor is la = 50 mm, the inner length of the outer conductor is lb = 60 mm, and the outer length of the outer conductor is lc = 62 mm. The structure shown in Fig. 5e has significant fringing capacitances at both ends. In order to compute the per-unit-length capacitance of the coaxial line (C '), we have to remove the effect of the fringing capacitances. In the middle zone of the structure, which is sufficiently far away from the ends, the structure of the electric field is practically the same as in an infinitely long line. If we increase the length of the structure for Dl (i.e., if we increase la, lb, and lc for Dl), without changing the gap widths, the fringing capacitances will remain the same. Hence, the corresponding increase in the mutual capacitance between the inner and the outer conductor can be attributed only to the increased capacitance of the middle zone. Following this reasoning, we compute the mutual capacitance for the original dimensions of the structure ( )1( 1 2 c ) and for the increased length ( )2( 12 c ). From these two results, ( 2) (1) 12 12 ( ) /C c c l = − D . Using 35 segments for the inner conductor and 93 for the outer conductor, for Dl = 2 mm, the computed per-unit-length capacitance is num C = 93.24421 pF/m. The theoretical per-unit-length capacitance is th r 0 2 / ln( / )C b ae e = = 93.25647 pF/m. The relative error between numC and t hC is 0.00013. 3.6. Run Time The run time of the program is primarily influenced by the number of unknowns. The program is not parallelized, i.e., it uses only one core. With 100 unknowns, the run time is less than 1 s on a desktop computer with Intel Core i7-3770 @ 3.4 GHz, 32 GB RAM, and 64-bit Windows operating system. 4. MEASUREMENTS USING COAXIAL CHAMBER In this section we implement the technique for the BoR analysis, described in Section 2, to the coaxial chamber shown in Fig. 1. We describe the model of the chamber and the calibration procedure, and present some measurement results. 152 A. DJORDJEVIĆ, D. OLĆAN, J. PETROVIĆ, N. OBRADOVIĆ, S. FILIPOVIĆ 4.1. BoR Model of Measurement Setup The segmented model of the chamber looks as in Fig. 6a. The plot shows the generatrix. The numbers of segments were chosen by an educated guess and numerical experiments (i.e., convergence tests) so to provide a good accuracy at a reasonable run time. Obviously, there are several differences between the model and the actual structure shown in Fig. 1: the generatrix in Fig. 6a does not completely follow the contours of the actual device. When analyzing antennas and various microwave circuits, the structure must have ports and it is excited at those ports [11]. This is the same situation as in actual measurements. However, in our electrostatic model, the excitation is “virtual”: the conductors are assumed to be at a certain potential with respect to the reference point. No interconnections are provided between the conductors and the surroundings. In Fig. 1, which shows the actual device, the inner conductor of the chamber extends all the way to the VNA reference plane at the bottom of the SMA connector. In measurements, the inner and the outer conductors of the SMA connector further extend into the VNA connector. However, in the electrostatic model, the conductors must be floating. Hence, in Fig. 6a, the inner conductor of the coaxial line is left open-circuited inside the SMA connector. The outer conductor of the chamber in Fig. 1 has an opening at the mouth of the SMA connector (where it extends to the mating SMA connector of the VNA), whereas in Fig. 6a there is no such opening. The structure shown in Fig. 6a is completely shielded so that there is no electric field outside. Hence, the shape of the outer surface of the outer conductor is irrelevant. For simplicity, we have taken a spherical shape. Shifted reference plane Shifted reference plane 0.01 0.005 0 −0.005 −0.01 −0.015 −0.02 0.006 0.004 0.002 0 −0.002 −0.004 −0.006 0.005 0.01 0.0150 0.02 0.025 0.03 0.004 0.0060 0.008 0.01 0.012 0.0140.002 SMA connector Teflon Dielectric sample Plunger Plate (a) (b) Fig. 6 Segmented model of (a) chamber shown in Fig. 1 and (b) coaxial-line section and its positive image; red segments are for inner conductor, blue segments are for outer conductor, and green segments are for dielectric-to-dielectric interfaces; coordinates are in meters High-accuracy Quasistatic Model for Bodies of Revolution 153 In Subsection 4.2 we present numerical results of the electrostatic BoR analysis. In Subsection 4.3 we describe a theoretically rigorous calibration procedure that relates the actual setup with the electrostatic model. The aim of the calibration is to obtain a unique and measurable result for the chamber capacitance as seen looking upwards from the shifted reference plane in Fig. 6a. 4.2. Numerical Results for Empty Chamber We consider an empty chamber, without a sample, but with a gap of mm 2=h between the brass plate and the plunger. (Equivalently, the relative permittivity of the sample is 1.) The electrostatic analysis of the structure shown in Fig. 6a yields the following matrix of the electrostatic induction coefficients: pF 81495.6881959039973.83418344 13843.8341832985513.83418356 ][       − − =B . Note that the numerically obtained matrix ][B is almost perfectly symmetrical (up to 8 significant digits). The matrix of the partial capacitances is pF 67651.8540126139973.83418344 13843.834183294550.00000012 ][       =C . The outer surface of the chamber in Fig. 6a is a sphere, which we approximate in the same way as described in Subsection 3.1. The theoretical capacitance, Csphere = 1.8543089243 pF, agrees with the computed c22 within the first four digits. 4.3. Calibration Referring to the previous subsection, the mutual capacitance Cmodel = c21  3.8342 pF is the capacitance between the inner conductor of the chamber and the outer conductor. The modeled structure (Fig. 6a) includes the inner conductor of a section of the coaxial line (within the zone of the SMA connector) whose length is 3 mm. This conductor is left open-circuited, but it contributes to Cmodel. Hence, its influence must be calibrated-out. There exists a strong fringing effect at the open end of the coaxial line. This is a similar situation as described in Subsection 3.5. There also exists a discontinuity at the transition between the coaxial line and the chamber. Hence, we cannot assume that the field structure along the whole line is the same as in an infinitely long line. (Note that the field in an infinitely long line corresponds to the electric field of the guided TEM wave.) We consider this coaxial-line section and its positive image in the shifted reference plane (Fig. 6b). This structure has two identical fringing zones. The computed capacitance between the inner conductor and the outer conductor is Ccoax_double = 0.7016 pF. One half of it can be ascribed to the coaxial line in Fig. 6a, assuming that the TEM field exists all the way up to the shifted reference plane (although this is not true). Hence, the apparent capacitance of the chamber, looking from the shifted reference plane upwards, is Cchamber = Cmodel − Ccoax_double / 2. Note that, theoretically, we cannot uniquely define Cchamber because it depends on the presence of the inner coaxial-line conductor, which affects the fringing field in the vicinity of the shifted reference plane. However, the described procedure of evaluating the apparent capacitance is essentially the same as used in actual measurements, where 154 A. DJORDJEVIĆ, D. OLĆAN, J. PETROVIĆ, N. OBRADOVIĆ, S. FILIPOVIĆ we measure the input admittance at the VNA calibration plane, then shift the reference plane, and calculate the new admittance. In this procedure, it is assumed that a pure TEM wave exists in the coaxial line all the way up to the shifted reference plane. On the other hand, from manufacturer’s data, we know the geometrical dimensions of the SMA connector and that its dielectric is Teflon. Hence, the actual length of the coaxial line, from the VNA calibration plane in Fig. 1 up to the shifted reference plane, is lcoax = 11.75 mm. As in Subsection 3.5, the per-unit-length capacitance of the coaxial line is calculated to be C ' = 96.045 pF/m, so that the capacitance of this section (assuming that the electric field has the same structure as in an infinitely long line) is Ccoax = lcoax C ' = 1.1285 pF. The apparent capacitance transformed back from the shifted reference plane to the VNA calibration plane is thus Cat VNA reference plane = Cchamber + Ccoax = Cmodel + 0.7777 pF = 4.6064 pF. This is physically the same result as evaluated by measurements at the VNA calibration plane, looking towards the chamber. This capacitance, measured at f = 30MHz, is Cmeasured = (4.60  0.01) pF and it agrees with Cat VNA reference plane within the measurement uncertainty. 4.4. Examples of Measurement of Dielectric Parameters In this subsection, we implement the complete measurement setup (VNA, coaxial chamber, and software) to evaluate parameters of various dielectric samples. The general procedure is to measure the reflection coefficient of the chamber and the dielectric sample (at the VNA reference plane) and, hence, calculate the corresponding complex admittance. From the admittance, we evaluate the complex capacitance of the chamber. Thereafter, we use the quasistatic model of the chamber with the sample. In that model, we vary (optimize) the relative complex permittivity of the sample so to obtain the same complex capacitance as measured. The procedure can be simplified because, in a reasonably wide range of permittivities, the capacitance is practically a linear function of the permittivity (i.e., Cat VNA reference plane = er + , where α and  are constants). Hence, it is sufficient to implement a linear fit in the complex domain between two capacitances computed for two assumed permittivities, which are selected, e.g., based on an educated guess. If the sample is small (i.e., d and h are sufficiently smaller than the diameter of the plate shown in Fig. 1), the electric field in the whole dielectric sample is practically homogeneous. In that case, the measurement procedure is simple. First, the sample is inserted, fixed by the plunger, and the complex capacitance C is measured. Second, the sample is removed, the plunger is positioned at the same elevation h as when the sample was present, and the capacitance C0 of the empty chamber is measured. This is an elementary situation in electrostatics, for which C − C0 = ((er − 1) e0d 2 / 4h). Hence, er can easily be calculated. In order to illustrate the applications of the coaxial chamber shown in Fig. 1, we present here results for three measured samples. Two samples are printed-circuit board (PCB) substrates, measured at f = 100 MHz. The first substrate is Taconic 602-250. The measured relative permittivity was er = 2.55, which agrees well with the manufacturer’s data (er = 2.50). The measured loss tangent was tan  < 0.001 (below the resolution of our measurement system). For the second substrate, FR-4, we obtained er = 4.49 and tan  = 0.025, which agrees well with the data from [7]. High-accuracy Quasistatic Model for Bodies of Revolution 155 The third example is a ceramic material – alumina (Al2O3) doped with nickel oxide (NiO), mechanically activated by ball-milling for 60 minutes and sintered at 1400 °C [27]. Fig. 7 shows the relative permittivity and loss tangent of the material in the frequency range from 1 MHz to 500 MHz. The material is lossy and, hence, the relative permittivity significantly decreases with frequency. Mathematically, this decay follows from the causality conditions [7]. The measurement uncertainty at the lowest frequencies (around 1 MHz) is large because the input admittance of the chamber is very small (i.e., the chamber behaves almost like an open circuit). Hence, very small measurement errors of the reflection coefficient cause huge errors of the input admittance. The accuracy at frequencies in the range from 10 MHz to 100 MHz is much better. The accuracy at higher frequencies decreases because the field in the chamber cannot be considered to be quasistatic anymore. For these frequencies, the estimation of the relative permittivity requires a full-wave model of the chamber. Fig. 7 Measured relative permittivity and loss tangent of alumina doped with nickel oxide 5. CONCLUSION A high-precision and efficient quasistatic numerical method for the analysis of arbitrary metallo-dielectric bodies of revolution was presented. The method has been developed for measurements of dielectric parameters of small disk-shaped samples, for frequencies up to several hundred MHz. For higher frequencies, up to around 10 GHz, a full-wave (dynamic) solver is under development. 156 A. DJORDJEVIĆ, D. OLĆAN, J. PETROVIĆ, N. OBRADOVIĆ, S. FILIPOVIĆ Acknowledgement: This paper was funded in part by the Project F133 of the Serbian Academy of Sciences and Arts, and in part by the Ministry of Education, Science and Technological Development of the Republic of Serbia. REFERENCES [1] Basics of Measuring the Dielectric Properties of Materials, Application Note, Keysight Technologies, Document available at: https://www.keysight.com/zz/en/assets/7018-01284/application-notes/5989-2589.pdf. [2] O. V. Tereshchenko, F. J. K. Buesink and F. B. J. Leferink, "An overview of the techniques for measuring the dielectric properties of materials", In Proceedings of the XXXth URSI Gen. Ass. Sci. Symp., vol. 1320, Istanbul, Turkey, 2011, pp. 1–4. [3] T. P. Marsland and S. Evans, "Dielectric measurements with an open-ended coaxial probe", IEE Proc. Microw., Antennas Propag., vol. 134, no. 4, pp. 341–349, 1987. [4] B. Sanadiki and M. Mostafavi, "Inversion of inhomogeneous continuously varying dielectric profiles using open-ended waveguides", IEEE Trans. Antennas Propag., vol. 39, no. 2, pp. 158–163, Feb. 1991. [5] D. K. Ghodgaonkar and V. V. Varadan, "A free-space method for measurement of dielectric constants and loss tangents at microwave frequencies", IEEE Trans. Instrum. Meas., vol. 37, pp. 789–793, 1989. [6] A. M. Nicolson and G. F. Ross, "Measurement of the intrinsic properties of materials by time domain techniques", IEEE Trans. Instrum. Meas., vol. IM–19, pp. 377–382, 1970. [7] A. R. Djordjević, R. M. Biljić, V. D. Likar-Smiljanić and T. K. Sarkar, "Wideband frequency-domain characterization of FR-4 and time-domain causality", IEEE Trans. Electromagn. Compat., vol. 43, pp. 662–667, 2001. [8] P. Dankov, B. Hadjistamov, I. Arestova and V. Levcheva, "Measurement of dielectric anisotropy of microwave substrates by two-resonator method with different pairs of resonators", PIERS Online, vol 5, pp. 501–505, Oct. 2009. [9] A. Đorđević, J. Dinkić, M. Stevanović, D. Olćan, S. Filipović and N. Obradović, "Measurement of permittivity of solid and liquid dielectrics in coaxial chambers", Microw. Rev., vol. 22, pp. 3–9, Dec. 2016. [10] R. F. Harrington, Time-Harmonic Electromagnetic Fields, Hoboken, NJ: Wiley-IEEE Press, 2001, Chapter 1. [11] B. M. Kolundžija and A. R. Djordjević, Electromagnetic Modeling of Composite Metallic and Dielectric Structures, Norwood, MA: Artech House, 2002, pp. 6–8. [12] R. F. Harrington, Field Computation by Moment Methods, Hoboken, NJ: Wiley-IEEE Press; 1993. [13] M. Salazar-Palma, T. K. Sarkar, L.-E. Garcia-Castillo, T. Roy and A.R. Djordjević, Iterative and Self- Adaptive Finite-Elements in Electromagnetic Modeling, Norwood, MA: Artech House, 1998. [14] J. V. Surutka and D. M. Veličković: "Some improvements of the charge simulation method for computing electrostatic fields", Bull. Serb. Acad. Sci. Arts, Class Sci. Techn., no. 15, pp. 27–44, 1981. [15] D. M. Veličković and A. Milovanović, "Electrostatic field of cube electrodes", Serbian J. Electr. Eng., vol. 1, pp. 187–198, June 2004. [16] A. R. Djordjević, M. B. Baždar, R. F. Harrington and T. K. Sarkar, LINPAR for Windows: Matrix Parameters for Multiconductor Transmission Lines, Norwood, MA: Artech House, 1999. [17] CST Studio Suite, Available at: https://www.3ds.com/products-services/simulia/products/cst-studio-suite/. [18] M. M. Nikolić, A. R. Djordjević and M. M. Nikolić, ES3D: Electrostatic Field Solver for Multilayer Circuits, Norwood, MA: Artech House, 2007. [19] R. A. Handelsman and J. B. Keller, "The electrostatic field around a slender conducting body of revolution", SIAM J. Appl. Math., vol. 15, pp. 824–841, July 1967. [20] R. Barshinger, "The electrostatic field about a thin oblate dielectric body of revolution", SIAM J. Appl. Math., vol. 52, pp. 651–675, May 1991. [21] O. Ciftja, A. Babineaux and N. Hafeez N, "On the electrostatic potential of a uniformly charged ring", Eur. J. Phy., vol. 30, 623–627, May 2009. [22] A. R. Djordjević, Electromagnetics, Belgrade, Serbia: Academic Mind, 2008, Section 2.5. [23] D. Olćan, Diakoptic Analysis of Electromagnetic Systems, Ph.D. Thesis, School of Electrical Engineering, University of Belgrade, Serbia, 2008, pp. 26–28. [24] IMSL Fortran and C Application Development Tools, Houston, TX: Visual Numerics, 1997. [25] J. Dinkić, D. Olćan, A. Djordjević, А. Zajić, "Design and optimization of nonuniform helical antennas with linearly varying geometrical parameters", IEEE Access, vol. 7, pp. 1–12, Oct. 2019. [26] A. R. Djordjević, Fundamentals of Electrical Engineering, Belgrade, Serbia: Academic Mind, 2016, Section 1.10.1. [27] S. Filipović, N. Obradović, S. Marković, A. Đorđević, I. Balać, A. Dapčević, J. Rogan, V. Pavlović, "Physical properties of sintered alumina doped with different oxides", Sci. Sinter., vol. 50, pp. 1–11, 2018.