Google ColabのJulia言語で1次元シュレーディンガー方程式の井戸型ポテンシャルによるトンネル効果を計算してアニメーションにしてみる。空間方向は普通の差分、時間は自前の13段8次のルンゲクッタ法を使用。

前回は波束の散乱だったが、今回はトンネル効果。と言ってもやったことは井戸型ポテンシャルの符号を変えただけ。

コードも省略でさっそく結果から。

高いポテンシャル障壁でも向こう側に染み出ている。

Google ColabのJulia言語で1次元シュレーディンガー方程式の井戸型ポテンシャルによる波束の散乱を計算してアニメーションにしてみる。空間方向は普通の差分、時間は自前の13段8次のルンゲクッタ法を使用。

今回はシュレーディンガー方程式をやってみる。

昔Excel VBAで計算して、Pythonでもやってみた1次元シュレーディンガー方程式のポテンシャルによる波束の散乱を計算しよう。
これはシッフの量子力学の例題で出ていたものを「計算物理」という本でfortranで書かれていたものの移植。

ただ、空間微分は普通に(φ(i+1) + φ(i-1) - 2 φ(i))/h²で差分化するが、時間方向は13段8次のルンゲクッタ法を使う。DifferentialEquations.jlを使えば一発だが、Google Colabで使うときは毎回Pkg.addとしないといけなくて、それが45分かかるのでもうやだ、ということで自前のもの（Excel VBA移植版）。

コードはこんな感じで。ルンゲクッタ部分の係数が長い…

using Plots

using Printf

function main()

    n = 2000 #空間分割数

    m = 5500 #時間分割数

    #波動関数の実部u、虚部v

    u = zeros(Float64, n + 2)

    v = zeros(Float64, n + 2)

    ud = zeros(Float64, n + 2)

    vd = zeros(Float64, n + 2)

    #ポテンシャル

    u_pot = zeros(Float64, n + 2)

    #物理パラメータ

    v_width = 0.064

    v0 = 0.6 * (70.7 * pi) ^ 2

    deltax = 0.035

    x0 = -0.3

    p0 = pi * 50.0

    t = 0.0

    dx = 0.001

    dt = dx * dx / 2.0

    xmax = dx * (n - 1) / 2.0

    xmin = -dx * (n - 1) / 2.0

    x = zeros(Float64, n + 2)

    #表示用

    results = []

    ndiv = 25

    #8次のルンゲクッタ法の係数

    a = zeros(Float64, 13)

    b = zeros(Float64, 13, 13)

    c = zeros(Float64, 13)

    ku = zeros(Float64, n + 1, 13)

    kv = zeros(Float64, n + 1, 13)

    a[1] = 0.0

    a[2] = 2.0 / 27.0

    a[3] = 1.0 / 9.0

    a[4] = 1.0 / 6.0

    a[5] = 5.0 / 12.0

    a[6] = 0.5

    a[7] = 5.0 / 6.0

    a[8] = 1.0 / 6.0

    a[9] = 2.0 / 3.0

    a[10] = 1.0 / 3.0

    a[11] = 1.0

    a[12] = 0.0

    a[13] = 1

    b[2, 1] = 2.0 / 27.0

    b[3, 1] = 1.0 / 36

    b[3, 2] = 1.0 / 12.0

    b[4, 1] = 1.0 / 24.0

    b[4, 3] = 1.0 / 8.0

    b[5, 1] = 5.0 / 12.0

    b[5, 3] = -25.0 / 16.0

    b[5, 4] = 25.0 / 16.0

    b[6, 1] = 1.0 / 20.0

    b[6, 4] = 1.0 / 4.0

    b[6, 5] = 1.0 / 5.0

    b[7, 1] = -25.0 / 108.0

    b[7, 4] = 125.0 / 108.0

    b[7, 5] = -65.0 / 27.0

    b[7, 6] = 125.0 / 54.0

    b[8, 1] = 31.0 / 300.0

    b[8, 5] = 61.0 / 225.0

    b[8, 6] = -2.0 / 9.0

    b[8, 7] = 13.0 / 900.0

    b[9, 1] = 2.0

    b[9, 4] = -53.0 / 6.0

    b[9, 5] = 704.0 / 45.0

    b[9, 6] = -107.0 / 9.0

    b[9, 7] = 67.0 / 90.0

    b[9, 8] = 3.0

    b[10, 1] = -91.0 / 108.0

    b[10, 4] = 23.0 / 108.0

    b[10, 5] = -976.0 / 135.0

    b[10, 6] = 311.0 / 54.0

    b[10, 7] = -19.0 / 60.0

    b[10, 8] = 17.0 / 6.0

    b[10, 9] = -1.0 / 12.0

    b[11, 1] = 2383.0 / 4100.0

    b[11, 4] = -341.0 / 164.0

    b[11, 5] = 4496.0 / 1025.0

    b[11, 6] = -301.0 / 82.0

    b[11, 7] = 2133.0 / 4100.0

    b[11, 8] = 45.0 / 82.0

    b[11, 9] = 45.0 / 164.0

    b[11, 10] = 18.0 / 41.0

    b[12, 1] = 3.0 / 205.0

    b[12, 6] = -6.0 / 41.0

    b[12, 7] = -3.0 / 205.0

    b[12, 8] = -3.0 / 41.0

    b[12, 9] = 3.0 / 41.0

    b[12, 10] = 6.0 / 41.0

    b[13, 1] = -1777.0 / 4100.0

    b[13, 4] = -341.0 / 164.0

    b[13, 5] = 4496.0 / 1025.0

    b[13, 6] = -289.0 / 82.0

    b[13, 7] = 2193.0 / 4100.0

    b[13, 8] = 51.0 / 82.0

    b[13, 9] = 33.0 / 164.0

    b[13, 10] = 12.0 / 41.0

    b[13, 12] = 1.0

    c[6] = 34.0 / 105.0

    c[7] = 9.0 / 35.0

    c[8] = 9.0 / 35.0

    c[9] = 9.0 / 280.0

    c[10] = 9.0 / 280.0

    c[12] = 41.0 / 840.0

    c[13] = 41.0 / 840.0

    #ポテンシャルと波動関数の初期値（ガウシアン波束）設定

    for i in 2:(n + 1)

        x[i] = xmin + (i - 1) * dx

        if x[i] >= -v_width / 2.0 && x[i] <= v_width / 2.0

            u_pot[i] = v0

        else

            u_pot[i] = 0.0

        end

        u[i] = exp(-((x[i] - x0) ^ 2) / (4 * deltax * deltax)) * cos(p0 * x[i]) / ((2 * pi * deltax * deltax) ^ 0.25)

        v[i] = exp(-((x[i] - x0) ^ 2) / (4 * deltax * deltax)) * sin(p0 * x[i]) / ((2 * pi * deltax * deltax) ^ 0.25)

    end

    u[1] = u[2]

    u[n + 2] = u[n + 1]

    v[1] = v[2]

    v[n + 2] = v[n + 1]

    x[1] = xmin - dx

    x[n + 2] = xmax + dx

    for i in 1:m

        # 8次のルンゲクッタ法計算

         @inbounds for ii in 2:(n + 1)

            ku[ii - 1, 1] = f(ii, t + dt, u, v, dx, u_pot, n) * dt

            kv[ii - 1, 1] = g(ii, t + dt, u, v, dx, u_pot, n) * dt

        end

        @inbounds for j in 2:13

            @inbounds @simd for ii in 2:(n + 1)

                ud[ii] = u[ii]

                vd[ii] = v[ii]

            end

            @inbounds for k in 1:(j - 1)

                @inbounds @simd for ii in 2:(n + 1)

                    ud[ii] = ud[ii] + b[j, k] * ku[ii - 1, k]

                    vd[ii] = vd[ii] + b[j, k] * kv[ii - 1, k]

                end

            end

            @inbounds for ii in 2:(n + 1)

                ku[ii - 1, j] = f(ii, t + a[j] * dt, ud, vd, dx, u_pot, n) * dt

                kv[ii - 1, j] = g(ii, t + a[j] * dt, ud, vd, dx, u_pot, n) * dt

            end

        end

        @inbounds @simd for j in 1:13

            for ii in 2:(n + 1)

                u[ii] = u[ii] + c[j] * ku[ii - 1, j]

                v[ii] = v[ii] + c[j] * kv[ii - 1, j]

            end

        end

        t = t + dt

        # ndivステップごとに結果(波動関数の大きさ)を配列に格納

        if i % ndiv == 0

            push!(results, sqrt.(u.^2 + v.^2))

        end

    end

    # 計算結果をアニメーションで表示

    anim = @animate for i in 1:length(results)

        plot(x, results[i], ylim = (0, 4), linewidth = 5,title="Tunnelling Effect", label="Wave function", xlabel="x", ylabel="phi", size=(900,500))

        plot!(x, u_pot, label = "Potential")

    end

    gif(anim, "Tunnel.gif", fps = 30)    

end

function f(i, t, u, v, dx, u_pot, n)

    v[1] = v[2]

    v[n + 2] = v[n + 1]

    d2 = (v[i + 1] + v[i - 1] - 2.0 * v[i]) / (dx * dx)

    return -d2 + u_pot[i] * v[i]

end

function g(i, t, u, v, dx, u_pot, n)

    u[1] = u[2]

    u[n + 2] = u[n + 1]

    d2 = (u[i + 1] + u[i - 1] - 2.0 * u[i]) / (dx * dx)

    return d2 - u_pot[i] * u[i]

end

実行結果。一部は通過して一部は反射している。

次はトンネル効果やってみよう。

PythonでFDTD法で電磁界シミュレーションできるopenEMSを使う(4)円筒座標のメッシュが使えるので円筒型のコンフォーマルパッチアンテナの計算をして電流分布を動画に、指向性などをグラフにする。

今回は例題に乗っていた曲がったパッチアンテナ。

https://docs.openems.de/python/openEMS/Tutorials/Bent_Patch_Antenna.html

https://wiki.openems.de/index.php/Tutorial:_Bent_Patch_Antenna.html

まずは必要なモジュールをインポート。

import os

import numpy as np

import matplotlib.pyplot as plt

from CSXCAD import CSXCAD

from openEMS.openEMS import openEMS

from openEMS.physical_constants import C0, EPS0

寸法や周波数の設定を行う。unitはmm単位にする。

Sim_Path = os.path.join("C:\\Users\\tomoh\\Documents\\Python Projects\\OpenEMS", 'Bent_Patch')

post_proc_only = False

unit = 1e-3 # all length in mm

f0 = 2.4e9 # center frequency, frequency of interest!

lambda0 = round(C0/f0/unit) # wavelength in mm

fc = 0.5e9 # 20 dB corner frequency

# patch width in alpha-direction

patch_width  = 32 # resonant length in alpha-direction

patch_radius = 50 # radius

patch_length = 40 # patch length in z-direction

#substrate setup

substrate_epsR   = 3.38

substrate_kappa  = 1e-3 * 2*np.pi*2.45e9 * EPS0*substrate_epsR

substrate_width  = 80

substrate_length = 90

substrate_thickness = 1.524

substrate_cells = 4

#setup feeding

feed_pos   = -5.5  #feeding position in x-direction

feed_width = 2     #feeding port width

feed_R     = 50    #feed resistance

# size of the simulation box

SimBox_rad    = 2*100

SimBox_height = 1.5*200

FDTDの設定。CoordSystem=1が円筒座標にする設定で、境界条件はすべてMurの吸収境界条件。

FDTD = openEMS(CoordSystem=1, EndCriteria=1e-4) # init a cylindrical FDTD

f0 = 2e9 # center frequency

fc = 1e9 # 20 dB corner frequency

FDTD.SetGaussExcite(f0, fc)

FDTD.SetBoundaryCond(['MUR', 'MUR', 'MUR', 'MUR', 'MUR', 'MUR']) # boundary conditions

# init a cylindrical mesh

CSX = CSXCAD.ContinuousStructure(CoordSystem=1)

FDTD.SetCSX(CSX)

mesh = CSX.GetGrid()

mesh.SetDeltaUnit(unit)

モデル作成。励振はLumpedポート。

# calculate some width as an angle in radiant

patch_ang_width = patch_width/(patch_radius+substrate_thickness)

substr_ang_width = substrate_width/patch_radius

feed_angle = feed_pos/patch_radius

# create patch

patch = CSX.AddMetal('patch') # create a perfect electric conductor (PEC)

start = [patch_radius+substrate_thickness, -patch_ang_width/2, -patch_length/2 ]

stop  = [patch_radius+substrate_thickness,  patch_ang_width/2,  patch_length/2 ]

patch.AddBox(priority=10, start=start, stop=stop) # add a box-primitive to the metal property 'patch'

FDTD.AddEdges2Grid(dirs='all', properties=patch)

# create substrate

substrate = CSX.AddMaterial('substrate', epsilon=substrate_epsR, kappa=substrate_kappa  )

start = [patch_radius                    , -substr_ang_width/2, -substrate_length/2]

stop  = [patch_radius+substrate_thickness,  substr_ang_width/2,  substrate_length/2]

substrate.AddBox(start=start, stop=stop)

FDTD.AddEdges2Grid(dirs='all', properties=substrate)

# save current density oon the patch

jt_patch = CSX.AddDump('Jt_patch', dump_type=3, file_type=0)

start = [patch_radius+substrate_thickness, -substr_ang_width/2, -substrate_length/2]

stop  = [patch_radius+substrate_thickness, +substr_ang_width/2,  substrate_length/2]

jt_patch.AddBox(start=start, stop=stop)

# create ground

gnd = CSX.AddMetal('gnd') # create a perfect electric conductor (PEC)

start = [patch_radius, -substr_ang_width/2, -substrate_length/2]

stop  = [patch_radius, +substr_ang_width/2, +substrate_length/2]

gnd.AddBox(priority=10, start=start, stop=stop)

FDTD.AddEdges2Grid(dirs='all', properties=gnd)

# apply the excitation & resist as a current source

start = [patch_radius                    ,  feed_angle, 0]

stop  = [patch_radius+substrate_thickness,  feed_angle, 0]

port = FDTD.AddLumpedPort(1 ,feed_R, start, stop, 'r', 1.0, priority=50, edges2grid='all')

メッシュの設定と、

# add the simulation domain size

mesh.AddLine('r', patch_radius+np.array([-20, SimBox_rad]))

mesh.AddLine('a', [-0.75*np.pi, 0.75*np.pi])

mesh.AddLine('z', [-SimBox_height/2, SimBox_height/2])

# add some lines for the substrate

mesh.AddLine('r', patch_radius+np.linspace(0,substrate_thickness,substrate_cells))

# generate a smooth mesh with max. cell size: lambda_min / 20

max_res = C0 / (f0+fc) / unit / 20

max_ang = max_res/(SimBox_rad+patch_radius) # max res in radiant

mesh.SmoothMeshLines(0, max_res, 1.4)

mesh.SmoothMeshLines(1, max_ang, 1.4)

mesh.SmoothMeshLines(2, max_res, 1.4)

近傍界から遠方界に変換する準備。

nf2ff = FDTD.CreateNF2FFBox()

そして計算。ここでモデルも確認しておく。

if 1:  # debugging only

    CSX_file = os.path.join(Sim_Path, 'bent_patch.xml')

    if not os.path.exists(Sim_Path):

        os.mkdir(Sim_Path)

    CSX.Write2XML(CSX_file)

    from CSXCAD import AppCSXCAD_BIN

    os.system(AppCSXCAD_BIN + ' "{}"'.format(CSX_file))

if not post_proc_only:

    FDTD.Run(Sim_Path, cleanup=True)

そしてSパラメータ（S11）、入力インピーダンス、指向性の図示。

f = np.linspace(max(1e9,f0-fc),f0+fc,401)

port.CalcPort(Sim_Path, f)

Zin = port.uf_tot / port.if_tot

s11 = port.uf_ref/port.uf_inc

s11_dB = 20.0*np.log10(np.abs(s11))

plt.figure()

plt.plot(f/1e9, s11_dB)

plt.grid()

plt.ylabel('s11 (dB)')

plt.xlabel('frequency (GHz)')

P_in = 0.5*np.real(port.uf_tot * np.conj(port.if_tot)) # antenna feed power

# plot feed point impedance

plt.figure()

plt.plot( f/1e6, np.real(Zin), 'k-', linewidth=2, label=r'$\Re(Z_{in})$' )

plt.grid()

plt.plot( f/1e6, np.imag(Zin), 'r--', linewidth=2, label=r'$\Im(Z_{in})$' )

plt.title( 'feed point impedance' )

plt.xlabel( 'frequency (MHz)' )

plt.ylabel( 'impedance ($\Omega$)' )

plt.legend( )

idx = np.where((s11_dB<-10) & (s11_dB==np.min(s11_dB)))[0]

if not len(idx)==1:

    print('No resonance frequency found for far-field calulation')

else:

    f_res = f[idx[0]]

    theta = np.arange(-180.0, 180.0, 2.0)

    print("Calculate NF2FF")

    nf2ff_res_phi0 = nf2ff.CalcNF2FF(Sim_Path, f_res, theta, 0, center=np.array([patch_radius+substrate_thickness, 0, 0])*unit, read_cached=True, outfile='nf2ff_xz.h5')

    plt.figure(figsize=(15, 7))

    ax = plt.subplot(121, polar=True)

    E_norm = 20.0*np.log10(nf2ff_res_phi0.E_norm/np.max(nf2ff_res_phi0.E_norm)) + nf2ff_res_phi0.Dmax

    ax.plot(np.deg2rad(theta), 10**(np.squeeze(E_norm)/20), linewidth=2, label='xz-plane')

    ax.grid(True)

    ax.set_xlabel('theta (deg)')

    ax.set_theta_zero_location('N')

    ax.set_theta_direction(-1)

    ax.legend(loc=3)

    phi = theta

    nf2ff_res_theta90 = nf2ff.CalcNF2FF(Sim_Path, f_res, 90, phi, center=np.array([patch_radius+substrate_thickness, 0, 0])*unit, read_cached=True, outfile='nf2ff_xy.h5')

    ax = plt.subplot(122, polar=True)

    E_norm = 20.0*np.log10(nf2ff_res_theta90.E_norm/np.max(nf2ff_res_theta90.E_norm)) + nf2ff_res_theta90.Dmax

    ax.plot(np.deg2rad(phi), 10**(np.squeeze(E_norm)/20), linewidth=2, label='xy-plane')

    ax.grid(True)

    ax.set_xlabel('phi (deg)')

    plt.suptitle('Bent Patch Anteanna Pattern\nFrequency: {} GHz'.format(f_res/1e9), fontsize=14)

    ax.legend(loc=3)

    print( 'radiated power: Prad = {:.2e} Watt'.format(nf2ff_res_theta90.Prad[0]))

    print( 'directivity:    Dmax = {:.1f} ({:.1f} dBi)'.format(nf2ff_res_theta90.Dmax[0], 10*np.log10(nf2ff_res_theta90.Dmax[0])))

    print( 'efficiency:   nu_rad = {:.1f} %'.format(100*nf2ff_res_theta90.Prad[0]/np.real(P_in[idx[0]])))

plt.show()

Calculate NF2FF
radiated power: Prad = 4.94e-26 Watt
directivity: Dmax = 5.5 (7.4 dBi)
efficiency: nu_rad = 89.9 %

となった。電流分布も動画にするとこうなる。

ものすごく時間がかかりデータも膨大になるが、

Et = CSX.AddDump('Et', dump_type=0, file_type=0)

Et.AddBox([patch_radius-20,-0.75*np.pi,-SimBox_height/2], [patch_radius+SimBox_rad,0.75*np.pi,SimBox_height/2])

を使えると電界分布も動画にできる。

Google ColabのJulia言語で2次元Swift-Hohenberg方程式（∂φ/∂t=φ - φ³ - (∇²+ k₀²)²φ、熱対流などを表す）を差分法で計算してパターンを動画にしてみる。

前回は複素TDGL方程式をやったので、今回は類似のSwift-Hohenberg方程式

∂φ/∂t = φ - φ³ - (∇²+ k₀²)²φ

をやってみよう。空間の4階微分が出てくる。コードはこんな感じで。

TDGLのプログラム流用なので無駄があり変数名がおかしいがまあ気にしないで。

using Plots

using Printf

using Random

function main()

    #パラメータ設定

    n = 128

    m = 1000

    ndiv = 2

    dt = 0.02

    dx = 1.0

    C1 = 1.0

    #初期設定

    X1 = zeros(n + 2, n + 2)

    Y1 = zeros(n + 2, n + 2)

    X2 = zeros(n + 2, n + 2)

    MT = MersenneTwister()

    Random.seed!(MT, 42)

    @inbounds for j in 2:(n + 1)

        @inbounds @simd for i in 2:(n + 1)

            X1[i, j] = 0.1 * (2.0 * rand(MT) - 1.0)

        end

    end

    # 結果を格納する配列

    results = []

    for t in 1:m

        #境界条件

        @inbounds for i in 2:(n + 1)

            X1[1, i] = X1[2, i]

            X1[n + 2, i] = X1[n + 1, i]

            X1[i, 1] = X1[i, 2]

            X1[i, n + 2] = X1[i, n + 1]

        end

        #一つ目のラプラシアン

        @inbounds for j in 2:(n + 1)

            @inbounds @simd for i in 2:(n + 1)

                lapX = (X1[i + 1, j] + X1[i - 1, j] + X1[i, j + 1] + X1[i, j - 1] - 4.0 * X1[i, j]) / (dx * dx)

                Y1[i, j] = lapX + C1 * C1 * X1[i, j]                                                                              

            end

        end

        #境界条件

        @inbounds for i in 2:(n + 1)    

            Y1[1, i] = Y1[2, i]

            Y1[n + 2, i] = Y1[n + 1, i]

            Y1[i, 1] = Y1[i, 2]

            Y1[i, n + 2] = Y1[i, n + 1]

        end

        #Swift-Hohenberg計算

        @inbounds for j in 2:(n + 1)

            @inbounds @simd for i in 2:(n + 1)

                lapY = (Y1[i + 1, j] + Y1[i - 1, j] + Y1[i, j + 1] + Y1[i, j - 1] - 4.0 * Y1[i, j]) / (dx * dx)

                X2[i, j] = X1[i, j] + dt * (X1[i, j] * (1.0 - X1[i, j] ^ 2) - (lapY + C1 * C1 * Y1[i, j]))                                                                            

            end

        end

        @inbounds for j in 2:(n + 1)

            @inbounds @simd for i in 2:(n + 1)

                X1[i, j] = X2[i, j]

            end

        end

        # ndivステップごとに結果を配列に格納

        if t % ndiv == 0

            push!(results, copy(X1))

        end

    end

    # 計算結果をアニメーションで表示

    anim = @animate for i in 1:length(results)

        heatmap(results[i], title="t = $(@sprintf("%.3f", i * ndiv * dt))", clim=(-1.2, 1.2), size = (950, 800))

    end

    gif(anim, "Swift-Hohenberg.gif", fps = 10)    

end

これでmain()とすると動画ができる。

いつの間にかWindows11のメモ帳がMarkdown記法が使えるようになっていた。使えるのは見出し、強調、斜体、リスト（数字付き、なし）、リンク。書式設定にすることもできる。Copilotもついてもう何が何だか…動画で動作を紹介。

この記事見た。

メモ帳がMarkdownに対応、Windows 11 Insider向けに順次展開

私もWindows11 Insiderに加入しているし、どうかなと思ったら使えるようになっていた。

動画で動作を紹介。

Copilotで要約やリライトもできるし、もう一体なにがメモ帳なのかわからない…

それの対局としてEditも入れてみている。

Microsoftの新しいCLIテキストエディター「Edit」に初めての更新、日本語メニューも完備

こっちのほうが落ち着く…

Google ColabのJulia言語で2次元複素TDGL方程式(Complex Time-Dependent Ginzburg Landau equation, ∂W/∂t = (1+iC₀)W + (1+iC₁)∇²W - (1+iC₂)|W|²W)を計算してスパイラルパターンをアニメーションにしてみる。

2次元複素TDGL方程式(Complex Time-Dependent Ginzburg Landau equation)は昔、Excel VBAで計算したり遊んでいた。

こんな形の方程式

∂W/∂t = (1+iC₀)W + (1+iC₁)∇²W - (1+iC₂)|W|²W

相転移に使われるGL方程式を複素数に拡張したもの。パラメータを選ぶとスパイラルパターン（らせん）が出ることはよく知られている。

今回はJuliaでやってみよう。VBAを移植するのはとても簡単。

コードはこんな感じで。

using Plots

using Printf

using Random

function main()

    #パラメータ設定

    n = 128

    m = 10000

    ndiv = 10

    dt = 0.1

    dx = 1.0

    C0 = 0.0

    C1 = 0.0

    C2 = 1.0

    #初期設定

    X1 = zeros(n + 2, n + 2)

    Y1 = zeros(n + 2, n + 2)

    X2 = zeros(n + 2, n + 2)

    Y2 = zeros(n + 2, n + 2)

    MT = MersenneTwister()

    Random.seed!(MT, 42)

    @inbounds for i in (Int(round(n * 0.4)) + 1):(Int(round(n * 0.6)) + 1)

        @simd for j in (Int(round(n * 0.4)) + 1):(Int(round(n * 0.6)) + 1)

            X1[i, j] = 0.1 * (2.0 * rand(MT) - 1.0)

            Y1[i, j] = 0.1 * (2.0 * rand(MT) - 1.0)

        end

    end

    # 結果を格納する配列

    results = []

    for t in 1:m

        # 境界条件

        @inbounds for i in 2:(n + 1)

            X1[1, i] = X1[2, i]

            X1[n + 2, i] = X1[n + 1, i]

            X1[i, 1] = X1[i, 2]

            X1[i, n + 2] = X1[i, n + 1]

            Y1[1, i] = Y1[2, i]

            Y1[n + 2, i] = Y1[n + 1, i]

            Y1[i, 1] = Y1[i, 2]

            Y1[i, n + 2] = Y1[i, n + 1]

        end

        # TDGLの計算

        @inbounds for i in 2:(n + 1)

            @simd for j in 2:(n + 1)

                R2 = X1[i, j] * X1[i, j] + Y1[i, j] * Y1[i, j]

                lapX = (X1[i + 1, j] + X1[i - 1, j] + X1[i, j + 1] + X1[i, j - 1] - 4.0 * X1[i, j]) / (dx * dx)

                lapY = (Y1[i + 1, j] + Y1[i - 1, j] + Y1[i, j + 1] + Y1[i, j - 1] - 4.0 * Y1[i, j]) / (dx * dx)

                X2[i, j] = X1[i, j] + dt * ((X1[i, j] - C0 * Y1[i, j]) + (lapX - C1 * lapY) - R2 * (X1[i, j] - C2 * Y1[i, j]))                                    

                Y2[i, j] = Y1[i, j] + dt * ((C0 * X1[i, j] + Y1[i, j]) + (C1 * lapX + lapY) - R2 * (C2 * X1[i, j] + Y1[i, j]))                                                                                

            end

        end

        @inbounds for i in 2:(n + 1)

            @simd for j in 2:(n + 1)

                X1[i, j] = X2[i, j]

                Y1[i, j] = Y2[i, j]

            end

        end

        # ndivステップごとに結果を配列に格納

        if t % ndiv == 0

            push!(results, copy(X1))

        end

    end

    # 計算結果をアニメーションで表示、heatmapのスケールを-1から1に固定

    anim = @animate for i in 300:length(results)

        heatmap(results[i], title="t = $(@sprintf("%.3f", i * ndiv * dt))", clim=(-1.0, 1.0), size = (950, 800))

    end

    gif(anim, "TDGLequation.gif", fps = 10)    

end

main()

とすれば数分で計算できる。

結果はこちら。

らせんを描いている！

PythonでFDTD法で電磁界シミュレーションできるopenEMSを使う(3)TE01モードの矩形導波管を計算してSパラメータ表示と電界分布を動画にする。

マイクロストリップラインノッチフィルタ、マイクロストリップラインパッチアンテナと来て今度は矩形導波管。

https://docs.openems.de/python/openEMS/Tutorials/Rect_Waveguide.html

https://wiki.openems.de/index.php/Tutorial:_Rectangular_Waveguide.html

まずはモジュールをインポート。

import os

import numpy as np

import matplotlib.pyplot as plt

from CSXCAD  import ContinuousStructure

from openEMS import openEMS

from openEMS.physical_constants import C0

次に寸法や周波数を設定。

Sim_Path = os.path.join("C:\\Users\\tomoh\\Documents\\Python Projects\\OpenEMS", 'Rect_WG')

post_proc_only = False

unit = 1e-6 #drawing unit in um

# waveguide dimensions

# WR42

a = 10700   #waveguide width

b = 4300    #waveguide height

length = 50000

# frequency range of interest

f_start = 20e9

f_0     = 24e9

f_stop  = 26e9

lambda0 = C0/f_0/unit

#waveguide TE-mode definition

TE_mode = 'TE10'

#targeted mesh resolution

mesh_res = lambda0/30

そして励振と境界条件を設定。0がPECで導波管の壁になる。3はPMLで終端に相当。

   0 = PEC      or  'PEC'

   1 = PMC      or  'PMC'

   2 = MUR-ABC  or  'MUR'

   3 = PML-ABC  or  'PML_x' with pml size x => 4..50

FDTD = openEMS(NrTS=1e4)

FDTD.SetGaussExcite(0.5*(f_start+f_stop),0.5*(f_stop-f_start))

# boundary conditions

FDTD.SetBoundaryCond([0, 0, 0, 0, 3, 3])

モデルを作る。導波管は直方体書くだけなので簡単。

CSX = ContinuousStructure()

FDTD.SetCSX(CSX)

mesh = CSX.GetGrid()

mesh.SetDeltaUnit(unit)

mesh.AddLine('x', [0, a])

mesh.AddLine('y', [0, b])

mesh.AddLine('z', [0, length])

次はポート設定。矩形導波管専用のポートがある。

ports = []

start=[0, 0, 10*mesh_res]

stop =[a, b, 15*mesh_res]

mesh.AddLine('z', [start[2], stop[2]])

ports.append(FDTD.AddRectWaveGuidePort( 0, start, stop, 'z', a*unit, b*unit, TE_mode, 1))

start=[0, 0, length-10*mesh_res]

stop =[a, b, length-15*mesh_res]

mesh.AddLine('z', [start[2], stop[2]])

ports.append(FDTD.AddRectWaveGuidePort( 1, start, stop, 'z', a*unit, b*unit, TE_mode))

mesh.SmoothMeshLines('all', mesh_res, ratio=1.4)

ここで電界分布をみるための設定をしておく。設定はこう。

• DumpType:
  
  • 0: for E-field time-domain dump (default)
  • 1: for H-field time-domain dump
  • 2: for electric current time-domain dump
  • 3: for total current density (rot(H)) time-domain dump
  • 10: for E-field frequency-domain dump
  • 11: for H-field frequency-domain dump
  • 12: for electric current frequency-domain dump
  • 13: for total current density (rot(H)) frequency-domain dump
  • 20: local SAR frequency-domain dump
  • 21: 1g averaging SAR frequency-domain dump
  • 22: 10g averaging SAR frequency-domain dump
  • 29: raw data needed for SAR calculations (electric field FD, cell volume, conductivity and density)

Et = CSX.AddDump('Et', file_type=0, sub_sampling=[2,2,2])

start = [0, 0, 0]

stop  = [a, b, length]

Et.AddBox(start, stop)

そして計算。一応モデルも確認。色がついているのは励振部分。

if 1:  # debugging only

    CSX_file = os.path.join(Sim_Path, 'rect_wg.xml')

    if not os.path.exists(Sim_Path):

        os.mkdir(Sim_Path)

    CSX.Write2XML(CSX_file)

    from CSXCAD import AppCSXCAD_BIN

    os.system(AppCSXCAD_BIN + ' "{}"'.format(CSX_file))

if not post_proc_only:

    FDTD.Run(Sim_Path, cleanup=False)

Sパラメータの図示。

freq = np.linspace(f_start,f_stop,201)

for port in ports:

    port.CalcPort(Sim_Path, freq)

s11 = ports[0].uf_ref / ports[0].uf_inc

s21 = ports[1].uf_ref / ports[0].uf_inc

ZL  = ports[0].uf_tot / ports[0].if_tot

ZL_a = ports[0].ZL # analytic waveguide impedance

plt.figure()

plt.plot(freq*1e-6,20*np.log10(abs(s11)),'k-',linewidth=2, label='$S_{11}$')

plt.grid()

plt.plot(freq*1e-6,20*np.log10(abs(s21)),'r--',linewidth=2, label='$S_{21}$')

plt.legend()

plt.ylabel('S-Parameter (dB)')

plt.xlabel(r'frequency (MHz) $\rightarrow$')

なるほど計算できてそうだ。PECなんでロスなしなのがいまいちだが。

電界分布もParaViewで見る。面でスライスしないと見えなかった。

Google ColabのJulia言語でFPUT問題(Fermi–Pasta–Ulam–Tsingou、非線形結合した振動子が最初に与えたモードに戻る再帰現象)をDifferentialEquations.jlの2階の常微分方程式ソルバーDynamicalODEProblemでシンプレクティック8次のKahanLi8で計算、振動子の動きも動画にしてみる。

FPU問題というのがある。最近ではFPUT問題（Fermi–Pasta–Ulam–Tsingou）ともいわれる（TはManiacでコードを書いた女性プログラマ）。

https://en.wikipedia.org/wiki/Fermi%E2%80%93Pasta%E2%80%93Ulam%E2%80%93Tsingou_problem

非線形結合した振動子が最初に与えたモードに戻る再帰現象を発見した。

これを確かめてみよう。

まず、必要なパッケージをインストールする。が45分かかった…ここだけがGoogle ColabでJuliaをやるときのネック（Pythonは一瞬）。

using Pkg

Pkg.add("DifferentialEquations")

Pkg.add("RecursiveArrayTools")

DifferentialEquations.jlには2階の常微分方程式専用のソルバーがある。

Dynamical, Hamiltonian, and 2nd Order ODE Solvers

DynamicalODEProblem{isinplace}(f1, f2, v0, u0, tspan, p = NullParameters(); kwargs...)
SecondOrderODEProblem{isinplace}(f, du0, u0, tspan, p = NullParameters(); kwargs...)
HamiltonianProblem{T}(H, p0, q0, tspan, p = NullParameters(); kwargs...)

おお、HamiltonianProblem使えば一瞬じゃないか、と思ったら…そんな関数はないとエラー、じゃあSecondOrderODEProblemは…これもエラー、DynamicalODEProblemは…これは使える。

まあシンプレクティック積分を使うだろうな、と思ったらいろいろある。その中で

”For symplectic methods, higher order algorithms are the most efficient when higher accuracy is needed, and when less accuracy is needed lower order methods do better. Optimized efficiency methods take more steps and thus have more force calculations for the same order, but have smaller error. Thus, the “optimized efficiency” algorithms are recommended if your force calculation is not too sufficiency large, while the other methods are recommended when force calculations are really large (for example, like in MD simulations VelocityVerlet is very popular since it only requires one force calculation per timestep). A good go-to method would be McAte5, and a good high order choice is KahanLi8.
”

ということなのでシンプレクティック8次のKahanLi8を使う。

コードはこんな感じで。RecursiveArrayToolsを使うとsolverから望むデータが取り出せる（とGeminiに教えてもらった）。

using DifferentialEquations

using Plots

using RecursiveArrayTools

using Printf

n = 32 #振動子の数

tmax = 12000.0 #最大時間

dt = 0.2 #時間刻み幅

function Velocity!(dx, x, v, p, t)

    dx .= v

end

#　非線形な結合を設定

function Force!(dv, x, v, p, t)

    dv[1] = 0.0

    dv[n + 2] = 0.0

    for i in 2:(n+1)

        dv[i] = (x[i + 1] + x[i - 1] - 2.0 * x[i]) + 0.25 * ((x[i + 1] - x[i])^2 - (x[i] - x[i - 1])^2)

    end

end

# 初期条件

x0 = zeros(Float64, n + 2)

v0 = zeros(Float64, n + 2)

for i in 2:(n+1)

    x0[i] = sin(pi * (i - 1) / (n + 1))

end

# 2階常微分方程式をシンプレクティック8次で計算

prob = DynamicalODEProblem(Velocity!, Force!, x0, v0, (0.0, tmax))

sol = solve(prob, KahanLi8(), dt = dt)

#モードごとのエネルギーを計算

m = length(sol.t)

a = zeros(Float64, n)

ad = zeros(Float64, n)

E = zeros(Float64, n, m)

for i in 1:m

    for k in 1:n

        a[k] = 0.0

        ad[k] = 0.0

        for j in 2:(n+1)

            a[k] += sol[i].x[1][j] * sin((j - 1) * k * pi / (n + 1))

            ad[k] += sol[i].x[2][j] * sin((j - 1) * k * pi / (n + 1))

        end

        E[k, i] = sqrt(2 / (n + 1)) * ((ad[k] ^ 2) / 2 + 2 * (a[k] ^ 2) * (sin(pi * k / 2 / (n + 1)) ^ 2))

    end

end

#5つのモードのエネルギーをプロット

plot(dt * (1:m), E[1,:], label = "mode 1", xlabel = "time", ylabel = "Energy", legend = :topright, legendfontsize=6, ylim=(0.0, 0.4))

plot!(dt * (1:m),E[2,:], label = "mode 2")

plot!(dt * (1:m),E[3,:], label = "mode 3")

plot!(dt * (1:m),E[4,:], label = "mode 4")

plot!(dt * (1:m),E[5,:], label = "mode 5")

ではプロットした結果。

確かにt=10000くらいで最初のモードに戻ってる！

この振動子の動きを動画にしてみた。

@gif for i in 1:10:m

    plot(sol[i].x[1], xlim=(1, n), ylim=(-1.2, 1.2), title = "FPU problem", label="t=$(@sprintf("%.0f", i*dt))",legend=:bottomleft, xlabel="Position", ylabel="Amplitude", size=(800,500))

end

で動画にできる。

結果はこちら。確かにt=10000くらいで最初の振動に戻っている。

PythonでFDTD法で電磁界シミュレーションできるopenEMSを使う(2)例題にあるマイクロストリップパッチアンテナ（MSA）を計算する。Sパラメータや入力インピーダンスだけでなく近傍界から遠方界の変換で指向性も計算できる。電流分布も動画で見る。給電はLumpedポートが使える。

前回はマイクロストリップラインのノッチフィルタを見てみた。

今回はマイクロストリップパッチアンテナ。

チュートリアルはPythonとMatlab/Octaveを両方見ないと分からない。

https://docs.openems.de/python/openEMS/Tutorials/Simple_Patch_Antenna.html

https://wiki.openems.de/index.php/Tutorial:_Simple_Patch_Antenna.html

まずはモジュールをインポート。

import os

import numpy as np

import matplotlib.pyplot as plt

from CSXCAD  import ContinuousStructure

from openEMS import openEMS

from openEMS.physical_constants import C0, EPS0

物理量を設定。Sim_Pathは各自の環境に合わせて変更。後で出ますが今回は寸法は㎜単位。

Sim_Path = os.path.join("C:\\Users\\tomoh\\Documents\\Python Projects\\OpenEMS", 'Simp_Patch')

post_proc_only = False

# patch width (resonant length) in x-direction

patch_width  = 32 #

# patch length in y-direction

patch_length = 40

#substrate setup

substrate_epsR   = 3.38

substrate_kappa  = 1e-3 * 2*np.pi*2.45e9 * EPS0*substrate_epsR

substrate_width  = 60

substrate_length = 60

substrate_thickness = 1.524

substrate_cells = 4

#setup feeding

feed_pos = -6 #feeding position in x-direction

feed_R = 50     #feed resistance

# size of the simulation box

SimBox = np.array([200, 200, 150])

# setup FDTD parameter & excitation function

f0 = 2e9 # center frequency

fc = 1e9 # 20 dB corner frequency

FDTDの設定。ここでパラメータの意味は

• NrTS – シミュレーションするタイムステップの最大数（例：デフォルト=1e9）

• EndCriteria – 終了基準（例：1e-5）。エネルギーがこの値まで減少するとシミュレーションが停止します（<1e-4 が推奨、デフォルト = 1e-5）。

で、境界条件は全方向でMurの吸収境界条件。SetDefaultUnitは1e-3なのでmm単位。メッシュサイズは波長の20分の1。

FDTD = openEMS(NrTS=30000, EndCriteria=1e-4)

FDTD.SetGaussExcite( f0, fc )

FDTD.SetBoundaryCond( ['MUR', 'MUR', 'MUR', 'MUR', 'MUR', 'MUR'] )

CSX = ContinuousStructure()

FDTD.SetCSX(CSX)

mesh = CSX.GetGrid()

mesh.SetDeltaUnit(1e-3)

mesh_res = C0/(f0+fc)/1e-3/20

Jt = CSX.AddDump('Jt', dump_type=3, file_type = 0)

Jt.AddBox([-substrate_width/2, -substrate_length/2, substrate_thickness], [substrate_width/2,  substrate_length/2, substrate_thickness])

パッチアンテナのモデルを作る。金属のパッチとGNDは端をグリッドにするとあるが一切説明ないがどうやらメッシュを見ると1/3ルールを自動でやってくれているようだ。基板の厚み部分は細かくメッシュ。給電はLumpedポート。遠方界のためのBoxも作る。

#initialize the mesh with the "air-box" dimensions

mesh.AddLine('x', [-SimBox[0]/2, SimBox[0]/2])

mesh.AddLine('y', [-SimBox[1]/2, SimBox[1]/2])

mesh.AddLine('z', [-SimBox[2]/3, SimBox[2]*2/3])

# create patch

patch = CSX.AddMetal( 'patch' ) # create a perfect electric conductor (PEC)

start = [-patch_width/2, -patch_length/2, substrate_thickness]

stop  = [ patch_width/2 , patch_length/2, substrate_thickness]

patch.AddBox(priority=10, start=start, stop=stop) # add a box-primitive to the metal property 'patch'

FDTD.AddEdges2Grid(dirs='xy', properties=patch, metal_edge_res=mesh_res/2)

# create substrate

substrate = CSX.AddMaterial( 'substrate', epsilon=substrate_epsR, kappa=substrate_kappa)

start = [-substrate_width/2, -substrate_length/2, 0]

stop  = [ substrate_width/2,  substrate_length/2, substrate_thickness]

substrate.AddBox( priority=0, start=start, stop=stop )

# add extra cells to discretize the substrate thickness

mesh.AddLine('z', np.linspace(0,substrate_thickness,substrate_cells+1))

# create ground (same size as substrate)

gnd = CSX.AddMetal( 'gnd' ) # create a perfect electric conductor (PEC)

start[2]=0

stop[2] =0

gnd.AddBox(start, stop, priority=10)

FDTD.AddEdges2Grid(dirs='xy', properties=gnd)

# apply the excitation & resist as a current source

start = [feed_pos, 0, 0]

stop  = [feed_pos, 0, substrate_thickness]

port = FDTD.AddLumpedPort(1, feed_R, start, stop, 'z', 1.0, priority=5, edges2grid='xy')

mesh.SmoothMeshLines('all', mesh_res, 1.4)

# Add the nf2ff recording box

nf2ff = FDTD.CreateNF2FFBox()

モデルを確認するとこんな感じ。

ここから実際の計算。

if 1:  # debugging only

    CSX_file = os.path.join(Sim_Path, 'simp_patch.xml')

    if not os.path.exists(Sim_Path):

        os.mkdir(Sim_Path)

    CSX.Write2XML(CSX_file)

    from CSXCAD import AppCSXCAD_BIN

    os.system(AppCSXCAD_BIN + ' "{}"'.format(CSX_file))

if not post_proc_only:

    FDTD.Run(Sim_Path, verbose=3, cleanup=False)

電流分布を見てみた。

そしてSパラメータ、入力インピーダンス、指向性の図示。

f = np.linspace(max(1e9,f0-fc),f0+fc,401)

port.CalcPort(Sim_Path, f)

s11 = port.uf_ref/port.uf_inc

s11_dB = 20.0*np.log10(np.abs(s11))

plt.figure()

plt.plot(f/1e9, s11_dB, 'k-', linewidth=2, label='$S_{11}$')

plt.grid()

plt.legend()

plt.ylabel('S-Parameter (dB)')

plt.xlabel('Frequency (GHz)')

idx = np.where((s11_dB<-10) & (s11_dB==np.min(s11_dB)))[0]

if not len(idx)==1:

    print('No resonance frequency found for far-field calulation')

else:

    f_res = f[idx[0]]

    theta = np.arange(-180.0, 180.0, 2.0)

    phi   = [0., 90.]

    nf2ff_res = nf2ff.CalcNF2FF(Sim_Path, f_res, theta, phi, center=[0,0,1e-3])

    plt.figure()

    E_norm = 20.0*np.log10(nf2ff_res.E_norm[0]/np.max(nf2ff_res.E_norm[0])) + nf2ff_res.Dmax[0]

    plt.plot(theta, np.squeeze(E_norm[:,0]), 'k-', linewidth=2, label='xz-plane')

    plt.plot(theta, np.squeeze(E_norm[:,1]), 'r--', linewidth=2, label='yz-plane')

    plt.grid()

    plt.ylabel('Directivity (dBi)')

    plt.xlabel('Theta (deg)')

    plt.title('Frequency: {} GHz'.format(f_res/1e9))

    plt.legend()

Zin = port.uf_tot/port.if_tot

plt.figure()

plt.plot(f/1e9, np.real(Zin), 'k-', linewidth=2, label='$\Re\{Z_{in}\}$')

plt.plot(f/1e9, np.imag(Zin), 'r--', linewidth=2, label='$\Im\{Z_{in}\}$')

plt.grid()

plt.legend()

plt.ylabel('Zin (Ohm)')

plt.xlabel('Frequency (GHz)')

plt.show()

だいぶ雰囲気は分かってきた。次は導波管系をやってみよう。

追記：

Et = CSX.AddDump('Et', dump_type=0, file_type = 0)

Et.AddBox([-SimBox[0]/2,-SimBox[1]/2,-SimBox[2]/3],[SimBox[0]/2,SimBox[1]/2,SimBox[2]*2/3])

を追加すると電界分布も見られる。

こちら。

PythonでFDTD法で電磁界シミュレーションできるopenEMSを使う(1)例題にあるマイクロストリップラインのノッチフィルタ（スタブ）を動かして電磁界分布を動画で見てみる。CSXCADでモデルは確認できるし、ParaViewで電磁界分布が見られる。Sパラメータも計算できる。

オープンソースの電磁界シミュレータはいろいろある。

Open-Source Electromagnetic Simulation: FDTD, FEM, MoM

今回はその中でEC-FDTE法のopenEMSをPythonから使ってみる。

まず公式サイトはこちら。

https://www.openems.de/

WindowsでPythonを使うときのインストールはこちら。

https://docs.openems.de/python/install.html#windows

Pythonのバージョンは3.10か3.11でないと動かないっぽいので3.11を仮想環境をvenvで作って使う。

以下をpipでインストールして、

openEMS-0.0.36-cp311-cp311-win_amd64.whl

CSXCAD-0.6.3-cp311-cp311-win_amd64.whl

環境変数を設定。

setx OPENEMS_INSTALL_PATH インストールしたフォルダ

これで動く。

問題はPythonで使うときのドキュメントがものすごく不親切（もともとMatlab/Octaveから使うもの）で、両方見ないと分からん…
とりあえず例題のマイクロストリップラインのノッチフィルタ（スタブ）やってみよう。
PythonとMatlabを見比べないと分からない部分が多くある。

https://docs.openems.de/python/openEMS/Tutorials/MSL_NotchFilter.html

https://wiki.openems.de/index.php/Tutorial:_Microstrip_Notch_Filter.html

まずは必要なモジュールをインポート。

import os

import numpy as np

import matplotlib.pyplot as plt

from CSXCAD  import ContinuousStructure

from openEMS import openEMS

from openEMS.physical_constants import C0

設定。ここでPathは各自保存したいところに変更。 

Sim_Path = os.path.join("C:\\Users\\tomoh\\Documents\\Python Projects\\OpenEMS", 'NotchFilter')

post_proc_only = False

unit = 1e-6 # specify everything in um

MSL_length = 50000

MSL_width = 600

substrate_thickness = 254

substrate_epr = 3.66

stub_length = 12e3

f_max = 7e9

FDTDの設定。ここで励振はガウシアンで、

• f0: the center frequency of the pulse
• fc : 20dB cutoff frequency --> bandwidth is 2*fc

また境界条件は [xmin xmax ymin ymax zmin zmax]の順で、PML-8は完全吸収層で8セル使うもの、MURはムーアの吸収境界条件、
PECは完全導体でGNDを表している。x方向がPML-8なのはポートがつく方向で、終端したいので。

FDTD = openEMS()

FDTD.SetGaussExcite( f_max/2, f_max/2 )

FDTD.SetBoundaryCond( ['PML_8', 'PML_8', 'MUR', 'MUR', 'PEC', 'MUR'] )

メッシュの基本設定。電流分布はAddDumpで設定する。

CSX = ContinuousStructure()

FDTD.SetCSX(CSX)

mesh = CSX.GetGrid()

mesh.SetDeltaUnit(unit)

resolution = C0/(f_max*np.sqrt(substrate_epr))/unit/50 # resolution of lambda/50

third_mesh = np.array([2*resolution/3, -resolution/3])/4

Jt = CSX.AddDump('Jt', dump_type=3, file_type = 0)

Jt.AddBox([-MSL_length, -MSL_length, substrate_thickness], [MSL_length, MSL_length, substrate_thickness])

でここが一番面倒なところで、メッシュを細かくする部分を手動で記述する。HFSSのようなアダプティブメッシュに慣れきっていると
これがわからないところかも。基本は1/3ルール（2D 金属のエッジでは電界が強くなるため、FDTD では正しく計算するのが困難になります。

つまり、例えばマイクロストリップライン（MSL）の場合、線路インピーダンスと波動伝播はシミュレーションと実測で多少異なります。そのため、精密計算では、金属の1/3を内側、2/3を外側に配置するメッシュラインを配置することが推奨されます。）を使う。メッシュは円筒形もあり。

mesh.AddLine('x', 0)

mesh.AddLine('x',  MSL_width/2+third_mesh)

mesh.AddLine('x', -MSL_width/2-third_mesh)

mesh.SmoothMeshLines('x', resolution/4)

mesh.AddLine('x', [-MSL_length, MSL_length])

mesh.SmoothMeshLines('x', resolution)

mesh.AddLine('y', 0)

mesh.AddLine('y',  MSL_width/2+third_mesh)

mesh.AddLine('y', -MSL_width/2-third_mesh)

mesh.SmoothMeshLines('y', resolution/4)

mesh.AddLine('y', [-15*MSL_width, 15*MSL_width+stub_length])

mesh.AddLine('y', (MSL_width/2+stub_length)+third_mesh)

mesh.SmoothMeshLines('y', resolution)

mesh.AddLine('z', np.linspace(0,substrate_thickness,5))

mesh.AddLine('z', 3000)

mesh.SmoothMeshLines('z', resolution)

基板の設定。Rogersを使うこと想定。

substrate = CSX.AddMaterial( 'RO4350B', epsilon=substrate_epr)

start = [-MSL_length, -15*MSL_width, 0]

stop  = [+MSL_length, +15*MSL_width+stub_length, substrate_thickness]

substrate.AddBox(start, stop )

マイクロストリップラインポートの設定。ポートをx=0で交差させると勝手にスルーラインができるそうだ。

port = [None, None]

pec = CSX.AddMetal( 'PEC' )

portstart = [ -MSL_length, -MSL_width/2, substrate_thickness]

portstop  = [ 0,  MSL_width/2, 0]

port[0] = FDTD.AddMSLPort( 1,  pec, portstart, portstop, 'x', 'z', excite=-1, FeedShift=10*resolution, MeasPlaneShift=MSL_length/3, priority=10)

portstart = [MSL_length, -MSL_width/2, substrate_thickness]

portstop  = [0         ,  MSL_width/2, 0]

port[1] = FDTD.AddMSLPort( 2, pec, portstart, portstop, 'x', 'z', MeasPlaneShift=MSL_length/3, priority=10 )

 スタブを追加。

start = [-MSL_width/2,  MSL_width/2, substrate_thickness]

stop  = [ MSL_width/2,  MSL_width/2+stub_length, substrate_thickness]

pec.AddBox(start, stop, priority=10 )

計算。確認のためにCSXCADでメッシュも見てみる。

if 1:  # debugging only

    CSX_file = os.path.join(Sim_Path, 'notch.xml')

    if not os.path.exists(Sim_Path):

        os.mkdir(Sim_Path)

    CSX.Write2XML(CSX_file)

    from CSXCAD import AppCSXCAD_BIN

    os.system(AppCSXCAD_BIN + ' "{}"'.format(CSX_file))

if not post_proc_only:

    FDTD.Run(Sim_Path, cleanup=False)

でSパラメータを計算して図示。

f = np.linspace( 1e6, f_max, 1601 )

for p in port:

    p.CalcPort( Sim_Path, f, ref_impedance = 50)

s11 = port[0].uf_ref / port[0].uf_inc

s21 = port[1].uf_ref / port[0].uf_inc

plt.plot(f/1e9,20*np.log10(abs(s11)),'k-',linewidth=2 , label='$S_{11}$')

plt.grid()

plt.plot(f/1e9,20*np.log10(abs(s21)),'r--',linewidth=2 , label='$S_{21}$')

plt.legend()

plt.ylabel('S-Parameter (dB)')

plt.xlabel('frequency (GHz)')

plt.ylim([-40, 2])

plt.show()

ParaViewで電磁界分布見てみる。ファイル開いてApply、リスケールするのと、表示をSolid ColorからRot-H Fieldに変更しないと見られない。
実はここにハマって時間消費した…

なるほどもっともらしい。メッシュを切る部分さえ慣れれば何でも計算できそうだ。

次はパッチアンテナやってみよう（続く）。