Appendix G
Clausing Factor Monte Carlo Calculation
Visual Basic Monte-Carlo calculation of Clausing Factor for thruster grids.
Inputs:
Clausing Factor Calculator Inputs thickScreen thickAccel rScreen rAccel gridSpace npart Radius (mm) 0.381 0.5 0.9525 0.5715 0.5 106 1.905 1.143 Diameter
Code:
Sub Clausing() thickScreen = Range("C4") thickAccel = Range("C5") rScreen = Range("C6") rAccel = Range("C7") gridSpace = Range("C8") npart = Range("C9") ' Monte Carlo Routine that calculates Clausing factor for CEX ' returns Clausing Factor and Downstream Correction factor Dim gone As Boolean Pi = 3.14159265358979 'assumes rTop = 1
483
484
Appendix G
rBottom = rScreen / rAccel lenBottom = (thickScreen + gridSpace) / rAccel lenTop = thickAccel / rAccel Length = lenTop + lenBottom iescape = 0 maxcount = 0 icount = 0 nlost = 0 vztot = 0# vz0tot = 0# For ipart = 1 To npart ' launch from bottom notgone = True r0 = rBottom * Sqr(Rnd) z0 = 0# costheta = Sqr(1# - Rnd) If (costheta > 0.99999) Then costheta = 0.99999 phi = 2 * Pi * Rnd sintheta = Sqr(1# - costheta ^ 2) vx = Cos(phi) * sintheta vy = Sin(phi) * sintheta vz = costheta rf = rBottom t = (vx * r0 + Sqr((vx ^ 2 + vy ^ 2) * rf ^ 2 - (vy * r0) ^ 2)) / (vx ^ 2 + vy ^ 2) z = z0 + vz * t vz0tot = vz0tot + vz icount = 0 Do While notgone icount = icount + 1 If (z < lenBottom) Then ' hit wall of bottom cylinder and is reemitted r0 = rBottom z0 = z costheta = Sqr(1# - Rnd) If (costheta > 0.99999) Then costheta = 0.99999 phi = 2 * Pi * Rnd sintheta = Sqr(1# - costheta ^ 2) vz = Cos(phi) * sintheta vy = Sin(phi) * sintheta vx = costheta rf = rBottom t = (vx * r0 + Sqr((vx ^ 2 + vy ^ 2) * rf ^ 2 - (vy * r0) ^ 2)) / (vx ^ 2 + vy ^ 2) z = z0 + t * vz End If ' bottom cylinder re-emission
Clausing Factor Monte Carlo Calculation
485
If ((z >= lenBottom) And (z0 < lenBottom)) Then ' emitted below but going up ' find radius at lenBottom t = (lenBottom - z0) / vz r = Sqr((r0 - vx * t) ^ 2 + (vy * t) ^ 2) If (r <= 1) Then ' continuing upward rf = 1# t = (vx * r0 + Sqr((vx ^ 2 + vy ^ 2) * rf ^ 2 - (vy * r0) ^ 2)) / (vx ^ 2 + vy ^ 2) z = z0 + vz * t Else ' hit the upstream side of the accel grid and is re-emitted downward r0 = r z0 = lenBottom costheta = Sqr(1# - Rnd) If (costheta > 0.99999) Then costheta = 0.99999 phi = 2 * Pi * Rnd sintheta = Sqr(1# - costheta ^ 2) vx = Cos(phi) * sintheta vy = Sin(phi) * sintheta vz = -costheta rf = rBottom t = (vx * r0 + Sqr((vx ^ 2 + vy ^ 2) * rf ^ 2 - (vy * r0) ^ 2)) / (vx ^ 2 + vy ^ 2) z = z0 + vz * t End If End If ' end upward If ((z >= lenBottom) And (z <= Length)) Then ' hit the upper cylinder wall and is re-emitted r0 = 1# z0 = z costheta = Sqr(1# - Rnd) If (costheta > 0.99999) Then costheta = 0.99999 phi = 2 * Pi * Rnd sintheta = Sqr(1# - costheta ^ 2) vz = Cos(phi) * sintheta vy = Sin(phi) * sintheta vx = costheta rf = 1# t = (vx * r0 + Sqr((vx ^ 2 + vy ^ 2) * rf ^ 2 - (vy * r0) ^ 2)) / (vx ^ 2 + vy ^ 2)
486
Appendix G
z = z0 + t * vz If (z < lenBottom) Then ' find z when particle bottom cylinder
hits
the
rf = rBottom If ((vx ^ 2 + vy ^ 2) * rf ^ 2 - (vy * r0) ^ 2 < 0#) Then t = (vx * r0) / (vx ^ 2 + vy ^ 2) 'if sqr arguement is less than 0 then set sqr term to 0 12 May 2004 Else t = (vx * r0 + Sqr((vx ^ 2 + vy ^ 2) * rf ^ 2 - (vy * r0) ^ 2)) / (vx ^ 2 + vy ^ 2) End If z = z0 + vz * t End If End If ' end upper cylinder emission If (z < 0#) Then notgone = False End If If (z > Length) Then iescape = iescape + 1 vztot = vztot + vz notgone = False End If If (icount > 1000) Then notgone = False icount = 0 nlost = nlost + 1 End If Loop ' while If (maxcount < icount) Then maxcount = icount Next ipart ' particles Range("C11") = (rBottom ^ 2) * iescape / npart Range("C12") = maxcount Range("C13") = nlost vz0av = vz0tot / npart vzav = vztot / iescape DenCor = vz0av / vzav ' Downstream correction factor End Sub ' Clausing
