Animating 1D Convection-Diffusion Equation to reach steady state

I'm trying to understand how the parameters ($$c,D$$) of the following equation:

$$\frac{{\partial x}}{{\partial t}} = D\frac{{{\partial ^2}x}}{{\partial {z^2}}} + c\frac{{\partial x}}{{\partial z}}$$

Affect the time it takes to get the system to steady-state.

The boundary conditions are fixed: $$x(0,t) = 17\,\,\,,\,\,\,x(6000,t) = 1$$ $$0 \le x \le 6000$$

I'm new to Mathematica, but it seems that a good strategy will be to use Manipulate to numerically solve the equation with sliders for the timestep, $$c$$ and $$D$$ values. I started with:

pde = D[x[z, t], t] == d*D[x[z, t], {z, 2}] + c*D[x[z, t], z]

bcs = {x[0, t] == 17, x[6000, t] == 1}

But I'm not sure how to continue and which solver should I use for this case. Or maybe there is a better way to approach this problem?

Thank you!!

You need I.C., so made one for you. Feel free to change it. I also changed your notation to make it more common. The dependent variable is $$u(x,t)$$ and the space variable is $$x$$.

In both solutions below, this initial conditions is used

To make it consistent with boundary conditions. Feel free to change it.

Solve for $$u(x,t)$$ with $$t>0, 0 and $$L_0=6000$$

$$\frac{{\partial u}}{{\partial t}} = d\frac{{{\partial ^2}u}}{{\partial {x^2}}} + c \frac{{\partial u}}{{\partial x}}$$

The boundary conditions are

\begin{align*} x(0,t) &= 17\\ x(L_0,t) &= 1 \end{align*}

And IC
$$u(x,0)=-\frac{171 x^3}{2 L_0^3}+\frac{162 x^2}{L_0^2}-\frac{185 x}{2 L_0}+17$$

This below all using V 12.1 on windows 10

NDSolve based solution

Manipulate[
Module[{d, c, solN, pars, L0, pde, ic, bc, x, u, t},
L0 = 6000;
pde = D[u[x, t], t] == d D[u[x, t], {x, 2}] + c D[u[x, t], x];
bc = {u[0, t] == 17, u[L0, t] == 1};
ic = u[x, 0] == -(171/(2 L0^3)) x^3 + 162/L0^2 x^2 - 185/(2 L0) x+ 17;(*made up IC*)
pars = {d -> d0, c -> c0};
solN = Quiet@NDSolve[Evaluate[{pde, ic, bc} /. pars],
u, {x, 0, L0}, {t, 0, t0}];

Quiet@Plot[Evaluate[u[x, t0] /. solN], {x, 0, L0},
PlotRange -> {Automatic, {-10, 17}},
GridLines -> Automatic, GridLinesStyle -> LightGray,
PlotStyle -> Red,
AxesLabel -> {"x", "u(x,t)"}, BaseStyle -> 12]
],
{{d0, 50, "D"}, 50, 5000, 10, Appearance -> "Labeled"},
{{c0, 0, "c"}, 0, 10, 0.1, Appearance -> "Labeled"},
{{t0, 0, "time"}, 0, maxTime, 0.1, Appearance -> "Labeled"},
{{maxTime, 200}, None},
TrackedSymbols :> {d0, c0, t0}]


DSolve based solution

I think DSolve solution could be wrong. When c is non zero, analytical solution no longer agrees the numerical solution. So I would use the numerical solution above for now.

I need to solve this by hand to compare with Mathematica solution to see why this happens.

ClearAll[u, x, t, d, c, d0, c0, n, L0];
pde = D[u[x, t], t] == d*D[u[x, t], {x, 2}] + c*D[u[x, t], x];
bc = {u[0, t] == 17, u[L0, t] == 1};
ic = u[x, 0] == -(171/(2 L0^3)) x^3 + 162/L0^2 x^2 - 185/(2 L0) x + 17;(*made up IC*)
sol = u[x, t] /. First@DSolve[{pde, ic, bc}, u[x, t], {x, t},
Assumptions -> {L0 > 0, 0 < x < L0, t > 0}];
sol = sol /. K[1] -> n


sol = sol /. Infinity -> 20;
sol = Activate[sol]

Manipulate[
L0 = 6000;
Quiet@Plot[Evaluate[sol /. {t -> t0, d -> d0, c -> c0}], {x, 0, L0},
PlotRange -> {Automatic, {-10, 17}},
GridLines -> Automatic, GridLinesStyle -> LightGray,
PlotStyle -> Red,
PlotLabel -> "Analytical solution",
AxesLabel -> {"x", "u(x,t)"}, BaseStyle -> 12
],
{{d0, 50, "D"}, 50, 5000, 10, Appearance -> "Labeled"},
{{c0, 0, "c"}, 0, 10, 0.1, Appearance -> "Labeled"},
{{t0, 0, "time"}, 0, maxTime, 0.01, Appearance -> "Labeled"},
{{maxTime, 200}, None}, TrackedSymbols :> {d0, c0, t0}]


• Thank you so much !!!! It's really helpful!! Can you shortly explain what these two lines are doing sol = u[x, t] /. First@DSolve[{pde, ic, bc}, u[x, t], {x, t}]; sol = Activate[sol /. Infinity -> 30];? – ValientProcess Mar 25 '20 at 1:16
• @ValientProcess The analytical solution is in terms of a series (Fourier). So the sum of terms is up to infinity. ofcourse, we can't sum to infinity, so I choose first 30 terms. Less terms will also work. More terms is more accurate but you really do not need more than few terms in these sums. The sum is returned as inactive so we have to activate it to use it. – Nasser Mar 25 '20 at 1:36
• Thank you !! Also when I tried to run it I get the error: ReplaceAll::rmix: Elements of {(u^(0,1))[x,t]==c (u^(1,0))[x,t]+d (u^(2,0))[x,t],u[x,0]==1,{u[0,t]==17,u[6000,t]==1}} are a mixture of lists and nonlists. and I don't see the plot. Is it a version error? I have version 12 – ValientProcess Mar 25 '20 at 1:45
• @ValientProcess You must have something else going on. Make sure you start with clean kernel. Initially there will be some delay to finish the DSolve. That is all. Using less terms speeds it up. I am using V 12.1. Make sure DSolve works on this in earlier version, May be it does not. I only tried it in V12.1. I can add a NDSolve version later tonight, but I have to go get some food first. – Nasser Mar 25 '20 at 1:48
• Thanks! I tried to restart the kernel and program and ran the code line by line. It seems that the line sol = u[x, t] /. First@DSolve[{pde, ic, bc}, u[x, t], {x, t}] produces the ReplaceAll error, and when I use manipulate I don't see any plot on the grid. – ValientProcess Mar 25 '20 at 1:57