1
$\begingroup$

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.

Static Electric Field Through Aperture

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:

r(t) for particles

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["NDSolve`FEM`"]

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]
$\endgroup$
7
  • 1
    $\begingroup$ 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? $\endgroup$ Jun 29 at 19:07
  • $\begingroup$ the initial velocities were set to 0 just to test, so yes. How does one make eforce a function of position? $\endgroup$
    – Zach
    Jun 29 at 19:14
  • 1
    $\begingroup$ 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 & $\endgroup$ Jun 29 at 19:53
  • $\begingroup$ @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. $\endgroup$
    – Zach
    Jun 29 at 20:50
  • $\begingroup$ 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 $\endgroup$ Jun 30 at 8:36

1 Answer 1

0
$\begingroup$

Another attempt trying to call individual

Clear["Global`*"];
Needs["NDSolve`FEM`"]

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]
```
$\endgroup$
7
  • $\begingroup$ 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]} $\endgroup$ Jul 1 at 8:20
  • $\begingroup$ @DanielHuber doesn't seem to make a difference. It is still a constant force. $\endgroup$
    – Zach
    Jul 1 at 13:55
  • $\begingroup$ 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. $\endgroup$
    – Zach
    Jul 1 at 14:17
  • $\begingroup$ 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 $\endgroup$ Jul 1 at 16:26
  • $\begingroup$ @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! $\endgroup$
    – Zach
    Jul 1 at 16:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.