# Surface Erosion

While it may not seem like a single particle will do much damage, when it is moving at 1 MeV and there are 10^20 of them hitting your spacecraft every second, then you might be having a bad day. Every time a particle impacts a surface there is a chance that it will knock off or “sputter” some material from the surface. Over time this can cause noticeable degradation in high flux areas such as Hall thruster channels, electrostatic acceleration grids, etc. Even structures far downstream of a thruster can erode over time if they are near the main plume. This article explains the basic algorithms involved in determining how a surface will erode when exposed to a given flux.

We can break the erosion algorithm down into three steps:

- Determine incoming flux to a surface element
- Calculate amount of sputtered material
- Move each node associated with the surface element

## Determine Incoming Flux to a Surface Element

Calculating the actual flux of something like a thruster is out of the scope of this article, but we can still do a *little* work. Given a certain uniform incoming flux into the domain, we need to get the cross sectional area of a surface element that sees that flux. You might be thinking that this is just a simple dot product, and you would be correct. We simply dot the surface normal **unit vector **with the flux **unit vector **to give us an area ratio. This ratio multiplied by the element area gives us the cross sectional area visible to the flux. We then get the overall flux to the surface element. I’m sure you are saying “if you had not made the flux a unit vector and just left it with a magnitude then you wouldn’t have to multiply by flux at the end, it would be just the dot product.” That is true, but it is a little easier to visualize if we take small steps!

$$\vec{F_{visible}}=\vec{F_{total}}\cdot \hat{N}$$

## Calculating Amount of Sputtered Material

*Sputter yield* is essentially the number of surface molecules kicked off for each incident molecule. There are a number of algorithms out there to calculate sputter yield and that I am not going to get into. What I *will* say is that most of them typically take the form of: one function to determine the sputter yield at normal incidence and another to adjust that for incidence angle.

$$Y=f(\text{materials}, \text{flow params})g(\theta, \text{materials})$$

What this tells us is that for a given input flux, F, we will be sputtering away X amount of material. While some models can be extremely complex, a simple cosine model goes a long way – choose a max energy the normal incidence, and a max angle (this is typically NOT at normal incidence, for a first cut 45 degrees is not a bad guess for the peak.) So for our simple model, we may have something like:

$$Y=\cos(E/E_{peak})\cos(\theta + \dfrac{\pi}{2})$$

The peak energy varies wildly by material, but just to try things out, maybe put a few thousand eV in there. If you are really cool, you can model where this sputtered material goes (*“**please don’t land on the optics!”)*, but that is a task for another day.

## Move each node associated with the surface element

Moving the nodes is the hard part. In reality, each surface interaction only knocks out a molecule or two of material, but we need our sim to run fast so we take off big chunks at a time. And of course by “take off big chunks” what we really mean is move the surface nodes to a new position. So how do we know where to move the nodes?

Let’s take our circle and introduce a uniform flux coming in from the top as shown in Figure 2. To illustrate the effects on the surface we are going to use a laughably coarse mesh, which you would never use in real life, right?

Typically we have small panels and small timesteps so it is OK to just loop through the panels and erode them one by one. For this example we want to erode the surface panel highlighted in green. We know how much material we want to remove, we know the density of the material, and we know the area of the panel, so it should be fairly straightforward to say “OK, my density is 1000 kg/m^3, I want to take out 1e-4 kg, and my area is 5e-2 m^2, soooooooo, if I move the panel in 2e-6 m then I should be all set!” The hard part is what do we mean when we say “move the panel in?” Figure 3 outlines some possible scenarios.

The first place most people would go is to just move along the surface normal. This is shown in Case A. For a very small amount of erosion this would probably be OK, but in our super exaggerated example you can see that the underside of the circle begins to lift up, implying that there is “erosion” on the bottom surface where this is no impingement.

So then you might think a beam of particles is sort of like a drill so maybe we should erode in the direction of the flux. This case is shown in Case B. This makes physical sense, but as shown here, if we move the nodes in the direction of the flux we can actually *add* material to the surface, which is right out.

OK, so maybe we can move surface nodes along the plane of a neighboring surface. This way we are removing material from *this* element without affecting any *other* elements. This technique is shown in Case C. At first this seems like a super deal, but what about the case where we have a flat plate? The nodes would simply slide around on the surface and no erosion would actually occur. HMMMMMMM.

What we really need is a combination of methods B and C. This is shown in Figure 4. We basically want to move a node in the flux direction until it hits a surface, and then move it along the surface. As with most problems in life, this one is again solved by dot products. If the surface normal has a component *opposite* the flux direction, then that surface is OK to move in the flux direction. But remember, we don’t generally move *surfaces*, we move *nodes*. So, when we are moving a node you have to check the normal of both surfaces connected to that node. If they both have a component opposite the flux, then we can move in the flux dir, such as in point A. For point B, however, the yellow arrow dots positively with the flux (has no opposite component), so we have to move point B down the relevant surface.

Now would be a good time to discuss the inherent conflict here between nodes and elements. It is easy to say a particle hits element X, which has an area of Z. It is also easy to move node A in direction B. As we saw above, however, it is a little tougher to have particles hit an *element* and then translate the *nodes*. One way to deal with this is to do everything relative to the nodes. What this means is that instead of looping over elements, we actually want to loop through nodes. We say that the “area” of the node, is just half the area of each neighboring element. Calculating both elements is really not that much more work, and it is nicer in a few ways:

- Each node only gets moved once. When we loop through panels, each node gets moved twice, which can get confusing quickly.
- Looping through nodes is much more natural as we already often loop through nodes for other reasons (potential for example)
- If we loop through
*elements*, the normal vector of the element remains ~ constant. This may seem like a good thing, but it restricts the surface from reshaping itself which can lead to some very large or small elements. Looping through*nodes*gives the mesh a little more freedom to adapt to the imposed shape changes.

Now that we know the direction to move the nodes how far do we move them? We can go back to that same simple thought from before about the density and the volume of material we want to remove. Let’s firm that up a bit:

$latex dx=dfrac{m}{rho Delta y Delta z}$

Our case here is 2D so we can throw out $latex Delta z$.

$latex dx=dfrac{m}{rho Delta y}$

And that is that….. if you are looping through surfaces. This is the one area where looping through the nodes is actually a little tricky as we are no longer just eroding out a rectangle, we are changing the shape of the mesh. Who knows how much volume we are taking out?

As it turns out, it is not as complicated as it might seem. We know the position of all the nodes so you can just look at the node of interest and the neighboring two nodes and treat those three points as a triangle. There are a few ways we can handle this.

- Actually calculate the area of the two triangles and solve for the distance to move. This is the most correct, but is a bit of a pain.
- Assume that, for calculating the distance, the erosion is normal to the triangle baseline. In this case we can just use the change in area of moving a node.$latex Delta A = dfrac{1}{2}bh_1 – dfrac{1}{2}bh_2$
- Since moving this node removes material from
*two*panels, we need to add their areas together when deciding how far to move the node.$latex Delta h = dfrac{2Delta A}{b_1+b_2}$ - If you would like to account for the fact that the erosion direction is not strictly along the “h” vector, that is easily done via the angle between the baseline (vector between baseline and the flux direction, as shown in Figure 5.)$latex Delta x = dfrac{Delta h}{sintheta}$

Now that we know how to handle a single erosion event, all we have to do is throw in a couple for loops (elements and time) and start eroding!

### Matlab Example

Here is a quick Matlab erosion routine you can try out yourself. It has some default values for all the material stuff and all you specify are the number of sides on the circle and the time. Let’s look at two example runs of this code.

First, we run with our laughably coarse 7 panel circle, run for 2e6 seconds. The result is shown in Figure 7. Note that the nodes all move in the flux direction except for the lower left node, which was forced to move along the surface. Other than this, however, it is tough to really get much out of this since so few panels were involved. It is difficult to say where the max erosion occurred or what shape the final rod will take.

Let’s look at Figure 8, which shows the same 2e6 seconds, but now with 50 panels. Here, the rod is much more symmetrical and now has a definite shape to it – we see that the rod has started to become less circular and more pointy. This makes sense because of our cosine angular profile. Putting in more complex flux distributions or geometries can yield some interesting results!

Here is an animation showing the progression of the erosion. Notice how the top of the circle tends to stay flat because this is a single element. This may be somewhat unrealistic, and could be alleviated by mesh refinement, or element splitting, which will be covered in a separate tutorial.

For more on this, checkout the new online sputter calculator.

Sputter Calculator (alpha version)

Detailed Surface Model

Modeling Diffuse Reflection (or How to Sample Cosine Distribution)

**Bio:**Alexander Barrie works for the heliophysics branch at NASA/GSFC. He has a masters degree from Virginia Tech in aerospace engineering with a thesis on differential spacecraft charging. He has extensive experience modeling plasma flows for electric propulsion applications with the Air Force Research Lab as well as charged particle and rarified gas modeling for NASA. He is fluent in C/C++, JAVA, Matlab, and Python and has developed algorithms and codes for spacecraft charging, PIC, surface erosion, DSMC, and material sputtering. He is currently working on the Fast Plasma Investigation for the Magnetosperic Multiscale mission modeling and testing the performance of the Dual Electron Spectrometer and the Instrument Data Processing Unit. In his spare time Alex enjoys rock climbing, mountain biking, and kayaking, among other outdoor adventures.

Thanks again PIC for another insightful article. I have two questions; here is the setup: when two adjacent overhead powerlines come close to each other (due to wind, for example), electricity will arc over from one powerline to the other powerline. My questions are: 1) can the arc be considered plasma? and 2) can the surface erosion on the powerlines due to the arc be calculated as in this article? Thank you kindly for your time.

Derek,

You could consider the immediate area of the arc a plasma since the atmosphere will likely be ionized, however because of the high density (1 atm) you can’t really model it with a particle approach and would need to use a fluid model. The plasma would also recombine very quickly. The surface erosion can’t really be modeled using this approach because an arc on a surface is burning away material, as opposed to knocking it off via particulate impacts.