How to solve for pressure using the velocity field?

22 views (last 30 days)
Orlando Rojas
Orlando Rojas on 25 Jan 2015
Edited: Youssef Khmou on 26 Jan 2015
Hi,
I'm trying to solve for the pressure field on a 120x120 grid. I have the velocity field, and I have to use the Pressure Poisson Equation to solve for the pressure. The known boundary condition is P=Patm along the right edge of the grid. Using this I can get all the boundary conditions. I've created my coefficient matrix and I've tried so many different ways to solve the problems but I keep getting and error saying "Out of memory". I've tried solving by the easiest way which is Ax=b, x=A\b. Also by using reduced row echelon form. I'm now trying to solve using lower triangular matrix and the other an upper triangular matrix, but when I try using the function lu on MATLAB, again I get the same "Out of memory" message. Any comments or ideas are appreciated.
Thank you,
Orlando R.
  5 Comments
Star Strider
Star Strider on 25 Jan 2015
OK. Thanks. I didn’t look at the .mat files yet.

Sign in to comment.

Answers (1)

Youssef  Khmou
Youssef Khmou on 26 Jan 2015
This problem is for non compressible fluid, first you have to explain the equation you want use, for this two dimensional case, the Poission Pressure equation is D(P)=d((1/Re).D(u)-u.d(u)), D is the Laplacian, d is divergence and Re denotes the Reynold number. What is not clear is that you need only one matrix for velocity u as field of x and y components, anyway i suggest that you use Jacobi method; for (i,j) component of P at iteration m+1 u have :
P(i,j,m+1) = ( P(i-1,j,m) + P(i+1,j,m) + P(i,j-1,m) + P(i,j+1,m) +C(i,j) )/4
C(i,j) is result of right hand side of equation, parallel to this approach try to use gradient function of velocity U to obtain C.
  2 Comments
Youssef  Khmou
Youssef Khmou on 26 Jan 2015
For discrete solution, you initialize the matrix P with zeros, you define the boundaries conditions numerically and the size of P(M,N), then you implement the loop:
C1=0.25*(Ux(x-1,y)+Ux(x+1,y)+Ux(x,y-1)+Ux(x,y+1));
C3=0.25*(Uy(x-1,y)+Uy(x+1,y)+Uy(x,y-1)+Uy(x,y+1));
C2=2*(0.5^2)*(Ux(x-1,y)+Ux(x+1,y))*(Uy(x,y-1)*Uy(x,y+1));
P(x,y)=0.25*(P(x-1,y)+P(x+1,y)+P(x,y-1)+P(x,y+1))+rho*(C1+C2+C3);

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!