openEMS Blows Up when Attempting to Simulate a 45-degree Turn in a Microstrip

[biergaizi]

I'm trying to simulate a 45-degree turn in a microstrip. In the code, I basically created a horizontal microstrip and then rotated it by 45 degrees with a CSTransform. However, my code and simulation setup cause openEMS to blow up with infinite energy without producing any meaningful output.

Without the rotation, the simulation runs without blowing up. I suspect the problem is caused by propagation_axis=Axis("x") in the Microstrip class. It can only be either horizontal or vertical. I suspect that rotating it by 45 degrees breaks some assumptions made in pyems, thus the code is incorrect, but I'm not knowledgeable enough to be sure.

What is the correct way to perform this simulation?

3D View

Code

• project.py

This file contains some board definitions and two utility functions for creating a rotated microstrip.

from math import sin, cos, radians
from CSXCAD.CSTransform import CSTransform
from pyems.pcb import PCBProperties
from pyems.material import Dielectric
from pyems.structure import SMDPassiveDimensions
from pyems.coordinate import Box2, Coordinate2

project_dielectrics = {
    "Generic-FR-4": Dielectric(epsr=[(1, 4.6)], loss=[(1, 0.03)]),
}

project_pcbs = {
    "layer1-to-layer3": PCBProperties(
        substrate=project_dielectrics["Generic-FR-4"],
        copper_thickness=[0.0356e-3, 0.0178e-3],
        substrate_thickness=[0.4825e-3],
        metal_conductivity=5.8e7,
        via_plating_thickness=0.020e-3,
    )
}


def BeginLengthWidth(
    begin: Coordinate2,
    length: float,
    width: float
):
    return (
        Box2(
            Coordinate2(begin.x, begin.y - width / 2),
            Coordinate2(begin.x + length, begin.y + width / 2),
        ),
        Coordinate2(begin.x + length, begin.y)
    )


def BeginLengthWidthRotation(
    begin: Coordinate2,
    length: float,
    width: float,
    degrees: int
):
    box = Box2(
            Coordinate2(begin.x, begin.y - width / 2),
            Coordinate2(begin.x + length, begin.y + width / 2),
    )
    
    tr = CSTransform()
    tr.AddTransform('RotateAxis', 'z', degrees)

    rotated_begin_x = (
                   (begin.x - box.center().x) * cos(radians(degrees)) -
                   (begin.y - box.center().y) * sin(radians(degrees))
               ) + box.center().x
    rotated_begin_y = (
                   (begin.x - box.center().x) * sin(radians(degrees)) +
                   (begin.y - box.center().y) * cos(radians(degrees))
               ) + box.center().y

    tr.AddTransform('Translate', [
        begin.x - rotated_begin_x,
        begin.y - rotated_begin_y,
        0
    ])

    end_x = begin.x + length
    end_y = begin.y
    rotated_end_x = (
                   (end_x - begin.x) * cos(radians(degrees)) -
                   (end_y - begin.y) * sin(radians(degrees))
               ) + begin.x
    rotated_end_y = (
                   (end_x - begin.x) * sin(radians(degrees)) +
                   (end_y - begin.y) * cos(radians(degrees))
               ) + begin.y

    return (
        box,
        Coordinate2(rotated_end_x, rotated_end_y),
        tr
    )

• main.py

This is the main simulation code.

import numpy as np
from CSXCAD.CSTransform import CSTransform
from pyems.pcb import common_pcbs
from pyems.simulation import Simulation
from pyems.field_dump import FieldDump
from pyems.utilities import print_table, mil_to_mm
from pyems.structure import (
    PCB,
    Microstrip,
    common_smd_passives,
    SMDPassive,
)
from pyems.coordinate import Box2, Box3, Coordinate2, Coordinate3, Axis
from pyems.mesh import Mesh
from project import project_pcbs, BeginLengthWidth, BeginLengthWidthRotation

unit = 1e-3
freq = np.arange(0, 18e9, 1e7)
sim = Simulation(freq=freq, unit=unit)
pcb_prop = project_pcbs["layer1-to-layer3"]
pcb_len = 9
pcb_width = 5
trace_width = 0.85
z0_ref = 50

pcb = PCB(
    sim=sim,
    pcb_prop=pcb_prop,
    length=pcb_len,
    width=pcb_width,
    layers=range(3),
    omit_copper=[0],
)

box, source_end = BeginLengthWidth(
    Coordinate2(-pcb_len / 2, 0), 6, trace_width
)
    
microstrip_source = Microstrip(
    pcb=pcb,
    position=box.center(),
    length=box.length(),
    width=box.width(),
    propagation_axis=Axis("x"),
    trace_layer=0,
    gnd_layer=1,
    port_number=1,
    excite=True,
)

box, straight_line_end, tr = BeginLengthWidthRotation(
    # changing the angle from 45 to 0 allows the simulation to run
    source_end, 1, trace_width + 0.06, 45
)
microstrip_corner = Microstrip(
    pcb=pcb,
    position=(box.center().x - 0.33, box.center().y - 0.11),
    length=box.length(),
    width=box.width(),
    transform=tr,
    propagation_axis=Axis("x"),
    trace_layer=0,
    gnd_layer=1,
)
straight_line_end.x -= 0.65
straight_line_end.y -= 0.21

box, _ = BeginLengthWidth(
    straight_line_end, pcb_len / 2 - straight_line_end.x, trace_width
)
    
microstrip_p2 = Microstrip(
    pcb=pcb,
    position=(box.center().x, box.center().y),
    length=box.length(),
    width=box.width(),
    propagation_axis=Axis("x"),
    trace_layer=0,
    gnd_layer=1,
    port_number=2,
)

dump = FieldDump(
    sim=sim,
    box=Box3(
        Coordinate3(-pcb_len / 2, -pcb_width / 2, 0),
        Coordinate3(pcb_len / 2, pcb_width / 2, 0),
    ),
)

mesh = Mesh(
    sim=sim,
    metal_res=1 / 120,
    nonmetal_res=1 / 40,
    smooth=(1.5, 1.5, 1.5),
    min_lines=5,
    expand_bounds=((0, 0), (8, 8), (8, 8)),
)

sim.run()
sim.view_field()

print_table(
    data=[sim.freq / 1e9, sim.s_param(1, 1), sim.s_param(2, 1)],
    col_names=["freq", "s11", "s21", "s31"],
    prec=[4, 4, 4],
)

[matthuszagh]

The fundamental problem here is that your structure is producing a very fine mesh in places. The x and y coordinates of the 45deg microstrip after the transformation don't line up perfectly with the corners of the microstrip trace corners. The mesher therefore sees miniscule features there that it tries to account for. It might be possible to adjust the transform to get the results to line up nicely with the microstrips, but I haven't tried.

A better approach would be to construct a polygon for the trace. There's a public-facing function for this, but you currently have to do a little manual work to get it to integrate with the PCB structure. It might be worth me adding a class that handles this. But, the current form really isn't so bad. In any event, take a look at the following code. I've butchered your example a bit, but it should be easy for you to fix that. The important point is how I created the polygon trace.

import numpy as np
from CSXCAD.CSTransform import CSTransform
from pyems.simulation import Simulation
from pyems.field_dump import FieldDump
from pyems.utilities import print_table
from pyems.structure import (
    PCB,
    Microstrip,
)
from pyems.coordinate import Box2, Box3, Coordinate2, Coordinate3, Axis
from pyems.mesh import Mesh
from math import sin, cos, radians
from pyems.pcb import PCBProperties
from pyems.material import Dielectric
from pyems.coordinate import Box2, Coordinate2
from pyems.csxcad import construct_polygon, add_conducting_sheet
from pyems.priority import priorities

project_dielectrics = {
    "Generic-FR-4": Dielectric(epsr=[(1, 4.6)], loss=[(1, 0.03)]),
}

project_pcbs = {
    "layer1-to-layer3": PCBProperties(
        substrate=project_dielectrics["Generic-FR-4"],
        copper_thickness=[0.0356e-3, 0.0178e-3],
        substrate_thickness=[0.4825e-3],
        metal_conductivity=5.8e7,
        via_plating_thickness=0.020e-3,
    )
}


def BeginLengthWidth(begin: Coordinate2, length: float, width: float):
    return (
        Box2(
            Coordinate2(begin.x, begin.y - width / 2),
            Coordinate2(begin.x + length, begin.y + width / 2),
        ),
        Coordinate2(begin.x + length, begin.y),
    )


def BeginLengthWidthRotation(
    begin: Coordinate2, length: float, width: float, degrees: int
):
    box = Box2(
        Coordinate2(begin.x, begin.y - width / 2),
        Coordinate2(begin.x + length, begin.y + width / 2),
    )

    tr = CSTransform()
    tr.AddTransform("RotateAxis", "z", degrees)

    rotated_begin_x = (
        (begin.x - box.center().x) * cos(radians(degrees))
        - (begin.y - box.center().y) * sin(radians(degrees))
    ) + box.center().x
    rotated_begin_y = (
        (begin.x - box.center().x) * sin(radians(degrees))
        + (begin.y - box.center().y) * cos(radians(degrees))
    ) + box.center().y

    tr.AddTransform(
        "Translate", [begin.x - rotated_begin_x, begin.y - rotated_begin_y, 0]
    )

    end_x = begin.x + length
    end_y = begin.y
    rotated_end_x = (
        (end_x - begin.x) * cos(radians(degrees))
        - (end_y - begin.y) * sin(radians(degrees))
    ) + begin.x
    rotated_end_y = (
        (end_x - begin.x) * sin(radians(degrees))
        + (end_y - begin.y) * cos(radians(degrees))
    ) + begin.y

    return (box, Coordinate2(rotated_end_x, rotated_end_y), tr)


unit = 1e-3
freq = np.arange(0, 18e9, 1e7)
sim = Simulation(freq=freq, unit=unit)
pcb_prop = project_pcbs["layer1-to-layer3"]
pcb_len = 9
pcb_width = 5
trace_width = 0.85
z0_ref = 50

pcb = PCB(
    sim=sim,
    pcb_prop=pcb_prop,
    length=pcb_len,
    width=pcb_width,
    layers=range(3),
    omit_copper=[0],
)

box, source_end = BeginLengthWidth(
    Coordinate2(-pcb_len / 2, 0), 6, trace_width
)

microstrip_source = Microstrip(
    pcb=pcb,
    position=box.center(),
    length=box.length(),
    width=box.width(),
    propagation_axis=Axis("x"),
    trace_layer=0,
    gnd_layer=1,
    port_number=1,
    excite=True,
)

_, straight_line_end, tr = BeginLengthWidthRotation(
    # changing the angle from 45 to 0 allows the simulation to run
    source_end,
    1,
    trace_width + 0.06,
    45,
)
straight_line_end.x -= 0.65
straight_line_end.y -= 0.21

new_box, _ = BeginLengthWidth(
    straight_line_end, pcb_len / 2 - straight_line_end.x, trace_width
)

microstrip_p2 = Microstrip(
    pcb=pcb,
    position=(new_box.center().x, new_box.center().y),
    length=new_box.length(),
    width=new_box.width(),
    propagation_axis=Axis("x"),
    trace_layer=0,
    gnd_layer=1,
    port_number=2,
)

trace_prop = add_conducting_sheet(
    csx=pcb.sim.csx,
    name="trace_poly",
    conductivity=pcb.pcb_prop.metal_conductivity(),
    thickness=pcb.pcb_prop.copper_thickness(0),
)

construct_polygon(
    prop=trace_prop,
    points=[
        Coordinate2(
            box.center().x + box.length() / 2,
            box.center().y - box.width() / 2,
        ),
        Coordinate2(
            box.center().x + box.length() / 2,
            box.center().y + box.width() / 2,
        ),
        Coordinate2(
            new_box.center().x - new_box.length() / 2,
            new_box.center().y + new_box.width() / 2,
        ),
        Coordinate2(
            new_box.center().x - new_box.length() / 2,
            new_box.center().y - new_box.width() / 2,
        ),
    ],
    normal=Axis("z"),
    elevation=pcb.copper_layer_elevation(0),
    priority=priorities["trace"],
)

dump = FieldDump(
    sim=sim,
    box=Box3(
        Coordinate3(-pcb_len / 2, -pcb_width / 2, 0),
        Coordinate3(pcb_len / 2, pcb_width / 2, 0),
    ),
)

mesh = Mesh(
    sim=sim,
    metal_res=1 / 120,
    nonmetal_res=1 / 40,
    smooth=(1.5, 1.5, 1.5),
    min_lines=5,
    expand_bounds=((0, 0), (8, 8), (8, 8)),
)

sim.run()
sim.view_field()

print_table(
    data=[sim.freq / 1e9, sim.s_param(1, 1), sim.s_param(2, 1)],
    col_names=["freq", "s11", "s21", "s31"],
    prec=[4, 4, 4],
)

Let me know if that's unclear.

A few other thoughts: the spacing between the x-axis microstrips is a bit small, which will increase your simulation time. There's no problem with that and I expect the penalty won't be that bad, but it's worth being aware of. If you intend to fabricate this exact structure and want to know the behavior then of course you want to simulate this as is. But, if all you really care about is the effect of the 45 degree bend in the trace it might be worth widening that gap a little to improve the simulation time.

It's also worth noting that while the openems blow-up in this case was probably caused by the really fine mesh, there are cases where openems will blow up from reasonable models. I was simulating a simple coax line the other day over a wide bandwidth to visualize non-TEM wave propagation and openems also saw infinite energy problems. The solution was to simulate narrower bandwidths in parts. The same is true for the horn antenna in the examples. This isn't something caused by pyems and I don't have a firm grasp on why it happens yet. Anyway, just something to be aware of.

Let me know if you run into any issues.