Bug fix for interface recombination (#21)

* Bug fix for interface recombination

The SRH interface recombination rates were incorrectly coded with n_i^2/k instead of n_i^2. This has been fixed and the rates have been rewritten in terms of the perovskite concentrations for clarity, meaning that the TL-related recombination parameters are now "effective" parameters rather than those that appear in the updated interface recombination rates Rl(n,p) and Rr(n,p).

* Fix bimolecular rates

Code/Common/nondimensionalise.m

@@ -116,11 +116,11 @@
     [gammaH, torH, torH3] = deal(0); % no perovskite/HTL interface recombination
 end
 % Total ETL/perovskite interface recombination rate
-Rl = @(nE,p) brateE*(nE.*p-ni2/kE) ... % bimolecular
-           + SRH(nE,p,gammaE,ni2/kE,torE,torE3); % SRH
+Rl = @(n,p) brateE*(n.*p-ni2) ... % bimolecular
+          + SRH(n,p,gammaE,ni2,torE,torE3); % SRH
 % Total perovskite/HTL interface recombination rate
-Rr = @(n,pH) brateH*(n.*pH-ni2/kH) ... % bimolecular
-           + SRH(pH,n,gammaH,ni2/kH,torH,torH3); % SRH
+Rr = @(n,p) brateH*(n.*p-ni2) ... % bimolecular
+          + SRH(p,n,gammaH,ni2,torH,torH3); % SRH
 
 % Spatial grid parameters (consistent with the choice of N above)
 X = 0.2; % percentage of grid points within an ionic Debye length of the interface

Code/Solver/AnJac.m

@@ -113,16 +113,16 @@
 JJJ(2*N+3:3*N+3,2*N+3:3*N+3) = JJJ(2*N+3:3*N+3,2*N+3:3*N+3) ...
     -1/2*gallery('tridiag',N+1, ...
     dx.*(RRn-RR)/del, ...
-    [(dx(1).*(RRn(1)-RR(1))+2*(Rl(nE(end)+del,p(1))-Rl(nE(end),p(1))))/del; ...
+    [(dx(1).*(RRn(1)-RR(1))+2*(Rl(n(1)+del,p(1))-Rl(n(1),p(1))))/del; ...
         (dx(2:N).*(RRn(2:N)-RR(2:N))+dx(1:N-1).*(RRn(1:N-1)-RR(1:N-1)))/del; ...
-        (dx(N)*(RRn(N)-RR(N))+2*(Rr(n(N+1)+del,pH(1))-Rr(n(N+1),pH(1))))/del], ...
+        (dx(N)*(RRn(N)-RR(N))+2*(Rr(n(N+1)+del,p(N+1))-Rr(n(N+1),p(N+1))))/del], ...
     dx.*(RRn-RR)/del);
 % n equation depends on p via GR, Rl and Rr
 JJJ(2*N+3:3*N+3,3*N+4:4*N+4) = -1/2*gallery('tridiag',N+1, ...
     dx.*(RRp-RR)/del, ...
-    [(dx(1).*(RRp(1)-RR(1))+2*(Rl(nE(end),p(1)+del)-Rl(nE(end),p(1))))/del; ...
+    [(dx(1).*(RRp(1)-RR(1))+2*(Rl(n(1),p(1)+del)-Rl(n(1),p(1))))/del; ...
         (dx(2:N).*(RRp(2:N)-RR(2:N))+dx(1:N-1).*(RRp(1:N-1)-RR(1:N-1)))/del; ...
-        (dx(N)*(RRp(N)-RR(N))+2*(Rr(n(N+1),pH(1)+del)-Rr(n(N+1),pH(1))))/del], ...
+        (dx(N)*(RRp(N)-RR(N))+2*(Rr(n(N+1),p(N+1)+del)-Rr(n(N+1),p(N+1))))/del], ...
     dx.*(RRp-RR)/del);
 % n equation depends on phiE
 JJJ(2*N+3,4*N+NE+4) = -KE*AvnE(end)./dxE(end);
@@ -147,9 +147,9 @@
 % p equation depends on n via GR, Rl and Rr
 JJJ(3*N+4:4*N+4,2*N+3:3*N+3) = -1/2*gallery('tridiag',N+1, ...
     dx.*(RRn-RR)/del, ...
-    [(dx(1)*(RRn(1)-RR(1))+2*(Rl(nE(end)+del,p(1))-Rl(nE(end),p(1))))/del; ...
+    [(dx(1)*(RRn(1)-RR(1))+2*(Rl(n(1)+del,p(1))-Rl(n(1),p(1))))/del; ...
         (dx(2:N).*(RRn(2:N)-RR(2:N))+dx(1:N-1).*(RRn(1:N-1)-RR(1:N-1)))/del; ...
-        (dx(N).*(RRn(N)-RR(N))+2*(Rr(n(N+1)+del,pH(1))-Rr(n(N+1),pH(1))))/del], ...
+        (dx(N).*(RRn(N)-RR(N))+2*(Rr(n(N+1)+del,p(N+1))-Rr(n(N+1),p(N+1))))/del], ...
     dx.*(RRn-RR)/del);
 % p equation depends on p
 JJJ(3*N+4:4*N+4,3*N+4:4*N+4) = Kp*gallery('tridiag',N+1, ...
@@ -161,9 +161,9 @@
 JJJ(3*N+4:4*N+4,3*N+4:4*N+4) = JJJ(3*N+4:4*N+4,3*N+4:4*N+4) ...
     -1/2*gallery('tridiag',N+1, ...
     dx.*(RRp-RR)/del, ...
-    [(dx(1)*(RRp(1)-RR(1))+2*(Rl(nE(end),p(1)+del)-Rl(nE(end),p(1))))/del; ...
+    [(dx(1)*(RRp(1)-RR(1))+2*(Rl(n(1),p(1)+del)-Rl(n(1),p(1))))/del; ...
         (dx(2:N).*(RRp(2:N)-RR(2:N))+dx(1:N-1).*(RRp(1:N-1)-RR(1:N-1)))/del; ...
-        (dx(N).*(RRp(N)-RR(N))+2*(Rr(n(N+1),pH(1)+del)-Rr(n(N+1),pH(1))))/del], ...
+        (dx(N).*(RRp(N)-RR(N))+2*(Rr(n(N+1),p(N+1)+del)-Rr(n(N+1),p(N+1))))/del], ...
     dx.*(RRp-RR)/del);
 % p equation does not depend on phiE
 % p equation does not depend on nE

Code/Solver/RHS.m

@@ -61,14 +61,14 @@
     -rH*dxH(1)*(1/2-pH(1,:)/3-pH(2,:)/6)/lamH2; % continuity
 
 % n equation
-dudt(2*N+3,:) = kE*fn(1,:)-fnE(end,:)+kE*(dx(1).*GR(1,:)/2-Rl(nE(end,:),p(1,:))); % continuity
+dudt(2*N+3,:) = kE*fn(1,:)-fnE(end,:)+kE*(dx(1).*GR(1,:)/2-Rl(n(1,:),p(1,:))); % continuity
 dudt(2*N+4:3*N+2,:) = fn(2:N,:)-fn(1:N-1,:)+(dx(2:N).*GR(2:N,:)+dx(1:N-1).*GR(1:N-1,:))/2;
-dudt(3*N+3,:) = -fn(N,:)+dx(N)*GR(end,:)/2-Rr(n(N+1,:),pH(1,:));
+dudt(3*N+3,:) = -fn(N,:)+dx(N)*GR(end,:)/2-Rr(n(N+1,:),p(N+1,:));
 
 % p equation
-dudt(3*N+4,:) = fp(1,:)+dx(1)*GR(1,:)/2-Rl(nE(end,:),p(1,:));
+dudt(3*N+4,:) = fp(1,:)+dx(1)*GR(1,:)/2-Rl(n(1,:),p(1,:));
 dudt(3*N+5:4*N+3,:) = fp(2:N,:)-fp(1:N-1,:)+(dx(2:N).*GR(2:N,:)+dx(1:N-1).*GR(1:N-1,:))/2;
-dudt(4*N+4,:) = fpH(1,:)-kH*fp(end,:)+kH*(dx(end)*GR(end,:)/2-Rr(n(N+1,:),pH(1,:))); % continuity
+dudt(4*N+4,:) = fpH(1,:)-kH*fp(end,:)+kH*(dx(end)*GR(end,:)/2-Rr(n(N+1,:),p(N+1,:))); % continuity
 
 % phiE equation
 dudt(4*N+5,:) = phiE(1,:)-psi(t);

parameters_template.m

@@ -90,12 +90,12 @@
 Augp  = 0;       % hole-dominated Auger recombination rate (m6s-1)
 
 % Interface recombination (max. velocity ~ 1e5)
-betaE = 0;       % ETL/perovskite bimolecular recombination rate (m3s-1)
-betaH = 0;       % perovskite/HTL bimolecular recombination rate (m3s-1)
-vnE   = 1e5;     % electron recombination velocity for SRH (ms-1)
+betaE = 0;       % effective ETL/perovskite bimolecular rate (m3s-1)
+betaH = 0;       % effective perovskite/HTL bimolecular rate (m3s-1)
+vnE   = 1e5;     % effective electron recombination velocity for SRH (ms-1)
 vpE   = 10;      % hole recombination velocity for SRH (ms-1)
 vnH   = 0.1;     % electron recombination velocity for SRH (ms-1)
-vpH   = 1e5;     % hole recombination velocity for SRH (ms-1)
+vpH   = 1e5;     % effective hole recombination velocity for SRH (ms-1)
 
 % Parasitic resistances (optional)
 % Rs    = 0;       % external series resistance (Ohm)