part2

import math


def get_r(x, t, xso2):
    R = 1.987
    Keff = 0

    if 693.15 <= t < 748.15:
        Keff = 7.6915 * (10 ** 18) * math.exp(-76062 / (R * t))
    if 748.15 <= t <= 873.15:
        Keff = 1.5128 * (10 ** 7) * math.exp(-35992 / (R * t))

    K = 2.3 * (10 ** (-8)) * math.exp(27200 / (R * t))
    Kp = 2.26203 * (10 ** (-5)) * math.exp(11295.3 / t)

    Pso2 = (xso2 - xso2 * x) / (1 - xso2 * x / 2)
    Pso3 = (xso2 * x) / (1 - xso2 * x / 2)
    Po2 = (0.17 - xso2 - xso2 * x / 2) / (1 - xso2 * x / 2)

    r1 = Po2 * Pso2 / Pso3
    r2 = Pso3 / (Pso2 * math.sqrt(Po2) * Kp)
    B = 48148 * math.exp(-7355.5 / t)

    r3 = math.sqrt(B + (B - 1) * (1 - x) / x) + math.sqrt(K * (1 - x) / x)

    r = Keff * K * r1 * (1 - r2 * r2) / (r3 * r3)

    return r


def get_y(x, t, xso2):
    h = 0.0001
    y = (get_r(x, t + h, xso2) - get_r(x, t - h, xso2)) / (2 * h)
    return y


def get_t(t, x0, x):
    H = -23135
    Cp = 254.9
    rou = 0.500
    c = 1.282

    l = -H * c / (rou * Cp)
    y = t + l * (x - x0)
    return y


def fun1(x, t, xso2):

    y = - get_y(x, t, xso2) / (get_r(x, t, xso2) * get_r(x, t, xso2))
    return y


def jifen(x, t0, xso2):
    sum = 0.0
    x1 = x
    t2 = 693.15
    h = 0.0001
    xout = 0

    t1 = get_t(t0, x, x1)
    x2 = x1 + h / 10
    t2 = get_t(t0, x, x2)

    sum += h * (fun1(x1, t1, xso2) + fun1(x2, t2, xso2)) / 20
    x1 = x2
    while sum < 0:
        t1 = get_t(t0, x, x1)
        x2 = x1 + h / 10
        t2 = get_t(t0, x, x2)

        if t2 > 873.15:
            xout = x1
            break
        sum += h * (fun1(x1, t1, xso2) + fun1(x2, t2, xso2)) / 20
        x1 = x2
        xout = x1 - h / 10
    return xout


def weijifen(xin, xou, tin, xso2):
    x1 = xin
    sum = 0.0
    h = 0.0001

    t1 = get_t(tin, xin, x1)
    x2 = x1 + h
    t2 = get_t(tin, xin, x2)

    sum += (1 / get_r(x1, t1, xso2) + 1 / get_r(x2, t2, xso2)) * h / 1000
    x1 = x2

    while x2 <= xou:
        t1 = get_t(tin, xin, x1)
        x2 = x1 + h
        t2 = get_t(tin, xin, x2)
        if t2 >= 873.15:
            break
        else:
            sum += (1 / get_r(x1, t1, xso2) + 1 / get_r(x2, t2, xso2)) * h / 1000
            x1 = x2
    wcat = sum * 131 * 1000 / 3600
    return wcat


def write_2_file(t, xso2):
    h = 0.0001
    wsum = 0.0
    x0 = 0.0001
    t0 = t

    with open('data-2.txt', 'w') as f:
        f.write("1   tin= %f    xin= %e    " % (t, x0))
        for i in range(4):
            xout = jifen(x0, t0, xso2)
            tout = get_t(t0, x0, xout)
            f.write("\n    tout=%f    xout=%f    " % (tout, xout))
            wcat = weijifen(x0, xout, t0, xso2)
            f.write("\n    Wcat=%f      \n" % wcat)
            wsum += wcat
            t1 = 693.15
            t1 = t1 + 0.01

            while (math.fabs((10 ** 5) * get_r(xout, t1, xso2) - (10 ** 5) * get_r(xout, tout, xso2))) > h:
                t1 = t1 + 0.01

            x0 = xout
            t0 = t1
            if i != 3:
                f.write("%d   tin= %f    xin= %f" % (i + 2, t0, xout))

        f.write("    Wsum=%f   \n" % wsum)
        t -= 0.1
        f.write("\n")

    if x0 <= 0.98:
        write_2_file(t, xso2)


def main():
    xso2 = 0.08
    t00 = 719
    write_2_file(t00, xso2)


if __name__ == '__main__':
    main()

u=1710533078,72194801&fm=26&gp=0.jpg

==============================
运行结果：
==============================

1   tin=717.800000    xin=1.000000e-04    
    tout=873.148833    xout=0.667660    
    Wcat=8349.667114      

2   tin= 723.600000    xin= 0.667660
    tout=778.519894    xout=0.903660    
    Wcat=11605.588273      

3   tin= 715.920000    xin= 0.903660
    tout=729.505692    xout=0.962040    
    Wcat=21557.219801      

4   tin= 693.190000    xin= 0.962040
    tout=697.388114    xout=0.980080    
    Wcat=50379.460102      

    Wsum=91891.935289   


Process finished with exit code 0