Gas Diffusion PDE

Question

Any help on how to get this started would be awesome, I have no idea where to being with this:

"Approximate the sensor output for the following model

$$u_t+k\cdot u_{xx}=0$$ $$u(x,0)=0$$ $$\pi\cdot a^2\cdot k\cdot u_x(0,t) = -r(t)$$ $$\pi\cdot a^2\cdot k\cdot u_x(l,t) =-E\cdot R\cdot u(l,t)$$ where: $k = 0.5*10^{-4}$; $a=4.0$; $l=25$; $E=0.9$; $R=1.5$; $r=1/11.2$. The desired output st (the sensor time trace) is the concentration at $x = 22$ in mg/L measured in minutes, which is given by: $$st(\tau)=21.3\cdot u(\rho,60\cdot \tau)$$ where $\tau$ is measured in minutes, $\rho=22$, and $c=20$ is the unit conversion factor $$2.0*10^4 mg/mol *10^{-3} m^3/L$$

Basically I am trying to find the gas concentration of some gas, 22m down a 25m tunnel. A sensor at 22m records the concentration once every minute. The gas is absorbed in a filter at the far end of the tunnel, at 90% of its concentration in the 1.5 m$^3$/s of air passing through the filter. I need to simulate the sensor output as a function of time from the instant the gas is released until the sensor reads 2 mg/L.

Answer 1

There is indeed a sign problem. Such an equation is very ill-posed. Forgetting the boundary condition and thinking of this equation in all $\mathbb R$, you can see that it is "intrinsecally" ill-posed by taking Fourier transform for instance : formally, you would have $\hat u(t,\xi) = Ce^{k\xi^2t}$. The function $e^{k\xi^2}$ is growing way too fast as $|\xi| \to +\infty$, so this does not make any sense, even in the Schwartz distribution sense : there are no solution, even in the Schwartz space.

If you change the sign of $k$ it is way better : you get the classic Gaussian function which makes a lot more sense.

Perhaps you have $k > 0$ but $u_x$ and not $u_{xx}$, this would model transport of gas instead of diffusion.

Answer 2

As others have said, $u_t + k u_{xx} = 0$ and $k > 0$ is not well-posed. If instead you consider $u_t - k u_{xx} = 0$ and $k > 0$, the problem is often well-posed (this still depends on the domain, the boundary conditions, etc.). For your particular problem (with domain $\Omega = [0, L]$ and mixed Dirichlet-Robin boundary conditions), you can prove well-posedness using Fourier Series (for a brief sketch of the proof, see this PDF).

The first step involved in numerically solving your problem is to discretize: Let us consider a finite difference discretization (but there are other discretizations, e.g. finite element and method of lines). For accuracy it would be best to use an implicit time stepping scheme (although implementing an explicit time stepping scheme might be considered easier).

The first implicit time-stepping schemes studied in a numerical analysis class are often Backwards Euler and Crank-Nicolson. The error due to the spatial discretization is $O(h)$ for Backwards Euler and $O(h^2)$ for Crank-Nicolson -- this means that as you decrease the step size $h$, you can expect your numerical solution to converge faster to the (unknown) true solution faster with Crank-Nicolson than with Backwards Euler. This does not automatically mean that you should use Crank-Nicolson. For some problem it is best to use one or the other... this requires a bit of background in numerical analysis to explain.

You also need to discretize the boundary conditions. The first one is fairly easy: if your grid of points is $x_0 = 0, x_1, \ldots, x_{N - 1}, x_N = L$, and if $u_0 \approx u(x_0) = u(0)$ then set up your problem so that $u_0 = 0$. The second boundary condition is not so easy. You can approximate the first derivative $u_x(l, t)$ using a "ghost node" or using a higher-order one-sided approximation of $u_x(l, t)$. Using a ghost node will typically cause your error to be $O(h)$ (even if you are using Crank-Nicolson), but it preserves the tridiagonal structure of the matrix involved in the discrete problem, which is nice in terms of solving the linear system of equations. If you are using Crank-Nicolson and want to preserve the $O(h^2)$ order of accuracy, then you should use a higher-order one-sided approximation, which means that your matrix will no longer be tridiagonal (on modern computers this likely will not be an issue). See this PDF for some of the details of finite difference approximation of elliptic problems, higher-order one-sided approximations of $\frac{\partial u}{\partial x}$, and Robin boundary conditions. Your problem is a parabolic problem, but some of the details are similar.