Blasius solution for a flat plate-MATLAB

Hello,I am here again with a new blog. which would discuss the topic of basic fluid mechanics. If we consider the flow over the flat plate & simplify momentum and continuity equations, there would be an x momentum equation(y momentum equation is inbuilt) and the continuity equation. For the flat plate, there would be no pressure gradient term in the x direction; which would turn out to be dp/dx=0. One should be aware that boundary layer equations are valid under i) Re>>1 ii) No boundary layer separation. So, we have two boundary layer equations, which are non-linear partial differential equations. To solve them it is required to write them in the form of ordinary differential equations, which can be done by doing stretching (similarity) transformation, so that equation would be a single variable function; lots of efforts were done by Blasius for this transformation. Blasius showed that velocity distribution downstream of the plate merges into a single curve, only if it has a solution like u/U = f(η); where U= freestream velocity, & η=y/δ(x). Here, we can write 1/δ(x) as g(x).

u du/dx + v du/dy= μ*d^2 u/dy^2

du/dx+dv/dy=0

As it, as shown above, we can write boundary layer equations in the function of η form and we can find the final solutions as d[f ''/f ']= -1/2*f*dη (above derivation can be found in the fluid mechanics book). If we integrate the f*dη, which would be similar to some function F.So, dF/dη =f, f '= F'',f '' = F '''. Here we know that f is like velocity, derivation of some special function(dF/dη =f) is f, so F has the physical meaning of like stream function. which would finally yield the equation as below.

F'''+1/2*F''*F=0.

After mathematical manipulation it is possible to derive the non-linear ordinary differential equation, we can solve it considering the boundary conditions below. (consider the x coordinate represents the length of the plate and y coordinate height)

@y=0,u=0 so η=0,f(0)=0,F'(0)=0

@y-->inf,u=U so η=1,f(1)=1,F'(1)=1,F''(1)=0

we have information about the velocity at the edge of the boundary layer becomes the free stream velocity but, we do not have information about how the velocity is changing from the base of the plate @y=0, we don't know what is f '(0) in other words F ''(0)=? we can consider the shooting method to solve this mystery. assume F=F1

we can consider,    dF1/dη = F2, 

                              dF2/dη = F3, 

                              dF3/dη= -0.5 * F1 *F3; three equations , which is the first order equations with some  conditions F1(0)=0,F1'(0)=0,F1'(η)=1, F1''(0)=?

so we know     F1(0)=0; 

                        F2(0)=F1'(0)=0; 

                        F3(0)=F2'(0)=F1''(0)=(?); so we can assume the value for F1''(0) such that F1'(η)=1

                                 = P (which is guessed value )  

 so, we need to guess some value of the P which we can use to find a condition such that as η--> inf, u=U, for first-time p was guessed as 1.0 for which we found that F1' achieves a value of 2(overshoot), which is beyond the targeted value of 1.

For the guessing value of 0.1 we found that as η--> inf, the value of F1P goes below(undershoot) the value of 1.0, which is targeted & which is shown as below in the figure.

 So, from the above results, it is required to find such a value that should lie between two guessed values 0.1(p1) & 1.0(p2), we require interpolation to find the value which gives the value of the F1'(η)=1.  with the guessing values of p1 and p2, the value for the F1'(η) were 0.4492 & 2.086, which is not our target value of 1 (it is better to mind it, F1'(η) represents the u/U). with some iterations, we can find the value of p such that F1'(η) becomes 1. it can be seen in the below graph.

In the above graph, the asterisk shows the convergence value of the p and the corresponding value of m which is 1. with iterations it was found that the value of p is 0.3320 which would shoot our target of F1'(η) = 1.  with the obtained value of p of 0.3320, it is possible to calculate the value of η for which u/U would be 1. if we plot the u/U vs η, we can find that for the value of η approximately 5.0 u reaches to freestream velocity. from the above equations, it is evident that η=y/δ(x)= y*g(x); considering the η = 5 & y=δ at the boundary layer edge, it is possible to find the simple relation for the flat plate δ/x=5 Re^-1/2, from this simple relation we can derive another variable. For reference MATLAB code with the Blasius function is also attached below.

link to download MATLAB file:

https://github.com/jayrajsolanki07/Matlab-CFD.git