Auto Flamelet V2计算过程解析

最新推荐文章于 2024-06-29 16:57:48 发布

pud_ding

最新推荐文章于 2024-06-29 16:57:48 发布

阅读量165

点赞数

分类专栏： CFD相关理论文章标签： python 开发语言

本文链接：https://blog.csdn.net/pud_ding/article/details/129753122

版权

CFD相关理论专栏收录该内容

1 篇文章 0 订阅

订阅专栏

代码来自于Github，用于湍流非预混火焰燃烧模拟中的火焰面模型中层流火焰面数据库的生成，我们对其计算过程进行梳理。

1. 输入文件的准备

给定需计算火焰面的个数，化学计量数处的混合物分数，初始流率，反应机理，对撞火焰入口之间的距离，求解对撞火焰面的初始网格节点数等信息。

Max Number of Flamelets:
60
Mixture Fraction at stoich:
0.351
Pressure:
101325
Temperature of Fuel:
300
Temperature of Oxidizer:
300
Initail Flow Rate:
0.1
Composition of Fuel(Mole Fraction):
CH4:0.25, O2:0.1575, N2:0.5925
Composition of Oxidizer(Mole Fraction):
O2:0.21, N2:0.79
Mechnism:
chem.yaml
Library Name for Output:
flamelet.d
Molecular Weights File Name:
mole_weights.d
One dimension width:
0.15
Number of Grid Point:
100
Number of Zeta:
10
===For FPV model===
Definition of progress variable (C):
CO2, H2O


/CH4:0.25, O2:0.1575, N2:0.5925, 0.351
/CH4:1, 0.0552
/mole fraction: X; mass fraction: Y
/CH4:0.221, H2:0.332, N2:0.447, 0.167

将相关信息读入后就可以进行最大、最小质量流率探测，将其保存在critical_rates数组中。

if switch:
    cdef = input_data[14]
if massread:
    critical_rates = np.load(file='mass_range.npy', allow_pickle=True)
else:
    critical_rates = massprobe(press, tempf, tempo, mdot, compf, compo, mech, width, npoint, air_multi)
if masssave:
    np.save(file='mass_range.npy', arr=critical_rates)

完成质量流率探测后，然后将质量流率在对数空间将其分割成给定的火焰面数量。

minrate = critical_rates[0]
maxrate = critical_rates[1]
exrate = critical_rates[2]
massflow = np.logspace(np.log10(minrate), np.log10(maxrate), flamenum, endpoint=False)

Question：将熄灭的那个火焰质量流率也放入到了质量流率数组中，后面还对该质量流率进行计算？是否是用来计算最大的可燃烧质量流率？

massflow = np.append(massflow, exrate)

2. 火焰面求解

对每一个质量流率的火焰面进行求解，保存每个火焰面的温度、密度分布、最高温度、组分分布信息、混合物分数分布。

并获取每个火焰面的应变率（strain rate），通过应变率计算标量耗散率。

tempset = flame.T
tmax = max(tempset)
rhoset = flame.density

a = flame.strain_rate('potential_flow_oxidizer')    # get the strain rate
dissrate = (2*a/np.pi)*np.exp(-2.0*(ss.erfcinv(2.0*zst)**2))

mixfracset = flame.mixture_fraction('Bilger')   # the Bilger mixture fraction
mixfracset = posfilter(mixfracset)   # filter

for isp in range(nsp):
    dataset.append(posfilter(list(reversed(flame.Y[isp]))))
dataset.append(list(reversed(tempset)))
dataset.append(list(reversed(rhoset)))
dataset.append(list(reversed(mixfracset)))

应变率计算方法有以下五种：

平均（mean）：整个计算域的平均轴向速度梯度

$a_{mean}=\left | \frac{\Delta u}{\Delta z} \right |$

最大值（max）：计算域中最大的速度梯度

$a_{max}=max \left( \left | \frac{d u}{d z} \right | \right )$

化学计量数（stoichiometric）：化学计量数处的速度梯度

$a_{stoichiometric}=\left| \frac{d u}{ d z}\mid_{\phi=1} \right|$

氧化剂/燃料势流入口：入口给定势流边界条件对应的应变率

$a_f=\sqrt{-\frac{\Lambda}{\rho_f}}$

$a_o=\sqrt{-\frac{\Lambda}{\rho_o}}$

对混合物分数进行标准化，使其在[0,1]区间均匀变化。只保存了最大温度后面一系列的火焰面。

    normzs = np.linspace(0, 1, npoint)
    tmax_marker = tmaxs.index(max(tmaxs))
    dissrate_stable = dissrates[tmax_marker]
    print('---Stable branch finished---')

    dissrates = dissrates[tmax_marker:]
    dataseries = dataseries[tmax_marker:]
    tmaxs = tmaxs[tmax_marker:]

对保存的一系列火焰面进行处理，利用混合物分数对温度、密度、组分信息进行插值，通过得到插值函数计算标准化后的混合物分数节点上的温度、密度组分信息。

for i in range(len(dissrates)):
    tempfunc = interpolate.interp1d((dataseries[i][nsp + 2]), (dataseries[i][nsp]))
    dataseries[i][nsp] = tempfunc(normzs)
    rhofunc = interpolate.interp1d((dataseries[i][nsp + 2]), (dataseries[i][nsp + 1]))
    dataseries[i][nsp + 1] = rhofunc(normzs)
    if switch:
        profunc = interpolate.interp1d((dataseries[i][nsp + 2]), (dataseries[i][nsp + 4]))
        varfunc = interpolate.interp1d((dataseries[i][nsp + 2]), (dataseries[i][nsp + 3]))
        dataseries[i][nsp + 4] = profunc(normzs)
        dataseries[i][nsp + 3] = varfunc(normzs)
        csts[i] = (max(dataseries[i][nsp + 4]))
    for isp in range(nsp):
        massfunc = interpolate.interp1d((dataseries[i][nsp + 2]), (dataseries[i][isp]))
        dataseries[i][isp] = massfunc(normzs)
    dataseries[i][nsp + 2] = normzs

3.火焰面数据库积分

程序中采用多线程同时对多个火焰面进行积分运算，我们只关心详细的积分运算过程：

混合物分数的方差zetas分布采用线性分布[0,0.11,0.22,0.33,...,0.99]。分别对温度、密度、组分信息积分。

def presume(package, switch):
    datas = package[0]
    nsp = package[1]
    nzeta = package[3]
    metadata = []
    zetas = np.linspace(0, 0.99, nzeta, endpoint=True)
    for z in zetas:
        intmassset = []
        inttemps = pdfinter(z, datas[nsp], datas[nsp + 2])        # temperature integration 
        if switch:
            intrhos = pdfinter(z, datas[nsp + 4], datas[nsp + 2])
            intvar = pdfinter(z, datas[nsp + 3], datas[nsp + 2])
        else:
            intrhos = pdfinter(z, datas[nsp + 1], datas[nsp + 2])   # density integration
        for isp in range(nsp):
            intmass = pdfinter(z, datas[isp], datas[nsp + 2])     # every species integration   
            intmassset.append(intmass)
        if switch:
            metadata.append([intmassset, inttemps, intrhos, datas[nsp + 2], intvar])
        else:
            metadata.append([intmassset, inttemps, intrhos, datas[nsp + 2]]) # save integration terms
    return metadata

积分过程：

zeta为归一化后的混合物分数方差；zs为对应火焰面的混合分数分布，该值是标准化后的值，对于每个火焰面是相同的。

varzs = zeta*(zs*(1.0 - zs))计算归一化之前的混合分数方差。

def pdfinter(zeta, sourcelist, zs):
    varzs = np.zeros(len(sourcelist))
    pdfs = np.zeros(len(sourcelist))
    palphas = np.zeros(len(sourcelist))
    pbetas = np.zeros(len(sourcelist))
    points = 250
    hzs = np.zeros(points)
    helpzs = np.zeros(points)
    hzs = np.zeros(points)
    deltas = np.zeros(points)
    datalist = sourcelist.copy()
    varzs = zeta*(zs*(1.0 - zs))

根据每一个混合分数均值zs和对应的混合分数方差varzs的值，来计算对应的beta函数中概率密度函数所用到的palphas，pbetas。

然后利用得到alpha，beta的值来计算混合分数均值对应的概率密度函数的值。