# Testing Turbulence Models

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 https://en.wikipedia.org/wiki/Spalart–Allmaras_turbulence_model

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;)];]

• bc1 is ill formulated. Please revise the given information. The MinPoints is in this form the limit. There is no real solution calculated. This reflects the conditions of the differential equation. It is really stiff. I can reproduce but no further. In general, bc1 is essential. NDSolve searches for a path to solve but does not solve. An evaluation monitor does not start during processing time. Search SE for the error type "The arguments should be ordered consistently" or if Dt for order must be lower than that of the differential equation to help yourself. Mar 28 '20 at 22:19
• bc1 works in v.12 with options, and in v.12.1 without options. Which version are you using? Mar 29 '20 at 10:02
• I use 12.0.0. I get results only without bc1. I get better results for the second model if I leave the MethodOfLines away. The third graphics looks rounded for small y. For larger t the values of W are lower. I tried Trace with NDSolveMethodData reports Mathematica suggests, prefers LSODA for this problem if allow to select. The neat for the stepsize correction for the spatial error appears later. Mar 29 '20 at 17:06
• In version 12.0, the code works with bc1, although longer than without bc1. But in version 12.1, it only works without options on a coarse grid of 15 points. Mar 30 '20 at 15:53
• "If we do not use bc1, then NDSolve[] automatically applies the homogeneous Neumann conditions at x=L. " It's not the case for TensorProductGrid, the zero Neumann value is automatically applied only for FiniteElement method. Please notice the bcart warning is a serious problem. A detailed discussion can be found here: mathematica.stackexchange.com/q/73961/1871 Apr 3 '20 at 3:16

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 NDSolveStateData[…]

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} = NDSolveProcessEquations[{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_ = NDSolveFiniteDifferenceDerivativeFunction[d_, __][h_]) :> (var ->
d@h), Infinity]

func /. rule /. HoldPattern[y = lst_] :> (y = Short@lst) 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,
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},
Map[sf # + D[#, t] &,
Flatten@{Map[del, ptoofunc@bc[[1 ;; 2]], {2}], ptoofunc@bc[], ptoofunc@bc1,
delL@ptoofunc@bcV[], 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] • Thank you for the detailed answer. Now I see that I did not spend time in vain to get code Update 2. This eliminates both problems with bc1 and Integrate. Apr 5 '20 at 11:50
• @alex One thing to supplement: V[x, L, t] == -W0 L^2/2 seems to be redundant. Apr 5 '20 at 12:38
• Yes, I agree, this condition can be removed. I wonder why with pdetoode code is always faster? Apr 5 '20 at 17:11
• @Alex If Method -> {"MethodOfLines", "DifferentiateBoundaryConditions" -> False, "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 100, "MaxPoints" -> 100, "DifferenceOrder" -> 2}} is added to your code, the timings become similar according to my test :) . (Mathematica 12.1, Win10 64bit. ) But indeed, there exist many cases that pdetoode makes the code faster (for example the dynamic Euler–Bernoulli beam equation linked above), whose reason remains a mystery to me. My guess is an explicit ODE system can be analysed better by NDSolve, but it's no more than a guess. Apr 6 '20 at 2:57
• It looks like in v.12.0 we have to use additional option sollst = NDSolveValue[{ode, odeic, odebc}, var, {t, 0, T}, Method -> {"EquationSimplification" -> "Residual"}] Apr 6 '20 at 20:47
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.

• At first glance, this does not solve the problem, since a grid of 15 points is used, and we want 150. Apr 2 '20 at 12:15
• Please, check Update 2 with a new system without Integrate. Apr 6 '20 at 20:52