Planar multilayer structure analysis: an educational approach

Planar multilayer structure analysis: an educational approach

D.B. Ferreira^I; A.F. Tinoco Salazar^II; I. Bianchi^II; J.C. da S. Lacava^II

^IFundação Centro de Pesquisa e Desenvolvimento em Telecomunicações – CPqD Rodovia Campinas – Mogi-Mirim, km 118,5 13086-902 Campinas – SP, Brazil [e-mail]

^IILaboratório de Antenas e Propagação, Instituto Tecnológico de Aeronáutica, Pr. Mal. Eduardo Gomes 50, 12228-900 São José dos Campos – SP, Brazil, [e-mail], [e-mail], [e-mail], [e-mail]

ABSTRACT

This paper discusses how symbolic computation combined with a circuit model can be used for analyzing planar multilayer structures, in a manner suitable for educational approach. Working in the Fourier domain, expressions for the transversal spectral Green's functions are evaluated in compact, closed form using the symbolic computation capability of the Mathematica package. Printed antennas were analyzed through the method of moments. Further validation was achieved with the IE3D and HFSS packages.

Index Terms: Planar multilayer structures, spectral fields, Green's functions, circuit model, method of moments, education.

I. INTRODUCTION

Sources positioned on the interfaces of planar multilayer structures can be used as a practical model in a wide variety of applications such as geophysical investigations, remote sensing, optoelectronics, microwaves, military surveillance and antenna theory [1], [2]. For accurate results, the analysis of these structures requires numerical methods such as the Method of Moments (MoM), Finite Elements or Finite Differences. It is presently widely accepted that MoM-based algorithms are suitable for rigorous numerical analysis of printed structures of small to medium sizes (in terms of wavelength) stacked up in layers [1]. However, for the application of this method, either in the spectral or in the spatial domain, the corresponding Green's functions need be derived. Due to their complex calculations, an essential question is posed: which is the most proper approach for determining the entire set of Green's functions, particularly for educational purposes?

In the spatial domain, Green's functions for multilayer media are traditionally represented by Sommerfeld's integrals. Due to the oscillatory nature of these integrals, their numerical computation is inefficient and time-consuming [1]. Consequently, this approach is not recommended for a first graduate course on applied electromagnetism. On the other hand, spectral Green's functions can be derived in closed form if appropriate techniques are utilized [3]-[5]. However, the calculations of spectral fields in multilayer structures are usually tedious and error-prone if done by hand, what is especially true in the analysis of structures with anisotropic and bi-isotropic materials [6], [7].

To overcome this limitation, this paper presents a new elegant procedure for calculating the electromagnetic fields in multilayer structures. Using the full-wave equivalent circuit, first established by Dreher in [8], and the symbolic computation capability of the Mathematica package [9], expressions for the transformed electromagnetic fields are derived in a straightforward, error-free way. Consequently, the corresponding spectral Green's functions can be determined in compact, closed form, with considerable reduction of the calculation time. Based on this procedure, a method-of-moments (MoM) algorithm was implemented in Mathematica with the purpose of analyzing printed antennas. The results confirm the usefulness of the new procedure, which can also be applied to more complex structures and substrates, like in the anisotropic and bi-isotropic cases. Furthermore, the implemented algorithm can be utilized in a microstrip antenna course for engineering students, since it allows for the fast evaluation of electrical and geometrical effects on the antenna parameters.

II. THEORY

A planar multilayer structure composed of N + 2 isotropic, linear and homogenous layers stacked up in the z direction is shown in Fig. 1. The layers are assumed to be unbounded along the x and y directions.

The lower layer, having complex permittivity εg and complex permeability µg, is denoted the ground layer and occupies the negative-z region. The next N layers are characterized by thickness ℓ_n, complex permittivity ε_nand complex permeability µn, where 1 ≤ n ≤ N. The planar interface z = d_Nseparates the N-th layer from free space (the upper layer). Metallic patches are printed at arbitrary positions on each one of the N + 1 interfaces of the structure. The surface current densities on these elements are given by J_Sξ(x, y)= x J_Sξ x(x, y)+ y J_Sξ y(x, y)-(ξ∈ {g, n}), where boldface letters represent vectors, and x and y are the unit vectors along the x-and the y-directions, respectively. These current densities are considered the virtual sources of the electromagnetic field within each layer. Since the layers are free of sources (the sources are located on the interfaces of the multilayer structure, as depicted in Fig. 1), the electric field in the spatial domain of a monochromatic wave can be written, according to spectral techniques, as a superposition of plane waves, as follows [10]

where the function E (k_x,k_y,z) is the spectral electric field, and k_xand k_yare the spectral variables.

For the n-th confined layer, located between the interfaces z = _dn-1 and z = d_n, the spectral fields E_n(k_x,k_y,z) and H_n(k_x,k_y,z) can be written as the superposition of two plane waves traveling in the ± z directions, that is

where e_nτ (k_x,k_y) and h_nτ (k_x,k_y) are the amplitudes of the spectral fields, k_nzτ =(−1)^τ(ω² µ_nε_n− u )^1/2 and u² = + . Since the propagation function is given by for the n-th confined layer then, for a plane wave traveling in the +z direction, the square root in the propagation constant k_nz1 needs to be adequately evaluated; in this case, Re{k_nz1} ≤ 0 and Im{k_nz1} ≥ 0. Similarly, for the plane wave traveling in the –z direction, the square root in k_nz2 is calculated considering that Re{k_nz2} ≥ 0 and Im{k_nz2} ≤ 0. With these restrictions on the signs of the real and imaginary parts of k_nzτ, Sommerfeld's radiation condition [14] is satisfied.

From Maxwell's equations and using the symbolic computation capability of the Mathematica package, the amplitudes of the spectral fields in the x and y directions (e_nxτ , e_nyτ , h_nxτ , and h_nyτ , with τ = 1 or 2) can be expressed in terms of the field amplitudes in the z direction, e_nzτ and h_nzτ , according to [5]. They are listed below and the Mathematica command lines utilized for their determination are given in the Appendix.

In order to establish the full-wave equivalent circuit for the n-th layer, firstly the x-and y-component of the spectral fields, given by equations (2)-(3), are specified at the z = dn interface, producing the following set of four equations in the unknown amplitudes e_nzτ and h_nzτ

Solving this system (equations (8)–(11)), expressions for e_nzτ and h_nzτ, in terms of E_nx(k_x,k_y,d_n), E_ny(k_x,k_y,d_n), H_nx(k_x,k_y,d_n), and H_ny(k_x,k_y,d_n) are obtained. Then, replacing these expressions in the same set of transversal components given by equations (2)-(3), but now specified at z = d_n-1, results in the following relationship between the transversal components of the spectral fields (t ∈ {x, y}) at the upper and the lower interfaces of the n-th layer

where

with Again, symbolic computation was performed in the Mathematica package. As an example, the calculation of the first row of the matrix is given in the Appendix. Consequently, for the n-th confined layer, a two-port full-wave equivalent circuit, as shown in Fig. 2, can be established.

An equivalent procedure is employed to evaluate the fields at the free space and the ground layer interfaces. The boundary conditions are then applied for the planar interfaces, resulting in the circuit elements shown in Table I. Then, the multilayer structure depicted in Fig. 1 can be modeled by the circuit illustrated in Fig. 3.

where,

Since the propagation functions for free space and the ground were chosen to be and , respectively, Sommerfeld's radiation condition is satisfied if the square roots in the propagation constants and k_gz are calculated as follows: Re{K_0z} ≥ 0, Im{K_0z} ≤ 0, Re{K_gz} ≥ 0 and Im{K_gz} ≤ 0.

Before finishing this section, it is important to point out that the theory developed here can also be applied to planar multilayer structures where the metallic patches are positioned inside the confined layers. For these cases, the metallic patch can simply be considered to be at the interface between two adjacent layers having the same electrical characteristics. Hence, the developed theory can be utilized without any change.

III. APPLICATIONS

As a first application, a printed Yagi-Uda antenna fed by a delta-gap generator was investigated. In this case, the multilayer structure, shown in Fig. 4, is composed of the free space–substrate–free space layers.

The radiator consists of two narrow flat dipoles printed on the dielectric substrate of thickness d₁, electric permittivity ε₁ and magnetic permeability µ₀. The length and width of the dipoles are denoted respectively by L₁ and 2w₁ (L₁ >> 2w₁), for the driven element, and L₂ and 2w₂(L₂ >> 2w₂), for the parasitic one. The distance between elements is a. The transversal Green's functions at the antenna interface z = d₁ are easily computed from the circuit model depicted in Fig. 5, as follows

where

η₀ is the intrinsic impedance of free space and ε_r1 is the relative electric permittivity of the substrate.

Notice that the transversal Green's functions in (19) to (22) are equivalent to the ones presented in [16], where they were evaluated analytically, i.e., without employing a full-wave equivalent circuit to represent the planar multilayer structure.

Once the Green's functions are derived, the next step is to determine the Yagi antenna's electric characteristics, such as its input impedance and radiation pattern. This involves solving integral equations constrained by the pertinent boundary conditions. In the present case, the following equations are obtained by enforcing the tangential electric field component being zero on the perfect conducting surface of the antenna elements,

where E_f (x, y) denotes the tangential electric field component created by the antenna feeder. Since L_v>> 2w_v, (v ∈ {1, 2}), the unknown consists of just the x-components of the surface current densities J_xv(x, y) excited on the dipoles.

Using the MoM procedure -the most widely utilized numerical technique for rigorous analysis of printed geometries on multilayer planar media - to solve the system of equations (25)-(26), the electric current densities J_xv(x, y) are expanded in a set of basis-functions as follows,

where J_{x 1 m} and J_x2nare the basis-functions, and I_x1m and I_x2n are coefficients to be determined.

After applying Galerkin method [11] and Parseval's theorem [12], the following linear system is established

where

J_xv(k_x, k_y) is the Fourier transform of J_xv(x, y), S1 is the driven element surface, p = 1, 2, … , M, and q = 1, 2, … , N.

Piecewise-linear sub-domain rooftop basis-functions (taking into account the edge condition) were used for expanding the surface current densities J_xv(x, y), so that

where

As a delta-gap generator is used to feed the driven element, then

where δ (x) is the Dirac's delta function. Consequently, the excitation equation (30) can be rewritten as

with M being an odd number.

Using this approach, a Mathematica-based CAD software capable of performing the analysis of printed antennas was implemented. For greater efficiency in the numerical calculation of the elements Z_pm, …, Z_qnof the impedance matrix [Z] (equation (29)), mathematical procedures -like parity analysis of the Green's functions, change of the coordinate system (rectangular to polar) and asymptotic extraction technique [13] -were applied. To illustrate their effect, after changing the coordinate system -(ρ, ϕ) are the polar coordinates -and undertaking the parity analysis, equation (30) assumes the simpler form [14]

For applying the asymptotic extraction technique, first the asymptotic term is calculated, resulting in

Consequently, Z_pmcan be evaluated from

The first double integral in (42) depends on the operating frequency, whereas the second one does not. These are the well-known Sommerfeld's integrals, whose integrands exhibit singularities in the form of branch points and poles, so that their numerical calculation requires careful attention. The poles, often complex, correspond to surface and leaky waves that can be excited in the dielectric layer. The number of poles depends on the thickness of the layer, its electrical permittivity and the wave number. An efficient way to calculate these integrals is the use of a deformed path of integration, with adaptive routines in the vicinity of the poles and branch points. Details on the numerical integration can be found in [15] – [17]. In this work, a triangular path was chosen to calculate the elements Z_pm, …, Z_qnof the impedance matrix. Finally, as the antenna is fed directly by a delta-gap generator, it is important to point out that the dielectric thickness poses no limitation on the application of the present theory.

Since M is an odd number, the antenna input impedance can be promptly calculated using the relation

where E₀ was selected as 1V.

To validate the developed CAD, a Yagi-Uda antenna with L₁=L₂=56.294mm, 2w₁=2w₂=3.0mm, a=28.174mm,d₁=3.048mm,ε_r1= 2.55and tgδ= 0.0022 (ε1=εr1ε0(1-itgδ)),was analyzed. Fig.6 shows the results obtained for the input impedance (M=N =17) compared to those simulated in two commercial software, HFSSandIE3D(aMoM-basedsoftware).

The radiation pattern was calculated using the followin gexpressions, derived from the stationary phase method [14]:

where J_xs (k_xe,k_ye) = J_x1 (k_xe,k_ye) + J_x2 (k_xe,k_ye) and the subscript e denotes the functions are calculated at the stationary phase point (kxe = k0 sinθ cosϕ; kye = k0sin θsinϕ), for 0 ≤ θ ≤ π/2 and 0 ≤ϕ< 2π.

Fig. 7 presents the results for theyz-planeradiation pattern calculated at 2.0 GHz.

The agreement between our results for the input impedance and radiation pattern with those simulated in HFSS and IE3D certifies the accuracy of the technique presented in this work.

As a second application, the three-layer structure (perfect ground plane – substrate – free space) shown in Fig. 8 was analyzed.

The radiator in question consists of a flat dipole printed atop the dielectric substrate of thickness d1, electric permittivity ε1 and magnetic permeability µ₀, such that 2w₁ denotes the dipole width and L₁ its length (L₁ >> 2w₁). The transversal Green's functions at the dipole interface z = d₁ are easily computed from the circuit model depicted in Fig. 9, as follows [14]

where

Once again, the transversal Green's functions evaluated through the full-wave equivalent circuit ((47) to (50)) are equivalent to the ones analytically calculated in [16], [17].

Graphics for the input impedance and radiation pattern of a flat dipole (L₁ = 53.134 mm and 2w₁ = 3.0 mm, printed on a substrate of ε_r1 = 2.55, d₁ = 3.048 mm, tg δ = 0.0022) are shown in Figs. 10 and 11, respectively.

Again, the results for the dipole input impedance and radiation pattern compared with simulations performed in HFSS and IE3D validate the accuracy of our CAD.

IV. FINAL COMMENTS

In this paper an elegant and efficient approach, based on a full-wave circuit model and employing the symbolic capability of the Mathematica package, to calculate the spectral Green's functions of planar multilayer structures was presented. Using this procedure and the method of moments, a Mathematicabased CAD software was implemented for the analysis of printed dipoles and Yagi-Uda antennas. Results obtained for the input impedance and the radiation pattern are in good agreement with the simulations performed in HFSS and IE3D, thus validating the proposed technique. It is noteworthy that the use of Mathematica accelerates both the Green's functions calculations and the program coding time, making this procedure suitable for a first graduate course on applied electromagnetism.

Received 26 Ago. 2011; for review 2 Sept. 2011; accepted 9 March 2012

APPENDIX

^Appendix

• [1] M. I. Aksun and G. Dural, "Clarification of issues on the closed-form Green's functions in stratified media," IEEE Transactions on Antennas and Propagation, vol. 53, no 11, pp. 3644-3653, November 2005.
• [2] F. Mesa, R. Marqués, and M. Horno, "On the computation of the complete spectral Green's dyadic for layered bianisotropic structures", IEEE Transactions on Microwave Theory and Techniques, vol. 46, no 8, pp. 1158-1164, August 1998.
• [3] T. Itoh, "Spectral domain immittance approach for dispersion characteristics of generalized printed transmission lines", IEEE Transactions on Microwave Theory and Techniques, vol. 28, no 7, pp. 733-736, July 1980.
• [4] N. K. Das and D. M. Pozar, "A generalized spectral-domain Green's function for multilayered dielectric substrate with applications to multilayer transmission lines," IEEE Transactions on Microwave Theory and Techniques, vol. 35, no 3, pp. 326335, March 1987.
• [5] J. C. S. Lacava and L. Cividanes, "A new technique for microstrip antenna analysis," 3Brazilian Microwave Symposium, Natal, pp. 258-266, July 1988 (In Portuguese).
• [6] J. C. S. Lacava, A. V. Proaño De la Torre, and L. Cividanes, "A dynamic model for printed apertures in anisotropic stripline structures," IEEE Transactions on Microwave Theory and Techniques, vol. 50, no 1, pp. 22-26, January 2002.
• [7] F. Lumini, "Analysis of chiral multilayer structures in spectral domain," Ph.D. Thesis, Instituto Tecnológico de Aeronáutica, São José dos Campos, Brazil, 2000 (In Portuguese).
• [8] A. Dreher, "A new approach to dyadic Green's function in spectral domain," IEEE Transactions on Antennas and Propagation, vol. 43, no 11, pp. 1297-1302, November 1995.
• [9] Mathematica, Wolfram Research. <http://www.wolfram.com/products/ mathematica/>
• [10] C. A. Balanis, Antenna Theory: Analysis and Design, 3rd ed., New York: John Wiley, 2005.
• [11] J. J. H. Wang, Generalized Moment Methods in Electromagnetics, New York: John Wiley, 1991.
• [12] R. F. Harrington, Time-Harmonic Electromagnetic Fields, New York: McGraw Hill, 1961.
• [13] L. Cividanes, "Analysis of microstrip antennas in dielectric uniaxial multilayer structures," Ph.D. Thesis, Instituto Tecnológico de Aeronáutica, São José dos Campos, Brazil, 1992 (In Portuguese).
• [14] D. B. Ferreira, "Spherical microstrip antennas," MsC. Thesis, Instituto Tecnológico de Aeronáutica, São José dos Campos, Brazil, 2011 (In Portuguese).
• [15] J. R. Mosig, Integral equation technique, Chapter 3, pp. 133-211, In: Numerical techniques for microwave and millimeterwave passive structures, T. Itoh (Ed.), New York: John Wiley, 1989.
• [16] S. J. S. Sant'Anna, J. C. S. Lacava, and D. Fernandes, From Maxwell's equations to polarimetric SAR images: a simulation approach. Sensors 2008; 8(11):7380-7409.
• [17] R. Garg, P. Bhartia, I. Bahl, and A. Ittipiboon, Microstrip antenna design handbook, Chapter 3 (Full-wave analysis of microstrip antennas), Norwood, Artech House, 2001.

Appendix

Publication Dates

• Publication in this collection
  16 July 2012
• Date of issue
  June 2012

History

• Received
  26 Aug 2011
• Accepted
  09 Mar 2012
• Reviewed
  02 Sept 2011