我参考cantera官网示例中的diffusion_flame_batch.py，设置一批应变率变化时的一维对冲火焰解。通过反复调用前一计算结果作为初始解，来逐渐增加或降低应变率，从而得到计算结果。当我按照代码设置应变率上升时（即和源代码中一样设置strain_factor = 1.25），边界结果无明显问题。但当我设置应变率下降时（将strain_factor改为0.8），在燃料侧的grid=0处，其结果中的组分组成随应变率下降而逐渐偏离边界条件。我的燃料设置为Y_O2:0.2，Y_NH3:0.1417，Y_N2:0.6583。取应变率下降循环第50次的结果而言，grid=0处的O2仅为0.15左右，NH3也偏差很大。这个可能是什么原因引起的？

示例中的源代码部分如下：

# PART 3: STRAIN RATE LOOP

# Compute counterflow diffusion flames at increasing strain rates at 1 bar
# The strain rate is assumed to increase by 25% in each step until the flame is
# extinguished
strain_factor = 1.25

# Exponents for the initial solution variation with changes in strain rate
# Taken from Fiala and Sattelmayer (2014)
exp_d_a = - 1. / 2.
exp_u_a = 1. / 2.
exp_V_a = 1.
exp_lam_a = 2.
exp_mdot_a = 1. / 2.

# Restore initial solution
file_name, entry = names("initial-solution")
f.restore(file_name, name=entry)

# Counter to identify the loop
n = 0
# Do the strain rate loop
while np.max(f.T) > temperature_limit_extinction:
    n += 1
    print('strain rate iteration', n)
    # Create an initial guess based on the previous solution
    # Update grid
    f.flame.grid *= strain_factor ** exp_d_a
    normalized_grid = f.grid / (f.grid[-1] - f.grid[0])
    # Update mass fluxes
    f.fuel_inlet.mdot *= strain_factor ** exp_mdot_a
    f.oxidizer_inlet.mdot *= strain_factor ** exp_mdot_a
    # Update velocities
    f.set_profile('velocity', normalized_grid,
                  f.velocity * strain_factor ** exp_u_a)
    f.set_profile('spread_rate', normalized_grid,
                  f.spread_rate * strain_factor ** exp_V_a)
    # Update pressure curvature
    f.set_profile('lambda', normalized_grid, f.L * strain_factor ** exp_lam_a)
    try:
        # Try solving the flame
        f.solve(loglevel=0)

CFD中文网

Rachel0096