https://cfd-online.com/W/index.php?title=Pfv_lid_driven_cavity&feed=atom&action=history
Pfv lid driven cavity - Revision history
2024-03-29T00:09:36Z
Revision history for this page on the wiki
MediaWiki 1.16.5
https://cfd-online.com/W/index.php?title=Pfv_lid_driven_cavity&diff=13018&oldid=prev
Jonas Holdeman: Created page with "Matlab script '''LDCW.m''' for lid-driven cavity. <pre> %LDCW LID-DRIVEN CAVITY % Finite element solution of the 2D Navier-Stokes equation using 4-node, 12 DOF, % ..."
2011-06-29T03:26:01Z
<p>Created page with "Matlab script '''LDCW.m''' for lid-driven cavity. <pre> %LDCW LID-DRIVEN CAVITY % Finite element solution of the 2D Navier-Stokes equation using 4-node, 12 DOF, % ..."</p>
<p><b>New page</b></p><div>Matlab script '''LDCW.m''' for lid-driven cavity. <br />
<br />
<pre><br />
%LDCW LID-DRIVEN CAVITY <br />
% Finite element solution of the 2D Navier-Stokes equation using 4-node, 12 DOF,<br />
% (3-DOF/node), simple-cubic-derived rectangular Hermite basis for <br />
% the Lid-Driven Cavity problem.<br />
%<br />
% This could also be characterized as a VELOCITY-STREAM FUNCTION or <br />
% STREAM FUNCTION-VELOCITY method.<br />
%<br />
% Reference: "A Hermite finite element method for incompressible fluid flow", <br />
% Int. J. Numer. Meth. Fluids, 64, P376-408 (2010). <br />
%<br />
% Simplified Wiki version <br />
% The rectangular problem domain is defined between Cartesian <br />
% coordinates Xmin & Xmax and Ymin & Ymax.<br />
% The computational grid has NumEx elements in the x-direction <br />
% and NumEy elements in the y-direction. <br />
% The nodes and elements are numbered column-wise from the <br />
% upper left corner to the lower right corner. <br />
%<br />
%This script calls the user-defined functions:<br />
% regrade - to regrade the mesh <br />
% DMatW - to evaluate element diffusion matrix <br />
% CMatW - to evaluate element convection matrix<br />
% GetPresW - to evaluate the pressure <br />
% ilu_gmres_with_EBC - to solve the system with essential/Dirichlet BCs <br />
%<br />
% Jonas Holdeman August 2007, revised June 2011<br />
<br />
clear all;<br />
disp('Lid-driven cavity');<br />
disp(' Four-node, 12 DOF, simple-cubic stream function basis.');<br />
<br />
% -------------------------------------------------------------<br />
nd = 3; nd2=nd*nd; % Number of DOF per node - do not change!!<br />
% -------------------------------------------------------------<br />
ETstart=clock;<br />
<br />
% Parameters for GMRES solver <br />
GMRES.Tolerance=1.e-14;<br />
GMRES.MaxIterates=20; <br />
GMRES.MaxRestarts=6;<br />
<br />
% Optimal relaxation parameters for given Reynolds number<br />
% (see IJNMF reference)<br />
% Re 100 1000 3200 5000 7500 10000 12500 <br />
% RelxFac: 1.04 1.11 .860 .830 .780 .778 .730 <br />
% ExpCR1 1.488 .524 .192 .0378 -- -- -- <br />
% ExpCRO 1.624 .596 .390 .331 .243 .163 .133<br />
% CritFac: 1.82 1.49 1.14 1.027 .942 .877 .804 <br />
<br />
% Define the problem geometry, set mesh bounds:<br />
Xmin = 0.0; Xmax = 1.0; Ymin = 0.0; Ymax = 1.0; <br />
<br />
% Set mesh grading parameters (set to 1 if no grading).<br />
% See below for explanation of use of parameters. <br />
xgrd = .75; ygrd=.75; % (xgrd = 1, ygrd=1 for uniform mesh) <br />
<br />
% Set " RefineBoundary=1 " for additional refinement at boundary, <br />
% i.e., split first element along boundary into two. <br />
RefineBoundary=1; <br />
<br />
% DEFINE THE MESH <br />
% Set number of elements in each direction<br />
NumEx = 16; NumEy = NumEx;<br />
<br />
% PLEASE CHANGE OR SET NUMBER OF ELEMENTS TO CHANGE/SET NUMBER OF NODES!<br />
NumNx=NumEx+1; NumNy=NumEy+1;<br />
<br />
% Define problem parameters: <br />
% Lid velocity<br />
Vlid=1.;<br />
<br />
% Reynolds number<br />
Re=1000.; <br />
<br />
% factor for under/over-relaxation starting at iteration RelxStrt <br />
RelxFac = 1.; % <br />
<br />
% Number of nonlinear iterations<br />
MaxNLit=10; %<br />
<br />
%--------------------------------------------------------<br />
<br />
% Viscosity for specified Reynolds number<br />
nu=Vlid*(Xmax-Xmin)/Re; <br />
<br />
% Grade the mesh spacing if desired, call regrade(x,agrd,e). <br />
% if e=0: refine both sides, 1: refine upper, 2: refine lower<br />
% if agrd=xgrd|ygrd is the parameter which controls grading, then<br />
% if agrd=1 then leave array unaltered.<br />
% if agrd<1 then refine (make finer) towards the ends<br />
% if agrd>1 then refine (make finer) towards the center.<br />
% <br />
% Generate equally-spaced nodal coordinates and refine if desired.<br />
if (RefineBoundary==1)<br />
XNc=linspace(Xmin,Xmax,NumNx-2); <br />
XNc=[XNc(1),(.62*XNc(1)+.38*XNc(2)),XNc(2:end-1),(.38*XNc(end-1)+.62*XNc(end)),XNc(end)];<br />
YNc=linspace(Ymax,Ymin,NumNy-2); <br />
YNc=[YNc(1),(.62*YNc(1)+.38*YNc(2)),YNc(2:end-1),(.38*YNc(end-1)+.62*YNc(end)),YNc(end)];<br />
else<br />
XNc=linspace(Xmin,Xmax,NumNx); <br />
YNc=linspace(Ymax,Ymin,NumNy); <br />
end<br />
if xgrd ~= 1 XNc=regrade(XNc,xgrd,0); end; % Refine mesh if desired<br />
if ygrd ~= 1 YNc=regrade(YNc,ygrd,0); end;<br />
[Xgrid,Ygrid]=meshgrid(XNc,YNc);% Generate the x- and y-coordinate meshes.<br />
<br />
% Allocate storage for fields <br />
psi0=zeros(NumNy,NumNx);<br />
u0=zeros(NumNy,NumNx);<br />
v0=zeros(NumNy,NumNx);<br />
<br />
%--------------------Begin grid plot-----------------------<br />
% ********************** FIGURE 1 *************************<br />
% Plot the grid <br />
figure(1);<br />
clf;<br />
orient portrait; orient tall; <br />
subplot(2,2,1);<br />
hold on;<br />
plot([Xmax;Xmin],[YNc;YNc],'k');<br />
plot([XNc;XNc],[Ymax;Ymin],'k');<br />
hold off;<br />
axis([Xmin,Xmax,Ymin,Ymax]);<br />
axis equal;<br />
axis image;<br />
title([num2str(NumNx) 'x' num2str(NumNy) ...<br />
' node mesh for Lid-driven cavity']);<br />
pause(.1);<br />
%-------------- End plotting Figure 1 ----------------------<br />
<br />
<br />
%Contour levels, Ghia, Ghia & Shin, Re=100, 400, 1000, 3200, ...<br />
clGGS=[-.1175,-.1150,-.11,-.1,-.09,-.07,-.05,-.03,-.01,-1.e-4,-1.e-5,-1.e-7,-1.e-10,...<br />
1.e-8,1.e-7,1.e-6,1.e-5,5.e-5,1.e-4,2.5e-4,5.e-4,1.e-3,1.5e-3,3.e-3];<br />
CL=clGGS; % Select contour level option<br />
if (Vlid<0) CL=-CL; end<br />
<br />
NumNod=NumNx*NumNy; % total number of nodes<br />
MaxDof=nd*NumNod; % maximum number of degrees of freedom<br />
EBC.Mxdof=nd*NumNod; % maximum number of degrees of freedom<br />
<br />
nn2nft=zeros(2,NumNod); % node number -> nf & nt<br />
NodNdx=zeros(2,NumNod);<br />
% Generate lists of active nodal indices, freedom number & type <br />
ni=0; nf=-nd+1; nt=1; % ________<br />
for nx=1:NumNx % | |<br />
for ny=1:NumNy % | |<br />
ni=ni+1; % |________|<br />
NodNdx(:,ni)=[nx;ny];<br />
nf=nf+nd; % all nodes have 4 dofs <br />
nn2nft(:,ni)=[nf;nt]; % dof number & type (all nodes type 1)<br />
end;<br />
end;<br />
%NumNod=ni; % total number of nodes<br />
nf2nnt=zeros(2,MaxDof); % (node, type) associated with dof<br />
ndof=0; dd=[1:nd];<br />
for n=1:NumNod<br />
for k=1:nd<br />
nf2nnt(:,ndof+k)=[n;k];<br />
end<br />
ndof=ndof+nd;<br />
end<br />
<br />
NumEl=NumEx*NumEy;<br />
<br />
% Generate element connectivity, from upper left to lower right. <br />
Elcon=zeros(4,NumEl);<br />
ne=0; LY=NumNy;<br />
for nx=1:NumEx<br />
for ny=1:NumEy<br />
ne=ne+1;<br />
Elcon(1,ne)=1+ny+(nx-1)*LY; <br />
Elcon(2,ne)=1+ny+nx*LY;<br />
Elcon(3,ne)=1+(ny-1)+nx*LY;<br />
Elcon(4,ne)=1+(ny-1)+(nx-1)*LY;<br />
end % loop on ny<br />
end % loop on nx<br />
<br />
% Begin essential boundary conditions, allocate space <br />
MaxEBC = nd*2*(NumNx+NumNy-2);<br />
EBC.dof=zeros(MaxEBC,1); % Degree-of-freedom index <br />
EBC.typ=zeros(MaxEBC,1); % Dof type (1,2,3)<br />
EBC.val=zeros(MaxEBC,1); % Dof value <br />
<br />
X1=XNc(2); X2=XNc(NumNx-1);<br />
nc=0;<br />
for nf=1:MaxDof<br />
ni=nf2nnt(1,nf);<br />
nx=NodNdx(1,ni);<br />
ny=NodNdx(2,ni);<br />
x=XNc(nx);<br />
y=YNc(ny); <br />
if(x==Xmin | x==Xmax | y==Ymin)<br />
nt=nf2nnt(2,nf);<br />
switch nt;<br />
case {1, 2, 3}<br />
nc=nc+1; EBC.typ(nc)=nt; EBC.dof(nc)=nf; EBC.val(nc)=0; % psi, u, v <br />
end % switch (type)<br />
elseif (y==Ymax)<br />
nt=nf2nnt(2,nf);<br />
switch nt;<br />
case {1, 3}<br />
nc=nc+1; EBC.typ(nc)=nt; EBC.dof(nc)=nf; EBC.val(nc)=0; % psi, v <br />
case 2<br />
nc=nc+1; EBC.typ(nc)=nt; EBC.dof(nc)=nf; EBC.val(nc)=Vlid; % u<br />
end % switch (type) <br />
end % if (boundary)<br />
end % for nf <br />
EBC.num=nc; <br />
<br />
if (size(EBC.typ,1)>nc) % Truncate arrays if necessary <br />
EBC.typ=EBC.typ(1:nc);<br />
EBC.dof=EBC.dof(1:nc);<br />
EBC.val=EBC.val(1:nc);<br />
end % End ESSENTIAL (Dirichlet) boundary conditions<br />
<br />
% partion out essential (Dirichlet) dofs<br />
p_vec = [1:EBC.Mxdof]'; % List of all dofs<br />
EBC.p_vec_undo = zeros(1,EBC.Mxdof);<br />
% form a list of non-diri dofs<br />
EBC.ndro = p_vec(~ismember(p_vec, EBC.dof)); % list of non-diri dofs<br />
% calculate p_vec_undo to restore Q to the original dof ordering<br />
EBC.p_vec_undo([EBC.ndro;EBC.dof]) = [1:EBC.Mxdof]; %p_vec';<br />
<br />
Q=zeros(MaxDof,1); % Allocate space for solution (dof) vector<br />
<br />
% Initialize fields to boundary conditions<br />
for k=1:EBC.num<br />
Q(EBC.dof(k))=EBC.val(k); <br />
end;<br />
<br />
errpsi=zeros(NumNy,NumNx); % error correct for iteration<br />
<br />
MxNL=max(1,MaxNLit);<br />
np0=zeros(1,MxNL); % Arrays for convergence info<br />
nv0=zeros(1,MxNL);<br />
<br />
Qs=[];<br />
<br />
Dmat = spalloc(MaxDof,MaxDof,36*MaxDof); % to save the diffusion matrix<br />
Vdof=zeros(nd,4);<br />
Xe=zeros(2,4); % coordinates of element corners <br />
<br />
NLitr=0; ND=1:nd;<br />
while (NLitr<MaxNLit), NLitr=NLitr+1; % <<< BEGIN NONLINEAR ITERATION <br />
<br />
tclock=clock; % Start assembly time <<<<<<<<<<br />
% Generate and assemble element matrices<br />
Mat=spalloc(MaxDof,MaxDof,36*MaxDof);<br />
RHS=spalloc(MaxDof,1,MaxDof);<br />
%RHS = zeros(MaxDof,1);<br />
Emat=zeros(1,16*nd2); % Values 144=4*4*(nd*nd) <br />
<br />
% BEGIN GLOBAL MATRIX ASSEMBLY<br />
for ne=1:NumEl <br />
<br />
for k=1:4<br />
ki=NodNdx(:,Elcon(k,ne));<br />
Xe(:,k)=[XNc(ki(1));YNc(ki(2))]; <br />
end<br />
<br />
if NLitr == 1 <br />
% Fluid element diffusion matrix, save on first iteration <br />
[DEmat,Rndx,Cndx] = DMatW(Xe,Elcon(:,ne),nn2nft);<br />
Dmat=Dmat+sparse(Rndx,Cndx,DEmat,MaxDof,MaxDof); % Global diffusion mat <br />
end <br />
<br />
if (NLitr>1) <br />
% Get stream function and velocities<br />
for n=1:4 <br />
Vdof(ND,n)=Q((nn2nft(1,Elcon(n,ne))-1)+ND); % Loop over local element nodes<br />
end<br />
% Fluid element convection matrix, first iteration uses Stokes equation. <br />
[Emat,Rndx,Cndx] = CMatW(Xe,Elcon(:,ne),nn2nft,Vdof); <br />
Mat=Mat+sparse(Rndx,Cndx,-Emat,MaxDof,MaxDof); % Global convection assembly <br />
end<br />
<br />
end; % loop ne over elements <br />
% END GLOBAL MATRIX ASSEMBLY<br />
<br />
Mat = Mat -nu*Dmat; % Add in cached/saved global diffusion matrix <br />
<br />
disp(['(' num2str(NLitr) ') Matrix assembly complete, elapsed time = '...<br />
num2str(etime(clock,tclock)) ' sec']); % Assembly time <<<<<<<<<<<<br />
pause(1);<br />
<br />
Q0 = Q; % Save dof values <br />
<br />
% Solve system<br />
tclock=clock; %disp('start solution'); % Start solution time <<<<<<<<<<<<<<<br />
<br />
RHSr=RHS(EBC.ndro)-Mat(EBC.ndro,EBC.dof)*EBC.val;<br />
Matr=Mat(EBC.ndro,EBC.ndro);<br />
Qs=Q(EBC.ndro);<br />
<br />
Qr=ilu_gmres_with_EBC(Matr,RHSr,[],GMRES,Qs);<br />
<br />
Q=[Qr;EBC.val]; % Augment active dofs with esential (Dirichlet) dofs<br />
Q=Q(EBC.p_vec_undo); % Restore natural order<br />
<br />
stime=etime(clock,tclock); % Solution time <<<<<<<<<<<<<<<br />
<br />
% ****** APPLY RELAXATION FACTOR *********************<br />
if(NLitr>1) Q=RelxFac*Q+(1-RelxFac)*Q0; end<br />
% ****************************************************<br />
<br />
% Compute change and copy dofs to field arrays<br />
dsqp=0; dsqv=0;<br />
for k=1:MaxDof<br />
ni=nf2nnt(1,k); nx=NodNdx(1,ni); ny=NodNdx(2,ni);<br />
switch nf2nnt(2,k) % switch on dof type <br />
case 1<br />
dsqp=dsqp+(Q(k)-Q0(k))^2; psi0(ny,nx)=Q(k);<br />
errpsi(ny,nx)=Q0(k)-Q(k); <br />
case 2<br />
dsqv=dsqv+(Q(k)-Q0(k))^2; u0(ny,nx)=Q(k);<br />
case 3<br />
dsqv=dsqv+(Q(k)-Q0(k))^2; v0(ny,nx)=Q(k);<br />
end % switch on dof type <br />
end % for <br />
np0(NLitr)=sqrt(dsqp); <br />
nv0(NLitr)=sqrt(dsqv); <br />
<br />
if (np0(NLitr)<=1e-15|nv0(NLitr)<=1e-15) <br />
MaxNLit=NLitr; np0=np0(1:MaxNLit); nv0=nv0(1:MaxNLit); end;<br />
disp(['(' num2str(NLitr) ') Solution time for linear system = '...<br />
num2str(etime(clock,tclock)) ' sec']); % Solution time <<<<<<<<<<<<<br />
<br />
%---------- Begin plot of intermediate results ----------<br />
% ********************** FIGURE 2 *************************<br />
figure(1);<br />
<br />
% Stream function (intermediate) <br />
subplot(2,2,3);<br />
contour(Xgrid,Ygrid,psi0,8,'k'); % Plot contours (trajectories)<br />
axis([Xmin,Xmax,Ymin,Ymax]);<br />
title(['Lid-driven cavity, Re=' num2str(Re)]);<br />
axis equal; axis image;<br />
<br />
% Plot convergence info <br />
subplot(2,2,2);<br />
semilogy(1:NLitr,nv0(1:NLitr),'k-+',1:NLitr,np0(1:NLitr),'k-o');<br />
xlabel('Nonlinear iteration number');<br />
ylabel('Nonlinear correction');<br />
axis square; <br />
title(['Iteration conv., Re=' num2str(Re)]);<br />
legend('U','Psi');<br />
<br />
% Plot nonlinear iteration correction contours <br />
subplot(2,2,4);<br />
contour(Xgrid,Ygrid,errpsi,8,'k'); % Plot contours (trajectories)<br />
axis([Xmin,Xmax,Ymin,Ymax]);<br />
axis equal; axis image;<br />
title(['Iteration correction']);<br />
pause(1);<br />
% ********************** END FIGURE 2 *************************<br />
%---------- End plot of intermediate results ---------<br />
<br />
if (nv0(NLitr)<1e-15) break; end % Terminate iteration if non-significant <br />
<br />
end; % <<< (while) END NONLINEAR ITERATION<br />
<br />
format short g;<br />
disp('Convergence results by iteration: velocity, stream function');<br />
disp(['nv0: ' num2str(nv0)]); disp(['np0: ' num2str(np0)]); <br />
<br />
% >>>>>>>>>>>>>> BEGIN PRESSURE RECOVERY <<<<<<<<<<<<<<<<<<<br />
% Essential pressure boundary condition <br />
% Index of node to apply pressure BC, value at node<br />
PBCnx=fix((NumNx+1)/2); % Apply at center of mesh<br />
PBCny=fix((NumNy+1)/2);<br />
PBCnod=0;<br />
for k=1:NumNod<br />
if (NodNdx(1,k)==PBCnx & NodNdx(2,k)==PBCny) PBCnod=k; break; end<br />
end<br />
if (PBCnod==0) error('Pressure BC node not found'); <br />
else<br />
EBCp.nodn = [PBCnod]; % Pressure BC node number<br />
EBCp.val = [0]; % set P = 0.<br />
end<br />
% Cubic pressure <br />
[P,Px,Py] = GetPresW(NumNod,NodNdx,Elcon,nn2nft,Xgrid,Ygrid,Q,EBCp,nu);<br />
% ******************** END PRESSURE RECOVERY *********************<br />
<br />
% ********************** CONTINUE FIGURE 1 *************************<br />
figure(1);<br />
<br />
% Stream function (final)<br />
subplot(2,2,3);<br />
[CT,hn]=contour(Xgrid,Ygrid,psi0,CL,'k'); % Plot contours (trajectories)<br />
clabel(CT,hn,CL([1,3,5,7,9,10,11,19,23]));<br />
hold on;<br />
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');<br />
hold off;<br />
axis([Xmin,Xmax,Ymin,Ymax]);<br />
axis equal; axis image;<br />
title(['Stream lines, ' num2str(NumNx) 'x' num2str(NumNy) ...<br />
' mesh, Re=' num2str(Re)]);<br />
<br />
% Plot pressure contours (final)<br />
subplot(2,2,4);<br />
CPL=[-.002,0,.02,.05,.07,.09,.11,.12,.17,.3];<br />
[CT,hn]=contour(Xgrid,Ygrid,P,CPL,'k'); % Plot pressure contours<br />
clabel(CT,hn,CPL([3,5,7,10]));<br />
hold on;<br />
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');<br />
hold off;<br />
axis([Xmin,Xmax,Ymin,Ymax]);<br />
axis equal; axis image;<br />
title(['Simple cubic pressure contours, Re=' num2str(Re)]);<br />
% ********************* END FIGURE 1 *************************<br />
<br />
disp(['Total elapsed time = '...<br />
num2str(etime(clock,ETstart)/60) ' min']); % Elapsed time from start <<<<br />
</pre></div>
Jonas Holdeman