パソコン・インターネット

2025年6月24日 (火)

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

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

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

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

 

2025年6月16日 (月)

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

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

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

 

2025年6月13日 (金)

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)

Openemsbendpatch1

そして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()

Openemsbendpatch3

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])

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

 

2025年6月 9日 (月)

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

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

∂φ/∂t = φ - φ3 - (∇2+ k02)2φ

をやってみよう。空間の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()とすると動画ができる。

 

 

2025年6月 6日 (金)

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

この記事見た。

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

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

動画で動作を紹介。

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

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

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

Windowsedit

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

2025年6月 2日 (月)

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()

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

結果はこちら。

 

らせんを描いている!

 

 

2025年5月28日 (水)

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)

 

Openemswaveguide1

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$')

Openemswaveguide2

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

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

 

 

 

2025年5月22日 (木)

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

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

Juliafpu3

これを確かめてみよう。

まず、必要なパッケージをインストールする。が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")

ではプロットした結果。

Fputjulia1

確かに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くらいで最初の振動に戻っている。

 

 

 

 

 

2025年5月19日 (月)

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()

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

Openemspatch2

ここから実際の計算。


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()

Openemspatch3

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

追記:

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])
を追加すると電界分布も見られる。
こちら。

2025年5月14日 (水)

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)
Openemsnotch5

で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()

 

Openemsnotch2

 

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

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

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

より以前の記事一覧

最近の記事

2025年6月
1 2 3 4 5 6 7
8 9 10 11 12 13 14
15 16 17 18 19 20 21
22 23 24 25 26 27 28
29 30          

最近のコメント

無料ブログはココログ
フォト