% This file can be downloaded at % http://www.tfd.chalmers.se/~lada/projects/proind.html % choose "Synthetic Inlet Boundary Conditions" on the left side % % % exemple of 1d Channel flow with a PDH k-omegaa model [1]. Re=u_tau*h/nu=8000 (h=half % channel height). % % Discretization described in detail in % http://www.tfd.chalmers.se/~lada/comp_fluid_dynamics/ % % [1} S-H Peng and L. Davidson and S. Holmberg", % "A Modified Low-{R}eynolds-Number k-omega Model for Recirculating Flows", % Journal of Fluids Engng, vol. 119, pp. 867-875, 1997 close all clear % max number of iterations niter=30000; % friction velocity u_*=1 % half channel width=1 % % create the grid nj=98; njm1=nj-1; yfac=1.15; dy=0.1; yc(1)=0; for j=2:nj/2 yc(j)=yc(j-1)+dy; dy=yfac*dy; end ymax=yc(nj/2); % cell faces for j=2:nj/2 yc(j)=yc(j)/ymax; yc(nj+1-j)=2-yc(j-1); end yc(nj)=0; % cell centres for j=2:nj-1 yp(j)=0.5*(yc(j)+yc(j-1)); end yp(nj)=yc(njm1); % viscosity viscos=1/8000; % under-relaxation urf=0.5; % plot k for each iteration at node jmon jmon=10; % turbulent constants sigma_k=0.8; sigma_om=1.35; betas=0.09; gama0 = 0.42; beta = 0.075; cr1=0.75 small=1e-10; great=1e10; % init u, k & omega k(2:nj)=1; u(2:njm1)=1; om(2:nj)=10; vist(1:nj)=0; % do a loop over all nodes (except the boundary nodes) for j=2:nj-1 % compute dy_s dy_s=yp(j)-yp(j-1); % compute dy_n dy_n=yp(j+1)-yp(j); % compute deltay delta_y(j)=yc(j)-yc(j-1); dn(j)=1/dy_n; ds(j)=1/dy_s; % interpolation factor del1=yc(j)-yp(j); del2=yp(j+1)-yc(j); fy(j)=del1/(del1+del2); end u(1)=0; vist(1)=0; vist(nj)=0; % do max. niter iterations for n=1:niter for j=2:nj-1 % compute turbulent viscosity vist_old=vist(j); ret=k(j)/(viscos*om(j)); ret=max(ret,1.e-5); fmu=0.025+(1-exp(-(ret/10.0)^0.75))*(0.975+(1.e-3/ret)*exp(-(ret/200.)^2)); vist(j)=urf*fmu*k(j)/om(j)+(1-urf)*vist_old; end % solve u for j=2:nj-1 % driving pressure gradient su(j)=delta_y(j); sp(j)=0; % interpolate turbulent viscosity to faces vist_n=fy(j)*vist(j+1)+(1-fy(j))*vist(j); vist_s=fy(j-1)*vist(j)+(1-fy(j-1))*vist(j-1); % compute an & as an(j)=(vist_n+viscos)*dn(j); as(j)=(vist_s+viscos)*ds(j); end % boundary conditions for u u(1)=0; u(nj)=0; for j=2:nj-1 % compute ap ap(j)=an(j)+as(j)-sp(j); % under-relaxate ap(j)= ap(j)/urf; su(j)= su(j)+(1-urf)*ap(j)*u(j); % use Gauss-Seidel u(j)=(an(j)*u(j+1)+as(j)*u(j-1)+su(j))/ap(j); end % use TDrMmA % u=tdma(ap,an,as,su,u,nj); % monitor the development of u_tau in node jmon tau_w(n)=viscos*u(2)/yp(2); % print iteration info disp(['Iteration number: ',num2str(n)]) % check for convergence (when converged, the wall shear stress must be one) if abs(tau_w(n)-1) < 0.001 % do at least 1000 iter if n > 1000 disp(' ') disp('Converged!') break end end % solve k dudy=gradient(u,yp); dudy2=dudy.^2; for j=2:nj-1 % production term su(j)=vist(j)*dudy2(j)*delta_y(j); % dissipation term ret=k(j)/(viscos*om(j)); ret=max(ret,1.e-5); fk=1-0.722*exp(-(ret/10.)^4); sp(j)=-fk*betas*om(j)*delta_y(j); % compute an & as vist_n=fy(j)*vist(j+1)+(1-fy(j))*vist(j); an(j)=(vist_n/sigma_k+viscos)*dn(j); vist_s=fy(j-1)*vist(j)+(1-fy(j-1))*vist(j-1); as(j)=(vist_s/sigma_k+viscos)*ds(j); end % boundary conditions for k k(1)=0; k(nj)=0; for j=2:nj-1 % compute ap ap(j)=an(j)+as(j)-sp(j); % under-relaxate ap(j)= ap(j)/urf; su(j)= su(j)+(1-urf)*ap(j)*k(j); % use Gauss-Seidel k(j)=(an(j)*k(j+1)+as(j)*k(j-1)+su(j))/ap(j); end % use TDmA % k=tdma(ap,an,as,su,k,nj); % monitor the development of k in node jmon k_iter(n)=k(jmon); %****** solve om-eq. dkdy=gradient(k,yp); domdy=gradient(om,yp); for j=2:nj-1 % compute an & as vist_n=fy(j)*vist(j+1)+(1-fy(j))*vist(j); an(j)=(vist_n/sigma_om+viscos)*dn(j); vist_s=fy(j-1)*vist(j)+(1-fy(j-1))*vist(j-1); as(j)=(vist_s/sigma_om+viscos)*ds(j); ret=k(j)/(viscos*om(j)); ret=max(ret,1.e-5); fmu=0.025+(1-exp(-(ret/10.0)^0.75))*(0.975+(1.e-3/ret)*exp(-(ret/200.)^2)); fomega=1+4.3*exp(-(ret/1.5)^0.5); gama = gama0*fomega; % production term su(j)=gama*fmu*dudy2(j)*delta_y(j); % dissipation term sp(j)=-beta*om(j)*delta_y(j); % cross diffusion term crosv=dkdy(j)*domdy(j); crost=cr1*vist(j)*crosv/k(j)*delta_y(j); su(j)=su(j)+max(crost,0); sp(j)=sp(j)+min(crost,0)/om(j); end % b.c. south wall dy=yp(2); omega=6.*viscos/0.075/dy^2; sp(2)=-great; su(2)=great*omega; % b.c. north wall dy=yp(nj)-yp(njm1); omega=6.*viscos/0.075/dy^2; sp(njm1)=-great; su(njm1)=great*omega; for j=2:nj-1 % compute ap ap(j)=an(j)+as(j)-sp(j); % under-relaxate ap(j)= ap(j)/urf; su(j)= su(j)+(1-urf)*ap(j)*om(j); % use Gauss-Seidel % om(j)=(an(j)*om(j+1)+as(j)*om(j-1)+su(j))/ap(j); end % use TDmA om=tdma(ap,an,as,su,om,nj); om_iter(n)=om(jmon); % flush the buffer (in octave) % fflush(stdout) end % compute shear stress uv=-vist.*dudy; % plot u %load ~/calc_channel/PDH-re8000/u_8000_PDH.dat figure(1) tau_turb=viscos*u(2)/yp(2); plot(u,yp,'linewidth',2) hold plot(u_8000_PDH(:,2),yp(2:njm1),'o') xlabel('u') ylabel('y') legend('PDH','Matlab') print u_8000.ps -deps % plot k %load ~/calc_channel/PDH-re8000/k_8000_PDH.dat figure(2) plot(k,yp,'linewidth',2) hold plot(k_8000_PDH,yp(2:njm1),'o') xlabel('k') ylabel('y') legend('PDH','Matlab') print k_8000.ps -deps % plot tau_w versus iteration number figure(3) plot(tau_w,'linewidth',2) title('wall shear stress') xlabel('Iteration number') ylabel('tauw') print tauw.ps -deps % plot k(jmon) versus iteration number figure(4) plot(k_iter,'linewidth',2) title('k in node jmon') xlabel('Iteration number') ylabel('k') print k_iter.ps -deps % plot om(jmon) versus iteration number figure(5) plot(om_iter,'linewidth',2) title('omega in node jmon') xlabel('Iteration number') ylabel('omega') print om_iter.ps -deps % save data ff(:,1)=yp; ff(:,2)=u; ff(:,3)=k; ff(:,4)=om; ff(:,5)=vist; ff(:,6)=uv; ff(:,7)=yc; %save yp_u_k_om_vist_uv_yc_PDH_8000.dat ff -ASCII save yp_u_k_om_vist_uv_yc_PDH_8000.dat ff -ascii