Help modifying 2D Capacitor example

1

I am trying to get into TCAD for a project I’ve been working on, but I’m having some trouble understanding how this works. I decided to try to modify the 2D capacitor example to put a dielectric slab halfway into the capacitor - a classic physics problem with a well-known solution - but no matter what I do I can’t get the right behavior.

I modified the 2D capacitor example to include a new region. I defined some coordinate limits for the new region and then made it:

add_2d_region(mesh=device, material='gas',   region='dielectric' , yl=yl_di  , yh=yh_di  , xl=xl_di  , xh=xh_di)

I made the interface:

add_2d_interface(mesh=device, name='air-dielectric', region0='dielectric', region1='air', yl=yl_di  , yh=yh_di  , xl=xl_di  , xh=xh_di)

I added a new permittivity, right now just using the same permittivity as for air, since that should leave everything unchanged:

set_parameter(device=device, region='dielectric',  name='Permittivity', value=3.9*8.85e-14)

I made a node_solution for the Potential and edge_models for the Electric and Displacement fields in the same way the example made those for the air region.

Finally, the only thing I could come up with that doesn’t give a DEVSIM_FATAL was:

interface_model( device=device, interface='air-dielectric', name='continuousD',  equation='(Permittivity@r0*Potential@r0 - Permittivity@r1*Potential@r1)')
interface_model( device=device, interface='air-dielectric', name='continuousD:Potential@r0', equation= 'Permittivity@r0')
interface_model( device=device, interface='air-dielectric', name='continuousD:Potential@r1', equation='-Permittivity@r1')
interface_equation(device=device, interface='air-dielectric',  name='PotentialEquation', interface_model='continuousD', type='continuous')

The problem is that the printout for get_contact_charge only changes out to the hundredth place when I change the permittivity in the dielectric region, by 1-2 orders of magnitude, simulating a high-k dielectric. Does anyone have any tips for what I should try to do, or maybe can point out some conceptual misunderstanding I have?

2

Hi @tpo

I would think that the proper boundary condition is the continous potential:

interface_model( device=device, interface='air-dielectric', name='continuousP',  equation='(Potential@r0 - Potential@r1)')
interface_model( device=device, interface='air-dielectric', name='continuousP:Potential@r0', equation= '1')
interface_model( device=device, interface='air-dielectric', name='continuousP:Potential@r1', equation='-1')
interface_equation(device=device, interface='air-dielectric',  name='PotentialEquation', interface_model='continuousP', type='continuous')

This interface equation will then replace the equation in region 1 at the nodes in the interface. Region 0 will then have the original equation from Region 1 permutated into its equation.

So at region interface nodes:

Eq0: \epsilon_0 E_0 + \epsilon_1 E_1\\ Eq1: \psi_0 - \psi_1

where the continuous boundary condtion has moved the bulk equation in region 1 to the equation in region 0.

Please feel free to post your full example.

3

I apologize for the imprecise form of the equations, but I hope you can imagine in 1D that the equations would be correct. In 2D there is a surface integration occurring.

\sum_i \epsilon_0 E_{0,i} \text{EdgeCouple}_{0,i} + \sum_j \epsilon_1 E_{1,j} \text{EdgeCouple}_{1,j}

where E_{0,i} is the electric field on an edge between the surface node and the other adjacent nodes in region 0. The EdgeCouple is a discretized form of a the surface integral and is depicted here
https://devsim.net/models.html#id4

4

image
image1270×1127 28.2 KB

5

Hi @Juan

Thanks for your response. Your answers to some of the other questions in this forum have been really helpful for me so far.

However, I don’t understand why the potential being continuous is the right boundary condition, and using that doesn’t give me the expected answer either when I change the relative permittivity in the dielectric from 3.9 to something much larger, like 50. The real boundary conditions at the interface of two dielectrics are really always that the component of the displacement field that is normal to the interface is continuous, (\vec{D}_1-\vec{D}_2)\cdot \hat{n}=0, and that the components of the electric field that are tangential to the interface are continous, (\vec{E}_1-\vec{E}_2)\times \hat{n}=0. See for example Interface conditions for electromagnetic fields - Wikipedia.

Is there a way to enforce those boundary conditions?

Full code 
from devsim import *
device="MyDevice"
region="MyRegion"

    
def cap2d():
    device="MyDevice"

    xmin, ymin  = -100, -50 # window bottom left
    xmax, ymax  = 100,  50  # window top right
    xlb, ylb    = -10,  -0.6 # bottom capacitor plate
    xhb, yhb    = 10,   -0.5
    xlt, ylt    = -10,  0.5  # top capacitor plate 
    xht, yht    = 10,   0.6
    yl_di, yh_di, xl_di, xh_di = yhb+0.2, ylt-0.2, 0*xlb, xhb #corners of dielectric layer

    create_2d_mesh(mesh=device)
    
    add_2d_mesh_line(mesh=device, dir="x", pos=xmin  , ps=5)
    add_2d_mesh_line(mesh=device, dir="x", pos=xmin/4, ps=0.1)
    add_2d_mesh_line(mesh=device, dir="x", pos=xmax/4, ps=0.1)
    add_2d_mesh_line(mesh=device, dir="x", pos=xmax  , ps=5)

    add_2d_mesh_line(mesh=device, dir="y", pos=ymin  , ps=5)
    add_2d_mesh_line(mesh=device, dir="y", pos=ymin/4, ps=0.1)
    add_2d_mesh_line(mesh=device, dir="y", pos=ymax/4, ps=0.1)
    add_2d_mesh_line(mesh=device, dir="y", pos=ymax  , ps=5)

    add_2d_region(mesh=device, material="gas"  , region="air", yl=ymin, yh=ymax, xl=xmin, xh=xmax) #metal sits on the air
    add_2d_region(mesh=device, material="metal", region="mb" , yl=ylb  , yh=yhb  , xl=xlb  , xh=xhb)
    add_2d_region(mesh=device, material="metal", region="mt" , yl=ylt  , yh=yht  , xl=xlt  , xh=xht)
    add_2d_region(mesh=device, material='gas',   region='dielectric' , yl=yl_di  , yh=yh_di  , xl=xl_di  , xh=xh_di)

    # must be air since contacts don't have any equations
    add_2d_contact(mesh=device, name="bot", region="air", material="metal", yl=ylb  , yh=yhb  , xl=xlb  , xh=xhb)
    add_2d_contact(mesh=device, name="top", region="air", material="metal", yl=ylt  , yh=yht  , xl=xlt  , xh=xht)
    add_2d_interface(mesh=device, name='air-dielectric', region0='dielectric', region1='air', yl=yl_di  , yh=yh_di  , xl=xl_di  , xh=xh_di)

    finalize_mesh( mesh=device)
    create_device( mesh=device, device=device)

    ### Set parameters on the region
    set_parameter(device=device, region='air', name="Permittivity", value=3.9*8.85e-14) #SiO2?
    set_parameter(device=device, region='dielectric',  name="Permittivity", value=50*8.85e-14)

    for region in ('air', 'dielectric'):
        ### Create the Potential solution variable
        node_solution(device=device, region=region, name="Potential")

        ### Creates the Potential@n0 and Potential@n1 edge model
        edge_from_node_model(device=device, region=region, node_model="Potential")

        ### Electric field on each edge, as well as its derivatives with respect to the potential at each node
        edge_model(device=device, region=region, name="ElectricField", equation="(Potential@n0 - Potential@n1)*EdgeInverseLength")
        edge_model(device=device, region=region, name="ElectricField:Potential@n0", equation="+EdgeInverseLength")
        edge_model(device=device, region=region, name="ElectricField:Potential@n1", equation="-EdgeInverseLength")


        ### Model the D Field
        edge_model(device=device, region=region, name="DField",              equation="Permittivity*ElectricField")
        edge_model(device=device, region=region, name="DField:Potential@n0", equation="diff(Permittivity*ElectricField, Potential@n0)")
        edge_model(device=device, region=region, name="DField:Potential@n1", equation="-DField:Potential@n0")

        element_from_edge_model(device=device, region=region, edge_model="DField")

        ### Create the bulk equation
        equation(device=device, region=region, name="PotentialEquation",
                variable_name="Potential", edge_model="DField", variable_update="default")
    
    interface_model(   device=device, interface='air-dielectric', name='continuousP', equation='Potential@r0 - Potential@r1')
    interface_model(   device=device, interface='air-dielectric', name='continuousP:Potential@r0', equation= '1')
    interface_model(   device=device, interface='air-dielectric', name='continuousP:Potential@r1', equation='-1')
    interface_equation(device=device, interface='air-dielectric', name='PotentialEquation', interface_model='continuousP', type='continuous')

    ### Contact models and equations
    for c in ("top", "bot"):
        contact_node_model(device=device, contact=c, name=f"{c}_bc", equation=f"Potential - {c}_bias")
        contact_node_model(device=device, contact=c, name=f"{c}_bc:Potential", equation="1")
        contact_equation(  device=device, contact=c, name="PotentialEquation", 
                           node_model=f"{c}_bc", edge_charge_model="DField")

    ### Set the contact
    set_parameter(device=device, name="top_bias", value=1.0e-0)
    set_parameter(device=device, name="bot_bias", value=-1.0e-0)

    edge_model(device=device, region="mb", name="ElectricField", equation="0") #E field is 0 in metal
    edge_model(device=device, region="mt", name="ElectricField", equation="0")
    node_model(device=device, region="mb", name="Potential", equation="bot_bias")
    node_model(device=device, region="mt", name="Potential", equation="top_bias")

    solve(type="dc", absolute_error=1.0, relative_error=1e-6, maximum_iterations=30, solver_type="direct")

    element_from_edge_model(edge_model="ElectricField", device=device, region='air')
    element_from_edge_model(edge_model="DField", device=device, region='air')
    element_from_edge_model(edge_model="ElectricField", device=device, region='dielectric')
    element_from_edge_model(edge_model="DField", device=device, region='dielectric')
    print(get_contact_charge(device=device, contact="top", equation="PotentialEquation"))
    print(get_contact_charge(device=device, contact="bot", equation="PotentialEquation"))

    write_devices(file="cap2d.msh", type="devsim")
    write_devices(file="cap2d.dat", type="vtk")

cap2d()
6

To clarify what I mean by the “expected result”, the unmodified example is just a parallel plate capacitor, so we have generally C=\epsilon_r \epsilon_0 \frac{A}{d}. When a half slab of dielectric material is inserted in between the plates, this becomes two capacitors in parallel, with the same area A and distance d, so the total capacitance is
C=\epsilon_r^{(1)} \epsilon_0 \frac{A}{d} + \epsilon_r^{(2)} \epsilon_0 \frac{A}{d}=(\epsilon_0 \frac{A}{d})(\epsilon_r^{(1)}+\epsilon_r^{(2)}).

Now, we know the relationship between the charge on the plates Q, the capacitance and the voltage difference V, is Q=CV. So when I modify C while keeping V the same, I can determine what I expect the charge to be reported will be. Namely, from the example \epsilon_r^{(1)}=3.9, and if I take initially \epsilon_r^{(2)}=3.9 I get some value Q_i. Then, if I increase \epsilon_r^{(2)}\to39, I expect Q\to(\frac{3.9+39}{3.9+3.9})Q_i=5.5Q_i.

Doing what I just described, what I saw initially was barely any increase in the reported charge from get_contact_charge(...). When I try your suggestion, I see just over a factor of 2 increase. I must be missing something here.

7

If you imagine the Potential to be continuous along the interface. Then the electric field is by the same on both sides. For the normal electric field, it will come just by applying the Poisson equation on both sides of the interface. I remember a simple derivation in a E&M book I no longer have.

8

Please experiment with adding additional space above and below the capacitor. To me, it appears the the field lines are going all the way to the top and the bottom, and that can possibly interfere with the solution.

image
image1461×1177 91.4 KB

9

I can also see the flattening of the field lines to the left and right. Imagine if you would this would be equivalent to having identical capacitors to the left, right, above, below this one.

10

Thank you for taking so much time to help me figure this out!

The potential can be continuous with a discontinuous derivative (meaning a discontinuous electric field), for example at a cusp. This is something that happens all the time in electrostatics.

Also, the electric field line flattening will always happen at the edge of the simulation window, won’t it? The potential is basically zero there though, and so is the electric field, so it should not matter. I tried moving the simulation boundary out from 100 to 150 in all directions, but the charge does not change. It still only increases by ~a factor of 2 when it should be increasing by a factor of 5.5.

11

Can you please make sure you are using the units you expect? You have the permittivity at 8.85e-14 F/cm. The width of your plates are each 20 cm. Is that the scale of the structure you are working with?

12

This is independent of the units because I put it in terms of a ratio, but yes I expect permittivity is in 8.9e-14 F/cm, the plates are 20 cm long with a 1cm gap between them, the voltage difference on the plates is 2 V, and the charge and capacitance would be in C/cm and F/cm because this is a 2D simulation.

Plugging in all the numbers, the charge I should get is ~1.4e-11 C/cm when the extra dielectric region has the same relative permittivity of the air region around, which is 3.9, unaltered from the example. The simulation gives me 1.5e-11 C/cm.

That suggests this simulation is working when there is effectively no dielectric slab, which is good. But when I put in a dielectric slab with a relative permittivity of 39, the answer I should get is 8.25e-11 C/cm. What I actually see is about 3.3e-11 C/cm.

13

I strongly recommend putting meshlines at positions where the regions and contacts begin and end:

    add_2d_mesh_line(mesh=device, dir="x", pos=xmin  , ps=5)
    add_2d_mesh_line(mesh=device, dir="x", pos=xmax  , ps=5)
    add_2d_mesh_line(mesh=device, dir="x", pos=xlt, ps=1)
    add_2d_mesh_line(mesh=device, dir="x", pos=xht, ps=1)
    add_2d_mesh_line(mesh=device, dir="x", pos=xlb, ps=1)
    add_2d_mesh_line(mesh=device, dir="x", pos=xhb, ps=1)
    add_2d_mesh_line(mesh=device, dir="x", pos=xl_di, ps=.1)
    add_2d_mesh_line(mesh=device, dir="x", pos=xh_di, ps=.1)

    add_2d_mesh_line(mesh=device, dir="y", pos=ymin  , ps=5)
    add_2d_mesh_line(mesh=device, dir="y", pos=ymax  , ps=5)
    add_2d_mesh_line(mesh=device, dir="y", pos=ylt, ps=1)
    add_2d_mesh_line(mesh=device, dir="y", pos=yht, ps=1)
    add_2d_mesh_line(mesh=device, dir="y", pos=ylb, ps=1)
    add_2d_mesh_line(mesh=device, dir="y", pos=yhb, ps=1)
    add_2d_mesh_line(mesh=device, dir="y", pos=yl_di, ps=.1)
    add_2d_mesh_line(mesh=device, dir="y", pos=yh_di, ps=.1)

otherwise the structure will not be the same as you are expecting. As far as the calculations are concerned, I would recommend spreading the plates apart in case there is any fringing effect. Please include the equation for the expected value at the end of the script.

14

Hi Juan,

I’ve implemented the mesh as you described, and added a line that prints my expected answer after the answer calculated by devsim. The code is below.

Full Code v2 
device="MyDevice"

xmin, ymin  = -100, -100 # window bottom left
xmax, ymax  = 100,  100  # window top right
xlb, ylb    = -20,  -0.6 # bottom capacitor plate
xhb, yhb    = 20,   -0.5
xlt, ylt    = -20,  0.5  # top capacitor plate 
xht, yht    = 20,   0.6
yl_di, yh_di, xl_di, xh_di = yhb+0.1, ylt-0.1, xlb, xhb #corners of dielectric layer

create_2d_mesh(mesh=device)

add_2d_mesh_line(mesh=device, dir="x", pos=xmin  , ps=5)
add_2d_mesh_line(mesh=device, dir="x", pos=xmax  , ps=5)
add_2d_mesh_line(mesh=device, dir="x", pos=xlt, ps=0.1)
add_2d_mesh_line(mesh=device, dir="x", pos=xht, ps=0.1)
add_2d_mesh_line(mesh=device, dir="x", pos=xlb, ps=0.1)
add_2d_mesh_line(mesh=device, dir="x", pos=xhb, ps=0.1)
add_2d_mesh_line(mesh=device, dir="x", pos=xl_di, ps=0.05)
add_2d_mesh_line(mesh=device, dir="x", pos=xh_di, ps=0.05)

add_2d_mesh_line(mesh=device, dir="y", pos=ymin  , ps=5)
add_2d_mesh_line(mesh=device, dir="y", pos=ymax  , ps=5)
add_2d_mesh_line(mesh=device, dir="y", pos=ylt, ps=0.1)
add_2d_mesh_line(mesh=device, dir="y", pos=yht, ps=0.1)
add_2d_mesh_line(mesh=device, dir="y", pos=ylb, ps=0.1)
add_2d_mesh_line(mesh=device, dir="y", pos=yhb, ps=0.1)
add_2d_mesh_line(mesh=device, dir="y", pos=yl_di, ps=0.05)
add_2d_mesh_line(mesh=device, dir="y", pos=yh_di, ps=0.05)

add_2d_region(mesh=device, material="gas"  , region="air", yl=ymin, yh=ymax, xl=xmin, xh=xmax) #metal sits on the air
add_2d_region(mesh=device, material="metal", region="mb" , yl=ylb  , yh=yhb  , xl=xlb  , xh=xhb)
add_2d_region(mesh=device, material="metal", region="mt" , yl=ylt  , yh=yht  , xl=xlt  , xh=xht)
add_2d_region(mesh=device, material='gas',   region='dielectric' , yl=yl_di  , yh=yh_di  , xl=xl_di  , xh=xh_di)

# must be air since contacts don't have any equations
add_2d_contact(mesh=device, name="bot", region="air", material="metal", yl=ylb  , yh=yhb  , xl=xlb  , xh=xhb)
add_2d_contact(mesh=device, name="top", region="air", material="metal", yl=ylt  , yh=yht  , xl=xlt  , xh=xht)
add_2d_interface(mesh=device, name='air-dielectric', region0='dielectric', region1='air', yl=yl_di  , yh=yh_di  , xl=xl_di  , xh=xh_di)

finalize_mesh( mesh=device)
create_device( mesh=device, device=device)

### Set parameters on the region
set_parameter(device=device, region='air', name="Permittivity", value=3.9*8.85e-14) #SiO2?
set_parameter(device=device, region='dielectric',  name="Permittivity", value=3.9*8.85e-14)

for region in ('air', 'dielectric'):
    ### Create the Potential solution variable
    node_solution(device=device, region=region, name="Potential")

    ### Creates the Potential@n0 and Potential@n1 edge model
    edge_from_node_model(device=device, region=region, node_model="Potential")

    ### Electric field on each edge, as well as its derivatives with respect to the potential at each node
    edge_model(device=device, region=region, name="ElectricField", equation="(Potential@n0 - Potential@n1)*EdgeInverseLength")
    edge_model(device=device, region=region, name="ElectricField:Potential@n0", equation="+EdgeInverseLength")
    edge_model(device=device, region=region, name="ElectricField:Potential@n1", equation="-EdgeInverseLength")


    ### Model the D Field
    edge_model(device=device, region=region, name="DField",              equation="Permittivity*ElectricField")
    edge_model(device=device, region=region, name="DField:Potential@n0", equation="diff(Permittivity*ElectricField, Potential@n0)")
    edge_model(device=device, region=region, name="DField:Potential@n1", equation="-DField:Potential@n0")

    element_from_edge_model(device=device, region=region, edge_model="DField")

    ### Create the bulk equation
    equation(device=device, region=region, name="PotentialEquation",
            variable_name="Potential", edge_model="DField", variable_update="default")

# interface_normal_model(device=device, region=region, interface='air-dielectric')
interface_model(   device=device, interface='air-dielectric', name='continuousP', equation='Potential@r0 - Potential@r1')
interface_model(   device=device, interface='air-dielectric', name='continuousP:Potential@r0', equation='1')
interface_model(   device=device, interface='air-dielectric', name='continuousP:Potential@r1', equation='-1')
interface_equation(device=device, interface='air-dielectric', name='PotentialEquation', interface_model='continuousP', type='continuous')

### Contact models and equations
for c in ("top", "bot"):
    contact_node_model(device=device, contact=c, name=f"{c}_bc", equation=f"Potential - {c}_bias")
    contact_node_model(device=device, contact=c, name=f"{c}_bc:Potential", equation="1")
    contact_equation(  device=device, contact=c, name="PotentialEquation", 
                        node_model=f"{c}_bc", edge_charge_model="DField")

### Set the contact
set_parameter(device=device, name="top_bias", value=1e-0)
set_parameter(device=device, name="bot_bias", value=-1e-0)

edge_model(device=device, region="mb", name="ElectricField", equation="0") #E field is 0 in metal
edge_model(device=device, region="mt", name="ElectricField", equation="0")
node_model(device=device, region="mb", name="Potential", equation="bot_bias")
node_model(device=device, region="mt", name="Potential", equation="top_bias")

solve(type="dc", absolute_error=1.0, relative_error=1e-6, maximum_iterations=30, solver_type="direct")

element_from_edge_model(edge_model="ElectricField", device=device, region='air')
element_from_edge_model(edge_model="DField", device=device, region='air')
element_from_edge_model(edge_model="ElectricField", device=device, region='dielectric')
element_from_edge_model(edge_model="DField", device=device, region='dielectric')
print()
print(f'Calculated charge on top plate: {get_contact_charge(device=device, contact="top", equation="PotentialEquation")} C/cm')

V = (get_parameter(device=device, name="top_bias")-get_parameter(device=device, name="bot_bias"))
C1 = get_parameter(device=device, region="air", name="Permittivity")*(xht-xlt)/(ylt-yhb)/2
C2 = get_parameter(device=device, region="dielectric", name="Permittivity")*(xht-xlt)/(ylt-yhb)/2
print(f'Expected charge on top plate: {(C1+C2)*V} C/cm') # C = \epsilon A/d, Q=CV
# print(get_contact_charge(device=device, contact="bot", equation="PotentialEquation"))

The output when the dielectric constants are the same, 3.9 and 3.9, is

Calculated charge on top plate: 2.9143489870920734e-11 C/cm
Expected charge on top plate: 2.7612e-11 C/cm

Which I’m willing to accept as similar enough (2.8 vs 2.9). The output when the dielectric constants are 3.9 and 39, however, is

Calculated charge on top plate: 1.0021476171746964e-10 C/cm
Expected charge on top plate: 1.5186599999999998e-10 C/cm

Which are pretty different (1 vs 1.5).

15

Hello @tpo

Looking at the mesh file in a viewer, I see that each contact is surrounded by air with a dielectric in the middle. I then calculated the expect charge for Capacitors using:

w=(xht-xlt)
d1=(ylt-yh_di) #between the top plate and dielectric
d2=(yl_di-yhb) #between the bottom plate and dielectric
d3=(yh_di-yl_di) #dielectric
e_a = get_parameter(device=device, region="air", name="Permittivity")
e_d = get_parameter(device=device, region="dielectric", name="Permittivity")
c = [e_a*w/d1, e_a*w/d2, e_d*w/d3]
ctot = 1.0 / sum([1/x for x in c])
#print(f'{d1}, {d2}, {d3}, {ylb}, {yhb}, {ylt}, {yht}, {yl_di}, {yh_di}')
print(ctot*V)

since there are 3 series capacitors, and they are therefore reciprocal summed to give a charge of:

9.861428571428573e-11

which seems close to the simulated value of:

Calculated charge on top plate: 1.0021476171746964e-10 C/cm

and I would be thinking that any additional charge in the simulation may be at least partially attributed to fringe effects.

16

To make it easier, you may consider using only 1 region for insulator. You can then set the permittivity for each edge, based on the midpoint of its x and y position

17

Hi Juan,

I don’t think that you’ve calculated the capacitance correctly - the structure I’ve defined has a half-slab of dielectric, so the full capacitor structure would be:
image

Before, I was assuming that the thin layers of air on the right do not contribute much compared to the high-K dielectric region, but that the other layer of air on the left was non-negligible. Including the thin layers does not fix my problem:

Calculated charge on top plate: 1.0021476171746964e-10 C/cm
Expected charge on top plate: 6.311314285714287e-11 C/cm
eps_die/eps_air=10.0
All capacitances in F/cm:
C_topright=6.903000000000002e-11, C_midright=8.628749999999999e-11, C_botright=6.903000000000002e-11
C_left=6.903e-12, C_right=2.4653571428571433e-11
C_total=3.1556571428571436e-11

Full Code 
from devsim import *
import numpy as np
import matplotlib
matplotlib.use('TkAgg')
from matplotlib import pyplot as plt

device="MyDevice"
region="MyRegion"

    
def cap2d():
    '''
    What is going on here? This forum post clears a lot up
    https://forum.devsim.org/t/potential-and-electrical-field-definition/107
    '''
    device="MyDevice"

    xmin, ymin  = -100, -100 # window bottom left
    xmax, ymax  = 100,  100  # window top right
    xlb, ylb    = -20,  -0.6 # bottom capacitor plate
    xhb, yhb    = 20,   -0.5
    xlt, ylt    = -20,  0.5  # top capacitor plate 
    xht, yht    = 20,   0.6
    yl_di, yh_di, xl_di, xh_di = yhb+0.1, ylt-0.1, xlb, xhb #corners of dielectric layer

    create_2d_mesh(mesh=device)

    add_2d_mesh_line(mesh=device, dir="x", pos=xmin  , ps=5)
    add_2d_mesh_line(mesh=device, dir="x", pos=xmax  , ps=5)
    add_2d_mesh_line(mesh=device, dir="x", pos=xlt, ps=0.1)
    add_2d_mesh_line(mesh=device, dir="x", pos=xht, ps=0.1)
    add_2d_mesh_line(mesh=device, dir="x", pos=xlb, ps=0.1)
    add_2d_mesh_line(mesh=device, dir="x", pos=xhb, ps=0.1)
    add_2d_mesh_line(mesh=device, dir="x", pos=xl_di, ps=0.05)
    add_2d_mesh_line(mesh=device, dir="x", pos=xh_di, ps=0.05)

    add_2d_mesh_line(mesh=device, dir="y", pos=ymin  , ps=5)
    add_2d_mesh_line(mesh=device, dir="y", pos=ymax  , ps=5)
    add_2d_mesh_line(mesh=device, dir="y", pos=ylt, ps=0.1)
    add_2d_mesh_line(mesh=device, dir="y", pos=yht, ps=0.1)
    add_2d_mesh_line(mesh=device, dir="y", pos=ylb, ps=0.1)
    add_2d_mesh_line(mesh=device, dir="y", pos=yhb, ps=0.1)
    add_2d_mesh_line(mesh=device, dir="y", pos=yl_di, ps=0.05)
    add_2d_mesh_line(mesh=device, dir="y", pos=yh_di, ps=0.05)

    add_2d_region(mesh=device, material="gas"  , region="air", yl=ymin, yh=ymax, xl=xmin, xh=xmax) #metal sits on the air
    add_2d_region(mesh=device, material="metal", region="mb" , yl=ylb  , yh=yhb  , xl=xlb  , xh=xhb)
    add_2d_region(mesh=device, material="metal", region="mt" , yl=ylt  , yh=yht  , xl=xlt  , xh=xht)
    add_2d_region(mesh=device, material='gas',   region='dielectric' , yl=yl_di  , yh=yh_di  , xl=xl_di  , xh=xh_di)

    # must be air since contacts don't have any equations
    add_2d_contact(mesh=device, name="bot", region="air", material="metal", yl=ylb  , yh=yhb  , xl=xlb  , xh=xhb)
    add_2d_contact(mesh=device, name="top", region="air", material="metal", yl=ylt  , yh=yht  , xl=xlt  , xh=xht)
    add_2d_interface(mesh=device, name='air-dielectric', region0='dielectric', region1='air', yl=yl_di  , yh=yh_di  , xl=xl_di  , xh=xh_di)

    finalize_mesh( mesh=device)
    create_device( mesh=device, device=device)

    ### Set parameters on the region
    set_parameter(device=device, region='air', name="Permittivity", value=3.9*8.85e-14) #SiO2?
    set_parameter(device=device, region='dielectric',  name="Permittivity", value=39*8.85e-14)

    for region in ('air', 'dielectric'):
        ### Create the Potential solution variable
        node_solution(device=device, region=region, name="Potential")

        ### Creates the Potential@n0 and Potential@n1 edge model
        edge_from_node_model(device=device, region=region, node_model="Potential")

        ### Electric field on each edge, as well as its derivatives with respect to the potential at each node
        edge_model(device=device, region=region, name="ElectricField", equation="(Potential@n0 - Potential@n1)*EdgeInverseLength")
        edge_model(device=device, region=region, name="ElectricField:Potential@n0", equation="+EdgeInverseLength")
        edge_model(device=device, region=region, name="ElectricField:Potential@n1", equation="-EdgeInverseLength")


        ### Model the D Field
        edge_model(device=device, region=region, name="DField",              equation="Permittivity*ElectricField")
        edge_model(device=device, region=region, name="DField:Potential@n0", equation="diff(Permittivity*ElectricField, Potential@n0)")
        edge_model(device=device, region=region, name="DField:Potential@n1", equation="-DField:Potential@n0")

        element_from_edge_model(device=device, region=region, edge_model="DField")

        ### Create the bulk equation
        equation(device=device, region=region, name="PotentialEquation",
                variable_name="Potential", edge_model="DField", variable_update="default")

    # interface_normal_model(device=device, region=region, interface='air-dielectric')
    interface_model(   device=device, interface='air-dielectric', name='continuousP', equation='Potential@r0 - Potential@r1')
    interface_model(   device=device, interface='air-dielectric', name='continuousP:Potential@r0', equation='1')
    interface_model(   device=device, interface='air-dielectric', name='continuousP:Potential@r1', equation='-1')
    interface_equation(device=device, interface='air-dielectric', name='PotentialEquation', interface_model='continuousP', type='continuous')

    ### Contact models and equations
    for c in ("top", "bot"):
        contact_node_model(device=device, contact=c, name=f"{c}_bc", equation=f"Potential - {c}_bias")
        contact_node_model(device=device, contact=c, name=f"{c}_bc:Potential", equation="1")
        contact_equation(  device=device, contact=c, name="PotentialEquation", 
                            node_model=f"{c}_bc", edge_charge_model="DField")

    ### Set the contact
    set_parameter(device=device, name="top_bias", value=1e-0)
    set_parameter(device=device, name="bot_bias", value=-1e-0)

    edge_model(device=device, region="mb", name="ElectricField", equation="0") #E field is 0 in metal
    edge_model(device=device, region="mt", name="ElectricField", equation="0")
    node_model(device=device, region="mb", name="Potential", equation="bot_bias")
    node_model(device=device, region="mt", name="Potential", equation="top_bias")

    solve(type="dc", absolute_error=1.0, relative_error=1e-6, maximum_iterations=30, solver_type="direct")

    element_from_edge_model(edge_model="ElectricField", device=device, region='air')
    element_from_edge_model(edge_model="DField", device=device, region='air')
    element_from_edge_model(edge_model="ElectricField", device=device, region='dielectric')
    element_from_edge_model(edge_model="DField", device=device, region='dielectric')
    print()
    print(f'Calculated charge on top plate: {get_contact_charge(device=device, contact="top", equation="PotentialEquation")} C/cm')

    V = (get_parameter(device=device, name="top_bias")-get_parameter(device=device, name="bot_bias"))
    eps_air = get_parameter(device=device, region="air", name="Permittivity")
    eps_die = get_parameter(device=device, region="dielectric", name="Permittivity")
    w = (xht - xlt)/2 #half-width of the full capacitor structure
    d_tb = ylt - yhb # length of gap between top and bottom plate
    d_td = ylt - yh_di # " " top plate and top of dielectric
    d_d  = yh_di - yl_di # " " top of dielectric and bottom of dielectric
    d_db = yl_di - yhb # " " bottom of dielectric and bottom plate

    C_left = eps_air*w/d_tb #left capacitor is just air
    C_topright = eps_air*w/d_td #top right capacitor, thin air layer
    C_midright = eps_die*w/d_d #middle right capacitor, dielectric layer
    C_botright = eps_air*w/d_db #bottom right capacitor, another thin air layer
    C_right = ( C_topright**-1 + C_midright**-1 + C_botright**-1 )**-1 # right capacitor is air-dielectric-air
    C_total = C_left + C_right #total capacitance

    print(f'Expected charge on top plate: {C_total*V} C/cm') # C = \epsilon A/d, Q=CV
    print(f'eps_die/eps_air={eps_die/eps_air}')
    print('All capacitances in F/cm:')
    print(f'C_topright={C_topright}, C_midright={C_midright}, C_botright={C_botright}')
    print(f'C_left={C_left}, C_right={C_right}')
    print(f'C_total={C_total}')
    # print(get_contact_charge(device=device, contact="bot", equation="PotentialEquation"))

cap2d()

I really feel that the boundary condition is not correct. Can you explain your reasoning? There should be a charge layer at the dielectric interfaces, so Poisson’s equation should have a charge layer there.

18

Hi @tpo

For the previous code you sent me:

What I am seeing is one dielectric slab sandwiched between the 2 plates, as shown in the plot figure.The dark blue being contacts, the light blue being air, and the red being dielectric.

image
image1698×1117 29.1 KB

I don’t see a half slab. Perhaps you did not mean to set:

xl_di = xlb
xh_di = xhb

If I set:

xl_di = 0

Then I see:

Calculated charge on top plate: 6.478672990772396e-11 C/cm
Expected charge on top plate: 6.311314285714287e-11 C/cm

If you want charge at the interfaces, you would need to add it to the Poisson’s equation.