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.61773310^-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.5210^-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_] = qGrad[sol[x, y, z], {x, y, z}][[1]];
eforceY[x_, y_, z_] = qGrad[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, .025tfin, tfin, .025tfin}];
video = ListAnimate[frames, SaveDefinitions -> True]
x[t], governed byx'[t] == vx[t]for some expressionvx[t], use a discrete variablex'[t] == a[t] * vx[t], with initial valuea[0] == 1. Then useWhenEvent[x[t] > upperboundary, a[t] -> 0]to stopx[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. – Michael E2 Jul 05 '22 at 16:30WhenEventfor each possible boundary crossing, potentially as many as two times the number of coordinates times the number of particles.) – Michael E2 Jul 05 '22 at 16:48eqs, I madez1[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]
– Zach Jul 05 '22 at 17:53... This throws an error"The function a[j] appears as the head of the expression a[j][t]"`.{a[j][t] -> 0, z1[j]'[t] -> 0}since it's a 2nd order equation. Also thej -> ishould go outside:Table[WhenEvent[ ...] /. {j -> i}, {i, 2}]– Michael E2 Jul 05 '22 at 18:20