\frac{\partial c}{\partial t} + \frac{1 - \epsilon}{\epsilon} \frac{\partial q}{\partial t} + u \frac{\partial c}{\partial x} = D \frac{\partial^2 c}{\partial x^2}
\frac{\partial q}{\partial t} = - k_f \left( q - \frac{q_m b c}{1 + b c} \right)
Initial conditions:
t = 0, c(x, 0) = q(x, 0) = 0, x > 0
Boundary conditions:
x = 0,
D \frac{\partial c}{\partial x} = u*(c - c0)
x = L,
D \frac{\partial c}{\partial x} = 0
Parameters:
L = 20 # Length of bed ,cm
ϵ = 0.4436 # bed porosity
c0 = 8.8 # Initial concentration ,mg/ml
qm = 76.85 # maximum adsorption capacity ,mg/ml
b = 0.0258 # Langmuir isotherm constant
u = 0.0424 # Velocity ,cm/s
k_f = 0.89 # mass transfer coefficient ,s-1
D = 2.69e-4 # Axial diffusion coefficient ,cm2/s
Discretize the spatial variables using the finite element orthogonal collection method, and divide the column into 30 finite elements, each with 2 configuration points, and solve the pde.
There are some tutorials using Matlab to solve the similar questions:
ERROR: ArgumentError: Upper boundary condition Differential(x)(c(t, 20.0)) ~ 0.0 on time variable is not supported, please use a change of variables t => -τ to make this an initial condition.
"
I think the condition for that error is erroneously including bcs in space that are on the upper end of the domain when they have time derivs, while it is only supposed to be thrown for final conditions in time.
…
Conditions like this one used to work before I changed the way initial conditions were parsed, and added this error message - the underlying discretization hasn’t changed though. This is a 10 line fix, I can get on this once I’m done debugging integrals
"
I’ve replicated this myself and get the same error.
I believe that the error is referring to the transformation τ=T−t, known as the “backward time” or “time reversal” technique, which is a common method used in solving PDEs with final time conditions like c(t=t_final). This is a strange request given that there are no final time conditions, but let’s lay out the details of this technique (as I understand it) incase I am missing something.
This change of variables effectively reverses the time direction, where:
τ is the new time variable
T is the final time
t is the original time variable
This means that:
When t = 0 (start of original time), τ = T
When t = T (end of original time), τ = 0
This converts a condition at the final time, t=T, into an initial condition, τ = 0.
Let’s consider how this transformation affects our PDE and its conditions:
PDE Equation:
The original PDE ∂u/∂t = … becomes -∂u/∂τ = …
(The negative sign appears due to the chain rule: ∂u/∂t = ∂u/∂τ * ∂τ/∂t = -∂u/∂τ)
Initial Condition:
The original final condition at t = T becomes the initial condition at τ = 0
Boundary Conditions:
Spatial boundary conditions remain the same, but are now expressed in terms of τ
This allows the use of standard time-marching schemes (forward in τ).
However, this is clearly not very useful in our case since we are concerned with a boundary condition in space, not time. Furthermore, our initial conditions q(t=0, x)=0 and c(t=0, x)=0 will both become final time conditions w.r.t. τ with this transformation.
Nonetheless, I implemented the change of variables. This time, I expected that the two transformed initial conditions (w.r.t. “t”) to cause the final time error, rather than the spatial boundary condition. To my surprise, I get exactly the same error triggered by the spatial boundary condition at the outlet.
ERROR: ArgumentError: Upper boundary condition Differential(x)(c(τ, 20.0)) ~ 0.0 on time variable is not supported, please use a change of variables `t => -τ` to make this an initial condition.
Stacktrace:
[1] (::MethodOfLines.var"#521#523"{…})(ic::PDEBase.UpperBoundary)
@ MethodOfLines ~/.julia/packages/MethodOfLines/iIYN1/src/discretization/generate_ic_defaults.jl:7
[2] _mapreduce(f::MethodOfLines.var"#521#523"{…}, op::typeof(vcat), ::IndexLinear, A::Vector{…})
@ Base ./reduce.jl:440
[3] _mapreduce_dim(f::Function, op::Function, ::Base._InitialValue, A::Vector{Any}, ::Colon)
@ Base ./reducedim.jl:365
[4] mapreduce
@ ./reducedim.jl:357 [inlined]
[5] generate_ic_defaults(tconds::Vector{…}, s::MethodOfLines.DiscreteSpace{…}, ::MOLFiniteDifference{…})
@ MethodOfLines ~/.julia/packages/MethodOfLines/iIYN1/src/discretization/generate_ic_defaults.jl:5
[6] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
@ PDEBase ~/.julia/packages/PDEBase/nzap9/src/symbolic_discretize.jl:84
[7] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…}; analytic::Nothing, kwargs::@Kwargs{})
@ PDEBase ~/.julia/packages/PDEBase/nzap9/src/discretization_state.jl:58
[8] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
@ PDEBase ~/.julia/packages/PDEBase/nzap9/src/discretization_state.jl:55
[9] top-level scope
@ ~/Dev/github_projects/julia/ModelingToolkit-practice/models/chromatography/change_vars.jl:62
Some type information was truncated. Use `show(err)` to see complete types.
This leads me to suspect @xtalax original hypothesis: “the condition for that error is erroneously including bcs in space that are on the upper end of the domain when they have time derivs, while it is only supposed to be thrown for final conditions in time.”
@Stephen I believe I caught a mistake in your inlet boundary condition, where:
∂x(c(t, 0)) ~ u * (c(t, x) - c0) / D,
should be
∂x(c(t, 0)) ~ u * (c(t, 0) - c0) / D,
but it will still yield the same error of course.
I happen to be modeling chromatography in Julia, too, including more advanced models and would be interested in collaborating. I’m especially interested in making these into MTK model components that can be connected in an equation oriented flowsheet. Send me a message if you’re interested and perhaps we can find a solution faster, together. https://www.linkedin.com/in/samuel-m-andersson/
Here is my refactor of @Stephen 's code for a change of variables incase it is helpful for debugging purposes. It does not include a re-transform back to the original forward time “t” coordinate system.