Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Particle simulations of colloid thruster beam
(USC Thesis Other)
Particle simulations of colloid thruster beam
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
PARTICLESIMULATIONSOFCOLLOIDTHRUSTERBEAM by Yinjian Zhao A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree of DOCTOR OF PHILOSOPHY (ASTRONAUTICAL ENGINEERING) December 2019 Copyright 2019 Yinjian Zhao ii EPIGRAPH Epigraph I doubt therefore I think, I think therefore I am. | Ren e Descartes ACKNOWLEDGMENTS iii Acknowledgments First, I would like to thank my parents for bringing me into this fantastic world, so I can sail forbidden seas of the mysterious and fascinating universe. I really appreciate their endless support and love throughout my entire life. Then, I must thank my beloved wife Yanan Zhang. Falling in love with her is God's best gift to me. She is like a light in the darkness, making me strong, courageous, and fearless, leading me to a colorful world full of happiness and joy. Next, I would like to thank my dear advisor Prof. Joseph J. Wang for his great eort put into helping me build my academic knowledge and enhance my research skills. His kindness, amiableness, and generousness also taught me how to become a better person. I appreciate all my labmates, Dr. Daoru Han, Dr. Scott Hughes, Dr. Kevin Chou, Dr. William Yu, Dr. Yuan Hu, Mr. Daniel Depew, Mr. Chen Cui, Mr. Robert Antypas, Mr. Jerey Asher, and Mr. Ziyu Huang. Without their help, I would not be able to overcome many obstacles occurred during my studies. I would like to thank Prof. Hideyuki Usui. I enjoyed our collaboration very much. I appreciate the advice and suggestion from Prof. Joseph Kunc and Prof. Aiichiro Nakano. I appreciate the kind help from Prof. Mike Gruntman and Prof. Daniel A. Erwin as well. I would like to say thanks to all the sta, Dell Causon, Linda Ly, Marlyn Lat, and Nicole Valdez, for their hard work to function the department. I would also like to thank all my friends who helped or accompanied me during my Ph.D. life. At last, I really appreciate the support from the Center for High-Performance Computing of USC for providing free resources to students. iv ABSTRACT Abstract Colloid thrusters have been widely utilized on small spacecrafts, because of their high thrust density, specic impulse, eciency, and accurate attitude control over long periods of time. Although the physics of an electried liquid surface has been strongly studied, the knowledge of the colloid thruster beam is very limited due to the lack of precise experimental measure- ments and numerical simulations. In this dissertation, state-of-the-art particle simulation models of the colloid thruster beam are developed: The Particle-Particle (PP) model aims at a high accuracy, such that its results can be treated as a baseline to provide guidance for other models; The Particle-Particle Particle-Mesh (PP-PM) hybrid model, which aims at both high accuracy and high simulation speed, is capable for simulating larger and more realistic cases. The commonly known computationally expensive PP model is rst successfully applied to simulate charged droplets in colloid thruster beam. Due to the high-accuracy nature of the PP model with respect to solving particle collisions, the most trustable results of droplet dynamics are obtained. When a relatively large number of droplets are involved, OpenMP or GPU CUDA parallelization can be easily implemented with high eciency on the PP model, because of its simple algorithm, directly applying Coulomb's law. The details of the PP model is described in the dissertation, some characteristic simulations considering droplet breakup are carried out. In addition, a novel concept of equivalent length is introduced to re ect the extent of particle short-range interactions (collisions). From the results of the PP model, the probability density function of the equivalent length can be obtained, which suggests a suciently minimum mesh size for a corresponding PM model that can capture all the necessary particle collisions. The space charge eects of the colloid thrusters under the parameter range simulated are studied. The PP algorithm and code are veried through a comparison with the Child- Langmuir law and a Particle-in-Cell (PIC) method. It is found that the Child-Langmuir law is inappropriate to be applied for the scenario of colloid thrusters. Possible reasons are given and discussed. In terms of the charge neutralization, the alternating current mode of a colloid thruster is studied using the PP model for the rst time. Dierent types of the alternating current functions are considered, and the transition time between every alternation is studied in detail. Also, a bipolar operation mode is simulated, in which two thrusters are placed side by side operating in an opposite polarity to avoid the charge accumulation on the spacecraft. ABSTRACT v Dierent parameters are considered to vary the space charge eects, and the simulation results indicate that the space charge eects are negligible, so the two owing beams are not aecting each other. Based on the development of the PP model, a PP-PM hybrid model is proposed. It applies the PP model on the upstream region of the colloid thruster beam, where the beam is dense and the scale of the particle interactions is short, to maintain the high accuracy of the simulation. And, it applies the PM model on the downstream region of the beam, where the beam is diuse and the scale of the particle interactions is long, to speed up the simulation. Due to the axially symmetric nature of the beam, the downstream PM model is developed to be 2-D cylindrical with nonuniform mesh. Comparisons between the hybrid model and the pure PP model indicate the validity of the hybrid model. Also, the suciently minimum mesh size found previously is veried using the hybrid model. At last, the hybrid model is compared with a multi-level Particle-in-Cell (PIC) model applied in a reference. It is found that the hybrid model produces a beam that expands more radially than that of the PIC model, due to the fact that it resolves the short-range particle interactions more accurately in the upstream region. Contents Epigraph ii Acknowledgments iii Abstract iv List of Figures x List of Tables xvii List of Symbols xviii Nomenclature xx 1 Introduction 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Basic Physics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.1 Starting Voltage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.2 Current and Charge . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.3 Limitations of Droplet Specic Charge . . . . . . . . . . . . . . . . . 4 1.2.4 Droplet Breakup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Motivations and Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 vi CONTENTS vii 1.4.1 Relevant Computational Works . . . . . . . . . . . . . . . . . . . . . 7 1.4.2 Relevant Charge Neutralization Works . . . . . . . . . . . . . . . . . 9 1.5 Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2 PP Simulation Model 12 2.1 Simulation Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 Input and Derived Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3 Initialization of droplets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.4 Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4.1 Particle-Particle Field . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4.2 Accelerating Electric Field . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4.3 Numerical Integration . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.5 Droplet Initial Breakup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.6 Simulation Procedure and Techniques . . . . . . . . . . . . . . . . . . . . . . 22 2.7 Diagnoses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.7.1 System Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.7.2 Equivalent length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.7.3 Droplet Macro-temperature . . . . . . . . . . . . . . . . . . . . . . . 26 2.7.4 Thrust and Specic Impulse . . . . . . . . . . . . . . . . . . . . . . . 27 3 PP Code Validation 28 3.1 Time Step Length Verication . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2 Space Charge Eect Validation . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.2.1 Child-Langmuir Law . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.2.2 The PIC Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.2.3 The PP Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4 Characteristic Simulation of Droplet Beam 37 viii CONTENTS 4.1 Chosen Parameters and Simulation Cases . . . . . . . . . . . . . . . . . . . . 37 4.2 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.2.1 Phase Space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.2.2 System Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.2.3 Number Density and Macro-temperature . . . . . . . . . . . . . . . . 44 4.2.4 Velocity Distribution Function . . . . . . . . . . . . . . . . . . . . . . 47 4.2.5 PDF of Equivalent Length . . . . . . . . . . . . . . . . . . . . . . . . 54 4.3 From PP to PM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.4 Initial Droplet Thermal Velocity . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.5 Space Charge Eects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.5.1 Calculation Based on Child-Langmuir Law . . . . . . . . . . . . . . . 64 4.5.2 Model of Two Parallel Plates . . . . . . . . . . . . . . . . . . . . . . 65 4.5.3 Model of Taylor Cone . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.5.4 The 1D Assumption of Child-Langmuir law . . . . . . . . . . . . . . 67 4.5.5 Conclusion for the Space Charge Eects . . . . . . . . . . . . . . . . 68 5 Colloid Thruster Neutralization of Droplet Beam 69 5.1 AC Mode Simulation Setup and Procedures . . . . . . . . . . . . . . . . . . 70 5.1.1 Type 1: Square Wave . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.1.2 Type 2: Step Function . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.2 Bipolar Operation Mode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.2.1 Parameter Calculation of Limitations . . . . . . . . . . . . . . . . . . 77 5.2.2 Simulations of Bipolar System . . . . . . . . . . . . . . . . . . . . . . 78 6 PP-PM Hybrid Model of Ion Beam 90 6.1 2-D Cylindrical Poisson Solver with Nonuniform Mesh . . . . . . . . . . . . . 90 6.1.1 Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 6.2 Charge Deposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 CONTENTS ix 6.2.1 Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 6.3 Particle Mapping and Push . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.4 PP-PM Hybrid Model Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 7 Hybrid Model Application of Ion Beam 103 7.1 Simulation Setup and Parameters . . . . . . . . . . . . . . . . . . . . . . . . 103 7.2 Accelerating Field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 7.3 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 7.4 Reproduction of Larger Expansion Angle . . . . . . . . . . . . . . . . . . . . 110 8 Summary and Future Work 113 Bibliography 117 List of Figures 2.1 Schematic diagram of a colloid thruster and the simulation domain. . . . . . 12 2.2 Injection volume. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3 The prolate spheroidal coordinate system. . . . . . . . . . . . . . . . . . . . 17 2.4 Schematic diagram of the needle and the Taylor cone. . . . . . . . . . . . . 19 2.5 Potential plot with R n =d = 0:2 and y = 0. . . . . . . . . . . . . . . . . . . . 19 2.6 Electric eld plots. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.7 Zoomed-in potential plot with electric eld vectors. . . . . . . . . . . . . . . 20 2.8 Simulation procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.9 Computation time comparison between using CPU and GPU. . . . . . . . . 23 2.10 Diagram for understanding the equivalent length. . . . . . . . . . . . . . . . 24 2.11 Probability density function of equivalent length. . . . . . . . . . . . . . . . 25 2.12 Components of the resultant force. . . . . . . . . . . . . . . . . . . . . . . . 26 3.1 Binary collision setup. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2 Simulated binary collision of Case 1. . . . . . . . . . . . . . . . . . . . . . . 29 3.3 The potential distribution (left), the electric eld distribution (right). . . . 31 3.4 The potential distribution (left), the electric eld distribution (right). . . . 31 3.5 The modied PP eld. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.6 The potential distribution (left), the electric eld distribution (right). . . . 33 3.7 The potential distribution (left), the electric eld distribution (right). . . . 34 x LIST OF FIGURES xi 3.8 The electric eld distribution of PIC (left), the electric eld distribution of PP (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.9 The potential distribution of PIC (left) and PP (right). . . . . . . . . . . . 35 3.10 The electric eld distribution of PIC (left) and PP (right). . . . . . . . . . . 35 4.1 Droplet spatial distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.2 Phase space plots of x-z of all A cases at time 2000= _ N d . . . . . . . . . . . . . 40 4.3 Phase space plots of z-v z of all A cases at time 2000= _ N d . . . . . . . . . . . . 41 4.4 Phase space plots of z-v xy of all A cases at time 2000= _ N d . . . . . . . . . . . . 41 4.5 Phase space plots of r xy -v xy of all A cases at time 2000= _ N d . . . . . . . . . . . 41 4.6 Phase space plots of r xy -v z of all A cases at time 2000= _ N d . . . . . . . . . . . 42 4.7 Phase space plots of v xy -v z of all A cases at time 2000= _ N d . . . . . . . . . . . 42 4.8 Time history of system energy. . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.9 Droplet number density n (m 3 ) of A 0 to A 4 ve cases, averaged from t _ N d = 1820 to 2000. The value of log 10 n is shown using colors. . . . . . . . . . . . 44 4.10 Droplet macro-temperature in axial direction T z (K) of A 0 to A 4 ve cases, averaged from t _ N d = 1820 to 2000. The value of log 10 (T z N s ) is shown using colors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.11 Droplet macro-temperature in radial direction T r (K) of A 0 to A 4 ve cases, averaged from t _ N d = 1820 to 2000. The value of log 10 (T r N s ) is shown using colors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.12 Droplet macro-temperature in tangential direction T t (K) of A 0 to A 4 ve cases, averaged from t _ N d = 1820 to 2000. The value of log 10 (T t N s ) is shown using colors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.13 Illustration of the three regions used for evaluating f v and f p . . . . . . . . . 47 4.14 Time history of f v (v z ) (top three panels), f v (v r ) (middle three panels), and f v (v t ) (bottom three panels) of case A 4 . The left, middle, and right columns show panels evaluated from region (1), (2), and (3), respectively. The values of f v have been normalized to range from 0 to 1. . . . . . . . . . . . . . . . 48 4.15 Time averaged plots of f v (v z ) (top), f v (v r ) (middle), and f v (v t ) (bottom) of case A 4 . l 1 , l 2 , and l 3 indicate the three regions evaluated. The values of f v have been normalized to range from 0 to 1. . . . . . . . . . . . . . . . . . . 49 xii LIST OF FIGURES 4.16 Time averaged plots off v (v z ) of all cases, evaluated in region (1), (2), and (3) shown in the top, middle, and bottom panel, respectively. The values of f v have been normalized to range from 0 to 1. . . . . . . . . . . . . . . . . . . 50 4.17 Time averaged plots off v (v r ) of all cases, evaluated in region (1), (2), and (3) shown in the top, middle, and bottom panel, respectively. CaseA 0 to caseA 4 corresponds to the line from the left to the right for all panels. The values of f v have been normalized to range from 0 to 1. . . . . . . . . . . . . . . . . . 51 4.18 Time averaged plots off v (v t ) of all cases, evaluated in region (1), (2), and (3) shown in the top, middle, and bottom panel, respectively. The values of f v have been normalized to range from 0 to 1. . . . . . . . . . . . . . . . . . . 52 4.19 Time history of f p ( ez ) (top three panels), f p ( er ) (middle three panels), and f p ( et ) (bottom three panels) of caseA 4 . The left, middle, and right columns show panels evaluated from region (1), (2), and (3), respectively. The values of f p have been normalized to range from 0 to 1. . . . . . . . . . . . . . . . 53 4.20 Time averaged plots of f p ( ez ) (top), f p ( er ) (middle), and f p ( et ) (bottom) of case A 4 . l 1 , l 2 , and l 3 indicate the three regions evaluated. The values of f p have been normalized to range from 0 to 1. . . . . . . . . . . . . . . . . . 54 4.21 Time averaged plots of f p ( ez ) of all cases, evaluated in region (1), (2), and (3) shown in the top, middle, and bottom panel, respectively. The values of f p have been normalized to range from 0 to 1. . . . . . . . . . . . . . . . . . 55 4.22 Time averaged plots of f p ( er ) of all cases, evaluated in region (1), (2), and (3) shown in the top, middle, and bottom panel, respectively. The values of f p have been normalized to range from 0 to 1. . . . . . . . . . . . . . . . . . 56 4.23 Time averaged plots of f p ( et ) of all cases, evaluated in region (1), (2), and (3) shown in the top, middle, and bottom panel, respectively. The values of f p have been normalized to range from 0 to 1. . . . . . . . . . . . . . . . . . 57 4.24 The mean interparticle distance of the case A 4 with t _ N = 1=(2N s ) in the unit of m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.25 Operation comparison between PP and PM. . . . . . . . . . . . . . . . . . . 59 4.26 An example of the operation comparison between PP and PM in the down- stream beam region. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.27 Time history of system energy for comparing dierent initial thermal veloci- ties. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.28 PDF of the equivalent length for comparing dierent initial thermal velocities. 61 4.29 Phase space x-z plot of the cases (A 0 ) with and without F pp (left), and an amplied view of the case withoutF pp (right). . . . . . . . . . . . . . . . . 62 LIST OF FIGURES xiii 4.30 Phase space z-v xy plot of the cases (A 0 ) with and without F pp (left), and an amplied view of the case withoutF pp (right). . . . . . . . . . . . . . . . . 63 4.31 Phase space z-v z plot of the cases (A 0 ) withoutF pp . . . . . . . . . . . . . . 63 4.32 Phase space plots of x-z for comparing dierent initial thermal velocities. A 1 cases are shown in the rst row; A 3 cases are shown in the second row. . . . 64 4.33 The electric eld distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.1 Phase spacex-z plots of the AC case: all droplets (top left); only recombined droplets (bottom left). The corresponding DC case: all droplets (top right); only recombined droplets (bottom right). Positive droplets are denoted by red dots; negative droplets are denoted by blue dots. The plotting simulation time is t _ N d = 16000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.2 Left: Phase space z-v z plot of the AC case. Right: The corresponding DC case. Positive droplets are denoted by red dots; negative droplets are denoted by blue dots. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.3 Time history of droplet exit speed of the square wave case (averaged in the region from0:1d 0 to 0:1d 0 ). . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.4 Time history of the total system energy. . . . . . . . . . . . . . . . . . . . . 72 5.5 Electric potential (in volts) caused by charged droplets. The AC case is shown on the left panel, while the corresponding DC case is shown on the right panel. 73 5.6 Time history of droplet exit speed of the step function case (averaged in the region from0:1d 0 to 0:1d 0 ). . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.7 Left panel: Time history of droplet exit speed (averaged in the region from 0:1d 0 to 0:1d 0 ) for cases with varied volume ow rates. The curves from the top to the bottom correspond to an increasing the volume ow rate. Right panel: Dependence of the transition time on the volume ow rate. . . . . . 75 5.8 Left panel: Time history of droplet exit speed (averaged in the region from 0:1d 0 to 0:1d 0 ) for cases with varied accelerating potentials. The curves from the top to the bottom correspond to an decreasing accelerating potential. Right panel: Dependence of the transition time on the accelerating potential. 75 5.9 Left panel: Time history of droplet exit speed (averaged in the region from 0:1d 0 to 0:1d 0 ) for cases with variedd 0 . The curves from the left to the right correspond to an increasing d 0 . Right panel: Dependence of the transition time on d 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 xiv LIST OF FIGURES 5.10 Thex-z phase space plot with separationx=d 0 = 0:4 of all droplets (left) and only recombined droplets (right). Red dots represent positive droplets. Blue dots represent negative droplets. . . . . . . . . . . . . . . . . . . . . . . . . 79 5.11 The y-z phase space plot (left) and the x-y phase space plot (right) with separation x=d 0 = 0:4 of all droplets. Red dots represent positive droplets. Blue dots represent negative droplets. . . . . . . . . . . . . . . . . . . . . . 79 5.12 The space charge potential of the beams with separation x=d 0 = 0:4. . . . . 80 5.13 The average velocity of v x along the z direction with separation x=d 0 = 0:4. 80 5.14 Thex-z phase space plot with separationx=d 0 = 0:1 of all droplets. Red dots represent positive droplets. Blue dots represent negative droplets. . . . . . . 81 5.15 The space charge potential of the beams with separation x=d 0 = 0:1. . . . . 81 5.16 The average velocity of v x along the z direction with separation x=d 0 = 0:1. 82 5.17 The x-z phase space plot with separation x=d 0 = 0:4 at time step 31000 and volume ow rate 10 _ V . Red dots represent positive droplets. Blue dots represent negative droplets. . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.18 The x-z phase space plot with separation x=d 0 = 0:4 at time step 27000 and volume ow rate 100 _ V . Red dots represent positive droplets. Blue dots represent negative droplets. . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.19 The space charge potential of the beams with separation x=d 0 = 0:4 and volume ow rate 10 _ V . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.20 The space charge potential of the beams with separation x=d 0 = 0:4 and volume ow rate 100 _ V . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.21 The average velocity of v x along the z direction with separation x=d 0 = 0:4 and volume ow rate 10 _ V . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.22 The average velocity of v x along the z direction with separation x=d 0 = 0:4 and volume ow rate 100 _ V . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.23 The x-z phase space plot with separation x=d 0 = 0:4 of all droplets at time step 26000. From the largest beam to the smallest beam, 0 changes from 2000 V, 3000 V, to 4000 V. . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.24 Left: The space charge potential of the beams with separationx=d 0 = 0:4 and 0 = 2000 V at time step 26000. Right: The average velocity of v x . . . . . . 86 5.25 Left: The space charge potential of the beams with separationx=d 0 = 0:4 and 0 = 3000 V at time step 26000. Right: The average velocity of v x . . . . . . 86 LIST OF FIGURES xv 5.26 Left: The space charge potential of the beams with separationx=d 0 = 0:4 and 0 = 4000 V at time step 26000. Right: The average velocity of v x . . . . . . 87 5.27 The x-z phase space plot with separation x=d 0 = 0:4 of all droplets at time step 26000 with d 0 = 2 mm, d 0 = 4 mm, and d 0 = 8 mm, respectively. . . . 87 5.28 Left: The space charge potential of the beams with separationx=d 0 = 0:4 and d 0 = 2 mm at time step 26000. Right: The average velocity of v x . . . . . . . 88 5.29 Left: The space charge potential of the beams with separationx=d 0 = 0:4 and d 0 = 4 mm at time step 26000. Right: The average velocity of v x . . . . . . . 88 5.30 Left: The space charge potential of the beams with separationx=d 0 = 0:4 and d 0 = 8 mm at time step 26000. Right: The average velocity of v x . . . . . . . 89 6.1 Left: The analytical solution of the potential (in volts). Right: The absolute dierence between the analytical potential and the simulated potential. . . . 93 6.2 Error over iteration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 6.3 Error over iteration for nonuniform meshes. . . . . . . . . . . . . . . . . . . 95 6.4 Left: Simulated potential. Right: Error over iteration. . . . . . . . . . . . . 97 6.5 Left: Comparison between simulated E z (labeled as S) and analytical E z (labeled as A) in the unit of V/m. Right: The absolute dierence. . . . . . 97 6.6 The absolute dierence of dierent cases. . . . . . . . . . . . . . . . . . . . 98 6.7 PP-PM hybrid test setup. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 6.8 Charge density distribution. The colorbar shows log 10 of volts. Electric po- tential distribution in volts. The case with uniform mesh is on the left, while the linear one is on the right. . . . . . . . . . . . . . . . . . . . . . . . . . . 99 6.9 Electric potential distribution in volts. The case with uniform mesh is on the left, while the linear one is on the right. . . . . . . . . . . . . . . . . . . . . 100 6.10 log 10 E r over z of the pure PP model. . . . . . . . . . . . . . . . . . . . . . 100 6.11 Relative errors of E r between the pure PP case and hybrid cases. . . . . . . 101 7.1 Simulation domain and setup. . . . . . . . . . . . . . . . . . . . . . . . . . 103 7.2 The accelerating potential distribution in volts with 500 500 mesh grids, and uniform mesh on the left, linear mesh on the right. . . . . . . . . . . . 104 7.3 E z of particles over z for dierent cases. . . . . . . . . . . . . . . . . . . . . 105 xvi LIST OF FIGURES 7.4 v z of particles over z for dierent cases. . . . . . . . . . . . . . . . . . . . . 105 7.5 Number density plots of case S1 (top left panel), S2 (top right panel), S3 (bottom left panel), and S4 (bottom right panel), in the unit of log 10 m 3 . . 106 7.6 Number density plots of case S3 with accelerating potential 1250 V (left panel) and 1750 V (right panel), in the unit of log 10 m 3 . . . . . . . . . . . . . . . 107 7.7 Phase space z-v z plots of case S1 (top left panel), S2 (top right panel), S3 (bottom left panel), and S4 (bottom right panel). . . . . . . . . . . . . . . . 108 7.8 Phase spacez-v z plots of case S3 with accelerating potentials 1750 V, 1500 V, and 1250 V from the top curve to the bottom curve, respectively. . . . . . . 109 7.9 Cuto applied in the PP part of the hybrid model. . . . . . . . . . . . . . . 110 7.10 Number density plots of case S3 with cutos in the unit of log 10 m 3 . . . . 111 7.11 Relation between r c and the expansion angle. . . . . . . . . . . . . . . . . . 111 7.12 Local view of the beam near the injection area. . . . . . . . . . . . . . . . . 112 List of Tables 2.1 Material parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Practical parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.1 Time step length verication. . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.1 Chosen material parameters (formamide). . . . . . . . . . . . . . . . . . . . 37 4.2 Chosen practical parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.3 Derived parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.4 Simulation cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.5 Percent error (%) of system energy of each case. . . . . . . . . . . . . . . . 44 4.6 Regions for evaluating f v and f p of all cases. . . . . . . . . . . . . . . . . . 47 6.1 Comparison of the elapsed simulation CPU time. . . . . . . . . . . . . . . . 102 xvii xviii LIST OF SYMBOLS List of Symbols " 0 Vacuum permittivity F T Thrust I sp Specic impulse 0 Accelerating potential d Distance between the Taylor cone tip and the accelerating plate d 0 Distance between the jet front and the accelerating plate K Conductivity " Relative permittivity Surface tension Mass density f(") Empirical parameter for calculating current : V Volume ow rate R n Needle radius I Current R d Radius of original droplet R j Radius of jet r Characteristic dimension of cone-jet transition region V d Volume of original droplet : N d Emission rate of original droplet Q d Charge of original droplet M d Mass of original droplet LIST OF SYMBOLS xix v d Initial axial speed of original droplet t Simulation time step length H Height of injection volume N in Number of injected original droplets per injection F pp Force due to other droplets E pp Electric eld due to other droplets F acc Force due to accelerating potential E acc Electric eld due to accelerating potential acc Accelerating potential a Half focal length of hyperboloids in prolate spheroidal coordinate system 0 A coordinate indicating liquid surface in prolate spheroidal coordinate system R c Radius of curvature T Taylor cone angle N s Number of droplets broken from original droplet E tot Average total energy of each droplet E p Potential energy E k Kinetic energy e Equivalent length f e ( e ) Probability density function of equivalent length k B Boltzmann constant T Macro-temperature of droplets g 0 Gravitational acceleration at the ground level r m Mean interparticle distance v th Thermal velocity t p Period of the AC mode % Charge density r c Cuto interparticle distance xx NOMENCLATURE Nomenclature FEEP Field emission electric propulsion ILIS Ionic liquid ion sources PM The Particle-Mesh model PIC The Particle-in-Cell method PP The Particle-Particle model PDF Probability density function DoD Drop-on-demand electrospray PP-PM The Particle-Particle Particle-Mesh hybrid model Chapter 1 1 Chapter 1 Introduction 1.1 Background Colloid thrusters are a subset of electrospray propulsion devices, which electrostatically accelerate charged particles that are produced from electried liquid surfaces to generate thrust. The idea of utilizing electrospray emissions for space propulsion can be traced back to the early 1960's from the work of Krohn [45, 41], in which liquid metals and viscous organic liquids (such as glycerol) were rst considered. Electrospray propulsion was then intensively studied till around 1975. The development of electrospray thrusters in those days mainly focused on applications that require relatively large thrust densities, because the emitted particles in electrospray thrusters usually have larger mass-to-charge ratio m=q than that of other electric propulsion devices (such as ion thrusters): a largerm=q corresponds to a larger accelerating voltage according to =mc 2 =(2q), wherec denotes the nal ow speed, and an increasing results in an increasing thrust density according to F T =A = 8" 0 2 =(9d 2 ), where d denotes the grid spacing and " 0 is the vacuum permittivity. Also, higher increases the eciency=(+ loss ), since the voltage used for producing charged particles, loss , becomes less signicant [51]. However, since the generated droplet mass-to-charge ratio, applying the technology back then, was so large that the corresponding accelerating voltage ranges from 10 to 100 kV (for a typical specic impulse I sp 1000s), which is too high and results in dicult insulation and packaging problems. In addition, the desired high thrust density requires arrays of a large number of individual liquid-dispensing needles (each providing a thrust of the order of 1 N), which leads to high complexity of the design and manufacture, thus further makes electrospray thrusters unattractive. Because of the above, the development of electrospray propulsion was stuck for nearly three decades. Fortunately, electrospray thrusters start to draw attention again recently due to the facts that [51, 49, 15]: (a) The demand on miniaturization of spacecraft. It is now desirable to have small thrust per emitter, which allows designs with precise controllability and high 2 Chapter 1 performance; (b) The advances of electrospray technology made in other elds, such as biological science, in which electrospray is used to extract charged biological macro-molecules from liquid samples. These advances would overcome previous limitations on the specic charge q=m of droplets, thus reduce the accelerating voltage to a more comfortable range (1 to 5 kV); (c) The advances of micro-manufacturing technologies now allow for ecient clustering of a large number of emitter tips on a small surface. Electrospray thrusters can be divided into three basic types: (a) Colloid thrusters, which apply solvents such as doped glycerol and formamide as propellant to produce and acceler- ate charged droplets (and also ions under some circumstances); (b) Field emission electric propulsion (FEEP), which uses liquid metals such as Cs and In to produce and accelerate charged metallic ions; (c) Ionic liquid ion sources (ILIS), which apply room-temperature molten salts to produce salt ion beams, or mixtures of ions and droplets. The advantages and disadvantages of each type depend on dierent missions [61], but colloid thrusters dif- fer from FEEP thrusters in a way that colloid thrusters can eject charged droplets, ions or a mixture of both. This additional exibility makes colloid thrusters more attractive for missions requiring either high performance (high ion fractions with high I sp and current) or high resolution in thrust determination (droplets only with low I sp and current). Therefore, rather than the other two types of electrospray thrusters, colloid thrusters are chosen to be the main focus in this work. 1.2 Basic Physics 1.2.1 Starting Voltage When a strong normal electric eld is applied on a at liquid surface, if the liquid contains free ions, those of the attracted polarity will concentrate on the surface [51, 49]. If the liquid surface deforms, the electric eld becomes stronger on the protruding parts, and more charge concentrates there. Then, a Taylor cone forms with a cone angle T = 49:29 [70, 81, 26]. In order for the liquid to overcome the surface tension forces and start owing, one needs to have 1 2 " 0 E 2 tip > 2 R c (1.1) where " 0 is the vacuum permittivity, is the surface tension, R c is the radius of curvature, and E tip is the electric eld at the Taylor cone tip. E tip = 2=R c ln(4d=R c ) : (1.2) This leads to a \starting voltage" [51, 75, 57] start = r R c " 0 ln 4d R c (1.3) Chapter 1 3 where d is the tip-to-plane distance. 1.2.2 Current and Charge The current transported by cone-jet electrosprays veried experimentally is [26, 37, 6] I =f(")( K _ V " ) 1=2 ; (1.4) where K is the conductivity, _ V is the volume ow rate, " is the relative permittivity, and f(") is an empirical parameter and has a value of about 18 for " > 40 [25]. Note that I is independent of the accelerating voltage, the electrode shape, and the liquid viscosity [14, 31, 36, 59]. The jet radius R j and the droplet radius R d emitted from the jet have an relation, and can be calculated from [35, 64, 24] R j = R d 1:89 = r 1:89 [ 3 f(") ] 2=3 ; (1.5) where r = ("" 0 _ V=K) 1=3 is the \characteristic dimension" of the cone-jet transition region [51, 66]. Note that the factor 3 is used here instead of 6 in Eq. (1.5), because experiments have shown that electrosprays produce streams of droplets charged to about half of their Rayleigh limit [51]. The Rayleigh limit indicates the maximum charge q R a droplet can hold, q R = 8 p " 0 R 3=2 d (1.6) where R d denotes the droplet radius. Based on the Rayleigh limit, and if we assume the droplet has a spherical shape, its mass is m = 4R 3 d =3, so the maximum specic charge carried by droplets would be q m max = 6 p " 0 R 3=2 d : (1.7) The droplet specic charge can be estimated from [51, 38, 58] q m = r K " _ V f(") ; (1.8) where _ V is the volume ow rate, is the mass density of the liquid. 4 Chapter 1 1.2.3 Limitations of Droplet Specic Charge As we have seen, high specic charge q=m can be obtained by increasing the liquid conduc- tivityK or decreasing the volume ow rate _ V . As _ V=K is reduced, the jet becomes thinner based on the relation r / ( _ V=K) 1=3 , as well as the droplet size, and the specic charge increases as ( K= _ V ) 1=2 . However, there are two phenomena that limit the increase of q=m. (a) Taylor cone instability Experiments have shown that the Taylor cone becomes unstable when = (K _ V= "" 0 ) 1=2 is less than some lower limit (e.g. of the order of 0.5 in conventional capillary tubes) [51, 73, 32]. The physics of this instability has not been fully understood, but one ex- planation is that the limit is the space charge that would lead to full separation of the positive and negative ions of the liquid used, q V max = q m= max =F 1000c d (C/m 3 ); (1.9) where F = 96500 C/mol is Faraday's constant, c d is the dissociated part of the solution's equivalent normality (mol/l). c d is linearly related toK through a \mobility parameter" 0 , c d =K= 0 . From Eq. 1.8, we can obtain f(") K " _ V 1=2 = 1000F 0 K , min = r " 0 0 1000F f(") " (1.10) Because 0 would be expected to depend on viscosity, so the argument is incomplete. One can use _ V min = "" 0 K 2 min (1.11) and Eq. 1.8, to obtain q m max = f(") " min K p " 0 (1.12) which reduces back to Eq. 1.9 using Eq. 1.10. Therefore, we can use Eq. 1.12 to estimate the maximum specic charge, if the Taylor cone instability is the main consideration. (b) Ion emission The normal electric eld increases towards the tip of the cone. We can estimate the maximum eld at the tip, E = 1:87 10 7 [f(")] 1=3 1=2 K " _ V 1=6 : (1.13) Chapter 1 5 If we evaluate it at the lowest stable ow rate, _ V min = "" 0 2 min K ; (1.14) then the maximum eld becomes E max = 1:3 10 9 1=6 1=3 min f(") K " 1=3 : (1.15) It is known from experiments that at normal elds in the range 1 - 2 V/nm, individual ions begin to be extracted from the liquid by eld evaporation [49, 51, 28]. 1.2.4 Droplet Breakup For the external electric eld E ex acting on a droplet, which could be due to the applied accelerating eld or other charged droplets, if it is over the surface tension of the droplet, the droplet could break up. We can thus estimate a maximum external electric eld that could result in the droplet breakup from 1 2 " 0 E 2 max = 2 R d ) E max = r 4 " 0 R d : (1.16) Therefore, besides the Rayleigh limit, this E max could be used as an criteria to break up droplets. Also, the internal electric eld E in can be added with the external electric eld, E in = q 4" 0 R 2 d ; (1.17) to judge if a droplet should break up. Namely, ifE in +E ex >E max , the droplet should break up. 1.3 Motivations and Objectives Due to the complex electrohydrodynamics, the multiscale features, and the lack of precise experimental data, the detailed physics of colloid thrusters is far from being fully under- stood. Thus, developing numerical models and implementing accurate simulations for colloid thrusters are necessary, which can provide insights to reveal more unknown physics. How- ever, there are currently rare simulation works on colloid thrusters, especially for the beam region. Therefore, one motivation of this research is to develop a fully kinetic particle model to simulate the colloid thruster beam, which is required to possess a high accuracy such that a complete and precise description of the droplet dynamics can be obtained. In particle simulations, Particle-Mesh (PM) models [39] have been widely used due to 6 Chapter 1 their relatively high accuracy and eciency. One well developed PM model for simulating electric propulsion devices (such as ion thrusters [40, 78], Hall thrusters [48], and some other types of thrusters [67, 88]) is the Particle-in-Cell (PIC) model [7], which applies a mesh system and the macro-particle assumption. If one attempts to simulate the beam of a colloid thruster using PIC, the required mesh resolution would change signicantly from the region near the jet to the downstream region of the beam, because the scale of the jet is much smaller than the expanded beam size. Also, since the colloid thruster beam can contain big droplets instead of small ions, the validity of using the macro-particle assumption on droplets is unknown. Therefore, the diculties in applying particle-mesh models motivate us to seek for another approach. The approach found to be suitable for the purpose of this research is the Particle-Particle (PP) model [39, 85, 89] (using real particles rather than macro-particles). The PP model does not have a mesh system, all interparticle forces are calculated using particle positions and the considered force law. Thus, applying the PP model can completely discard the mesh resolution issue discussed before. Knowing the disadvantage of the PP model very well, namely the large computation which scales as the square of the number of particles, its advantages tend to be ignored. For example, if we limit ourselves in classical electro- statics, which is believed to be the dominant physics in the beam of colloid thrusters, some advantages of the PP model are as follows. First, the PP model directly applies Coulomb's law to calculate all the interparticle forces, thus if the simulation time step length is small enough, it can describe the precisest particle movements. Second, the PP model contains the least numerical errors compared to the PM model, which includes errors (such as truncation errors and interpolation errors) caused by solving Poisson's equation on mesh points and interpolations between particle positions and mesh points. Third, because the simulation algorithm of the PP model is simple and straightforward, when using shared-memory paral- lel computing techniques, the parallel eciency can be extremely high. Especially with the development of GPU computing [1], the PP model is becoming more and more practicable than it used to be [89, 84, 86, 87, 85]. Another motivation of this work is the fact that the PP model can resolve Coulomb collisions with the highest accuracy, and the Coulomb collisions between charged droplets occurred in the beam of colloid thrusters are intrinsically important but have not been well studied. For PIC or other PM models, the interparticle forces within the range of a mesh are either weaken non-physically or ignored, thus these models cannot resolve Coulomb collisions unless the mesh size is chosen to be suciently small, which however would make the PM model even more computational expensive than the corresponding PP model. In addition, the results of the PP model, especially with respect to Coulomb collisions, can be utilized as a baseline for developing, validating, and optimizing other models. For example, the PP model could give a relation between the mesh size of a PM model and the amount or extent of Coulomb collisions that model can capture. The objectives of this work are therefore as follows: (a) Develop a PP model that can simulate the colloid thruster beam with high accuracy. (b) Apply the PP model to obtain a complete and precise description of the droplet dynamics and reveal more unknown physics. Chapter 1 7 (c) Study in detail the Coulomb collisions between droplets and its variations in dierent circumstances as well as its dependence on other parameters. (d) Analyze results of the PP model to give criteria and guidance for developing other models, such as a PP-PM hybrid model. 1.4 Literature Review Colloid thrusters were rst implemented in space through the Space Technology 7 Mission [90], which launched on December 3, 2015. The Space Technology 7 Mission aimed to detect and measure gravitational waves [5]. It tested a Disturbance Reduction System that was planned for launch on the LISA Pathnder mission in 2013 [92, 29, 20, 21]. The colloid thrusters used in the Space Technology 7 Mission were two clusters of four colloid micro- newton thruster systems. Each of eight thruster systems can provide thrust between 5 and 30 N with resolution0.1 N and thrust noise0.1 N/ p Hz. In 2017, the colloid microthruster ight performance results from the space technology 7 disturbance reduction system was reported [91, 27]. The system performed as expected and required, meeting all the requirements based on the obtained on-orbit data and analysis. The microthrusters operated for more than 2400 hours in ight during commissioning activities, a 90-day experiment and the extended mission. 1.4.1 Relevant Computational Works Overall, simulations of colloid thruster beam are very rare. Here, we discuss several previous works that are relevant to the focus of our research. Rather than colloid thrusters, M. Tajmar and J. Wang [68] in 2000 applied 3-D PIC simulations on a FEEP thruster, and focused on the topic of beam neutralization. They found that to neutralize the beam, one must either operate the thruster in an ambient plasma with a density comparable with that of the beam, or place a neutralizer at a location suciently far from the thruster emitter. Later in 2001, M. Tajmar and J. Wang [69] studied the back ow contamination of FEEP thruster using the previous PIC code plus a Monte Carlo collision model. They found that the back ow ion current is on the order of 0.01% of the total emitting current, and operating at lower neutral ux, lower emitter current, or higher accelerating potentials helps to reduce the back ow contamination. The group of Dr. Deborah A. Levin at Pennsylvania State University applied PIC to simulate colloid thruster beam in recent years. The following are three works done by this group. In 2012, Kumar et al. [43] applied an electrostatic 3-D PIC model without the consideration of particle collisions to simulate a colloid thruster beam with EMIM-BF 4 [23] as the propellant. They simulated dierent cases using ions of EMIM + or BF 4 rather than large droplets, and used results acquired from previous molecular dynamics simulations of Taylor cone [10] as initial particle conditions. They obtained the particle velocity distributions, 8 Chapter 1 thrust, specic impulse, etc, and found that the beam under their simulation remains fairly collimated. Also in 2012, Wang et al. [79] parallelized the above PIC code using domain decomposition such that larger cases can be simulated. They presented some new results such as distributions of the number density of the beam. In 2013, Wang et al. [80] further improved their PIC code with a multi-level grid, such that the computation can be reduced, because larger meshes are applied for the downstream region that is far away from the injection source, while smaller meshes are applied near the injection source. They considered both ions and droplets, and found the simulation results agree with experimental data. One remark is that, when using a ner mesh near the injection source, the beam expands more radially. It is because using a ner mesh includes more short-range interparticle forces, i.e. Coulomb collisions, as discussed in Sec. 1.3. In 2015, Borner et al [9] did a study using both molecular dynamics and a 3-D Poisson solver. They focused on the importance of parameters relative to the grid used to solve Poisson's equation, such as the grid cell size and the boundary conditions. It was found that the boundary conditions have an important impact on the potential and electric eld. In recent years, the group of Dr. Deborah A. Levin also did some research on dierent topics of electrospray via molecular dynamics simulations. In 2017, they did molecular dynamics electrospray simulations of coarse-grained Ethylammonium Nitrate (EAN) and 1-Ethyl-3-Methylimidazolium Tetra uoroborate (EMIM-BF 4 ), and found that the emission modes depend on the electro-chemical properties of the ionic liquid propellant [52]. For example, EMIM-BF 4 -based colloid thrusters operate in ion mode and EAN-based devices operate in the droplet mode. Also in 2017, they compared the behaviors of two protic ionic liquid in the presence of an electric eld using molecular dynamics [53]. In 2018, they studied the sensitivity of electrospray molecular dynamics simulations to long-range Coulomb interaction models [54], in which three models are compared: the direct Coulomb (DC) model, the shifted force Coulomb sum (SFCS) model, and the particle-particle particle-mesh (PPPM) long-range Coulomb model. In 2019, they developed a new octree-based Coulomb interaction model to simulate the electrospray of ionic liquids (ILs) [55]. They compared this model with computationally more expensive models, such as the direct Coulomb model. It is found that the octree-based method is capable of reproducing Coulomb energy in agreement with the direct Coulomb model. There are also some other works they have done [13, 12, 11] In addition, recently M. Rahmanpour and R. Ebrahimi [60] applied the idea of a PP model to simulate the droplet breakup in a 2-D electrospray beam. Their model included the PP forces, the accelerating eld forces, and the air drag forces, as well as the droplet evaporation. The simulation results of droplet size and velocity distribution are obtained, and found to be in good agreement with other theoretical and experimental studies. As can be seen, none of the above research studied the detail of Coulomb collisions occurred in colloid thruster beam, and when PIC is used, the eects caused by using dierent cell sizes have not been well studied. In this work, we will thus use the PP model to research on and attempt to answer these unsolved questions. Chapter 1 9 1.4.2 Relevant Charge Neutralization Works The charge neutralization of colloid thrusters is one topic that will be studied in the dis- sertation, thus here we review some relevant literatures. Electrospray thrusters do not emit electrons inherently, but usually are operated at only positive polarity with a neutralizing electron source placed outside of the thruster [22, 65, 83, 90, 27]. Since electrospray thrusters can emit either positive or negative ions of similar mass, resulting in either a positive or neg- ative ion beam, we can operate one thruster in positive polarity and another in negative polarity to let the total emitted currents be zero, i.e. a bipolar system, thus the spacecraft should not charge [56]. Without the charge neutralization, a potential dierence between the spacecraft and its environment could build up quickly, causing damaging arcs or discharges, disruption of the emitted beam, corruption of scientic measurements, etc [44, 34]. Mier-Hicks and Lozano at MIT in 2017 studied the bipolar mode of electrospray thrusters [56]. They developed theoretical and experimental methods to investigate the charging characteristics expected to be observed on spacecraft during the operation of electrospray thrusters. Their experimental results demonstrate that the neutralization from the bipolar mode is indeed possible. The thrusters are able to be operated for long periods of time, inducing bounded spacecraft charging in the range from -400 to +600 V, when emitting cur- rents of about 20A. The thruster arrays they used are housed in a 13 12 2.4 mm silicon frame. 480 laser-ablated tips are held in a porous glass emitter substrate, which is contained in the silicon frame. A gold-coated silicon extractor grid is above these tips. The thruster array produces thrust about 0.1 N per microampere of emitted current. They operated the arrays with current around 150 A and thrust 15 N per thruster. The propellant used was EMI-BF4 [42]. When operating, the polarity of each thruster is periodically alternated about every 30 s to avoid electrochemical reactions that could degrade the emitter tips [50]. Another study of neutralization from applying bipolar pairs of electrospray thrusters was done by Courtney et al. again in 2017 [22]. Each thruster has a porous boroslicate substrate with a diameter 1 cm, thickness 3 mm, and nine 7.5 mm long triangular prisms. Each prism has a sharp edge with a radius of a few 10's m. An extracting grid is roughly 100 m above the prisms. A accelerating potential of about 2 kV is applied, facilitating ion beams of several hundredA, and is alternated under a 1Hz square wave. Again, the propellant EMI- BF4 was used. They evaluated experimentally an active current-balance control circuit and a passive self-balancing conguration. The passive conguration produced an unbalanced charge collection, indicating the unsuitability of that architecture. The active conguration produced eective charge neutralization, during stable periods of emission, but insucient neutralization during potential alternations. And they provided methods to improve the congurations for neutralization. 10 Chapter 1 1.5 Organization The dissertation considers two dierent types of beams: one is droplet only, and the other is ion only. A beam of a mixture of droplets and ions are not considered. The detail of the basic frame of the PP simulation model is given in chapter 2, including the simulation domain setup, input and derived parameters, droplet initialization, equations of motion, and some important diagnoses. The model assumes that the beam consists of only droplets and the droplets in dierent simulations are all of the same size. In chapter 3, the PP model is veried via a time step length verication using the binary collision theory, and via a space charge eect study using Child-Langmuir law, as well as by comparing with a 3D electrostatic PIC model. In chapter 4, reasonable parameters are chosen and characteristic simulations are ex- ecuted, in which the time step length and the number of broken droplets are varied to investigate the PP model and the eects of number of droplets. The simulation results are analyzed and shown, such as the time history of system energy, droplet phase space plots, droplet number density, macro-temperature (dened in Sec. 2.7.3), velocity distribution function, and the probability density function of equivalent length (dened in Sec. 2.7.2). A corresponding PM model is discussed based on the PP results, and the initial droplet thermal velocity is also studied for completeness. The space charge eects of the colloid thrusters under the parameter range simulated are studied. It is found that the Child-Langmuir law is inappropriate to be applied for the scenario of colloid thrusters, and possible reasons are given and discussed. In chapter 5, the AC mode of a colloid thruster is studied. Two types of alternation are considered. The rst type is a square wave which helps to study the droplet dynamics when both positively and negatively charged droplets are involved. The second type is a step function which helps to study the dependence of the transition time, the time that the thruster changes back to a steady state after one alternation. In addition, a bipolar operation mode is studied, in which two thrusters with dierent polarity are placed side by side to avoid the charge accumulation on the spacecraft. Space charge eects are studied by varying the separation between the two thrusters, the accelerating potential, the tip-plane distance, and the volume ow rate. In chapter 6, a PP-PM hybrid model is developed, which applies 3-D PP for the upstream beam region and 2-D cylindrical PM with nonuniform mesh for the downstream beam region, so that the simulation can be faster while still maintain the high accuracy. A test simulation indicates the validity of this hybrid model and veries that the PP model can provide the suciently minimum mesh size for a corresponding PM model. In chapter 7, the proposed PP-PM hybrid model is applied to simulate some cases studied in the reference [80]. The results of the hybrid model are compared with those in the reference, including the number density contours, particle exit speed plots, etc. Chapter 1 11 At the end in chapter 8, a summary of the PP model and PP-PM hybrid model is given, and some future works are discussed. 12 Chapter 2 Chapter 2 PP Simulation Model The detail of the PP model is given in this chapter, including the simulation domain, in- put and derived parameters, initialization of droplets, equations of motion, droplet initial breakup, simulation procedure, and diagnoses. The model assumes that the beam consists of only droplets and the droplets in dierent simulations are all of the same size. 2.1 Simulation Domain Figure 2.1: Schematic diagram of a colloid thruster and the simulation domain. A schematic diagram of a colloid thruster is shown in Figure 2.1, where the needle/capillary and the accelerating plate/grid (shown in gray) are connected by a power supply with a potential dierence 0 , and there is liquid in the needle (shown in blue). Due to the potential Chapter 2 13 dierence, i.e. the electric eld acting on the liquid surface, a Taylor cone is formed at the end of the needle, and a jet is emitted from the tip of the cone. At some point, the jet breaks up into droplets, and large droplets further break up into smaller ones, all of which are accelerated toward the plate. The simulation domain is chosen to be from the location where the jet breaks to the accelerating plate, as the green shaded region ABDC shown in Figure 2.1. The origin of the simulation domain is chosen to be the point O labeled in the gure; the z-axis is perpendicular to the accelerating plate pointing into the needle along the symmetric axis of the needle; the x-axis is along the plate; and the y-axis points into the paper (not shown in the gure). Therefore the simulation domain is a three dimensional cube containing all the droplets. d is used to denote the distance from the tip of the Taylor cone to the accelerating plate, and d 0 is used to denote the distance from the location where the jet breaks to the accelerating plate. However, one should note that since the PP model does not have a mesh system, the simulation domain given above is just for the purpose of indicating the main simulation region. In the simulation, droplets are allowed to move radially to be anywhere far away from the z-axis, namely there are no boundaries of the simulation domain in the x and y directions. For the z direction, we will choose a simulation time such that the beam front arrives at the location we desire. 2.2 Input and Derived Parameters Parameters for the simulation model are divided into three categories: the parameters of the liquid material (labeled as material parameters), the parameters of the setup and operation of the colloid thruster (labeled as practical parameters), and parameters of the simulation algorithm (labeled as numerical parameters). The material parameters are listed in Tab. 2.1, the practical parameters are listed in Tab. 2.2, while the numerical parameters will be given later when needed. Parameter name Symbol SI unit Conductivity K S m 1 (or kg 1 m 3 s C 2 ) Relative permittivity " Dimensionless Surface tension N m 1 (or kg s 2 ) Mass density kg m 3 Empirical current factor f(") Dimensionless Table 2.1: Material parameters. 14 Chapter 2 Parameter name Symbol SI unit Volume ow rate _ V m 3 s 1 Distance between the liquid tip and the plate d m Distance between the jet front and the plate d 0 m The accelerating voltage 0 V (or kg m 2 C 1 s 2 ) Needle radius R n m Table 2.2: Practical parameters. If the input parameters listed in Tab. 2.1 and Tab. 2.2 are specied, all other required parameters in the simulation can be obtained analytically: (a) The ow current I can be obtained from (shown in Eq. 1.4) I =f(")( K _ V " ) 1=2 : (2.1) (b) The radius of the original dropletR d and the radius of the jetR j can be calculated from (shown in Eq. 1.5) R j = R d 1:89 = r 1:89 [ 3 f(") ] 2=3 : (2.2) (c) The shape of all droplets are assumed to be spherical, thus the volume of the original droplet is V d = 4R 3 d =3: (2.3) (d) The time rate of the number of the original droplets emitted from the jet is _ N d = _ V=V d : (2.4) (e) The charge of the original droplet is Q d =I= _ N d : (2.5) (f) The mass of the original droplet is M d =V d : (2.6) (g) If we assume the original droplet has the same speed as that of the jet, the initial axial droplet speed v d is (along the negative z-axis direction) v d = _ V R 2 j : (2.7) Chapter 2 15 2.3 Initialization of droplets The simulation starts when the jet begins to break up, thus the applied initialization is to place the appropriate number of original droplets in a corresponding volume at the jet front. If the simulation time step length t is specied, an injection volume can be determined. Figure 2.2 shows the injection volume, which is a right cylinder with radius R j and a height H. The point A is located at (0; 0;d 0 ) in the simulation domain, and the symmetric axis of the cylinder is the z-axis. The number of the original droplets that should be injected is N in = _ N d t: (2.8) If N in 1, N in droplets will be injected at each time step. If N in < 1, injection may not happen every time step. The height of the injection volume is H = v d t. Of course, it is convenient if we choose certain values of t such that N in are integers or unit fractions. Figure 2.2: Injection volume. The positions of the injected droplets are set to be random in the injection volume, namely the initial coordinates of each droplet are obtained from x =R j p U sin (2U 0 ); y =R j p U cos (2U 0 ); z =d 0 U 00 H; (2.9) where U, U 0 , and U 00 are three independent random numbers uniformly distributed in the interval (0; 1). For the initial velocity of each injected droplet, it is chosen to be v =v d ^ z = (0; 0;v d ); (2.10) where ^ z is the unit vector along z-axis, namely only the axial initial velocity is considered, the radial initial velocity is ignored. 16 Chapter 2 2.4 Equations of Motion For any charged droplet, including the original droplets as well as those smaller droplets broken from the original droplets (through Coulomb ssion [47] in the reality), the equations of motion can be described as m d 2 r dt 2 =m dv dt =F =F pp +F acc =q(E pp +E acc ); (2.11) where m denotes the mass, q denotes the charge, r denotes the position vector, v denotes the velocity vector, and F denotes the resultant force. Two types of forces are considered: F pp is the electrostatic force acting on this droplet due to all other droplets, E pp is the corresponding electrostatic eld;F acc is the electric force caused by the applied accelerating potential,E acc is the corresponding accelerating electric eld. 2.4.1 Particle-Particle Field The electric eld at the location of the ith droplet caused by all other droplets can be calculated using Eq. (2.12), where r ij is the distance vector pointing from the jth droplet to the ith droplet, N denotes the total number of simulated droplets at that time. E (i) pp = 1 4" 0 N X j=1;j6=i q j jr ij j 3 r ij (2.12) One should note that there are two types of systems with respect to charge. For systems with unlike charges, i.e. both positive and negative charges are involved, since positive charges attract negative ones,jr ij j could approach 0 physically, which results in an innitely large eclectic eld if the recombination of these two charges is ignored in the simulation, unless a cuto is introduced [84, 86, 87]. For systems with like charges, i.e. charges are either all positive or all negative,jr ij j should not approach 0 physically, since the interparticle forces are always repulsive. In addition, the time step length used to numerically integrate Eq. (2.11) should always be small enough. Because if the time step length is too large, one droplet could suddenly move to a position that is very close to another droplet after one integration, which causes nonphysical large forces in the next time step, no matter what charges these two droplets have. 2.4.2 Accelerating Electric Field The accelerating electric eld is given using an approximate analytical solution [51] described as follows. The prolate spheroidal coordinates (;;') are used for convenience, as shown Chapter 2 17 in Figure 2.3, = (r 1 +r 2 )=a; = (r 1 r 2 )=a; (2.13) r 1 = [x 2 +y 2 + (z +a=2) 2 ] 1=2 ; r 2 = [x 2 +y 2 + (za=2) 2 ] 1=2 ; (2.14) where a is the half focal length of the hyperboloids, and ' is the azimuthal angle about the lineFF 0 , which is not shown in the gure, as well as they-axis which points into the paper. Lines of constant are confocal ellipsoids with the same foci (F and F 0 ), while lines of constant are confocal hyperboloids. In the simulation model, the surface = 0 represents the accelerating plate, the surface = 0 represents the liquid surface, i.e. the Taylor cone, thus the vertical distance between the liquid surface tip and the accelerating plate is d. Figure 2.3: The prolate spheroidal coordinate system. Since the potential acc on the liquid surface ( = 0 ) can be assumed to be a constant 0 , and zero on the accelerating plate ( = 0), the solution for acc in the whole space will only depend on . The part of Laplace's equation is then @ @ [(1 2 ) @ acc @ ] = 0; (2.15) which can be integrated easily using R (1 2 ) 1 d = [ln (1 +) ln (1)]=2 = tanh 1 with the boundary conditions at = 0 and = 0 to obtain acc () = 0 tanh 1 tanh 1 0 : (2.16) Knowing the potential, the electric eld can be solved from E acc =r acc = ( @ acc @ @ @x ; @ acc @ @ @y ; @ acc @ @ @z ) (2.17) 18 Chapter 2 where @ acc @ = 0 tanh 1 0 1 1 2 ; (2.18) @ @x = x a ( 1 r 1 1 r 2 ); @ @y = y a ( 1 r 1 1 r 2 ); @ @z = 1 a ( z +a=2 r 1 za=2 r 2 ): (2.19) Now the values ofa and 0 need to be determined. Letr xy denotes the cylindrical radius, i.e. r 2 xy =x 2 +y 2 , from = (r 1 r 2 )=a, the (r xy ,z) relation for a constant hyperboloid can be described as a = [r 2 xy + (z +a=2) 2 ] 1=2 [r 2 xy + (za=2) 2 ] 1=2 ; (2.20) which can be rearranged to (for z > 0) z =( a 2 4 + r 2 xy 1 2 ) 1=2 ; (2.21) by squaring the equation rst and then moving the terms with square roots to one side and squaring the equation again. The radius of curvature of this surface can thus be calculated, R c = (1 +z 02 ) 3=2 z 00 ; (2.22) wherez 0 andz 00 denote the rst and the second derivative ofz with respect tor xy , respectively, based on Eq. (2.21). Substituting z 0 = dz dr xy = 1 2 ( a 2 4r 2 xy + 1 1 2 ) 1=2 ; z 00 = d 2 z dr 2 xy = 1 2 a 2 4 r 3 xy ( a 2 4r 2 xy + 1 1 2 ) 3=2 (2.23) into Eq. (2.22) yields R c = a(1 2 ) 2 [1 + 4r 2 xy a 2 (1 2 ) 2 ] 3=2 : (2.24) Chapter 2 19 Figure 2.4: Schematic diagram of the needle and the Taylor cone. At the tip of the liquid surface, x = 0, y = 0, and z = d, therefore = (r 1 r 2 )=a = 2d=a = 0 , here r 1 = d +a=2 and r 2 = a=2d. Substitute = 0 = 2d=a into Eq. (2.24), and notice that r xy = 0 when x = 0 and y = 0, one can nd that a = 2d(1 +R c =d) 1=2 ; 0 = (1 +R c =d) 1=2 : (2.25) 0.1 0.1 0.2 0.2 0.3 0.3 0.4 0.4 0.5 0.5 0.6 0.6 0.7 0.8 -1 0 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Figure 2.5: Potential plot with R n =d = 0:2 and y = 0. Therefore, if R c and d are specied, a and 0 can be determined. While the value of d depends on dierent applications, the value of R c can be obtained if the needle radius R n is given. Figure 2.4 schematically shows the needle (the cylinder) and the liquid surface (the Taylor cone). Since the Taylor cone angle is T = 49:29 [70], R c =R n = cos T , where R n is the radius of the needle. 20 Chapter 2 An example of the potential obtained from Eq. (2.16) choosingR n =d = 0:2 at the plane y = 0 is shown in Figure 2.5. The corresponding E x and E z obtained from Eq. (2.17), Eq. (2.18), and Eq. (2.19) are plotted in Figure 2.6 with the values normalized by the maximum electric eld in the plane y = 0, i.e. E max = max( p j^ xE acc j 2 +j^ zE acc j 2 ). At last, a zoomed-in view near the liquid surface tip is shown in Figure 2.7 with electric eld vectors on the potential contour. -0.15 -0.1 -0.05 0 0 0 0.05 0.1 0.15 0.2 -1 0 1 0 0.2 0.4 0.6 0.8 1 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 -0.6 -0.5 -0.4 -0.3 -0.3 -0.2 -0.2 -0.2 -0.2 -1 0 1 0 0.2 0.4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0 Figure 2.6: Electric eld plots. -0.4 -0.2 0 0.2 0.4 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Figure 2.7: Zoomed-in potential plot with electric eld vectors. 2.4.3 Numerical Integration The velocity Verlet algorithm [77] is used, as shown in Eq. (2.26), to numerically integrate the equations of motion. At time t, if the positionr(t), the velocityv(t) and the forceF (t) are known, the position at time t + t, i.e. r(t + t), can be calculated according to Eq. Chapter 2 21 (2.26). From r(t + t), the forces F pp (t + t) and F acc (t + t) can be obtained, thus the velocityv(t + t) can be updated. r(t + t) =r(t) +v(t)t + F (t) 2m t 2 v(t + t) =v(t) + F (t) +F (t + t) 2m t (2.26) 2.5 Droplet Initial Breakup The droplet breakup is considered as follows. A numerical parameter N s is introduced here, which denotes the number of droplets broken from the original droplet. If N s is set to be 1, the original droplets do not break. If N s is set to be other integers greater than 1, N s times of N in droplets will be injected during injection following the method introduced in Sec. 2.3, with the same amount of chargeq =Q d =N s and massm =M d =N s . In other words, no droplets break up during the simulation, all the droplet breakups are addressed in the initialization. Of course, this droplet breakup is not realistic, but is very handy to study the collisional eects of the colloid thruster beam as a rst attempt. In addition, one may want to consider the mechanism given in Sec. 1.2.4 to judge if a droplet should break. Here, we use the parameters of case A 0 that will be given in Sec. 4.1 to do a calculation. The maximum electric eld that balances the surface tension is E max = r 4 " 0 R d 2:189 V/nm: The internal electric eld is E in = q 4" 0 R 2 d 1:094 V/nm: The accelerating electric eld near the Taylor cone tip is about 0:0851 V/nm, and the particle-particle eld may be about 0:0415 V/nm (see Sec. 4.5.3 for detail). Therefore E ex +E in = 0:0851 + 0:0415 + 1:094 (V/nm) = 1:2206 (V/nm)<E max ; so the droplet should not break, if the electric eld and the surface tension are the only factors considered. Also note that, the subdroplets after breakup have higher stability, and the external eld in the downstream region of the beam is much less than that in the upstream. Thus, if a droplet should not break near the tip, it has less and less possibilities to break in the future. Therefore, we consider this articial droplet breakup to study the dynamics of subdroplets. 22 Chapter 2 2.6 Simulation Procedure and Techniques Start Droplet Injection Force Calculation Droplet Push End Output Figure 2.8: Simulation procedure. The simple simulation procedure is shown in Fig. 2.8. During each time step, new droplets are injected rst (Droplet Injection), then the force of each droplet is solved (Force Calcula- tion), at last the position and velocity of each droplet are updated (Droplet Push), and some simulation results may be outputted (Output). The simulation ceases when the simulation time reaches our needs. Note that if some droplets pass the accelerating plate atz = 0, they will be deleted in the simulation. Two types of parallel computing techniques are considered: (a) OpenMP [4] paralleliza- tion using CPU (central processing unit); (b) CUDA [1] parallelization using GPU (graphics processing unit). Both types apply the parallelization among the outer loop of the force calculation. For example, the OpenMP code would look like: !$OMP PARALLEL PRIVATE( i , j , . . . ) !$OMP DO do i =1,N do j =1,N . . . end do end do !$OMP END DO !$OMP END PARALLEL where the outer index i is parallelized among CPU threads. The CUDA code would look like: i =(blockidx%x1)∗blockdim%x+threadidx%x do j =1,N . . . end do where the outer index i is parallelized among the GPU threads. The simulation code is written in Fortran 90 [2]. Chapter 2 23 10 0 10 1 10 2 10 3 10 4 8 9 10 11 12 13 14 15 32768 16384 8192 4096 2048 1024 A typical run for 1000 time steps Time (s) log 2 N p CPU-1 CPU-16 GPU Figure 2.9: Computation time comparison between using CPU and GPU. A computation time comparison between using CPU (OpenMP) and GPU (CUDA) of a typical PP model is shown in Fig. 2.9, whereN p denotes the number of simulated particles, the GPU used is k20 NVIDIA with 2496 cores, the CPU used is Intel Xeon of 2.4 GHz. As can be seen, using GPU is much faster than using CPU without OpenMP parallelization (the CPU-1 line shown in the gure). When using OpenMP with 16 threads (the CPU-16 line shown in the gure), GPU is faster than CPU for N p > 2048. Therefore, GPU CUDA parallelization will be mainly used in this work. 2.7 Diagnoses 2.7.1 System Energy E tot =E tot =N = (E p +E k )=N E p = N X i=1 [q i acc (r i )] + 1 2 N X i=1 [q i N X j=1;j6=i ( 1 4" 0 q j jr ij j )] E k = N X i=1 (m i v 2 i =2) (2.27) System energy would be the most basic diagnosis for any simulations. In our system, the total energy includes the kinetic energy E k and the electric eld potential energy E p . Since we keep injecting droplets at a constant rate, the total energy of the system is increasing linearly. For convenience, we divide the total energy by the total number of droplets, thus we obtain the averaged energy of each droplet, which does not change much over time, as shown in Eq. (2.27). Note that the potential energy has two components, one is due to the accelerating eld, and the other is due to the Coulomb potential of other droplets. The diagnosis of E tot will be used as the rst validation of the simulation, especially to verify 24 Chapter 2 if the chosen time step t is small enough. If t is not small enough, E tot will not be conserved. 2.7.2 Equivalent length One quantity named equivalent length denoted as e is introduced here, its denition can be given as follows: the equivalent length of a test particle is the distance from this test particle to an imagined eld particle with a specied charge, which is located at a unique position, such that the Coulomb force acting on this test particle due to the imagined eld particle is equivalent to the resultant Coulomb force acting on this test particle due to all other real particles. For example, imagine there are three particles, P , P 1 , and P 2 , as shown in Fig. 2.10. According to Coulomb's law, P feels forces F 1 and F 2 caused by P 1 and P 2 , respectively. The resultant force acting on P isF =F 1 +F 2 . SinceF can also be obtained equivalently by placing a single particle P at an appropriate and unique location, if the charge of P is specied, the equivalent length e is thus dened as the distance between P and P . Since we can nd e for each particle, a probability density function (PDF) f p ( e ) can be obtained from these e , which statistically re ects how collisional the system is. Fig. 2.11 gives an example of f p ( e ) of two cases A and B. As can be seen, A has a higher possibility for smaller e than B, because smaller e results in larger force, we can say that the system A is more collisional than B, if A and B have the same temperature. Figure 2.10: Diagram for understanding the equivalent length. By analogy, for the system of charged droplets, if the charge of the imagined droplet is specied, we can nd e for each droplet, thus f p ( e ) can be obtained to re ect how collisional the system is. However, the resultant forceF of each droplet in our system is not only due to short-range Coulomb collisions, but includes the accelerating eld as well as the long-range eld caused by the charged beam body. One way to only consider the short-range Coulomb collisions is as follows. Chapter 2 25 f p ( e ) e A B Figure 2.11: Probability density function of equivalent length. First, we neglect thez component of the resultant force vector and the position vector, F (F x ;F y ;F z )!F xy (F x ;F y ; 0); r(r x ;r y ;r z )!r xy (r x ;r y ; 0); (2.28) as shown in Fig. 2.12 (left), since thez component is mainly caused by the accelerating eld. Then, we splitF xy into two components, one is parallel tor xy , the other is perpendicular to r xy , F xy =F r +F t jF r j =jF xy ^ r xy j jF t j =jF xy (^ r xy ^ z)j; (2.29) where ^ r xy = r xy =jr xy j, as shown in Fig. 2.12 (right). F r is due to both the short-range Coulomb collisions and the long-range repulsion caused by other droplets if the evaluated droplet is not on the z axis, namely the further the droplet away from the symmetric axis of the beam cone, the more F r is due to long-range repulsion rather than the short-range Coulomb collisions. Fortunately, F t is believed to be only contributed by the short-range Coulomb collisions, thus we useF t to solve for the equivalent length, et = ( q 2 4" 0 jF t j ) 1=2 ; (2.30) where we have chosen the charge of the imagined droplet to be the same as the evaluated droplet. Of course, it may also be instructive to see the equivalent length obtained fromF z and F r , but the components caused by the accelerating eld (E acc ) should be subtracted before 26 Chapter 2 evaluation: ez = ( q 2 4" 0 jF z qE acc ^ zj ) 1=2 ; (2.31) Figure 2.12: Components of the resultant force. er = ( q 2 4" 0 j[F xy q(E acc ^ x)^ xq(E acc ^ y)^ y] ^ r xy j ) 1=2 : (2.32) Also note that the projection ofE acc onF t is zero. Therefore, based on et , ez and er of each droplet, we can obtain f p ( et ), f p ( ez ), andf p ( er ), among whichf p ( et ) only re ects the short-range collisions, f p ( ez ) andf p ( er ) re ect the superposition of short-range collisions and long-range forces. 2.7.3 Droplet Macro-temperature In kinetic theory, the temperature T of particles in equilibrium is dened as 1 2 m<v 2 > = 3 2 k B T; (2.33) where k B is the Boltzmann constant, the angle brackets <> denote the statistical average among particles, v 2 =v 2 x +v 2 y +v 2 z , andT may be split into three directions T = (T x +T y + T z )=3. Although in our beam system the state is not equilibrium and the droplets are much larger than those small particles, it is still instructive to diagnose the \temperature" using Chapter 2 27 the same calculation as Eq. 2.33 and treat each droplet as a particle. Obviously, one should note that the calculated T indicates the thermal energy of droplets of the beam system, instead of the liquid temperature of each droplet, which re ects the thermal energy of the molecules within each droplet. Under our circumstances, instead of considering x, y, and z components, it is better to split T into T z , T r , and T t , namely the axial, the radial, and the tangential components. To evaluate the thermal motion only, we need to subtract the collective motion < v i >. Therefore the three components of T can be calculated from T i = m k B (<v 2 i ><v i > 2 ); i =z;r;t: (2.34) where v r and v t can be obtained using similar procedures as that shown in Sec. 2.7.2, v r =jv ^ r xy j; v t =jv (^ r xy ^ z)j: (2.35) In sum, T z , T r , and T t re ect the thermal energy of the random movements of droplets in z, r, and t directions, respectively, and T = (T z +T r +T t )=3. To distinguish T from the normal temperature, we call it macro-temperature. 2.7.4 Thrust and Specic Impulse The mean droplet exit speed v z (z = 0) can be obtained from the simulation, thus the thrust can be obtained from F T = _ V v z : (2.36) The specic impulse can be obtained from I sp = v z =g 0 ; (2.37) where g 0 = 9:8 ms 2 is the acceleration due to gravity at the ground level. 28 Chapter 3 Chapter 3 PP Code Validation In this chapter, we verify the PP model via a time step length verication using the binary collision theory, and via a space charge eect study using Child-Langmuir law. 3.1 Time Step Length Verication Figure 3.1: Binary collision setup. To verify if a chosen time step length is small enough, one way we consider is to apply the PP model on the scenario of the binary collision theory. The binary collision theory [8] gives us the scattering angle , tan 2 = b 0 b ; b 0 = 1 4" 0 qq 0 g 2 ; (3.1) where b is the impact parameter, g is the relative speed of two particles having charges q, q 0 , and masses m, m 0 , =mm 0 =(m +m 0 ) is the reduced mass. The verication of the time Chapter 3 29 step length is thus a comparison between the simulated scattering angle with this analytical solution. The setup of the PP model is shown in Fig. 3.1, where two charged particles q and q 0 move toward each other with relative speed g, separation L in x direction, b in y direction. Since L cannot be innitely large in the simulation, we set L = 1000b for all cases. Case q m g b t Error 1 3.77E-18 C 8.15E-22 kg 257 m/s 2.94 nm 5.15E-11 s 2.51 % 2 2.36E-19 C 5.09E-23 kg 257 m/s 2.94 nm 3.22E-22 s 0.05 % 3 1.60E-19 C 1.85E-25 kg 751 m/s 0.92 nm 1.28E-13 s 0.07 % Table 3.1: Time step length verication. Three cases are considered with parameters given in Tab. 3.1. These parameters are chosen to be the same as some simulation cases that will be presented in later chapters. The Error is the relative error of the simulated scattering angle with respect to the analytical scattering angle. The particle trajectories of the case 1 are given in Fig. 3.2 as an example. As we can see from Tab. 3.1, all errors are small, thus the t used should be small enough. Since the simulated cases produce the largest scattering angles among the simulation cases that will be presented later, we can say that t should be small enough for all simulation cases involved in this work. 0:1 0:05 0 0:05 0:1 0:1 0:05 0 0:05 0:1 y (m) x (m) Figure 3.2: Simulated binary collision of Case 1. 30 Chapter 3 3.2 Space Charge Eect Validation 3.2.1 Child-Langmuir Law The simplest law of space charge limit would be the Child-Langmuir law [18, 46, 74]. It describes the scenario of two innitely long parallel plates, one serves as an anode, the other serves as a cathode which has an innite supply of free electrons at rest. If the distance between the two plates isd, the potential dierence is 0 , the maximum current density that can be drawn across this gap is thus the space-charge-limited current density, J SCL = 4 9 " 0 r 2e m 3=2 0 d 2 ; (3.2) where " 0 is the vacuum permittivity, e is the unit charge, m is the electron mass. If we assume the ow is along the z direction from z = 0 to z =d, and the potential at z = 0 is (0) = 0, we have the solution of the potential distribution, (z) = 3 2 4=3 J SCL " 0 2=3 m 2e 1=3 z 4=3 : (3.3) Since E z = d dz , we have the solution of the axial electric eld distribution, E z (z) = 4 3 3 2 4=3 J SCL " 0 2=3 m 2e 1=3 z 1=3 : (3.4) In the following sections, we will simulate the scenario of electrons with a zero initial velocity, and ions with an nonzero initial velocity. But since the nonzero initial velocity is much smaller than the nial particle velocity, Child-Langmuir law can still be applied approximately. 3.2.2 The PIC Model First, we use the 3D PIC to simulate a cubic domain. The xed potential boundary condition (Dirichlet) is applied on z = 0 and z = d, (0) = 0, (d) = 0 . Those particles that travel out of the domain in thez direction are deleted. The periodic boundary condition is applied on x and y directions for both the eld and particles. The parameters are chosen to be as follows. The potential dierence is 0 = 1 V. The separation of the two plates is d = 1 m, which is the length of the simulation domain in z direction. So, the constant electric eld without space charge eect is E 0 = 0 =d = 1 V/m. The length of the simulation domain in x and y directions are also set to be d. Thus the simulation domain is a cubic box. The simulated particles are electrons with chargee and mass m, e = 1:6 10 19 C, m = 9:1 10 31 kg. The resultant space-charge-limited current Chapter 3 31 density is J SCL 2:332 A/m 2 . 0 0:2 0:4 0:6 0:8 1 0 0.2 0.4 0.6 0.8 1.0 (z) (V) z (m) Child PIC 1:4 1:2 1 0:8 0:6 0:4 0:2 0 0 0.2 0.4 0.6 0.8 1.0 E z (z) (V/m) z (m) Child PIC Figure 3.3: The potential distribution (left), the electric eld distribution (right). 5000 4000 3000 2000 1000 0 1000 0 0.2 0.4 0.6 0.8 1.0 (z) (V) z (mm) 1 10 6 0 1 10 6 2 10 6 3 10 6 4 10 6 5 10 6 6 10 6 7 10 6 0 0.2 0.4 0.6 0.8 1.0 E z (z) (V/m) z (mm) Child 0:1J SCL 0:5J SCL 1:0J SCL 2:0J SCL 5:0J SCL Figure 3.4: The potential distribution (left), the electric eld distribution (right). We rst let the emitted current density to be J =J SCL . Thus the current is I = Jd 2 =2:332A. we choose the injection period to be t = 0:001t 0 , i.e. we inject electrons every t, where t 0 = p 2dm=eE 0 3:373 s is the single electron travel time between the two plates under the constant electric eld E 0 . There should be _ N R =jI=ej 1:458 10 13 32 Chapter 3 real electrons injected per second. We choose to inject two simulated electrons per injection, N in = 2, so the macro-particle weight is N W = _ N R t=N in 24583, meaning that each simulated macro-particle represents N W real electrons. Within each injection period, we solve the Poisson's equation and push particles N loop = 2 times, so the time step length used to integrate the equations of motions is t 0 = t=N loop . The number of grid points are 50 50 50, which results in a cell size x = y = z = 0:02 m. The maximum speed a single electron can gain is v max = p 2 0 e=m 5:930 10 5 m/s. So, the time step length should be small enough to satisfyv max t 0 = 0:001(m)< x. The maximum simulation time is chosen to be t max = 10t 0 , we found a steady state achieved after about 5t 0 . The simulation results of the potential distribution and the electric eld distribution averaged using all grid points in x and y along z are given in Fig. 3.3. We can see that the PIC results match very well with the Child-Langmuir law, thus the PIC model is valid. Then, we use the parameters of our case A 0 that will be given in the following chapter. q = 3:769 10 18 C, m = 8:149 10 22 kg, 0 =5000 V, d = 1 mm, v 0 = 257:1 m/s. Note that the particles have an initial axial velocity v 0 this time. The value of v 0 = 257:1 m/s is much smaller than the nal accelerated velocity of p 2q 0 =m 6800 m/s. We vary the injected current density to be 0.1, 0.5, 1.0, 2.0, and 5.0 times of J SCL by increasing the number density while keep the emission velocity at v 0 . The results are shown in Fig. 3.4. Since the particle emission velocity is much smaller than the nal accelerated velocity, the potential prole for J >J SCL in this test still approximately matches the monotonic prole at the space charge limit. 3.2.3 The PP Model Now, we have a valid PIC, so we will mainly compare the PP model with the PIC model. First, for the PIC model, we still inject electrons from the whole surface at z = 0 which ranges from x = 0 to 1 m, y = 0 to 1 m, but change the periodic boundary condition to Neumann, i.e. the normal electric eld E n = 0 on the corresponding surface. For particles, we still apply the periodic boundary condition when they travel out of the domain from x or y direction. Then, we simply replace the Poisson eld solver by the PP eld solver (Coulomb's law) to obtain the PP model. One should note that in this PP model, macro-particles are used and a cuto is applied when solving the electric eld from Coulomb's law. Coulomb's eld becomes innity as the inter-particle distancer approaches zero. In PIC, the eld within the range of a cell is reduced to zero, thus we articially reduce the eld of PP as well. Ideally, we should apply a eld function that is the same as the one produced by PIC. Here, for simplicity we apply a cuto and let the eld reduce to zero linearly as shown in Fig. 3.5, because in the PIC model, we use linear interpolation. This cuto model has been used previously in a work to simulate ion beam neutralization [89]. Here, we choose r cf = x. The boundary condition for particles are the same as that of PIC. No boundary condition for the eld is used, meaning that a constant electric eldE 0 is added when pushing particles, Chapter 3 33 and no periodic boundary condition is used in x and y directions. 0 kq=r 2 cf E 0 r cf r Figure 3.5: The modied PP eld. 0 0:2 0:4 0:6 0:8 1 0 0.2 0.4 0.6 0.8 1.0 (z) (V) z (m) Child PP PIC 1:4 1:2 1 0:8 0:6 0:4 0:2 0 0 0.2 0.4 0.6 0.8 1.0 E z (z) (V/m) z (m) Child PP PIC Figure 3.6: The potential distribution (left), the electric eld distribution (right). The simulation results are shown in Fig. 3.6, which are the results along the beam center, not averaged as Fig. 3.3. We can see that, the PIC results still match well with Child-Langmuir law, the PP results also match with Child-Langmuir law, but with a small deviation. The small deviation is believed due to the fact that the eld function of PP we 34 Chapter 3 applied does not match well with that of PIC. However, since our focus is to apply PP to simulate real particles using the full Coulomb's law, this small deviation here does not matter much. 0 0:2 0:4 0:6 0:8 1 0 0.2 0.4 0.6 0.8 1.0 (z) (V) z (m) Child PP PIC 1:4 1:2 1 0:8 0:6 0:4 0:2 0 0 0.2 0.4 0.6 0.8 1.0 E z (z) (V/m) z (m) Child PP PIC Figure 3.7: The potential distribution (left), the electric eld distribution (right). 0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 x (m) z (m) E z (V/m) 1:1 1:05 1 0:95 0:9 0:85 0:8 0:75 0:7 0:65 0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 x (m) z (m) E z (V/m) 1:1 1:05 1 0:95 0:9 0:85 0:8 0:75 0:7 0:65 Figure 3.8: The electric eld distribution of PIC (left), the electric eld distribution of PP (right). Then, for the same PIC and PP model, we reduce the injection surface to be a circle centered at (0:5; 0:5; 0) m on thez = 0 plane, and with a radiusR = 0:25 m. We still choose Chapter 3 35 the emitted current density J =J SCL , thus the current changes to I = JR 2 0:458 A. Other parameters are still the same, except those that are related to the current will be changed accordingly. 5000 4000 3000 2000 1000 0 1000 0 0.2 0.4 0.6 0.8 1.0 (z) (V) z (mm) 5000 4000 3000 2000 1000 0 1000 0 0.2 0.4 0.6 0.8 1.0 (z) (V) z (mm) Figure 3.9: The potential distribution of PIC (left) and PP (right). 2 10 6 0 2 10 6 4 10 6 6 10 6 0 0.2 0.4 0.6 0.8 1.0 E z (z) (V/m) z (mm) Child 0:1J SCL 1:0J SCL 5:0J SCL 2 10 6 0 2 10 6 4 10 6 6 10 6 0 0.2 0.4 0.6 0.8 1.0 E z (z) (V/m) z (mm) Child 0:1J SCL 1:0J SCL 5:0J SCL Figure 3.10: The electric eld distribution of PIC (left) and PP (right). The simulation results of the potential and electric eld distributions are shown in Fig. 3.7. We can see that this time both PIC and PP deviate from the Child-Langmuir law. The 36 Chapter 3 main reason is believed to be that the electron beam is allowed to expand radially, which will weaken the axial space charge eects. However, since PP and PIC still match with each other, we can say PP is valid. Also, we plotted the 2D electric eld color map in the y = 0:5 m plane. We can see that, overall the PP results match with that of PIC, but the electron beam of PP expands slightly more radially, which is again because the eld function of PP is not exactly the same as that of PIC, and it produces slightly larger forces. Then, we use the parameters of our case A 0 that will be given in the following chapter. q = 3:769 10 18 C, m = 8:149 10 22 kg, 0 =5000 V, d = 1 mm, v 0 = 257:1 m/s. We reduce the injection surface to be R = 0:25 mm centered in the z = 0 plane. We vary the injected current density to be 0.1, 1.0, and 5.0 times of J SCL . Again both PIC and PP are tested. The results are shown in Fig. 3.9 and Fig. 3.10. We can see that again both PP and PIC deviate from Child's law, but PP still matches with PIC. The larger deviation between PP and PIC of the case with J = 5J SCL is believed due to the dierent force functions applied, which is enhanced for a larger current. Chapter 4 37 Chapter 4 Characteristic Simulation of Droplet Beam 4.1 Chosen Parameters and Simulation Cases In this chapter, some reasonable parameters are chosen, and several characteristic simula- tions are executed. The chosen material parameters for all simulations in this chapter are given in Tab. 4.1, which correspond to those of formamide (CH 3 NO). The chosen practical parameters are given in Tab. 4.2. We choose d = d 0 for simplicity, meaning that the jet breaks immediately after emitted from the Taylor cone. This assumption is reasonable due to the fact that the characteristic lengthL of the cone-jet transition region can be calculated from L 2R 0 _ V= _ V 0 24:22 nmd 0 , where 2R 0 = [ " 2 0 =(K 2 )] 1=3 , _ V 0 = " 0 =(K). Based on Tab. 4.1 and Tab. 4.2, other relevant parameters can be calculated according to Sec. 2.2 and are listed in Tab. 4.3. The simulation only needs the following input parameters: Q d , M d , R j , v d , d 0 , 0 , a, and 0 . Parameter name Symbol Value SI unit Conductivity K 1 S m 1 Relative permittivity " 100 Dimensionless Surface tension 0.059 N m 1 Mass density 1130 kg m 3 Empirical current factor f(") 18 Dimensionless Table 4.1: Chosen material parameters (formamide). Some numerical parameters need to be specied, such as t and N s . t will be varied to test the conservation of the system energy. N s will also be varied to study the eects of droplet breakup. Five simulation cases are listed in Tab. 4.4, where N s is varied from 1 to 38 Chapter 4 16, and t is varied three times for each N s . All the simulations were run using k20 NVIDIA GPU CUDA parallelization at the Center for High Performance Computing of the University of Southern California [3]. The computation time varies for dierent cases, among which the caseA 4 with t _ N d = 1=48 has the longest computation time which is about 24 hours. Parameter name Symbol Value SI unit Volume ow rate _ V 7 10 15 m 3 s 1 Distance between the liquid tip and the plate d 1 mm Distance between the jet front and the plate d 0 1 mm The accelerating voltage 0 5 kV Needle radius R n 15 m Table 4.2: Chosen practical parameters. Parameter name Symbol Value SI unit Flow current I 36.58 nA Characteristic dimension r 18.37 nm Radius of jet R j 2.94 nm Radius of original droplet R d 5.56 nm Mass of original droplet M d 8.15 10 22 kg Charge of original droplet Q d 3.77 10 18 C Initial axial droplet speed v d 257.18 m s 1 Emission rate of original droplet _ N d 9.71 10 9 s 1 Radius of curvature R c 19.79 m Half focal length a 2.02 mm coordinate of liquid surface 0 0.99 Dimensionless Table 4.3: Derived parameters. Case N s t _ N d =N in A 0 1 1/1, 1/2, 1/3 A 1 2 1/2, 1/4, 1/6 A 2 4 1/4, 1/8, 1/12 A 3 8 1/8, 1/16, 1/24 A 4 16 1/16, 1/32, 1/48 Table 4.4: Simulation cases. Chapter 4 39 4.2 Simulation Results 4.2.1 Phase Space 0:1 0:05 0 0:05 0:1 0:1 0:05 0 0:05 0:1 0 1 z=d 0 x=d 0 y=d 0 z=d 0 Figure 4.1: Droplet spatial distribution. First, Fig. 4.1 gives a general idea of the 3-D spatial distribution of the colloid thruster beam, where droplets are injected from the injection volume located at (0; 0;d 0 ), accelerated along z direction and expand in the radial direction. More detailed and quantitative comparisons between cases are given in 2-D plots as follows. Fig. 4.2 shows thex-z phase space plots of the ve cases fromA 0 toA 4 at time 2000= _ N d , respectively, and each plotted dot represents a droplet. In each row of Fig. 4.2, the left panel has the corresponding largest t, and we can see that some droplets move out of the main beam region, while the middle and the right panels do not show outsiders of droplets. The reason is that the t used for the cases of the left panels are not small enough, thus some droplets gain non-physical large energies due to coarse numerical integration. Fig. 4.2 also shows that the beam cone angle increases as N s increases, thus the beam is expected to expand more in the radial direction, when more smaller droplets are broken from the original large droplets. However, one should note that the beam radius (measured at z = 0) increases with a decreasing rate as N s increases. Also, based on the beam size and shape obtained, the x and y components of the accelerating eld are negligible, which has been tested using other simulations not shown here. Therefore, all radial expansion of the beam is due to droplet Coulomb interactions. The phase space plot ofz-v z at time 2000= _ N d is shown in Fig. 4.3, which is the same for allA cases with all dierent t. We can see that all the injected droplets are accelerated in thez direction from v d to about 25v d . Therefore, the thrust obtained from the simulation of one needle is about F T = 25v d _ V 5:1 10 8 N, and the specic impulse is about I sp = 25v d =g 0 656 s. 40 Chapter 4 0 0:2 0:4 0:6 0:8 1 -0.1 0 0.1 z=d 0 x=d 0 t _ N d = 1/1 0 0:2 0:4 0:6 0:8 1 -0.1 0 0.1 z=d 0 x=d 0 t _ N d = 1/2 0 0:2 0:4 0:6 0:8 1 -0.1 0 0.1 z=d 0 x=d 0 t _ N d = 1/3 0 0:2 0:4 0:6 0:8 1 -0.1 0 0.1 z=d 0 x=d 0 t _ N d = 1/2 0 0:2 0:4 0:6 0:8 1 -0.1 0 0.1 z=d 0 x=d 0 t _ N d = 1/4 0 0:2 0:4 0:6 0:8 1 -0.1 0 0.1 z=d 0 x=d 0 t _ N d = 1/6 0 0:2 0:4 0:6 0:8 1 -0.1 0 0.1 z=d 0 x=d 0 t _ N d = 1/4 0 0:2 0:4 0:6 0:8 1 -0.1 0 0.1 z=d 0 x=d 0 t _ N d = 1/8 0 0:2 0:4 0:6 0:8 1 -0.1 0 0.1 z=d 0 x=d 0 t _ N d = 1/12 0 0:2 0:4 0:6 0:8 1 -0.1 0 0.1 z=d 0 x=d 0 t _ N d = 1/8 0 0:2 0:4 0:6 0:8 1 -0.1 0 0.1 z=d 0 x=d 0 t _ N d = 1/16 0 0:2 0:4 0:6 0:8 1 -0.1 0 0.1 z=d 0 x=d 0 t _ N d = 1/24 0 0:2 0:4 0:6 0:8 1 -0.1 0 0.1 z=d 0 x=d 0 t _ N d = 1/16 0 0:2 0:4 0:6 0:8 1 -0.1 0 0.1 z=d 0 x=d 0 t _ N d = 1/32 0 0:2 0:4 0:6 0:8 1 -0.1 0 0.1 z=d 0 x=d 0 t _ N d = 1/48 Figure 4.2: Phase space plots of x-z of all A cases at time 2000= _ N d . Chapter 4 41 0 5 10 15 20 25 0 0:2 0:4 0:6 0:8 1 jv z j=v d z=d 0 Figure 4.3: Phase space plots of z-v z of all A cases at time 2000= _ N d . 0 0:5 1 1:5 2 0 0:2 0:4 0:6 0:8 1 v xy =v d z=d 0 A 0 0 0:5 1 1:5 2 0 0:2 0:4 0:6 0:8 1 v xy =v d z=d 0 A 1 0 0:5 1 1:5 2 0 0:2 0:4 0:6 0:8 1 v xy =v d z=d 0 A 2 0 0:5 1 1:5 2 0 0:2 0:4 0:6 0:8 1 v xy =v d z=d 0 A 3 0 0:5 1 1:5 2 0 0:2 0:4 0:6 0:8 1 v xy =v d z=d 0 A 4 Figure 4.4: Phase space plots of z-v xy of all A cases at time 2000= _ N d . 0 0:5 1 1:5 2 0 0.05 0.1 v xy =v d r xy =d 0 A 0 0 0:5 1 1:5 2 0 0.05 0.1 v xy =v d r xy =d 0 A 1 0 0:5 1 1:5 2 0 0.05 0.1 v xy =v d r xy =d 0 A 2 0 0:5 1 1:5 2 0 0.05 0.1 v xy =v d r xy =d 0 A 3 0 0:5 1 1:5 2 0 0.05 0.1 v xy =v d r xy =d 0 A 4 Figure 4.5: Phase space plots of r xy -v xy of all A cases at time 2000= _ N d . 42 Chapter 4 0 5 10 15 20 25 0 0.05 0.1 jv z j=v d r xy =d 0 A 0 0 5 10 15 20 25 0 0.05 0.1 jv z j=v d r xy =d 0 A 1 0 5 10 15 20 25 0 0.05 0.1 jv z j=v d r xy =d 0 A 2 0 5 10 15 20 25 0 0.05 0.1 jv z j=v d r xy =d 0 A 3 0 5 10 15 20 25 0 0.05 0.1 jv z j=v d r xy =d 0 A 4 Figure 4.6: Phase space plots of r xy -v z of all A cases at time 2000= _ N d . 0 5 10 15 20 25 0 0:5 1 1:5 2 jv z j=v d v xy =v d A 0 0 5 10 15 20 25 0 0:5 1 1:5 2 jv z j=v d v xy =v d A 1 0 5 10 15 20 25 0 0:5 1 1:5 2 jv z j=v d v xy =v d A 2 0 5 10 15 20 25 0 0:5 1 1:5 2 jv z j=v d v xy =v d A 3 0 5 10 15 20 25 0 0:5 1 1:5 2 jv z j=v d v xy =v d A 4 Figure 4.7: Phase space plots of v xy -v z of all A cases at time 2000= _ N d . Fig. 4.4 shows the z-v xy plots at time 2000= _ N d of all A cases with t _ N d = 1=(2N s ), where v xy = (v 2 x +v 2 y ) 1=2 . As can be seen, the increase of N s results in larger v xy , which agrees with the greater radial expansion shown before. In addition, in the region near the accelerating plate, for example z=d 0 2 [0; 0:4], v xy seems to be approaching a steady state, namely the radial accelerations of droplets become smaller and smaller as z decreases, so droplets move with almost constant radial velocities eventually. Fig. 4.5 shows the r xy -v xy plots at time 2000= _ N d of all A cases with t _ N d = 1=(2N s ), where r xy = (x 2 +y 2 ) 1=2 . A straight slant line can be seen, which indicates that only those droplets with large v xy can arrive at large r xy , and only those droplets with small v xy can stay near the axis r xy = 0. At last, the r xy -v z plots at time 2000= _ N d are given in Fig. 4.6 and the v xy -v z plots at time 2000= _ N d are given in Fig. 4.7 for completeness. Chapter 4 43 4.2.2 System Energy The system energy is evaluated using Eq. 2.27. The top left panel of Fig. 4.8 shows the averaged total energy per droplet E tot , the averaged potential energy per droplet E p , and the averaged kinetic energy per droplet E k of the caseA 0 with t _ N d = 1=2. As can be seen, E p drops and E k raises such that E tot conserves over time. Since the results of all other cases are similar to this panel, we do not present them redundantly. The other ve panels of Fig. 4.8 show E tot of the ve cases with three curves plotted in each panel indicating the results obtained using dierent t, where E 0 equals the initial output value of E tot of each corresponding case. The panel of A 0 (top right) shows the E tot obtained using t _ N d = 1=1, t _ N d = 1=2, and t _ N d = 1=3, respectively, from the top curve to the bottom curve. As can be seen, the total system energy conserves better when using smaller t. The panels of A 1 , A 2 , A 3 , and A 4 show similar results as that of A 0 , but curves of t _ N d = 1=(2N s ) and t _ N d = 1=(3N s ) of these four cases are so close that they overlap each other in these plots. 0 0:2 0:4 0:6 0:8 1 0 400 800 1200 1600 2000 E tot , E p , E k t _ N d A 0 (t _ N d = 1=2) 0:99996 1:00012 1:00028 1:00044 0 400 800 1200 1600 2000 E tot = E 0 t _ N d A 0 1 1:0004 1:0008 1:0012 0 400 800 1200 1600 2000 E tot = E 0 t _ N d A 1 1 1:0004 1:0008 1:0012 0 400 800 1200 1600 2000 E tot = E 0 t _ N d A 2 0:996 0:997 0:998 0:999 1 0 400 800 1200 1600 2000 E tot = E 0 t _ N d A 3 0:96 0:97 0:98 0:99 1 0 400 800 1200 1600 2000 E tot = E 0 t _ N d A 4 Figure 4.8: Time history of system energy. To present the energy results more quantitatively, Tab. 4.5 is given, in which the values are the percent error with respect to E 0 of each case evaluated at time t _ N d = 2000, i.e. 44 Chapter 4 [ E tot (t _ N d = 2000)= E 0 1] 100%. Therefore, considering both Tab. 4.5 and the outsiders of droplets shown in Fig. 4.1, we can say that when t satises t _ N d = 1=(2N s ), the simulation results are valid with respect to system energy conservation. Also, as shown in Sec. 2 Tab. 3.1, the case 1 corresponds to the case A 0 with t _ N d = 1=2, since the relative error is 2.51 %, we can say that the integration of the equations of motion Eq. (2.11) is accurate. A 0 A 1 A 2 A 3 A 4 t _ N d = 1=(1N s ) 0.0328 0.101 0.0610 0.358 3.97 t _ N d = 1=(2N s ) 0.0175 0.0229 0.0201 0.0161 0.0161 t _ N d = 1=(3N s ) 0.0207 0.0240 0.0203 0.0195 0.0212 Table 4.5: Percent error (%) of system energy of each case. 4.2.3 Number Density and Macro-temperature A 0 0 0:01 0:02 0:03 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 14 15 16 17 18 A 1 0 0:02 0:04 0:06 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 14 15 16 17 18 A 2 0 0:04 0:08 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 14 15 16 17 18 A 3 0 0:03 0:06 0:09 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 14 15 16 17 18 A 4 0 0:05 0:1 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 14 15 16 17 18 Figure 4.9: Droplet number densityn (m 3 ) ofA 0 toA 4 ve cases, averaged fromt _ N d = 1820 to 2000. The value of log 10 n is shown using colors. The droplet number density n of A 0 to A 4 ve cases, averaged from t _ N d = 1820 to 2000, are plotted in Fig. 4.9. Note that the horizontal axis r xy varies from case to case, because the beam size is dierent, and in the time period used for calculating the average, n does Chapter 4 45 not change much (not shown here). Similar ways of calculating averages will also be carried out in later diagnoses of macro-temperatures. As can be seen, ve cases produce similar n, n decreases signicantly along the downstream axial direction, and increases along the outer radial direction. This result is reasonable, since the beam expand signicantly, and all droplets are charged positively, thus those droplets in the center of the beam are not stable and will be repelled toward the outer radial direction. For all ve cases,n ranges from about 10 18 to 10 15 m 3 . A 0 0 0:01 0:02 0:03 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 1 2 3 4 5 6 A 1 0 0:02 0:04 0:06 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 1 2 3 4 5 6 A 2 0 0:04 0:08 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 1 2 3 4 5 6 A 3 0 0:03 0:06 0:09 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 1 2 3 4 5 6 A 4 0 0:05 0:1 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 1 2 3 4 5 6 Figure 4.10: Droplet macro-temperature in axial direction T z (K) of A 0 to A 4 ve cases, averaged from t _ N d = 1820 to 2000. The value of log 10 (T z N s ) is shown using colors. The droplet macro-temperature in the axial directionT z ofA 0 toA 4 ve cases, averaged from t _ N d = 1820 to 2000, are plotted in Fig. 4.10. As can be seen, T z is high near the injection volume, because droplets are accelerated signicantly there, and E acc has a large gradient along z direction, thus two droplets separated with a small distance will have a large velocity dierence. From A 0 to A 4 , the mass of droplet reduces with a factor of N s , since T/ m, we plotted T z N s to only re ect the thermal speed of v z . Except large noises occurred inA 0 andA 1 due to less number of droplets, dierent cases produce similar T z . T z drops along the downstream direction from about 50000 to 5000 K. The droplet macro-temperature in the radial directionT r ofA 0 toA 4 ve cases, averaged from t _ N d = 1820 to 2000, are plotted in Fig. 4.11. Again, we plotted T r N s to only re ect the thermal speed of v r . As can be seen, T r is higher in the upstream region than that in the downstream region, because the larger radial acceleration in the upstream region results in a larger thermal speed of v r . From case A 0 to A 4 , it is found that T r increases. 46 Chapter 4 A 0 0 0:01 0:02 0:03 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 1:5 2 2:5 3 3:5 4 4:5 5 A 1 0 0:02 0:04 0:06 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 1:5 2 2:5 3 3:5 4 4:5 5 A 2 0 0:04 0:08 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 1:5 2 2:5 3 3:5 4 4:5 5 A 3 0 0:03 0:06 0:09 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 1:5 2 2:5 3 3:5 4 4:5 5 A 4 0 0:05 0:1 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 1:5 2 2:5 3 3:5 4 4:5 5 Figure 4.11: Droplet macro-temperature in radial direction T r (K) of A 0 to A 4 ve cases, averaged from t _ N d = 1820 to 2000. The value of log 10 (T r N s ) is shown using colors. A 0 0 0:01 0:02 0:03 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 1 0:5 0 0:5 1 1:5 2 A 1 0 0:02 0:04 0:06 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 1 0:5 0 0:5 1 1:5 2 A 2 0 0:04 0:08 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 1 0:5 0 0:5 1 1:5 2 A 3 0 0:03 0:06 0:09 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 1 0:5 0 0:5 1 1:5 2 A 4 0 0:05 0:1 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 1 0:5 0 0:5 1 1:5 2 Figure 4.12: Droplet macro-temperature in tangential direction T t (K) of A 0 to A 4 ve cases, averaged from t _ N d = 1820 to 2000. The value of log 10 (T t N s ) is shown using colors. Chapter 4 47 The droplet macro-temperature in the tangential direction T t of A 0 to A 4 ve cases, averaged from t _ N d = 1820 to 2000, are plotted in Fig. 4.12. Again, we plotted T t N s to only re ect the thermal speed of v t . As can be seen, the absolute values of T t are extremely low for all cases, and T t decreases from A 0 to A 4 , meaning that when more droplets are broken, the thermal speed of v t decreases. 4.2.4 Velocity Distribution Function Since the beam has a cone shape, it is very non-uniform spatially. We divide the beam into three regions for evaluating f v as well as f p . Tab. 4.6 shows the three regions of each case, in which the intervals of r xy vary due to the fact that the beam radial sizes are dierent from case to case. The three regions in the beam are illustrated in Fig. 4.13. Overall, region (1) represents the upstream region, region (2) represents the central downstream region, and region (3) represents the outer downstream region. Region (1) Region (2) Region (3) z=d 0 r xy =d 0 z=d 0 r xy =d 0 z=d 0 r xy =d 0 A 0 (0:5; 1) (0; 0:015) (0; 0:5) (0; 0:015) (0; 0:5) (0:015;1) A 1 (0:5; 1) (0; 0:030) (0; 0:5) (0; 0:030) (0; 0:5) (0:030;1) A 2 (0:5; 1) (0; 0:040) (0; 0:5) (0; 0:040) (0; 0:5) (0:040;1) A 3 (0:5; 1) (0; 0:045) (0; 0:5) (0; 0:045) (0; 0:5) (0:045;1) A 4 (0:5; 1) (0; 0:045) (0; 0:5) (0; 0:045) (0; 0:5) (0:045;1) Table 4.6: Regions for evaluating f v and f p of all cases. Figure 4.13: Illustration of the three regions used for evaluating f v and f p . 48 Chapter 4 0 0:4 0:8 1:2 1:6 2 0 5 10 15 20 25 30 t _ N d =1000 v z =v d 0 0:2 0:4 0:6 0:8 1 0 0:4 0:8 1:2 1:6 2 0 5 10 15 20 25 30 t _ N d =1000 v z =v d 0 0:2 0:4 0:6 0:8 1 0 0:4 0:8 1:2 1:6 2 0 5 10 15 20 25 30 t _ N d =1000 v z =v d 0 0:2 0:4 0:6 0:8 1 0 0:4 0:8 1:2 1:6 2 0 0:5 1 1:5 2 t _ N d =1000 v r =v d 0 0:2 0:4 0:6 0:8 1 0 0:4 0:8 1:2 1:6 2 0 0:5 1 1:5 2 t _ N d =1000 v r =v d 0 0:2 0:4 0:6 0:8 1 0 0:4 0:8 1:2 1:6 2 0 0:5 1 1:5 2 t _ N d =1000 v r =v d 0 0:2 0:4 0:6 0:8 1 0 0:4 0:8 1:2 1:6 2 0:003 0 0:003 t _ N d =1000 v t =v d 0 0:2 0:4 0:6 0:8 1 0 0:4 0:8 1:2 1:6 2 0:003 0 0:003 t _ N d =1000 v t =v d 0 0:2 0:4 0:6 0:8 1 0 0:4 0:8 1:2 1:6 2 0:003 0 0:003 t _ N d =1000 v t =v d 0 0:2 0:4 0:6 0:8 1 Figure 4.14: Time history of f v (v z ) (top three panels), f v (v r ) (middle three panels), and f v (v t ) (bottom three panels) of case A 4 . The left, middle, and right columns show panels evaluated from region (1), (2), and (3), respectively. The values of f v have been normalized to range from 0 to 1. The time history of the velocity distribution function f v (v z ), f v (v r ), and f v (v t ) of the case A 4 are plotted in Fig. 4.14. Note that here we neglect the negative sign of v z . The top left panel shows thatf v (v z ) moves fromv z = 0 tov z =v d 23 as time goes by in the upstream region, indicating the acceleration of droplets along z direction. After t _ N d = 1000, f v (v z ) peaks at about v z =v d = 23. The top middle panel shows the central downstream region of f v (v z ), in which f v (v z ) starts to have non-zero values after about t _ N d = 1000 when droplets reach this region, and peaks at about v z =v d = 25. The top right panel shows the outer downstream region of f v (v z ), in which f v (v z ) peaks at about v z =v d = 25 after t _ N d = 1000 like the top middle panel. Chapter 4 49 The left panel in the middle row shows thatf v (v r ) moves fromv r = 0 tov r =v d 1:7 as time goes by in the upstream region, indicating the acceleration of droplets along the radial direction. f v (v r ) reaches a steady state after about t _ N d = 1000. The middle panel in the middle row shows that f v (v r ) starts to have non-zero values after t _ N d = 1000 and becomes stable soon in the central downstream region. The right panel in the middle row shows that f v (v r ) starts to have non-zero values after t _ N d = 1000, and f v (v r ) changes from peaking at 2v d to gradually spread as time goes by, because initially only those droplets having large v r can reach the outer downstream region, and then droplets with smaller v r move into this region. 0 0:2 0:4 0:6 0:8 1 0 5 10 15 20 25 30 f v (v z ) v z =v d l 1 l 2 l 3 0 0:2 0:4 0:6 0:8 1 0 0:5 1 1:5 2 f v (v r ) v r =v d l 1 l 2 l 3 0 0:2 0:4 0:6 0:8 1 0:003 0:002 0:001 0 0:001 0:002 0:003 f v (v t ) v t =v d l 1 l 2 l 3 Figure 4.15: Time averaged plots of f v (v z ) (top), f v (v r ) (middle), and f v (v t ) (bottom) of case A 4 . l 1 , l 2 , and l 3 indicate the three regions evaluated. The values of f v have been normalized to range from 0 to 1. 50 Chapter 4 The left bottom panel shows that f v (v t ) becomes more and more like a Maxwellian distribution as time goes by, and becomes steady aftert _ N d = 1000. The middle bottom and the right bottom panels are similar, showing a Maxwellian-like distribution but colder than that of the left bottom panel. For presenting more detail of f v , time averaged plots are shown in Fig. 4.15 using the data presented in Fig. 4.14. l 1 of f v (v z ) and f v (v r ) are obtained using the results between t _ N d = 1220 andt _ N d = 2000;l 1 off v (v t ) is obtained using the results betweent _ N d = 1020 and t _ N d = 2000;l 2 andl 3 off v (v z ) andf v (v r ) are obtained using the results betweent _ N d = 1820 and t _ N d = 2000; l 2 and l 3 of f v (v t ) are obtained using the results between t _ N d = 1620 and t _ N d = 2000. 0 0:2 0:4 0:6 0:8 1 0 3 6 9 12 15 18 21 24 27 30 f v (v z ) v z =v d A 4 A 3 A 2 A 1 A 0 0 0:2 0:4 0:6 0:8 1 0 3 6 9 12 15 18 21 24 27 30 f v (v z ) v z =v d A 4 A 3 A 2 A 1 A 0 0 0:2 0:4 0:6 0:8 1 0 3 6 9 12 15 18 21 24 27 30 f v (v z ) v z =v d A 4 A 3 A 2 A 1 A 0 Figure 4.16: Time averaged plots of f v (v z ) of all cases, evaluated in region (1), (2), and (3) shown in the top, middle, and bottom panel, respectively. The values of f v have been normalized to range from 0 to 1. Chapter 4 51 In the top panel of Fig. 4.15, the linel 1 includes those droplets that are in the upstream with relatively small v z , and those droplet that have been accelerated to about v z =v d = 23. Both linel 2 andl 3 are centered at aboutv z =v d = 25, but the peak value ofl 2 is less than that of l 3 , because the region (2) is a rectangular in the r xy -z plot, which includes more droplets with relatively small v z in the midstream region, while the region (3) is triangular in the r xy -z plot, which includes less droplets in the midstream region. 0 0:2 0:4 0:6 0:8 1 0 0:25 0:5 0:75 1 1:25 1:5 1:75 2 f v (v r ) v r =v d 0 0:2 0:4 0:6 0:8 1 0 0:25 0:5 0:75 1 1:25 1:5 1:75 2 f v (v r ) v r =v d 0 0:2 0:4 0:6 0:8 1 0 0:25 0:5 0:75 1 1:25 1:5 1:75 2 f v (v r ) v r =v d Figure 4.17: Time averaged plots off v (v r ) of all cases, evaluated in region (1), (2), and (3) shown in the top, middle, and bottom panel, respectively. Case A 0 to case A 4 corresponds to the line from the left to the right for all panels. The values off v have been normalized to range from 0 to 1. In the middle panel of Fig. 4.15,l 1 andl 2 are similar, butl 1 peaks at av r that is slightly greater than that of l 2 , since l 1 includes more droplets that have large enough v r to pass 52 Chapter 4 region (2) and enter region (3). l 3 peaks at a relatively large v r compared with l 1 and l 2 , because only those droplets having high v r can reach region (3). In the bottom panel of Fig. 4.15, all lines have a Maxwellian shape, l 1 is more thermal than l 2 and l 3 , and l 2 is almost the same as l 3 . This is because l 1 includes droplets in the upstream region where droplets are denser thus interact with each other more signicantly to spread v t . 0 0:2 0:4 0:6 0:8 1 0:01 0:005 0 0:005 0:01 f v (v t ) v t =v d A 4 A 3 A 2 A 1 A 0 0 0:2 0:4 0:6 0:8 1 0:01 0:005 0 0:005 0:01 f v (v t ) v t =v d A 4 A 3 A 2 A 1 A 0 0 0:2 0:4 0:6 0:8 1 0:01 0:005 0 0:005 0:01 f v (v t ) v t =v d A 4 A 3 A 2 A 1 A 0 Figure 4.18: Time averaged plots of f v (v t ) of all cases, evaluated in region (1), (2), and (3) shown in the top, middle, and bottom panel, respectively. The values of f v have been normalized to range from 0 to 1. Next we start to comparef v of dierent cases, using the same time averaged diagnoses. Fig. 4.16 shows f v (v z ) of all cases, evaluated in region (1), (2), and (3), shown in the top, Chapter 4 53 middle, and bottom panel, respectively. As can be seen, dierent cases produce similar f v (v z ), except A 0 is a little bit o due to the large statistical noise. Fig. 4.17 shows f v (v r ) of all cases, evaluated in region (1), (2), and (3), shown in the top, middle, and bottom panel, respectively. As can be seen, from case A 0 to A 4 , f v (v r ) shifts from small v r to large v r , as expected from Fig. 4.4 and Fig. 4.5. 0 0:4 0:8 1:2 1:6 2 0 4 8 12 16 20 t _ N d =1000 ez (m) 0 0:2 0:4 0:6 0:8 1 0 0:4 0:8 1:2 1:6 2 0 4 8 12 16 20 t _ N d =1000 ez (m) 0 0:2 0:4 0:6 0:8 1 0 0:4 0:8 1:2 1:6 2 0 4 8 12 16 20 t _ N d =1000 ez (m) 0 0:2 0:4 0:6 0:8 1 0 0:4 0:8 1:2 1:6 2 0 2 4 6 8 t _ N d =1000 er (m) 0 0:2 0:4 0:6 0:8 1 0 0:4 0:8 1:2 1:6 2 0 2 4 6 8 t _ N d =1000 er (m) 0 0:2 0:4 0:6 0:8 1 0 0:4 0:8 1:2 1:6 2 0 2 4 6 8 t _ N d =1000 er (m) 0 0:2 0:4 0:6 0:8 1 0 0:4 0:8 1:2 1:6 2 0 4 8 12 16 20 t _ N d =1000 et (m) 0 0:2 0:4 0:6 0:8 1 0 0:4 0:8 1:2 1:6 2 0 4 8 12 16 20 t _ N d =1000 et (m) 0 0:2 0:4 0:6 0:8 1 0 0:4 0:8 1:2 1:6 2 0 4 8 12 16 20 t _ N d =1000 et (m) 0 0:2 0:4 0:6 0:8 1 Figure 4.19: Time history of f p ( ez ) (top three panels), f p ( er ) (middle three panels), and f p ( et ) (bottom three panels) of case A 4 . The left, middle, and right columns show panels evaluated from region (1), (2), and (3), respectively. The values of f p have been normalized to range from 0 to 1. Fig. 4.18 shows f v (v t ) of all cases, evaluated in region (1), (2), and (3), shown in the top, middle, and bottom panel, respectively. As can be seen, from case A 0 to A 4 , f v (v t ) 54 Chapter 4 changes from wider to narrower, namely from more thermal to colder, indicating that the thermal speed of v t decreases as more droplets are broken. 4.2.5 PDF of Equivalent Length 0 0:2 0:4 0:6 0:8 1 0 2 4 6 8 10 12 14 16 18 20 f p ( ez ) ez (m) l 1 l 2 l 3 0 0:2 0:4 0:6 0:8 1 0 1 2 3 4 5 6 7 8 f p ( er ) er (m) 0 0:2 0:4 0:6 0:8 1 0 2 4 6 8 10 12 14 16 18 20 f p ( et ) et (m) Figure 4.20: Time averaged plots of f p ( ez ) (top), f p ( er ) (middle), and f p ( et ) (bottom) of case A 4 . l 1 , l 2 , and l 3 indicate the three regions evaluated. The values of f p have been normalized to range from 0 to 1. The time history of the probability density function of the equivalent length,f p ( ez ),f p ( er ), andf p ( et ), of the caseA 4 are plotted in Fig. 4.19. f p ( ez ) does not include the eects caused by the accelerating eld, but does include the long-range Coulomb eld due to droplets in the Chapter 4 55 z direction. The top left panel in Fig. 4.19 showsf p ( ez ) in the upstream region. We can see that, initially f p ( ez ) has values at very small ez , because droplets are closely located near the injection volume. Then,f p ( ez ) expands and gets a steady state after aboutt _ N d = 1200. The top middle and top right panels show f p ( ez ) of the central downstream region and the outer downstream, both of which start to have values after t _ N d = 1000. 0 0:2 0:4 0:6 0:8 1 0 2 4 6 8 10 12 14 16 18 20 f p ( ez ) ez (m) A 4 A 3 A 2 A 1 A 0 0 0:2 0:4 0:6 0:8 1 0 2 4 6 8 10 12 14 16 18 20 f p ( ez ) ez (m) A 4 A 3 A 2 A 1 A 0 0 0:2 0:4 0:6 0:8 1 0 2 4 6 8 10 12 14 16 18 20 f p ( ez ) ez (m) A 4 A 3 A 2 A 1 A 0 Figure 4.21: Time averaged plots of f p ( ez ) of all cases, evaluated in region (1), (2), and (3) shown in the top, middle, and bottom panel, respectively. The values of f p have been normalized to range from 0 to 1. f p ( er ) does not include the eects caused by the accelerating eld, but does include the long-range Coulomb force due to droplets in the radial direction. Note that the range of er plotted in the middle row of Fig. 4.19 is decreased to from 0 to 8m, which is dierent from those of ez and et , because f p ( er ) peaks at smaller er values, indicating that the radial 56 Chapter 4 droplet interactions are stronger. Overall, the patterns of the three panels in the middle row are similar to those in the top row, but er of region (2) spreads more. f p ( et ) only re ects the short-range Coulomb collisions. Compare the last row in Fig. 4.19 with the other two rows, we can see that f p ( et ) locates at larger values of et , due to the fact that the tangential interparticle force is relatively small, and f p ( et ) spreads much more than f p ( ez ) and f p ( er ). 0 0:2 0:4 0:6 0:8 1 0 1 2 3 4 5 6 7 8 f p ( er ) er (m) A 4 A 3 A 2 A 1 A 0 0 0:2 0:4 0:6 0:8 1 0 1 2 3 4 5 6 7 8 f p ( er ) er (m) A 4 A 3 A 2 A 1 A 0 0 0:2 0:4 0:6 0:8 1 0 1 2 3 4 5 6 7 8 f p ( er ) er (m) A 4 A 3 A 2 A 1 A 0 Figure 4.22: Time averaged plots of f p ( er ) of all cases, evaluated in region (1), (2), and (3) shown in the top, middle, and bottom panel, respectively. The values of f p have been normalized to range from 0 to 1. For presenting more detail of f p , time averaged plots are shown in Fig. 4.20 using the data presented in Fig. 4.19 following the same averaging method described in Sec. 4.2.4. Chapter 4 57 Fig. 4.20 (top) shows that the axial forces in the upstream region are much greater than those in the downstream region as mentioned before, and the values off p ( ez ) in the central downstream region and in the outer downstream region are similar. Fig. 4.20 (middle) shows that the radial forces in the upstream region are much greater than those in the downstream region, and l 2 has a portion at larger er than that of l 3 , because those droplets in the central region feel less radial forces, while those in the outer region feel stronger outward repulsive forces. Fig. 4.20 (bottom) shows that the tangential forces (short-range collisions) in the upstream region are much greater than those in the downstream region, andl 3 locates slightly right tol 2 , indicating that the collisions are more signicant in the outer downstream than in the central downstream. 0 0:2 0:4 0:6 0:8 1 0 2 4 6 8 10 12 14 16 18 20 f p ( et ) et (m) A 4 A 3 A 2 A 1 A 0 0 0:2 0:4 0:6 0:8 1 0 2 4 6 8 10 12 14 16 18 20 f p ( et ) et (m) A 4 A 3 A 2 A 1 A 0 0 0:2 0:4 0:6 0:8 1 0 2 4 6 8 10 12 14 16 18 20 f p ( et ) et (m) A 4 A 3 A 2 A 1 A 0 Figure 4.23: Time averaged plots of f p ( et ) of all cases, evaluated in region (1), (2), and (3) shown in the top, middle, and bottom panel, respectively. The values of f p have been normalized to range from 0 to 1. 58 Chapter 4 Next we start to comparef p of dierent cases, using the same time averaged diagnoses. Fig. 4.21 shows that for all regions, having more droplets results in stronger droplet interac- tions (including both the long-range and short-range) in thez direction, sincef p ( ez ) moves toward the left asN s increases. Similarly, Fig. 4.22 shows that having more droplets results in stronger droplet interactions (including both the long-range and short-range) in the radial direction, since f p ( er ) moves toward the left as N s increases. Again, Fig. 4.22 shows that having more droplets results in stronger droplet short-range collisions, since f p ( et ) moves toward the left as N s increases. Therefore, in terms of short-range Coulomb collisions, the beam is more collisional in the upstream region than in the downstream region, slightly more collisional in the central region than the outer region, and more collisional as the original droplets break into more sub-droplets. 4.3 From PP to PM 0 0:02 0:04 0:06 0:08 r xy =d 0 0 0:2 0:4 0:6 0:8 1 z=d 0 0 2 4 6 8 10 1 ○ 1.4527 2 ○ 5.6238 3 ○ 5.8566 (0:045; 0:5) Figure 4.24: The mean interparticle distance of the case A 4 with t _ N = 1=(2N s ) in the unit of m. In this section, we discuss what else we can obtain from the simulation results of f p . If we want to apply a PM model to simulate the colloid thruster beam without missing any physics due to both the long-range and short-range droplet interactions, the mesh size of the PM model should be small enough to include all interactions captured by f p . For example, if we want to develop a valid PM model to simulate the downstream beam region z=d 0 2 (0; 0:5) using the parameters of case A 4 , the PP results of f p suggest a suciently minimum mesh size. As we can see from Fig. 4.20 (lines of l 2 or l 3 ), f p ( ez ) and f p ( er ) drop to zero approximately when e < 0:9 m, and f p ( et ) drop to zero approximately when e < 1:5 m. Because f p ( et ) purely re ects the short-range collisions, and f p ( ez ) and f p ( er ) have Chapter 4 59 values at smaller e due to the addition of long-range interactions, the suciently minimum mesh size should be about 1.5 m. Meaning that if we let the mesh size of a PM model to be about 1.5m, it should be able to capture all the necessary short-range particle collisions in the downstream region of the beam. In addition, 1.5 m is approximately the suciently minimum mesh size for the cases A 0 , A 1 , A 2 , and A 3 as well, as can be seen from Fig. 4.23. One instructive length quantity is the mean interparticle distance r m = n 1=3 . We evaluated r m of the case A 4 from the droplet number density plot shown in Fig. 4.9, and plotted the spatial distribution of r m in Fig. 4.24. From Fig. 4.24, we obtained three averaged values ofr m in the three evaluated regions according to r m = P (nr m )= P n, where the sum is over the bins in both r xy and z directions of the plot in Fig. 4.24. These three values of r m are labeled in Fig. 4.24 and drawn in Fig. 4.20 as black vertical dashed lines for convenience. 2 4 6 8 10 12 0 1 2 3 4 5 6 N m = 1024 256 64 16 4 PP PP is faster PM is faster log 10 N o log 10 N p Figure 4.25: Operation comparison between PP and PM. 0 5 10 15 20 25 30 0 2 4 6 8 10 12 14 16 N o =10 8 N p =1000 PP PM PP/PM 1.66 Figure 4.26: An example of the operation comparison between PP and PM in the down- stream beam region. 60 Chapter 4 The concept of PM model was originally designed to speed up the simulation for systems with a large number of particles. A simple comparison of the computation between PP model and PM model is given here. According to the operation count Hockney derived in his book [39], a typical PP model has 10N 2 p N p operations, and a typical PM model has 20N p + 5N 3 m log 2 N 3 m operations, where N p denotes the number of particles and N 3 m denotes the number of mesh points. Based on these two formulas, Fig. 4.25 is plotted, in which N o denotes the number of operations. As can be seen, the PM model is more ecient in the region with small N m and large N p . In our case, the number of droplets is increasing as time goes by, and the PM model may not be more advantageous than PP if we apply it from the jet front, where the required mesh size is so small. However, it would be advantageous if we develop a hybrid PP and PM model, which applies PP on the upstream region of the beam and PM on the downstream region. For example, if we consider the downstream regionz=d 0 2 (0; 0:5) and use the parameters of the caseA 4 , a PM model with a cylindrical domain will haveN 3 m (100=1:5) 2 (500=1:5) 167 3 , if the mesh size is chosen to be the suciently minimum mesh size 1.5m. Using 10N 2 p N p for PP and 20N p + 5N 3 m log 2 N 3 m for PM to solve for operations as N p increases from 1 to 16000, we obtain Fig. 4.26. Integrate the curve with respect to N p , we can approximately obtain the number of cumulative operations, which represents the computation time. It is found that PM is about 1.66 times faster than PP under this consideration. 0:9995 1 1:0005 1:001 1:0015 1:002 1:0025 0 400 800 1200 1600 2000 E tot = E 0 t _ N d A 1 :v th =v d = 1:0 A 3 :v th =v d = 1:0 Figure 4.27: Time history of system energy for comparing dierent initial thermal velocities. Chapter 4 61 4.4 Initial Droplet Thermal Velocity 0 0:2 0:4 0:6 0:8 1 0 4 8 12 16 20 24 28 32 36 40 44 48 f p ( et ) et (m) A 1 :v th =v d = 0:0 A 1 :v th =v d = 0:1 A 1 :v th =v d = 1:0 0 0:2 0:4 0:6 0:8 1 0 4 8 12 16 20 24 28 32 36 40 44 48 f p ( et ) et (m) A 3 :v th =v d = 0:0 A 3 :v th =v d = 0:1 A 3 :v th =v d = 1:0 Figure 4.28: PDF of the equivalent length for comparing dierent initial thermal velocities. Here, we study the eects caused by the initial droplet macro-temperature, i.e. the initial thermal velocity. Since in the reality, when droplets are formed at the jet front, they may experience Coulomb ssion, which must result in thermal velocities in all directions. There- fore, here we add thermal velocities when injecting droplets. The injected droplet velocity distributions in x and y directions are chosen to be Maxwellian; in z direction, it is also chosen to be a Maxwellian distribution but shifted to be centered atv d . Two thermal velocities are considered,v th =v d = 0:1 and 1.0, for all directions and for the casesA 1 andA 3 with t _ N d = 1=(2N s ). First, the system energy is checked, shown in Fig. 4.27, in which only the two cases with v th =v d = 1:0 are plotted distinguishable, since they have relatively large ranges of E tot = E 0 . The system energy of the A 1 case with v th =v d = 1:0 conserves within about 0.25%; the system energy of the A 3 case with v th =v d = 1:0 conserves within about 0.03%. Therefore, energy conservation is satised for all cases. Second, the phase space plots ofx-z are plotted in Fig. 4.32. As can be seen, adding an initial thermal velocity leads to a more diusive beam in the radial direction. However, in the axial direction, it is found that (results not shown here) the eects of the initial thermal velocity is negligible for the chosen v th . 62 Chapter 4 Last, the PDF of the equivalent lengthf p ( et ) is compared, shown in Fig. 4.28. Here, we consider the whole downstream region of the beam, i.e. the region 2 plus the region 3 dened in Sec. 4.2.4, and we average the data using the same way as the one used in Sec. 4.2.5. Fig. 4.28 shows that f p ( et ) moves to the right for a larger v th , i.e. a higher T , thus considering an initial thermal velocity results in a more collisionless system, because an initial thermal velocity leads to a more diuse beam. Therefore, since the initial thermal velocity (or the initial droplet macro-temperature) aects f p , if we want to choose an appropriate mesh size for a PM model based on the results of f p , the appropriate initial thermal velocity should be taken into account. 4.5 Space Charge Eects To further study the space charge eects captured from F pp , we simulate the case A 0 with t _ N d = 1=2 without considering F pp , meaning that droplets only feel the accelerating eld forceF acc . Fig. 4.29 shows the x-z phase space plot of the two cases with and without F pp (left), and an amplied view of the case without F pp (right). We can see that, the beam expands radially about 1000 times less whenF pp is not considered. 0 0:2 0:4 0:6 0:8 1 0:03 0:015 0 0:015 0:03 z=d 0 x=d 0 0 0:2 0:4 0:6 0:8 1 0:1 0:05 0 0:05 0:1 z=d 0 x=(0:001d 0 ) Figure 4.29: Phase space x-z plot of the cases (A 0 ) with and without F pp (left), and an amplied view of the case withoutF pp (right). Chapter 4 63 0 0:1 0:2 0:3 0:4 0:5 0:6 0:7 0 0:2 0:4 0:6 0:8 1 v xy =v d z=d 0 0 0:4 0:8 1:2 1:6 2 0 0:2 0:4 0:6 0:8 1 v xy =(0:001v d ) z=d 0 Figure 4.30: Phase space z-v xy plot of the cases (A 0 ) with and without F pp (left), and an amplied view of the case withoutF pp (right). Fig. 4.30 shows thez-v xy phase space plot of the two cases with and withoutF pp (left), and an amplied view of the case without F pp (right). We can see that, the radial droplet velocity is also about 1000 times less whenF pp is not considered. 0 5 10 15 20 25 0 0:2 0:4 0:6 0:8 1 v z =v d z=d 0 Figure 4.31: Phase space z-v z plot of the cases (A 0 ) withoutF pp . 64 Chapter 4 0 0:2 0:4 0:6 0:8 1 -0.3 -0.15 0 0.15 0.3 z=d 0 x=d 0 v th =v d = 0:0 0 0:2 0:4 0:6 0:8 1 -0.3 -0.15 0 0.15 0.3 z=d 0 x=d 0 v th =v d = 0:1 0 0:2 0:4 0:6 0:8 1 -0.3 -0.15 0 0.15 0.3 z=d 0 x=d 0 v th =v d = 1:0 0 0:2 0:4 0:6 0:8 1 -0.3 -0.15 0 0.15 0.3 z=d 0 x=d 0 v th =v d = 0:0 0 0:2 0:4 0:6 0:8 1 -0.3 -0.15 0 0.15 0.3 z=d 0 x=d 0 v th =v d = 0:1 0 0:2 0:4 0:6 0:8 1 -0.3 -0.15 0 0.15 0.3 z=d 0 x=d 0 v th =v d = 1:0 Figure 4.32: Phase space plots of x-z for comparing dierent initial thermal velocities. A 1 cases are shown in the rst row; A 3 cases are shown in the second row. Fig. 4.31 shows the z-v z phase space plot of the case withoutF pp . If one compare this gure to the result obtained with F pp shown in Fig. 4.3, the dierence between the two curves is negligible. Therefore, the radial space charge eect is signicant and the beam radial expansion is mainly due to F pp rather than F acc . The axial space charge eect is negligible, since the axial component of F pp is much less than the axial component of F acc . 4.5.1 Calculation Based on Child-Langmuir Law Let us rst use the Child-Langmuir law to estimate the space-charge-limited current density for the case A 0 studied in the previous chapter. The droplet charge is Q d = 3:77 10 18 C, mass is M d = 8:15 10 22 kg. The accelerating potential is 0 = 5000 V, the tip-to-plate distance is d 0 = 1 mm. Therefore, J SCL = 4 9 " 0 r 2Q d M d 3=2 0 d 2 0 133:8 A/m 2 : However, our emitted current is I = 36:58 nA, the jet radius is R j = 2:94 nm, thus the emitted current density is J = I=(R 2 j ) 1:344 10 9 A/m 2 , which is seven orders of magnitude larger thanJ SCL , but we still can obtain an accelerated owing beam of droplets with negligible space charge eects in the simulation. Because we have already veried the PP model, the only possible answer to this con ict is that the Child-Langmuir law is not Chapter 4 65 appropriate to be applied for the colloid thruster droplet emission, at least within the range of the parameters we simulated. To understand the reasons that cause the inappropriateness of the Child-Langmuir law, we rst do some calculations shown in the following section. 4.5.2 Model of Two Parallel Plates First, let us consider the 1D scenario of the Child-Langmuir law. The accelerating potential is 0 . The distance between the two plates is d 0 . So, the uniform constant electric eld between the two plates is E 0 = 0 =d 0 . The original droplet charge is Q d , mass is M d , and radius is R d . The droplet emitting rate is _ N d , so every t 0 = 1= _ N d (seconds), an original droplet should be released. Let us use the parameters of the case A 0 : 0 = 5000 V and d 0 = 0:001 m) E 0 = 0 =d 0 = 5 10 6 V/m. Q d = 3:769 10 18 C, M d = 8:149 10 22 kg, R d = 5:563 10 9 m. _ N d = 9:707 10 9 s 1 , t 0 = 1= _ N d = 1:030 10 10 s. Now, let the emission start. At time t = 0, the rst original droplet D 1 is emitted at z 1 (t = 0) = 0 m with velocity v 1 (t = 0) = 0 m/s. D 1 is then accelerated by E 0 = 5 10 6 V/m. At time t = t 0 , v 1 (t = t 0 ) = v 1 (t = 0) + Q d E 0 M d t 0 2:382 m/s z 1 (t = t 0 ) = z 1 (t = 0) +v 1 (t = 0)t 0 + 1 2 Q d E 0 M d (t 0 ) 2 1:227 10 10 m Still at time t = t 0 , the second original droplet D 2 is emitted at z 2 (t = t 0 ) = 0 m, with velocity v 2 (t = t 0 ) = 0 m/s. The Coulomb electric eld on D 2 due to D 1 is E C21 = 1 4" 0 Q d [z 1 (t = t 0 )z 2 (t = t 0 )] 2 2:249 10 12 V/m So, E C21 E 0 , and E C21 is pointing at the opposite direction of E 0 . Thus the second original droplet D 2 will not be accelerated toward the desired direction. Also, note that jz 1 (t = t 0 )z 2 (t = t 0 )j< 2R d , so the two droplets should actually be recombined. Therefore, if we consider this scenario, the space charge limit has already been reached right after the rst droplet is emitted. 66 Chapter 4 4.5.3 Model of Taylor Cone Rather than a uniform constant electric eld E 0 , the eld used in the simulation is not uniform, it has a larger magnitude near the Taylor cone tip (or the jet front). The electric eld at the Taylor cone tip is about E T = 8:506 10 7 V/m. Then, we will have z 1 (t = t 0 ) = 1 2 Q d E T M d (t 0 ) 2 2:087 10 9 m E C21 = 1 4" 0 Q d [z 1 (t = t 0 )] 2 7:776 10 9 V/m Again, we obtain E C21 >E T . However, the original droplets in the simulation have an initial axial velocity, which equals the ow speed of the jet, v d = 257:1 m/s. Considering v 1 (t = 0) =v d , we can obtain z 1 (t = t 0 ) = v d t 0 + 1 2 Q d E T M d (t 0 ) 2 2:858 10 8 m E C21 = 1 4" 0 Q d [z 1 (t = t 0 )] 2 4:145 10 7 V/m Now we have E C21 <E T , so the second original droplet can be accelerated toward the desired direction. Further more, in the simulation, the original droplets are not emitted at the same point (0; 0; 0). Instead, they are randomly distributed in the injection volume, which has a radius R j = 2:943 10 9 m, and a height H = v d t 0 = 2:650 10 8 m, to simulate the random formation of the original droplets from the jet. Thus, it further increases the possible distances between the new droplet being injected and previous droplets that have been injected. Therefore, nearly all droplets in the simulation can be accelerated as desired. In addition, for the beam of case A 0 in the nal steady state. The electric eld due to all droplets at the location of the Taylor cone tip is about 4:1 10 7 V/m, which is about half ofE T = 8:506 10 7 V/m. Thus we can further see that the space charge eect has not been reached. If we compare the nal beam size with the jet radius, we can see that the beam expands a lot radially. The nal beam radius at z = 0 is about 0.04 mm for case A 0 . The jet radius is R j = 2:94 nm = 0:00000294 mm. Chapter 4 67 4.5.4 The 1D Assumption of Child-Langmuir law Another obvious reason is that Child-Langmuir law is derived based on the one-dimension assumption, namely the two parallel plates are innitely large and there is an innite number of charged particles in the radial direction (o-axis). In terms of the space charge eects, a particle that is emitted on the axis will feel a signicant amount of repulsive forces that are due to those previously injected particles that are o-axis. In our case, however, we only have a thin beam, these repulsive forces thus do not exist. In addition, although the beam mainly ows along the axial direction, it also expands in the radial direction. An expanded beam front will have a much less impact on the newly injected particles at the beam end. 0:2 0:15 0:1 0:05 0 0 0.2 0.4 0.6 0.8 1.0 E z (z) (mV/m) z (m) Child n = 5 n = 0 Figure 4.33: The electric eld distribution. One way to evaluate the dierence between 1D and 3D is to apply the image charges in the 3D PP code. When calculating the Coulomb eld, some image charges in the x and y directions are also involved. The electric eld at the location of the ith particle E (i) pp can be evaluated from E (i) pp = 1 4" 0 X n N X j=1 0q j (r i r j +nL) jr i r j +nLj 3 ; (4.1) where L is the domain length in the x and y direction, nL = n 1 L^ x +n 2 L^ y denotes an arbitrary repeat vector,n 1 andn 2 are arbitrary integers, the 0 symbol is to exclude the term j =i if and only ifn = 0. This calculation refers to the basics of Ewald summation [30, 82]. The parameters for the PP model are chosen to be the same as those described in Sec. 3.2.3, except the potential which is reduced to be 0.1 mV, so the macro-particle weight is reduced to N W 2:46. And, the full Coulomb eld is used without cuto. Since we are now simulating nearly real particles, we use the PP algorithm to evaluate the electric eld, 68 Chapter 4 rather than solving the Poisson's equation on mesh points. We evaluate the electric eld E z along the z-axis, on 100 uniformly distributed points from z = 0 to z = 1 m. The electric eld distribution of two runs are plotted in Fig. 4.33. The curve of n = 0 is a run without image particles, while the curve of n = 5 is a run with ve image particles in each of the x,x, y, andy directions. Namely, for every simulated particle, there are 10 10 corresponding image particles. From Fig. 4.33, we can see that a larger n leads to a E z distribution that is closer to the Child-Langmuir law. When n = 0, however, i.e. the 3D scenario, the space charge eects are much less than that of the 1D scenario. 4.5.5 Conclusion for the Space Charge Eects The conclusions are as follows: The PP algorithm and code are veried through the com- parison with the Child-Langmuir law and the PIC method. The Child-Langmuir law is inappropriate to be applied for the scenario of colloid thrusters at least within the parame- ters we simulated. The reasons that the Child-Langmuir law is inappropriate are as follows: (1) In our case, the applied electric eld is not a constant as the one in the Child- Langmuir law produced from two parallel plates. Instead, it has a much larger magnitude near the Taylor cone tip, and drops signicantly toward the downstream region. (2) Our droplets have initial velocities that are determined by the jet speed, which is dierent from the zero initial velocity assumption in the Child-Langmuir law. (3) We randomly inject droplets in an injection volume to simulate the real breakup of droplets from the jet front, which increases the distances between those particles that are just injected and those that are about to be injected. (4) The Child-Langmuir law assumes a 1D scenario, whose particles shall be innite large charged sheets. In our case, we simulate a real thin beam of real charged droplets along the beam axis. Thus, there is not an innite number of other charged droplets that are located o-axis and produce repulsive forces on those on-axis real droplets. (5) The emitted droplets under Coulomb's law expand radially signicantly. An ex- panded beam front will have a much less impact on the newly injected particles at the beam end. Chapter 5 69 Chapter 5 Colloid Thruster Neutralization of Droplet Beam Sometimes, an alternating current mode may be applied to generate an electrospray. The main advantage of the AC mode is that there will be no charge accumulation on the space- craft, compared to the DC mode, thus we do not need an extra neutralizer placed outside the thruster to neutralize the beam. Depending on the ratio of the characteristic time of the electrical actuation to the characteristic time of the dynamical response of the electri- ed liquid surface, one may distinguish two limiting regimes [33]: drop-on-demand (DoD) or forced dripping, when the ratio is large, and AC electrospray, when the ratio is small. For the application of colloid thrusters, the ratio is required to be large, since the goal is to generate thrust. For the DoD electrospray, the behavior of the liquid surface has been studied previously [17, 62, 63, 71, 72], but the droplet behavior during the electric potential transition of a colloid thruster has not been studied. In this chapter, the detailed droplet dynamics during the electric potential transition of the AC mode of colloid thruster is studied using the PP model. Two AC modes are considered, one is a square wave function, the other is a step function. Since both positively charged droplets and negatively charged droplets are contained in the beam, the recombi- nation between droplets is considered. More details of the simulation setup and procedures are given in the following section. After the AC mode, a bipolar system is studied, as introduced in 1.4.2, which contains two thrusters operating in opposite polarity to achieve the charge neutralization on the spacecraft. Here, we focus on the beams of the bipolar system, and study the space charge eects due to the two charged beams. 70 Chapter 5 5.1 AC Mode Simulation Setup and Procedures 10 8 6 4 2 0 2 0:30:20:1 0 0:1 0:2 0:3 z=d 0 x=d 0 10 8 6 4 2 0 2 0:30:20:1 0 0:1 0:2 0:3 z=d 0 x=d 0 10 8 6 4 2 0 2 0:30:20:1 0 0:1 0:2 0:3 z=d 0 x=d 0 10 8 6 4 2 0 2 0:30:20:1 0 0:1 0:2 0:3 z=d 0 x=d 0 Figure 5.1: Phase space x-z plots of the AC case: all droplets (top left); only recombined droplets (bottom left). The corresponding DC case: all droplets (top right); only recombined droplets (bottom right). Positive droplets are denoted by red dots; negative droplets are denoted by blue dots. The plotting simulation time is t _ N d = 16000. Two types of AC mode are considered: Type (1) is a square wave, namely from timet = 0 to t p =2, the accelerating voltage is set to be positive 0 , positively charged droplets are injected; from timet =t p =2 tot p , the accelerating voltage is set to be negative 0 , negatively charged droplets are injected, wheret p denotes the alternating period. Then, the accelerating voltage Chapter 5 71 alternates in the same manner. Type (2) is a step function, namely the accelerating voltage changes only once at a specic time that will be mentioned later. Since both positively charged and negatively charged droplets are involved in the system, recombinations should be considered. The recombination of any two droplets is considered as follows. If droplet 1 and droplet 2, denoted using subscripts 1 and 2, are located closer than the sum of their radii, i.e.jr 12 j<R 1 +R 2 , these two droplets recombine into a new droplet. This new droplet has massm =m 1 +m 2 , chargeq =q 1 +q 2 , radiusR = (R 3 1 +R 3 2 ) 1=3 , and is located atr = (r 1 +r 2 )=2. For the velocity of this new droplet, we only consider momentum conservation, thusv = (m 1 v 1 +m 2 v 2 )=m. Note that if the new droplet is over the Rayleigh limit [66],q> 8(" 0 ) 1=2 R 3=2 , it is not stable and should break up again, thus we do not let these two droplets recombine. Previously in Sec. 2.3, the height of the injection volume is set to be H = v d t. If a t is chosen such that H < 2R d , those droplets that are just injected should recombine immediately, which is not what we desire. Therefore, here we x H = v d = _ N d for all t to avoid the recombination occurred during injection, where 1= _ N d represents the time duration in which one original droplet should be injected. 5.1.1 Type 1: Square Wave 30 20 10 0 10 20 12108 6 4 2 0 2 v z =v d z=d 0 30 20 10 0 10 20 12108 6 4 2 0 2 v z =v d z=d 0 Figure 5.2: Left: Phase spacez-v z plot of the AC case. Right: The corresponding DC case. Positive droplets are denoted by red dots; negative droplets are denoted by blue dots. The parameters of the case A 0 listed in Sec. 4.1 are used here, namely we only consider the original droplets without any breakups. The alternating period is chosen to be t p = 72 Chapter 5 4000 _ N d 0:412 s. Note that t p =2 is about the droplet travel time from the injection volume to the accelerating plate. The simulation time step is chosen to be t _ N d = 1/2. 0 5 10 15 20 25 0 5000 10000 15000 20000 j v z j=v d t _ N d AC DC Figure 5.3: Time history of droplet exit speed of the square wave case (averaged in the region from0:1d 0 to 0:1d 0 ). 0 0:5 1 1:5 2 2:5 3 3:5 4 0 5000 10000 15000 20000 E tot =10 10 J t _ N d AC DC Figure 5.4: Time history of the total system energy. First, the x-z phase space plots are shown in Fig. 5.1. For the AC case, we can see that, from the top left panel, positive and negative droplets are mixed to dierent extent at Chapter 5 73 dierent locations of the beam. The beam shape remains the same when compared with the DC case, shown in the top right panel. We also plotted those droplets which have electric charges over the original value, i.e.jqj>Q d , in the bottom two panels. We found that those droplets occur in both AC and DC cases, and are all formed near the injection volume. This is due to the fact that the droplet number density is the highest near the injection volume, and the random injection mechanism applied. It is found that 16000 original droplets have been injected since the beginning of the simulation tillt _ N d = 16000. 816 droplets are vanished due to ying back to the Taylor cone. 983 droplets are vanished due to recombination. Therefore, there are 14201 droplets left at that time, out of which the number of original droplets is 13512. The numbers of those droplets having mass M d to 21M d are 13512, 659, 2, 3, 7, 3, 5, 3, 0, 1, 1, 1, 2, 0, 0, 0, 1, 0, 0, 0, and 1, respectively. The numbers of those droplets having charge Q d to 6Q d are 6702, 319, 5, 4, 2, and 1, respectively. The numbers of those droplets having chargeQ d to6Q d are 6813, 336, 4, 2, 5, and 0, respectively. The numbers of those droplets having charge 0 are 8. Therefore, we found that recombinations between like charges occur but only due to the simulation injection mechanism applied, and recombinations between unlike charges seldom occur. As expected, thez-v z plot of the AC case, Fig. 5.2 (left), shows oscillations ofv z , rather than a constant v z as the DC case, Fig. 5.2 (right). An averaged droplet exit speed,j v z j, is shown in Fig. 5.3, which indicates that the produced thrust from the AC case oscillates around 15v d _ V . 0:30:20:1 0 0:1 0:2 0:3 x=d 0 9 7 5 3 1 1 z=d 0 0:6 0:4 0:2 0 0:2 0:4 0:6 0:30:20:1 0 0:1 0:2 0:3 x=d 0 9 7 5 3 1 1 z=d 0 0:2 0:3 0:4 0:5 0:6 0:7 0:8 Figure 5.5: Electric potential (in volts) caused by charged droplets. The AC case is shown on the left panel, while the corresponding DC case is shown on the right panel. Under this consideration, the time history of the total system energy E tot is plotted in 74 Chapter 5 Fig. 5.4. Note that the averaged energy E tot diagnosed before in Sec. 4.2.2 is not used here, since the total number of droplets keeps changing during the simulation. We can see from Fig. 5.4 that E tot of the DC case increases linearly as droplets are continuously injected into the system with a constant rate. However, due to the alternation of the accelerating potential and the droplet charge,E tot of the AC case experiences a drop at every alternation. The electric potential caused only by charged droplets is shown in Fig. 5.5, where the potential at innity is assumed to be zero. The AC case is shown on the left panel, while the corresponding DC case is shown on the right panel. As can be seen, the dierence between AC and DC is obvious. 5.1.2 Type 2: Step Function One meaningful study is to nd the transition time t x after which the thrust changes back to its stable value after one alternation. Therefore, two other AC cases are tested using the step function, which are the same as previous AC cases in Sec. 5.1.1, except that AC1 stopped alternating after t _ N d = 10000, AC2 only alternates once at t _ N d = 10000. The results of the averaged droplet exit speed are shown in Fig. 5.6. For both AC1 and AC2 cases, the transition time is about 6000= _ N d , if we count at the time of the start of the last/one alternation. Since the type of AC2 is closer to the reality, we will focus on using this type to study the transition time. 0 5 10 15 20 25 0 5000 10000 15000 20000 j v z j=v d t _ N d AC1 AC2 DC Figure 5.6: Time history of droplet exit speed of the step function case (averaged in the region from0:1d 0 to 0:1d 0 ). Chapter 5 75 0 5 10 15 20 25 0 5000 10000 15000 20000 j v z j=v d t _ N d 5000 7000 9000 11000 13000 15000 1 2 4 8 16 32 t x _ N d _ V= _ V 0 t x / _ V 0:27 _ V 1=4 Figure 5.7: Left panel: Time history of droplet exit speed (averaged in the region from 0:1d 0 to 0:1d 0 ) for cases with varied volume ow rates. The curves from the top to the bottom correspond to an increasing the volume ow rate. Right panel: Dependence of the transition time on the volume ow rate. 0 5 10 15 20 25 0 5000 10000 15000 20000 j v z j=v d t _ N d 5500 6000 6500 7000 7500 8000 8500 9000 9500 2000 3000 4000 5000 t x _ N d 0 (V) t x / 0:52 0 1=2 0 Figure 5.8: Left panel: Time history of droplet exit speed (averaged in the region from 0:1d 0 to 0:1d 0 ) for cases with varied accelerating potentials. The curves from the top to the bottom correspond to an decreasing accelerating potential. Right panel: Dependence of the transition time on the accelerating potential. Next, we study the dependence of the transition time on the volume ow rate _ V , the 76 Chapter 5 accelerating potential 0 , and the accelerating plate spacing d 0 . Based on the case A 0 , we alternate the potential once at the time when 1000 droplets pass the accelerating plate. Fig. 5.7 shows the results of the averaged droplet exit speed obtained from cases with 1, 2, 4, 8, 16, and 32 times of the volume ow rate of the case A 0 (with _ V = 7 10 15 m 3 s 1 ). We can see that a lager _ V leads to a longer transition time. If we do a power t, we obtain t x / _ V 1=4 . Fig. 5.8 shows the results of the averaged droplet exit speed obtained from cases with dierent accelerating potentials 5000, 4500, 4000, 3500, 3000, 2500, and 2000 V. We can see that a higher 0 leads to a shorter transition time. If we do a power t, we obtaint x / 1=2 0 . Fig. 5.9 shows the results of the averaged droplet exit speed obtained from cases with dierent accelerating plate spacings from 0.5 mm to 15 mm. We can see that a lagerd 0 leads to a longer transition time. If we do a power t, we obtain t x /d 0 . 0 5 10 15 20 25 0 5000 10000 15000 20000 j v z j=v d t _ N d 2000 3000 4000 5000 6000 7000 8000 9000 4 6 8 10 12 14 16 t x _ N d d 0 =10 4 m t x /d 0:994 0 d 0 Figure 5.9: Left panel: Time history of droplet exit speed (averaged in the region from 0:1d 0 to 0:1d 0 ) for cases with varied d 0 . The curves from the left to the right correspond to an increasing d 0 . Right panel: Dependence of the transition time on d 0 . The reasoning of the above relations is as follows. Since the accelerating potential drops signicantly only near the region of the injection volume, we can ignore the acceleration pro- cess and assume all droplets gain their nal kinetic energy immediately after being injected. Based on this assumption, we can say that the transition time should be proportional to the time that one droplet travels from the injection volume to the accelerating plate. Note that we ignore the eects of the space charges, since it is much smaller compared to the eects of the accelerating eld. Therefore, from the energy point of viewM d v 2 =2 =Q d 0 , we can obtain the nal speed Chapter 5 77 v = (2Q d 0 =M d ) 1=2 , thus the transition time t x / d 0 =v = d 0 =(2Q d 0 =M d ) 1=2 . From the equations given in Sec. 2.2, we can obtain the relations: Q d / I= _ N d , M d / V d , I/ _ V 1=2 , V d /R 3 d / (r ) 3 / _ V , _ N d / _ V=V d . Thus, we found thatt x /d 0 1=2 0 _ V 1=4 , which agrees with the simulation results. 5.2 Bipolar Operation Mode As introduced in section 1.4, another way to avoid the charge accumulation on a spacecraft is to apply two thrusters side by side with dual polarity, which can be called a bipolar system. Here, we simulate a bipolar system with each thruster has only one emitter/nee- dle. One thruster is operated using a positive accelerating potential and emits positively charged droplets, while the other is operated using a negative accelerating potential and emits negatively charged droplets. In the previous chapter, we chose some reasonable parameters to test the particle-particle model for the droplet acceleration process. Here, we rst study in more detail about the operation parameters. 5.2.1 Parameter Calculation of Limitations Based on the case A 0 used in the previous chapter, we calculate the following parameters. (a) The starting voltage start = r R c " 0 ln 4d R c 1928 (V); (5.1) which indicates the minimum reasonable accelerating potential. (b) The Rayleigh limit q R = 8 p " 0 R 3=2 d 7:537 10 18 (C); (5.2) and we chose Q d = 3:769 10 18 C of the case A 0 to be half of q R , to mimic the maxi- mum charge observed in experiments which is about half of the Rayleigh limit. For those subdroplets broken from the original droplet in A 0 , their charges are more away from their corresponding Rayleigh limit. For example, the Rayleigh limit charge of the case A 4 is 1:884 10 18 C, and the charge of the subdroplets in the case A 4 is 2:355 10 19 C, which is 0.125 times of the the Rayleigh limit. (c) The minimum volume ow rate If we consider the ion emission criterion, and choose the electric eld to be 1.5 V/nm, 78 Chapter 5 we can obtain a minimum volume ow rate from E = 1:5 10 9 (V=m) = 1:87 10 7 [f(")] 1=3 1=2 K " _ V min 1=6 ) _ V min = 2:498 10 15 (m 3 /s): (5.3) Thus our volume ow rate of the case A 0 = 7:000 10 15 m 3 /s is higher than this _ V min . If we consider the Taylor cone instability, and choose the parameter min = 0:5, we can obtain another minimum volume ow rate, _ V min = "" 0 K 2 min 1:156 10 14 (m 3 /s): (5.4) However, the volume ow rate of the case A 0 is lower than _ V min . Therefore, for further simulations, it is more reasonable to only only increase the volume ow rate. Change of the Space Charge Eects If we believe the scale of the Child-Langmuir law still applies, and if we want to enhance the space charge eects, we may want to reduce J SCL . For xed 0 and d 0 , we can reduce Q d =M d . To reduce Q d =M d , we can only increase _ V , because Q d M d = f(") K " _ V 1=2 For xedQ d =M d andd 0 , we can reduce 0 , and make sure 0 is greater than the starting voltage start . For xed Q d =M d and 0 , we can increase d 0 , however, increasing d 0 leads to a longer simulation time. 5.2.2 Simulations of Bipolar System With respect to the simulation, we still choose the caseA 0 used in Sec. 4.1 with t _ N d = 1=2. The separation between the two thrusters is rst set to bex=d 0 = 0:4 only in thex direction. Two analytical accelerating elds are set up for two thrusters with opposite polarity. Thex-z phase space plot of all droplets at time step 38000 is shown in Fig. 5.10 (left), we can see that the two beams ow away without a signicant attraction from each other. The recombined droplets are shown in Fig. 5.10 (right). It is found again, the recombination only occurs near the injection volume of each thruster, and considering or not the recombination in the simulation does not aect much the results. The y-z phase space plot of all droplets Chapter 5 79 is shown in Fig. 5.11 (left), and thex-y phase space plot of all droplets is shown in Fig. 5.11 (right). 14 12 10 8 6 4 2 0 2 -0.5 0 0.5 1.0 z=d 0 x=d 0 14 12 10 8 6 4 2 0 2 -0.5 0 0.5 1.0 z=d 0 x=d 0 Figure 5.10: The x-z phase space plot with separation x=d 0 = 0:4 of all droplets (left) and only recombined droplets (right). Red dots represent positive droplets. Blue dots represent negative droplets. 14 12 10 8 6 4 2 0 2 -0.5 -0.25 0 0.25 0.5 z=d 0 y=d 0 -0.5 0 0.5 -0.5 0 0.5 1.0 y=d 0 x=d 0 Figure 5.11: The y-z phase space plot (left) and the x-y phase space plot (right) with separation x=d 0 = 0:4 of all droplets. Red dots represent positive droplets. Blue dots represent negative droplets. 80 Chapter 5 -0.5 0 0.5 1 -10 -5 0 -1 -0.5 0 0.5 1 Figure 5.12: The space charge potential of the beams with separation x=d 0 = 0:4. -2 -1 0 1 2 -10 -5 0 + - Figure 5.13: The average velocity of v x along the z direction with separation x=d 0 = 0:4. The space charge potential of the two beams is shown in Fig. 5.12. We can see that again, the highest potential is only about 1 volt. To see the extent of the attraction of the two beams, we analyze the average droplet velocity of v x along the z direction, which is shown in Fig. 5.13. We can see that positive droplets are slightly accelerated toward the x direction, and negative droplets are slightly accelerated toward thex direction, andj v x j reaches about 1 m/s at around z=d 0 = 12. Chapter 5 81 14 12 10 8 6 4 2 0 2 -0.4 -0.2 0 0.2 0.4 z=d 0 x=d 0 Figure 5.14: The x-z phase space plot with separation x=d 0 = 0:1 of all droplets. Red dots represent positive droplets. Blue dots represent negative droplets. -0.4 -0.2 0 0.2 0.4 -10 -5 0 -1 -0.5 0 0.5 1 Figure 5.15: The space charge potential of the beams with separation x=d 0 = 0:1. If we calculate the electric eld due to the negative beam at the location of the Taylor cone tip of the positive beam, i.e. (0; 0; 1) mm, we can see that it equals about 169.4 V/m in thex direction, which is much smaller that the axial accelerating eld at the Taylor cone tip, which is about E T = 8:506 10 7 V/m. Thus the attraction between the two beams is not signicant. 82 Chapter 5 Reduce the Separation Distance Then, a smaller separation is tested, x=d 0 = 0:1. We can see from the phase space x-z plot at time step 38000 in Fig. 5.14, the two beams start to overlap in the downstream beam region. The potential is plotted in Fig. 5.15, we can see that the overlapping region has an approximately zero potential. -4 -2 0 2 4 -10 -5 0 + - Figure 5.16: The average velocity of v x along the z direction with separation x=d 0 = 0:1. 6 5 4 3 2 1 0 1 2 -0.5 0 0.5 1.0 z=d 0 x=d 0 Figure 5.17: The x-z phase space plot with separation x=d 0 = 0:4 at time step 31000 and volume ow rate 10 _ V . Red dots represent positive droplets. Blue dots represent negative droplets. The average velocity is plotted in Fig. 5.16, it shows a larger value ofj v x j about 3 m/s Chapter 5 83 atz=d 0 = 12 as expected, since the attraction between the two beams is more signicant for a smaller separation. 3 2 1 0 1 2 -0.5 0 0.5 1.0 z=d 0 x=d 0 Figure 5.18: The x-z phase space plot with separation x=d 0 = 0:4 at time step 27000 and volume ow rate 100 _ V . Red dots represent positive droplets. Blue dots represent negative droplets. -0.5 0 0.5 1 -6 -4 -2 0 2 -3 -2 -1 0 1 2 3 Figure 5.19: The space charge potential of the beams with separation x=d 0 = 0:4 and volume ow rate 10 _ V . Enhance the Space Charge Eects (a) Varied volume ow rate 84 Chapter 5 To enhance the eect of space charge, we increase the volume ow rate to be 10 and 100 times of that of the case A 0 , and the separation between the two thrusters is still kept at x=d 0 = 0:4. The results of the x-z phase space are plotted in Fig. 5.17 and Fig. 5.18, respectively. And the potential plots are shown in Fig. 5.19 and Fig. 5.20. We can see from the potential plots that increasing the volume ow rate leads to the increase of the potential due to space charge, because a higher volume ow rate results in a bigger droplet size, thus a larger droplet charge. -0.5 0 0.5 1 -3 -2 -1 0 1 2 -15 -10 -5 0 5 10 15 Figure 5.20: The space charge potential of the beams with separation x=d 0 = 0:4 and volume ow rate 100 _ V . -2 -1 0 1 2 -6 -4 -2 0 2 + - Figure 5.21: The average velocity of v x along the z direction with separation x=d 0 = 0:4 and volume ow rate 10 _ V . The average velocity of v x are plotted in Fig. 5.21 and Fig. 5.22. They only show a slight increase of v x compared to the v x shown in Fig. 5.13. Chapter 5 85 -2 -1 0 1 2 -3 -2 -1 0 1 2 + - Figure 5.22: The average velocity of v x along the z direction with separation x=d 0 = 0:4 and volume ow rate 100 _ V . (b) Varied accelerating potential 8 7 6 5 4 3 2 1 0 1 -0.25 0 0.25 0.5 0.75 z=d 0 x=d 0 Figure 5.23: The x-z phase space plot with separation x=d 0 = 0:4 of all droplets at time step 26000. From the largest beam to the smallest beam, 0 changes from 2000 V, 3000 V, to 4000 V. The accelerating potential 0 of values 2000 V, 3000 V, and 4000 V are tested. The x-z phase space plots are shown in Fig. 5.23. We can see that a smaller 0 leads to a wider beam, since droplets are accelerated to a smaller nal axial speed, thus they have more time to repel each other to form a more radially expanded beam. 86 Chapter 5 -0.2 0 0.2 0.4 0.6 -6 -4 -2 0 2 -1 -0.5 0 0.5 1 -2 -1 0 1 2 -8 -6 -4 -2 0 2 + - Figure 5.24: Left: The space charge potential of the beams with separation x=d 0 = 0:4 and 0 = 2000 V at time step 26000. Right: The average velocity of v x . -0.2 0 0.2 0.4 0.6 -6 -4 -2 0 2 -1 -0.5 0 0.5 1 -2 -1 0 1 2 -8 -6 -4 -2 0 2 + - Figure 5.25: Left: The space charge potential of the beams with separation x=d 0 = 0:4 and 0 = 3000 V at time step 26000. Right: The average velocity of v x . The potential and the averagev x plots for the cases with the three dierent 0 are shown in Fig. 5.24, Fig. 5.25, and Fig. 5.26, respectively. We can see that as 0 decreases, the space charge potential slightly increases, but the dierence of the average v x is not obvious. Chapter 5 87 -0.2 0 0.2 0.4 0.6 -6 -4 -2 0 2 -1 -0.5 0 0.5 1 -2 -1 0 1 2 -8 -6 -4 -2 0 2 + - Figure 5.26: Left: The space charge potential of the beams with separation x=d 0 = 0:4 and 0 = 4000 V at time step 26000. Right: The average velocity of v x . (c) Varied tip-plate distance The tip-plate distance d 0 is increased to be 2 mm, 4 mm, and 8 mm. The x-z phase space plots are shown in Fig. 5.27. We can see that the beam expand more radially when d 0 is larger, since the acceleration process is longer in time and distance, droplets can have a chance to repel each other more radially. 6 4 2 0 2 4 6 8 -0.25 0 0.25 0.5 0.75 z=d 0 x=d 0 Figure 5.27: The x-z phase space plot with separation x=d 0 = 0:4 of all droplets at time step 26000 with d 0 = 2 mm, d 0 = 4 mm, and d 0 = 8 mm, respectively. 88 Chapter 5 -0.2 0 0.2 0.4 0.6 -6 -4 -2 0 2 -1 -0.5 0 0.5 1 -2 -1 0 1 2 -6 -4 -2 0 2 + - Figure 5.28: Left: The space charge potential of the beams with separation x=d 0 = 0:4 and d 0 = 2 mm at time step 26000. Right: The average velocity of v x . -0.2 0 0.2 0.4 0.6 -4 -2 0 2 4 -1 -0.5 0 0.5 1 -2 -1 0 1 2 -4 -2 0 2 4 + - Figure 5.29: Left: The space charge potential of the beams with separation x=d 0 = 0:4 and d 0 = 4 mm at time step 26000. Right: The average velocity of v x . The potential and the average v x plots are shown in Fig. 5.28 for d 0 = 2 mm, Fig. 5.29 for d 0 = 4 mm, and Fig. 5.30 for d 0 = 8 mm. We can see that the space charge potentials are similar. The average v x for the case of d 0 = 8 mm is around zero for both positive and negative beams, because they are still within the accelerating region, where they do not feel each other. Chapter 5 89 -0.2 0 0.2 0.4 0.6 0 2 4 6 8 -1 -0.5 0 0.5 1 -2 -1 0 1 2 0 2 4 6 8 + - Figure 5.30: Left: The space charge potential of the beams with separation x=d 0 = 0:4 and d 0 = 8 mm at time step 26000. Right: The average velocity of v x . Therefore, none of our simulations with dierent parameters (varied thruster separations and varied space charge eects) leads to a mixing of the two beams. And usually the separation between the two bipolar thrusters are much more than the separation studied in our simulations. So, the space charge of the two beams has negligible eects to cause a beam bending or mixing. 90 Chapter 6 Chapter 6 PP-PM Hybrid Model of Ion Beam As mentioned in Sec. 4.3, a PM model can be used to speed up the simulation in the downstream region. In this chapter, a PP-PM hybrid model will be developed and tested. 6.1 2-D Cylindrical Poisson Solver with Nonuniform Mesh We start at the 3-D cylindrical Poisson's equation (6.1) in SI units, where r, , and z denote the radial coordinate, the azimuth coordinate, and the axial coordinate, respectively, denotes the electric potential,% denotes the charge density," 0 is the vacuum permittivity. If we have axial symmetry, @=@ =@ 2 =@ 2 = 0, Eq. (6.1) reduces to Eq. (6.2), which is the 2-D cylindrical Poisson's equation in (r;z) coordinates. r 2 = @ 2 @r 2 + 1 r @ @r + 1 r 2 @ 2 @ 2 + @ 2 @z 2 = % " 0 (6.1) @ 2 @r 2 + 1 r @ @r + @ 2 @z 2 = % " 0 (6.2) Now, let's apply Taylor series to approximate the partial dierential terms in 1-D dis- crete direction. The Taylor series of at i + i and i i1 can be expressed as Eq. (6.3), ignoring higher-order terms, where the subscript i denotes the mesh grid index, i denotes the distance between i and i+1 , i1 denotes the distance between i1 and i . Chapter 6 91 Using subscripts to express the evaluating locations, we obtain Eq. (6.4). ( i + i ) =( i ) + 0 ( i ) i + 00 ( i ) 2 i =2 ( i i1 ) =( i ) 0 ( i ) i1 + 00 ( i ) 2 i1 =2 (6.3) i+1 = i + 0 i i + 00 i 2 i =2 i1 = i 0 i i1 + 00 i 2 i1 =2 (6.4) Solving 0 i and 00 i from Eq. (6.4), we obtain Eq. (6.5), where ! , shown in Eq. (6.6), is used to simplify the equations. Note that, if i = i1 = , Eq. (6.5) reduces to the nite dierence formulas for uniform mesh, i.e. 0 i = ( i+1 i1 )=(2) and 00 i = ( i+1 + i1 2 i )= 2 . 0 i = i+1 2 i1 i1 2 i i ( 2 i1 2 i ) i i1 ( i + i1 ) = i+1 ! 2 i1 i1 ! 2 i i ! [ 2 i1 2 i ] 00 i = i+1 i1 + i1 i i ( i + i1 ) i i1 ( i + i1 )=2 = 2 i+1 ! i1 + 2 i1 ! i 2 i ! ( i + i1 ) (6.5) ! = [ i i1 ( i + i1 )] 1 (6.6) Now, consider the 2-D (r;z) coordinates, we discretizer andz intor i andz j ,i andj are varied from 1 to N r and N z , respectively. The mesh sizes are given in Eq. (6.7). Therefore, we have N r N z mesh grids, (N r 1) (N z 1) meshes. Letting = r or = z, we substitute Eq. (6.5) into Eq. (6.2) to obtain the solution formula Eq. (6.8) with coecients listed in Eq. (6.9), where ! becomes ! r and ! z shown in Eq. (6.10). r i =jr i+1 r i j; i = 1; 2; 3;:::;N r 1 z j =jz j+1 z j j; j = 1; 2; 3;:::;N z 1 (6.7) i+1;j C i+1;j + i1;j C i1;j + i;j C i;j + i;j+1 C i;j+1 + i;j1 C i;j1 =% i;j =" 0 (6.8) 92 Chapter 6 C i+1;j =! r r i1 (2 + r i1 =r i ) C i1;j =! r r i (2 r i =r i ) C i;j =[2! r (r i + r i1 ) +! r (r 2 i1 r 2 i )=r i + 2! z (z j + z j1 )] C i;j+1 = 2! z z j1 C i;j1 = 2! z z j (6.9) ! r = [r i r i1 (r i + r i1 )] 1 ! z = [z j z j1 (z j + z j1 )] 1 (6.10) If mesh grids on the symmetry axis are also considered, i.e. r 1 = 0, a problem occurs when i = 1, because 1=r 1 !1. Using L'Hospital's rule, lim r!0 (@=@r)=r = @ 2 =@r 2 . Therefore, as r! 0, Eq. (6.2) becomes Eq. (6.11). Since @=@r = 0 as r! 0, set 0 1 = 0 in Eq. (6.5), and note that r 0 = r 1 due to symmetry, we obtain the condition 0 = 2 . Substituting Eq. (6.5) into Eq. (6.11) and using 0 = 2 , we can obtain the solution formula on the axis Eq. (6.12) with extra coecients listed in Eq. (6.13). 2 @ 2 @r 2 + @ 2 @z 2 = % " 0 (6.11) 2;j D 2;j + 1;j D 1;j + 1;j+1 C 1;j+1 + 1;j1 C 1;j1 =% 1;j =" 0 (6.12) D 2;j = 8! r r 1 D 1;j =[8! r r 1 + 2! z (z j + z j1 )] (6.13) To calculate the electric eld from the potential, we start at the relation Eq. (6.14). Substitute Eq. (6.5) into Eq. (6.14), we obtain Eq. (6.15). For the on-axis grids, E r = 0 due to symmetry. E r = @ @r E z = @ @z (6.14) Chapter 6 93 E r (i;j) =! r [ i+1;j r 2 i1 i1;j r 2 i i;j (r 2 i1 r 2 i )] E z (i;j) =! z [ i;j+1 z 2 j1 i;j1 z 2 j i;j (z 2 j1 z 2 j )] (6.15) 6.1.1 Test One test case is that of a vacuum potential, i.e. % = 0. The analytical solution is given in Eq. (6.16), where I 0 (x) is the modied Bessel function of order zero. The simulation domain is that r2 [0;L r ] andz2 [0;L z ], whereL r = 0:02m, L z = 0:08m, and 0 is set to be 1V. This test case was also used by Chao et al. [16]. Here, the three boundary conditions applied are all Dirichlet, namely (L r ;z), (r; 0), and (r;L z ) are given and xed. For simplicity, the Jacobi method is used to solve the Poisson's equation. Namely, the solution formulas Eq. (6.8) and Eq. (6.12) are directly used to solve for i;j in each iteration. (r;z) = 0 cos (2z=L z )I 0 (2r=L z ) I 0 (2L r =L z ) (6.16) 4 8 12 16 20 24 28 32 z 4 8 12 16 20 24 28 32 r 1 0:8 0:6 0:4 0:2 0 0:2 0:4 0:6 0:8 1 4 8 12 16 20 24 28 32 z 4 8 12 16 20 24 28 32 r 0 0:0002 0:0004 0:0006 0:0008 0:001 0:0012 Figure 6.1: Left: The analytical solution of the potential (in volts). Right: The absolute dierence between the analytical potential and the simulated potential. First, constant r and z are used, r = Lr=(N r 1) and z = Lz=(N z 1). N r and N z are both chosen to be 32. The analytical solution of the potential, Eq. (6.16), is plotted in Fig. 6.1 (left), in which the values of the colorbar are in the unit of volts. The 94 Chapter 6 absolute dierence between the analytical potential ( a ) and the simulated potential ( s ), i.e. j s a j, is plotted in Fig. 6.1 (right). The values of z and r in Fig. 6.1 indicate the indices of the mesh grid. The error evaluated using Eq. (6.17) is plotted over iteration in Fig. 6.2 (the solid line). We can see that the computation converges after about 3000 iterations, and Error3:5 eventually. In addition, Fig. 6.2 also shows the error of a case with N r N z = 100 100. We can see the computation converges after about 30000 iterations, and Error4:5 eventually. Error = log 10 [ Nr X i=1 Nz X j=1 j s (i;j) a (i;j)j=(N r N z )] (6.17) 5 4 3 2 1 0 0 10000 20000 30000 40000 Error Iteration 3232 100100 Figure 6.2: Error over iteration. Next, we test dierent nonuniform meshes. Linear mesh is shown in Eq. (6.18). Square- root mesh is shown in Eq. (6.19), where r 1 = 0:0001699948693166944 forN r N z = 3232 and r 1 = 0:00003023600957467768 for N r N z = 100 100 are used, and z 1 = 4r 1 . N 1 X i=1 i = N 1 X i=1 i 1 =L ) 1 = 2L N (N 1) (6.18) Chapter 6 95 N 1 X i=1 i = N 1 X i=1 p i 1 =L (6.19) 5 4 3 2 1 0 0 10000 20000 30000 40000 Error Iteration L: 32 L: 100 S: 32 S: 100 Figure 6.3: Error over iteration for nonuniform meshes. The Error plots of both linear and square-root tests are shown in Fig. 6.3, where L denotes linear, S denotes the square-root, and the numbers in the labels indicate the mesh sizes, i.e. N r , N z = 32 or 100. As we can see, the Poisson solver produces good accuracy when using nonuniform meshes. 6.2 Charge Deposition Suppose we have a particle at r p andz p in the 2-D (r;z) cylindrical domain. We rst locate the mesh that the particle is in. For example, if r p r ip and r p <r ip+1 , and if z p z jp and z p < z jp+1 , the particle is located in the mesh surrounded by the four mesh grids: (ip;jp), (ip + 1;jp), (ip;jp + 1), and (ip + 1;jp + 1). Now, we follow Verboncoeur's charge deposition method [76], which is capable to accurately compute the charge density in nonuniform meshes using arbitrary particle-mesh interpolation. The formulas with linear interpolation are given in Eq. (6.20), where Q = 2r p q , q is the linear charge density in the unit of C/m, the 96 Chapter 6 corrected volumes V i are given in Eq. (6.21). % ip;jp = Q V ip jr ip+1 r p j r ip jz jp+1 z p j z jp % ip+1;jp = Q V ip+1 jr p r ip j r ip jz jp+1 z p j z jp % ip;jp+1 = Q V ip jr ip+1 r p j r ip jz p z jp j z jp % ip+1;jp+1 = Q V ip+1 jr p r ip j r ip jz p z jp j z jp (6.20) V i =z jp [r i+1 (r i +r i+1 )r i1 (r i1 +r i )]=3; 1<i<N r V i =z jp (r i+1 r i )(2r i +r i+1 )=3; i = 1 V i =z jp (r i r i1 )(r i1 + 2r i )=3; i =N r (6.21) 6.2.1 Test The test case is the axial electric eld due to a ring of charge, which is equivalent to placing a single particle with linear charge density q in the r-z simulation domain, r2 [0;L r ] and z2 [0;L z ]. If the particle is placed at (r p ;z p ), the analytical axial electric eld on thez axis is Eq. (6.22). E z = r p q (zz p ) 2" 0 [(zz p ) 2 +r 2 p ] 3=2 (6.22) Here we choose L r = 1 m, L z = 5 m, N r = 201, N z = 301, q = 8:85 10 12 C/m, the iteration times to be 20000, and use the Poisson solver to obtain the simulated axial electric eld. For simplicity, the three boundary conditions applied are all Dirichlet, namely (L r ;z), (r; 0), and (r;L z ) are xed to be 0 V. First, we test the uniform mesh. When r p = 0:1 m and z p = 2:5 m, the particle is located exactly on the mesh grid (21; 105). The potential is shown in Fig. 6.4 (left), the colorbar is in the unit of log 10 Volts. The Error, Eq. (6.17), over iteration is plotted in Fig. 6.4 (right). The comparison plot of E z is shown in Fig. 6.5 (left), where the label S represents the simulated results, the label A represents the analytical solution. The absolute dierence between the simulated results and the analytical solution is shown in Fig. 6.5 (right). Chapter 6 97 0 50 100 150 200 250 300 z 0 25 50 75 100 125 150 175 200 r 8 7 6 5 4 3 2 1 0 7 6:5 6 5:5 5 0 5000 10000 15000 20000 Error Iteration Figure 6.4: Left: Simulated potential. Right: Error over iteration. 2 1:5 1 0:5 0 0:5 1 1:5 2 0 50 100 150 200 250 300 E z z S A 0 0:005 0:01 0:015 0:02 0 50 100 150 200 250 300 Dierence z Figure 6.5: Left: Comparison between simulated E z (labeled as S) and analytical E z (labeled as A) in the unit of V/m. Right: The absolute dierence. Next, we change the particle position to r p = 0:1 + r=2 m, z p = 2:5 + z=2 m, and r p = 0:1 + r=4 m, z p = 2:5 + z=4 m, and label these two cases as 2 and 4, respectively. Also, we test the linear mesh using r p = 0:1 m, z p = 2:5 m, and label this case as L. The absolute dierence between the simulated results and the analytical solution of these three cases are all plotted in Fig. 6.6. As we can see, all cases produce correct results with an 98 Chapter 6 acceptable accuracy. 0 0:005 0:01 0:015 0:02 0:025 0:03 0:035 0 50 100 150 200 250 300 Dierence z 2 4 L Figure 6.6: The absolute dierence of dierent cases. 6.3 Particle Mapping and Push As the pure PP model, the hybrid model still keeps the 3-D particle positions and velocities. Given the position of one particle, x, y, and z, we map the particle into the 2-D (r;z) coordinates as shown in Eq. (6.23), and map the solved radial eld E r back into the 3-D (x;y;z) coordinates as shown in Eq. (6.24), where ^ r xy = r xy =jr xy j and r xy is the radial position vector in thex-y plane as used in Eq. (2.28). KnowingE x ,E y , andE z , the particle push algorithm explained in Sec. 2.4.3 is still used. Therefore, the PP-PM hybrid model is really a 3-D PP and 2-D PM hybrid model. r p = p x 2 +y 2 z p =z (6.23) E x =E r ^ r xy ^ x E y =E r ^ r xy ^ y E z =E z (6.24) Chapter 6 99 6.4 PP-PM Hybrid Model Test Figure 6.7: PP-PM hybrid test setup. 0 0:025 0:05 0:075 0:1 r=d 0 -2 -1 0 1 2 z=d 0 7 6 5 4 3 2 1 0 1 0 0.006 0.025 0.055 0.1 r=d 0 -2 -1 0 1 2 z=d 0 7 6 5 4 3 2 1 0 1 Figure 6.8: Charge density distribution. The colorbar shows log 10 of volts. Electric potential distribution in volts. The case with uniform mesh is on the left, while the linear one is on the right. 100 Chapter 6 0 0:025 0:05 0:075 0:1 r=d 0 -2 -1 0 1 2 z=d 0 0 0:1 0:2 0:3 0:4 0:5 0:6 0:7 0:8 0 0.006 0.025 0.055 0.1 r=d 0 -2 -1 0 1 2 z=d 0 0 0:1 0:2 0:3 0:4 0:5 0:6 0:7 0:8 Figure 6.9: Electric potential distribution in volts. The case with uniform mesh is on the left, while the linear one is on the right. 1 2 3 4 5 6 7 1 0:750:50:25 0 0:25 0:5 0:75 1 log 10 E r z=d 0 Figure 6.10: log 10 E r over z of the pure PP model. In this section, a PP-PM hybrid model is presented. As shown in Fig. 6.7, the PP model is used in the upstream region z d 0 =2, the PM model is used in the downstream region z <d 0 =2. In more detail, for those droplets in the upstream region, their forces are calculated using the PP algorithm considering the PP forces due to droplets in both regions, thus the computation operation scales as N up N, where N up denotes the number of droplets in the Chapter 6 101 upstream region, N denotes the total number of droplets. When the rst droplet passes z =d 0 =2, the PM model starts to work. The simulation domain of the PM model is shown using the red lines in Fig. 6.7, which is a rectangular region in the 2-D (r;z) coordinates with z2 [2d 0 ; 2d 0 ] andr2 [0; 0:1d 0 ]. The electric potential in the whole simulation domain due to all droplets is solved using the PM algorithm, but only those droplets in the downstream region are using the electric eld information computed from the PM model. For simplicity, the top, bottom, and right boundary conditions of the PM model is set to be Dirichlet with = 0 V. Note that the analytical solution of the accelerating eld (Sec. 2.4.2) is still used for all droplets. The chosen material, practical, and simulation parameters are as the same as those used in the case A 0 with t _ N d = 1=2 in Sec. 4.1. 0 0:05 0:1 0:15 0:2 0:25 0:3 1 0:750:50:25 0 0:25 0:5 0:75 1 Relative Error of E r z=d 0 6 11 21 67 101 0 0:05 0:1 0:15 0:2 0:25 0:3 1 0:750:50:25 0 0:25 0:5 0:75 1 Relative Error of E r z=d 0 L6 L11 L21 L67 L101 Figure 6.11: Relative errors of E r between the pure PP case and hybrid cases. 102 Chapter 6 First, consider the numbers of grid points N r = 101 and N z = 101. One case applies both uniform meshes inr andz directions, which results in r=d 0 = 0:001 and z=d 0 = 0:04. Another case applies uniform mesh inz direction and linear mesh inr direction, which results in r 1 =d 0 = 2 10 5 . When applying the Poisson solver, the Successive Overrelaxation (SOR) method is used with an overrelaxation factor 1.95. The iteration ceases when the maximum value ofj 0 i;j i;j j=j 0 i;j + i;j j is less than 1 10 5 , where the prime denotes the new updated value. The simulations run for 6500 time steps, at the end of which the beam front reaches z=d 0 =1. The charge density on the 2-D mesh grids is shown in Fig. 6.8. The corresponding electric potential is shown in Fig. 6.9. Both are at the simulation time step 6500. To compare the hybrid model with the pure PP model, we plotted log 10 E r over z for the pure PP model, shown in Fig. 6.10, because E r is sensitive to the eld solver. In Fig. 6.10, each dot denotes a droplet, and the black line is an averaged result. This averaged result of the hybrid model is also diagnosed and the relative error between it with that of the pure PP model is shown in Fig. 6.11 (the uniform case on the top, the line labeled 101; the linear case on the bottom, the line labeled L101), namelyj(E p r E h r )=E h r j, where E p r denotes the PP result, E h r denotes the hybrid result. As we can see, the hybrid results match very well with those of the pure PP model for both uniform mesh and linear mesh. Then, dierent numbers of mesh grid inr direction are tested, and the results are shown in Fig. 6.11 as well, labeled using the corresponding numbers of the grid. The direct result is that smaller numbers of mesh grid have larger errors. One signicant nding is that, for the uniform case, when the mesh size is smaller than the case with 67, the errors reach a smallest value. The mesh size of the case 67 is about 1.5 m, which is the sucient minimum mesh size we discussed in Sec. 4.3. Case PP 6 11 22 67 101 L6 L11 L22 L67 L101 Time (s) 755 447 803 1554 4673 7605 450 779 1434 3987 7521 Table 6.1: Comparison of the elapsed simulation CPU time. At last, the elapsed simulation CPU time is shown in Tab. 6.1. We can see that for this test, most cases of the hybrid model are slower than the PP model, due to the fact that the parameters we chose result in a small number of droplets, as well as the early time at which we stop the simulation. Chapter 7 103 Chapter 7 Hybrid Model Application of Ion Beam In this chapter, the developed PP-PM hybrid model is applied to simulate a more realistic case. We refer the PIC simulation of a colloid thruster done by the group of Prof. Deborah Levin [80]. The same simulation setup and parameters as those studied in [80] will be applied here using the PP-PM hybrid model, thus our simulation results can be compared with theirs. 7.1 Simulation Setup and Parameters Figure 7.1: Simulation domain and setup. 104 Chapter 7 Following [80], the ionic liquid 1-ethyl-3-methylimidazolium/tetra ouroborate (EMIM-BF 4 ) is used. EMIM-BF 4 has a chemical formula C 6 H 11 BF 4 N 2 . A positive potential will be set on the needle (or capillary), and the accelerating plate (or grid) will be grounded to have zero volts. We thus consider the injected particles are all positive ions, EMIM + (C 6 H 11 N 2 ), which has an ion mass about 1:846 10 25 kg and charge 1:602 10 19 C. If EMIM-BF 4 has a mass density 1240 kg/m 3 , a relative permittivity 14, a conductivity 1.31 S/m, and we still use f(") = 18, when the mass ow rate is 2:47 10 12 kg/s (the same as one case in [80]), we can calculate the radius of the jet R j 0:919 nm using the formulas given in Sec. 2.2. The initial axial velocity of each particle can thus be determined, v d 751:16 m/s, which equals the volume ow rate over the cross section area of the jet. The simulation domain and setup are shown in Fig. 7.1. The length of OD is set to be 0.3 mm, OF is 0.15 mm, the rectangular region surrounded by GOA simulates the needle with size OG= 2:5 m, OA= 5 m, and the potential on it is xed to be 1500 V. The boundary conditions of GF and FE are Neumann, while DE is xed at zero volts. Particles are injected using the same algorithm described in Sec. 2.3. The injection volume starts at the surface corresponding to A. The radius of the volume is R j , the height of it is v d t. Particle positions are randomly distributed in the injection volume. The axial velocity is xed to be v d , and the radial velocity is ignored. Note that the injection applied here is dierent from the one used in [80], since they used the data from Molecular Dynamics simulations. 0 100 200 300 400 500 0 100 200 300 400 500 r z 0 300 600 900 1200 1500 0 100 200 300 400 500 0 100 200 300 400 500 r z 0 300 600 900 1200 1500 Figure 7.2: The accelerating potential distribution in volts with 500 500 mesh grids, and uniform mesh on the left, linear mesh on the right. Before the start of the simulation, the background accelerating eld labeled as 0 , E r0 , and E z0 are prepared in advance, which may be obtained using a ner mesh to provide a more accurate acceleration for particles. Now, let's denote the elds that are only due to those particles in the upstream region (between A and B in Fig. 7.1) E r1 (r;z) andE z1 (r;z) Chapter 7 105 with potential 1 , and the elds that are only due to those particles in the downstream region (between B and C in Fig. 7.1) E r2 (r;z) and E z2 (r;z) with potential 2 . Therefore, the elds at the positions of those particles in the upstream region are the superpositions E r0 +E r2 +E rpp andE z0 +E z2 +E zpp , where the subscriptpp denotes the elds solved using the PP algorithm in the upstream region. The elds at the positions of those particles in the downstream region are the superpositions E r0 +E r1 +E r2 and E z0 +E z1 +E z2 . 1:2 1:4 1:6 1:8 2 2:2 2:4 2:6 2:8 3 5:5 6 6:5 7 7:5 8 E z =10 8 (V/m) z (nm) Figure 7.3: E z of particles over z for dierent cases. 0 5 10 15 20 25 30 35 5:5 6 6:5 7 7:5 8 v z (km/s) z (nm) Figure 7.4: v z of particles over z for dierent cases. 106 Chapter 7 When solving 0 , E r0 , and E z0 , the boundary conditions shown in Fig. 7.1 are used. When solving 1 , 2 ,E r1 ,E z1 ,E r2 , andE z2 , the rectangular needle region GOA is removed, and the boundary condition on OF is all @=@z = 0. Other conditions are still the same as those shown in Fig. 7.1. Under this treatment, the computation of the PP part of the hybrid model scales asN up N up , whereN up denotes the number of particles in the upstream region. The Poisson's equation needs to be solved twice in each simulation loop for the upstream region and the downstream region, respectively. 0 0.1 0.2 0.3 z (mm) 0 0.05 0.1 0.15 r (mm) 14 15 16 17 18 19 0 0.1 0.2 0.3 z (mm) 0 0.05 0.1 0.15 r (mm) 14 15 16 17 18 19 0 0.1 0.2 0.3 z (mm) 0 0.05 0.1 0.15 r (mm) 14 15 16 17 18 19 0 0.1 0.2 0.3 z (mm) 0 0.05 0.1 0.15 r (mm) 14 15 16 17 18 19 Figure 7.5: Number density plots of case S1 (top left panel), S2 (top right panel), S3 (bottom left panel), and S4 (bottom right panel), in the unit of log 10 m 3 . Chapter 7 107 7.2 Accelerating Field Here, the accelerating elds 0 obtained using both uniform and nonuniform meshes with dierent number of mesh grids are tested and compared. Fig. 7.2 shows 0 obtained using 500 500 mesh grids. The left panel used a uniform mesh, while the right panel used a linear mesh in bothr andz directions. A pure PP model is applied using these accelerating elds with those parameters described in Sec. 7.1. The obtained plots of E z over z are compared. As shown in Fig. 7.3, the lines from the top to the bottom correspond to the cases with 1000 1000 linear mesh grids, 500 500 linear mesh grids, 1000 1000 uniform mesh grids, 500 500 uniform mesh grids, 200 200 linear mesh grids, and 200 200 uniform mesh grids, respectively. The lines do not match with each other very well, because we inject particles on the locations of specic grids that are close to z = 5:4 m. Thus the actual injection locations are dierent from case to case. Here, we do not inject particles at the exact location z = 5 = OA m to avoid letting those injected particles be trapped in the potential well (the GOA region shown in Fig. 7.1). Therefore, the cases with 200 200 linear mesh grids and 200 200 uniform mesh grids have more obviously smaller E z than other cases, because particles are injected at locations of larger z, which have smaller E z . Other than the dierence caused by the injection location, the line shapes of all cases are similar. 0 0.1 0.2 0.3 z (mm) 0 0.05 0.1 0.15 r (mm) 14 15 16 17 18 19 0 0.1 0.2 0.3 z (mm) 0 0.05 0.1 0.15 r (mm) 14 15 16 17 18 19 Figure 7.6: Number density plots of case S3 with accelerating potential 1250 V (left panel) and 1750 V (right panel), in the unit of log 10 m 3 . The resultant axial velocity distribution over z is plotted in Fig. 7.4 by shifting all the injection locations to 5.4 m for the cases with 1000 1000 linear mesh grids, 500 500 linear mesh grids, 1000 1000 uniform mesh grids, and 500 500 uniform mesh grids. As 108 Chapter 7 we can see that they match with each other good enough. Since the linear cases have ner mesh near the needle, where the major potential drop appears, we will use the accelerating eld produced by the 1000 1000 linear mesh grids in all later simulations. The injection location is at the grid that corresponds to z = 5:42158 m. 0 10000 20000 30000 40000 50000 60000 0 0:05 0:1 0:15 0:2 0:25 0:3 v z (m/s) z (mm) 0 10000 20000 30000 40000 50000 60000 0 0:05 0:1 0:15 0:2 0:25 0:3 v z (m/s) z (mm) 0 10000 20000 30000 40000 50000 60000 0 0:05 0:1 0:15 0:2 0:25 0:3 v z (m/s) z (mm) 0 10000 20000 30000 40000 50000 60000 0 0:05 0:1 0:15 0:2 0:25 0:3 v z (m/s) z (mm) Figure 7.7: Phase space z-v z plots of case S1 (top left panel), S2 (top right panel), S3 (bottom left panel), and S4 (bottom right panel). Chapter 7 109 7.3 Simulation Results 0 10000 20000 30000 40000 50000 60000 0 0:05 0:1 0:15 0:2 0:25 0:3 v z (m/s) z (mm) Figure 7.8: Phase space z-v z plots of case S3 with accelerating potentials 1750 V, 1500 V, and 1250 V from the top curve to the bottom curve, respectively. The mass ow rate is varied to be 0.72 10 12 kg/s, 1.38 10 12 kg/s, 2.47 10 12 kg/s, and 3.51 10 12 kg/s, following the cases labeled as S1, S2, S3, and S4 in [80]. And for the case S3, the accelerating potential is varied to be 1250 V and 1750 V, besides 1500 V. For all cases, the simulation time step length is chosen to be t = 1=2 _ N d . t = 1=5 _ N d is also tested and the results (not shown here) indicate that t = 1=2 _ N d is small enough. The separation between the PP model and the PM model (the location of B shown in Fig. 7.1) is atz = 7m, which is large enough to cover all the short-range interactions that can only be captured by PP. The number density plots of the cases S1 to S4 are shown in Fig. 7.5. We can see that a larger mass ow rate leads to a larger beam expansion in the radial direction, which agrees with the results found in [80]. The beam radii at z = 0:3 mm are measured from Fig. 7.5. They are about 0.06 mm, 0.075 mm, 0.09 mm, and 0.105 mm, respectively, for cases from S1 to S4. However, those shown in the reference [80] are about 0.014 mm, 0.018 mm, 0.024 mm, and 0.03 mm, respectively, all of which are smaller than our results. It is believed that the multi-scale PIC model utilized in [80] still does not accurately resolve the short-range particle interactions in the beam upstream region. Because the interactions are underestimated, the beam expands less in the radial direction. The number density plots of the case S3 with dierent accelerating potentials are shown in Fig. 7.6. We can see that a larger potential leads to a smaller beam expansion in the radial direction, because the beam has a smaller axial velocity when the potential is smaller, 110 Chapter 7 which provides more time for the beam to expand radially. Again, this trend agrees with that found in [80], but the beam radii shown in [80] are smaller than ours. The phase-space plotsz-v z are shown in Fig. 7.7 for the cases S1 to S4, and Fig. 7.8 for the case S3 with dierent accelerating potentials. Overall these results are similar to those obtained in [80]. For all the cases, the CPU time is about 18.98 hours with about 85400 ions at the end of the simulation. If OpenMP is used with 16 threads, the elapsed simulating time would be about just 1.19 hours. However, the elapsed simulation time is not mentioned in the reference [80], so we cannot compare our results with theirs. 7.4 Reproduction of Larger Expansion Angle 0 r c Force Interparticle Distance Figure 7.9: Cuto applied in the PP part of the hybrid model. To further verify the accuracy of the PP-PM results, we apply a cuto (weakened force) on the PP part of the hybrid model. The force function of the PP model is modied as shown in Fig. 7.9, where at the cuto interparticle distance r c , the Coulomb force starts to be weakened linearly till zero. This articial force is analog to the weakened force of a PM model within a mesh using linear interpolation, and was applied previously [89]. Fig. 7.10 shows the number density plots using r c = 30R j 28 nm (left panel) and r c = 42R j 39 nm (right panel) for the case S3. We can see that a larger r c leads to a smaller radial expansion angle of the beam. Fig. 7.10 (right panel) indicates that we can reproduce the beam size 0.03 mm shown in the paper [80] using r c = 42R j . Chapter 7 111 0 0.1 0.2 0.3 z (mm) 0 0.05 0.1 0.15 r (mm) 14 15 16 17 18 19 0 0.1 0.2 0.3 z (mm) 0 0.05 0.1 0.15 r (mm) 14 15 16 17 18 19 Figure 7.10: Number density plots of case S3 with cutos in the unit of log 10 m 3 . 5 10 15 20 25 30 35 1 10 100 1000 Expansion Cone Angle (Degree) r c =R j Figure 7.11: Relation between r c and the expansion angle. In addition, more values of r c are tested, the relation between r c and the expansion angle is obtained in Fig. 7.11. We can see that the angle deviates from the correct value (the black dashed line) signicantly when r c > 10R j . At last, a local view of the beam near the injection area is shown in Fig. 7.12, which is the case S3 without a cuto. A quadratic function can t well the shape of the beam surface 112 Chapter 7 near the injection area. When z > 7:5 m, the shape of the beam surface becomes linear, indicating that the radial particle interactions become less and less signicant, thus most particles eventually move with almost constant velocities in the radial direction. Therefore, we can say that the maximum necessary PP region is about z < 7:5 m for this simulation case. For the minimum PP region, however, it depends on the applied PM model and the mesh size, so it is hard to give a specic value. 4 3 2 1 0 1 2 3 4 6 7 8 9 10 11 12 13 14 15 x (m) z (m) Figure 7.12: Local view of the beam near the injection area. Chapter 8 113 Chapter 8 Summary and Future Work At the end, we give a summary of all the works that have been done in this dissertation and some possible future works. The main contributions of this dissertation are listed as follows. A new particle simulation model, the Particle-Particle (PP) model is developed for simulating the beam of colloid thrusters with high accuracy. The physics of the charge neutralization of colloid thrusters is investigated using the PP model, including an alternating current (AC) mode and a bipolar mode. A new PP-PM hybrid model is developed to speed up the simulations for cases with a larger number of particles, while still maintaining the high accuracy. The detailed summary is given as follows. In chapter 1, a background of electrospray thrusters and colloid thrusters is given. The motivations of this work include the facts that many unknowns remain with respect to the physics of colloid thrusters, simulations on the colloid thruster beam is rare, the commonly used particle-in-cell (PIC) method has its diculties for simulating the colloid thruster beam, and a PP model which does not apply a mesh system can resolve short-range interactions accurately, thus it is satiable for simulating the colloid thruster beam. The objectives of this work are therefore: (a) Develop a PP model that can simulate the colloid thruster beam with high accuracy. (b) Apply the PP model to obtain a complete and precise description of the droplet dynamics and reveal more unknown physics. (c) Study in detail the Coulomb collisions between droplets and its variations in dierent circumstances as well as its dependence on other parameters. (d) Analyze results of the PP model to give criteria and guidance for developing other models. Then, a literature review is given, including the most relevant simulation works on electrospray thruster beam, among which the paper [80] about PIC simulations on a colloid thruster beam is of considerable referential importance, since in chapter 7, the simulation results using a proposed PP-PM hybrid model are compared with those in [80]. 114 Chapter 8 In chapter 2, the development of the PP simulation model is described in detail, in- cluding: the simulation domain and setup; the input and derived parameters based on the available theories of colloid thrusters; the reasonable simulation initialization and assump- tions on droplet breakup; the particle-particle eld, the analytical accelerating electric eld, and the integration of the equations of motion; the simulation procedure and parallelization using OpenMP and CUDA; diagnoses such as system energy, droplet macro-temperature, thrust and specic impulse, and equivalent length. In chapter 3, the PP code is veried through a comparison with the binary collision theory and a space charge eect study. The binary collision comparison indicates that the PP model is cable to produce accurate particle-particle Coulomb collisions, and the time step length used in all simulation cases are small enough. The space charge eects of the colloid thrusters under the parameter range simulated are studied. The PP algorithm and code are veried through a comparison with the Child-Langmuir law and a Particle-in-Cell (PIC) method. It is found that the Child-Langmuir law is inappropriate to be applied for the scenario of colloid thrusters, and possible reasons are given and discussed: (1) In our case, the applied electric eld is not a constant as the one in the Child-Langmuir law produced from two parallel plates. Instead, it has a much larger magnitude near the Taylor cone tip, and drops signicantly toward the downstream region. (2) Our droplets have initial velocities that are determined by the jet speed, which is dierent from the zero initial assumption in the Child-Langmuir law. (3) We randomly inject droplets in an injection volume to simulate the real breakup of droplets from the jet front, which increases the distances between those particles that are just injected and those that are about to be injected. (4) The Child- Langmuir law assumes a 1D scenario, whose particles shall be innite large charged sheets. In our case, we simulate a real thin beam of real charged droplets along the beam axis. Thus, there is not an innite number of other charged droplets that are located o-axis and produce repulsive forces on those on-axis real droplets. (5) The emitted droplets under Coulomb's law expand radially signicantly. An expanded beam front will have a much less impact on the newly injected particles at the beam end. In chapter 4, some reasonable parameters are chosen and several characteristic simula- tions are executed. Simulation results of phase space plots, system energy, droplet number density, macro-temperature, velocity distribution function, and the probability density func- tion (PDF) of the equivalent length are obtained, which indicate that the PP model can be successfully applied to simulate the colloid thruster beam and all information of droplet dynamics can be obtained from the simulation results. Also, from the PDF of the equivalent length, a suciently minimum mesh size of a corresponding PM model, which can capture all the necessary short-range collisions, can be obtained. In addition, the initial droplet thermal velocity is studied at the end of chapter 4, and found that it aects the PDF of the equivalent length, thus one may need to consider a reasonable initial droplet thermal velocity for more practical applications. The space charge eects are also discussed in more details based on the simulation parameters. In chapter 5, the PP model is applied to investigate the alternating current (AC) mode of colloid thrusters, which is used to avoid charge accumulation. Since, both positively and Chapter 8 115 negatively charged droplets are involved, the recombination of droplets are considered. Two types of alternating functions are studied, one is a square wave, which shows the droplet dynamics under the alternating potential, and the fact that recombination only occurs near the injection volume under the assumptions of the simulation; the other is a step function, meaning that the potential only alternates once, the transition time (t x ), during which the thrust changes back to its stable value after one alternation, is studied under dierent volume ow rates ( _ V ), accelerating potentials ( 0 ), and the accelerating plate spacing (d 0 ). The simulation results indicate that t x / d 0 1=2 0 _ V 1=4 , and it is found that this relation agrees with the reasoning that t x is proportional to the travel time of one droplet from the injection volume to the accelerating plate under the assumption that it gains the nal exit speed very soon after being injected. In addition, a bipolar operation mode is studied, where two thrusters are placed side by side with dual polarity to avoid the charge accumulation on the spacecraft. Space charge eects are studied by varying the separation between the two thrusters, the accelerating potential, the tip-plane distance, and the volume ow rate. The simulation results indicate that the space charge eects are negligible for the bipolar system, so the two owing beams are not aecting each other. In chapter 6, a PP-PM hybrid model is proposed to speed up the simulation for cases with larger number of droplets. Since the colloid thruster beam has a symmetric nature, the PM part of the model is chosen to be 2-D cylindrical. A 2-D cylindrical Poisson solver with nonuniform mesh is introduced. The nonuniform mesh is meaningful and important, because the number density of the colloid thruster beam drops signicantly from the upstream region to the downstream region. The charge deposition for the 2-D cylindrical nonuniform mesh system is described as well referring Verboncoeur's method [76]. A test case that is the same as a previous pure PP case is simulated using the PP-PM hybrid model. The results show the validity of the hybrid model, and show that the previously found suciently minimum mesh size (1.5m) is really the suciently minimum mesh size, since further decreasing the mesh size does not further contribute to obtain a smaller relative error, when comparing the results from the hybrid model with those from the pure PP model. In chapter 7, the developed PP-PM hybrid model is applied to simulate a more realistic case, the one that is studied by the group of Prof. Deborah Levin [80]. The simulation setup and the material used (EMIM-BF 4 ) all follow those in the reference [80]. Instead of large droplets, small EMIM + ions are used, which results in a large number of particles, so the pure PP model is not capable any more. The trend of the simulation results agree with those in [80], but it is found that the beam expands more radially than that in [80], because their multi-level PIC still use a mesh size that is not small enough near the injection volume, so the short-range interactions there can not be fully resolved, while the PP-PM hybrid model applies the accurate PP model near the injection volume, which captures all short-range interactions. The limitations of the dissertation are listed as follows. No droplet breakup and deformation are included during the simulation. Only the Coulomb force is included in the droplet interaction. 116 Chapter 8 Mixtures of droplets of dierent size, dierent q=m, or mixtures of droplet and ions are not considered. The actual accelerator grid is not considered to simulate the droplet extraction from the thruster. For the PP model, boundary conditions are not easy to be implemented as those in the PM model. Possible future works are listed as follows. More realistic droplet breakup can be considered based on the calculation of surface tension, electrostatic forces, and evaporation. And droplet breakup can be considered during the beam ow. Mixtures of droplets of dierent size, dierent q=m, or mixtures of droplet and ions can be considered. The actual accelerator grid can be considered to simulate the droplet extraction from the thruster. More convenient approaches to implement boundary conditions for the PP model can be investigated. Further optimization and parallelization for both the PP model and the PP-PM hybrid model can be worked on. More realistic cases can be simulated using the developed models, and compare with experimental results. Bibliography [1] URL https://developer.nvidia.com/cuda-zone/. [2] URL http://www.fortran.com/. [3] URL http://hpcc.usc.edu/. [4] URL https://www.openmp.org/. [5] B. P. Abbott and et al. Observation of gravitational waves from a binary black hole merger. Phys. Rev. Lett., 116:061102, Feb 2016. doi:10.1103/PhysRevLett.116.061102. [6] P. K. Bhattacharjee, T. M. Schneider, M. P. Brenner, G. H. McKinley, and G. C. Rutledge. On the measured current in electrospinning. Journal of Applied Physics, 107 (4):044306, 2010. doi:10.1063/1.3277018. [7] C. Birdsall and A. Langdon. Plasma physics via computer simulation. McGraw-Hill, New York, 1985. [8] J. A. Bittencourt. Fundamentals of Plasma Physics. Springer, 2004. [9] Arnaud Borner and Deborah A. Levin. Coupled molecular dynamics|3-d poisson sim- ulations of ionic liquid electrospray thrusters. IEEE Transactions on Plasma Science, 43:295, 2015. doi:10.1109/TPS.2014.2327913. [10] Arnaud Borner, Zheng Li, and Deborah A. Levin. Modeling of an ionic liquid electro- spray using molecular dynamics with constraints. The Journal of Chemical Physics, 136:124507, 2012. doi:10.1063/1.3696006. [11] Arnaud Borner, Zheng Li, and Deborah A. Levin. Modeling of an ionic liquid electro- spray using molecular dynamics with constraints. The Journal of Chemical Physics, 123:124507, 2013. doi:10.1063/1.3696006. [12] Arnaud Borner, Zheng Li, and Deborah A. Levin. Prediction of fundamental properties of ionic liquid electrospray thrusters using molecular dynamics. The Journal of Chemical Physics, 117:6768, 2013. doi:10.1021/jp402092e. [13] Arnaud Borner, Pengxiang Wang, and Deborah A. Levin. In uence of electrical bound- ary conditions on molecular dynamics simulations of ionic liquid electrosprays. The Journal of Chemical Physics, 90:063303, 2015. doi:10.1103/PhysRevE.90.063303. 117 118 BIBLIOGRAPHY [14] Jorge Carretero and Manuel Mart nez-S anchez. Numerical simulation of a colloidal thruster in the droplet regime. Computer Physics Communications, 164(1):202 { 208, 2004. ISSN 0010-4655. doi:10.1016/j.cpc.2004.06.074. Proceedings of the 18th Interna- tional Conferene on the Numerical Simulation of Plasmas. [15] Jorge Alejandro Carretero-Benignos. Numerical simulation of a single emitter colloid thruster in pure droplet cone-jet mode. Ph.D. dissertation, Massachusetts Institute of Technology, Department of Aeronautics and Astronautics, 2005. [16] Edward H. Chao, Stephen F. Paul, Ronald C. Davidson, and Kevin S. Fine. Di- rect numerical solution of poisson's equation in cylindrical (r,z) coordinates. Prince- ton Univ., Princeton Plasma Physics Lab., NJ (United States), PPPL-3251, 1997. doi:10.2172/304205. [17] Paul R. Chiarot, Sergey I. Gubarenko, Ridha Ben Mrad, and Pierre E. Sullivan. On the pulsed and transitional behavior of an electried uid interface. Journal of Fluids Engineering, 131:091202, 2009. doi:10.1115/1.3203203. [18] C. D. Child. Discharge from hot cao. Phys. Rev., 32:492{511, 1911. doi:10.1103/PhysRevSeriesI.32.492. [19] Roland Clift, John R. Grace, Martin E. Weber, and Martin F. Weber. Bubbles, drops, and particles. Academic Press, 1978. [20] Daniel G. Courtney, Nereo Alvarez, and Nathaniel R. Demmons. Electrospray Thrusters for Small Spacecraft Control: Pulsed and Steady State Operation. doi:10.2514/6.2018- 4654. [21] Daniel G. Courtney, Simon Dandavino, and Herbert Shea. Comparing direct and in- direct thrust measurements from passively fed ionic electrospray thrusters. Journal of Propulsion and Power, 32(2):392{407, 2016. doi:10.2514/1.B35836. [22] Daniel G. Courtney, Herbert Shea, Kathe Dannenmayer, and Alexandra Bulit. Charge neutralization and direct thrust measurements from bipolar pairs of ionic electrospray thrusters. Journal of Spacecraft and Rockets, 55(1):54{65, 2018. doi:10.2514/1.A33863. [23] John W. Daily and Michael M. Micci. Ionic velocities in an ionic liquid under high electric elds using all-atom and coarse-grained force eld molecular dynamics. The Journal of Chemical Physics, 131:094501, 2009. doi:10.1063/1.3197850. [24] L. de Juan and J.Fern andez de la Mora. Charge and size distributions of electrospray drops. Journal of Colloid and Interface Science, 186(2):280 { 293, 1997. ISSN 0021-9797. doi:10.1006/jcis.1996.4654. [25] J. Fern andez de la Mora and I. G. Loscertales. The current emitted by highly conducting taylor cones. Journal of Fluid Mechanics, 260:155{184, 1994. doi:10.1017/S0022112094003472. BIBLIOGRAPHY 119 [26] J. Fern andez De La Mora and I. G. Loscertales. The current emitted by highly conducting taylor cones. Journal of Fluid Mechanics, 260:155{184, 1994. doi:10.1017/S0022112094003472. [27] Nathaniel R. Demmons, Nick Lamarre, John K. Ziemer, Morgan Parker, and Douglas Spence. Electrospray Thruster Propellant Feedsystem for a Gravity Wave Observatory Mission. doi:10.2514/6.2016-4739. [28] Weiwei Deng and Alessandro Gomez. In uence of space charge on the scale-up of multiplexed electrosprays. Journal of Aerosol Science, 38(10):1062 { 1078, 2007. ISSN 0021-8502. doi:10.1016/j.jaerosci.2007.08.005. [29] M. Armano et al. Sub-femto-g free fall for space-based gravitational wave ob- servatories: Lisa pathnder results. Phys. Rev. Lett., 116:231101, Jun 2016. doi:10.1103/PhysRevLett.116.231101. [30] P. P. Ewald. The calculation of optical and electrostatic grid potential. ANNALEN DER PHYSIK, 64:253{287, 1921. [31] Alfonso M. Ga~ n an Calvo. Cone-jet analytical extension of taylor's electrostatic solution and the asymptotic universal scaling laws in electrospraying. Phys. Rev. Lett., 79:217{ 220, Jul 1997. doi:10.1103/PhysRevLett.79.217. [32] A M Ga~ n an-Calvo, N Rebollo-Mu~ noz, and J M Montanero. The minimum or nat- ural rate of ow and droplet size ejected by taylor cone{jets: physical symmetries and scaling laws. New Journal of Physics, 15(3):033035, mar 2013. doi:10.1088/1367- 2630/15/3/033035. [33] Alfonso M. Ga~ n an-Calvo, Jos e M. L opez-Herrera, Miguel A. Herrada, Antonio Ramos, and Jos e M. Montanero. Review on the physics of electrospray: From electrokinetics to the operating conditions of single and coaxial taylor cone-jets, and ac electrospray. Journal of Aerosol Science, 125:32{56, 2018. doi:10.1016/j.jaerosci.2018.05.002. [34] Henry Berry Garrett. The charging of spacecraft surfaces. Reviews of Geophysics, 19 (4):577{616, 1981. doi:10.1029/RG019i004p00577. [35] A. Gomez and W. Deng. Fundamentals of Cone-Jet Electrospray. 2011. doi:10.1002/9781118001684.ch20. cited By 7. [36] V. R. Gundabala, N. Vilanova, and A. Fern andez-Nieves. Current-voltage characteris- tic of electrospray processes in micro uidics. Phys. Rev. Lett., 105:154503, Oct 2010. doi:10.1103/PhysRevLett.105.154503. [37] Ji-Huan He, Yu-Qin Wan, and Jian-Yong Yu. Scaling law in electrospinning: relation- ship between electric current and solution ow rate. Polymer, 46(8):2799 { 2801, 2005. ISSN 0032-3861. doi:10.1016/j.polymer.2005.01.065. [38] C. D. HENDRICKS and J. HOGAN. Investigation of the charge-to-mass ratio of elec- trically sprayed liquid particles. AIAA Journal, 3(2):296{301, 1965. doi:10.2514/3.2845. 120 BIBLIOGRAPHY [39] R. W. Hockney and J. W. Eastwood. Computer simulation using particles. IOP Pub- lishing Ltd, 1988. [40] Yuan Hu and Joseph Wang. Fully kinetic simulations of collisionless mesothermal plasma emission: macroscopic plume structure and microscopic electron characteristics. Physics of Plasmas, 24:035510, 2017. doi:10.1063/1.4978484. [41] Victor E. Krohn Jr. Glycerol droplets for electrostatic propulsion. Electric Propul- sion Development, Progress in Astronautics and Aeronautics, pages 435{440, 1962. doi:10.2514/5.9781600864865.0435.0440. [42] Ha-Na Kown, Su-Jin Jang, Yun Chan Kang, and Kwang Chul Roh. The eect of ils as co-salts in electrolytes for high voltage supercapacitors. Scientic Reports, 9(1), 2019. doi:10.1038/s41598-018-37322-y. [43] Rakesh Kumar, Arnaud Bauner, Zheng Li, and Deborah A. Levin. Electrospray simu- lation in a colloid thruster using particle-in-cell (pic) approach. 50th AIAA Aerospace Sciences Meeting, AIAA 2012-0788, 2012. doi:10.2514/6.2012-788. [44] S. T. Lai. Spacecraft Charging. 2011. doi:10.2514/4.868375. [45] David B. Langmuir, Ernst Stuhlinger, and Jr. J. M. Sellen. Liquid metal droplets for heavy particle propulsion. Electrostatic Propulsion, Progress in Astronautics and Aeronautics, pages 73{80, 1961. doi:10.2514/5.9781600864797.0073.0080. [46] I. Langmuir. The eect of space charge and residual gases on thermionic currents in high vacuum. Phys. Rev., 2:450{486, 1913. doi:10.1103/PhysRev.2.450. [47] L. Liao and R. J. A. Hill. Shapes and ssility of highly charged and rapidly rotating levitated liquid drops. Physical Review Letters, 119:114501, 2017. doi:10.1103/PhysRevLett.119.114501. [48] Hui Liu, Boying Wu, Daren Yu, Yong Cao, and Ping Duan. Particle-in-cell simu- lation of a hall thruster. Journal of Physics D: Applied Physics, 43:165202, 2010. doi:10.1088/0022-3727/43/16/165202. [49] Paulo Lozano. Studies on the ion-droplet mixed regime in colloid thrusters. Ph.D. dis- sertation, Massachusetts Institute of Technology, Department of Aeronautics and As- tronautics, 2003. [50] Paulo Lozano and Manuel Mart nez-S anchez. Ionic liquid ion sources: suppression of electrochemical reactions using voltage alternation. Journal of Colloid and Interface Science, 280(1):149 { 154, 2004. ISSN 0021-9797. doi:10.1016/j.jcis.2004.07.037. [51] Manuel Martinez-Sanchez and Paulo Lozano. 16.522 space propulsion. Massachusetts Institute of Technology: MIT OpenCourseWare, Spring 2015. URL https://ocw.mit.edu. License: Creative Commons BY-NC-SA. BIBLIOGRAPHY 121 [52] Neil A. Mehta and Deborah A. Levin. Molecular dynamics electrospray simulations of coarse-grained ethylammonium nitrate (ean) and 1-ethyl-3-methylimidazolium tetra u- oroborate (emim-bf 4 ). Aerospace, 5:1, 2017. doi:10.3390/aerospace5010001. [53] Neil A. Mehta and Deborah A. Levin. Comparison of two protic ionic liquid behaviors in the presence of an electric eld using molecular dynamics. The Journal of Chemical Physics, 147:234505, 2017. doi:10.1063/1.5001827. [54] Neil A. Mehta and Deborah A. Levin. Sensitivity of electrospray molecular dynamics simulations to long-range coulomb interaction models. Physical Review E, 97:033306, 2018. doi:10.1103/PhysRevE.97.033306. [55] Neil A. Mehta and Deborah A. Levin. Electrospray molecular dynamics simulations using an octree-based coulomb interaction method. Physical Review E, 99:033302, 2019. doi:10.1103/PhysRevE.99.033302. [56] Fernando Mier-Hicks and Paulo C. Lozano. Spacecraft-charging characteristics induced by the operation of electrospray thrusters. Journal of Propulsion and Power, 33(2): 456{467, 2017. doi:10.2514/1.B36292. [57] M. Paine and S. Gabriel. A micro-fabricated colloidal thruster array. doi:10.2514/6.2001- 3329. [58] Robert J. Pfeifer and Charles D. Hendricks. Charge-to-mass relationships for electrohy- drodynamically sprayed liquid droplets. The Physics of Fluids, 10(10):2149{2154, 1967. doi:10.1063/1.1762011. [59] Boguslaw P. Pozniak and Richard B. Cole. Current measurements within the electro- spray emitter. Journal of the American Society for Mass Spectrometry, 18(4):737{748, Apr 2007. ISSN 1879-1123. doi:10.1016/j.jasms.2006.11.012. [60] M. Rahmanpour and R. Ebrahimi. Numerical simulation of highly charged droplets dynamics with primary break-up in electrospray at sub rayleigh limit. Journal of Applied Fluid Mechanics, 10:541{550, 2017. doi:10.18869/acadpub.jafm.73.238.26186. [61] Jerey Reichbach. Micropropulsion system selection for precision formation ying satel- lites. M.S. dissertation, Massachusetts Institute of Technology, Department of Aeronau- tics and Astronautics, 2001. [62] S. N. Reznik, A. L. Yarin, A. Theron, and E. Zussman. Transient and steady shapes of droplets attached to a surface in a strong electric eld. Journal of Fluid Mechanics, 516:349{377, 2004. doi:10.1017/S0022112004000679. [63] S. N. Reznik, A. L. Yarin, E. Zussman, and L. Bercovici. Evolution of a compound droplet attached to a core-shell nozzle under the action of a strong electric eld. Physics of Fluids, 18:062101, 2006. doi:10.1063/1.2206747. 122 BIBLIOGRAPHY [64] J. Rosell-Llompart, J. Grifoll, and I.G. Loscertales. Electrosprays in the cone-jet mode: From taylor cone formation to spray development. Journal of Aerosol Science, 125: 2{31, 2018. doi:10.1016/j.jaerosci.2018.04.008. cited By 12. [65] J.-F Roussel, Thomas Tondu, J.-C Mateo-Velez, Enrico Chesta, Stephane D'Escrivan, and Laurent Perraud. Modeling of feep plume eects on microscope space- craft. Plasma Science, IEEE Transactions on, 36:2378 { 2386, 11 2008. doi:10.1109/TPS.2008.2002541. [66] Lord Rayleigh F. R. S. Xx. on the equilibrium of liquid conducting masses charged with electricity. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 14:184{186, 1882. doi:10.1080/14786448208628425. [67] F. Taccogna, S. Longo, M. Capitelli, and R. Schneider. Particle-in-cell simulation of stationary plasma thruster. Contributions to plasma physics, 47:635{656, 2007. doi:10.1002/ctpp.200710074. [68] M. Tajmar and J. Wang. Three-dimensional numerical simulation of eld-emission- electric-propulsion neutralization. Journal of Propulsion and Power, 16:536{544, 2000. doi:10.2514/2.5602. [69] M. Tajmar and J. Wang. Three-dimensional numerical simulation of eld-emission- electric-propulsion back ow contamination. Journal of Propulsion and Power, 38:69{78, 2001. doi:10.2514/2.3656. [70] Sir Georey Taylor. Disintegration of water drops in an electric eld. Proceedings of the Royal Society A, 280:383{397, 1964. doi:10.1098/rspa.1964.0151. [71] Si Bui Quang Tran, Doyoung Byun, Vu Dat Nguyen, and Tae Sam Kang. Liquid menis- cus oscillation and drop ejection by ac voltage, pulsed dc voltage, and superimposing dc to ac voltages. Physical Review E, 80:026318, 2009. doi:10.1103/PhysRevE.80.026318. [72] Si Bui Quang Tran, Doyoung Byun, Hadi Teguh Yudistira, and Je Hoon Oh. Semi- analytical study of hemispherical meniscus oscillation with an anchored edge on a conductive at plate under an ac electric eld. Physics of Fluids, 23:022006, 2011. doi:10.1063/1.3556610. [73] R. J. Turnbull. On the instability of an electrostatically sprayed liquid jet. In Conference Record of the 1990 IEEE Industry Applications Society Annual Meeting, pages 783{788 vol.1, Oct 1990. doi:10.1109/IAS.1990.152275. [74] R. J. Umstattd, C. G. Carr, C. L. Frenzen, J. W. Luginsland, and Y. Y. Lau. A simple physical derivation of child{langmuir space-charge-limited emission using vacuum capacitance. American Journal of Physics, 73(2):160{163, 2005. doi:10.1119/1.1781664. [75] L. F. Velasquez-Garcia, A. I. Akinwande, and M. Martinez-Sanchez. A planar array of micro-fabricated electrospray emitters for thruster applications. Journal of Microelectromechanical Systems, 15(5):1272{1280, Oct 2006. ISSN 1057-7157. doi:10.1109/JMEMS.2006.879710. BIBLIOGRAPHY 123 [76] J. P. Verboncoeur. Symmetric spline weighting for charge and current den- sity in particle simulation. Journal of Computational Physics, 174:421{427, 2001. doi:10.1006/jcph.2001.6923. [77] Loup Verlet. Computer \experiments" on classical uids. i. thermodynamical properties of lennard-jones molecules. Physical Review, 159:98, 1967. doi:10.1103/PhysRev.159.98. [78] Joseph Wang, Ouliang Chang, and Yong Cao. Electron-ion coupling in mesother- mal plasma beam emission. IEEE Transactions on Plasma Science, 40:230{236, 2012. doi:10.1109/TPS.2011.2179066. [79] Pengxiang Wang, Arnaud Borner, Zheng Li, and Deborah Levin. An advanced particle- in-cell (pic) approach for electrospray simulation in a colloid thruster using molecular dynamics simulation results. 43rd AIAA Thermophysics Conference, AIAA 2012-2993, 2012. doi:10.2514/6.2012-2993. [80] Pengxiang Wang, Arnaud Borner, Burak Korkut, Zheng Li, and Deborah Levin. Sim- ulations of electrospray in a colloid thruster with a high resolution three-dimensional particle-in-cell model. 44th AIAA Plasmadynamics and Lasers Conference, AIAA 2013- 2629, 2013. doi:10.2514/6.2013-2629. [81] A. L. Yarin, S. Koombhongse, and D. H. Reneker. Taylor cone and jetting from liquid droplets in electrospinning of nanobers. Journal of Applied Physics, 90(9):4836{4846, 2001. doi:10.1063/1.1408260. [82] Yasushige Yonezawa. A long-range electrostatic potential based on the wolf method charge-neutral condition. The Journal of Chemical Physics, 136(24):244103, 2012. doi:10.1063/1.4729748. [83] Sidney Zafran, John C. Beynon, Philip W. Kidd, Haywood Shelton, and Frederick A. Jackson. One-mlb colloid thruster system development. Journal of Spacecraft and Rockets - J SPACECRAFT ROCKET, 10, 08 1973. doi:10.2514/3.61920. [84] Yinjian Zhao. Investigation of eective impact parameters in electron-ion temperature relaxation via particle-particle coulombic molecular dynamics. Physics Letters A, 381: 2944{2948, 2017. doi:10.1016/j.physleta.2017.05.066. [85] Yinjian Zhao. A binary collision monte carlo model for temperature relaxation in mul- ticomponent plasmas. AIP Advances, 8:075016, 2018. doi:10.1063/1.5035362. [86] Yinjian Zhao. A binary collision monte carlo model for electron-ion temperature relax- ation. Plysics of Plasmas, 25:032707, 2018. doi:10.1063/1.5025581. [87] Yinjian Zhao. Three-dimensional particle-particle simulations: Dependence of relaxation time on plasma parameter. Plysics of Plasmas, 25:052112, 2018. doi:10.1063/1.5025431. 124 BIBLIOGRAPHY [88] Yinjian Zhao, Hui Liu, Daren Yu, Peng Hu, and Huan Wu. Particle-in-cell simulations for the eect of magnetic eld strength on a cusped eld thruster. Journal of Physics D: Applied Physics, 47:045201, 2014. doi:10.1088/0022-3727/47/4/045201. [89] Yinjian Zhao, Joseph Wang, and Hideyuki Usui. Simulations of ion thruster beam neutralization using a particle{particle model. Journal of Propulsion and Power, 34: 1109{1115, 2018. doi:10.2514/1.B36770. [90] J. K. Ziemer, T. M. Randolph, G. W. Franklin, V. Hruby, D. Spence, N. Demmons, T. Roy, E. Ehrbar, J. Zwahlen, R. Martin, and W. Connolly. Colloid micro-newton thrusters for the space technology 7 mission. In 2010 IEEE Aerospace Conference, pages 1{19, 2010. doi:10.1109/AERO.2010.5446760. [91] John Ziemer, Colleen Marrese-Reading, Charley Dunn, Andrew Romero-Wolf, Curt Cutler, Shahram Javidnia, Thanh Li, Irena Li, Garth Franklin, and Phil Barela. Colloid microthruster ight performance results from space technology 7 disturbance reduction system. In International Electric Propulsion Conference (IEPC), number 20170010216, 2017. [92] John K. Ziemer, Thomas Randolph, Vlad Hruby, Douglas Spence, Nathaniel Demmons, Tom Roy, William Connolly, Eric Ehrbar, Jurg Zwahlen, and Roy Martin. Colloid mi- crothrust propulsion for the space technology 7 (st7) and lisa missions. AIP Conference Proceedings, 873(1):548{555, 2006. doi:10.1063/1.2405097.
Abstract (if available)
Abstract
Colloid thrusters have been widely utilized on small spacecrafts, because of their high thrust density, specific impulse, efficiency, and accurate attitude control over long periods of time. Although the physics of an electrified liquid surface has been strongly studied, the knowledge of the colloid thruster beam is very limited due to the lack of precise experimental measurements and numerical simulations. In this dissertation, state-of-the-art particle simulation models of the colloid thruster beam are developed: The Particle-Particle (PP) model aims at a high accuracy, such that its results can be treated as a baseline to provide guidance for other models
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Numerical and experimental investigations of ionic electrospray thruster plume
PDF
The cathode plasma simulation
PDF
Kinetic studies of collisionless mesothermal plasma flow dynamics
PDF
Spacecraft trajectory optimization: multiple-impulse to time-optimal finite-burn conversion
PDF
Three-dimensional kinetic simulations of whistler turbulence in solar wind on parallel supercomputers
PDF
Numerical and experimental investigations of dust-plasma-asteroid interactions
PDF
The space environment near the Sun and its effects on Parker Solar Probe spacecraft and FIELDS instrumentation
PDF
Development of the two-stage micro pulsed plasma thruster
PDF
Grid-based Vlasov method for kinetic plasma simulations
PDF
An investigation into magnetically induced extractor-less electrospray propulsion devices
PDF
Stability and folding rate of proteins and identification of their inhibitors
PDF
A study of ignition effects on thruster performance of a multi-electrode capillary discharge using visible emission spectroscopy diagnostics
PDF
Optimization of relative motion rendezvous via modified particle swarm and primer vector theory
PDF
An experimental investigation of cathode erosion in high current magnetoplasmadynamic arc discharges
PDF
Particle-in-cell simulations of kinetic-scale turbulence in the solar wind
PDF
Mission design study of an RTG powered, ion engine equipped interstellar spacecraft
PDF
The development of an autonomous subsystem reconfiguration algorithm for the guidance, navigation, and control of aggregated multi-satellite systems
PDF
A molecular dynamics study of interactions between the enzyme lysozyme and the photo-responsive surfactant azobenzene trimethylammonium bromide (azoTAB)
PDF
Dynamic modeling and simulation of flapping-wing micro air vehicles
PDF
Novel processing of liquid substrates via initiated chemical vapor depostion
Asset Metadata
Creator
Zhao, Yinjian
(author)
Core Title
Particle simulations of colloid thruster beam
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Astronautical Engineering
Publication Date
10/10/2019
Defense Date
08/30/2019
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
colloid thruster,OAI-PMH Harvest,particle simulation,particle-mesh model,particle-particle model
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Wang, Joseph (
committee chair
), Kunc, Joseph (
committee member
), Nakano, Aiichiro (
committee member
)
Creator Email
yin.work.0@gmail.com,yinjianz@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c89-223626
Unique identifier
UC11674727
Identifier
etd-ZhaoYinjia-7852.pdf (filename),usctheses-c89-223626 (legacy record id)
Legacy Identifier
etd-ZhaoYinjia-7852.pdf
Dmrecord
223626
Document Type
Dissertation
Rights
Zhao, Yinjian
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
colloid thruster
particle simulation
particle-mesh model
particle-particle model