File:Navier Stokes Laminar.svg

from __future__ import divisionfrom numpy import arange, meshgrid, sqrt, zeros, sumimport matplotlib.pyplot as pltfrom mpl_toolkits.mplot3d import Axes3Dfrom matplotlib.ticker import ScalarFormatterfrom matplotlib import rcParams rcParams['font.family'] = 'serif'rcParams['font.size'] = 16 # the layout of the duct laminarx_max = 5 # duct lengthy_max = 1 # duct width# draw the frames, including the angles and labelsax = Axes3D(plt.figure(figsize=(10, 8)), azim=20, elev=20)ax.set_xlabel(r"$x$", fontsize=20)ax.set_ylabel(r"$y$", fontsize=20)ax.zaxis.set_rotate_label(False)ax.set_zlabel(r"$v_x$", fontsize=20, rotation='horizontal')formatter = ScalarFormatter(useMathText=True)formatter = ScalarFormatter()formatter.set_scientific(True)formatter.set_powerlimits((-2,2))ax.w_zaxis.set_major_formatter(formatter)ax.set_xlim([0, x_max])ax.set_ylim([0, y_max])# initial speed of the airini_v = 3e-3mu = 1e-5rho = 1.3# the acceptable difference when terminationaccept_diff = 1e-5# time intervaltime_delta = 1.0# coordinate intervaldelta = 1e-2;X = arange(0, x_max + delta, delta)Y = arange(0, y_max + delta, delta)# number of coordinate pointsx_size = len(X) - 1y_size = len(Y) - 1Vx = zeros((len(X), len(Y)))Vy = zeros((len(X), len(Y)))new_Vx = zeros((len(X), len(Y)))new_Vy = zeros((len(X), len(Y)))# initial conditionsVx[1: x_size - 1, 2:y_size - 1] = ini_v# start evolution and computationres = 1 + accept_diffrounds = 0alpha = mu/(rho * delta**2)while (res>accept_diff and rounds<100): """ The iterations here are based on the velocity change rate, which is given by  \frac{\partial v}{\partial t} = \alpha\nabla^2 v - v \cdot \nabla v  with \alpha = \mu/\rho. """ new_Vx[2:-2, 2:-2] = Vx[2:-2, 2:-2] + time_delta*(alpha*(Vx[3:-1, 2:-2] + Vx[2:-2, 3:-1] - 4*Vx[2:-2, 2:-2] + Vx[2:-2, 1:-3] + Vx[1:-3, 2:-2]) - 0.5/delta * (Vx[2:-2, 2:-2] * (Vx[3:-1, 2:-2] - Vx[1:-3, 2:-2]) + Vy[2:-2, 2:-2]*(Vx[2:-2, 3:-1] - Vx[2:-2, 1:-3]))) new_Vy[2:-2, 2:-2] = Vy[2:-2, 2:-2] + time_delta*(alpha*(Vy[3:-1, 2:-2] + Vy[2:-2, 3:-1] - 4*Vy[2:-2, 2:-2] + Vy[2:-2, 1:-3] + Vy[1:-3, 2:-2]) - 0.5/delta * (Vy[2:-2, 2:-2] * (Vy[2:-2, 3:-1] - Vy[2:-2, 3:-1]) + Vx[2:-2, 2:-2]*(Vy[3:-1, 2:-2] - Vy[1:-3, 2:-2])))  rounds = rounds + 1  # copy the new values Vx[2:-2, 2:-2] = new_Vx[2:-2, 2:-2] Vy[2:-2, 2:-2] = new_Vy[2:-2, 2:-2] # set free boundary conditions: dv_x/dx = dv_y/dx = 0. Vx[-1, 1:-1] = Vx[-3, 1:-1] Vx[-2, 1:-1] = Vx[-3, 1:-1] Vy[-1, 1:-1] = Vy[-3, 1:-1] Vy[-2, 1:-1] = Vy[-3, 1:-1] # there exists a still object in the plane Vx[x_size//3:x_size//1.5, y_size//2.0] = 0 Vy[x_size//3:x_size//1.5, y_size//2.0] = 0 # calculate the residual of Vx res = (Vx[3:-1, 2:-2] + Vx[2:-2, 3:-1] - Vx[1:-3, 2:-2] - Vx[2:-2, 1:-3])**2 res = sum(res)/(4 * delta**2 * x_size * y_size)# prepare the plot dataZ = sqrt(Vx**2)# refine the region boundaryZ[0, 1:-2] = Z[1, 1:-2]Z[-2, 1:-2] = Z[-3, 1:-2]Z[-1, 1:-2] = Z[-3, 1:-2]Y, X = meshgrid(Y, X);ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap="summer", lw=0.1, edgecolors="k")plt.savefig("Navier_Stokes_Laminar.svg")