function [u_x, f] = constant_rho_left_right(Nx, Ny, f, u_x, rho_left, rho_right) % This function should be used in conjuction with stream_2 or similar. % It sents the boundary conditions on the left and right boundaries. % Both boundaries are constant-pressure (constant-density) boundaries. % rho_left is the density on the left boundary. % rho_right is the density on the right boundary. % To avoid compressibility (mass balance) problems, these should not be % very far apart from each other. I have used rho_left = 1.01 and % rho_right = 0.99 and there was no problem with those values. i = 1; % left boundary for j = 1:Ny u_x(i,j) = 1 - (f(9,i,j) + f(2,i,j) + f(4,i,j) + 2*(f(3,i,j) + f(6,i,j) + f(7,i,j)))/rho_left; f(1,i,j) = f(3,i,j) + (2/3)*rho_left*u_x(i,j); f(5,i,j) = f(7,i,j) - (1/2)*(f(2,i,j) - f(4,i,j)) + (1/6)*rho_left*u_x(i,j); f(8,i,j) = f(6,i,j) + (1/2)*(f(2,i,j) - f(4,i,j)) + (1/6)*rho_left*u_x(i,j); end i = Nx; % right boundary for j = 1:Ny u_x(i,j) = -1 + (f(9,i,j) + f(2,i,j) + f(4,i,j) + 2*(f(1,i,j) + f(5,i,j) + f(8,i,j))) / rho_right; f(3,i,j) = f(1,i,j) - (2/3)*rho_right*u_x(i,j); f(6,i,j) = f(8,i,j) - (1/6)*rho_right*u_x(i,j) - (1/2)*(f(2,i,j) - f(4,i,j)); f(7,i,j) = f(5,i,j) - (1/6)*rho_right*u_x(i,j) + (1/2)*(f(2,i,j) - f(4,i,j)); end end