3
$\begingroup$

I am attempting to implement an event into my simulation which makes any moving particle in an electric field completely stop at the instant they pass some specified threshold.

I was trying to utilize DiscreteVariable, as presented in other threads to accomplish this:

Using WhenEvent to Bound Variables

Using WhenEvent to limit the derivative,

but can't seem to work it out in my particular application.

The boundary I am trying to implement is seen here under event:

eqs = Table[{x1[j]''[t] == -(1/mass)*eforceX[x, y, z] /. {x -> 
       x1[j][t], y -> y1[j][t], z -> z1[j][t]},
    y1[j]''[t] == -(1/mass)*eforceY[x, y, z] /. {x -> x1[j][t], 
      y -> y1[j][t], z -> z1[j][t]},
    z1[j]''[t] == -(1/mass)*eforceZ[x, y, z] /. {x -> x1[j][t], 
      y -> y1[j][t], 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}];

partRad = 0.025/2;

vars = Flatten[Table[{x1[j][t], y1[j][t], z1[j][t]}, {j, numbodies}]];

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

sol1 = NDSolve[{eqs, event}, Head /@ vars, {t, 0, tfin}]

How would one go about setting up an actual boundary condition, such that if one of the N particles passes the boundary z1[j][t] >= -partRadius, they completely stop?

The full program is shown here:

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, -6.5}, {2, 2, 0.5}], M], 
    Cuboid[{0, 0, 0.4}, {2, 2, 0.5}]], 
   Cuboid[{0, 0, -6.5}, {2, 2, -6.4}]];
ToElementMesh[r]["Wireframe"];
pol = -1;

V0 = 1000;
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, (-6.5 <= z <= -6.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, PlotLegends -> Automatic], 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, PlotLegends -> Automatic],
    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];


eforceX[x_, y_, z_] = q*Grad[sol[x, y, z], {x, y, z}][[1]];
eforceY[x_, y_, z_] = q*Grad[sol[x, y, z], {x, y, z}][[2]];
eforceZ[x_, y_, z_] = q*Grad[sol[x, y, z], {x, y, z}][[3]];


tfin = 400;

eqs = Table[{x1[j]''[t] == -(1/mass)*eforceX[x, y, z] /. {x -> 
       x1[j][t], y -> y1[j][t], z -> z1[j][t]},
    y1[j]''[t] == -(1/mass)*eforceY[x, y, z] /. {x -> x1[j][t], 
      y -> y1[j][t], z -> z1[j][t]},
    z1[j]''[t] == -(1/mass)*eforceZ[x, y, z] /. {x -> x1[j][t], 
      y -> y1[j][t], 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}];

partRad = 0.025/2;

vars = Flatten[Table[{x1[j][t], y1[j][t], z1[j][t]}, {j, numbodies}]];

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

sol1 = NDSolve[{eqs, event}, Head /@ vars, {t, 0, tfin}]


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, .025*tfin, tfin, .025*tfin}];
video = ListAnimate[frames, SaveDefinitions -> True]
$\endgroup$
6
  • 2
    $\begingroup$ Just a quick hint, which I hope is accurate: To stop a variable, say x[t], governed by x'[t] == vx[t] for some expression vx[t], use a discrete variable x'[t] == a[t] * vx[t], with initial value a[0] == 1. Then use WhenEvent[x[t] > upperboundary, a[t] -> 0] to stop x[t]. You'll need a discrete variable for each particle, and you'll need to use it to stop each coordinate of a given particle. $\endgroup$
    – Michael E2
    Jul 5, 2022 at 16:30
  • 1
    $\begingroup$ (You'll also need a WhenEvent for each possible boundary crossing, potentially as many as two times the number of coordinates times the number of particles.) $\endgroup$
    – Michael E2
    Jul 5, 2022 at 16:48
  • $\begingroup$ @MichaelE2 That makes sense to me, but I suppose implementation is what is confusing. For instance, in my set of eqs, I made z1[j]''[t] == a[j][t]*(-1/mass)*eforceZ[x, y, z] and event = Table[ WhenEvent[ z1[j][t] >= -partRad , {a[j][t] -> 0} /. {j -> i}], {i, numbodies}]; sol1 = NDSolve[{eqs, event}, Head /@ vars, {t, 0, tfin}, DiscreteVariables -> a]... This throws an error "The function a[j] appears as the head of the expression a[j][t]". $\endgroup$
    – Zach
    Jul 5, 2022 at 17:53
  • 1
    $\begingroup$ You need to set both {a[j][t] -> 0, z1[j]'[t] -> 0} since it's a 2nd order equation. Also the j -> i should go outside: Table[WhenEvent[ ...] /. {j -> i}, {i, 2}] $\endgroup$
    – Michael E2
    Jul 5, 2022 at 18:20
  • $\begingroup$ @MichaelE2 The answer is shown below. Thank you for your help! $\endgroup$
    – Zach
    Jul 5, 2022 at 19:41

1 Answer 1

3
$\begingroup$

Here is the solution. Thanks to @Michael E2 in original question comments.

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"];*)

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, -6.5}, {2, 2, 0.5}], M], 
    Cuboid[{0, 0, 0.4}, {2, 2, 0.5}]], 
   Cuboid[{0, 0, -6.5}, {2, 2, -6.4}]];
ToElementMesh[r]["Wireframe"];
pol = -1;

V0 = 1000;
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, (-6.5 <= z <= -6.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, PlotLegends -> Automatic], 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, PlotLegends -> Automatic],
    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];
discreteStart = Table[Partition[{1, 1, 1}, 1], numbodies];

eforceX[x_, y_, z_] = q*Grad[sol[x, y, z], {x, y, z}][[1]];
eforceY[x_, y_, z_] = q*Grad[sol[x, y, z], {x, y, z}][[2]];
eforceZ[x_, y_, z_] = q*Grad[sol[x, y, z], {x, y, z}][[3]];


tfin = 600;

eqs = Table[{x1[j]''[t] == -a[t]*(1/mass)*eforceX[x, y, z] /. {x -> 
       x1[j][t], y -> y1[j][t], z -> z1[j][t]},
    y1[j]''[t] == -a[t]*(1/mass)*eforceY[x, y, z] /. {x -> x1[j][t], 
      y -> y1[j][t], z -> z1[j][t]},
    z1[j]''[t] == -a[t]*(1/mass)*eforceZ[x, y, z] /. {x -> x1[j][t], 
      y -> y1[j][t], 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]],
    a[0] == 1},
   {j, numbodies}];

partRad = 0.0125;

vars = Flatten[Table[{x1[j][t], y1[j][t], z1[j][t]}, {j, numbodies}]];

(*event=Table[{WhenEvent[ z1[j][t]\[GreaterEqual] -partRad \
,{x1[j]'[t]\[Rule] 0,y1[j]'[t]\[Rule] 0,z1[j]'[t]\[Rule] 0}]}/.{j\
\[Rule]i},{i,numbodies}];*)

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

sol1 = NDSolve[{eqs, event}, Head /@ vars, {t, 0, tfin}, 
  DiscreteVariables -> a]


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, .015*tfin, tfin, .015*tfin}];
video = ListAnimate[frames, SaveDefinitions -> True]
```
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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