# Lid driven flow in shallow rectangular domain issues, MATLAB

 Register Blogs Members List Search Today's Posts Mark Forums Read

 November 10, 2011, 21:41 Lid driven flow in shallow rectangular domain issues, MATLAB #1 Member   Peter Join Date: Oct 2011 Posts: 52 Rep Power: 5 So i'm modifying the code made by Benjamin Seibold to deal with shallow rectangular cavities for lid driven flows. Original code is located here and is 6.7 under MATLAB codes along with associated documentation. http://www-math.mit.edu/cse/ My problem is that the pressure correction step is not enforcing a divergence free velocity field. To change to a rectangular domain, I have non dimensionalized using the following: which gives the N-S equations Here is the modified code that I have so far. It plots the stream lines, the u and v velocity at different vertical slices of the rectangle and prints the maximum divergence in the field. if you feel like running it, just copy paste to matlab m file, save and hit F5 and everything should work. By adjusting AR at the top, you can change the size of the rectangle where AR = H/L If the AR is 1, code runs fine and quantitative values have been compared with ghia and do not show any significant errors. Any one know whats wrong?????? Thanks Code: function []=debug() clear all close all clc %% Paramters AR=0.5; % Aspect Ratio Pe=500; % Peclet Number Re = 100; % Reynolds number tf = 40e-0; % final time dt = 0.01; % time step lx = 1; % width of non-D domain ly = 1; % height of nond-D domain hxmax=0.025; % Maximum x grid size hymax=0.025; % Maximum y grid size ny=ceil(1/hymax); % Number of y grid points nx=ceil(1/AR/hxmax); % Number of x grid points hy=1/ny; % Actual y grid size hx=1/nx; % Actual x grid size % Dimensional Variables H=1; % Height of droplet L=H/AR; % Length of droplet dy=H/ny; % Width of y grid dx=L/nx; % Width of x grid dimx(1:nx+1)=(0:nx)*dx; % X coordinates dimy(1:ny+1)=(0:ny)*dy; % Y coordinates % Fluid Properties rho=1000; % Density mu=0.86e-3; % Viscosity c=4184; % Specific heat k=0.6; % Thermal conductivity ubulk=(Re*mu)/(rho*H); % Bulk Velocity nsteps = 10; % number of steps with graphic output %----------------------------------------------------------------------- nt = ceil(tf/dt); dt = tf/nt; % number of time iterations x = linspace(0,lx,nx+1); % non-D x coordinates y = linspace(0,ly,ny+1); % non-D y coordinates [X,Y] = meshgrid(y,x); %----------------------------------------------------------------------- % initial conditions U = zeros(nx-1,ny); V = zeros(nx,ny-1); % boundary conditions uN = x*0+1; uN = [0 uN(2:end-1) 0]; uS = x*0; uS = [0 uS(2:end-1) 0]; uW = avg(y)*0; uE = avg(y)*0; vN = avg(x)*0; vS = avg(x)*0; vW = y*0; vE = y*0; %% Time stepping Ubc = dt/Re*([2*uS(2:end-1)' zeros(nx-1,ny-2) 2*uN(2:end-1)']/hy^2+... [uW;zeros(nx-3,ny);uE]/hx^2); Vbc = dt/Re*([vS' zeros(nx,ny-3) vN']/hy^2+... [2*vW(2:end-1);zeros(nx-2,ny-1);2*vE(2:end-1)]/hx^2); fprintf('initialization') Lp = kron(speye(ny),AR*K1(nx,hx,1))+kron((1/AR)*K1(ny,hy,1),speye(nx)); Lp(1,1) = 3/2*Lp(1,1); perp = symamd(Lp); Rp = chol(Lp(perp,perp)); Rpt = Rp'; Lu = speye((nx-1)*ny)+dt/Re*(kron(speye(ny),AR^2*K1(nx-1,hx,2))+... kron(K1(ny,hy,3),speye(nx-1))); peru = symamd(Lu); Ru = chol(Lu(peru,peru)); Rut = Ru'; Lv = speye(nx*(ny-1))+dt/Re*(kron(speye(ny-1),AR^2*K1(nx,hx,3))+... kron(K1(ny-1,hy,2),speye(nx))); perv = symamd(Lv); Rv = chol(Lv(perv,perv)); Rvt = Rv'; Lq = kron(speye(ny-1),K1(nx-1,hx,2))+kron(K1(ny-1,hy,2),speye(nx-1)); perq = symamd(Lq); Rq = chol(Lq(perq,perq)); Rqt = Rq'; fprintf(', time loop\n--20%%--40%%--60%%--80%%-100%%\n') for k = 1:nt % treat nonlinear terms gamma = min(1.2*dt*max(max(max(abs(U)))/hx,max(max(abs(V)))/hy),1); Ue = [uW;U;uE]; Ue = [2*uS'-Ue(:,1) Ue 2*uN'-Ue(:,end)]; Ve = [vS' V vN']; Ve = [2*vW-Ve(1,:);Ve;2*vE-Ve(end,:)]; Ua = avg(Ue')'; Ud = diff(Ue')'/2; Va = avg(Ve); Vd = diff(Ve)/2; UVx = diff(Ua.*Va-gamma*abs(Ua).*Vd)/hx; UVy = diff((Ua.*Va-gamma*Ud.*abs(Va))')'/hy; Ua = avg(Ue(:,2:end-1)); Ud = diff(Ue(:,2:end-1))/2; Va = avg(Ve(2:end-1,:)')'; Vd = diff(Ve(2:end-1,:)')'/2; U2x = diff(Ua.^2-gamma*abs(Ua).*Ud)/hx; V2y = diff((Va.^2-gamma*abs(Va).*Vd)')'/hy; U = U-dt*AR*(UVy(2:end-1,:)+U2x); V = V-dt*AR*(UVx(:,2:end-1)+V2y); % implicit viscosity rhs = reshape(U+Ubc,[],1); u(peru) = Ru\(Rut\rhs(peru)); U = reshape(u,nx-1,ny); rhs = reshape(V+Vbc,[],1); v(perv) = Rv\(Rvt\rhs(perv)); V = reshape(v,nx,ny-1); % pressure correction % rhsv = reshape((diff([uW;U;uE])/hx+diff([vS' V vN']')'/hy),[],1); % rhsv = rhsv.*AR^2; % rhsu = reshape((diff([uW;U;uE])/hx+diff([vS' V vN']')'/hy),[],1); % pu(perp) = -Rp\(Rpt\rhsu(perp)); % pv(perp) = -Rp\(Rpt\rhsv(perp)); % Pu = reshape(pu,nx,ny); % Pv = reshape(pv,nx,ny); % U = U-AR*diff(Pu)/hx; % V = V-(1/AR)*diff(Pv')'/hy; rhs = reshape(diff([uW;U;uE])/hx+diff([vS' V vN']')'/hy,[],1); p(perp) = -Rp\(Rpt\rhs(perp)); P = reshape(p,nx,ny); U = U-diff(P)/hx; V = V-diff(P')'/hy; % visualization if floor(25*k/nt)>floor(25*(k-1)/nt), fprintf('.'), end if k==1|floor(nsteps*k/nt)>floor(nsteps*(k-1)/nt) % stream function rhs = reshape(diff(U')'/hy-diff(V)/hx,[],1); q(perq) = Rq\(Rqt\rhs(perq)); Q = zeros(nx+1,ny+1); Q(2:end-1,2:end-1) = reshape(q,nx-1,ny-1); % clf, contourf(avg(x),avg(y),P',20,'w-'), hold on contourf(dimx,dimy,Q',20,'w-'); Ue = [uS' avg([uW;U;uE]')' uN']; Ve = [vW;avg([vS' V vN']);vE]; Len = sqrt(Ue.^2+Ve.^2+eps); % quiver(x,y,(Ue./Len)',(Ve./Len)',.4,'k-') hold off, axis equal, axis([0 L 0 H]) % p = sort(p); caxis(p([8 end-7])) title(sprintf('Re = %0.1g t = %0.2g',Re,k*dt)) drawnow end end fprintf('\n') ufield=Ue'; vfield=Ve'; ulab=abs(ufield-1); mid=ceil((nx+1)/2); quart=ceil(1*(nx+1)/100); if quart==1 quart=2; end third=ceil(5*(nx+1)/100); next=ceil(1*(nx+1)/5); figure(3) plot(ulab(:,quart),y,ulab(:,third),y,ulab(:,next),y,ulab(:,mid),y)%,ulab(:,end-quart),y) legend('0.02H','0.1H','0.4L','H','other','location','West') title('U profiles') figure(4) plot(vfield(:,quart),y,vfield(:,third),y,vfield(:,next),y,vfield(:,mid),y)%,vfield(:,end-quart),y) legend('0.02H','0.1H','0.4L','H','other','location','West') title('V profiles') %% Divergence Test div=divergence(ufield,vfield); max_div = max(max(div(2:end-1,2:end-1))) function B = avg(A,k) if nargin<2, k = 1; end if size(A,1)==1, A = A'; end if k<2, B = (A(2:end,:)+A(1:end-1,:))/2; else, B = avg(A,k-1); end if size(A,2)==1, B = B'; end function A = K1(n,h,a11) % a11: Neumann=1, Dirichlet=2, Dirichlet mid=3; A = spdiags([-1 a11 0;ones(n-2,1)*[-1 2 -1];0 a11 -1],-1:1,n,n)'/h^2;