我无法从出版物中实施模型。

; KL 黄。霍尔森，TM；塞尔曼，JR 工业工程师。化学。水库。2003, 42, 15, 3620–3625

scihub 链接：https ://sci-hub.se/10.1021/ie030109q

我想使用我自己的细胞参数模拟硫酸氢盐通过 nafion 膜的扩散。该系统是两个 10 L 罐，由~200 微米的 Nafion 膜隔开。一个罐比另一个罐中的硫酸氢盐浓度更高，并且会通过膜发生扩散。

作者给出的方程和边界条件是

$$C_{m} = \text{Concentration within nafion membrane} \\ L = \text{thickness of the nafion membrane}\\ D_{e} = \text{Diffusion Constant} \\ V_{c/d} = \text{Volume of the concentrated/dilute compartment}\\ C_{c/d} = \text{Concentration of the concentrated/dilute compartment}\\ A_{m} = \text{Area of nafion membrane connecting the two tanks}\\ K = \text{bisulfite partition between nafion and solution } (0 \leq K \leq 1)$$

我的python代码如下，包括我正在使用的参数：

D = (36e-3 * 10**-4) * np.exp(-26.3e3 / (8.3144 * (273 + 31))) # m^2 / s, Diffusion Constant
A = 14e3 * 10**-4 # m^2, area of the membrane
V = 0.01 # m^3, the volume in each bath (10 L)
Cc_0 = 60e3 # ppm (unitless), the initial concentration in the cathode compartment
Cd_0 = .3 * Cc_0 # ppm, initial concentration of the anode compartment

membrane_thickness = 183e-6 # m, the thickness of the nafion membrane in meters
Vc = Vd = V # Setting the volume of each tank, tank volumes are both 10 L

因为不依赖于，所以我不需要计算隔间中的浓度梯度，我只包括了一个网格一个网格点，C_$\frac{dC_{c/d}}{dt}$$\frac{dC_{c/d}}{dx}$$C_{c}$$C_{d}$

L = membrane_thickness
n_x = 1e3
dx = L / n_x
fresh = np.zeros((int(n_x)+2, 1)) # added 2 to provide a mesh point for the concentrated and dilute compartments 
fresh[0] = Cc_0 # initialize the concentrated compartment (eq 14b in the paper)
fresh[-1] = Cd_0 # initialize the dilute compartment (eq 14c)

@jit
def solve(n_t):
    cell = fresh.copy()
    cells=[cell]
    kappa_c = (A * D)  / Vc
    kappa_d = -(A * D)  / Vd
    for k in range(1,n_t):
        cell_iplus = np.roll(cell, -1)
        cell_iminus = np.roll(cell, 1)
        
        cell_update = cells[k-1] + ((2 * D * dt) / dx**2) * (cell_iplus + cell_iminus - 2 * cell) # eq 13a
        cell_update[1] = k_partition * cell[0] # eq 14d, boundary condition 1
        cell_update[-2]= k_partition * cell[-1]# eq 14e, boundary condition 2
        
        dCm_dx = np.diff(cell, axis=0) / dx # d(Cm)/dx
        
        dCm_dx_c = dCm_dx[1] # for eq 13b, index 1 corresponds to x=0 in the membrane
        dCm_dx_d = dCm_dx[-2]# for eq 13c, index -2 corresponds to x=L in the membrane
        
        Cc = (kappa_c * dCm_dx_c * 2 * dt) + cells[k-1][0] # eq 13b, scaler value to update the first mesh point
        Cd = (kappa_d * dCm_dx_d * 2 * dt) + cells[k-1][-1]# eq 13c, scaler value to update the final mesh point 
        
        cell_update[0] = Cc
        cell_update[-1]= Cd
        cell = cell_update
        cells.append(cell_update)
    return cells

这段代码似乎给出了合理的结果，但我不确定。问题是我的 dx 值非常小，因为 nafion 宽度在微米量级，但我想建模到 60 天。我发现我必须将 dt 的值保持在 dx 附近，否则结果会变得不稳定（我认为这就是正在发生的事情）。如果，则结果很快变为 NaN。$dt = dx * 10^3$

这是由于我的代码中的错误还是解决方案真的那么不稳定？如果我想建模到 60 天，我希望 dt 更大。该出版物称“模型方程是使用正交搭配方法和使用 Gear 方法的数值积分的组合来求解的。” 他们提供长达 50 天的建模数据。

作为参考，这里是正在建模的物理系统的图表。

更新

在继续工作后，我更新了我的代码。我使用了作者的系统参数（V、A、C_0 等）来确保我可以重新创建他们的图作为检查。

我现在使用具有自适应时间步长的 Runge-Kutta。自适应时间步长似乎是关键，因为在离膜平衡很远的低时刻，dt 必须非常小，因为膜中有非常大的梯度（非常非常大）。但是随后必须允许 dt 随着扩散的发生而增加（并且减少），这样一秒以上的建模就不需要一个小时。$\frac{dC_m}{dx}$$\frac{dC_m}{dx}$

我通过取边界条件的时间导数合并了边界条件 1 和 2 (eqs 14d/e)，因此然后我可以使用 eq 13b 更新$\frac{dC_m}{dt} = K\cdot \frac{dC_c}{dt}$$C_m$并对 x = L 和 eq 13c 执行类似的操作

使用作者参数的单元设置：

D = (36e-3 * 10**-4) * np.exp(-26.3e3 / (8.314 * (273 + 25))) # m^2 / s, Diffusion Constant
A = .001 # m^2, area of the membrane
V = 0.001 # m^3, the volume in each bath
Cc_0 = 9.606e4 # ppm (unitless), the initial concentration in the cathode compartment
Ca_0 = 0 #ppm, initial concentration of the anode compartment

membrane_thickness = 183e-6 # m, the thickness of the nafion membrane in meters
Vc = Va = V # Setting volume of each tank, tank volumes are both 1 L

C_c = Cc_0 # Scaling constant for concentration in the concentrated tank set to inital concentration
Ca_c = C_c * (Vc / Va) # scaling constant for concentration in the dilute tank

k_partition = .23

求解方程的函数：

def f(mesh, dx_, x_c = 1, t_c = 1): #f = dCi / dt at each mesh point, performs calculations for membrane and both tanks
    # x_c and t_c are scaling constants to make equations dimensionless, default is 1 (not scaled)
    kappa_13b = (A * D * t_c)/(Vc * x_c) 
    kappa_13c = -(A * D * t_c)/(Va * x_c)
    kappa_13a = ( D * t_c ) / x_c**2
    
    dCm_dx = np.gradient(mesh[padding:-padding], dx_, axis = 0) # for eq 13b/c, only use the membrane region in the gradient
    
    dx0_ = dCm_dx[0] # for eq 13b. dCm/dx at x = 0. A scaler value
    dxl_ = dCm_dx[-1] #for eq 13c. dCm/dx at x = L. A scaler value
    
    # for eq 13a
    mesh_plus = np.roll(mesh, -1) # gives mesh value at (x + 1) for each mesh point. A vector 
    mesh_minus = np.roll(mesh, 1) # gives mesh value at (x - 1) for each mesh point. A vector
    # The first and last entry of mesh_plus and mesh_minus are invalid values. But are not used in eq 13a, so no error.
    
    second_deriv = ((mesh_plus+mesh_minus - (2*mesh)) / dx_**2)[1:-1] # for eq 13a. use entire mesh to calculate second derivative, keep only membrane region
    # The first and last values of the 2nd derivative within the membrane region are incorrect because they use values from the tanks during their calculation
    # However, these two values will be replaced with boundary conditions so no error
    
    #put calculations together to form the RHS of eqns 13a, 13b, and 13c. A single vector is returned
    fun = np.empty(mesh.shape) # the return mesh which will contain dCi/dt at each mesh point
     
    # f(t,x) for the cathode tank
    fun[:padding] = kappa_13b * dx0_ # eq 13b, f = dC/dt for the cathode tank (concentrated tank)
    
    #f(t,x) for the membrane region
    fun[padding:-padding] = kappa_13a * second_deriv # eq 13a, f = dC/dt to update interior of membrane
    
    #f(t,x) for the anode tank
    fun[-padding:] = kappa_13c * dxl_ # eq 13c, f = dC/dt for the anode tank (dilute tank)
    
    # Boundary Conditions
    # f(t,x) = dC/dt at x = 0 and x = L  set by the boundary conditions BC1 and BC2 from the paper
    fun[padding] = k_partition * kappa_13b * dx0_ # time derivative of BC1 (eq 14d) gives dCm/dt = k_partition * dCc/dt at x=0
    fun[-(padding+1)] = k_partition * kappa_13c * dxl_ # time derivative of BC2 (eq 14e) gives dCm/dt = k_partition * dCa/dt at x=L
    
    return fun

@jit
def rk4_step(f, mesh, dx_, dt_, x_c, t_c):
    
    k1 = dt_ * f(mesh, dx_, x_c, t_c)
    k2 = dt_ * f(mesh + .5*k1, dx_, x_c, t_c)
    k3 = dt_ * f(mesh + .5*k2, dx_, x_c, t_c)
    k4 = dt_ * f(mesh + k3, dx_, x_c, t_c)
    
    return (k1 + 2*k2 + 2*k3 + k4)/6

@jit
def solver(mesh, T, dt, dx, t_c=1, x_c=1, target_accuracy = 1e-5):
    cell = mesh.copy()
    
    # t_c and x_c are scaling constants
    # default scaling constant is 1 (no scaling)
    dt_ = dt / t_c
    dx_ = dx / x_c
    
    T_ = T / t_c
    
    t_ = 0
    cells = [(t_c*t_, t_c*dt_, cell)] # (current time, timestep for stepping to next point, cell state)
    while t_ < T_:
        # adaptive timestep method
        #     first test_step starts at t and steps twice using dt_ to get to test_step_1
        single_step  = cell + rk4_step(f, cell, dx_, dt_, x_c, t_c)
        test_step_1 = single_step + rk4_step(f, single_step, dx_, dt_, x_c, t_c)
        
        #     second test_step starts at t and steps once using 2*dt_ to get to test_step_2
        test_step_2 = cell + rk4_step(f, cell, dx_, 2*dt_, x_c, t_c)
        
        difference = np.abs(test_step_2 - test_step_1) # a vector of differences for each mesh point
        largest_diff = np.max(difference) # want to look at the worst discrepency between the two test steps
        
        if largest_diff != 0: # if largest_diff == 0, then dt_ does not need to be updated
            rho = (30 * dt_ * target_accuracy) / largest_diff # ratio of target accuracy and actual accuracy
            if rho > 1: # accuracy is better at every mesh point than required, keep data and update dt_ to a larger step size
                t_ = t_ + dt_
            dt_ = np.power(rho, .25) * dt_ # calculate correct time step based on desired accuracy using the worst rho value
            #dt_ = min(dt_calc, 2*dt_) # do no allow the new dt_ be more than twice the original dt_
            
            if rho < 1: # dt_ was too large for target accuracy, calculate step again using new dt_
                single_step = cell + rk4_step(f, cell, dx_, dt_, x_c, t_c)
                t_ = t_+dt_
        else:
            t_ = t_ + dt_
        
        new_cell = single_step
        cells.append((t_c*t_, t_c*dt_, new_cell))        
        cell = new_cell
    
    return cells

这是我用来实现和使用上述代码的 jupyter 单元：

x_c = membrane_thickness# scaling constant for position
t_c = np.power(x_c, 2) / D # scaling constant for time

membrane_mesh_size = 15 # found to strongly affect the solving time

dx = membrane_thickness / membrane_mesh_size
dt = 1e-12 # starting timestep, the adaptive timestep in the solver will correct this to a more appropriate value

padding = 1 # the number of extra mesh points padding the start and end of the membrane mesh
            # holds values for the tanks

n_x =int( membrane_mesh_size + 2*padding)

fresh = np.zeros((n_x, 1))
fresh[:padding] = Cc_0 #/ C_c # pad mesh points set to tank starting concentration
fresh[-padding:]= Ca_0 #/ Ca_c# terminal pad mesh points set to anode tank starting concentration

fresh[padding] = k_partition * fresh[0] # starting membrane conditions at x=0
fresh[-(padding + 1)] = k_partition * fresh[-1] # starting membrane conditions at x=L

hours = 24 * 5
T = hours * 3600

t_c = x_c = 1

%time _ = solver(fresh, T, 1e-6, dx, t_c, x_c, target_accuracy=1e-6)

使用此代码，我可以复制出版物中的图 4。

稀释罐中硫酸氢盐的比例：浓。我建模的坦克与他们得到的比率一致，所以我对我的代码感觉很好。我认为我实现自适应时间步长的方式可能存在错误，我认为我可以加快速度。在没有缩放常数 (x_c = t_c = 1) 的情况下，在我的机器上模拟 5 天的扩散需要 9 分钟。