Lorentz Wideband Dielectric Model, Problems with big numbers?

[toammann]

Hello,
currently I am trying to implement the wideband Djordjevic-Sarkar model for openEMS. The my chosen workflow is:

1. Calculate a wideband model from the Djordjevic-Sarkar from a single epsR and trand Value
2. Fit a wideband debye model using mulitple terms to the result of 1)
   
3. Fit the openEMS Lortentz model to the wideband Debye Model

The fit to the Debye equation in Step 2) is necessary because the coefficents can be calulated from measured data (or the result of the Djordjevic-Sarkar) model easily. Debye model fits are quite common in commercial solvers

The Debye model is a a simplification of the more general Lorentz model. However, it is not possible to calculate the Lorentz model prameters from the Debye ones directly beause of the f^2 term in the demoninator if the Lorentz model.

I therefore made an approximation by choosing:

In principle this seems to work well. For a tand=0.02 and Er=4.3 result of the fit looks promising:

fp =

       933.5751e+09
       892.2961e+09
       894.3945e+09
       892.0289e+09
       899.3542e+09
       875.9840e+09
       948.0429e+09
       608.2609e+09

    fl = 5.0329e+012
    taup =

       100.0000e-24
         1.0000e-21
        10.0000e-21
       100.0000e-21
         1.0000e-18
        10.0000e-18
       100.0000e-18
         1.0000e-15

    epsInf = 4.1317e+000
    sigma = 80.0000e-012

The curves of the Debye and openEMS Lorentz model are on top of each other. I have created a second function which creates the multi term material based on the calculated poles.

---

I am currently testing the models in openEMS and this is where I am facing thre first problems. The numbers for fp are very big. The values for taup are very small. The higher the upper operating frequency is chosen the higher and lower these numbers become. If, in openEMS fp/taup is calculated the number will explode :)

In the equation these numbers are squared what makes the problem even worse. For corner frequencies >200GHz (a typical value) the openEMS simulation does not converge.The energy in the model goes up to infinity.

Create FDTD operator (compressed SSE + multi-threading)
FDTD simulation size: 107x77x44 --> 362516 FDTD cells
FDTD timestep is: 7.63218e-14 s; Nyquist rate: 262 timesteps @2.50046e+10 Hz
Excitation signal length is: 3003 timesteps (2.29194e-10s)
Max. number of timesteps: 1000000000 ( --> 333000 * Excitation signal length)
Create FDTD engine (compressed SSE + multi-threading)
Running FDTD engine... this may take a while... grab a cup of coffee?!?
[@        5s] Timestep:          130 || Speed:    8.3 MC/s (4.367e-02 s/TS) || Energy: ~6.75e-22 (- 0.00dB)
[@        9s] Timestep:          260 || Speed:   11.0 MC/s (3.302e-02 s/TS) || Energy: ~5.33e-21 (- 0.00dB)
.
.
[@     1m12s] Timestep:         2535 || Speed:   12.6 MC/s (2.888e-02 s/TS) || Energy: ~1.90e+20 (- 0.00dB)
[@     1m18s] Timestep:         2730 || Speed:   12.7 MC/s (2.854e-02 s/TS) || Energy: ~2.64e+24 (- 0.00dB)
[@     1m23s] Timestep:         2925 || Speed:   12.9 MC/s (2.807e-02 s/TS) || Energy: ~   inf (-  nandB)

I suppose that in openEMS a floating point overlow occures which results in corrupted model paramers. When choosing a less good approximation I can scale the numbers down by one or two number of magnitudes. This can be adjusted choosing another value for kLor in line 166 of "calcDjordjevicSarkarApprox.m."

The Lorentz parameters are in this case:

    fp =

        93.3575e+09
        89.2296e+09
        89.4395e+09
        89.2029e+09
        89.9354e+09
        87.5984e+09
        94.8043e+09
        60.8261e+09

    fl = 503.2921e+009
    taup =

        10.0000e-21
       100.0000e-21
         1.0000e-18
        10.0000e-18
       100.0000e-18
         1.0000e-15
        10.0000e-15
       100.0000e-15

    epsInf = 4.1317e+000
    sigma = 80.0000e-012

The simulation runs just fine in this case and the simulation converges. The results are also correct. Here is an example of a microstrip line.

I don´t know if this should be considered a "bug" in openEMS or it is problem with my "approximation" of the Lorentz models. The numbers for the poles fit well into the frame of realMin=22.2507e-309, realMax = 179.7693e+306. On their own I would not consider them to be too big.

The cleanest way would be to implement the Debye model equation in openEMS. The debye model can also be fitted to measured data very easily. This is not the case for the Lorentz model. However I suspect that this is a signifcant effort.

The usage of my proposed function is very easy. The wideband material can be created with a single line:

 calcDjordjevicSarkarApprox(4.3, 0.02, 10e9, 200e9, 'f1', 10^4/(2*pi), 'sigma',...
 80e-12, 'nTermsPerDec', 1, 'plotEn')

Maybe I am completly wrong here and somethig else is causing the simulation instablity. However the magnitudes of the numbers is the only differnce I can spot here...

Regards
Tobias

permitttivity_model.zip

[gadiLhv]

Hey Tobias

I really need to check what you built there a lot more carefully, but I'll try to add some data to the mix.

In general, all these calculations for the dielectric models happen in parallel to this capacitor, in the EC-FDTD YEE cell

It's value is usually set by the volumne of the cell, and the dielectric constant of the material there with the highest priority.

While developing some engine extension, I notice that if this capacitance value is set, for some reason, to be very small (e.g., open circuit), these misconvergences you saw start to happen.

I'm not taking off the table that it could be a floating point precision issue, but this could also be a point to look about.

Cheers

[toammann]

Hey gadiLhv,
thank you very much for your comment. Do you know a frequency dependent permittivy is considered here? In this case epsilon is not a single number.

I suppose that the floating point issue (if there is one) happens before the simulation starts while openEMS evaluates the Lorentz equation by using the user input data. The value of epsilon is almost identical in both cases (The factor kLor specifies how closely the Lorentz model will follow the Debye model for higher frequencies).

The only difference is that the input numbers for the Lorentz model are scaled by a vactor of 10 (sqrt(kLor)).

Debye-Lorentz fit scaling factor	First Term plasma frequency	epsR @ 10GHz
kLor = 1/1000	fp(1) = 933.5751e+09	4.3001e+00 - 91.5526e-03i
kLor = 1/10	fp(1) = 93.3575e+09	4.3002e+00 - 91.5928e-03i

So, the output of the Lorentz equation is not diffrent. But the magnitudes of the input numbers into the model are.

edit: Here is a plot showing epsR over a broad frequency range. The resuluts of the Lorentz model equation are neraly identical for the two values of kLor above. Some small deviations are visible around 1000GHz.

Regards,
Tobias

[toammann]

Hey,
I have just noiced that I forgot to add the functions to my first post. Sorry :)

Here is the octave code:
permitttivity_model.zip

Regards,
Tobias

[gadiLhv]

Hey @toammann

Thank you for supplying the codes, for reference. This way I can export the XML file and debug the C++ code more easily.

I didn't really have the time yet to debug your request until now, hopefully tomorrow.
This is where I would put my effort:

The Lorentz model's coefficients are calculated as such (I'm attaching the article here):

Numerical instablity can happen while calculating these (It happened to me while using a bad formulation). I need to debug the numbers slowly though.

Will let you know

Equivalent_Circuit_EC_FDTD_Method_for_the_Modeling.pdf

[gadiLhv]

Hey @toammann

I briefly checked out the numbers. All the values produced so far seem to be Ok. Here is one example:

One thing that didn't make sense is that you chose the corner\cutoff frequency to be the same as the center frequency. DC calculations are a bit problematic, but let's forget that for now.

Can you supply me with an example script for something that causes the instability you noticed? This will help me see where the issues are.

[toammann]

Hello gadiLhV,
thanks for looking into this. The model fails in example for kLor = 1/1000 (line 166, calcDjordjevicSarkarApprox.m).

I have changed the value in the attached script below.

One thing that didn't make sense is that you chose the corner\cutoff frequency to be the same as the center frequency. DC calculations are a bit problematic, but let's forget that for now.

I agree, these values come from the MSL_losses example and remained unchanged .

Regards,
Tobias

permittivity_model_energy_inf.zip