Project 2 Numerical Methods I: MIE334 H1S
Department of Mechanical and Industrial Engineering University of Toronto
Winter 2021
Introduction:
Project 2
MIE334: Numerical Methods I Due Date: April 16, 2021, 11:59 pm
In fluid mechanics, a boundary layer is a thin layer of a flowing gas or liquid that comes into contact with a surface such as an airplane or the inside of a pipe. Shearing forces act on the fluid in the boundary layer. If the fluid is in contact with the surface, a range of velocities occurs within the thin boundary layer, varying from a maximum to zero. The boundary layers of an aircraft wing are thinner at the leading edge and thicker at the trailing edge. In such boundary layers, the flow is usually laminar at the leading or upstream end and turbulent at the trailing or downstream end, as shown in Figure 1.
Figure 1. Boundary layer of an aircraft wing (from: tec-science.com).
As long as the contour of a body has smooth transitions and the fluid viscosity is large enough to resist the inertial forces, a flow can follow the contour. A typical hydrodynamic boundary layer develops around the body, the thickness of which is largely influenced by the viscosity. (see Figure 2)
Figure 2. Hydrodynamic boundary layer smooth development along the profile of a wing.
In 1905 Ludwig Prandtl, an engineer by profession and therefore motivated to find realistic fields near bodies of various shapes, first hypothesized about the boundary layers, and relation between the fluid viscosity and the thickness of these boundary layers. The equations of fluid motion within the boundary layer can be simplified because of the layer’s thinness, and exact or approximate solutions can be obtained in many cases. In addition, boundary-layer phenomena provide explanations for the lift and drag characteristics of bodies of various shapes in high Reynolds number flows, including turbulent flows.
1
A fundamental assumption of boundary-layer theory is that the layer is thin compared to other length scales such as the length of the surface or its local radius of curvature. Across this thin layer, which can exist only in high Reynolds number flows, the velocity varies rapidly enough for viscous effects to be dominant. This is shown in Figure 3, where the boundary-layer thickness is greatly exaggerated. (On a typical airplane wing, which may have a chord (width) of several meters, the boundary-layer thickness is of order one centimeter.) However, thin viscous layers exist not only next to solid walls but also in the form of jets, wakes, and shear layers if the Reynolds number is sufficiently high. So, to be specific, we shall first consider the boundary layer contiguous to a solid surface, adopting a curvilinear coordinate system that conforms to the surface where x increases along the surface and y increases normal to it. Here the surface may be curved but the radius of curvature of the surface must be much larger than the boundary-layer thickness. We shall refer to the solution of the irrotational flow outside the boundary layer as the outer problem and that of the boundary-layer flow as the inner problem.
Figure 3. A boundary layer forms when a viscous fluid moves over a solid surface. Only the boundary layer on the top surface of the foil is depicted in the figure and its thickness, 𝛿, is greatly exaggerated. Here, 𝑈∞ is the
oncoming free-stream velocity and 𝑈𝑒 is the velocity at the edge of the boundary layer. (From: Kundu, P. K., Cohen, I. M., & Hu, H. H. (2004). Fluid mechanics.)
Problem Statement
A researcher wants to determine how the boundary layer thickness 𝛿 changes along the solid surface, which is crucial to determine the drag and lift force applied on the solid body by the fluid flow. Calculating boundary layer thickness at each point along the 𝑥 direction involves velocity profile measurement in 𝑦 direction, which can be rather expensive given the length and shape of the solid surface profile. To simplify the problem, they first try find a cheaper and simpler solution to find 𝛿 for a laminar boundary layer on a long-enough flat plate, as shown in Figure 4.
Figure 4. Development of a boundary layer over a flat plate. (From: aerospaceengineeringblog.com) 2
An experimental setup for measurement of the velocity profile at a single location 𝑥 is also shown in Figure 4. (See Martin et al. “Micro-PIV Measurements and Challenges Towards Characterizing Active Flow Control Actuators for Rotor Applications.” (2009). for more information).
Figure 5. Schematic of the experimental setup for flat plate boundary layer measurements.
Using the above experimental setup, the velocity profile depicted below Figure 6 was successfully measured by
the researcher at an arbitrary position on the flat plate (data points are given in the table).
𝒚/𝜹
𝒖/𝒖𝟎
0.00000
0.00000
0.04613
0.07241
0.12175
0.19009
0.21479
0.33103
0.32132
0.48357
0.43915
0.63640
0.56685
0.77733
0.70339
0.89341
0.84798
0.97162
1.00000
1.00000
Figure 6. Schematic Velocity profile an arbitrary position on the flat plate
3
Since the researcher is not familiar with any numerical methods or boundary layer theory, they have asked you to help them with a method to find the boundary layer thickness as a function of position, i.e. 𝛿(𝑥), with this single measurement without need to repeat the same experiment for every position along the boundary layer.
Solution Approach
One of the earliest and, until recently, most widely used approximate methods for the solution of the boundary layer equation is that developed by Pohlhausen. This method is based on the momentum equation of von Karman (a student of Prandtl), which is obtained by integrating the boundary layer equation across the layer. (The theory behind this method is beyond the scope of this project, but you can find more information here). In the case of steady flow over an impermeable flat plate without any pressure gradient, the momentum equation reduces to:
𝜏𝑤 = 𝑑 𝜃
𝜌𝑢2 𝑑𝑥 0
(1)
In which 𝜏𝑤 is the wall shear stress and: 𝜏𝑤=𝜇(𝜕𝑢)|
In the Pohlhausen’s method, which is an approximate method, a form for the velocity profile 𝑢(𝑥, 𝑦) is sought which satisfies the momentum equation and some of the boundary conditions It is hoped that this form will approximate the exact profile, which satisfies all the conditions. In this method, the following parameters are defined as:
𝑢 =𝑓(𝜂)with𝜂= 𝑦 𝑢0 𝛿(𝑥)
It can be shown that 𝑓(𝜂) in this problem is independent of 𝑥 (self-similarity). By substituting the defined parameters, the reduced momentum equation takes the form of:
𝜕𝑦 𝑦=0
and𝜃=∫1(𝑢)(1−𝑢)𝑑(𝑦 )
0 𝑢0 𝑢0 𝛿(𝑥) (2)
𝜈′𝑑1
𝑢 𝛿 𝑓 (0) = 𝑑𝑥 {𝛿(𝑥) ∫ 𝑓 × (1 − 𝑓)𝑑𝜂}
00
(3)
Having known the non-dimensional velocity profile, 𝑓(𝜂), either by guessing and enforcing boundary conditions (Pohlhausen’s method), or through experimental measurement of the velocity profile at one single position across the boundary layer, let us to calculate the boundary layer thickness 𝛿(𝑥) as a function of the position along the flat plate. Data points related to the measurement of 𝑓(𝜂) is given to you in Figure 6.
To obtain 𝛿(𝑥), you need to solve the equation 3. Therefore, first you need to fit a polynomial curve to the data- points given to you to have a smoother function for 𝑓(𝜂) for numerical calculations. You can use Lagrange polynomials for polynomial interpolation.
𝑛
𝑓(𝜂) = ∑𝑓𝐿 (𝜂) (4) 𝑖𝑖
𝑖=0
which is a linear combination the Lagrange basis polynomials 4
𝑛 𝜂−𝜂𝑗
𝐿𝑖(𝜂)=∏𝜂 −𝜂 (5)
𝑗=0𝑖 𝑗
After fitting the curve and data interpolation (see later on how many interpolation points to use), a high resolution relation between 𝑓(𝜂) and 𝜂 is obtained that can now for numerical differentiation and integration of 𝑓(𝜂) that you need to solve Equation 3. Ensure to use as many data as needed to increase the accuracy of calculations.
Subsequently, the following integral should be numerically calculated (see later for details) using the data you generated:
1
𝜃=∫ 𝑓×(1−𝑓)𝑑𝜂 (6) 0
After calculating the value of 𝜃 numerically, equation 3 reduces to:
𝜈 𝑓′(0)=𝜃 𝑑 𝛿(𝑥)
(7) Finally, if you solve equation 7, then rearrange to get the analytical answer in the following form:
𝑢0𝛿(𝑥) 𝑑𝑥
𝛿(𝑥)= 𝛽 𝑥 √𝑅𝑒𝑥
(8)
Where 𝑅𝑒𝑥 = 𝑢0𝑥. 𝜈
You need to calculate 𝛽.
Project Requirements: Online Report Submissions:
You should upload your report in the designated place on Quercus. You have to submit a single *.pdf file containing all parts.
1. First, you need to develop a function in MATLAB capable of fitting the Lagrange polynomial to the data-points. In this procedure, your function should first implement the coefficients 𝐿0,𝐿1,…,𝐿𝑘 for the polynomial as in Equation 5. Then the function should use the equation 4 to interpolate 20 new- data points (equally spaced). This way, instead of only 10 data-points, you will have 20 points to increase the accuracy of calculations.
Submission: Plot the results of polynomial curve ( 𝜂 𝑣𝑒𝑟𝑠𝑢𝑠 𝑓(𝜂) ) together with the original data- points. The original data-points should be shown as discrete symbols, while in the same figure, the polynomial curve should be plotted as a solid line. The figures should also be complete in terms of axes labels, legends, and title.
2. Estimate 𝑓′(0) using interpolated points. You can do this using second order forward difference at 𝜂 = 0.
Hint: Second order forward difference:
−3𝑓+4𝑓 −𝑓
𝑓′ = 𝑖 𝑖
𝑖+1 𝑖+2 + 𝑂(h2) 2h
5
3.
4.
Submission: Report your 𝑓′(0) .
Use Trapezoidal rule to calculate the integral of Equation 6 using interpolated points.
Submission: Report the equation that you used to calculate the integral and report the obtained estimate of 𝜃.
Derive 𝛿(𝑥) as a function of f’(0), 𝜃 and 𝑅𝑒𝑥 using Equation 7. Then derive an equation for beta 𝛽 by comparing your derived 𝛿(𝑥) from Equation 7 and Equation 8. Calculate the value of beta 𝛽 using the interpolated data.
𝜈 𝑢0𝛿(𝑥)
Hint: This is how you should start with Equation 7 to derive 𝛿(𝑥) :
𝑓′(0)=𝜃 𝑑 𝛿(𝑥) → 𝜈𝑓′(0) 𝑑𝑥=𝛿(𝑥)𝑑𝛿(𝑥) → ∫𝑥𝜈𝑓′(0) 𝑑𝑥′=∫𝛿(𝑥)𝛿(𝑥′)𝑑𝛿(𝑥′)→⋯
5.
Submission: Show your work to derive 𝛿(𝑥) and β. Then, report the obtained estimate of β. Also, plot 𝛿(𝑥) vs x for 0 < 𝑥 < 𝐿 where L is the length of the flat plate and:
L=1, 𝜈 = 0.001, 𝑎𝑛𝑑 𝑢0 = 1 .This shows the development of the boundary layer thickness over a plate.
Hint:
To plot 𝛿(𝑥) versus x, you consider x as x=linspace(0.0001,L,1000).
Recalculate beta using various numbers of interpolated data points. Consider number of data points
from 20 to 1000 with the increment size of 50 (20:50:1000). Calculate beta for each case and report
the error for each two subsequent n.
e.g.: Error(n=70) =|beta(n=70)−beta(n=20)| beta(n=70)
Submission: Plot beta vs n (n=20:50:1000). Comment on your results. Plot the error versus n (n=70:50:1000). Report the smallest n at which the error becomes less than 0.01 %.
𝑑𝑥 𝑢0𝜃 0 𝑢0𝜃 0
Electronic CODE Submissions:
You have to submit a single *.zip (or *.rar) file containing all of your m-files. The name of the *.zip (or *.rar) file should be:
lastname_studentID_project2.zip
Forexample: einstein_9999999999_project2.zip
1. A function with the name of “lagrange” that calculated the value of fitted curve for an arbitrary function using Lagrange method. The inputs of this function are original data point x, original data point f, and the new_x (the x values of the fitted curve equally spaced). The output of this function is the new f (the y value of the fitted curve).
2. A function with the name of “forwardifference” that calculates derivative of an arbitrary function at a particular point using the second order forward difference. The inputs of this function are data points X , data points Y, and i (the index that the user would like to calculate the derivative) and the output is the derivative of f at that point.
6
3. A function with the name of “integration” that calculates the integral of an arbitrary function using Trapezoidal rule. The inputs of this function are data points X, data points Y, the index of the lower limit, and the index of the upper limit of the integral. The output in the calculated integral.
4. The last function is the main function that should use all of the previous functions in it. The inputs of this function are original data points etha, original data points f, number of the interpolated data points that users would like to have in their fitted curve (ni), length of the plate (L), the velocity u0 and nu. The outputs of this function are:
• One figure that shows etha vs f(etha) for both original data points and the fitted curve. The original data points should be shown by discrete markers (e.g. circles). The figure should be complete in terms of axis labels, legends, and title.
• One figure that shows delta vs x . The figure should be complete in terms of axis labels, and title.
• Calculated f’(0), theta , beta, and delta.
To test your program, the TAs will call your main function in the Command Window in the following form:
eta = [0.00000, 0.04613, 0.12175, 0.21479, 0.32132, 0.43915,
0.56685, 0.70339, 0.84798, 1.00000];
f = [0.00000, 0.07241, 0.19009, 0.33103, 0.48357, 0.63640, 0.77733,
0.89341, 0.97162, 1.00000];
ni = 20;
L = 1;
u_0 = 1;
nu = 0.001;
[df, theta, beta, delta] =
einstein_9999999999_project2(eta,f,ni,L,u_0,nu)
where df is the 𝑓′(0), theta is the 𝜃, beta is the 𝛽, eta and f are the original data-points, and ni is the number of interpolated points L,u_0,nu are length of the flat plate, upstream velocity, and kinematic viscosity of the fluid, respectively. You should define your main function to make possible such an action in MATLAB. Once called, your program should print out the 𝑓′(0), 𝜃, 𝛽, a figure where 𝑓(𝜂) is plotted against 𝜂 (with the conditions required for Report Submission, Part 2), and finally a figure where 𝛿(𝑥) is plotted against 𝑥. It is strongly recommended that you check your program in the same manner to avoid any problems. It is your responsibility to make sure that your code works properly without any further adjustments by the TAs.
Remember not to capitalize your last name at any point. The commented header at the beginning of each file which describes your function must be completed.
7