Hello everyone, I am trying to calculate the Joule heat only and its influence to temperature distribution in the MOSFET. For now, there are two problems bothering me. The equation is :
First, the heat conduction equation has same form with Potential equation. So I modified the Poisson equation to get the heat conduction equation. For simplicity, I use the uniform heat source to test the code. But the results shows that the heat source didn’t work. My modified code is shown below：
The gmsh_mos2d.py :
# Copyright 2013 Devsim LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
from devsim.python_packages.simple_physics import *
from devsim.python_packages.ramp import *
import gmsh_mos2d_create
device = "mos2d"
silicon_regions=("gate", "bulk")
oxide_regions=("oxide",)
regions = ("gate", "bulk", "oxide")
interfaces = ("bulk_oxide", "gate_oxide")
for i in regions:
CreateSolution(device, i, "Potential")
for i in silicon_regions:
SetSiliconParameters(device, i, 300)
CreateSiliconPotentialOnly(device, i)
for i in oxide_regions:
SetOxideParameters(device, i, 300)
CreateOxidePotentialOnly(device, i, "log_damp")
### Set up contacts
contacts = get_contact_list(device=device)
for i in contacts:
tmp = get_region_list(device=device, contact=i)
r = tmp[0]
print("%s %s" % (r, i))
CreateSiliconPotentialOnlyContact(device, r, i)
# set_parameter(device=device, name=GetContactBiasName(i), value=0.0)
set_parameter(device=device, name="drain_bias", value=0.0)
set_parameter(device=device, name="source_bias", value=0.0)
set_parameter(device=device, name="gate_bias", value=0.0)
set_parameter(device=device, name="body_bias", value=0.0)
for i in interfaces:
CreateSiliconOxideInterface(device, i)
# for d in get_device_list():
# for gn in get_parameter_list():
# print("{0} {1}").format(gn, get_parameter(device=d, name=gn))
# for gn in get_parameter_list(device=d):
# print("{0} {1} {2}").format(d, gn, get_parameter(device=d, nasme=gn))
# for r in get_region_list(device=d):
# for gn in get_parameter_list(device=d, region=r):
# print("{0} {1} {2} {3}").format(d, r, gn, get_parameter(device=d, region=r, name=gn))
# write_devices(file="foo.msh", type="devsim")
#
CreateSolution(device, "bulk", "Temperature")
CreateTemperaturefiled(device,"bulk")
set_node_values(device=device, region="bulk", name="Temperature", init_from="init_temperature")
CreateThermalcontact(device,"bulk","body")
CreateThermalcontact(device,"bulk","drain")
CreateThermalcontact(device,"bulk","source")
set_parameter(device=device, name="drain_biasT", value=300)
set_parameter(device=device, name="source_biasT", value=300)
set_parameter(device=device, name="body_biasT", value=300)
solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=30)
# solve(type="dc", absolute_error=1.0e-13, relative_error=1e-12, maximum_iterations=30)
#
# write_devices -file gmsh_mos2d_potentialonly.flps -type floops
# write_devices(file="gmsh_mos2d_potentialonly", type="vtk")
for i in silicon_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 silicon_regions:
node_model(device=device, region=r, name="logElectrons", equation="log(Electrons)/log(10)")
##write_devices -file gmsh_mos2d_dd.flps -type floops
##write_devices -file gmsh_mos2d_dd -type vtk
##write_devices -file gmsh_mos2d_dd.msh -type devsim_data
for r in silicon_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)
#
# printAllCurrents(device)
# rampbias(device, "drain", 1, 0.01, 0.001, 100, 1e-10, 1e30, printAllCurrents) ##current unit is A/cm
# rampbias(device, "gate", 1, 0.01, 0.001, 100, 1e-10, 1e30, printAllCurrents)
#
write_devices(file="gmsh_mos2d_dd.dat", type="tecplot")
The simple_physics.py:
# Copyright 2013 Devsim LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
from .simple_dd import *
from devsim import *
#TODO: make this a class so that paramters can be changed
contactcharge_node="contactcharge_node"
contactcharge_edge="contactcharge_edge"
ece_name="ElectronContinuityEquation"
hce_name="HoleContinuityEquation"
celec_model = "(1e-10 + 0.5*abs(NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"
chole_model = "(1e-10 + 0.5*abs(-NetDoping+(NetDoping^2 + 4 * n_i^2)^(0.5)))"
q = 1.6e-19 # coul
k = 1.3806503e-23 # J/K
eps_0 = 8.85e-14 # F/cm^2
#T = 300 # K
eps_si = 11.1
eps_ox = 3.9
# TODO: make this temperature dependent
n_i = 1.0e10 # #/cm^3
# constant in our approximation
mu_n = 400
mu_p = 200
Sithermalconducitivity = 19100 ##W/K*cm
heatsource = 10000 ##w/cm^2
init_T = 300
boundarytemperature = 250
def GetContactBiasName(contact):
return "{0}_bias".format(contact)
def GetContactTemperatureName(contact):
return "{0}_biasT".format(contact)
def GetContactNodeModelName(contact):
return "{0}nodemodel".format(contact)
def GetContactTemperatureNodeModelName(contact):
return "{0}Tnodemodel".format(contact)
def PrintCurrents(device, contact):
'''
print out contact currents
'''
# TODO add charge
contact_bias_name = GetContactBiasName(contact)
electron_current= get_contact_current(device=device, contact=contact, equation=ece_name)
hole_current = get_contact_current(device=device, contact=contact, equation=hce_name)
total_current = electron_current + hole_current
voltage = get_parameter(device=device, name=GetContactBiasName(contact))
print("{0}\t{1}\t{2}\t{3}\t{4}".format(contact, voltage, electron_current, hole_current, total_current))
#in the future, worry about workfunction
def CreateOxideContact(device, region, contact):
conteq="Permittivity*ElectricField"
contact_bias_name = GetContactBiasName(contact)
contact_model_name = GetContactNodeModelName(contact)
eq = "Potential - {0}".format(contact_bias_name)
CreateContactNodeModel(device, contact, contact_model_name, eq)
CreateContactNodeModelDerivative(device, contact, contact_model_name, eq, "Potential")
#TODO: make everyone use dfield
if not InEdgeModelList(device, region, contactcharge_edge):
CreateEdgeModel(device, region, contactcharge_edge, "Permittivity*ElectricField")
CreateEdgeModelDerivatives(device, region, contactcharge_edge, "Permittivity*ElectricField", "Potential")
contact_equation(device=device, contact=contact, name="PotentialEquation",
node_model=contact_model_name, edge_charge_model= contactcharge_edge)
#####
##### Constants
#####
def SetOxideParameters(device, region, T):
'''
Sets physical parameters
'''
set_parameter(device=device, region=region, name="Permittivity", value=eps_ox * eps_0)
set_parameter(device=device, region=region, name="ElectronCharge", value=q)
def SetSiliconParameters(device, region, T):
'''
Sets physical parameters assuming constants
'''
#### TODO: make T a free parameter and T dependent parameters as models
set_parameter(device=device, region=region, name="Permittivity", value=eps_si * eps_0)
set_parameter(device=device, region=region, name="ElectronCharge", value=q)
set_parameter(device=device, region=region, name="n_i", value=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=mu_n)
set_parameter(device=device, region=region, name="mu_p", value=mu_p)
set_parameter(device=device, region=region, name="thermalconducitivity", value=Sithermalconducitivity)
set_parameter(device=device, region=region, name="heatsource", value=heatsource)
set_parameter(device=device, region=region, name="initialTemperature", value=init_T)
set_parameter(device=device, region=region, name="boundaryTemperature", value=boundarytemperature)
#default SRH parameters
set_parameter(device=device, region=region, name="n1", value=n_i)
set_parameter(device=device, region=region, name="p1", value=n_i)
set_parameter(device=device, region=region, name="taun", value=1e-5)
set_parameter(device=device, region=region, name="taup", value=1e-5)
def CreateSiliconPotentialOnly(device, region):
'''
Creates the physical models for a Silicon region
'''
if not InNodeModelList(device, region, "Potential"):
print("Creating Node Solution Potential")
CreateSolution(device, region, "Potential")
elec_i = "n_i*exp(Potential/V_t)"
hole_i = "n_i^2/IntrinsicElectrons"
charge_i = "kahan3(IntrinsicHoles, -IntrinsicElectrons, NetDoping)"
pcharge_i = "-ElectronCharge * IntrinsicCharge"
# require NetDoping
for i in (
("IntrinsicElectrons", elec_i),
("IntrinsicHoles", hole_i),
("IntrinsicCharge", charge_i),
("PotentialIntrinsicCharge", pcharge_i)
):
n = i[0]
e = i[1]
CreateNodeModel(device, region, n, e)
CreateNodeModelDerivative(device, region, n, e, "Potential")
### TODO: Edge Average Model
for i in (
("ElectricField", "(Potential@n0-Potential@n1)*EdgeInverseLength"),
("PotentialEdgeFlux", "Permittivity * ElectricField")
):
n = i[0]
e = i[1]
CreateEdgeModel(device, region, n, e)
CreateEdgeModelDerivatives(device, region, n, e, "Potential")
equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
node_model="PotentialIntrinsicCharge", edge_model="PotentialEdgeFlux", variable_update="log_damp")
def CreateSiliconPotentialOnlyContact(device, region, contact, is_circuit=False):
'''
Creates the potential equation at the contact
if is_circuit is true, than use node given by GetContactBiasName
'''
# Means of determining contact charge
# Same for all contacts
if not InNodeModelList(device, region, "contactcharge_node"):
CreateNodeModel(device, region, "contactcharge_node", "ElectronCharge*IntrinsicCharge")
#### TODO: This is the same as D-Field
if not InEdgeModelList(device, region, "contactcharge_edge"):
CreateEdgeModel(device, region, "contactcharge_edge", "Permittivity*ElectricField")
CreateEdgeModelDerivatives(device, region, "contactcharge_edge", "Permittivity*ElectricField", "Potential")
# set_parameter(device=device, region=region, name=GetContactBiasName(contact), value=0.0)
contact_model = "Potential -{0} + ifelse(NetDoping > 0, \
-V_t*log({1}/n_i), \
V_t*log({2}/n_i))".format(GetContactBiasName(contact), celec_model, chole_model)
contact_model_name = GetContactNodeModelName(contact)
CreateContactNodeModel(device, contact, contact_model_name, contact_model)
# Simplify it too complicated
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,"Potential"), "1")
if is_circuit:
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,GetContactBiasName(contact)), "-1")
if is_circuit:
contact_equation(device=device, contact=contact, name="PotentialEquation",
node_model=contact_model_name, edge_model="",
node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
node_current_model="", edge_current_model="", circuit_node=GetContactBiasName(contact))
else:
contact_equation(device=device, contact=contact, name="PotentialEquation",
node_model=contact_model_name, edge_model="",
node_charge_model="contactcharge_node", edge_charge_model="contactcharge_edge",
node_current_model="", edge_current_model="")
def CreateSRH(device, region):
USRH="(Electrons*Holes - n_i^2)/(taup*(Electrons + n1) + taun*(Holes + p1))"
Gn = "-ElectronCharge * USRH"
Gp = "+ElectronCharge * USRH"
CreateNodeModel(device, region, "USRH", USRH)
CreateNodeModel(device, region, "ElectronGeneration", Gn)
CreateNodeModel(device, region, "HoleGeneration", Gp)
for i in ("Electrons", "Holes"):
CreateNodeModelDerivative(device, region, "USRH", USRH, i)
CreateNodeModelDerivative(device, region, "ElectronGeneration", Gn, i)
CreateNodeModelDerivative(device, region, "HoleGeneration", Gp, i)
def CreateECE(device, region, mu_n):
CreateElectronCurrent(device, region, mu_n)
NCharge = "-ElectronCharge * Electrons"
CreateNodeModel(device, region, "NCharge", NCharge)
CreateNodeModelDerivative(device, region, "NCharge", NCharge, "Electrons")
equation(device=device, region=region, name="ElectronContinuityEquation", variable_name="Electrons",
time_node_model = "NCharge",
edge_model="ElectronCurrent", variable_update="positive", node_model="ElectronGeneration")
def CreateHCE(device, region, mu_p):
CreateHoleCurrent(device, region, mu_p)
PCharge = "ElectronCharge * Holes"
CreateNodeModel(device, region, "PCharge", PCharge)
CreateNodeModelDerivative(device, region, "PCharge", PCharge, "Holes")
equation(device=device, region=region, name="HoleContinuityEquation", variable_name="Holes",
time_node_model = "PCharge",
edge_model="HoleCurrent", variable_update="positive", node_model="HoleGeneration")
def CreatePE(device, region):
pne = "-ElectronCharge*kahan3(Holes, -Electrons, NetDoping)"
CreateNodeModel(device, region, "PotentialNodeCharge", pne)
CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Electrons")
CreateNodeModelDerivative(device, region, "PotentialNodeCharge", pne, "Holes")
equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
node_model="PotentialNodeCharge", edge_model="PotentialEdgeFlux",
time_node_model="", variable_update="log_damp")
def CreateSiliconDriftDiffusion(device, region, mu_n="mu_n", mu_p="mu_p"):
CreatePE(device, region)
CreateBernoulli(device, region)
CreateSRH(device, region)
CreateECE(device, region, mu_n)
CreateHCE(device, region, mu_p)
def CreateSiliconDriftDiffusionAtContact(device, region, contact, is_circuit=False):
'''
Restrict electrons and holes to their equilibrium values
Integrates current into circuit
'''
contact_electrons_model = "Electrons - ifelse(NetDoping > 0, {0}, n_i^2/{1})".format(celec_model, chole_model)
contact_holes_model = "Holes - ifelse(NetDoping < 0, +{1}, +n_i^2/{0})".format(celec_model, chole_model)
contact_electrons_name = "{0}nodeelectrons".format(contact)
contact_holes_name = "{0}nodeholes".format(contact)
CreateContactNodeModel(device, contact, contact_electrons_name, contact_electrons_model)
#TODO: The simplification of the ifelse statement is time consuming
# CreateContactNodeModelDerivative(device, contact, contact_electrons_name, contact_electrons_model, "Electrons")
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_electrons_name, "Electrons"), "1")
CreateContactNodeModel(device, contact, contact_holes_name, contact_holes_model)
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_holes_name, "Holes"), "1")
#TODO: keyword args
if is_circuit:
contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
node_model=contact_electrons_name,
edge_current_model="ElectronCurrent", circuit_node=GetContactBiasName(contact))
contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
node_model=contact_holes_name,
edge_current_model="HoleCurrent", circuit_node=GetContactBiasName(contact))
else:
contact_equation(device=device, contact=contact, name="ElectronContinuityEquation",
node_model=contact_electrons_name,
edge_current_model="ElectronCurrent")
contact_equation(device=device, contact=contact, name="HoleContinuityEquation",
node_model=contact_holes_name,
edge_current_model="HoleCurrent")
def CreateOxidePotentialOnly(device, region, update_type="default"):
'''
Create electric field model in oxide
Creates Potential solution variable if not available
'''
if not InNodeModelList(device, region, "Potential"):
print("Creating Node Solution Potential")
CreateSolution(device, region, "Potential")
efield="(Potential@n0 - Potential@n1)*EdgeInverseLength"
# this needs to remove derivatives w.r.t. independents
CreateEdgeModel(device, region, "ElectricField", efield)
CreateEdgeModelDerivatives(device, region, "ElectricField", efield, "Potential")
dfield="Permittivity*ElectricField"
CreateEdgeModel(device, region, "PotentialEdgeFlux", dfield)
CreateEdgeModelDerivatives(device, region, "PotentialEdgeFlux", dfield, "Potential")
equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
edge_model="PotentialEdgeFlux", variable_update=update_type)
def CreateSiliconOxideInterface(device, interface):
'''
continuous potential at interface
'''
model_name = CreateContinuousInterfaceModel(device, interface, "Potential")
interface_equation(device=device, interface=interface, name="PotentialEquation", interface_model=model_name, type="continuous")
#
##TODO: similar model for silicon/silicon interface
## should use quasi-fermi potential
def CreateSiliconSiliconInterface(device, interface):
'''
Enforces potential, electron, and hole continuity across the interface
'''
CreateSiliconOxideInterface(device, interface)
ename = CreateContinuousInterfaceModel(device, interface, "Electrons")
interface_equation(device=device, interface=interface, name="ElectronContinuityEquation", interface_model=ename, type="continuous")
hname = CreateContinuousInterfaceModel(device, interface, "Holes")
interface_equation(device=device, interface=interface, name="HoleContinuityEquation", interface_model=hname, type="continuous")
def CreateTemperaturefiled(device, region):
'''
Creates a thermal Equation in a region
'''
if not InNodeModelList(device, region, "Temperature"):
print("Creating Node Solution Temperature")
CreateSolution(device, region, "Temperature")
for i in (
("Source", "heatsource"),
("init_temperature", "initialTemperature"),
):
n = i[0]
e = i[1]
CreateNodeModel(device, region, n, e)
CreateNodeModelDerivative(device, region, n, e, "Temperature")
### TODO: Edge Average Model
for i in (
("Temperaturefield", "(Temperature@n0-Temperature@n1)*EdgeInverseLength"),
("ThermalEdgeflux", "thermalconducitivity * Temperaturefield")
):
n = i[0]
e = i[1]
CreateEdgeModel(device, region, n, e)
CreateEdgeModelDerivatives(device, region, n, e, "Temperature")
equation(device=device, region=region, name="ThermalEquation", variable_name="Temperature",
node_model="Source", edge_model="ThermalEdgeflux", variable_update="log_damp")
def CreateThermalcontact(device, region, contact): ##constant temperature boundary
# Boundary_name_temperature="%sbiasT" % contact
# format_dict = { "contact" : contact}
# set_parameter(device=device, region=region, name=Boundary_name_temperature, value=300)
# for name, equation in (
# ("%(contact)snodemodel", "Temperature-%(contact)sbiasT"),
# ("%(contact)snodemodel:Potential", "1"),
# ):
# name_sub = name % format_dict
# equation_sub = equation % format_dict
# contact_node_model(device=device, contact=contact, name=name_sub, equation=equation_sub)
# contact_equation(device=device, contact=contact, name="ThermalEquation",
# node_model="%snodemodel" % contact)
if not InNodeModelList(device, region, "contacttemperature_node"):
CreateNodeModel(device, region, "contacttemperature_node", "-heatsource")
#### TODO: This is the same as D-Field
if not InEdgeModelList(device, region, "contacttemperature_edge"):
CreateEdgeModel(device, region, "contacttemperature_edge", "thermalconducitivity*Temperaturefield")
CreateEdgeModelDerivatives(device, region, "contacttemperature_edge", "thermalconducitivity*Temperaturefield", "Temperature")
contact_model = "Temperature -{0}".format(GetContactTemperatureName(contact))
contact_model_name = GetContactTemperatureNodeModelName(contact)
CreateContactNodeModel(device, contact, contact_model_name, contact_model)
# Simplify it too complicated
CreateContactNodeModel(device, contact, "{0}:{1}".format(contact_model_name,"Temperature"), "1")
contact_equation(device=device, contact=contact, name="ThermalEquation",
node_model=contact_model_name, edge_model="",
node_charge_model="contacttemperature_node", edge_charge_model="contacttemperature_edge",
node_current_model="", edge_current_model="")
I didn’t change the simple_dd.py and gmsh_mos2d_create.py. And the temperature boundary I set just for test, which didn’t have actual physical background. The code works fine, but it turns out that the heat source isn’t doing the trick. Could anyone tell me where I am wrong?
The second problem is that the current and electric field both belongs to the edge model. But to calculate the heat generation, it seem to I need transfer them to a node model. Has anyone done something similar ？Please give me some advice.
