CST仿真色散图流程

One of the most important reports in metamaterial topics is “Dispersion curve”, which determines the band-gap behavior of the band-gap structures. In this article two examples that could clearly show the right way of getting this report including two extra Matlab code scripts, are presented.

I hope this article could be helpful for the students who are love Electromagnetic !

“Thanks to Dr Miguel Navarro -Cía Imperial College London for his valuable comment”.

Dispersion Diagram Example 1: Slow Wave Helix General Description

This example demonstrates an eigenmode calculation using periodic boundaries in z-direction and a very meaningful example for showing how we can get the dispersion diagram in CST MWS for periodic structures. This example already is available in example folder of CST MWS in path of:

“C:\Program Files (x86)\CST STUDIO SUITE 2013\Examples\MWS\Eigenmode\TET\Slow Wave\helix.cst ”

Structure description:

Inner conductor (Fig.1(a)) is a helical metal ribbon which has been hold by three dielectric bricks (d) inside a PEC cylindrical hole (b) as outer conductor.

Fig.1 (a) inner conductor

(b) outer conductor

(c) vacuum

(d) dielectric holder

Frequency Range

The frequency range starts at DC and ends at 20 GHz.

Background Material:

Fig.2 Background material setting.

第一个例子

Mesh setting:

Click on the objects in navigation tree and R. Click on the object; then choose the “Local mesh properties” and set it as the following suggestions:

Vacuum: Meshgroup1

The automesh option is disabled for the two inlets which are not parallel to the coordinate axes. Other elements: Meshgroup2 The automesh option is enabled. (Inner&outer conductor, holders)

Fig.3 mesh setting of the model components.

Boundary Conditions

The boundary conditions for x and y planes are set to “electric” or Et=0, which is equivalent to PEC , while for the two boundary conditions in z-direction have been defined as “periodic” in order to model an infinite extent of the proposed helical inner conductor coax unit ell.

A parameter “phase” is assigned to the periodic boundary, so that the phase shift can be used in a parameter sweep, for measuring the frequency dispersive behavior of the present structure.

Fig.4 Boundary conditions and phase shift setting of the one dimensional periodic structure.

*A userdefined watch is already available and may serve as a template for building custom watches. The userdefined parameter sweep watch adds the group velocity, phase velocity and dispersion plots of the first mode to the navigation tree in the folder '1D Results'. Also, the power

flow and Pierce Impedance can be seen in the folder '1D Results'. In parameter sweep window, select "Edit" to view or modify the source code.

>

Fig.5 the user-defined VBA codes of the important reports of the simulation as dispersion

curve, phase & group velocity and Pierce impedance.

But I preferred that the following method for setting the dispersion diagram.

Post Processing:

The dispersion diagram, which is a graph based on “frequency” vs. spatial phase variations of the predefined parameter of “phase”, could obtained by the following step by step process:

1)Recall the TBP and choose the “2D and 3D Field Results>3D Eigenmode Result”

Fig.6 the outline of defining the postprocessing setup of dispersion plot in CST

MWS.

2)For obtaining the dispersion data of the modes start one by one from the first mode, then

after finishing simulation operation export the data in “txt” format(available in “Table”

folder in Navigation tree) for our future postprocessing in Matlab environment.

3)Evaluate

Solver Setup

Whenever the “Eigenmode Solver” is started, a specific number of the lowest resonance frequencies of the structure are calculated. Since only the fundamental mode is of interest, the number of modes is reduced to “1”, and the JDM eigenmode solver is chosen, which is faster for the given example.

Fig.7 Eigenmode Solver setting for calculating the dispersion characteristics of the proposed

periodic structure.

After the parameter sweep has been selected from the Eigenmode Solver's dialog, a new sequence is added, and the parameter "phase" is chosen to be swept from 0 to 180 degrees in 19 steps for an equal steps.

Fig.8 Parametric sweep setup for phase shift parameter of “phase” in the space periodicity

direction.

Start:

If you set lower sweep limit at zero, the machine temporary would interrupt the process and alert we the parameter is zero, don’t care about it and click Ok to continue the operation.

>>

Fig.10

Results

During the run process the real time dispersion results is accessible from the Table folder, as a plot based on frequency vs. phase parameter, as following:

(a)

(b)

Fig. 11 The dispersion diagram based on (a) TBP result template, and (b) VBA

userdefined.

As shown in Fig. 11, all the data are the same but in Fig.11(a) the horizontal axis is the phase.

Export the data of the Mode “1”

>>

>>

Fig. 12 extraction of table dispersion data

Open via excel

Open a new Excel file and open from it your txt file of “Mode1.txt”:

>

>

>>

Fig. 13 The preparation of the dispersion text data file for producing the

compatible *.csv file

We should keep just the phase and value column and eliminate others, then edit the title of the value column as “Frequency(Mode1)”, finally save as a CSV file format, as shown in Fig. 13. We should repeat the simulation for remained modes frequencies as presented in Fig. 14, and extract the favorite data for updating the previous “mode1.csv” data file.

Fig. 15 How we can update the calculation for the higher mode?

Fig. 16 all the dispersion data gathered in one csv file for plotting the final

dispersion diagram for first forth modes of the proposed 1D-periodic structure.

Mtalab Code for plotting the dispersion diagram

For publishing all of your plot result, especially in IEEE publications, we should import our data into Matlab environment, as a perfect engineering tool to sketch various plots in a unique plot in a high resolution view. Here we present a brief quick guide about it.

1)Run Matlab

2) import the dispersion CSV file.

>

4)disable the header and text box.

>>

5)Run the plot m.file(see the appendix)

Final result:

Fig. the final sketch of the dispersion diagram obtained by a MATLAB plotting

program (Appendix).

Phase(Deg)

F r e q u e n c y (

G

H z )

Example2. 2D-periodic structure dispersion diagram

In this section, dispersion characterization of surface waves propagating on a mushroom-like periodic structure is investigated using CST MWS. Mushroom-like EBG structures

Wavenumber k is an important parameter to describe the propagation property of electromagnetic waves. In a lossless case , the phase constant is β = k .

Usually, β is a function of frequency ω. Once the phase constant is obtained, the phase velocity (v p) and group velocity (v g) can be derived:

??=

?

?

??? ??=

????

(2-1)

Furthermore, the field distribution can also be determined, such as the field variation in a

transverse direction. For a plane wave in free space, the relation between β and ω is a linear function:

????=?= ??????

(2-2)

For surface waves propagating in an EBG structure, it is usually difficult to give an explicit expression for the wavenumber k . One has to either solve an eigen-value equation or perform a full wave simulation to determine the wavenumber. It is important to point out that the solution of an eigen-value equation may not be unique.

In another words, there may exist several different propagation constants at the same frequency. Each one is known as a specific mode with its own phase velocity, group velocity, and field distribution. The relation between β and ω is often plotted out and referred to as the dispersion diagram .

For a periodic structure such as the EBG, the field distribution of a surface wave is also periodic with a proper phase delay determined by the wavenumber ? ???β? and periodicity p . Thus, each surface wave mode can be decomposed into an infinite series of space harmonic waves.

?????,?,??=∑??????,???

??β?? ????? , β????? = β? ??? + ? (2-3)

Here, we assume the periodic and propagation direction is the x direction. Although these space

harmonics have different phase velocities, they share the same group velocity. Furthermore, these space harmonics cannot exist individually because each single harmonic does not satisfy the boundary conditions of the periodic structure. Only their summation satisfies the boundary conditions. Thus, they are considered to be the same mode.

Another important observation from (2-3) is that the dispersion curve β? ??? is periodic along the β- axis with a periodicity of 2π/p. Therefore, we only need to plot the dispersion relation within one single period, namely, 0 ≤ β?? ≤ 2?/??, which is known as the Brillouin zone. This concept can be easily extended to two-dimensional periodic structures, where the Brillouin zone becomes a two-dimensional square area:

0 ≤ β?? ≤ 2?/?? , 0 ≤ β?? ≤ 2?/??

Figure 2-1 shows the dispersion diagram of the mushroom-like EBG structure.

The vertical axis shows the frequency and the horizontal axis represents the

values of the transverse wavenumbers (βx, βy). Three specific points are: Γ, X

and M.

Let us consider a plane wave impinging on a periodic surface. For surface waves, the angle of incidence is equal to 90° (parallel with the surface). For this value, the wave propagation on a periodic structure cannot be investigated using the plane wave response. The exact position of pass bands and band gaps in the frequency spectrum can be obtained only by the dispersion relation of surface waves (i.e. calculating the resonant frequencies of eigenmodes) along the contour of the irreducible Brillouin zone.

Brillouin Zone in CST Microwave Studio

In dispersion diagram of a 2D-periodic structure, we are faced with an special definition of the frequency dispersion effects of the surface wave phase variations in three main directions usually labeled as Γ to X, X to M and M to Γ, which we know as “Brillouin Zone”, as shown in Fig..

Fig.2-2. Brillouin zone definition of a 2D-periodic structure (in a full period of

p x=p y=d, we have 2π radian phase change).

Now we should describe three directions phase constant variations with a 2D-plot. In this plot the independent variable of β-axis is a set of ordered pairs of:

??????_?

?????_?

Then, the first set has been called as “Γ to X”:

??????_?:0 ?? 180 ???.

?????_?=0

And the second set, “X to M”:

??????_?=180 ???.

?????_?: 0 ?? 180 ???

Finally, the last set is “M to Γ”:

??????_?:0 ?? 180 ???.

?????_?:0 ?? 180 ???;"???? ??? ????????? ???????? ??? ????????"

Fig.2-3. example of “dispersion diagram” of a 2D-periodic structure, based on

Brillouin Zone definition (in this example the step of sweep is 30deg).

This simulator uses finite integration technique (FIT). The unit cell of the structure under investigation has to be drawn with periodic boundary conditions applied in the appropriate directions. The phase shift is changed along the boundary of the irreducible Brillouin zone and frequencies of eigenmodes are obtained in each step. Because of slow-wave behavior of surface

waves, dispersion curves are calculated only in the region under the light line, as shown in Fig.

2.1. Band gaps occur in frequency intervals, where no dispersion curves in the slow-wave region are present.

Model and Dimension

The most popular mushroom-like EBG structure which is used here as an example of 2D-periodic structures, is composed of a grounded FR4 substrate with periodic square patch above, which the middle of the patch is connected with a via to the ground plane.

Parameter List of our example:

Substrate:(xs,ys) 15×15mm2, FR4(εr=4.5, loss tan. 0.025,

h s=1.6mm)

Patch size: (xs-gap, ys-gap) 14×14mm2 gap=1mm,

Via: radius Rp=0.5mm

Conductor: PEC, thickness(tm):0.035mm

Fig.2-4.

Fig. 2.5 shows the unit cell setup of a periodic structure with the proposed mushroom-like square lattice, for dispersion analysis in CST MWS. The main difference between the computational models in CST MWS and HFSS consists in the following fact. HFSS uses a perfectly matched layer (PML) to represent an infinite air layer above the unit cell, while in CST MWS, open boundaries are not allowed in combination with periodic walls, but only perfect electric conductor (PEC) or perfect magnetic conductor (PMC) boundary conditions can be applied.

(a)

(b)

Fig.2-5. The setup of (a) Background material and (b) six boundaries of the unit cell.

The Eigenmode Solver of CST Microwave Studio does not support open boundaries. In this case, the boundary condition of side walls of the unit cell can be considered periodic.

After many computer simulations with this program and comparing the obtained results with analytical and experimental considerations, a practical rule for the correct surface wave dispersion diagram computation was stated for background material and boundary setting.

Therefore, the key parameter in the computer model is then the height of the air space above of the dielectric substrate to emulate the free space over the structure (in “Background Setting” this space is “Upper Z-distance”), an airbox with the height of about many times of the dielectric slab thickness has to be placed over the unit cell. Based on numerous computer simulations and by comparing the obtained results with analytical and experimental considerations, the correct choice of approximately ten times the substrate thickness was established, “Upper Z-distance ?10·hs”[1].

The boundaries at the top and the bottom of the model should be defined as electric conductor or PEC (Et=0).

If the EBG surfaces are used as ground plane of antenna, like monopole this Figure 3(b) shows the frequency response of transmission coefficient S21, both TM and TE surface waves measured by using a pair of small monopole antennas oriented normally (TM mode) and parallel (TE mode) to the EBG surface.

Frequency Setting: 0-8GHz

Mesh setting:

For not complicated structures like this poject mesh setting is not required but if you want more and more precise results you can pay time for more dense mesh size!

Phase Shift Setting

“Phase_x” and “phase_y” are the phase variation of the slow surface wave on the periodic structure, in x- and y-directions, respectively. We know those parameters have a obvious relation with correspondent wave numbers as the following:

?????_?≡???????? ??????????? ????? ??? x-axis

?????_?≡?????

?????_?≡???????? ??????????? ????? ??? x-axis

?????_?≡?????

For symmetry with respect to the x- and y-axis’s, we define both of xs and ys with a new parameter of “d”.

Γ to X plot:

Sweep the “phase_x (βx d)” from 0 t0 180, with sample: 19 and “phase_y (βy d)=0”

X to M plot:

Sweep the “phase_y (βy d)” from 0 t0 180, with sample: 19 and “phase_x (βx d)=180”

Eigenmode Solver Setting

Start:

Don’t care about alert zero value for phase_x parameter.

Γ to X plot:

Sweep the “phase_x (βx d)” from 0t0 180, with sample: 19and “phase_y (βy d)=0” , the simulation takes not much more time, about 25 minutes by a typical laptop.

The figure shows phase variation linearity below 2.2GHz for Mode#1 and for surface wave propagation phase constant with respect to x-direction.

Export the data and save as “Mode1_Gamma2X.txt”:

The last value of phase_x is 180, that’s perfect! So go to Eigenmode>para sweep and delete phase_x and define phase_y from 0 to 180, then open “Result Template…” and do “Evaluate” for reset for new plot, “X to M”, then click the “Start”:

During the new parameter simulation, still the previous plot would be shown, don’t care about it! The real time new results are available in table format:

Updating the phase_y sweep data:

Recall the table plot, set the phase_y as, “abscissa” and phase_x at 180 as fixed x direction phase which is required for X to M plot,

It’s obvious that the upper limit of Frequency Band Gap region is about 2.13GHz.

Export:

Nav. Tree>R.Click on Table and choose the table. Then change phase_x t0 180 and apply, finally export:

>>

M to Γ Plot:

For convenience to extract this step data we first clean all the stored data of previous sweep by drawing for example a Brick then apply it to message of deletion would be shown, click Ok. In the next step delete that object and now the project is free of any data.

The last plot of Brillouin zone is based on diagonal propagation phase shift which is equivalent to k x=k y or phase_x=phase_y. So we can satisfy this by setting value of phase_y with phase_x, in parameters list window.

>>

Now recall Eigenmode solver and sweep par. And delete the phase_y and define phae_x: 180 to 0deg, with 19 samples, this constraint the calculations performed with descent variation, which is required for M to Γ plot.

Start the Sweep!

Results:

Export the results of the last parametric calculations as “M2Gamma_Mode1.txt”.

Preparing the final Dispersion diagram CSV file

In this section we should create the CSV file format of three correspondent text files of Γto X, X to M and M to Γ, respectively.

Γto X X to M M to Γ

In this step we should prepared a new excel book which contains in a serial format all the three calculated data set, compatible with the Brillouin zone definition, as following:

(a)

(b)

For plotting the desirable Brillouin zone diagram with MATLAB, We insert a new sequential numbering column, as shown in Fig. (b). It is important to determine boarder row numbers of the CSV data file, this issue has shown in the following table.

Plot Name Row range in main CSV data file

Γ to X 1-19 X to M 19-38 M to Γ

38-57