speth / ember

Ember: unsteady strained flame solver
MIT License
47 stars 25 forks source link

jSoret = 0 for all transport models #16

Closed marina8888 closed 2 years ago

marina8888 commented 2 years ago

Taking a gas susceptible to thermal diffusion (H2), the following code generates a thermal flux profile jSoret = 0 for both mixture averaged and multicomponent formulation.

Solution: If this is an error (rather than something for documentation), i think cantera has soret_enabled = False as default. Would adding this feature to the chemistry0D.cpp, void CanteraGas::setOptions function help?

To recreate: Cantera = 2.6.0 Ember = dev (main) branch

def main():
    conf = Config(
        Paths(outputDir='src/solutions'),
        InitialCondition(fuel='CH4:1.0',
                         equivalenceRatio=0.70),
        StrainParameters(initial=400,
                         final=400),
        Chemistry(mechanismFile="src/mech/gri-3.0.cti", transportModel = 'Approx'),
        TerminationCondition(tEnd=0.010))
    conf.run()`

    `filename = 'src/solutions/prof000007.h5'
    s = utils.load(filename)
    gas_obj = ct.Solution("src/mech/gri-3.0.cti")

    print(f' testing species {gas_obj.species_names[0]}')
    utils.expandProfile(s, gas_obj, diffusion=True, reaction_rates=True)

    plt.title('Fickian Diffusion Flux')
    plt.plot(s.x, s.jFick[0,:])
    plt.show()

    plt.title('Soret Diffusion Flux')
    plt.plot(s.x, s.jSoret[0,:])
    plt.show()

if __name__ == '__main__':
    main()
marina8888 commented 2 years ago

Just to follow up on this with some more information. I believe various transport models are not being activated (for example the code below will create the same mixture fraction vs grid profile for all transport models (although the profile should be flat for UnityLewis).

def main():
    conf = Config(
        Paths(outputDir='src/solutions'),
        InitialCondition(fuel='CH4:1.0',
                         equivalenceRatio=0.70),
        StrainParameters(initial=400,
                         final=400),
        Chemistry(mechanismFile="src/mech/gri-3.0.cti", transportModel = 'UnityLewis'),
        TerminationCondition(tEnd=0.010))
    conf.run()

    filename = 'src/solutions/prof000007.h5'
    s = utils.load(filename)
    gas_obj = ct.Solution("src/mech/gri-3.0.cti")

    print(f' testing species {gas_obj.species_names[0]}')
    utils.expandProfile(s, gas_obj, diffusion=True, reaction_rates=True)

    Z = np.zeros(s.T.shape)

    # fill mixture fraction:
    for n in range(len(s.x)):
        print(n)
        gas_obj.TPY = s.T[n], s.P, s.Y[:, n]
        Z[n] = gas_obj.mixture_fraction(fuel = 'CH4:1.0', oxidizer='O2:0.21, N2:0.79')

    plt.plot(s.x, Z)
    plt.show()

if __name__ == '__main__':
    main()
speth commented 2 years ago

I assume that in your first example, the argument should be transportModel = 'Multi' rather than transportModel = 'Approx', correct? In this case, I believe the actual calculations of the flame are being done correctly.

The problem is that the gas_obj being passed into utils.expandProfile doesn't have the same transport model as the one used to solve the flame -- it just has the default for the gri30.cti input file, which is the mixture-averaged transport model. If you correct this by adding the line

    gas_obj.transport_model = 'multicomponent'

just after creating gas_obj, I think you'll see the expected result.

In the second case, while the resulting plot looks like there is a non-flat profile for the unity-Le transport model, the default axes used by Matplotlib are a bit misleading here, and the actual variation in the profile is only on the order of the 5th significant digit.

marina8888 commented 2 years ago

@speth Thank you for your response gas_obj.transport_model = 'multicomponent' Did the trick, thank you very much!