PISO algorithm
2 Attachment(s)
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!!! |
Quote:
|
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, |
matlab PISO algorithm
Quote:
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 |
|
Quote:
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. |
All times are GMT -4. The time now is 09:47. |