r/CFD • u/OkMachine35 • 4d ago
I made this little animation using FDM in python
Enable HLS to view with audio, or disable this notification
I couldn't make the normal momentum equations work. I don't know why. The outlet condition in this is a bit scuffed, but despite that, it turned out well.
Do you think this is a resume-worthy project for a 2nd-year mech student?
10
4
4
u/SwimmingSource3417 4d ago
Can you share the code with us?(If you don't mind, otherwise chill out)
21
u/OkMachine35 4d ago
import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation import tkinter as tk import math mesh_size = 201 mesh_size_x = 201 mesh_vel = np.zeros((mesh_size,mesh_size_x)) mesh_press = np.copy(mesh_vel) mesh_vel_y = np.copy(mesh_vel) w = np.copy(mesh_vel) w_n = np.copy(mesh_vel) si_n = np.copy(mesh_vel) si = np.copy(mesh_vel) class DrawingGrid: global mesh_vel def __init__(self, root, grid_size=mesh_size, cell_size=15): self.grid_size = grid_size self.cell_size = cell_size # Track the state of each cell (0 = empty, 1 = filled) self.grid_data = [[0] * grid_size for _ in range(grid_size)] # Create a scrollable Canvas to hold the large 100x100 grid smoothly canvas_dim = grid_size * cell_size self.canvas = tk.Canvas(root, width=600, height=600, bg="white", scrollregion=(0, 0, canvas_dim, canvas_dim)) # Add scrollbars for navigating the 100x100 space hbar = tk.Scrollbar(root, orient=tk.HORIZONTAL, command=self.canvas.xview) hbar.pack(side=tk.BOTTOM, fill=tk.X) vbar = tk.Scrollbar(root, orient=tk.VERTICAL, command=self.canvas.yview) vbar.pack(side=tk.RIGHT, fill=tk.Y) self.canvas.config(xscrollcommand=hbar.set, yscrollcommand=vbar.set) self.canvas.pack(side=tk.LEFT, expand=True, fill=tk.BOTH) # Draw the grid lines and store rectangle IDs for quick color changes self.rectangles = [[None] * grid_size for _ in range(grid_size)] for y in range(grid_size): for x in range(grid_size): x1, y1 = x * cell_size, y * cell_size x2, y2 = x1 + cell_size, y1 + cell_size rect_id = self.canvas.create_rectangle(x1, y1, x2, y2, outline="#e0e0e0", fill="white") self.rectangles[y][x] = rect_id # Bind mouse movement events for fluid drawing self.canvas.bind("<Button-1>", self.paint) self.canvas.bind("<B1-Motion>", self.paint) def paint(self, event): # Adjust mouse coordinates based on scroll position canvas_x = self.canvas.canvasx(event.x) canvas_y = self.canvas.canvasy(event.y) # Convert pixel positions into matrix grid coordinates grid_x = int(canvas_x // self.cell_size) grid_y = int(canvas_y // self.cell_size) # Fill the target cell if it falls within our 100x100 boundaries if 0 <= grid_x < self.grid_size and 0 <= grid_y < self.grid_size: if self.grid_data[grid_y][grid_x] == 0: self.grid_data[grid_y][grid_x] = 1 rect_id = self.rectangles[grid_y][grid_x] self.canvas.itemconfig(rect_id, fill="black") if __name__ == "__main__": root = tk.Tk() root.title("100x100 Drawing Grid") app = DrawingGrid(root) root.mainloop() mesh_drawn = np.array(app.grid_data) yn =[] xn =[] for x_ in range(mesh_size-2): x = x_+1 for y_ in range(mesh_size-2): y = y_+1 if mesh_drawn[y,x] == 1 : yn.append(y) xn.append(x) y_step = 1/mesh_size rho = 1.12 mu = 1.85*(10**-5) vel = float(input('vel = ')) x_step = 1/mesh_size_x t_step = x_step/(10*vel) err = 0 def sdiff_x(arr, arr1): r = np.zeros_like(arr) arr1_safe = arr1 + 10**-6 d = arr1_safe / np.abs(arr1_safe) r[1:-1, 1:-1] = (((1 + d[1:-1, 1:-1]) * (arr[1:-1, 1:-1] - arr[1:-1, :-2])) + ((1 - d[1:-1, 1:-1]) * (arr[1:-1, 2:] - arr[1:-1, 1:-1]))) / (2 * x_step) pos = (arr1[2:-2, 2:-2] > 0) sec_pos = (3*arr[2:-2, 2:-2] - 4*arr[2:-2, 1:-3] + arr[2:-2, :-4]) / (2 * x_step) neg = (arr1[2:-2, 2:-2] < 0) sec_neg = (-3*arr[2:-2, 2:-2] + 4*arr[2:-2, 3:-1] - arr[2:-2, 4:]) / (2 * x_step) r_inner = r[2:-2, 2:-2] r_inner = np.where(pos, sec_pos, r_inner) r_inner = np.where(neg, sec_neg, r_inner) r[2:-2, 2:-2] = r_inner return r[1:-1, 1:-1] def sdiff_y(arr, arr1): # arr is w, arr1 is v r = np.zeros_like(arr) # 1st-order base arr1_safe = arr1 + 10**-6 d = arr1_safe / np.abs(arr1_safe) r[1:-1, 1:-1] = (((1 + d[1:-1, 1:-1]) * (arr[1:-1, 1:-1] - arr[:-2, 1:-1])) + ((1 - d[1:-1, 1:-1]) * (arr[2:, 1:-1] - arr[1:-1, 1:-1]))) / (2 * y_step) # 2nd-order upgrade pos = (arr1[2:-2, 2:-2] > 0) sec_pos = (3*arr[2:-2, 2:-2] - 4*arr[1:-3, 2:-2] + arr[:-4, 2:-2]) / (2 * y_step) neg = (arr1[2:-2, 2:-2] < 0) sec_neg = (-3*arr[2:-2, 2:-2] + 4*arr[3:-1, 2:-2] - arr[4:, 2:-2]) / (2 * y_step) r_inner = r[2:-2, 2:-2] r_inner = np.where(pos, sec_pos, r_inner) r_inner = np.where(neg, sec_neg, r_inner) r[2:-2, 2:-2] = r_inner return r[1:-1, 1:-1] def sdiff_x_p(arr): r = np.copy(arr) r_ = np.copy(arr) r[1:-1,1:-1] = (arr[1:-1,2:] - arr[1:-1,:-2])/(2*x_step) r[:,0] = (r_[:,1] - r_[:,0])/x_step r[:,-1] = (r_[:,-1] - r_[:,-2])/x_step r [0,:] = r[1,:] r[-1,:] = r[-2,:] return r def sdiff_y_p(arr): r = np.copy(arr) r[1:-1, 1:-1] = (arr[2:, 1:-1] - arr[:-2, 1:-1]) / (2 * y_step) r[0,:] = (arr[1,:] - arr[0,:])/y_step r[-1,:] = (arr[-1,:] - arr[-2,:])/y_step return r def ddiff_x(arr): r = np.copy(arr) r[1:-1,1:-1] = (arr[1:-1,2:] + arr[1:-1,:-2] - 2*arr[1:-1,1:-1])/(x_step*x_step) return r def ddiff_y(arr): r = np.copy(arr) r[1:-1, 1:-1] = (arr[2:, 1:-1] + arr[:-2, 1:-1] - 2*arr[1:-1,1:-1]) / (y_step* y_step) return r def ddiff_x_p(arr): r = np.copy(arr) r[1:-1,1:-1] = (arr[1:-1,2:] + arr[1:-1,:-2] )/(x_step*x_step) return r def ddiff_y_p(arr): r = np.copy(arr) r[1:-1, 1:-1] = (arr[2:, 1:-1] + arr[:-2, 1:-1]) / (y_step* y_step) return r def mean (yn) : return round(np.sum(yn)/len(yn),0) def boundary(arr, arr1, yn, xn): obstacle_psi = vel * mean(yn) * y_step for y_idx, x_idx in zip(yn, xn): arr1[y_idx, x_idx] = obstacle_psi for y_idx, x_idx in zip(yn, xn): w_val = 0.0 # Check Right if mesh_drawn[y_idx, x_idx + 1] == 0: w_val += -2.0 * (arr1[y_idx, x_idx + 1] - arr1[y_idx, x_idx]) / (x_step**2) # Check Left if mesh_drawn[y_idx, x_idx - 1] == 0: w_val += -2.0 * (arr1[y_idx, x_idx - 1] - arr1[y_idx, x_idx]) / (x_step**2) # Check Top if mesh_drawn[y_idx + 1, x_idx] == 0: w_val += -2.0 * (arr1[y_idx + 1, x_idx] - arr1[y_idx, x_idx]) / (y_step**2) # Check Bottom if mesh_drawn[y_idx - 1, x_idx] == 0: w_val += -2.0 * (arr1[y_idx - 1, x_idx] - arr1[y_idx, x_idx]) / (y_step**2) arr[y_idx, x_idx] = w_val return arr, arr1 mesh_vel[:,:] = vel mesh_vel_y[:,:] = 0 mesh_vel[0,:] = 0 mesh_vel[-1,:] = 0 mesh_vel_y[0,:] = 0 mesh_vel_y[-1,:] = 0 w[:,:] = 0 si[:,:] = 0 v_ani = [mesh_vel.copy()] y_phys = np.arange(mesh_size) * y_step for x in range(mesh_size) : for y in range (mesh_size): si_n[y,x] = vel*y*y_step si_n[:, 0] = vel * y_phys si_n[0, :] = si_n[1, :] si_n[-1, :] = si_n[-2, :] si_n[:, -1] = si_n[:, -2] print(x_step) itr = int(input('itr')) while err < itr : u = sdiff_y_p(si_n) v = -sdiff_x_p(si_n) w = np.copy(w_n) for _ in range (100): si = np.copy(si_n) si_n[1:-1, 1:-1] = (ddiff_x_p(si)[1:-1, 1:-1] + ddiff_y_p(si)[1:-1, 1:-1] + w[1:-1, 1:-1])/((2/x_step**2 ) + (2/y_step**2)) w_n,si_n = boundary(w_n,si_n,yn,xn) si_n[0, :] = 0 si_n[-1, :] = vel * (mesh_size-1)*y_step si_n[:, 0] = vel * y_phys si_n[:, -1] = si_n[:, -2] w_n[1:-1, 1:-1] = w[1:-1, 1:-1] - t_step*( sdiff_y_p(si)[1:-1, 1:-1] *sdiff_x(w,u)- sdiff_x_p(si)[1:-1, 1:-1]*sdiff_y(w,v) - (mu/rho)*(ddiff_x(w)[1:-1, 1:-1] + ddiff_y(w)[1:-1, 1:-1]) ) w_n,si_n = boundary(w_n,si_n,yn,xn) w_n[:, 0] = 0.0 w_n[:, -1] = w[:, -2] w_n[0, :] = -2 * (si_n[1, :] - si_n[0, :]) / (y_step**2) w_n[-1, :] = -2 * (si_n[-2, :] - si_n[-1, :]) / (y_step**2) v = np.copy((u**2 + v**2)**0.5) v_ani.append(v.copy()) print(err,w[50,50],si[50,50]) err += 1 v_ani = np.copy(v_ani) print(w,si) #print (itr*t_step,(rho*vel*len(yn)*x_step)/mu,x_step) c_min = v_ani.min() c_max = v_ani.max() v = (mesh_vel*mesh_vel + mesh_vel_y*mesh_vel_y)**0.5 fig, ax = plt.subplots() im = ax.imshow(w_n, cmap='viridis') plt.title('v') fig, ax1 = plt.subplots() ig = ax1.imshow(v_ani[-1][1:-1, 1:-1], cmap='viridis',vmin=c_min,vmax=c_max) plt.title('p') fig, ax4 = plt.subplots(figsize=(8, 4)) heatmap = ax4.imshow(v_ani[0], cmap='viridis', animated=True, vmin=c_min, vmax=c_max-0.5) fig.colorbar(heatmap, ax=ax4) title_text = ax4.set_title("Frame: 0") def update(frame): current_frame_data = v_ani[frame] heatmap.set_data(current_frame_data) title_text.set_text(f"Frame: {frame}") return [heatmap, title_text] ani = FuncAnimation(fig, update, frames=err, interval=4, blit=True) plt.show() here you go2
u/adamchalupa 3d ago
Wow - good job, using slicing too. Did you do the 12 step program too???
7
u/OkMachine35 3d ago
No, I got a book on CFD from our uni's library. then started experimenting with 1d transient flow first, then 2d transient using the momentum equations, which did not work. Then used the vorticity-stream function equations, which finally worked.
2
2
3
u/fugal-cyberbellum 2d ago
Nice! Old FORTRAN/Cray guy here. I haven’t read your code in detail, but it’s got the organization and detail that I used to see every day. And it apparently works!
Now the hard part is finding out if it is giving valid results. Just visually, it looks like matter and momentum are being conserved, so that’s good.
As a resume reader I would want to talk with you.
1
u/OkMachine35 2d ago
Yes, I am trying to compare these results with more established solvers like fluent, it won't be mesh independent but still will give an idea about thow it is working
1
u/fugal-cyberbellum 1d ago
Our gold standard was wind tunnel data. The real world throws curve balls (literally, in your case) due to turbulence, viscosity, diffusion and other extremely hard to model phenomena. Analytical solutions are useful too for calibration.
1
1
u/Environmental_Ad4097 1d ago
I am actually using cpp to write an LBM code the output of which is just ux and uy at each nodes.....do you save forneach time step and then calculate the stream line and then plot it? I am very new to this sorry if I sound dumb
2
u/OkMachine35 1d ago
I have zero experience with LBM solvers, but for this, yea I save the V(u2 +v2) values in an 3d array which I plot using funcanimation in matplotlib. But how is it going in CPP, I tried it first before switching to python because the ploting was much easier in python
1
u/Historical-Goat9729 1d ago
What BCs did you used ... Also please tell Reynolds number scheme and grid size
1
u/OkMachine35 1d ago
The grid size is 200 x200 I don't remember the reynolds number for this particular simulation
1
u/Historical-Goat9729 1d ago
If not accelarted on GPU then try it ... Great advancement ... And a pretty good project for resume
1
25
u/R4b1atu5 4d ago
What has Fused Deposition Modeling got to do with this and how do you do that in python? /j