Animating 1D Convection-Diffusion Equation to reach steady state

Question

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!!

Answer 1

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<x<L_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}]