SHARE
TWEET

Untitled

a guest Mar 19th, 2019 50 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))
  62. fig.subplots_adjust(right=0.85)
  63. cbar_ax = fig.add_axes([0.9, 0.15, 0.03, 0.7])
  64. cbar_ax.set_xlabel('$T$ / K', labelpad=20)
  65. fig.colorbar(im, cax=cbar_ax)
  66. plt.show()
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top