ENGR20005
Numerical Methods in Engineering
Workshop 3
Part A: MATLAB Livescripts
3.1 The livescript ENGR20005 Workshop3p1.mlx runs through the use of functions in MATLAB.
(a) Read through the livescript and make sure you understand what each line of code does.
(b) Convert your .m files from workshop 2 to functions that output the root xr.
3.2 The livescript ENGR20005 Workshop3p2.mlx shows you how to analyse the function
f(x)=x3 −x (3.1) and implements the Bisection method to solve Eq. (3.1)
(a) Read through the livescript and make sure you understand what each of the commands are doing.
(b) Write a MATLAB code (a .m file) to analyse and use the Bisection method find the roots of the function
f(x)=e−x −x (3.2)
(c) Check that your answer is correct by finding the roots of the above equation using
the fzero() function in MATLAB.
Part B: Problems
3.3 Starting from the code snippets in the livescript ENGR20005 Workshop3p2.mlx, write a MATLAB code (a .m file) that uses the Method of False Position to solve Eq. (3.1). Check your answer with the answer given using the Bisection Method.
1
3.4 *When we turn on a tap we expect that water flows out. Such flows are the result of the balance between pressure and the resistance of motion (shear stresses) at the pipe wall.
To get an idea of how these forces act on the fluid: (a) Show that the force balance is given by
πD2
(p1−p2) 4 =τwLπD
(3.3)
(b) Determine the pressure gradient required to sustain the flow.
Therefore, to analyse these systems, we need a way to characterise the frictional resis-
tance at the wall. It can be shown that the shear stress, τw, is given by
τw = 81fρU2 (3.4)
with ρ being the fluid density, U the fluid velocity, and f is a friction factor (much like the drag coefficient that you might have seen in ESD2) that may be modelled using the Colebrook equation
ε1 ε2.51
g f;Re,D =√f+2.0log10 3.7D+Re√f =0 (3.5)
where ε is a parameter that describes the roughness of the pipe, and Re is a parameter called the Reynolds number, which is given by
Re = ρUD (3.6) μ
where μ is the dynamic viscosity, and D is the pipe diameter.
You don’t have to know what these symbols mean, that’s for next year. But for context, the Reynolds number characterises the balance of inertial forces and viscous forces. For pipe flow, it is used to characterise the different flow regimes.
Our goal for this section is to plot the Colebrook equation for a variety of Reynolds numbers. For simplicity sake, we’ll assume that the pipe is perfectly smooth, i.e. ε = 0, but you should be able to generalise your code to find the friction factor for an arbitrary roughness ratio, ε/D.
2
(c) To begin with assume a Reynolds number of Re = 104, plot g(f;Re = 104,ε/D = 0) to obtain a good guess of the root.
(d) Write a MATLAB function (.m file) that returns the root of Eq. (3.5).
(e) Repeat parts (c) and (d) for Reynolds numbers of Re = {105, 106, 107, 108} and
plot the results on a log-log scale.
(f) Instead of graphically estimating the root every time, the Haaland equation is an approximate explicit formula given by
√1 ≈−1.8log 6.9+ε/D1.11 (3.7) f 10Re3.7
Convince yourself that this is in fact a good guess of the root. Then modify your function in part (d) so that you don’t have to provide an initial guess.
However, the Colebrook equation is only valid for turbulent flow (Re 4000), for laminar flow (Re 2000), f takes the form
f = 64 (3.8) Re
(f) Generalise your function from part (c) to account for these two cases.
In the region 2000 Re 4000, a transitional regime exists, where both laminar and turbulent characteristics exist.
(g) Have your function return a null value in the transition region.
3