前几日检索到了由韦伯州立大学物理系的Dan Schroeder教授编写的利用格子玻尔兹曼方法进行流体模拟的程序（原网站：Fluid Dynamics Simulation (weber.edu)）。网站给出的Python程序仅对一个线状截面进行了模拟，我心血来潮在其基础上使用Python编写了GUI程序，通过调节滑块实现机翼截面角度的变化，实现了对其周围流体进行模拟。
但正如Dan Schroeder教授所说，此模拟用于定性和半定量教育演示，而非用于严肃的工程用途，明显的局限性在于，它仅模拟了二维流体而非三维流体。（This simulation is intended for qualitative and semi-quantitative educational demonstrations—not for serious engineering use. One obvious limitation is that it simulates a fluid in only two dimensions rather than three.）因此不能用于科学研究。
如果你也有过跟我相同的想法，欢迎与我交流。
我编写的程序如下所示，所用到的库有Tkinter、Matplotlib以及Numpy。另外，此代码使用 NACA 2412 翼型坐标作为示例生成一个简单的翼型图：
import tkinter as tk
from tkinter import ttk
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import numpy as np
import time
class AirfoilApp:
def __init__(self, root):
self.root = root
self.root.title("Fluid Simulation At Different Airfoil Angles")
self.create_widgets()
def create_widgets(self):
self.frame = ttk.Frame(self.root)
self.frame.pack(fill=tk.BOTH, expand=1)
# Create matplotlib figure and axis
self.fig, (self.ax1, self.ax2) = plt.subplots(1, 2, figsize=(10, 5))
#self.fig, (self.ax2) = plt.subplots(1, figsize=(10, 5))
self.canvas = FigureCanvasTkAgg(self.fig, master=self.frame)
self.canvas.get_tk_widget().pack(fill=tk.BOTH, expand=1)
# Slider for rotation
self.slider = ttk.Scale(self.root, from_=0, to=-360, orient=tk.HORIZONTAL, command=self.update_plot)
self.slider.pack(fill=tk.X)
self.initialize_simulation()
self.draw_airfoil()
self.ani = animation.FuncAnimation(self.fig, self.update_simulation, interval=50, blit=True,cache_frame_data=False)
def initialize_simulation(self):
self.height = 200
self.width = 200
self.viscosity = 0.1
self.omega = 1 / (3 * self.viscosity + 0.5)
self.u0 = 0.2
self.four9ths = 4.0 / 9.0
self.one9th = 1.0 / 9.0
self.one36th = 1.0 / 36.0
performanceData = False
self.n0 = self.four9ths * (np.ones((self.height, self.width)) - 1.5 * self.u0 ** 2)
self.nN = self.one9th * (np.ones((self.height, self.width)) - 1.5 * self.u0 ** 2)
self.nS = self.one9th * (np.ones((self.height, self.width)) - 1.5 * self.u0 ** 2)
self.nE = self.one9th * (np.ones((self.height, self.width)) + 3 * self.u0 + 4.5 * self.u0 ** 2 - 1.5 * self.u0 ** 2)
self.nW = self.one9th * (np.ones((self.height, self.width)) - 3 * self.u0 + 4.5 * self.u0 ** 2 - 1.5 * self.u0 ** 2)
self.nNE = self.one36th * (np.ones((self.height, self.width)) + 3 * self.u0 + 4.5 * self.u0 ** 2 - 1.5 * self.u0 ** 2)
self.nSE = self.one36th * (np.ones((self.height, self.width)) + 3 * self.u0 + 4.5 * self.u0 ** 2 - 1.5 * self.u0 ** 2)
self.nNW = self.one36th * (np.ones((self.height, self.width)) - 3 * self.u0 + 4.5 * self.u0 ** 2 - 1.5 * self.u0 ** 2)
self.nSW = self.one36th * (np.ones((self.height, self.width)) - 3 * self.u0 + 4.5 * self.u0 ** 2 - 1.5 * self.u0 ** 2)
self.rho = self.n0 + self.nN + self.nS + self.nE + self.nW + self.nNE + self.nSE + self.nNW + self.nSW
self.ux = (self.nE + self.nNE + self.nSE - self.nW - self.nNW - self.nSW) / self.rho
self.uy = (self.nN + self.nNE + self.nNW - self.nS - self.nSE - self.nSW) / self.rho
self.create_airfoil_barrier(0)
def create_airfoil_barrier(self, angle):
self.barrier = np.zeros((self.height, self.width), bool)
x, y = self.generate_airfoil()
x_rot, y_rot = self.rotate_airfoil(x, y, angle)
for i in range(len(x_rot)):
xi = int((x_rot[i] + 1) * self.width / 2)
yi = int((y_rot[i] + 1) * self.height / 2)
if 0 <= xi < self.width and 0 <= yi < self.height:
self.barrier[yi//2+50, xi//2+50] = True
self.update_barrier_neighbors()
def update_barrier_neighbors(self):
self.barrierN = np.roll(self.barrier, 1, axis=0)
self.barrierS = np.roll(self.barrier, -1, axis=0)
self.barrierE = np.roll(self.barrier, 1, axis=1)
self.barrierW = np.roll(self.barrier, -1, axis=1)
self.barrierNE = np.roll(self.barrierN, 1, axis=1)
self.barrierNW = np.roll(self.barrierN, -1, axis=1)
self.barrierSE = np.roll(self.barrierS, 1, axis=1)
self.barrierSW = np.roll(self.barrierS, -1, axis=1)
def generate_airfoil(self):
t = 0.12 # Maximum thickness as a fraction of the chord
c = 1.0 # Chord length
x = np.linspace(0, c, 100)
yt = 5 * t * c * (0.2969 * np.sqrt(x/c) - 0.1260 * (x/c) - 0.3516 * (x/c)**2 + 0.2843 * (x/c)**3 - 0.1015 * (x/c)**4)
return np.concatenate([x, x[::-1]]), np.concatenate([yt, -yt[::-1]])
def rotate_airfoil(self, x, y, angle):
angle_rad = np.radians(angle)
x_rot = x * np.cos(angle_rad) - y * np.sin(angle_rad)
y_rot = x * np.sin(angle_rad) + y * np.cos(angle_rad)
return x_rot, y_rot
def draw_airfoil(self, angle=0):
self.ax1.clear()
x, y = self.generate_airfoil()
x_rot, y_rot = self.rotate_airfoil(x, y, angle)
self.ax1.plot(x_rot, y_rot)
self.ax1.set_aspect('equal')
self.ax1.grid(True)
self.ax1.set_xlim(-1.1,1.1)
self.ax1.set_ylim(-1.1,1.1)
self.create_airfoil_barrier(angle)
self.canvas.draw()
def update_plot(self, angle):
self.draw_airfoil(float(angle))
def stream(self):
self.nN = np.roll(self.nN, 1, axis=0)
self.nNE = np.roll(self.nNE, 1, axis=0)
self.nNW = np.roll(self.nNW, 1, axis=0)
self.nS = np.roll(self.nS, -1, axis=0)
self.nSE = np.roll(self.nSE, -1, axis=0)
self.nSW = np.roll(self.nSW, -1, axis=0)
self.nE = np.roll(self.nE, 1, axis=1)
self.nNE = np.roll(self.nNE, 1, axis=1)
self.nSE = np.roll(self.nSE, 1, axis=1)
self.nW = np.roll(self.nW, -1, axis=1)
self.nNW = np.roll(self.nNW, -1, axis=1)
self.nSW = np.roll(self.nSW, -1, axis=1)
self.nN[self.barrierN] = self.nS[self.barrier]
self.nS[self.barrierS] = self.nN[self.barrier]
self.nE[self.barrierE] = self.nW[self.barrier]
self.nW[self.barrierW] = self.nE[self.barrier]
self.nNE[self.barrierNE] = self.nSW[self.barrier]
self.nNW[self.barrierNW] = self.nSE[self.barrier]
self.nSE[self.barrierSE] = self.nNW[self.barrier]
self.nSW[self.barrierSW] = self.nNE[self.barrier]
def collide(self):
self.rho = self.n0 + self.nN + self.nS + self.nE + self.nW + self.nNE + self.nSE + self.nNW + self.nSW
self.ux = (self.nE + self.nNE + self.nSE - self.nW - self.nNW - self.nSW) / self.rho
self.uy = (self.nN + self.nNE + self.nNW - self.nS - self.nSE - self.nSW) / self.rho
ux3 = self.ux ** 2
uy3 = self.uy ** 2
u2 = ux3 + uy3
omu215 = 1 - 1.5 * u2
uxuy = self.ux * self.uy
self.n0 = self.n0 * (1 - self.omega) + self.omega * self.four9ths * omu215
self.nN = self.nN * (1 - self.omega) + self.omega * self.one9th * (omu215 + 3 * self.uy + 4.5 * uy3)
self.nS = self.nS * (1 - self.omega) + self.omega * self.one9th * (omu215 - 3 * self.uy + 4.5 * uy3)
self.nE = self.nE * (1 - self.omega) + self.omega * self.one9th * (omu215 + 3 * self.ux + 4.5 * ux3)
self.nW = self.nW * (1 - self.omega) + self.omega * self.one9th * (omu215 - 3 * self.ux + 4.5 * ux3)
self.nNE = self.nNE * (1 - self.omega) + self.omega * self.one36th * (omu215 + 3 * (self.ux + self.uy) + 4.5 * (u2 + 2 * uxuy))
self.nNW = self.nNW * (1 - self.omega) + self.omega * self.one36th * (omu215 + 3 * (-self.ux + self.uy) + 4.5 * (u2 - 2 * uxuy))
self.nSE = self.nSE * (1 - self.omega) + self.omega * self.one36th * (omu215 + 3 * (self.ux - self.uy) + 4.5 * (u2 - 2 * uxuy))
self.nSW = self.nSW * (1 - self.omega) + self.omega * self.one36th * (omu215 + 3 * (-self.ux - self.uy) + 4.5 * (u2 + 2 * uxuy))
def curl(self,ux,uy):
return np.roll(uy, -1, axis=1) - np.roll(uy, 1, axis=1) - np.roll(ux, -1, axis=0) + np.roll(ux, 1, axis=0)
def update_simulation(self, frame):
self.stream()
self.collide()
#im = self.ax2.imshow(self.rho, cmap='jet', animated=True)
im = self.ax2.imshow(self.curl(self.ux, self.uy), origin='lower', norm=plt.Normalize(-.1, .1),
cmap=plt.get_cmap('jet'), interpolation='none')
bImageArray = np.zeros((self.height, self.width, 4), np.uint8) # an RGBA image
bImageArray[self.barrier, 3] = 255 # set alpha=255 only at barrier sites
barrierImage = self.ax2.imshow(bImageArray, origin='lower', interpolation='none')
return im,barrierImage
root = tk.Tk()
app = AirfoilApp(root)
root.mainloop()
3052
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
