I'm trying to solve numerically a non-linear problem in order to determine the velocity field ($U$) and the film thickness ($h$) of a non-Newtonian fluid over an inclined plane. The equations are,

$$\Biggl(I+(1-I)\Biggl[1+\Biggl(L\frac{dU(y)}{dy}\Biggl)^2\Biggl]^{(n- 1)/2}\Biggl)\frac{dU(y)}{dy}= -y\,$$

with the boundary conditions,

$$U(h)=0\,$$ $$\int_0^h U(y)dy = 1\,$$

These equation are dimensionless with domain of $y \in [0,h]$. The parameters $I,L $ and $n$ are known. I tried to overcome the difficulty of an unknown domain choosing the transformation $z=\dfrac{y}{h}$. With this transformation the domain $y \in [0,h]$ becomes $z \in [0,1]$. Using $\dfrac{dz}{dy}=\dfrac{1}{h}$ the equations become,

$$\Biggl(I+(1-I)\Biggl[1+\Biggl(\frac{L}{h}\frac{dU(z)}{dz}\Biggl)^2\Biggl]^{(n-1)/2}\Biggl)\frac{1}{h}\frac{dU(z)}{dz}= -zh\,$$

with the boundary conditions,

$$U(1)=0\,$$ $$\int_0^1 U(z)hdz = 1\,$$

I suppose that (in this new system) $h$ becomes a unknown parameter that I must find together with $U(z)$. Therefore; I'm trying to solve this system numerically with NDSolve. I'm following a hint that I found in the Wolfram documentation in how to handle unknown parameters using NDSolve in this way (sadly without any success),


diff = (i + (1 - i)*(1 + ((L/h[z])*(f'[z]))^2)^((n - 1)/2))*(1/h[z])*
    f'[z] == -(z*h[z])

sol = NDSolve[{diff, h'[z] == 0, f[1] == 0, 
      WhenEvent[NIntegrate[f[z]*h[z],z, 0, 1}] == 1,"StopIntegration"]}, {f, h}, {z, 0, 1}]

I believe that my problem is how to implement the integral boundary properly using NDSolve. Any help or suggestion will be much appreciated. Thank you for your attention.

PS: Reference to the original article that I found this problem.


This problem can be solved as follows. First, solve diff with h[z] treated as a parameter and the constants as listed in the queston to obtain an explicit expression for f'[z].

diff /. {h[z] -> h, n -> 1/2, i -> 1/2, L -> 2/5}
(* (f'[z] (1/2 + 1/(2 (1 + (4 f'[z]^2)/(25 h^2))^(1/4))))/h == -h z *)
Simplify[(h # - i f'[z]) & /@ %];
solf = Equal @@ Simplify[Solve[#^4 & /@ %, f'[z], Reals], z > 0][[2, 1]]
(* (f'[z] == Root[100 h^10 z^4 + 200 h^8 z^3 #1 + (150 h^6 z^2 + 16 h^8 z^4) #1^2 + 
   (50 h^4 z + 32 h^6 z^3) #1^3 + 24 h^4 z^2 #1^4 + 8 h^2 z #1^5 + #1^6 &, 2] *)

At this point, the first order ODE could be solved for f as a function of h by means of ParametricNDSolve, the solution integrated over {z, 0, 1}, and FindRoot used to vary h until the integral equals 1/h. However, it is easier to add an auxiliary dependent variable, g[z],

g'[z] == h f[z]

which is equal to the integral over h f[z], if g[0] == 0.

 solp = ParametricNDSolveValue[{solf, g'[z] == h f[z], f[1] == 0, g[0] == 0}, 
    {g, f}, {z, 0, 1}, {h}];

The integral now can be plotted as a function of h.

Plot[First[solp[h]][1], {h, .1, 2}, ImageSize -> Large, AxesLabel -> {h, g[1]}, 
    LabelStyle -> {15, Black, Bold}]

enter image description here

The desired value of h is determined by

solh = h /. FindRoot[First[solp[h]][1] == 1, {h, 1}, Evaluated -> False]
(* 1.43127 *)

and f[z] plotted.

Plot[Last[solp[solh]][z], {z, 0, 1}, ImageSize -> Large, AxesLabel -> {z, f}, 
    LabelStyle -> {15, Black, Bold}]

enter image description here

Addendum: Alternative Solution Technique

This problem also can be solved using the method described here. To do so, replace h by h[z] and include h'[z] == 0.

NDSolveValue[{solf /. h -> h[z], g'[z] == h[z] f[z], h'[z] == 0, g[0] == 0, f[1] == 0, 
    g[1] == 1} /. h[z] -> Max[h[z], .1], {f[z], h[1]}, {z, 0, 1}, Method -> {"Shooting", 
    "StartingInitialConditions" -> {f[1] == 0, g[1] == 1, h[1] == 1}}]

Note that /. h[z] -> Max[h[z], .1] also is necessary to prevent the search for h[1] conducted by NDSolve from trying h[1] == 0, where diff becomes singular. I do not understand why this additional precaution should be necessary. In any case, the code just given reproduces the results in the first part of this answer.

Addendum: More General Solution

The two preceding solutions require that n be a rational number, so that diff can be transformed into a polynomial, preferably of low order, in f'[z]. NDSolve also can integrate ODEs that cannot be solved explicitly for the leading derivative by automatically using the {"EquationSimplification" -> "Residual"} Method to convert the ODE into a DAE. Unfortunately, the DAE solver cannot handle boundary value problems, forcing the following approach: Solve the original equation in the question, for which f is a function of y and the integral constraint is over the range {y, 0, h}, with h the value of y for f[y] == 0. Then vary f[0] until the integral equals 1.

s = ParametricNDSolveValue[{diff /. {n -> 1/2, i -> 1/2, L -> 2/5}, g'[y] == f[y], 
    f[0] == f0, g[0] == 0, WhenEvent[f[y] == 0, h = y; "StopIntegration"]}, 
    {g, f}, {y, 0, 10}, {f0}];
sf = FindRoot[First[s[f00]][h] - 1, {f00, 1}, Evaluated -> False];
Through[(s[f00] /. %)[h]]
(* {1., -1.73472*10^-17} *)
(* 1.43127 *)
Plot[Evaluate@Through[(s[f00] /. sf)[y]], {y, 0, h}, ImageSize -> Large, 
    AxesLabel -> {y, "g,f"}, LabelStyle -> {15, Black, Bold}]

enter image description here

With y == h z taken into account, the y curve is indistinguishable to the eye from the second plot above. Note, however, that h computed here is about 0.7% smaller than that computed in the preceding approaches. This small discrepancy perhaps is due to the difference in integration methods used for ODEs and DAEs.

To test this approach for irrational n and strong nonlinearity, repeat the computation with {n -> Pi, i -> 1/2, L -> 2} to yield h = 1.80483 and

enter image description here

  • $\begingroup$ Thank you so much for such a clear and helpful answer. If you allow me one more question, I would like to know if your code segment (solg = Equal @@ Simplify[Solve[#^4 & /@ %, f'[z], Reals], z > 0][[2, 1]]) took a few hours to run. Again, thank you for your answer and your time. $\endgroup$ – brunopc1 Jan 18 '20 at 1:12
  • $\begingroup$ @brunopc1 Just a few seconds. To be sure that I had copied my code correctly to my answer, I copied it from my answer to a new notebook and ran it. No problem. By the way, I plan to modify my answer a bit tonight. Results will not change, but the explanation my be more understandable. $\endgroup$ – bbgodfrey Jan 18 '20 at 2:46
  • $\begingroup$ @brunopc1 Just added an additional approach that works for n that are irrational as well as rational numbers. $\endgroup$ – bbgodfrey Jan 18 '20 at 21:05
  • $\begingroup$ Thank you so much for all your help. About the time problem that I mentioned before, it was a small mistake in my code. Everything is running properly now. $\endgroup$ – brunopc1 Jan 19 '20 at 13:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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