Solving Blasius Equation with Fortran 77

[carpediem85]

hi I am trying to write a Fortran 77 code to solve the Blasius equation numerically:

blasius equation: F''' + (1/2)*F*F'' = 0

boundary conditions: f(0)=0.0, F'(0)=0, and limit of F'(eta) as eta approaches infinity is 1.0.

This differential equation represents the velocity profile for an incompressible and laminar flow over a flat plate.

So far I know this about how to do it:

1. I need to take this 3rd order ode and build asystem of of 1st order odes. I do not how many, probably 3 or 2, I am not sure at this points.

2. i need to use the shooting method to approximate the last boudnary condition and find f''(0) that I can use to continue to implement Runge-Kutta march.


If anyone please knows or has already solved this problem , i would really appreciate that you write to me or post here your advice. i am confused and cannt get tarted..once i know the differential equations and the initial conditions properly, i can write the runge kutta algorithm in no time.

thanks.

 

[Clausius2]

I need to take this 3rd order ode and build asystem of of 1st order odes. I do not how many, probably 3 or 2, I am not sure at this points.

What?? Don't you know how many first order equations can you write with a third order DE? Just three, man.

2. i need to use the shooting method to approximate the last boudnary condition and find f''(0) that I can use to continue to implement Runge-Kutta march.

Yeah, you're right. The only funny thing here is that the equation is non linear, and so the shooting method you should employ. In total, you should end up with three first order equations coming from the third order original equation plus three more first order equations which helps you to calculate the next value of the shooting parameter using a Newton-Raphson approximation.

 

[kasman]

Hallo! I also have a question. In the blasius equation, since you don't know the value of f(η) in every point, how can you compute the value for the point f(x+h/2,y+k1/2)
it's needed in order to compute the vlaue o k2 and so on... Any suggestions are welcome!

 

[gato_]

hi I am trying to write a Fortran 77 code to solve the Blasius equation numerically:

blasius equation: F''' + (1/2)*F*F'' = 0

boundary conditions: f(0)=0.0, F'(0)=0, and limit of F'(eta) as eta approaches infinity is 1.0.

This differential equation represents the velocity profile for an incompressible and laminar flow over a flat plate.

So far I know this about how to do it:

1. I need to take this 3rd order ode and build asystem of of 1st order odes. I do not how many, probably 3 or 2, I am not sure at this points.

2. i need to use the shooting method to approximate the last boudnary condition and find f''(0) that I can use to continue to implement Runge-Kutta march.If anyone please knows or has already solved this problem , i would really appreciate that you write to me or post here your advice. i am confused and cannt get tarted..once i know the differential equations and the initial conditions properly, i can write the runge kutta algorithm in no time.

thanks.

This is how it goes. First, consider that F is a function of an aditional parameter, representing F'' at 0, that is, $$F(x,\alpha)$$. The system, written for variables $$\vec{x}\equiv [F,F',F'',F_{,\alpha},F'_{,\alpha},F''_{\alpha}]$$ reads:
$$x'=[x_{2},x_{3},-1/2x_{1}{x_{3},x_{4},x_{5},-1/2(x_{1}x_{6}+x_{4}x_{3})]$$
with boundary conditions:
$$x_{0}\equiv [0,0,\alpha,0,0,0]$$
why including $$F_{,\alpha}$$? because at each step, your error is $$F(L,\alpha)-1$$, so that next alpha is calculated via a Newton Raphson method:
$$\Delta \alpha= -F(L,\alpha)/F_{,\alpha}(L,\alpha)$$

 

[kasman]

thanks gato! quite helpful!:)

 

[nanunath]

This is how it goes. First, consider that F is a function of an aditional parameter, representing F'' at 0, that is, $$F(x,\alpha)$$. The system, written for variables $$\vec{x}\equiv [F,F',F'',F_{,\alpha},F'_{,\alpha},F''_{\alpha}]$$ reads:
$$x'=[x_{2},x_{3},-1/2x_{1}{x_{3},x_{4},x_{5},-1/2(x_{1}x_{6}+x_{4}x_{3})]$$
with boundary conditions:
$$x_{0}\equiv [0,0,\alpha,0,0,0]$$
why including $$F_{,\alpha}$$? because at each step, your error is $$F(L,\alpha)-1$$, so that next alpha is calculated via a Newton Raphson method:
$$\Delta \alpha= -F(L,\alpha)/F_{,\alpha}(L,\alpha)$$

Hello...can anyone make me clear what exactly and how the 'alpha' factor here helps in calculating the value of f''(0) for next iteration..
I also don't get how to solve this equation by making the marching process (f''(0)) automatically ..
I know what Newton Raphson method is and used ..but not why F(x,alpha) is assumed ..and how to solve those 6 equation..
also what does [PLAIN]https://www.physicsforums.com/latex_images/24/2464600-4.png
can anyone please xplain to me that..
thanks in advance...

 

[gato_]

The $$\alpha$$ parameter is included just because you are solving an initial value problem of third order without the information required (in this case, the value $$F''(0)$$). The procedure involves making some guess ($$F''(0)=\alpha$$), then advancing the solution form x=0 to a reasonable length (with a Runge Kutta scheme, for instance), then taking $$F(L)-1\equiv F(\infty)-1$$ as the error of your calculation. You use a Newton Raphson method to iteratively minimize that error as much as you need. The variable $$F_{,\alpha}$$ is included so that you have the derivative of the error for the Newton Raphson. Hope this helps.

 

[jackmell]

Hello guys. Been following this problem. Not very good though with boundary-value problems other than trial-and error to find the initial conditions. However, I just learned that boundary-value problems can be solved directly in Mathematica:

http://reference.wolfram.com/mathematica/tutorial/NDSolveBVP.html

I realize though you want to write the code directly and not use Mathematica so perhaps just use this as a reference for future work on BVPs you may wish to solve. Also, the reference cited above does describe the method so perhaps you can review it to gain further insight into implementing it in FORTRAN.

Note below I use the method "shooting" and I only went as far as x=20. Not sure how far out I could really go but it does compute a value of f''(0) very close to the actual answer which is approx 0.332 (see http://en.wikipedia.org/wiki/Blasius_boundary_layer).

In[53]:=
sol = NDSolve[{Derivative[3][f][x] + 
      (1/2)*f[x]*Derivative[2][f][x] == 
     0, f[0] == 0, Derivative[1][f][0] == 
     0, Derivative[1][f][20] == 1}, f, x, 
   Method -> "Shooting"]
N[D[f[x] /. sol, x] /. x -> 20]
N[D[f[x] /. sol, {x, 2}] /. x -> 0]

Out[53]=
{{f -> InterpolatingFunction[]}}

Out[54]=
{0.9999999644878674}

Out[55]=
{0.33205733312444885}

 

[nanunath]

Thanks gato and jackwell,

But gato dear, I don't understand how you get the three more 1st order ODE's
And what does "F _,alpha" actually mean..I haven't seen this notation before...

I hope I was able to highlight what my problem is...

Btw..I solved blausius equation using automatic marching.using a very simple method (for a layman like me) as follows:
Enter a guess value of f''(0), increment it by .1 (just like that) in next iteration...and interpolate the line passing thru 2 these points to the point where f'(0)-1=0, which gives the next iteration value of f''(0)...5-6 iterations are sufficient...

But I want to know what the method ur talking of is all about..a layman like me feels a strong need to understand what u guys are talking about..please explain to me ..
thanks..bye