Equations of TLs illuminated by incident EM wave

As shown in Fig. 1, the coupling of the external EM wave to a TL is modeled using two additional induced sources in the per-unit-length equivalent circuit model of the line. Therefore, by considering the quasi-TEM approximation, the time-domain equations of a TL excited by the EM wave illumination can be written as follows [Reference Agrawal, Price and Gurbaxani12]:

(1)$$\displaystyle{{\partial V( {x, \;t} ) } \over {\partial x}} + l\displaystyle{{\partial I( {x, \;t} ) } \over {\partial t}} + rI( {x, \;t} ) = V_E( {x, \;t} ) , \;$$

(2)$$\displaystyle{{\partial I( {x, \;t} ) } \over {\partial x}} + c\displaystyle{{\partial V( {x, \;t} ) } \over {\partial t}} + gV( {x, \;t} ) = I_E( x, \;t) , \;$$

where l, r, c, and g denote the per-unit-length parameters for the TL. Also, the per-unit-length induced sources due to the EM excitation field (V _E and I _E) can be obtained from the incident magnetic ($\vec{B}_{inc}$) and electric ($\vec{E}_{inc}^t$) fields as follows:

(3)$$V_E( {x, \;t} ) = \displaystyle{\partial \over {\partial t}}\int_0^d {{\vec{B}}_{inc}\cdot {\vec{a}}_ndl} $$

(4)$$I_E( {x, \;t} ) = {-}g\int_0^d {{\vec{E}}_{inc}.\overrightarrow {dl} -c\displaystyle{\partial \over {\partial t}}} \int_0^d {{\vec{E}}_{inc}\cdot \overrightarrow {dl} } , \;$$

where $\vec{a}_n$ is the unit vector normal to the plane created by the x-axis and the connecting line between the centers of the two conductors of the TL (the x–y plane as shown in Fig. 1). For the uniform plane wave excitation, the electric field can be represented by:

(5)$$\vec{E}_{inc}( {x, \;y, \;z, \;t} ) = \quad\quad\quad\vec{P}f_W\left({t + \sqrt {\mu \varepsilon } ( {\sin \theta_E\sin \varphi_Ex+ \cos \theta_Ey}+ \sin \theta_E\cos \varphi_Ez} ) \right), $$

where θ _E and φ_E are the incidence angles, and function f _W denotes the waveform of the incident wave. According to Fig. 1, the polarization vector ($\vec{P}$) of the incident field can be written as

(6)$$\vec{P} = \left[{\matrix{ {P_x} \cr {P_y} \cr {P_z} \cr } } \right] = \left[{\matrix{ {-\sin \phi \cos \theta_E\sin \varphi_E + \cos \phi \cos \varphi_E} \cr {\sin \phi \sin \theta_E} \cr {-\sin \phi \cos \theta_E\cos \varphi_E-\cos \phi \sin \varphi_E} \cr } } \right], \;$$

where ϕ denotes the polarization angle of the electric field. Furthermore, in the case of a TL with a conductor at height d placed above a ground plane and illuminated by the uniform plane wave, the induced sources can be formulated as follows:

(7)$$V_E( {x, \;t} ) = {-}2d\sqrt {\mu \varepsilon } ( \cos \theta _EP_x + \sin \theta _E\sin \varphi _EP_y) \displaystyle{\partial \over {\partial t}}f_W\left({t + \sqrt {\mu \varepsilon } \sin \theta_E\sin \varphi_Ex} \right), \;$$

(8)$$\eqalign{I_E( {x, \;t} ) &= {-}2dP_y\left[{gf_W\left({t + \sqrt {\mu \varepsilon } \sin \theta_E\sin \varphi_Ex} \right)}\right. \cr &\left. {+ c\displaystyle{\partial \over {\partial t}}f_W\left({t + \sqrt {\mu \varepsilon } \sin \theta_E\sin \varphi_Ex} \right)} \right].}$$

It should be noted that these equations are derived by considering the incident filed and reflected field from the ground plane [Reference Paul5].

Fig. 1.

Modeling of a TL with a conductor at height d placed above a ground plane, which is illuminated by the external incident uniform EM wave with the propagation vector of $\vec{k}$, incidence angles of θ _E and φ_E, and the polarization angle of ϕ. The per-unit-length inductance, resistance, capacitance, and conductance of the line are denoted by l, r, c, and g, respectively. Also, V _E and I _E are the per-unit-length induced sources due to the EM excitation field, and V(x, t) and I(x, t) denote the line voltage and current, respectively, as a function of position (x) and time (t).

Generalized FDTD solution method for TLs with linear termination

In this subsection, for simplicity, the development of the proposed generalized FDTD approach is presented for a two-conductor TL as shown in Fig. 2. In the proposed method, the forward difference approximation of the spatial derivative is first used to express (1) in discrete space as follows:

(9)$$V_{i + 1}-V_i + r\Delta xI_i + l\Delta x\displaystyle{{\partial I( {x, \;t} ) } \over {\partial t}} = \Delta xV_E i , \;\quad \;i = 1\colon N_x$$

where V _i ≡ V((i − 1)Δx, t) and I _i ≡ I((i − 0.5)Δx, t). Moreover, it is assumed that the line is discretized into N _x cells of length Δx. Similarly, using the backward difference approximation of the spatial derivative, equation (2) can be discretized as follows:

(10)$$I_i-I_{i-1} + g\Delta xV_i + c\Delta x\displaystyle{{\partial V( {x, \;t} ) } \over {\partial t}} = \Delta xI_E i , \;\;\quad i = 2\colon N_x.$$

Fig. 2.

Illumination of the incident EM wave to a TL terminated with the linear load.

According to Fig. 2, by using the current–voltage relationships in the source (I _S = −V ₁/R _S + V _S/R _S) and load ($I_L = V_{N_x + 1}/R_L-V_L/R_L$), and replacing Δx ⇒ 0.5Δx in equation (10), the terminal relations (i = 0, i =  N _x + 1) can be obtained as follows:

(11)$$I_1 + ( {1/R_S + 0.5g\Delta x} ) V_1 + 0.5c\Delta x\displaystyle{{\partial V( {x, \;t} ) } \over {\partial t}} = 0.5\Delta xI_E 1 + V_S/R_S, \;$$

(12)$$-I_{N_x} + ( {1/R_L + 0.5g\Delta x} ) V_{N_x + 1} + 0.5c\Delta x\displaystyle{{\partial V( {x, \;t} ) } \over {\partial t}} = 0.5\Delta xI_E {N_x + 1} + V_L/R_L.$$

By defining the vectors

(13)$$\bar{{\bf V}}\equiv \left[{\matrix{ {V_1} & \cdots & {V_{N_x}} & {V_{N_x + 1}} \cr } } \right]^T_{( {N_x + 1} ) \times 1} , \;$$

(14)$$\bar{{\bf I}}\equiv \left[{\matrix{ {I_1} & \cdots & {I_{N_x-1}} & {I_{N_x}} \cr } } \right]^T_{{\rm N}_x \times 1} , \;$$

(15)$$\overline {{\bf V}_{\bf S}} \equiv \Delta x.\left[{\matrix{ {V_E 1 } & \cdots & {V_E {N_x} } & {V_E {N_x + 1} } \cr } } \right]^T_{{\rm N}_x \times 1} , \;$$

equation (9) can be written in the matrix form as follows:

(16)$${ \bar{\bf M}\bar{V}} + { \bar{\bf R}\bar{I}} + \bar{{\bf L}}\displaystyle{\partial \over {\partial t}}\bar{{\bf I}} = \overline {{\bf V}_{\bf S}} , \;$$

wherein:

(17)$$\bar{{\bf {\rm M}}}\equiv \left[{\matrix{ {-1} & 1 & 0 & \cdots & 0 \cr 0 & {-1} & 1 & \cdots & 0 \cr 0 & 0 & \ddots & \ddots & \vdots \cr \vdots & \ddots & \ddots & 1 & 0 \cr 0 & \cdots & 0 & {-1} & 1 \cr } } \right]_{N_x \times ( {N_x + 1} ) }, \;$$

(18)$$\bar{{\bf R}}\equiv \left[{\matrix{ {r\Delta x} & 0 & 0 & 0 & 0 \cr 0 & {r\Delta x} & \ddots & \vdots & 0 \cr 0 & 0 & \ddots & 0 & \vdots \cr \vdots & \vdots & \ddots & {r\Delta x} & 0 \cr 0 & 0 & \cdots & 0 & {r\Delta x} \cr } } \right]_{N_x \times N_x}, \;$$

(19)$$\bar{{\bf L}}\equiv \left[{\matrix{ {l\Delta x} & 0 & 0 & 0 & 0 \cr 0 & {l\Delta x} & \ddots & \vdots & 0 \cr 0 & 0 & \ddots & 0 & \vdots \cr \vdots & \vdots & \ddots & {l\Delta x} & 0 \cr 0 & 0 & \cdots & 0 & {l\Delta x} \cr } } \right]_{N_x \times N_x}.$$

Equations (10)–(12) can be similarly expressed in the matrix form as follows:

(20)$${ \bar{\bf N}\bar{I}} + { \bar{\bf G}\bar{V}} + \bar{{\bf C}}\displaystyle{\partial \over {\partial {\rm t}}}\bar{{\bf V}} = \overline {{\bf I}_{\bf S}} , \;$$

wherein:

(21)$$\overline {{\bf I}_{\bf S}} \equiv \left[{\matrix{ {0.5\Delta xI_E 1 + V_S/R_S} \cr {\Delta xI_E 2 } \cr \vdots \cr {\Delta xI_E {N_x} } \cr {0.5\Delta xI_F {N_x + 1} + V_L/R_L} \cr } } \right]_{( {N_x + 1} ) \times 1}, \;$$

(22)$$\bar{{\bf N}}\equiv \left[{\matrix{ 1 & 0 & 0 & \cdots & 0 \cr {-1} & 1 & 0 & \cdots & 0 \cr 0 & {-1} & 1 & \ddots & \vdots \cr \vdots & \ddots & \ddots & 1 & 0 \cr 0 & \cdots & \cdots & 0 & {-1} \cr } } \right]_{( {{\boldsymbol N}_{\boldsymbol x} + 1} ) \times {\boldsymbol N}_{\boldsymbol x}}, \;$$

(23)$$\bar{{\bf C}}\equiv \left[{\matrix{ {0.5c\Delta x} & 0 & 0 & 0 & 0 \cr 0 & {c\Delta x} & \ddots & \vdots & 0 \cr 0 & 0 & \ddots & 0 & \vdots \cr \vdots & \vdots & \ddots & {c\Delta x} & 0 \cr 0 & 0 & \cdots & 0 & {0.5c\Delta x} \cr } } \right]_{( {N_x + 1} ) \times ( {N_x + 1} ) }, \;$$

(24)$$\bar{{\bf G}}\equiv \left[{\matrix{ {1/R_S + 0.5g\Delta x} & 0 & 0 & 0 & 0 \cr 0 & {g\Delta x} & \ddots & \vdots & 0 \cr 0 & 0 & \ddots & 0 & \vdots \cr \vdots & \vdots & \ddots & {g\Delta x} & 0 \cr 0 & 0 & \cdots & 0 & {1/R_L + 0.5g\Delta x} \cr } } \right]_{( {N_x + 1} ) \times ( {N_x + 1} ) }.$$

The source matrices $\overline {{\bf V}_{\bf S}}$ and $\overline {{\bf I}_{\bf S}}$ include the line excitation sources and field-to-line coupled sources. Applying the forward difference approximation in the time domain to equations (16) and (20) results in:

(25)$${ \bar{\bf M}\bar{V}} + { \bar{\bf R}\bar{I}} + \Delta t^{{-}1}\bar{{\bf L}}( {{\bar{{\bf I}}}^{{\boldsymbol n} + 3/2}-{\bar{{\bf I}}}^{{\boldsymbol n} + 1/2}} ) = \overline {{\bf V}_{\bf S}} , \;$$

(26)$${ \bar{\bf N}\bar{I}} + { \bar{\bf G}\bar{V}} + \Delta t^{{-}1}\bar{{\bf C}}( {{\bar{{\bf V}}}^{{\boldsymbol n} + 1}-{\bar{{\bf V}}}^{\boldsymbol n}} ) = \overline {{\bf I}_{\bf S}} , \;$$

where superscript n denotes the time step. In the proposed method, voltage and current vectors are defined as the weighted average over the time step as follows:

(27)$$\bar{{\bf V}} = \alpha \bar{{\bf V}}^{{\boldsymbol n} + 1} + ( {1-\alpha } ) \bar{{\bf V}}^{\boldsymbol n}, \;$$

(28)$$\bar{{\bf I}} = \alpha \bar{{\bf I}}^{{\boldsymbol n} + 3/2} + ( {1-\alpha } ) \bar{{\bf I}}^{{\boldsymbol n} + 1/2}, \;$$

(29)$$\overline {{\bf V}_{\boldsymbol S}} = \alpha \overline {{\bf V}_{\boldsymbol S}} ^{{\boldsymbol n} + 1} + ( {1-\alpha } ) \overline {{\bf V}_{\boldsymbol S}} ^{\boldsymbol n}, \;$$

(30)$$\overline {{\boldsymbol I}_{\boldsymbol S}} = \alpha \overline {{\boldsymbol I}_{\boldsymbol S}} ^{{\boldsymbol n} + 1} + ( {1-\alpha } ) \overline {{\boldsymbol I}_{\boldsymbol S}} ^{\boldsymbol n}, \;$$

where α is the weighting factor (0 ≤ α ≤ 1). It can be observed that the proposed generalized FDTD method with specific values of α leads to the conventional FDTD approaches such as the explicit scheme (α = 0), the Crank–Nicholson scheme (α = 0.5), and the fully implicit scheme (α = 1). Finally, substituting the average vectors (equations (27)–(30)) into equations (25) and (26) gives the final iteration formulas based on proposed generalized FDTD for the voltages and currents of the line shown in Fig. 2, in the matrix form as follows:

(31)$$\left[{\matrix{ {\bar{{\bf V}}} \cr {\bar{{\bf I}}} \cr } } \right]^{n + 1} = \bar{{\bf \Phi }}\left[{\matrix{ {\bar{{\bf V}}} \cr {\bar{{\bf I}}} \cr } } \right]^n + \bar{{\bf \Psi }}\left({{\boldsymbol \alpha }{\left[{\matrix{ {\overline {{\boldsymbol I}_{\boldsymbol S}} } \cr {\overline {{\bf V}_{\boldsymbol S}} } \cr } } \right]}^{n + 1} + ( {1-\alpha } ) {\left[{\matrix{ {\overline {{\boldsymbol I}_{\boldsymbol S}} } \cr {\overline {{\bf V}_{\boldsymbol S}} } \cr } } \right]}^n} \right), \;$$

wherein:

(32)$$\bar{{\bf \Psi }} = \left[{\matrix{ {\Delta t^{{-}1}\bar{{\bf C}} + \alpha \bar{{\bf G}}} & {\alpha \bar{{\bf N}}} \cr {\alpha \bar{{\bf M}}} & {\Delta t^{{-}1}\bar{{\bf L}} + \alpha \bar{{\bf R}}} \cr } } \right]^{{-}1}, \;$$

(33)$$\eqalign{\bar{{\bf \Phi }} = \left[{\matrix{ {\Delta t^{{-}1}\bar{{\bf C}} + \alpha \bar{{\bf G}}} & {\alpha \bar{{\bf N}}} \cr {\alpha \bar{{\bf M}}} & {\Delta t^{{-}1}\bar{{\bf L}} + \alpha \bar{{\bf R}}} \cr } } \right]^{{-}1}\cr\left[{\matrix{ {\Delta t^{{-}1}\bar{{\bf C}}-( {1-\alpha } ) \bar{{\bf G}}} & {-( {1-\alpha } ) \bar{{\bf N}}} \cr {-( {1-\alpha } ) \bar{{\bf M}}} & {\Delta t^{{-}1}\bar{{\bf L}}-( {1-\alpha } ) \bar{{\bf R}}} \cr } } \right]}.$$

In order to guarantee the stability of the proposed method, the spectral radius of the square matrix $\bar{{\bf \Phi }}$ (that is the largest absolute value of its eigenvalues) must be less than one ($\rho ( {\bar{{\bf \Phi }}} ) \le 1$) [Reference Kantartzis and Tsiboukis13]. It will be shown in the next section that this condition is satisfied for 0.5 ≤ α ≤ 1, and consequently, the proposed method is unconditionally stable even for arbitrarily large time steps.

Generalized FDTD solution method for nonlinearly loaded TLs

As shown in Fig. 3(a), consider a TL that is terminated with a Schottky diode as a nonlinear load and excited by the external EM wave illumination. For incorporating the nonlinear terminal conditions in the proposed FDTD algorithm, an accurate large-circuit model of the nonlinear device must be used. Thus, the nonlinear lumped-circuit model of the Schottky diode investigated in this paper, which includes the intrinsic part with the nonlinear elements and the extrinsic parasitic elements, is presented in Fig. 3(b). In the PN-junction diode model [Reference Antognetti and Massobrio14], which is used for modeling of the intrinsic part, the nonlinear diode current and capacitor are given by:

(34)$$\eqalign{& I_d( {V_{dj}} ) \cr & = \left\{{\matrix{ {I_S( {e^{{-}10}-1} ) + ( {I_Se^{{-}10}/( {N\cdot V_t} ) } ) ( {V_{dj}-10N\cdot V_t} ) } & {V_{dj} < -10N\cdot V_t} \cr {I_S( {e^{V_{dj}/N.V_t}-1} ) } & \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! {-10N\cdot V_t \le V_{dj} \le V_{max}} \cr {I_{max} + ( {V_{dj}-V_{max}} ) ( {I_E + I_S} ) /( {N\cdot V_t} ) } & {V_{dj} > V_{max}} \cr } } \right.,} \;$$

(35)$$C_{dj}( {V_{dj}} ) = \left\{{\matrix{ {Area\cdot C_{\,jo}{( {1-V_{dj}/V_j} ) }^{{-}M}} & {V_{dj} \le F_CV_j} \cr {Area\cdot \displaystyle{{C_{\,jo}} \over {1-F_C^M }}\left[{1 + \displaystyle{{M( V_{dj}-F_CV_j) } \over {V_j( 1-F_C) }}} \right]} & {V_{dj} > F_CV_j} \cr } } \right., \;$$

wherein V _max = N ⋅ V _tLn(I _E/I _S + 1). Also, I _S, I _E, N, V _t, I _max, Area, C _jo, V _j, M, and F _C are the model parameters of the diode, and V _dj denotes the internal junction voltage. The descriptions and values of model parameters for the Schottky diode used in this paper are presented in Table 1. Therefore, applying Kirchhoff's law at the nonlinear load termination network shown in Fig. 3 gives the relations for the load terminal as:

(36)$$-V_{N_x + 1}-V_D + 2L_p\displaystyle{{\partial I_L} \over {\partial t}} = 0, \;$$

(37)$$V_{dj}-V_D + R_DI_{LD} + L_D\displaystyle{{\partial I_{LD}} \over {\partial t}} = 0, \;$$

(38)$$I_L + I_{LD} + C_D\displaystyle{{\partial V_D} \over {\partial t}} = 0, \;$$

(39)$$-I_{LD} + C_{dj}( {V_{dj}} ) \displaystyle{{\partial V_{dj}} \over {\partial t}} = {-}I_d( {V_{dj}} ) .$$

Fig. 3.

(a) Illumination of the incident EM wave to a TL terminated with a diode. (b) The schematic diagram of the large-signal circuit model of the Schottky diode.

Table 1.

Model parameters of the diode SMS7621-001LF.

Moreover, substituting Δx ↠ 0.5Δx and i =  N _x + 1 into equation (10) gives:

(40)$$I_L-I_{N_x} + 0.5g\Delta xV_{N_x + 1} + 0.5c\Delta x\displaystyle{{\partial V( {x, \;t} ) } \over {\partial t}} = 0.5\Delta xI_E {N_x + 1} .$$

As in the previous subsection, representing equation (10) and the constraint equations, including the nonlinear load termination circuit equations (38)–(40) and source terminal equation (11) together in the matrix form results in:

(41)$${ \bar{\bf N}\bar{\bf I}} + { \bar{\bf G}\bar{\bf V}} + \bar{{\bf C}}\displaystyle{\partial \over {\partial {\rm t}}}\bar{{\bf V}} = \overline {{\bf I}_{\bf S}} , \;$$

wherein:

(42)$$\bar{{\bf V}}\equiv \left[{\matrix{ {V_1} & \cdots & {V_{N_x}} & {V_{N_x + 1}} & {V_D} & {V_{dj}} \cr } } \right]^T_{( {N_x + 3} ) \times 1} , \;$$

(43)$$\bar{{\bf I}}\equiv \left[{\matrix{ {I_1} & \cdots & {I_{N_x-1}} & {I_{N_x}} & {I_L} & {I_{LD}} \cr } } \right]^T_{( {N_x + 2} ) \times 1} , \;$$

(44)$$\overline {{\bf I}_{\bf S}} \equiv \left[{\matrix{ {0.5\Delta xI_E 1 + V_S/R_S} \cr {\Delta xI_E 2 } \cr \vdots \cr {\Delta xI_E {N_x} } \cr {0.5\Delta xI_E {N_x + 1} } \cr 0 \cr {-I_d( {V_{dj}} ) } \cr } } \right]_{( {N_x + 3} ) \times 1}, \;$$

(45)$$\bar{{\bf N}}\equiv \left[{\matrix{ 1 & 0 & \cdots & \cdots & \cdots & 0 \cr {-1} & 1 & \ddots & \ddots & \ddots & \vdots \cr 0 & \ddots & \ddots & \ddots & \ddots & \vdots \cr \vdots & \ddots & {-1} & 1 & \ddots & \vdots \cr \vdots & \ddots & \ddots & {-1} & 1 & 0 \cr \vdots & \ddots & \ddots & \ddots & 1 & 1 \cr 0 & \cdots & \cdots & \cdots & 0 & {-1} \cr } } \right]_{( {N_x + 3} ) \times ( {N_x + 2} ) }, \;$$

(46)$$\bar{{\bf C}}\equiv \left[{\matrix{ {0.5c\Delta x} & 0 & \cdots & \cdots & \cdots & \cdots & 0 \cr 0 & {c\Delta x} & \ddots & \ddots & \ddots & \ddots & \vdots \cr \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \cr \vdots & \ddots & \ddots & {c\Delta x} & \ddots & \ddots & \vdots \cr \vdots & \ddots & \ddots & \ddots & {0.5c\Delta x} & \ddots & \vdots \cr \vdots & \ddots & \ddots & \ddots & \ddots & {C_D} & 0 \cr 0 & \cdots & \cdots & \cdots & \cdots & 0 & {C_{dj}( {V_{dj}} ) } \cr } } \right]_{( {N_x + 3} ) \times ( {N_x + 3} ) }, \;$$

(47)$$\bar{{\bf G}}\equiv \left[{\matrix{ {1/R_S + 0.5g\Delta x} & 0 & \cdots & \cdots & \cdots & \cdots & 0 \cr 0 & {g\Delta x} & \ddots & \ddots & \ddots & \ddots & \vdots \cr \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \cr \vdots & \ddots & \ddots & {g\Delta x} & \ddots & \ddots & \vdots \cr \vdots & \ddots & \ddots & \ddots & {0.5g\Delta x} & \ddots & \vdots \cr \vdots & \ddots & \ddots & \ddots & \ddots & 0 & 0 \cr 0 & \cdots & \cdots & \cdots & \cdots & 0 & 0 \cr } } \right]_{( {N_x + 3} ) \times ( {N_x + 3} ) }.$$

Moreover, equations (9), (36), and (37) can be similarly written in the matrix form as follows:

(48)$${ \bar{\bf M}\bar{\bf V}} + { \bar{\bf R}\bar{\bf I}} + \bar{{\bf L}}\displaystyle{\partial \over {\partial t}}\bar{{\bf I}} = \overline {{\bf V}_{\bf S}} , \;$$

wherein:

(49)$$\bar{{\bf {\rm M}}}\equiv \left[{\matrix{ {-1} & 1 & 0 & \cdots & \cdots & \cdots & 0 \cr 0 & {-1} & 1 & \ddots & \ddots & \ddots & \vdots \cr \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \cr \vdots & \ddots & \ddots & {-1} & 1 & \ddots & \vdots \cr \vdots & \ddots & \ddots & \ddots & {-1} & {-1} & 0 \cr 0 & \cdots & \cdots & \cdots & 0 & {-1} & 1 \cr } } \right]_{( {N_x + 2} ) \times ( {N_x + 1} ) }, \;$$

(50)$$\overline {{\bf V}_{\bf S}} \equiv \Delta x.\left[{\matrix{ {V_E 1 } & \cdots & {V_E {N_x} } & {V_E {N_x + 1} } & 0 & 0 \cr } } \right]^T_{( {N_x + 2} ) \times 1} , \;$$

(51)$$\bar{{\bf R}}\equiv \left[{\matrix{ {r\Delta x} & 0 & \cdots & \cdots & 0 \cr 0 & \ddots & \ddots & \ddots & \vdots \cr \vdots & \ddots & {r\Delta x} & \ddots & \vdots \cr \vdots & \ddots & \ddots & 0 & 0 \cr 0 & \cdots & \cdots & 0 & {R_D} \cr } } \right]_{( {N_x + 2} ) \times ( {N_x + 2} ) }, \;$$

(52)$$\bar{{\bf L}}\equiv \left[{\matrix{ {l\Delta x} & 0 & \cdots & \cdots & 0 \cr 0 & \ddots & \ddots & \ddots & \vdots \cr \vdots & \ddots & {l\Delta x} & \ddots & \vdots \cr \vdots & \ddots & \ddots & {2L_p} & 0 \cr 0 & \cdots & \cdots & 0 & {L_D} \cr } } \right]_{( {N_x + 2} ) \times ( {N_x + 2} ) }.$$

Finally, using the forward difference approximation and substituting the proposed average vectors (equations (27)–(30)) into the time-domain equations (41) and (48), the matrix form of the final recursion relation for the nonlinearly loaded line can be determined as:

(53)$${\boldsymbol F}( {\bar{{\boldsymbol S}}} ) = \bar{{\boldsymbol S}}^{n + 1}-\bar{{\bf \Phi }}\bar{{\boldsymbol S}}^n-\bar{{\bf \Psi }}( {{\boldsymbol \alpha }{\bar{{\boldsymbol U}}}^{n + 1} + ( {1-\alpha } ) {\bar{{\boldsymbol U}}}^n} ) = 0$$

where $\bar{{\boldsymbol S}} = \left[{\matrix{ {\bar{{\bf V}}} \cr {\bar{{\bf I}}} \cr } } \right]$ and $\bar{{\boldsymbol U}} = \left[{\matrix{ {\overline {{\bf I}_{\boldsymbol S}} } \cr {\overline {{\bf V}_{\boldsymbol S}} } \cr } } \right]$. Also, matrices $\bar{{\bf \Psi }}$ and $\bar{{\bf \Phi }}$ are given in (32) and (33), respectively. An iterative root-finding method can be used to solve the nonlinear equation (53). Since locally convergent methods are not converged for any arbitrary initial conditions in the highly nonlinear equations, a globally convergent approach is utilized to ensure the stability and convergence of the presented algorithm. Applying the global Newton method [Reference Kelley15] to solve the nonlinear equation (53) results in:

(54)$$\bar{{\boldsymbol S}}_{k + 1}^{n + 1} \Delta \;\bar{{\boldsymbol S}}_k^{n + 1} -\delta _k^\ast {\boldsymbol J}( {\bar{{\boldsymbol S}}_k^{n + 1} } ) ^{{-}1}{\boldsymbol F}( {\bar{{\boldsymbol S}}_k^{n + 1} } ) ^T, \;$$

where the Jacobian matrix of F (${\boldsymbol J} = \partial {\boldsymbol F}( {\bar{{\boldsymbol S}}} ) /\partial \bar{{\boldsymbol S}}$) can be obtained as follows:

(55)$$\boldsymbol{J}\left({\overline{\boldsymbol{S}}}^{\boldsymbol{n}\boldsymbol{+}\boldsymbol{1}}\right)=\overline{\boldsymbol{I}}\boldsymbol{-}{\overline{\boldsymbol{\mathrm{\Phi }}}}^{\boldsymbol{'}}{\overline{\boldsymbol{S}}}^n\boldsymbol{+}\overline{\boldsymbol{\mathrm{\Phi }}}\boldsymbol{-}\overline{\boldsymbol{\mathrm{\Psi }}}\left(\alpha {\overline{\boldsymbol{\mathrm{U}}}}^{\boldsymbol{n}\boldsymbol{+}\boldsymbol{1}}\boldsymbol{'+}\left(1-\alpha \right){\overline{\boldsymbol{\mathrm{U}}}}^{\boldsymbol{n}}\boldsymbol{'}\right)-\overline{\boldsymbol{\mathrm{\Psi }}}\boldsymbol{'}\left(\boldsymbol{\alpha }{\overline{\boldsymbol{U}}}^{n+1}+\left(\mathrm{1-}\alpha \right){\overline{\boldsymbol{U}}}^n\right).$$

The parameter $\delta _k^{\rm \ast }$ is called the Newton step and must be obtained in each iteration to ensure the convergence of the proposed method (${\boldsymbol F}( {{\bar{{\boldsymbol S}}}_{{\boldsymbol k} + 1}} ) \ll {\boldsymbol F}( {{\bar{{\boldsymbol S}}}_{\boldsymbol k}} )$).

As mentioned above, the proposed algorithm is applicable in the frequency range that the quasi-TEM mode assumption is satisfied. For the microstrip circuits, the quasi-TEM propagation mode is limited to the cutoff frequency, which is given by [Reference Leone and Singer16]:

(56)$$f_{QTEM} = 0.0213/\left({\sqrt {\varepsilon_r + 1} ( W + 2d) } \right), \;$$

where d and ε _r are the thickness and relative permittivity of the substrate, respectively, and W denotes the trace width. According to (56), the cutoff frequency range for typical microstrip lines (d/5 ≤ W ≤ 3d, 1 ≤ ɛ_r ≤ 6) supporting the quasi-TEM mode is 2–10 GHz. Moreover, it is worth noting that the large-signal lumped-circuit models for the nonlinear components are not valid at high frequency since the dimensions of the active devices could not be ignored compared to the wavelength of the incident field. Under these conditions, the distributed TL models [Reference Seifi, Abdipour and Mirzavand17] for nonlinear devices could be used in the presented algorithm.

EM coupling of a single microstrip line with linear termination

As a first example, a single microstrip TL on the RO4003 substrate (ɛ_r = 3.38 and d = 0.8 mm) is investigated to validate the proposed method (see Fig. 4(a)). The line has a width of 1.8 mm and a length of 30 mm and is terminated by linear loads. The generalized FDTD method is first applied to analyze the microstrip line without the incident wave, and the terminal conditions are R _s = 50 Ω, R _L = 50 Ω, and V _s = sin(2π10⁹t). The stability of the proposed FDTD method for different values of α is investigated based on the spectral radius of the amplification matrix. The values of $\rho ( {\bar{{\bf \Phi }}} )$ versus the time step are obtained from equation (33), and the results are shown in Fig. 6. It can be observed that the proposed FDTD method for 0.5 ≤ α ≤ 1 is unconditionally stable since, in this case, the spectral radius of the amplification matrix for any time step size is less than or equal to one ($\rho ( {\bar{{\bf \Phi }}} ) \le 1$). Furthermore, as Fig. 6 shows, $\rho ( {\bar{{\bf \Phi }}} )$ is equal to one for the Crank–Nicholson scheme (α = 0.5). As previously stated, unconditional stability does not mean that the solution is accurate and physically realistic for any large step size, and the accuracy of the numerical method must also be considered. The maximum percentage error between the proposed method results and the exact solution for the linearly loaded line is shown in Fig. 7.

Fig. 6.

The spectral radius of the amplification matrix of the proposed FDTD method for different values of α.

Fig. 7.

The maximum percentage error between the proposed method results and the exact solution for the linearly loaded line versus the time step for different values of α.

It is worth noting that for no incident field illumination, the exact closed-form solution for the time-domain voltage and current on each point of the line, which is shown in Fig. 2, are given by [Reference Paul5]:

(57)$$V( {x, \;t} ) = V^ + e^{-\alpha x}\cos ( {\omega t-\beta x + \theta^ + } ) + V^- e^{-\alpha x}\cos ( {\omega t + \beta x + \theta^-} ) , \;$$

(58)$$I( {x, \;t} ) = \displaystyle{{V^ + } \over {Z_C}}e^{-\alpha x}\cos ( {\omega t-\beta x + \theta^ + {-}\theta_{ZC}} ) -\displaystyle{{V^-} \over {Z_C}}e^{\alpha x}\cos ( {\omega t + \beta x + \theta^-{-}\theta_{ZC}} ) , \;$$

wherein:

(59)$$\bar{Z}_C\equiv Z_Ce^{\,j\theta _{ZC}} = \sqrt {\displaystyle{{r + j\omega l} \over {g + j\omega c}}} , \;$$

(60)$$\gamma = \alpha + j\beta = ( {r + j\omega l} ) ( {g + j\omega c} ) , \;$$

(61)$$\left[{\matrix{ {{\bar{V}}^ + } \cr {{\bar{V}}^-} \cr } } \right]\equiv \left[{\matrix{ {V^ + e^{\,j\theta^ + }} \cr {V^- e^{\,j\theta^-}} \cr } } \right] = \left[{\matrix{ {{\bar{Z}}_C + R_S} & {{\bar{Z}}_C-{\bar{Z}}_S} \cr {( {{\bar{Z}}_C-R_L} ) e^{-\gamma L}} & {( {{\bar{Z}}_C + R_L} ) e^{\gamma L}} \cr } } \right]^{{-}1}\left[{\matrix{ {{\bar{Z}}_C{\bar{V}}_S} \cr {{\bar{Z}}_C{\bar{V}}_L} \cr } } \right].$$

In the above equations, $\bar{Z}_C$ and γ are the characteristic impedance and the propagation constant of the line, respectively. Since the fully implicit FDTD scheme (α = 1) is an O(Δt ²) perturbation to the Crank–Nicholson method (α = 0.5), resulting in second-order truncation error terms, the accuracy of this approach degrades quickly in the time step above the CFL limit [Reference Crawford and Workman18]. Therefore, as can be observed in Fig. 7, the proposed FDTD method with the smaller weighting factor (α) gives more accurate solutions for large time steps. In other words, in problems that require large time steps to reduce the simulation time, the proposed method with α = 0.5, while ensuring stability, gives lower errors than the other schemes with larger values of α. Moreover, as shown in Fig. 7, when α becomes large for small time steps, the error of the proposed method decreases. It can be concluded that by selecting the optimum value of weighting factor, the proposed FDTD method, in addition to ensuring the stability of the solution for a wide range of time step sizes, presents good accuracy.

As mentioned in the introduction, in the conventional explicit FDTD method, the time step size must be limited to the CFL number (Δt _CFL = Δx/v _max) to avoid the instability of the solution. For evaluating the efficiency of the proposed method, the spectral radius of the square matrix $\bar{{\bf \Phi }}$ and computational time of the conventional explicit FDTD approach and the proposed algorithm for different time step sizes are presented in Table 2. As can be observed from this table, the computational time is reduced significantly in the proposed algorithm compared to the conventional explicit FDTD method, while good accuracy is obtained. Furthermore, Table 2 shows that the spectral radius $\rho ( \bar{{\bf \Phi }})$ in the proposed algorithm is less than or equal to 1 for all time step sizes, while in the conventional explicit FDTD algorithm, the spectral radius $\rho ( \bar{{\bf \Phi }})$ becomes larger than 1 in time step sizes beyond Δt _CFL, i.e., the algorithm becomes unstable.

Table 2.

Comparison of the spectral radius of the square matrix $\bar{{\bf \Phi }}$ and computational time of the conventional explicit FDTD approach and the proposed algorithm.

In Fig. 8, the proposed method is used to compute the magnitude of the induced voltage and current distributions on the line endfire illumination ($\theta _E = {-}{\rm \;}\varphi _E = \phi = 90^\circ$) with a frequency of 1.3 GHz and an electric-field amplitude of 100 V/m. According to (5) and (6), in the endfire excitation, the incident electric field is polarized in the +y direction, and its amplitude depends only on the x coordinate ($\vec{E}_{inc}( {x, \;t} ) = \vec{a}_yf_W\left({t + \sqrt {\mu \varepsilon } x} \right)$). Therefore, since the induced sources due to the incident field (equations (7) and (8)) are symmetric to the center of the line, the magnitude of the induced voltage and current at the line terminals, which are terminated to matched load, will be equal, as shown in Fig. 8. Moreover, as expected, maximum voltage and maximum current are coupled to the line in the open-circuited and short-circuited terminals, respectively. Further, as Fig. 8 shows, the simulated results are in excellent agreement with analytical reference data [Reference Leone and Singer16]. Furthermore, Fig. 9 shows the measured and simulated results for the induced output terminal voltage of the single microstrip line versus the incident electric field amplitude at 1.3 GHz in the case of broadside excitation ($\theta _E = {\rm \;}\phi = 90^\circ , \;{\rm \;}\varphi _E = 0^\circ$). The line is terminated with the matched load. As expected, the increase in the incident EM wave amplitude has led to a rise in the coupled voltage of the TL. Besides, a good agreement between the measured and analytical results is observed in Fig. 9, which confirms the validity of the proposed algorithm.

Fig. 8.

The magnitude of the induced (a) voltage and (b) current distributions on the single microstrip line for endfire illumination at 1.3 GHz and electric-field amplitude 100 V/m.

Fig. 9.

Measured and simulated results for the induced output terminal voltage of a single microstrip line versus the incident electric field amplitude for broadside illumination at 1.3 GHz.

Diode power detector excited by uniform EM wave illumination

In this subsection, the proposed method is applied to the power detector circuit that generates DC voltage by rectifying the RF/microwave signal through a diode. Furthermore, for verification of the paper method, the designed detector circuit is implemented on a two-layer RO4003 PCB, as shown in Fig. 5(b). The parameters of the circuit are also tabulated in Table 3. The schematic and the PCB layout of the fabricated power detector circuit are also presented in Fig. 10. It should be noted that since the dimensions of the TLs of the detector circuit are much smaller than the wavelength of the incident field at low microwave frequencies, the quasi-TEM assumption of the proposed method is valid.

Fig. 10.

(a) The schematic diagram and (b) the PCB layout of the fabricated power detector circuit.

Table 3.

Parameters of the fabricated power detector circuit.

Transient analysis of the detector circuit, which is excited by −1 dBm CW source at 1 GHz, is carried out using the proposed method and ADS Agilent software. Figure 11 shows the time-domain waveforms of the output voltage (V _out), diode current (I _d), and diode internal junction voltage (V _dj) of the circuit, which indicates a good agreement between the two methods. The proposed method results of the detector circuit, which is subjected to 300 V/m incident EM field illumination at 1.1 GHz, are also presented in Fig. 11. In this section, the detector is subjected to endfire illumination by the uniform sinusoidal plane wave. As shown in Fig. 11(a), the incident field coupling to the detector circuit has distorted and increased the output voltage. Further, the mixing of the desired detector input signal and the induced signal due to the nonlinear operation of the diode can be observed in the voltage and current waveforms in Figs 11(b) and 11(c).

Fig. 11.

Time-domain waveforms of (a) the output voltage, (b) diode current, and (c) diode internal junction voltage of the power detector circuit, which is excited by −1 dBm CW source at 1 GHz as well as the incident plane wave illumination at 1.1 GHz.

Besides, Fig. 12 shows the measured and proposed method results for the detected DC output voltage versus the input power of the detector circuit subjected to the external uniform plane wave with an amplitude of 500 V/m at 1.1 GHz. The good agreement between the results confirms the accuracy of the proposed method. As can be seen, the high power incident field illumination on the detector circuit has increased the output voltage and also forced the detector to saturate with lower input power. Furthermore, the simulation time and error of the proposed FDTD method for analyzing the diode detector circuit versus the weighting factor parameter (α) are shown in Fig. 13. It is apparent from this figure that good accuracy can be achieved in the shortest possible simulation time by selecting the optimal value for the parameter α in the proposed method.

Fig. 12.

The measured and proposed method results for the detected DC output voltage versus the input power of the detector circuit, which is excited by the CW source at 1 GHz and the incident plane wave illumination at 1.1 GHz.

Fig. 13.

The simulation time and error of the proposed FDTD method for analyzing the diode detector circuit versus the parameter α.