Open - Source coded 3D Navier–Stokes equations in C++

     Here are 3D Navier–Stokes equations configured for lid-driven cavity flow. The syntax is C++. The results from this code are shown in Fig. 1. From the pressure-Poisson equation, I removed mixed derivative terms to improve solution stability. 😆

p[i][j][k] = ((p[i+1][j][k] + p[i-1][j][k]) * h * h + (p[i][j+1][k] + p[i][j-1][k]) * h * h + (p[i][j][k+1] + p[i][j][k-1]) * h * h) / (2 * (h * h + h * h + h * h)) - h * h * h * h * h * h / (2 * (h * h + h * h + h * h)) * (rho * (1 / dt * ((u(i+1, j, k) - u(i-1, j, k)) / (2 * h) + (v(i, j+1, k) - v(i, j-1, k)) / (2 * h) + (w(i, j, k+1) - w(i, j, k-1)) / (2 * h))));

p[0][j][k] = p[1][j][k];

p[num_i - 1][j][k] = p[num_i - 2][j][k];

p[i][0][k] = p[i][1][k];

p[i][num_j - 1][k] = p[i][num_j - 2][k];

p[i][j][0] = p[i][j][1];

p[i][j][num_k - 1] = 0.0;

u[i][j][k] = u[i][j][k] - u[i][j][k] * dt / (2 * h) * (u[i+1][j][k] - u[i-1][j][k]) - v[i][j][k] * dt / (2 * h) * (u[i][j+1][k] - u[i][j-1][k]) - w[i][j][k] * dt / (2 * h) * (u[i][j][k+1] - u[i][j][k-1]) - dt / (2 * rho * h) * (p[i+1][j][k] - p[i-1][j][k]) + nu * (dt / (h * h) * (u[i+1][j][k] - 2 * u[i][j][k] + u[i-1][j][k]) + dt / (h * h) * (u[i][j+1] [k] - 2 * u[i][j][k] + u[i][j-1][k]) + dt / (h * h) * (u[i][j][k+1] - 2 * u[i][j][k] + u[i][j][k-1]));

u[0][j][k] = 0.0;

u[num_i - 1][j][k] = 0.0;

u[i][0][k] = 0.0;

u[i][num_j - 1][k] = 0.0;

u[i][j][0] = 0.0;

u[i][j][num_k - 1] = -Uinf;

v[i][j][k] = v[i][j][k] - u[i][j][k] * dt / (2 * h) * (v[i+1][j][k] - v[i-1][j][k]) - v[i][j][k] * dt / (2 * h) * (v[i][j+1][k] - v[i][j-1][k]) - w[i][j][k] * dt / (2 * h) * (v[i][j][k+1] - v[i][j][k-1]) - dt / (2 * rho * h) * (p[i][j+1][k] - p[i][j-1][k]) + nu * (dt / (h * h) * (v[i+1][j][k] - 2 * v[i][j][k] + v[i-1][j][k]) + dt / (h * h) * (v[i][j+1][k] - 2 * v[i][j][k] + v[i][j-1][k]) + dt / (h * h) * (v[i][j][k+1] - 2 * v[i][j][k] + v[i][j][k-1]));

v[0][j][k] = 0.0;

v[num_i - 1][j][k] = 0.0;

v[i][0][k] = 0.0;

v[i][num_j - 1][k] = 0.0;

v[i][j][0] = 0.0;

v[i][j][num_k - 1] = 0.0;

w[i][j][k] = w[i][j][k] - u[i][j][k] * dt / (2 * h) * (w[i+1][j][k] - w[i-1][j][k]) - v[i][j][k] * dt / (2 * h) * (w[i][j+1][k] - w[i][j-1][k]) - w[i][j][k] * dt / (2 * h) * (w[i][j][k+1] - w[i][j][k-1]) - dt / (2 * rho * h) * (p[i][j][k+1] - p[i][j][k-1]) + nu * (dt / (h * h) * (w[i+1][j][k] - 2 * w[i][j][k] + w[i-1][j][k]) + dt / (h * h) * (w[i][j+1][k] - 2 * w[i][j][k] + w[i][j-1][k]) + dt / (h * h) * (w[i][j][k+1] - 2 * w[i][j][k] + w[i][j][k-1]));

w[0][j][k] = 0.0;

w[num_i - 1][j][k] = 0.0;

w[i][0][k] = 0.0;

w[i][num_j - 1][k] = 0.0;

w[i][j][0] = 0.0;

w[i][j][num_k - 1] = 0.0;

Fig. 1, Velocity and pressure iso-surfaces

     Of course constants need to be defined, such as grid spacing in space and time, density, kinematic viscosity. These equations have been validated, as you might have read and here already! Happy coding!

     If you want to hire me as your PhD student in the research projects related to turbo-machinery, aerodynamics, renewable energy, please reach out. Thank you very much for reading.