0
$\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.

When it comes time to plot the motion of the particles and the electric field/aperture (seen above), I am running into some problems in the last few lines of the code

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

which produces these errors:

Errors

Here is the full code:

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

q = 1.60217733*10^-19*10;(*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 = 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]

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];





eforce1 = q*electricField1[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^-8;(*particle mass in kg/m^3*)
numbodies = 3;
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, 2}], RandomReal[{0, 2}], 
    RandomReal[{0.00001, 1}]}, 1], numbodies]

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

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


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


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 -> Automatic, 
     Axes -> False], 
    Graphics3D[
     Table[{Hue[.35], 
       Sphere[{x1[j][tf], y1[j][tf], z1[j][tf]} /. sol1, 0.0005]}, {j,
        numbodies}]]], {tf, 0.01 tfin, tfin, .01 tfin}];
ListAnimate[frames]
```
$\endgroup$
4
  • 1
    $\begingroup$ Your initial conditions are lists with 1 element, but numbers are expected. E.g. instead of x1[j][0] == pos0[[j, 1, ]] write x1[j][0] == pos0[[j, 1, 1]] $\endgroup$ Jun 28 at 18:17
  • $\begingroup$ @DanielHuber That was it! And of course it was something simple that I missed lol. Thank you. $\endgroup$
    – Zach
    Jun 28 at 18:38
  • $\begingroup$ @DanielHuber How do I mark a comment as the answer? I am new and apparently can't figure it out. $\endgroup$
    – Zach
    Jun 28 at 18:50
  • $\begingroup$ You can not. I will write an answer. $\endgroup$ Jun 28 at 18:52

1 Answer 1

1
$\begingroup$

Your initial conditions are lists with 1 element, but numbers are expected. E.g. instead of x1[j][0] == pos0[[j, 1, ]] write x1[j][0] == pos0[[j, 1, 1]]

$\endgroup$

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.