I am trying to model a chromatography column using ModelingToolkit and MethodOfLines, and I can't get past this error message.
From my testing the problematic term is - (1 - Īµ_column)*(3/rp)*Deff*dz(cā(t, z, rp)) - commenting it out discretizes and solves the model. To be clear this should be the the partial derivative evaluated at rp, I'm not sure if this is how it should be coded.
I don't understand how to fix it - previous reports of the same error were due to mixed derivatives, but I don't think I have them do I?
To be clear: the numbers in the Any[...] vector belong to the z domain I think, changing zL changes them too.
I'd be very grateful for some help!
Minimal Reproducible Example š
@parameters t, z, r
@variables c(..), cā(..)
zL = 0.1 # [m] - length of the column
dp = 85e-6 # [m] diameter of resin particle
rp = dp/2 # [m] radius of resin particle
Īµ_column = 0.36 # [-] void fraction of column
Īµ_particle = 0.88 # [-] void fraction of resin particle
kfit = 0.083 # [-] fitting parameter taking resin particle pore diameter, antibody molecule diameter and pore tortuosity into account
qmax = 5.47 # [mol m-3] maximum binding capacity of affinity resin
kd = 7.7e-4 # [mol m-3] equilibrium desorption constant
Dm = 5.2154833353911946e-11 # [m2 s-1] Molecular diffusivity
Dc = 0.008 # [m] Internal diameter of the column
Lc = 0.1 # [m] Length of the column
Ac = pi * (Dc/2)^2 # [m2] Cross sectional area of the column
M = 145000 # [g mol-1] Molecular weight of mAb
cin = 2.27586207e-02 # [mol m-3]
Qin = 1.05000000e-08 # [m3 s-1] inlet volumetric flowrate RT: 8 0.63 mL/min
u = Qin/(Ac*Īµ_column) # [m s-1] Liquid interstitial velocity
Re = 997*dp*u/0.001002 # [-] Reynold's number of the flow (using water dynamic viscosity and density)
Dax = dp*u*Īµ_column/(0.339 + 0.033*Re^0.48) # [m2 s-1] Axial dispersion coefficient
Deff = kfit*Īµ_particle*Dm # [m2 s-1] Effective pore diffusion coefficient
k_film = 3/(rp) # [m-1] film mass transfer coefficient
dt = Differential(t)
dz = Differential(z)
dr = Differential(r)
dz_sq = Differential(z)^2
dr_sq = Differential(r)^2
eqs = [
Īµ_column*dt(c(t, z)) ~ -u*Īµ_column* dz(c(t, z)) + Īµ_column*Dax*dz_sq(c(t, z)) - (1 - Īµ_column)*(3/rp)*Deff*dz(cā(t, z, rp)), #### this last term is the problem
Īµ_particle*dt(cā(t, z, r)) ~ Deff*(dr_sq(cā(t, z, r)) + 2/(r)*dr(cā(t, z, r)))- (1 - Īµ_particle)*kd*qmax*c(t, z)/(kd + c(t, z))^2,
]
bcs = [
# boundary conditions
c(t,0) ~ cin + (Dax/u) * dz(c(t,0)),
dz(c(t,zL)) ~ 0,
cā(t, z, rp) ~ c(t, z),
dr(cā(t,z, 0)) ~ 0,
# Initial conditions
c(0, z) ~ 0,
cā(0, z, r) ~ 0,
]
domain = [z ā (0.0, zL), r ā (0.0, rp), t ā (0.0, 400)]
@named HPLC = PDESystem(eqs, bcs, domain, [t, z, r], [c(t, z), cā(t, z, r)])
# prob = discretize(HPLC, MOLFiniteDifference([z=>40, r=>40], t), advection_scheme = WENOScheme())
prob = discretize(HPLC, MOLFiniteDifference([z=>40, r=>40], t), advection_scheme = UpwindScheme())
sol = solve(prob, Tsit5(), saveat = 1)
Describe the bug š
I am trying to model a chromatography column using ModelingToolkit and MethodOfLines, and I can't get past this error message.
From my testing the problematic term is
- (1 - Īµ_column)*(3/rp)*Deff*dz(cā(t, z, rp))- commenting it out discretizes and solves the model. To be clear this should be the the partial derivative evaluated at rp, I'm not sure if this is how it should be coded.I don't understand how to fix it - previous reports of the same error were due to mixed derivatives, but I don't think I have them do I?
To be clear: the numbers in the Any[...] vector belong to the z domain I think, changing zL changes them too.
I'd be very grateful for some help!
Minimal Reproducible Example š
Error & Stacktrace ā ļø
Environment (please complete the following information):
using Pkg; Pkg.status()using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)versioninfo()Additional context
Add any other context about the problem here.