Add functionality to find fusion reaction rates, cross sections, and energy generation rates #110
Comments
namurphy
added
plasmapy.particles
|
A review article: Solar fusion cross sections |
|
Is there any reason why the input is a string that is being parsed through here? Further development would be much more convenient if the function took in two lists: reactants and products (e.g. nuclear_reaction_energy("D + T --> alpha + n") would become nuclear_reaction_energy(reactants=[D, T], products=[alpha, n]). That way, if you were checking against a database of reactants to products to automate this process, you could easy check by lengths of the arrays => possible reactions. |
|
Also, I'm not an expert on this BY ANY MEANS but as far as the initial temperature you mentioned in the description - is that the ignition temperature? If so, shouldn't it default to that temperature depending on the reaction? I haven't taken a physics class in years but was reading more into this to tackle this issue: https://www.ems.psu.edu/~radovic/Chapter14.pdf |
@jessicalostinspace - This is a great idea! You're right that we should have this be able to take a list of reactants and a list of products. I put this as a string to maximize readability, but I was operating under the (incorrect!) assumption that people would only care about two or three reactions at a time. This does ignore nucleosynthesis in massive stars and supernovae which have to deal with a lot of reactions. Thank you for the code review! I wonder if the best way to do this would be to allow both the string reaction input (which would be most readable to those of us who care only about a couple of reactions) as well as your idea for the |
|
Oof, it took me a few minutes to figure out how to add an image in one of these comments! The temperature is not the ignition temperature but rather the temperature of the plasma at the reaction is occuring in, so the goal is to figure out how much each reaction is happening at different temperatures as in this plot. [Image description: Fusion reactivity as functions of temperature for deuterium-deuterium, deuterium-tritium, and deuterium-helium-3 reactions. The reactivity for each reaction increases with temperature until it reaches a maximum and then gradually declines. Licensed under CC BY 2.5 under the associated link.] |
|
@namurphy would you suggest opening up a ticket for the refactor? or just branching off on this ticket? It should probably be refactored before implementing the logic for reaction rates, etc. They should also be in two different PR's. There won't really be a need to parse a string so a lot of that code can be cleaned up and more succinct. Thank you for your feedback :) |
|
Branching off this ticket sounds fine to me, and I agree that two separate PRs are needed. I'm wondering if we should require that Again, thank you! |
Originally, the function nuclear_reaction_energy required a string containing the nuclear reaction as the sole argument. In issue PlasmaPy#110, @jessicalostinspace suggested that the inputs be lists of the reactants and products in the reaction which will make it easier to iterate over different reactions. This pull request implements this suggestion. I did require that reactants and products be keyword-only arguments because there is potential for the two lists to get mixed up, and the sign of the reaction is important. I decided to keep the old API which has the argument being a string representing the function, but with some ambivalence. The advantage of this is that it is much more human readable if we are doing a small number of reactions. The disadvantage is that it makes the implementation a bit more complex. I broke the code into different functions so that the flow of the main part of the function is smoother. Electrons and positrons are now allowed in the reactions since they are important for beta decays.
|
Another resource is SkyNet which has nucleosynthesis data for astrophysics (see Lippuner & Roberts 2017). This uses a BSD 3-clause license as well. |
namurphy
removed
the
plasmapy.particles
|
I've written a python library, https://github.com/swurzel/lawson-criterion-paper which contains all the code needed to calculate fusion cross sections and reactivities for the main fusion reactions (D-T, both D-D branches, D-3He, and p-11B), among other things (like making plots of achieved Lawson parameters, temperatures, and fusion triple products. At the APS-DPP meeting in Spokane I met a few PlasmaPy folks who suggested I open an issue to see if someone might be interested in incorporating this library into PlasmaPy. I found this existing issue from 2017 which is basically requesting the same functionality. I wrote this code in order to generate the plots for this paper "Progress toward fusion energy breakeven and gain as measured against the Lawson criterion" which is open access and available here, Below are are a few plots from that paper generated by this code. If anyone is interested in taking on this project I'd be happy to advise! |
|
@swurzel — that's awesome! It looks like the main changes would be incorporating I noticed in the repo that you mention that it's not optimized for speed and that it'd help performance to cache the intermediate results. A related possibility would be to use something like Cython to get closer to compiled speeds. Also, have you considered posting the current version of that repository to Zenodo? That way, you'd be able to get a DOI for it, which would in turn make it possible to cite it from PlasmaPy's documentation. Since Zenodo records have a metadata field for references, you could also reference and connect to the Lawson criterion paper. Thank you for following up on this! |
|
Yes, those sound like the right steps to be taken. Maybe even the first step is to incorporating Using Cython sounds like a good idea too. As for Zenodo, I had not heard of it but I'd be happy to post it there so we could get a DOI for it. I can take care of that shortly. |
|
I would be happy to help take care of this, I've been itching to get a feature in rather than just do maintenance work :) |
|
@StanczakDominik That's great! I'm at your disposal to answer any and all questions about the code. |
|
@namurphy Question I probably should know the answer to.. what do the primes denote in nuclear reactions? In ENDF, this reaction must be written But this reaction must be written So only some of the primes you have included above... Best answer I see is here: https://physics.stackexchange.com/questions/580154/what-does-the-prime-symbol-indicate-in-a-nuclear-reaction So is including them somewhat qualitative? |
|
It looks to me that those are single quotation marks and are artifacts of copying any pasting from python strings. The reactions listed inside the single quotation marks are the reactions themselves. |
Ah, I think you're right in Nick's post! However the primes as written in my comment above seem to be different - they are required in ENDF's notation in some places like this D + Li-6 --> 2He-4' So at least when searching ENDF: https://www-nds.iaea.org/exfor/endf.htm This reaction is valid: But this is not |
My best guess is that the prime here indicates that the alpha particles are in an excited state rather than the ground state. I don't have a good understanding beyond that, though. |
I think this is right - slightly unfortunate since we don't have a way to represent that energy easily in the notation we use for particles in PlasmaPy. When PlasmaPy users look for this reaction, they'll probably need to look for "D + Li-6 --> 2He-4" |
In laboratory fusion devices and in stellar astrophysics, it's important to be able to calculate quantities associated with different fusion reactions. We will need a way to calculate the following quantities as functions of temperatures and densities:
We have some of the capability already in the atomic package, including
nuclear_binding_energy,energy_from_nuclear_reaction(which might be worth renamingnuclear_reaction_energy), andhalf_life. In particular,energy_from_nuclear_reactionhas code that parses reactions. Some of the reactions will be radioactive decay, so we will need the capability to calculate rates for that too.It will probably be reasonable for the first pull request to be able to do only
'D + T --> He-4 + n'since getting that will be tricky enough!There is a huge number of reactions that happen, so it will probably be worth working with applications in this order:
Reactions that are important in the laboratory include:
'D + T --> He-4 + n''D + D --> T + p''D + D --> T + p''D + He-3 --> He-4 + p''T + T --> He-4 + 2n''He-3 + He-3 --> He-4 + 2p''He-3 + T --> He-4 + p + n''He-3 + T --> He-3 + D''D + Li-6 --> 2He-4''D + Li-6 --> He-3 + He-4 + n''D + Li-6 --> Li-7 + p''D + Li-6 --> Be-7 + n''He-3 + Li-6 --> He-4 + He-3''p + B-11 --> 3He-4'I got these from the Wikipedia article on nuclear fusion, it turns out. The articles on the CNO cycle and proton-proton chain should also be helpful. There's even more in the article on stellar nucleosynthesis.
We'll probably want to first assume cases with a Maxwellian particle distribution and have the corresponding input be a temperature. Things will get more complicated if we have a neutron/particle beam so that we can't assume a particular temperature...but we can probably ignore that for now.
There's a chance that we could learn from or use other packages (like PyNE, though I think that may be more on the nuclear engineering side of things which is kind of different).
This will need documentation and tests as well. Reactions that aren't implemented should raise a
NotImplementedError.My sincere thanks!
-Nick
The text was updated successfully, but these errors were encountered: