4
$\begingroup$

streamlinespressure as a funtion o psi eq 14graph I need to make 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}]      
$\endgroup$
9
  • $\begingroup$ is this an ode or pde? what is $\psi$ a function of? just $y$ ? then why the book uses partial derivatives? $\endgroup$
    – Nasser
    Nov 20, 2021 at 17:55
  • $\begingroup$ 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. $\endgroup$ Nov 27, 2021 at 7:27
  • 1
    $\begingroup$ It is not clear in what coordinates stream plot in Figure 9 has been created. $\endgroup$ Dec 12, 2021 at 10:50
  • $\begingroup$ For completeness, please give a link to the paper you have taken a screenshot of (apparently, it's from World Scientific). $\endgroup$ Dec 13, 2021 at 10:35
  • $\begingroup$ worldscientific.com/doi/10.1142/S0219887819501391 $\endgroup$
    – umar
    Dec 13, 2021 at 10:43

2 Answers 2

6
$\begingroup$

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}]

Figure 1

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

Figure 2

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]

Figure 3

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]

Figure 4

$\endgroup$
23
  • $\begingroup$ thank you.....can you please make this graph I attached. $\endgroup$
    – umar
    Nov 21, 2021 at 15:08
  • $\begingroup$ Yes, we can make any graph. But can you show expression for pressure as a function of $\psi$? $\endgroup$ Nov 21, 2021 at 17:56
  • 1
    $\begingroup$ kindly help me. $\endgroup$
    – umar
    Dec 24, 2021 at 11:18
  • 1
    $\begingroup$ @umar Could you open this paper for me? I have no access to this page worldscientific.com/doi/epdf/10.1142/S0219887819501391 $\endgroup$ Dec 25, 2021 at 5:04
  • 1
    $\begingroup$ @umar I got the paper, thank you. Actually they solved this equation in a form of series by $We$ powers. I can reproduce solution. $\endgroup$ Dec 26, 2021 at 15:53
1
$\begingroup$

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}]

  
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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