0
$\begingroup$

I have created a 2D axisymmetric simulation of the motion of charged particles in an electric field. The electric field is created by a needle attached to a top plate both of which are at -3800V with a grounded plate close below the needle. The needle is positioned above a hole in the plate which will allow particles to travel to the area below the grounded plate. The particles have an initial velocity, and in addition to the electric field, the forces acting on the particles include gravity, drag (approximated with Stokes Law), and the Coulomb force between particles.

It seems I cannot have more than two particles without it being too computationally intensive for my computer to solve (i.e. when there are 3 particles, it takes hours to solve for sol1 in the code below). Does anyone have any suggestions to improve the code and reduce computational time? Thanks!

ClearAll["Global`*"]
Needs["NDSolve`FEM`"]
q = -1.60217733*10^-19*10;(*particle charge*)
voltage = -3800;(*needle and top plate voltage*)
r1 = 0.0065; (*hole radius*)
r2 = 0.0365; (*domain radius*)
r3 = 0.00015; (*needle radius*)
z1 = 0.07; (*height of domain*)
z2 = 0.065; (*height of tip of needle*)
l = 0.005; (*distance between tip of needle and hole in middle plate*)
z3 = z2 - l; (*height of middle plate top surface*)
z4 = z3 - 0.0016; (*height of middle plate bottom surface*)
reg1 = ImplicitRegion[True, {{r, r3, r2}, {z, z2, z1}}]; (*region to the right of needle, above middle plate top surface*)
reg2 = ImplicitRegion[True, {{r, 0, r2}, {z, z3, z2}}]; (*region between tip of needle and middle plate top surface*)
reg3 = ImplicitRegion[True, {{r, 0, r1}, {z, z4, z3}}]; (*region to the left of middle plate*)
reg4 = ImplicitRegion[True, {{r, 0, r2}, {z, 0, z4}}]; (*region below middle plate bottom surface*)
region = RegionUnion[reg1, reg2, reg3, reg4]; (*merge all the previously created regions, space left out represents needle and middle plate*)
meshRefine[vertices_, area_] := area > 0.0000001;
mesh = ToElementMesh[DiscretizeRegion[region], MeshRefinementFunction -> meshRefine]; 
bc1 = {DirichletCondition[phi[r, z] == voltage, (z == z2 && 0 <= r <= r3)],
DirichletCondition[phi[r, z] == voltage, (z == z1 && r3 <= r <= r2)],
DirichletCondition[phi[r, z] == voltage, (r == r3 && z2 <= z <= z1)],
DirichletCondition[phi[r, z] == 0, (z == z3 && r1 <= r <= r2)],
DirichletCondition[phi[r, z] == 0, (z == z4 && r1 <= r <= r2)],
DirichletCondition[phi[r, z] == 0, (r == r1 && z4 <=z <=z3)]};  (*boundary conditions*)
sol = NDSolveValue[{1/r*D[r*D[phi[r, z], r], r] + 
      D[phi[r, z], z, z] == 0, bc1}, phi, {r, z} \[Element] mesh];
electricField = -Evaluate[Grad[sol[r, z], {r, z}]];
eforce = q*electricField;
dreal = 25*10^-9; (*real particle diameter in m*)
d = 7.5*10^-5; (*particle diameter in m (for simulation/visual purposes)*)
rho = 998; (*real particle density in kg/m^3*)
mu = 1.81*10^-5; (*air dynamic viscosity in Pa*s *)
volume = N[4/3*Pi*(dreal/2)^3]; (*real particle volume in m^3*)
mass = volume*rho;(*real particle mass in kg*)
VFR = 10/60/60*10^-6*0.001; (*volumetric flow rate through the needle in m^3/s*)
ivel = VFR/(Pi*r3^2); (*particle initial velocity out of the needle in m/s*)
numbodies = Round[r3/d]; (*number of particles that fit across the radius of the needle*)
gravity = -9.81; (*acceleration due to gravity in m/s^2*)
ke = 9*10^9; (*Coulomb constant*)
vel0 = Table[{0, -ivel}, numbodies];
pos0 = Riffle[Range[0, r3, r3/numbodies], z2]~Partition~2;
eqs = Table[{x[j]''[t] == 1/mass*Sum[((x[j][t] - x[i][t])/Sqrt[(x[j][t] - x[i][t])^2 
+ (y[j][t] - y[i][t])^2])*(ke*q^2/Sqrt[(x[j][t] - x[i][t])^2 + (y[j][t] - y[i][t])^2]),
 {i, Delete[Range[numbodies], j]}] + 1/mass*6.0*Pi*mu*dreal/2*- x[j]'[t]
 + 1/mass*eforce[[1]] /. {r -> x[j][t], z -> y[j][t]}, 
y[j]''[t] == gravity + 1/mass*Sum[((x[j][t] - x[i][t])/Sqrt[(x[j][t] - x[i][t])^2 
+ (y[j][t] - y[i][t])^2])*(ke*q^2/Sqrt[(x[j][t] - x[i][t])^2 + (y[j][t] - y[i][t])^2]), 
{i, Delete[Range[numbodies], j]}] + 1/mass*6.0*Pi*mu*dreal/2*-y[j]'[t] 
+ 1/mass*eforce[[2]] /. {r -> x[j][t], z -> y[j][t]}, x[j][0] == pos0[[j, 1]], 
y[j][0] == pos0[[j, 2]], x[j]'[0] == vel0[[j, 1]], y[j]'[0] == vel0[[j, 2]]}, {j, numbodies}];
vars = Flatten[Table[{x[j], y[j]}, {j, numbodies}]];
event = Table[{WhenEvent[x[j][t] == 0, {x[j]'[t] -> 0, y[j]'[t] -> 0}],
WhenEvent[x[j][t] == 0, {x[j]'[t] -> 0, y[j]'[t] -> 0}], 
WhenEvent[x[j][t] == r2, {x[j]'[t] -> 0, y[j]'[t] -> 0}], 
WhenEvent[y[j][t] == z1, {x[j]'[t] -> 0, y[j]'[t] -> 0}], 
WhenEvent[y[j][t] == z3 && r1 <= x[j][t] <= r2, {x[j]'[t] -> 0, y[j]'[t] -> 0}], 
WhenEvent[y[j][t] == z4 && r1 <= x[j][t] <= r2, {x[j]'[t] -> 0, y[j]'[t] -> 0}], 
WhenEvent[x[j][t] == r1 && z4 <= y[j][t] <= z3, {x[j]'[t] -> 0, y[j]'[t] -> 0}]} /. j -> i, {i, numbodies}];
tfin = 25;
sol1 = NDSolve[{eqs, event}, vars, {t, 0, tfin}][[1]];
dp = DensityPlot[sol[r, z], {r, z} \[Element] mesh, ColorFunction -> "Rainbow", 
PlotLegends -> Automatic, PlotRange -> All, Frame -> False, AspectRatio -> Automatic];
frames = Table[Show[dp, ParametricPlot[Table[{x[j][t], y[j][t]} /. sol1, {j, numbodies}], {t, 0, tf}, 
PlotRange -> {{0, r2}, {0, z1}}, Axes -> False], Graphics[Table[{Hue[.35], 
Disk[{x[j][tf], y[j][tf]} /. sol1, d]}, {j, numbodies}]]], {tf, 0.01 tfin, tfin, .01 tfin}];
ListAnimate[frames]
$\endgroup$
  • $\begingroup$ This is fast on my computer. What do I need to do to see the slowness? Also, there is no sol1 only sol in the code you provided. Could you clarify? $\endgroup$ – user21 Jan 8 at 14:39
  • $\begingroup$ @user21 On mine (MMA 12, Win10-64) I let the code run for a few minutes then aborted it as nothing was forthcoming. The bottleneck was NDSolve, so I don't know about the plot stuff after that. $\endgroup$ – MarcoB Jan 8 at 15:08
  • $\begingroup$ @MarcoB, hm, what happens when use this: mesh = ToElementMesh[region, "MaxCellMeasure" -> 0.0001]; and make the MCM smaller? $\endgroup$ – user21 Jan 8 at 15:22
  • 2
    $\begingroup$ @kjcole: Please simplify your code to highlight the problem!! The difficult point is the second NDSolve, but the structure of your problem is completely lost in the morass of code. Split it up in chunks, add verbal explanations of your purpose, remove the plotting stuff that is irrelevant, etc etc. $\endgroup$ – MarcoB Jan 8 at 15:45
  • 1
    $\begingroup$ @kjcole, you absolutely need to simplify your code. $\endgroup$ – user21 Jan 8 at 16:16

Your Answer

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

Browse other questions tagged or ask your own question.