A Simulation-Assisted Non-destructive Approach for Permittivity Measurement Using an Open-Ended Microwave Waveguide

A new convenient and non-destructive permittivity measurement method is presented. No physical cut of specimens is needed here for material characterisation. In the setup, the material under test is placed in the near-field region of a microwave open-ended waveguide. An electromagnetic model of the setup is built in the Computer Simulation Technology simulation software. Employing optimisation, the permittivity is obtained from the measured reflection coefficients S11\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {S}_{11}$$\end{document}. Using the same technique, the effect of the model size is investigated that could reduce the modelling effort for large structures. The efficiency of a traditional method (i.e., Newton) and an intelligent algorithm (i.e. particle swarm optimisation) for permittivity calculation is thoroughly studied and compared. The proposed methodology is validated by experimental data. It is demonstrated that the proposed method can provide more accurate permittivity results than the intrusive in-waveguide measurement. The proposed methodology can contribute to electromagnetic analysis, thickness measurement and non-destructive evaluation.


Introduction
Microwaves are commonly used for telecommunications and food processing, while extensive attention has been given to their potential in material characterisation (e.g., porosity evaluation [1,2] and moisture measurement [3]) and non-destructive testing (e.g., corrosion detection [4], crack detection [5][6][7], thickness variation [8], delamination detection [9,10] and impact damage inspection [11][12][13][14]). Microwaves propagate in air and dielectric materials with low attenuation. The electric permittivity is an intrinsic parameter of a material that describes the interaction with the electromagnetic field, and microwaves are highly sensitive to its variation. With the knowledge of the permittivity, the quality or condition of the material can be readily assessed. A number of microwave techniques have been used for permittivity measurement: resonance methods, transmission line technique, free space methods and open-ended rectangular waveguide/coaxial probe techniques. The resonance approach is inherently narrowband, and it requires careful sample preparation (spheres or cylinders required) and calibration [15]. For the transmission line technique, samples with the same inner dimensions of the waveguide need to be cut and perfectly fit within the waveguide. And the resolution of loss tangent measurements provided by this method is limited (typically ± 0.01) [16]. In the free space method, the sample and the sample holder should be placed between two horn antennae, and special attention must be paid to the sample geometry and location. Hence, in these three methods, cutting and machining of test samples is required. Considering the time for contacting the mechanical workshop, task scheduling, actual machining work and delivery, the total time could be a couple of days. This kind of mechanical work is labour intensive and sample cutting from the examined component may not be permissible due to coatings [17] and could undermine structural integrity. Therefore, the open-ended method is well-suited for industrial non-destructive measurements. Compared with the coaxial open-ended probe, the open-ended rectangular waveguide is more capable for measurements of low-permittivity materials and anisotropic materials [16]. And the field region around the waveguide is larger than that around the coaxial probe, where electromagnetic fringe fields closely surround the tip of the probe [18].
For the open-ended waveguide method, permittivity calculation requires rigorous mathematical formulation of the electromagnetic field outside the waveguide flange. As the exact solution of Maxwell's equations for the field is not available, approximate models have been reported [19][20][21]. In these models, a 2D problem is considered where the material space is commonly assumed infinite, and the analysis is limited to layered structures. Hence, the finite element method can be employed to provide more accurate description of the experimental arrangements and the electromagnetic field.
An appropriate optimisation method is also needed for permittivity extraction from one-port measurement data. The permittivity value used in the finite element modelling should to be optimised in order to produce the simulated reflection coefficients S 11 approximate the measured data. The optimisation techniques can be classified into two groups: iterative methods (e.g., Gauss-Newton, quasi-Newton, gradient and conjugate gradient) and heuristic methods (e.g., genetic algorithm (GA), particle swarm optimisation (PSO) and ant colony optimisation (ACO)). For the iterative methods, the gradient of the objective function usually needs to be calculated at each iteration, and the accuracy of the results depends on the selection of the initial estimation. However, for the heuristic algorithms, the gradient of the problem being optimised is not required. In addition, they have better global search abilities than traditional iterative methods.
In this paper, the non-destructive measurement of complex permittivity at microwave frequencies using an open-ended rectangular waveguide is presented. The experimental setup is modelled in the Computer Simulation Technology (CST) simulation software [22], and the input permittivity value is optimised until a converged solution is found. The calculation procedure is demonstrated by an example of a thin ceramic coated plate. The S 11 obtained from the forward simulation is used as the input for the inverse calculation. The iterative quasi-Newton method and the heuristic PSO method are employed, and the accuracy and efficiency of both methods are compared. The effect of the model size is investigated for large workpieces. Sample thickness is also evaluated using the same methodology. An experiment is conducted to measure the permittivity of the dielectric layer of a printed circuit board (PCB) with one-sided copper cladding. The experimental results obtained are compared with the manufacturer's data and results given by the microwave transmission line technique.

Definition of Permittivity
The complex permittivity ε can be described as: where ε 0 = 8.8542×10 −12 F/m is the permittivity of free space, and ε r is the relative permittivity. The real part ε r , or dielectric constant, characterises the ability of a material to store the electric field energy. The imaginary part ε r , or dielectric loss factor, reflects the ability of the material to dissipate the energy in the form of heat. ε r is positive due to energy conservation. In practice, the electrical behaviour is described by dielectric constant and loss tangent (tanδ= ε r /ε r ).

Non-destructive Measurement of Permittivity
The schematic diagram of the proposed setup for complex permittivity measurement is presented in Fig. 1. The signal response of the material is measured with a Vector Network Analyser (VNA). The material test sample can be of an arbitrary shape, with or without a metal backplate and of a multi-layered dielectric construction with a layer of unknown permittivity or thickness. In the test, the sample should be placed within the near-field region of the waveguide, at a standoff distance d of [23]: where a is the broad inner dimension of the waveguide, and λ is the wavelength of the incident electromagnetic field. For a frequency range, the minimum frequency is used for the estimation of the allowable standoff distance, so that the requirement can be met for all the frequencies used.
The flowchart that describes the permittivity calculation procedure is shown in Fig. 2. The reflection coefficients S 11 measured are imported to MATLAB ® software. A 3D model of the measurement setup is built in the CST simulation tool Fig. 2 Flowchart of the non-destructive permittivity calculation procedure using CST simulation and optimisation using the parametric modelling technique. The MATLAB-CST interaction technique reported in [9,10] is adopted here. MATLAB is utilised to control the simulation process via visual basic for application (VBA), which is the programming language in the CST environment. MATLAB calls the VBA programme to modify the permittivity, update the model and perform the electromagnetic calculation. The error function used for the evaluation of the solution is defined as: where N f is the number of the frequency sampling points. S 11,meas and S 11,C ST are the measured and simulated S 11 , respectively. |S 11 | and S 11 are the magnitude and phase of S 11 in decibels (dB) and degrees (deg), respectively. If the convergence is not satisfied, a new solution is generated using an appropriate optimisation algorithm. In this work, two different methods are employed and discussed: the quasi-Newton method and the particle swarm optimisation.

(1) Quasi-Newton method
The quasi-Newton method is iterative, involving a series of line searches [24]. At each iteration, the function value and the derivative are computed. The main procedure is as follows: (i) Set the iteration counter k = 1 and an estimated permittivity x 1 = [ε r tan δ]. Initialise a matrix H 1 = I (2×2 identity matrix). Calculate the derivative of F err (x k ), ∇ F err (x k ), using finite-difference approximation.
Here the bisection method is used for the 1D search of α. The initial points of the 1D search are carefully selected to guarantee that all the solutions to be produced are physically acceptable (ε r ≥ 1 and tanδ ≥ 0 for the present case). Repeat until the local maximum iteration k max is met.
Here the Broyden, Fletcher, Goldfarb and Shanno (BFGS) update expression of U k [25] is employed.
Step ii with k = k+1. Repeat until the convergence criteria or the maximum iteration k max is met.

(2) Particle swarm optimisation (PSO)
The PSO is a population-based algorithm, which was inspired by the social and cognitive behaviour of animals like fish schooling and bird flocking, adapting to the environment to find food sources [25]. A population (called a swarm) of candidate solutions (called particles) is used in the algorithm. The movements of these particles are guided by their own best-known positions as well as the knowledge gained by the swarm. For the permittivity calculation, the position of each particle with ε r and tanδ is in a 2D search space. The fundamental steps of the PSO are addressed: Step 1 Initialisation: an initial particle swarm is randomly distributed within the search space and random "velocities" are assigned.
Step 2 Evaluation: the location of each particle of the swarm (e.g., x k i =[ε r tanδ], the location of Particle i at Iteration k) is evaluated using a fitness function, which is defined as a measure of the quality of the solutions produced at each iteration: where ψ is a small value to avoid singularity. This objective function needs to be maximised, and the limit is 1/ψ.
Step 3 Velocity vector updating: the velocity vector used to update the current position of each particle is calculated by [26]: where w is the inertia weight factor. v k i and v k+1 i are the velocities of Particle i at Iterations kand k+1, respectively. r 1 and r 2 are random numbers between 0 and 1 used to maintain the diversity of the swarm. p k i and p k g are the best position of Particle i (personal best) and the global best position in the swarm within k iterations, respectively. Constant c 1 represents the cognitive factor that pulls the particle to its own best position. Constant c 2 is the 'social' factor that pushes the swarm to converge to the current global best position.
The inertia weight factor w k is adjusted dynamically throughout the optimisation process [27]: where w min and w max are the minimum and maximum inertia weight factors, respectively. k max is the maximum iteration number.
Step 4 Position updating: the position of Particle i at Iteration k+1 is updated by: where v k+1 i is the corresponding velocity vector, and tis the time step value that is set to 1 in the present work.
Step 5 Check convergence: return to Step 2 with k = k+1. Repeat until the maximum iteration k max is met.

Forward Calculation
Here a CST model is built and the values of the reflection coefficient calculated are used as an input for the subsequent inverse calculation of permittivity. As shown in Fig. 3, the model is made up of a waveguide, a dielectric material layer and a metal plate. The X-band microwave frequency range (8)(9)(10)(11)(12) with 401 sampling points is used. The dielectric is a 96% glazed alumina ceramic, ε r = 7.2 and tanδ =  Table 1.
For the given waveguide size and operating frequency range, the allowable maximum standoff distance d is 3.4 cm. In the present case, d is set to zero for simplicity. The inside dimensions of the waveguide are denoted by a × b. The ceramic and aluminium layers are of the same planar dimensions, which are set to ξ a × ξ b. Here the size factor ξ is defined as the ratio of the model size to the inner dimensions of the waveguide. In the present case, a size factor of 10 is employed. The time domain solver is used to calculate the magnitude and phase of S 11 , which are shown in Fig. 4.

Inverse Calculation
Main parameters of the quasi-Newton method and PSO used for the inverse calculation are listed in Tables 2 and 3, respectively. A high-performance computer (HPC) with an Intel Xeon Central Processing Unit (CPU) and 96 GB memory were used. The simulation results provided by the two optimisation methods are presented in Table 4. The errors of ε r given by both methods are within 0.1 %, while the error of tanδ predicted by both methods is almost the same (3.75 %). The history of the objective functions, F err and F PSO , is presented in Fig. 5. It is seen that the quasi-Newton method converges rapidly requiring only 10 iterations, while for PSO a sharp increase of the fitness function value occurs at the end of the optimisation process, 50 iterations. The maximum fitness function value is F PSO = 12.15, which corresponds to  the error function value F err = 0.07 according to Equation (5). However, the minimum error function value obtained by the quasi-Newton method is 0.11. The PSO method does not converge as expected, more accurate results being obtained through the global search. As shown in Fig. 6, in the quasi-Newton case ε r starts to stabilise after 11 iterations, and tanδ levels off after 13 iterations. This indicates that it is slightly faster to find an approximate value for ε r than tanδ. Compared with the Newton method, the solutions by the PSO method are closer to the real values in the first 10 iterations. In the PSO case ε r starts to stabilise within 20 iterations, while the tanδ curve varies considerably in the first 40 iterations. The computation time of the PSO was nearly 55.58 h, while it took approximately 33.35 h to complete the quasi-Newton code.

Model Size Effect
From the simulation point of view, for a small test piece it is straightforward to create an exact model. However, for a large structure, it is practical to use a relatively small model to speed up the computation without affecting the accuracy. Hence, in order to investigate the effect of the model size, the responses in four other cases with the same permittivity but varied size factors (ξ = 1, 2, 5, 15 and 20) are studied. As shown in Fig. 7, the signal converges with increasing ξ , and little difference can be found between the cases with ξ ≥10. The curve for ξ = 2 is slightly different from the other curves, which could be caused by electromagnetic resonance.
The size of the sample modelled at ξ = 2 (i.e., 45.72 mm × 20.32 mm) is close to the wavelengths of the incident waves (37.5-25 mm for 8-12 GHz). The computation time of a single simulation run for ξ = 10, 15 and 20 is approximately 2.2, 4 and 6.5 min, respectively. This suggests the CST model with ξ = 10 can describe the electromagnetic behaviour of a much larger size structure, hence reducing the calculation time.
It should be noted that the size factor is associated with the inner dimensions of the waveguide, which also relates to the available operating frequency range. With a smaller wavelength (i.e., a higher frequency), a smaller mesh grid should be adopted, which would increase the computation time but improve the computational accuracy. Hence, there is a compromise between the frequency and the efficiency.

Thickness Calculation
With the known permittivity, the thickness of the sample layer can be calculated using the CST based methodology. Here 1D Newton-Raphson method is used, since there is only one unknown parameter. The update equation for the thickness parameter x is: The process is repeated until the absolute value of the derivative F err is less than a predefined tolerance. The initial estimation of the thickness for the Newton-Raphson method  Tables 2 and 3, respectively. The accurate thickness of 2 mm is obtained by both methods. The histories of the objective function values and the thickness variable in both cases are presented in Figs. 8 and 9, respectively. In Fig. 8b, the fitness function value, F P SO , increases as would be expected, since a better solution is generated as it iterates. Similar to the permittivity calculation, the convergence of the Newton-Raphson method is quicker than that of the PSO, and in both cases the solution approaches the actual value after the second iteration. However, to achieve the exact thickness value of 2 mm, the computation time of the PSO programme was roughly 35.10 h for 30 iterations, while for the Newton-Raphson code it took around 56 min to complete 10 iterations. In this example, the Newton-Raphson method is faster and preferred for estimating the sample thickness.

Test and Permittivity Calculation for a Rogers PCB Laminate
The permittivity of the dielectric layer in a Rogers ® RO4350B PCB laminate with one-sided copper cladding was evaluated. The dielectric was made of a glass woven fabric-reinforced thermoset composite filled with ceramic particles. The length and width of the panel were 304.8 mm and 228.6 mm, respectively. The thickness of the composite and that of the copper backplate was 1.52 mm and 35 µm, respectively. The permittivity of the composite is available in the datasheet provided by the PCB laminate supplier. It was mea-  Variations of the calculated sample thickness with respect to the iteration using the two optimisation methods sured using the IPC-TM-650 2.5.5.5 standard [29], which is widely adopted in industry. This is a narrowband mea- surement, as the permittivity is given at a specific resonance frequency, here of 10 GHz [30].
Two microwave permittivity measurement methods were adopted here: the transmission line broadband technique and the proposed methodology. Like in the IPC test, the material anisotropy due to the woven fibre architecture with respect to the incident electromagnetic field is not considered here but can be accounted for [31].
The setup for the transmission line technique is schematically illustrated in Fig. 10. The VNA was calibrated before the test using the thru-reflect-line (TRL) standard [32]. Calibration accuracy was checked using a flush short (a 6.24 mm thick copper sheet). Over the operating frequency range, the attenuation was less than 0.1 dB and the phase shift was 180 • ± 0.5 • . After the calibration, the S-parameter measurement was carried out at the indicated reference planes, Fig. 10. A personal computer was connected to the VNA by an IEEE-488 cable. A MATLAB programme was developed for data acquisition and permittivity computation. A sample with the inner dimensions of the rectangular waveguide was cut and mounted on the waveguide flange. The transmission coefficients (S 21 ) were retrieved for permittivity calculation, since this is generally more accurate than using the reflection coefficients (S 11 ) [33]. The Nicolson-Ross-Weir algorithm [34] was used for the permittivity calculation. The calculated dielectric constants and loss tangents over 8-12 GHz are presented in Fig. 11. The real part of the permittivity remains relatively stable over the frequency range. Considering the effect of the calibration errors, the uncertainty for the dielectric constant is less than ± 0.5 %, while for the loss tangent the uncertainty is approximately ± 40 %.
For the methodology proposed in this work, a single port calibration of the HP8510C VNA used was carried out before the measurement. The measured S 11 is shown in Fig. 12. The quasi-Newton method fails to find a reasonable solution, while the PSO method successfully obtains its best solution with a fitness value of 100 at the 8th iteration.  The permittivity comparison between the manufacturer data [35] and the measured results is demonstrated in Table 5. There is good agreement in the prediction of the dielectric constant. Due to the low-loss characteristic of the material, the loss tangent provided by the transmission line technique is nearly double the manufacturer value. However, the simulation-assisted method proposed offers a more accurate estimation of the loss tangent (within 6 %).

Concluding Remarks
A novel non-destructive approach for complex permittivity determination with the use of optimisation and electromagnetic software has been proposed. Using an open-ended waveguide facilitates the measurement process without the need for sample machining, which is relatively operator friendly. In combination with electromagnetic simulation and optimisation, the proposed methodology is readily applied to obtain the electric permittivity ε r from the measured reflection coefficient S 11 data. Two representative optimisation methods (i.e., Newton's method and PSO) have been used. From the simulation results of a ceramic coated plate, it is seen that both methods can offer satisfactory results. The convergence of the Newton method is better than that of the PSO. The PSO can provide approximate solutions without the careful selection of initial guess points required by the Newton method, while its computation process is relatively time-consuming, several hours. Therefore, when choosing a suitable optimisation method, there is a compromise between the accuracy and computation time.
For relatively large test pieces, in order to reduce the computational time, the effect of the CST model size on the microwave signal has been studied. For the same material, little variation of the magnitude and phase is observed when the model is 10 times larger than the inner dimensions of the open-ended rectangular waveguide. In addition, it has been demonstrated that the thickness can also be evaluated using the proposed method when the permittivity is known.
Measurement on a dielectric layer, a woven glass fibrereinforced ceramic filled thermoset material, with a copper cladding was carried out. It has been revealed that the proposed method with the use of PSO can provide more accurate permittivity results than the intrusive in-waveguide measurement. The unsuccessful implementation of the quasi-Newton method is probably due to the computational complexity induced by the two sharp resonance dips, which are observed in the magnitude of S 11 in Fig. 12.
The computational time is associated with the implementation of the optimisation algorithm, model size, settings of the electromagnetic simulation (e.g., mesh size and solver parameters), performance of the computer used and the material properties to be evaluated. The primary reason for the long computational time demonstrated here is the low loss tangent value of the material used in the example (in the order of 10 −3 ). It is known that it is more difficult to obtain an accurate measurement of the loss tangent than that of the real part of the permittivity. However, if the results obtained after a few iterations are sufficient for applications, the computational time can be greatly reduced. For the thickness calculation, it is seen that the code can be executed within 1 h with an acceptable value. Future work is required to improve the computational efficiency, e.g., use of an efficient optimisation method and a high-speed computer.