How to correctly set Neumann type boundary condition?

scripting
1

Hello, I am new to Devsim. I tried to make a test problem simulating current flow across 1d device with given conductivity. On one side (contact1), I would like to set the current density entering the device. On the other side (contact2), I set the potential to zero. So, it is a problem with one Neumann and one Dirichlet boundary condition. I get the correct linear function for the potential, but the current density is one half of the required value (set by the contact1 equation). I might have done something wrong in the problem setup, but I cannot figure out what, although I have read the manual and examples thoroughly. Can anybody help me, please? I attach my script.

    dev = 'device'
    reg = 'region'
    msh = 'mesh1'
    length = 1.0
    resistivity = 1.0

    ds.create_1d_mesh(mesh=msh)
    ds.add_1d_mesh_line(mesh=msh, pos=0.0, ps=0.1, tag='contact1')
    ds.add_1d_mesh_line(mesh=msh, pos=length, ps=0.1, tag='contact2')
    ds.add_1d_contact(mesh=msh, name='contact1', tag='contact1', material='metal')
    ds.add_1d_contact(mesh=msh, name='contact2', tag='contact2', material='metal')
    ds.add_1d_region(mesh=msh, material='Si', region=reg, tag1='contact1', tag2='contact2')
    ds.finalize_mesh(mesh=msh)
    ds.create_device(mesh=msh, device=dev)

    ds.set_parameter(device=dev, region=reg, name='sigma', value=1./resistivity)

    # potential variable
    ds.node_solution(device=dev, region=reg, name='V')
    ds.edge_from_node_model(device=dev, region=reg, node_model='V')
    # electric field
    ds.edge_model(device=dev, region=reg, name='E', equation='(V@n0 - V@n1)*EdgeInverseLength')
    ds.edge_model(device=dev, region=reg, name='E:V@n0', equation='EdgeInverseLength')
    ds.edge_model(device=dev, region=reg, name='E:V@n1', equation='-EdgeInverseLength')
    # current density
    ds.edge_model(device=dev, region=reg, name='j', equation='sigma*E')
    ds.edge_model(device=dev, region=reg, name='j:V@n0', equation='sigma*EdgeInverseLength')
    ds.edge_model(device=dev, region=reg, name='j:V@n1', equation='-sigma*EdgeInverseLength')

    # main equation: div(j) = div(-sigma*grad(phi)) = 0   (j = current density, sigma = conductivity, phi = potential)
    ds.equation(device=dev, region=reg, name='PotentialEquation', variable_name='V', edge_model='j',
                variable_update='default')

    # left boundary condition, set incoming current density equal to 1
    ds.set_parameter(device=dev, region=reg, name='current', value=1.0)
    ds.contact_edge_model(device=dev, contact='contact1', name='contact1_bc', equation='j-current')
    ds.contact_edge_model(device=dev, contact='contact1', name='contact1_bc:V@n0', equation='sigma*EdgeInverseLength')
    ds.contact_edge_model(device=dev, contact='contact1', name='contact1_bc:V@n1', equation='-sigma*EdgeInverseLength')

    # right boundary condition, set potential to 0
    ds.contact_node_model(device=dev, contact='contact2', name='contact2_bc', equation='V')
    ds.contact_node_model(device=dev, contact='contact2', name='contact2_bc:V', equation='1')

    # contact equations
    ds.contact_equation(device=dev, contact='contact1', name='PotentialEquation', edge_model='contact1_bc',
                        edge_current_model='j')
    ds.contact_equation(device=dev, contact='contact2', name='PotentialEquation', node_model='contact2_bc',
                        edge_current_model='j')

    ds.solve(type='dc', absolute_error=1.0e-6, relative_error=1.0e-10)

    print(ds.get_node_model_values(device=dev, region=reg, name='V'))
    print(ds.get_edge_model_values(device=dev, region=reg, name='j'))
    print(ds.get_edge_model_values(device=dev, region=reg, name='contact1_bc'))
    print(ds.get_node_model_values(device=dev, region=reg, name='contact2_bc'))

Results:
V = [0.5, 0.45, 0.4, 0.35, 0.3, 0.25, 0.2, 0.15000000000000002, 0.10000000000000003, 0.050000000000000044, 0.0]
j = [0.4999999999999999, 0.4999999999999999, 0.5000000000000002, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]
contact1_bc = [-0.5000000000000001, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
contact2_bc = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

2

Hi @Tomas_Kozak

Thanks for trying the software and sharing your script. The contact boundary conditions should be more like this:

ds.contact_node_model(device=dev, contact='contact1', name='contact1_nbc', equation='-current/NodeVolume')
ds.contact_equation(device=dev, contact='contact1', name='PotentialEquation', edge_model='j', node_model='contact1_nbc',
                    edge_current_model='j')

This is so that the current boundary condition is applied once at each contact node, and is integrated properly. This condition can be tricky at higher dimensions (2d and 3d), and it may be more approriate to incorporate the contact surface area:

ds.contact_node_model(device=dev, contact='contact1', name='contact1_nbc', equation='-ContactSurfaceArea*current/NodeVolume')
3

Hi Juan,
thank you very much for the quick solution. I did not realize that I can combine the ‘j’ edge model and the BC node model in the contact equation. Now it makes perfect sense. I really like the flexibility of Devsim.

Can I have one more question? How would you setup the models and equation to solve the current continuity equation, but with tensor conductivity, i.e. j_i = sigma_ij * E_j. My idea is to solve a Hall effect problem, where an applied magnetic field causes that the current density vector is not parallel with the E field vector. Can this be done in Devsim? I suppose that the Element Edge Models must be used?

Tomas

4

Hi @Tomas,

I have not considered the tensor conductivity problem, but it sounds like the element edge formulation may help. In this paper:
https://ieeexplore.ieee.org/document/9491929
my coauthor considered the ferroelectric effect with a polarization vector which was not aligned with the applied electric field. Please email me at info@devsim.com if you would like a copy of the paper.

Converting the electric field to its vector components, it should be possible to calculate the current.

For your case, it is an interesting problem to consider how the current should be integrated as it flows off of the element edge.

It may be possible to solve the problem with only edge models.
Is the tensor diagonal? Then perhaps you only need the edge ElectricField, and to employ the unitx, unity, unitz edge models to get its vector components directed along each edge.

J = \left( \sigma_x \cdot \textrm{unitx} + \sigma_y \cdot \textrm{unity} + \sigma_z \cdot \textrm{unitz} \right) \cdot \textrm{ElectricField}

where \sigma_x, \sigma_y, \sigma_z are calculated from the magnetic field.