2D_Matalb_code_twophase_flow.rar_2d flow_FVM;matlab_cfd matlab_m

%======================================================================== % A very simple Navier-Stokes solver for a drop falling in a rectangular % box. A forward in time, centered in space discretization is used. % The density is advected by a front tracking scheme. % This version uses a stretched grid and ENO for the advection %======================================================================== %domain size and physical variables Lx=1.0;Ly=1.0;gx=0.0;gy=100.0; rho1=1.0; rho2=2.0; m1=0.01; m2=0.05; sigma=10; rro=rho1; un=0;us=0;ve=0;vw=0;time=0.0; rad=0.15;xc=0.5;yc=0.7; % Initial drop size and location % Numerical variables nx=32;ny=32;dt=0.00125;nstep=300;maxit=200;maxError=0.001;beta=1.2; % Zero various arrys u=zeros(nx+1,ny+2); v=zeros(nx+2,ny+1); p=zeros(nx+2,ny+2); ut=zeros(nx+1,ny+2); vt=zeros(nx+2,ny+1); tmp1=zeros(nx+2,ny+2); uu=zeros(nx+1,ny+1); vv=zeros(nx+1,ny+1); tmp2=zeros(nx+2,ny+2); fx=zeros(nx+2,ny+2); fy=zeros(nx+2,ny+2); dk=zeros(1,nx+2); dl=zeros(1,ny+2); dkh=zeros(1,nx+1); dlh=zeros(1,ny+1); ux=zeros(nx+2,ny+2); uy=zeros(nx+2,ny+2); vx=zeros(nx+2,ny+2); vy=zeros(nx+2,ny+2); % Set the stretched grid for i=1:nx+2, s=Lx*(i-1.5)/(nx); x(i)=1.5*s*(0.5-s)*(1-s)+s; end; for j=1:ny+2, s=Ly*(j-1.5)/(ny); y(j)=1.5*s*(0.5-s)*(1-s)+s; end; for i=1:nx+1, xh(i)=0.5*(x(i+1)+x(i)); end; for j=1:ny+1, yh(j)=0.5*(y(j+1)+y(j)); end for i=1:nx+1, dkh(i)=x(i+1)-x(i); end; for j=1:ny+1, dlh(j)=y(j+1)-y(j); end; for i=2:nx+1, dk(i)=xh(i)-xh(i-1); end; for j=2:ny+1, dl(j)=yh(j)-yh(j-1); end; dk(1)=dk(2);dk(nx+2)=dk(nx+1); dl(1)=dl(2);dl(ny+2)=dl(ny+1); % Set density and viscosity in the domain and the drop r=zeros(nx+2,ny+2)+rho1;m=zeros(nx+2,ny+2)+m1; for i=2:nx+1,for j=2:ny+1; if ( (x(i)-xc)^2+(y(j)-yc)^2 < rad^2), r(i,j)=rho2;m(i,j)=m2; end, end,end %================== SETUP THE FRONT =================== Nf=100; xf=zeros(1,Nf+2);yf=zeros(1,Nf+2); uf=zeros(1,Nf+2);vf=zeros(1,Nf+2); tx=zeros(1,Nf+2);ty=zeros(1,Nf+2); for l=1:Nf+2, xf(l)=xc-rad*sin(2.0*pi*(l-1)/(Nf)); yf(l)=yc+rad*cos(2.0*pi*(l-1)/(Nf));end %------------find the mapped front coordinates-------- sxf=zeros(1,Nf+2);syf=zeros(1,Nf+2); for l=1:Nf+2, for i=1:nx+1 if xf(l) > x(i+1), % DO NOTHING else sxf(l)=i*(x(i+1)-xf(l))/dkh(i)+(i+1)*(xf(l)-x(i))/dkh(i);break end end for j=1:ny+1 if yf(l) >= y(j+1), % DO NOTHING else syf(l)=j*(y(j+1)-yf(l))/dlh(j)+(j+1)*(yf(l)-y(j))/dlh(j);break end end end %================== START TIME LOOP====================================== for is=1:nstep,is %------------------ FIND SURFACE TENSION -------------- fx=zeros(nx+2,ny+2);fy=zeros(nx+2,ny+2); % Set fx & fy to zero for l=1:Nf+1, ds=sqrt((xf(l+1)-xf(l))^2+(yf(l+1)-yf(l))^2); tx(l)=(xf(l+1)-xf(l))/ds; ty(l)=(yf(l+1)-yf(l))/ds; % Tangent vectors end tx(Nf+2)=tx(2);ty(Nf+2)=ty(2); for l=2:Nf+1 nfx=sigma*(tx(l)-tx(l-1));nfy=sigma*(ty(l)-ty(l-1)); ip=floor(sxf(l)-0.5);jp=floor(syf(l)); ax=sxf(l)-0.5-ip;ay=syf(l)-jp; fx(ip,jp) =fx(ip,jp)+(1.0-ax)*(1.0-ay)*nfx/dkh(ip)/dl(jp); fx(ip+1,jp) =fx(ip+1,jp)+ax*(1.0-ay)*nfx/dkh(ip+1)/dl(jp); fx(ip,jp+1) =fx(ip,jp+1)+(1.0-ax)*ay*nfx/dkh(ip)/dl(jp+1); fx(ip+1,jp+1)=fx(ip+1,jp+1)+ax*ay*nfx/dkh(ip+1)/dl(jp+1); ip=floor(sxf(l));jp=floor(syf(l)-0.5); ax=sxf(l)-ip;ay=syf(l)-0.5-jp; fy(ip,jp) =fy(ip,jp)+(1.0-ax)*(1.0-ay)*nfy/dk(ip)/dlh(jp); fy(ip+1,jp) =fy(ip+1,jp)+ax*(1.0-ay)*nfy/dk(ip+1)/dlh(jp); fy(ip,jp+1) =fy(ip,jp+1)+(1.0-ax)*ay*nfy/dk(ip)/dlh(jp+1); fy(ip+1,jp+1)=fy(ip+1,jp+1)+ax*ay*nfy/dk(ip+1)/dlh(jp+1); end %--------------------------------------------------------- % tangential velocity at boundaries u(1:nx+1,1)=2*us-u(1:nx+1,2);u(1:nx+1,ny+2)=2*un-u(1:nx+1,ny+1); v(1,1:ny+1)=2*vw-v(2,1:ny+1);v(nx+2,1:ny+1)=2*ve-v(nx+1,1:ny+1); for i=2:nx-1,for j=2:ny+1; % finding the ENO values for the face velocities ss=0.5*sign(u(i+1,j)+u(i,j)); ux(i+1,j)=(0.5+ss)*(u(i,j)+0.5*min(abs(u(i+1,j)-u(i,j)),abs(u(i,j)-u(i-1,j))) )... +(0.5-ss)*(u(i+1,j)-0.5*min(abs(u(i+2,j)-u(i+1,j)),abs(u(i+1,j)-u(i,j))) ); end,end ux(2,2:ny+1)=0.5*u(2,2:ny+1)+u(1,2:ny+1); ux(nx+1,2:ny+1)=0.5*u(nx,2:ny+1)+u(nx+1,2:ny+1); %Bdry for i=2:nx,for j=3:ny; ss=0.5*sign(v(i+1,j)+v(i,j)); uy(i,j+1)=(0.5+ss)*(u(i,j)+0.5*min(abs(u(i,j+1)-u(i,j)),abs(u(i,j)-u(i,j-1))) )... +(0.5-ss)*(u(i,j+1)-0.5*min(abs(u(i,j+2)-u(i,j+1)),abs(u(i,j+1)-u(i,j))) ); end,end uy(1:nx+1,2)=0.0;uy(1:nx+1,ny+1)=0.0; % Bottom and Top for i=2:nx,for j=2:ny+1 % TEMPORARY u-velocity-ADVECTION ut(i,j)=u(i,j)+dt*(-(u(i,j)*(ux(i+1,j)-ux(i,j))/dkh(i)+... 0.25*(v(i,j-1)+v(i+1,j-1)+v(i+1,j)+v(i,j))*(uy(i,j+1)-uy(i,j))/dl(j) )+... fx(i,j)/(0.5*(r(i+1,j)+r(i,j))) ... - (1.0 -rro/(0.5*(r(i+1,j)+r(i,j))) )*gx ); end,end for i=3:nx,for j=2:ny; % finding the ENO values for the face velocities ss=0.5*sign(u(i,j+1)+u(i,j)); vx(i+1,j)=(0.5+ss)*(v(i,j)+0.5*min(abs(v(i+1,j)-v(i,j)),abs(v(i,j)-v(i-1,j))) )... +(0.5-ss)*(v(i+1,j)-0.5*min(abs(v(i+2,j)-v(i+1,j)),abs(v(i+1,j)-v(i,j))) ); end,end vx(2,1:ny+1)=0.0;vx(nx+1,1:ny+1)=0.0; % Left and Right for i=2:nx+1,for j=2:ny-1; ss=0.5*sign(v(i,j+1)+v(i,j)); vy(i,j+1)=(0.5+ss)*(v(i,j)+0.5*min(abs(v(i,j+1)-v(i,j)),abs(v(i,j)-v(i,j-1))) )... +(0.5-ss)*(v(i,j+1)-0.5*min(abs(v(i,j+2)-v(i,j+1)),abs(v(i,j+1)-v(i,j))) ); end,end vy(2:nx+1,2)=0.5*v(2:nx+1,2)+v(2:nx+1,1); vy(2:nx+1,ny+1)=0.5*v(2:nx+1,ny)+v(2:nx+1,ny+1); %Bdry for i=2:nx+1,for j=2:ny % TEMPORARY v-velocity-ADVECTION vt(i,j)=v(i,j)+dt*(-(0.25*(u(i-1,j)+u(i,j)+u(i,j+1)+... u(i-1,j+1))*(vx(i+1,j)-vx(i,j))/dk(i)+ ... v(i,j)*(vy(i,j+1)-vy(i,j))/dlh(j))+ ... fy(i,j)/(0.5*(r(i,j+1)+r(i,j))) ... - (1.0 -rro/(0.5*(r(i,j+1)+r(i,j))) )*gy ); end,end for i=2:nx,for j=2:ny+1 % TEMPORARY u-velocity-DIFFUSION ut(i,j)=ut(i,j)+dt*(... (1./dkh(i))*2.*(m(i+1,j)*(1./dk(i+1))*(u(i+1,j)-u(i,j)) - ... m(i,j) *(1./dk(i))*(u(i,j)-u(i-1,j)) ) ... +(1./dl(j))*( 0.25*(m(i,j)+m(i+1,j)+m(i+1,j+1)+m(i,j+1))* ... ((1./dlh(j))*(u(i,j+1)-u(i,j)) + (1./dkh(i))*(v(i+1,j)-v(i,j)) ) - ... 0.25*(m(i,j)+m(i+1,j)+m(i+1,j-1)+m(i,j-1))* ... ((1./dlh(j-1))*(u(i,j)-u(i,j-1))+ (1./dkh(i-1))*(v(i+1,j-1)- v(i,j-1))) )... )/(0.5*(r(i+1,j)+r(i,j)) ); end,end for i=2:nx+1,for j=2:ny % TEMPORARY v-velocity-DIFFUSION vt(i,j)=vt(i,j)+dt*(... (1./dk(i))*( 0.25*(m(i,j)+m(i+1,j)+m(i+1,j+1)+m(i,j+1))* ... ((1./dlh(j))*(u(i,j+1)-u(i,j)) + (1./dkh(i))*(v(i+1,j)-v(i,j)) ) - ... 0.25*(m(i,j)+m(i,j+1)+m(i-1,j+1)+m(i-1,j))* ... ((1./dlh(j))*(u(i-1,j+1)-u(i-1,j))+ (1./dkh(i-1))*(v(i,j)- v(i-1,j))) )... +(1./dlh(j))*2.*(m(i,j+1)*(1./dl(j+1))*(v(i,j+1)-v(i,j)) - ... m(i,j) *(1./dl(j))*(v(i,j)-v(i,j-1)) ) ... )/(0.5*(r(i,j+1)+r(i,j)) ); end,end %======================================================================== % Compute source term and the coefficient for p(i,j) rt=r; lrg=1000; rt(1:nx+2,1)=lrg;rt(1:nx+2,ny+2)=lrg; rt(1,1:ny+2)=lrg;rt(nx+2,1:ny+2)=lrg; for i=2:nx+1,for j=2:ny+1 tmp1(i,j)= (0.5/dt)*( (ut(i,j)-ut(i-1,j))/dk(i)+(vt(i,j)-vt(i,j-1))/dl(j) ); tmp2(i,j)=1.0/( (1./dk(i))*( 1./(dkh(i