Testing Turbulence Models

Question

In the 90s, several turbulence models were proposed. Using Mathematica vv 10-12.1 we tested two models:

1. Spalart-Allmaras turbulence model (SA) for aerodynamic applications, published in AIAA Paper 92 - 0439, 1992, see also wikipedia

2. turbulence model for a nonisothermal atmospheric boundary layer with roughness, buoyancy and flow rotation (TT), published on https://www.witpress.com/Secure/elibrary/papers/AIR96/AIR96028FU.pdf https://www.witpress.com/Secure/elibrary/papers/AIR97/AIR97011FU.pdf

Test. In 2D, at the initial moment of time, a linear velocity profile is set. Define a profile at a later time. Model SA depends on 10 parameters:

sigma = 2/3; kap = .41; cb1 = .1355; cb2 = .622; eps = 10^-6; d = 
 Sqrt[y^2 + eps^2]; cw1 = 
 cb1/kap^2 + (1 + cb2)/
   sigma; cw2 = .3; cw3 = 2; cv1 = 7.1; ct1 = 1; ct2 = 2; ct3 = 1.2; \
ct4 = .5; mu = 1.711 10^-5;

om = Sqrt[ (D[u[t, x, y], y] - D[u[t, x, y], x])^2];
S = om + (1 - (nu[t, x, y]/
         mu)/(1 + (nu[t, x, y]/
             mu)^4/(cv1^3 + (nu[t, x, y]/mu)^3))) nu[t, x, 
      y]/(kap d)^2;
r = nu[t, x, y]/S/(kap d)^2;
fw = (r + 
     cw2 (r^6 - r)) ((1 + cw3^6/(cw3^6 + (r + cw2 (r^6 - r))^6)))^(1/
      6);
ft2 = ct3 Exp[-ct4 (nu[t, x, y]/mu)^2];
nut = nu[t, x, y] (nu[t, x, y]/mu)^3/(cv1^3 + (nu[t, x, y]/mu)^3);
(*d=y for the wall*)

L = 10^4; t0 = 15; px = 0;

eq = {D[nu[t, x, y], t] + 
     u[t, x, y] D[nu[t, x, y], x] == (cb1 (1 - ft2) S  nu[t, x, y] + 
       1/sigma (D[(nut + mu) D[nu[t, x, y], y], y] + 
          D[(nut + mu) D[nu[t, x, y], x], x] + 
          cb2 (D[nu[t, x, y], x]^2 + D[nu[t, x, y], y]^2)) - (cw1 fw -
           cb1/kap^2 ft2) (nu[t, x, y]/d)^2 )/mu, 
   D[u[t, x, y], t] + u[t, x, y] D[u[t, x, y], x] + px == 
    D[(nut/mu + 1) D[u[t, x, y], y], y] + D[u[t, x, y], x, x]};
bc = {nu[t, x, 0] == 0, u[t, x, 0] == 0, u[t, x, L] == 1, 
   nu[t, x, L] == 0.1, nu[t, 0, y] == 0.1 y/L, u[t, 0, y] == y/L};
ic = {nu[0, x, y] == 0.1 y/L, u[0, x, y] == y/L};
{nU, U} = 
  NDSolveValue[{eq, ic, bc}, {nu, u}, {t, 0, t0}, {x, 0, L}, {y, 0, 
    L}];

Visualisation

{Plot3D[U[t, L/2, y], {t, 0, t0}, {y, 0, L}, Mesh -> None, 
  ColorFunction -> "Rainbow", AxesLabel -> Automatic], 
 Plot3D[nU[t, L/2, y], {t, 0, t0}, {y, 0, L}, Mesh -> None, 
  ColorFunction -> "Rainbow", AxesLabel -> Automatic]}

{Plot3D[U[t0, x, y], {x, 0, L}, {y, 0, L}, PlotRange -> All, 
  AxesLabel -> {x, y, ""}, PlotLabel -> "U", Mesh -> None, 
  ColorFunction -> "Rainbow"], 
 Plot3D[nU[t0, x, y], {x, 0, L}, {y, 0, L}, PlotRange -> All, 
  AxesLabel -> {x, y, ""}, PlotLabel -> "nU", Mesh -> None, 
  ColorFunction -> "Rainbow"]}

Model TT depends on 2 parameters:

T = 3; L = 1000; W0 = .00002; U0 = 1; R = 13.22; px = 0.13; b = \
.0001; eq = {D[W[x, y, t], t] + u[x, y, t] D[W[x, y, t], x] + 
    Integrate[W[x, y, t], {y, 0, y}]*D[W[x, y, t], y] - 
    2*R*y*D[W[x, y, t], y] - R*(1 + y^2)*D[W[x, y, t], y, y] - 
    D[W[x, y, t], x, x] == (y/(1 + y^2))*
     Integrate[W[x, y, t], {y, 0, y}] + b*y/(1 + y^2), 
  D[u[x, y, t], t] + u[x, y, t] D[u[x, y, t], x] + 
    Integrate[W[x, y, t], {y, 0, y}]*D[u[x, y, t], y] - 
    R*y*D[u[x, y, t], y] - R*(1 + y^2)*D[u[x, y, t], y, y] - 
    D[u[x, y, t], x, x] + px == 0};
ic = {W[x, y, 0] == W0*(y - L), 
  u[x, y, 0] == U0*y/L}; bc = {W[x, L, t] == 0, W[x, 0, t] == -W0*L, 
  W[0, y, t] == W0*(y - L), u[x, 0, t] == 0, u[x, L, t] == U0, 
  u[0, y, t] == U0*y/L}; bc1 = {Derivative[1, 0, 0][u][L, y, t] == 0, 
  Derivative[1, 0, 0][W][L, y, t] == 0};
sol = NDSolve[{eq, ic, bc}, {W, u}, {x, 0, L}, {y, 0, L}, {t, 0, T}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 150}}];

Visualisation

{Plot3D[W[L/2, y, t] /. First[sol], {y, 0, L}, {t, 0, T}, 
  PlotRange -> All, AxesLabel -> Automatic, PlotLabel -> W, 
  Mesh -> None, ColorFunction -> "Rainbow"], 
 Plot3D[u[L/2, y, t] /. First[sol], {y, 0, L}, {t, 0, T}, 
  PlotRange -> All, AxesLabel -> Automatic, PlotLabel -> U, 
  Mesh -> None, ColorFunction -> "Rainbow"]}

{Plot3D[W[x, y, T] /. First[sol], {x, 0, L}, {y, 0, L}, 
  PlotRange -> All, AxesLabel -> {x, y, ""}, PlotLabel -> W, 
  Mesh -> None, ColorFunction -> "Rainbow"], 
 Plot3D[u[x, y, T] /. First[sol], {x, 0, L}, {y, 0, L}, 
  PlotRange -> All, AxesLabel -> {x, y, ""}, PlotLabel -> U, 
  Mesh -> None, ColorFunction -> "Rainbow"]}

We see that both models describe the transition from a linear to a logarithmic profile in a finite time. We did not use bc1 in the last code, as this will increase the computing time several times with the same result. Why is time increasing?

#Update 1

I found a method to equalize computer time with bc1 and without bc1. It also allow us to answer the question. Code without bc1

 T = 3; L = 1000; W0 = .00002; U0 = 1;

 R = 13.22; px = 0.13; b = .0001; eq = \
{D[W[x, y, t], t] + u[x, y, t] D[W[x, y, t], x] + 
    Integrate[W[x, y, t], {y, 0, y}]*D[W[x, y, t], y] - 
    2*R*y*D[W[x, y, t], y] - R*(1 + y^2)*D[W[x, y, t], y, y] - 
    D[W[x, y, t], x, x] == (y/(1 + y^2))*
     Integrate[W[x, y, t], {y, 0, y}] + b*y/(1 + y^2), 
  D[u[x, y, t], t] + u[x, y, t] D[u[x, y, t], x] + 
    Integrate[W[x, y, t], {y, 0, y}]*D[u[x, y, t], y] - 
    R*y*D[u[x, y, t], y] - R*(1 + y^2)*D[u[x, y, t], y, y] - 
    D[u[x, y, t], x, x] + px == 0};
ic = {W[x, y, 0] == W0*(y - L), 
  u[x, y, 0] == U0*y/L}; bc = {W[x, L, t] == 0, W[x, 0, t] == -W0*L, 
  W[0, y, t] == W0*(y - L), u[x, 0, t] == 0, u[x, L, t] == U0, 
  u[0, y, t] == U0*y/L}; bc1 = {Derivative[1, 0, 0][u][L, y, t] == 0, 
  Derivative[1, 0, 0][W][L, y, t] == 0};
 Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[
 sol = NDSolve[{eq, ic, bc}, {W, u}, {x, 0, L}, {y, 0, L}, {t, 0, T}, 
    Method -> {"MethodOfLines", 
      "DifferentiateBoundaryConditions" -> False, 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MinPoints" -> 150}}, 
    EvaluationMonitor :> (currentTime = t;)];]

Computational time {172.642, Null}, and we have a message

NDSolve::bcart: Warning: an insufficient number of boundary conditions have been specified for the direction of independent variable x. Artificial boundary effects may be present in the solution. 

If we use homogeneous Neumann conditions at x=L= bc1, then the code is

Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[
 sol1 = NDSolve[{eq, ic, bc, bc1}, {W, u}, {x, 0, L}, {y, 0, L}, {t, 
     0, T}, Method -> {"MethodOfLines", 
      "DifferentiateBoundaryConditions" -> False, 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MinPoints" -> 150}}, 
    EvaluationMonitor :> (currentTime = t;)];]

In this case we have computational time {167.784, Null}, and there are no messages. Now we must explain why, without option "DifferentiateBoundaryConditions" -> False, time increases several times. Here we solve a system of integro-differential equations. With option, the problem boils down to solving a system of algebraic equations. Without this option, the system is transformed into a system of differential equations that are solved together.

#Update 2

As suggested by @xzczd, we test code without integrals. Integrate[W[x, y, t], {y, 0, y}] is transformed into new function D[V[x,y,t],y]==W[x,y,t], and from it we derive an equation D[V[x, y, t], y, t] + D[V[x, y, t], y, x] == D[W[x, y, t], t] + D[W[x, y, t], x]. With this equation we solve a problem and compare with the code from Update 1 containing integrals. There is not much difference between the two solutions. New code:

T = 3; L = 1000; W0 = .00002; U0 = 1; R = 13.22; px = 0.13; b = \
.0001; eq = {D[W[x, y, t], t] + u[x, y, t] D[W[x, y, t], x] + 
V[x, y, t]*D[W[x, y, t], y] - 2*R*y*D[W[x, y, t], y] - 
    R*(1 + y^2)*D[W[x, y, t], y, y] - 
    D[W[x, y, t], x, x] == (y/(1 + y^2))*V[x, y, t] + b*y/(1 + y^2), 
  D[u[x, y, t], t] + u[x, y, t] D[u[x, y, t], x] + 
    V[x, y, t]*D[u[x, y, t], y] - R*y*D[u[x, y, t], y] - 
    R*(1 + y^2)*D[u[x, y, t], y, y] - D[u[x, y, t], x, x] + px == 0, 
  D[V[x, y, t], y, t] + D[V[x, y, t], y, x] == 
   D[W[x, y, t], t] + D[W[x, y, t], x]};
ic = {W[x, y, 0] == W0*(y - L), u[x, y, 0] == U0*y/L, 
  V[x, y, 0] == W0 (y^2/2 - L y)}; bc = {W[x, L, t] == 0, 
  W[x, 0, t] == -W0*L, W[0, y, t] == W0*(y - L), u[x, 0, t] == 0, 
  u[x, L, t] == U0, u[0, y, t] == U0*y/L};
bcV = {V[x, 0, t] == 0, 
  V[0, y, t] == 
   W0 (y^2/2 - L y)}; bc1 = {Derivative[1, 0, 0][u][L, y, t] == 0, 
  Derivative[1, 0, 0][W][L, y, t] == 0};
 Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[
 sol = NDSolve[{eq, ic, bc, bc1, bcV}, {W, V, u}, {x, 0, L}, {y, 0, 
     L}, {t, 0, T}, 
    Method -> {"MethodOfLines", 
      "DifferentiateBoundaryConditions" -> False, 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MinPoints" -> 150}}, 
    EvaluationMonitor :> (currentTime = t;)];]  

Update 3 Thanks to user21 remarks about SA model, I have modified bc for SA as

 bc = {nu[t, x, 0] == 0, u[t, x, 0] == 0, u[t, x, L] == 1, 
   nu[t, x, L] == 0.1, nu[t, 0, y] == 0.1 y/L, u[t, 0, y] == y/L, 
   Derivative[0, 1, 0][nu][t, L, y] == 0, 
   Derivative[0, 1, 0][u][t, L, y] == 0};

With the Neumann boundary condition at x==L we have numerical solution without message from NDSolve. New solution for nu looks different from shown above

Answer 1

Seems that Alex has solved the problem himself, but I still want to compete for the bounty think it's still better to elaborate my points in the comments.

It should be noted that, NDSolve actually doesn't parse the Integrate[…] correctly. This can be verified by checking NDSolve`StateData[…]：

With[{W = W[x, y, t], u = u[x, y, t]}, 
  eq = {D[W, t] + u D[W, x] + Integrate[W, {y, 0, y}] D[W, y] - 2 R y D[W, y] - 
      R (1 + y^2) D[W, y, y] - D[W, x, x] == 
     y/(1 + y^2) Integrate[W, {y, 0, y}] + (b y)/(1 + y^2), 
    D[u, t] + u D[u, x] + Integrate[W, {y, 0, y}] D[u, y] - R y D[u, y] - 
      R (1 + y^2) D[u, y, y] - D[u, x, x] + px == 0};
  ic = {W == W0*(y - L), u == U0*y/L} /. t -> 0; 
  bc = {{W == 0, u == U0} /. y -> L, {W == -W0*L, u == 0} /. 
     y -> 0, {W == W0*(y - L), u == U0*y/L} /. x -> 0}; 
  bc1 = {D[u, x] == 0, D[W, x] == 0} /. x -> L];
{state} = NDSolve`ProcessEquations[{eq, ic, bc, bc1}, {W, u}, {x, 0, L}, {y, 0, L}, {t, 
    0, T}];
func = state["NumericalFunction"]["FunctionExpression"];

The output involves messy intermediate variables and long list, we make some replacements to make it easier to read:

rule = Cases[func, 
  HoldPattern@(var_ = NDSolve`FiniteDifferenceDerivativeFunction[d_, __][h_]) :> (var -> 
     d@h), Infinity]

func /. rule /. HoldPattern[y = lst_] :> (y = Short@lst) 

Pictured by Simon Wood's shadow.

Comparing with the original system, it's not hard to notice $\int_0^y W(x,y,t) d y$ becomes $W(x,y,t) y$ inside NDSolve, probably because the integrand has been treated as constant.

This can be further verified by replacing Integrate[W, {y, 0, y}] with W y in eq and comparing the outputs of NDSolve.

As to bcart warning, I still don't think it's a good idea to bear with it, despite it doesn't seem to cause significant problem in this case. Readers are interested in the topic please check this post.

To resolve the issue, I think the approach in update 2 is the right way to go, the b.c. V[x, L, t] == -W0 L^2/2 seems to be redundant though.

BTW, it's good to see the DAE solver is improved in recent versions. The code in update 2 just crashes the kernel in v9.0.1.

---

To make this answer more interesting, I'd like to add a solution that also works in v9.0.1. pdetoode is used to discretize the PDE system to an ODE system:

T = 3; L = 1000; W0 = .00002; U0 = 1; R = 13.22; px = 0.13; b = .0001;
With[{W = W[x, y, t], u = u[x, y, t], V = V[x, y, t]}, 
  eq = {D[W, t] + u D[W, x] + V D[W, y] - 2 R y D[W, y] - R (1 + y^2)*D[W, y, y] - 
      D[W, x, x] == y/(1 + y^2) V + b y/(1 + y^2), 
    D[u, t] + u D[u, x] + V D[u, y] - R y D[u, y] - R (1 + y^2) D[u, y, y] - D[u, x, x] +
       px == 0, D[V, y, t] + D[V, y, x] == D[W, t] + D[W, x]};
  ic = {W == W0*(y - L), u == U0*y/L, V == W0 (y^2/2 - L y)} /. t -> 0;
  bc = {{W == 0, u == U0} /. y -> L, {W == -W0*L, u == 0} /. 
     y -> 0, {W == W0*(y - L), u == U0*y/L} /. x -> 0};
  bcV = {(*V\[Equal]-W0 L^2/2/.y\[Rule]L,*)V == 0 /. y -> 0, 
    V == W0 (y^2/2 - L y) /. x -> 0};
  bc1 = {D[u, x] == 0, D[W, x] == 0} /. x -> L];

points@x = points@y = 100; difforder = 2;
domain@x = domain@y = {0, L};
(grid@# = Array[# &, points@#, domain@#]) & /@ {x, y};

(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
ptoofunc = pdetoode[{W, u, V}[x, y, t], t, grid /@ {x, y}, difforder];
del = #[[2 ;; -2]] &;
delL = Rest;
ode = {del /@ del@# & /@ ptoofunc@eq[[1 ;; 2]], delL /@ delL@ptoofunc@eq[[-1]]};
odeic = ptoofunc@ic;
odebc = With[{sf = 1}, 
   diffbc[t, sf]@{Map[del, ptoofunc@bc[[1 ;; 2]], {2}], ptoofunc@bc[[3]], 
     ptoofunc@bc1, delL@ptoofunc@bcV[[1]], ptoofunc@bcV[[2]]}];
var = Outer[#[#2, #3] &, {W, u, V}, grid@x, grid@y, 1];
sollst = NDSolveValue[{ode, odeic, odebc}, var, {t, 0, T}, 
                      SolveDelayed -> True]; // AbsoluteTiming
(* {54.518346, Null} *)
sol = {W, u, V} -> (rebuild[#, {grid@x, grid@y}, 3] & /@ sollst) // Thread

The option SolveDelayed is red, but don't worry. Alternatively you can use Method -> {"EquationSimplification" -> "Residual"}.

Limited by the RAM of my laptop, I only use 100 points for each dimension, but the result is already good:

plot[expr_, rangex_, rangey_, label_] := 
 Plot3D[expr, rangex, rangey, PlotRange -> All, AxesLabel -> Automatic, 
  PlotLabel -> label, ColorFunction -> "AvocadoColors", Mesh -> None]

GraphicsGrid[
 {{plot[W[L/2, y, t] /. sol, {y, 0, L}, {t, 0, T}, W], 
   plot[u[L/2, y, t] /. sol, {y, 0, L}, {t, 0, T}, U]},
  {plot[W[x, y, T] /. sol, {x, 0, L}, {y, 0, L}, W], 
   plot[u[x, y, T] /. sol, {x, 0, L}, {y, 0, L}, U]}}, ImageSize -> Large]

Answer 2

T = 3; L = 1000; W0 = .00002; U0 = 1; R = 13.22; px = 0.13; b = \
.0001; eq = {D[W[x, y, t], t] + u[x, y, t] D[W[x, y, t], x] + 
    Integrate[W[x, y, t], {y, 0, y}]*D[W[x, y, t], y] - 
    2*R*y*D[W[x, y, t], y] - R*(1 + y^2)*D[W[x, y, t], y, y] - 
    D[W[x, y, t], x, x] == (y/(1 + y^2))*
     Integrate[W[x, y, t], {y, 0, y}] + b*y/(1 + y^2), 
  D[u[x, y, t], t] + u[x, y, t] D[u[x, y, t], x] + 
    Integrate[W[x, y, t], {y, 0, y}]*D[u[x, y, t], y] - 
    R*y*D[u[x, y, t], y] - R*(1 + y^2)*D[u[x, y, t], y, y] - 
    D[u[x, y, t], x, x] + px == 0};
ic = {W[x, y, 0] == W0*(y - L), 
  u[x, y, 0] == U0*y/L}; bc = {W[x, L, t] == 0, W[x, 0, t] == -W0*L, 
  W[0, y, t] == W0*(y - L), u[x, 0, t] == 0, u[x, L, t] == U0, 
  u[0, y, t] == U0*y/L}; bc1 = {Derivative[1, 0, 0][u][L, y, t] == 0, 
  Derivative[1, 0, 0][W][L, y, t] == 0};
sol = NDSolve[{eq, ic, bc, bc1}, {W, u}, {x, 0, L}, {y, 0, L}, {t, 0, 
    T}, Method -> "StiffnessSwitching"];

{Plot3D[W[L/2, y, t] /. First[sol], {y, 0, L}, {t, 0, T}, 
  PlotRange -> All, AxesLabel -> Automatic, PlotLabel -> W, 
  Mesh -> None, ColorFunction -> "Rainbow"], 
 Plot3D[u[L/2, y, t] /. First[sol], {y, 0, L}, {t, 0, T}, 
  PlotRange -> All, AxesLabel -> Automatic, PlotLabel -> U, 
  Mesh -> None, ColorFunction -> "Rainbow"]}

{Plot3D[W[x, y, T] /. First[sol], {x, 0, L}, {y, 0, L}, 
  PlotRange -> All, AxesLabel -> {x, y, ""}, PlotLabel -> W, 
  Mesh -> None, ColorFunction -> "Rainbow"], 
 Plot3D[u[x, y, T] /. First[sol], {x, 0, L}, {y, 0, L}, 
  PlotRange -> All, AxesLabel -> {x, y, ""}, PlotLabel -> U, 
  Mesh -> None, ColorFunction -> "Rainbow"]}

Alternative:

sol = NDSolve[{eq, ic, bc, bc1}, {W, u}, {x, 0, L}, {y, 0, L}, {t, 0, 
    T}, Method -> {"ExplicitRungeKutta", "StiffnessTest" -> False}];

sol = NDSolve[{eq, ic, bc, bc1}, {W, u}, {x, 0, L}, {y, 0, L}, {t, 0, 
    T}, Method -> {"ExplicitRungeKutta", "StiffnessTest" -> True}];

NDSolve[{eq, ic, bc, bc1}, {W, u}, {x, 0, L}, {y, 0, L}, {t, 0, T}, 
 Method -> {"ExplicitRungeKutta", 
   "StiffnessTest" -> {True, "MaxRepetitions" -> {1, 1}, 
     "SafetyFactor" -> 1} }]

Using other sets of methods can bring the system of equation closer to the better one. The results shed some light on the critics of the results of the given better solution but did not really make the solution better on all points of interest.

The options are taken from an example in Stiffnesstest.

With this:

Tt = 3; L = 0.1; W0 = 2/1000000(*.00002*); U0 = 1; R = 
 1322/10000 (*13.11*); px = 0(*0.13*); b = 
 1/100000; eq = {D[W[x, y, t], t] + u[x, y, t] D[W[x, y, t], x] + 
    Integrate[W[x, y, t], {y, 0, y}]*D[W[x, y, t], y] - 
    2*R*y*D[W[x, y, t], y] - R*(1 + y^2)*D[W[x, y, t], y, y] - 
    D[W[x, y, t], x, x] == (y/(1 + y^2))*
     Integrate[W[x, y, t], {y, 0, y}] + b*y/(1 + y^2), 
  D[u[x, y, t], t] + u[x, y, t] D[u[x, y, t], x] + 
    Integrate[W[x, y, t], {y, 0, y}]*D[u[x, y, t], y] - 
    R*y*D[u[x, y, t], y] - R*(1 + y^2)*D[u[x, y, t], y, y] - 
    D[u[x, y, t], x, x] + px == 0};
ic = {W[x, y, 0] == W0*(y - L), 
  u[x, y, 0] == U0*y/L}; bc = {W[x, L, t] == 0, W[x, 0, t] == -W0*L, 
  W[0, y, t] == W0*(y - L), u[x, 0, t] == 0, u[x, L, t] == U0, 
  u[0, y, t] == U0*y/L}; bc1 = {Derivative[1, 0, 0][u][L, y, t] == 0, 
  Derivative[1, 0, 0][W][L, y, t] == 0};
sol = NDSolve[{eq, ic, bc, bc1}, {W, u}, {x, 0, L}, {y, 0, L}, {t, 0, 
    Tt}, Method -> "StiffnessSwitching"];

I get different results:

{Plot3D[W[L/2, y, t] /. First[sol], {y, 0, L}, {t, 0, Tt}, 
  PlotRange -> All, AxesLabel -> Automatic, PlotLabel -> W, 
  Mesh -> None, ColorFunction -> "Rainbow"], 
 Plot3D[u[L/2, y, t] /. First[sol], {y, 0, L}, {t, 0, Tt}, 
  PlotRange -> All, AxesLabel -> Automatic, PlotLabel -> U, 
  Mesh -> None, ColorFunction -> "Rainbow"]}

{Plot3D[W[x, y, Tt] /. First[sol], {x, 0, L}, {y, 0, L}, 
  PlotRange -> All, AxesLabel -> {x, y, ""}, PlotLabel -> W, 
  Mesh -> None, ColorFunction -> "Rainbow"], 
 Plot3D[u[x, y, Tt] /. First[sol], {x, 0, L}, {y, 0, L}, 
  PlotRange -> All, AxesLabel -> {x, y, ""}, PlotLabel -> U, 
  Mesh -> None, ColorFunction -> "Rainbow"]}

It seems to be a direct path into chaos immanent in nonlinear differential equations. The calmed solution is not very good.

It shows the ripples from the calmed solution much more at the ground floor level than atop the great buckle. It provides more similar solutions than the calmed one for all of the four plots. It is not so a great domain than the other ones. The domain is now {-500000,500000} and {-10^7,10^7}. It is not all positive as might be physical, but it is plain in most parts of the defined domain for {t,x,y}.

It first attempted to make the domain smaller. That failed and only proved, time is passing by for the system.

Second I altered the parameters since this seems to gain more insight into the behaviour of the systems under consideration. The did the trick. It resembled, on the other hand, the most important critics from the scientific community for fluid dynamics on the model under consideration. The calming might be due to implicit change in the parameters. That too is possibly still under the regime of chaos in this system.

Nevertheless, it still has potential that chaos is introduced by the methods used to solve the problem and not just the parameters in use. The results presented here are chosen due to physical consideration. As far as I know, this is the first time such results are published for the problem presented here. This is not critics of the methods in NDSolve as offered at present by Wolfram Research.

The computational experiment shows up the fine power of StiffnessSwitching on a very stiff problem with immense singularities of not so point-like character.