Add Humid Air Capability

[jdgratz10]

Added files to allow for humid air at flowstart, and made a slight modification to the Combustor element to accomodate the changes.

[askprash]

Hi @jbergeson I was trying to work with humidity as well, and ran into convergence issues in chem_eq . I see you've defined a new thermo set here, does that fix it?

Interestingly I found the that I run into issues only if my element definitionHUMID_AIR_MIX = {'N2': 78.2, 'O2': 20.78, 'H2O':1.0, 'CO2':.01, 'Ar':0.01} has both H and C.

I was able to fix this by changing the residual weighting in chem_eq.py as follows:

def _resid_weighting(n):
        np.seterr(under='ignore')
        return (1 / (1 + np.exp(-1e4 * n)) - .5) * 2

instead of,

def _resid_weighting(n):
        np.seterr(under='ignore')
        return (1 / (1 + np.exp(-1e5 * n)) - .5) * 2

Note change form -1e5 to -1e4 However the convergence of the default AIR_MIX with janaf now takes roughly twice as long to converge :/

[askprash]

Just saw your note in wet_simple_turbojet.py! Any idea why this happens, is this due to missing species that help drive the chemical potential down? Seems like the Janaf definition in NPSS has CH3OH as well, I thought maybe adding that might help the chemical convergence but alas, no luck.

[eshendricks]

@askprash You found this PR as we (mostly @jbergeson ) is working on trying to address your issue with humidity. We are still reviewing a few changes but we think we have figured out the answer. I was planning to message you once we had this fully sorted out and merged in.

The challenge is really that when water (specifically hydrogen) is added into the mix, by default the janaf thermo file include all the products of combustion. To achieve chemical equilibrium upstream of the combustor, many of these products end up being trace species with very low concentrations. These trace species present convergence challenges for the chemical equilibrium solver. You were on the right track by changing the scalar in the _resid_weighting above, but that doesn't entirely fix the problem. The best solution seems to be to completely remove those trace species from the analysis upstream of the combustor resulting in a more stable solution.

[askprash]

@eshendricks and @jbergeson ah thank you so much for looking into this! Sorry I didn't mean to jump the gun on this just got excited when I saw this PR haha! 😄

[jdgratz10]

@askprash thank you for your patience! I added an extra component that allows you to specify a water to air ratio (WAR) option in the flight conditions element. I added the implementation of this in the wet_propulsor.py and wet_simple_turbojet.py examples so you can take a look. If you'd like to go ahead and give it a try feel free, and let us know if you run into any issues with it!

[JustinSGray]

@askprash we're going to leave this PR open for a few days, to give you a chance to test it out. @jbergeson added some examples cases that show you how to use the WAR calculation. Let us know if this works for you.

The key fix was creating a reduced set of species that included H2O and OH but none of the hydrocarbon species. The details of why that is needed related to numerical challenges with trace species in the chemical equilibrium calculations. You can read more about those details in this paper if you like.

Note from the examples that were included how you have to update the thermo_data for all the elements before the combustor to handle the humid air case.

You might wonder why we needed to do this for humid air, but not for dry air (for which the examples use the full Janaf thermo-data). The answer lies in the details of how we process that data. When you specify the elements, pyCycle process all the species and removes any that require elements that are missing. In dry air, there is no hydrogen, so all species with any H in them get discarded. When you switch to wet air, obviously that trick doesn't work anymore, so you need to be a bit more selective as to what species you let the chem_eq consider

[askprash]

Perfect! I'll test it out! Thanks a ton @jbergeson, @eshendricks and @JustinSGray. Thanks for the explanation as well! It makes sense, as soon we have hydrogen and carbon, a whole swath of HC compounds need to be accounted for even though they will be driven to trace concentrations upstream of the combustor. One question I had was, why does simply setting remove_trace_species = True inchem_eq.py` not do the trick? From the paper you linked I'm guessing it is more that we do not want to use this so we have well defined gradients? I'll go through the paper in more detail! Thanks!

[JustinSGray]

the remove_trace_species stuff is a numerical trick that I developed. Its very add_hoc, but the code itself flips the flag as the last step of its residual evaluation:

        if np.linalg.norm(resids['n']) < 1e-4:
            self.remove_trace_species = True
        else:
            self.remove_trace_species = False

So you manually changing it won't have a permanent effect. The goal of this hack is to prevent trace species from limiting convergence. In a sense, once something has become a trace species you don't care about its residual any more. However if you're just going to drop it from the solution (which is what remove_trace_species effectively does) you have to take care that it doesn't become trace to early and then you get stuck.

So the residual damping helps with that. If the overall convergence is poor, then the residual damping helps to shrink the value of the residual for anything that is currently trace. However, it still provides a path for the solver to choose to increase the concentration back up later on.

Once the remove_trace_species has been triggered though, the way I modify the derivatives effectively blocks that species from ever coming back in. All the interaction partial derivatives go to 0, so the newton solver would never choose to add it back in. That is what this code does:

        # Replace J for tiny values of n with identity
        if self.remove_trace_species:
            n = outputs['n']
            for j in range(num_prod):
                if n[j] <= 1.0e-10:
                    dRdy[j, :] = 0.0
                    dRdy[j, j] = -1.0

This is a numerically smooth equivalent operation to what the original CEA code did: If something was trace toward the end of convergence, then it removed that species from the list and re-solves a new chem_eq problem with a small set of possible chemical species. In OpenMDAO, we need to keep the sizes of arrays constant, so I numerically remove it from the solve without actually resizing the arrays.

I just wait till I get below 1e-4 tolerance before activating that, to ensure sufficient time for the solution to find a good neighboorhood to be in. The 1e-4 is admittedly ad-hoc. It is what I found to work well in the cases I played with. It hasn't failed me yet though.