# NDSolve mathematica fluid mechanics

I have made Mathematica code for equation 16 from this paper. How can we make these streamlines I attached? I am using DensityPlot command but can't get it right and if I have to make plot of dp/dx and x. What changes should I do? I am just replacing q values with x. Is it OK?

M == 0.5
\[Alpha] == 0.2
n == 0.5
\[Theta] == 0.5
We == 0.02
\[Epsilon] == 0.1
\[Delta] == 0.0006
h == 1 + \[Epsilon]Cos[2 \[Pi] x]
F == \!$$\*SubsuperscriptBox[\(\[Integral]$$, $$0$$, $$h$$]$$\*SubscriptBox[\(\[PartialD]$$, $$y$$]w \[DifferentialD]y\)\)

sol = NDSolve[{((w^
'')[y]) D[1 + n (We ((w^
'')[y]) - 1), {y, 2}] - M^2 (Cos[\[Theta]]^2)
\!$$\*SuperscriptBox[\(w$$, $$''$$,
MultilineFunction->None]\)[y] == 0 ,
\!$$\*SuperscriptBox[\(w$$, $$''$$,
MultilineFunction->None]\)[0] == 0, w[0] == 0,
w'[h] == -1 - 2 \[Pi]\[Epsilon]\[Alpha]\[Delta]Cos[2 \[Pi] x],
w[h] == F}, w[y], {y, 0, h}]

• is this an ode or pde? what is $\psi$ a function of? just $y$ ? then why the book uses partial derivatives? Nov 20, 2021 at 17:55
• It could be better to define hydrodynamic problem as it is, without method of solution described in the book. Then we can solve it with FEM or FDM and compare with method from the book. See update to my answer. Nov 27, 2021 at 7:27
• It is not clear in what coordinates stream plot in Figure 9 has been created. Dec 12, 2021 at 10:50
• For completeness, please give a link to the paper you have taken a screenshot of (apparently, it's from World Scientific). Dec 13, 2021 at 10:35
• worldscientific.com/doi/10.1142/S0219887819501391
– umar
Dec 13, 2021 at 10:43

We can reproduce equation (16) and numerical solution as follows

M = 0.5;
\[Alpha] = 0.2;
n = 0.5;
\[Theta] = 0.5; m = M^2 Cos[\[Theta]]^2;
We = 0.02;
\[Epsilon] = 0.1;
\[Delta] = 0.0006; F = 1; x = 0;
h = 1 + \[Epsilon] Cos[2 \[Pi] x];
eq = {psi''[y] == u[y],
D[(1 + n (We u[y] - 1)) u[y], y, y] - m u[y] == 0}; bc = {psi[0] ==
0, u[0] == 0,
psi'[h] == -1 - 2 Pi \[Epsilon] \[Alpha] \[Delta] Cos[2 Pi x],
psi[h] == F};
sol = NDSolve[Join[eq, bc], {psi, u}, {y, 0, h}]


Visualization of velocity profile $$\psi '(y)$$

Plot[psi'[y] /. sol[[1]], {y, 0, h}]


Using Module we can reproduce Figure 1 from paper cited

M = 0.5;
\[Alpha] = 0.2;
n = 0.5;
\[Theta] = 0.5; m = M^2 Cos[\[Theta]]^2;
We = 0.02;
\[Delta] = 0.0006; x = 2/3;
p[eps_, qu_] := Module[{e = eps, q = qu}, h = 1 + e Cos[2 \[Pi] x];
eq = {psi''[y] == u[y],
D[(1 + n (We u[y] - 1)) u[y], y, y] - m u[y] == 0};
bc = {psi[0] == 0, u[0] == 0,
psi'[h] == -1 - 2 Pi e \[Alpha] \[Delta] Cos[2 Pi x],
psi[h] == q - 1};
px = NDSolveValue[
Join[eq,
bc], (D[(1 + n (We u[y] - 1)) u[y], y] - m psi'[y]) /. {y ->
h/2}, {y, 0, h}]; px]

Plot[Evaluate[Table[p[eps, k], {eps, {.1, .3, .5, .6}}]], {k, -2, 2},
PlotLegends ->
Table[Row[{"\[Epsilon] = ", eps}], {eps, {0.1, 0.3, 0.5, 0.6}}],
Frame -> True, FrameLabel -> { "Q","\[CapitalDelta]p"}] // Quiet


Update 1. To plot stream lines as shown in Figure 9, we can use code

e = .7; q = .9; M = 0.5;
\[Alpha] = 0.92;
n = 0.99;
\[Theta] = 0.9; m = M^2 Cos[\[Theta]]^2;
We = 0.0001;
\[Delta] = 0.0006;
p[xe_, ye_] := Module[{x = xe, yc = ye}, h = 1 + e Cos[2 \[Pi] x];
eq = {psi''[y] == u[y],
D[(1 + n (We u[y] - 1)) u[y], y, y] - m u[y] == 0};
bc = {psi[0] == 0, u[0] == 0,
psi'[h] == -1 - 2 Pi e \[Alpha] \[Delta] Cos[2 Pi x],
psi[h] == q - 1};
Ps = NDSolveValue[Join[eq, bc], psi, {y, 0, h},
Method -> {"StiffnessSwitching",
Method -> {"ExplicitRungeKutta", Automatic}}, AccuracyGoal -> 5,
PrecisionGoal -> 4]; Ps[yc]]

ContourPlot[p[xe, ye], {xe, -.3, .3}, {ye, 0., 1 + e}, Contours -> 30,
ColorFunction -> "Rainbow", PlotLegends -> Automatic]


Update 2. We can also solve this problem analytically with using series $$\psi =\psi_0+We \psi_1+...$$ as in a paper cited. We have

Clear["Global*"]

sol0 = DSolve[{D[psi[y], {y, 4}] - k D[psi[y], {y, 2}] == 0,
psi[0] == 0, psi''[0] == 0, psi'[h] == f0[x], psi[h] == F0},
psi[y], y];

i = D[(D[psi[y] /. sol0[[1]], {y, 2}])^2, {y, 2}] // Simplify;

sol1 = DSolve[{D[psi[y], {y, 4}] - k D[psi[y], {y, 2}] == i n/(n - 1),
psi[0] == 0, psi''[0] == 0, psi'[h] == 0, psi[h] == F1}, psi[y],
y];
f0[x_] := -1 - 2 Pi e a d Cos[2 Pi x]; h = 1 + e Cos[2 Pi x]; k =
M^2 Cos[tet]^2/(1 - n);

Psi0 = psi[y] /. sol0[[1]]; Psi1 = psi[y] /. sol1[[1]];


Visualization (Figure 9)

{e = .7, a = 0.92, tet = 0.9, M = .5, We = .0005, F0 = .9 - 1, F1 = 0,
n = 0.99, d = .0006};
ContourPlot[Evaluate[Psi0 + We Psi1], {x, -.3, .3}, {y, 0, 1 + e},
Contours -> 20, ColorFunction -> "Rainbow", PlotLegends -> Automatic]


• thank you.....can you please make this graph I attached.
– umar
Nov 21, 2021 at 15:08
• Yes, we can make any graph. But can you show expression for pressure as a function of $\psi$? Nov 21, 2021 at 17:56
• kindly help me.
– umar
Dec 24, 2021 at 11:18
• @umar I got the paper, thank you. Actually they solved this equation in a form of series by $We$ powers. I can reproduce solution. Dec 26, 2021 at 15:53

Here is your code with typos corrected. But there remains a big mistake. You give 2 ODE', namely w''[y] D[1 + n We (w''[y] - 1), {y, 2}] - M^2 Cos[θ]^2 , w''[y] == 0  instead of one. Presumably  w''[y] == 0 is wrong.:

M == 0.5
α == 0.2
n == 0.5
θ == 0.5
We == 0.02
ϵ == 0.1
δ == 0.0006
h == 1 + ϵ  Cos[2 π x]
F == \!$$\*SubsuperscriptBox[\(∫$$, $$0$$, $$h$$]$$\*SubscriptBox[\(∂$$, $$y$$]w \[DifferentialD]y\)\)

sol = NDSolve[{w''[y] D[1 + n  We (w''[y] - 1), {y, 2}] -
M^2 Cos[θ]^2  ,
w''[y] == 0 ,
w''[0] == 0, w[0] == 0,
w'[h] == -1 - 2 π ϵ α δ Cos[2 π x],
w[h] == F}, w[y], {y, 0, h}]

`