10
$\begingroup$

So I am trying to solve the movement in space and time of a spreading gravity current. The interface satisfies the following PDE:

$ \frac{\partial h}{\partial t} = \frac{\partial}{\partial x}\left(h^3 \frac{\partial h}{\partial x}\right) $.

The nose, $h(x,t) = 0$, of the current can move forward in space. If we say the initial condition is $h(x,0) = 1-x$, then the nose, $x_N$ is initially at $x = 1$. The base of the current can be fixed at $h(0,t) = 1$. So I have one IC and one BC, I need another BC to close the system. Now, as time moves forward the nose (fixed at h = 0) propagates according to the expression:

$\int_0^{x_N} h \, dx = \frac{1}{2} + t$

or kinematically $\dot{x_N} = -h^2 h_x$.

Clearly the second boundary condition in $x$ needs to come from the nose condition, but I have no idea how to input an integral boundary condition or the kinematic condition. Secondly, I have no idea how to deal with xmax in NDSolve as clearly my xmax is moving forward with each time step.

Here is my attempt, which gives a "There are fewer dependent variables error"

NDSolve[{D[h[x, t], t] == D[h[x, t]^3 D[h[x, t], x], x], h[0, t] == 1,
h[x, 0] == 1 - x (D[h[x, t], t] /. x -> 1) == -h[x, t]^2 (D[h[x, t], x] /. 
  x -> 1)}, h, {x, 0, 1}, {t, 0, 1}]

EDIT

So I solved this problem by amending the code in this Stefan problem post. I modified the statement of the problem slightly to work more closely to @ybeltukov excellent solution, and actually solved the following system:

$ \frac{\partial h}{\partial t} = \frac{\partial}{\partial x}\left(h^3 \frac{\partial h}{\partial x}\right) \\ \dot{s} = -h^2 h_x \\ s(0) = 0 \\ h(x,0) = \begin{cases}1 \quad \text{for } x = 0 \\ 0 \quad \text{otherwise}\end{cases} \\ h(s(t),t) = 0 \\ h(0,t) = 1 $.

I made the same change to a normalised variable for my PDE and then amended finite difference code:

n = 100;
\[Delta]\[Xi] = 1./n;

ClearAll[dv, t];
dv[v_List] := 
With[{s = First@v, u = Rest@v}, 
With[{ds = u[[-1]]^3/(s \[Delta]\[Xi]), \[Xi] = N@Range[n - 1]/n, 
 d1 = ListCorrelate[{-0.5, 0, 0.5}/\[Delta]\[Xi], #] &, 
 d2 = ListCorrelate[{1, -2, 1}/\[Delta]\[Xi]^2, #] &}, 
Prepend[3 u^2 d1[#]^2/s^2 + u^3 d2[#]/s^2 + \[Xi] ds d1[#]/s &@
  Join[{1}, u, {0.}], ds]]];
 s0 = 0.001;
 v0 = Flatten@Prepend[ConstantArray[0.001, n - 2], {s0, 1.}];
 sol = NDSolve[{v'[t] == dv[v[t]], v[0] == v0}, v, {t, 0, 50}][[1, 1, 
2]];

Returning from the normalised variable is identical code.

You can then make a nice time series with this code:

ListAnimate[
Table[Plot[u[t, x], {x, 0, 10}, PlotRange -> {{0, 10}, {0, 3}}], {t, 
 0, 50, 0.1}]]

Which returns a spreading current as desired.

Spreading viscous gravity current

$\endgroup$
  • $\begingroup$ Could you explain what do you mean by the kinematic condition $\dot{x}_N=-h^2h_x$? Do I understand correctly that the right-hand-side should be evaluated at $x=0$? or at $x=x_N$ (for which it vanishes)? $\endgroup$ – yohbs May 12 '16 at 19:30
  • $\begingroup$ Also, your initial condition is inconsistent with the integral constraint because at $t=0$ you have $$\int_0^{x_N} h\, dx=\int_0^1(1-x)\,dx=\frac{1}{2}\ne t=0$$ $\endgroup$ – yohbs May 12 '16 at 19:38
  • $\begingroup$ Similar to mathematica.stackexchange.com/questions/109201/… $\endgroup$ – Dr. Wolfgang Hintze May 12 '16 at 20:28
  • $\begingroup$ @ojlm I'm not sure if your statements about the moving nose are correct. You might wish to consult my solution, and particularly the steady state solution. $\endgroup$ – Dr. Wolfgang Hintze May 12 '16 at 20:35
  • $\begingroup$ @Dr.WolfgangHintze The nose needs to propagate downstream as this is a model for a fluid mechanical process. Pinning the interface to the bottom boundary at x=1 is not physically what happens. I have given a link to a paper in my other comment. Governing equations are on page 3. I have simplified the propagation flux condition as I am only considering a one layered spreading current in this model. $\endgroup$ – M Horsley May 12 '16 at 20:44
3
$\begingroup$

The PDE can be easily solved numerically if a simple additional boundary condition is imposed.

Adding a second bundary condition h[1,t] = 0 we find numerically, without any error message

hh[x_, t_] = 
  h[x, t] /. 
   NDSolve[{D[h[x, t], t] == D[h[x, t]^3 D[h[x, t], x], x], h[0, t] == 1, 
      h[1, t] == 0, h[x, 0] == (1 - x)}, h[x, t], {x, 0, 1}, {t, 0, 5}][[1]];

Plot3D[hh[x, t], {x, 0, 1}, {t, 0, 3}, PlotRange -> {0, 1}, 
 PlotLabel -> "Solution of a PDE", AxesLabel -> {"x", "t", "h[x,t]"}]

enter image description here

Executing NDSolve without a secod boundary condition leads to about the same graph and the error message that a second boundary condition is missing.

The time dependence at various positions is given by the graph

Plot[{hh[0.2, t], hh[0.5, t], hh[0.8, t], hh[0.95, t]}, {t, 0, 3}, 
 PlotRange -> {0, 1}, 
 PlotLabel -> 
  "Solution of a PDE\nTime dependence for various locations x = 1 - h[x,0]", 
 AxesLabel -> {"t", "h[x,t]"}]

enter image description here

The asymptotic function for large t is given by

ha[x_] = (1 - x)^(1/4)

which solves the steady state ODL derived from the original PDE.

We could study the approch of the solution to the equlibirum by linearizing the equation, i.e. letting

h = hL[x] + ha[x]

and dropping all non linear terms in hL.

Partial solutions are given by Bessel functions. The interested reader will surely be able to follow this path on his own.

Moving structures

Let us have a look at moving structures allowed by the PDE.

Letting

h[x,t] = g[x - v t]

The PDE reduces to an ODE which can be solved

gg[x_] = g[x] /. 
   DSolve[g'[x] v == 3 g[x]^2 g'[x]^2 + g[x]^3 g''[x], g[x], x][[1]] /. 
  C[2] -> 0

(* Out[341]= InverseFunction[-((C[1]^3 Log[C[1] + v #1])/v^4) + (C[1]^2 #1)/v^3 - (
    C[1] #1^2)/(2 v^2) + #1^3/(3 v) &][x] *)

Defining

h1[s_] = (gg[x][[0, 1]] /. {C[1] -> -10, v -> 1} // FullSimplify)@s

(* Out[393]= 100 s + 5 s^2 + s^3/3 + 1000 Log[-10 + s] *)

Gives

Plot[Re[h1[s]], {s, 0, 10}, PlotRange -> {0, 3540}, 
 PlotLabel -> "'Nose' profile", AxesLabel -> {"x - t", "g"}]

enter image description here

This section is just a sketch which must be worked out more thoroughly.

$\endgroup$
  • $\begingroup$ Thanks for the response. I know it's possible to solve the PDE with h[1,t] = 0 boundary condition. However, since this is the propagation of a current, the nose needs to move over the bottom boundary. I want to be able to to create a time series such that the current moves along down the x axis. If you see this paper, you can see how the current propgates down the x-axis as time progresses. Having the h[1,t] = 0 boundary condition does not allow me this. repository.cam.ac.uk/bitstream/handle/1810/247602/… $\endgroup$ – M Horsley May 12 '16 at 20:41
  • $\begingroup$ I appreciate your further analysis about travelling wave solutions. This problem is well known and the best way to study it analytically is via a similarity solution. However, the point is I am using this model in an illustrative way. I am studying a much harder system, that has an interface moving transiently along a boundary. The key thing I am interested in do is to create a time series of the spreading thinning profile as seen in the paper I sent. Would you be able to help me in achieving this? I want to be have the profile spreading and thinning with accordance to the volume BC $\endgroup$ – M Horsley May 13 '16 at 8:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.