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: 7 Sponsored Links 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;

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Perumal Main CFD Forum 1 January 22, 2007 14:27 ani Main CFD Forum 3 April 2, 2006 13:27 diaw Main CFD Forum 104 February 16, 2006 06:44 Axel Rohde Main CFD Forum 1 November 19, 2001 13:19 ram Main CFD Forum 5 June 17, 2000 21:31