Debye Dispersion
Ampere’s law is written as
(2-9)
The total current is given by
(2-10)
where
is the induction current,
is the electric displacement in the high frequency limit (in the sense that the excitation duration is smaller than the shortest charge relaxation time).
The charge relaxation is modeled using the polarization current (Ref. 1)
where each term obeys an equation of the following form:
(2-11)
where τm is the relaxation time, and Δεm is the relative permittivity contribution. Note that both quantities can be taken as diagonal matrices to cover the anisotropic case.
Introduce new variables, auxiliary electric field from
Equation 2-11 can be integrated once with respect to time assuming zero initial fields
(2-12)
The polarization current can be expressed the as either
(2-13)
or
By applying the divergence operator, Equation 2-9 yields the current conservation equation:
Assume that the electric conductivity is small, so that the induction current can be neglected. Using Equation 2-13, the current conservation equation can be written as
which can be integrated once with respect to time assuming zero initial fields to give the charge conservation equation
(2-14)
Equation 2-14 can be solved using the scalar electric potential V as a dependent variable determine the electric field as , and it needs to be solved together with N vector Equation 2-12 for N vector dependent variables em. In time domain, these auxiliary field variables can be treated as the local state variables, and each of N corresponding ODE equations can be integrated locally.
Frequency domain
In frequency domain, Equation 2-13 for the auxiliary variables can be solved analytically, so that Equation 2-14 gives
where ω is the angular frequency. By separating the real and imaginary parts, the equation can be further rewritten as
where
and
Each term in the sum in the above equations represents a pole on the complex w plane. Thus, the model is referred to as Multipole Debye Dispersion model. To determine the number of poles needed to approximate the lossy material behavior can be a difficult task. One possible approach is presented in Constant loss tangent section.
Eigenfrequency study
Direct use of the complex frequency-dependent permittivity would lead to a nonlinear eigenvalue problem, for which the solution would depend on the eigenvalue linearization point (transform point frequency, specified on the solver). The solution would be accurate only for the eigenfrequency that is close to the transform point frequency. Thus, the eigenfrequencies would need to be computed one by one using iterative updates of the transform point frequency value.
Instead, one can solve N algebraic equations for the auxiliary variables em
(2-15)
together with Equation 2-15. Even though this approach requires extra degrees of freedom, it will produce frequency independent contributions to the damping and stiffness matrices. The corresponding linear damped eigenvalue problem can be solved using the default solver in a single run for any specified number of eigenfrequencies which will be computed exactly.
Constant loss tangent
In frequency domain, the loss tangent is defined as
(2-16)
and it is often used to characterize lossy materials based on experimental data.
In Ref. 2, a special parameterization was suggested for the multipole Debye model based on the original result derived in Ref. 3 for equivalent electric circuits. It allows to fit the model parameters for Multipole Debye Model so that the loss tangent is nearly constant at a certain frequency range.
The input parameters are the loss tangent η(fc) together with the center frequency fc, at which is has been measured, and the model bandwidth nd that defines a frequency interval (in decades) centered at fc, in which the loss tangent will be approximately constant and equal to η(fc).
The relaxation times are computed as , where the corresponding frequencies are equidistantly spaced in the logarithmic space as
where N is the number of Debye poles, and .
The corresponding relative permittivity contribution is given by
where k > 1 is the spacing parameter. Note that the relative electric permittivity in the low and high frequency limits are related simply as .
COMSOL Multiphysics software will automatically deduce the necessary number of poles N together with the values of the relaxation times tm and relative permittivity contributions Δεm, which will be used in computations to maintain the requested bandwidth and accuracy.
The accuracy can be selected as normal to give
for
which requires poles, and the spacing parameter with the loss angle d computed from η(fc).
Alternatively, the accuracy can be selected as high to get
in the same interval, which will require poles, and the spacing parameter .
Thus computed parameters will be used automatically in a Time Dependent study using Equation 2-12, and in an Eigenfrequency study using Equation 2-16.
Temperature Effects
For many dielectric materials such as polymers, the charge relaxation properties have a strong dependence on the temperature. A common assumption is that a change in the temperature can be transformed directly into a change in the time scale. The reduced time is defined as
where αT(T) is a temperature-dependent shift function.
The implication is that the problem can be solved using the original material data, provided that the time is transformed into the reduced time. Think of the shift function αT(T) as a multiplier to the relaxation time in Equation 2-12 and Equation 2-16.
Vogel-Fulcher Shift
The shift function is given by
here, a base-e logarithm is used, Q is the activation energy (SI unit: J/mol), and R is the universal gas constant. Temperature T0 is the temperature at which the relaxation time would becomes infinite.
Arrhenius Shift
The shift function is assumed to follow an Arrhenius law
here, a base-e logarithm is used, Q is the activation energy (SI unit: J/mol), and R is the universal gas constant.
Williams–Landel–Ferry Shift
One commonly used shift function is defined by the WLF (Williams–Landel–Ferry) equation:
where a base-10 logarithm is used. This shift is only valid over a certain temperature range, typically around the glass transition temperature.
The first step to compute the shift function αT(T) consists of building a master curve based on experimental data. To do this, the curves of the dielectric properties versus time or frequency are measured at a reference temperature Tref. Then, the same properties are measured at different temperatures.
The shift value of each curve, with respect to the master curve obtained at the temperature Tref, defines the shift factor αT(T). The constants C1 and C2 are material dependent and are calculated after plotting log(αT) versus T − Tref.
The shift factor at the reference temperature equals αT(Tref) = 1, so that Tref is the temperature at which the master curve is given. If the temperature T drops below Tref − C2, the WLF equation is no longer valid.
Since the master curve is measured at an arbitrary reference temperature Tref, the shift factor αT(T) can be derived with respect to any temperature, and it is commonly taken as the shift with respect to the glass transition temperature.
Tool–Narayanaswamy–Moynihan Shift
Charge relaxation in glass can be modeled using the so-called Tool–Narayanaswamy–Moynihan shift factor which is given as
here, a base-e logarithm is used, Q is the activation energy (SI unit: J/mol), R is the universal gas constant, T is the current temperature, Tref is a reference temperature, χ is a dimensionless activation energy fraction, and Tf is the so-called fictive temperature. The fictive temperature is given as the weighted average of partial fictive temperatures.
with
Here, wi are the weights and Tfi are the partial fictive temperatures. The partial fictive temperatures are determined from a system of coupled ordinary differential equations (ODEs) which follow Tool’s equation
here, λ0i is a relaxation time parameter.
References for debye dispersion
1. N.S. Stoykov, T.A. Kuiken, M.M. Lowery, and A. Taflove, “Finite-element time-domain algorithms for modeling linear Debye and Lorentz dielectric dispersions at low frequencies,” IEEE Transactions on Biomedical Engineering, vol. 50, no. 9, pp. 1100–1107, 2003.
2. A.E. Engin, I. Ndip, K. Lang, and J. Aguirre, “Closed-Form Multipole Debye Model for Time-Domain Modeling of Lossy Dielectrics,” IEEE Trans. on Electromagnetic Compatibility, vol. 61, no. 3, pp. 966–968, 2019.
3. R. Morrison, “RC constant-argument driving-point admittances,” IRE Trans. Circuit Theory, vol. 6, no. 3, pp. 310–317, 1959.
4. J.A. Diego, J. Sellares, A. Aragoneses, M. Mudarra, J.C. Canadas, and J. Belana, “TSDC study of the glass transition: correlation with calorimetric data,” Journal of Physics D: Applied Physics, vol. 40, pp.  138–1145, 2007.