Detailed Description

Method of Moments with Interpolative Closure (MOMIC) model.

Definition at line 15 of file sootModel_MOMIC.h.

#include <sootModel_MOMIC.h>

Public Member Functions
virtual void	setSourceTerms (state &state)

virtual double	pahSootCollisionRatePerDimer (const state &state, const double mDimer)

	sootModel_MOMIC (size_t nsoot_, nucleationModel nucl_, growthModel grow_, oxidationModel oxid_, coagulationModel coag_)

	sootModel_MOMIC (size_t nsoot_, nucleationMech Nmech, growthMech Gmech, oxidationMech Omech, coagulationMech Cmech)

virtual	~sootModel_MOMIC ()

Public Member Functions inherited from soot::sootModel
virtual void	setSourceTerms (state &state)=0

void	checkSpec ()

virtual double	pahSootCollisionRatePerDimer (const double mDimer) const

virtual double	get_M0_sectional (const state &state)

virtual double	get_M1_sectional (const state &state)

	sootModel (size_t nsoot_, nucleationModel nucl_, growthModel grow_, oxidationModel oxid_, coagulationModel coag_)

	sootModel (size_t nsoot_, nucleationMech Nmech, growthMech Gmech, oxidationMech Omech, coagulationMech Cmech)

virtual	~sootModel ()

Private Member Functions
void	downselectIfNeeded (std::vector< double > &M)

double	f_grid (int x, int y)

double	g_grid (int y)

std::vector< double >	MOMICCoagulationRates (const state &state, std::vector< double > &M)

double	Mr (const double r)

void	set_diffTable (const std::vector< double > &M)

void	set_fractional_moments_Mp6_Mq6 ()

void	set_mDimerPowers ()

Private Attributes
size_t	Nmom
	number of soot moments (may decrese from downselectIfNeeded) More...

std::vector< double >	Mp6
	arrays holding fractional moments M_{p/6} where p is posiitive (Continuum) More...

std::vector< double >	Mq6
	arrays holding fractional moments M_{q/6} where q is negative (FM) More...

std::vector< size_t >	np
	# of Mp6 entries needed based on Nmom; if Nmom=3, then np[3] is # Mp6 entries to fill More...

std::vector< size_t >	nq
	# of Mq6 entries needed based on Nmom; if Nmom=3, then nq[3] is # Mq6 entries to fill More...

std::vector< std::vector< double > >	diffTable
	set in set_diffTable, used in Mr More...

double	mDn46
	mDimer^{-4/6}, etc.; for continuum More...

double	mDn26

double	mD26

double	mDn36
	mDimer^{-3/6}, etc.; for free molecular More...

double	mDn16

double	mD16

double	mD36

double	mD56

double	mD76

double	mD96

double	mD116

double	mD136

double	mD156

double	mD176

double	mD196

Additional Inherited Members
Public Attributes inherited from soot::sootModel
size_t	nsoot
	# of soot variables: moments or sections More...

nucleationModel *	nucl
	pointer to nucleation mechanism More...

growthModel *	grow
	pointer to growth mechanism More...

oxidationModel *	oxid
	pointer to oxidation mechanism More...

coagulationModel *	coag
	pointer to coagulation mechanism More...

bool	mechsNewedHere
	flag to delete "new" objects More...

psdMech	psdMechType
	one of MONO, LOGN, QMOM, MOMIC, SECT, etc. More...

std::vector< double >	mBins
	mass in sections for the sectional model More...

sourceTerms	sources
	struct containing soot, gas, and pah source terms vectors More...

Constructor & Destructor Documentation

◆ sootModel_MOMIC() [1/2]

sootModel_MOMIC::sootModel_MOMIC	(	size_t	nsoot_,
		nucleationModel *	nucl_,
		growthModel *	grow_,
		oxidationModel *	oxid_,
		coagulationModel *	coag_
	)

Constructor taking pointers to chemistry models as input. User creates these pointers nominally by "new-ing" them.

Parameters

nsoot_	input: number of soot moments (3-8).
nucl_	input: pointer to nucleation model.
grow_	input: pointer to growth model.
oxid_	input: pointer to oxidation model.
coag_	input: pointer to coagulation model.

Definition at line 21 of file sootModel_MOMIC.cc.

◆ sootModel_MOMIC() [2/2]

sootModel_MOMIC::sootModel_MOMIC	(	size_t	nsoot_,
		nucleationMech	Nmech,
		growthMech	Gmech,
		oxidationMech	Omech,
		coagulationMech	Cmech
	)

Constructor taking enumerations names as input. Chemistry pointers are created (new-ed) here based on those enumerations.

Parameters

nsoot_	input: number of soot moments (3-8).
Nmech	input: one of enum class nucleationMech in sootDefs.h
Gmech	input: one of enum class growthMech in sootDefs.h
Omech	input: one of enum class oxidationMech in sootDefs.h
Cmech	input: one of enum class coagulationMech in sootDefs.h

Definition at line 64 of file sootModel_MOMIC.cc.

◆ ~sootModel_MOMIC()

virtual soot::sootModel_MOMIC::~sootModel_MOMIC ( )

inlinevirtual

Definition at line 70 of file sootModel_MOMIC.h.

Member Function Documentation

◆ setSourceTerms()

void sootModel_MOMIC::setSourceTerms ( state & state )

virtual

Primary user interface.

Parameters

state input: gas and soot state, set by user.

sets sources.sootSources vector sets sources.gasSources vector sets sources.pahSources vector

Assumes mDn36, etc. (mDimer^powers) have been set in set_mDimerPowers().

Implements soot::sootModel.

Definition at line 107 of file sootModel_MOMIC.cc.

◆ pahSootCollisionRatePerDimer()

virtual double soot::sootModel_MOMIC::pahSootCollisionRatePerDimer	(	const state &	state,
		const double	mDimer
	)

virtual

◆ downselectIfNeeded()

void sootModel_MOMIC::downselectIfNeeded ( std::vector< double > & M )

private

Reduces the number of moments if needed to avoid invalid inversion.

Parameters

M	input: vector of moment values

Definition at line 238 of file sootModel_MOMIC.cc.

◆ f_grid()

double sootModel_MOMIC::f_grid	(	int	x,
		int	y
	)

private

Calculates the grid function described in Frenklach 2002 MOMIC paper using polynomial interpolation between whole order moments.

Parameters

x	input: x grid point
y	input: y grid point

Returns: f_{1/2}

/// Mq6: (-3 -1 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51 53 55
/// indx:  0  1 2 3 4 5 6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
/// f_l^(x,y) is composed of M_(x+7/6)*M_(y+3/6), etc. where fractions vary
///      All M_(frac) are in Mq6. All frac are odd sixth fractions. 
///      If x = 1, then, e.g., x+7/6 --> (6+7)/6=13/6; index 7/6 is 5, shift by 3 --> 13/6 is ind 8
///      If x = 2, then shift by 6; x = 3, then shift by 9, etc.
///

Definition at line 294 of file sootModel_MOMIC.cc.

◆ g_grid()

double soot::sootModel_MOMIC::g_grid ( int y )

private

◆ MOMICCoagulationRates()

std::vector< double > soot::sootModel_MOMIC::MOMICCoagulationRates	(	const state &	state,
		std::vector< double > &	M
	)

private

◆ Mr()

double sootModel_MOMIC::Mr ( const double r )

private

Compute fractional moment r by Newton forward polynomial.

/// x points are 0,   1,  2,  3, ... (hard-coded this way)
/// y points are M0, M1, M2, M3, ... (log10 of those, or what is set in set_diffTable)
///     but returned value is 10^l10Mr here.
///

Assumes diffTable is set. Mr = M0 + r*dM0 + r(r-1)*ddM0/2!+ r(r-1)(r-2)*dddM0/3! + ... Positive and negative moments get different treatement. Positive interpolates among all moments. Negative extrapolates from M0, M1, M2 (log of those). Code verified by comparison to np.polyfit, np.polyval.

Parameters

r	input: fractional moment r

Returns: fractional moment Mr

Definition at line 645 of file sootModel_MOMIC.cc.

◆ set_diffTable()

void sootModel_MOMIC::set_diffTable ( const std::vector< double > & M )

private

Calculates the grid function described in Frenklach 2002 MOMIC paper using polynomial interpolation between whole order moments

Parameters

y	input: y grid point

Returns: g_grid value.

Assumes mDn36, etc. (mDimer^powers) have been set in set_mDimerPowers()

/// Mq6: (-3 -1 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51 53 55
/// indx:  0  1 2 3 4 5 6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
/// f_l^(x,y) is composed of M_(x+7/6)*M_(y+3/6), etc. where fractions vary
///      All M_(frac) are in Mq6. All frac are odd sixth fractions. 
///      If x = 1, then, e.g., x+7/6 --> (6+7)/6=13/6; index 7/6 is 5, shift by 3 --> 13/6 is ind 8
///      If x = 2, then shift by 6; x = 3, then shift by 9, etc.
/// \end verbatim
///
////////////////////////////////////////////////////////////////////////////////

double sootModel_MOMIC::g_grid(int y) {

    int yi = (y-1)*3;       // add this index

    double g0, g1, g2, g3;

    //----------- 

    g0 =            mD16   * Mq6[yi+0] +
                 2.*mDn16  * Mq6[yi+1] + 
                    mDn36  * Mq6[yi+2];

    g1 =            mD76   * Mq6[yi+0] +  
                 2.*mD56   * Mq6[yi+1] +
                    mD36   * Mq6[yi+2] +
                    mD16   * Mq6[yi+3] +
                 2.*mDn16  * Mq6[yi+4] +
                    mDn36  * Mq6[yi+5] ;

    if (y<=3)    
        g2 =        mD136  * Mq6[yi+0] +  
                 2.*mD116  * Mq6[yi+1] +
                    mD96   * Mq6[yi+2] +
                 2.*mD76   * Mq6[yi+3] +
                 4.*mD56   * Mq6[yi+4] +
                 2.*mD36   * Mq6[yi+5] +
                    mD16   * Mq6[yi+6] +
                 2.*mDn16  * Mq6[yi+7] +
                    mDn36  * Mq6[yi+8] ;

    if (y<=2)    
        g3 =        mD196 * Mq6[yi+0] +  
                 2.*mD176 * Mq6[yi+1] +
                    mD156 * Mq6[yi+2]  +
                 3.*mD136 * Mq6[yi+3]  +
                 6.*mD116 * Mq6[yi+4]  +
                 3.*mD96  * Mq6[yi+5]  +
                 3.*mD76  * Mq6[yi+6]  +
                 6.*mD56  * Mq6[yi+7]  +
                 3.*mD36  * Mq6[yi+8]  +
                    mD16  * Mq6[yi+9]  +
                 2.*mDn16 * Mq6[yi+10]  +
                    mDn36 * Mq6[yi+11]  ;

    //----------- Lagrange interpolation to point 1/2 using some or all of points 0, 1, 2, 3

    double ghalf;
    if (y<=2)       // 4 point interpolant
        ghalf = pow(g0,5./16.) * pow(g1,15./16.) * pow(g2, -5./16.) * pow(g3, 1./16.);
    else if (y<=3)  // 3 point interpolant
        ghalf = pow(g0, 3./8.) * pow(g1,  3./4.) * pow(g2,  -1./8.);
    else            // h point interpolant
        ghalf = pow(g0, 1./2.) * pow(g1,  1./2.);

    return ghalf;
}

////////////////////////////////////////////////////////////////////////////////
///
/// Get coagulation rates for moment k
/// 
/// @param state \input        gas and soot state, set by user.
/// @param M     \input        vector of soot moments
/// @return vector of coagulation rates
///
////////////////////////////////////////////////////////////////////////////////

vector<double> sootModel_MOMIC::MOMICCoagulationRates(const state& state, vector<double>& M){

    if (M[0] <= 0.0) return vector<double>(Nmom, 0.0);

    vector<double> Rates_C(Nmom, 0.0);
    vector<double> Rates_FM(Nmom, 0.0);
    vector<double> Rates(Nmom, 0.0);

    //----------- free-molecular regime

    if (coag->mechType == coagulationMech::FM || 
        coag->mechType == coagulationMech::HM) {

        const double Kfm = coag->getKfm(state);

        for (size_t r=0; r<Nmom; r++) {
            if (r==1) continue; 
            if (r==0)
                Rates_FM[r] = -0.5*Kfm*f_grid(0,0);
            else {
                for (size_t k=1; k<=r-1; k++)
                    Rates_FM[r] += binomial_coefficient(r,k) * f_grid(k, r-k);
                Rates_FM[r] *= 0.5*Kfm;
            }
        }
    }

    //----------- continuum regime

    if (coag->mechType == coagulationMech::CONTINUUM || 
             coag->mechType == coagulationMech::HM) {

        const double Kc  = coag->getKc( state);
        const double Kcp = coag->getKcp(state);


        for (size_t r=0; r<Nmom; r++) {
            if (r==1) continue; 
            if (r==0)
                Rates_C[r] = -Kc*( Mp6[2]*Mp6[2] + Mp6[1]*Mp6[3] +            // M_0*M_0 + M_{-2/6}*M_{2/6}
                                   Kcp*( Mp6[1]*Mp6[2] + Mp6[0]*Mp6[3] ) );   // M_{-2/6}*M_0 + M_{-4/6}*M_{2/6}
            else {
                size_t kk, rk;        // index shifters
                for (size_t k=1; k<=r-1; k++) {
                    kk = k*3;
                    rk = (r-k)*3;
                    Rates_C[r] = binomial_coefficient(r,k) * (
                                       Mp6[kk+1]*Mp6[rk+3] +                  // M_{k-2/6}*M_{r-k+2/6}
                                 2.0*  Mp6[kk+2]*Mp6[rk+2] +                  // M_{k+0/6}*M_{r-k+0/6}
                                       Mp6[kk+3]*Mp6[rk+1] +                  // M_{k+2/6}*M_{r-k-2/6}
                                 Kcp*( Mp6[kk+0]*Mp6[rk+3] +                  // M_{k-4/6}*M_{r-k+2/6}
                                       Mp6[kk+1]*Mp6[rk+2] +                  // M_{k-2/6}*M_{r-k+0/6}
                                       Mp6[kk+2]*Mp6[rk+1] +                  // M_{k+0/6}*M_{r-k-2/6}
                                       Mp6[kk+3]*Mp6[rk+0] ) );               // M_{k+2/6}*M_{r-k-4/6}
                }
                Rates_C[r] *= 0.5*Kc;
            }
        }
    }

    //----------- finalize and return

    for (size_t r=0; r<Nmom; r++) {
        if (r==1) continue;
        if (coag->mechType == coagulationMech::FM)
            Rates[r] = Rates_FM[r];
        else if (coag->mechType == coagulationMech::CONTINUUM)
            Rates[r] = Rates_C[r];
        else      // harmonic mean
            Rates[r] = Rates_FM[r]*Rates_C[r] / (Rates_FM[r]+Rates_C[r]);
    }

    return Rates;
}

////////////////////////////////////////////////////////////////////////////////
///
/// Compute PAH condensation terms for MOMIC model.
/// Function split out from setSourceTerms so that it can be called in nucleationModel_PAH
/// for computing the pah dimer concentration.
///
/// Function only called if nucleationMech::PAH.
/// Function called by nucleationModel_PAH::getNucleationSootRate
///
/// @param state  \input gas and soot state, set by user.
/// @param mDimer \input dimer mass (kg)
/// @return pah/soot sollision rate per dimer. Call it I. I*mDimer*nDimer = Cnd1 (=) kg/m3*s
///
////////////////////////////////////////////////////////////////////////////////

double sootModel_MOMIC::pahSootCollisionRatePerDimer(const state &state, const double mDimer) {

    if (nucl->mechType != nucleationMech::PAH)
        return 0.0;

    set_mDimerPowers();

    //------ continuum

    const double Kc  = coag->getKc( state);
    const double Kcp = coag->getKcp(state);

    double Ic1 = Kc * (  
                  mD26*Mp6[1] +            //  mD26*M_{-2/6}
                    2.*Mp6[2] +            //    2.*M_{ 0/6}
                 mDn26*Mp6[3] +            // mDn26*M_{ 2/6}
                 Kcp*( mD26*Mp6[0] +       //  mD26*M_{-4/6}
                            Mp6[1] +       //       M_{-2/6}
                      mDn26*Mp6[2] +       // mDn26*M_{ 0/6}
                      mDn46*Mp6[3] ) );    // mDn46*M_{ 2/6}

    //------ free molecular

    const double Kfm = coag->getKfm(state);

    double Ifm1 = Kfm*g_grid(1);

    //------ harmonic mean

    return Ifm1*Ic1/(Ifm1 + Ic1);
}

///////////////////////////////////////////////////////////////////////////////
/// 
/// Compute the difference table for a Newton polynomial.
/// \verbatim
/// x points are 0,   1,  2,  3, ... (hard-coded this way).
/// y points are M0, M1, M2, M3, ... (or log10 of those, or whatever is passed to l10M)
/// M0 dM0=M1-M0 ddM0=dM1-dM0 dddM0=ddM1-ddM0
/// M1 dM1=M2-M1 ddM1=dM2-dM1
/// M2 dM2=M3-M2           
/// M3
///

(Interpolating among log10(M) here; nomenclature in code reflects that but there is no change other than the var names.)

Function used by get_Mr. Function only needs to be set when the set of moments change, so call once at the top of the rates routine. Code verified by comparison to np.polyfit, np.polyval.

Parameters

M	input: vector of Moments, converted to log10(Moments).

Definition at line 611 of file sootModel_MOMIC.cc.

◆ set_fractional_moments_Mp6_Mq6()

void sootModel_MOMIC::set_fractional_moments_Mp6_Mq6 ( )

private

Compute arrays of fractional moments. Used to compute coagulation.

/// Mp6 = M_{p/6}          evens (Continuum)
///     Mp6 = (-4 -2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38) / 6
///     indx:   0  1 2 3 4 5 6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21
///     indx = (p+4)/2
/// Mq6 = M_{q/6}          odds (FM)
///     Mq6 = (-3 -1 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51 53 55) / 6
///     indx:   0  1 2 3 4 5 6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
///     indx = (q+3)/2
///

Definition at line 676 of file sootModel_MOMIC.cc.

◆ set_mDimerPowers()

void sootModel_MOMIC::set_mDimerPowers ( )

private

Set powers of mDimer used in other routines.

Definition at line 694 of file sootModel_MOMIC.cc.

Member Data Documentation

◆ Nmom

size_t soot::sootModel_MOMIC::Nmom

private

number of soot moments (may decrese from downselectIfNeeded)

Definition at line 20 of file sootModel_MOMIC.h.

◆ Mp6

std::vector<double> soot::sootModel_MOMIC::Mp6

private

arrays holding fractional moments M_{p/6} where p is posiitive (Continuum)

Definition at line 22 of file sootModel_MOMIC.h.

◆ Mq6

std::vector<double> soot::sootModel_MOMIC::Mq6

private

arrays holding fractional moments M_{q/6} where q is negative (FM)

Definition at line 23 of file sootModel_MOMIC.h.

◆ np

std::vector<size_t> soot::sootModel_MOMIC::np

private

# of Mp6 entries needed based on Nmom; if Nmom=3, then np[3] is # Mp6 entries to fill

Definition at line 24 of file sootModel_MOMIC.h.

◆ nq

std::vector<size_t> soot::sootModel_MOMIC::nq

private

# of Mq6 entries needed based on Nmom; if Nmom=3, then nq[3] is # Mq6 entries to fill

Definition at line 25 of file sootModel_MOMIC.h.

◆ diffTable

std::vector<std::vector<double> > soot::sootModel_MOMIC::diffTable

private

set in set_diffTable, used in Mr

Definition at line 27 of file sootModel_MOMIC.h.

◆ mDn46

double soot::sootModel_MOMIC::mDn46

private

mDimer^{-4/6}, etc.; for continuum

Definition at line 29 of file sootModel_MOMIC.h.

◆ mDn26

double soot::sootModel_MOMIC::mDn26

private

Definition at line 30 of file sootModel_MOMIC.h.

◆ mD26

double soot::sootModel_MOMIC::mD26

private

Definition at line 30 of file sootModel_MOMIC.h.

◆ mDn36

double soot::sootModel_MOMIC::mDn36

private

mDimer^{-3/6}, etc.; for free molecular

Definition at line 31 of file sootModel_MOMIC.h.

◆ mDn16

double soot::sootModel_MOMIC::mDn16

private

Definition at line 32 of file sootModel_MOMIC.h.

◆ mD16

double soot::sootModel_MOMIC::mD16

private

Definition at line 32 of file sootModel_MOMIC.h.

◆ mD36

double soot::sootModel_MOMIC::mD36

private

Definition at line 32 of file sootModel_MOMIC.h.

◆ mD56

double soot::sootModel_MOMIC::mD56

private

Definition at line 32 of file sootModel_MOMIC.h.

◆ mD76

double soot::sootModel_MOMIC::mD76

private

Definition at line 32 of file sootModel_MOMIC.h.

◆ mD96

double soot::sootModel_MOMIC::mD96

private

Definition at line 33 of file sootModel_MOMIC.h.

◆ mD116

double soot::sootModel_MOMIC::mD116

private

Definition at line 33 of file sootModel_MOMIC.h.

◆ mD136

double soot::sootModel_MOMIC::mD136

private

Definition at line 33 of file sootModel_MOMIC.h.

◆ mD156

double soot::sootModel_MOMIC::mD156

private

Definition at line 33 of file sootModel_MOMIC.h.

◆ mD176

double soot::sootModel_MOMIC::mD176

private

Definition at line 33 of file sootModel_MOMIC.h.

◆ mD196

double soot::sootModel_MOMIC::mD196

private

Definition at line 33 of file sootModel_MOMIC.h.

The documentation for this class was generated from the following files:

/Users/dol4/Desktop/sootlib/src/sootModel_MOMIC.h
/Users/dol4/Desktop/sootlib/src/sootModel_MOMIC.cc

Detailed Description

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Constructor & Destructor Documentation

◆ sootModel_MOMIC() [1/2]

◆ sootModel_MOMIC() [2/2]

◆ ~sootModel_MOMIC()

Member Function Documentation

◆ setSourceTerms()

◆ pahSootCollisionRatePerDimer()

◆ downselectIfNeeded()

◆ f_grid()

◆ g_grid()

◆ MOMICCoagulationRates()

◆ Mr()

◆ set_diffTable()

◆ set_fractional_moments_Mp6_Mq6()

◆ set_mDimerPowers()

Member Data Documentation

◆ Nmom

◆ Mp6

◆ Mq6

◆ np

◆ nq

◆ diffTable

◆ mDn46

◆ mDn26

◆ mD26

◆ mDn36

◆ mDn16

◆ mD16

◆ mD36

◆ mD56

◆ mD76

◆ mD96

◆ mD116

◆ mD136

◆ mD156

◆ mD176

◆ mD196