15 October 2013

Generating the Blasius boundary layer profile

I've been working on some boundary layer edge detection routines for Suzerain and wanted to use the Blasius boundary layer as a test case. The Blasius boundary layer produces a slick ODE that can only be solved numerically. Ganapol has a nice article on it and the more general Falkner—Skan equation. He also presents some nice tabulated data which I'd fit with an Akima, but I've found that (surprise surprise) the Akima natural end conditions disturb the second derivative of my fits sufficiently that it makes my test cases lousier than I'd like. Since Ganapol was the nicest tabulated Blasius data that I'd seen (meaning lots of digits, large similarity coordinate eta) and it didn't meet my needs, I set off to generate more (comparable number of digits, even larger similarity coordinate) using the GNU Scientific Library's ODEIV2 capabilities:

This code produces output matching Ganapol's everywhere to a minimum of eight digits. So though I'm reporting more than eight here, do not trust more than eight. The results are:

                   eta                       f                      fp                     fpp
0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  3.3205733621519629e-01
2.0000000000000001e-01  6.6409997145970585e-03  6.6407792096250848e-02  3.3198383711462798e-01
4.0000000000000002e-01  2.6559884017994733e-02  1.3276416076102221e-01  3.3146984420145503e-01
6.0000000000000009e-01  5.9734637498038382e-02  1.9893725242221899e-01  3.3007912757429625e-01
8.0000000000000004e-01  1.0610822083904069e-01  2.6470913872311741e-01  3.2738927014925229e-01
1.0000000000000000e+00  1.6557172578927976e-01  3.2978003124966698e-01  3.2300711668694282e-01
1.2000000000000002e+00  2.3794871728892703e-01  3.9377610443395650e-01  3.1658919106111544e-01
1.4000000000000001e+00  3.2298157382953580e-01  4.5626176470513380e-01  3.0786539179016786e-01
1.6000000000000001e+00  4.2032076550162573e-01  5.1675678442261486e-01  2.9666346145571931e-01
1.8000000000000000e+00  5.2951803774381023e-01  5.7475814388945690e-01  2.8293101725975589e-01
2.0000000000000000e+00  6.5002436993528900e-01  6.2976573650238588e-01  2.6675154569727827e-01
2.2000000000000002e+00  7.8119333701057347e-01  6.8131037723602750e-01  2.4835091319037131e-01
2.4000000000000004e+00  9.2229012563454160e-01  7.2898193506257503e-01  2.2809176068668360e-01
2.6000000000000001e+00  1.0725059767878351e+00  7.7245502114856535e-01  2.0645462679942550e-01
2.8000000000000003e+00  1.2309773023143091e+00  8.1150962319910902e-01  1.8400659386536070e-01
3.0000000000000000e+00  1.3968082308703464e+00  8.4604444365799358e-01  1.6136031954087834e-01
3.2000000000000002e+00  1.5690949600067619e+00  8.7608145517247749e-01  1.3912805557242577e-01
3.4000000000000004e+00  1.7469500939488432e+00  9.0176122143676474e-01  1.1787624608752340e-01
3.6000000000000001e+00  1.9295251697605302e+00  9.2332966588164878e-01  9.8086278784280084e-02
3.8000000000000003e+00  2.1160298171720697e+00  9.4111799672931340e-01  8.0125918139197061e-02
4.0000000000000000e+00  2.3057464184620726e+00  9.5551822981069434e-01  6.4234121091690577e-02
4.2000000000000002e+00  2.4980396627069701e+00  9.6695707375360584e-01  5.0519747486645686e-02
4.4000000000000004e+00  2.6923609375430866e+00  9.7587083213647785e-01  3.8972610853629498e-02
4.6000000000000005e+00  2.8882479900178160e+00  9.8268350076070032e-01  2.9483772011648642e-02
4.8000000000000007e+00  3.0853206551774823e+00  9.8778952620078742e-01  2.1871186347443238e-02
5.0000000000000000e+00  3.2832736651563144e+00  9.9154190016439347e-01  1.5906798685318153e-02
5.2000000000000002e+00  3.4818676115086635e+00  9.9424553535992810e-01  1.1341788968929135e-02
5.4000000000000004e+00  3.6809190628975084e+00  9.9615530396275309e-01  7.9276598147061551e-03
5.6000000000000005e+00  3.8802906776333757e+00  9.9747776821291501e-01  5.4319576799275147e-03
5.8000000000000007e+00  4.0798819392396917e+00  9.9837549365019085e-01  3.6484136667473462e-03
6.0000000000000000e+00  4.2796209225138426e+00  9.9897287243586053e-01  2.4020398437572996e-03
6.2000000000000002e+00  4.4794572972826767e+00  9.9936254171909067e-01  1.5501706906563828e-03
6.4000000000000004e+00  4.6793566154314172e+00  9.9961170171767588e-01  9.8061511700977232e-04
6.6000000000000005e+00  4.8792958110602429e+00  9.9976787021003244e-01  6.0804426478449074e-04
6.8000000000000007e+00  5.0792597724490030e+00  9.9986381903703803e-01  3.6956257014031428e-04
7.0000000000000000e+00  5.2792388110290958e+00  9.9992160414794995e-01  2.2016895527113460e-04
7.2000000000000002e+00  5.4792268473431696e+00  9.9995571727792665e-01  1.2856980723515674e-04
7.4000000000000004e+00  5.6792201473338295e+00  9.9997545768489782e-01  7.3592983389223794e-05
7.6000000000000005e+00  5.8792164658040509e+00  9.9998665513912843e-01  4.1290311113698544e-05
7.8000000000000007e+00  6.0792144810719346e+00  9.9999288116610119e-01  2.2707751402803766e-05
8.0000000000000000e+00  6.2792134313460615e+00  9.9999627453530060e-01  1.2240926243249920e-05
8.2000000000000011e+00  6.4792128866784848e+00  9.9999808745923757e-01  6.4679786108453658e-06
8.4000000000000004e+00  6.6792126094414330e+00  9.9999903687484326e-01  3.3499397531974408e-06
8.5999999999999996e+00  6.8792124710150580e+00  9.9999952424780447e-01  1.7006679885717248e-06
8.8000000000000007e+00  7.0792124032167214e+00  9.9999976948975089e-01  8.4628412130585306e-07
9.0000000000000000e+00  7.2792123706452294e+00  9.9999989045379900e-01  4.1278790156475232e-07
9.2000000000000011e+00  7.4792123552968919e+00  9.9999994893901312e-01  1.9735668380412326e-07
9.4000000000000004e+00  7.6792123482031105e+00  9.9999997665715090e-01  9.2489158677487331e-08
9.6000000000000014e+00  7.8792123449874136e+00  9.9999998953401714e-01  4.2485812620375726e-08
9.8000000000000007e+00  8.0792123435577192e+00  9.9999999539788897e-01  1.9129831304430804e-08
1.0000000000000000e+01  8.2792123429343132e+00  9.9999999801539019e-01  8.4429158651022967e-09
1.0200000000000001e+01  8.4792123426677222e+00  9.9999999916068538e-01  3.6524803913092198e-09
1.0400000000000000e+01  8.6792123425559158e+00  9.9999999965190545e-01  1.5488074704694815e-09
1.0600000000000001e+01  8.8792123425099305e+00  9.9999999985842569e-01  6.4375569794427733e-10
1.0800000000000001e+01  9.0792123424913829e+00  9.9999999994353506e-01  2.6227617757128940e-10
1.1000000000000000e+01  9.2792123424840440e+00  9.9999999997791622e-01  1.0473955128094568e-10
1.1200000000000001e+01  9.4792123424811994e+00  9.9999999999153044e-01  4.0999322474265197e-11
1.1400000000000000e+01  9.6792123424801169e+00  9.9999999999681477e-01  1.5731015212516414e-11
1.1600000000000001e+01  9.8792123424797147e+00  9.9999999999882527e-01  5.9163099772618016e-12
1.1800000000000001e+01  1.0079212342479567e+01  9.9999999999957512e-01  2.1810177063481187e-12
1.2000000000000000e+01  1.0279212342479513e+01  9.9999999999984923e-01  7.8810042675859475e-13
1.2200000000000001e+01  1.0479212342479494e+01  9.9999999999994749e-01  2.7913740190928823e-13
1.2400000000000000e+01  1.0679212342479486e+01  9.9999999999998201e-01  9.6910000670083129e-14
1.2600000000000001e+01  1.0879212342479486e+01  9.9999999999999389e-01  3.2978678842745192e-14
1.2800000000000001e+01  1.1079212342479485e+01  9.9999999999999789e-01  1.1000489193601443e-14
1.3000000000000000e+01  1.1279212342479484e+01  9.9999999999999922e-01  3.5967050786090774e-15
1.3200000000000001e+01  1.1479212342479485e+01  9.9999999999999967e-01  1.1526879056075572e-15
1.3400000000000000e+01  1.1679212342479484e+01  9.9999999999999978e-01  3.6210349563133137e-16
1.3600000000000000e+01  1.1879212342479484e+01  9.9999999999999978e-01  1.1149817624085130e-16

If you look at the last two rows of fp you'll see that double precision has most likely run out of steam by eta = 13.4. Using gsl_odeiv2_step_rk4imp or some other choice instead of gsl_odeiv2_step_rk8pd will tend to cause the solver to fail somewhere before here rather than to simply stop producing updated values in fp. Interestingly, gsl_odeiv2_step_rk8pd will run out to about eta = 38.6 producing fpp = 1.35e-155 before fpp becomes non-decreasing (indicating the physics is glaringly b0rked as opposed to just wholly b0rked).

1 comment:

Anonymous said...

Thanks for making the code and solution available.
It's been very useful for me (I will use your approach to generate boundary conditions for a CFD solver).

Subscribe Subscribe to The Return of Agent Zlerich