-
Notifications
You must be signed in to change notification settings - Fork 0
/
StripPlanar_SolvePoisson2D.m
180 lines (162 loc) · 7.93 KB
/
StripPlanar_SolvePoisson2D.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve Poisson equation to compute the potential %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Bulk = Bulk thickness [um]
% PitchX = Pitch along X [um]
% BiasV = Sensor backplane voltage [V] == 0 ? compute weighting field
% epsR = Relative permittivity
% rho = Charge density in the bulk [(Coulomb/um^3)]
function [pdem,Potential,DecomposedGeom,BulkStart,BulkStop,VolumeHeight] =...
StripPlanar_SolvePoisson2D(Bulk,PitchX,BiasV,epsR,rho)
TStart = cputime; % CPU time at start
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Variable initialization %
%%%%%%%%%%%%%%%%%%%%%%%%%%%
NStripLayers = 1; % Number of strip layers
eps0 = 8.85e-18; % Vacuum permittivity [F/um]
VolumeHeight = 3; % Volume height [units of bulk thickness]
MeshMax = 5; % Maximum mesh edge length [um]
MetalThick = 5; % Metalization thickness [um]
MetalWidth = PitchX-20; % Metalization width [um]
BulkStart = 0; % Bulk start coordinate [um]
BulkStop = Bulk+(NStripLayers-1)*MetalThick; % Bulk stop coordinate [um]
NStrips = 13; % Total number of strips
BiasW = 0; % Bias to compute weighting potential
if BiasV == 0
BiasW = 1;
end
%%%%%%%%%%%%%%%%%%%%
% Create PDE model %
%%%%%%%%%%%%%%%%%%%%
fprintf('@@@ I''m solving Poisson equation in 2D to calculate the potential @@@\n');
pdem = createpde('electromagnetic','electrostatic');
%%%%%%%%%%%%%%%%%%%%%%
% Create 2D geometry %
%%%%%%%%%%%%%%%%%%%%%%
gd = zeros(10,NStripLayers*NStrips+2);
bulkStrips = ' (RSV - (';
sf = 'RWV - (';
ns = zeros(3,NStripLayers*NStrips+2);
for l = 1:NStripLayers
% Central strips
indx = 1+(l-1)*NStrips;
gd(:,indx) = [ 3 4 -MetalWidth/2 MetalWidth/2 MetalWidth/2 -MetalWidth/2 ...
BulkStart+l*Bulk/NStripLayers+MetalThick BulkStart+l*Bulk/NStripLayers+MetalThick BulkStart+l*Bulk/NStripLayers BulkStart+l*Bulk/NStripLayers ]';
sf = strcat(sf,sprintf('R%02d+',indx));
ns(:,indx) = sprintf('R%02d',indx);
bulkStrips = strcat(bulkStrips,sprintf('R%02d+',indx));
% Positive strips
for s = 1:(NStrips-1)/2
indx = s+1+(l-1)*NStrips;
gd(:,indx) = [ 3 4 s*PitchX-MetalWidth/2 s*PitchX+MetalWidth/2 s*PitchX+MetalWidth/2 s*PitchX-MetalWidth/2 ...
BulkStart+l*Bulk/NStripLayers+MetalThick BulkStart+l*Bulk/NStripLayers+MetalThick BulkStart+l*Bulk/NStripLayers BulkStart+l*Bulk/NStripLayers ]';
sf = strcat(sf,sprintf('R%02d+',indx));
ns(:,indx) = sprintf('R%02d',indx);
bulkStrips = strcat(bulkStrips,sprintf('R%02d+',indx));
end
% Negative strips
for s = 1:(NStrips-1)/2
indx = s+1+(NStrips-1)/2+(l-1)*NStrips;
gd(:,s+1+(NStrips-1)/2+(l-1)*NStrips) = [ 3 4 -(s*PitchX-MetalWidth/2) -(s*PitchX+MetalWidth/2) -(s*PitchX+MetalWidth/2) -(s*PitchX-MetalWidth/2) ...
BulkStart+l*Bulk/NStripLayers+MetalThick BulkStart+l*Bulk/NStripLayers+MetalThick BulkStart+l*Bulk/NStripLayers BulkStart+l*Bulk/NStripLayers ]';
sf = strcat(sf,sprintf('R%02d+',indx));
ns(:,indx) = sprintf('R%02d',indx);
bulkStrips = strcat(bulkStrips,sprintf('R%02d+',indx));
end
end
% Sensor volume RSV
gd(:,NStrips*NStripLayers+1) = [ 3 4 -((NStrips-1)/2*PitchX+PitchX/2) ((NStrips-1)/2*PitchX+PitchX/2) ((NStrips-1)/2*PitchX+PitchX/2) -((NStrips-1)/2*PitchX+PitchX/2) ...
BulkStop BulkStop BulkStart BulkStart ]';
% Whole volume RWV
gd(:,NStrips*NStripLayers+2) = [ 3 4 -((NStrips-1)/2*PitchX+PitchX/2) ((NStrips-1)/2*PitchX+PitchX/2) ((NStrips-1)/2*PitchX+PitchX/2) -((NStrips-1)/2*PitchX+PitchX/2) ...
Bulk*VolumeHeight Bulk*VolumeHeight BulkStart BulkStart ]';
bulkStrips = strcat(bulkStrips(1:end-1),'))');
sf = strcat(sf(1:end-1),strcat(') +',bulkStrips));
ns(:,NStrips*NStripLayers+1) = 'RSV';
ns(:,NStrips*NStripLayers+2) = 'RWV';
ns = char(ns);
DecomposedGeom = decsg(gd,sf,ns);
geometryFromEdges(pdem,DecomposedGeom);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Apply boundary conditions (only on conductors) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if NStripLayers == 1
% Central strip
electromagneticBC(pdem,'Voltage',BiasW,'edge',1:2);
% Positive strips
electromagneticBC(pdem,'Voltage',0,'edge',3:14);
% Negative strips
electromagneticBC(pdem,'Voltage',0,'edge',15:26);
% Top all strips
electromagneticBC(pdem,'Voltage',0,'edge',28:33);
electromagneticBC(pdem,'Voltage',BiasW,'edge',34);
electromagneticBC(pdem,'Voltage',0,'edge',35:40);
% Bottom all strips
electromagneticBC(pdem,'Voltage',0,'edge',linspace(42,52,6));
electromagneticBC(pdem,'Voltage',BiasW,'edge',54);
electromagneticBC(pdem,'Voltage',0,'edge',linspace(56,66,6));
% Top edge
electromagneticBC(pdem,'Voltage',0,'edge',27);
% Right edge sensor
electromagneticBC(pdem,'SurfaceCurrentDensity',0,'edge',68);
% Right edge air
electromagneticBC(pdem,'SurfaceCurrentDensity',0,'edge',69);
% Bottom edge
electromagneticBC(pdem,'Voltage',BiasV,'edge',70);
% Left edge sensor
electromagneticBC(pdem,'SurfaceCurrentDensity',0,'edge',71);
% Left edge air
electromagneticBC(pdem,'SurfaceCurrentDensity',0,'edge',72);
else
% All internal strips
top = [149 166 214];
bottom = [179 196 2];
for i = 0:NStripLayers-2
electromagneticBC(pdem,'Voltage',BiasV/NStripLayers*(NStripLayers-i-1),'edge',top(i+1):top(i+1)+NStrips-1);
electromagneticBC(pdem,'Voltage',BiasV/NStripLayers*(NStripLayers-i-1),'edge',bottom(i+1):bottom(i+1)+NStrips-1);
electromagneticBC(pdem,'Voltage',BiasV/NStripLayers*(NStripLayers-i-1),'edge',linspace(61+i,109+i,(NStrips-1)/2+1));
electromagneticBC(pdem,'Voltage',BiasV/NStripLayers*(NStripLayers-i-1),'edge',linspace(57+i,105+i,(NStrips-1)/2+1));
electromagneticBC(pdem,'Voltage',BiasV/NStripLayers*(NStripLayers-i-1),'edge',linspace(113+i,145+i,(NStrips-1)/2-1));
electromagneticBC(pdem,'Voltage',BiasV/NStripLayers*(NStripLayers-i-1),'edge',linspace(117+i,141+i,(NStrips-1)/2-2));
electromagneticBC(pdem,'Voltage',BiasV/NStripLayers*(NStripLayers-i-1),'edge',192+i);
electromagneticBC(pdem,'Voltage',BiasV/NStripLayers*(NStripLayers-i-1),'edge',210+i);
electromagneticBC(pdem,'Voltage',BiasV/NStripLayers*(NStripLayers-i-1),'edge',162+i);
end
% Top strips
electromagneticBC(pdem,'Voltage',0,'edge',15:15+NStrips-1);
electromagneticBC(pdem,'Voltage',0,'edge',linspace(29,29+(NStrips-1)*2,NStrips));
electromagneticBC(pdem,'Voltage',0,'edge',[213,195,165,148,144,140,136,132,...
128,124,120,116,72,68,80,76,88,84,96,92,104,100,112,108]);
% Central top strip
electromagneticBC(pdem,'Voltage',BiasW,'edge',21);
electromagneticBC(pdem,'Voltage',BiasW,'edge',41);
electromagneticBC(pdem,'Voltage',BiasW,'edge',60);
electromagneticBC(pdem,'Voltage',BiasW,'edge',64);
% Top edge
electromagneticBC(pdem,'Voltage',0,'edge',209);
% Right edge sensor
electromagneticBC(pdem,'SurfaceCurrentDensity',0,'edge',55);
% Right edge air
electromagneticBC(pdem,'SurfaceCurrentDensity',0,'edge',56);
% Left edge sensor
electromagneticBC(pdem,'SurfaceCurrentDensity',0,'edge',227);
% Left edge air
electromagneticBC(pdem,'SurfaceCurrentDensity',0,'edge',228);
% Bottom edge
electromagneticBC(pdem,'Voltage',BiasV,'edge',1);
end
%%%%%%%%%%%%%%%%%
% Generate mesh %
%%%%%%%%%%%%%%%%%
msh = generateMesh(pdem,'Hmax',MeshMax,'GeometricOrder','quadratic');
%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve Poisson equation %
%%%%%%%%%%%%%%%%%%%%%%%%%%
pdem.VacuumPermittivity = eps0;
electromagneticSource(pdem,'face',1,'ChargeDensity',0); % Air
electromagneticProperties(pdem,'RelativePermittivity',1,'face',1); % Air
electromagneticSource(pdem,'face',2,'ChargeDensity',rho); % Sensor
electromagneticProperties(pdem,'RelativePermittivity',epsR,'face',2); % Sensor
Potential = solve(pdem);
fprintf('CPU time --> %.2f [min]\n\n',(cputime-TStart)/60);
end