van Roosbroeck system
In both of the following examples, we solve the van Roosbroeck equations, a system of partial differential equations which describe current flow in a bipolar multi layer device:
\[\begin{aligned} + - \nabla \cdot (\varepsilon_s \nabla \psi) &= q \Big( (p(\psi, \varphi_p) - C_p ) - (n(\psi, \varphi_n) - C_n) \Big),\\ + q \partial_t n(\psi, \varphi_n) -\nabla \cdot \mathbf{j}_n &= -qR(n,p), \\ + q \partial_t p(\psi, \varphi_p) + \nabla \cdot \mathbf{j}_p &= -qR(n,p). +\end{aligned}\]
Ohmic contacts will be used as boundary conditions. We will proceed as follows
Step 1: Initialize grid
Step 2: Initialize physical model
Step 3: Solve the problem in equilibrium
Step 4: Solve the problem for an applied bias
Example 1: Stationary 1D problem (region doping)
We consider a three-layer GaAs p-i-n device in one dimension. We will explain the PIN example in greater detail.
Step 1: Initialize grid
We have three layers and two external boundaries. We would like to solve the van Roosbroeck system on a uniform mesh with local grid refinement. We declare subregions and external boundaries.
## region numbers
+regionAcceptor = 1 # p doped region
+regionIntrinsic = 2 # intrinsic region
+regionDonor = 3 # n doped region
+regions = [regionAcceptor, regionIntrinsic, regionDonor]
+numberOfRegions = length(regions)
+
+## boundary region numbers
+# Note that by convention we have 1 for the left boundary and 2 for the right boundary. If
+# adding additional interior boundaries, continue with 3, 4, ...
+bregionAcceptor = 1
+bregionDonor = 2
+bregionJunction1 = 3
+bregionJunction2 = 4
+
+## grid
+refinementfactor = 2^(n-1)
+h_pdoping = 2.0 * μm
+h_intrinsic = 2.0 * μm
+h_ndoping = 2.0 * μm
+h_total = h_pdoping + h_intrinsic + h_ndoping
+w_device = 0.5 * μm # width of device
+z_device = 1.0e-4 * cm # depth of device
+coord = initialize_pin_grid(refinementfactor,
+ h_pdoping,
+ h_intrinsic,
+ h_ndoping)
+
+grid = simplexgrid(coord)
+
+## cellmask! for defining the subregions and assigning region number
+cellmask!(grid, [0.0 * μm], [h_pdoping], regionAcceptor) # p-doped region = 1
+cellmask!(grid, [h_pdoping], [h_pdoping + h_intrinsic], regionIntrinsic) # intrinsic region = 2
+cellmask!(grid, [h_pdoping + h_intrinsic], [h_pdoping + h_intrinsic + h_ndoping], regionDonor) # n-doped region = 3
+
+## bfacemask! for setting different boundary regions. At exterior boundaries they are automatically set by
+## ExtendableGridsjl. Thus, there the following two lines are actually unneccesarry, but are only written for completeness.
+bfacemask!(grid, [0.0], [0.0], bregionAcceptor) # outer left boundary
+bfacemask!(grid, [h_total], [h_total], bregionDonor) # outer right boundary
+bfacemask!(grid, [h_pdoping], [h_pdoping], bregionJunction1) # first inner interface
+bfacemask!(grid, [h_pdoping + h_intrinsic], [h_pdoping + h_intrinsic], bregionJunction2) # second inner interface
Step 2: Initialize physical model
Next, we choose relevant physical models such as the underlying statistics function or the recombination model. Additional options are stated in the comments. Furthermore, we define the charge carrier indices. The index for the electrostatic potential is set automatically to numberOfCarriers + 1
.
# Set indices for the quasi Fermi potentials
+iphin = 1 # electrons
+iphip = 2 # holes
+numberOfCarriers = 2
+
+# Initialize Data instance
+data = Data(grid, numberOfCarriers)
+
+# Solve the stationary problem instead of the transient one
+data.modelType = Stationary
+
+# Choose statistical relation between density and qF potential
+# options: Boltzmann, FermiDiracOneHalfBednarczyk,
+# FermiDiracOneHalfTeSCA FermiDiracMinusOne, Blakemore
+data.F .= Boltzmann
+
+# Enable/Disable recombination processes, the default is stationary SRH recombination.
+data.bulkRecombination = set_bulk_recombination(;iphin = iphin, iphip = iphip,
+ bulk_recomb_Auger = true,
+ bulk_recomb_radiative = true,
+ bulk_recomb_SRH = true)
+
+# choose boundary models
+# exterior boundaries: OhmicContact and SchottkyContact
+# interior boundaries: InterfaceModelNone, InterfaceModelSurfaceReco.
+data.boundaryType[bregionAcceptor] = OhmicContact
+data.boundaryType[bregionDonor] = OhmicContact
+
+# choose flux discretization scheme: ScharfetterGummel ScharfetterGummelGraded,
+# ExcessChemicalPotential, ExcessChemicalPotentialGraded, DiffusionEnhanced, GeneralizedSG
+data.fluxApproximation .= ExcessChemicalPotential
Next, we fill in pre-defined or externally read in parameter values.
# params contains all necessary physical parameters
+params = Params(grid, numberOfCarriers)
+params.temperature = T
+params.UT = (kB * params.temperature) / q
+params.chargeNumbers[iphin] = -1
+params.chargeNumbers[iphip] = 1
+
+for ireg in 1:numberOfRegions # region data
+
+ params.dielectricConstant[ireg] = εr * ε0
+
+ # effective DOS, band-edge energy and mobilities
+ params.densityOfStates[iphin, ireg] = Nc
+ params.densityOfStates[iphip, ireg] = Nv
+ params.bandEdgeEnergy[iphin, ireg] = Ec
+ params.bandEdgeEnergy[iphip, ireg] = Ev
+ params.mobility[iphin, ireg] = mun
+ params.mobility[iphip, ireg] = mup
+
+ # recombination parameters
+ params.recombinationRadiative[ireg] = Radiative
+ params.recombinationSRHLifetime[iphin, ireg] = SRH_LifeTime
+ params.recombinationSRHLifetime[iphip, ireg] = SRH_LifeTime
+ params.recombinationSRHTrapDensity[iphin, ireg] = SRH_TrapDensity
+ params.recombinationSRHTrapDensity[iphip, ireg] = SRH_TrapDensity
+ params.recombinationAuger[iphin, ireg] = Auger
+ params.recombinationAuger[iphip, ireg] = Auger
+
+end
+
+# doping
+params.doping[iphin, regionDonor] = Nd
+params.doping[iphin, regionIntrinsic] = ni
+params.doping[iphip, regionIntrinsic] = 0.0
+params.doping[iphip, regionAcceptor] = Na
+
+# Initialize a ChargeTransport struct
+data.params = params
+ctsys = System(grid, data, unknown_storage=unknown_storage)
Step 3: Solve the problem in equilibrium
Solve the equilibrium. Note that control
refers to the SolverControl parameters given in VoronoiFVM
.
solution = equilibrium_solve!(ctsys, control = control)
+inival = solution
Step 4: Solve the problem for an applied bias
Starting from the equilibrium solution, we increase the applied voltage.
maxBias = voltageAcceptor # bias at acceptor boundary
+biasValues = range(0, stop = maxBias, length = 32)
+
+for Δu in biasValues
+ set_contact!(ctsys, bregionAcceptor, Δu = Δu) # non equilibrium bc
+ solution = solve(ctsys; inival = inival, control = control)
+ inival .= solution
+end
To be consistent with the latest changes of VoronoiFVM, please do not use the solve!() function anymore. Otherwise, you will get deprecation warnings.
Step 5: Postprocessing
By adding the following line to the previous loop
current = get_current_val(ctsys, solution)
we have the possibility to calculate the total current.
Moreover, there are several different plotting routines, see ct_plotting.jl.
Example 2: Stationary 1D problem (nodal doping)
Now, instead of using regionwise doping it is possible to apply a nodal doping. (This is indeed also possible for other physical parameters, see the description of ParamsNodal.) For this, go to previous Step 2, where you build your parameter set and adjust the doping initialization (code snippet is from this example)
paramsnodal = ParamsNodal(grid, numberOfCarriers)
+
+# initialize the space dependent doping
+NDoping = 1.0e17 / cm^3; κ = 500.0
+for icoord = 1:numberOfNodes
+ t1 = tanh( (0.1 - coord[icoord]/μm) *κ )
+ t2 = 1.0 + tanh( (coord[icoord]/μm - 0.2) * κ )
+ paramsnodal.doping[icoord] = NDoping * 0.5 * ( 1.0 + t1 - t2 )
+end
+
+data.paramsnodal = paramsnodal