-
Notifications
You must be signed in to change notification settings - Fork 0
/
rankinevortexsourcesink.py
144 lines (120 loc) · 4.86 KB
/
rankinevortexsourcesink.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
import numpy
import math
from matplotlib import pyplot
import sys
def main():
N = 200 # Number of points in each direction
x_start, x_end = -4.0, 4.0 # x-direction boundaries
y_start, y_end = -2.0, 2.0 # y-direction boundaries
x = numpy.linspace(x_start, x_end, N) # 1D-array for x
y = numpy.linspace(y_start, y_end, N) # 1D-array for y
X, Y = numpy.meshgrid(x, y) # generates a mesh grid
u_inf = 1.0 # freestream speed
# compute the freestream velocity field
u_freestream = u_inf * numpy.ones((N, N), dtype=float)
v_freestream = numpy.zeros((N, N), dtype=float)
# compute the stream-function
psi_freestream = u_inf * Y
strength_source = 5.0 # strength of the source
x_source, y_source = -1.0, 0.0 # location of the source
# compute the velocity field
u_source, v_source = get_velocity(strength_source, x_source, y_source, X, Y)
# compute the stream-function
psi_source = get_stream_function(strength_source, x_source, y_source, X, Y)
global u,v
u = u_freestream + u_source
v = v_freestream + v_source
psi = psi_freestream + psi_source
# plot the streamlines
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyplot.figure(figsize=(width, height))
pyplot.grid(True)
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
pyplot.streamplot(X, Y, u, v, density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyplot.scatter(x_source, y_source, color='#CD2305', s=80, marker='o')
# calculate the stagnation point
x_stagnation = x_source - strength_source / (2 * numpy.pi * u_inf)
y_stagnation = y_source
# display the stagnation point
pyplot.scatter(x_stagnation, y_stagnation, color='g', s=80, marker='o')
# display the dividing streamline
pyplot.contour(X, Y, psi,
levels=[-strength_source / 2, strength_source / 2],
colors='#CD2305', linewidths=2, linestyles='solid');
pyplot.show()
strength_sink = -5.0 # strength of the sink
x_sink, y_sink = 1.0, 0.0 # location of the sink
# compute the velocity field on the mesh grid
u_sink, v_sink = get_velocity(strength_sink, x_sink, y_sink, X, Y)
# compute the stream-function on the grid mesh
psi_sink = get_stream_function(strength_sink, x_sink, y_sink, X, Y)
u = u_freestream + u_source + u_sink
v = v_freestream + v_source + v_sink
psi = psi_freestream + psi_source + psi_sink
# plot the streamlines
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyplot.figure(figsize=(width, height))
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
pyplot.streamplot(X, Y, u, v,
density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyplot.scatter([x_source, x_sink], [y_source, y_sink],
color='#CD2305', s=80, marker='o')
pyplot.contour(X, Y, psi,
levels=[0.], colors='#CD2305', linewidths=2, linestyles='solid');
pyplot.show()
def get_velocity(strength, xs, ys, X, Y):
"""
Returns the velocity field generated by a source/sink.
Parameters
----------
strength: float
Strength of the source/sink.
xs: float
x-coordinate of the source (or sink).
ys: float
y-coordinate of the source (or sink).
X: 2D Numpy array of floats
x-coordinate of the mesh points.
Y: 2D Numpy array of floats
y-coordinate of the mesh points.
Returns
-------
u: 2D Numpy array of floats
x-component of the velocity vector field.
v: 2D Numpy array of floats
y-component of the velocity vector field.
"""
u = strength / (2 * numpy.pi) * (X - xs) / ((X - xs)**2 + (Y - ys)**2)
v = strength / (2 * numpy.pi) * (Y - ys) / ((X - xs)**2 + (Y - ys)**2)
return u, v
def get_stream_function(strength, xs, ys, X, Y):
"""
Returns the stream-function generated by a source/sink.
Parameters
----------
strength: float
Strength of the source/sink.
xs: float
x-coordinate of the source (or sink).
ys: float
y-coordinate of the source (or sink).
X: 2D Numpy array of floats
x-coordinate of the mesh points.
Y: 2D Numpy array of floats
y-coordinate of the mesh points.
Returns
-------
psi: 2D Numpy array of floats
The stream-function.
"""
psi = strength / (2 * numpy.pi) * numpy.arctan2((Y - ys), (X - xs))
return psi
main()