9
$\begingroup$

I am having trouble implementing a model from a publication.

Huang, K-L.; Holsen, T.M.; Selman, J.R. Ind. Eng. Chem. Res. 2003, 42, 15, 3620–3625

scihub link: https://sci-hub.se/10.1021/ie030109q

I want to model the diffusion of bisulfate through a nafion membrane using my own cell parameters. The system is two 10 L tanks separated by a ~200 micrometer nafion membrane. One tank is more concentrated in bisulfate than the other, and diffusion will occur through the membrane.

The equations and boundary conditions given by the authors are

The equations given by the authors to model the system

$$ 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) $$

My python code is below, including the parameters I am using:

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

Because $\frac{dC_{c/d}}{dt}$ do not depend on $\frac{dC_{c/d}}{dx}$, so I wont need to calculate the concentration gradient in the compartments, I only included one mesh point for $C_{c}$ and one mesh point for $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

This code seems to give reasonable results but I am not sure. The problem is that my dx value is very small because the nafion width is on the order of microns but I want to model out to 60 days.I have found that I must keep the value of dt near dx or the results become unstable (I think that's whats happening). If $dt = dx * 10^3$, then the results quickly become NaN.

Is this due to an error in my code or is the solution really that unstable? If I want to model out to 60 days, I would like dt to be larger. The publication said "The model equations were solved using a combination of the orthogonal collocation method and numerical integration using Gear's method." They present modeled data out to 50 days. enter image description here

For reference here is a diagram of the physical system being modeled.

enter image description here

Update

After continuing to work on this I have updated my code. I used the authors system parameters (V, A, C_0 etc.) to ensure that I can recreate their plots as a check.

I am now using Runge-Kutta with adaptive time steps. The adaptive timesteps seem to be key, because at low times, very far from membrane equilibration, dt must be very small because there are very large gradients in the membrane ($\frac{dC_m}{dx}$ very very large). But then dt must be allowed to increase as diffusion occurs (and $\frac{dC_m}{dx}$ decreases) so that modeling longer than a second does not take an hour.

I incorporated boundary conditions 1 and 2 (eqs 14d/e) by taking the time derivative of the boundary conditions so then $\frac{dC_m}{dt} = K\cdot \frac{dC_c}{dt}$ and then I can use eq 13b to update $C_m$ at x=0 and perform a similar operation for x = L and eq 13c

Cell setup using the author parameters:

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

Functions for solving the equations:

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

And here is the jupyter cell I use to implement and play around with the above code:

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)

With this code, I can replicate Figure 4 from the publication.

enter image description here

The ratio of bisulfate in the dilute tank:conc. tank that I model agrees with the ratios they got, so I feel okay about my code. I think there may be an error in how I implemented the adaptive time step and I think I can speed it up. It takes 9 minutes to model 5 days of diffusion on my machine with no scaling constant (x_c = t_c = 1).

$\endgroup$
10
$\begingroup$

I have found that I must keep the value of dt near dx or the results become unstable

This behavior you have noticed is known as the Courant–Friedrichs–Lewy (CFL) condition. Indeed, there are timestep restrictions that are related to spatial refinement and choice of time integration method. Due to the diffusion term in 13a, the timestep must satisfy an inequality of the form $$ dt \leq C dx^2. $$ The exponent of $2$ comes from the second derivative in space. In other words, the stiffness of the problem grows quadratically with the spatial refinement. The constant C depends on the diffusion rate and stability region of the time integration method.

Is this due to an error in my code or is the solution really that unstable?

I haven't reviewed you code carefully enough to says it is a correct implementation of the PDE model, but the instability for large stepsize is expected. The main issue is the use of a suboptimal numerical method. From the line of code

cell_update = cells[k-1] + ((2 * D * dt) / dx**2) * (cell_iplus + cell_iminus - 2 * cell) # eq 13a

it looks like you are using a forward Euler discretization in time (although I not sure why there is a 2 multiplying D*dt). This is one of the simplest but least stable methods. It leads to a relatively small $C$.

The model equations were solved using a combination of the orthogonal collocation method and numerical integration using Gear's method.

Gear's method (aka backward differentiation formula) is an example of an implicit time integration. At the cost of having to solve (non)linear systems of equations, implicit methods can bypass the CFL stepsize restriction. These are ideal for solving stiff problems where explicit methods require an impractically large number of timesteps. You should look into implementing implicit methods like backward Euler, implicit Runge-Kutta, BDF, Rosenbrock, etc. If you want to stick with explicit schemes, use a higher order Runge-Kutta with better stability or a Runge-Kutta-Chebyshev method.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.