Solving a challenging diffusion problem

143 Views Asked by At

I have the given problem:

A very long pipe is understood simplistically as semi-infinite ( $0<x<\infty$ ). The tube initially contains clean water and is closed. It is lowered into a container with polluted water. At $t=0$ it is opened at the end at $x=0$. Contaminated water then diffuses into the pipe. The concentration of the polluted substance at $x=0$ is then assumed to be $q_0$. a) Set up a mathematical model for the concentration q(x,t) in the tube of the substance as a function of position and time. b) Solve this math problem.

This problem is clearly the diffusion equation. The answer to a) is:

\begin{equation} \alpha u_{xx}=u_t \ \ \ \ \ \ 0\leq x<\infty \\ u(0,t)=q_0\\ u(x,0)=0 \end{equation}

The answer to b), can be solved by using convolution formula:

\begin{equation} u(x,t)=\int_{-\infty}^\infty S(x-y,t)\phi(y)dy=\frac{1}{\sqrt{4\alpha\pi t}}\int_{-\infty}^\infty e^{\frac{-(x-y)^2}{4\alpha t}}\phi(y)dy \end{equation}

Since the convolution formula is for the infinite domain, we define $u(x)$ on an infinite domain using the reflection method:

\begin{equation} \phi(y)=\begin{cases} u(x)\ \ \ \ \ x>0 \\ u(-x) \ \ \ \ x<0 \end{cases} \end{equation}

Then we evaluate $u(x,t)=\phi(x,t)_{-\infty}^\infty$ and modify accordingly to include only the semi-finite domain. However, this gives:

\begin{equation} u(x,t)=\frac{1}{2\sqrt{\pi t}}\int_{-\infty}^\infty e^{-(x-y)^2/4t}q_0dy \end{equation}

Trying to solve this using the substitution $p=\frac{-(x-y)}{2\sqrt{t}}$, $dy=-2\sqrt{t}dp$ we get:

\begin{equation} u(x,t)=-\frac{q_02\sqrt{t}}{2\sqrt{\pi t}}\int_{-\infty}^\infty e^{p^2}dp\\ u(x,t)=-\frac{q_0}{\sqrt{\pi}}\int_{-\infty}^\infty e^{p^2}dp\\ u(x,t)=-\frac{q_0}{\sqrt{\pi}}\bigg[\int_{-\infty}^0 e^{p^2}dp+\int_{0}^\infty e^{p^2}dp\bigg]\\\\ \end{equation}

Since the original function is defined only on the latter domain, we can omit the integral to the left, and get:

\begin{equation} u(x,t)=-\frac{q_0}{\sqrt{\pi}}\bigg[\int_{0}^\infty e^{p^2}dp\bigg]\\ u(x,t)=-\frac{q_0}{2} \end{equation}

But this is clearly not right.

Is there a way to use this approach and get the correct answer?

Thanks

UPDATE:

Following K.defaoite's suggestion we use:

\begin{equation} \int_{0}^\infty \Phi(x-\xi,t)f(\xi)\mathrm d\xi+\int_0^t\int_{0}^\infty\Phi(x-\xi,t-\tau)h(\xi,\tau)\mathrm d\xi\mathrm d\tau \tag{3} \end{equation}

which gives with $f(\xi)=q_0$ and with $p=\frac{x-\xi}{2\sqrt{t}}$:

\begin{equation} \int_{0}^\infty \Phi(x-\xi,t)f(\xi)\mathrm d\xi=\int_0^\infty e^{-p^2}dp=\frac{\sqrt{\pi}}{2} \end{equation}

Then for the second integral, we set $p=\frac{x-\xi}{\sqrt{4(t-\tau)}}$, and still $f(\xi)=q_0$:

\begin{equation} \int_0^t\int_{0}^\infty\Phi(x-\xi,t-\tau)h(\xi,\tau)\mathrm d\xi\mathrm d\tau =-\int_0^t \frac{\sqrt{4(t-\tau)}}{\sqrt{4\pi t}} \int_{0}^\infty e^{-p^2}dpd\tau \end{equation}

since $\int_{0}^\infty e^{-p^2}dp=\frac{\sqrt{\pi}}{2}$, we get:

\begin{equation} -\frac{\sqrt{\pi}}{2}\int_0^t \frac{\sqrt{4(t-\tau)}}{\sqrt{4\pi t}} d\tau=-\frac{t}{3} \end{equation}

Summing up with the first integral above, we get:

\begin{equation} u(x,t)=\frac{\sqrt{\pi}}{2}-\frac{t}{3} \end{equation} However, this is also wrong, as it can be inserted in the PDE and gives an incorrect answer.

Where is the error here?

Thanks

2

There are 2 best solutions below

4
On BEST ANSWER

Full solution.

Referring to my other answer, we let $p(t)=q_0$ and let $u(x,0)=\phi(x)=0$. We then let $$U(x,t)=u(x,t)-p(t)=u(x,t)-q_0$$ Which, you can check, satisfies $$\begin{cases}\partial_t U-k\partial_x^2U=0 \\ U(0,t)=0 \\ U(x,0)=-q_0\equiv \varphi(x)\end{cases}$$ We now extend this to the whole real line $x\in (-\infty,\infty)$ by constructing the odd extension of $\varphi(x)$: $$\varphi_{\text{ext}}(x)=\begin{cases}\varphi(x)=-q_0 & x\geq 0 \\ -\varphi(-x)=q_0 & x<0\end{cases}$$ We then solve $$\partial_tU-k\partial_x^2U=0 ~~~~~x\in\Bbb R,t\in[0,\infty) \\ U(x,0)=\varphi_{\text{ext}}(x)$$ We convolve with the heat kernel: $$U(x,t)=\int_{\Bbb R}\frac{\exp\left(\frac{-(x-\xi)^2}{4kt}\right)}{\sqrt{4\pi kt}}\varphi_{\text{ext}}(\xi)\mathrm d\xi$$ Finally, $$u(x,t)=U(x,t)+q_0=q_0+\int_{\Bbb R}\frac{\exp\left(\frac{-(x-\xi)^2}{4kt}\right)}{\sqrt{4\pi kt}}\varphi_{\text{ext}}(\xi)\mathrm d\xi$$

Shown below is a few snapshots of the solution at $t=10^{-1},10^0,10^1,10^2$, taking $k=q_0=1$, plotted in $x\in[0,10]$.

0.1

1

10

100

As we expect the heat traverses down the pipe.

In terms of constructing the solution analytically, this boils down finding a closed form for the integral $$\int_{\Bbb R}\Phi(x-\xi,t)\varphi_{\text{ext}}(\xi)\mathrm d\xi =\int_{\Bbb R}\Phi(\xi,t)\varphi_{\text{ext}}(x-\xi)\mathrm d\xi \\ $$ We can write this as $$\int_{-\infty}^x \Phi(\xi,t)\varphi_{\text{ext}}(x-\xi)\mathrm d\xi + \int_{x}^\infty \Phi(\xi,t)\varphi_{\text{ext}}(x-\xi)\mathrm d\xi \\ =-q_0\int_{-\infty}^x \Phi(\xi,t)\mathrm d\xi +q_0\int_{x}^\infty \Phi(\xi,t)\mathrm d\xi \\ =q_0\left(\int_x^\infty \Phi(\xi,t)\mathrm d\xi -\int_{-\infty}^x \Phi(\xi,t)\mathrm d\xi\right) \\ =-q_0\operatorname{erf}\left(\frac{x}{\sqrt{4kt}}\right)$$ Hence $$\boxed{u(x,t)=q_0\left(1-\operatorname{erf}\left(\frac{x}{\sqrt{4kt}}\right)\right)}$$

If you want to play with this solution you can do so here.

1
On

Here's a step-by-step approach to solving the heat (or diffusion, whatever you want to call it) equation on the half line. Consider the Dirichlet problem

$$\begin{cases}\partial_tu-k\partial_x^2 u=0 & x\in[0,\infty)~,~t\in[0,\infty) \\ u(0,t)=p(t) \\ u(x,0)=\phi(x)\end{cases}\tag{1}$$ Let $U(x,t)=u(x,t)-p(t)$. We find that $U$ satisfies the forced heat equation with homogeneous Dirichlet BCs: $$\begin{cases}\partial_tU-k \partial_x^2 U=-p'(t) & x\in\Bbb R~,~t\in[0,\infty) \\ U(0,t)=0\\ U(x,0)=\phi(x)-p(0)\equiv \varphi(x) \end{cases}\tag{2}$$

Now, since we have a homogeneous Dirichlet condition at the origin, we can make use of the odd extension of $\varphi$: $$\varphi_{\text{ext}}(x):=\begin{cases}\varphi(x)&x\geq 0 \\ -\varphi(-x) & x<0\end{cases}$$


Now recall the forced heat equation on the whole line: $$\begin{cases}\partial_tu-k\partial_x^2u=h(x,t)& t\in[0,\infty)~,~x\in\Bbb R \\ u(x,0)=f(x) \\ \int_{\Bbb R}u(x,t)^2 \mathrm dx<\infty& \forall t \end{cases}$$ $$u(x,t)=\int_{\Bbb R} \Phi(x-\xi,t)f(\xi)\mathrm d\xi+\int_0^t\int_{\Bbb R}\Phi(x-\xi,t-\tau)h(\xi,\tau)\mathrm d\xi\mathrm d\tau \tag{3}$$

(Source: Introduction to PDEs by Peter J Olver.)

Where $$\Phi(x,t)=\frac{1}{\sqrt{4\pi k t}}\exp\left(\frac{-x^2}{4kt}\right)$$


So, use formula $(3)$ to solve problem $(2)$, taking $h(x,t)=-p'(t)$ and $f(x)=\varphi_{\text{ext}}(x)$, to construct $U(x,t)$ and then simply take $u(x,t)=U(x,t)+p(t)$.