I have tried to calculate the curl of the velocity in cylindrical coordinates, following the equations from https://en.wikipedia.org/wiki/Del_in_cylindrical_and_spherical_coordinates
I have assumed finite differences for my derivatives. This is the curl calculation
def curl(grid,u):
# Define the spacing between grid points in the θ direction
delta_theta =2* np.pi/2 / 30 # Replace with the actual spacing between grid points in the θ direction
delta_z = 2*np.pi/1.7 / 42 # Replace with the actual spacing between grid points in the z direction
delta_r = np.diff(grid)
# Create arrays to store the derivatives
du_theta_dr = np.zeros_like(u[..., 0])
du_r_dtheta = np.zeros_like(u[..., 1])
du_z_dr = np.zeros_like(u[..., 0])
du_r_dz = np.zeros_like(u[..., 2])
du_z_dtheta = np.zeros_like(u[..., 1])
du_theta_dz = np.zeros_like(u[..., 2])
# Calculate the derivatives using finite differences
for i in range(1, u.shape[0]-1):
for j in range(1, u.shape[1]-1):
for k in range(1, u.shape[2]-1):
dr_plus = grid[i+1] - grid[i]
dr_minus = grid[i] - grid[i-1]
du_theta_dr[i, j, k] = (u[i+1, j, k, 1] - u[i-1, j, k, 1]) / (dr_plus + dr_minus)
du_r_dtheta[i, j, k] = (u[i, j+1, k, 0] - u[i, j-1, k, 0]) / (2 * delta_theta)
du_z_dr[i, j, k] = (u[i+1, j, k, 2] - u[i-1, j, k, 2]) / (dr_plus + dr_minus)
du_r_dz[i, j, k] = (u[i, j, k+1, 0] - u[i, j, k-1, 0]) / (2 * delta_z)
du_z_dtheta[i, j, k] = (u[i, j+1, k, 2] - u[i, j-1, k, 2]) / (2 * delta_theta)
du_theta_dz[i, j, k] = (u[i, j, k+1, 1] - u[i, j, k-1, 1]) * grid[i]/ (2 * delta_z)
curl_u_r = 1/grid[:, np.newaxis, np.newaxis] * du_z_dtheta - du_theta_dz
curl_u_theta = (du_r_dz - du_z_dr)
curl_u_z = (1/grid[:, np.newaxis, np.newaxis]) * (du_theta_dr - du_r_dtheta)
curl_u = np.stack((curl_u_r, curl_u_theta, curl_u_z), axis=-1)
return curl_u
but I am not sure whether it is correct.
