# How to create an animation of charged particles' motion in an electric field

I am trying to create an animation of 3 particles' motion in an electric field. The electric field is created by a needle attached to a top plate both of which are at -3800V with a grounded plate close below the needle. The needle is positioned above a hole in the plate which will allow particles to travel to the area below the grounded plate. I have created a 2-D axisymmetric simulation of the electric field without issue. The code for which is below.

ClearAll["Global*"]
Needs["NDSolveFEM"]
q = -1.60217733*10^-19*10;(*particle charge*)
voltage = -3800;(*needle and top plate voltage*)
r1 = 0.0065; (*hole radius*)
r2 = 0.0365; (*domain radius*)
r3 = 0.00015; (*needle radius*)
z1 = 0.07; (*height of domain*)
z2 = 0.065; (*height of tip of needle*)
l = 0.005; (*distance between tip of needle and hole in middle plate*)
z3 = z2 - l; (*height of middle plate top surface*)
z4 = z3 - 0.0016; (*height of middle plate bottom surface*)
reg1 = ImplicitRegion[True, {{r, r3, r2}, {z, z2, z1}}]; (*region to the right of needle, above middle plate top surface*)
reg2 = ImplicitRegion[True, {{r, 0, r2}, {z, z3, z2}}]; (*region between tip of needle and middle plate top surface*)
reg3 = ImplicitRegion[True, {{r, 0, r1}, {z, z4, z3}}]; (*region to the left of middle plate*)
reg4 = ImplicitRegion[True, {{r, 0, r2}, {z, 0, z4}}]; (*region below middle plate bottom surface*)
region = RegionUnion[reg1, reg2, reg3, reg4]; (*merge all the previously created regions, space left out represents needle and middle plate*)
meshRefine[vertices_, area_] := area > 0.0000001;
mesh = ToElementMesh[DiscretizeRegion[region], MeshRefinementFunction -> meshRefine];
bc1 = {DirichletCondition[phi[r, z] == voltage, (z == z2 && 0 <= r <= r3)],
DirichletCondition[phi[r, z] == voltage, (z == z1 && r3 <= r <= r2)],
DirichletCondition[phi[r, z] == voltage, (r == r3 && z2 <= z <= z1)],
DirichletCondition[phi[r, z] == 0, (z == z3 && r1 <= r <= r2)],
DirichletCondition[phi[r, z] == 0, (z == z4 && r1 <= r <= r2)],
DirichletCondition[phi[r, z] == 0, (r == r1 && z4 <=z <=z3)]};  (*boundary conditions*)
sol = NDSolveValue[{1/r*D[r*D[phi[r, z], r], r] + D[phi[r, z], z, z] ==0, bc1}, phi, {r, z} \[Element] mesh];
electricField[r_,z_] := -Grad[sol[r, z], {r,z}];
eforce[r_,z_]:= q*electricField[r,z];


I am having trouble using the steady electric field to find the position of each particle in time. (Eventually, I would like to add the force of gravity as well once I have the simulation working with just the electric field.) The particles originate at the tip of the needle and initial velocity is zero. When a particle reaches a boundary I would like it to stop and stay there. Below is the code I have developed so far.

d = 5*10^-5; (*particle diameter in m*)
mass = 6.52*10^-8; (*particle mass in kg/m^3*)
numbodies = 3;
vel0 = Table[{0, 0}, numbodies];
pos0 = Riffle[Range[0, r3, r3/numbodies], z2]~Partition~2;
force[j_, t_] := eforce[p[[j, 1]][t], p[[j, 2]][t]];
odesys = Table[{p[j]''[t] == 1/mass*force[j][t], p[j] == pos0[[j]], p[j]' == vel0[[j]],
WhenEvent[p[[j, 2]][t] == 0, Norm[p[j]'[t]] -> 0], WhenEvent[p[[j, 1]][t] == 0, Norm[p[j]'[t]] -> 0],
WhenEvent[p[[j, 1]][t] == r2, Norm[p[j]'[t]] -> 0], WhenEvent[p[[j, 2]][t] == z1, Norm[p[j]'[t]] -> 0],
WhenEvent[p[[j, 2]][t] == z3 && r1 <= p[[j, 1]][t] <= r2, Norm[p[j]'[t]] -> 0],
WhenEvent[p[[j, 2]][t] == z4 && r1 <= p[[j, 1]][t] <= r2, Norm[p[j]'[t]] -> 0],
WhenEvent[p[[j, 1]][t] == r1 && z4 <= p[[j, 2]][t] <= z3, Norm[p[j]'[t]] -> 0]}, {j, numbodies}];
depvars = Flatten[Table[{p[j]}, {j, numbodies}]];
tfin = 1;
sol2 = NDSolve[odesys, depvars, {t, 0, tfin}][];
pos = Array[p, {numbodies}] /. sol2;
Animate[Show[ParametricPlot[#[t] & /@ sol2, {t, Max[0, tfin], tfin}, ImageSize -> 400, Frame -> True, PlotRange -> {0, 0.07}], Graphics[MapIndexed[{Hue[.35], Disk[#1[tfin], 0.5]} &, pos]]], {tfin, 0.1, tfin, 1}]


This code produces a couple of error messages; the first being in NDSolve saying there are more dependent variables than equations so the system is underdetermined, and the second being in the ReplaceAll in the line below the NDSolve saying what I am trying to replace is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. The code was pieced together from examples I was able to find of people doing something similar to what I am looking to do, but I am very much a beginner with Mathematica, so I don't always understand the functions and syntax used by the examples I pulled from or why I am getting errors when I use the same. If anyone has any suggestions to help me achieve what I am looking to do, or documentation that could help me, I would greatly appreciate it!

• The ReplaceAll error message is likely linked to the NDSolve one; Since NDSolve does not find a solution, then sol2 contains the unevaluated call, rather than a list of rules, so ReplaceAll complains as well. – MarcoB Dec 20 '19 at 21:31

## 1 Answer

Here we can use the coordinates {x[j][t],y[j][t]} to describe the particle trajectories. The purpose of the large number WhenEvent is not entirely clear. If this is a condition for the absorption of particles at the boundaries of the region, then this is not true. Here is an example of working code.

ClearAll["Global*"]
Needs["NDSolveFEM"]
q = -1.60217733*10^-19*10;(*particle charge*)voltage = -3800;(*needle \
and top plate voltage*)r1 = 0.0065;(*hole radius*)r2 = \
0.0365;(*domain radius*)r3 = 0.00015;(*needle radius*)z1 = \
0.07;(*height of domain*)z2 = 0.065;(*height of tip of needle*)l = \
0.005;(*distance between tip of needle and hole in middle plate*)z3 =
z2 - l;(*height of middle plate top surface*)z4 =
z3 - 0.0016;(*height of middle plate bottom surface*)reg1 =
ImplicitRegion[
True, {{r, r3, r2}, {z, z2,
z1}}];(*region to the right of needle,above middle plate top \
surface*)reg2 =
ImplicitRegion[
True, {{r, 0, r2}, {z, z3,
z2}}];(*region between tip of needle and middle plate top \
surface*)reg3 =
ImplicitRegion[
True, {{r, 0, r1}, {z, z4,
z3}}];(*region to the left of middle plate*)reg4 =
ImplicitRegion[
True, {{r, 0, r2}, {z, 0,
z4}}];(*region below middle plate bottom surface*)region =
RegionUnion[reg1, reg2, reg3,
reg4];(*merge all the previously created regions,space left out \
represents needle and middle plate*)
meshRefine[vertices_, area_] := area > 0.0000001;
mesh = ToElementMesh[DiscretizeRegion[region],
MeshRefinementFunction -> meshRefine]; mesh["Wireframe"]
bc1 = {DirichletCondition[
phi[r, z] == voltage, (z == z2 && 0 <= r <= r3)],
DirichletCondition[
phi[r, z] == voltage, (z == z1 && r3 <= r <= r2)],
DirichletCondition[
phi[r, z] == voltage, (r == r3 && z2 <= z <= z1)],
DirichletCondition[phi[r, z] == 0, (z == z3 && r1 <= r <= r2)],
DirichletCondition[phi[r, z] == 0, (z == z4 && r1 <= r <= r2)],
DirichletCondition[
phi[r, z] ==
0, (r == r1 && z4 <= z <= z3)]};(*boundary conditions*)sol =
NDSolveValue[{1/r*D[r*D[phi[r, z], r], r] + D[phi[r, z], z, z] == 0,
bc1}, phi, {r, z} \[Element] mesh];

electricField = -Evaluate[Grad[sol[r, z], {r, z}]];
eforce = q*electricField;

StreamPlot[eforce, {r, z} \[Element] mesh]

DensityPlot[sol[r, z], {r, z} \[Element] mesh,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotRange -> All] Particle trajectories and animation

d = 5*10^-5;(*particle diameter in m*)mass =
6.52*10^-8;(*particle mass in kg/m^3*)numbodies = 3;
vel0 = Table[{0, 0}, numbodies];
pos0 = Riffle[Range[0, r3, r3/numbodies], z2]~Partition~2;
eqs = Table[{x[j]''[t] == 1/mass*eforce[] /. {r -> x[j][t],
z -> y[j][t]},
y[j]''[t] == 1/mass*eforce[] /. {r -> x[j][t], z -> y[j][t]},
x[j] == pos0[[j, 1]], y[j] == pos0[[j, 2]],
x[j]' == vel0[[j, 1]], y[j]' == vel0[[j, 2]]}, {j,
numbodies}];
vars = Flatten[Table[{x[j], y[j]}, {j, numbodies}]];

event = Table[{WhenEvent[
x[j][t] == 0, {x[j]'[t] -> 0, y[j]'[t] -> 0}],
WhenEvent[x[j][t] == 0, {x[j]'[t] -> 0, y[j]'[t] -> 0}],
WhenEvent[x[j][t] == r2, {x[j]'[t] -> 0, y[j]'[t] -> 0}],
WhenEvent[y[j][t] == z1, {x[j]'[t] -> 0, y[j]'[t] -> 0}],
WhenEvent[
y[j][t] == z3 && r1 <= x[j][t] <= r2, {x[j]'[t] -> 0,
y[j]'[t] -> 0}],
WhenEvent[
y[j][t] == z4 && r1 <= x[j][t] <= r2, {x[j]'[t] -> 0,
y[j]'[t] -> 0}],
WhenEvent[
x[j][t] == r1 && z4 <= y[j][t] <= z3, {x[j]'[t] -> 0,
y[j]'[t] -> 0}]} /. j -> i, {i, numbodies}];

tfin = 150; sol1 = NDSolve[{eqs, event}, vars, {t, 0, tfin}][]

dp = DensityPlot[sol[r, z], {r, z} \[Element] mesh,
ColorFunction -> "Rainbow", PlotRange -> All, Frame -> False,
AspectRatio -> Automatic];

frames = Table[
Show[dp,
ParametricPlot[
Table[{x[j][t], y[j][t]} /. sol1, {j, numbodies}], {t, 0, tf},
PlotRange -> {{0, r2}, {0, z1}}, Axes -> False],
Graphics[
Table[{Hue[.35], Disk[{x[j][tf], y[j][tf]} /. sol1, 0.0005]}, {j,
numbodies}]]], {tf, 0.01 tfin, tfin, .01 tfin}];
ListAnimate[frames] • To remove the jitter in the plot you could try to set the ImageSize explicitly. Either that or expand the PlotRange to add more padding on the left and right. – b3m2a1 Dec 22 '19 at 0:46