FACTA UNIVERSITATIS Series: Electronics and Energetics Vol. 33, N o 1, March 2020, pp. 73-82 https://doi.org/10.2298/FUEE2001073K © 2020 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND NUMERICAL COMPACT MODELING APPROACH OF DISPERSIVE MAGNETOELECTRIC MEDIA BASED ON SCATTERING PARAMETERS  Miloš Kostić 1 , Nebojša Dončov 2 , Zoran Stanković 2 , John Paul 3 1 Innovation Center of Advanced Technologies, Niš, Serbia 2 Faculty of Electronic Engineering, University of Niš, Serbia 3 Electromagnetics scientist, Nottingham, United Kingdom Abstract: Z-TLM based compact modeling approach for dispersive media exhibiting magnetoelectric coupling is presented in this paper. Scattering parameters based representation of considered medium is created in a form of compact model by extracting effective electromagnetic parameters using a retrieval method, and implementing them into a non-uniform TLM grid. Proposed approach is illustrated here on the example of dispersive isotropic chiral medium modeling. Key words: Dispersive media, Compact models, Scattering parameters, Z-TLM method, Retrieval method, Non-uniform mesh. 1. INTRODUCTION Numerical modeling techniques nowadays represent important tools in a research process of complex materials especially when it is not possible or worthwhile to solve a problem with analytical approach. Two most used discrete time domain numerical techniques, the Finite Difference Time Domain (FD-TD) method [1] and Transmission Line Matrix (TLM) method [2] are very suitable for solving problems of electromagnetic (EM) wave propagation through complex structures and media. Even though the FD-TD method is often favored by researchers, the TLM method offers in some cases a more straightforward approach for describing and modeling different discontinuities, internal boundaries, propagation in dispersive media etc. This is a result of TLM feature that both electric and magnetic field components are solved in center of the TLM cell simultaneously without a need for temporal and spatial averaging. Modifying and extending TLM method with Z transformation techniques create valuable means for an efficient time domain modeling of linear isotropic and anisotropic, bi-isotropic, nonlinear, quantum, chiral materials and metamaterials [3-7]. This so-called Received March 21, 2019; received in revised form June 18, 2019 Corresponding author: Miloš Kostić Innovation Center of Advanced Technologies, 18000 Niš, Serbia (E-mail: r.i.p.romeo@gmail.com)  74 M. KOSTIĆ, N. DONĈOV, Z. STANKOVIĆ, J. PAUL Z-TLM method supports implementation of the Debye, Drude, Lorentz and other dispersion models along with specific methods which allow for describing materials with complex frequency dependencies. Compact models allow for complex structure, artificial or multilayered material to be represented as one effective material block via scattering parameters which to some degree simplifies numerical analysis and modeling process. In addition, compact models can be also used to reduce computational and time costs of the simulation by using a much coarser mesh for modeling of thin material panel, where instead of direct modeling by a fine mesh the material is replaced with single interface between two TLM cells [8,9]. In this paper, a formulation based on nonlinear constitutive relations and discretization of Maxwell’s equations which allows implementation of most general properties of dispersive and anisotropic materials into the Z-TLM non-linear grid, is described. Procedures for applying TLM method in modeling of dispersive and general anisotropic media inside of non-uniform mesh are given in [9-11]. Z-TLM based approach presented in [12,13] is here expanded to allow modeling of dispersive materials with magnetoelectric coupling characteristics while preserving the advantage of including the arbitrary frequency dependencies of modeled material EM parameters, i.e. these dependences do not have to necessarily follow some of the known dispersion models. Effective EM parameters which are used to characterize materials with magnetoelectric coupling are extracted from S parameters through retrieval procedure [14-16], approximated through the vector fitting (VF) method [17-19] and then used to form a compact model after applying the bilinear Z transforms, which is later included into the TLM scattering algorithm. Created model efficiently describes studied material based on provided S parameters and enables analyzing and observing EM field propagation throughout the medium. Proposed approach is demonstrated by modeling dispersive isotropic chiral material slab exhibiting magnetoelectric coupling [20]. 2. FORMULATION OF MAXWELL’S EQUATIONS FOR NON-UNIFORM MESH Non-uniform mesh generally enables a proper resource handling especially while modeling materials and structures with nonlinear characteristics. It allows for a usage of smaller size cells in areas which are physically small but have a greater EM importance and also bigger cells in less complex or less EM important areas of the structure which improves an overall efficiency of the simulations. In a non-uniform TLM cell (see Fig. 1) one or more directional space steps (x, z and z) are not the equal. Relations between incident V i and reflected V r voltage wave components of the twelve-port non uniform hybrid TLM cell [2] can be presented as: 1 1 2 2 2 1 3 3 10 4 4 9 5 5 6 6 6 5 r i z z r i z z r i x y r i y x r i z y r i z y V V I Z V V V I Z V V V I Z V V V I Z V V V I Z V V V I Z V                   , 7 7 8 8 8 7 9 9 4 10 10 3 11 11 12 12 12 11 r i z x r i z x r i y x r i x y r i y z r i y z V V I Z V V V I Z V V V I Z V V V I Z V V V I Z V V V I Z V                   (1) Dispersive Media Numerical Compact Modeling Method Based on Scattering Parameters 75 where impedance of each link line can be calculated based on link line inductance, time and space steps: / , k k Z L i t   (1, 2,...,12)k  (2) Fig. 1 Non-uniform TLM cell Using notations for fields, current and flux densities Maxwell’s curl equation can be written as , ef mf H J D t E J B t           (3) After given constitutive relations for electric and magnetic current and flux densities , , , ef mf J J D B , (3) can be written as: 0 0 0 0 ( / ), ( / ) ef e e r mf m r m H J E E E c H t t E J H H c E H t t                                     (4) where sign  denotes time domain convolution, ,e m  and ,e m  are electric and magnetic conductivity and susceptibility matrices respectively, 0, 0 are free space permittivity and permeability, respectively, and ,r r  are dimensionless matrices describing magnetoelectric coupling coefficients given as 1 r rr r r r r rr xyxx xz yx yy yz r zyzx zz c                        , 1 r rr r r r r rr xyxx xz yx yy yz r zyzx zz c                        (5) Formulation (4) is extended further with transformations and introduced additional compact notations to create normalized form of Maxwell’s equations (6): 76 M. KOSTIĆ, N. DONĈOV, Z. STANKOVIĆ, J. PAUL 1 1 2 2 2 ( ) , 2 2 2 ( ) t ef b e e r t mf b m m r C i i V V g V A V i T T T C V V i i r i A i V T T T                                                            (6) with 1 , , , b A C   representing background susceptibility matrix, matrix of inverse cell areas, normalized curl matrix and matrix of inverse cell length, respectively [9]. Previous form can be further transformed into the traveling wave format by using following relation where vector i V represents sum of appropriate incident wave intensities [9]: 2 2 4 , 2 2 4 TLM i TLM i C i V V V T C V i V i T              (7) This leads to: 1 1 2 4 2 ( ) , 2 4 2 ( ) i t ef e b e r i t mf m b m r V i V g V V A V i T V V i r i i A i V T                                                    (8) If reflected voltage wave r V and reflected free current r i are introduced as in [9], and matrix P is defined: 2 , 2 r i ef r i mf V V i i V V      (9) 1tP A       (10) while 1 b P   , (8) can be written as: 2 4 2 ( ( )), 2 4 2 ( ( )) r e b e r r m b m r V V g V V P V i T i i r i i P i V t                                (11) 3. MATERIAL MODELING BASED ON EXTRACTED EFFECTIVE PARAMETERS Modeling process based on proposed approach (see Fig. 2) begins by applying appropriate retrieval method [14-16] on the scattering matrix parameters of considered material, obtained either analytically, numerically or experimentally, in order to acquire effective permittivity, permeability and magnetoelectric coupling coefficients in the frequency range of interest. Dispersive Media Numerical Compact Modeling Method Based on Scattering Parameters 77 In order to calculate effective parameters from S parameters the index of refraction n and characteristic impedance of considered medium zred need to be first determined [15]. Based on calculated n and zred effective parameters of the material can be determined as: red n z   , rednz  , n   (12) After straightforward conversion of effective permittivity, permeability and magnetoelectric coupling coefficients to appropriate effective susceptibilities, (electric–e, magnetic–m and magnetoelectric–r which is related to (5)), they are then further approximated using the VF method [17-19] in order to represent each susceptibility in the form of rational function: [ , , ] 1 [ , , ] [ , , ] [ , , ]0 ( ) e m r e m r i e m r e m r pii C s s s       NP (13) In (13) NP[e,m,r] stands for the number of poles, s[e,m,r]pi represents the set of complex pole frequencies, and C[e,m,r]i are the pole residues of VF approximated susceptibilities. Next, the bilinear Z-transform is applied: 1 1 2 1 1 Z z s t z             (14) to obtain a discrete-time model: [ , , ] [ , , ] [ , , ] 0 [ , , ] [ , , ] 0 ( ) e m r e m r i e m r i i e m r i e m r i i B z z A z         NP NP (15) where A[e,m,r]i and B[e,m,r]i are real coefficients and z is time-shift operator. Fig. 2 Z-TLM compact modelling approach A and B coefficients are used to define compact model which is incorporated into TLM scattering procedure where electric and magnetic fields are calculated as: 78 M. KOSTIĆ, N. DONĈOV, Z. STANKOVIĆ, J. PAUL 1 1 1 1 (2 ) ( 2 ) 2 2 r e e r m m re e e rm rm m m re T V z S V i T i z S S V k V S S S i k i S S                                       (16) where 1 0 0 0 (4 4 4 ) e e e r T g        , 1 0 0 0 (4 4 4 ) m m m r T r        , 0 1 (4 4 ) e e e k g     , 0 1 (4 4 ) m m m k r     and 1 1 , , e m re rm S S S S represent additional material accumulator vectors. For simplicity, it is assumed that electric conductivity and magnetic resistivity are not dependent on frequency. 4. NUMERICAL RESULTS Proposed approach is illustrated here for an efficient modeling of dispersive isotropic d = 200 mm wide chiral material slab placed between two isotropic free-spaces [20]. In this case due to isotropic nature of the material slab, magnetoelectric susceptibilities have the form 0 0 1 0 0 0 0 r r r r c               , r r   (17) By using the retrieval method, effective susceptibilities are first extracted from S parameters analytically obtained in [20] (Figs. 3 - 6, marked with triangle). The retrieved electric, magnetic and magnetoelectric susceptibilities are shown in Figs. 7, 8 and 9, respectively (marked with red triangle). Fig. 3 Magnitude of scattering parameters S11 and S21 Dispersive Media Numerical Compact Modeling Method Based on Scattering Parameters 79 Fig. 4 Phase of scattering parameters S11 and S21 Fig. 5 Magnitude of scattering parameters S12 and S22 Fig. 6 Phase of scattering parameters S12 and S22 80 M. KOSTIĆ, N. DONĈOV, Z. STANKOVIĆ, J. PAUL Retrieved susceptibilities are approximated with the VF method [17-19] with fourth order rational function (NP = 4) (12). Obtained A’s and B’s coefficients for electric and magnetic susceptibility as well as magnetoelectric susceptibility are presented in Table 1. Accuracy of approximation is confirmed through perfectly matched comparison of retrieved parameters values (red triangles on the graph) and values calculated based on A and B coefficients (solid line) (see Figs. 7-9). Discrete time-models described by (14) are further incorporated into the Z-TLM scattering algorithm in (15) taking into consideration magnetoelectric coupling characteristic of isotropic chiral slab. Total of 800 TLM cells are used in x direction, while chiral slab itself is modeled with 200 cells and simulation is executed within 2000 time steps. Model was excited with initial z polarized Gaussian pulse which propagates along +x direction. Fig. 7 Real and imaginary parts of effective electric susceptibility of isotropic chiral medium Fig. 8 Real and imaginary parts of effective magnetic susceptibility of isotropic chiral medium Dispersive Media Numerical Compact Modeling Method Based on Scattering Parameters 81 Fig. 9 Real and imaginary parts of effective magnetoelectric susceptibility of isotropic chiral medium Two simulations are performed, first where the chiral medium was not considered present inside of the mesh in order to obtain incident field values, and second simulation with chiral medium included in the mesh in order to obtain the total field at the first interface, free space-chiral slab, and transmitted fields at the second interface, chiral slab- free space. Reflected field was calculated by deducting incident field from total field values. Accuracy of the approach is confirmed through a comparison of magnitudes and phases of scattering parameters from [20] and simulated scattering parameters shown with solid lines in Figs. 3 - 6. 4. CONCLUSION In this paper, an approach for compact modeling of dispersive media exhibiting magnetoelectric coupling, described with effective parameters extracted from scattering matrix, is presented. Since there is a wide variety of dispersive media exhibiting magnetoelectric coupling the compact models presented here provide efficient tools for characterization and inside view of EM waves propagation through these media. In future research this approach will be used in the design process of devices based on these media because it allows optimization and fine tuning of their characteristics in the frequency range of interest. Table 1 Coefficients of discrete-time model in (15) Coefficient e m r = r A 1 1 1 -3.994898e+00 -1.971846e+00 -3.994897e+00 5.985618e+00 -2.760262e-02 5.985617e+00 -3.986539e+00 1.971845e+00 -3.986539e+00 9.958199e-01 -9.723764e-01 9.958197e-01 B 5.471056e-05 -4.988997-01 5.224473e-03 2.651185e-08 9.860555e-01 -1.044641e-02 -1.093681e-04 1.174398e-02 -3.898947e-09 2.649748e-08 -9.860555e-01 1.044641e-02 5.471051e-05 4.871557e-01 -5.224469e-03 82 M. KOSTIĆ, N. DONĈOV, Z. STANKOVIĆ, J. PAUL Acknowledgement: This work has been supported by the Ministry of Education, Science and Technological Development of Serbia, project number TR32024. REFERENCES [1] K.S. Kunz and R.J. Luebbers, The Finite Difference Time Domain Method for Electromagnetics, CRC Press, 1993 [2] C. Christopoulos, The Transmission-Line Modelling (TLM) Method, IEEE/OUP Press, 1995. [3] J. Paul, C. Christopoulos and D.W.P. Thomas, "Generalized material models in TLM – Part I: Materials with frequency-dependent properties", IEEE Trans. Antennas and Propagation, vol. 47, no. 10, pp. 1528- 1534, 1999. [4] J. Paul, C. Christopoulos and D.W.P. Thomas, "Generalized material models in TLM - Part 2: Materials with anisotropic properties", IEEE Trans. Antennas and Propagation, vol. 47, no. 10, pp. 1535-1542, 1999. [5] J. Paul, C. Christopoulos and D.W.P. Thomas, "Time-domain modelling of electromagnetic wave propagation in complex materials", Electromagnetics, vol. 19, no. 6, pp. 527-546, 1999. [6] J. Paul, C. Christopoulos and D.W.P. Thomas, "Generalized material models in TLM – Part III: Materials with nonlinear properties", IEEE Trans. Antennas and Propagation, vol. 50, no. 7, pp. 997-1004, 2002. [7] J. Paul, C. Christopoulos, and D.W.P. Thomas, "Time-domain simulation of electromagnetic wave propagation in two-level dielectrics", International Journal of Numerical Modelling: Electronic Networks, Devices and Fields, vol. 22, no. 2, pp. 129-141, 2009. [8] N. Donĉov, M. Kostić, Z. Stanković, "Compact numerical models for efficient representation of EM field propagation through dispersive and anisotropic media", In Proceedings of the 5th International Conference IcETRAN 2018, Palić, Serbia, 2018, pp. 582-587. [9] M. Kostić, N. Donĉov, Z. Stanković, J. Paul, "Efficient TLM-based approach for compact modeling of anisotropic materials and composites", Applied Computational Electromagnetics Society (ACES) Journal, vol. 34, no. 1, pp. 1-10, 2019. [10] P. Saguet, H. Louzani and F. Ndagijimana, "The use of Z-transform in TLM with non-uniform meshes", In Proceedings of the Workshop on Computational Electromagnetics in Time-domain, CEM-TD, Atlanta, USA, 2005, pp. 68-71. [11] A. L. Farhat, S. Le Maguer, P. Queffelec and M. Ney, "TLM Extension to Electromagnetic Field Analysis of Anisotropic and Dispersive Media: A Unified Field Equation", IEEE Transactions on Microwave Theory and Techniques, vol. 60, no. 8, pp. 2339-2351, 2012. [12] T. Asenov, M. Kostić, N. Donĉov and B. Milovanović, “Z-TLM Method Simulation of Left-Handed Metamaterials Based on Retrieved Effective Parameters”, In Proceedings of the 2nd International Conference IcETRAN 2015, Silver Lake, Serbia, pp. MTI1.7.1-5, 2015. [13] M. Kostić, N. Donĉov, Z. Stanković and T. Asenov, "3D Z-TLM Modeling of Dispersive Lossy Metamaterial Structures Described by Scattering Parameters", Proceedings of the 3nd International Conference on Electrical, Electronic and Computing Engineering, IcETRAN 2016, Zlatibor, Serbia, pp. MTI2.7.1-4, 2016. [14] F. J. Hsieh and W. C. Wang, "Full extraction methods to retrieve effective refractive index and parameters of a anisotropic metamaterial based on material dispersion models", Journal of Applied Physics, vol. 112, 064907, September 2012. [15] D. R. Smith, D. C. Vier, Th. Koschny and C. M. Soukoulis, "Electromagnetic parameter retrieval from inhomogeneous metamaterials", Physical review E, 71, 036617, 2005. [16] V. Milosević, B. Jokanović and R. Bojanić, "Effective Electromagnetic Parameters of Metamaterial Transmission Line Loaded With Asymmetric Unit Cells", IEEE Transactions on Microwave Theory and Techniques, vol. 61, no. 8, pp. 2761-2772, Aug. 2013. [17] B. Gustavsen and A. Semlyen, "Rational approximation of frequency domain responses by vector fitting", IEEE Transactions on Power Delivery, vol. 14, no. 3, pp. 1052-1061, 1999. [18] B. Gustavsen, "Improving the pole relocating properties of vector fitting", IEEE Transactions on Power Delivery, vol. 21, no. 3, pp. 1587-1592, 2006. [19] D. Deschrijver, M. Mrozowski, T. Dhaene, and D. De Zutter, "Macromodeling of multiport systems using a fast implementation of the vector fitting method", IEEE Microwave and Wireless Components Letters, vol. 18, no. 6, pp. 383-385, 2008. [20] J. Paul, C. Christopoulos and D. W. P. Thomas (1999) “Time-Domain Modeling of Electromagnetic Wave Propagation in Complex Materials”, Electromagnetics, vol. 19, no. 6, pp. 527-546, 1999.