Work Description

INTRODUCTION:
This deposity serves as the record of data and analysis for our manuscript "Robust Design of a Porous Electrospray Array Thruster". Included in the deposit are all of the raw data used in the study and every piece of analysis code necessary to reproduce the published results, as well as additional materials (further plots, trade studies, etc.) than could be accommodated in the publication. These extra plots should be included in the supplementary materials of linked publications, but are obviously also available upon request.

We prioritize this level of transparency and reproducibility in our results, and hope you do, too. Designing robust electrospray array thrusters under uncertainty constitutes a principal challenge in enabling this potentially game-changing technology.

Regards,
Collin B. Whittaker
[email protected]

Last updated: 24 Feb 2025

RESEARCH OVERVIEW:
This study follows after work conducted first for my dissertation and is presently being prepared for journal submission. The goal of our analysis was to analyze a small design space for an electrospray array thruster---varying the geometry of its emitters, the size of its extractor apertures, and its operating voltage---to determine designs robust to uncertainty. That is, we use a model for array performance whose input parameters we treat as uncertain (stemming from approximations to higher-order physics, manufacturing tolerances in fabricating a thruster, and so on). Making these predictions as a function of design, then, we can identify configurations that are performant robust to this uncertainty (i.e., still meet required performance specifications with high confidence).

The data which inform this trade study are taken pricipally from our pending manuscript "Emitter Model Inference from Electrospray Array Thruster Tests", and from my thesis, "Designing Porous Electrospray Array Thrusters Under Uncertainty" (linked to the dataset as published).
The analysis was conducted in January and February of 2025.
This work was supported by a NASA Space Technology Graduate Research Opportunity (80NSSC21K1247).
This research was also supported in part through computational resources and services provided by Advanced Research Computing, a division of Information and Technology Services at the University of Michigan, Ann Arbor.

METHODS:
At a high level, analysis is straightforward: We consider an array of different possible designs, compute probabilistic predictions for each design that account for different sources of uncertainty, and then we compute the optimal design over the space for various objective functions (i.e., maximize this, minimize that, and so on). Practically, though, there are many fiddly bits, and hence this repository.

The first component of analysis is our predictive model itself, which we derive primarily from:
C. B. Whittaker and B. A. Jorns, Journal of Applied Physics 134, 083301 (2023).
This model represents the emission current of a porous conical type electrospray emitter as the sum over individual menisci on the surface of that emitter which may deform into Taylor cone structures and begin extracting charged particles from the electrified fluid. To model the interception of this beam by the extraction electrode for the emitter, we extend this framework to predict divergence in the beam accounting for the position-dependent formation of emission sites on the emitter surface. Consequently, we can compute the performance of an array of emitters as the sum over its constituent elements. Crucially, we model the operational lifetime of an array is inversely proportional to the maximum current intercepted by any of its emitters (because this emitter is most likely to fail first as a result of a short induced by propellant accumulation).

Having formulated this model, we then itemize and quantify sources of uncertainty in the fabrication, assembly, and operation of each would-be design (e.g., uncertainty in the temperature of the propellant, manufacturing tolerances in the emitters, etc.). This uncertainty in the inputs to the model propagates to uncertainty in its predictions over the array thruster's performance. Because the performance its uncertain (i.e., described by a probability distribution), objectives over the performance take the form of statistical statements over this distribution (e.g., maximize the average thrust over the distribution, minimize the standard deviation in specific impulse over the distribution, etc.). As such, we compute which of the discrete designs are optimal over several of these objectives.

ANALYSIS FLOW:
Below is an overview of the analysis procedure, including the order the various scripts should be run to reproduce the results.

As a note, I have included only the data necessary to initiate analysis. Over the course of evaluating the scripts in this analysis flow, many GB of data will be generated, to which effect I have inserted comments below. If you want some intermediate or final product of the analysis for comparison, please reach out.
-----
Data/geoms.data
- starting place for analysis
- this file is generated out of another code base I maintain oriented toward actually manufacturing the emitters and extractors (i.e., procedurally generating machine G-code); access to this code base can be granted given proper licensure and upon reasonable request
- it contains info about the position and number of emitters on each design of the chip, which for the most part is only dictated by the aperture radius
- however, the structure of the file (groups over base radius nested within groups over aperture radius) will be propagated through analysis

VVVVV

Code/superGeomDistCalc.m
- this MATLAB script records the heuristic methods by which we arrived at the hierarchical distribution over emitter geometry distributions (i.e., the distribution over tools)
- it simply prints the results of analysis to terminal, which have been mechanically transcribed to "Code/geometry_sampling.py", so running this script isn't strictly necessary
- it does serve as a record of analysis, however

VVVVV

Code/geom_sampling.py
- called as is (run in terminal)
- Python script (refer to Python spec file) that encodes the hierarchical information from "superGeomDistCalc.m" and samples over the geometry for each design for each realization of the geometry
- adds its results into the "geoms.data" file

VVVVV

Code/extractorGeomCompute.m
- run as is
- MATLAB script that interpolates samples over the reference sites for the extractor geoemtry across the rest of the array for each realization
- apologies for the transfers back and forth from Python and MATLAB; some things are easier in one vs. the other
- adds its results to "geoms.data"

VVVVV

Code/interception_intersection_compute.py
- called as is
- once the extractor geometry is determined, we can then compute for each realization whether there is an intersection between the emitter and extractor and what the critical position (normal to emitter surface) is where interception first occurs
- adds its results to "geoms.data"

VVVVV

Code/fieldCompute.m
- called over a SLURM cluster as a job array (see also "Code/HPC_fieldCompute.sh")
- for each design (different base radius/aperture radius pair) for each realization, runs an electrostatic FEM model to compute the field strength of that emitter as a function of voltage
- exports to "fields_X.data" data structures depending on design
- requires "emitterSE.m" to be on the MATLAB path

VVVVV

Code/field_consol.py
- called with target directory as an argument (command line compatible)
- consolidates individual "fields_X.data" files into a singular "fields.data" file with all the field profiles

VVVVV

Code/paramssurfssubspropsbeams_sampling.py
- called as is
- creates several different HDF5 files that contain samples of other input parameters to the model (i.e., those defined at the array level and not at the individual emitter level)
- samples are drawn from a prior distribution specified over the parameters
- for scaling coefficients that need to be trained from data, drawn from previous inference over the model from data on an actual system (see links)

VVVVV

Code/opt_predicts.py
- called over a SLURM cluster as a job array (see also "Code/HPC_opt_predicts.sh")
- this is the main script for analysis
- for evaluates the model for each design for each sample for each realization for each emitter, some trillions of model calls
- exports to "predicts_X.data" HDF5 files for each different design

VVVVV

Code/predict_consol.py
- called with target directory/file as an argument (command line)
- consolidates individual "predicts_X.data" files into a singular "predicts.data" file needed by plotting scripts, etc.

VVVVV

Code/time_predicts.py
- called on a SLURM cluster (see also "Code/HPC_time_predicts.sh")
- for a single design (maximum thrust), computes time-dependent performance predictions assuming emitters fail individually rather than array failing collectively

VVVVV

Code/Plotting
- called as is
- various plotting scripts to render results
- congrats, analysis completed
-----
FILE INVENTORY:
Below appears a description of all the files included here.
--------------------
Python-spec-file.txt

A spec file for conda that can be used to pull the necessary packages in a new conda environment by running:

conda create --name myenv --file Python-spec-file.txt

This includes all dependencies, the most important of which is that since the code relies on MPI to parallelize jobs, you need a parallel build of HDF5.
For MATLAB codes, using whatever version is called for in the respective batch script should suffice.
--------------------
README.txt

This readme file.
--------------------
Data/

Directory including raw data files necessary to initiate analysis and serving as the target for exporting much data.
--------------------
Data/geoms.data

HDF5 file serving as a starting point for analysis over the emitter geometries of the chip. Additional datasets are added to this chip by various other scripts.
--------------------
Data/EMIM_IM_*.csv

These files are simple tables of temperature-dependent propellant properties for EMIM-IM. Used to generate samples for the properties subject to an uncertainty on the temperature.
--------------------
Data/geoms_old.data

HDF5 file containing information from previous inference done over emitter geometry for a couple chips. Used to inform distribution over tools for producing emitter geometry. That is, this contains training data for a heuristic inference over the tool distribution.
--------------------
Data/onset_params_old.data

HDF5 file containing samples of the prior over parameters of the onset criterion model for menisci. Based on prior inference based on simulations of individual menisci dynamics. Refer to "paramssurfssubspropsbeams_sampling.py".
--------------------
Data/params_old.data

HDF5 file containing samples of the prior over the current model parameters learned from previous studies. Refer to "paramssurfssubspropsbeams_sampling.py". The system these were trained on is located within the middle of the design space.
--------------------
Data/ExtractorFrame_ChipCut_new.png

Binary image file (saved as .png) used to mask areas where the extractor chip superstructure occludes emitters that would be underneath. Not used in any scripts aside from plotting ones.
--------------------
Code/

This directory contains all of the analysis code used in the study.
--------------------
Code/EmissionModel/current_model.py

This module implements the multi-site current model at the individual emitter level.
--------------------
Code/EmissionModel/__init__.py

This is just so that Python recognizes the directory 'Code/EmissionModel' as a package. I have a broader code base with some other utilities/models in here for other projects (e.g., fitting geometry), so this has pulled only the relevant code.
--------------------
Code/Sampling/MHDRAM.py

A module with a slew of helpful functions for performing inference, generating samples, etc. Used frequently throughout the study.
--------------------
Code/Sampling/__init__.py

Same story as for the EmissionModel package.
--------------------
Code/superGeomDistCalc.m

MATLAB script that computes hierarchical distributions over emitter geometry distribution parameters (i.e., the mean and standard deviations in the means and standard deviations over emitter geometry). Prints results to terminal. Important to note that we have exercised some scientific judgement in which datasets to count for which parameters, on the basis of conditions of the chips, which we detail in the documentation.
--------------------
Code/geom_sampling.py

Python script that samples from the hierarchical distributions over emitter geometries to populate each possible design of the array (different choice of aperture radius and average base radius) with potential realizations of the geometry (samples) which we evaluate the model over.
--------------------
Code/extractorGeomCompute.m

MATLAB script that completes sampling over the emitter and extractor geometry by interpolating all the sites in the array from "reference sites" for the extractor (individual sites within the array that we model as controlling the otherwise smooth variation in extractor geometry). Transfers into MATLAB to utilize their implementation of a natural neighbor method (smooth indeed).
--------------------
Code/interception_intersection_compute.py

Python script that computes properties of interception and intersection over all emitters of all realizations of all designs. In particular, computes the critical arc position (position down the cross section of the emitter) below which all spray is predicted to intercept the extractor. Also computes whether there is an intersection between an emitter geometry with its extractor (i.e., they overlap and would be unphysical), indicating an infeasible design.
--------------------
Code/fieldCompute.m

MATLAB function that computes all the field simulations for each emitter for each realization of a single design. Represents a key bottleneck in the study and would benefit from surrogate modeling. Requires "emitterSE.m" to be on the MATLAB path. Designed to be run over a SLURM job array to exploit high performance computing resources. Functionality counts on the presence of the job array to divide by design and detect which dataset to analyze.
--------------------
Code/emitterSE.m

MATLAB function that encapsulates the process of computing the field profile for an individual emitter. Contains some code toward an abortive attempt to motivate an alternative divergence model that does individual particle tracing assuming particles are emitted from the surface of the emitter. Could be used to capture collimation of the beam by the externally applied field, which we have not considered using our simple line-of-sight model (though for many emitter geometries the discrepancy is slim).
--------------------
Code/HPC_field_compute.sh

An example batch script to run "fieldCompute.m" on a SLURM cluster. YMMV on your HPC system. It's always an adventure getting things set up to run on a different platform. At a minimum, your cluster configuration is likely to vary and hence your resource allocations, so should be modified to something usable.
--------------------
Code/field_consol.py

Command-line compatible Python script designed to coalesce "fieldCompute.m" calls over all designs made over a job array into a single data file. Recommended to call on your HPC resource then migrate down rather than vice versa (pretty disk write hungry and it's nice to have some fast scratch/temp storage to work on).
--------------------
Code/paramssurfssubspropsbeams_sampling.py

Python script to draw samples from prior distribution over array-level model parameters to store and prepare for evaluations. Notably, draws from previous inference work to train the emitter model on real array thruster tests (see references).
--------------------
Code/opt_predicts.py

Python script that represents the main analysis. Also designed to be called over a job array (see "Code/HPC_opt_predicts.py") where each job in the array targets a different combination of aperture radius and emitter base radius. Note that this code relies on MPI (specifically mpi4py's mpi4py.futures module), and is---frankly---likely to be quite buggy in migrating to a different architecture (i.e., your cluster may have a different distribution/version of MPI available, may not interact with the intercommunicator functionality in the same way, etc.).

Generates a bunch of different individual prediction files of from "predicts_X.data" which contain the important results. Note that we make more predictions here than are actually plotted/included in the manuscript, including some predictions over individual emitter properties (variance in emitter performance as a function of position in the array, quantiles over the population of emitters for each sample for each realization, etc.).

Importantly, this script utilizes temporary files and has hard-coded directories to which to put out data, which will need to be changed (unless you also have access to a "/scratch/bjorns_root/bjorns_0/cbwhitt" directory on your HPC cluster, which I doubt. Due to the size (many GBs) of data and the nature of the computations (over many emitters and so on), it is necessary to dump many of the intermediate results to disk before then reconsolidating at the end. This is partially a concession to memory limits (i.e., our cluster has a limit of 5 GB requested per CPU).
--------------------
Code/HPC_opt_predicts.py

SLURM batch script implementation to call "opt_predicts.py". Module loading statements will need to be changed based on your local configuration and whatnot to make sure that your HPC environment has access to the code you need to run it.
--------------------
Code/predict_consol.py

Python script to consolidate individual predictions, much like "field_consol.py".
--------------------
Code/time_predicts.py

Python script designed to be run on a SLURM cluster that does time-variable predictions for performance for a single design (in the included example, the highest-thrust case of Ra=175 um, rb = 150 um). Assumes that emitters fail individually instead of an array failing all at once. The need to make time-dependent emitter by emitter predictions over all samples and realizations defining our probabilistic predictions is the single most memory intensive thing we do in the study (end up with temp files a couple hundred GB in size, so be prepared to need that much remote storage).

As a side note, we sort of bungled the temp file implementation and neglected to chunk the data storage on disk to match the access pattern in the HDF5 file we would need, so this job is extremely disk hungry. This should be reimplemented if you don't want to waste money spinning CPUs while only one of them reads and writes to disk at the end (and waste time). We include the code in the state it is because that's how we ran the analysis (by the time we noticed that it was going too slowly and should have been chunked, it was cheaper to just let it keep running rather than reconfigure and restart the job. C'est la vie.
--------------------
Code/HPC_time_predicts.py

Example SLURM HPC batch script to run "time_predicts.py". See comments for earlier batch scripts.
--------------------
Code/Plotting/

The scripts in this directory are those that generate the plots used in the paper.
--------------------
Code/Plotting/emitterPackingPlot.m

Generates plots that show how emitters are arranged on the chip for each choice of aperture radius. Though only one is shown in the paper (for Ra=225), one for each different choice of pitch is generated here.
--------------------
Code/emitterPlots.m

Renders figure that shows example emitter geometries for each pair of aperture radius and mean base radius over the design space. Relies on the helper function "emitterPlot.m" to render a single cross section.
--------------------
Code/emitterPlot.m

A help function in MATLAB that generates a line plot showing an individual emitter subject to our spherically-capped cone model. I put in the effort to make it pretty adaptable (by which I mean it inherits functionality from the underlying line plot).
--------------------
Code/predictPlots.m

Generates prediction plots for a slew of different performance parameters in the gridded style presented in the manuscript. Generates prediction plots for many more parameters than included in the manuscript. If you want these or other additional plots without actually having to spend the money and time to regenerate them, I intend to include them in the supplementary materials of the linked publications (or you can just email me and I'll send them to you).
--------------------
Code/failedEmitterPlots.m

Generates plots showing credible intervals over the tip radius of which emitter fails first over the entire design space. Again, only one choice of aperture radius and base radius are included in the manuscript, but the remainder are generated here.
--------------------
Code/optimPlots.m

Computes the optimal designs over trade space and renders the plots that illustrate them. As before, there are more objectives considered here than we decided to focus on in the manuscript. An important note is that in documenting results we often refer to the idea of infinite "within machine precision". We were getting some round-off errors in integrating current curves to estimate transmitted and intercepted current, so some of the raw predictions are tens of orders of magnitude larger than otherwise predicted (because the computer is dividing terms of order machine precision, which for double floating point is about 1e-16).
--------------------
Code/timePredictPlots.m

Computes credible intervals over time-dependent performance and renders as the gridded plot that appears in the manuscript.