Heat convection differential equations from 1952 - Mathematica "fails to converge" - Mathematica Stack Exchange most recent 30 from mathematica.stackexchange.com 2019-09-21T18:55:48Z https://mathematica.stackexchange.com/feeds/question/104170 https://creativecommons.org/licenses/by-sa/4.0/rdf https://mathematica.stackexchange.com/q/104170 16 Heat convection differential equations from 1952 - Mathematica "fails to converge" dearN https://mathematica.stackexchange.com/users/204 2016-01-15T18:51:02Z 2017-09-20T02:05:44Z <p>I am trying to solve a fundamental problem in analytical convective heat transfer: laminar free convection flow and heat transfer from a flat plate parallel to the direction of the generating body force.</p> <p><strong>Brief History of the problem</strong></p> <p>Effectively: a flat plate is vertical and parallel to the direction of gravity vector. The plate is hot and the ambient is not. Heat transfer occurs from the plate to the ambient through natural convection due to density stratification. </p> <p><a href="https://en.wikipedia.org/wiki/Simon_Ostrach" rel="nofollow noreferrer">Simon Ostrach</a>, a distinguished scientist in the field of microgravity science <a href="https://dl.dropboxusercontent.com/u/13223318/Pohlhausen1952a.pdf" rel="nofollow noreferrer">solved this problem through a coupled set of equations</a>. <strong>In Ostrach's work, these equations were solved by an <a href="https://www-03.ibm.com/ibm/history/exhibits/vintage/vintage_4506VV2198.html" rel="nofollow noreferrer">IBM Card Programmed Electronic Calculator</a></strong></p> <p>$$F''' + 3 FF'' - 2 (F')^2 + H = 0$$ $$H'' + 2 \text{Pr} F H' = 0$$</p> <p>The Boundary conditions are: $$F'(0) = F(0) = 0$$ $$H(0) = 1$$ $$F'(\infty) = H(\infty) = 0$$</p> <p>Here, $F$ provides the hydrodynamic solution while $H$ provides the thermal solution with Pr being the Prandtl number which is a property of the fluid that the plate is "immersed" in.</p> <p><strong>My Mathematica code ... it runs selectively</strong></p> <pre><code>Clear[max, Pr, T, f, η, p]; max = 50; Pr = 0.72; pohl = NDSolve[{f'''[η] + 3 f[η] f''[η] - 2 (f'[η])^2 + T[η] == 0, T''[η] + 2 Pr f[η] T'[η] == 0, f == f' == 0, f'[max] == 0, T == 1, T[max] == 0}, {f, T}, {η, max}] p4 = Plot[{Evaluate[f'[η] /. pohl]}, {η, 0, max}, PlotRange -&gt; All, PlotLabel -&gt; Style[Framed["Hydrodynamic development is depicted on this plot"], 10, Blue, Background -&gt; Lighter[Yellow]], ImageSize -&gt; Large, BaseStyle -&gt; {FontWeight -&gt; "Bold", FontSize -&gt; 18}, AxesLabel -&gt; {"η", "f'[η]"}, PlotLegends -&gt; "Expressions"] </code></pre> <p>For a Prandtl number of 0.72 (Air) I get a velocity profile ($F'$) as suggested by the Ostrach in his pivotal report. However, for, many Prandtl numbers, the following warning message is sometimes flashed and I get incorrect velocity profiles (negative velocities) per the publication. For instance try Pr=6.</p> <blockquote> <p>FindRoot::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations. >></p> <p>NDSolve::berr: The scaled boundary value residual error of 2.9035865095898766*^7 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found.</p> </blockquote> <p>I have experimented with the <code>LSODA</code> <a href="https://mathematica.stackexchange.com/q/11630/204">method because this system of diff eqs is stiff</a> and LSODA has proven to be a 'magic wand' in the past. What gives? How do I select a method for this problem? I wonder if this is a problem with the method of choice (or default method with no options) or my definition of the "free stream limit" $\infty$.</p> <p><strong>Pr=0.01</strong></p> <p><a href="https://i.stack.imgur.com/FfJZM.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/FfJZM.png" alt="Pr=0.01"></a></p> <p><strong>Pr=0.72</strong></p> <p><a href="https://i.stack.imgur.com/ShFmq.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/ShFmq.png" alt="Pr=0.72"></a></p> <p><strong>Pr=0.6 (what went wrong? Warning message was displayed too...)</strong> <a href="https://i.stack.imgur.com/2MxSH.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/2MxSH.png" alt="Pr=0.6"></a></p> https://mathematica.stackexchange.com/questions/104170/heat-convection-differential-equations-from-1952-mathematica-fails-to-converg/118911#118911 12 Answer by Michael E2 for Heat convection differential equations from 1952 - Mathematica "fails to converge" Michael E2 https://mathematica.stackexchange.com/users/4999 2016-06-20T19:34:54Z 2016-06-21T12:40:08Z <p>The problem is with the default starting initial conditions used by the shooting method in <code>NDSolve</code>. The shooting method is where <code>FindRoot</code> is being used internally, so the OP's error message is a strong hint that this is the problem. Getting convergence in a nonlinear system can depend greatly on the starting conditions. </p> <p>Having luckily solved the system for <code>Pr = 0.72</code>, we can use its initial conditions as starting values for <code>Pr = 0.6</code>. We hope that it will be suitably close. (If not, we could have tried solving for, say, <code>Pr = 0.66</code> and edged our way bit by bit to <code>0.6</code>, hoping that the dependence on <code>Pr</code> is continuous.)</p> <pre><code>Pr = 0.72; pohl72 = NDSolve[{f'''[η] + 3 f[η] f''[η] - 2 (f'[η])^2 + T[η] == 0, T''[η] + 2 Pr f[η] T'[η] == 0, f == f' == 0, f'[max] == 0, T == 1, T[max] == 0}, {f, T}, {η, max}]; Pr = 0.6; pohl = NDSolve[{f'''[η] + 3 f[η] f''[η] - 2 (f'[η])^2 + T[η] == 0, T''[η] + 2 Pr f[η] T'[η] == 0, f == f' == 0, f'[max] == 0, T == 1, T[max] == 0}, {f, T}, {η, max}, Method -&gt; {"Shooting", "StartingInitialConditions" -&gt; Thread[{f, f', f'', T, T'} == ({f, f', f'', T, T'} /. First@pohl72)]}] </code></pre> <p>Plot:</p> <pre><code>Plot[{Evaluate[f'[η] /. pohl]}, {η, 0, max}, PlotRange -&gt; All, PlotLabel -&gt; Style[Framed["Hydrodynamic development is depicted on this plot"], 10, Blue, Background -&gt; Lighter[Yellow]], ImageSize -&gt; Large, BaseStyle -&gt; {FontWeight -&gt; "Bold", FontSize -&gt; 18}, AxesLabel -&gt; {"η", "f'[η]"}, PlotLegends -&gt; "Expressions"] </code></pre> <p><img src="https://i.stack.imgur.com/qDhNT.png" alt="Mathematica graphics"></p> https://mathematica.stackexchange.com/questions/104170/heat-convection-differential-equations-from-1952-mathematica-fails-to-converg/118994#118994 2 Answer by george2079 for Heat convection differential equations from 1952 - Mathematica "fails to converge" george2079 https://mathematica.stackexchange.com/users/2079 2016-06-21T15:12:35Z 2016-06-21T15:12:35Z <p>The problem with your system of equation is that the "zero" far field solution is unstable, hence extremely sensitive to the initial conditions. Posing the problem as an initial value problem, with the "known" conditions from the successful solution:</p> <pre><code>max = 50; Pr = .72; a = 0.7172594734816521 b = -0.4344414944896132` pohl = NDSolve[{f'''[\[Eta]] + 3 f[\[Eta]] f''[\[Eta]] - 2 (f'[\[Eta]])^2 + T[\[Eta]] == 0, T''[\[Eta]] + 2 Pr f[\[Eta]] T'[\[Eta]] == 0, f == f' == 0, f'' == a, T == 1, T' == b}, {f, T}, {\[Eta], max}] </code></pre> <p><a href="https://i.stack.imgur.com/MFfi4.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/MFfi4.png" alt="enter image description here"></a></p> <p>now change the <code>T'</code> initial value by 1/2%:</p> <pre><code>b=.995b </code></pre> <p>and you see you get the essential solution but it eventually blows up:</p> <p><a href="https://i.stack.imgur.com/IAtkK.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/IAtkK.png" alt="enter image description here"></a></p> <p>so the shooting method needs an exceptionally good initial guess to work.</p>