Behavior of neutral density and question on config parameters

[bboyeleven]

Hi, I am still learning the code and physics of Hall Thruster (I wrote an issue with the tutorial not long ago) and I have new questions !

1. I did try to study the sensitivity of simulation to the different parameters of the config, and I'm a bit concern regarding the impact of ion_temperature and neutral_temperature. After running several simulations and plotted them all together, it appears that the value change of plasma properties, performances and frequencies is very low. For the neutral_temperature, I expected at least a variation of the neutral density at the inlet. Here are the plots :

And it is the same for the ion_temperature, do you have an idea ?

2. I plotted many different plasma properties compared to data I get from
   2007_Hofer_HPHall2-mobility-model_AIAA-2007-5267.pdf and tried to match the better I could with HallThruster.jl. It is quite accurate for electron temperature, plasma potential, ion density, ion velocity (from another dataset) and electron mobility. But I have an issue with the neutral density, I can't make it match Hofer's data and apart from modifying the neutral_velocity at the inlet (which modify also the electron temperature) I have no other parameters to modify that properties. What do you think the inlet neutral density is so far from Hofer's calculation ?
   
   
   
   
   
   

3. I would like to understand the physics behind the parameter wall_loss_model, the default value in the code is wall_loss_model = WallSheath(HallThruster.BNSiO2, 1.0) for BNSiO2 it is the material, so I understand you entered ϵ_star = 40 and σ₀ = 0.54 which are physical properties of the material. But what of the second parameter ? In the source code it is defined as α and intervenes in the formula of the electron-wall collision frequency, but you don't specify what is it and what it changes when you modify it ? Is it simply a proportional value ?

4. Almost the same question for the parameter electron_plume_loss_scale = 1.0, where in the equations does it intervene ? And what should happen in the behavior of the plasma properties if I modify it ?

5. And last question (I have many 😅), I didn't find in the source code nor in the documentation how the discharge_current was calculated from the simulation. Because once again, the one from my simulation is far from Hofer's one. 7 to 4.5 A.

For information, here is my configuration (for the comparison) :

config = HallThruster.Config(
ncharge = 1,
discharge_voltage = 300u"V",
thruster = HallThruster.SPT_100,
domain = (0.0, 0.05),
anode_mass_flow_rate = 5.0u"mg/s",
propellant = Xenon,
neutral_velocity = 225.0u"m/s",
ion_temperature = 500.0u"K",
neutral_temperature = 750.0u"K",
cathode_Te = 2.0u"eV",
anode_Te = 2.5u"eV",
anom_model = HallThruster.TwoZoneBohm(0.035*(1/16), 1.0*(1/16)), # calibrated with Hofer document
background_pressure = 1.7e-5u"Torr",
wall_loss_model = HallThruster.WallSheath(HallThruster.BNSiO2, 1.5),
electron_plume_loss_scale = 1.0,
ion_wall_losses = true
)

Thanks in advance for your answer 🙏.

[archermarx]

Hey, thanks for your questions.

I think the neutral temperature doesn't do what you're expecting it to do. You're trying to modify the anode temperature in order to adjust the neutral velocity, correct? To do that, you would need to modify the neutral_velocity flag in the config file. That said, at present, the neutral temp option does nothing really, so I think I can modify it to work as you expect, as that seems sensible.

To your third and fourth question, the parameters alpha and electron_plume_loss_scale scales the magnitude of the wall losses. Effectively, they modify the assumed edge-to-center density ratio. By default, we assume that the density at the walls is ~1/2 of the channel-averaged density. When alpha is not one, the density ratio is then

density_at_walls = alpha * 0.5 * channel_density

So a higher value of alpha increases wall losses.

The electron plume loss scale is applied on top of this value, as in certain configurations you might want to adjust the plume losses separately from those in the channel. In general, 1.0 is a good default. In the plume, then, the effective "density at the walls" is

density_at_walls = alpha * electron_plume_loss_scale * 0.5 * center_density

Lastly, the discharge current is computed here, in update_electrons.jl

function _discharge_current(params)
    (;cache, Δz_edge, ϕ_L, ϕ_R, ncells, iteration) = params
    (;∇pe, μ, ne, ji, Vs, channel_area) = cache

    int1 = 0.0
    int2 = 0.0

    apply_drag = false & !params.config.LANDMARK & (iteration[] > 5)

    if (apply_drag) 
        (;νei, νen, νan, ui) = cache
    end

    @inbounds for i in 1:ncells-1
        Δz = Δz_edge[i]


        int1_1 = (ji[i]   / e / μ[i]   + ∇pe[i])   / ne[i]
        int1_2 = (ji[i+1] / e / μ[i+1] + ∇pe[i+1]) / ne[i+1]
        
        if (apply_drag)
            ion_drag_1 = ui[1, i  ] * (νei[i  ] + νan[i  ]) * me / e
            ion_drag_2 = ui[1, i+1] * (νei[i+1] + νan[i+1]) * me / e
            neutral_drag_1 = params.config.neutral_velocity * νen[i  ] * me / e
            neutral_drag_2 = params.config.neutral_velocity * νen[i+1] * me / e
            int1_1 -= ion_drag_1 + neutral_drag_1
            int1_2 -= ion_drag_2 + neutral_drag_2
        end

        int1 += 0.5 * Δz * (int1_1 + int1_2)

        int2_1 = inv(e * ne[i] * μ[i] * channel_area[i])
        int2_2 = inv(e * ne[i+1] * μ[i+1] * channel_area[i+1])

        int2 += 0.5 * Δz * (int2_1 + int2_2)
    end

    ΔV = ϕ_L - ϕ_R + Vs[]

    I = (ΔV + int1) / int2

    return I
end

Here, we integrate Ohm's law axially to express the discharge current in terms of the discharge voltage and other electron parameters. This allows us to ensure that it is constant axially. I have also seen that for SPT-100 simulations, we tend to over-predict the discharge current compared to 2D sims, especially using two-zone Bohm profiles. I don't have a good idea why, but I don't think it's this code specifically, as we match the Landmark benchmark discharge currents, all of which are 7-8 A for an SPT-100-like case. (table 1 here)

[archermarx]

I am working on modifying the neutral velocity to work as you expect and to update the documentation to better describe the wall loss scaling. Bear with me and I should have a change pushed shortly!

[archermarx]

Anothing thing I noticed. We prescribe a transition length between the inner and outer anomalous transport zones. By default, this is 1/5 of a channel length, but you may want to shorten it for your purposes. This parameter is "transition_length_m" in the config file. I might actually change this parameter to 1/10 by default since 1/5 seems a bit too long.

[bboyeleven]

Thank you for the quick answer ! (Unfortunately I was not available this weekend to answer and look at the modification in detail, sorry).
I will check everything you did implement and keep you in touch.

Nota bene : when I was looking at the changes, I noticed that wall_collision_freq was defined in config.md but never implemented in config.jl.

[archermarx]

Good catch! Fixed that.

[bboyeleven]

It seems that there is a unity problem with the change in the configuration.jl code, probably it comes from the kB constant used in the neutral velocity thermal formula which must be incorrectly defined, or the program does not understand where it comes from (from the package Unitful in the module HallThruster.jl).
When I run my code without the neutral_velocity implemented and only the neutral_temperature I get this error :
ERROR: DimensionError: m s^-1 and K^1/2 are not dimensionally compatible.
with ref to : @ HallThruster C:\Users\.julia\packages\HallThruster\muv2R\src\simulation\configuration.jl:128 and @ HallThruster C:\Users\.julia\packages\HallThruster\muv2R\src\simulation\configuration.jl:178
or
ERROR: Unitful.DimensionError(Unitful.FreeUnits{(Unitful.Unit{:Meter, Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}()}(0, 1//1), Unitful.Unit{:Second, Unitful.Dimensions{(Unitful.Dimension{:Time}(1//1),)}()}(0, -1//1)), Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1), Unitful.Dimension{:Time}(-1//1))}(), nothing}(), Unitful.FreeUnits{(Unitful.Unit{:Kelvin, Unitful.Dimensions{(Unitful.Dimension{:Temperature}(1//1),)}()}(0, 1//2),), Unitful.Dimensions{(Unitful.Dimension{:Temperature}(1//2),)}(), nothing}())
with the same ref in the code.

I tried using Unitful.kB but didn't work and tried to replace by neutral_velocity = ustrip(u"m/s", neutral_velocity) but didn't work as well. Any ideas ? I'm not familiar with the structure of creating package on Julia.

[archermarx]

Ok, I see. I'm working on a fix right now.

[archermarx]

New version (0.16.2) should fix things!

[bboyeleven]

Hello Thomas. Questions about how the neutral density is computed inside the package in initialization.jl source file and the relation to ionization rate.

Implementation

1. From which model does the nn_0 come from (inlet neutral density)? Now it is a function of ion_density, ion_velocity and neutral_velocity at the inlet as well as anode_mass_flow_rate.
   

2. In the code, the neutral density and ionization rate are computed independently. The ionization rate is extracted from external data (probably mesured from experimental data I guess) and the neutral density is a smoothed between the inlet neutral function and 1 % of it on a distance of L_ch / 6. Is it a correct approximation ? Since ionization rate and neutral density are related by the formula (with ζ the ionization coefficient) :
   
   Why not have implemented related neutral density and ionization rate ?

Remark :

In the landmark_rates.csv file, the second column is said to be the ionization rate in [1/m³.s] but after plotting it and comparing it to 2008_Goebel__foundamentals_electric propulsion_Ion_HET.pdf formula it appears to be the ionization coefficient (in [m³/s]) rather than the ionization rate :

[archermarx]

The code you've referenced is just the initial neutral velocity at timestep zero. During the remainder of the code, the neutral density is computed using the continuity equation:

$$ \frac{\partial \rho_n A}{\partial t} + u_n \frac{\partial \rho_n A}{\partial t} = -\sum_{j=1}^3 \rho_n n_{i,j} k_{iz,j}(T_e) $$

This is handled in simulation/update_heavy_species.jl and sourceterms.jl

At the left boundary, the neutral density is just

$$ \rho_n(z = 0) = \frac{\dot{m}_a}{u_n A(0)} -{\left. \frac{\rho_i u_i}{u_n}\right\vert}_0 $$

This is computed in boundaryconditions.jl

You're right about the units in landmark_rates.csv. I will update that file in the next version.

[archermarx]

Hey, just checking in to see how things are going. have your issues been resolved?

[bboyeleven]

Hey, thank you for the interest. Unfortunately, not really 😅 I tried many different things to match the ionization rate, but the neutral density was still too low and plasma density too high.

I switched to another thruster, the P5 (developed in your laboratory, if I well remember) because I found experimental data on the PEPL website (especially in 2001_haas_thesis.pdf). And now I am trying to match the data on the center of the channel for the 1.6kW case (R=12.5 mm) with HallThruster.jl, but I have issues with the original position of the cathode (not said in the document), the discharge current which is far too high compared to the data (between 7 and 15 A and is supposed to be 5.4 A for the 1.6kW case), the position of the peak of electron temperature, an early drop in the plasma potential resulting in the accumulation of plasma density ner the anode and the correct value to use for the anomalous model (I use a TwoZoneBohm).

Have you ever tried to simulate it in HallThruster.jl ?

EDIT : I saw in the plotting.jl file the definition of current_spectrum function and was wondering what it was ? Is it a function to calculate the breathing mode frequency or something like that ? It seems not finished since the package FFTW is not even in the dependencies of HallThruster.jl.

[archermarx]

I have not tried to simulate the P5 before. The current_spectrum function is obsolete at this point and will be removed/replaced in a future update.

Can you share some plots of your problems? I can see if there's anything that sticks out to me.

[bboyeleven]

Sure ! For the anomalous transport I used TwoZoneBohm(k_in*(1/16), k_out*(1/16)) with k_in = 0.5 and k_out = 1.0. discharge_voltage is set to 300 V and cathode_potential to -21 V.

• I adjusted the magnetic field to the one used in 2001_haas_thesis.pdf
• For the ion density, by setting the right parameters I can match the order of magnitude, but the peak is too far to the left of that of the blue curve. And I can't get that double peak that they get.
• For the electron temperature, it is the same problem but on the right side. I don't know which parameters to adjust to control the electron temperature (regarding the physics).
• Since the measure given in the document start 0.06 m from the anode, we don't know the exact position (distance from the anode) of the cathode. It is written in the article that they used -21 V for the cathode_potential, so as a first guess we could admit that the potential is supposed to decrease to -21 V while reaching the cathode, and it is why I supposed the position of the cathode to be at 0.0875 m from the anode.
• I also have issues with the beginning of the plasma potential curve. It is supposed to be more constant inside the discharge channel, and here it is not.

I think that most of my problems come from the anomalous modelization, I never know which function and values to use. I always have to try different values to adapt to the case I have. (I also didn't really understand how to use the MultiLogBohm function).

[archermarx]

Thanks for sharing. Looking through the Haas thesis, it looks like you're actually plotting measurements instead of simulation results. There are some simulation results in section 6.6, beginning at page 186. Those simulation results don't show a double peak.

For the electron temperature, it is the same problem but on the right side. I don't know which parameters to adjust to control the electron temperature (regarding the physics).

This will largely be controlled by the anomalous transport

I think that most of my problems come from the anomalous modelization, I never know which function and values to use. I always have to try different values to adapt to the case I have.

That's certainly true. This is the most challenging part of Hall thruster modeling, and the anomalous transport must always be tuned.

To match this data, I would recommend a MultiLogBohm model with three nodes Perhaps ([0.015, 0.02, 0.025], [0.5, 0.05, 1] * 1/16) or something like that. You'll need to iterate on that to match the data, though.

However, i would caution that probes like those employed in this work are known to be highly perturbative, and may cause the spatial location and magnitude of these variables to shift around, so matching these data exactly may be difficult. For example the double-peaked number density seems pretty suspect to me.

[bboyeleven]

That's certainly true. This is the most challenging part of Hall thruster modeling, and the anomalous transport must always be tuned.

Yes, I understood that since the anomalous phenomenon is not well describe by physics literature.

The curves are better after adjusting the anomalous model to a MultiLogBohm as you advised :

However, i would caution that probes like those employed in this work are known to be highly perturbative, and may cause the spatial location and magnitude of these variables to shift around, so matching these data exactly may be difficult.

I see what you mean because after reading section 6.6 of Haas thesis and looking at the plot, the region of peak electron temperature seems to always follow the region of potential drop. But the blue curves are not following this relation.

[bboyeleven]

Hello. Regarding the original issue, I have not managed to solve it and there are behavior of plasma properties that I don't understand yet :

While I think the anomalous model is not perfectly adjust to the problem, there are still behaviors that perhaps are not related to it, and more to other parameters. Especially for the plasma density and neutral density.
For the plasma density, the ionization rate is far too high while the neutral density is well-adjusted (with this angle on the curve but with good order of magnitude).

I can match the order of magnitude of the plasma density by changing the wall_loss_model and switch $\alpha$ to 6.5, but I'm not sure regarding physical reality of that assumption. And furthermore, the neutral density is acting weird then.

Is it, when trying to match data, mandatory to adjust the wall_loss_model ? This parameter is making the curves go crazy while varying $\alpha$. When you increase $\alpha$ the plasma density decrease at first but sometimes re-increase suddenly or drop to 0.

[archermarx]

Can you share plots both with and without adjusting the wall loss coefficient? to 18. heinäk. 2024 klo 10.44 ap. bboyeleven ***@***.***> kirjoitti:

…

Hello. Regarding the original issue, I have not managed to solve it and there are behavior of plasma properties that I don't understand yet : Figure_1.png (view on web) <https://github.com/user-attachments/assets/c5fb862c-5744-45fd-b7e3-b66522ddeac3> While I think the anomalous model is not perfectly adjust to the problem, there are still behaviors that perhaps are not related to it, and more to other parameters. Especially for the plasma density and neutral density. For the plasma density, the ionization rate is far too high while the neutral density is well-adjusted (with this angle on the curve but with good order of magnitude). I can match the order of magnitude of the plasma density by changing the wall_loss_model and switch $\alpha$ to 6.5, but I'm not sure regarding physical reality of that assumption. And furthermore, the neutral density is acting weird then. [image: Figure_1] <https://github.com/user-attachments/assets/c6acd5c5-fa9d-47da-b62c-3a57148cece8> Is it, when trying to match data, mandatory to adjust the wall_loss_model ? This parameter is making the curves go crazy while varying $\alpha$. When you increase $\alpha$ the plasma density decrease at first but sometimes re-increase suddenly or drop to 0. — Reply to this email directly, view it on GitHub <#136 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AOP3QGXNHNVPDVUK4K74NGDZM7IDHAVCNFSM6AAAAABKNQVHKWVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDEMZWG43DCMRQGY> . You are receiving this because you commented.Message ID: ***@***.***>

[bboyeleven]

It is already the case, between the first and second series of plots, only the wall_loss_coefficient $\alpha$ change. From 1.0 (default) to 6.5.
Or perhaps you want plots with several changes of $\alpha$ on the same graphic ?

[archermarx]

No that's fine! I'll give it some thought.

[bboyeleven]

Hello @archermarx. I have questions regarding the issue.
Does the wall_loss_model have an impact on ion_density and neutral_density ? It should have a direct impact on ion_density but does it have a direct or indirect influence on neutral_density ?

Wall temperature is not adjustable in HallThruster.jl ? It could be perhaps a clue to the issue (neutrals interaction with the wall).

Also for the discharge current we have those results :

• First series of plasma properties ($\alpha$ = 1.0) : discharge current = 6.8 A
• Second series : Wall loss model adaptation ($\alpha$ = 6.5) : discharge current = 9.38 A

The discharge current is always higher than it is expected. It looks strange to me that with the wall loss function adaption, the computed current almost doubles the expected value.

[archermarx]

Apologies for the delay on the response. I totally missed this.

1. Wall loss model does impact neutral density if ion wall losses are enabled. This is because ions are assumed to recombine at the walls.
2. Wall temperature is not modelled in the code because ions and neutrals are assumed isothermal, so it is not able to affect things.
3. That's strange, I would definitely expect increasing the wall losses to decrease the discharge current

Let me know if you've had any further issues!