I am trying to solve drift-diffusion equatons (Poisson's equation, continuity equations for electrons and holes and Kirchhoff's equation) for a reversely biased diode in a stationary state (no time dependences therefore). Recombination, impact ionization and drift velocity nonlinearity are taken into account. However, I need to implement four boundary conditions for ohmic contacts (two for each) instead of initial conditions. Commercial TCADs do it quite well, but I need Mathematica to solve it. My code is
q = 16*10^-20; t0 =
1 10^-6; A = 10^9; \[Epsilon] = 12; Subscript[\[Epsilon], 0] =
882*10^-16; W = 10^-2; R = 50; Vs = 10^7;
Esp = 23200;
Esn = 8000;
bn = 11*10^5;
bp = 22*10^5;
\[Alpha]ns = 74*10^4;
\[Alpha]ps = 725*10^3;
S = 10^-2; Er = 5/W;
Ur = 999;
\[Tau]n = 10*10^-6; \[Tau]p = 10*10^-6;
ni = 58*10^8;
T = 300;
k = 138*10^-25;
large = 100000;
\[Delta] = 10^-4;
SetSystemOptions["CheckMachineUnderflow" -> False];
Nd[x_] =
10^18*1/(1 + Exp[2*large*x]) - 10^18*1/(1 + Exp[2*large*(W - x)]);
Vn[x_] = Vs*F[x]/(Esn + F[x]); Vp[x_] = Vs*F[x]/(Esp + F[x]);
Difn[x_] = (k*T)/q*Vs*1/(Esn + F[x]);
Difp[x_] = (k*T)/q*Vs*1/(Esp + F[x]); \[Alpha]n[x_] = \[Alpha]ns*
Exp[-bn/(F[x])]; \[Alpha]p[x_] = \[Alpha]ps*Exp[-bp/(F[x])];
equ = {
0 == \[Alpha]n[x]*Vn[x]*n[x] + \[Alpha]p[x]*Vp[x]*p[x] +
D[Difp[x]*D[p[x], x], x] - (
p[x]*n[x] - ni^2)/(\[Tau]p*(n[x] + ni) + \[Tau]n*(p[x] + ni)) -
D[Vp[x]*p[x], x],
D[F[x], x] ==
q/(\[Epsilon]*Subscript[\[Epsilon], 0])*(Nd[x] + p[x] - n[x]), (
A*t0 - Ur)/(q*R*S) ==
Vp[x]*p[x] + Vn[x]*n[x] - Difp[x]*D[p[x], x] + Difn[x]*D[n[x], x],
n[0 - \[Delta]] == (Sqrt[(Nd[0 - \[Delta]])^2/4 + ni^2] +
Nd[0 - \[Delta]]/2),
p[0 - \[Delta]] ==
ni^2/(Sqrt[(Nd[0 - \[Delta]])^2/4 + ni^2] + Nd[0 - \[Delta]]/2),
n[W + \[Delta]] ==
ni^2/(Sqrt[(Nd[W + \[Delta]])^2/4 + ni^2] - Nd[W + \[Delta]]/2),
p[W + \[Delta]] == (Sqrt[(Nd[W + \[Delta]])^2/4 + ni^2] -
Nd[W + \[Delta]]/2)};
soln = NDSolve[equ, {p, F, n}, {x, 0 - \[Delta], W + \[Delta]}];
Which gives an error:
NDSolve::ntdv: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{"EquationSimplification"->"Residual"}.
The same error occurs when I try shooting method.
When I do what Mathematica suggested
soln = NDSolve[equ, {p, F, n}, {x, 0 - \[Delta], W + \[Delta]},
Method -> {"EquationSimplification" -> "Residual"}];
the error becomes:
NDSolve::bvdae: Differential-algebraic equations must be given as initial value problems.
Is there any way to solve my problem?
EDIT #1
When I use
soln = NDSolve[equ, {p, F, n}, {x, 0 - \[Delta], W + \[Delta]},
Method -> {"EquationSimplification" -> "Solve"}];
as xzczd suggested, I get a better result (at least something):
Power::infy: Infinite expression 1/0. encountered.
Infinity::indet: Indeterminate expression E^ComplexInfinity encountered.
Actually, the electric field F(x) must be positive everywhere, so I don't know what's wrong here.
Method -> {"EquationSimplification" -> "Solve"}
orSolveDelayed->False
to forceNDSolve
to find the explicit formula,NDSolve
is capable of doing that in your case. Then you'll see some other warning, perhaps something is wrong with your model, perhaps you need to tackle the initial guess ofShooting
method carefully, I haven't looked into it so I'm not sure, but at least you're in the correct direction now. $\endgroup$