Engineering Computation in Dr. Kyle Horne Department of Mechanical - - PowerPoint PPT Presentation

Engineering Computation in

• Dr. Kyle Horne

Department of Mechanical Engineering Spring, 2018

What is It?

20 40 60 80 100 120 140 Position x [px] 20 40 60 80 100 120 140 Position y [px]

20 25 30 35 40 45 Radial Position r [um] 22.5 23.0 23.5 24.0 24.5 Temperature T [C]

Fit Data

Explicit Implicit

15 10 5 5 10 15 Position x [px] 15 10 5 5 10 15 Position y [px] 225 230 235 240 Probe Laser r =18.80 [µm] 225 230 235 240 Position y [µm] r =32.30 [µm] 170 180 190 200 210 220 230 Position x [µm] 225 230 235 240 r =44.80 [µm]

Engineering Computation in Python; Dr. Kyle Horne 2 / 44

Where to Learn More?

Python Language Numpy (Array Support) Scipy (Numerical Methods) Sympy (Symbolic Math) Matplotlib (Plotting)

Engineering Computation in Python; Dr. Kyle Horne 3 / 44

Basic Usage

Procedural Execution, Comments and Text Formatting

g = 9.81 h = 100.0 # h = 0.5*g*t**2 -> t = sqrt(2*h/g) import math t = math.sqrt(2*h/g) print('t = %f [s]'%t)

Results

t = 4.515236

Engineering Computation in Python; Dr. Kyle Horne 4 / 44

Basic Usage

Lists and Arrays

import pylab as pl A = [1.0,2.0,3.0] A.append(4.0) B = pl.array(A) C = 2*B print(A) print(B) print(C)

Results

[1.0 , 2.0 , 3.0 , 4.0] [ 1. 2. 3. 4.] [ 2. 4. 6. 8.]

Engineering Computation in Python; Dr. Kyle Horne 5 / 44

Flow Control

If-Then-Else

a = 1 if a==1: print('a=1') elif a==2: print('a=2') else: print('a!=1')

Results

a=1

Engineering Computation in Python; Dr. Kyle Horne 6 / 44

Flow Control

For-Loops

s = 0 for k in range(4): x = 2.0*k-1.0 s = s+x print('k = %5d; x = %10f'%(k,x)) print('s = %f'%s)

Results

k = 0; x =

• 1.000000

k = 1; x = 1.000000 k = 2; x = 3.000000 k = 3; x = 5.000000 s = 8.000000

Engineering Computation in Python; Dr. Kyle Horne 7 / 44

Flow Control

Functions

def f(x): y = x**2-1.0 return y g = lambda x: x**2+1.0 x = 2.0 y = f(x) print('x = %f; y = %f; g(x) = %f'%(x,y,g(x)))

Results

x = 2.000000; y = 3.000000; g(x) = 5.000000

Engineering Computation in Python; Dr. Kyle Horne 8 / 44

Plotting

Line Plots and Boilerplate

import pylab as pl N = 100 x = pl.linspace(0,1,N) f = lambda x: x**2-1 y = f(x) fig = pl.figure( tight_layout=True, figsize=(5,4)) ax = fig.add_subplot(1,1,1) ax.plot(x,y,'b-',lw=2) ax.plot(x[::10],y[::10],'rs',mew=0) ax.set_xlabel('Abscissa $x$ [-]') ax.set_ylabel('Ordinate $y$ [-]') pl.savefig('examplePlot.pdf')

0.0 0.2 0.4 0.6 0.8 1.0 Abscissa x [-] 1.0 0.8 0.6 0.4 0.2 0.0 Ordinate y [-]

Engineering Computation in Python; Dr. Kyle Horne 9 / 44

Calculation of π (Series)

Theory

Leibniz’s series formula: π = 4

∞

∑

k=0

• (−1)k

2k + 1

• 10

20 30 40 50 Term Number k 10

1

100 Term Magnitude vk 10 20 30 40 50 Term Count N 2.6 3.0 3.4 3.8 Terms Sum vk 3.121595; N = 50

Engineering Computation in Python; Dr. Kyle Horne 10 / 44

Calculation of π (Series) I

Example Solution

1

import numpy

2

import pylab

3 4

N = 50

5 6

terms = numpy.empty(N)

7

for k in range(N):

8

terms[k] = 4.0*(-1)**k/(2.0*k+1.0)

9 10

sums = numpy.empty(N)

11

for k in range(N):

12

sums[k] = terms[:k+1].sum()

13 14

numbers = numpy.linspace(0,N-1,N)

15 16

fig = pylab.figure(tight_layout=True,figsize=(5,4))

17

ax = fig.add_subplot(1,1,1)

18

ax.plot(numbers,abs(terms),'-',lw=2)

19

ax.set_yscale('log')

20

ax.set_xlabel('Term Number $k$') Engineering Computation in Python; Dr. Kyle Horne 11 / 44

Calculation of π (Series) II

Example Solution

21

ax.set_ylabel('Term Magnitude $v_k$')

22

ax.set_title(r' ')

23

fig.savefig('piTermsPlot.pdf')

24 25

fig = pylab.figure(tight_layout=True,figsize=(5,4))

26

ax = fig.add_subplot(1,1,1)

27

ax.axhline(numpy.pi,color='k',linestyle='--',linewidth=2)

28

ax.plot(numbers,sums,'.',lw=2)

29

ax.set_yticks([2.6,3.0,numpy.pi,3.4,3.8])

30

ax.set_yticklabels(['2.6','3.0','$\pi$','3.4','3.8'])

31

ax.set_xlabel(r'Term Count $N$')

32

ax.set_ylabel(r'Terms Sum $\sum v_k$')

33

ax.text(5,3.8,r'$\pi \approx %f;\,N = %d$'%(terms.sum(),N))

34

fig.savefig('piSumsPlot.pdf') Engineering Computation in Python; Dr. Kyle Horne 12 / 44

Calculation of π (Monte Carlo)

Ac = πR2 As = (2R)2 Ac As = π 4 ri =

• x2

i + y2 i

π 4 = lim

N→∞

  count

i∈[1,N] (ri < 1)

N  

0.0 0.2 0.4 0.6 0.8 1.0 x-Coordinate 0.0 0.2 0.4 0.6 0.8 1.0 y-Coordinate

3.07200; N = 500

Engineering Computation in Python; Dr. Kyle Horne 13 / 44

Calculation of π (Monte Carlo) I

Example Solution

1

import pylab

2

import numpy

3

from numpy.random import rand

4 5

def calculatePI(N):

6

x = rand(N,2)

7

r = numpy.sqrt( (x**2).sum(1) )

8

c = (r<1.0).sum()

9 10

pi = 4.0*c/float(N)

11 12

th = numpy.linspace(0,numpy.pi/2,100)

13

xt = numpy.cos(th)

14

yt = numpy.sin(th)

15 16

xp = x[:,0]

17

yp = x[:,1]

18 19

fig = pylab.figure(tight_layout=True,figsize=(4,4))

20

ax = fig.add_subplot(1,1,1,aspect=1.0) Engineering Computation in Python; Dr. Kyle Horne 14 / 44

Calculation of π (Monte Carlo) II

Example Solution

21

ax.plot(xt,yt,'k--',lw=2)

22

ax.plot(xp[r<1],yp[r<1],'x',ms=7,mew=1)

23

ax.plot(xp[r>1],yp[r>1],'+',ms=7,mew=1)

24

ax.set_xlabel('$x$-Coordinate')

25

ax.set_ylabel('$y$-Coordinate')

26

ax.set_title( r'$\pi \approx %.5f; N=%d$'%(pi,N) )

27

fig.savefig('randomPI.pdf')

28 29

calculatePI(500) Engineering Computation in Python; Dr. Kyle Horne 15 / 44

Sieve of Eratosthenes

Theory

Sift the Twos and sift the Threes, The Sieve of Eratosthenes! When the multiples sublime, The numbers that remain are Prime.

Algorithm

1 P = {1, 2, . . . , N} 2 For i ∈ [2, N] 1 if P [i] is 0, then skip i 2 For j ∈ [2i, N] 1

P [j] = 0

101 102 103 104 Number n [#] 101 102 103 Number of Primes [#] 0.8 0.9 1.0 1.1 1.2 1.3 Ratio [1] (n) f(n) = n/log(n) f(n)/ (n)

Engineering Computation in Python; Dr. Kyle Horne 16 / 44

Sieve of Eratosthenes I

Example Solution

1

import numpy

2

import pylab

3 4

N = 10000+1

5

P = numpy.empty(N,dtype=int)

6 7

for i in range(N):

8

P[i] = i

9

for i in range(2,N):

10

if P[i]==0:

11

continue

12

for j in range(2*i,N,P[i]):

13

P[j] = 0

14 15

s = 2

16

px = (P[P!=0])[s:]

17

py = numpy.linspace(s,px.size,px.size)

18

fy = px/numpy.log(px)

19 20

fig = pylab.figure(tight_layout=True,figsize=(5,4)) Engineering Computation in Python; Dr. Kyle Horne 17 / 44

Sieve of Eratosthenes II

Example Solution

21

ax = fig.add_subplot(1,1,1)

22

ax.set_xscale('log')

23

ax.set_yscale('log')

24

lp, = ax.plot(px,py,'.')

25

lf, = ax.plot(px,fy,'-',lw=2)

26

ax.set_xlabel('Number $n$ [#]')

27

ax.set_ylabel('Number of Primes $\pi$ [#]')

28

ax2 = ax.twinx()

29

lr, = ax2.plot(px,fy/py,'C2_',mew=2)

30

ax2.set_ylabel(r'Ratio [1]')

31

lines = [lp,lf,lr]

32

labels = [r'$\pi(n)$',r'$f(n)=n/log(n)$',r'$f(n)/\pi(n)$']

33

ax2.legend(lines,labels,loc='upper center')

34

fig.savefig('seivePlot.pdf') Engineering Computation in Python; Dr. Kyle Horne 18 / 44

Orbital Simulation

Theory

1 au

Engineering Computation in Python; Dr. Kyle Horne 19 / 44

Orbital Simulation

Solver

Compute the acceleration on the Earth using Newton’s law of

• gravitation. Use this to update velocity at each time step. Use

velocity to update position at each time step.

• x (t)

= x0 +

t

• v (τ) dτ
• v (t)

=

• v0 +

t

• a (τ) dτ
• a (

x) = −G · M x3 x

• x0

= r0ˆ i

• v0

= ˆ j

• G · M

r0 r0 = 1 au The integrals can be converted into much simpler expressions by using Euler’s explicit method: f (x + ∆x) = f (x) +

x+∆x

x

f ′ (ξ) dξ ⋍ f (x) + ∆x · f ′ (x)

Engineering Computation in Python; Dr. Kyle Horne 20 / 44

Orbital Simulation

Algorithm

1

x ∈ R2 = (r0)ˆ i

2

v ∈ R2 =

• G·M

r0

• ˆ

j

3 t ∈ R = 0 4 While t < tend: 1

x = x + ∆t · v

2 r = |

x|

3

a = − G·M

r3

x

4

v = v + ∆t · a

5 t = t + ∆t

1.0 0.5 0.0 0.5 1.0 Position x [au] 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00 Position y [au] Orbit Sun Earth

Engineering Computation in Python; Dr. Kyle Horne 21 / 44

Orbital Simulation I

Example Solution

1

import pylab

2

import numpy

3 4

MS = 1.98892e30

5

au = 149.598e9

6

G = 6.67E-11

7 8

N = 365

9 10

U0 = numpy.sqrt(G*MS/au)

11

Y0 = 31556926.0

12

dt = Y0/365.0

13 14

xh = numpy.empty( (N,2) )

15 16

x = numpy.array([au,0.0])

17

v = numpy.array([0.0,U0])

18 19

xh[0,:] = x

20

Engineering Computation in Python; Dr. Kyle Horne 22 / 44

Orbital Simulation II

Example Solution

21

for k in range(1,N):

22

x = x+v*dt

23

r = numpy.sqrt( (x**2).sum() )

24

a = -G*MS*x/r**3

25

v = v+a*dt

26 27

xh[k] = x

28 29

fig = pylab.figure(figsize=(4,4),

30

tight_layout={'rect':(0,0,1,0.85)})

31

ax = fig.add_subplot(1,1,1,aspect=1.0)

32

ax.plot(xh[:,0]/au,xh[:,1]/au,'k--',lw=2,label='Orbit')

33

ax.plot([0],[0],'o',color='orange',

34

ms=30,mew=1.5,mec='yellow',label='Sun')

35

ax.plot([1],[0],'o',color='green',

36

ms=10,mew=1.5,mec='cyan',label='Earth')

37

ax.set_xlim(-1.1,1.1)

38

ax.set_ylim(-1.1,1.1)

39

ax.set_xlabel('Position $x$ [au]')

40

ax.set_ylabel('Position $y$ [au]') Engineering Computation in Python; Dr. Kyle Horne 23 / 44

Orbital Simulation III

Example Solution

41

ax.legend(frameon=False,ncol=3,numpoints=1,handletextpad=0.5,

42

bbox_to_anchor=(0,1.05),loc='lower left')

43

fig.savefig('orbitPlot.pdf') Engineering Computation in Python; Dr. Kyle Horne 24 / 44

Curve Fitting

Theory

Find Γ, α such that R (Γ, α) is minimized Vθ (r) = Γ r exp

• −αr2

R (Γ, α) =

• ∑

i

• (Vi − Vθ (ri))2

0.0 0.2 0.4 0.6 0.8 Position r [m] 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 Tangential Velocity V [m/s] Fit Experiment

from scipy.optimize import minimize

Engineering Computation in Python; Dr. Kyle Horne 25 / 44

Curve Fitting I

Example Solver

1

import numpy

2

import pylab

3

from scipy.optimize import minimize

4 5

re,Ve = numpy.loadtxt('vortex.txt',unpack=True)

6 7

def candidate(G,a):

8

g = lambda r: G/r*(1-numpy.exp(-a*r**2))

9

return g

10 11

def residual(x):

12

G,a = x

13

g = candidate(G,a)

14

R = numpy.sqrt( ( (Ve-g(re))**2 ).sum() )

15

return R

16 17

x0 = [1.0,200.0]

18

• pt = minimize(residual,x0)

19

G,a = opt.x

20

print(G,a) Engineering Computation in Python; Dr. Kyle Horne 26 / 44

Curve Fitting II

Example Solver

21 22

N = 100

23

rf = numpy.linspace(re.min(),re.max(),N)

24

gf = candidate(G,a)

25

Vf = gf(rf)

26 27

fig = pylab.figure(tight_layout=True,figsize=(5,4))

28

ax = fig.add_subplot(1,1,1)

29

ax.plot(rf,Vf,'--',lw=2,label='Fit')

30

ax.plot(re,Ve,'.',label='Experiment')

31

ax.set_xlabel(r'Position $r$ [m]')

32

ax.set_ylabel(r'Tangential Velocity $V_\theta$ [m/s]')

33

ax.legend()

34

fig.savefig('vortexPlot.pdf') Engineering Computation in Python; Dr. Kyle Horne 27 / 44

Artillery Fire

Theory

x (t) = v0t cos (θ) t (x) = x v0 cos (θ) y (x) = v0t (x) sin (θ) − 1 2gt (x)2 g (x) = x 50 2 g (xh) = y (xh)

100 200 300 400 500 Position x [m] 100 200 300 400 500 Position y [m]

from scipy.optimize import newton

Engineering Computation in Python; Dr. Kyle Horne 28 / 44

Artillery Fire I

Example Solver

1

import pylab

2

import numpy

3

from scipy.optimize import newton

4 5

def py(x,th,v0):

6

g = 9.81

7

t = x/(v0*numpy.cos(th))

8

y = t*(v0*numpy.sin(th)-0.5*g*t)

9

return y

10

def gy(x):

11

y = (x/50)**2

12

return y

13

def hit(th,v0):

14

f = lambda x: py(x,th,v0)-gy(x)

15

x = newton(f,1000.0)

16

return x

17 18

def plotTrajectory(th,v0):

19

N = 100

20

xh = hit(th,v0) Engineering Computation in Python; Dr. Kyle Horne 29 / 44

Artillery Fire II

Example Solver

21

yh = gy(xh)

22

xp = numpy.linspace(0,xh,N)

23

xg = numpy.linspace(-0.1*xh,1.2*xh,N)

24

p = py(xp,th,v0)

25

g = gy(xg)

26

fig = pylab.figure(tight_layout=True,figsize=(5,4))

27

ax = fig.add_subplot(1,1,1,aspect=1.0)

28

lg = ax.fill_between(xg,g,-30,color='g')

29

lp = ax.plot(xp,p,'--',lw=2)

30

lh = ax.plot([xh],[yh],'kx',ms=10,mew=4)

31

ax.set_xlabel('Position $x$ [m]')

32

ax.set_xlim(xg.min(),xg.max())

33

ax.set_ylim(-30,1.1*p.max())

34

ax.set_ylabel('Position $y$ [m]')

35

fig.savefig('artilleryPlot.pdf')

36 37

th = 75.0*numpy.pi/180.0

38

v0 = 100.0

39

plotTrajectory(th,v0) Engineering Computation in Python; Dr. Kyle Horne 30 / 44

Heated Wall

Theory

import sympy k d2T dx2 + Q = 0 k dT dx + Qx = C0 kT + Q 2 x2 = C0x + C1 ki Ti+1 − 2Ti + Ti−1 ∆x2 + Qi = 0

0.0 0.2 0.4 0.6 0.8 1.0 Position x [m] 26 28 30 32 34 36 38 Temperature T [C] Exact Numerical

ξ (i − 1) i (i + 1) N + 1 L x

Engineering Computation in Python; Dr. Kyle Horne 31 / 44

Heated Wall I

Example Solver

1

import pylab

2

import numpy

3

import scipy.sparse as sparse

4

import scipy.sparse.linalg as linalg

5

import sympy

6 7

Q = 5000.0 # W

8

k = 50.0 # W/m.K

9

L = 1.0 # m

10 11

T0 = 25.0 # K

12

TL = 25.0 # K

13 14

def plotSolution(xe,Te,xn,Tn):

15

fig = pylab.figure(tight_layout=True,figsize=(5,4))

16

ax = fig.add_subplot(1,1,1)

17

eline = ax.plot(xe,Te,'-',lw=2,label='Exact')

18

nline = ax.plot(xn,Tn,'v',lw=2,label='Numerical')

19

ax.legend(loc='best')

20

ax.set_xlabel('Position $x$ [m]') Engineering Computation in Python; Dr. Kyle Horne 32 / 44

Heated Wall II

Example Solver

21

ax.set_ylabel('Temperature $T$ [C]')

22

fig.savefig('wallPlot.pdf')

23 24

def finiteDifferences(N):

25

dx = L/(N+1)

26 27

DP = -2.0*numpy.ones(N)

28

DL = numpy.ones(N-1)

29

DU = numpy.ones(N-1)

30

A = sparse.diags([DL,DP,DU],(-1,0,1))

31

A = A.tocsc()

32 33

S = (-Q*dx**2/k)*numpy.ones(N)

34

S[0] = S[0]-T0

35

S[N-1] = S[N-1]-TL

36 37

T = linalg.spsolve(A,S)

38 39

x = numpy.linspace(0,L,N+2)

40

T = numpy.insert(T,0,T0) Engineering Computation in Python; Dr. Kyle Horne 33 / 44

Heated Wall III

Example Solver

41

T = numpy.append(T,TL)

42 43

return x,T

44 45

def differentialEquation(N):

46

x = sympy.symbols('x')

47

T = sympy.symbols('T',cls=sympy.Function)

48

deq = sympy.Eq( k*T(x).diff(x,x)+Q , 0 )

49

sol = sympy.dsolve(deq,T(x))

50

C1,C2 = sympy.symbols('C1 C2')

51 52

system = []

53

system.append( sol.subs(x,0).subs(T(0),T0) )

54

system.append( sol.subs(x,L).subs(T(L),TL) )

55

coeffs = sympy.solve(system,[C1,C2])

56

Ts = sympy.solve( sol.subs(coeffs) , T(x) )[0]

57

Tf = sympy.lambdify(x,Ts,'numpy')

58 59

x = numpy.linspace(0,L,N+2)

60

T = Tf(x) Engineering Computation in Python; Dr. Kyle Horne 34 / 44

Heated Wall IV

Example Solver

61 62

return x,T

63 64

xn,Tn = finiteDifferences(6)

65

xe,Te = differentialEquation(100)

66

plotSolution(xe,Te,xn,Tn) Engineering Computation in Python; Dr. Kyle Horne 35 / 44

Duct Flow

Theory

from scipy.integrate import ode ∂w ∂t = −1 ρ ∂P ∂z + ν

• ∂2w

∂x2 + ∂2w ∂y2

• 1

2 3 4 5 6 7 Time t [ms] 2 4 6 8 10 12 14 Max Velocity w [m/s] 4 2 2 4 Position x [cm] 4 2 2 4 Position y [cm] 0.0 1.6 3.2 4.8 6.4 8.0 9.6 11.2 12.8 14.4 Velocity w [m/s]

Engineering Computation in Python; Dr. Kyle Horne 36 / 44

Duct Flow I

Example Solver

1

import pylab

2

import numpy

3

from scipy.integrate import ode

4 5

#=============#

6

#= Constants =#

7

#=============#

8 9

dPdz = -10.0e3

10

mu = 0.1

11

rho = 1.0

12

nu = mu/rho

13

L = 10.0e-2

14

D = L/2.0

15 16

Dc = 4*(L**2-numpy.pi*D**2/4)/(4*L)

17

Vc = numpy.sqrt(2*abs(dPdz)*Dc/rho)

18

gRe = rho*Vc*Dc/mu

19 20

N = 75 Engineering Computation in Python; Dr. Kyle Horne 37 / 44

Duct Flow II

Example Solver

21

x = numpy.linspace(-L/2,L/2,N)

22

y = numpy.linspace(-L/2,L/2,N)

23

dx = L/(N-1)

24

dy = L/(N-1)

25 26

X,Y = numpy.meshgrid(x,y)

27

R = numpy.sqrt(X**2+Y**2)

28 29

mask = numpy.logical_and(Y<L/15,Y>-L/15)

30

mask = numpy.logical_and(X>0.0,mask)

31

mask = numpy.logical_or(R<D/2.0,mask)

32 33

#==============#

34

#= Derivative =#

35

#==============#

36 37

def dwdt(t,w):

38

wl = w.reshape( (N,N) )

39

• = numpy.zeros( (N,N) )

40

Engineering Computation in Python; Dr. Kyle Horne 38 / 44

Duct Flow III

Example Solver

41

• [1:-1,1:-1] = -dPdz/rho+nu*(

42

(wl[2:,1:-1]-2*wl[1:-1,1:-1]+wl[:-2,1:-1])/dx**2 +

43

(wl[1:-1,2:]-2*wl[1:-1,1:-1]+wl[1:-1,:-2])/dy**2

44

)

45 46

• [mask] = 0.0

47

return o.reshape(-1)

48 49

#==================#

50

#= Initialization =#

51

#==================#

52 53

w0 = numpy.zeros( (N,N) ).reshape(-1)

54

t = numpy.linspace(0.0,0.007,100)

55

wm = numpy.zeros(t.size)

56

wm[0] = w0.max()

57 58

I = ode( dwdt )

59

I.set_initial_value(w0,0.0)

60

Engineering Computation in Python; Dr. Kyle Horne 39 / 44

Duct Flow IV

Example Solver

61

#============#

62

#= Solution =#

63

#============#

64 65

for k in range(1,t.size):

66

I.integrate(t[k])

67

wm[k] = I.y.max()

68

print('%10d %20.5e %20.5e'%(k,t[k],wm[k]))

69 70

#============#

71

#= Plotting =#

72

#============#

73 74

w = I.y.reshape( (N,N) )

75 76

fig = pylab.figure(tight_layout=True,figsize=(5,4))

77

ax = fig.add_subplot(1,1,1)

78

ax.plot(1e3*t,wm,'-',lw=2)

79

ax.set_ylim(0,1.1*wm.max())

80

ax.set_xlabel('Time $t$ [ms]') Engineering Computation in Python; Dr. Kyle Horne 40 / 44

Duct Flow V

Example Solver

81

ax.set_ylabel('Max Velocity $w$ [m/s]')

82

fig.savefig('velocityMaxPlot.pdf')

83 84

fig = pylab.figure(tight_layout=True,figsize=(5,4))

85

ax = fig.add_subplot(1,1,1,aspect=1.0)

86

cf = ax.contourf(1e2*x,1e2*y,w,20,cmap=pylab.get_cmap('viridis'))

87

cl = ax.contour(1e2*x,1e2*y,w,20,cmap=pylab.get_cmap('viridis'))

88

ax.set_xlabel('Position $x$ [cm]')

89

ax.set_ylabel('Position $y$ [cm]')

90

fig.colorbar(cf,ax=ax,label='Velocity $w$ [m/s]')

91

fig.savefig('velocityProfilePlot.pdf')

92 93

V = w.max()

94

Re = rho*V*Dc/mu

95 96

print('')

97

print(gRe)

98

print( Re) Engineering Computation in Python; Dr. Kyle Horne 41 / 44

Group Work

Theory

from glog import glob z = −1 2g0t2 + ˙ z0t + z0 2h = g0t2

f

0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 Fall Time Squared, t2

f [s2]

2 4 6 8 10 Doubled Fall Height, 2 h [m] Data 1 Data 2 Data 3 Data 4 g0 = 9.27 m/s2

Data

# Height [m], Time [s] 0.057 0.168 0.125 0.239 0.501 0.404 0.613 0.449 2.396 0.781 2.404 0.752 2.539 0.798 2.710 0.780 3.068 0.867 3.481 0.862

Engineering Computation in Python; Dr. Kyle Horne 42 / 44

Group Work I

Example Solver

1

import pylab as pl

2

import numpy as np

3

from glob import glob

4 5

def getData():

6

fns = sorted(glob('*.txt'))

7

D = []

8

for fn in fns:

9

D.append( np.loadtxt(fn) )

10

return np.array(D)

11 12

def doPlot(D,C):

13

fig = pl.figure(figsize=(5,4),tight_layout=True)

14

ax = fig.add_subplot(1,1,1)

15 16

ms = 'os^v'

17

for k in range(D.shape[0]):

18

h = D[k,:,0]

19

t = D[k,:,1]

20

ax.plot(t**2,h,ms[k],label='Data %d'%(k+1)) Engineering Computation in Python; Dr. Kyle Horne 43 / 44

Group Work II

Example Solver

21 22

x = np.linspace((D[:,:,1]**2).min(),(D[:,:,1]**2).max())

23

y = np.polyval(C,x)/2.0

24

ax.plot(x,y,'k--',label='$g_0=%.2f$ $m/s^2$'%C[0],zorder=1)

25 26

ax.set_xlabel('Fall Time Squared, $t_f^2$ [$s^2$] ')

27

ax.set_ylabel('Doubled Fall Height, $2 \cdot h$ [$m$]')

28

ax.legend()

29 30

fig.savefig('Group.pdf')

31 32

def doFit(D):

33

x = D[:,:,1].flatten()**2

34

y = 2.0*D[:,:,0].flatten()

35

C = np.polyfit(x,y,1)

36

return C

37 38

D = getData()

39

C = doFit(D)

40

doPlot(D,C) Engineering Computation in Python; Dr. Kyle Horne 44 / 44