# Help with setting up coupled PDE system for adsorption column
**URL:** https://discourse.julialang.org/t/help-with-setting-up-coupled-pde-system-for-adsorption-column/130751
**Category:** Numerics
**Created:** [July 15, 2025, 2:28pm UTC](https://discourse.julialang.org/t/help-with-setting-up-coupled-pde-system-for-adsorption-column/130751 "2025-07-15T14:28:52Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@shisac](https://discourse.julialang.org/u/shisac)
#### Post date: [July 15, 2025, 2:28pm UTC](https://discourse.julialang.org/t/help-with-setting-up-coupled-pde-system-for-adsorption-column/130751/1 "2025-07-15T14:28:52Z")
Hi all,
I’m trying to model an adsorption column with a system of PDEs, which have been discretized in space using a WENO finite volume scheme. I’m new to julia and this is my first time simulating a complex system like this. I am currently unsure what the best framework is in which to write out the equations and use solvers from DifferentialEquations.jl. I’ve attempted to write the model both as an `ODEProblem()` and as a `DAEProblem()`, but the solvers I’ve tried either crash by reducing the time step to an extremely small value, or run very slowly.
Specifically the PDEs that I’m solving are (after removing most constants):
- Component Mass balance:
\frac{\partial y\_i}{\partial t} = \frac{T}{P} \frac{\partial}{\partial z} \left(\frac{P}{T} \frac{\partial y\_i}{\partial z} \right) - \frac{T}{P} \frac{\partial x\_i}{\partial t} - \frac{y}{P} \frac{\partial P}{\partial t} + \frac{y}{T} \frac{\partial T}{\partial t}
and (hopefully self-enforced) y\_1 + y\_2 + y\_3 = 1.
- Total Mass balance:
\frac{\partial P}{\partial t} = - T \frac{\partial }{\partial Z} (P v / T) - T \sum\_{i=1}^2 \frac{\partial x\_i}{\partial t} + \frac{P}{T} \frac{\partial T}{\partial t}
- Solid phase balance:
\frac{\partial x\_i}{\partial t} = \alpha\_i (x\_i^\*(T, P) - x\_i)
- Column energy balance:
\frac{\partial T}{\partial t} = \frac{\partial^2 T}{\partial Z^2} - \frac{\partial}{\partial Z} (v P) - T \sum\_{i=1}^2 \frac{\partial x\_i}{\partial t} - \Omega(x\_1, x\_2) \frac{dP}{dt}
- Pressure drop
- \frac{dP}{dZ} = a v + b v |v|
These equations and boundary/initial conditions are taken directly from the paper _Multiobjective Optimization of a Four-Step Adsorption Process for Postcombustion CO2 Capture Via Finite Volume Simulation_.
The only purely algebraic equation is the pressure drop relation, which in theory could be solved separately at each time step. However, because the system is so tightly coupled, I’m not sure whether I should be formulating this as an `ODEProblem` or a `DAEProblem`. I also tried writing the system in mass matrix form, but it seems that isn’t supported when the coefficients in front of the time derivatives depend on the state variables.
The other issue I’m having is that the simulations without the solid phase balance equation (dx\_i/dt = 0) are using up an enormous amount of memory, e.g. `45.900171 seconds (3.93 G allocations: 62.034 GiB, 12.26% gc time, 6.25% compilation time)`. Even though I’m using preallocated buffers for intermediate arrays, this seems excessive. I’m not sure what’s causing such high memory usage.
My project can be found at [GitHub - shivang57721/adsorption\_process: Numerical simulation of adsorption](https://github.com/shivang57721/adsorption_process). I apologize that it’s not a ‘minimal example’, but I’m not really sure how to best isolate the problem. I would greatly appreciate any advice on debugging memory/stability issues, or maybe more general suggestions for organizing coupled PDEs in julia. For example, I see that the package `ModelingToolkit.jl` might be useful, but I’m not really sure.
---
### Author:  [@j-fu](https://discourse.julialang.org/u/j-fu)
#### Post date: [July 15, 2025, 3:07pm UTC](https://discourse.julialang.org/t/help-with-setting-up-coupled-pde-system-for-adsorption-column/130751/2 "2025-07-15T15:07:55Z")
Interesting problem. For the mass matrix formulation you could use an internal energy formulation so that you have \frac{\partial U}{\partial t} instead of \frac{P}{T}\frac{\partial T}{\partial t} and introduce an algebraic equation for U(P,T). In principle then it should be possible to write the full system as in the paper as a DAE system and to use some stiff DAE capable solver.
That said, we tackled a slightly similar problem in this paper: [https://doi.org/10.1016/j.cej.2025.162027](https://doi.org/10.1016/j.cej.2025.162027) with VoronoiFVM.jl (Disclaimer: I am the main author). This could work for the upwind discretization, but probably not for WENO et al.
---
### Author:  [@shisac](https://discourse.julialang.org/u/shisac)
#### Post date: [July 16, 2025, 3:34pm UTC](https://discourse.julialang.org/t/help-with-setting-up-coupled-pde-system-for-adsorption-column/130751/3 "2025-07-16T15:34:43Z")
Thanks a lot for your reply. I’ve looked into VoronoiFVM.jl to see if I can model the system with the package with an upwind scheme for the flux. I also introduced U(P,T) = \log P - \log T so that the total mass balance equation (without any adsorption) simplifies to
\frac{dU}{dt} + \frac{T}{P} \frac{\partial}{\partial Z} \left( \frac{P}{T} v \right) = 0
I was struggling a bit with implementing this equation with VoronoiFVM.jl. The flux term has a coefficient \frac{T}{P}, but all the examples I’ve seen don’t have any coefficients in front of the flux term. Is there a way to encode this in the VoronoiFVM.jl format for the flux ?
Another tangential thing I was unsure about was why the diffusion flux for \Delta u is encoded as u\_r - u\_l in all the examples in the documentation, rather than (u\_r - u\_l) / \Delta x ?
---
### Author:  [@moyner](https://discourse.julialang.org/u/moyner)
#### Post date: [July 16, 2025, 7:20pm UTC](https://discourse.julialang.org/t/help-with-setting-up-coupled-pde-system-for-adsorption-column/130751/4 "2025-07-16T19:20:36Z")
Hi, we have a code for this specific type of system that we will open source after the summer once docs are finished. Feel free to drop me a message if you are interested.
---
### Author:  [@Vinicius\_Viena](https://discourse.julialang.org/u/Vinicius_Viena)
#### Post date: [September 3, 2025, 6:47am UTC](https://discourse.julialang.org/t/help-with-setting-up-coupled-pde-system-for-adsorption-column/130751/5 "2025-09-03T06:47:42Z")
Hi @shisac . I’ve solved PDEs in chromatography but for liquid phase using orthogonal collocation on finite elements with cubic hermite polynomials. [ude\_chromatography/UDE\_paper\_chromatography/PDE\_gradients\_lux\_mechanistic.jl at master · viniviena/ude\_chromatography · GitHub](https://github.com/viniviena/ude_chromatography/blob/master/UDE_paper_chromatography/PDE_gradients_lux_mechanistic.jl).
And it is an old code and would do many different things today. Also, if you need support for isotherms, I also maintain [GitHub - ClapeyronThermo/Langmuir.jl: Single and Multi-component Adsorption Equilibrium.](https://github.com/ClapeyronThermo/Langmuir.jl). I’d be glad to help making this problem work.
---
### Author:  [@Bruno\_De\_Jonckheere](https://discourse.julialang.org/u/Bruno_De_Jonckheere)
#### Post date: [September 17, 2025, 6:42pm UTC](https://discourse.julialang.org/t/help-with-setting-up-coupled-pde-system-for-adsorption-column/130751/6 "2025-09-17T18:42:00Z")
Moyner, I was wondering if you had made your publications.
---
### Author:  [@moyner](https://discourse.julialang.org/u/moyner)
#### Post date: [September 23, 2025, 8:49am UTC](https://discourse.julialang.org/t/help-with-setting-up-coupled-pde-system-for-adsorption-column/130751/7 "2025-09-23T08:49:11Z")
Hi,
The code is available here: [GitHub - sintefmath/Mocca.jl: CO2 capture modules in Julia. Name subject to change.](https://github.com/sintefmath/Mocca.jl)
It is still being made ready for the initial release (docs/examples) but it is possible to use with some care.
---
### Author:  [@Rui\_Matias1](https://discourse.julialang.org/u/Rui_Matias1)
#### Post date: [November 6, 2025, 8:56am UTC](https://discourse.julialang.org/t/help-with-setting-up-coupled-pde-system-for-adsorption-column/130751/8 "2025-11-06T08:56:40Z")
Hello @moyner ,
I’ve been working with your code over the past few days and wanted to ask if you have any idea when an initial “final” version might be released.
I noticed that the current implementation assumes an adiabatic system (`h_in = 0` and `h_out = 0`). When I modify these values to simulate non-adiabatic conditions, the code fails to run. Do you plan to extend the model to handle non-adiabatic systems in the future?
Also, I found a small issue that might be helpful to point out. When changing the initial gas composition in the column, the code still assumes that `q_CO2` is zero — meaning that CO₂ is never initially present in the column, even if specified. I tracked down the cause in the `init_adsorption_column.jl` file:
```julia-auto
qN2 = map(1:ncells) do i
qstar = compute_equilibrium(model.system, c[i,:], temperature_init[i])
qstar[2]
end
qCO2 = zeros(eltype(qN2), ncells)
q_init = hcat(qCO2, qN2)
```
This should instead be:
```julia-auto
qN2 = map(1:ncells) do i
qstar = compute_equilibrium(model.system, c[i,:], temperature_init[i])
qstar[2]
end
qCO2 = map(1:ncells) do i
qstar = compute_equilibrium(model.system, c[i,:], temperature_init[i])
qstar[1]
end
q_init = hcat(qCO2, qN2)
```
Just wanted to give you a heads-up in case you haven’t already updated that part.
---
### Author:  [@moyner](https://discourse.julialang.org/u/moyner)
#### Post date: [November 7, 2025, 2:46pm UTC](https://discourse.julialang.org/t/help-with-setting-up-coupled-pde-system-for-adsorption-column/130751/9 "2025-11-07T14:46:35Z")
Hi! Thanks for the bug report - happy to see that you were able to make use of the code. There are updates that should fix the non-adiabatic conditions as well as the initial condition.
For the full release, we still have a few todos, but we hope to get it out by the end of the year. We are working on more examples and making sure that the setup is as robust as possible for optimization usage.
---
### Author:  [@FrancescaWatson](https://discourse.julialang.org/u/FrancescaWatson)
#### Post date: [February 4, 2026, 7:53am UTC](https://discourse.julialang.org/t/help-with-setting-up-coupled-pde-system-for-adsorption-column/130751/10 "2026-02-04T07:53:19Z")
Hi @Rui_Matias1 we now have the initial release of the Mocca code registered as a package in Julia.
It is still a new code under active development so the interface will almost certainly change but we hope this version is reasonably easy to understand and use for basic simulation and optimisation.
Let us know if you find any issues with it or if you have any other feedback / suggestions 🙂