• API
• FAQ
• Tools
• Archive
daily pastebin goal
12%
SHARE
TWEET

# Untitled

a guest Mar 19th, 2019 48 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import numpy as np
2. import matplotlib.pyplot as plt
3.
4. # plate size
5. w = h = 10.
6. # intervals in x-, y- directions
7. dx = dy = 0.1
8. # Thermal diffusivity of steel
9. D = 4.
10.
11. Tcool, Thot = 300, 1700
12.
13. nx, ny = int(w/dx), int(h/dy)
14.
15. dx2, dy2 = dx*dx, dy*dy
16. dt = dx2 * dy2 / (2 * D * (dx2 + dy2))
17.
18. u0 = Tcool * np.ones((nx, ny))
19. u = np.empty((nx, ny))
20.
21. for i in range(1, nx-1):
22.     for j in range(1, ny-1):
23.         uxx = (u0[i+1,j] - 2*u0[i,j] + u0[i-1,j]) / dx2
24.         uyy = (u0[i,j+1] - 2*u0[i,j] + u0[i,j-1]) / dy2
25.         u[i,j] = u0[i,j] + dt * D * (uxx + uyy)
26.
27. # Initial conditions - ring of inner radius r, width dr centred at (cx,cy) (mm)
28. r, cx, cy = 2, 5, 5
29. r2 = r**2
30. def applyBC():
31.     for i in range(nx):
32.         for j in range(ny):
33.             p2 = (i*dx-cx)**2 + (j*dy-cy)**2
34.             if p2 < r2:
35.                 u0[i,j] = Thot
36.
37. def do_timestep(u0, u):
38.     applyBC()
39.     # Propagate with forward-difference in time, central-difference in space
40.     u[1:-1, 1:-1] = u0[1:-1, 1:-1] + D * dt * (
41.           (u0[2:, 1:-1] - 2*u0[1:-1, 1:-1] + u0[:-2, 1:-1])/dx2
42.           + (u0[1:-1, 2:] - 2*u0[1:-1, 1:-1] + u0[1:-1, :-2])/dy2 )
43.
44.     u0 = u.copy()
45.     return u0, u
46.
47. # Number of timesteps
48. nsteps = 2000
49. # Output 4 figures at these timesteps
50. mfig = [0, 10, 50, 1999]
51. fignum = 0
52. fig = plt.figure()
53. for m in range(nsteps):
54.     u0, u = do_timestep(u0, u)
55.     if m in mfig:
56.         fignum += 1
57.         print(m, fignum)
58.         ax = fig.add_subplot(220 + fignum)
59.         im = ax.imshow(u.copy(), cmap=plt.get_cmap('hot'), vmin=Tcool,vmax=Thot)
60.         ax.set_axis_off()
61.         ax.set_title('{:.1f} ms'.format(m*dt*1000))