Hi Hussein,
I understand that you are facing difficulty with constructing a matrix system, specifically (Ax = b), for simulating incompressible 2D Navier-Stokes equations (NSEs) using a dual-splitting scheme. Following is a breakdown of process into steps outlined and then constructing the matrix system (Ax = b) for each step. The dual-splitting scheme helps in decoupling the velocity and pressure fields, making the system easier to solve.
Step 1: Convective-Diffusion Equation
First, solve for an intermediate velocity field (U*) hat does not satisfy the incompressibility constraint. This step involves solving convective-diffusion equation:
For simplicity, assuming a uniform grid and using finite difference methods for spatial discretization and an explicit scheme for time discretization. Let (h) be the spatial step size and (Δt) the time step. The grid points are ((𝑥𝑖,𝑦𝑗)) with (i,j) indexing the grid in the (x) and (y) directions, respectively.
Matrix Formulation- For a 2D grid, vectorize the velocity fields (u) and (v) into long vectors and construct the corresponding matrices for the Laplacian and gradient operations.The system for (U*) becomes:
where,
encapsulates both the convection and diffusion effects and 
includes the forcing term (f) and contributions from the boundary conditions. Step 2: Poisson-Pressure Correction
Next, solve for the pressure correction by ensuring the divergence-free condition on the velocity field. This leads to a Poisson equation for the pressure:
Use central differences for the Laplacian of the pressure.
Matrix Formulation- Vectorize the pressure field (P) into a long vector. The matrix system for the pressure correction step becomes:
where,
is the discrete Laplacian operator, and 
is derived from the divergence of (U*) and the time step (Δt).MATLAB Implementation Example:
Here is a rough outline to implement these steps in MATLAB. This example does not include specific boundary conditions or initialization of matrices, which one will need to adapt based on problem setup.
% Parameters
nx = 50; % Number of grid points in x
ny = 50; % Number of grid points in y
dx = 1/(nx-1); % Spatial step size
dt = 0.01; % Time step
mu = 1; % Viscosity
% Initialize fields
U_star = zeros(nx*ny, 1); % Intermediate velocity field
P = zeros(nx*ny, 1); % Pressure field
f = @(t,x,y) sin(x)*cos(y); % Forcing function
% Construct A1 and A2 matrices for discretized operators
A1 = constructA1(nx, ny, dx, dt, mu); % You need to define this function
A2 = constructA2(nx, ny, dx); % You need to define this function
% Step 1: Solve for U_star
b1 = constructB1(f, nx, ny, dx, dt); % You need to define this function
U_star = A1 \ b1;
% Step 2: Solve for pressure correction
b2 = constructB2(U_star, nx, ny, dx, dt); % You need to define this function
P_new = A2 \ b2;
% Update velocity field based on pressure correction
% This step requires back-substitution to get the corrected velocity field
Note: constructA1, constructA2, constructB1, and constructB2 are placeholder functions which needs to be defined based on discretization of operators and specifics of problem setup. This example is simplified and assumes familiarity with constructing finite difference matrices in MATLAB. For non-linear terms and complex boundary conditions, the construction of (A) matrices and (b) vectors will be more involved.
Hope it helps!
Best Regards,
Simar
