Let’s say that we have a darcy flow with a coefficient \(D\) in a pipe of length \(L\) and cross-section area \(A\) with a pressure drop \(\Delta p\). The velocity can be calculated from a relation: \[\frac{\Delta p}{L} = Du\] Now let’s have a expotentialy decreasing concentration: \[\dot C = -\lambda C\] The solution of this equation is simpli expressed as \(C(t)=e^{-\lambda t}C(0)\). The length of time spent in the porous pipe is \(\frac{L}{u}\). The concentration lost in the pipe would be equal to the change in concetration \(C(0)-e^{-\lambda t}C(0)\) times the flux \(Q=uA\): \[ AC(0)\left(1-e^{-\lambda \frac{L}{u}}\right)u = AC(0)\left(1-e^{-\frac{\lambda L^2D}{\Delta p}}\right)\frac{\Delta p}{LD}\] It’s clear that \(A\) and \(C(0)\) is inconsequential so we will drop them. Let now: \[ f(L) = \left(1-e^{-\frac{\lambda L^2D}{\Delta p}}\right)\frac{\Delta p}{LD}\] Let’s see how this function look for example variables:

lambda = 1
dp = seq(0,2,len=5)
L = seq(0.1,5,len=300)
D = 1

f = outer(L,dp,function(L,dp)dp/(L*D)*(1-exp(-D*L*L*lambda/dp)))
matplot(L, f, type="l",lty=1)
legend("topright", lty=1, col=1:length(dp),legend = dp,title = "dp")

It’s clear that there is a maximum with a zero derivative: \[ f'(L) = -e^{-\frac{\lambda L^2D}{\Delta p}}\left(-\frac{\lambda 2LD}{\Delta p}\right)\frac{\Delta p}{LD}+\left(1-e^{-\frac{\lambda L^2D}{\Delta p}}\right)\frac{\Delta p}{D}\left(-\frac{1}{L^2}\right)\] \[ f'(L) = 2\lambda e^{-\frac{\lambda L^2D}{\Delta p}}-\left(1-e^{-\frac{\lambda L^2D}{\Delta p}}\right)\frac{\Delta p}{DL^2}\] Let’s set \(w=\frac{\lambda L^2D}{\Delta p}\): \[ f'(L) = 2\lambda e^{-w}-\left(1-e^{-w}\right)\frac{\lambda}{w}\] After equating it to \(0\) we get: \(1+2w - e^w = 0\). The root of this equation can be easily found (one root is \(0\) of course):

w = uniroot(function(w) 1+2*w-exp(w),interval = c(0.0001,10))$root
print(w)
## [1] 1.256432

Now we can calculate the optimal lengths for specific pressure drops:

optL = sqrt(w*dp/(lambda*D))
optf = dp/(optL*D)*(1-exp(-D*optL*optL*lambda/dp))
matplot(L, f, type="l",lty=1)
points(optL,optf)