How to solve continuity equations together with Poisson equation?

13 次查看（过去 30 天）
Sven B.
Sven B. 2016-12-22
回答： Sven B. 2016-12-30
As I'm working a lot with semiconductor phyics, I wonder if there is a way to solve the common continuity equations together with the Poisson equation. Perhaps this is a known issue? I already tried it with pdepe, but continue to receive following error message:
Error using pdepe (line 293)
Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial
derivative.
Error in ContinuityEquations (line 27)
sol = pdepe(0,pdefun,icfun,bcfun,x,t);
I tried the following two versions of pdefun (all physical constants are set to 1). First version doesn't work because of the zero in the 3rd component of the flux-term I guess. Does anyone know what is wrong with the 2nd version?
function [c,f,s] = pdefun(x,t,u,DuDx)
c = [1;1;0];
f = [DuDx(1); DuDx(2); 0];
s = [u(1)*DuDx(3)+u(3)*DuDx(1); ...
-u(2)*DuDx(3)-u(3)*DuDx(2); ...
-DuDx(3)+u(2)-u(1)];
end
function [c,f,s] = pdefun(x,t,u,DuDx)
c = [1;1;0];
f = [DuDx(1)-u(1)*DuDx(3); ...
DuDx(2)+u(2)*DuDx(3); ...
DuDx(3)];
s = [-DuDx(3)*DuDx(1)+DuDx(1)*DuDx(3); ...
DuDx(3)*DuDx(2)-DuDx(2)*DuDx(3); ...
u(2)-u(1)];
end

请先登录，再回答此问题。

回答（2 个）

Sharmila Raghu
Sharmila Raghu 2016-12-29
The above error might occur if the boundary conditions are ill-posed. Please verify the boundary conditions to see if they are ill-posed. The boundary conditions specified as "p" and "q", follow this relationship:
p + q*f = 0
If "pr" and "qr" (the parameters for the right boundary) are both zero, this becomes 0+0=0, which is an ill-posed problem. To resolve this, please try setting "qr" to anything besides zero. This is equivalent to "f=0", which was probably the intended result.
Sven B.
Sven B. 2016-12-30
Thanks for you answer! Indeed I had an ill-posed boundary condition, but unfortunately I still suffer from the same problem.
Find attached my complete code so far and the whole task in a better readable mathematical notation.
% Mesh
x = (-2:0.01:2)*1e-4; % cm
t = logspace(-16,-10,1000); % s
% Parameters
N_D = 1e17; % cm-3
N_A = 1e10; % cm-3
mu_n = 1000; % cm²/Vs
mu_p = 350; % cm²/Vs
epsilon_0 = 8.854187817e-14; % As/Vcm
epsilon_r = 11.9;
e = 1.602177e-19; % As
k_B = 8.617330e-5; % eV/K
% Starting conditions
gauss = @(x,A,mu,sigma) A*exp(-(x-mu).^2/(2*sigma^2));
n0 = @(x) N_D + gauss(x,1e18,0,1e-5); % cm-3
p0 = @(x) zeros(1,numel(x))+N_A; % cm-3
dEdx0 = @(x) e/(epsilon_0*epsilon_r)*(p0(x)-n0(x)+N_D-N_A); % V/cm2
E0 = @(x) int(x,dEdx0)-0.5*integral(dEdx0,x(1),x(end)); % V/cm (see below for 'int')
% Solution
pdefun = @(x,t,u,DuDx) pdefunH(x,t,u,DuDx,N_D,N_A,mu_n,mu_p,epsilon_0,epsilon_r,e,k_B);
icfun = @(x) icfunH(x,n0,p0,E0);
bcfun = @(xl,ul,xr,ur,t) bcfunH(xl,ul,xr,ur,t,E0,x);
sol = pdepe(0,pdefun,icfun,bcfun,x,t);
function [c,f,s] = pdefunH(x,t,u,DuDx,N_D,N_A,mu_n,mu_p,epsilon_0,epsilon_r,e,k_B)
c = [1;1;0];
D_n = k_B*300/e*mu_n; %%cm2/s
D_p = k_B*300/e*mu_p; %%cm2/s
f = [D_n*DuDx(1); ...
D_p*DuDx(2); ...
eps*DuDx(3)];
s = [mu_n*u(1)*DuDx(3)+mu_n*u(3)*DuDx(1); ...
-mu_p*u(2)*DuDx(3)-mu_p*u(3)*DuDx(2); ...
DuDx(3)-e/(epsilon_0*epsilon_r)*(u(2)-u(1)+N_D-N_A)];
end
function u = icfunH(x,n0,p0,E0)
u = [n0(x);p0(x);E0(x)];
end
function [pl,ql,pr,qr] = bcfunH(xl,ul,xr,ur,t,E0,x)
pl = [0;0;0];
ql = [1;1;1];
pr = [0;0;0];
qr = [1;1;1];
end
function F = int( x, f )
F = zeros(1,numel(x));
for n = 1:numel(x)
F(n) = integral(f,x(1),x(n));
end
end

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by