| new_at_this |
December 21, 2011 23:57 |
so this is code is driving me crazy. I'm going to post it up here and if someone has some time to look at it I would really appreciate it. The code is based off the lid driven cavity written by Benjamin Seibold. It has been discussed before on the forum so I'm hoping someone is familiar with it.
I have modified the code so that both the lid and the base have velocity which produces a symmetric flow. If you copy the whole code and then run "testcode(0)" in the command line it will compute the entire domain. If you run "testcode(1)" then it will run the code for only one half of the domain using a symmetry BC on the lower boundary.
Code:
function []=testcode(sym)
close all
clc
%% Paramters
% Input Parameters --------------------------------------------------------
AR=0.5; % Aspect Ratio
Pe=500; % Peclet Number
Re = 100; % Reynolds number
tf = 40e-0; % final time
gridnum=50; % Number of grid points
dt = 0.01; % time step
nsteps = 10; % number of steps with graphic output
mode=1; % 1: basic , 2: extended with more plots
% Run parameters ----------------------------------------------------------
lx = 1; % width of non-D domain
ly = 1; % height of nond-D domain
nx=gridnum; % number of x gridpoints
ny=gridnum; % number of y gridpoints
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
%-----------------------------------------------------------------------
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];
if sym == 0
uS = x*0+1; uS = [0 uS(2:end-1) 0];
end
uW = avg(y)*0;
uE = avg(y)*0;
vN = avg(x)*0;
vS = avg(x)*0;
vW = y*0;
vE = y*0;
%% Time stepping
if sym == 0
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);
end
fprintf('initialization')
if sym == 1
Lp = kron(speye(ny),AR*K1(nx,hx,1))+kron((1/AR)*K7(ny,hy,1),speye(nx));
Lp(1,1) = 3/2*Lp(1,1);
Lu = speye((nx-1)*ny)+dt/Re*(kron(speye(ny),AR^2*K1(nx-1,hx,2))+...
kron(K7(ny,hy,3),speye(nx-1)));
elseif sym == 0
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)));
end
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
if sym == 1
uS = [0 U(:,1)' 0];
Ubc = dt/Re*([0*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);
end
% 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
rhs = reshape(diff([uW;U;uE])/hx+diff([vS' V vN']')'/hy,[],1);
if sym == 0
p(perp) = -Rp\(Rpt\rhs(perp));
P = reshape(p,nx,ny);
elseif sym == 1
p2=Lp\rhs;
P2=reshape(p2,nx,ny);
P=-P2;
end
U = U-AR*diff(P)/hx;
V = V-(1/AR)*diff((P)')'/hy;
% Store U and V including ghost points
vghost=[vW; [vS' V vN'] ;vE];
ughost=[uS' [uW;U;uE] uN'];
Ue = [uS' avg([uW;U;uE]')' uN'];
Ve = [vW;avg([vS' V vN']);vE];
% 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-'), hold on
Ue = [uS' avg([uW;U;uE]')' uN'];
Ve = [vW;avg([vS' V vN']);vE];
Len = sqrt(Ue.^2+Ve.^2+eps);
quiver(dimx,dimy,(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);
%% Profiles
figure(1)
contour(dimx,dimy,Q',20,'Linewidth',1.2);
axis equal, axis([0 L 0 H])
set(gca,'FontSize',14)
xlabel('X','Fontsize',16)
ylabel('Y','Fontsize',16)
title('Droplet Streamlines (Re = 100, AR = 0.5)','Fontsize',20)
colorbar
mid=ceil((nx+1)/2);
cx1=ceil(1*(nx+1)/100);
if cx1==1
cx1=2;
end
cx2=ceil(5*(nx+1)/100);
cx3=ceil(1*(nx+1)/5);
if mode==2
figure(2)
plot(ufield(:,cx1),y,'k-',ufield(:,cx2),y,':',ufield(:,cx3),y,'m-.',ufield(:,mid),y,'--');
leg=legend('0.02H','0.1H','0.4L','H','other','location','East');
xlabel('U^*','Fontsize',16)
ylabel('Y','Fontsize',16)
title('U profiles Droplet Frame(Re = 100, AR = 0.5)','Fontsize',16)
set(leg,'FontSize',14);
end
figure(3)
uprof=plot(ulab(:,cx1),y,'k-',ulab(:,cx2),y,':',ulab(:,cx3),y,'m-.',ulab(:,mid),y,'--');
set(uprof,'LineWidth',1.5)
set(gca,'FontSize',14)
leg=legend('0.02H','0.1H','0.4L','H','other','location','West');
xlabel('U^*','Fontsize',16)
ylabel('Y','Fontsize',16)
title('U profiles Lab Frame (Re = 100, AR = 0.5)','Fontsize',16)
set(leg,'FontSize',14);
if mode==2
figure(4)
plot(vfield(:,cx1),x,vfield(:,cx2),x,vfield(:,cx3),x,vfield(:,mid),x)%,vfield(:,end-quart),y)
legend('0.02H','0.1H','0.4L','H','other','location','NorthWest')
title('V profiles along vertical cut')
figure(5)
plot(x,vfield(cx1,:),x,vfield(cx2,:),x,vfield(cx3,:),x,vfield(mid,:))%,vfield(:,end-quart),y)
legend('0.02H','0.1H','0.4L','H','other','location','NorthWest')
title('V profiles along horizontal cut')
end
%% Divergence Test
div=divergence(ufield,vfield);
max_div = max(max(div(2:end-1,2:end-1)))
%% Vorticity Test
if mode==2
[CURLZ, CAV] = curl(X,Y,ufield,vfield);
figure(6)
contourf(x,y,CURLZ,50)
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;
function A = K7(n,h,a11)
% a11: Neumann=1, Dirichlet=2, Dirichlet mid=3;
A = spdiags([-1 1 0;ones(n-2,1)*[-1 2 -1];0 a11 -1],-1:1,n,n)'/h^2;
function B = avg(A,k)
if nargin<2, k = 1; end
if size(A,1)==1, A = A'; end
if k<2, B = conv2(A,[1;1]/2,'valid'); else B = avg(A,k-1); end
if size(A,2)==1, B = B'; end
|
