Working around inconsistent Boundary Conditions for a PDE

Question

I have a question concerning how to solve a 1D PDE with inconsistent boundary conditions using mathematica. The function U is in x and t, with two derivatives in x and one in t. Here is what the PDE looks like (foam drainage): $$ c_1 U \frac{\delta U}{\delta x} + c_2 (\frac{1}{2\sqrt{U}} (\frac{\partial U}{\partial x})^{2} + \sqrt{U}\frac{\partial^2 U}{\partial x^2}) - \frac{\partial U}{\partial x} = \frac{\partial U}{\partial t} $$

Where C1 and C2 are constant parameters. (Yes, this could be simplified somewhat). We need three boundary conditions. For the IC

$$U[x,t=0] = z\:(a\:constant) $$

For the BCs there is first: $$U[x=0,t] = z\:(the\:same\:constant) $$

and the BC that is causing me to lose sleep:
$$U^{3/2}[x=L,t] = -\frac{c_2}{c_1}\frac{\partial U}{\partial x}[x=L,t] $$ that is, at U[x=L,t] the first derivative in x is non-zero (it is "z^(3/2)", in my code, z = 1), but U[x,t=0] is constant (and so a derivative in x for U[x=L,t=0] should be zero). Thus the inconsistent boundary conditions.

In order to clarify what I am doing, here is the code block I built up from the heat equation step by step. I didn't have any problems until adding the Neumann BC for this problem:

C1 = 0.1; 
C2 = 0.1; 
L = 1; 
tfinal34 = 1;  
GovEqFD = 
  C2 Sqrt[u[x, t]] D[u[x, t], {x, 2}] + 
  C2 (1/(2 Sqrt[ u[x, t]])) (D[u[x, t], x])^2 +  
  C1 (u[x, t] - 1) D[u[x, t], x] - D[u[x, t], t] == 0  
BC34 = 
  {u[0, t] == 1, 
   -C2/C1 (D[u[x, t], x] /. x -> L) == (u[L, t])^(3/2) /. x -> L)]}  
IC34 = u[x, 0] == 1  
EqsFD4 = Join[{GovEqFD}, BC34, {IC34}]

FDsol4 = NDSolve[EqsFD4, u[x, t], {x, 0, L}, {t, 0, tfinal34}]

In order to try to find a solution around the inconsistent boundary conditions I have looked into the Method of Lines, but to be honest, I haven't gone too far into trying to figure out how to code it yet. It's not clear to me from the documentation exactly how complicated it would be to use FDM on the above governing equation.

In attempting a work around, I have thought to use WhenEvent to define the boundary condition @ t= 0 and also for t>0. I didn't have any problem with using WhenEvent in terms of getting warnings.

I changed the list for BC34 to:

BC34 = 
  {u[0, t] == 1, 
   -C2/C1 (D[u[x, t], x] /. {x -> L,t -> 0}) == 0, 
   WhenEvent[
     (t > 0), (-C2/C1 (D[u[x, t], x] /. x -> L) == 
       (u[L, t])^(3/2) /. x -> L)]}

As I said there were no warnings, and I saw that beloved output for the interpolating function. However, in plotting I see the solution is the same unstable plane it shouldn't be for the inconsistent solution. It seems that the BC remains -C2/C1 (D[u[x, t], x] /. x -> L) == 0, for times greater than t==0. Can someone help me understand if I am using WhenEvent incorrectly, or how I might implement it in this case to resolve the inconsistency in the BC @ x=L and IC? Or perhaps suggest other solution strategies?

If there is anyway I can clarify my question, formatting, or the information I am providing, please let me know.

Answer 1

Do pretty much the same as what you have:

C1 = 1/10;
C2 = 1/10;
L = 1;
tfinal34 = 6;

GovEqFD = 
 C2 Sqrt[u[x, t]] D[u[x, t], {x, 2}] + 
   C2 (1/(2 Sqrt[u[x, t]])) (D[u[x, t], x])^2 + 
   C1 (u[x, t] - 1) D[u[x, t], x] - D[u[x, t], t] == 0;

BC34 = {u[0, t] == 1, (D[u[x, t], x] /. x -> L) == -C1/C2 ((u[L, t])^(3/2) /. x -> L)};

The trick is choosing an IC that closely matches the IC we want, but is not too inconsistent with our BC's. I generally try to tack on a parabola on the IC at the BC that is unmatched. Let the deviation from our original IC to a parabola occur at x = u1. In this case u0 is the constant value of the original IC. With the equation of a parabola with vertex at {u0,u1}, use a as the parameter to match the BC at x = L.

u0 = 1;
u1 = 99/100;

y[x_] = u0 - a (x - u1)^2

Find the value of a that matches the right slope at the boundary.

(y'[x] /. x -> L) == -y[L]^(3/2) // N

a = a /. Solve[%, a][[1]]
(*49.6283*)

Transition from the constant value u0 to the parabola with a UnitStep.

IC34 = u[x, 0] == u0 - (UnitStep[x - u1]) a (x - u1)^2

This is what our adjusted IC looks like at the end.

Plot[IC34[[2]], {x, .98, 1}]

The transition is fairly smooth and the deviation from the desired IC is small.

EqsFD4 = Join[{GovEqFD}, BC34, {IC34}];

solu = FDsol4 = 
  NDSolve[EqsFD4, u[x, t], {x, 0, L}, {t, 0, tfinal34}, 
   PrecisionGoal -> Infinity, MaxStepFraction -> 1/200, 
   MaxSteps -> {50000, Automatic}];

u[x_, t_] = u[x, t] /. %[[1]]

gifs = Table[
   Plot[u[x, t], {x, 0, 1}, PlotRange -> {.5, 1}, PlotLabel -> t "t", 
    AxesLabel -> {"x", "u"}], {t, 0, tfinal34, .1}];
ListAnimate[gifs]