Adjusting combustion component leads to no progress while solving

[HaSchneider]

Hi!

I am trying to include NO generation in the combustion component in my tespy fork but have some solving issues.

When I test it with:

from tespy.components import (
    Source, Sink, 
    SimpleHeatExchanger, 
    DiabaticCombustionChamber,
)
from tespy.tools.fluid_properties.wrappers import PyromatWrapper

from tespy.connections import Connection

from rva_classes import ModelTemplate

class my_PyromatWrapper(PyromatWrapper):
    def _set_constants(self):
        self._p_min, self._p_max = 100, 1000e5
        self._T_crit, self._p_crit = 180.15, 6400000 #self.AS.critical() not supported by pyromat
        self._T_min, self._T_max = self.AS.Tlim()
        self._molar_mass = self.AS.mw()



class create_solvent_combustion(ModelTemplate):
        
    def _create_network(self) -> None:
        super()._create_network()
        amb = Source('ambient air')
        sf = Source('fuel')
        fg = Sink('flue gas outlet')
        self.comb = DiabaticCombustionChamber('combustion chamber')
        self.hex = SimpleHeatExchanger('heat exchanger')
        
        
        self.amb_comb = Connection(amb, 'out1', self.comb, 'in1')
        self.sf_comb = Connection(sf, 'out1', self.comb, 'in2')
        self.comb_hex = Connection(self.comb, 'out1', self.hex, 'in1')
        self.hex_fg = Connection(self.hex, 'out1', fg, 'in1')
        self.nw.add_conns(self.sf_comb, self.amb_comb, self.comb_hex, self.hex_fg, 
        )
        self.comb.set_attr(ti=500000, pr=0.95, eta=1, lamb=4, f_nox= 0.00025)
        self.amb_comb.set_attr(
            p=1.2, T=20,
            fluid={'Ar': 0.0129, 'N2': 0.7553, 'CO2': 0.0004, 'O2': 0.2314, 'ig::NO':0},
            fluid_engines={"NO": my_PyromatWrapper}, 
        )
        self.sf_comb.set_attr(
            p=1.3, T=25, fluid={'CO2': 0.03, 'H2': 0.01, 'CH4': 0.96, 'ig::NO':0},
            fluid_engines={"NO": my_PyromatWrapper}, 
        )
        
        
        self.hex.set_attr(pr=0.95)

        #self.comb_hex.set_attr(T=1200)
        self.hex_fg.set_attr(T=200)

rva= create_solvent_combustion()

rva.nw.solve('design')

I get the following output:

 iter  | residual   | progress   | massflow   | pressure   | enthalpy   | fluid      | component  
-------+------------+------------+------------+------------+------------+------------+------------
 1     | 5.30e+07   | 0 %        | 1.41e+02   | 0.00e+00   | 1.15e+07   | 7.98e+01   | 0.00e+00   
 2     | 6.66e+08   | 0 %        | 1.42e+02   | 0.00e+00   | 1.47e+06   | 5.34e-03   | 0.00e+00   
 3     | 4.06e+06   | 0 %        | 1.29e-09   | 0.00e+00   | 5.83e+06   | 7.67e-01   | 0.00e+00   
 4     | 1.82e+04   | 19 %       | 6.37e-10   | 0.00e+00   | 3.49e+04   | 7.10e-01   | 0.00e+00   
 5     | 1.29e+04   | 20 %       | 3.29e-14   | 0.00e+00   | 2.62e+04   | 6.96e-01   | 0.00e+00   
 6     | 4.37e+03   | 26 %       | 4.69e-16   | 0.00e+00   | 8.85e+03   | 6.94e-01   | 0.00e+00   
...
 13    | 6.93e+00   | 57 %       | 3.49e-17   | 0.00e+00   | 1.40e+01   | 6.93e-01   | 0.00e+00   
 14    | 4.85e+00   | 59 %       | 3.49e-17   | 0.00e+00   | 9.91e+00   | 6.93e-01   | 0.00e+00   
 15    | 5.44e-01   | 69 %       | 3.50e-17   | 0.00e+00   | 6.46e-01   | 6.93e-01   | 0.00e+00   
 16    | 4.83e-01   | 70 %       | 3.50e-17   | 0.00e+00   | 1.09e-07   | 6.93e-01   | 0.00e+00   
 17    | 4.83e-01   | 70 %       | 3.50e-17   | 0.00e+00   | 1.37e-10   | 6.93e-01   | 0.00e+00   
 18    | 4.83e-01   | 70 %       | 3.50e-17   | 0.00e+00   | 1.37e-10   | 6.93e-01   | 0.00e+00   
 19    | 4.83e-01   | 70 %       | 3.50e-17   | 0.00e+00   | 1.37e-10   | 6.93e-01   | 0.00e+00   
The solver does not seem to make any progress, aborting calculation. Residual value is 4.83e-01. This frequently happens, if the solver pushes the fluid properties out of their feasible range.
 20    | 4.83e-01   | 70 %       | 3.50e-17   | 0.00e+00   | 1.37e-10   | 6.93e-01   | 0.00e+00   
...
 42    | 4.83e-01   | 70 %       | 3.50e-17   | 0.00e+00   | 1.37e-10   | 6.93e-01   | 0.00e+00   
Total iterations: 42, Calculation time: 0.37 s, Iterations per second: 112.82

I needed to adjust some other things to make use of pyromat fluid engine for NO.
It doesn't seem to calculate the NO content in the output properly. But I included it in the stoichiometry function.

What am I am missing here?

[fwitte]

Hi @HaSchneider,

great addition, thank you! I think there is two issues in the implementation right now:

First, you do need to subtract the mass of pure nitrogen vanishing from the mass balance due to conversion to NO. It seems (to me) in the first elif block that you were calculating the difference between what was available and what is converted. I think this would be the correct formulation (but please confirm):

            if self.lamb.val_SI < 1:
                dm=0
            elif n_nitrogen >= inl[1].m.val_SI*self.f_nox.val_SI /0.03001 *0.5:  #molar mass of NO: 0.03001 kg/mol
                dm = -inl[1].m.val_SI* self.f_nox.val_SI /0.03001 *inl[0].fluid.wrapper[self.n2]._molar_mass *0.5 #molar mass of NO: 0.03001 kg/mol
            else:
                dm= -n_nitrogen * inl[0].fluid.wrapper[self.n2]._molar_mass

Second, you also have to add a dm for the NO created in the process. That is missing:

        elif fluid == "NO":
            # nitrogen
            n_nitrogen = 0
            for i in inl:
                n_nitrogen += (
                    i.m.val_SI
                    * i.fluid.val.get(self.n2, 0)
                    / inl[0].fluid.wrapper[self.n2]._molar_mass
                )
            if self.lamb.val_SI < 1:
                dn=0
            elif n_nitrogen >= inl[1].m.val_SI*self.f_nox.val_SI /0.03001 *0.5:  #molar mass of NO: 0.03001 kg/mol
                dn = inl[1].m.val_SI* self.f_nox.val_SI /0.03001#molar mass of NO: 0.03001 kg/mol
            else:
                dn = n_nitrogen

            dm = dn * (inl[0].fluid.wrapper[self.n2]._molar_mass + inl[0].fluid.wrapper[self.o2]._molar_mass) / 2

Can you test if this gives the expected results?

Best

Francesco

[HaSchneider]

Yes, both things where right! And I found another wrong calculation of the nitrogen residuals.
It works now!