Problem with boundary values for a partial differential equation

Question

I have a PDE problem that I can't solve. and I wonder if you could help me with it.

I have a Terzaghi Consolidation equation (similar to the heat equation) that is written as follows:

$\frac{du}{dt}=c\frac{d^2u}{dt^2}$

$u$ is pore water pressure

$z$ is depth

$t$ is time

$c$ is a constant

pore water pressure is equal at all depth at $t = 0$, but it becomes $0$ at top and bottom directly after that.

u = u[z,t] ($z$ is depth and $t$ is time)

with boundary conditions as follows:

u[z, 0] = 100 ($u$ at all depth when $t = 0$)

u[2, t] = 0 ($u$ at top of the layer; in this case, $2$ is the height of the top layer)

u[0, t] = 0 ($u$ at bottom of the layer)

I have solved this kind of problem before using a finite difference technique in a spread sheet, I wanted to find out whether Mathematica can solve this problem directly, but unfortunately I failed.

This is what I entered into Mathematica:

NDSolve[{Derivative[0, 1][u][z ,t] == Derivative[2, 0][u][z, t],
 u[0, t > 0] == 0, u[2,t > 0] == 0, u[2 > z > 0, 0] == 100}, u, {z, 0, 2},{t, 0, 100}]

Unfortunately, evaluating this expression leads gives the error:

DSolve::litarg:

"To avoid possible ambiguity, the arguments of the dependent variable in u[0,t>0]==0; should literally match the independent variables"

I dont understand why I'm getting this message.

I found a heat problem similar to this. I tried to use the that equation but changed the boundary problem, but that did not work.

Can you help me detect which part of equation is wrong, and how I can fix the boundary values?

Answer 1

You can get rid of the warnings if you modify your initial condition close to the space boundaries :

limit = 3;
newIC[z_] = 100 (HeavisideTheta[z - 10^-limit] - HeavisideTheta[z - (2 - 10^-limit)])

Remove[sol]
sol[z_, t_] = NDSolve[{Derivative[0, 1][u][z, t] == Derivative[2, 0][u][z, t], 
  u[0, t] == 0, u[2, t] == 0, u[z, 0] == newIC[z]}, 
  u[z, t], {z, 0, 2}, {t, 0, 100}][[1, 1, 2]]

Plot[{newIC[z], sol[z, 0]}, {z, 0, 2}, PlotRange -> All, Exclusions -> None]

Answer 2

I need to look more into this, the warning about

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent.

As I've seen this many times. Are you sure the IC's you are using are correct? (ie. consistent) ? Please double check the textbook or from the source you obtained this.

But just wanted to say that you still get a solution and can animate it just fine. A different method solver might be needed like MethodOfLines. But need to look more into this. Here is the solution

ClearAll[u, z, t]
eq = Derivative[2, 0][u][z, t] == Derivative[0, 1][u][z, t];
ic = {u[0, t] == 0, u[2, t] == 0, u[z, 0] == 10};

sol = u /. First@NDSolve[Flatten[{eq, ic}], u, {z, 0, 2}, {t, 0, 10}]

Plot3D[sol[z, t], {z, 0., 2.}, {t, 0., 10}, PlotRange->All, AxesLabel->{z, t, u[z, t]}]