[Python×Ansys Lumerical FDTD]任意の３次元曲線の光導波路をLumerical FDTDに実装する方法

はじめに

この記事を書く動機として、光電融合デバイスにおける長尺の三次元光配線や機能素子をAnsys Cloud Burst ComputeでFDTDシミュレーションしてみたいと思ったからです。

Ansys Cloud Burst Computeとは?

Ansys Cloud Burst ComputeとはAnsysの主要ソルバーと統合されたオンデマンド型の高性能コンピューティング（HPC）サービスで、エンジニアが自前のハードウェア制限を超えてクラウド上でシミュレーションを実行できるサービスです。 セットアップやクラウドインフラ管理の専門知識が不要で、必要なときに柔軟にスケールアップ・スケールダウンして計算資源を利用できます。
https://www.ansys.com/ja-jp/products/cloud/ansys-cloud-burst-compute

近年、光電融合の実装が進むにつれて、ポリマー光配線の立体的な引き回しや、偏波ローテータ・偏波制御素子のように三次元的な形状を本質的に含む構造を扱う場面が増えています。一方で、これらの構造をそのまま FDTD で解析しようとすると、モデルの自由度や計算領域が大きくなり、計算コストが急激に増大してしまいます。

本記事では、こうした三次元の大規模 FDTD 問題を扱うために、Python を用いて長尺の三次元形状を自動生成し、STL として FDTD に取り込むワークフローを紹介します。

特に、曲線導波路や立体配線を柔軟に記述できる点は、ポリマー光配線の引き回し設計や、三次元偏波ローテータの形状検討において有効です。また、Lumerical FDTD の Ansys Cloud Burst Compute 機能を前提とすることで、ローカル環境では扱いにくい大規模モデルでも計算を回せることを想定しています。

スクリプト全体

今回はテストとしてx方向に進行しながらy 方向にのみサイン的に振動し
z 方向は単調に（あるいは一定値のまま）変化するだけで、周期的な振動は持たないパターンを作ってみます。
下記のスクリプトは、Lumerical FDTD 用に三次元の曲線導波路形状を Python から自動生成するものです。中心線は x 方向に進むサインカーブとして定義されますが、端部で不連続な折れが生じないよう、サイン波の振幅に包絡関数を掛けています。具体的には
$y(x)=A*w(x)\sin\left(\dfrac{2\pi n x}{L}\right)$
とし、包絡関数 $w(x)$ は
$w(0)=w(L)=0$ かつ $w'(0)=w'(L)=0$
を満たすように設計されています。この条件により、端点で
$\dfrac{dy}{dx}=0$
となり、導波路の入口・出口で中心線の接線が $+x$ 方向に揃います。

z 方向についても同様に、必要に応じて
$z(x)=z_{\mathrm{rise}},s(x/L)$
（$s'(0)=s'(1)=0$ を満たす smoothstep 関数）を用いることで、端点で
$\dfrac{dz}{dx}=0$
を実現できます。これにより、入口および出口の断面は幾何学的に x 方向に垂直になります。

得られた中心線に対しては、parallel transport frame を用いて断面のねじれを抑えながら楕円断面（$a=b$ で円）を掃引し、stlポリゴンメッシュとして立体形状を構築します。
最終的に、この形状を STL ファイルとして書き出し、stlimport() を用いて FDTD に直接読み込むことで、入口・出口の面向きと滑らかな曲線形状を両立した三次元導波路モデルを自動生成しています。

# 3d_sine_wg_ellipse_eased_full.py
# ------------------------------------------------------------
# 端部をイーズイン/イーズアウトして接線不連続（折れ）を解消しつつ、
# 入口/出口で dy/dx=0（さらにdz/dx=0も可能）を満たして
# 入口/出口断面が x に垂直になるようにする。
#
# 断面: 楕円（a,b） a=bで円
# STL -> FDTD stlimport()
# ------------------------------------------------------------

import os
import struct
import numpy as np
import lumapi


# ============================================================
# PARAMETERS
# ============================================================
L = 200e-6
A = 10e-6
N_PERIODS = 2.0
N = 900

# 端部の“滑らか化”長さ（全長に対する割合 or 絶対長）
EASE_LEN = 20e-6   # [m] 端からこの長さで包絡を0→1に立ち上げ/立ち下げ

# z の扱い
Z_MODE = "ease"    # "ease" or "linear"
Z_SLOPE = 0.2      # linearの場合: z = Z_SLOPE*x
Z_RISE = 40e-6     # easeの場合: zは 0→Z_RISE を端で滑らかに持ち上げ（例）

# 断面（楕円）
ELLIPSE_A = 0.30e-6
ELLIPSE_B = 0.15e-6   # a=bで円
ELLIPSE_M = 64

# STL / FDTD
STL_FILENAME = "wg_3d_sine_ellipse_eased_m_units.stl"
FDTD_OBJ_NAME = "wg_3d_sine_ellipse_eased"
FDTD_MATERIAL = "Si (Silicon) - Palik"
STLIMPORT_SCALING = 1.0


# ============================================================
# Utility
# ============================================================
def unit(v: np.ndarray, eps: float = 1e-15) -> np.ndarray:
    n = float(np.linalg.norm(v))
    if n < eps:
        return np.zeros_like(v)
    return v / n


def smoothstep01(u: np.ndarray) -> np.ndarray:
    """
    0..1 -> 0..1 の C1 スムーズステップ
      s(u)=3u^2-2u^3
    端で s'(0)=s'(1)=0
    """
    u = np.clip(u, 0.0, 1.0)
    return u * u * (3.0 - 2.0 * u)


def ease_window(x: np.ndarray, L: float, ease_len: float) -> np.ndarray:
    """
    両端で0、内部で1になる包絡 w(x)
    - 端から ease_len で 0->1 (smoothstep)
    - 端から ease_len で 1->0 (smoothstep)
    w'(端)=0 なので dy/dx の不連続が出にくい
    """
    if ease_len <= 0:
        return np.ones_like(x)

    # 左端
    uL = x / ease_len
    wL = smoothstep01(uL)

    # 右端
    uR = (L - x) / ease_len
    wR = smoothstep01(uR)

    return np.minimum(wL, wR)


def make_centerline_eased(L, A, n_periods, N, ease_len, z_mode, z_slope, z_rise) -> np.ndarray:
    """
    端で dy/dx=0 を実現する中心線。
    y = A * w(x) * sin(2π n x/L)
    w(0)=w(L)=0 かつ w'(0)=w'(L)=0

    z:
      - linear: z=z_slope*x（端でdz/dx≠0）
      - ease:   z=z_rise * s(x/L)（端でdz/dx=0）
    """
    x = np.linspace(0.0, L, int(N))
    w = ease_window(x, L, ease_len)

    y = A * w * np.sin(2.0 * np.pi * n_periods * x / L)

    if z_mode == "linear":
        z = z_slope * x
    elif z_mode == "ease":
        # 0..L を 0..1 に正規化して smoothstep で持ち上げ（端でdz/dx=0）
        s = smoothstep01(x / L)
        z = z_rise * s
    else:
        raise ValueError('Z_MODE must be "ease" or "linear".')

    return np.column_stack([x, y, z]).astype(np.float64)


# ============================================================
# parallel transport frame
# ============================================================
def parallel_transport_frames(P: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    Np = P.shape[0]
    if Np < 2:
        raise ValueError("中心線点列 P は2点以上必要です。")

    dP = np.zeros_like(P)
    dP[1:-1] = P[2:] - P[:-2]
    dP[0] = P[1] - P[0]
    dP[-1] = P[-1] - P[-2]
    t = np.array([unit(v) for v in dP], dtype=np.float64)

    ref = np.array([0.0, 0.0, 1.0], dtype=np.float64)
    if abs(float(np.dot(t[0], ref))) > 0.9:
        ref = np.array([0.0, 1.0, 0.0], dtype=np.float64)

    n0 = unit(np.cross(ref, t[0]))
    if np.linalg.norm(n0) < 1e-12:
        ref = np.array([1.0, 0.0, 0.0], dtype=np.float64)
        n0 = unit(np.cross(ref, t[0]))
    b0 = unit(np.cross(t[0], n0))

    n = np.zeros_like(P, dtype=np.float64)
    b = np.zeros_like(P, dtype=np.float64)
    n[0], b[0] = n0, b0

    for i in range(1, Np):
        t_prev = t[i - 1]
        t_curr = t[i]

        axis = np.cross(t_prev, t_curr)
        axis_norm = float(np.linalg.norm(axis))

        if axis_norm < 1e-12:
            ni = n[i - 1].copy()
        else:
            axis_u = axis / axis_norm
            c = float(np.clip(np.dot(t_prev, t_curr), -1.0, 1.0))
            ang = float(np.arccos(c))
            v = n[i - 1]
            ni = (v * np.cos(ang)
                  + np.cross(axis_u, v) * np.sin(ang)
                  + axis_u * np.dot(axis_u, v) * (1.0 - np.cos(ang)))

        ni = unit(ni - t_curr * np.dot(t_curr, ni))
        bi = unit(np.cross(t_curr, ni))
        n[i] = ni
        b[i] = bi

    return t, n, b


# ============================================================
# swept ellipse mesh
# ============================================================
def swept_ellipse_mesh(P: np.ndarray, a: float, b: float, M: int = 64) -> tuple[np.ndarray, np.ndarray]:
    if a <= 0 or b <= 0:
        raise ValueError("a,b は正である必要があります。")
    if M < 8:
        raise ValueError("M は 8 以上を推奨します。")

    _, n, bvec = parallel_transport_frames(P)
    Np = P.shape[0]

    th = np.linspace(0.0, 2.0 * np.pi, int(M), endpoint=False)
    ct = np.cos(th)
    st = np.sin(th)

    ring = np.zeros((Np, M, 3), dtype=np.float64)
    for i in range(Np):
        c = P[i]
        ni = n[i]
        bi = bvec[i]
        ring[i, :, :] = c + (a * ct)[:, None] * ni[None, :] + (b * st)[:, None] * bi[None, :]

    vertices = ring.reshape(-1, 3).astype(np.float32)

    def vid(i, j):
        return i * M + j

    faces = []

    for i in range(Np - 1):
        for j in range(M):
            j2 = (j + 1) % M
            a0 = vid(i, j)
            a1 = vid(i, j2)
            b1 = vid(i + 1, j2)
            b0 = vid(i + 1, j)
            faces.append([a0, a1, b1])
            faces.append([a0, b1, b0])

    # caps
    start_center_id = vertices.shape[0]
    end_center_id = vertices.shape[0] + 1
    vertices = np.vstack([vertices, P[0].astype(np.float32)[None, :], P[-1].astype(np.float32)[None, :]])

    for j in range(M):
        j2 = (j + 1) % M
        faces.append([start_center_id, vid(0, j2), vid(0, j)])

    for j in range(M):
        j2 = (j + 1) % M
        faces.append([end_center_id, vid(Np - 1, j), vid(Np - 1, j2)])

    return vertices, np.array(faces, dtype=np.int32)


# ============================================================
# STL writer
# ============================================================
def write_binary_stl(path: str, vertices: np.ndarray, faces: np.ndarray, header: bytes = b"wg_stl"):
    header = header[:80].ljust(80, b"\0")
    tri_count = int(faces.shape[0])

    with open(path, "wb") as f:
        f.write(header)
        f.write(struct.pack("<I", tri_count))

        for tri in faces:
            v0 = vertices[int(tri[0])]
            v1 = vertices[int(tri[1])]
            v2 = vertices[int(tri[2])]

            nrm = np.cross(v1 - v0, v2 - v0)
            nrm = unit(nrm).astype(np.float32)

            f.write(struct.pack("<3f", float(nrm[0]), float(nrm[1]), float(nrm[2])))
            f.write(struct.pack("<3f", float(v0[0]), float(v0[1]), float(v0[2])))
            f.write(struct.pack("<3f", float(v1[0]), float(v1[1]), float(v1[2])))
            f.write(struct.pack("<3f", float(v2[0]), float(v2[1]), float(v2[2])))
            f.write(struct.pack("<H", 0))


# ============================================================
# FDTD import via stlimport()
# ============================================================
def import_stl_into_fdtd_via_stlimport(fdtd, stl_path: str, obj_name: str, material: str, scaling_factor: float = 1.0):
    stl_path = os.path.abspath(stl_path).replace("\\", "/")
    cmd = f"""
    stlimport("{stl_path}", {scaling_factor});
    set("name","{obj_name}");
    set("material","{material}");
    """
    fdtd.eval(cmd)


# ============================================================
# MAIN
# ============================================================
if __name__ == "__main__":
    P = make_centerline_eased(
        L=L, A=A, n_periods=N_PERIODS, N=N,
        ease_len=EASE_LEN,
        z_mode=Z_MODE, z_slope=Z_SLOPE, z_rise=Z_RISE
    )

    vertices, faces = swept_ellipse_mesh(P, a=ELLIPSE_A, b=ELLIPSE_B, M=ELLIPSE_M)

    stl_path = os.path.abspath(STL_FILENAME)
    write_binary_stl(stl_path, vertices, faces, header=b"3d_sine_ellipse_wg_eased_m_units")

    fdtd = lumapi.FDTD()
    import_stl_into_fdtd_via_stlimport(
        fdtd,
        stl_path=stl_path,
        obj_name=FDTD_OBJ_NAME,
        material=FDTD_MATERIAL,
        scaling_factor=STLIMPORT_SCALING
    )

出力結果

こちらが出力結果になります。

sineカーブがLumerical FDTDにインポートされていることが確認できます。
以下が導波路の左端

続いて、右端

作製した導波路の端面がいずれもx軸に対して垂直になっていることが確認できます。

こうすることで端面の端部分の導波路モードを設定すれば、FDTDシミュレーションできるようになるでしょう。

※本記事は筆者個人の見解であり、所属組織の公式見解を示すものではありません。

問い合わせ

光学シミュレーションソフトの導入や技術相談、
設計解析委託をお考えの方はサイバネットシステムにお問合せください。

光学ソリューションサイトについては以下の公式サイトを参照：
👉 光学ソリューションサイト(サイバネット)

光学分野のエンジニアリングサービスについては以下の公式サイトを参照：
👉 光学エンジニアリングサービス(サイバネット)