Thanks very much for your suggestion. I appreciate the time you have taken to analyze the question I posed. I do not fully understand your code but it appears that you do 10 iterations starting from a linear relation between $u$ and $x$. The first plot shows nice convergence for $u$ vs. $x$ as does the second $u$ vs. $t$. I am puzzled as to why the $u$ vs. $t$ plots for the two methods are quite different. Similarly the 3D plots for the two methods look different.

Thank you again. Regarding -- Warning: boundary and initial conditions are inconsistent., Tim Laska uses an exponential ramp to circumvent this issue - see his posts above. The inconsistencies in the two results is a problem - which is correct?

I'm too much an amateur to post a sol'n.

But my book says it's often the case measurements are made (with sensors) and later (if found useful) looked back upon in retrospect with pde. The example given was temperature of curing concrete in a dam - that sensors are used.

Rutton,

At this year's Wolfram Conference, Wolfram developer, Yi-Lin Chiu demonstrated how to handle PDE's for acoustics. While your PDE's will be a bit different, I think this talk will be helpful for generally solving PDE's relating space and time with various boundary conditions in Mathematica.

On this page- search for Modeling Acoustics with Differential Equations

You can use his notebook as a tutorial. I did not compare it to your example, I'd be curious if you find it helpful.

Regards,

Neil

Hi Rutton,

I would say that the answer to your question (1) is probably yes there is a mistake in your formulation at least at this time because your derivative term evaluated at the boundary is infinite at $t=0$. The condition that you imposed would look like UnitStep[-x] at $t=0$.

Mathematica says this is indeterminate.

(D[UnitStep[-x], x] /. x -> 0) (*  Indeterminate *)

When you execute your code, you get the warning

NDSolveValue::delpde: Delay partial differential equations are not currently supported by NDSolve.

Perhaps, that is what you need to get beyond your infinite term at $t=0$, but it is not implement yet. What is interesting, is that since the derivative at the boundary is proportional to molar velocity, and it would be physically impossible to have infinite velocity.

Best regards,

Tim

I will add the following for completeness. If you read the The Numerical Method of Lines Tutorial, it shows you how you can roll your own discretization scheme using NDSolve`FiniteDifferenceDerivative to accomplish what you need.

Use equations 19.1-12 and 19.1-13 from BSL to reduce $\mathbf{v}^*=f(t\ only)$ as shown:

$$\frac{{\partial c}}{{\partial t}} = - \left( {\nabla \cdot c{\mathbf{v}^*}} \right) + \sum\limits_{\alpha = 1}^N {{R_\alpha }} \qquad (19.1-12) $$

$$ \left( {\nabla \cdot {\mathbf{v}^*}} \right) =\frac{1}{c} \sum\limits_{\alpha = 1}^N {{R_\alpha }} \qquad (19.1-13) $$

$$\left( {\nabla \cdot {{\mathbf{v}}^*}} \right) = 0$$

Since the divergence of $\mathbf{v}^*$ is zero and this is a 1-d problem it means this value is a constant with respect to $z$. Equation (19.1-17) simplifies to the boxed form for a non-reactive system.

$$c\left( {\frac{{\partial {x_\alpha }}}{{\partial t}} + \left( {{{\mathbf{v}}^*} \cdot \nabla {x_\alpha }} \right)} \right) = c{\mathcal{D}_{AB}}{\nabla ^2}{x_A} + \left( {{x_B}{R_A} - {x_A}{R_B}} \right) \qquad (19.1-17) $$

$$\boxed{\left( {\frac{{\partial {x_\alpha }}}{{\partial t}} + \left( {{{\mathbf{v}}^*} \cdot \nabla {x_\alpha }} \right)} \right) = {\mathcal{D}_{AB}}{\nabla ^2}{x_A}}$$

Now, we need to solve for $\mathbf{v}^*$, which is defined by

$$\mathbf{v}^*=\frac{N_A+N_B}{c}$$

Since $\frac{N_A+N_B}{c}$ does not depend on $z$, we choose the interface condition where $N_B$ equals 0 similar to the situation of steady-state and we obtain the same relation for $N_A$.

$$N_{A z}=-\frac{c \mathcal{D}_{A B}}{1-x_A}\frac{\partial x_A}{\partial z}$$

By rolling our own scheme, we have access to the derivative at the boundary. The following shows the Mathematica implementation and evaluation of the boxed equation evaluated up to $\frac{1}{4}$ length.

molMod[x_, ramp_, tmax_] := 
 Module[{n, h , sf, U, bc0, bc0d , bc1, bc1d, delconc, laplacianconc, 
   vstar, eqns, ics, system, lines, f, xsat, m},
  n = 50; Subscript[h, n] = 1/n; sf = 1;
  xsat = Rationalize[x];
  m = Rationalize[ramp];
  U[t_] = Table[Subscript[u, i][t], {i, 0, n}];
  bc0 = Subscript[u, 0][t] == xsat (1 - Exp[-m t]);
  bc0d = Map[D[#, t] + sf # &, bc0];
  bc1 = Subscript[u, n][t] == 0;
  bc1d = Map[D[#, t] + sf # &, bc1];
  delconc = 
   NDSolve`FiniteDifferenceDerivative[1, Subscript[h, n] Range[0, n], 
    U[t]];
  laplacianconc = 
   NDSolve`FiniteDifferenceDerivative[2, Subscript[h, n] Range[0, n], 
    U[t]];
  vstar = -1/(1 - Subscript[u, 0][t]) delconc[[1]];
  eqns = Thread[D[U[t], t] + vstar delconc -  laplacianconc == 0];
  eqns[[1]] = bc0d;
  eqns[[-1]] = bc1d;
  ics = Thread[U[0] == Table[0, {n + 1}]];
  system = Join[eqns, ics];
  lines = NDSolve[system, U[t], {t, 0, 1}];
  f = Interpolation[
    Flatten[Table[{{t, i Subscript[h, n]}, 
       First[Subscript[u, i][t] /. lines]}, {i, 0, n}, {t, 0, tmax, 
       tmax/500}], 1], InterpolationOrder -> 2];
  f]
f75 = molMod[0.75, 20000, 0.01];
f50 = molMod[0.5, 20000, 0.01];
f25 = molMod[0.25, 20000, 0.01];
inv[f_, s_] := Function[{t}, s /. FindRoot[f - t, {s, 1}]]
xa0[phi_] := 1/(1 + 1/(Sqrt[Pi] (1 + Erf[phi]) phi Exp[phi^2]))
phiinv = inv[xa0[phi], phi];
bigX[Z_, x_] := (1 - Erf[Z - phiinv[x]])/(1 + Erf[phiinv[x]])
frame[legend_] := 
 Framed[legend, FrameStyle -> Red, RoundingRadius -> 10, 
  FrameMargins -> 0, Background -> Lighter[LightYellow, 0.25]]
legend = LineLegend[
   Automatic, {"3/4 \[Infinity]", "3/4", "1/2 \[Infinity]", "1/2", 
    "1/4 \[Infinity]", "1/4", "0 \[Infinity]"}, 
   LegendFunction -> frame];
Plot[{bigX[z, 0.75], 
  Callout[f75[1/256, 2 Sqrt[1/256] z]/0.75, "3/4", 1, 
   CalloutMarker -> "CirclePoint"], bigX[z, 0.5], 
  Callout[f50[1/256, 2 Sqrt[1/256] z]/0.5, "1/2", 0.9, 
   CalloutMarker -> "CirclePoint"], bigX[z, 0.25], 
  Callout[f25[1/256, 2 Sqrt[1/256] z]/0.25, "1/4", 0.8, 
   CalloutMarker -> "CirclePoint"], bigX[z, 0]}, {z, 0, 2}, 
 PlotRange -> {0, 1}, GridLines -> Automatic, 
 Filling -> {{1 -> {2}}, {3 -> {4}}, {5 -> {6}}}, 
 PlotStyle -> {Dashed, Thick, Dashed, Thick, Dashed, Thick, Dashed}, 
 ImageSize -> Large, Frame -> True, 
 FrameLabel -> {Style["Z", 18], Style["X", 18], 
   Style[StringForm["Dimensionless Concentration Profiles"], 18]}, 
 PlotLegends -> Placed[legend, {Right, Center}]]

Very tight agreement with the semi-infinite case. The following plot shows how different the new custom formulation (dashed lines) is versus my initial quick and dirty model (solid lines).

There looks to be about a 10% error that exists at dimensionless time scales between 0.05-0.15. In the end, one needs to operate within the constraints of the currently available solver. Mathematica 11.3 has the capability to solve this case, but you probably will need to do a deep dive into the MOL tutorial and start with a different form of the balance equation.

By rolling your own Method Of Lines (MOL), you get more flexibility on what you can do with the derivatives. If you look at "vstar", then you see it was only evaluated at the interface $u_0$ (note that part 1 of delconc is also at the interface).

vstar = -1/(1 - Subscript[u, 0][t]) delconc[[1]];

So, vstar is evaluated at the interface only where the flux of B is zero and it is only a function of $t$ because its divergence is zero. Therefore, it should satisfy your conditions of only being evaluated at the interface and not throughout the domain.

Regarding the figure 20.1-1, if the characteristic length scale is, $L$, and the characteristic time scale is $\frac{L^2}{D_{AB}}$, then Z using dimensionless units becomes

$$Z=\frac{z}{\sqrt{4t}}$$

Given that figure 20.1-1 has a Z value that goes from 0 to 2 and we know that z is between 0 and 1, we can solve for t in terms of z.

Solve[z/(2*Sqrt[t]) == 2, t]
(* {{t -> z^2/16}} *)

So, I chose $z=\frac{1}{4}<1$ and that meant the time when $Z=2$ was $t=\frac{1}{256}$. Then, I solved for $z$ in terms of $Z$ using.

Solve[bigZ == z/(2*Sqrt[t]), z]
(* {{z -> 2*bigZ*Sqrt[t]}} *)

Taken together, that means the following is the numerically derived $X$ for $x_{A0}=0.5$ @ $z=\frac{1}{4}$ (note that z in the Mathematica implementation is actually $Z$ to avoid beginning with capitals).

f50[1/256, 2 Sqrt[1/256] z]/0.5

If we want to look at the case where $z=1$, then $t=\frac{1}{16}$. Now, instead of following the infinite tube, we expect the numerical result to satisfy the $x_{A0}=0$ boundary condition.

Plot[{bigX[z, 0.75], 
  Callout[f75[1/16, 2 Sqrt[1/16] z]/0.75, "3/4", 1, 
   CalloutMarker -> "CirclePoint"], bigX[z, 0.5], 
  Callout[f50[1/16, 2 Sqrt[1/16] z]/0.5, "1/2", 0.9, 
   CalloutMarker -> "CirclePoint"], bigX[z, 0.25], 
  Callout[f25[1/16, 2 Sqrt[1/16] z]/0.25, "1/4", 0.8, 
   CalloutMarker -> "CirclePoint"], bigX[z, 0]}, {z, 0, 2}, 
 PlotRange -> {0, 1}, GridLines -> Automatic, 
 Filling -> {{1 -> {2}}, {3 -> {4}}, {5 -> {6}}}, 
 PlotStyle -> {Dashed, Thick, Dashed, Thick, Dashed, Thick, Dashed}, 
 ImageSize -> Large, Frame -> True, 
 FrameLabel -> {Style["Z", 18], Style["X", 18], 
   Style[StringForm["Dimensionless Concentration Profiles"], 18]}, 
 PlotLegends -> Placed[legend, {Right, Center}]]

It appears that the BC is satisfied.

Now, we can test your assertion that my model does not allow for mass conservation by simply evaluating $N_{Bz}$ at $z=1$. You can use the relations stated previously to derive an expression for the molar flux of B. Start with the dimensionless molar velocity.

$$\mathbf{v}^*=N_A+N_B$$

Use the expression for $N_{B z}$ in terms of $x_A$. $$N_{B z} = \left (1-x_A \right ) \left(N_{A z}+N_{B z}\right)+\frac{\partial x_A}{\partial z} $$

Substitute $\mathbf{v}^*$ for $N_A+N_B$ $$N_{B z} = \left (1-x_A \right ) \mathbf{v}^*+\frac{\partial x_A}{\partial z}$$

Evaluate at the exit ( $z=1$).

$${N_{Bz}} = {\left. {\left( {1 - {x_A}} \right){{\mathbf{v}}^*} + \frac{{\partial {x_A}}}{{\partial z}}} \right|_{z = 1}}$$

Evaluate $\mathbf{v}^*$ at the interface $${N_{Bz}} = {\left. {\left( {\left( {1 - {x_A}} \right)\left( {{{\left. {\frac{1}{{1 - {x_{A,sat}}}}\frac{{\partial {x_A}}}{{\partial z}}} \right|}_{z = 0}}} \right) + \frac{{\partial {x_A}}}{{\partial z}}} \right)} \right|_{z = 1}}$$

We can implement and plot the last equation in Mathematica code for the $x_{A,sat}=0.5$ case.

Plot[((1 - f50[t, 1]) (-D[f50[t, x], x]/(1 - 0.5)) /. {t -> tp, 
     x -> 0}) + ((D[f50[t, x], x]) /. {t -> tp, x -> 1}), {tp, 0.0002,
   0.1}, PlotRange -> {-1, 40}, ImageSize -> Large, Frame -> True, 
 FrameLabel -> {Style["Dimensionless Time", 18], 
   Style["\!\(\*SubscriptBox[\(N\), \(Bz\)]\)", 18], 
   Style[StringForm["Flux of B at z=1"], 18]}]

Now, we did not achieve infinite flux because this is a numerical model, but the flux of B at the exit starts out very high and asymptotically approaches 0 at large t. So, the assertion of the model having zero B flux along the axis is false.

Thank you Rutton!

Since $Z$ is a combined variable of z and t, I think you need to multiply the $\phi$ expression by $\sqrt{t}$ to make it equivalent to dimensionless $\mathbf{v}^*$. If we choose the $t=\frac{1}{256}$ again, then we can plot the numerically derived result versus the analytical solution using:

Plot[phiinv[x], {x, 0.01, 0.99}, ImageSize -> Large, Frame -> True, 
 FrameLabel -> {Style["xa0", 18], Style["\[Phi]", 18], 
   Style[StringForm["Dimensionless Interface Flux"], 18]}, 
 Epilog -> {PointSize[0.02], 
   Point[{0.25, ((-Sqrt[t] D[f25[t, x], x])/(1 - 0.25))}] /. {t -> 
      1/256, x -> 0}, 
   Point[{0.5, ((-Sqrt[t] D[f50[t, x], x])/(1 - 0.5))}] /. {t -> 
      1/256, x -> 0}, 
   Point[{0.75, ((-Sqrt[t] D[f75[t, x], x])/(1 - 0.75))}] /. {t -> 
      1/256, x -> 0}}]

You can see very tight agreement between the analytical our numerical solution for the evaluated points.

I have both editions, but it has been a couple or more decades since I have read those sections. I will revisit it when I have a few moments.

Personally, I have not seen cases where there is something that looks like a boundary condition expressed in the main equation in NDSolve. That was some cleverness that BSL applied in the 60s before the advent of Computer Algebra Systems and probably is not appropriate for a general purpose PDE solver like NDSolve.

You may be able to express your problem as a system of PDEs $x_A, N_A, N_B$ since you know initial and boundary concentrations and you know that $N_B$ equals zero at $z=0$.

I will take another look when I have a spare moment if nobody else chimes in.