Assignment # 2
ME-831
Computational Fluid Dynamics
Code for explicit time marching technique
Google colab was used to develop the code for assignment. It is highlighted that code
can be just copied and pasted online on google colab
(https://colab.research.google.com/) without use of python to check its authenticity and
functioning.
It is also highlighted that online solvers for solution of one-dimensional heat equation were
available on google colab, however, they were not used and code was written watching
tutorials and reference material on youtube and google colab.
Loops for different values of alpha is already made part of the code to avoid repetitive
running of code.
import numpy as np
import matplotlib.pyplot as plt
h = 0.1
k = 0.1
x = np.arange(0, 1 + h, h)
t = np.arange(0, 3, k)
n = len(x)
m = len(t)
boundaryConditions = [300, 370]
initialConditions = [300]
# Different values of alpha to be tested
alpha_values = [0.01, 0.03, 0.05, 0.08, 0.1,]
# Create subplots
fig, axs = plt.subplots(len(alpha_values), 1, figsize=(8, 6 *
len(alpha_values)))
# Loop over each alpha value
for idx, alpha in enumerate(alpha_values):
T = np.ones((n, m)) * 300
factor = alpha * (k / h**2)
for j in range(0, m-1):
for i in range(1, n-1):
T[i, j + 1] = T[i, j] + factor * (T[i + 1, j] - 2 * T[i, j] +
T[i - 1, j])
# Apply boundary conditions
T[0, j + 1] = boundaryConditions[0]
T[-1, j + 1] = boundaryConditions[1]
# Apply initial conditions (optional)
T[:, 0] = initialConditions[0]
T = np.round(T, 3)
# Plot temperature profiles at selected times for the current alpha
axs[idx].plot(x, T)
#axs[idx].plot(x, T[:, 0], label=f'Time = 0s')
#axs[idx].plot(x, T[:, int(1 / k)], label=f'Time = 1s')
#axs[idx].plot(x, T[:, int(2 / k)], label=f'Time = 2s')
#axs[idx].plot(x, T[:, -1], label=f'Time = 3s')
axs[idx].set_xlabel('Distance [m]')
axs[idx].set_ylabel('Temperature [K]')
axs[idx].set_title(f'Temperature Profiles for alpha = {alpha}')
#axs[idx].legend()
# Adjust layout for better appearance
plt.tight_layout()
plt.show()
Code for Implicit time marching technique
import numpy as np
import matplotlib.pyplot as plt
# Parameters
L = 1.0 # Length of the rod
T = 3.0 # Total simulation time
Nx = 11 # Number of spatial grid points
Nt = 30 # Number of time steps
initial_temperature = 300.0
left_boundary_temperature = 300.0
right_boundary_temperature = 370.0
# Spatial and Temporal Discretization
dx = L / (Nx - 1)
dt = T / Nt
# Create spatial grid
x = np.linspace(0, L, Nx)
# Initialize temperature field
T = np.ones((Nx, Nt+1)) * initial_temperature
# Set initial and boundary conditions
T[:, 0] = initial_temperature
T[0, :] = left_boundary_temperature
T[-1, :] = right_boundary_temperature
# Coefficients for the Crank-Nicolson method
alpha = 0.04 * dt / (dx**2 * 2)
A = np.diag(1 + 2 * alpha * np.ones(Nx - 2)) + np.diag(-alpha * np.ones(Nx
- 3), -1) + np.diag(-alpha * np.ones(Nx - 3), 1)
B = np.diag(1 - 2 * alpha * np.ones(Nx - 2)) + np.diag(alpha * np.ones(Nx
- 3), -1) + np.diag(alpha * np.ones(Nx - 3), 1)
# Time-stepping loop
for j in range(0, Nt):
# Right-hand side vector
b = B @ T[1:-1, j]
b[0] += alpha * (T[0, j] + T[0, j+1])
b[-1] += alpha * (T[-1, j] + T[-1, j+1])
# Solve the system of equations using NumPy's linalg.solve
solution = np.linalg.solve(A, b)
# Update temperature for the next time step
T[1:-1, j+1] = solution
# Plot temperature profiles at selected times
plt.plot(x, T[:, 0], label='Time = 0s')
plt.plot(x, T[:, int(Nt/6)], label='Time = 0.5s')
plt.plot(x, T[:, int(Nt/3)], label='Time = 1s')
plt.plot(x, T[:, int(Nt/2)], label='Time = 1.5s')
plt.plot(x, T[:, int(2*Nt/3)], label='Time = 2s')
plt.plot(x, T[:, int(5*Nt/6)], label='Time = 2.5s')
plt.plot(x, T[:, -1], label='Time = 3s')
plt.xlabel('Distance [m]')
plt.ylabel('Temperature [K]')
plt.title('Temperature Profile using Crank-Nicolson Method')
plt.legend()
plt.show()
