# Simulating the Response of Charged Particles to Apertures with Inhomogeneous Electric Fields Passing through them

I have produced a simulation which allows me to observe a static electric field produced by a charged plate as it passes through some complex apertures in a different, oppositely charged metal plate.

I am curious what one might do to insert N charged "particles" which are moving randomly (with radius r) into a particular region of my simulation (the negative-z space of my "airbox") and observe their motion over some time (final result is fine). Another caveat is that the particles must stop at the aperture opening at z=0 (if they are pulled toward it). Furthermore, the particles will not be allowed to occupy the same volume simultaneously (must stack on top of one another if there is no room at aperture boundary, and can collide).

I assume the kinematic equations would follow some form like r''[t] + a*r'[t] - (q/m)*E[r] == 0, where "a" is a constant and E[r] is the electric field strength in the region the particle occupies... initial conditions would give the particles slightly random r'[0] and r[0] and slightly different charge/mass ratios.

Here is the .obj file:

https://www.dropbox.com/s/dds8rm3odg2m7gu/largeAp.obj?dl=0

and here is the code.

I will continue to experiment and figure this out as well, and will post a solution if I ever find one. Thank you!

Clear["Global*"];
Needs["NDSolveFEM"]

R = Import["C:\\Users\\zaust\\OneDrive\\Desktop\\largeAp.obj"];
RegionDimension[R];
M = BoundaryMeshRegion[MeshCoordinates[R], MeshCells[R, 2]];
RegionDimension[M];
Volume[M];

r = RegionDifference[
RegionDifference[
RegionDifference[Cuboid[{0, 0, -0.5}, {2, 2, 0.5}], M],
Cuboid[{0, 0, 0.4}, {2, 2, 0.5}]],
Cuboid[{0, 0, -0.5}, {2, 2, -0.4}]];
mesh = ToElementMesh[r]["Wireframe"];

pol = -1;

V0 = 4000;
sol = NDSolveValue[{Laplacian[V[x, y, z], {x, y, z}] == 0,
DirichletCondition[
V[x, y, z] == -pol*
V0/2, (0.4 <= z <= 0.5) && (0 <= y <= 2) && (0 <= x <= 2)],
DirichletCondition[
V[x, y, z] ==
pol*V0/2, (0.0071 <= z <= 0.0072) && (0 <= y <= 2) && (0 <= x <=
2)], DirichletCondition[
V[x, y, z] ==
0, (0 <= z <= 0.0070) && (0 <= y <= 2) && (0 <= x <= 2)],
DirichletCondition[
V[x, y, z] ==
0, (-0.5 <= z <= -0.4) && (0 <= y <= 2) && (0 <= x <= 2)]},
V, {x, y, z} \[Element] r]

plotXZ = ContourPlot[sol[x, 0.75, z], {x, 0, 2}, {z, -0.4, 0.1},
ContourShading -> Automatic, ColorFunction -> "Rainbow",
Contours -> 100]
plotYZ = ContourPlot[sol[0.75, y, z], {y, 0, 2}, {z, -0.4, 0.1},
ContourShading -> Automatic, ColorFunction -> "Rainbow",
Contours -> 100]
plotXY = ContourPlot[sol[x, y, -0.05], {x, 0.5, 1.5}, {y, 0.5, 1.5},
ColorFunction -> "Rainbow", Contours -> 50, PlotLegends -> Automatic]

SliceContourPlot3D[sol[x, y, z],
"CenterPlanes", {x, 0.5, 1}, {y, 0.5, 1}, {z, 0, -0.1},
ColorFunction -> "Rainbow", Contours -> 200];
DensityPlot3D[
sol[x, y, z], {x, 0.5, 1}, {y, 0.5, 1}, {z, -0.4, -0.00001},
PlotTheme -> "Detailed", ColorFunction -> "Rainbow",
PerformanceGoal -> "Quality"];

d = ContourPlot3D[
sol[x, y, z], {x, 0.5, 1}, {y, 0.5, 1}, {z, -0.25, -0.00001},
PlotTheme -> "Detailed", ColorFunction -> "Rainbow",
PerformanceGoal -> "Quality", Contours -> 1, Mesh -> None];

electricField[x_, y_, z_] = -Grad[sol[x, y, z], {x, y, z}];
v = Show[VectorPlot3D[
electricField[x, y, z], {x, 0.5, 1}, {y, 0.5,
1}, {z, -0.25, -0.00001}, PlotTheme -> "Detailed",
ColorFunction -> "Rainbow", PerformanceGoal -> "Quality",
VectorScale -> 0.05, VectorPoints -> 7], M]
$$$$
New contributor
zak720 is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.