[ari92]

Hello everybody!

I have to find the solution for the velocity and pressure using PISO algorithm in Matlab, but I don't know how how to obtain b'' for the boundary and next to-boundary pressure nodes. could someone help me??
thank you very much!!!

[lcarasik]

Originally Posted by ari92

Hello everybody!

I have to find the solution for the velocity and pressure using PISO algorithm in Matlab, but I don't know how how to obtain b'' for the boundary and next to-boundary pressure nodes. could someone help me??
thank you very much!!!

What have you worked on so far? It is easier for forum members to help if you can show your thinking.

[Roger9318]

Hi, ari92
right now I´m in the same situation and I´m wondering if you could you solve your dobut?

could you share with me?
best regards,

[ari92]

Originally Posted by Roger9318

Hi, ari92
right now I´m in the same situation and I´m wondering if you could you solve your dobut?

could you share with me?
best regards,

Hi Roger9318,
I can send you the copy of my matlab code..you can take the information you need from it..hoping that is good for you:

while (r_cont>tol_cont || r_u>tol_u)
%%% Solve discretised momentum equations %%%
a=zeros(n,n);
S=zeros(n,1);
u_star=zeros(n,1);
d=zeros(n,1);
%node 1
a(1,1)=rho*Ap(2)*(u_guess(1)+u_guess(2))/2 + 1/2*rho*u_guess(1)*Av(1)*(Av(1)/Ap(1))^2;
S(1)=(p0-p_guess(2))*Av(1) + rho*u_guess(1)*Av(1)*(Av(1)/Ap(1))*u_guess(1);
d(1)=Av(1)/a(1,1);
%node i
for i=2:n-1
a(i,i)=rho*Ap(i+1)*(u_guess(i)+u_guess(i+1))/2;
a(i,i-1)=-(rho*Ap(i)*(u_guess(i-1)+u_guess(i))/2);
S(i)=(p_guess(i)-p_guess(i+1))*Av(i);
d(i)=Av(i)/a(i,i);
end
%node n
a(n,n)=rho*Av(n)*u_guess(n);
a(n,n-1)=-(rho*Ap(N-1)*(u_guess(n-1)+u_guess(n))/2);
S(n)=(p_guess(N-1)-p_out)*Av(n);
d(n)=Av(n)/a(n,n);
u_star=a\S; % first guess velocity (u*)
%%% FIRST CORRECTION STEP: Solve pressure correction equation %%%
a1=zeros(N,N);
p_corr=zeros(N,1);
b1=zeros(N,1);
p_corr(1)=0;
p_corr(N)=0;
%node A
a1(1,1)=1;
a1(1,2)=0;
b1(1)=0;
%node i
for i=2:N-1
a1(i,i-1)=-(rho*d(i-1)*Av(i-1));
a1(i,i+1)=-(rho*d(i)*Av(i));
a1(i,i)=rho*d(i-1)*Av(i-1) + rho*d(i)*Av(i);
b1(i)=rho*Av(i-1)*u_star(i-1) - rho*Av(i)*u_star(i);
end
%node N
a1(N,N)=1;
a1(N,N-1)=0;
b1(N)=0;
p_corr=a1\b1; % first pressure correction (p')
%%% Calculate u** and p** %%%
u_2star=zeros(n,1); %second guess velocity (u**)
for i=1:n
u_2star(i)=u_star(i)+d(i)*(p_corr(i)-p_corr(i+1));
end
p_2star=zeros(N,1); % second guess pressure (p**)
p_2star(1)=p0-(1/2)*rho*u_2star(1)^2*(Av(1)/Ap(1))^2;
p_2star(n)=p_out;
for i=2:N-1
p_2star(i)=p_guess(i)+p_corr(i);
end
%%% Second corrector step - Solve second pressure correction equation %%%
a2=zeros(N,N);
p_2corr=zeros(N,1);
b2=zeros(N,1);
p_2corr(1)=0;
p_2corr(N)=0;
%node A
a2(1,1)=1;
a2(1,2)=0;
b2(1)=0;
% node B
a2(2,1)=-rho*d(1)*Av(1);
a2(2,3)=-rho*d(2)*Av(2);
a2(2,2)=rho*d(1)*Av(1) + rho*d(2)*Av(2);
b2(2)=-(rho*Av(1)/a(1,1))*(a(1,2)*(u_2star(2)-u_star(2)))+ (rho*Av(2)/a(2,2))*(a(2,1)*(u_2star(1)-u_star(1))+a(2,3)*(u_2star(3)-u_star(3)));
%node i
for i=3:N-2
a2(i,i-1)=-(rho*d(i-1)*Av(i-1));
a2(i,i+1)=-(rho*d(i)*Av(i));
a2(i,i)=rho*d(i-1)*Av(i-1) + rho*d(i)*Av(i);
b2(i)=-(rho*Av(i-1)/a(i-1,i-1))*(a(i-1,i-2)*(u_2star(i-2)-u_star(i-2))+a(i-1,i)*(u_2star(i)-u_star(i)))+ (rho*Av(i)/a(i,i))*(a(i,i-1)*(u_2star(i-1)-u_star(i-1))+a(i,i+1)*(u_2star(i+1)-u_star(i+1)));
end
%node N
a2(N,N)=1;
a2(N,N-1)=0;
b2(N)=0;
%node N-1
a2(N-1,N-2)=-(rho*d(n-1)*Av(n-1));
a2(N-1,N)=-(rho*d(n)*Av(n));
a2(N-1,N-1)=rho*d(n-1)*Av(n-1)+rho*d(n)*Av(n);
b2(N-1)=-(rho*Av(n-1)/a(n-1,n-1))*(a(n-1,n-2)*(u_2star(n-2)-u_star(n-2))+a(n-2,n)*(u_2star(n)-u_star(n)))+ (rho*Av(n)/a(n,n))*(a(n,n-1)*(u_2star(n-1)-u_star(n-1)));
p_2corr=a2\b2; % second pressure correction (p'')
%%% correct pressure and velocity %%%
U=zeros(n,1);
U(1)=u_2star(1)+(1/a(1,1))*(a(1,2)*(u_2star(2)-u_star(2))) +d(1)*(p_2corr(1)-p_2corr(2));
U(n)=u_2star(n)+(1/a(n,n))*(a(n,n-1)*(u_2star(n-1)-u_star(n-1))) +d(n)*(p_2corr(n)-p_2corr(n+1));
for i=2:n-1
U(i)=u_2star(i)+(1/a(i,i))*((a(i,i-1)*(u_2star(i-1)-u_star(i-1))+a(i,i+1)*(u_2star(i+1)-u_star(i+1)))) +d(i)*(p_2corr(i)-p_2corr(i+1));
end
P=zeros(N,1);
P(1)=p0-(1/2)*rho*U(1)^2*(Av(1)/Ap(1))^2;
P(n)=p_out;
for i=2:N-1
P(i)=p_2star(i)+p_2corr(i);
end
%%%% Residuals for momentum equation
r_u=norm(a*U-S)./norm(diag(a).*U);
%%%% Residuals for continuity
r_cont=sum(abs(b2));
residual_u(end+1)=r_u;
residual_cont(end+1)=r_cont;
%%% Under-relaxation
alfa_u=0.5;
alfa_p=0.5;
U_rel=(1-alfa_u)*u_guess+alfa_u*U;
P_rel=(1-alfa_p)*p_guess+alfa_p*P;
u_guess=U_rel;
p_guess=P_rel;
k=k+1;
end

[Roger9318]

ari92, thank you for share your code, now I'm going to look at it.



best regards!!!

[midoronald]

Originally Posted by ari92

Hi Roger9318,
I can send you the copy of my matlab code..you can take the information you need from it..hoping that is good for you:



%%% Under-relaxation
alfa_u=0.5;
alfa_p=0.5;
U_rel=(1-alfa_u)*u_guess+alfa_u*U;
P_rel=(1-alfa_p)*p_guess+alfa_p*P;
u_guess=U_rel;
p_guess=P_rel;
k=k+1;
end

I'm wondering why you are using under-relaxation factor while using PISO algorithm. It is supposed to be a non-iterative algorithm and if you want to reach the steady state case, you just have to wait until the flow field doesn't change any more at t=T_steady_state.

Also, you don't have to add your P_rel for each time step. The algorithm is already dealing with the pressure differences and the pressure correction equation at the boundaries are supposed to be 0 at the outflow and to be calculated in the inflow, in case of fixing the total pressure at the inflow. For some cases, the velocities are explicitly given at the inflow e.g., the classical open channel flow. In that case the pressure needs to be calculated at the inflow boundary and approximately equal to zero just like the SIMPLE algorithm boundary condition for the pressure Poisson equation by declaring aw=0, assuming the inflow is on the left hand side of the channel.