Simulation of simple 2D SiC diode

Hi Juan and everyone!

I have to say that I’m a big fan of Devsim. It has allowed me to start learning about TCAD. Currently, I’m working on simulating a simple 2D SiC diode and I’ve encountered some issues.

I need to simulate the diode under reverse and high voltage conditions, but I’m encountering a devsim_py3.error: Convergence failure! when ramping up the voltage.

Here is my code:

from devsim import (
    create_gmsh_mesh,
    add_gmsh_region,
    add_gmsh_contact,
    add_gmsh_interface,
    finalize_mesh,
    create_device,
    write_devices,
    node_model,
    get_region_list,
    get_contact_list,
    set_parameter,
    get_interface_list,
    solve,
    element_from_edge_model,
    set_node_values,
    get_contact_current,
    reset_devsim
    )

from devsim.python_packages.model_create import (
    CreateSolution,
    )

from devsim.python_packages.simple_physics import(
    CreateSiliconPotentialOnly,
    CreateSiliconPotentialOnlyContact,
    GetContactBiasName,
    CreateSiliconOxideInterface,
    CreateSiliconDriftDiffusion,
    CreateSiliconDriftDiffusionAtContact,
    PrintCurrents,
    )

import matplotlib.pyplot as plt

def ExtendSolver(x):
    if x is True:
        set_parameter(name = "extended_solver", value=True)
        set_parameter(name = "extended_model", value=True)
        set_parameter(name = "extended_equation", value=True)
        print("Solver Extended: True")
    else:
        print("Solver Extended: False")

ExtendSolver(True)

name="JBS_Diode_Lite"

def CreateMesh(name):
    create_gmsh_mesh(file=f"{name}.msh", mesh=name)

    add_gmsh_region(mesh=name, gmsh_name="P", region="P", material= "SiC")
    add_gmsh_region(mesh=name, gmsh_name="N", region= "N", material="SiC")

    add_gmsh_contact(mesh=name, gmsh_name= "CP", region="P", name="TOP", material= "Metal")
    add_gmsh_contact(mesh=name, gmsh_name= "CN", region="N", name="BOT", material= "Metal")

    add_gmsh_interface(mesh=name, gmsh_name="I1", region0= "P", region1= "N", name="I1")

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

CreateMesh(name)

P_doping= "1e16"
N_doping= "2e16"

def Doping(P_doping, N_doping):
    node_model(name="Acceptors", device=name, region="P", equation=P_doping)
    node_model(name="Donors", device=name, region="N", equation=N_doping)

    node_model(name="NetDoping", device=name, region="P",equation="-Acceptors")
    node_model(name="NetDoping", device=name, region="N", equation="Donors")

Doping(P_doping, N_doping)

Si_dic={ 
        "n_i": 1.18e-18,
         "mu_n": 950,
         "mu_p": 125,
         "taun": 2.5e-6 ,
         "taup": 0.5e-6 ,
         "eps": 9.66
         }

def SetSiCParameters(device, region, dictonary):
    """
    Sets physical parameters assuming constants
    """
    x=dictonary
    eps_0 = 8.85e-14
    k = 1.3806503e-23
    T=300 
    q= 1.6e-19 
    
    set_parameter(device=device, region=region, name="Permittivity", value=x["eps"] * eps_0)
    set_parameter(device=device, region=region, name="ElectronCharge", value=q)
    set_parameter(device=device, region=region, name="n_i", value=x["n_i"])
    set_parameter(device=device, region=region, name="T", value=T)
    set_parameter(device=device, region=region, name="kT", value=k * T)
    set_parameter(device=device, region=region, name="V_t", value=k * T / q)
    set_parameter(device=device, region=region, name="mu_n", value=x["mu_n"])
    set_parameter(device=device, region=region, name="mu_p", value=x["mu_p"])
    set_parameter(device=device, region=region, name="n1", value=x["n_i"])
    set_parameter(device=device, region=region, name="p1", value=x["n_i"])
    set_parameter(device=device, region=region, name="taun", value=x["taun"])
    set_parameter(device=device, region=region, name="taup", value=x["taup"])

regions= get_region_list(device=name)
contacts= get_contact_list(device=name)
interfaces= get_interface_list(device=name)

def InitialSolution(name):
    
    for i in regions:
        CreateSolution(name, i, "Potential")
        SetSiCParameters(name, i, Si_dic)
        CreateSiliconPotentialOnly(name, i)
                  
    for j in contacts:
        tmp=get_region_list(device=name, contact= j)
        r= tmp[0]
        CreateSiliconPotentialOnlyContact(name, r, j)
        set_parameter(device=name, name=GetContactBiasName(j), value=0)
  
    for k in interfaces:
        CreateSiliconOxideInterface(name, k)

    solve(type="dc", absolute_error=1.0e30, relative_error=1e-5, maximum_iterations=30)
    write_devices(file=f"{name}.msh", type="vtk")

def DriftDiff(device):
    
    regions= get_region_list(device=name)
    contacts= get_contact_list(device=name)
    interfaces= get_interface_list(device=name)
    
    for i in regions:
        CreateSolution(device, i, "Electrons")
        CreateSolution(device, i, "Holes")
        set_node_values(
            device=device, region=i, name="Electrons", init_from="IntrinsicElectrons"
        )
        set_node_values(device=device, region=i, name="Holes", init_from="IntrinsicHoles")
        CreateSiliconDriftDiffusion(device, i, "mu_n", "mu_p")

    for c in contacts:
        tmp = get_region_list(device=device, contact=c)
        r = tmp[0]
        CreateSiliconDriftDiffusionAtContact(device, r, c)

    solve(type="dc", absolute_error=1.0e30, relative_error=1e-5, maximum_iterations=30)

    for r in regions:
        node_model(
            device=device, region=r, name="logElectrons", equation="log(Electrons)/log(10)")

    for r in regions:
        element_from_edge_model(edge_model="ElectricField", device=device, region=r)
        element_from_edge_model(edge_model="ElectronCurrent", device=device, region=r)
        element_from_edge_model(edge_model="HoleCurrent", device=device, region=r)

    solve(type="dc", absolute_error=1.0e30, relative_error=1e-5, maximum_iterations=30)    

    write_devices(file=f"{name}.msh", type="vtk")

import keyboard

def IV(name, contact, voltage, step):
    
    v=[]
    I=[]
    
    while True:
    
        set_parameter(device=name, name= GetContactBiasName(contact), value=0-voltage)
        solve(type="dc", absolute_error=1.0e30, relative_error=1e-5, maximum_iterations=30)
        PrintCurrents(name, contact)
        
        I_electron= get_contact_current(device=name, contact=contact, equation="ElectronContinuityEquation")
        I_hole= get_contact_current(device=name, contact=contact, equation="HoleContinuityEquation")
        I_all=I_electron+I_hole
        
        v.append(voltage)
        I.append(I_all)
        print(voltage)
        if abs(I_all) > 1e-3:  # Przykładowy próg przebicia
            print(f"Napięcie przebicia: {v} V")
            break

        voltage+=step
        
    plt.figure()
    plt.title(contact)
    plt.xlabel("Voltage (V)")
    plt.ylabel("Current (A)")
    plt.plot(v,I)
    
    plt.savefig(contact, format= "png" ) 

InitialSolution(name)
DriftDiff(name)
IV(name, "TOP", 100, 0.1)
plt.show()
reset_devsim()

I’m new to TCAD and Devsim software, so any tips would be greatly appreciated.

Hi @Kacper_J,

Thank you for trying the software. For a Silicon simulation, 100mV is what I think would be a very large voltage step. For SiC simulation I am not sure how the situation would be much different.

Please try smaller steps, such as 0.005, 0.010, or 0.025 to see if that improves the how far you get in your simulation.

There is a ramp_bias function, which is catches convergence failure, and attempts to reduce the voltage step. It is used in some of the devsim examples, that I should be able to help you with further if the previous experiments helps in your bias sweep.

I am curious about your use of CreateSiliconOxideInterface. Once you switch to the Drift Diffusion Simulation, you would need to make sure that the interface equations are handled properly for electrons and holes.

Please consider CreateSiliconSiliconInterface function as an example showing how to allow Electrons and Holes to flow across the material interface in your example.

Alternatively, you may want to consider using a single region for both the p and n doped areas of the device. Then you can control the doping profile just using step or other function.