Falkner – Skan equation is a third order non-linear ordinary differential equation which arises in the laminar boundary layer flow past wedge-like objects. (More details here). The equation reads,

where *m* is a constant representing the pressure gradient parameter. Our objective is to solve this differential equation for , for a given value of *m* using the boundary conditions,

Generally, initial value problems (IVP) are preferred for ODEs. In the case of IVP, We will be given where to start and which direction to proceed. Then we use the differential equation to progress in that direction in a step by step manner. But here we have a boundary value problem. Shooting method is used in situations where a boundary value problem has to be solved using initial value methods. The method is described as follows.

## Shooting method:

- Guess two values for .
- Solve the FS equation using RK4 method with initial conditions .
- Solve the FS equation using RK4 method with initial conditions .
- Find out the resulting boundary value from both these solutions.
- If the boundary value is different from the required value , find a better initial guess using Secant method.
- Solve the FS equation by RK4 method using the new initial value guess (obtained from secant method).
- Repeat the process until the required boundary value is obtained.

A Fortran program for solving the Falkner-Skan equation implementing the above algorithm is provided below. (Click on any part of the code and use right arrow key to scroll to right).

program falkner_skan ! this program solves the Falkner-Skan equation using shooting method and fourth order Runge-Kutta method. ! the FS equation is given by, ! f'''+((m+1)/2)ff''+ m*(1-f'^2) = 0 ! m,the pressure gradient parameter, as in u_inf=u0*x^m. ! boundary conditions: f(0)=0, f'(0)=0, f'(1)=1 ! the third order ode is converted into a system of three first order odes where the ! dependent variables are f, u (=f'), v (=u'=f''). ! we solve the initial value problem f(0)=0,f'(0)=0, f''(0)=guess and shoot for solutions such that f'(1)=1. implicit none integer :: i,j,np,niter real,allocatable,dimension(:) :: f,u,v,eta real :: m,mh,re,ue,f0,u0,v0,deta,x,eps,cnu,slope,unp,unp1,vnp,vnp1 real,allocatable,dimension(:) :: vvel,y,vvel1 real :: delta_star_temp,theta_temp,deltas,theta niter=1000 ! no. of iterations for the shooting method deta=0.001 ! non-dimensional spacing in the wall-normal direction np=8001 ! no. of points in the wall-normal direction eps=1.0e-16 ! error margin used in shooting method iteration m=0.0 ! falkner-skan pressure gradient parameter mh=0.5*(m+1.0) allocate(f(np),u(np),v(np),eta(np),vvel(np),vvel1(np),y(np)) ! form the grid in the wall-normal direction eta(1)=0.0 do j=2,np eta(j)=eta(j-1)+deta enddo ! initial values set 1 f0=0.0 u0=0.0 v0=0.335 call rk4_fs(f,u,v,deta,mh,m,np,f0,u0,v0) vnp=v(1) unp=u(np) ! initial values set 2 f0=0.0 u0=0.0 v0=0.33 call rk4_fs(f,u,v,deta,mh,m,np,f0,u0,v0) vnp1=v(1) unp1=u(np) ! loop for the shooting method ******************************************************************** do i=1,niter if (abs(unp1-1.0) .ge. eps) then slope=(vnp1-vnp)/(unp1-unp) v0=vnp1+slope*(1.0-unp1) unp=unp1 vnp=vnp1 call rk4_fs(f,u,v,deta,mh,m,np,f0,u0,v0) vnp1=v(1) unp1=u(np) else write(*,*) 'iteration converged' write(*,*) 'v(1)=',v(1) exit endif if(i .eq. niter) then write(*,*) 'maximum number of iterations reached' endif enddo ! end of shooting method ************************************************************************** do j=1,np write(38,*) eta(j),f(j),u(j),v(j) ! these are transformed variables - input for box method enddo ! calculation of normal velocity component profile x=0.1 u0=1.0 ue=u0*(x**m) cnu=15.0e-06 re=ue*x/cnu do j=1,np y(j)=eta(j)*sqrt(cnu*x/ue) ! distance from the wall in meters vvel(j)=-((m-1.0)*eta(j)*u(j)+(m+1.0)*f(j))*(0.5/sqrt(re)) ! vvel=v/u_inf write(48,*) eta(j),u(j),vvel(j),y(j),u(j)*ue,vvel(j)*sqrt(re) enddo ! calculation of boundary layer parameters theta_temp = 0.0 delta_star_temp = 0.0 do i = 1,(np-1) theta_temp = theta_temp + 0.5*(eta(i+1)-eta(i))*(u(i+1)*(1.0-u(i+1)) + u(i)*(1.0-u(i)) ) delta_star_temp = delta_star_temp + 0.5*(eta(i+1)-eta(i))*((1.0-u(i+1))+(1.0-u(i))) enddo deltas=delta_star_temp*x/sqrt(re) theta=theta_temp*x/sqrt(re) write(*,*) '************************************************************************' write(*,*) 'Solution of the Falkner-Skan equation' write(*,*) '************************************************************************' write(*,'(a,f8.5)') 'pressure gradient parameter, m=',m write(*,'(a,f8.5,a)') 'distance from the leading edge, x=',x,' m' write(*,'(a,f8.2)') 'Re = u_inf*x/cnu = ',re write(*,'(a,f8.5,a)') 'displacement thickness, delta*=',deltas,' m' write(*,'(a,f8.5,a)') 'momentum thickness, theta=',theta,' m' write(*,*) '************************************************************************' call system('gnuplot -p FS.plt') return end program falkner_skan subroutine rk4_fs(f,u,v,deta,mh,m,np,f0,u0,v0) ! subroutine for solving a system of three first order odes with fourth order Runge-Kutta method ! f0,u0,v0 are the initial values ! f,u,v arrays give the solution ! initial values are propagated to np steps using step spacing deta. integer :: j real :: kf(4),ku(4),kv(4) real,intent(in) :: m,mh,deta,f0,u0,v0 integer,intent(in) :: np real :: f(np),u(np),v(np) f(1)=f0 u(1)=u0 v(1)=v0 do j=2,np ku(1)=v(j-1) kf(1)=u(j-1) kv(1)=-mh*f(j-1)*v(j-1)-m*(1.0-u(j-1)*u(j-1)) ku(2)=v(j-1)+0.5*kv(1)*deta kf(2)=u(j-1)+0.5*ku(1)*deta kv(2)=-mh*(f(j-1)+0.5*kf(1)*deta)*(v(j-1)+0.5*kv(1)*deta)-m*(1.0-(u(j-1)+0.5*ku(1)*deta)**2) ku(3)=v(j-1)+0.5*kv(2)*deta kf(3)=u(j-1)+0.5*ku(2)*deta kv(3)=-mh*(f(j-1)+0.5*kf(2)*deta)*(v(j-1)+0.5*kv(2)*deta)-m*(1.0-(u(j-1)+0.5*ku(2)*deta)**2) ku(4)=v(j-1)+kv(3)*deta kf(4)=u(j-1)+ku(3)*deta kv(4)=-mh*(f(j-1)+kf(3)*deta)*(v(j-1)+kv(3)*deta)-m*(1.0-(u(j-1)+ku(3)*deta)**2) f(j)=f(j-1)+(1.0/6.0)*deta*(kf(1)+2.0*kf(2)+2.0*kf(3)+kf(4)) u(j)=u(j-1)+(1.0/6.0)*deta*(ku(1)+2.0*ku(2)+2.0*ku(3)+ku(4)) v(j)=v(j-1)+(1.0/6.0)*deta*(kv(1)+2.0*kv(2)+2.0*kv(3)+kv(4)) enddo end subroutine rk4_fs

Line 114 in the above code calls the *gnuplot* program FS.plt to plot the solution. The most relevant plots for fluid dynamics in this problem are the streamwise velocity profile and the wall-normal velocity profile given by,

Sample results for the case m=0 (Blasius boundary layer) are shown below.

the solution for m=1 is not true and lead to overflow. please help!

For the m = 1 case, run the program with the following initial conditions:

Intitial values set 1:

F0=0.0

U0=0.0

V0=1.1

Intitial values set 2:

F0=0.0

U0=0.0

V0=1.3

yes. its true. thank you so much.

Hello,

I notice that, once I change NP (everything else remains unchanged, m = 1 in my case), I receive ‘NAN’ error. I tried to increase NP to 72001 instead of original 8001 it did not work. I also tried 9001, it ‘partially’ worked, but ‘V(1)’ changed significantly. Would you know I shall I fix it?

I also have a quick question, if you do not mind. If I need to change m to some other value, how could I quickly set two reasonable initial values for V0?

Thanks for the help!

Hi,

I have noticed such behavior as well. From what I understand, it is related to the approximate boundary layer thickness used for solving the problem. For example, by simply increasing NP, we are increasing the domain of solution. But if this domain length (NP*DETA) far exceeds 8*delta (where delta is the usual scale of boundary layer thickness that can be estimated from laminar boundary layer theory), the numerical approach seem to fail and blows up or falls short of the freestream value. It is theoretically true that this domain length should not affect the final solution since the freestream boundary condition is applied at infinity. But the numerical approach seems to have a problem with that.

The value of the initial guess for the second derivative is by trial and error. We normally expect the shooting method to work with two good positive values (since this corresponds to the wall shear stress). And you can estimate a good guess from the known value for the m=0 case problem. Note that changing ‘m’ corresponds to changing the acceleration in the flow which helps to determine whether the wall shear stress increases or decreases. But, this equation is non-linear and its behavior drastically changes under different parametric spaces. It is not guaranteed that a particular guess will work correctly and the freestream velocity will be smoothly approached by the solution. There are other, better, numerical approaches to deal with tricky situations. But, since the shooting method is simple to understand and apply, it is presented here. You can check out the book ‘Fundamentals of engineering numerical analysis’ by Parviz Moin which discusses some of these methods.

Thank you for the reply. I will think about it.