I need to develop a code to solve Stokes Equations in 3D in cubic geometries (structured grid, uniform mesh spacing).
My code needs to take a pressure gradient in one direction as a BC (pinlet=p1, poutlet=p0) and no-slip BC's in the other 4 faces of the cubic geometry.
I started to develop this code in Python, using a staggered grid FV approach; pressure at the center of the FV cells, ux velocities at the faces normal to the x direction and so on.
My code gives strange results, which has made me wonder if I am considering the correct approach to solve such system.
I am currently assembling a system of equations that considers [p, ux, uy] as unknowns and solving it in "one go". That is, I build a matrix and solve the linear system, which gives me the solution. The degrees of freedom are coupled through Stokes equations and continuity equation, but the coupling is "solved" by assembling the coefficients to the matrix.
However, whenever I research numerical codes and books on Navier-Stokes equations, they usually give a couple of strategies to treat the pressure-velocity coupling, typically "guess and correct" approaches, where you guess a pressure field, and correct it iteratively to satisfy continuity (e.g. SIMPLE, SIMPLEC, PISO). These approaches are usually formulated in the Navier-Stokes context, but Stokes equations have the same type of pressure-velocity coupling; they just don't have the nonlinear inertia term. So I assume these strategies would be applicable to Stokes equations as well.
So my question to the community is: Is there something "wrong" with my approach to solve Stokes equations in one go, implicitly, by solving one linear system? Or do I need to consider iterative methods to treat the pressure-velocity coupling?
Thank you! Rafael March.