# Properly Updating Particle Motion in 3D Electric Field Simulation

I have a code (very bottom of post) which plots a static electric field as it passed through some metal sheet with an aperture in it.

I seek to observe the motion of some N charged particles (point charges) randomly dispersed near the aperture-field (seen above). I have created a set of equations for the N particles in x[t], y[t] and z[t] as they sit in the electric field, and can solve the set of differential equations which exist for each particle.

I can plot the particles with my aperture and show SOME motion, but it isn't the correct motion as they all share the same vector. Their trajectories also remain constant in time, as seen here:

I suspect I have chosen some point in the electric field with some magnitude and have passed it to all the particles, but have failed to tell them to evaluate the electric field strength at each point in time/space in my set of equations:

eqs = Table[{x1[j]''[t] == -(1/mass)*eforce[[1, 1]] /. {x -> x1[j][t]},
y1[j]''[t] == -(1/mass)*eforce[[2, 1]] /. {y -> y1[j][t]},
z1[j]''[t] == -(1/mass)*eforce[[3, 1]] /. {z -> z1[j][t]},
x1[j][0] == pos0[[j, 1, 1]],
y1[j][0] == pos0[[j, 2, 1]],
z1[j][0] == pos0[[j, 3, 1]],
x1[j]'[0] == vel0[[j, 1, 1]],
y1[j]'[0] == vel0[[j, 2, 1]],
z1[j]'[0] == vel0[[j, 3, 1]]},
{j, numbodies}];


Here is the full code:

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

q = 1.617733*10^-18;(*Net ion Charge*)

R= Import["https://www.dropbox.com/s/dds8rm3odg2m7gu/largeAp.obj?dl=\
1"];

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}]];
ToElementMesh[r]["Wireframe"];
pol = -1;

V0 = 10000;
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];

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.5, 0.1},
PlotTheme -> "Detailed", ColorFunction -> "Rainbow",
PerformanceGoal -> "Quality", VectorScale -> 0.05,
VectorPoints -> 7], M];

eforce = -q*Grad[sol[x, y, z], {x, y, z}];

vecForce =
Show[VectorPlot3D[
q*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];

mass = 6.52*10^-10;(*particle mass in kg/m^3*)
numbodies = 3;

vel0 = Table[Partition[{0, 0, 0}, 1], numbodies];
(*vel0=Table[Partition[{RandomReal[{-0.0001,0.0001}],RandomReal[{-0.\
0001,0.0001}],RandomReal[{-0.0001,0.0001}]},1],numbodies]*)
pos0 = Table[
Partition[{RandomReal[{0.5, 1}], RandomReal[{0.5, 1}],
RandomReal[{-0.4, -0.00001}]}, 1], numbodies];

eqs = Table[{x1[j]''[t] == -(1/mass)*eforce[[1, 1]] /. {x -> x1[j][t]},
y1[j]''[t] == -(1/mass)*eforce[[2, 1]] /. {y -> y1[j][t]},
z1[j]''[t] == -(1/mass)*eforce[[3, 1]] /. {z -> z1[j][t]},
x1[j][0] == pos0[[j, 1, 1]],
y1[j][0] == pos0[[j, 2, 1]],
z1[j][0] == pos0[[j, 3, 1]],
x1[j]'[0] == vel0[[j, 1, 1]],
y1[j]'[0] == vel0[[j, 2, 1]],
z1[j]'[0] == vel0[[j, 3, 1]]},
{j, numbodies}];
vars = Flatten[Table[{x1[j], y1[j], z1[j]}, {j, numbodies}]];

event = Table[{WhenEvent[
z1[j][t] == 0, {x1[j]'[t] -> 0, y1[j]'[t] -> 0,
z1[j]'[t] -> 0}]} /. j -> i, {i, numbodies}];

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

plotXZ = ContourPlot[sol[x, 0.75, z], {x, 0, 2}, {z, -0.4, 0.1},
ContourShading -> Automatic, ColorFunction -> "Rainbow",
Contours -> 100];

frames = Table[
Show[v, ParametricPlot3D[
Table[{x1[j][t], y1[j][t], z1[j][t]} /. sol1, {j,
numbodies}], {t, 0, tf}, PlotRange -> All, Axes -> Off],
Graphics3D[
Table[{Hue[.35],
Sphere[{x1[j][tf], y1[j][tf], z1[j][tf]} /. sol1, 0.03]}, {j,
numbodies}]]], {tf, 0.01 tfin, tfin, .01 tfin}];
video = ListAnimate[frames]

• eforce[1,1],eforce[2,1],eforce[3,1] are constants instead of functions of position. Further vel0, all velocities are zero, is this intended? Jun 29, 2022 at 19:07
• the initial velocities were set to 0 just to test, so yes. How does one make eforce a function of position?
– Zach
Jun 29, 2022 at 19:14
• eforce[[1]] is a function of x,y,z. But eforce[[1,1]] is just the coefficient. You will have to write eforce[[1]] /. {x->x -> x1[j][t],y->y -> y1[j][t],z->z -> z1[j][t]}. But it will be easier if you define eforce as as functions like e.g. eforce [[1]]= Function[{x,y,z}, ....] or shorter: eforce[[1]]= expression of #1,#2,#3 & Jun 29, 2022 at 19:53
• @DanielHuber Thanks. I see what you mean, but I am a little lost regarding how I would make each of my equations:{x1,y1,z1[j]''[t] == -(1/mass)*eforce[[1]], independently functions of x, y and z by doing so? To clarify, I need x1[j][t] to only look at x components of field, y1[j][t] to only look at y components of the field, etc.
– Zach
Jun 29, 2022 at 20:50
• If  eforce[[1]] /. {x->x -> x1[j][t],y->y -> y1[j][t],z->z -> z1[j][t]} is a vector, then ( eforce[[1]] /. {x->x -> x1[j][t],y->y -> y1[j][t],z->z -> z1[j][t]})[[1]] is its x-component Jun 30, 2022 at 8:36

Another attempt trying to call individual

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

q = 1.617733*10^-18;(*Net ion Charge*)

R= Import["https://www.dropbox.com/s/dds8rm3odg2m7gu/largeAp.obj?dl=\
1"];

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}]];
ToElementMesh[r]["Wireframe"];
pol = -1;

V0 = 15000;
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];

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.5, 0.1},
PlotTheme -> "Detailed", ColorFunction -> "Rainbow",
PerformanceGoal -> "Quality", VectorScale -> 0.05,
VectorPoints -> 7], M];

vecForce =
Show[VectorPlot3D[
q*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];

numbodies = 3;
mass = 6.52*10^-11;(*particle mass in kg/m^3*)

vel0 = Table[Partition[{0, 0, 0}, 1], numbodies];
(*vel0=Table[Partition[{RandomReal[{-0.0001,0.0001}],RandomReal[{-0.\
0001,0.0001}],RandomReal[{-0.0001,0.0001}]},1],numbodies]*)
pos0 = Table[
Partition[{RandomReal[{0.5, 1}], RandomReal[{0.5, 1}],
RandomReal[{-0.4, -0.04}]}, 1], numbodies];

eforce1 = -q*(Grad[sol[x, y, z], {x, y, z}]);

eforceX = (eforce1[[1]] /. {x -> x -> x1[j][t], y -> y -> y1[j][t],
z -> z -> z1[j][t]}[[1]])
eforceY = (eforce1[[1]] /. {x -> x -> x1[j][t], y -> y -> y1[j][t],
z -> z -> z1[j][t]}[[2]])
eforceZ = (eforce1[[1]] /. {x -> x -> x1[j][t], y -> y -> y1[j][t],
z -> z -> z1[j][t]}[[3]])

eqs = Table[{x1[j]''[t] == -(1/mass)*eforceX[[1]],
y1[j]''[t] == -(1/mass)*eforceY[[1]],
z1[j]''[t] == -(1/mass)*eforceZ[[1]],
x1[j][0] == pos0[[j, 1, 1]],
y1[j][0] == pos0[[j, 2, 1]],
z1[j][0] == pos0[[j, 3, 1]],
x1[j]'[0] == vel0[[j, 1, 1]],
y1[j]'[0] == vel0[[j, 2, 1]],
z1[j]'[0] == vel0[[j, 3, 1]]},
{j, numbodies}];
vars = Flatten[Table[{x1[j], y1[j], z1[j]}, {j, numbodies}]];

Print[eforce1, {j, numbodies}, {t, 0, tf}]

event = Table[{WhenEvent[
z1[j][t] == -0.025, {x1[j]'[t] -> 0, y1[j]'[t] -> 0,
z1[j]'[t] -> 0}]} /. j -> i, {i, numbodies}];

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

plotXZ = ContourPlot[sol[x, 0.75, z], {x, 0, 2}, {z, -0.4, 0.1},
ContourShading -> Automatic, ColorFunction -> "Rainbow",
Contours -> 100];

frames = Table[
Show[v, ParametricPlot3D[
Table[{x1[j][t], y1[j][t], z1[j][t]} /. sol1, {j,
numbodies}], {t, 0, tf}, PlotRange -> All, Axes -> Off],
Graphics3D[
Table[{Hue[.35],
Sphere[{x1[j][tf], y1[j][tf], z1[j][tf]} /. sol1, 0.025]}, {j,
numbodies}]]], {tf, 0.025 tfin, tfin, .025 tfin}];
video = ListAnimate[frames]
$$$$

• Print eforceX and you will see that is still a constant. But it should be a function of t. You should write. Then "x -> x ->..." is wrong. Try: eforceX = eforce1[[1]] /. {x -> x1[j][t], y -> y1[j][t], z -> z1[j][t]} Jul 1, 2022 at 8:20
• @DanielHuber doesn't seem to make a difference. It is still a constant force.
– Zach
Jul 1, 2022 at 13:55
• I changed vars = Flatten[ Table[{x1[j], y1[j], z1[j]}, {j, numbodies}]]; to be vars = Flatten[ Table[{x1[j][t], y1[j][t], z1[j][t]}, {t, 0, tfin}, {j, numbodies}]]; and it produced a table of x1,y1,z1[j][t] values for t: 0->tfin, which is promising. However, when it comes time to plot in the last few lines there must be some issue, as it is producing a red 3d plot. Error says: Coordinate {\$CellContext x1[1][200.],...} should be a triple of number, or a Scaled form.
– Zach
Jul 1, 2022 at 14:17
• If I use eforceX = eforce1[[1]] /. {x -> x1[j][t], y -> y1[j][t], z -> z1[j][t]} and then print eforceX, I get: -1.61773*10^-18 -16 InterpolatingFunction[{{0., 2.}, {-1.11022 10 , 2.}, {-0.4, 0.4}}, <>][x1[j][t], y1[j][t], z1[j][t]]` what is a function of "t", NOT a constant Jul 1, 2022 at 16:26
• @DanielHuber I guess I meant that although eforceX became a function of t, I was still observing a constant force, which is most likely due to my sol1 definition. So I changed vars and now this new problem has arose. But it does seem like a positive change in that vars spits out a table of all the positions of the particles at each time step!
– Zach
Jul 1, 2022 at 16:49