The complete Python program CDsphere.py used to plot the drag
coefficient in the example above is listed below. The program uses a
function cd_sphere which results from a curve fit to the data of
Evett and Liu [3]. In our setting we will use this
function for two purposes, namely to demonstrate how functions and modules are
implemented in Python and finally use these functions in the solution of the ODE in Equations (13) and (14).
# chapter1/programs_and_modules/CDsphere.py
from numpy import logspace, zeros
# Define the function cd_sphere
def cd_sphere(Re):
"This function computes the drag coefficient of a sphere as a function of the Reynolds number Re."
# Curve fitted after fig . A -56 in Evett and Liu: "Fluid Mechanics and Hydraulics"
from numpy import log10, array, polyval
if Re <= 0.0:
CD = 0.0
elif Re > 8.0e6:
CD = 0.2
elif Re > 0.0 and Re <= 0.5:
CD = 24.0/Re
elif Re > 0.5 and Re <= 100.0:
p = array([4.22, -14.05, 34.87, 0.658])
CD = polyval(p, 1.0/Re)
elif Re > 100.0 and Re <= 1.0e4:
p = array([-30.41, 43.72, -17.08, 2.41])
CD = polyval(p, 1.0/log10(Re))
elif Re > 1.0e4 and Re <= 3.35e5:
p = array([-0.1584, 2.031, -8.472, 11.932])
CD = polyval(p, log10(Re))
elif Re > 3.35e5 and Re <= 5.0e5:
x1 = log10(Re/4.5e5)
CD = 91.08*x1**4 + 0.0764
else:
p = array([-0.06338, 1.1905, -7.332, 14.93])
CD = polyval(p, log10(Re))
return CD
# Calculate drag coefficient
Npts = 500
Re = logspace(-1, 7, Npts, True, 10)
CD = zeros(Npts)
i_list = range(0, Npts-1)
for i in i_list:
CD[i] = cd_sphere(Re[i])
# Make plot
from matplotlib import pyplot
pyplot.plot(Re, CD, '-b')
font = {'size' : 16}
pyplot.rc('font', **font)
pyplot.yscale('log')
pyplot.xscale('log')
pyplot.xlabel('$Re$')
pyplot.ylabel('$C_D$')
pyplot.grid('on', 'both', 'both')
pyplot.savefig('example_sphere.png', transparent=True)
pyplot.show()
In the following, we will break up the program and explain the different parts. In the first code line,
from numpy import logspace, zeros
the functions logspace and zeros are imported from the package numpy. The numpy package (NumPy is an abbreviation for Numerical Python) enables the use of array objects. Using numpy a wide range of mathematical operations can be done directly on complete arrays, thereby removing the need for loops over array elements. This is commonly called vectorization and may cause a dramatic increase in computational speed of Python programs. The function logspace works on a logarithmic scale just as the function linspace works on a regular scale. The function zeros creates arrays of a certain size filled with zeros. Several comprehensive guides to the numpy package may be found at http://www.numpy.org.
In CDsphere.py a function cd_sphere was defined as follows:
def cd_sphere(Re):
"This function computes the drag coefficient of a sphere as a function of the Reynolds number Re."
# Curve fitted after fig . A -56 in Evett and Liu: "Fluid Mechanics and Hydraulics"
from numpy import log10, array, polyval
if Re <= 0.0:
CD = 0.0
elif Re > 8.0e6:
CD = 0.2
elif Re > 0.0 and Re <= 0.5:
CD = 24.0/Re
elif Re > 0.5 and Re <= 100.0:
p = array([4.22, -14.05, 34.87, 0.658])
CD = polyval(p, 1.0/Re)
elif Re > 100.0 and Re <= 1.0e4:
p = array([-30.41, 43.72, -17.08, 2.41])
CD = polyval(p, 1.0/log10(Re))
elif Re > 1.0e4 and Re <= 3.35e5:
p = array([-0.1584, 2.031, -8.472, 11.932])
CD = polyval(p, log10(Re))
elif Re > 3.35e5 and Re <= 5.0e5:
x1 = log10(Re/4.5e5)
CD = 91.08*x1**4 + 0.0764
else:
p = array([-0.06338, 1.1905, -7.332, 14.93])
CD = polyval(p, log10(Re))
return CD
The function takes Re as an argument and returns the value CD. All Python functions begin with def, followed by the function name, and then inside parentheses a comma-separated list of function arguments, ended with a colon. Here we have only one argument, Re. This argument acts as a standard variable inside the function. The statements to perform inside the function must be indented. At the end of a function it is common to use the return statement to return the value of the function.
Variables defined inside a function, such as p and x1 above, are local variables that cannot be accessed outside the function. Variables defined outside functions, in the "main program", are global variables and may be accessed anywhere, also inside functions.
Three more functions from the numpy package are imported in the function. They are not used outside the function and are therefore chosen to be imported only if the function is called from the main program. We refer to the documentation of NumPy for details about the different functions.
The function above contains an example of the use of the if-elif-else block. The block begins with if and a boolean expression. If the boolean expression evaluates to true the indented statements following the if statement are carried out. If not, the boolean expression following the elif is evaluated. If none of the conditions are evaluated to true the statements following the else are carried out.
In the code block
Npts = 500
Re = logspace(-1, 7, Npts, True, 10)
CD = zeros(Npts)
i_list = range(0, Npts-1)
for i in i_list:
CD[i] = cd_sphere(Re[i])
the function cd_sphere is called. First, the number of data points to be calculated are stored in the integer variable Npts. Using the logspace function imported earlier, Re is assigned an array object which has float elements with values ranging from \( 10^{-1} \) to \( 10^7 \). The values are uniformly distributed along a 10-logarithmic scale. CD is first defined as an array with Npts zero elements, using the zero function. Then, for each element in Re, the drag coefficient is calculated using our own defined function cd_sphere, in a for loop, which is explained in the following.
The function range is a built-in function that generates a list containing arithmetic progressions. The for i in i_list construct creates a loop over all elements in i_list. In each pass of the loop, the variable i refers to an element in the list, starting with i_list[0] (0 in this case) and ending with the last element i_list[Npts-1] (499 in this case). Note that element indices start at 0 in Python. After the colon comes a block of statements which does something useful with the current element; in this case, the return of the function call cd_sphere(Re[i]) is assigned to CD[i]. Each statement in the block must be indented.
Lastly, the drag coefficient is plotted and the figure generated:
from matplotlib import pyplot
pyplot.plot(Re, CD, '-b')
font = {'size' : 16}
pyplot.rc('font', **font)
pyplot.yscale('log')
pyplot.xscale('log')
pyplot.xlabel('$Re$')
pyplot.ylabel('$C_D$')
pyplot.grid('on', 'both', 'both')
pyplot.savefig('example_sphere.png', transparent=True)
pyplot.show()
To generate the plot, the package matplotlib is used. matplotlib is the standard package for curve plotting in Python. For simple plotting the matplotlib.pyplot interface provides a Matlab-like interface, which has been used here. For documentation and explanation of this package, we refer to http://www.matplotlib.org.
First, the curve is generated using the function plot, which takes the x-values and y-values as arguments (Re and CD in this case), as well as a string specifying the line style, like in Matlab. Then changes are made to the figure in order to make it more readable, very similarly to how it is done in Matlab. For instance, in this case it makes sense to use logarithmic scales. A png version of the figure is saved using the savefig function. Lastly, the figure is showed on the screen with the show function.
To change the font size the function rc is used. This function takes in the object font, which is a dictionary object. Roughly speaking, a dictionary is a list where the index can be a text (in lists the index must be an integer). It is best to think of a dictionary as an unordered set of key:value pairs, with the requirement that the keys are unique (within one dictionary). A pair of braces creates an empty dictionary: {}. Placing a comma-separated list of key:value pairs within the braces adds initial key:value pairs to the dictionary. In this case the dictionary font contains one key:value pair, namely 'size' : 16.
Descriptions and explanations of all functions available in pyplot may be found here.