Python coding help

Run steady-state simulations in python. Repeat as transient simulations.

How long would it take to reach steady-state in this problem? (For transient use initial temperatures of 20C and 500C respectively).

A thin-film 500-W electrical heater of negligible thickness and dimensions 100 mm 100 mm is sandwiched between two slabs made of copper alloy (k = 120 W/mK) and stainless steel (15 W/mK). The two plates, copper and stainless steel with thickness 50 mm and 50 mm, are surrounded by air with heat transfer coefficients of 75 W/m2 K and 25 W/m2 K, respectively. The plates are subjected to a constant temperature of 50C and 30C at the top and the bottom surface, respectively, as shown in Fig. P583. The temperature of the surrounding air is maintained at 20C. Determine (a) the finite difference formulation at each node and (b) the location of maximum temperature.

import numpy as np import matplotlib.pyplot as plt print("2D heat equation solver") x_length = 0.5 y_length = 0.3 z_length = 5 max_iters = 100 #maximum number of iterations on solver kappa = 23 # W/m-K delta_x = 0.1 # 10 cm xnorm = x_length/delta_x + 1 #node on each end need to add 1 ynorm = y_length/delta_x + 1 xnorm = round(xnorm) ynorm = round(ynorm) print("nodes in x: ", xnorm) print("nodes in y: ", ynorm) # Initialize solution: the grid of u(k, i, j) T = np.empty((ynorm, xnorm)) # Initial condition everywhere inside the grid T_initial = 273 # Boundary conditions # need to set up qdot for top #T_top = 100 qtop = 8000/(z_length*x_length) # W of heat added to top, 8000W / (5 m * 0.5m) T_left = 273 #temperature in K #T_bottom = 0.0 #bottom is insulated BC so we won't use this T_right = 273 #Set the initial condition T.fill(T_initial) #Set the boundary conditions #top is heat in so: #qtop*deltax + k*deltay/2*(Tleft-T0)/deltax + k*deltay/2*(Tright-T0)/deltax + #k*deltax*(Tbot-T0)/deltax == rho*Cp * (T0(t+1) - T0(t))/deltat #because dx = dy this simplifies to: #T0(t+1)= alpha*deltat* [qtop*deltax +k*(Tleft/2+Tright/2+Tbot-2*T0)] +T0 #This equation will be used in the for loop, for now just set to initial condition #T[(x_length-1):, :] = T_top #this would be for fixed T BC #Fixed temperatures for left and right: T[:, :1] = T_left T[:, (xnorm-1):] = T_right #Insulated BC #T[:1, 1:] = T_bottom maxTchange = 1 tempT = np.empty((ynorm, xnorm)) iter_to_solve = 0 while(maxTchange>=0.001): #for k in range(max_iters): iter_to_solve = iter_to_solve + 1 tempT[:, :] = T.copy() for i in range(0, ynorm, 1): for j in range(1, xnorm-1, 1): if i == ynorm - 1: # top row of matrix has qdot 8000... use equation #defined above, indexing starts at 0 so need ynorm - 1 index for top T[i][j] = (T[i - 1][j] + T[i][j + 1] / 2 + T[i][j - 1] / 2) / 2 + qtop / kappa * delta_x / 2 elif i == 0: # insulated BC on bottom, so no T(i-1,j) T[i][j] = (T[i + 1][j] + T[i][j + 1] / 2 + T[i][j - 1] / 2) / 2 else: T[i][j] = (T[i + 1][j] + T[i - 1][j] + T[i][j + 1] + T[i][j - 1]) / 4 deltaT = np.subtract(T, tempT[:, :]) deltaT = abs(deltaT) maxTchange = np.max(deltaT) print("Max delta T: ", maxTchange) print("Maximum Temperature : ", np.max(T)) print(f"Solution took {iter_to_solve:.0f} iterations to complete") # convert to C from K T_convert = np.empty((ynorm, xnorm)) T_convert.fill(273) T = np.subtract(T, T_convert) plt.clf() plt.title("Temperature (K)") plt.xlabel("x") plt.ylabel("y") # This is to plot u_k (u at time-step k) xrange = np.mgrid[slice(0, x_length + delta_x, delta_x)] yrange = np.mgrid[slice(0, y_length + delta_x, delta_x)] plt.pcolormesh(xrange, yrange, T, cmap=plt.cm.jet, vmin=round(np.min(T)), vmax=round(np.max(T))) plt.colorbar() print("Done!")