PDES+ODEs with ModelingToolkit

New to Julia
1

I am simulating an electrochemical for Li+ purification using membranes. The problem includes several PDEs modeling the diffusion and migration in membrane coupled with ODEs representing the balances at chambers. The code currently is: using ModelingToolkit, DifferentialEquations, MethodOfLines, PDEBase, DomainSets
using SciMLBase, LinearAlgebra, Plots, Sundials

----------------------

1. Physical and System Constants

----------------------

@parameters F D_Li D_OH RT z_Li z_OH k_ext V_anodo V_catodo A
@parameters F_LiCl_rec LiCl_feed_rec F_LiOH_rec Li_OH_rec F_Wat_Cat kf kr dPhi L
@parameters C_LiAIn C_LiCIn C_OHAIn C_OHCIn

----------------------

1. Physical and System Constants (numeric)

----------------------

const F_val = 96485.3329 # Faraday constant (C/mol)
const R_gas_val = 8.3145 # Ideal gas constant (J/mol/K)
const Temp_val = 341.15 # Temperature (K)
const D_Li_val = 9.0e-10 # Diffusion coefficient of Li+ (m^2/s)
const D_OH_val = 2.7e-10 # Diffusion coefficient of OH- (m^2/s)
const z_Li_val = 1 # Charge of Li+ ion
const z_OH_val = -1 # Charge of OH- ion
const L_val = 1.5e-3 # Length of the domain (m)
const N_val = 20 # Number of spatial discretization points
const k_ext_val = 2.0 # External mass transfer coefficient (m/s)

----------------------

2. Process Input/Flow Constants (numeric)

----------------------

const V_anodo_val = 15.0
const V_catodo_val = 7.0
const A_val = 1.5
const kf_val = 1.0
const kr_val = 0.97
const F_LiCl_rec_val = 192/3600
const F_LiCl_fre_val = 0.15/3600
const F_Wat_And_val = 0.17/3600
const F_LiOH_rec_val = 192/3600
const F_Wat_Cat_val = 0.60/3600
const LiCl_feed_rec_val = 3.2
const LiCl_feed_fre_val = 11.83
const Li_OH_rec_val = 1.98

----------------------

3. Boundary/Initial Concentrations (numeric)

----------------------

const C_LiAIn_val = 0.25 # Li+ concentration at Anode Inlet
const C_LiCIn_val = 0.15 # Li+ concentration at Cathode Inlet
const C_OHAIn_val = 1e-4 # OH- concentration at Anode Inlet
const C_OHCIn_val = 0.15 # OH- concentration at Cathode Inlet

----------------------

4. Derived Constants (numeric)

----------------------

const dx_val = L_val/(N_val-1)
const RT_val = R_gas_val * Temp_val
const RT_over_F_val = RT_val / F_val

Electric potentials

const Phi_AIn_val = RT_over_F_val * log(C_LiAIn_val / C_LiCIn_val)
const Phi_CIn_val = 0.0

Gradient of potential along x (numeric)

const dPhi_dx_val = (Phi_CIn_val - Phi_AIn_val)/L_val

Helper functions for Initial Conditions

C_Li_init_func(x_val) = C_LiAIn + (C_LiCIn-C_LiAIn)*x_val/L
C_OH_init_func(x_val) = C_OHAIn + (C_OHCIn-C_OHAIn)*x_val/L

----------------------

5. Symbolic PDEs Definition

----------------------

@variables x
@independent_variables t

PDE variables

@variables C_Li(..) C_OH(..)

ODE variables

@variables C_LiA(..) C_LiC(..) C_LiOHC(..) C_OHC(..)

Fluxes and reaction rate (functions of time)

@variables flux_Li_anod(..) flux_Li_cathod(..) flux_OH_cathod(..) reaction_rate_cathode(..)

Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

PDEs

pde_eqs = [
Dt(C_Li(x,t)) ~ D_LiDx(D_x(C_Li(x,t))) + D_Li(F/RT)dPhiDx(C_Li(x,t)),
Dt(C_OH(x,t)) ~ D_OHDx(D_x(C_OH(x,t))) - D_OH(F/RT)dPhiDx(C_OH(x,t))
]

Flux/reaction equations

anode_flux_eq_Li = flux_Li_anod(t) ~ -D_Li * Dx(C_Li(0, t)) - (z_Li * D_Li * F / RT) * C_Li(0, t) * dPhi
cathod_flux_eq_Li = flux_Li_cathod(t) ~ -D_Li * Dx(C_Li(L, t)) - (z_Li * D_Li * F / RT) * C_Li(L, t) * dPhi
cathod_flux_eq_OH = flux_OH_cathod(t) ~ -D_OH * Dx(C_OH(L, t)) - (z_OH * D_OH * F / RT) * C_OH(L, t) * dPhi
reaction_rate_eq = reaction_rate_cathode(t) ~ kf * C_LiC(t) * C_OHC(t) - kr * C_LiOHC(t)

ODEs

ode_eqs = [
Dt(C_LiA(t)) ~ (F_LiCl_rec/V_anodo)(LiCl_feed_rec - C_LiA(t)) - (A/V_anodo) * flux_Li_anod(t),
Dt(C_LiC(t)) ~ (F_LiOH_rec/V_catodo)(Li_OH_rec - C_LiC(t)) + (A/V_catodo) * flux_Li_cathod(t) - reaction_rate_cathode(t),
Dt(C_LiOHC(t)) ~ (F_LiOH_rec/V_catodo) * Li_OH_rec - ((F_LiOH_rec + F_Wat_Cat) / V_catodo) * C_LiOHC(t) + reaction_rate_cathode(t),
Dt(C_OHC(t)) ~ (F_LiOH_rec/V_catodo) * Li_OH_rec - ((F_LiOH_rec + F_Wat_Cat) / V_catodo) * C_OHC(t) + (A/V_catodo) * flux_OH_cathod(t) - reaction_rate_cathode(t)
]

Combine all equations

eqs = vcat(pde_eqs, [anode_flux_eq_Li, cathod_flux_eq_Li, cathod_flux_eq_OH, reaction_rate_eq], ode_eqs)

bcs = [
# Initial Conditions (t=0) for PDEs
C_Li(x, 0) ~ C_Li_init_func(x),
C_OH(x, 0) ~ C_OH_init_func(x),

# Initial Conditions for the ODE variables
C_LiA(0) ~ C_LiAIn,
C_LiC(0) ~ C_LiCIn,
C_OHC(0) ~ C_OHCIn,
C_LiOHC(0) ~ 0.0,

# Boundary Flux Conditions (x=0, Anode side)
Dx(C_Li(0,t)) + (z_Li*F/RT)*C_Li(0,t)*dPhi ~ -k_ext/D_Li*(C_Li(0,t) - C_LiA(t)),
Dx(C_OH(0,t)) + (z_OH*F/RT)*C_OH(0,t)*dPhi ~ -k_ext/D_OH*(C_OH(0,t) - C_OHAIn),

# Boundary Flux Conditions (x=L, Cathode side)
Dx(C_Li(L,t)) + (z_Li*F/RT)*C_Li(L,t)*dPhi ~ -k_ext/D_Li*(C_Li(L,t) - C_LiC(t)),
Dx(C_OH(L,t)) + (z_OH*F/RT)*C_OH(L,t)*dPhi ~ -k_ext/D_OH*(C_OH(L,t) - C_OHC(t))

]

----------------------

6. System Construction

----------------------

domains = [x ∈ Interval(0.0, L_val), t ∈ Interval(0.0, 100.0)]

Collect all parameters needed in the defaults dictionary

all_params = Dict(
F => F_val,
D_Li => D_Li_val,
D_OH => D_OH_val,
RT => RT_val,
z_Li => z_Li_val,
z_OH => z_OH_val,
k_ext => k_ext_val,
V_anodo => V_anodo_val,
V_catodo => V_catodo_val,
A => A_val,
F_LiCl_rec => F_LiCl_rec_val,
LiCl_feed_rec => LiCl_feed_rec_val,
F_LiOH_rec => F_LiOH_rec_val,
Li_OH_rec => Li_OH_rec_val,
F_Wat_Cat => F_Wat_Cat_val,
C_LiAIn => C_LiAIn_val,
C_LiCIn => C_LiCIn_val,
C_OHAIn => C_OHAIn_val,
C_OHCIn => C_OHCIn_val,
kf => kf_val,
kr => kr_val,
dPhi => dPhi_dx_val,
L => L_val,
)

@named pde_system = PDESystem(
eqs,
bcs,
domains,
[x, t],
[C_Li(x,t), C_OH(x,t), C_LiA(t), C_LiC(t), C_OHC(t), C_LiOHC(t),
flux_Li_anod(t), flux_Li_cathod(t), flux_OH_cathod(t), reaction_rate_cathode(t)],
collect(keys(all_params));
defaults=all_params
)

----------------------

7. Discretization and Solution

----------------------

@show dx_val
xspan = 0.0:Float64(dx_val):Float64(L_val)
discretization = MOLFiniteDifference([x => dx_val], t; advection_scheme=UpwindScheme())
prob_ode = discretize(pde_system, discretization)

Solve the problem

sol = solve(prob_ode, TRBDF2(), saveat=1.0, reltol=1e-6, abstol=1e-8)

----------------------

8. Results Visualization

----------------------

Extract discrete x values

x_vals = range(0, stop=L, length=N)

Extract concentration matrices

C_Li_mat = sol[C_Li(x, t)]
C_OH_mat = sol[C_OH(x, t)]

Target times for plotting

target_times = [0.0, 25.0, 50.0, 75.0, 100.0]

Create plots for each time

p = plot(layout=(length(target_times), 1), size=(800, 1000))

for (i, t_target) in enumerate(target_times)
# Find closest time index
t_idx = argmin(abs.(sol.t .- t_target))

# Extract profiles at that time
Li_profile = C_Li_mat[:, t_idx]
OH_profile = C_OH_mat[:, t_idx]

# Plot both species
plot!(
    p[i],
    x_vals, Li_profile,
    label="Li⁺",
    lw=2,
    color=:blue,
    title="t = $(t_target) s",
    legend=:topright
)
plot!(
    p[i],
    x_vals, OH_profile,
    label="OH⁻",
    lw=2,
    color=:red,
    legend=:topright
)

# Format axes
xlabel!(p[i], "Position (m)")
ylabel!(p[i], "Concentration (mol/m³)")

end

Save the figure

savefig(p, “concentration_profiles.png”)
but I got an error: ERROR: LoadError: ArgumentError: invalid index: nothing of type Nothing Stacktrace: [1] to_index(i::Nothing) @ Base .\indices.jl:315 [2] to_index(A::Vector{Symbolics.VarDomainPairing}, i::Nothing) @ Base .\indices.jl:292 [3] to_indices @ .\indices.jl:368 [inlined] [4] to_indices @ .\indices.jl:360 [inlined] [5] getindex @ .\abstractarray.jl:1312 [inlined] [6] (::PDEBase.var"#54#62"{Vector{Symbolics.VarDomainPairing}})(x::SymbolicUtils.BasicSymbolic{Real}) The error is in prob_ode = discretize(pde_system, discretization) line. Previously I solved the problem including only the membrane. But when the balances are included the BC’s of the PDEs are intelinked with them. If someone has ideas how to handle the problem, I will be very grateful.

2

Small remark:
If you share code, please enclose it in triple backticks, like ```

That makes it much easier to read your code, but also to copy-and-paste it for people who want to reproduce your problem.

1 Like
3 
using SciMLBase, LinearAlgebra, Plots, Sundials

# ----------------------
# 1. Physical and System Constants
# ----------------------
@parameters F D_Li D_OH RT z_Li z_OH k_ext V_anodo V_catodo A 
@parameters F_LiCl_rec LiCl_feed_rec F_LiOH_rec Li_OH_rec F_Wat_Cat kf kr dPhi L
@parameters C_LiAIn C_LiCIn C_OHAIn C_OHCIn

# ----------------------
# 1. Physical and System Constants (numeric)
# ----------------------
const F_val = 96485.3329          # Faraday constant (C/mol)
const R_gas_val = 8.3145          # Ideal gas constant (J/mol/K)
const Temp_val = 341.15            # Temperature (K)
const D_Li_val = 9.0e-10           # Diffusion coefficient of Li+ (m^2/s)
const D_OH_val = 2.7e-10           # Diffusion coefficient of OH- (m^2/s)
const z_Li_val = 1                  # Charge of Li+ ion
const z_OH_val = -1                 # Charge of OH- ion
const L_val = 1.5e-3                # Length of the domain (m)
const N_val = 20                     # Number of spatial discretization points
const k_ext_val = 2.0               # External mass transfer coefficient (m/s)

# ----------------------
# 2. Process Input/Flow Constants (numeric)
# ----------------------
const V_anodo_val = 15.0
const V_catodo_val = 7.0
const A_val = 1.5
const kf_val = 1.0
const kr_val = 0.97
const F_LiCl_rec_val = 192/3600
const F_LiCl_fre_val = 0.15/3600
const F_Wat_And_val = 0.17/3600
const F_LiOH_rec_val = 192/3600
const F_Wat_Cat_val = 0.60/3600
const LiCl_feed_rec_val = 3.2
const LiCl_feed_fre_val = 11.83
const Li_OH_rec_val = 1.98

# ----------------------
# 3. Boundary/Initial Concentrations (numeric)
# ----------------------
const C_LiAIn_val = 0.25           # Li+ concentration at Anode Inlet
const C_LiCIn_val = 0.15           # Li+ concentration at Cathode Inlet
const C_OHAIn_val = 1e-4           # OH- concentration at Anode Inlet
const C_OHCIn_val = 0.15           # OH- concentration at Cathode Inlet

# ----------------------
# 4. Derived Constants (numeric)
# ----------------------
const dx_val = L_val/(N_val-1)
const RT_val = R_gas_val * Temp_val
const RT_over_F_val = RT_val / F_val

# Electric potentials
const Phi_AIn_val = RT_over_F_val * log(C_LiAIn_val / C_LiCIn_val)
const Phi_CIn_val = 0.0

# Gradient of potential along x (numeric)
const dPhi_dx_val = (Phi_CIn_val - Phi_AIn_val)/L_val

# Helper functions for Initial Conditions
C_Li_init_func(x_val) = C_LiAIn + (C_LiCIn-C_LiAIn)*x_val/L
C_OH_init_func(x_val) = C_OHAIn + (C_OHCIn-C_OHAIn)*x_val/L

# ----------------------
# 5. Symbolic PDEs Definition
# ----------------------

@variables x
@independent_variables t
# PDE variables
@variables C_Li(..) C_OH(..)

# ODE variables
@variables C_LiA(..) C_LiC(..) C_LiOHC(..) C_OHC(..)

# Fluxes and reaction rate (functions of time)
@variables flux_Li_anod(..) flux_Li_cathod(..) flux_OH_cathod(..) reaction_rate_cathode(..)

Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

# PDEs
pde_eqs = [
    Dt(C_Li(x,t)) ~ D_Li*Dx(D_x(C_Li(x,t))) + D_Li*(F/RT)*dPhi*Dx(C_Li(x,t)),
    Dt(C_OH(x,t)) ~ D_OH*Dx(D_x(C_OH(x,t))) - D_OH*(F/RT)*dPhi*Dx(C_OH(x,t))
]

# Flux/reaction equations
anode_flux_eq_Li = flux_Li_anod(t) ~ -D_Li * Dx(C_Li(0, t)) - (z_Li * D_Li * F / RT) * C_Li(0, t) * dPhi
cathod_flux_eq_Li = flux_Li_cathod(t) ~ -D_Li * Dx(C_Li(L, t)) - (z_Li * D_Li * F / RT) * C_Li(L, t) * dPhi
cathod_flux_eq_OH = flux_OH_cathod(t) ~ -D_OH * Dx(C_OH(L, t)) - (z_OH * D_OH * F / RT) * C_OH(L, t) * dPhi
reaction_rate_eq = reaction_rate_cathode(t) ~ kf * C_LiC(t) * C_OHC(t) - kr * C_LiOHC(t)

# ODEs
ode_eqs = [
    Dt(C_LiA(t)) ~ (F_LiCl_rec/V_anodo)*(LiCl_feed_rec - C_LiA(t)) - (A/V_anodo) * flux_Li_anod(t),
    Dt(C_LiC(t)) ~ (F_LiOH_rec/V_catodo)*(Li_OH_rec - C_LiC(t)) + (A/V_catodo) * flux_Li_cathod(t) - reaction_rate_cathode(t),
    Dt(C_LiOHC(t)) ~ (F_LiOH_rec/V_catodo) * Li_OH_rec - ((F_LiOH_rec + F_Wat_Cat) / V_catodo) * C_LiOHC(t) + reaction_rate_cathode(t),
    Dt(C_OHC(t)) ~ (F_LiOH_rec/V_catodo) * Li_OH_rec - ((F_LiOH_rec + F_Wat_Cat) / V_catodo) * C_OHC(t) + (A/V_catodo) * flux_OH_cathod(t) - reaction_rate_cathode(t)
]

# Combine all equations
eqs = vcat(pde_eqs, [anode_flux_eq_Li, cathod_flux_eq_Li, cathod_flux_eq_OH, reaction_rate_eq], ode_eqs)

bcs = [
    # Initial Conditions (t=0) for PDEs
    C_Li(x, 0) ~ C_Li_init_func(x),
    C_OH(x, 0) ~ C_OH_init_func(x),
    
    # Initial Conditions for the ODE variables
    C_LiA(0) ~ C_LiAIn,
    C_LiC(0) ~ C_LiCIn,
    C_OHC(0) ~ C_OHCIn,
    C_LiOHC(0) ~ 0.0,
    
    # Boundary Flux Conditions (x=0, Anode side)
    Dx(C_Li(0,t)) + (z_Li*F/RT)*C_Li(0,t)*dPhi ~ -k_ext/D_Li*(C_Li(0,t) - C_LiA(t)),
    Dx(C_OH(0,t)) + (z_OH*F/RT)*C_OH(0,t)*dPhi ~ -k_ext/D_OH*(C_OH(0,t) - C_OHAIn),
    
    # Boundary Flux Conditions (x=L, Cathode side)
    Dx(C_Li(L,t)) + (z_Li*F/RT)*C_Li(L,t)*dPhi ~ -k_ext/D_Li*(C_Li(L,t) - C_LiC(t)),
    Dx(C_OH(L,t)) + (z_OH*F/RT)*C_OH(L,t)*dPhi ~ -k_ext/D_OH*(C_OH(L,t) - C_OHC(t))
]
# ----------------------
# 6. System Construction
# ----------------------
domains = [x ∈ Interval(0.0, L_val), t ∈ Interval(0.0, 100.0)]

# Collect all parameters needed in the defaults dictionary
all_params = Dict(
    F => F_val,
    D_Li => D_Li_val,
    D_OH => D_OH_val,
    RT => RT_val,
    z_Li => z_Li_val,
    z_OH => z_OH_val,
    k_ext => k_ext_val,
    V_anodo => V_anodo_val,
    V_catodo => V_catodo_val,
    A => A_val,
    F_LiCl_rec => F_LiCl_rec_val,
    LiCl_feed_rec => LiCl_feed_rec_val,
    F_LiOH_rec => F_LiOH_rec_val,
    Li_OH_rec => Li_OH_rec_val,
    F_Wat_Cat => F_Wat_Cat_val,
    C_LiAIn => C_LiAIn_val, 
    C_LiCIn => C_LiCIn_val,
    C_OHAIn => C_OHAIn_val, 
    C_OHCIn => C_OHCIn_val,
    kf => kf_val,
    kr => kr_val,
    dPhi => dPhi_dx_val,
    L => L_val,
)

@named pde_system = PDESystem(
    eqs,
    bcs,
    domains,
    [x, t],
    [C_Li(x,t), C_OH(x,t), C_LiA(t), C_LiC(t), C_OHC(t), C_LiOHC(t), 
    flux_Li_anod(t), flux_Li_cathod(t), flux_OH_cathod(t), reaction_rate_cathode(t)],
    collect(keys(all_params));
    defaults=all_params
)

# ----------------------
# 7. Discretization and Solution
# ----------------------
@show dx_val
xspan = 0.0:Float64(dx_val):Float64(L_val)
discretization = MOLFiniteDifference([x => dx_val], t; advection_scheme=UpwindScheme())
prob_ode = discretize(pde_system, discretization)

# Solve the problem
sol = solve(prob_ode, TRBDF2(), saveat=1.0, reltol=1e-6, abstol=1e-8)

# ----------------------
# 8. Results Visualization
# ----------------------
# Extract discrete x values
x_vals = range(0, stop=L, length=N)

# Extract concentration matrices
C_Li_mat = sol[C_Li(x, t)]
C_OH_mat = sol[C_OH(x, t)]

# Target times for plotting
target_times = [0.0, 25.0, 50.0, 75.0, 100.0]

# Create plots for each time
p = plot(layout=(length(target_times), 1), size=(800, 1000))

for (i, t_target) in enumerate(target_times)
    # Find closest time index
    t_idx = argmin(abs.(sol.t .- t_target))
    
    # Extract profiles at that time
    Li_profile = C_Li_mat[:, t_idx]
    OH_profile = C_OH_mat[:, t_idx]
    
    # Plot both species
    plot!(
        p[i],
        x_vals, Li_profile,
        label="Li⁺",
        lw=2,
        color=:blue,
        title="t = $(t_target) s",
        legend=:topright
    )
    plot!(
        p[i],
        x_vals, OH_profile,
        label="OH⁻",
        lw=2,
        color=:red,
        legend=:topright
    )
    
    # Format axes
    xlabel!(p[i], "Position (m)")
    ylabel!(p[i], "Concentration (mol/m³)")
end

# Save the figure
savefig(p, "concentration_profiles.png")
1 Like